text:
function Base.findfirst(node::Leaf{T,N}, reg::Region{T,N}, id::Any) where {T,N}
    for (i, el) in enumerate(children(node))
        if isequal_rtree(el, reg, id)
            return (node, i)
        end
    end
    return nothing
end

function Base.findfirst(node::Branch{T,N,V}, reg::Region{T,N}, id::Any) where {T,N,V}
    for child in children(node)
        if contains(mbr(child), reg)
            res = findfirst(child, reg, id)
            if res !== nothing
                return res::Tuple{Leaf{T,N,V}, Int}
            end
        end
    end
    return nothing
end

"""
    findfirst(tree::RTree{T,N}, reg::Region{T,N}, [id]) where {T,N}

Find the element in the `tree` by its region (`reg`) and, optionally, its `id`.
The region (MBR) of the element and `reg` should match exactly.

Returns the tuple of `Leaf` and position of the element or `nothing`.
"""
Base.findfirst(tree::RTree{T,N}, reg::Region{T,N}, id::Any = nothing) where {T,N} =
    findfirst(tree.root, reg, id)

# FIXME: currently isempty() doesn't allow specifying how
#        to treat overlapping elements (inside or not), currently treated as outside
"""
    isempty(tree::RTree, region::Region)

Check if there are `tree` elements inside `region`.
"""
Base.isempty(tree::RTree{T,N}, region::Region{T,N}) where {T,N} =
    _isempty(tree.root, region)

function _isempty(node::Node, region::Region)
    isempty(node) && return true
    nodebr = mbr(node)
    if in(nodebr, region) # there are elements inside rect
        return false
    elseif intersects(nodebr, region) # there could be node elements inside region
        for child in children(node)
            if node isa Branch # should be optimized out at compile time
                _isempty(child, region) || return false
            elseif node isa Leaf
                in(mbr(child), region) && return false
            end
        end
    end
    return true
end

# the RTreeIterator state
struct RTreeIteratorState{T,N,V}
    leaf::Leaf{T,N,V}       # current leaf node
    indices::Vector{Int}    # indices of the nodes (in their parents) in the current subtree
end

# get the current data element pointed by `RTreeIteratorState`
Base.get(state::RTreeIteratorState) = @inbounds(state.leaf[state.indices[1]])

# iterate all R-tree data elements
function Base.iterate(tree::RTree)
    isempty(tree) && return nothing
    node = tree.root
    indices = fill(1, height(tree))
    # get the first leaf
    while level(node) > 1
        node = node[1]
    end
    state = RTreeIteratorState(node, indices)
    return get(state), state # first element of the first leaf
end

function Base.iterate(tree::RTree, state::RTreeIteratorState)
    @inbounds if state.indices[1] < length(state.leaf) # fast branch: next data element in the same leaf
        state.indices[1] += 1
        return get(state), state
    end
    # leaf iterations is done, go up until the first non-visited branch
    node = state.leaf
    while state.indices[level(node)] >= length(node)
        hasparent(node) || return nothing # returned to root, iteration finished
        node = parent(node)
    end
    # go down into the first leaf of the new subtree
    ix = state.indices[level(node)] += 1
    node = node[ix]
    @inbounds while true
        state.indices[level(node)] = 1
        level(node) == 1 && break
        node = node[1]
    end
    new_state = RTreeIteratorState(node, state.indices)
    return get(new_state), new_state
end

# iterates R-tree data elements matching `Q` query w.r.t `region`
struct RTreeRegionQueryIterator{T,N,V,Q,TT,R} <: SpatialQueryIterator{T,N,V,Q}
    tree::TT
    region::R

    function RTreeRegionQueryIterator{T,N}(kind::QueryKind, tree::TT, region::R) where
            {T, N, V, TT <: RTree{T,N,V}, R <: Region{T,N}}
        new{T,N,V,kind,TT,R}(tree, region)
    end
end

# the RTreeRegionQueryIterator state
struct RTreeQueryIteratorState{T,N,V}
    leaf::Leaf{T,N,V}       # current leaf node
    indices::Vector{Int}    # indices of the nodes (in their parents) in the current subtree
    needtests::BitVector    # whether children MBRs should be tested or not (because they automatically satisfy query)
end

# get the current data element pointed by `RTreeIteratorState`
Base.get(state::RTreeQueryIteratorState) = @inbounds(state.leaf[state.indices[1]])

function Base.iterate(iter::RTreeRegionQueryIterator)
    isempty(iter.tree) && return nothing
    # no data or doesn't intersect at all
    root_match = should_visit(iter.tree.root, iter)
    root_match == QueryNoMatch && return nothing
    return _iterate(iter, iter.tree.root, fill(1, height(iter.tree)),
                    root_match == QueryMatchComplete ? falses(height(iter.tree)) : trues(height(iter.tree)))
end

function Base.iterate(iter::RTreeRegionQueryIterator,
                      state::RTreeQueryIteratorState)
    @inbounds ix = state.indices[1] = _nextchild(state.leaf, state.indices[1] + 1, state.needtests[1], iter)[1]
    if ix <= length(state.leaf) # fast branch: next data element in the same leaf
        return get(state), state
    else
        return _iterate(iter, state.leaf, state.indices, state.needtests)
    end
end

"""
    contained_in(index::SpatialIndex, region::Region)

Get iterator for `index` elements contained in `region`.
"""
contained_in(tree::RTree{T,N}, region::Region{T,N}) where {T,N} =
    RTreeRegionQueryIterator{T,N}(QueryContainedIn, tree, region)

"""
    intersects_with(index::SpatialIndex, region::Region)

Get iterator for `index` elements intersecting with `region`.
"""
intersects_with(tree::RTree{T,N}, region::Region{T,N}) where {T,N} =
    RTreeRegionQueryIterator{T,N}(QueryIntersectsWith, tree, region)

# whether the R-tree node/data element should be visited (i.e. its children examined)
# by the region iterator
function should_visit(node::Node, iter::RTreeRegionQueryIterator)
    kind = querykind(iter)
    if kind == QueryContainedIn || kind == QueryIntersectsWith
        contains(iter.region, mbr(node)) && return QueryMatchComplete # node (and all its children) fully contained in query region
        intersects(iter.region, mbr(node)) && return QueryMatchPartial # node (and maybe its children) intersects query region
        return QueryNoMatch
    else
        throw(ArgumentError("Unknown spatial query kind: $kind"))
    end
end

should_visit(el::Any, iter::RTreeRegionQueryIterator) =
    ((querykind(iter) == QueryContainedIn) && contains(iter.region, mbr(el))) ||
    ((querykind(iter) == QueryIntersectsWith) && intersects(iter.region, mbr(el))) ?
    QueryMatchComplete : QueryNoMatch
    # FIXME update for NotContainedIn etc

# get the index of the first child of `node` starting from `pos` (including)
# that satifies `iter` query (or length(node) + 1 if not found)
@inline function _nextchild(node::Node, pos::Integer, needtests::Bool, iter::RTreeRegionQueryIterator)
    res = needtests ? QueryNoMatch : QueryMatchComplete # if tests not needed, all node subtrees are considered fully matching
    while pos <= length(node) && needtests &&
            ((res = should_visit(@inbounds(node[pos]), iter)) == QueryNoMatch)
        pos += 1
    end
    return pos, res
end

# do depth-first search starting from the `node` subtree and return the
# `RTreeIteratorState` for the first leaf that satisfies `iter` query or
# `nothing` if no such leaf in the R-tree.
# The method modifies `indicies` array and uses it for the returned iteration state
function _iterate(iter::RTreeRegionQueryIterator, node::Node,
                  indices::AbstractVector{Int}, needtests::AbstractVector{Bool})
    #@debug"_iterate(): enter lev=$(level(node)) indices=$indices"
    @assert length(indices) == height(iter.tree)
    ix = @inbounds(indices[level(node)])
    while true
        ix_new, queryres = _nextchild(node, ix, needtests[level(node)], iter)
        #@debug "node=$(Int(Base.pointer_from_objref(node))) lev=$(level(node)) ix_new=$ix_new"
        if ix_new > length(node) # all node subtrees visited, go up one level
            while ix_new > length(node)
                if !hasparent(node)
                    #@debug "_iterate(): finished lev=$(level(node)) indices=$indices ix_new=$ix_new"
                    return nothing # returned to root, iteration finished
                end
                #@debug "_iterate(): up lev=$(level(node)) indices=$indices ix_new=$ix_new"
                node = parent(node)
                @inbounds ix_new = indices[level(node)] += 1 # next subtree
            end
            ix = ix_new
            #@debug "_iterate(): next subtree lev=$(level(node)) indices=$indices ix_new=$ix_new"
        else # subtree found
            ix_new > ix && @inbounds(indices[level(node)] = ix_new)
            if node isa Branch
                # go down into the first child
                node = node[ix_new]
                indices[level(node)] = ix = 1
                needtests[level(node)] = queryres != QueryMatchComplete
                #@debug "_iterate(): down lev=$(level(node)) indices=$indices"
            else # Leaf
                #@debug "_iterate(): return lev=$(level(node)) indices=$indices"
                state = RTreeQueryIteratorState(node, indices, needtests)
                return get(state), state
            end
        end
    end
end


text:
"Conditioning"
module Cond
export cond, condf, condp, |ₚ
using ..Omega
using DocStringExtensions

export observe

include("pred.jl")        # Conditioning on general predicates
include("observe.jl")     # Condtion on observations
include("sugar.jl")       # Syntactic sugar

end

ext: jl


text:
using eFEMpart, JLD

using LinearAlgebra # for condition number calculation

function main()
  function computeNorms(N)
    start = time()
    # load mesh
    mesh = squareMeshFluid([-2,2,-1,1],N)

    μ0     = 2.3
    μ(x,y) = μ0 + exp(x*y)
    xm = [i.x for i in mesh.xy]
    ym = [i.y for i in mesh.xy]
    μarr = [μ(xm[i],ym[i]) for i=1:length(mesh.xy)]
    param  = FluidParam(μarr)

    OperatorType = :Stokes2D

    dNodes = Dirichlet(:left,:top,:bottom,:right)
    Nodes = [dNodes]

    dUBCf = Dirichlet((x,y) -> -x^4*y^2)
    dVBCf = Dirichlet((x,y) ->  4*x^3*y^3/3)

    dudx(x,y) = -4*x^3*y^2; dudy(x,y) = -2*x^4*y
    dvdx(x,y) = 4*x^2*y^3;  dvdy(x,y) = 4*x^3*y^2
    cross(x,y) = dudy(x,y) + dvdx(x,y)

    d2dudx(x,y) = -12*x^2*y^2; d2dudy(x,y) = -2*x^4
    d2dvdx(x,y) = 8*x*y^3;     d2dvdy(x,y) = 8*x^3*y
    lapu(x,y) = d2dudx(x,y) + d2dudy(x,y)
    lapv(x,y) = d2dvdx(x,y) + d2dvdy(x,y)

    dpdx(x,y) = 3*x^2*y^3; dpdy(x,y) = 3*x^3*y^2

    dμdx(x,y) = y*exp(x*y)
    dμdy(x,y) = x*exp(x*y)

    Fx(x,y) = -μ(x,y)*lapu(x,y) - (2*dμdx(x,y)*dudx(x,y) + dμdy(x,y)*cross(x,y)) + dpdx(x,y)
    Fy(x,y) = -μ(x,y)*lapv(x,y) - (2*dμdy(x,y)*dvdy(x,y) + dμdx(x,y)*cross(x,y)) + dpdy(x,y)

    ffx   = Forcing(Fx)
    ffy   = Forcing(Fy)
    bcfun = [dUBCf,dVBCf,ffx,ffy]

    prob = Problem(mesh,Nodes,bcfun,OperatorType)
    sol = solve(prob,mesh,param)

    elapsed = time() - start

    # compute condition number of operator matrix
    #LinOp = GenerateSystem(mesh,prob,param)
    #ApplyBC!(LinOp,mesh,prob,param,OperatorType)
    #κ = cond(Array(LinOp.Op),2)
    κ = 0.5

    # compute mesh h
    h = hCalc(mesh)

    # generate exact solution
    Uexact(x,y) = dUBCf.f(x,y)
    Vexact(x,y) = dVBCf.f(x,y)
    Pexact(x,y) = x^3*y^3
    UexactArr = [Uexact.(mesh.xy[i].x,mesh.xy[i].y) for i=1:length(mesh.xy)]
    VexactArr = [Vexact.(mesh.xy[i].x,mesh.xy[i].y) for i=1:length(mesh.xy)]
    PexactArr = [Pexact.(mesh.xyp[i].x,mesh.xyp[i].y) for i=1:length(mesh.xyp)]

    # compute difference
    velErr = sqrt.((sol.u - UexactArr).^2 + (sol.v -VexactArr).^2)
    presErr = sol.p - PexactArr

    # compute Li norm
    L1v   = DomainNorm(mesh.xy,mesh.cm,velErr;normID="1")
    L2v   = DomainNorm(mesh.xy,mesh.cm,velErr;normID="2")
    Linfv = DomainNorm(mesh.xy,mesh.cm,velErr;normID="Inf")
    L1p   = DomainNorm(mesh.xyp,mesh.cmp,presErr;normID="1")
    L2p   = DomainNorm(mesh.xyp,mesh.cmp,presErr;normID="2")
    Linfp = DomainNorm(mesh.xyp,mesh.cmp,presErr;normID="Inf")

    # plot solution
    sD = ScalarData(sol.p)
    sN = ScalarNames("pressure")
    vD = VectorData([sol.u,sol.v])
    vN = VectorNames("velocity")
    fn = Path("solution_var")
    vtksave(mesh,sD,sN,vD,vN,fn)

    # plot solution
    sD = ScalarData(PexactArr)
    sN = ScalarNames("pressure")
    vD = VectorData([UexactArr,VexactArr])
    vN = VectorNames("velocity")
    fn = Path("exact_var")
    vtksave(mesh,sD,sN,vD,vN,fn)

    return κ,h,L1v,L2v,Linfv,L1p,L2p,Linfp,elapsed
  end

  Narr     = [4,4,8,16,32]#,64,96]
  N        = length(Narr)
  κarr     = zeros(N)
  harr     = zeros(N)
  L1arrV   = zeros(N)
  L2arrV   = zeros(N)
  LInfarrV = zeros(N)
  L1arrP   = zeros(N)
  L2arrP   = zeros(N)
  LInfarrP = zeros(N)
  timearr  = zeros(N)

  for i=1:N
    n = Narr[i]
    κarr[i],harr[i],L1arrV[i],L2arrV[i],LInfarrV[i],
      L1arrP[i],L2arrP[i],LInfarrP[i],timearr[i] =
      computeNorms(n)
      println("completed N=$(n)")
  end

  κarr     = κarr[2:end]
  harr     = harr[2:end]
  L1arrV   = L1arrV[2:end]
  L2arrV   = L2arrV[2:end]
  LInfarrV = LInfarrV[2:end]
  L1arrP   = L1arrP[2:end]
  L2arrP   = L2arrP[2:end]
  LInfarrP = LInfarrP[2:end]
  timearr  = timearr[2:end]

  # export output to *.jld file
  run(`mkdir -p data`)

  save("data/TEMP_stokes2D_var_validation.jld", "harr",harr,
                                          "VL1arr", L1arrV,
                                          "VL2arr", L2arrV,
                                          "VLInfarr",LInfarrV,
                                          "PL1arr", L1arrP,
                                          "PL2arr", L2arrP,
                                          "PLInfarr",LInfarrP,
                                          "timearr",timearr,
                                          "κarr",κarr)

  nothing
