text:
"""
    IndicatorSOS1Bridge{T, BC <: MOI.AbstractScalarSet}

The `IndicatorSOS1Bridge` replaces an indicator constraint of the following form:
``z \\in \\mathbb{B}, z == 1 \\implies f(x) \\leq b`` with a SOS1 constraint:
``z \\in \\mathbb{B}, w \\leq 0, f(x) + w \\leq b, SOS1(w, z)``.
`GreaterThan` constraints are handled in a symmetric way:
``z \\in \\mathbb{B}, z == 1 \\implies f(x) \\geq b`` is reformulated as:
``z \\in \\mathbb{B}, w \\geq 0, f(x) + w \\geq b, SOS1(w, z)``.
Other scalar sets are handled without a bound constraint:
``z \\in \\mathbb{B}, z == 1 \\implies f(x) == b`` is reformulated as:
``z \\in \\mathbb{B}, w \\text{ free}, f(x) + w == b, SOS1(w, z)``.

If `BC !<: Union{LessThan, GreaterThan}`, `bound_constraint_index` is `nothing`.
"""
struct IndicatorSOS1Bridge{T, BC <: MOI.AbstractScalarSet, MaybeBC <: Union{MOI.ConstraintIndex{MOI.SingleVariable, BC}, Nothing}} <: AbstractBridge
    w_variable_index::MOI.VariableIndex
    z_variable_index::MOI.VariableIndex
    affine_func::MOI.ScalarAffineFunction{T}
    bound_constraint_index::MaybeBC
    sos_constraint_index::MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.SOS1{T}}
    linear_constraint_index::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, BC}
end

function bridge_constraint(::Type{IndicatorSOS1Bridge{T,BC,MaybeBC}}, model::MOI.ModelLike, f::MOI.VectorAffineFunction{T}, s::MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE, BC}) where {T <: Real, BC, MaybeBC}
    f_scalars = MOIU.eachscalar(f)
    (w, bound_constraint) = _add_bound_constraint!(model, BC)
    z = convert(MOI.SingleVariable, f_scalars[1]).variable
    sos_vector = MOI.VectorOfVariables([w, z])
    sos_constraint = MOI.add_constraint(model, sos_vector, MOI.SOS1{T}([0.4, 0.6]))
    affine_func = f_scalars[2]
    affine_expr = MOIU.operate(+, T, affine_func, MOI.SingleVariable(w))
    linear_constraint = MOI.add_constraint(model, affine_expr, s.set)
    return IndicatorSOS1Bridge{T,BC,MaybeBC}(w, z, affine_func, bound_constraint, sos_constraint, linear_constraint)
end

function _add_bound_constraint!(model::MOI.ModelLike, ::Type{BC}) where {T <: Real, BC <: Union{MOI.LessThan{T}, MOI.GreaterThan{T}}}
    return MOI.add_constrained_variable(model, BC(zero(T)))
end

function _add_bound_constraint!(model::MOI.ModelLike, ::Type{<:MOI.AbstractScalarSet})
    return (MOI.add_variable(model), nothing)
end

function MOI.supports_constraint(::Type{<:IndicatorSOS1Bridge},
                                 ::Type{<:MOI.AbstractVectorFunction},
                                 ::Type{<:MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE, <:MOI.AbstractScalarSet}})
    return true
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintSet,
                 b::IndicatorSOS1Bridge)
    return MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(
        MOI.get(model, attr, b.linear_constraint_index),
    )
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 b::IndicatorSOS1Bridge{T}) where {T}
    z = b.z_variable_index
    terms = [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(one(T), z))]
    for affine_term in b.affine_func.terms
        push!(terms, MOI.VectorAffineTerm(2, affine_term))
    end
    return MOI.VectorAffineFunction(
        terms, [zero(T), b.affine_func.constant]
    )
end

function MOI.delete(model::MOI.ModelLike, bridge::IndicatorSOS1Bridge)
    if bridge.bound_constraint_index !== nothing
        MOI.delete(model, bridge.bound_constraint_index)
    end
    MOI.delete(model, bridge.sos_constraint_index)
    MOI.delete(model, bridge.linear_constraint_index)
    MOI.delete(model, bridge.w_variable_index)
    return
end

function MOIB.added_constrained_variable_types(::Type{<:IndicatorSOS1Bridge{T, BC}}) where {T, BC <: Union{MOI.LessThan{T}, MOI.GreaterThan{T}}}
    return [(BC,)]
end

function MOIB.added_constrained_variable_types(::Type{<:IndicatorSOS1Bridge{T, BC}}) where {T, BC}
    return []
end

function MOIB.added_constraint_types(::Type{<:IndicatorSOS1Bridge{T, BC}}) where {T, BC <: Union{MOI.LessThan{T}, MOI.GreaterThan{T}}}
    return [(MOI.VectorOfVariables, MOI.SOS1{T}),
            (MOI.ScalarAffineFunction{T}, BC),
            ]
end

function MOIB.added_constraint_types(::Type{<:IndicatorSOS1Bridge{T, S}}) where {T, S <: MOI.AbstractScalarSet}
    return [(MOI.VectorOfVariables, MOI.SOS1{T}),
            (MOI.ScalarAffineFunction{T}, S),
           ]
end

function concrete_bridge_type(::Type{<:IndicatorSOS1Bridge{T}},
                              ::Type{<:MOI.AbstractVectorFunction},
                              ::Type{MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE, S}}) where {T, S<:Union{MOI.LessThan, MOI.GreaterThan}}
    return IndicatorSOS1Bridge{T, S, MOI.ConstraintIndex{MOI.SingleVariable, S}}
end

function concrete_bridge_type(::Type{<:IndicatorSOS1Bridge{T}},
                              ::Type{<:MOI.AbstractVectorFunction},
                              ::Type{MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE, S}}) where {T, S <: MOI.AbstractScalarSet}
    return IndicatorSOS1Bridge{T, S, Nothing}
end

# Attributes, Bridge acting as a model

function MOI.get(::IndicatorSOS1Bridge, ::MOI.NumberOfVariables)
    return 1
end

function MOI.get(b::IndicatorSOS1Bridge, ::MOI.ListOfVariableIndices)
    return [b.w_variable_index]
end

function MOI.get(::IndicatorSOS1Bridge{T, BC, Nothing}, ::MOI.NumberOfConstraints{MOI.SingleVariable, BC}) where {T, BC}
    return 0
end

function MOI.get(::IndicatorSOS1Bridge{T, BC, CI}, ::MOI.NumberOfConstraints{MOI.SingleVariable, BC}) where {T, BC, CI <: MOI.ConstraintIndex{MOI.SingleVariable, BC}}
    return 1
end

function MOI.get(::IndicatorSOS1Bridge, ::MOI.NumberOfConstraints{MOI.VectorOfVariables, <:MOI.SOS1})
    return 1
end

function MOI.get(::IndicatorSOS1Bridge{T, BC}, ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, BC}) where {T, BC, CI <: MOI.ConstraintIndex{MOI.SingleVariable, BC}}
    return 1
end

function MOI.get(b::IndicatorSOS1Bridge{T, BC, CI}, ::MOI.ListOfConstraintIndices{MOI.SingleVariable, BC}) where {T, BC, CI <: MOI.ConstraintIndex}
    return [b.bound_constraint_index]
end

function MOI.get(::IndicatorSOS1Bridge{T, BC, Nothing}, ::MOI.ListOfConstraintIndices{MOI.SingleVariable, BC}) where {T, BC}
    return MOI.ConstraintIndex{MOI.SingleVariable, BC}[]
end

function MOI.get(b::IndicatorSOS1Bridge{T}, ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables, <:MOI.SOS1}) where {T}
    return [b.sos_constraint_index]
end

function MOI.get(b::IndicatorSOS1Bridge{T, BC}, ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, BC}) where {T, BC}
    return [b.linear_constraint_index]
end

function MOI.supports(
    ::MOI.ModelLike,
    ::MOI.ConstraintPrimalStart,
    ::Type{<:IndicatorSOS1Bridge})
    return true
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal, bridge::IndicatorSOS1Bridge)
    zvalue = MOI.get(model, MOI.VariablePrimal(attr.N), bridge.z_variable_index)
    wvalue = MOI.get(model, MOI.VariablePrimal(attr.N), bridge.w_variable_index)
    lin_primal_start = MOI.get(model, attr, bridge.linear_constraint_index)
    return [zvalue, lin_primal_start - wvalue]
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimalStart, bridge::IndicatorSOS1Bridge)
    zstart = MOI.get(model, MOI.VariablePrimalStart(), bridge.z_variable_index)
    wstart = MOI.get(model, MOI.VariablePrimalStart(), bridge.w_variable_index)
    lin_primal_start = MOI.get(model, attr, bridge.linear_constraint_index)
    return [zstart, lin_primal_start - wstart]
end

function MOI.set(model::MOI.ModelLike, attr::MOI.ConstraintPrimalStart,
                 bridge::IndicatorSOS1Bridge{T}, value) where {T}
    zvalue = value[1]
    lin_start = value[2]
    MOI.set(model, MOI.VariablePrimalStart(), bridge.z_variable_index, zvalue)
    wstart = MOI.get(model, MOI.VariablePrimalStart(), bridge.w_variable_index)
    wstart = wstart === nothing ? zero(T) : wstart
    MOI.set(model, attr, bridge.linear_constraint_index, lin_start + wstart)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    LearnedModel(𝓂, θ)

An object that stores a learning model `m`
along with its learned parameters `θ`.
"""
struct LearnedModel{ℳ,Θ}
  𝓂::ℳ
  θ::Θ
end

"""
    learn(𝒯, 𝒟, 𝓂)

Learn the task `𝒯` with geospatial data `𝒟`
using a learning model `𝓂`.
"""
function learn(𝒯::LearningTask, 𝒟, 𝓂)
  # retrieve table of values
  table = values(𝒟)

  # learn model with table
  if issupervised(𝒯)
    X = TO.select(table, features(𝒯)...)
    y = Tables.getcolumn(table, label(𝒯))
    θ, _, __ = MI.fit(𝓂, 0, X, y)
  else
    X = TO.select(table, features(𝒯)...)
    θ, _, __ = MI.fit(𝓂, 0, X)
  end

  # return learned model
  LearnedModel(𝓂, θ)
end

"""
    perform(𝒯, 𝒟, 𝓂̂)

Perform the task `𝒯` with geospatial data `𝒟` using
a *learned* model `𝓂̂`.
"""
function perform(𝒯::LearningTask, 𝒟, 𝓂̂)
  # unpack model and learned parameters
  𝓂, θ = 𝓂̂.𝓂, 𝓂̂.θ

  # retrieve table of values
  table = values(𝒟)

  # apply model to the data
  X = TO.select(table, features(𝒯)...)
  ŷ = MI.predict(𝓂, θ, X)

  # post-process result
  var = outputvars(𝒯)[1]
  val = if issupervised(𝒯)
    isprobabilistic(𝓂) ? mode.(ŷ) : ŷ
  else
    ŷ
  end

  ctor = constructor(typeof(𝒟))
  dom  = domain(𝒟)
  tab  = (; var=>val)
  dat  = Dict(paramdim(dom) => tab)

  ctor(dom, dat)
end

ext: jl


text:
using ProgressMeter

"""
    +(self, table[, N])

Computes with N (default 1,000) Monte Carlo simulations the odds of the hand
winning the showdown given the cards shown on the table

"""
function calculate_odds(self::Hand, table::Hand; N = 1_000)
    wins = 0
    for i = 1:N
        deck = fresh_deck - self - table
        opp, deck = deal(deck, 2)
        extra_cards, deck = deal(deck, 5 - length(table))
        if opp + table + extra_cards < self + table + extra_cards
            wins += 1
        end
    end
    return wins / N
end

"""
    get_preflop_odds([N])

Computes with N (default 1,000) Monte Carlo simulations the odds of the each
starting hand winning the showdown given no cards shown on the table yet.
This table is optionally used to speed up game calculations.

"""
function get_preflop_odds(; N = 1000)
    preflop_odds = zeros(52, 52)
    @showprogress for i = 1:52, j = 1:52
        if i != j
            preflop_odds[i, j] = calculate_odds(
                Hand([fresh_deck_cards[i], fresh_deck_cards[j]]),
                Hand([]), N = N)
        end
    end
    return preflop_odds
end


text:
module ComoniconTargetExpr

using ExproniconLite
using ComoniconTypes
using ComoniconOptions

export emit_expr

include("emit.jl")

end


text:
#-------------------------------------------------------------------------------
# TSeries struct and constructors
#-------------------------------------------------------------------------------

"""
    TSeries

Data structure representing a time-series vector. The following
operations are allowed:

 - indexing using `MIT` (aka "moment-in-time") and `UnitRange{MIT}`
 - assignment using `MIT` and `UnitRange{MIT}`

In addition, most of the operations available to Julia vectors (+, -, *, etc.)
are supported by `TSeries` as well.

### Examples
 - Create `TSeries`
```julia-repl
julia> x = TSeries(qq(2020, 1), ones(4))
TSeries{Quarterly} of length 4
2020Q1: 1.0
2020Q2: 1.0
2020Q3: 1.0
2020Q4: 1.0
```

 - Index into `TSeries`
```julia-repl
julia> x[2000Q1]
1.0

julia> x[qq(2020, 1):qq(2020, 2)]
TSeries{Quarterly} of length 2
2020Q1: 1.0
2020Q2: 1.0
```

- Assignment using `MIT`
```julia-repl
julia> x[qq(2020, 1)] = 100; x
TSeries{Quarterly} of length 4
2020Q1: 100.0
2020Q2: 1.0
2020Q3: 1.0
2020Q4: 1.0

julia> x[qq(2020, 1):qq(2020, 2)] = 100; x
TSeries{Quarterly} of length 4
2020Q1: 100.0
2020Q2: 100.0
2020Q3: 1.0
2020Q4: 1.0
```

 - Arithmetic Operations on `TSeries`
```julia-repl
julia> x = TSeries(qq(2020, 1), ones(4))
TSeries{Quarterly} of length 4
2020Q1: 1.0
2020Q2: 1.0
2020Q3: 1.0
2020Q4: 1.0

julia> 2*x + 98
TSeries{Quarterly} of length 4
2020Q1: 100.0
2020Q2: 100.0
2020Q3: 100.0
2020Q4: 100.0

julia> log(exp(x))
TSeries{Quarterly} of length 4
2020Q1: 1.0
2020Q2: 1.0
2020Q3: 1.0
2020Q4: 1.0
```
"""
mutable struct TSeries{T <: Frequency, C <: AbstractVector{Float64}} <: AbstractVector{Float64}
    firstdate::MIT{T}
    values::C
end

# We work only with Float64. All other numbers are converted to it.
# Might be inefficient
TSeries(fd::MIT, v::AbstractVector{<:Number}) = TSeries(fd, Vector{Float64}(v))

TSeries(v::AbstractVector{<:Number}) = TSeries(1U, v)

TSeries(I::AbstractUnitRange{<:MIT}) = TSeries(first(I), Vector{Float64}(undef, length(I)))
TSeries(I::AbstractUnitRange{<:MIT}, ::UndefInitializer) = TSeries(I)

# TSeries constructor with a range and data
function TSeries(I::AbstractUnitRange{<:MIT}, V::AbstractVector{<:Number})
    if length(I) ≠ length(V)
        throw(ArgumentError("Range and data lengths don't match."))
    end
    TSeries(first(I), V)
end

# TSeries constructor with a range and a single value
function TSeries(I::AbstractUnitRange{<:MIT}, v::Number)
    # use the inner constructor
    TSeries(first(I), fill(Float64(v), length(I)))
end

#-------------------------------------------------------------------------------
# Base.show
#-------------------------------------------------------------------------------

Base.summary(io::IO, t::TSeries) = isempty(t) ?
        print(io, "Empty ", frequencyof(t), " TSeries") :
        print(IOContext(io, :compact=>true), length(t.values), "-element ", frequencyof(t), " TSeries from ", t.firstdate)

Base.show(io::IO, ::MIME"text/plain", t::TSeries) = show(io, t)
function Base.show(io::IO, t::TSeries)
    summary(io, t)
    isempty(t) && return
    print(io, ":")
    limit = get(io, :limit, true)
    nval = length(t.values)
    from = t.firstdate
    nrow, ncol = displaysize(io)
    if limit && nval > nrow - 5
        top = div(nrow - 5, 2)
        bot = nval - nrow + 6 + top
        for i = 1:top
            print(io, "\n", lpad(from+(i-1), 8), " : ", t.values[i])
        end
        print(io, "\n    ⋮")
        for i = bot:nval
            print(io, "\n", lpad(from+(i-1), 8), " : ", t.values[i])
        end
    else
        for i = 1:nval
            print(io, "\n", lpad(from+(i-1), 8), " : ", t.values[i])
        end
    end
end

#-------------------------------------------------------------------------------
# Base.getindex and Base.setindex
#-------------------------------------------------------------------------------

# """
# Since `TSeries` is a subtype of `AbstractVector`, we have to provide implementations for `size`, `getindex`, and `setindex!`
# """
Base.:(==)(x::TSeries{T}, y::TSeries{T}) where T <: Frequency = (x.firstdate == y.firstdate && isequal(x.values, y.values))
Base.size(ts::TSeries) = size(ts.values)

# Indexing with plain integers simply indexes within the values
Base.getindex(ts::TSeries, i::Int64) = getindex(ts.values, i)
Base.getindex(ts::TSeries, I::AbstractUnitRange{Int}) = TSeries(firstdate(ts) - 1 .+ I, getindex(ts.values, I))
Base.getindex(ts::TSeries, I::AbstractVector{Int}) = getindex(ts.values, I)
Base.setindex!(ts::TSeries, v::Number, i::Int64) = setindex!(ts.values, v, i)
Base.setindex!(ts::TSeries, v::Number, I::AbstractUnitRange{Int}) = ts.values[I] .= v
Base.setindex!(ts::TSeries, v::Number, I::AbstractVector{Int}) = ts.values[I] .= v
Base.setindex!(ts::TSeries{F}, v::TSeries{F}, I::AbstractUnitRange{Int}) where F <: Frequency = ts.values[I] .= v.values
Base.setindex!(ts::TSeries{F}, v::TSeries{F}, I::AbstractVector{Int}) where F <: Frequency = ts.values[I] .= v.values
Base.setindex!(ts::TSeries, v, I::AbstractUnitRange{Int}) = setindex!(ts.values, v, I)
Base.setindex!(ts::TSeries, v, I::AbstractVector{Int}) = setindex!(ts.values, v, I)

Base.axes(t::TSeries) = (mitrange(t),)
Base.axes1(t::TSeries) = mitrange(t)
Base.axes(r::AbstractUnitRange{<:MIT}) = (Base.OneTo(length(r)),)
Base.axes1(r::AbstractUnitRange{<:MIT}) = Base.OneTo(length(r))
Base.getindex(r::AbstractUnitRange{<:MIT}, I::AbstractUnitRange{Int}) = r[first(I)]:r[last(I)]
Base.getindex(r::AbstractUnitRange{<:MIT}, I::AbstractVector{Int}) = [r[i] for i in I]

function Base.view(t::TSeries, I::AbstractUnitRange{<:Integer})
    if !<:(eltype(I), MIT)
        I = firstdate(t) - 1 .+ I
    end
    @boundscheck  checkbounds(t, I)
    TSeries(I, @inbounds view(t.values, I .- t.firstdate .+ 1))
end

Base.similar(t::TSeries) = TSeries(firstdate(t), similar(getfield(t, :values)))
Base.dataids(t::TSeries) = Base.dataids(getfield(t, :values))
Base.IndexStyle(::TSeries) = IndexLinear()

"""
`getindex` using `MIT`
"""
Base.getindex(ts::TSeries{T}, i::MIT{T}) where T <: Frequency = begin
    # return `nothing` if accessing outside of the `ts` range
    if i < ts.firstdate || i > ts.firstdate + length(ts) - 1
        return nothing
    end

    i_int = i - ts.firstdate + 1
    ts[i_int]
end

"""
`getindex` using `Vector{MIT}`. Note the difference between `Vector{MIT}` and `UnitRange{MIT}`
"""
Base.getindex(s::TSeries{T},v::AbstractVector{MIT{T}}) where T <: Frequency = begin
    [s[i] for i in v]
end

"""
`getindex` using `UnitRange{MIT}`
"""
Base.getindex(ts::TSeries{T}, I::AbstractUnitRange{MIT{T}}) where T <: Frequency = begin

    I_int = I .- ts.firstdate .+ 1

    I_common = intersect(I_int, 1:length(ts))

    if I_common.start > I_common.stop
        println("Warning: $I is outside of TSeries bounds $(mitrange(ts)).")
        return nothing
    end

    firstdate_new = MIT{T}(Int64(ts.firstdate) + I_common.start - 1)

    return TSeries(firstdate_new, ts.values[I_common])
end

"""
`setindex` a value using `MIT`
"""
Base.setindex!(ts::TSeries{T}, v::Number, mit::MIT{T}) where T <: Frequency = begin
    # Step 1: find the `distance` between
    # - mit::MIT{T} and
    # - ts.firstdate::MIT{T}
    distance = mit - ts.firstdate + 1

    # Step 2: vectorize v::IntOrFloat to append in place to ts.values::Vector{Float64}
    # in Case 2 and 3
    v_singleton = Vector{Float64}([v])

    if 1 <= distance <= length(ts)          # Case 1: place v in an array
        ts.values[distance] = Float64(v)
    elseif distance < 1                     # Case 2: extend ts.values on the left side
        val_nan_vector = append!(v_singleton, fill(NaN, abs(distance)))
        prepend!(ts.values, val_nan_vector)
        ts.firstdate = mit
    else # length(ts) < distance            # Case 3: extend ts.values on the right side
        nan_val_vector = append!(fill(NaN, abs(distance) - length(ts) - 1), v_singleton)
        append!(ts.values, nan_val_vector)
    end
end

"""
`setindex` a value using `UnitRange{MIT}`
"""
Base.setindex!(ts::TSeries{T}, v::Number, I::AbstractUnitRange{MIT{T}}) where T <: Frequency = begin
    for i in I
        ts[i] = v
    end
end

"""
`setindex` a vector using `UnitRange{MIT}`
"""
Base.setindex!(ts::TSeries{T}, v::AbstractVector{<:Number}, I::AbstractUnitRange{MIT{T}}) where T <: Frequency = begin
    for (i, val) in zip(I, v)
        ts[i] = val
    end
end

# Base.setindex!(ts::TSeries{T}, v::Vector{Int64}, I::UnitRange{MIT{T}}) where T <: Frequency = begin
#     for (i, val) in zip(I, v)
#         ts[i] = val
#     end
# end

"""
`setindex` values from other `TSeries` using `MIT`
"""
Base.setindex!(ts::TSeries{T}, v::TSeries{T}, mit::MIT{T}) where T <: Frequency = begin
    mit in ts.firstdate:lastdate(ts) || error("Given date:$mit is not in ", ts)
    mit in v.firstdate:lastdate(v) || error("Given date:$mit is not in ", v)

    ts[mit] = v[mit].values[1]
end

"""
`setindex` values from other `TSeries` using `UnitRange{MIT}`
"""
Base.setindex!(ts::TSeries{T}, v::TSeries{T}, I::AbstractUnitRange{MIT{T}}) where T <: Frequency = begin
    commonrange = intersect(I, mitrange(v))
    ts[commonrange] = v[commonrange].values
end

# """
# `setindex` values from other `Vector` using `MIT`
# """
# Base.setindex!(ts::TSeries{T}, v::AbstractVector{<:Number}, mit::MIT{T}) where T <: Frequency = begin
#     length(v) == 1 || error(v, " can contain only one element.")
#     ts[mit] = v[mit][1]
# end

#-------------------------------------------------------------------------------
# Operations
#-------------------------------------------------------------------------------

"""
    firstdate(x::TSeries)

Return an `MIT` indicating the first date in the TSeries.

### Examples
```julia-repl
julia> firstdate(TSeries(qq(2020, 1), ones(10)))
2020Q1
```
"""
firstdate(s::TSeries) = s.firstdate

"""
    lastdate(x::TSeries)

Return an `MIT` indicating the last date in the TSeries.

### Examples
```julia-repl
julia> lastdate(TSeries(qq(2020, 1), ones(10)))
2022Q2
```
"""
lastdate(s::TSeries) = (s.firstdate + length(s) - 1)

#
# Base.range(t::TSeries) = firstdate(t):lastdate(t)

"""
    Horizonatal Concatenation of `TSeries`

### Examples
```julia-repl
julia> a = TSeries(ii(1), ones(3));
julia> b = TSeries(ii(2), ones(3));
julia> [a b]
4×2 Array{Float64,2}:
   1.0  NaN
   1.0    1.0
   1.0    1.0
 NaN      1.0

```
"""
function Base.hcat(tuple_of_ts::Vararg{TSeries{T}, N}) where T <: Frequency where N
    firstdate = [i.firstdate for i in tuple_of_ts] |> minimum
    lastdate  = [TimeSeriesEcon.lastdate(i) for i in tuple_of_ts] |> maximum

    holder = Array{Float64, 1}()

    for ts in tuple_of_ts
        for date in firstdate:lastdate
            if ts[date] == nothing
                v = [NaN]
            else
                v = ts[date]
            end

            append!(holder, v)
        end
    end

    reshaped_array = reshape(holder, (length(firstdate:lastdate), length(tuple_of_ts)))

    return  reshaped_array
end

"""
    mitrange(x::TSeries)

Return an `UnitRange{MIT{<:Frequency}}` associated with `x`.

### Examples
```julia-repl
julia> mitrange(TSeries(qq(2020, 1), ones(4)))
2020Q1:2020Q4
```
"""
mitrange(ts::TSeries) = firstdate(ts):lastdate(ts)

function Base.:(+)(x::TSeries{T}, y::TSeries{T}) where T <: Frequency
    # I = intersect(x.firstdate:lastdate(x), y.firstdate:lastdate(y))
    I = intersect(mitrange(x), mitrange(y))

    I.start <= I.stop || error("There are no dates in common, operation can't be performed.")
    TSeries(I.start, x[I].values + y[I].values)
end

function Base.:(-)(x::TSeries{T}, y::TSeries{T}) where T <: Frequency
    # I = intersect(x.firstdate:lastdate(x), y.firstdate:lastdate(y))
    I = intersect(mitrange(x), mitrange(y))
    I.start <= I.stop || error("There are no dates in common, operation can't be performed.")
    TSeries(I.start, x[I].values - y[I].values)
end

function Base.:(-)(x::Number, y::TSeries)
    TSeries(y.firstdate, Float64(x) .- y.values)
end

Base.:(-)(s::TSeries{T}) where T <: Frequency = TSeries(s.firstdate, -s.values)

Base.log(ts::TSeries{T}) where T <: Frequency = TSeries(ts.firstdate, log.(ts.values))
Base.exp(ts::TSeries{T}) where T <: Frequency = TSeries(ts.firstdate, exp.(ts.values))
Base.:(+)(ts::TSeries{T}, a::Number) where T <: Frequency = TSeries(ts.firstdate, ts.values .+ a)
Base.:(+)(a::Number, ts::TSeries{T}) where T <: Frequency = TSeries(ts.firstdate, ts.values .+ a)

Base.:(-)(ts::TSeries{T}, a::Number) where T <: Frequency = TSeries(ts.firstdate, ts.values .- a)

"""
    diff(x::TSeries, k::Int64 = -1)

Same as Iris implementation of diff
"""
function Base.diff(ts::TSeries{T}, k::Int64 = -1) where T <: Frequency
    y = deepcopy(ts);

    y.firstdate = y.firstdate - k

    return ts - y
end

"""
    3 * ts
    ts * 3 returns timeseries with every element multiplied by 3
"""
function Base.:(*)(ts::TSeries{T}, s::Number) where T <: Frequency
    TSeries(ts.firstdate, ts.values .* Float64(s))
end

function Base.:(*)(s::Number, ts::TSeries{T}) where T <: Frequency
    TSeries(ts.firstdate, ts.values .* Float64(s))
end

"""
    ts / 3 returns timeseries with every element divided by 3
"""
function Base.:(/)(ts::TSeries{T}, s::Number) where T <: Frequency
    TSeries(ts.firstdate, ts.values ./ Float64(s))
end

function Base.:(/)(s::Number, ts::TSeries{T}) where T <: Frequency
    TSeries(ts.firstdate, Float64(s) ./ ts.values)
end

function Base.:(^)(ts::TSeries{T}, s::Number) where T <: Frequency
    TSeries(ts.firstdate, ts.values .^ s)
end