end

main()

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
Base.promote_rule(::Type{VaxInt16}, ::Type{T}) where T <: Union{Int8,Int16} = Int16
Base.promote_rule(::Type{VaxInt16}, ::Type{T}) where T <: Union{Int32,Int64,Int128} = T
Base.promote_rule(::Type{VaxInt16}, ::Type{T}) where T <: IEEEFloat = T

Base.promote_rule(::Type{VaxInt32}, ::Type{T}) where T <: Union{Int8,Int16,Int32,VaxInt16} = Int32
Base.promote_rule(::Type{VaxInt32}, ::Type{T}) where T <: Union{Int64,Int128,IEEEFloat} = T
Base.promote_rule(::Type{VaxInt32}, ::Type{T}) where T <: IEEEFloat = T

Base.promote_rule(::Type{VaxFloatF}, ::Type{T}) where T <: Union{Int8,Int16,Int32,Int64,Int128} = Float32
Base.promote_rule(::Type{VaxFloatF}, ::Type{T}) where T <: Union{Float16,Float32} = Float32
Base.promote_rule(::Type{VaxFloatF}, ::Type{T}) where T <: VaxInt = Float32
Base.promote_rule(::Type{VaxFloatF}, ::Type{Float64}) = Float64

Base.promote_rule(::Type{VaxFloatD}, ::Type{T}) where T <: Union{Int8,Int16,Int32,Int64,Int128} = Float64
Base.promote_rule(::Type{VaxFloatD}, ::Type{T}) where T <: Union{VaxFloatF,VaxFloatG} = Float64
Base.promote_rule(::Type{VaxFloatD}, ::Type{T}) where T <: VaxInt = Float64
Base.promote_rule(::Type{VaxFloatD}, ::Type{T}) where T <: IEEEFloat = Float64

Base.promote_rule(::Type{VaxFloatG}, ::Type{T}) where T <: Union{Int8,Int16,Int32,Int64,Int128} = Float64
Base.promote_rule(::Type{VaxFloatG}, ::Type{T}) where T <: Union{VaxFloatF,VaxFloatG} = Float64
Base.promote_rule(::Type{VaxFloatG}, ::Type{T}) where T <: VaxInt = Float64
Base.promote_rule(::Type{VaxFloatG}, ::Type{T}) where T <: IEEEFloat = Float64

Base.promote_rule(::Type{BigFloat}, ::Type{<:VaxFloat}) = BigFloat
Base.promote_rule(::Type{BigFloat}, ::Type{<:VaxInt}) = BigFloat

Base.promote_rule(::Type{BigInt}, ::Type{<:VaxFloat}) = BigFloat
Base.promote_rule(::Type{BigInt}, ::Type{<:VaxInt}) = BigFloat

Base.promote_type(::Type{VaxInt16}, ::Type{VaxInt16}) = Int16
Base.promote_type(::Type{VaxInt32}, ::Type{VaxInt32}) = Int32

Base.promote_type(::Type{VaxFloatF}, ::Type{VaxFloatF}) = Float32
Base.promote_type(::Type{VaxFloatD}, ::Type{VaxFloatD}) = Float64
Base.promote_type(::Type{VaxFloatG}, ::Type{VaxFloatG}) = Float64


text:
module NPFinancial

export fv, pmt, nper, ipmt, ppmt, pv, rate, irr, npv, mirr, eir

using Polynomials: Poly, roots

# Internal functions

const _when_to_num = Dict(:begin => 1, :end => 0)

function _convert_when(when)
    return _when_to_num[when]
end

function _rbl(rate, per, pmt, pv, when)
    return fv(rate, (per - 1), pmt, pv, when)
end

"""
    fv(rate::Real, nper::Integer, pmt::Real, pv::Real, when = :end)

Compute the future value given the present value `pv`, an interest rate `rate` that
is compounded once per period, over `nper` number of periods. A fixed payment `pmt`
may be specified in the `when` argument, which is paid either at the beginning of
each period (`:begin`) or `:end` of the period.

# Examples
```
julia> fv(0.05, 1, 0, -100)
105.0

julia> fv(0.05, 2, 0, -100)
110.25
```
"""
function fv(rate::Real, nper::Integer, pmt::Real, pv::Real, when = :end)
    when = _convert_when(when)
    temp = (1+rate)^nper
    fact = rate == 0 ? nper : (1 + rate*when)*(temp - 1)/rate
    return -(pv * temp + pmt * fact)
end

"""
    pmt(rate::Real, nper::Integer, pv::Real, fv = 0.0, when = :end)

Compute the payment given the present value `pv`, an interest rate `rate` that
is compounded once per period, over `nper` number of periods, such that at the
end of the last period the value becomes `fv`. The payment is expected to be paid
at the beginning of each period (`:begin`) or `:end` of the period, as specified
in the `when` argument.

# Examples

Let's say I have need to repay a mortage loan with amount 300000 (present value)
with monthly interst rate of 4.25% over the next 30 years.  What would be my
monthly payment?
```
julia> pmt(0.0425/12, 30*12, 300000)
-1475.8196732384283
```
"""
function pmt(rate::Real, nper::Integer, pv::Real, fv = 0.0, when = :end)
    when = _convert_when(when)
    temp = (1 + rate).^nper
    mask = (rate == 0)
    masked_rate = mask ? 1 : rate
    fact = mask != 0 ? nper : (1 + masked_rate * when) * (temp - 1)/masked_rate
    return -(fv + pv * temp) / fact
end

"""
    nper(rate::Real, pmt::Real, pv::Real, fv = 0.0, when = :end)

Compute how many periods the present value `pv` may accrue/repaid till the future
value `fv` given a specific interest rate `rate` and a fixed payment `pmt`.
The payment is expected to be paid
at the beginning of each period (`:begin`) or `:end` of the period, as specified
in the `when` argument.
"""
function nper(rate::Real, pmt::Real, pv::Real, fv = 0.0, when = :end)
    when = _convert_when(when)
    if rate == 0
        return -(fv + pv) / pmt
    else
        z = pmt * (1 + rate * when) / rate
        return log((-fv+z) / (pv+z))/log(1+rate)
    end
end

"""
    ipmt(rate::Real, per::Integer, nper::Integer, pv::Real, fv = 0.0, when = :end)

Compute the interest component of the periodic payment.  This useful for any
loan that has a repayment schedule.

# Examples

Let's say I have need to repay a mortage loan with amount 300000 (present value)
with monthly interst rate of 4.25% over the next 30 years.  What would be the
interest component of my monthly payment?  Initially, the interest component is
a large part of the payment but towards the end, it would be a smaller portion
as illustrated below for periods 1, 2, 3, 358, 359, and 360:

```
julia> pmt(0.0425/12, 30*12, 300000)
-1475.8196732384283

julia> ipmt.(0.0425/12, [1,2,3,358,359,360], 30*12, 300000)
6-element Array{Float64,1}:
 -1062.5
 -1061.04
 -1059.57
   -15.5702
   -10.3984
    -5.20841
```
"""
function ipmt(rate::Real, per::Integer, nper::Integer, pv::Real, fv = 0.0, when = :end)
    #when = _convert_when(when)
    #rate, per, nper, pv, fv, when = np.broadcast_arrays(rate, per, nper, pv, fv, when)
    total_pmt = pmt(rate, nper, pv, fv, when)
    ipmt = _rbl(rate, per, total_pmt, pv, when)*rate
    ipmt = when == 1 ? ipmt/(1 + rate) : ipmt
    ipmt = (when == 1) && (per == 1) ? 0 : ipmt
    return ipmt
end

"""
    ppmt(rate, per, nper, pv, fv = 0.0, when = :end)

Compute the principal component of the periodic payment.  This useful for any
loan that has a repayment schedule.

# Examples

```
julia> pmt(0.0425/12, 30*12, 300000)
-1475.8196732384283

julia> ppmt.(0.0425/12, [1,2,3,358,359,360], 30*12, 300000)
6-element Array{Float64,1}:
  -413.32
  -414.784
  -416.253
 -1460.25
 -1465.42
 -1470.61
```
"""
function ppmt(rate::Real, per::Integer, nper::Integer, pv::Real, fv = 0.0, when = :end)
    total = pmt(rate, nper, pv, fv, when)
    return total - ipmt(rate, per, nper, pv, fv, when)
end

"""
    pv(rate::Real, nper::Integer, pmt::Real, fv = 0.0, when = :end)

Compute the present value given the future value `fv`, an interest rate `rate` and
a fixed periodic payment `pmt` over a number of periods `nper`.  The payment is
expected to be paid at the beginning of each period (`:begin`) or `:end` of the period,
as specified in the `when` argument.
"""
function pv(rate::Real, nper::Integer, pmt::Real, fv = 0.0, when = :end)
    when = _convert_when(when)
    temp = (1+rate)^nper
    fact = rate == 0 ? nper : (1+rate*when)*(temp-1)/rate
    return -(fv + pmt*fact)/temp
end

function _g_div_gp(r, n, p, x, y, w)
    t1 = (r+1)^n
    t2 = (r+1)^(n-1)
    return ((y + t1*x + p*(t1 - 1)*(r*w + 1)/r) /
                (n*t2*x - p*(t1 - 1)*(r*w + 1)/(r^2) + n*p*t2*(r*w + 1)/r +
                p*(t1 - 1)*w/r))
end

"""
    rate(nper::Integer, pmt::Real, pv::Real, fv::Real, when=:end, guess=0.1, tol=1e-6, maxiter=100)

Compute interest rate given present value `pv`, future value `fv`, and fixed periodic
payment `pmt` over a number of periods `nper`.

This implementation uses Newton's iteration until the change is less than 1e-6.
Newton's rule is

``r_{n+1} = r_{n} - \\frac{g(r_n)}{g'(r_n)}``

where
* g(r) is the formula
* g'(r) is the derivative with respect to r.

# Examples
```
julia> rate(1, 0, -100, 101)
0.010000000000000155
```
"""
function rate(nper::Integer, pmt::Real, pv::Real, fv::Real,
              when=:end, guess=0.1, tol=1e-6, maxiter=100)
    when = _convert_when(when)
    rn = guess
    iterator = 0
    close = false
    while (iterator < maxiter) && !close
        rnp1 = rn - _g_div_gp(rn, nper, pmt, pv, fv, when)
        diff = abs(rnp1-rn)
        close = diff < tol
        iterator += 1
        rn = rnp1
    end
    if !close
        # Return nan's in array of the same shape as rn
        return NaN + rn
    else
        return rn
    end
end

"""
     irr(values::Real)

Calculate internal rate of return given an array of cash flow `values`
(nearest one first)

# Examples
```
julia> irr([-100, 101])
0.010000000000000009
```
"""
function irr(values)
    res = roots(Poly(values))
    mask = map(x -> isa(x, Complex) ? x.im == 0 && x.re > 0 : x > 0, res)
    if !any(mask)
        return NaN
    end
    res = map(x -> isa(x, Complex) ? x.re : x, res[mask])
    rate = @. 1 / res - 1
    # TODO not sure why we can just do min(abs(rate))... seems convoluted
    rate = rate[findmin(abs.(rate))[2]]
    return rate
end

"""
    npv(rate::Real, values::AbstractVector{<:Real})

Compute the Net Present Value (NPV) of a cash flow series `values` given
an internal rate of return `rate`.

The (fixed) time interval between cash flow "events" must be the same as that for
which `rate` is given (i.e., if `rate` is per year, then precisely
a year is understood to elapse between each cash flow event).  By
convention, investments or "deposits" are negative, income or
"withdrawals" are positive; `values` must begin with the initial
investment, thus `values[1]` will typically be negative.

# Examples
```
julia> npv(0.281,[-100, 39, 59, 55, 20])
-0.00847859163845488
```

# Reference
* L. J. Gitman, "Principles of Managerial Finance, Brief," 3rd ed.,
    Addison-Wesley, 2003, pg. 346.
"""
function npv(rate::Real, values::AbstractVector{<:Real})
    return sum(values ./ (1+rate).^(0:length(values)-1))
end

"""
    mirr(values::AbstractVector{<:Real}, finance_rate::Real, reinvest_rate::Real)

Compute the modified internal rate of return (MIRR) given a series of cash flows,
a `finance_rate` (interest rate paid on the cash flows) and `reinvest_rate` (interest
rate received on the cash flows upon reinvestment).
"""
function mirr(values::AbstractVector{<:Real}, finance_rate::Real, reinvest_rate::Real)
    n = length(values)
    pos = values .> 0
    neg = values .< 0
    if !(any(pos) && any(neg))
        return NaN
    end
    numer = abs(npv(reinvest_rate, values .* pos))
    denom = abs(npv(finance_rate, values .* neg))
    return (numer/denom)^(1/(n - 1))*(1 + reinvest_rate) - 1