"""
    shift(x::TSeries, n::Int64)

Shift dates of `x` back by `k` periods.
__Note:__ The implementation of is similar to IRIS `ts{1}`.

Examples
```julia-repl
julia> shift(TSeries(qq(2020, 1), ones(4)), 1)
TSeries{Quarterly} of length 4
2019Q4: 1.0
2020Q1: 1.0
2020Q2: 1.0
2020Q3: 1.0

julia> shift(TSeries(qq(2020, 1), ones(4)), -1)
TSeries{Quarterly} of length 4
2020Q2: 1.0
2020Q3: 1.0
2020Q4: 1.0
2021Q1: 1.0
```

See also: [`shift!`](@ref)
"""
function shift(ts::TSeries{T}, k::Int64) where T <: Frequency
    return TSeries(ts.firstdate - k, ts.values)
end

"""
    shift!(x::TSeries, n::Int64)

Shift dates of `x` back by `k` periods, in-place.
__Note:__ The implementation of is similar to IRIS `ts{1}`.

Examples
```julia-repl
julia> x = TSeries(qq(2020, 1), ones(4));

julia> shift!(x, 1);

julia> x
TSeries{Quarterly} of length 4
2019Q4: 1.0
2020Q1: 1.0
2020Q2: 1.0
2020Q3: 1.0
```

See also: [`shift`](@ref)
"""
function shift!(ts::TSeries{T}, k::Int64) where T <: Frequency
    ts.firstdate = ts.firstdate - k
    return ts
end

"""
    ppy(::MIT)
    ppy(::Type{MIT})

When applied to an [`MIT`](@ref) instance or type, return the `ppy` of its frequency.
"""
ppy(::MIT{F}) where F <: Frequency = ppy(F)
ppy(::Type{MIT{F}}) where F <: Frequency = ppy(F)

"""
    ppy(::TSeries)
    ppy(::Type{TSeries})

When applied to a [`TSeries`](@ref) instance or type, return the `ppy` of its frequency.
"""
ppy(::TSeries{T}) where T <: Frequency = ppy(T)
ppy(::Type{TSeries{T}}) where T <: Frequency = ppy(T)

function Base.:(/)(x::TSeries{T}, y::TSeries{T}) where T <: Frequency
    I = intersect(mitrange(x), mitrange(y))

    firstdate = I.start;
    values = x[I].values./y[I].values;

    TSeries(firstdate, values)
end

function Base.:(*)(x::TSeries{T}, y::TSeries{T}) where T <: Frequency
    I = intersect(mitrange(x), mitrange(y))

    firstdate = I.start;
    values = x[I].values.*y[I].values;

    TSeries(firstdate, values)
end

"""
    pct(x::TSeries, shift_value::Int64, islog::Bool)

Calculate percentage growth in `x` given a `shift_value`.

__Note:__ The implementation is similar to IRIS.

Examples
```julia-repl
julia> x = TSeries(yy(2000), Vector(1:4));

julia> pct(x, -1)
TSeries{Yearly} of length 3
2001Y: 100.0
2002Y: 50.0
2003Y: 33.33333333333333
```
See also: [`apct`](@ref)
"""
function pct(ts::TSeries, shift_value::Int64; islog::Bool = false)
    if islog
        a = exp(ts);
        b = shift(exp(ts), shift_value);
    else
        a = ts;
        b = shift(ts, shift_value);
    end

    result = ( (a - b)/b ) * 100

    TSeries(result.firstdate, result.values)
end

Base.round(ts::TSeries; digits = 2) = begin
    values = (x -> round(x, digits=digits)).(ts.values)
    TSeries(ts.firstdate, values)
end

"""
    apct(x::TSeries, islog::Bool)

Calculate annualised percent rate of change in `x`.

__Note:__ The implementation is similar to IRIS.

Examples
```julia-repl
julia> x = TSeries(qq(2018, 1), Vector(1:8));

julia> apct(x)
TSeries{Quarterly} of length 7
2018Q2: 1500.0
2018Q3: 406.25
2018Q4: 216.04938271604937
2019Q1: 144.140625
2019Q2: 107.35999999999999
2019Q3: 85.26234567901243
2019Q4: 70.59558517284461
```

See also: [`pct`](@ref)
"""
function apct(ts::TSeries, islog::Bool = false)
    if islog
        a = exp(ts);
        b = shift(exp(ts), - 1);
    else
        a = ts;
        b = shift(ts, -1);
    end

    values_change = a/b
    firstdate = values_change.firstdate

    values = (values_change.^ppy(ts) .- 1) * 100

    TSeries( (a/b).firstdate, values)
end

function Base.cumsum(s::TSeries)
    TSeries(s.firstdate, cumsum(s.values))
end

function Base.cumsum!(s::TSeries)
    s.values = cumsum(s.values)
    return s
end

"""
    leftcropnan!(x::TSeries)

Remove `NaN` values from starting at the beginning of `x`, in-place.

__Note__: an internal function.
"""
function leftcropnan!(s::TSeries)
    while isequal(s[firstdate(s)], NaN)
        popfirst!(s.values)
        s.firstdate = s.firstdate + 1
    end
    return s
end

"""
rightcropnan!(x::TSeries)

Remove `NaN` values from the end of `x`

__Note__: an internal function.
"""
function rightcropnan!(s::TSeries)
    while isequal(s[lastdate(s)], NaN)
        pop!(s.values)
    end
    return s
end

"""
    nanrm!(s::TSeries, type::Symbol)

Remove `NaN` values that are either at the beginning of the `s` and/or end of `x`.

Examples
```
julia> s = TSeries(yy(2018), [NaN, NaN, 1, 2, NaN]);

julia> nanrm!(s);

julia> s
TSeries{Yearly} of length 2
2020Y: 1.0
2021Y: 2.0
```
"""
function nanrm!(s::TSeries, type::Symbol=:both)
    if type == :left
        leftcropnan!(s)
    elseif type == :right
        rightcropnan!(s)
    elseif type == :both
        leftcropnan!(s)
        rightcropnan!(s)
    else
        error("Please select between :left, :right, or :both.")
    end
    return s
end

"""
Can be applied to a [`TSeries`](@ref) instance or a range of MIT to return its [`Frequency`](@ref).
"""
frequencyof(::TSeries{T}) where T <: Frequency = T
frequencyof(::Type{TSeries{T}}) where T <: Frequency = T
frequencyof(::AbstractUnitRange{MIT{T}}) where T <: Frequency = T
frequencyof(::Type{<:AbstractUnitRange{MIT{T}}}) where T <: Frequency = T

# -----------------------------------------------------
# Broadcasting Interface: `BroadcastStyle` and `similar`
# https://docs.julialang.org/en/v1.0.5/manual/interfaces/#man-interfaces-broadcasting-1
# -----------------------------------------------------
Base.BroadcastStyle(::Type{<:TSeries}) = Broadcast.ArrayStyle{TSeries}()

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{TSeries}}, ::Type{ElType}) where ElType
    # Scan the inputs for the TSeries:
    ts = find_tseries(bc)

    similar(ts)
end

"""
    find_tseries

Return the first TSeries among the arguments.

__Note:__ An internal function used for broadcasting support.
"""
find_tseries(bc::Base.Broadcast.Broadcasted) = find_tseries(bc.args)
find_tseries(args::Tuple) = find_tseries(find_tseries(args[1]), Base.tail(args))
find_tseries(x) = x
find_tseries(a::TSeries, rest) = a
find_tseries(::Any, rest) = find_tseries(rest)


text:
const _available_compressors = Dict{Symbol,Compressor}()
const _available_filters = Dict{Symbol, FilterDef}()

function __init__()
    fill_available_compressors!(_available_compressors)
    global _default_compressor_name = _get_default_compressor_name()

    fill_buildin_filters()
end


text:
export onehot, onehot!

using DataFrames: AbstractDataFrame

"""
    onehot(df, col, cate = sort(unique(df[!, col])); outnames = Symbol.(:ohe_, cate))

    onehot!(df, col, cate = sort(unique(df[!, col])); outnames = Symbol.(:ohe_, cate))

One-hot encode a column and create columns. The output columns will be overwritten WITHOUT warning

Arguments:

    df   -   The DataFrame
    col  -   The column to onehot encode
    cate -  The categories

"""
function onehot!(df::AbstractDataFrame, col, cate = sort(unique(df[!, col])); outnames = Symbol.(:ohe_, cate))
    transform!(df, @. col => ByRow(isequal(cate)) .=> outnames)
end

function onehot(df::AbstractDataFrame, col, cate = sort(unique(df[!, col])); outnames = Symbol.(:ohe_, cate))
    transform(df, @. col => ByRow(isequal(cate)) .=> outnames)
end

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
@testset "K8sClusterManager" begin
    @testset "pods not found" begin
        try
            K8sClusterManager(1)
        catch ex
            # Show the original stacktrace if an unexpected error occurred.
            ex isa KubeError || rethrow()

            @test ex isa KubeError
            @test length(Base.catch_stack()) == 1
        end
    end
end

@testset "TimeoutException" begin
    e = K8sClusterManagers.TimeoutException("time out!")
    @test sprint(showerror, e) == "TimeoutException: time out!"
end


text:
module aero105

#using Plots
#using DataInterpolations
using Interpolations

Mach_1 = [0, 0.875, 0.925, 0.965, 0.990, 1.025, 1.085, 1.19, 1.35, 1.80, 2.0, 2.5]
CD0 = [0.124, 0.124, 0.15, 0.2, 0.35, 0.375, 0.415, 0.415, 0.385, 0.335, 0.318, 0.276]
CD0_inter = LinearInterpolation(Mach_1, CD0, extrapolation_bc=Flat())

Mach_2 = [0, 0.88, 0.97, 0.99, 1.15, 1.25, 1.3, 2.5]
CDδ2 = [3.2, 3.2, 6.3, 4.0, 5.0, 5.4, 5.5, 5.5]
CDδ2_inter = LinearInterpolation(Mach_2, CDδ2, extrapolation_bc=Flat())

Mach_3 = [0, 0.43, 0.70, 0.91, 1.4, 1.75, 2.1, 2.5]
Cl_p = [-0.0178, -0.0149, -0.0135, -0.0126, -0.0110, -0.0101, -0.0094, -0.0087]
Cl_p_inter = LinearInterpolation(Mach_3, Cl_p, extrapolation_bc=Flat())

Mach_4 = [0, 0.4, 0.7, 0.89, 0.99, 1.09, 1.5, 2.0, 2.5]
Cl_α0 = [1.63, 1.63, 1.41, 1.22, 1.73, 1.57, 1.97, 2.25, 2.5]
Cl_α0_inter = LinearInterpolation(Mach_4, Cl_α0, extrapolation_bc=Flat())

Mach_5 = [0, 0.2, 0.6, 0.8, 0.985, 1.09, 1.3, 1.5, 2.0, 2.5]
Cl_α2 = [0.1, 0.1, 3.5, 6.6, 9.2, 8.8, 12.0, 13.7, 16.0, 17.0]
Cl_α2_inter = LinearInterpolation(Mach_5, Cl_α2, extrapolation_bc=Flat())

Mach_6 = [0, 0.46, 0.61, 0.78, 0.87, 0.925, 0.97, 1.09, 1.5, 2.5]
CM_α0 = [3.55, 3.55, 3.76, 3.92, 3.96, 4.85, 4.0, 3.83, 3.75, 3.75]
CM_α0_inter = LinearInterpolation(Mach_6, CM_α0, extrapolation_bc=Flat())

Mach_7 = [0, 0.4, 0.45, 0.65, 0.78, 0.885, 0.98, 1.075, 1.25, 1.5, 2.0, 2.5]
CM_α2 = [-2.9, -2.9, -3.1, -4.4, -3.45, -1.78, -3.0, -2.1, -3.325, -4.45, -4.6, -4.6]
CM_α2_inter = LinearInterpolation(Mach_7, CM_α2, extrapolation_bc=Flat())

Mach_8 = [0, 0.79, 1.15, 1.55]
CM_q_plus_CM_α = [-3.15, -3.15, -9.1, -9.5]
CM_q_plus_CM_α_inter = LinearInterpolation(Mach_8, CM_q_plus_CM_α, extrapolation_bc=Flat())

Mach_9 = [0, 0.22, 0.31, 0.48, 0.999, 1.001, 1.55]
α_t_sq = [0, 21.4, 348.5, 364.5, 632, 638, 706, 908, 1316]
CN_pα_Mach_inter = Vector{}(undef,length(Mach_9))

CN_pα_Mach1 = [-0.34, -0.91, -1.42, -2.63]
α_t_sq_Mach1 = [0, 632, 908, 1316]

CN_pα_Mach_inter[1] = LinearInterpolation(α_t_sq_Mach1, CN_pα_Mach1, extrapolation_bc=Flat())

CN_pα_Mach2 = [-0.34, -0.91, -1.42, -2.63]
α_t_sq_Mach2 = [0, 632, 908, 1316]

CN_pα_Mach_inter[2] = LinearInterpolation(α_t_sq_Mach2, CN_pα_Mach2, extrapolation_bc=Flat())

CN_pα_Mach3 = [-0.125, -0.465, -0.503, -1.015, -2.92]
α_t_sq_Mach3 = [0, 21.4, 364.5, 638, 1316]

CN_pα_Mach_inter[3] = LinearInterpolation(α_t_sq_Mach3, CN_pα_Mach3, extrapolation_bc=Flat())

CN_pα_Mach4 = [-0.34, -0.591, -2.45]
α_t_sq_Mach4 = [0, 348.5, 1316]

CN_pα_Mach_inter[4] = LinearInterpolation(α_t_sq_Mach4, CN_pα_Mach4, extrapolation_bc=Flat())

CN_pα_Mach5 = [-0.34, -0.591, -2.45]
α_t_sq_Mach5 = [0, 348.5, 1316]

CN_pα_Mach_inter[5] = LinearInterpolation(α_t_sq_Mach5, CN_pα_Mach5, extrapolation_bc=Flat())

CN_pα_Mach6 = [-0.36, -1.68]
α_t_sq_Mach6 = [0, 706]

CN_pα_Mach_inter[6] = LinearInterpolation(α_t_sq_Mach6, CN_pα_Mach6, extrapolation_bc=Flat())

CN_pα_Mach7 = [-0.36, -1.68]
α_t_sq_Mach7 = [0, 706]

CN_pα_Mach_inter[7] = LinearInterpolation(α_t_sq_Mach7, CN_pα_Mach7, extrapolation_bc=Flat())