end

"""
    eir(rate::Real, nper::Real)

Computes the effective interest rate `ieff` given
* `rate` is the nominal (annual) rate
* `nper` is the number of compounding periods
"""

function eir(rate::Real, nper::Real)
    return (((1 + (rate / nper))^nper)-1)
end

end


text:
"""
Geometry struct containing parameters and arrays describing the Gaussian grid
and the vertical levels. NF is the number format used for the precomputed constants.
"""
struct Geometry{NF<:AbstractFloat}      # NF: Number format

    # GRID-POINT SPACE
    nlon::Int           # Number of longitudes
    nlat::Int           # Number of latitudes
    nlev::Int           # Number of vertical levels
    nlat_half::Int      # Number of latitudes in one hemisphere
    nlon_half::Int      # Half the number of longitudes

    dlon::NF            # grid spacing in longitude
    dlat::NF            # average grid spacing in latitude
    lon::Array{NF,1}    # array of longitudes
    lat::Array{NF,1}    # array of latitudes

    # VERTICAL SIGMA COORDINATE σ = p/p0 (fraction of surface pressure)
    n_stratosphere_levels::Int      # number of upper levels for stratosphere
    σ_levels_half::Array{NF,1}      # σ at half levels
    σ_levels_full::Array{NF,1}      # σ at full levels
    σ_levels_thick::Array{NF,1}     # σ level thicknesses
    σ_levels_half⁻¹_2::Array{NF,1}  # 1/(2σ_levels_full)       #TODO rename?
    σ_f::Array{NF,1}                # akap/(2σ_levels_thick)   #TODO rename?

    # SINES AND COSINES OF LATITUDE
    sinlat::Array{NF,1}         # sin of latitudes
    coslat::Array{NF,1}         # cos of latitudes
    sinlat_NH::Array{NF,1}      # only northern hemisphere
    coslat_NH::Array{NF,1}      # only northern hemisphere
    radang::Array{NF,1}         # radians of latitudes TODO rename to radlat?
    cosg::Array{NF,1}           # this should be sinlat?!
    cosg⁻¹::Array{NF,1}         # rename to sinlat⁻¹?
    cosg⁻²::Array{NF,1}         # rename to sinlat⁻²?

    # CORIOLIS FREQUENCY
    f::Array{NF,1}              # = 2Ω*sin(lat)

    # GEOPOTENTIAL CALCULATION WORK ARRAYS
    xgeop1::Array{NF,1}                  # ?
    xgeop2::Array{NF,1}                  # ?
    lapserate_correction::Array{NF,1}    # ?
end

"""
Defines the geometry.
"""
function Geometry(P::Parameters)

    @unpack nlon,nlat,nlev,trunc = P
    @unpack R,Ω,akap = P
    @unpack n_stratosphere_levels = P

    nlat_half = nlat ÷ 2
    nlon_half = nlon ÷ 2

    # GRID SPACE ARRAYS GAUSSIAN GRID - lon is equi-spaced, lat is not!
    dlon = 360 / nlon                       # grid spacing in longitude
    dlat = 180 / nlat                       # average grid spacing in latitude
    lon  = Array(0:dlon:360-dlon)           # array of longitudes
    # array of latitudes (North to South) corresponding to the zeros
    # of the (unassociated) legendre polynomial order nlat
    lat  = reverse(asind.(gausslegendre(nlat)[1]))

    # VERTICAL SIGMA COORDINATE
    # σ = p/p0 (fraction of surface pressure)
    # sorted such that σ_levels_half[end] is at the planetary boundary
    σ_levels_half = vertical_coordinates(P)
    σ_levels_full = 0.5*(σ_levels_half[2:end] + σ_levels_half[1:end-1])
    σ_levels_thick = σ_levels_half[2:end] - σ_levels_half[1:end-1]
    σ_levels_half⁻¹_2 = 1 ./ (2σ_levels_thick)
    σ_f = akap ./ (2σ_levels_full)

    # SINES AND COSINES OF LATITUDE
    sinlat = sind.(lat)
    coslat = cosd.(lat)
    sinlat_NH = sinlat[1:nlat_half] # sinlat only for northern hemisphere = NH
    coslat_NH = coslat[1:nlat_half]
    radang = asin.(sinlat)
    cosg   = sinlat                 # inconsistent here due to the sin/cos swap
    cosg⁻¹ = 1 ./ cosg
    cosg⁻² = 1 ./ cosg.^2

    # CORIOLIS FREQUENCY
    f = 2Ω*sinlat

    # GEOPOTENTIAL coefficients to calculate geopotential (TODO reference)
    xgeop1 = zeros(nlev)
    xgeop2 = zeros(nlev)
    for k in 1:nlev
        xgeop1[k] = R*log(σ_levels_half[k+1]/σ_levels_half[k])
        if k != nlev
            xgeop2[k+1] = R*log(σ_levels_full[k+1]/σ_levels_half[k+1])
        end
    end

    # LAPSE RATE correction (TODO reference)
    lapserate_correction = zeros(nlev-2)
    for k in 2:nlev-1
        lapserate_correction[k-1] = 0.5*xgeop1[k]*
                    log(σ_levels_half[k+1]/σ_levels_full[k]) / log(σ_levels_full[k+1]/σ_levels_full[k-1])
    end

    # conversion to number format NF happens here
    Geometry{P.NF}( nlon,nlat,nlev,nlat_half,nlon_half,
                    dlon,dlat,lon,lat,
                    n_stratosphere_levels,
                    σ_levels_half,σ_levels_full,σ_levels_thick,σ_levels_half⁻¹_2,σ_f,
                    sinlat,coslat,sinlat_NH,coslat_NH,radang,
                    cosg,cosg⁻¹,cosg⁻²,f,xgeop1,xgeop2,lapserate_correction)
end

"""Vertical sigma coordinates defined by their nlev+1 half levels `σ_levels_half`. Sigma coordinates are
fraction of surface pressure (p/p0) and are sorted from top (stratosphere) to bottom (surface).
The first half level is at 0 the last at 1. Evaluate a generalised logistic function with
coefficients in `P` for the distribution of values in between. Default coefficients follow
the L31 configuration historically used at ECMWF."""
function vertical_coordinates(P::Parameters)
    @unpack nlev,GLcoefs = P

    halflevels_normalised = range(0,1,nlev+1)   # normalised = level/nlev
    σ_levels_half = generalised_logistic(halflevels_normalised,GLcoefs)
    σ_levels_half[1] = 0           # topmost half-level is at 0 pressure
    σ_levels_half[end] = 1         # lowermost half-level is at p=p_surface

    @assert isincreasing(σ_levels_half) "Vertical coordinates are not increasing."
    return σ_levels_half
end

"""Generalised logistic function based on the coefficients in `coefs`."""
function generalised_logistic(x,coefs::GenLogisticCoefs)
    @unpack A,K,C,Q,B,M,ν = coefs
    return @. A + (K-A)/(C+Q*exp(-B*(x-M)))^(1/ν)
end


text:
# This file is a part of TextClassification.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export KnnClassifierConfig, KnnClassifierConfigSpace, KnnClassifier

using StatsBase: counts
import StatsBase: predict

@with_kw struct KnnClassifierConfig
    k::Int32 = 1
    keeptop::Float32 = 1.0
end

create(config::KnnClassifierConfig, train_X, train_y, dim) = KnnClassifier(config, train_X, train_y)

@with_kw struct KnnClassifierConfigSpace <: AbstractSolutionSpace
    k=[1, 5]
    keeptop=[1.0]
    scale_k = (lower=1, s=1.5, upper=100)
    scale_keeptop = (lower=0.1, s=1.5, upper=1.0)
end

Base.eltype(::KnnClassifierConfigSpace) = KnnClassifierConfig

function Base.rand(space::KnnClassifierConfigSpace)
    KnnClassifierConfig(rand(space.k), rand(space.keeptop))
end

function combine(a::KnnClassifierConfig, b::KnnClassifierConfig)
    k = div(a.k + b.k, 2)
    KnnClassifierConfig(k, a.keeptop)
end

function mutate(space::AbstractSolutionSpace, a::KnnClassifierConfig, iter)
    k = space.scale_k === nothing ? a.k : SearchModels.scale(a.k; space.scale_k...)
    keeptop = space.scale_keeptop === nothing ? a.keeptop : SearchModels.scale(a.keeptop; space.scale_keeptop...)
    KnnClassifierConfig(k, keeptop)
end

struct KnnClassifier
    config::KnnClassifierConfig
    index::InvIndex
    labels::CategoricalArray
end

function KnnClassifier(config::KnnClassifierConfig, X, y::CategoricalArray)
    invindex = InvIndex(X)
    invindex = config.keeptop < typemax(Int) ? prune(invindex, config.keeptop) : invindex
    KnnClassifier(config, invindex, y)
end

function predict(cls::KnnClassifier, vec::SVEC)
    res = search(cls.index, vec, cls.config.k)
    knn_most_frequent_label(cls, res)
end

function knn_most_frequent_label(knn::KnnClassifier, res::KnnResult)
    c = counts([knn.labels.refs[id] for (id, dist) in res], 1:length(levels(knn.labels)))
    findmax(c)[end]
end


text:
start,rules = split(read("input.txt", String), "\n\n")
re = r"(\w+) -> (\w)"
rules = Dict(k => v for (k,v) = ([match(re, rule) for rule = eachline(IOBuffer(rules))]))
freq(X) = Dict([(x, count(==(x), X)) for x = unique(X)])
repeatf(f, x, n) = n > 1 ? f(repeatf(f, x, n-1)) : f(x)
reaction(str,c) = str * (haskey(rules, string(str[end], c)) ? rules[string(str[end], c)] : "") * c

str = repeatf(s -> foldl(reaction, s), start, 10)
elements = sort(collect(freq(str)), by=x->x[2])
println(last(elements).second - first(elements).second)


text:
module Models

using ..Scruff
using ..Scruff.Operators
using ..Scruff.SFuncs

import ..Scruff: make_initial, make_transition

include("models/instantmodel.jl")
include("models/timelessinstantmodel.jl")
include("models/simplemodel.jl")
include("models/fixedtimemodel.jl")
include("models/timelessfixedtimemodel.jl")
include("models/homogeneousmodel.jl")
include("models/variabletimemodel.jl")
include("models/staticmodel.jl")

end


text:
module Regularismus

import Base: ==, hash, cat, int
export @wx_str, rx
export bots, tops, alt, star, plus
export botr, topr
export botm, topm, botd, topd
export mac

typealias WordType ASCIIString
const ABET = Set("abcd")
const NULL = 'e'

@assert typeof(NULL) == eltype(ABET)
@assert !(NULL in ABET)

include("Ex.jl")

include("mac.jl")

end


text:
function eos(model::SAFT, z, v, T)
    return N_A*k_B*sum(z) * T*(a_ideal(model,z,v,T)+a_res(model,z,v,T))
end

"""
The EoS is extensible to other types of equations.

julia > function Eos(model::SRK, p, T, z)
            return some_function(p, T, z)
        end

Then create some get_properties functions for other equations.
"""


text:
function schedule_day(employees, env)
    ok, message = validate_day(employees, env)
    if !ok
        Logging.info( message)
        return false, {}
    end

    ok, schedule = calculate_day(employees, env)
    if !ok
        return false, {}
    end

    # again a boolean "ok" followed by message
    ok, message = validate_day_schedule(employees, env, schedule)
    if !ok
        Logging.info( message)
        return false, {}
    end

    # boolean "ok" followed by data
    return true, schedule
end