CN_pα = Array{Float64,2}(undef, length(Mach_9), length(α_t_sq))
for i=1:length(Mach_9)

    for j=1:length(α_t_sq)
    CN_pα[i,j] = CN_pα_Mach_inter[i](α_t_sq[j])

    #println(CN_pα[i,j])
end
end

CN_pα_inter = LinearInterpolation((Mach_9,α_t_sq), CN_pα)

Mach_10 = [0, 0.22, 0.31, 0.48, 0.81, 0.87, 0.92, 0.96, 0.995, 1.02, 1.1, 1.21, 1.28, 1.46, 1.55]
α_t_sq_10 = [0, 27.5, 315.3, 322.2, 375.2, 403.6, 410.8, 630.2, 637.7, 705.7, 743.9, 915.9, 1316]
CM_pα_Mach_inter = Vector{}(undef,length(Mach_10))

CM_pα_Mach1 = [0.10, 0.173, 0.345, 2.35]
α_t_sq_10_Mach1 = [0, 403.6, 630.2, 1316]

CM_pα_Mach_inter[1] = LinearInterpolation(α_t_sq_10_Mach1, CM_pα_Mach1, extrapolation_bc=Flat())

CM_pα_Mach2 = [0.10, 0.173, 0.345, 2.35]
α_t_sq_10_Mach2 = [0, 403.6, 630.2, 1316]

CM_pα_Mach_inter[2] = LinearInterpolation(α_t_sq_10_Mach2, CM_pα_Mach2, extrapolation_bc=Flat())

CM_pα_Mach3 = [0.10, 0.133, 0.471, 1.276, 2.35]
α_t_sq_10_Mach3 = [0, 410.8, 637.7, 915.9, 1316]

CM_pα_Mach_inter[3] = LinearInterpolation(α_t_sq_10_Mach3, CM_pα_Mach3, extrapolation_bc=Flat())

CM_pα_Mach4 = [-0.46, 0.08, 0.022, 0.94]
α_t_sq_10_Mach4 = [0, 27.5, 375.2, 1316]

CM_pα_Mach_inter[4] = LinearInterpolation(α_t_sq_10_Mach4, CM_pα_Mach4, extrapolation_bc=Flat())

CM_pα_Mach5 = [-0.46, 0.08, 0.22, 0.94]
α_t_sq_10_Mach5 = [0, 27.5, 375.2, 1316]

CM_pα_Mach_inter[5] = LinearInterpolation(α_t_sq_10_Mach5, CM_pα_Mach5, extrapolation_bc=Flat())

CM_pα_Mach6 = [0.4175, 0.053, 0.285]
α_t_sq_10_Mach6 = [0, 315.3, 743.9]

CM_pα_Mach_inter[6] = LinearInterpolation(α_t_sq_10_Mach6, CM_pα_Mach6, extrapolation_bc=Flat())

CM_pα_Mach7 = [0.4175, 0.053, 0.285]
α_t_sq_10_Mach7 = [0, 315.3, 743.9]

CM_pα_Mach_inter[7] = LinearInterpolation(α_t_sq_10_Mach7, CM_pα_Mach7, extrapolation_bc=Flat())

CM_pα_Mach8 = [0.3747, 0.05, 0.665]
α_t_sq_10_Mach8 = [0, 322.2, 1316]

CM_pα_Mach_inter[8] = LinearInterpolation(α_t_sq_10_Mach8, CM_pα_Mach8, extrapolation_bc=Flat())

CM_pα_Mach9 = [0.3747, 0.05, 0.665]
α_t_sq_10_Mach9 = [0, 322.2, 1316]

CM_pα_Mach_inter[9] = LinearInterpolation(α_t_sq_10_Mach9, CM_pα_Mach9, extrapolation_bc=Flat())

CM_pα_Mach10 = [0.2, 0.301]
α_t_sq_10_Mach10 = [0, 375.2]

CM_pα_Mach_inter[10] = LinearInterpolation(α_t_sq_10_Mach10, CM_pα_Mach10, extrapolation_bc=Flat())

CM_pα_Mach11 = [0.2, 0.301]
α_t_sq_10_Mach11 = [0, 375.2]

CM_pα_Mach_inter[11] = LinearInterpolation(α_t_sq_10_Mach11, CM_pα_Mach11, extrapolation_bc=Flat())

CM_pα_Mach12 = [0.193, 0.5, 0.445]
α_t_sq_10_Mach12 = [0, 403.6, 705.7]

CM_pα_Mach_inter[12] = LinearInterpolation(α_t_sq_10_Mach12, CM_pα_Mach12, extrapolation_bc=Flat())

CM_pα_Mach13 = [0.193, 0.5, 0.445]
α_t_sq_10_Mach13 = [0, 403.6, 705.7]

CM_pα_Mach_inter[13] = LinearInterpolation(α_t_sq_10_Mach13, CM_pα_Mach13, extrapolation_bc=Flat())

CM_pα_Mach14 = [0.215, 0.495]
α_t_sq_10_Mach14 = [0, 410.8]

CM_pα_Mach_inter[14] = LinearInterpolation(α_t_sq_10_Mach14, CM_pα_Mach14, extrapolation_bc=Flat())

CM_pα_Mach15 = [0.215, 0.495]
α_t_sq_10_Mach15 = [0, 410.8]

CM_pα_Mach_inter[15] = LinearInterpolation(α_t_sq_10_Mach15, CM_pα_Mach15, extrapolation_bc=Flat())

CM_pα = Array{Float64,2}(undef, length(Mach_10), length(α_t_sq_10))
for i=1:length(Mach_10)

    for j=1:length(α_t_sq_10)
    CM_pα[i,j] = CM_pα_Mach_inter[i](α_t_sq_10[j])

    #println(CM_pα[i,j])
end
end

CM_pα_inter = LinearInterpolation((Mach_10,α_t_sq_10), CM_pα)

end


text:
# geothermal energy simulation; see http://fvt.simulkade.com/posts/2017-04-07-heat-transport-in-porous-media.html
# written by Ali A. Eftekhari, while on paternity leave!
# 14 June 2017; Lyngby public library
# output:
# pressure gradient (pressure difference between the injection and the production well)
# outlet temperature: temperature at the middle of the right boundary
# vs simulation time
#

# You can install CoolProp by:
# Pkg.clone("https://github.com/vimalaad/CoolProp.jl")
# Pkg.build("CoolProp")
using CoolProp, DataFrames, PyPlot, Polynomials, JFVM
import GR
# a table of density values
g_gravity = 9.81 # [m/s^2]
depth_reservoir = 2500         # [m]
T_res = 80 + 273.15            # [K]
T_inj = 40 + 273.15            # [K]
p_surf= 101325.0               # [Pa]
p_res = PropsSI("D", "T", T_res, "P", p_surf, "H2O")*g_gravity*depth_reservoir # [Pa]
T     = collect(linspace(T_inj, T_res, 10))
rho_water = zeros(length(T))
mu_water  = zeros(length(T))
for i in eachindex(T)
    rho_water[i] = PropsSI("D", "T", T[i], "P", p_res, "water")           # [kg/m^3]
    mu_water[i]  = PropsSI("viscosity", "T", T[i], "P", p_res, "water")   # [Pa.s]
end

mu_fit  = polyfit(T, mu_water, 2)
rho_fit = polyfit(T, rho_water, 2)
figure()
subplot(2,1,1)
scatter(T, mu_water, label = "μ data")
plot(T, polyval(mu_fit, T))
xlabel("T [K]")
ylabel("Viscosity [Pa.s]")
legend()
axis([minimum(T), maximum(T), minimum(mu_water), maximum(mu_water)])

subplot(2,1,2)
scatter(T, rho_water, label = "ρ data")
plot(T, polyval(rho_fit, T))
xlabel("T [K]")
ylabel("Density [kg/m^3]")
legend()

flow_m3_h = 100        # [m3/h]
flow_prod  = flow_m3_h/3600 # m3/s
flow_inj   = flow_m3_h/3600*PropsSI("D", "T", T_res, "P", p_surf, "H2O")/PropsSI("D", "T", T_inj, "P", p_surf, "H2O")
perm_D  = 0.02        # Darcy
perm    = perm_D*1e-12 # [m^2]
L_res    = 500            # [m] will be changed later in a loop
L_left = 500           # [m] distance between the left boundary and injection well
L_right = 500          # [m] distance between the right boundary and injection well
L = L_res + L_left + L_right # [m] domain length
W     = 500            # [m]
thick_res = 100        # [m]
poros = 0.20           # [-]
T0    = 25.0 + 273.15  # [K] reference temperature
V_dp  = 0.2            # Dykstra-Parsons coef
clx   = 0.1            # Correlation length in x direction
cly   = 0.1            # Correlation length in y direction
T_inj = 40 + 273.15    # [K]
theta_inj = T_inj - T0 # [K]
T_res = 80 + 273.15    # [K]
theta_res = T_res - T0 # [K]
λ_water  = PropsSI("conductivity", "T", T_res, "P", p_res, "water") # W/m/K
λ_rock   = 0.16        # [W/m/K] rock conductivity (for sandstone from my thesis)
λ_eff    = λ_water^poros*λ_rock^(1-poros) # effective conductivity
cp_water = PropsSI("CPMASS", "T", T_res, "P", p_res, "water") # J/kg/K
cp_rock  = 837.0       # [J/kg/K] rock heat capacity (for sandstone)
rho_rock = 2650.0      # [kg/m^3] rock density (for sandstone)

well_diameter = 7*0.0254 # [m]; converted from a well diameter of 7 inch
nx_left = 20
nx_right = 20
nx_well           = 10
ny_well           = 10
nx                = 40
ny                = 15
x_range =
  collect([linspace(0, L_left-(nx_well+1)*well_diameter, nx_left);
           linspace(L_left-nx_well*well_diameter, L_left+nx_well*well_diameter, 2*nx_well);
           linspace(L_left+(nx_well+1)*well_diameter, L_left+L_res-(nx_well+1)*well_diameter, nx);
           linspace(L_left+L_res-nx_well*well_diameter, L_left+L_res+nx_well*well_diameter, 2*nx_well);
           linspace(L_left+L_res+(nx_well+1)*well_diameter, L, nx_right)])
y_range = collect([linspace(0, 0.5*W-(ny_well+1)*well_diameter, ny);
                   linspace(0.5*W-ny_well*well_diameter, 0.5*W+ny_well*well_diameter, 2*ny_well);
                   linspace(0.5*W+(ny_well+1)*well_diameter, W, ny)])
ind_inj  = [nx_left+nx_well+1, ny+ny_well]
ind_prod = [nx_left+nx_well+nx_well+nx+nx_well, ny+ny_well]
m   = createMesh2D(x_range, y_range)  # 2D Cartesian grid
v_cell = thick_res*cellVolume(m)
Nx  = length(x_range)-1
Ny  = length(y_range)-1

perm_val   = perm # permfieldlogrnde(Nx, Ny, perm, V_dp, clx, cly)
perm_field = createCellVariable(m, perm_val)

left_range  = ny+ny_well:ny+ny_well
right_range = ny+ny_well:ny+ny_well
u_inj = flow_inj/(m.cellsize.y[ny+ny_well+1]*thick_res)
BCp = createBC(m)                 # pressure boundary
BCp.left.a[:] = 0.0
BCp.left.b[:] = 1.0
BCp.left.c[:] = p_res
BCp.right.a[:] = 0.0
BCp.right.b[:] = 1.0
BCp.right.c[:] = p_res

# BCp.left.a[left_range]   =
#     perm_field.value[2, left_range]/polyval(mu_fit, T_inj)
# BCp.left.b[left_range]   = 0.0
# BCp.left.c[left_range]   = -u_inj
# BCp.right.a[right_range] = 0.0
# BCp.right.b[right_range] = 1.0
# BCp.right.c[right_range] = p_res

BCt = createBC(m)                 # temperature boundary
# Danckwertz with assumptions (density of water is almost constant)
rho_water_inj  = polyval(rho_fit, T_inj)
# BCt.left.a[left_range]  = λ_eff
# BCt.left.b[left_range]  = -rho_water_inj*cp_water*u_inj
# BCt.left.c[left_range]  = -rho_water_inj*cp_water*u_inj*theta_inj

# figure(figsize=(8,2))
# visualizeCells(perm_field)
# title("Permeability [m^2]", fontsize = 10)
# colorbar()

theta_init = createCellVariable(m, theta_res, BCt)
theta_val  = createCellVariable(m, theta_res, BCt)
p_init     = createCellVariable(m, p_res, BCp)
p_val      = createCellVariable(m, p_res, BCp)
rho_init   = createCellVariable(m,
    polyval(rho_fit, theta_init.value+T0))
rho_val    = copyCell(rho_init)
mu_init    = createCellVariable(m,
    polyval(mu_fit, theta_init.value+T0))
mu_val     = copyCell(mu_init)

dt_init    = (L*W*thick_res*poros)/flow_inj/50 # (L/Nx)/(u_inj/poros)   # [s] time step
final_time = 100*dt_init              # [s]
t_step     = collect([linspace(0, 5*dt_init, 30);
                      linspace(5.5*dt_init, final_time, 70)])
n_steps = length(t_step) # number of time steps
# variables to be stored:
dp_res = zeros(n_steps)
T_out  = zeros(n_steps)
p_inj  = zeros(n_steps)
p_out  = zeros(n_steps)
T_out[1] = T_res

# discretization
M_BCp, RHS_BCp = boundaryConditionTerm(BCp)
M_BCt, RHS_BCt = boundaryConditionTerm(BCt)
M_conductivity = diffusionTerm(-harmonicMean(
        createCellVariable(m, λ_eff)))
q_inj = createCellVariable(m, 0)
q_inj.value[ind_inj+1 ...] = flow_inj/v_cell.value[ind_inj+1 ...]
RHSinjt = constantSourceTerm(rho_water_inj*cp_water*(T_inj-T0)*q_inj)

q_prod = createCellVariable(m, 0)
q_prod.value[ind_prod+1 ...] = -flow_prod/v_cell.value[ind_prod+1 ...]

q_cont = createCellVariable(m, 0)