function calculate_day(employees, env)
    total_time = length(env["coverage"])
    num_employees = length(employees)

    sum_coverage = sum(env["coverage"])

    if GUROBI
        m = Model(
            solver=GurobiSolver(
                LogToConsole=0,
                TimeLimit=CALCULATION_TIMEOUT,
                OutputFlag=0,
            )
        )
    else
        m = Model(
            solver=CbcSolver(
                threads=Base.CPU_CORES,
            )
        )
    end

    ###########################################################################
    # SHIFTS
    ###########################################################################

    # Build variable
    shift_start_max = Uint8[]
    length_max = Uint8[]
    length_min = Uint8[]
    worker_keys = Any[]
    for e in keys(employees) # e is a string
        push!(shift_start_max, total_time + 1 - employees[e]["min"])
        push!(length_max, employees[e]["max"])
        push!(length_min, employees[e]["min"])
        push!(worker_keys, e)
    end

    # Shift start time
    @defVar(m, 1 <= shift_start[1:num_employees] <= total_time, Int)
    for e in 1:num_employees # e is an int
        setUpper(shift_start[e], shift_start_max[e])
    end
    # Shift length - max & min
    @defVar(m, shift_length[e=1:num_employees], Int)
    for e in 1:num_employees # e is an int
        setLower(shift_length[e], length_min[e])
        setUpper(shift_length[e], length_max[e])

        # The "delivery" constraint:
        # prevent unassigned shifts from starting after a certain time
        if "no_shifts_after" in keys(env)
            @addConstraint(m,
                shift_start[e] <= env["no_shifts_after"],
            )
        end
    end

    ###########################################################################
    # BINARY CONVERSION
    ###########################################################################
    #=
    Set up binary variable for every employee for every week to determine
    whether somebody is working. This makes it easy to sum over the week and
    see how many people are working and whether s comply with availability.
    =#

    # build some helper variables to make my life more sane
    bin_total = total_time * num_employees
    i_to_t = Int[]
    i_to_e = Uint8[] # basically what is the equivalent s index
    i_employee = String[]
    #(could have used a modulus of cycle?)

    # we are stacking cycles together, offset by total_time
    # so all t=1 first, then all t=2, etc
    for t in 1:total_time
        e_index = 0
        for e in keys(employees)
            e_index+=1
            push!(i_to_t, t)
            push!(i_to_e, e_index)
            push!(i_employee, e)
        end
    end

    # This whole section is about converting a shift start time and length to
    # a binary array of whether it has started at time t and whether it is
    # active at time t

    # define variables
    @defVar(m, started[1:bin_total], Bin)
    @defVar(m, active[1:bin_total], Bin)

    # Need to split this into its constituent parts
    # (allows us to determine whether it is equal to zero, meaning
    # that the shift started)
    @defVar(m, 0 <= start_min_now_pos[1:bin_total] <= total_time, Int)
    @defVar(m, 0 <= start_min_now_neg[1:bin_total] <= total_time, Int)
    # Is it positive?
    @defVar(m, start_min_now_helper[1:bin_total], Bin)

    # Constraints
    for i in 1:bin_total
        # instead of SOS constraint on the start_min variables:
        # Use a big M approach
        M = total_time + 1
        @addConstraints m begin
            start_min_now_pos[i] <= M * start_min_now_helper[i]
            start_min_now_neg[i] <= M * (1 - start_min_now_helper[i])
            start_min_now_pos[i] - start_min_now_neg[i] == shift_start[i_to_e[i]] - i_to_t[i]
            1 <= start_min_now_pos[i] + start_min_now_neg[i] + started[i]
        end

        # Create one contiguous chunk of ones
        if i_to_t[i] == 1
            @addConstraint(m,
                started[i] == active[i],
            )
        else # not first time
            @addConstraint(m,
                # basically "t-1" for same shift
                started[i] >= active[i] - active[i-num_employees],
            )
        end
    end

    for e in 1:num_employees
        @addConstraints m begin
            # each shift can only have one start time
            1 == sum{started[i], i=e:num_employees:bin_total}
            # Shift binary array must sum to length variable
            shift_length[e] == sum{active[i], i=e:num_employees:bin_total}
        end
    end

    ###########################################################################
    # AVAILABILITY CONSTRAINTS
    ###########################################################################

    for i in 1:bin_total
        # Employees dont' work when they are unavailable
        if employees[i_employee[i]]["availability"][i_to_t[i]] == 0
            @addConstraint(m, active[i] == 0)
        end
    end

    ###########################################################################
    # ENVIRONMENT CONSTRAINTS
    ###########################################################################

    # Make sure that the organization has the proper number of people working
    for t in 1:total_time
        # No coverage == nobody working
        if env["coverage"][t] > 0
            @addConstraint(m,
                env["coverage"][t] <= sum{active[i], i=1:bin_total; i_to_t[i] == t}
            )
        else
            # nobody needed, so nobody should be working
            @addConstraint(m,
                0 == sum{active[i], i=1:bin_total; i_to_t[i] == t}
            )
        end
    end

    ###########################################################################
    # OPTIMIZATION
    ###########################################################################

    # Set the objective
    @setObjective(m, Min, sum{shift_length[e], e=1:num_employees})

    # optimize
    status = solve(m)
    sumHours = 0
    if status == :Optimal
        hours = getObjectiveValue(m)
        Logging.debug("Day scheduler: scheduled $hours / min $sum_coverage")
        e_index = 0
        schedule = Dict()
        for e in keys(employees)
            e_index += 1
            # NOTE - we round here. If we're debugging - don't round. Might get weird floats.
            start = convert(Int, round(getValue(shift_start[e_index])))
            length = convert(Int, round(getValue(shift_length[e_index])))
            sumHours += length
            schedule[e] = {
                "start" => start,
                "length" => length,
            }
        end
        # TODO - build return function
    else
        Logging.debug( "No optimal solution found")
        gc()
        return false, {}
    end
    # ok
    gc()
    return true, schedule
end

function validate_day(employees, env)
    total_time = length(env["coverage"])
    min_hours = 0
    max_hours = 0
    coverage_total = sum(env["coverage"])

    for e in keys(employees)
        if length(employees[e]["availability"]) != total_time
            return false, string("Incorrect availability length for ", e)
        end

        longest = longest_consecutive(employees[e]["availability"])

        if longest < employees[e]["min"]
            return false, string("employee ", e, " not available long enough")
        end

        if employees[e]["min"] > total_time
            return false, string("employee ", e, " has too long a min length")
        end

        min_hours += employees[e]["min"]
        max_hours += employees[e]["max"]
    end

    Logging.debug("Coverage $coverage_total / Min hours: $min_hours / Max hours: $max_hours")
    if coverage_total > max_hours
        return false, "coverage exceeds max employee hours"
    end

    # check that enough employees are available to work each shift
    for t in 1:total_time
        availability_sum = 0
        for e in keys(employees)
            availability_sum += employees[e]["availability"][t]
        end

        if availability_sum < env["coverage"][t]
            return false, string("not enough employees available at time ", t)
        end
    end

    return true, "42"
end

function validate_day_schedule(employees, env, schedule)
    total_time = length(env["coverage"])

    # helper variable for later
    coverage_check = zeros(Int, total_time)

    for e in keys(employees)
        start = schedule[e]["start"]
        length = schedule[e]["length"]

        # shift start time >= 1
        if start < 1
            return false, string("employee ", e, " start time out of bounds")
        end

        if length <= 0
            return false, string("employee ", e, " shift length invalid")
        end

        # end  time in bounds
        if (start + length - 1) > total_time # fucking non-zero-indices
            return false, string("employee ", e, " length out of bounds")
        end

        for t in start:(start + length - 1)
            # make sure employee is available then
            if employees[e]["availability"][t] == 0
                return false, string("employee ", e, " scheduled when unavailable")
            end
            # also add to coverage tab for checking later
            coverage_check[t] += 1
        end
    end

    for t in 1:total_time
        # make sure we have enough peeps working
        if coverage_check[t] < env["coverage"][t]
            return false, string("not enough people working at time ", t)
        end
    end
    return true, ""
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
type AdamSolver <: Solver
  timestep::Int64
  alpha::Float64
  beta1::Float64
  beta2::Float64
  epsilon::Float64

  mts::Array{Any}
  vts::Array{Any}

  AdamSolver() = new(1, 0.001, 0.9, 0.999, 1e-8, [], [])
end

type AdamSolverParams <: SolverParams
    clipval::Float64
end

function step!(solver::AdamSolver, model::Model, sparams::AdamSolverParams)

  # perform parameter update

  # All of the matrices used by the model
  modelMatices = model.matrices

  # init
  if length(solver.mts) == 0
    for m in modelMatices
      push!(solver.mts, zeros(m.n, m.d))
    end
  end

  if length(solver.vts) == 0
    for m in modelMatices
      push!(solver.vts, zeros(m.n, m.d))
    end
  end
  alpha = solver.alpha;
  beta1 = solver.beta1;
  beta2 = solver.beta2;
  at = alpha*(sqrt(1-(beta2^solver.timestep)))/(1 - beta1^solver.timestep);

  for k = 1:length(modelMatices)
    @inbounds A = modelMatices[k] # mat ref
    @inbounds mt = solver.mts[k]
    @inbounds vt = solver.vts[k]
    for i = 1:A.n
      for j = 1:A.d

        # get the gradients
        @inbounds gt = A.dw[i,j]

        # gradient clip
        if gt > sparams.clipval || gt < -sparams.clipval
          gt = sparams.clipval*sign(gt);
        end

        # update biased 1st moment estimate
        @inbounds mt[i,j] = beta1*mt[i,j] + (1-beta1)*gt;
        @inbounds vt[i,j] = beta2*vt[i,j] + (1-beta2)*(gt.*gt);

        # update
        @inbounds A.w[i,j] += -at * mt[i,j]/(sqrt(vt[i,j]) + solver.epsilon);
        @inbounds A.dw[i,j] = 0. # reset gradients for next iteration

      end
    end
  end
  solver.timestep = solver.timestep + 1;
  return
end

max_forks_repo_name: haroldsoh/RecurrentNN


text:
using MPI
using JuliaChem.JERI

#== initialize JuliaChem runtime ==#
function initialize()
  if(!MPI.Initialized())
    MPI.Init()
    JERI.initialize()
  else
    println("JuliaChem has already been initialized!")
  end

  #==set up scratch directory==#

end
export initialize

#== finalize JuliaChem runtime ==#
function finalize()

  #== finalize MPI ==#
  if(!MPI.Finalized())
    MPI.Finalize()
    JERI.finalize()
  else
    println("JuliaChem has already been finalized!")
  end

  #==clean scratch directory==#
end
export finalize


text:
module MarkableFloats

export MarkableFloat, MarkFloat64, Marked, Unmarked, mark, unmark,
    square, cube

import Base: AbstractFloat, Float64, float

import Base: @pure, promote_type, promote_rule, convert,
    mark, unmark, ismarked,
    typemax, typemin, realmax, realmin,
    string, show,
    (<=), (<), (==), (!=), (>=), (>), isless, isequal,
    (~), (&), (|), (⊻),
    findall

import Base.Math: zero, one, iszero, isone, isinteger,
    isodd, iseven, sign, signbit, abs, copysign, flipsign,
    (+), (-), (*), (/), (\), (^),
    round, trunc, ceil, floor,
    div, fld, cld, mod, rem,
    sqrt, cbrt,
    log, log2, log10, log1p, exp, expm1,
    sin, cos, tan, sec, csc, cot, sinpi, cospi,
    asin, acos, atan, asec, acsc, acot, atan2,
    sinh, cosh, tanh, sech, csch, coth,
    asinh, acosh, atanh, asech, acsch, acoth

include("type.jl")
include("math.jl")

end # MarkableFloats


text:
import Test
import PredictMD

table1 = [20 5; 10 15]
Test.@test(
    isapprox(PredictMD.cohen_kappa(table1), 0.4; atol = 0.00000000001)
    )

table2 = [45 15; 25 15]
Test.@test(
    isapprox(PredictMD.cohen_kappa(table2), 0.1304; atol = 0.0001)
    )

table3 = [25 35; 5 35]
Test.@test(
    isapprox(PredictMD.cohen_kappa(table3), 0.2593; atol = 0.0001)
    )

table4 = [9 3 1; 4 8 2 ; 2 1 6]
Test.@test(
    isapprox(PredictMD.cohen_kappa(table4), 0.45; atol = 0.001)
    )


text:
export lbgk
using .LB: Lattice
using .Units: dimensionless, LBUnits

# Sets up calculation with standard lbgk
# All units should be in SI:
#  - time -> s
#  - length -> m
#  - pressure -> Pa
#  - viscosity -> Pa * s
function lbgk{N, I <: Integer}(
        lattice::Lattice, dimensions::NTuple{N, I},
        δt::Number,                       # s
        δx::Number;                       # m
        viscosity::Number=1e-3,           # Pa*s
        μ₀::Vector=[0., 0],               # m/s,
        p₀::Number=80.0 * 133.3223684211, # mmHG * x = Pa
        Δp::Number=0,                     # Pa
        kwargs...)

    const T = typeof(lattice).parameters[1]
    const ρ₀=1e3                    # kg*m^-3
    const δm = ρ₀ * δx^3
    const units = LBUnits{T}(δt, δx, δm)

    const μ = dimensionless(units, :velocity, μ₀)
    const p = dimensionless(units, :pressure, p₀)
    const c_s² = 1./3.
    const c_s = sqrt(c_s²)
    const τ = 1. + dimensionless(units, :viscosity, viscosity) / c_s²

    density(p) = 1. + dimensionless(units, :pressure, p) / c_s²

    # Set up main item
    result = SandBox(lattice, dimensions; ρ₀=density(p₀ - 0.5Δp), μ₀=μ, kwargs...)
    # Set up kernel and streamers with some default values
    streamers = Dict{Playground.Feature, LB.Streaming}(
        Playground.INLET  => NashZeroOrderPressure{T}(T[0, 1], density(p₀)),
        Playground.OUTLET => NashZeroOrderPressure{T}(T[0, 1], convert(T, density(p₀ - Δp))),
        Playground.FLUID  => FluidStreaming(),
        Playground.SOLID  => HalfWayBounceBack()
    )
    result.kernels = Dict(
        Playground.FLUID =>
            FluidKernel(SingleRelaxationTime{T}(1.0/τ), streamers))

    # Set up playground as box with inflow on left, outflow on right
    if length(dimensions) == 2
        Playground.initialize(result.playground) do i, j
            if j == 1 || j == dimensions[2]
                return Playground.SOLID
            elseif i == 1
                return Playground.INLET
            elseif i == dimensions[1]
                return Playground.OUTLET
            else
                return Playground.FLUID
            end
        end
    else
        Playground.initialize(result.playground) do i, j, k
            if j == 1 || j == dimensions[2] || k == 1 || k == dimensions[3]
                return Playground.SOLID
            elseif i == 1
                return Playground.INLET
            elseif i == dimensions[1]
                return Playground.OUTLET
            else
                return Playground.FLUID
            end
        end
    end
    result
end

function lbgk{N, I <: Integer}(lattice::Symbol, dimensions::NTuple{N, I},
                               δt::Number, δx::Number; kwargs...)
    lbgk(getfield(LB, lattice)::LB.Lattice, dimensions, δt, δx; kwargs...)
end


text:
# Exercise "Histogram of the English alphabets"
# Julia code generate the histogram
using Plots,DelimitedFiles
f = readdlm("ch3_data_english.txt")
bar(f/100,xticks = (1:26,'a':'z'),yticks = 0.0:0.02:0.1,label = false)

# Exercise "Histogram of the throwing a die"
# Julia code generate the histogram
using Plots,Random
q = rand(1:6,100)
histogram(q,bins = 6,normalize = true,label = "N = 100")

# Exercise "Histogram of an exponential random variable"
# Julia code used to generate the plots
using Plots,Random,Distributions
λ = 1
X = rand(Exponential(1/λ),1_000)
histogram(X,bins = 200,label = "K = 200")

# Exercise "Cross validation loss"
using Plots,Random,Distributions,StatsBase
λ = 1
n = 1_000
X = rand(Exponential(1/λ),n)
m = 6:200
J = zeros(195)
for i = 1:195
    hist = fit(Histogram,X,range(minimum(X),maximum(X),length = m[i]+1))
    h = n/m[i]
    J[i] = 2/((n-1)*h)-((n+1)/((n-1)*h))*sum( (hist.weights/n).^2 )
end

plot(m,J,label = false)

# Exercise "Expectation E[X] where X is uniform from [0,1]"
using Random,Statistics
X = rand(10_000)
mean(X)

# Exercise "Mean from PMF"
# Julia code to compute the expectation
p = [0.25,0.5,0.25]
x = [0,1,2]
EX = sum(p .* x)