for t_ind in 2:length(t_step)
  dt = t_step[t_ind]-t_step[t_ind-1]
  t  = t_step[t_ind]
  println(t_ind)
  for i in 1:3 # internal loop
    # solve pressure equation
    RHS_ddt_p   = constantSourceTerm(poros/dt*(rho_val - rho_init))
    water_mobil = harmonicMean(perm_field./mu_val)
    rho_face = arithmeticMean(rho_val)
    M_diff_p    = diffusionTerm(-rho_face.*water_mobil)
    q_cont.value[ind_inj+1 ...] = rho_val.value[ind_inj+1 ...]*flow_inj/v_cell.value[ind_inj+1 ...]
    q_cont.value[ind_prod+1 ...] = -rho_val.value[ind_prod+1 ...]*flow_prod/v_cell.value[ind_prod+1 ...]
    RHScont = constantSourceTerm(q_cont)
    p_val       = solveLinearPDE(m, M_diff_p + M_BCp,
        RHS_BCp - RHS_ddt_p + RHScont)

    # velocity vector
    u = -water_mobil.*gradientTerm(p_val)

    # solve heat equation
    α                  = poros*rho_val*cp_water+(1-poros)*rho_rock*cp_rock
    M_trans, RHS_trans = transientTerm(theta_init, dt, α)
    RHS_ddt_t          = constantSourceTerm(cp_water*poros*theta_val.*(rho_val - rho_init)/dt)
    M_conv             = convectionUpwindTerm(cp_water*rho_face.*u)
    # RHSprodt = constantSourceTerm(-rho_val.value[ind_prod+1 ...]*cp_water*T0*q_prod)
    q_inj.value[ind_inj+1 ...] = flow_inj/v_cell.value[ind_inj+1 ...]
    RHSinjt = constantSourceTerm(rho_water_inj*cp_water*(T_inj-T0)*q_inj)
    Mprodt   = linearSourceTerm(rho_val.value[ind_prod+1 ...]*cp_water*q_prod)
    theta_val          = solveLinearPDE(m,
        M_BCt + M_conv + M_trans + M_conductivity - Mprodt,
        RHS_BCt + RHS_trans + RHSinjt - RHS_ddt_t)

    # update density and viscosity values
    rho_val.value[:] = polyval(rho_fit, theta_val.value[:] + T0)
    mu_val.value[:]  = polyval(mu_fit, theta_val.value[:] + T0)
    flow_inj   = flow_prod*rho_val.value[ind_prod+1 ...]/rho_water_inj
  end # end of inner loop
  # GR.imshow(internalCells(theta_val)+T0)
  T_out[t_ind]  = theta_val.value[ind_prod+1 ...]+T0
  p_inj[t_ind]  = p_val.value[ind_inj+1 ...]
  p_out[t_ind]  = p_val.value[ind_prod+1 ...]
  dp_res[t_ind] = p_val.value[ind_inj+1 ...]-p_val.value[ind_prod+1 ...]
                #  0.5*(p_val.value[end, 1+ny+ny_well]+p_val.value[end-1, 1+ny+ny_well])
  rho_init   = copyCell(rho_val)
  p_init     = copyCell(p_val) # not necessary
  theta_init = copyCell(theta_val)
end
dp_res[1] = dp_res[2] # pressure grad at time zero; assume that process has already started
p_inj[1] = p_inj[2]
p_out[1] = p_out[2]
df1 = DataFrame(T_s = t_step, dp_Pa = dp_res, T_K = T_out)
# DataFrames.writetable("Q-$flow_m3_h-L-$flow_m3_h-k$perm_D.csv", df1)
figure(figsize=(8,2))
visualizeCells(theta_val+T0)
title("temperature profile")
colorbar()
figure()
plot(t_step/(365*24*3600), T_out, ".")
xlabel("time [year]")
ylabel("Water temperature [K]")
figure()
plot(t_step/(365*24*3600), dp_res/1e5)
xlabel("time [year]")
ylabel("Pressure difference [bar]")
figure()
plot(t_step/(365*24*3600), p_inj/1e5)
plot(t_step/(365*24*3600), p_out/1e5)
xlabel("time [year]")
ylabel("Pressure [bar]")

max_stars_repo_name: simulkade/geothermal

max_issues_repo_name: simulkade/geothermal

max_forks_repo_name: simulkade/geothermal


text:
using Test

for submodule in ["Grids",
                  "StateVecs",
                  "StateVecFuncs",
                  ]

  println("Testing $submodule")
  include(joinpath("../test",submodule*".jl"))
end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
module HadoopBlocks

using Base.Order
using URIParser
using Compat
using Elly
using Blocks

import Blocks: Block

import Elly: status

import Base: close, eof, read, write, readbytes, peek, readall, flush, nb_available, position, filesize, seek, seekend, seekstart, skip,
             start, done, next, wait

export  close, eof, read, write, readbytes, peek, readall, flush, nb_available, position, filesize, seek, seekend, seekstart, skip,
        # from hdfs_jobs.jl
        dmap, dmapreduce, results, status, unload, wait, times, JobId,
        # from hdfs_reader.jl
        MRInput, BlockIO,
        MRMapInput, MRHdfsFileInput, MRFsFileInput,
        MapResultReader, HdfsBlockReader, FsBlockReader

export Block

abstract MapInputReader
abstract MapStreamInputReader <: MapInputReader
abstract MRInput

function close(r::MapStreamInputReader)
    try
        close(get_stream(r))
    catch
        # do nothing
    end
end

function is_loaded(x::Symbol)
    try
        return (Module == typeof(eval(Main,x)))
    catch
        return false
    end
end

# enable logging only during debugging
using Logging
const logger = Logging.configure(level=DEBUG)
#const logger = Logging.configure(filename="/tmp/elly$(getpid()).log", level=DEBUG)
logmsg(s) = debug(s)
#logmsg(s) = nothing

typealias JobId     Int64
typealias FuncNone  Union{Function,Void}

include("scheduler.jl")
using HadoopBlocks.Scheduler
Scheduler.prep_remotes()

include("hdfs_blocks.jl")
include("fs_reader.jl")
include("hdfs_reader.jl")
include("map_result_reader.jl")

include("hdfs_jobs.jl")
include("hdfs_mrutils.jl")

#if(is_loaded(:DataFrames))
#    include("hdfs_mrutils_dataframes.jl")
#end

end # module

max_stars_repo_name: tanmaykm/HadoopBlocks

max_issues_repo_name: tanmaykm/HadoopBlocks

max_forks_repo_name: tanmaykm/HadoopBlocks


text:
cubic(x::Real) = x * x * x

function apy2(x::Real, y::Real)
    isnan(x) && return x
    isnan(y) && return y
    # general case
    xabs = abs(x)
    yabs = abs(y)
    w = max(xabs, yabs)
    z = min(xabs, yabs)
    iszero(z) && return w
    return w * sqrt(1.0 + abs2(z / w))
end

function sqr2(x::Real, y::Real)
    isnan(x) && return x
    isnan(y) && return y
    # general case
    xabs = abs(x)
    yabs = abs(y)
    w = max(xabs, yabs)
    z = min(xabs, yabs)
    iszero(z) && return abs2(w)
    return abs2(w) * (1.0 + abs2(z / w))
end

function unsafe_copy!(des::AbstractArray, src::AbstractArray)
    @simd for i in eachindex(des)
        @inbounds des[i] = src[i]
    end
    return nothing
end

nrm2(x::VecI{Tx}) where Tx<:Real = sqrt(dot(x, x))

function nrm2(x::VecI{Tx}, y::VecI{Ty}, b::VecB{Tb}) where {Tx<:Real,Ty<:Real,Tb<:Real} # @code_warntype ✓
    @simd for i in eachindex(b)
        @inbounds b[i] = abs2(x[i] - y[i])
    end
    return sqrt(sum(b))
end

function swap!(v::VecIO, i::Int, j::Int) # @code_warntype ✓
    @inbounds temp = v[i]
    @inbounds v[i] = v[j]
    @inbounds v[j] = temp
    return nothing
end

function swap!(m::MatI, i1::Int, j1::Int, i2::Int, j2::Int) # @code_warntype ✓
    @inbounds temp     = m[i1,j1]
    @inbounds m[i1,j1] = m[i2,j2]
    @inbounds m[i2,j2] = temp
    return nothing
end

function checksquare(A::MatI)
    m, n = size(A)
    m == n || error("tr: A should be a square matrix")
    return n
end

function diagmax(A::MatI, n::Int)
    r = -Inf
    @inbounds for i in eachindex(1:n)
        r = max(r, A[i,i])
    end
    return r
end

diagmax(A::MatI) = (n = checksquare(A); return diagmax(A, n))

function tr(A::MatI, n::Int)
    ret = 0.0
    @inbounds for i in eachindex(1:n)
        ret += A[i,i]
    end
    return ret
end

"""
    tr(A)
Matrix trace. Sums the diagonal elements of `A`.
`A` must be a square matrix.
"""
tr(A::MatI) = (n = checksquare(A); return tr(A, n))

function dot(x::VecI, y::VecI, n::Int)
    r = 0.0
    m = mod(n, 5)
    if m ≠ 0
        @inbounds for i in 1:m
            r += x[i] * y[i]
        end
        n < 5 && return r
    end
    m += 1
    @inbounds for i in m:5:n
        r += x[i] * y[i] + x[i+1] * y[i+1] + x[i+2] * y[i+2] + x[i+3] * y[i+3] + x[i+4] * y[i+4]
    end
    return r
end

function dot(x::VecI, m::Int, A::MatI, y::VecI, n::Int)
    ret = 0.0
    for j in eachindex(1:n)
        @inbounds yj = y[j]
        if !iszero(yj)
            tmp = 0.0
            @inbounds for i in eachindex(1:m)
                tmp += A[i,j] * x[i]
            end
            ret += tmp * yj
        end
    end
    return ret
end

"""
    dot(x, y)
Compute the dot product `x' * y` between two vectors `x` and `y`.
"""
function dot(x::VecI, y::VecI)
    n = length(x)
    n == length(y) || error("dot: length(y) ≠ length(x) = $n.")
    return dot(x, y, n)
end

"""
    dot(x, A, y)
Compute the dot product `x' * A * y` between two vectors `x` and `y`.
"""
function dot(x::VecI, A::MatI, y::VecI)
    m, n = size(A)
    m == length(x) || error("dot: length(x) ≠ size(A, 1) = $m.")
    n == length(y) || error("dot: length(y) ≠ size(A, 2) = $n.")
    return dot(x, m, A, y, n)
end

function logdet(A::MatI, n::Int)
    r = 0.0
    @inbounds for i in eachindex(1:n)
        r += log(abs(A[i,i]))
    end
    return r
end

"""
    logdet(A)
Compute the log-determinant `log(det(A))` where `A` is a triangular matrix.
"""
logdet(A::MatI) = (n = checksquare(A); return logdet(A, n))


text:
"""
$(TYPEDEF)

# Fields
$(TYPEDFIELDS)
"""
struct PermutedBlock <: MultiArmRandomisationModel
    target::AbstractVector{<:AbstractFloat}
    blocksize::Integer
    weights::AbstractVector{<:Integer}
    W::Integer
    sequence::AbstractVector{<:Integer}
    dist::AbstractVector{<:Integer}
    PermutedBlock(
        target::AbstractVector{<:AbstractFloat},
        blocksize::Integer,
        weights::AbstractVector{<:Integer},
        W::Integer,
        sequence::AbstractVector{<:Integer},
        dist::AbstractVector{<:Integer},
    ) = new(target, blocksize, weights, W, sequence, dist)
end

"""
$(TYPEDSIGNATURES)
"""
function PermutedBlock(weights::AbstractVector{<:Integer}, blocksize::Int)
    if !all(weights .≥ zero(eltype(weights)))
        throw(DomainError(weights, "All values must be ≥ 0"))
    end
    target = normalize(weights, 1)
    W = sum(weights)
    if mod(blocksize, W) != 0
        throw(DomainError(blocksize, "blocksize not integer multiple of sum(weights)."))
    end
    PermutedBlock(target, blocksize, weights, W, zeros(Int, 0), zeros(Int, length(target)))
end

"""
$(TYPEDSIGNATURES)
"""
function PermutedBlock(target::AbstractVector{<:AbstractFloat}, blocksize::Int)
    if !isposvec(target)
        throw(DomainError(target, "All values must be ≥ 0"))
    end
    if !isnormvec(target)
        target = target ./ sum(target)
    end
    weights = round.(Int, target .* blocksize)
    W = sum(weights)
    if mod(blocksize, W) != 0
        errstr = "blocksize not integer multiple of sum(target*blocksize)."
        throw(DomainError(blocksize, errstr))
    end
    PermutedBlock(target, blocksize, weights, W, zeros(Int, 0), zeros(Int, length(target)))
end

Base.deepcopy(m::PermutedBlock) = PermutedBlock(deepcopy(m.target), deepcopy(m.blocksize))

function show(io::IO, PB::PermutedBlock)
    k = length(PB.target)
    n = length(PB.sequence)
    print(io, "PermutedBlock(arms = $k, blocksize = $(blocksize(PB)), n = $n)")
end

blocksize(PB::PermutedBlock) = PB.blocksize
weights(PB::PermutedBlock) = PB.weights
target(PB::PermutedBlock) = PB.target

"""
$(TYPEDSIGNATURES)
"""
function prob(PB::PermutedBlock)
    b = blocksize(PB)
    w = weights(PB)
    W = PB.W
    λ = b ÷ W
    n = dist(PB)
    i = sum(n) + 1
    k = floor((i - 1) ÷ b) # completed blocks
    return (λ * w + λ * w * k - n) / (b + b * k - (i - 1))
end

"""
$(TYPEDSIGNATURES)
"""
function update!(PB::PermutedBlock, y::Int)
    push!(PB.sequence, y)
    PB.dist[y] += 1
    return nothing
end


text:
module PlotUtils

# import ..TabularTMazeUtils
using FileIO
using JLD2

function build_tabular_tmaze_heatmap(itm)
    values = load_result(itm, :value_map)
    map_set = @load "../src/data/TTMazeValueSet.jld2" ValueSet

    num_steps = size(values)[1]

    value_heatmap = zeros(num_steps, 9, 7)

    # position_mapping = TabularTMazeUtils.valid_state_mask()
    function valid_state_mask()
        """
        return a mask of valid states that is 9x7
        """
        world = [["G1", "0", "0", "0", "0", "0", "G3"],
                          ["1", "0", "0", "0", "0", "0", "1"],
                          ["1", "1", "1", "1", "1", "1", "1"],
                          ["1", "0", "0", "1", "0", "0", "1"],
                          ["G2", "0", "0", "1", "0", "0", "G4"],
                          ["0", "0", "0", "1", "0", "0", "0"],
                          ["0", "0", "0", "1", "0", "0", "0"],
                          ["0", "0", "0", "1", "0", "0", "0"],
                          ["0", "0", "0", "1", "0", "0", "0"]]

        world = permutedims(hcat(world...))
        valid_states = findall(x-> x!="0", world)
        return valid_states
    end
    position_mapping = valid_state_mask()

    for i in 1:num_steps

        for (set_ind,obs) in enumerate(ValueSet["states"])
            heatmap_ind = position_mapping[Int(obs[1])]
            value_heatmap[i,heatmap_ind] = values[i,set_ind,1]
        end
    end
    return value_heatmap