# Exercise "Mean of geometric random variable"
k = 1:100
p = 0.5.^k
EX = sum(p .* k)

# Exercise "Bernoulli Random Variables"
using Plots,Distributions
p = 0.5
n = 1
X = rand(Binomial(n,p),1_000)
histogram(X,bins = 2,label = false)

# Julia code to generate Erdos Renyi Graph
using Plots,Random,LinearAlgebra,GraphRecipes
A = rand(40,40) .< 0.1
A = triu(A,1)
A = A + A'
graphplot(A)

# Exercise "Binomial Random Variables" 2
using Plots,Distributions
p = 0.5
n = 10
X = rand(Binomial(n,p),5_000)
histogram(X,label = false)

# Exercise "Binomial CDF"
using Plots,Distributions
x = 0:10
p = 0.5
n = 10
F = cdf.(Binomial(n,p),x)
plot(x,F,linetype = :steppost,label = false,markershape = :circle)

# Exercise "Poisson CDF"
## Code for the Poisson PMF
using Plots,Distributions
λ_set = [1,4,10]
p = zeros(21,3)
k = 0:1:20
for i = 1:3
    p[:,i] = pdf.(Poisson(λ_set[i]),k)
end
plot(k,p[:,1],line = :stem,markershape = :circle,label = "λ = 1")
plot!(k,p[:,2],line = :stem,markershape = :circle,label = "λ = 4")
plot!(k,p[:,3],line = :stem,markershape = :circle,label = "λ = 10")

## Code for Poisson CDF
using Plots,Distributions
λ_set = [1,4,10]
p = zeros(21,3)
k = 0:1:20
for i = 1:3
    p[:,i] = cdf.(Poisson(λ_set[i]),k)
end
plot(k,p[:,1],line = :steppost,markershape = :circle,label = "λ = 1")
plot!(k,p[:,2],line = :steppost,markershape = :circle,label = "λ = 4")
plot!(k,p[:,3],line = :steppost,markershape = :circle,label = "λ = 10",legend = :bottomright)

# Exercise "Poisson-Binomial approximation"
using Plots,Distributions
n = 5_000
p = 0.01
λ = n*p
x = 0:120
y = pdf.(Binomial(n,p),x)
z = pdf.(Poisson(λ),x)
plot(x,y,line = :stem,markershape = :circle,label = "Binomial,n = 5000,p = 0.01")
plot!(x,z,label = "Poisson,λ = 50")

# Photon Shot Noise
using Distributions,Images,ImageView,Random
x = Float64.(load("cameraman.tif"))

α = 10
λ = α*x
X1 = rand.(Poisson.(λ))
imshow(X1)

α = 100
λ = α*x
X1 = rand.(Poisson.(λ))
imshow(X1)

α = 1_000
λ = α*x
X1 = rand.(Poisson.(λ))
imshow(X1)


text:
"""
Body of a request for creating or updating a data element.
"""
mutable struct BigDataElementRequest
    id::String
    sourceName::String
    description::String
    data::String
    mimeType::String
    BigDataElementRequest(id::String, sourceName::String, description::String, data::String, mimeType::String="application/octet-stream") = new(id, sourceName, description, data, mimeType)
end

function show(io::IO, c::BigDataElementRequest)
    println(io, "Graff Data Element Request: ID = $(c.id), MIME type = $(c.mimeType), length = $(length(c.data)), description = $(c.description)")
end

"""
Summary of data entry returned from request.
"""
mutable struct BigDataEntryResponse
    id::String
    sourceId::String
    neoNodeId::Int
    sourceName::String
    sourceParams::Dict{String, Any}
    description::String
    mimeType::String
    lastSavedTimestamp::String
    links::Union{Nothing, Dict{String, Any}}
end

function show(io::IO, c::BigDataEntryResponse)
    println(io, "Graff Data Entry Response: ID = $(c.id), MIME type = $(c.mimeType), description = $(c.description)")
end

"""
Complete data element response.
"""
mutable struct BigDataElementResponse
    id::String
    mimeType::String
    data::String
end

function show(io::IO, c::BigDataElementResponse)
    println(io, "Graff Data Element Response: ID = $(c.id), MIME type = $(c.mimeType), length = $(length(c.data))")
end


text:
include("core.jl")
include("types.jl")
#include("dev.jl")


text:
using Abaco
using Test

Abaco.loginit()

# My time zone: Saturday, October 30, 2021 3:43:35 PM GMT+02:00 DST
metric_ts = 1635601415
metric_sn = "Marmolada"

width = -1
ages = 4

expected_triggers = 4
actual_triggers = 0

values = Dict(
    "x" => 2.0,
    "y" => 4.0,
    "z" => 8.0
    )

function onresult(ts, ne, name, value, inputs)
    global actual_triggers
    @info "age [$ts]: scope: [$ne] $name = $value"
    @test ts in [metric_ts, metric_ts+5]
    @test ne == metric_sn
    @test value == 14.0
    actual_triggers += 1
end

abaco = abaco_init(onresult, interval=width, emitone=false)

# add a formula that is not evaluated because variable t is not feeded
formula(abaco, "r = (x + y) / exp(t)")

# this gets evaluated (2 + 4 + 8 = 14)
formula(abaco, "w = x + y + z")

ingest(abaco, metric_ts, metric_sn, values)

ingest(abaco, metric_ts+5, metric_sn, values)

@test actual_triggers === expected_triggers


text:
# ==========================================================================

immutable SVC{TModel <: SVM, TEnc <: BinaryClassEncoding}
  model::TModel
  encoding::TEnc
end

function SVC(model::SVM, encoding::BinaryClassEncoding)
  SVC{typeof(model),typeof(encoding)}(model, encoding)
end

# ==========================================================================

immutable SVR{TModel <: SVM}
  model::TModel
end

SVR(model::SVM) = SVR{typeof(model)}(model)

# ==========================================================================

for kind = (:SVC, :SVR)
  for op = (:details, :isconverged, :iterations, :features, :params, :intercept,
            :nsv, :svindex, :xsv, :predmodel, :bias, :minimum, :minimizer,
            :decision_function, :nobs, :model_response, :coef, :predict, :accuracy)
    @eval begin
      ($op)(svm::$kind) = ($op)(svm.model)
    end
  end
  @eval (predict)(svm::$kind, X) = (predict)(svm.model, X)
end

targets(svc::SVC) = labeldecode(svc.encoding, targets(svc.model))
ysv(svc::SVC) = labeldecode(svc.encoding, ysv(svc.model))

targets(svc::SVR) = targets(svc.model)
ysv(svc::SVR) = ysv(svc.model)

# ==========================================================================

labels(svc::SVC) = labels(svc.encoding)

classify(svc::SVC) = classify(svc, features(svc.model))

function classify(svc::SVC, X)
  t = classify(svc.model, X)
  labeldecode(svc.encoding, t)
end

function accuracy(svc::SVC, X, y)
  accuracy(svc.model, X, labelencode(svc.encoding, y))
end

# ==========================================================================

function scatterplot(scv::SVC, args...; nargs...)
  scatterplot(scv.model, args...; lbl = map(string, labels(scv)), nargs...)
end

# ==========================================================================
# Base.show

function show(io::IO, fit::SVC)
  println(io, typeof(fit), "\n")
  if typeof(fit.model) <: PrimalSVM
    _showprimal(io, fit)
  else
    _showdual(io, fit)
  end
end


text:
show_location(X, Y, Z) = string("(", string(typeof(X())), ", ",
                                     string(typeof(Y())), ", ",
                                     string(typeof(Z())), ")")

show_location(field::AbstractLocatedField{X, Y, Z}) where {X, Y, Z} = show_location(X, Y, Z)

short_show(a) = string(typeof(a))
shortname(a::Array) = string(typeof(a).name.wrapper)

show(io::IO, field::Field) =
    print(io,
          short_show(field), '\n',
          "├── data: ", typeof(field.data), '\n',
          "└── grid: ", typeof(field.grid), '\n',
          "    ├── size: ", size(field.grid), '\n',
          "    └── domain: ", show_domain(field.grid), '\n')

short_show(field::AbstractLocatedField) = string("Field at ", show_location(field))


text:
#=
ambiguity.jl
Avoiding method ambiguities
This may be a bit of a whack-a-mole™ exercise...
=#

import LinearAlgebra

const Trans = LinearAlgebra.Transpose{<: Any, <: LinearAlgebra.RealHermSymComplexSym}
const Adjoi = LinearAlgebra.Adjoint{<: Any, <: LinearAlgebra.RealHermSymComplexHerm}

Base.:(*)(A::LinearMapAM, D::LinearAlgebra.Diagonal) = AM_M(A, D)
Base.:(*)(D::LinearAlgebra.Diagonal, B::LinearMapAM,) = M_AM(D, B)

Base.:(*)(A::LinearMapAM, B::LinearAlgebra.AbstractTriangular) = AM_M(A, B)
Base.:(*)(A::LinearAlgebra.AbstractTriangular, B::LinearMapAM) = M_AM(A, B)

Base.:(*)(A::LinearMapAM, B::Trans) = AM_M(A, B)
Base.:(*)(A::Trans, B::LinearMapAM) = M_AM(A, B)

Base.:(*)(A::LinearMapAM, B::Adjoi) = AM_M(A, B)
Base.:(*)(A::Adjoi, B::LinearMapAM) = M_AM(A, B)

# see https://github.com/Jutho/LinearMaps.jl/issues/118
Base.:(*)(A::LinearMapAM,
   B::LinearAlgebra.Adjoint{<: Any, <: LinearAlgebra.AbstractRotation}) =
   throw("AbstractRotation lacks size so * is unsupported")
#Base.:(*)(A::Given, B::LinearMapAM) = M_AM(A, B)

Base.:(*)(x::LinearAlgebra.AdjointAbsVec, A::LinearMapAM) = (A' * x')'
Base.:(*)(x::LinearAlgebra.TransposeAbsVec, A::LinearMapAM) =
    transpose(transpose(A) * transpose(x))


text:
include("LDA_VWN.jl")

function test_main()
    rhoe = 1.5*ones(5)
    epsxc = excVWN( rhoe )
    depsxc = excpVWN( rhoe )
    println(epsxc)
    println(depsxc)
end

test_main()


text:
module HiQGA
module transD_GP
    include("GP.jl")
    include("TransD_GP_MCMC.jl")
    include("AbstractOperator.jl")
    include("MCMC_Driver.jl")
    include("GradientInversion.jl")
    include("GeophysOperator.jl")
end
export transD_GP
end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
function F₀(x::AbstractArray{T}, α::T, λ::T) where T<:Real
    out = zero(T)
    @inbounds for i in 1:length(x)
        δᵢ = max(abs(x[i])-λ,0)
        out += δᵢ
    end
    out -= α+λ
    out
end

function F₁(x::AbstractArray{T}, α::T, λ::T) where T<:Real
    if λ==0
        return sum(abs,diff(x))
    else
        x_prox = fit(FusedLasso,x,λ).β
        return sum(abs,diff(x_prox))-α-λ
    end
end

function absepi_projection!(θ_proj::AbstractArray{T}, skip::Int, θ::AbstractArray{T}, α::T) where T<:Real
    copyto!(θ_proj,θ)
    n = length(θ)
    @views if sum(abs, θ[(skip+1):n]) <= α
        return α
    else
        λstar = find_zero(λ -> F₀(θ[(skip+1):n],α,λ),(0,maximum(abs,θ[(skip+1):n])+max(0,-α)),Bisection())
        @inbounds for i in 1:(n-skip)
            θ_proj[skip+i] = sign(θ[skip+i])*max(abs(θ[skip+i])-λstar,0)
        end
        return α + λstar
    end
end

function diffepi_projection!(θ_proj::AbstractArray{T}, skip::Int, θ::AbstractArray{T}, α::T, inv_DDᵀ_D::AbstractArray{T}) where T<:Real
    copyto!(θ_proj,θ)
    n = length(θ)
    if sum(abs, diff(θ[(skip+1):n])) <= α
        return α
    else
        λₘ = maximum(abs, inv_DDᵀ_D*θ[(skip+1):n])
        λstar = find_zero(λ -> F₁(θ[(skip+1):n],α,λ),(0,λₘ+max(0,-α)),Bisection())
        θ_proj[(skip+1):n] = fit(FusedLasso,θ[(skip+1):n],λstar).β
        return α + λstar
    end
end


text:
root_path = raw"D:\AppliRTE\repo\scopf-quanti"
push!(LOAD_PATH, root_path);
cd(root_path);

include(joinpath(root_path, "AmplTxt.jl"));
include(joinpath(root_path, "SCOPF.jl"));

# Compute ptdf and export them as node branch value file

# test_name = "RC_PPE2035_S0_deboucle"
test_name = "5buses_wind"
test_name = "2buses"

amplTxt = AmplTxt.read(test_name);
test_path = joinpath(root_path, test_name)
network = SCOPF.Network(amplTxt);

ref_bus = 1 ;
B = SCOPF.get_B(network, 1e-6);
Binv = SCOPF.get_B_inv(B, ref_bus);
PTDF = SCOPF.get_PTDF(network, Binv, ref_bus);

PTDF_distributed = SCOPF.distribute_slack(PTDF, 1);

PTDF_TRIMMER = 1e-6;

# println(PTDF);
file_path =  joinpath(test_path, "ptdf.txt")
SCOPF.write_PTDF(file_path, network, ref_bus, PTDF, PTDF_TRIMMER)
file_path =  joinpath(test_path, "pscopf_ptdf.txt")
SCOPF.write_PTDF(file_path, network, ref_bus, PTDF_distributed, PTDF_TRIMMER)

max_stars_repo_name: azized/pscopf

max_issues_repo_name: azized/pscopf

max_forks_repo_name: azized/pscopf


text:
struct Ctx1 end
struct Ctx2 end
@context Ctx1
@context Ctx2

@testset "inline contextual" begin

    @testset "@contextualized" begin
        function get_contextualized()::Int64
            if @contextualized(Ctx1)
                return 1
            elseif @contextualized(Ctx2)
                return 2
            else
                return 0
            end
        end

        @test get_contextualized() == 0
        @test TinyCassette.execute(nothing, get_contextualized) == 0
        @test TinyCassette.execute(Ctx1(), get_contextualized) == 1
        @test TinyCassette.execute(Ctx2(), get_contextualized) == 2
    end

    @testset "@contextual" begin
        function get_ctx1()::Ctx1
            return @contextual(Ctx1)
        end

        function get_ctx2()::Ctx2
            return @contextual(Ctx2)
        end

        @test_throws MethodError get_ctx1()
        @test_throws MethodError get_ctx2()
        @test isa(TinyCassette.execute(Ctx1(), get_ctx1), Ctx1)
        @test isa(TinyCassette.execute(Ctx2(), get_ctx2), Ctx2)
    end
end


text:
"""
`Resize <: ImageTransformation`

Description
============

Transforms the image intoto a fixed specified pixel size. This
operation does not take any measures to preserve aspect ratio
of the source image.  Instead, the original image will simply be
resized to the given dimensions. This is useful when one needs a
set of images to all be of the exact same size.

Usage
======

    Resize()

    Resize(width = 64, height = 64)

Arguments
==========

- **`width`** : The width in pixel of the output image

- **`height`** : The height in pixel of the output image

Methods
========

- **`transform`** : Applies the transformation to the given Image
or set of images

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
========

    using Augmentor
    using TestImages

    # load an example image
    img = testimage("lena")

    # Resize the image to exactly 32x32 pixel
    img_new = transform(Resize(32, 32), img)

see also
=========

`ImageTransformation`, `ProbableOperation`, `transform`
"""
immutable Resize <: ImageTransformation
    width::Int
    height::Int

    function Resize(width::Int, height::Int)
        width >= 1 || throw(ArgumentError("width has to be greater than 0"))
        height >= 1 || throw(ArgumentError("height has to be greater than 0"))
        new(width, height)
    end
end

Resize(;width::Int = 64, height::Int = 64) = Resize(width, height)

Base.showcompact(io::IO, op::Resize) = print(io, "Resize to $(op.width)x$(op.height).")
multiplier(::Resize) = 1

function transform{T<:AbstractImage}(op::Resize, img::T)
    result = if isxfirst(img)
        copyproperties(img, Images.imresize(img, (op.width, op.height)))
    else
        copyproperties(img, Images.imresize(img, (op.height, op.width)))
    end::T
    _log_operation!(result, op)
end


text:
partial = false
include("lorenz.jl")
include("../extra/makie.jl")
using Makie, GLVisualize
using GeometryTypes, Colors
using Bridge: _viridis
skippoints = 1;
extra = 80
scene = Scene(resolution = (800, 800))
eyeposition = Float32[-37.2262, 65.8173, 37.4411]
sphere = Sphere(Point3f0(0,0,0), 1.0f0)
circle = Sphere(Point2f0(0,0), 1.0f0)
#RGB(Bridge._viridis[20]...)
maxviri = 200

en = 1000

#for i in indices(Paths)[1]
#    viri = _viridis[round.(Int,linspace(1,maxviri, length(Paths[i])))]
#    lines(Paths[i], linewidth = 0.1, color = map(x->RGBA(Float32.(x)..., 0.5f0), viri) )
#end
viri = _viridis[round.(Int,range(1, stop=maxviri, length=length(Xtrue)))]
lines(( Xtrue.yy[1:skippoints:en]), linewidth = 0.4, color = map(x->RGBA(Float32.(x)..., 1.0f0), viri[1:skippoints:en]) )
#lines(Xtrue.yy, linewidth = 0.2, color = :black )
scatter(Xtrue.yy[1:1], marker=circle,  markersize = 0.4, color= map(x->RGBA(Float32.(x)..., 1.0f0),viri[1:1]))

viri = _viridis[round.(Int,range(1, stop=maxviri, length=length(Xmean)))]
#lines((@view Xmean[1:extra*skippoints:end]), linewidth = 0.2, color = map(x->RGBA(Float32.(x)..., 0.5f0), viri[1:extra*skippoints:end]) )
scatter( Xmean[1:extra*skippoints:en], marker = circle, markersize=0.2, color = map(x->RGBA(Float32.(x)..., 0.5f0), viri[1:extra*skippoints:en]) )

viri = _viridis[round.(Int,range(1, stop=maxviri, length=length(Xscal)))]
visualize_uncertainty(scene, (Xmean[1:en], Xrot[1:en], Xscal[1:en]), extra*skippoints; color = map(x->RGBA(Float32.(x)..., 0.05f0), viri[1:extra*skippoints:en]) )

ma = Float32[10,10,10]
axis(map(x -> UnitRange(0, round.(Int,x)), ma)..., showaxis = (true,true,true), showgrid = (false,false,false))
ma2 = ma .+ 2
for (x, c) in zip(ℝ{3}[(ma2[1], 0, 0), (0, ma2[2], 0),(0, 0, ma2[3]), (0,0,0)], ['x', 'y', 'z', '•'])
    scatter([x], marker = c, markersize = 1., color = :black)
end

center!(scene)
println(Makie.getscreen(scene).cameras[:perspective].view.value)

push!(Makie.getscreen(scene).cameras[:perspective].view,
    @SArray Float32[0.974148 0.225913 2.23517f-8 -5.723; -0.181213 0.781398 -0.597143 5.67828; -0.134902 0.581705 0.802135 -59.5837; 0.0 0.0 0.0 1.0]
)

Makie.save(joinpath("output", simname, "lorenz6n$simid.png"), scene);


text:
module DataStructures

    import Base: length, isempty, start, next, done,
                 show, dump, empty!, getindex, setindex!, get, get!,
                 in, haskey, keys, merge, copy, cat,
                 push!, pop!, shift!, unshift!,
                 union!, delete!, similar, sizehint,
                 isequal, hash,
                 map, reverse

    export Deque, Stack, Queue
    export deque, enqueue!, dequeue!, update!
    export capacity, num_blocks, front, back, top, sizehint

    export Accumulator, counter
    export ClassifiedCollections
    export classified_lists, classified_sets, classified_counters

    export IntDisjointSets, DisjointSets, num_groups, find_root, in_same_set
    export push!

    export AbstractHeap, compare, extract_all!
    export BinaryHeap, binary_minheap, binary_maxheap
    export MutableBinaryHeap, mutable_binary_minheap, mutable_binary_maxheap

    export OrderedDict, OrderedSet
    export DefaultDict, DefaultOrderedDict
    export Trie, subtrie, keys_with_prefix, path

    export LinkedList, Nil, Cons, nil, cons, head, tail, list, filter, cat,
           reverse

    include("delegate.jl")

    include("deque.jl")
    include("stack.jl")
    include("queue.jl")
    include("accumulator.jl")
    include("classifiedcollections.jl")
    include("disjoint_set.jl")
    include("heaps.jl")

    include("hashdict.jl")
    include("ordereddict.jl")
    include("orderedset.jl")
    include("defaultdict.jl")
    include("trie.jl")

    include("list.jl")

    @deprecate stack Stack
    @deprecate queue Queue
    @deprecate add! push!
end


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    C() {
    }
    int xTESTOUTPUT = 7;
    D d = new D();
    void m() {
	int y = ({L /this.d.c H}) 42;

      d.foo();

      this.xTESTOUTPUT = y; // At this point, the level of y should be H, since this.d.c has been set.

    }

    public static void main(String[] args) {
	new C().m();
    }
}