end

"""
    Given an ItemCollection, returns all logs at the logger_key
"""
function load_results(ic, logger_key; return_type = "tensor")
    num_results = length(ic)
    results = []
    for itm in ic.items
        data = FileIO.load(joinpath(itm.folder_str, "results.jld2"))["results"]
        push!(results,data[logger_key])
    end

    if return_type == "tensor"
        return cat(results..., dims = 3)
    elseif return_type == "array"
        return results
    end
end

function load_result(itm, logger_key)
    data = FileIO.load(joinpath(itm.folder_str, "results.jld2"))["results"]
    return data[logger_key]
end

end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using FactCheck, LibBSON, Mongo, DataStructures

facts("Mongo") do
    client = MongoClient()
    collection = MongoCollection(client, "foo", "bar")
    oid = BSONOID()

    context("variables") do
      @fact collection.client --> client
      @fact collection.db --> "foo"
      @fact collection.name --> "bar"
    end

    context("insert") do
        insert(collection, ("_id" => oid, "hello" => "before"))
        @fact count(collection, ("_id" => oid)) --> 1
        @fact count(collection) --> 1
        for item in find(collection, ("_id" => oid), ("_id" => false, "hello" => true))
            @fact dict(item) --> Dict("hello" => "before")
        end
    end

    context("update") do
        update(
            collection,
            ("_id" => oid),
            set("hello" => "after")
            )
        @fact count(collection, ("_id" => oid)) --> 1
        for item in find(collection, ("_id" => oid), ("_id" => false, "hello" => true))
            @fact dict(item) --> Dict("hello" => "after")
        end
    end

    context("command_simple") do
        reply = command_simple(
            client,
            "foo",
            OrderedDict(
               "count" => "bar",
               "query" => Dict("_id" => oid))
            )
        @fact reply["n"] --> 1
    end

    context("delete") do
        delete(
            collection,
            ("_id" => oid)
            )
        @fact count(collection, ("_id" => oid)) --> 0
        @fact count(collection) --> 0
    end
end

facts("Mongo: bad host/port") do
    client = MongoClient("bad-host-name", 9999)
    collection = MongoCollection(client, "foo", "bar")
    @fact_throws insert(collection, ("foo" => "bar"))
end

facts("Query building helpers") do
    client = MongoClient()
    ppl = MongoCollection(client, "foo", "ppl")
    person(name, age) = insert(ppl, ("name" => name, "age" => age))
    person("Tim", 25)
    person("Jason", 21)
    person("Jim", 87)
    context("orderby") do
        @fact first(find(ppl, (query(), orderby("age" => -1))))["name"] --> "Jim"
        @fact first(find(ppl, (query(), orderby("age" => 1))))["name"] --> "Jason"
    end
    context("gt and lt") do
        @fact first(find(ppl, query("age" => lt(25))))["name"] --> "Jason"
        @fact first(find(ppl, query("age" => gt(50))))["name"] --> "Jim"
    end
    context("in and nin") do
        @fact first(find(ppl, query("age" => in([21]))))["name"] --> "Jason"
        @fact first(find(ppl, query("age" => nin([21,25]))))["name"] --> "Jim"
    end
    context("eq and ne") do
        @fact first(find(ppl, query("age" => eq(21))))["name"] --> "Jason"
        @fact first(find(ppl, query("age" => ne(87))))["name"] == "Jim" --> false
    end
    context("update with operator") do
        update(ppl, ("age" => 87), set("age" => 88))
        @fact first(find(ppl, query("name" => "Jim")))["age"] --> 88
    end
    delete(ppl, ())
end


text:
using SamplingHS3

## SMALL TEST:

G_test         = [eye(2);-eye(2)]
h_test         = [10.;10;10;10]
f_test         = x1 -> sqrt(sum(x1.^2,1))[1]
threshold_test = 4.0
x              = zeros(2)
Sigma_test     = 4.0*eye(2)
X_Collected    = zeros(2,0)

mh_tracking_data = SamplingHS3.MHSamplingTrackingData()

for zi in 1:10_000
    X_Collected = hcat( X_Collected , x )
    mhstep_result =  SamplingHS3.mh_step(f_test,threshold_test,G_test,h_test[:],x[:],Sigma_test,10,mh_tracking_data)
    x = mhstep_result[2]
end

if(false)
    using Plots
    Plots.plot(X_Collected[1,:],X_Collected[2,:])
end


text:
function getvalidpermnos(df::DataFrame)::Vector{Int}
  permnosym = :permno ∈ names(df) ? :permno : :lpermno
  validpermnos::Vector{Int} = unique(df[!, permnosym])

  return validpermnos
end

function getvalidpermnos()::Vector{Int}
  local ccm::DataFrame = loadccm()
  ccm = preprocessccm!(ccm)

  return getvalidpermnos(ccm)
end

#do some pre-filtering to reduce the processing time
function prefiltercrsp!(crsp::DataFrame, validpermnos::Vector{Int})::Nothing
  #rename columns to lower case

  crsp.tokeep = trues(size(crsp,1))

  crspgrp::GroupedDataFrame = groupby(crsp, :permno)
  @mpar for i ∈ 1:length(crspgrp)
    scrsp::SubDataFrame = crspgrp[i]

    if scrsp.permno[1] ∉ validpermnos
      scrsp.tokeep .= false
    end
  end

  filter!(r::DataFrameRow->r.tokeep, crsp)
  select!(crsp, Not([:tokeep]))

  return nothing

end

@inline normalizecrspnames!(crsp::DataFrame) = names!(crsp, (s::Symbol->Symbol(lowercase(string(s)))).(names(crsp)))

function prepcrsp(op!::Function = (noop(x)=x);
  crsptype = :not_set,
  refreshcrsp=true,
  datapath=DATA_PATH,
  incsvstream::Function = IN_CSV_STREAM,
  injlsstream::Function = IN_JLS_STREAM,
  outjlsstream::Function = OUT_JLS_STREAM,
  validpermnos::Union{Nothing, Vector{Int}} = getvalidpermnos(),
  prefix::String = "",
  crspcolumns::Union{Nothing,Vector{Symbol}} = [:permno, :date, :ret, :lret],
  testoutput::Bool = TEST_OUTPUT,
  archive::Bool=true,
  yearrange::UnitRange{Int} = YEAR_RANGE[])::DataFrame

  crsptype ∈ [:monthly, :daily, :dailyexpanded] || error("Set crsptype to :monthly or :daily instead of $crsptype")

  local crspname
  local crspparts

  if crsptype == :monthly
    crspname = CRSPM_NAME
    crspparts = CRSPM_PARTS
  elseif crsptype==:daily
    crspname = CRSPD_NAME
    crspparts = CRSPD_PARTS
  elseif crsptype==:dailyexpanded
    crspname = CRSPD_EXPANDED_NAME
    crspparts = CRSPD_EXPANDED_PARTS
  else
    @assert false
  end

  local p::String = "$datapath\\$(prefix)$crspname"
  local crsp::DataFrame

  #makes a serialization object so we don't have to keep parsing the CSV
  if refreshcrsp
    #do the below so we can pre-filter

    pin::String = "$datapath\\$(crspparts[1])"
    crsp = incsvstream("$pin.csv.lz4") |> CSV.read |> DataFrame
    normalizecrspnames!(crsp)
    (!isnothing(validpermnos)) && (prefiltercrsp!(crsp, validpermnos))

    if length(crspparts) > 1
      #try to do this is in a space efficient way
      for crsppart ∈ crspparts[2:end]
        pin = "$datapath\\$crsppart"
        local addedcrsp::DataFrame = incsvstream("$pin.csv.lz4") |> CSV.read |> DataFrame
        normalizecrspnames!(addedcrsp)
         (!isnothing(validpermnos)) && (prefiltercrsp!(addedcrsp, validpermnos))
        append!(crsp, addedcrsp)
      end
    end

    #filter by year first
    crsp.date = (s->parsecrsp(Date, s)).(crsp.date)
    filter!(r::DataFrameRow->year(r.date) ∈ yearrange, crsp) #restrict to desired years

    crsp = preprocesscrsp(crsp, crsptype)
    op!(crsp)
    select!(crsp, crspcolumns)

    if testoutput
      d = Date(2015,12,31)
      testoutput && CSV.write("output\\$(prefix)_$(crspname)_$d.csv", crsp[(cd->cd==d).(crsp.date), :])
    end

    local istream  = outjlsstream("$p.jls.lz4")
    local wrotefile::Bool = false
    #=try
      serialize(istream, crsp)
      wrotefile = true
    finally
      flush(istream)
      close(istream)
      wrotefile || error("Write to file failed $err")
    end=#

    outjlsstream("$p.jls.lz4") do s
      serialize(s, crsp)
    end

    archive && archivefile("$(prefix)$crspname.jls.lz4") #never hurts to have a backup

  else
    try
      injlsstream("$p.jls.lz4") do s
        crsp = deserialize(s)
      end
    catch err
      @warn "Could not load $(prefix)$crspname. Error: $err"
      @info "Attempting to restore back-up and load from archive..."
      unarchivefile("$(prefix)$crspname.jls.lz4")
      injlsstream("$p.jls.lz4") do s
        crsp = deserialize(s)
      end
      @info "Successfully restored from archive."
    end
  end

  @info "CRSP data loaded and/or saved into file $(prefix)$crspname.jls.lz4."

  return crsp
end

function computecrspreturns!(crsp::DataFrame)
  crsp.ret = (s::MString->parsecrsp(Float64, s)).(crsp.ret)
  crsp.dlret= (s::MString->parsecrsp(Float64, s)).(crsp.dlret)

  rename!(crsp, :ret=>:retlisted)
  crsp.retlisted .= (r->
    !ismissing(r) && (r > -1.0) ? r : missing).(crsp.retlisted)
  crsp.dlret .= (r->
    !ismissing(r) && (r > -1.0) ? r : missing).(crsp.dlret)

  crsp.ret = ((retlisted,dlret)->
    ismissing(retlisted) && ismissing(dlret) ? missing :
      ismissing(dlret) ? retlisted :
      ismissing(retlisted) ? dlret :
      (1+retlisted)*(1+dlret) - 1.
    ).(crsp.retlisted,crsp.dlret)

  #winsorizekeydata!(crsp, fields2winsorize = [:ret])

  #get log returns
  crsp.lret = (log).(1.0 .+ crsp.ret)
  deleterows!(crsp, (ismissing).(crsp.ret))

  return crsp
end

function filtercrsp!(crsp::DataFrame; filternonfinancials::Bool = false)

    #NOTE: SHRCD==10,11 : US Common Shares
    #NOTE:  :EXCHCD ∈ [1,2,3] : NYSE, Amex, NASDAQ
    deleterows!(crsp, (!).(completecases(crsp[:,[:shrcd, :exchcd]])))
    local shrcdvalid::Set{Int} = Set{Int}((10,11))
    local exchcdvalid::Set{Int} = Set{Int}((1,2,3))
    filter!(r::DataFrameRow->
      (r.shrcd ∈ shrcdvalid) && (r[:exchcd] ∈ exchcdvalid), crsp)

    #NOTE: r[:SICCD] ∉ 6000:6999 : No financial companies

    if filternonfinancials
        error("Need to get the data again- sic code no longer incldued in crsp file")
        crsp.siccd = (s::MString->parsecrsp(Int, s)).(crsp.siccd)
        local siccdinvalid::Set{Int} = Set{Int}(6000:6999)
        filter!(r::DataFrameRow->
          (!ismissing(r.siccd)) && (r.siccd ∉ siccdinvalid), crsp)
    end

    return crsp
end

function computecrspprices!(crsp::DataFrame, crsptype::Symbol)
  #Now work on the price
  #in this case, use the mean of the two values
  #price cannot be zero
  crsp.prc .= (
    prc->coalesce(iszero(prc),true) ? missing : prc).(crsp.prc)
  crsp.dlprc .= (
    dlprc->coalesce(iszero(dlprc),true) ? missing : dlprc).(crsp.dlprc)

  crsp.price = ((prc,dlprc)->
    ismissing(prc) && ismissing(dlprc) ? missing :
      ismissing(dlprc) ? prc :
      ismissing(prc) ? dlprc :
      (prc+dlprc)/2
  ).(crsp.prc, crsp.dlprc)

  #add in the alternative price if nothing else is available
  if crsptype == :monthly
    crsp.altprc .= (
      altprc->coalesce(iszero(altprc),true) ? missing : altprc).(crsp.altprc)
    crsp.altprc .= (abs).(crsp.altprc)
    crsp.price .= ((price,altprc)->
      ismissing(price) ? altprc : price).(crsp.price,crsp.altprc)
  end

  #negative values computed via bid-ask spread method
  crsp.price .= (abs).(crsp.price)

  return crsp
end

function computecrspmktcap!(crsp::DataFrame)

  #adjust the units and compute the market cap
  crsp.mktcap = Vector{MInt}((round).(crsp.price .* crsp.shrout .* 1000))
  crsp.mktcap .= (i::MInt->coalesce(iszero(i),true) ? missing : i).(crsp.mktcap)

  deleterows!(crsp, (ismissing).(crsp.mktcap))
  return crsp
end

function computenetreturns!(crsp::DataFrame)

  #look up the benchmark against an index
  #crsp doesn't handle irregular intervals very well
  mergeewindd!(crsp)
  crsp.lnet = Vector{MFloat64}(undef, size(crsp,1))

  sort!(crsp, [:permno, :date])

  scrsps::GroupedDataFrame = groupby(crsp, :permno)
  print("time to check (net ret formation): ")
  @time @mpar for i ∈ 1:length(scrsps)
    scrsp::SubDataFrame = scrsps[i]

    scrsp[2:end, :lnet] .= scrsp[2:end, :lret] .- (scrsp[2:end, :lewindd] .- scrsp[1:(end-1), :lewindd])

  end

  #compute the logged versions
  crsp.net = (exp).(crsp.lnet) .- 1.0

  #=if TEST_OUTPUT
    CSV.write("output\\crspnetseries.csv", crsp[(end-10_000):end, :])
    @info "crspnetseries.csv written."
  end=#

  return nothing