class D {
    static final Condition{L} c = new Condition();
    void foo() {
	c.set();
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
using FluxExtensions
using Base.Test
using Flux
using ForwardDiff
using CoordinateTransformations

@testset "Flux extensions" begin
    srand(1)
    models = [
        Chain(
            AffineMap(randn(1, 1), randn(1)),
            Dense(1, 1, elu),
            Dense(1, 1, elu),
        ),
        Chain(
            Dense(1, 10, elu),
            Dense(10, 1, elu),
            AffineMap(randn(1, 1), randn(1)),
        ),
        Chain(
            Dense(1, 10, elu),
            AffineMap(randn(10, 10), randn(10)),
            Dense(10, 1, elu),
        ),
    ]
    for m in models
        mp = FluxExtensions.TangentPropagator(m)
        p = FluxExtensions.untrack(m)

        for i in 1:100
            x = randn(1)
            y = m(x)
            y2, J = mp(x)
            @test FluxExtensions.value(y) ≈ FluxExtensions.value(y2)
            @test p(x) ≈ FluxExtensions.value(y)
            @test FluxExtensions.value(J) ≈ ForwardDiff.jacobian(p, x)
        end

        lf = (x, y, J) -> begin
            ŷ, Ĵ = mp(x)
            Flux.mse(ŷ, y) + Flux.mse(Ĵ, J)
        end
        train_data = [
            ([1.0], [1.2], [1.5])
        ]
        opt = Flux.Optimise.Momentum(params(mp))
        for i in 1:1000
            Flux.train!(lf, train_data, opt)
        end
        ŷ, Ĵ = mp([1.0])
        @test FluxExtensions.value(ŷ) ≈ [1.2]
        @test FluxExtensions.value(Ĵ) ≈ [1.5]

    end
end


text:
using AntennaExpts, PIFA_GE

pifa_ge()


text:
"""
    package ModiaMedia

Media models for use with Modia and other Julia packages.
It is planned that this package contains media property models to be used in package
[Modia](https://github.com/ModiaSim/Modia.jl), but also in other Julia packages.
The initial goal is to achieve a similar functionality as
[Modelica.Media](https://doc.modelica.org/Modelica%203.2.3/Resources/helpDymola/Modelica_Media.html#Modelica.Media),
the standard media library for Modelica models, but with improvements based on Julia's features
such as multiple dispatch.

A medium model is basically a struct that has the following structure:

```
mutable struct MediumName <: AbstractMedium
    infos::ModiaMedium.FluidInfos
    fluidConstants::SVector{1,ModiaMedium.AbstractFluidConstants}
    fluidLimits::ModiaMedium.FluidLimits
    data::MediumSpecificData
end
```

and functions operating on instances of such a struct.

This package is currently under development.
"""
module ModiaMedia

const path    = dirname(dirname(@__FILE__))          # Absolute path of package directory
const Version = "0.1.0-dev"
const Date    = "2019-02-25"

println(" \nImporting ModiaMedia Version $Version ($Date)")

# Export Abstract types
export AbstractMedium
export PureSubstance
export MixtureMedium
export CondensingGases
export ThermodynamicState
export MixtureThermodynamicState
export AbstractFluidConstants

# Export types used in a Medium
export FluidInfos
export FluidLimits
export BasicFluidConstants
export IdealGasFluidConstants

# Export State structs
export SimpleMediumState
export SingleGasNasaState
export SimpleIdealGasMediumState
export MoistAirState

# Export general functions
export getMedium, listMedia, standardCharacteristics, standardPlot

# Export set state functions
export setState_pTX, setState_pTX!, setState_pT, setState_pT!
export setState_phX, setState_phX!, setState_ph, setState_ph!
export setState_psX, setState_psX!, setState_ps, setState_ps!
export setState_dTX, setState_dTX!, setState_dT, setState_dT!
export isenthalpicState, isenthalpicState!

# Export thermodynamic property functions for all media
export density, density_phX, density_pTX, density_der_1, density_pT, density_pT_der_1, density_pT_der_2, density_pT_der_3
export specificInternalEnergy_T, specificInternalEnergy_T_der_1, specificInternalEnergy_T_der_2
export temperature, temperature_phX, temperature_ph
export pressure, pressure_dT
export specificEnthalpy, specificEnthalpy_pTX, specificEnthalpy_dT, specificEnthalpy_T
export specificInternalEnergy, specificHeatCapacityCp
export dynamicViscosity
export gasConstant

# Export PureSubstance thermodynamic property  functions
export density_ph, temperature_ph, pressure_dT, specificEnthalpy_dT

# Export MoistAir specific functions
export saturationPressureOfLiquidWater
export sublimationPressureIce
export saturationPressureOfWater

# Enumeration values
export IndependentVariables, IndependentVariables_T, IndependentVariables_pT, IndependentVariables_ph
export IndependentVariables_phX, IndependentVariables_pT, IndependentVariables_dTX
export ReferenceEnthalpy, ReferenceEnthalpy_ZeroAt0K, ReferenceEnthalpy_ZeroAt25C, ReferenceEnthalpy_UserDefined
export ReferenceEntropy , ReferenceEntropy_ZeroAt0K , ReferenceEntropy_ZeroAt0C  , ReferenceEntropy_UserDefined

### Abstract types -------------------------------------------------------------------------------
# The following structures and names are identical to Modelica.Media.Interfaces
# with the only exception, that "Partial" is replaced by "Abstract"

"""
    abstract type AbstractMedium

Abstract type of all media.
"""
abstract type AbstractMedium end

"""
    abstract type PureSubstance <: AbstractMedium

Abstract type of all media consisting of a pure substance.
"""
abstract type PureSubstance <: AbstractMedium end

"""
    abstract type MixtureMedium <: AbstractMedium

Abstract type of all media consisting of a mixture of media.
"""
abstract type MixtureMedium <: AbstractMedium end

"""
    abstract type CondensingGases <: AbstractMedium

Abstract type of all media consisting of condensing media.
"""
abstract type CondensingGases <: MixtureMedium end

"""
   abstract type ThermodynamicState

Abstract type of all media states.
"""
abstract type ThermodynamicState end

"""
    abstract type MixtureThermodynamicState <: ThermodynamicState

Abstract type of the states of all media consisting of a mixture of media.
"""
abstract type MixtureThermodynamicState <: ThermodynamicState end

"""
    abstract type AbstractFluidConstants

Abstract type of all FluidConstants structures.
"""
abstract type AbstractFluidConstants end

### Importing packages -------------------------------------------------------------------------------
using  JSON
using  StaticArrays
using  Unitful
import DataFrames
import ModiaMath
import Serialization

### Including files for the ModiaMedia module --------------------------------------------------------
include("Interfaces/PartialMedium.jl")
include("Interfaces/PartialPureSubstance.jl")
include("Interfaces/PartialMixtureMedium.jl")
include("Interfaces/PartialCondensingGases.jl")

include("Media/SimpleMedium.jl")
include("Media/SimpleIdealGasMedium.jl")
include("Media/SingleGasNasa.jl")
include("Media/MoistAir.jl")

const mediumDictFile   = "$path/src/Media/media.julia_serializer"
const generateDictFile = "$path/dict/GenerateMediumDict.jl"
const mediumDict       = Array{Any}(nothing,1)

### Load medium dictionary from file
function loadMediumDictFile()
    global mediumDictFile
    println("... Read media dictionary from file:\n",
            "    ", mediumDictFile)
    f = open(mediumDictFile)
    dict = Serialization.deserialize(f)
    close(f)
    return dict
end

function loadMediumDict()
    global generateDictFile

    if isfile(mediumDictFile)
        try
            dict = loadMediumDictFile()
        catch
            println("... Regenerate media dictionary")
            include(generateDictFile)
            dict = loadMediumDictFile()
        end
    else
        include(generateDictFile)
        dict = loadMediumDictFile()
    end

    return dict
end

"""
    Medium = getMedium(name::AbstractString)

Return `Medium` object from medium `name`.
Possible values of argument `name` can be inquired via `listMedia()`
([Available media](@ref)).

# Examples
```julia
medium = getMedium("SimpleAir")
```

"""
function getMedium(name::AbstractString)::AbstractMedium
    global mediumDict
    if typeof(mediumDict[1]) == Nothing
        mediumDict[1] = loadMediumDict()
    end
    return mediumDict[1][name]
end

"""
    listMedia()

List available media of ModiaMedia.
"""
function listMedia()::Nothing
    global mediumDict
    if typeof(mediumDict[1]) == Nothing
        mediumDict[1] = loadMediumDict()
    end

    dict = mediumDict[1]
    media_table = DataFrames.DataFrame(name=AbstractString[], type=Symbol[])

    for key in sort(collect(keys(dict)))
        push!(media_table, [key, Symbol(typeof(dict[key]))])
    end

    println("\nMedia available in ModiaMedia:\n")
    show(media_table,allrows=true,allcols=true,summary=false)
    return nothing
end

# Import packages that are used in examples and tests
# (in order that there are no requirements on the environment
#  in which the examples and tests are executed).
import JSON
import ModiaMath
import Serialization
import StaticArrays
import Unitful
import Test

end # module


text:
@testset "Torsion" begin
  # Trivial torsion module

  A = diagonal_matrix(fmpq[1, 1])
  T = Hecke.TorQuadMod(A)
  @test order(T) == 1

  D4_gram = matrix(ZZ, [[2, 0, 0, -1],
                        [0, 2, 0, -1],
                        [0, 0, 2, -1],
                        [-1 ,-1 ,-1 ,2]])

  L = Zlattice(gram = D4_gram)
  T = @inferred discriminant_group(L)
  @test elem_type(T) == typeof(gens(T)[1])
  @test order(T) == 4
  @test elementary_divisors(T) == fmpz[2, 2]

  q1 = discriminant_group(root_lattice(:D,4))
  q2 = discriminant_group(Zlattice(gram=ZZ[0 2; 2 0]))
  @test Hecke.gram_matrix_quadratic(q1) != Hecke.gram_matrix_quadratic(q2)
  @test Hecke.gram_matrix_bilinear(q1) == Hecke.gram_matrix_bilinear(q2)

  @test sprint(show, q1) isa String
  @test sprint(show, "text/plain", q1) isa String
  @test sprint(show, gens(q1)[1]) isa String
  a,b = gens(q1)
  @test lift(inner_product(a,b)) == 1//2
  @test order(a) == 2
  @test order(0*a) == 1
  set_attribute!(q1, :name, "q1")
  f = hom(q1,q1, ZZ[2 0; 0 1])
  @test sprint(show, f) isa String

  @test b == preimage(f,b)
  @test_throws ErrorException preimage(f,a)
  @test !isbijective(f)

  T, i = primary_part(q1, 3)
  @test order(T) == 1

  Ld = dual(L)

  T = @inferred torsion_quadratic_module(Ld, L, snf = false)
  @test ngens(T) == 4
  G = gens(T)
  for i in 1:4
    @test T(lift(G[i])) == G[i]
  end

  TT, mTT = @inferred sub(T, [T([1, 1//2, 1//2, 1])])
  @test order(TT) == 2

  L = Zlattice(ZZ[1 0; 0 1])
  M = lattice(ambient_space(L), ZZ[2 0; 0 3])
  T = @inferred torsion_quadratic_module(L, M, gens = [[1, 1]])
  @test order(T) == 6

  @test_throws ArgumentError torsion_quadratic_module(L, M, gens = [[0, 0]])
  @test_throws ArgumentError torsion_quadratic_module(L, M, gens = [[1//2, 0]])
  @test_throws ArgumentError torsion_quadratic_module(L, M, gens = [[1, 1, 2]])
  @test_throws ArgumentError torsion_quadratic_module(L, lattice(ambient_space(L), QQ[1//2 0; 0 0]))

  #primary part of a TorQuadMod
  L = Zlattice(matrix(ZZ, [[2,0,0],[0,2,0],[0,0,2]]))
  T = Hecke.discriminant_group(L)
  @test basis_matrix(Hecke.cover(Hecke.primary_part(T,fmpz(2))[1])) == matrix(QQ, 3, 3, [1//2, 0, 0, 0, 1//2, 0, 0, 0, 1//2])
  L1 = Zlattice(identity_matrix(ZZ, 3))
  T1 = torsion_quadratic_module((1//6)*L1, L1)
  @test gram_matrix(Hecke.cover(Hecke.primary_part(T1,fmpz(2))[1])) == matrix(QQ, 3, 3, [1//4, 0, 0, 0, 1//4, 0, 0, 0, 1//4])
  @test ambient_space(Hecke.cover(Hecke.primary_part(T1, exponent(T1))[1]))==ambient_space(Hecke.cover(T1))
  @test Hecke.cover(Hecke.primary_part(T1, exponent(T1))[1]) == Hecke.cover(T1)
  @test gram_matrix(Hecke.cover(Hecke.primary_part(T1, exponent(T1))[1])) == gram_matrix(Hecke.cover(T1))

  #orthogonal submodule to a TorQuadMod
  L = Zlattice(matrix(ZZ, [[2,0,0],[0,2,0],[0,0,2]]))
  T = Hecke.discriminant_group(L)
  S, i = sub(T, gens(T))
  @test all([preimage(i,i(s))==s for s in gens(S)])
  @test basis_matrix(Hecke.cover(Hecke.orthogonal_submodule_to(T, S)[1])) == basis_matrix(L)
  L1 = Zlattice(identity_matrix(ZZ,10))
  T1 = torsion_quadratic_module(L1, 3*L1)
  S1, _ = sub(T1, gens(T1)[1:5])
  @test ambient_space(Hecke.cover(Hecke.orthogonal_submodule_to(T1, S1)[1])) == ambient_space(L1)

  #checks if a TorQuadMod is degenerate
  @test Hecke.isdegenerate(T) == false
  t = Hecke.TorQuadMod(matrix(QQ,1,1,[1//27]))
  d = sub(t, gens(t)*3)[1]
  @test Hecke.isdegenerate(d) == true

  #test for rescaled torsion quadratic module
  @test Hecke.gram_matrix_quadratic(Hecke.rescale(T, -1)) == matrix(QQ, 3, 3, [7//4,0,0,0,7//4,0,0,0,7//4])
  t = Hecke.TorQuadMod(QQ[1//3 0; 0 1//9])
  @test Hecke.gram_matrix_quadratic(Hecke.rescale(t, -1)) == matrix(QQ, 2, 2, [2//3,0,0,8//9])
  #This form is defined modulo `2`
  @test Hecke.gram_matrix_quadratic(Hecke.rescale(t, 2)) == matrix(QQ, 2, 2, [2//3,0,0,2//9])
  #The next form is defined modulo `3`
  @test Hecke.gram_matrix_quadratic(Hecke.rescale(t, 3)) == matrix(QQ, 2, 2, [1,0,0,1//3])
  #The next form is defined modulo `4`
  @test Hecke.gram_matrix_quadratic(Hecke.rescale(t, 4)) == matrix(QQ, 2, 2, [4//3,0,0,4//9])

  #test for normal form
  L1 = Zlattice(gram=matrix(ZZ, [[-2,0,0],[0,1,0],[0,0,4]]))
  T1 = Hecke.discriminant_group(L1)
  @test Hecke.gram_matrix_quadratic(normal_form(T1)[1]) == matrix(QQ, 2, 2, [1//2,0,0,1//4])

  L1 = Zlattice(gram=ZZ[-4 0 0; 0 4 0; 0 0 -2])
  AL1 = discriminant_group(L1)
  L2 = Zlattice(gram=ZZ[-4 0 0; 0 -4 0; 0 0 2])
  AL2 = discriminant_group(L2)
  n1 = normal_form(AL1)[1]
  g1 = QQ[1//2   0   0;
            0 1//4   0;
            0   0 5//4]
  @test Hecke.gram_matrix_quadratic(n1) == g1
  n2 = normal_form(AL2)[1]
  @test Hecke.gram_matrix_quadratic(n2) == g1
  L3 = Zlattice(gram=matrix(ZZ, [[2,0,0,-1],[0,2,0,-1],[0,0,2,-1],[-1,-1,-1,2]]))
  T=torsion_quadratic_module((1//6)*dual(L3), L3)
  n3 = normal_form(T)[1]
  g3 = QQ[1//6 1//12      0     0     0     0     0     0;
          1//12 1//6      0     0     0     0     0     0;
            0     0   1//12 1//24     0     0     0     0;
            0     0   1//24 1//12     0     0     0     0;
            0     0       0     0   2//9    0     0     0;
            0     0       0     0     0    2//9   0     0;
            0     0       0     0     0     0    2//9   0;
            0     0       0     0     0     0     0   2//9]
  @test Hecke.gram_matrix_quadratic(n3) == g3
  T2 = torsion_quadratic_module((1//6)*dual(L3), L3, modulus=fmpq(1//36))
  n4 = normal_form(T2)[1]
  g4 = QQ[1//36 1//72;
          1//72 1//36]
  @test Hecke.gram_matrix_quadratic(n4) == g4

  #test for brown invariant
  L1 = Zlattice(gram=matrix(ZZ, [[2,-1,0,0],[-1,2,-1,-1],[0,-1,2,0],[0,-1,0,2]]))
  T1 = discriminant_group(L1)
  @test Hecke.brown_invariant(T1) == 4
  L2 = Zlattice(matrix(ZZ, 2,2,[4,2,2,4]))
  T2 = Hecke.discriminant_group(L2)
  @test Hecke.brown_invariant(T2) == 2
  L3 = Zlattice(gram=matrix(ZZ, [[1,0,0],[0,1,0],[0,0,1]]))
  T3 = torsion_quadratic_module((1//10)*L3, L3)
  @test_throws ArgumentError Hecke.brown_invariant(T3)

  #test for genus
  L = Zlattice(gram=diagonal_matrix(fmpz[1,2,3,4]))
  D = discriminant_group(L)
  @test genus(D, (4,0)) == genus(L)
  L1 = Zlattice(gram=matrix(ZZ, [[2, -1, 0, 0, 0, 0],[-1, 2, -1, -1, 0, 0],[0, -1, 2, 0, 0, 0],[0, -1, 0, 2, 0, 0],[0, 0, 0, 0, 6, 3],[0, 0, 0, 0, 3, 6]]))
  T1 = discriminant_group(L1)
  @test genus(T1, (6,0)) == genus(L1)

  #test for isgenus
  L = Zlattice(gram=diagonal_matrix(fmpz[1,2,3,4]))
  D = discriminant_group(L)
  @test isgenus(D, (4,0))
  L1 = Zlattice(gram=matrix(ZZ, [[2, -1, 0, 0, 0, 0],[-1, 2, -1, -1, 0, 0],[0, -1, 2, 0, 0, 0],[0, -1, 0, 2, 0, 0],[0, 0, 0, 0, 6, 3],[0, 0, 0, 0, 3, 6]]))
  T1 = discriminant_group(L1)
  @test isgenus(T1, (6,0)) == true
  @test isgenus(T1, (4,2)) == false
  @test isgenus(T1, (16,2)) == true
  @test isgenus(T1, (5,1)) == false
  G = genus(diagonal_matrix(fmpz[2, 6, 6]))
  D = discriminant_group(G)
  @test isgenus(D, (2,0)) == false
  @test isgenus(D, (3,0)) == true

  N, i = normal_form(D)
  @test N === normal_form(N)[1]
  j = inv(i)
  @test all(g == i(j(g)) for g in gens(N))
end


text:
"""
$(SIGNATURES)

Write, to the current working directory, a .tex file with metadata (table name and description, line numbers and descriptions, and table notes) for the NIPA table contained in `bnt`.

Arguments
---------
* `bnt` -- a [`BeaNipaTable`](@ref) object

"""
function table_metadata_tex(bnt::BeaNipaTable)
    tablenum = bnt.tablenum
    tableid = bnt.tableid
    tabledesc = bnt.tabledesc
    linedesc = bnt.linedesc
    tabnotes = bnt.tablenotes

    tablename = DataFrames.makeidentifier(string("nipa", tablenum))
    fname = string(tablename, ".tex")

    f = open(fname, "w") do f
        write(f, "\\begin{longtable}{lp{5in}} \n")
        write(f, "\\toprule \n")
        write(f, string("Table & ", tablenum, " \\\\ \n"))
        write(f, string("Table ID & ", tableid, " \\\\ \n"))
        write(f, string("Description & ", tabledesc, " \\\\ \n"))
        write(f, "\\midrule \n")
        write(f, "Line No. & Description \\\\ \n")
        write(f, "\\midrule \\endhead \n")
        for (k, v) in zip(keys(linedesc), values(linedesc))
            write(f, string(k, " & ", v, " \\\\ \n"))
        end
        if tabnotes != ""
            write(f, "& \\\\ \n")
            write(f, "\\midrule \n")
            write(f, "Note & Description \\\\ \n")
            write(f, "\\midrule \n")
            for (k, v) in zip(keys(tabnotes), values(tabnotes))
                write(f, string(k, " & ", v, " \\\\ \n"))
            end
        end
        write(f, "\\bottomrule \n")
        write(f, "\\end{longtable} \n\n")
    end
end


text:
layers(::Type{<:CHELSA{<:Future{BioClim}}}) = layers(BioClim)
layerkeys(T::Type{<:CHELSA{<:Future{BioClim}}}, args...) = layerkeys(BioClim, args...)

layers(::Type{<:CHELSA{<:Future{Climate}}}) = (:prec, :temp, :tmin, :tmax)

date_step(::Type{<:CHELSA{<:Future{Climate,CMIP5}}}) = Year(20)
date_step(::Type{<:CHELSA{<:Future{Climate,CMIP6}}}) = Year(30)

# A modified key is used in the file name, while the key is used as-is in the path
const CHELSAKEY = (prec="pr", temp="tas", tmin="tasmin", tmax="tasmax", bio="bio")

"""
    getraster(T::Type{CHELSA{Future{BioClim}}}, [layer]; date) => String

Download CHELSA [`BioClim`](@ref) data, choosing layers from: `$(layers(CHELSA{BioClim}))`.

See the docs for [`Future`](@ref) for model choices.

Without a layer argument, all layers will be downloaded, and a `NamedTuple` of paths
returned.

## Keywords

- `date`: a `Date` or `DateTime` object, a Vector, or Tuple of start/end dates.
    Note that CHELSA CMIP5 only has two datasets, for the periods 2041-2060 and
    2061-2080. Dates must fall within these ranges.
"""
function getraster(
    T::Type{<:CHELSA{<:Future{BioClim}}}, layers::Union{Tuple,Int,Symbol}; date
)
    _getraster(T, layers, date)
end

"""
    getraster(T::Type{CHELSA{Future{Climate}}}, [layer]; date, month) => String

Download CHELSA [`Climate`](@ref) data, choosing layers from: `$(layers(CHELSA{BioClim}))`.

See the docs for [`Future`](@ref) for model choices.

Without a layer argument, all layers will be downloaded, and a `NamedTuple` of paths returned.

## Keywords

- `date`: a `Date` or `DateTime` object, a Vector, or Tuple of start/end dates.
    Note that CHELSA CMIP5 only has two datasets, for the periods 2041-2060 and
    2061-2080. Dates must fall within these ranges.
- `month`: the month of the year, from 1 to 12, or a array or range of months like `1:12`.
"""
function getraster(
    T::Type{<:CHELSA{<:Future{Climate}}}, layers::Union{Tuple,Symbol}; date, month
)
    _getraster(T, layers, date, month)
end

function _getraster(T::Type{<:CHELSA{<:Future{Climate}}}, layers, date, months::AbstractArray)
    map(month -> _getraster(T, layers, date, month), months)
end
function _getraster(
    T::Type{<:CHELSA{<:Future{Climate}}}, layers, dates::AbstractArray, months::AbstractArray
)
    map(date -> _getraster(T, layers, date, months), dates)
end
function _getraster(T::Type{<:CHELSA{<:Future{Climate}}}, layers, dates::AbstractArray, month)
    map(date -> _getraster(T, layers; date, month), dates)
end
function _getraster(T::Type{<:CHELSA{<:Future{Climate}}}, layers, dates::Tuple, months::AbstractArray)
    _getraster(T, layers, date_sequence(T, dates), months)
end
function _getraster(T::Type{<:CHELSA{<:Future{Climate}}}, layers, dates::Tuple, month)
    _getraster(T, layers, date_sequence(T, dates), month)
end
function _getraster(T::Type{<:CHELSA{<:Future{Climate}}}, layers, date, month)
    _getraster(T, layers; date, month)
end
function _getraster(T::Type{<:CHELSA{<:Future{BioClim}}}, layers, dates::AbstractArray)
    map(date -> _getraster(T, layers, date), dates)
end
function _getraster(T::Type{<:CHELSA{<:Future{BioClim}}}, layers, date::TimeType)
    _getraster(T, layers; date)
end

# We have the extra args as keywords again to generalise rasterpath/rasterurl
function _getraster(T::Type{<:CHELSA{<:Future}}, layers::Tuple; kw...)
    _map_layers(T, layers; kw...)
end
function _getraster(T::Type{<:CHELSA{<:Future}}, layer::Union{Symbol,Integer}; kw...)
    _check_layer(T, layer)
    path = rasterpath(T, layer; kw...)
    url = rasterurl(T, layer; kw...)
    return _maybe_download(url, path)
end

function rastername(T::Type{<:CHELSA{<:Future}}, layer; kw...)
    _rastername(_phase(T), T, layer; kw...)
end

function _rastername(
    ::Type{CMIP5}, T::Type{<:CHELSA{<:Future{BioClim}}}, layer::Integer; date
)
    date_string = _date_string(_phase(T), date)
    mod = _format(CHELSA, _model(T))
    scen = _format(CHELSA, _scenario(T))
    return "CHELSA_bio_mon_$(mod)_$(scen)_r1i1p1_g025.nc_$(layer)_$(date_string)_V1.2.tif"
end
function _rastername(
    ::Type{CMIP5}, T::Type{<:CHELSA{<:Future{Climate}}}, layer::Symbol; date, month
)
    date_string = _date_string(_phase(T), date)
    mod = _format(CHELSA, _model(T))
    scen = _format(CHELSA, _scenario(T))
    key = CHELSAKEY[layer]
    suffix = layer === :prec ? "" : "_V1.2" # prec filenames dont end in _V1.2
    return "CHELSA_$(key)_mon_$(mod)_$(scen)_r1i1p1_g025.nc_$(month)_$(date_string)$(suffix).tif"
end
function _rastername(::Type{CMIP6}, T::Type{<:CHELSA{<:Future{BioClim}}}, layer::Integer; date)
    date_string = _date_string(_phase(T), date)
    mod = _format(CHELSA, _model(T))
    scen = _format(CHELSA, _scenario(T))
    return "CHELSA_bio$(layer)_$(date_string)_$(mod)_$(scen)_V.2.1.tif"
end
function _rastername(
    ::Type{CMIP6}, T::Type{<:CHELSA{<:Future{Climate}}}, layer::Symbol; date, month
)
    # CMIP6 Climate uses an underscore in the date string, of course
    date_string = replace(_date_string(_phase(T), date), "-" => "_")
    mod = _format(CHELSA, _model(T))
    scen = _format(CHELSA, _scenario(T))
    key = CHELSAKEY[layer]
    mon = lpad(month, 2, '0')
    return "CHELSA_$(mod)_r1i1p1f1_w5e5_$(scen)_$(key)_$(mon)_$(date_string)_norm.tif"
end

function rasterpath(T::Type{<:CHELSA{<:Future}})
    joinpath(rasterpath(CHELSA), "Future", string(_dataset(T)), string(_scenario(T)), string(_model(T)))
end
function rasterpath(T::Type{<:CHELSA{<:Future}}, layer; kw...)
    joinpath(rasterpath(T), rastername(T, layer; kw...))
end

function rasterurl(T::Type{<:CHELSA{<:Future}}, layer; date, kw...)
    date_str = _date_string(_phase(T), date)
    key = _chelsa_layer(_dataset(T), layer)
    path = _urlpath(_phase(T), T::Type{<:CHELSA{<:Future}}, key, date_str)
    joinpath(rasterurl(CHELSA), path, rastername(T, layer; date, kw...))
end

_chelsa_layer(::Type{<:BioClim}, layer) = :bio
_chelsa_layer(::Type{<:Climate}, layer) = layer

function _urlpath(::Type{CMIP5}, T::Type{<:CHELSA{<:Future}}, name, date_str)
    return "chelsa_V1/cmip5/$date_str/$name/"
end
function _urlpath(::Type{CMIP6}, T::Type{<:CHELSA{<:Future}}, name, date_str)
    # The model is in uppercase in the URL for CMIP6
    mod = uppercase(_format(CHELSA, _model(T)))
    scen = _format(CHELSA, _scenario(T))
    key = CHELSAKEY[name]
    return "chelsa_V2/GLOBAL/climatologies/$date_str/$mod/$scen/$key/"
end

function _date_string(::Type{CMIP5}, date)
    if date < DateTime(2041)
        _cmip5_date_error(date)
    elseif date < DateTime(2061)
        "2041-2060"
    elseif date < DateTime(2081)
        "2061-2080"
    else
        _cmip5_date_error(date)
    end
end

function _date_string(::Type{CMIP6}, date)
    if date < DateTime(1981)
        _cmip6_date_error(date)
    elseif date < DateTime(2011)
        "1981-2010"
    elseif date < DateTime(2041)
        "2011-2040"
    elseif date < DateTime(2071)
        "2041-2070"
    elseif date < DateTime(2101)
        "2071-2100"
    else
        _cmip6_date_error(date)
    end
end

_cmip5_date_error(date) = error("CMIP5 covers the period from 2041-2080, not including $date")
_cmip6_date_error(date) = error("CMIP6 covers the period from 1981-2100, not including $date")

_dataset(::Type{<:CHELSA{F}}) where F<:Future = _dataset(F)
_phase(::Type{<:CHELSA{F}}) where F<:Future = _phase(F)
_model(::Type{<:CHELSA{F}}) where F<:Future = _model(F)
_scenario(::Type{<:CHELSA{F}}) where F<:Future = _scenario(F)

# Climate model string formatters for CHELSA Future

# CMIP5
_format(::Type{CHELSA}, ::Type{ACCESS1}) = "ACCESS1-0"
_format(::Type{CHELSA}, ::Type{BNUESM}) = "BNU-ESM"
_format(::Type{CHELSA}, ::Type{CCSM4}) = "CCSM4"
_format(::Type{CHELSA}, ::Type{CESM1BGC}) = "CESM1-BGC"
_format(::Type{CHELSA}, ::Type{CESM1CAM5}) = "CESM1-CAM5"
_format(::Type{CHELSA}, ::Type{CMCCCMS}) = "CMCC-CMS"
_format(::Type{CHELSA}, ::Type{CMCCCM}) = "CMCC-CM"
_format(::Type{CHELSA}, ::Type{CNRMCM5}) = "CNRM-CM5"
_format(::Type{CHELSA}, ::Type{CSIROMk3}) = "CSIRO-Mk3"
_format(::Type{CHELSA}, ::Type{CanESM2}) = "CanESM2"
_format(::Type{CHELSA}, ::Type{FGOALS}) = "FGOALS-g2"
_format(::Type{CHELSA}, ::Type{FIOESM}) = "FIO-ESM"
_format(::Type{CHELSA}, ::Type{GFDLCM3}) = "GFDL-CM3"
_format(::Type{CHELSA}, ::Type{GFDLESM2G}) = "GFDL-ESM2G"
_format(::Type{CHELSA}, ::Type{GFDLESM2M}) = "GFDL-ESM2M"
_format(::Type{CHELSA}, ::Type{GISSE2HCC}) = "GISS-E2-H-CC"
_format(::Type{CHELSA}, ::Type{GISSE2H}) = "GISS-E2-H"
_format(::Type{CHELSA}, ::Type{GISSE2RCC}) = "GISS-E2-R-CC"
_format(::Type{CHELSA}, ::Type{GISSE2R}) = "GISS-E2-R"
_format(::Type{CHELSA}, ::Type{HadGEM2AO}) = "HadGEM2-AO"
_format(::Type{CHELSA}, ::Type{HadGEM2CC}) = "HadGEM2-CC"
_format(::Type{CHELSA}, ::Type{IPSLCM5ALR}) = "IPSL-CM5A-LR"
_format(::Type{CHELSA}, ::Type{IPSLCM5AMR}) = "IPSL-CM5A-MR"
_format(::Type{CHELSA}, ::Type{MIROCESMCHEM}) = "MIROC-ESM-CHEM"
_format(::Type{CHELSA}, ::Type{MIROCESM}) = "MIROC-ESM"
_format(::Type{CHELSA}, ::Type{MIROC5}) = "MIROC5"
_format(::Type{CHELSA}, ::Type{MPIESMLR}) = "MPI-ESM-LR"
_format(::Type{CHELSA}, ::Type{MPIESMMR}) = "MPI-ESM-MR"
_format(::Type{CHELSA}, ::Type{MRICGCM3}) = "MRI-CGCM3"
_format(::Type{CHELSA}, ::Type{MRIESM1}) = "MRI-ESM1"
_format(::Type{CHELSA}, ::Type{NorESM1M}) = "NorESM1-M"
_format(::Type{CHELSA}, ::Type{BCCCSM1}) = "bcc-csm-1"
_format(::Type{CHELSA}, ::Type{Inmcm4}) = "inmcm4"

# CMIP6
_format(::Type{CHELSA}, ::Type{GFDLESM4}) = "gfdl-esm4"
_format(::Type{CHELSA}, ::Type{IPSLCM6ALR}) = "ipsl-cm6a-lr"
_format(::Type{CHELSA}, ::Type{MPIESMHR}) = "mpi-esm1-2-hr"
_format(::Type{CHELSA}, ::Type{MRIESM2}) = "mri-esm2-0"
_format(::Type{CHELSA}, ::Type{UKESM}) = "ukesm1-0-ll"

# Format scenarios
_format(::Type{CHELSA}, ::Type{RCP26}) = "rcp26"
_format(::Type{CHELSA}, ::Type{RCP45}) = "rcp45"
_format(::Type{CHELSA}, ::Type{RCP60}) = "rcp60"
_format(::Type{CHELSA}, ::Type{RCP85}) = "rcp85"

_format(::Type{CHELSA}, ::Type{SSP126}) = "ssp126"
_format(::Type{CHELSA}, ::Type{SSP245}) = "ssp245"
_format(::Type{CHELSA}, ::Type{SSP370}) = "ssp370"
_format(::Type{CHELSA}, ::Type{SSP585}) = "ssp585"

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
include("/Users/guilhermebodin/Documents/PUC/MIPTests_Branch_Guilherme/runmiptests.jl")
include("/Users/guilhermebodin/Documents/PUC/MIPSolver2017_04.jl/src/branch_and_bound.jl")

runtests(SolveMIP, GurobiSolver())