end

function preprocesscrsp(crsp::DataFrame, crsptype::Symbol;
  crspdateformat::DateFormat = CRSP_DATE_FORMAT,
  testoutput::Bool = TEST_OUTPUT,
  usedcrspfields::Vector{Symbol} = USED_CRSP_FIELDS,
  months2stale::Int = MONTHS_2_STALE_SHORT)::DataFrame where T<:DatePeriod

  local scrsps::GroupedDataFrame

  filtercrsp!(crsp)

  #Doing this before the dedup eliminates rows without returns,
  #prices and shares outstanding first. Implicitly requires that
  #these fields exist
  computecrspreturns!(crsp)
  computecrspprices!(crsp, crsptype)

  #now dedup crsp
  dedupctr::Int = 0
  crsp.keep = trues(size(crsp,1))

  scrsps = groupby(crsp, [:permno, :date])
  @mpar for i ∈ 1:length(scrsps)
    local scrsp::SubDataFrame = scrsps[i]
    if size(scrsp,1) > 1
      #first try to pick records by number missing
      nmissing::Vector{Int} = zeros(Int, size(scrsp,1))
      for f ∈ usedcrspfields
        nmissing .+= (ismissing).(scrsp[!,f])
      end
      minmissing::Int = minimum(nmissing)
      scrsp.keep .= nmissing .== minmissing

      #keep the top record iff it is identical to the others
      if sum(scrsp.keep) > 1
        sscrsp = view(scrsp, scrsp.keep, :)
        nchecks::Int = size(sscrsp,1)
        keep1::Bool = true
        for i ∈ 2:nchecks
          for f ∈ usedcrspfields
            #checks if the given used field is equal to the top row
            (coalesce(sscrsp[i,f],-99.0) ≠ coalesce(sscrsp[1,f], -99.0)) && (keep1 = false)
          end
        end

        sscrsp.keep .= false
        sscrsp.keep[1] = keep1
      end
    end
  end

  deleterows!(crsp, (!).(crsp.keep))

  return crsp
end

#applies a function to a vector of trailing variables
function computetrailing!(F::Function,
  crsp::DataFrame, periods::Int,  Ftarget::Symbol,   Fnew::Symbol;
  offset::Int=0, months2stale::Int = MONTHS_2_STALE_LONG,
  sorted::Bool = false)::Nothing

  periods2stale = Month(months2stale)

  N::Int = size(crsp,1)

  #allocate spcae for the new fields
  crsp[!,Fnew] = Vector{MFloat64}(undef, N)

  #need to sort
  (!sorted) && (sort!(crsp, [:permno, :date]))

  #avoid missing data in the returns for simplicity
  scrspnomissing::SubDataFrame = view(crsp, (!ismissing).(crsp[!,Ftarget]),:)
  scrsps::GroupedDataFrame = groupby(scrspnomissing, :permno)

  #now compute the trailing returns
  @mpar for i ∈ 1:length(scrsps)
    local scrsp = scrsps[i]

    Nscrsp = size(scrsp, 1)
    if periods ≤ Nscrsp #if there are enough  entries to compute the trailing returns
      for r ∈ periods:Nscrsp
        r₀::Int = r-periods+1-offset #this si the first row of the period
        rₜ::Int = r - offset
        staledate::Date = firstdayofmonth(scrsp[r,:date]) - periods2stale
        mindate::Date = scrsp[r₀,:date]

        #if the returns are not stale, record them
        if staledate ≤ mindate
          scrsp[r, Fnew] = F(scrsp[r₀:rₜ , Ftarget])
        end
      end
    end
  end

  return nothing
end

#computes the trailing returns from the logged returns
function computetrailingreturns!(crsp::DataFrame, periods::Int;
  months2stale = MONTHS_2_STALE_LONG,
  Flret::Symbol = :lret,
  Fnewret::Symbol = Symbol(:t,:ret,period),
  offset::Int = 0,
  sorted::Bool = false
  )::Nothing

  local Flnewret::Symbol = Symbol(:l, Fnewret)
  computetrailing!(sum, crsp, periods, Flret, Flnewret,
    months2stale=months2stale, offset=offset, sorted=sorted)

  #create the unlogged version
  crsp[!,Fnewret] = (exp).(crsp[!,Flnewret]) .- 1.0

  return nothing
end

#computes the trailing returns from the logged returns
function computetrailingvol!(crsp::DataFrame, periods::Int;
  months2stale = MONTHS_2_STALE_LONG,
  Fret::Symbol = :ret,
  Fvol::Symbol = Symbol(:t,:vol,periods),
  offset::Int = 0,
  sorted::Bool = false
  )::Nothing

  computetrailing!(std, crsp, periods, Fret, Fvol, months2stale=months2stale, offset=offset)

  if TEST_OUTPUT
    TEST_OUTPUT && CSV.write(
      "output\\vol_p$(periods)_m2s$(months2stale)_off$(offset).csv",
        crsp[(end-20000):end, :])
  end

  return nothing
end

function computevolumes!(crsp::DataFrame)::Nothing
  checkmissing(f::Union{Real,Missing}) = (
    ismissing(f) || (f == 99) || (f < 0) ? missing : f)

  crsp.vol .= (checkmissing).(crsp.vol)
  crsp.numtrd .= (checkmissing).(crsp.numtrd)

  return nothing
end

#make tryparse handle missing values
#Base.tryparse(::Type{<:Any}, v::Missing)::Nothing = nothing

#parse a value where the type is known: Generic case
function parsecrsp(::Type{T}, s::String;
    parsedmissings = CRSP_PARSED_MISSINGS) where T

  #if it doesn't parse to the right type, set to missing
  v::Union{T,Missing} = something(tryparse(T, s), missing)

  #check against the list of missing codes
  (!ismissing(v)) && (v ∈ parsedmissings) && (v=missing)

  return v
end

#the date case
parsecrsp(::Type{Date}, s::String;
    crspdateformat::DateFormat = CRSP_DATE_FORMAT) = Dates.Date(s, crspdateformat)

#helper methods and special cases
parsecrsp(::Type{<:Any}, v::Missing) = missing
parsecrsp(::Type{Date}, i::Int) = parsecrsp(Date, "$i")

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using LinearAlgebra, Plots, Statistics, TaylorSeries, Printf, JLD2, FileIO
# if !(pwd() in LOAD_PATH) push!(LOAD_PATH, pwd()) end
# using TaylorBigF

path1=realpath(dirname(@__FILE__)*"/..")
include(string(path1,"/src/TaylorBigF.jl"))

p=2000;setprecision(p); setrounding(BigFloat, RoundUp);
n=BigInt(20);# zigos lliws i prwti grammi tou Vall einai =0 sto x=0, (to idio k sta ODES)
nb=BigInt(2)*n+BigInt(2)*(n-BigInt(2)) # boundary points
dx=BigInt(1)/BigInt(99)
L=(n-BigInt(1))*dx
x0=zeros(BigFloat,n).-L/BigFloat("2")
x=Array{BigFloat}(undef,0)
for i=BigInt(1):n
    global x0,x,dx
    x=[x;x0.+(i-BigInt(1))*dx]
end
y=convert(Array{BigFloat},range(-L/BigFloat("2"),L/BigFloat("2"),length=n))
y=repeat(y,n)
# scatter(x,y)

V=TaylorBigF.calc_Vandermonde_2D(x,y,n*n,n*n+nb)
Vxxxx=TaylorBigF.calc_Vandermonde_xxxx_2D(x,y,n*n,n*n+nb)
Vyyyy=TaylorBigF.calc_Vandermonde_yyyy_2D(x,y,n*n,n*n+nb)
Vxxyy=TaylorBigF.calc_Vandermonde_xxyy_2D(x,y,n*n,n*n+nb)

Vall0=Vxxxx+2Vxxyy+Vyyyy
# boundaries
i_keep1=x.==-L/BigFloat("2")
i_keep2=x.==L/BigFloat("2")
i_keep3=y.==-L/BigFloat("2")
i_keep4=y.==L/BigFloat("2")
i_keep=i_keep1+i_keep2+i_keep3+i_keep4
i_b=i_keep.>0
# scatter!(x[i_b],y[i_b])
BigInt(sum(i_b))==nb
Vall1=TaylorBigF.calc_Vandermonde_2D(x[i_b],y[i_b],nb,n*n+nb)

Vall=[Vall0;Vall1]
q=[-ones(n*n);zeros(BigFloat,nb)]
@time a=Vall\q
maximum(abs.(Vall*a-q))
# scatter(x,y,zcolor=V*a)
# minimum(V*a)
# scatter3d(x,y,V*a,label=string("deflection (w), p=",p),legend=:topleft)

Vxxx=TaylorBigF.calc_Vandermonde_xxx_2D(x[i_b],y[i_b],nb,n*n+nb)
Vyyy=TaylorBigF.calc_Vandermonde_yyy_2D(x[i_b],y[i_b],nb,n*n+nb)
Vxyy=TaylorBigF.calc_Vandermonde_xyy_2D(x[i_b],y[i_b],nb,n*n+nb)
Vyxx=TaylorBigF.calc_Vandermonde_yxx_2D(x[i_b],y[i_b],nb,n*n+nb)
Qx=-(Vxxx+Vxyy)*a
Qy=-(Vyxx+Vyyy)*a
# scatter(x[i_b],y[i_b],zcolor=Qx)

i1=x.==-L/BigFloat("2")
i2=x.==L/BigFloat("2")
i12=i1+i2
i12=i12.>0
Vxxx=TaylorBigF.calc_Vandermonde_xxx_2D(x[i12],y[i12],2n,n*n+nb)
Vxyy=TaylorBigF.calc_Vandermonde_xyy_2D(x[i12],y[i12],2n,n*n+nb)
Qx=-(Vxxx+Vxyy)*a
# scatter(x[i12],y[i12],zcolor=Qx)
i3=y.==-L/BigFloat("2")
i4=y.==L/BigFloat("2")
i34=i3+i4
i34=i34.>0
Vyyy=TaylorBigF.calc_Vandermonde_yyy_2D(x[i34],y[i34],2n,n*n+nb)
Vyxx=TaylorBigF.calc_Vandermonde_yxx_2D(x[i34],y[i34],2n,n*n+nb)
Qy=-(Vyxx+Vyyy)*a
# scatter!(x[i34],y[i34],zcolor=Qy)

mean(abs.(Qx))*2*L+mean(abs.(Qy))*2*L-1*L*L

# point load
q=zeros(BigFloat,n*n+nb)
dd=zeros(BigFloat,n*n)
for i=1:n*n
    dd[i]=sqrt(x[i].^2+y[i].^2)
end
iso=sortperm(dd)[1:4]
# scatter(x[iso],y[iso])
q[iso].=-BigFloat("1")

@time a=Vall\q
maximum(abs.(Vall*a-q))
# scatter(x,y,zcolor=V*a)
minimum(V*a)
# scatter3d(x,y,V*a,label=string("deflection (w), p=",p),legend=:topleft)
# savefig(string("pdes_p_",p,".pdf"))

Qx=-(Vxxx+Vxyy)*a
Qy=-(Vyxx+Vyyy)*a
# scatter(x[i34],y[i34],zcolor=Qy)

mean(abs.(Qx))*2*L+mean(abs.(Qy))*2*L-1*dx*dx


text:
# Autogenerated wrapper script for SHOT_jll for i686-linux-gnu-cxx03
export amplexe, libshotsolver

using ASL_jll
using Cbc_jll
using Ipopt_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("SHOT")
JLLWrappers.@declare_library_product(libshotsolver, "libSHOTSolver.so")
JLLWrappers.@declare_executable_product(amplexe)
function __init__()
    JLLWrappers.@generate_init_header(ASL_jll, Cbc_jll, Ipopt_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libshotsolver,
        "lib/libSHOTSolver.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        amplexe,
        "bin/SHOT",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ReachabilityAnalysis
using Plots
using LaTeXStrings

# IVP parameters
A = [0.0 1; -4π 0]
X0 = BallInf([0.0, 10.0], 0.1)
X = Universe(2)
P = @ivp(x' = A*x, x ∈ X, x(0) ∈ X0)

# discretization parameters
δ = 0.1
alg = Forward()
k = 10

# result for original time step
Pd = discretize(P, δ, alg)
Ω0 = Pd.x0

# result for shrunk time step
δ_small = δ / k
Pd_small = discretize(P, δ_small, alg)
Ωs = ConvexHullArray()
push!(array(Ωs), Pd_small.x0)
Φ_small_pow = Pd_small.s.A
for i in 1:k-1
    push!(array(Ωs), Φ_small_pow * Pd_small.x0)
    global Φ_small_pow *= Pd_small.s.A
end
Ω0_small = concretize(Ωs)

# plot sets
plot(leg=:bottomleft, legendfontsize=20, tickfontsize=20)
plot!(Ω0, alpha=0, color=:blue, linecolor=:blue, linealpha=1, lw=2, lab=L"\Omega_0\ \delta = %$δ")
plot!(Ω0_small, alpha=0.3, color=:green, linecolor=:green, linealpha=1, lw=2, lab=L"\Omega_0\ \gamma = %$δ_small")
for i in eachindex(array(Ωs))
    plot!(Ωs[i], color=:gray, alpha=0.2, lab="")
end
plot!(X0, alpha=1, color="yellow", lab=L"X_0")
x0s = vcat(sample(X0, 0, include_vertices=true),
           sample(X0, 30, sampler=LazySets.FaceSampler(1)),
           sample(X0, 30))
for x0 in x0s
    local P = @ivp(x' = A*x, x ∈ X, x(0) ∈ x0)
    sol = solve(P, δ, ORBIT(δ=δ/10))
    plot!(sol, vars=(1, 2), c=:magenta, seriestype=:path, marker=:none)
end

savefig("shrink_delta.pdf")


text:
import FiniteDiff
import Calculus
import ForwardDiff

# TODO: derivative, second_derivative

include("../test_functions.jl")
include("../gradient_functions.jl")
include("../hessian_functions.jl")
include("../max_error.jl")

io = open("errors.tsv", "w")

@printf(io, "timestamp\tf\tx\tmethod\terr\tlatency\n")

for test_func in test_funcs
    for x1 in test_func.x1s
        for x2 in test_func.x2s
            x = [x1, x2]

            g_true = ForwardDiff.gradient(test_func.f, x)
            h_true = ForwardDiff.hessian(test_func.f, x)

            # TODO: Make these latencies less crazy for the first run of each
            # function.
            for grad_func in grad_funcs
                latency = @elapsed g_approx = grad_func.f(test_func.f, x)
                @printf(
                    io,
                    "%s\t%s\t%s\t%s\t%s\t%s\n",
                    Dates.now(),
                    test_func.name,
                    x,
                    grad_func.name,
                    max_error(g_true, g_approx),
                    latency,
                )
            end

            for hess_func in hess_funcs
                latency = @elapsed h_approx = hess_func.f(test_func.f, x)
                @printf(
                    io,
                    "%s\t%s\t%s\t%s\t%s\t%s\n",
                    Dates.now(),
                    test_func.name,
                    x,
                    hess_func.name,
                    max_error(h_true, h_approx),
                    latency,
                )
            end
        end
    end
end

close(io)


text:
mutable struct SPClamp{T<:VariateType} <: SumProductRule{Clamp{T}} end
outboundType(::Type{SPClamp{T}}) where T<:VariateType = Message{PointMass, T}
isApplicable(::Type{SPClamp{T}}, input_types::Vector{<:Type}) where T<:VariateType = true


text:
using HYT939Logger
using Test

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


text:
# This code presents an example to illustrate how NLboxsolve.jl can be used

using NLboxsolve

function test_nlboxsolve_outplace()

    function test3(x)

        f = similar(x)

        f[1] = 0.5*x[1] + x[2] - x[3] + x[1]*x[2] - 1.5
        f[2] = x[1] - 0.5*x[2] + x[3] + x[2]^2 - 2.5
        f[3] = -x[1] + x[2] + 0.5*x[3] + x[2]*x[3] - 1.5
        f[4] = x[1]^2 - x[2]^2 + x[3]^2 + x[4]^2 - 2.0

        return f
    end

    # Solve test3 model

    x0 = [1.5, 1.5, 1.5, 1.5]
    l = [-Inf,-Inf,-Inf,-Inf]
    u = [Inf,Inf,Inf,Inf]

    soln_newton     = constrained_newton(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_lm         = constrained_levenberg_marquardt(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_lm_kyf     = constrained_levenberg_marquardt_kyf(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_lm_ar      = constrained_levenberg_marquardt_ar(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_dogleg     = constrained_dogleg_solver(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_dogleg_bmp = constrained_dogleg_bmp_solver(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_nk         = constrained_newton_krylov(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_jfnk       = constrained_jacobian_free_newton_krylov(test3,x0,l,u,xtol=1e-15,ftol=1e-15)

    soln_newton_sparse     = constrained_newton_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_sparse         = constrained_levenberg_marquardt_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_kyf_sparse     = constrained_levenberg_marquardt_kyf_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_ar_sparse      = constrained_levenberg_marquardt_ar_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_sparse     = constrained_dogleg_solver_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_bmp_sparse = constrained_dogleg_solver_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_nk_sparse         = constrained_newton_krylov_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)

    soln1 = nlboxsolve(test3,x0,l,u,method = :dogleg,xtol=1e-15,ftol=1e-15)
    soln2 = nlboxsolve(test3,x0,l,u,xtol=1e-15,ftol=1e-15)

    test_one   = maximum(abs,soln_newton.zero       - soln2.zero) < 1e-14
    test_two   = maximum(abs,soln_lm.zero           - soln2.zero) < 1e-14
    test_three = maximum(abs,soln_lm_kyf.zero       - soln2.zero) < 1e-14
    test_four  = maximum(abs,soln_lm_ar.zero        - soln2.zero) < 1e-14
    test_five  = maximum(abs,soln_dogleg.zero       - soln2.zero) < 1e-14
    test_six   = maximum(abs,soln_dogleg_bmp.zero   - soln2.zero) < 1e-14
    test_seven = maximum(abs,soln_nk.zero           - soln2.zero) < 1e-14
    test_eight = maximum(abs,soln_jfnk.zero         - soln2.zero) < 1e-14

    test_nine     = maximum(abs,soln_newton_sparse.zero       - soln2.zero) < 1e-14
    test_ten      = maximum(abs,soln_lm_sparse.zero           - soln2.zero) < 1e-14
    test_eleven   = maximum(abs,soln_lm_kyf_sparse.zero       - soln2.zero) < 1e-14
    test_twelve   = maximum(abs,soln_lm_ar_sparse.zero        - soln2.zero) < 1e-14
    test_thirteen = maximum(abs,soln_dogleg_sparse.zero       - soln2.zero) < 1e-14
    test_fourteen = maximum(abs,soln_dogleg_bmp_sparse.zero   - soln2.zero) < 1e-14
    test_fifteen  = maximum(abs,soln_nk_sparse.zero           - soln2.zero) < 1e-14

    return test_one, test_two, test_three, test_four, test_five, test_six, test_seven, test_eight, test_nine, test_ten, test_eleven,
           test_twelve, test_thirteen, test_fourteen, test_fifteen

end

function test_nlboxsolve_inplace()

    function test3(f,x)

        f[1] = 0.5*x[1] + x[2] - x[3] + x[1]*x[2] - 1.5
        f[2] = x[1] - 0.5*x[2] + x[3] + x[2]^2 - 2.5
        f[3] = -x[1] + x[2] + 0.5*x[3] + x[2]*x[3] - 1.5
        f[4] = x[1]^2 - x[2]^2 + x[3]^2 + x[4]^2 - 2.0

    end

    # Solve test3 model

    x0 = [1.5, 1.5, 1.5, 1.5]
    l = [-Inf,-Inf,-Inf,-Inf]
    u = [Inf,Inf,Inf,Inf]

    soln_newton     = constrained_newton(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_lm         = constrained_levenberg_marquardt(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_lm_kyf     = constrained_levenberg_marquardt_kyf(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_lm_ar      = constrained_levenberg_marquardt_ar(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_dogleg     = constrained_dogleg_solver(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_dogleg_bmp = constrained_dogleg_bmp_solver(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_nk         = constrained_newton_krylov(test3,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_jfnk       = constrained_jacobian_free_newton_krylov(test3,x0,l,u,xtol=1e-15,ftol=1e-15)

    soln_newton_sparse     = constrained_newton_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_sparse         = constrained_levenberg_marquardt_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_kyf_sparse     = constrained_levenberg_marquardt_kyf_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_ar_sparse      = constrained_levenberg_marquardt_ar_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_sparse     = constrained_dogleg_solver_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_bmp_sparse = constrained_dogleg_solver_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_nk_sparse         = constrained_newton_krylov_sparse(test3,x0,l,u,ftol=1e-15,xtol=1e-15)

    soln1 = nlboxsolve(test3,x0,l,u,method = :dogleg,xtol=1e-15,ftol=1e-15)
    soln2 = nlboxsolve(test3,x0,l,u,xtol=1e-15,ftol=1e-15)

    test_one   = maximum(abs,soln_newton.zero       - soln2.zero) < 1e-14
    test_two   = maximum(abs,soln_lm.zero           - soln2.zero) < 1e-14
    test_three = maximum(abs,soln_lm_kyf.zero       - soln2.zero) < 1e-14
    test_four  = maximum(abs,soln_lm_ar.zero        - soln2.zero) < 1e-14
    test_five  = maximum(abs,soln_dogleg.zero       - soln2.zero) < 1e-14
    test_six   = maximum(abs,soln_dogleg_bmp.zero   - soln2.zero) < 1e-14
    test_seven = maximum(abs,soln_nk.zero           - soln2.zero) < 1e-14
    test_eight = maximum(abs,soln_jfnk.zero         - soln2.zero) < 1e-14

    test_nine     = maximum(abs,soln_newton_sparse.zero       - soln2.zero) < 1e-14
    test_ten      = maximum(abs,soln_lm_sparse.zero           - soln2.zero) < 1e-14
    test_eleven   = maximum(abs,soln_lm_kyf_sparse.zero       - soln2.zero) < 1e-14
    test_twelve   = maximum(abs,soln_lm_ar_sparse.zero        - soln2.zero) < 1e-14
    test_thirteen = maximum(abs,soln_dogleg_sparse.zero       - soln2.zero) < 1e-14
    test_fourteen = maximum(abs,soln_dogleg_bmp_sparse.zero   - soln2.zero) < 1e-14
    test_fifteen  = maximum(abs,soln_nk_sparse.zero           - soln2.zero) < 1e-14

    return test_one, test_two, test_three, test_four, test_five, test_six, test_seven, test_eight, test_nine, test_ten, test_eleven,
           test_twelve, test_thirteen, test_fourteen, test_fifteen

end

test_nlboxsolve_outplace()
test_nlboxsolve_inplace()


text:
#!/usr/bin/env julia

printstyled("Starting QuickMag\n",bold=:true)
println()

printstyled("Updating dependencies...\n",bold=:true)
include(joinpath(".","src","updateDependencies.jl"))
println()

using JuliaDB
using Dates
using REPL.TerminalMenus

printstyled("Checking database status...\n",bold=:true)
ExistingDatabase=isfile(joinpath(".","HostFiles","WhiteList.jldb"));
if ExistingDatabase
	CurrentTime=Dates.now()
	WhiteListTable=WhiteListTable=load(joinpath(".","HostFiles","WhiteList.jldb"));
	UpdateTime=WhiteListTable[1].TimeStamp
	DayFrac=round(Millisecond(CurrentTime-UpdateTime)/Millisecond(Day(1)),digits=2); #Get number of days since last update
	if DayFrac<1.0
		printstyled("Database has recently been updated\n")			#Skip option to update if less than 24 hours
		println()
	else
		printstyled("Database last updated $DayFrac days ago\n")		#Offer update if >24 hours and databases exist
		options = ["Yes", "No"]
		menu = RadioMenu(options, pagesize=2)
		choice = request("Would you like to get fresh data?", menu)
		println()
		if choice==1
			printstyled("Rebuilding database of BOINC hosts:\n",bold=:true)
			include(joinpath(".","src","updateDatabase.jl"))
			println()
		else
			printstyled("Continuing with current database of BOINC hosts\n")
			println()
		end

	end

else											#Force update if there is no database files
	printstyled("No database found. Building new database of BOINC hosts:\n",bold=:true)
	include(joinpath(".","src","updateDatabase.jl"))
	println()
end
include(joinpath(".","src","QM_EvalFunctions.jl"))

printstyled("QuickMag is ready:\n",bold=:true)

while(true)

	RunOptions = ["Estimate CPU Performance", "Estimate GPU Performance", "BOINC CPU Survey", "BOINC GPU Survey", "Exit"]
	RunMenu = RadioMenu(RunOptions, pagesize=5)
	RunChoice = request("\nWhat would you like to do?", RunMenu)

	if RunChoice == 1		#Run QuickMagCPU(CPUid,Number)
		println()
		try
		println("Enter CPUid string to search for: (e.g. 'i7-6700 CPU', 'ryzen 9 3900X', ...)")
		LocCPUid=uppercase(readline());
		println("Max number of results per project? (Default 5)")
		LocNumber=readline();
		if isempty(LocNumber)
			LocNumber="5"
		end
		LocNumber = parse(Int32, LocNumber)

		println("Calculating...")
		println()
		QuickMagCPU(LocCPUid,LocNumber)
		println()

		catch LocError
			println("Failed: $LocError")
			println()
		end
	end
	if RunChoice == 2		#Run QuickMagGPU(GPUid,Number)
		println()
		try
		println("Enter GPUid string to search for: (e.g. 'RTX 3060', 'RX 5700 XT', ...)")
		LocGPUid=uppercase(readline());
		println("Max number of results per project? (Default 5)")
		LocNumber=readline();
		if isempty(LocNumber)
			LocNumber="5"
		end
		LocNumber = parse(Int32, LocNumber)

		println("Calculating...")
		println()
		QuickMagGPU(LocGPUid,LocNumber)
		println()

		catch LocError
			println("Failed: $LocError")
			println()
		end
	end

	if RunChoice == 3		#Run CPU survey
		println()
		try
		println("Max number of CPU models to return? (Default 100)")
		LocNumber=readline();
		if isempty(LocNumber)
			LocNumber="100"
		end
		LocNumber = parse(Int32, LocNumber)

		println("Calculating...")
		println()
		CPUsurvey(LocNumber)
		println()

		catch LocError
			println("Failed: $LocError")
			println()
		end
	end
	if RunChoice == 4		#Run GPU survey
		println()
		try
		println("Max number of GPU models to return? (Default 100)")
		LocNumber=readline();
		if isempty(LocNumber)
			LocNumber="100"
		end
		LocNumber = parse(Int32, LocNumber)

		println("Calculating...")
		println()
		GPUsurvey(LocNumber)
		println()

		catch LocError
			println("Failed: $LocError")
			println()
		end
	end
	if RunChoice == 5		#Exit Program
		break
	end
end

exit()


text:
try Pkg.clone("https://github.com/zsunberg/VDPTag.jl") end
using POMDPs
Pkg.build("VDPTag");        using VDPTag
Pkg.add("Reel");            using Reel
Pkg.add("Plots");           using Plots
Pkg.add("GR");
Pkg.add("ParticleFilters"); using ParticleFilters
Pkg.add("ProgressMeter");   using ProgressMeter
Pkg.add("POMDPToolbox");    using POMDPToolbox
Pkg.add("Distributions");   using Distributions

pomdp = VDPTagPOMDP()
filter = SIRParticleFilter(pomdp, 1000, rng=MersenneTwister(100))

hist = sim(pomdp, updater=filter, max_steps=100, rng=MersenneTwister(1)) do b
    # Policy: move towards predicted target position; if uncertainty area > 0.01, take measurement
    agent = first(particles(b)).agent
    target_particles = Array{Float64}(2, n_particles(b))
    for (i, s) in enumerate(particles(b))
        target_particles[:,i] = s.target
    end
    normal_dist = fit(MvNormal, target_particles)
    angle = action(ToNextML(mdp(pomdp)), TagState(agent, mean(normal_dist)))
    return TagAction(sqrt(det(cov(normal_dist))) > 0.01, angle)
end

gr()
frames = Frames(MIME("image/png"), fps=2)
@showprogress "Creating gif..." for i in 1:length(hist)
    push!(frames, plot(pomdp, view(hist, 1:i)))
end
write("out.gif", frames)