text:
"""
    struct MeshData{Dim, Tv, Ti}

MeshData: contains info for a high order piecewise polynomial discretization on an
unstructured mesh.

Example:
```julia
N, K1D = 3, 2
rd = RefElemData(Tri(), N)
VXY, EToV = uniform_mesh(Tri(), K1D)
md = MeshElemData(VXY, EToV, rd)
@unpack x, y = md
```
"""
Base.@kwdef struct MeshData{Dim, Tv, Ti, DimTimesDim}

    # num_elements::Ti            # number of elements
    VXYZ::NTuple{Dim, Vector{Tv}}  # vertex coordinates
    EToV::Matrix{Ti}              # mesh vertex array
    FToF::Matrix{Ti}              # face connectivity

    xyz::NTuple{Dim, Matrix{Tv}}   # physical points
    xyzf::NTuple{Dim, Matrix{Tv}}  # face nodes
    xyzq::NTuple{Dim, Matrix{Tv}}  # phys quad points, Jacobian-scaled weights
    wJq::Matrix{Tv}

    # arrays of connectivity indices between face nodes
    mapM::Matrix{Ti}
    mapP::Matrix{Ti}
    mapB::Vector{Ti}

    # volume geofacs Gij = dx_i/dxhat_j
    rstxyzJ::SMatrix{Dim, Dim, Matrix{Tv}, DimTimesDim}
    J::Matrix{Tv}

    # surface geofacs
    nxyzJ::NTuple{Dim, Matrix{Tv}}
    Jf::Matrix{Tv}

    is_periodic::NTuple{Dim, Bool}
end

function Base.show(io::IO, md::MeshData{DIM}) where {DIM}
    @nospecialize md
    print(io,"MeshData{$DIM}")
end
function Base.show(io::IO, ::MIME"text/plain", md::MeshData{DIM}) where {DIM}
    @nospecialize md
    print(io,"MeshData of dimension $DIM with $(md.K) elements")
end

# enable use of @set and setproperties(...) for MeshData
ConstructionBase.constructorof(::Type{MeshData{A, B, C}}) where {A, B, C} = MeshData{A, B, C}

function Base.propertynames(x::MeshData{1}, private::Bool = false)
    return (fieldnames(MeshData)...,
            :num_elements, :VX, :x, :xq, :xf, :nxJ, :rxJ)
end
function Base.propertynames(x::MeshData{2}, private::Bool = false)
    return (fieldnames(MeshData)...,
            :num_elements, :VX, :VY, :x, :y, :xq, :yq, :xf, :yf,
            :nxJ, :nyJ, :rxJ, :sxJ, :ryJ, :syJ)
end
function Base.propertynames(x::MeshData{3}, private::Bool = false)
    return (fieldnames(MeshData)...,
            :num_elements, :VX, :VY, :VZ, :x, :y, :z, :xq, :yq, :zq, :xf, :yf, :zf,
            :nxJ, :nyJ, :nzJ, :rxJ, :sxJ, :txJ, :ryJ, :syJ, :tyJ, :rzJ, :szJ, :tzJ)
end

# convenience routines for unpacking individual tuple entries
function Base.getproperty(x::MeshData, s::Symbol)

    if s==:VX
        return getfield(x, :VXYZ)[1]
    elseif s==:VY
        return getfield(x, :VXYZ)[2]
    elseif s==:VZ
        return getfield(x, :VXYZ)[3]

    elseif s==:x
        return getfield(x, :xyz)[1]
    elseif s==:y
        return getfield(x, :xyz)[2]
    elseif s==:z
        return getfield(x, :xyz)[3]

    elseif s==:xq
        return getfield(x, :xyzq)[1]
    elseif s==:yq
        return getfield(x, :xyzq)[2]
    elseif s==:zq
        return getfield(x, :xyzq)[3]

    elseif s==:xf
        return getfield(x, :xyzf)[1]
    elseif s==:yf
        return getfield(x, :xyzf)[2]
    elseif s==:zf
        return getfield(x, :xyzf)[3]

    elseif s==:nxJ
        return getfield(x, :nxyzJ)[1]
    elseif s==:nyJ
        return getfield(x, :nxyzJ)[2]
    elseif s==:nzJ
        return getfield(x, :nxyzJ)[3]

    elseif s==:rxJ
        return getfield(x, :rstxyzJ)[1,1]
    elseif s==:sxJ
        return getfield(x, :rstxyzJ)[1,2]
    elseif s==:txJ
        return getfield(x, :rstxyzJ)[1,3]
    elseif s==:ryJ
        return getfield(x, :rstxyzJ)[2,1]
    elseif s==:syJ
        return getfield(x, :rstxyzJ)[2,2]
    elseif s==:tyJ
        return getfield(x, :rstxyzJ)[2,3]
    elseif s==:rzJ
        return getfield(x, :rstxyzJ)[3,1]
    elseif s==:szJ
        return getfield(x, :rstxyzJ)[3,2]
    elseif s==:tzJ
        return getfield(x, :rstxyzJ)[3,3]
    elseif s==:K || s==:num_elements # old behavior where K = num_elements
        return size(getfield(x, :EToV), 1)

    elseif s==:sJ
        return getfield(x, :Jf)

    # return getfield(x,:num_elements) # num rows in EToV = num elements
    else
        return getfield(x, s)
    end
end

"""
    MeshData(VXYZ, EToV, rd::RefElemData)

Returns a MeshData struct with high order DG mesh information from the unstructured
mesh information (VXYZ..., EToV).

    MeshData(md::MeshData, rd::RefElemData, xyz...)

Given new nodal positions `xyz...` (e.g., from mesh curving), recomputes geometric terms
and outputs a new MeshData struct. Only fields modified are the coordinate-dependent terms
    `xyz`, `xyzf`, `xyzq`, `rstxyzJ`, `J`, `nxyzJ`, `sJ`.
"""

# splats VXYZ
MeshData(VXYZ::T, EToV, rd::RefElemData{NDIMS}) where {NDIMS, T <: NTuple{NDIMS}} =
    MeshData(VXYZ..., EToV, rd)

function MeshData(VX::AbstractVector{Tv}, EToV, rd::RefElemData{1}) where {Tv}

    # Construct global coordinates
    @unpack V1 = rd
    x = V1*VX[transpose(EToV)]
    K = size(EToV,1)
    Nfaces = 2

    FToF = zeros(Int,Nfaces,K)
    sk = 1
    for e = 1:K
        l = 2*e-1
        r = 2*e
        FToF[1:2,e] .= [l-1; r+1]
        sk += 1
    end
    FToF[1,1] = 1
    FToF[Nfaces,K] = Nfaces*K

    # Connectivity maps
    @unpack Vf = rd
    xf = Vf*x
    mapM = reshape(1:2*K,2,K)
    mapP = copy(mapM)
    mapP[1,2:end] .= mapM[2,1:end-1]
    mapP[2,1:end-1] .= mapM[1,2:end]
    mapB = findall(@. mapM[:]==mapP[:])

    # Geometric factors and surface normals
    J = repeat(transpose(diff(VX)/2),length(rd.r),1)
    rxJ = one.(J)
    nxJ = repeat([-1.0; 1.0],1,K)
    sJ = abs.(nxJ)

    @unpack Vq,wq = rd
    xq = Vq*x
    wJq = diagm(wq)*(Vq*J)

    is_periodic = (false,)
    return MeshData(tuple(VX),EToV,FToF,
                    tuple(x),tuple(xf),tuple(xq),wJq,
                    collect(mapM),mapP,mapB,
                    SMatrix{1,1}(tuple(rxJ)),J,
                    tuple(nxJ),sJ,
                    is_periodic)

end

function MeshData(VX, VY, EToV, rd::RefElemData{2})

    @unpack fv = rd
    FToF = connect_mesh(EToV,fv)
    Nfaces,K = size(FToF)

    #Construct global coordinates
    @unpack V1 = rd
    x = V1*VX[transpose(EToV)]
    y = V1*VY[transpose(EToV)]

    #Compute connectivity maps: uP = exterior value used in DG numerical fluxes
    @unpack Vf = rd
    xf = Vf*x
    yf = Vf*y
    mapM,mapP,mapB = build_node_maps(FToF,xf,yf)
    Nfp = size(Vf,1) ÷ Nfaces
    mapM = reshape(mapM,Nfp*Nfaces,K)
    mapP = reshape(mapP,Nfp*Nfaces,K)

    #Compute geometric factors and surface normals
    @unpack Dr,Ds = rd
    rxJ, sxJ, ryJ, syJ, J = geometric_factors(x, y, Dr, Ds)
    rstxyzJ = SMatrix{2,2}(rxJ,ryJ,sxJ,syJ)

    @unpack Vq,wq = rd
    xq,yq = (x->Vq*x).((x,y))
    wJq = diagm(wq)*(Vq*J)

    nxJ,nyJ,sJ = compute_normals(rstxyzJ,rd.Vf,rd.nrstJ...)

    is_periodic = (false,false)
    return MeshData(tuple(VX,VY),EToV,FToF,
                    tuple(x,y),tuple(xf,yf),tuple(xq,yq),wJq,
                    mapM,mapP,mapB,
                    SMatrix{2,2}(tuple(rxJ,ryJ,sxJ,syJ)),J,
                    tuple(nxJ,nyJ),sJ,
                    is_periodic)

end

function MeshData(VX, VY, VZ, EToV, rd::RefElemData{3})

    @unpack fv = rd
    FToF = connect_mesh(EToV,fv)
    Nfaces,K = size(FToF)

    #Construct global coordinates
    @unpack V1 = rd
    x,y,z = (x->V1*x[transpose(EToV)]).((VX,VY,VZ))

    #Compute connectivity maps: uP = exterior value used in DG numerical fluxes
    @unpack r,s,t,Vf = rd
    xf,yf,zf = (x->Vf*x).((x,y,z))
    mapM,mapP,mapB = build_node_maps(FToF,xf,yf,zf)
    Nfp = convert(Int,size(Vf,1)/Nfaces)
    mapM = reshape(mapM,Nfp*Nfaces,K)
    mapP = reshape(mapP,Nfp*Nfaces,K)

    #Compute geometric factors and surface normals
    @unpack Dr,Ds,Dt = rd
    rxJ,sxJ,txJ,ryJ,syJ,tyJ,rzJ,szJ,tzJ,J = geometric_factors(x,y,z,Dr,Ds,Dt)
    rstxyzJ = SMatrix{3,3}(rxJ,ryJ,rzJ,sxJ,syJ,szJ,txJ,tyJ,tzJ)

    @unpack Vq,wq = rd
    xq,yq,zq = (x->Vq*x).((x,y,z))
    wJq = diagm(wq)*(Vq*J)

    nxJ,nyJ,nzJ,sJ = compute_normals(rstxyzJ,rd.Vf,rd.nrstJ...)

    is_periodic = (false,false,false)
    return MeshData(tuple(VX,VY,VZ),EToV,FToF,
                    tuple(x,y,z),tuple(xf,yf,zf),tuple(xq,yq,zq),wJq,
                    mapM,mapP,mapB,
                    rstxyzJ,J,tuple(nxJ,nyJ,nzJ),sJ,
                    is_periodic)
end

MeshData(md::MeshData, rd::RefElemData, xyz...) = MeshData(rd, md, xyz...)

function MeshData(rd::RefElemData, md::MeshData{Dim}, xyz...) where {Dim}

    # compute new quad and plotting points
    xyzf = map(x->rd.Vf*x,xyz)
    xyzq = map(x->rd.Vq*x,xyz)

    #Compute geometric factors and surface normals
    geo = geometric_factors(xyz...,rd.Drst...)
    if Dim==1
        rstxyzJ = SMatrix{Dim,Dim}(geo[1])
    elseif Dim==2
        rstxyzJ = SMatrix{Dim,Dim}(geo[1],geo[3],
                                   geo[2],geo[4])
    elseif Dim==3
        rstxyzJ = SMatrix{Dim,Dim}(geo[1],geo[4],geo[7],
                                   geo[2],geo[5],geo[8],
                                   geo[3],geo[6],geo[9])
    end
    geof = compute_normals(rstxyzJ,rd.Vf,rd.nrstJ...)

    setproperties(md,(xyz=xyz,xyzq=xyzq,xyzf=xyzf,
                  rstxyzJ=rstxyzJ,J=last(geo),
                  nxyzJ=geof[1:Dim],Jf=last(geof)))
end

# physical normals are computed via G*nhatJ, G = matrix of geometric terms
function compute_normals(geo::SMatrix{Dim, Dim}, Vf, nrstJ...) where {Dim}
    nxyzJ = ntuple(x->zeros(size(Vf,1),size(first(geo),2)),Dim)
    for i = 1:Dim, j = 1:Dim
        nxyzJ[i] .+= (Vf*geo[i,j]).*nrstJ[j]
    end
    sJ = sqrt.(sum(map(x->x.^2,nxyzJ)))
    return nxyzJ...,sJ
end


text:
using InformationSensitivityFunctions
using Test

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

ext: jl


text:
function Base.show(io::IO, x::CLRObject)
    h = gethandle(x)
    if h == 0
        print(io, "null")
        return
    end
    type = clrtypeof(x)
    if isassignable(T"System.Type", type)
        print(io, string(x))
    else
        print(io, string(type))
        print(io, "(")
        unboxed = unbox(x)
        if typeof(unboxed) == CLRObject
            print(io, repr(string(x)))
        else
            print(io, repr(unboxed))
        end
        print(io, ")")
    end
end

function makegenericmethod(mi::CLRObject, args...)
    invokemember(T"System.Reflection.MethodInfo", mi, :MakeGenericMethod, args...)
end

function makegenerictype(ty::CLRObject, args...)
    invokemember(T"System.Type", ty, :MakeGenericType, args...)
end

function getmember(type::CLRObject, name)
    ret = invokemember(type, :GetMember, typeof(name) == Symbol ? string(name) : name)
    collect(ret)
end

function getmethod(type::CLRObject, name, argtypes...)
    arr = arrayof(T"System.Type", length(argtypes))
    for (i, t) in enumerate(argtypes)
        arraystore(arr, i - 1, t)
    end
    invokemember(type, :GetMethod, name, arr)
end


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

"""
    LocalKriging(var₁=>param₁, var₂=>param₂, ...)

LocalKriging estimation solver, where `var::Symbol` is the variable name
and `param` is a `NamedTuple` containing the parameters below:

## Parameters

* `variogram` - Reference variogram model
* `mean`      - Simple Kriging mean
* `method`    - LocalKriging method. :MovingWindows or :KernelConvolution
  (default to :MovingWindows)
* `localaniso`    - Local parameters of the domain
* `localanisohd`  - Local parameters of the samples. Only necessary for
  :KernelConvolution method. They are automatically passed via NN from
  `localaniso` if not informed.
* `minneighbors` - Minimum number of neighbors (default to 1)
* `maxneighbors` - Maximum number of neighbors (default to 20)
* `neighborhood` - Search neighborhood (default to `nothing`)
* `distance`     - Distance used to find nearest neighbors (default to `Euclidean()`)
"""
@estimsolver LocalKriging begin
  @param variogram = ExponentialVariogram()
  @param mean = nothing #0.0 or [1.0,0.8,....]
  @param method = :MovingWindows
  @param localaniso = nothing
  @param localanisohd = nothing
  @param minneighbors = 1
  @param maxneighbors = 20
  @param neighborhood = nothing
  @param distance = Euclidean()
end


text:
module Solve_TF_module
    using LinearAlgebra

    struct Solve_TF_param
        ELE_TOTAL::Int64
        INTEGTABLENUM::Int64
        NODE_TOTAL::Int64
        USETHREAD::Bool
        X_MAX::Float64
        X_MIN::Float64
        Y0::Float64
        YMAX::Float64
    end

    mutable struct Solve_TF_variables
        length::Array{Float64, 1}
        mat_A_glo::SymTridiagonal{Float64,Array{Float64,1}}
        node_num_seg::Array{Int64, 2}
        node_x_ele::Array{Float64, 2}
        node_x_glo::Array{Float64, 1}
        tmp::Array{Float64, 1}
        ug::Array{Float64, 1}
        vec_b_ele::Array{Float64, 2}
        vec_b_glo::Array{Float64, 1}
        w::Array{Float64, 1}
        x::Array{Float64, 1}
    end
end


text:
# Michelle Medina
# June 29, 2020

module SparseGrids

include("Types.jl")
include("Interpolation.jl")
include("1DBasis.jl")
include("Utilities.jl")
include("2DBasis.jl")
include("Derivatives.jl")
include("FiniteDifferencing.jl")

end # module


text:
abstract type AbstractSamples{T} <: AbstractVector{T} end

struct WeightedSamples{T} <: AbstractSamples{T}
    distances::Vector{T}
    weights::Vector{T}
end

struct UniformSamples{T} <: AbstractSamples{T}
    distances::Vector{T}
end

function WeightedSamples{T}(::UndefInitializer, N) where T
    WeightedSamples(Vector{T}(undef, N),Vector{T}(undef, N))
end
function WeightedSamples{T}(::UndefInitializer, N, nthreads) where T
    samples = Vector{WeightedSamples{T}}(undef, nthreads)
    Threads.@threads for thread ∈ 1:nthreads
        samples[thread] = WeightedSamples{T}(undef, N)
    end
    samples
end
function UniformSamples{T}(::UndefInitializer, N) where T
    UniformSamples(Vector{T}(undef, N))
end
function UniformSamples{T}(::UndefInitializer, N, nthreads) where T
    samples = Vector{UniformSamples{T}}(undef, nthreads)
    Threads.@threads for thread ∈ 1:nthreads
        samples[thread] = UniformSamples{T}(undef, N)
    end
    samples
end

Base.size(s::AbstractSamples) = size(s.distances)
Base.length(s::AbstractSamples) = length(s.distances)
Base.getindex(ws::WeightedSamples, i) = (ws.distances[i], ws.weights[i])
Base.getindex(us::UniformSamples, i) = us.distances[i]
function Base.setindex!(ws::WeightedSamples, (d,w), i)
    ws.distances[i] = d
    ws.weights[i] = w
    (d,w)
end
Base.setindex!(us::UniformSamples, d, i) = us.distances[i] = d
Base.IndexStyle(::Type{<:AbstractSamples}) = IndexLinear()
function Base.resize!(s::WeightedSamples, N)
    resize!(s.distances, N)
    resize!(s.weights)
    s
end
function Base.resize!(s::UniformSamples, N)
    resize!(s.distances, N)
    s
end

function KernelDensityDistributionEsimates.KDE(distances::UniformSamples)
    KDE(KernelDensityDistributionEsimates.kde(distances.distances,npoints=2048))
end
function KernelDensityDistributionEsimates.KDE(distances::WeightedSamples)
    KDE(KernelDensityDistributionEsimates.kde(distances.distances, weights = distances.weights,npoints=2048))
end

function silverman(data, α = 0.9)
    N = length(data)
    N <= 1 && return α
    var_width = std(data)
    q25, q75 = quantile(data, (0.25, 0.75))
    quant_width = (q75 - q25) * 0.7462686567164178

    width = min(var_width, quant_width)

    α * width * N^(-0.2)
end

function KernelDensity.kde(distances::UniformSamples)
    bw = silverman(distances.distances)
    if bw == 0
        x = distances.distances[1]
        return UnivariateKDE(x:1.0:x,[1.0])
    else
        return KernelDensity.kde( distances.distances, bandwidth = bw, npoints=2048 )
    end
end
function KernelDensity.kde(distances::WeightedSamples)
    bw = silverman(distances.distances)
    if bw == 0
        x = distances.distances[1]
        return UnivariateKDE(x:1.0:x,[1.0])
    else
        return KernelDensity.kde( distances.distances, bandwidth = bw, weights = distances.weights, npoints=2048 )
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

module LocallyWeightedRegression

using GeoStatsBase

using LinearAlgebra
using NearestNeighbors
using StaticArrays
using Distances
using Variography

import GeoStatsBase: solve

export LocalWeightRegress

"""
    LocalWeightRegress(var₁=>param₁, var₂=>param₂, ...)

Locally weighted regression (LOESS) estimation solver.

## Parameters

* `neighbors` - Number of neighbors (default to all data locations)
* `variogram` - A variogram defined in Variography.jl (default to `ExponentialVariogram()`)
* `distance`  - A distance defined in Distances.jl (default to `Euclidean()`)

### References

Cleveland 1979. *Robust Locally Weighted Regression and Smoothing Scatterplots*
"""
@estimsolver LocalWeightRegress begin
  @param neighbors = nothing
  @param variogram = ExponentialVariogram()
  @param distance = Euclidean()
end

function solve(problem::EstimationProblem, solver::LocalWeightRegress)
  # retrieve problem info
  pdata = data(problem)
  pdomain = domain(problem)

  # result for each variable
  μs = []; σs = []

  for (var,V) in variables(problem)
    # get user parameters
    if var ∈ keys(solver.params)
      varparams = solver.params[var]
    else
      varparams = LocalWeightRegressParam()
    end

    # get valid data for variable
    X, z = valid(pdata, var)

    # number of data points for variable
    ndata = length(z)

    @assert ndata > 0 "estimation requires data"

    # allocate memory
    varμ = Vector{V}(undef, npoints(pdomain))
    varσ = Vector{V}(undef, npoints(pdomain))

    # fit search tree
    kdtree = KDTree(X, varparams.distance)

    # determine number of nearest neighbors to use
    k = varparams.neighbors == nothing ? ndata : varparams.neighbors

    @assert k ≤ ndata "number of neighbors must be smaller or equal to number of data points"

    # determine kernel (or weight) function
    γ = varparams.variogram
    kern(x, y) = sill(γ) - γ(x, y)

    # pre-allocate memory for coordinates
    x = MVector{ndims(pdomain),coordtype(pdomain)}(undef)

    # estimation loop
    for location in LinearPath(pdomain)
      coordinates!(x, pdomain, location)

      inds, dists = knn(kdtree, x, k)

      Xₗ = [ones(eltype(X), k) X[:,inds]']
      zₗ = view(z, inds)

      Wₗ = Diagonal([kern(x, view(X,:,j)) for j in inds])

      # weighted least-squares
      θₗ = Xₗ'*Wₗ*Xₗ \ Xₗ'*Wₗ*zₗ

      # add intercept term to estimation location
      xₗ = [one(eltype(x)), x...]

      # linear combination of response values
      rₗ = Wₗ*Xₗ*(Xₗ'*Wₗ*Xₗ\xₗ)

      varμ[location] = θₗ ⋅ xₗ
      varσ[location] = norm(rₗ)
    end

    push!(μs, var => varμ)
    push!(σs, var => varσ)
  end

  EstimationSolution(pdomain, Dict(μs), Dict(σs))
end

end


text:
"""
Provider module for the `biber` executable.
"""
module Biber

export biber

using Pkg.Artifacts

binary() = joinpath(artifact"biber_bin", "biber$(Sys.iswindows() ? ".exe" : "")")

"""
    biber(f)

Call function `f` with the `biber` command available on `PATH` for the duration of `f`.

```julia
biber() do bin
    run(`\$bin file`)
end
```
"""
function biber(f)
    dir, bin = splitdir(binary())
    path = ENV["PATH"]
    pathsep = Sys.iswindows() ? ";" : ":"
    withenv("PATH" => "$dir$pathsep$path") do
        f(bin)
    end
end

version() = VersionNumber(lstrip(!isnumeric, biber(bin -> readchomp(`$bin --version`))))

end


text:
#=
prompt.jl
prompt user to hit a key to continue
=#

export prompt

import Plots # isplotnull, plot!
import REPL

prompt_state = :prompt # global prompt state

"""
    prompt( ; gui::Bool=true)

Prompt user to hit any key to continue, after `gui()`.
Some keys have special actions: `[q]uit [d]raw [n]odraw`.
Call `prompt(:prompt)` to revert to default.
"""
function prompt( ;
    gui::Bool=true,
    io_in::IO = stdin,
    kwargs...,
)

    global prompt_state

    prompt_state === :nodraw && return nothing

    gui && !Plots.isplotnull() && display(plot!()) # Plots.gui()

    (prompt_state === :draw) && return nothing

    # call wait_for_key only if interactive or testing
    if isinteractive() || io_in != stdin
        c = wait_for_key( ; io_in, kwargs...)
    else
        c = 'd'
    end

    (c == 'd') && (prompt_state = :draw)
    (c == 'n') && (prompt_state = :nodraw)
    (c == 'q') && throw("quit")

    return nothing
end

"""
    function wait_for_key( ; io_in = stdin, io_out = stdout, prompt=?)

From [discourse](https://discourse.julialang.org/t/wait-for-a-keypress/20218).
"""
function wait_for_key( ;
    io_in::IO = stdin,
    io_out::IO = stdout,
    prompt::String = "press any key [d]raw [n]odraw [q]uit : ",
)

    print(io_out, prompt)

    t = REPL.TerminalMenus.terminal
    REPL.Terminals.raw!(t, true)
    char = read(io_in, Char)
    REPL.Terminals.raw!(t, false)

    write(io_out, char)
    write(io_out, "\n")

    return char
end

"""
    prompt(key::symbol)

Set prompt state to one of:
- `:prompt` call `gui()` if possible, then prompt user.
- `:draw` call `gui()` if possible, then continue.
- `:nodraw` do not call `gui()`, just continue.

Use `prompt(:state)` to query current state.

Actually it calls `display(plot!())` instead of `gui()`.
"""
function prompt(key::Symbol)
    global prompt_state

    if key === :state
        return prompt_state
    end

    key ∉ (:prompt, :draw, :nodraw) && throw("prompt $prompt")
    prompt_state = key
end


text:
using SenseHat
using Test

using ColorTypes

const LED = led_matrix()

LED[:] = SenseHat.JULIA_LOGO
sleep(3)
led_clear()


text:
using Documenter, ScenTreesMakie

push!(LOAD_PATH,"../src/")
makedocs(
	sitename =  "ScenTreesMakie.jl",
	authors = "Ruben Schlotter, Alois Pichler, Kirui Kipngeno",
	format = Documenter.HTML(prettyurls = false),
	pages = ["Home" => "index.md",
		"Tutorials" => Any["tutorial/tutorial1.md",
				    "tutorial/tutorial2.md",
				    "tutorial/tutorial3.md",
				    "tutorial/tutorial31.md",
				    "tutorial/tutorial4.md",
				    "tutorial/tutorial41.md",
				    "tutorial/tutorial5.md"
				],
		"Function Library" => "library.md",
		]
)

deploydocs(
	repo = "github.com/rubsc/ScenTreesMakie.jl.git",
	devbranch = "main"
	#target = "build",
	#versions = ["stable" => "v^", "v#.#", "dev" => "master"]
)


text:
using Documenter, DocumenterMarkdown, LaurentPolynomials

makedocs(sitename="Laurent polynomials documentation",format=Markdown(),modules=[LaurentPolynomials])


text:
# Dependency Graph is a predictor based on Variable order Markov Chains

"""
    DG{SymbolType}(win_length::Int)

Creates a Dependency Graph predictor for `SymbolType` with lookahead window of `win_length`.
`SymbolType` can be any valid type including `Char`, `Int` etc.,

## Examples
```julia-repl
jjulia> p = DG{Char}( 4 )
DiscretePredictors.DG{Char}([*] (0)
, Char[], 4)

julia> p = DG{Int}( 3 )
DiscretePredictors.DG{Int}([*] (0)
, Int[], 3)
```

Reference:
Padmanabhan, Venkata N., and Jeffrey C. Mogul. "Using predictive prefetching to improve world wide web latency." ACM SIGCOMM Computer Communication Review 26.3 (1996): 22-36.
"""
mutable struct DG{T} <: BasePredictor{T}
    model::Trie{T,Int}
    window::Vector{T}
    win_length::Int

    DG{T}( _c::Int ) where {T} = new( Trie{T,Int}(), Vector{T}(), _c )
    # (::Type{DG{T}}){T}( _c::Int )   = new{T}( Trie{T,Int}(), Vector{T}(), _c );
end

function add!( p::DG{T}, sym::T ) where {T}
    p.model.value += 1;
    # Check if this symbol is new
    if sym ∉ keys(p.model.children)
        p.model[[sym]]  = 0;    # Start with a support of 0
    end

    for s ∈ p.window
        if haskey( p.model, [s;sym] )
            p.model[[s;sym]] += 1   # Increment occurence
            p.model[[s]] += 1       # Increment support
        else
            p.model[[s;sym]] = 1    # Create new node
            p.model[[s]] += 1        # Increment support
        end
    end

    push!( p.window, sym )   # Add this symbol to window
    if length(p.window) > p.win_length    # Trim window
        popfirst!( p.window )     # Discard oldest symbol
    end

    nothing
end

function predict( p::DG{T} ) where {T}
    # Create a dictionary with symbols
    symbols = Dict( k => 0.0 for k ∈ keys(children(p.model,Vector{T}())) )
    if isempty(symbols)
        # Just return all symbols with equal probability
        prob    = 1.0 / length( keys(symbols) )
        for s ∈ keys(symbols)
            symbols[s]  = prob
        end
    else
        support = p.model[p.window[end:end]]
        if support == 0
            prob    = 1.0 / length( keys(symbols) )
            for s ∈ keys(symbols)
                symbols[s]  = prob
            end
        end
        for (symbol,count) in children(p.model,p.window[end:end])
            symbols[symbol] = count.value / support
        end
    end
    return symbols
end

function info_string( p::DG{T} ) where {T}
    return @sprintf( "DG(%d)", p.win_length )
end

function unique_string( p::DG{T} ) where {T}
    return @sprintf( "DG_%02d", p.win_length )
end


text:
# JIT compilation of Julia code to PTX

const to = Ref(TimerOutput())

function timings!(new=TimerOutput())
  global to
  to[] = new
  return
end

timings() = (TimerOutputs.print_timer(to[]); println())

include("compiler/common.jl")
include("compiler/irgen.jl")
include("compiler/optim.jl")
include("compiler/validation.jl")
include("compiler/rtlib.jl")
include("compiler/mcgen.jl")
include("compiler/debug.jl")
include("compiler/driver.jl")

function __init_compiler__()
    # enable generation of FMA instructions to mimic behavior of nvcc
    LLVM.clopts("-nvptx-fma-level=1")

    timings!()
end


text:
import Pkg; Pkg.add("StaticArrays")
import Pkg; Pkg.add("HypergeometricFunctions")
import Pkg; Pkg.add("BenchmarkTools")
import Pkg; Pkg.add("HDF5")
import Pkg; Pkg.add("DelimitedFiles")
import Pkg; Pkg.add("Interpolations")
import Pkg; Pkg.add("ArgParse")
import Pkg; Pkg.add("SpecialFunctions")
import Pkg; Pkg.add("Distributions")
import Pkg; Pkg.add("Plots")


text:
"""
    owens_t(x, a)

Owen's T function

Ported from Matlab implementation of John Burkardt, which is distributed under the GNU LGPL license.

Reference:
    MA Porter, DJ Winstanley, Remark AS R30: A Remark on Algorithm AS76: An Integral Useful in
    Calculating Noncentral T and Bivariate Normal Probabilities, Applied Statistics, Volume 28,
    Number 1, 1979, page 113.

    JC Young, Christoph Minder, Algorithm AS 76: An Algorithm Useful in Calculating Non-Central T and
    Bivariate Normal Distributions, Applied Statistics, Volume 23, Number 3, 1974, pages 455-457.
"""
function owens_t(x::Real, a::Real)
    ng = 5
    r = [0.1477621, 0.1346334, 0.1095432, 0.0747257, 0.0333357]
    tp = 0.159155
    tv1 = 1.0E-35
    tv2 = 15.0
    tv3 = 15.0
    tv4 = 1.0E-05
    u = [0.0744372, 0.2166977, 0.3397048, 0.4325317, 0.4869533]
    # Test for X near zero.
    if abs(x) < tv1
        return tp * atan(a)
    end
    # Test for large values of abs(X).
    if tv2 < abs(x)
        return 0.
    end
    # Test for `a` near zero.
    if abs(a) < tv1
        return 0.
    end
    # Test whether abs(`a`) is so large that it must be truncated.
    xs = -0.5 * x * x
    x2 = a
    fxs = a * a
    #  Computation of truncation point by Newton iteration.
    if tv3 <= log(1.0 + fxs) - xs * fxs
        x1 = 0.5 * a
        fxs = 0.25 * fxs
        while true
            rt = fxs + 1.
            x2 = x1 + (xs * fxs + tv3 - log(rt)) / (2. * x1 * ( 1. / rt - xs))
            fxs = x2 * x2
            if abs(x2 - x1) < tv4
                break
            end
            x1 = x2
        end
    end
    # Gaussian quadrature.
    rt = 0.0
    for i in 1:ng
        r1 = 1.0 + fxs * (0.5 + u[i])^2
        r2 = 1.0 + fxs * (0.5 - u[i])^2
        rt = rt + r[i] * (exp(xs * r1) / r1 + exp(xs * r2) / r2)
    end
    return rt * x2 * tp
end


text:
include("check_manifold.jl")

@testset "PlaneParabola manifold" begin
    M = PlaneParabola()

    check_manifold(M)

    comparison = check_e2r_gradient_hessian(M)
    for (k, slope_theoretical) in [(:frgrad, 2), (:frhess, 3)]
        regressiondata = remove_small_functionvals(comparison[k])

        # Some values may remain for t close to 1
        if length(regressiondata) > 5
            slope_exp, residual = get_sloperesidual(regressiondata)
            @test isapprox(slope_exp, slope_theoretical; atol = 1e-1)
        end
    end
end

@testset "maxquadAL regularizer" begin
    x = Float64[1, 1]
    reg = maxquadAL()

    @test g(reg, x) == 0.0

    # Upper parabola point
    y, M = prox_αg(reg, [0.0, 4.0], 0.25)
    @test y[1] == 0.0
    @test 0 <= y[2] < 4.0
    @test M == Euclidean(2)

    # Lower parabola point
    y, M = prox_αg(reg, [0.0, -4.0], 0.25)
    @test y[1] == 0.0
    @test -4.0 < y[2] < 0.0
    @test M == Euclidean(2)

    # Parabola point
    y, M = prox_αg(reg, x, 0.25)
    @test y == [1.0, 1.0]
    @test M == PlaneParabola()

    y, M = prox_αg(reg, x, 0.25)
    @test y == [1.0, 1.0]
    @test M == PlaneParabola()

    @test ∇M_g(reg, M, y) == zeros(2)
    @test norm(∇²M_g_ξ(reg, M, y, rand(2))) == 0.0

    comparison = check_regularizer_gradient_hessian(M, reg)
    for (k, slope) in [(:frgrad, 2), (:frhess, 3)]
        regressiondata = remove_small_functionvals(comparison[k])

        # Some vaues may remain for t close to 1
        if length(regressiondata) > 5
            slope, residual = get_sloperesidual(regressiondata)
            @test isapprox(slope, slope; atol = 1e-1)
            @test isapprox(residual, 0; atol = 5e-2)
        end
    end
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Compute the satellite position.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export satellite_position_i

"""
    satellite_position_i(a::Number, e::Number, i::Number, RAAN::Number, w::Number, f::Number)

Compute the satellite position in the Earth-Centered Inertial (ECI) reference
frame given the orbital elements `a`, `e`, `i`, `RAAN`, `w`, and `f`.

Notice that the ECI frame used will be the same as the frame of the orbital
elements.

# Args

* `a`: Semi-major axis.
* `e`: Eccentricity.
* `i`: Inclination [rad].
* `RAAN`: Right ascension of the ascending node [rad].
* `w`: Argument of perigee [rad].
* `f`: True anomaly [rad].

# Returns

* The satellite position vector represented in the ECI reference frame.
* The unit vector perpendicular to the satellite position vector that lies on
  the orbit plane represented in the ECI reference frame.

# Remarks

The satellite position vector will have the same unit of the semi-major axis.

"""
function satellite_position_i(a::Number, e::Number, i::Number, RAAN::Number,
                              w::Number, f::Number)
    # Compute the radius from the focus.
    norm_r = a*(1-e^2)/(1+e*cos(f))

    # Let s be the coordinate system in which:
    #     - The X axis points towards the satellite;
    #     - The Z axis is normal to the orbit plane (right-hand direction);
    #     - The Y axis completes a right-hand coordinate system.
    #
    # Thus, the satellite vector represented in the s coordinate frame is:
    r_s = SVector{3}(1,0,0)*norm_r

    # rt is the versor perpendicular to the r vector that lies on the orbit
    # plane.
    rt_s = SVector{3}(0,1,0)

    # Compute the matrix that rotates from the S coordinate frame to the
    # Inertial coordinate Frame.
    Dis = angle_to_dcm(RAAN, i, w+f, :ZXZ)'

    # Compute the satellite vector represented in the Inertial coordinate
    # frame.
    r_i = Dis*r_s

    # Compute unit vector `rt` represented in the Inertial coordinate frame.
    rt_i = Dis*rt_s

    # Return.
    return r_i, rt_i
end


text:
function dijkstra_test()
    (A,xy,labels) = load_matrix_network_metadata("airports")
    A = -A; # fix funny encoding of airport data
    lax = 247; rst = 355
    (d,pred) = dijkstra(A,lax)

    if maximum(d) != 540
        error("dijkstra failed")
    end

    if pred[end] != lax
        error("dijkstra failed")
    end
    return true
end


text:
using Clang
using Clang: LibClang

const dir = @__DIR__

"get file and line number"
function expand_location(cursor)
    loc = location(cursor)
    line = Cint[0]
    column = Cint[0]
    offset = Cint[0]
    LibClang.clang_getExpansionLocation(loc, C_NULL, line, column, offset)
    return filename(cursor), line[] #, column[], offset[]
end

function extract_comment(fn, line)
    lines = readlines(fn)
    commentlines = []
    if line - 1 < 1
        return ""
    end
    if strip(lines[line-1]) == "*/"
        for idx in line-2:-1:1
            !startswith(strip(lines[idx]), "*") && break
            push!(commentlines, replace(lines[idx], r"^\s*\* ? ?" => ""))
        end
        reverse!(commentlines)
    end
    return join(commentlines, "\n")
end

function doxygen2julia(str)
    str = replace(str, r"\\c ([^ ]+)" => s"`\1`")
    str = replace(str, "\\param " => "param: ")
    str = replace(str, "\\brief " => "")
    str = replace(str, "\\return " => "Return ")
    str = replace(str, r"\\sa (.*)" => s" - `\1`")
    str = replace(str, "\\code" => "```C")
    str = replace(str, "\\endcode" => "```")
    str = replace(str, "\\note " => "!!!note\n\n    ")
    str = replace(str, "\\warning " => "!!!warning\n\n    ")
    str = replace(str, "\\li " => " - ")
    str = replace(str, "\\\\" => Char(0))
    str = replace(str, r"\\([^\\])" => s":\1")
    str = replace(str, "\0" => "\\\\")
    str
end

function docstring(cursor)
    filename(cursor) == "" && return nothing
    fn, line = expand_location(cursor)
    body = doxygen2julia(extract_comment(fn, line))
    body == "" && return nothing
    # TODO: get the real signature of the thing, if it's a function.
    # rewrite can get that, probably.
    sig = replace(repr(cursor), r"CLCursor [^ ]* " => "")
    return """
        $sig

    $body
    """
end

function main()
    sdl2_version = "2.0.12"

    LIBSDL_INCLUDE = "SDL2-$sdl2_version/include/"

    # readdir() sorts the names, and Clang seems to generate the bindinds
    # sequentially in this order, so that should make our diffs easy to read
    # between versions
    LIBSDL_HEADERS = [joinpath(LIBSDL_INCLUDE, header) for header in readdir(LIBSDL_INCLUDE) if endswith(header, ".h")]

    # I think Jonathan excluded these, too.
    filter!(!fn -> occursin("test", fn), LIBSDL_HEADERS)
    filter!(!fn -> occursin("opengl", fn), LIBSDL_HEADERS)
    filter!(!fn -> occursin("config", fn), LIBSDL_HEADERS)
    filter!(!fn -> occursin("SDL_log", fn), LIBSDL_HEADERS)

    shitlist = [
        "SDL_FILE",
        "SDL_LINE",
        "BYTECAP",
        "OUT_Z_CAP",
        "OUT_CAP",
        "BYTEORDER",
       ]

    docstrings = Dict()

    # create a work context
    ctx = DefaultContext()

    # parse headers
    parse_headers!(ctx, LIBSDL_HEADERS,
                   args=["-I", joinpath(LIBSDL_INCLUDE, ".."), map(x->"-I"*x, find_std_headers())...],
                   includes=vcat(LIBSDL_INCLUDE, CLANG_INCLUDE),
                   )

    # settings
    ctx.libname = "libsdl2"
    ctx.options["is_function_strictly_typed"] = false
    ctx.options["is_struct_mutable"] = true

    # write output
    api_file = joinpath(@__DIR__, "new_bindings", "libsdl2_api.jl")
    api_stream = open(api_file, "w")

    trans_unit = first(ctx.trans_units)

    for trans_unit in ctx.trans_units
        root_cursor = getcursor(trans_unit)
        push!(ctx.cursor_stack, root_cursor)
        header = spelling(root_cursor)
        @info "wrapping header: $header ..."
        # loop over all of the child cursors and wrap them, if appropriate.
        ctx.children = children(root_cursor)
        for (i, child) in enumerate(ctx.children)
            child_name = name(child)
            child_header = filename(child)
            ctx.children_index = i
            # choose which cursor to wrap
            startswith(child_name, "__") && continue  # skip compiler definitions
            child_name in keys(ctx.common_buffer) && continue  # already wrapped
            child_header != header && continue  # skip if cursor filename is not in the headers to be wrapped

            if any(occursin.(shitlist, child_name))
                @info "Ignoring $child_name"
                continue
            end

            docname = replace(child_name, r"\(.*\)" => "")
            if !isempty(docname)
                doc = docstring(child)
                if !isnothing(doc)
                    docstrings[docname] = doc
                end
            end

            wrap!(ctx, child)
        end
        @info "writing $(api_file)"
        println(api_stream, "# Julia wrapper for header: $(basename(header))")
        println(api_stream, "# Automatically generated using Clang.jl\n")
        #= ctx.api_buffer = rewrite(ctx.api_buffer) =#
        print_buffer(api_stream, ctx.api_buffer)
        empty!(ctx.api_buffer)  # clean up api_buffer for the next header
    end
    close(api_stream)

    # write "common" definitions: types, typealiases, etc.
    common_file = "$dir/new_bindings/libsdl2_types.jl"
    open(common_file, "w") do f
        println(f, "# Automatically generated using Clang.jl\n")
        print_buffer(f, dump_to_buffer(ctx.common_buffer))
    end

    docsfile = "$dir/new_bindings/libsdl2_docs.jl"
    open(docsfile, "w") do f
        println(f, "# Automatically generated using Clang.jl\n")
        for (name, doc) in docstrings
            # This is pretty ugly, but it gets the job done.
            println(f, """if isdefined(@__MODULE__, Symbol("$name"))
                        @eval @doc \"""
                    $doc
                    \$(
                        if $name === Nothing
                            "Treat as a void pointer\\n"
                        else
                            replace(string(@doc $name), r"^No documentation found.\\n" => "")
                        end
                    )
                    \""" $name
                    end
                    """)
        end
    end

    # remove SDL_ and copy over.
    function rm_SDL_(fn)
        str = String(read(fn))
        str = replace(str, r"(::|[^:])SDL_" => s"\1");
        write(fn, str)
    end
    rm_SDL_(docsfile)
    rm_SDL_(api_file)
    rm_SDL_(common_file)

    return nothing
end

# uncomment the following code to generate dependency and template files
# copydeps(dirname(api_file))
# print_template(joinpath(dirname(api_file), "LibTemplate.jl"))

#=
# Seems to work.
for c in ctx.children
    filename(c) == "" && continue
    fn, line = expand_location(c)
    body = doxygen2julia(extract_comment(fn, line))
    body == "" && continue
    # TODO: get the real signature of the thing, if it's a function.
    # rewrite can get that, probably.
    sig = replace(repr(c), r"CLCursor [^ ]* " => "")
    println("""
        $sig

    $body
    """)
end
=#

#= function rewrite!(e::Expr) =#
#=     # Add deprecated warning to some functions =#
#=     #   ref: [Remove unused CompilationDatabase::MappedSources](https://reviews.llvm.org/D32351) =#
#=     fname = e.args[1].args[1] =#
#=     deprecated_func = [ =#
#=         # :clang_CompileCommand_getNumMappedSources, # not export =#
#=         :clang_CompileCommand_getMappedSourcePath, =#
#=         :clang_CompileCommand_getMappedSourceContent, =#
#=     ] =#

#=     if e.head == :function && fname in deprecated_func =#
#=         msg = """ =#
#=         `$fname` Left here for backward compatibility. =#
#=         No mapped sources exists in the C++ backend anymore. =#
#=         This function just return Null `CXString`. =#
#=         See: =#
#=         - [Remove unused CompilationDatabase::MappedSources](https://reviews.llvm.org/D32351) =#
#=         """ =#
#=         insert!(e.args[2].args, 1, Expr(:macrocall, Symbol("@warn"), :Base, msg)) =#
#=     end =#
#=     e =#
#= end =#
#= rewrite!(x) = x =#
#= rewrite(v::Vector) = map(rewrite!, v) =#


text:
using CUDA
using GemmKernels

M = parse(Int, ARGS[1])
N = parse(Int, ARGS[2])
K = parse(Int, ARGS[3])

transpose_a = (ARGS[4] == "n" ? false : ARGS[4] == "t" ? true : error("Invalid memory layout for A: $(ARGS[4])"))
transpose_b = (ARGS[5] == "n" ? false : ARGS[5] == "t" ? true : error("Invalid memory layout for B: $(ARGS[5])"))

function benchmark_matmul(a, b, c)
    alpha = rand(Float32)
    beta = rand(Float32)

    CUDA.@sync begin
        GemmKernels.BLAS.gemmEx!(!transpose_a ? 'N' : 'T', !transpose_b ? 'N' : 'T', alpha, a, b, beta, c)
    end
end

a_h = rand(Float16, (M, K)) / sqrt(Float16(K))
b_h = rand(Float16, (K, N)) / sqrt(Float16(K))
c_h = rand(Float32, (M, N))

# Transpose input if necessary
a_h = transpose_a ? transpose(a_h) : a_h
b_h = transpose_b ? transpose(b_h) : b_h

a   = CuArray(a_h)
b   = CuArray(b_h)
c   = CuArray(c_h)

# warmup
benchmark_matmul(a, b, c)

# profile
for i = 1 : 10
    CUDA.@profile benchmark_matmul(a, b, c)
end


text:
# TODO: Make more tests with Settings
@testset "Operation Model kwargs" begin
    template = get_thermal_dispatch_template_network()
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys5_re = PSB.build_system(PSITestSystems, "c_sys5_re")

    @test_throws MethodError OperationsProblem(template, c_sys5; bad_kwarg = 10)

    op_problem = OperationsProblem(
        template,
        c_sys5;
        use_forecast_data = false,
        optimizer = GLPK_optimizer,
    )
    @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT
    @test PSI.get_use_forecast_data(
        PSI.get_settings(PSI.get_optimization_container(op_problem)),
    ) == false

    op_problem = OperationsProblem(
        MockOperationProblem,
        get_thermal_dispatch_template_network(),
        c_sys5_re;
        optimizer = GLPK_optimizer,
        balance_slack_variables = true,
    )
    op_problem = OperationsProblem(
        get_thermal_dispatch_template_network(),
        c_sys5;
        use_forecast_data = false,
        optimizer = GLPK_optimizer,
    )
    @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT
    @test PSI.get_use_forecast_data(
        PSI.get_settings(PSI.get_optimization_container(op_problem)),
    ) == false

    #"Test passing custom JuMP model"
    my_model = JuMP.Model()
    my_model.ext[:PSI_Testing] = 1
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    op_problem = OperationsProblem(
        get_thermal_dispatch_template_network(),
        c_sys5,
        my_model;
        optimizer = GLPK_optimizer,
        use_parameters = true,
    )
    build!(op_problem; output_dir = mktempdir(cleanup = true)) == PSI.BuildStatus.BUILT
    @test haskey(PSI.get_optimization_container(op_problem).JuMPmodel.ext, :PSI_Testing)
    @test (
        :ParameterJuMP in keys(PSI.get_optimization_container(op_problem).JuMPmodel.ext)
    ) == true
end

@testset "Set optimizer at solve call" begin
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    template = get_thermal_standard_uc_template()
    set_service_model!(template, ServiceModel(VariableReserve{ReserveUp}, RangeReserve))
    UC = OperationsProblem(template, c_sys5)
    output_dir = mktempdir(cleanup = true)
    @test build!(UC; output_dir = output_dir) == PSI.BuildStatus.BUILT
    @test solve!(UC; optimizer = GLPK_optimizer) == RunStatus.SUCCESSFUL
    res = ProblemResults(UC)
    @test isapprox(get_objective_value(res), 340000.0; atol = 100000.0)
    vars = res.variable_values
    @test :P__ThermalStandard in keys(vars)
    export_results(res)
    results_dir = joinpath(output_dir, "results")
    @test isfile(joinpath(results_dir, "optimizer_stats.csv"))
    variables_dir = joinpath(results_dir, "variables")
    @test isfile(joinpath(variables_dir, "P__ThermalStandard.csv"))
end

@testset "Test optimization debugging functions" begin
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    template = get_thermal_standard_uc_template()
    set_service_model!(template, ServiceModel(VariableReserve{ReserveUp}, RangeReserve))
    op_problem = OperationsProblem(template, c_sys5; optimizer = GLPK_optimizer)
    @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT
    optimization_container = PSI.get_optimization_container(op_problem)
    MOIU.attach_optimizer(optimization_container.JuMPmodel)
    constraint_indices = get_all_constraint_index(op_problem)
    for (key, index, moi_index) in constraint_indices
        val1 = get_con_index(op_problem, moi_index)
        val2 = optimization_container.constraints[key].data[index]
        @test val1 == val2
    end
    @test isnothing(get_con_index(op_problem, length(constraint_indices) + 1))

    var_indices = get_all_var_index(op_problem)
    for (key, index, moi_index) in var_indices
        val1 = get_var_index(op_problem, moi_index)
        val2 = optimization_container.variables[key].data[index]
        @test val1 == val2
    end
    @test isnothing(get_var_index(op_problem, length(var_indices) + 1))
end

# @testset "Test print methods" begin
#     template = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
#     c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
#     op_problem = OperationsProblem(
#         MockOperationProblem,
#         template,
#         c_sys5;
#         optimizer = GLPK_optimizer,
#         use_parameters = true,
#     )
#     list = [template, op_problem, op_problem.optimization_container, services]
#     _test_plain_print_methods(list)
#     list = [services]
#     _test_html_print_methods(list)
# end

@testset "Operation Model Solve with Slacks" begin
    c_sys5_re = PSB.build_system(PSITestSystems, "c_sys5_re")
    networks = [StandardPTDFModel, DCPPowerModel, ACPPowerModel]
    for network in networks
        template = get_thermal_dispatch_template_network(network)
        op_problem = OperationsProblem(
            template,
            c_sys5_re;
            balance_slack_variables = true,
            optimizer = ipopt_optimizer,
            PTDF = PTDF(c_sys5_re),
        )
        @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        @test solve!(op_problem) == RunStatus.SUCCESSFUL
    end
end

@testset "Default Operations Constructors" begin
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    op_problem_ed = EconomicDispatchProblem(c_sys5; output_dir = mktempdir())
    moi_tests(op_problem_ed, false, 120, 0, 120, 120, 24, false)
    op_problem_uc = UnitCommitmentProblem(c_sys5; output_dir = mktempdir())
    moi_tests(op_problem_uc, false, 480, 0, 240, 120, 144, true)
    ED_output = run_economic_dispatch(
        c_sys5;
        output_dir = mktempdir(),
        optimizer = fast_lp_optimizer,
    )
    UC_output =
        run_unit_commitment(c_sys5; output_dir = mktempdir(), optimizer = fast_lp_optimizer)
    @test ED_output == RunStatus.SUCCESSFUL
    @test UC_output == RunStatus.SUCCESSFUL
end

@testset "Test Locational Marginal Prices between DC lossless with PowerModels vs StandardPTDFModel" begin
    networks = [DCPPowerModel, StandardPTDFModel]
    sys = PSB.build_system(PSITestSystems, "c_sys5")
    parameters = [true, false]
    ptdf = PTDF(sys)
    # These are the duals of interest for the test
    dual_constraint =
        [[:nodal_balance_active__Bus], [:CopperPlateBalance, :network_flow__Line]]
    LMPs = []
    for (ix, network) in enumerate(networks), p in parameters
        template = get_template_dispatch_with_network(network)
        op_problem = OperationsProblem(
            template,
            sys;
            optimizer = OSQP_optimizer,
            use_parameters = p,
            PTDF = ptdf,
            constraint_duals = dual_constraint[ix],
        )
        @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        @test solve!(op_problem) == RunStatus.SUCCESSFUL
        res = ProblemResults(op_problem)

        # These tests require results to be working
        if network == StandardPTDFModel
            push!(LMPs, abs.(psi_ptdf_lmps(res, ptdf)))
        else
            duals = res.dual_values[:nodal_balance_active__Bus]
            duals = abs.(duals[:, propertynames(duals) .!== :DateTime])
            push!(LMPs, duals[!, sort(propertynames(duals))])
        end
    end
    # TODO: the above calculation executes, but the following test does not pass
    #@test isapprox(convert(Array, LMPs[1]), convert(Array, LMPs[2]), atol = 100.0)
end

@testset "Test ProblemResults interfaces" begin
    sys = PSB.build_system(PSITestSystems, "c_sys5_re")
    template = get_template_dispatch_with_network(CopperPlatePowerModel)
    op_problem = OperationsProblem(
        template,
        sys;
        optimizer = OSQP_optimizer,
        use_parameters = true,
        constraint_duals = [:CopperPlateBalance],
    )
    @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT
    @test solve!(op_problem) == RunStatus.SUCCESSFUL

    optimization_container = PSI.get_optimization_container(op_problem)
    constraints = PSI.get_constraints(optimization_container)[:CopperPlateBalance]
    dual_results = read_duals(optimization_container, [:CopperPlateBalance])
    for i in axes(constraints)[1]
        dual = JuMP.dual(constraints[i])
        @test isapprox(dual, dual_results[:CopperPlateBalance][i, 1])
    end

    system = PSI.get_system(op_problem)
    params = PSI.get_parameters(optimization_container)[:P__max_active_power__PowerLoad]
    param_vals = PSI.axis_array_to_dataframe(params.parameter_array)
    param_mult = PSI.axis_array_to_dataframe(params.multiplier_array)
    for load in get_components(PowerLoad, system)
        name = get_name(load)
        vals = get_time_series_values(Deterministic, load, "max_active_power")
        @test all([-1 * x == get_max_active_power(load) for x in param_mult[!, name]])
        @test all(vals .== param_vals[!, name])
    end

    res = ProblemResults(op_problem)
    @test length(get_existing_variables(res)) == 1
    @test length(get_existing_parameters(res)) == 1
    @test length(get_existing_duals(res)) == 1
    @test get_model_base_power(res) == 100.0
    @test isa(get_objective_value(res), Float64)
    @test isa(get_variables(res), Dict{Symbol, DataFrames.DataFrame})
    @test isa(get_total_cost(res), Float64)
    @test isa(get_optimizer_stats(res), PSI.OptimizerStats)
    @test isa(get_duals(res), Dict{Symbol, DataFrames.DataFrame})
    @test isa(get_parameters(res), Dict{Symbol, DataFrames.DataFrame})
    @test isa(get_resolution(res), Dates.TimePeriod)
    @test isa(get_system(res), PSY.System)
    @test length(get_timestamps(res)) == 24
end

@testset "Test Serialization, deserialization and write optimizer problem" begin
    path = mktempdir(cleanup = true)
    sys = PSB.build_system(PSITestSystems, "c_sys5_re")
    template = get_template_dispatch_with_network(CopperPlatePowerModel)
    op_problem = OperationsProblem(
        template,
        sys;
        optimizer = OSQP_optimizer,
        use_parameters = true,
        constraint_duals = [:CopperPlateBalance],
    )
    @test build!(op_problem; output_dir = path) == PSI.BuildStatus.BUILT
    @test solve!(op_problem) == RunStatus.SUCCESSFUL

    file_list = sort!(collect(readdir(path)))
    @test "OptimizationModel.json" in file_list
    @test "OperationProblem.bin" in file_list
    filename = joinpath(path, "OperationProblem.bin")
    ED2 = OperationsProblem(filename, optimizer = OSQP_optimizer)
    build!(ED2, output_dir = path)
    psi_checksolve_test(ED2, [MOI.OPTIMAL], 240000.0, 10000)

    path2 = mktempdir(cleanup = true)
    op_problem_no_sys = OperationsProblem(
        template,
        sys;
        optimizer = OSQP_optimizer,
        use_parameters = true,
        system_to_file = false,
        constraint_duals = [:CopperPlateBalance],
    )

    @test build!(op_problem_no_sys; output_dir = path2) == PSI.BuildStatus.BUILT
    @test solve!(op_problem) == RunStatus.SUCCESSFUL

    file_list = sort!(collect(readdir(path2)))
    @test .!all(occursin.(r".h5", file_list))
    filename = joinpath(path2, "OperationProblem.bin")
    ED3 = OperationsProblem(filename; system = sys, optimizer = OSQP_optimizer)
    build!(ED3, output_dir = path2)
    psi_checksolve_test(ED3, [MOI.OPTIMAL], 240000.0, 10000)
end


text:
"""This script is intended to do the post processing part of our speciation study
The end goal is to be able to split the .fa files according to timestep using the inds file and
the seqs file. It automatically takes the first half of the genome and concatenates it into a
"chromosome" for further analysis.
These are then loaded into a further tool (eg. FasTree)

this is meant to be parallelised via an sbatch command in slurm

No warranties for anything. Use at your own risk.
Robin Rölz, 21/09/2021"""

results_path = "gemm/results/"
n_genes = 22

"""
    searchdir(path, key)

short helper function to find files in a directory matching specifications
"""

function searchdir(path, key)
    item = filter(x->occursin(key,x), readdir(path))
    return item[1]
end

"""
    concat_genome(genome)

takes the genome of one individual in .fa format and takes the first half of
genes and then concatenates those into a single .fa entry such that each gene
is represented once
"""
function concat_chrs(genome)
    chrs_len = Int(length(genome)/2) #get number of genes in chromosome

    label = split(genome[1], "_") #split first line to get label
    label = join(label[1:2], "_") #label consists of species and id
    chrs = label*"\n" #the chromosome that is going to be printed

    i = 1
    for l in genome[1:chrs_len]
        if i%2 == 0
            chrs = chrs*l
        end
        i+=1
    end

    return chrs*"\n"
end

"""
    count_inds(folderpath)

counts the individuals in an ind file
returns a tuple with the timestep and the number of individuals
"""

function count_inds(folderpath)
    #individuals in each timestep
    #count the amount of individuals in each timestep
    n_inds = []
    i=0
    first = true
    t_index = String
    t_all = []
    for line in eachline(joinpath(folderpath, searchdir(folderpath, r"ind.*\.tsv")))
        if first == true
            first = false
            continue
        end
        t = split(line, "\t")[1] #this gets me the time
        if t != t_index #is this a new time index?
            t_index = t #set new indexed timestep
            push!(t_all, t_index) #collect all t
            i+=1 # increment counter
            push!(n_inds, 0) #make a new space at the end of the n_vector
        end
        n_inds[i] +=1
    end
    return hcat(t_all, n_inds)
end

"""
    get_results(results_path)

grab all results folders and returns one according to the SLURM_ARRAY_TASK_ID
"""

function get_results(results_path)
    folders = readdir(results_path)
    task_id = Base.parse(Int, ENV["SLURM_ARRAY_TASK_ID"])
    taskfolder = folders[task_id]
    return joinpath(results_path, taskfolder)
end

"""
    split_fa(folder)

split the seqs.fa file in the folder according to the amount of individuals
in single fasta files to allow for post processing.
"""

function split_fa(folderpath, n_genes)

    fa_length = n_genes*2 #number of lines the genome takes in the fasta file

    indtable = count_inds(folderpath)

    open(joinpath(folderpath, searchdir(folderpath, r"seqs.*\.fa"))) do seqs
        #iterating over the timesteps
        chrsdir = joinpath(folderpath, "chrs")
        mkpath(chrsdir)
        for n in 1:size(indtable, 1)
            outseqs = ""
            #this is all the genes from one time step
            for m in 1:indtable[n,2]
                genome = []
                for o in 1:fa_length
                    push!(genome, readline(seqs))
                end
                outseqs = outseqs*concat_chrs(genome)
            end

            filename_seqs = "chrs_"*string(indtable[n,1])*".fa" #labels files after the timestep
            filename_seqs = joinpath(chrsdir, filename_seqs)
            open(filename_seqs, "w") do io
                write(io, outseqs)
            end
        end
    end
end

println("Starting Analysis")
folderpath = get_results(results_path)
println("Splitting .fa in "*folderpath)
split_fa(folderpath, n_genes)
println("Finished splitting .fa in "*folderpath)

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
function [ varp ] = rayleighTop(varz,varztop,varp,varb,raydmpcoef,raydmpz,trigpi,i,k)
%UNTITLED2 Summary of this function goes here
%   Set varb to 0 for w,th,and pic

  % Compute Local Damping Coefficient
coef = raydmpcoef*0.5*(1-cos(trigpi*(varz-raydmpz)/(varztop-raydmpz)));

% Apply Sponge to Slowly Remove Perturbations
varp = varp - coef*(varp-varb);

end

max_stars_repo_name: jeffreycunn/toy-atmospheric-model

max_issues_repo_name: jeffreycunn/toy-atmospheric-model

max_forks_repo_name: jeffreycunn/toy-atmospheric-model


text:
# Find the Most Frequent Words in a String

function most_frequent_k_mers(text, k)
    counts = Dict{AbstractString, Int64}()
    for ix in 1:(length(text)-k+1)
        k_mer = text[ix:ix+k-1]
        count = get(counts, k_mer, 0)
        counts[k_mer] = count + 1
    end
    maxcount = maximum(values(counts))
    [k_mer for (k_mer, count) in counts if count == maxcount]
end

function main()
    text = readline()
    k = parse(Int64, readline())
    result = most_frequent_k_mers(text, k)
    println(join(result, " "))
end

main()


text:
testobject(val) = Dict("tree" => Dict(
    "val" => val
))

@info "deleteenvkey!"
@test !haskey(ENV, "FOOTEST")
@test_nowarn Configs.deleteenvkey!("FOOTEST")
ENV["FOOTEST"] = "foo"
@test_nowarn Configs.deleteenvkey!("FOOTEST")
@test !haskey(ENV, "FOOTEST")

@info "parseenvkey"
@test Configs.parseenvkey("FOOTEST", "nope") === "nope"
ENV["FOOTEST"] = "yep"
@test Configs.parseenvkey("FOOTEST", "nope") === "yep"
Configs.deleteenvkey!("FOOTEST")

@info "parsecustomenv!"
foo = testobject("FOOTEST")
@test_nowarn Configs.parsecustomenv!(foo)
@test isequal(foo, Dict("tree" => Dict()))
foo = testobject("FOOTEST")
ENV["FOOTEST"] = "yep"
@test_nowarn Configs.parsecustomenv!(foo)
@test isequal(foo, testobject("yep"))
Configs.deleteenvkey!("FOOTEST")

@info "makeimmutable"
isimmutable = @test_nowarn Configs.makeimmutable(testobject(["one", Dict("one" => 1)]))
@test isequal(isimmutable, (tree = (val = ("one", (one = 1,)),),))

@info "override!"
base = Dict()
add = testobject("test")
@test_nowarn Configs.override!(base, add)
@test isequal(base, add)
add = Dict("foo" => ["one", "two"], "bar" => 2, "tree" => Dict("added" => 2))
@test_nowarn Configs.override!(base, add)
@test isequal(base, Dict{Any,Any}("bar" => 2,"tree" => Dict{Any,Any}("val" => "test","added" => 2),"foo" => ["one", "two"]))

@info "getconffiles"
@test_throws Configs.Configserror Configs.getconffiles("foobar")
@test_nowarn Configs.getconffiles("configs")
@test_nowarn Configs.getconffiles(joinpath(pwd(), "configs"))

@info "readconffile"
@test_nowarn Configs.readconffile(joinpath(pwd(), "configs"), "default.yml")
file = Configs.readconffile(joinpath(pwd(), "configs"), "default.yml")
@test file isa Dict
@test file["database"]["connection"]["port"] === 3600

@info "parseconfigs"
@test_nowarn Configs.parseconfigs("dev", "configs")
@test_throws Configs.Configserror Configs.parseconfigs("dev", "badconfigs")
conf = Configs.parseconfigs("DEPLOYMENT", "configs")
@test isequal(conf.files, ["default.yml", "staging.jl", "custom-environment-variables.json"])

@info "pathtodict"
@test_nowarn Configs.configpathtodict("one,two.three", "test")
dict = Configs.configpathtodict("one.two.three", "test")
@test dict isa Dict{String, Any}
@test dict["one"]["two"]["three"] === "test"

@info "parseconfigpath"
testtuple = (; one=(; two = "test"))
@test_nowarn Configs.parseconfigpath(testtuple, "one.two")
@test_throws Configs.Configserror Configs.parseconfigpath(testtuple, "one.two.foo")
@test Configs.parseconfigpath(testtuple, "one.two") === "test"

max_stars_repo_name: citkane/Config

max_issues_repo_name: citkane/Config

max_forks_repo_name: citkane/Config


text:
#include("common.jl")
#include("algorithm_parameters.jl")
#include("dynamical_systems.jl")

"""
`TMState` is a named tuple of the state variables of a temporal memory.

- `α`: active neurons
- `Π`: predictive neurons
- `WN`: winning neurons
- `Πₛ`: predictive dendritic segments (to calculate winning segments)
- `Mₛ`: matching dendritic segments (didn't receive enough input to activate, but enough to learn)
- `ovp_Mₛ`:
"""
mutable struct TMState
  state::NamedTuple{
    (:α, :Π, :WN, :Πₛ, :Mₛ, :ovp_Mₛ),  # Ncell, Ncell, Ncell, Nseg, Nseg, Nseg
    Tuple{CellActivity, CellActivity, CellActivity,
          BitArray{1}, BitArray{1}, Vector{Int}}
  }
  init_Nseg::Int
end
TMState(Nₙ::Int, Nseg)= TMState((
  α=falses(Nₙ), Π=falses(Nₙ), WN=falses(Nₙ),
  Πₛ=falses(Nseg), Mₛ=falses(Nseg),
  ovp_Mₛ=zeros(Nseg)
), Nseg)
Base.getproperty(s::TMState, name::Symbol)= name === :state ?
    getfield(s,:state) :
    name === :init_Nseg ? getfield(s, :init_Nseg) :
      getproperty(getfield(s,:state),name)
update_TMState!(s::TMState; Nseg,α,Π,WN,Πₛ,Mₛ,ovp_Mₛ)=
    s.state= (α=α, Π= Π, WN= WN, Πₛ= padfalse(Πₛ,Nseg),
              Mₛ= padfalse(Mₛ,Nseg), ovp_Mₛ= padfalse(ovp_Mₛ,Nseg))
reset!(s::TMState)= begin
  Nₙ= length(s.α);
  update_TMState!(s, Nseg=s.init_Nseg, α=falses(Nₙ), Π=falses(Nₙ), WN=falses(Nₙ),
    Πₛ=falses(s.init_Nseg), Mₛ=falses(s.init_Nseg), ovp_Mₛ=zeros(s.init_Nseg))
end

"""
    TemporalMemory(params::TMParams; recurrent= true, distal_input_size=0)

`TemporalMemory` learns to predict sequences of input Sparse Distributed Representations (SDRs),
usually generated by a [`SpatialPooler`](@ref).
It learns to represent each input symbol in the temporal context of the symbols that come before it in the sequence,
using the individual neurons of each minicolumn and their distal synapses.

When considering a neuron layer with proximal and distal synapses, the spatial pooler is a way to
activate and learn the proximal synapses, while the temporal memory is a way to activate and learn
the distal synapses.

Parameters:

- `params`: [`TMParams`](@ref) defining the layer size and algorithm constants
- `recurrent`: presence of recurrent distal connections between the layer's neurons
- `distal_input_size`: the number of presynaptic neurons for each layer with incoming distal connections (apart from same layer's neurons)

## High-order predictions and Ambiguity

The [neurons-thousand-synapses paper](https://www.frontiersin.org/articles/10.3389/fncir.2016.00023/full) describes the
Temporal Memory's properties, and especially the ability to
- make "high-order" predictions, based on previous inputs potentially going far back in time
- represent ambiguity in the form of simultaneous predictions

For more information see figures 2,3 in the paper.

# TM activation

Overview of the temporal memory's process:

1. Activate neurons (fire, proximal synapses)
2. Predict neurons (depolarize, distal/apical synapses)
3. Learn distal/apical synapses:
    - adapt existing synapses
    - create new synapses/dendrites

## Activation

An SDR input activates some minicolumns of the neuron layer.
If some neurons in the minicolum were predicted at the previous step, they activate faster than the rest and inhibit them.
If no neuron was predicted, all the neurons fire (minicolumn bursting).

TODO

---
See also: [`TMParams`](@ref) for parameter and symbol description, [`DistalSynapses`](@ref), [`TMState`](@ref)
"""
struct TemporalMemory
  params::TMParams
  distalSynapses::DistalSynapses
  previous::TMState
  recurrent::Bool
end
function TemporalMemory(params::TMParams; recurrent= true, distal_input_size=0)
  @unpack Nₙ, Nc, k = params
  distal_input_size= Int(distal_input_size)
  @assert recurrent || distal_input_size > 0 "At least 1 of recurrent or external input connections are needed"

  # NOTE: it's not possible to have independent distal synapses for each presynaptic input,
  # because we need a single address space for postsynaptic segments. Otherwise,
  # it would be impossible to mix signals from multiple sources.
  N_presynaptic= recurrent ? Nₙ + distal_input_size : distal_input_size
  Nseg_init= 0

  TemporalMemory(params,
      DistalSynapses(N_presynaptic, Nc, k;
        Nseg_init= Nseg_init, params= DistalSynapseParams(params)),
      TMState(Nₙ,Nseg_init),
      recurrent)
end

Nₙ(tm::TemporalMemory)= tm.params.Nₙ
reset!(tm::TemporalMemory)= begin
  reset!(tm.distalSynapses)
  reset!(tm.previous)
end

"""
`step!(tm::TemporalMemory, c, distal_input=[])` evolves the Temporal Memory to the next timestep
given the minicolumn activations and any distal input.
Returns the state of the region's neurons: active, predictive, bursting (minicolumns)

See also: [`tm_activate`](@ref), [`tm_predict`](@ref)
"""
# Given a column activation pattern `c` (SP output), step the TM
function step!(tm::TemporalMemory, c, distal_input=falses(0))
  s= tm.distalSynapses; p= tm.previous

  α, B, WN= tm_activate(tm, c, p.Π)
  # If recurrent, concatenate distal input activity with this layer's activity
  distal_input= tm.recurrent ? [α;distal_input] : distal_input
  Π, Πₛ, Mₛ, ovp_Mₛ= tm_predict(tm, distal_input)

  # Learn
  WS, WS_burst= calculate_WS!(s, p.Πₛ,p.ovp_Mₛ,α,B)
  # Update winner neurons with entries from bursting columns
  WN[NS(s)*WS_burst .> 0].= true
  step!(s, p.WN,WS, α, p.α,p.Mₛ,p.ovp_Mₛ)
  update_TMState!(p, Nseg=Nₛ(s),
                  α=α, Π=Π, WN=WN, Πₛ=Πₛ, Mₛ=Mₛ, ovp_Mₛ=ovp_Mₛ)
  return (
    active= α,
    predictive= Π,
    bursting= B
  )
end

"""
Active and predictive cells given the minicolumn activations and any distal input.

See also: [`tm_activate`](@ref), [`tm_predict`](@ref)
"""
(tm::TemporalMemory)(c, distal_input=falses(0))= begin
  α= tm_activate(tm,c,tm.previous.Π)[1]
  # If recurrent, concatenate distal input activity with this layer's activity
  distal_input= tm.recurrent ? [α;distal_input] : distal_input
  (
    active= α,
    predictive= tm_predict(tm,distal_input)[1]
  )
end

"""
`tm_activate(tm::TemporalMemory, c, Π)` calculates

1. which minicolumns burst,
2. which neurons in the layer are activated
3. which become predictive

for minicolumn activation `c` (size `Nc`) given by the [`SpatialPooler`](@ref)
and the previously predictive neurons `Π` (size `Nₙ`).

# Returns

1. `a`: neuron activation (`Nₙ`)
2. `B`: bursting minicolumns (`Nc`)
3. `WN`: "winning" neurons (`Nₙ`)

See also: [`tm_predict`](@ref), [`DistalSynapses`](@ref)
"""
function tm_activate(tm::TemporalMemory, c, Π)
  @unpack Nc, k = tm.params

  # bursting minicolumns (Nc)
  burst(c,Π)= c .& .!@percolumn(any, Π, k)
  # minicolumns with a predictive neuron (k × Nc)
  predicted(c,Π)= @percolumn(&,Π,c, k)
  # activation for the whole layer given the predicted/bursting neurons
  activate(A_pred, B)= (A_pred .| B')|> vec
  # cache parts of the output because they're used 2 times
  B= burst(c,Π)           # Nc
  A_pred= predicted(c,Π)  # k × Nc
  return activate(A_pred,B), B, A_pred|>vec
end

"""
`tm_predict(tm::TemporalMemory, α)` calculates which neurons will be predictive at the next step
given the currently active neurons + distal input `α`.
`size(a)` must match the presynaptic length of `tm.distalSynapses`.

# Returns

1. `Π`: predictive neurons (`Nₙ`)
2. `Πₛ`: predictive dendritic segments ('Nₛ') (caching)
3. `Mₛ`: matching dendritic segments (`Nₛ`) (learning)
4. `ovp_Mₛ`: subthreshold-matching dendritic segments (`Nₛ`) (learning)
"""
function tm_predict(tm::TemporalMemory, α::CellActivity)
  @unpack θ_stimulus_activate, θ_stimulus_learn = tm.params
  distal= tm.distalSynapses

  # Segment depolarization (prediction)
  Πₛ= Wd(distal)'α .> θ_stimulus_activate
  # Neuron depolarization (prediction)
  Π(Πₛ)= NS(distal)*Πₛ .> 0  # NOTE: params.θ_segment_act instead of 0
  # Sub-threshold segment stimulation sufficient for learning
  ovp_Mₛ= distal.Dd'α
  Mₛ= ovp_Mₛ .> θ_stimulus_learn
  return Π(Πₛ),Πₛ, Mₛ,ovp_Mₛ
end


text:
#module DynamicCorporateModels

using PyPlot
using Optim
#using BlackBoxOptim
using LinearAlgebra, Statistics
using Parameters, QuantEcon, DataFrames, Plots, Random
using Interpolations
using BenchmarkTools

#Here we start putting the Neoclassical models and add extensions
function NeoClassical(;
#Flavours
    cash::Bool=false,
    costly_equity::Bool=true,
    asymmetric_investment::Bool=false,
    quadratic_investment::Bool=true,
    risk_neutral::Bool=true,
    systematic::Bool=true,
    idiosyncratic::Bool=false,
#Parameters
    α::Float64=0.3, #Elasticity of capital
    r::Float64=0.01, #Risk free rate
    f::Float64=0.0, #Fixed cost
    τ::Float64=0.3, # Tax rate
    δ::Float64=0.01, # Depreciation rate
    x_bar::Float64=2.0, #Unconditional mean of the systematic shock
    ρ_x::Float64=0.9, #Persistence of Systematic Shock
    σ_x::Float64=0.1, #Volatility of systematic shock innovations
    ρ_z::Float64=0.99, #Persistence of idiosyncratic shock
    σ_z::Float64=0.1, #Volatility of idiosyncratic shock innovations
    λ::Float64=0.001, #Equity issuance cost

# Sizes
    n_k::Int64=30,
    n_x::Int64=10,
#Adjustment costs
    η::Float64=1.0, # Map i×k -> real)
    ) # End of arguments
# Here we create the matrices with the value function and policy function
    V::Matrix=zeros(n_k, n_x)
    i_opt::Matrix=zeros(n_k, n_x)
    mc = tauchen(n_x, ρ_x, σ_x, x_bar);
    Π_x=mc.p;
    x_grid = collect(range(x_bar-3*σ_x/sqrt(1-ρ_x^2), x_bar+3*σ_x/sqrt(1-ρ_x^2), length = n_x));
    k_grid = collect(range(1.0, 400.0, length=n_k))

    # Dividend payout in the simplest case
    Ψ(i,k)=k*0.5*η*(i/k)^2
    e(x, k, i)=(1-τ)*(exp(x)*k^α-f-δ*k)-i-Ψ(i,k)
    function d(x, k, i)
        earnings=e(x,k,i)
        return earnings+λ*abs(earnings)*(earnings<0)
    end

    function TV!(V, i_opt)
        V_func=interpolate((k_grid, x_grid), V, Gridded(Linear()));
        V_func=extrapolate(V_func, (Interpolations.Flat(),  Interpolations.Flat()));
        for (i_k, k) in enumerate(k_grid)
            for (i_x, x) in enumerate(x_grid)
                obj(i)=-d(x, k, i)-(1.0/(1+r))*sum(Π_x[i_x, j]*V_func(k*(1-δ)+i, x_grid[j]) for j in 1:n_x)
                # Here we get the optimal one, we can use one dimensional optimization
                res=Optim.optimize(obj, -maximum(k_grid), maximum(k_grid))
                V[i_k, i_x]=-res.minimum
                i_opt[i_k, i_x]=res.minimizer
            end
        end
    end

    function solve_VFI(V, i_opt; max_iter::Int64=1000, max_tol::Float64=1e-6)
        it=1
        tol=Inf
        while tol>max_tol && it<=max_iter
            V_old=copy(V)
            TV!(V, i_opt)
            tol=maximum(abs(x - y) for (x, y) in zip(V_old, V))
            it=it+1
        end
        return V, i_opt
    end

    function simulate_model(V_opt::Matrix, i_opt::Matrix, T::Int64=100)
        #Solves the model
        V_func=interpolate((k_grid, x_grid), V, Gridded(Linear()));
        V_func=extrapolate(V_func, (Interpolations.Flat(),  Interpolations.Flat()));

        I_func=interpolate((k_grid, x_grid), i_opt, Gridded(Linear()));
        I_func=extrapolate(I_func, (Interpolations.Flat(),  Interpolations.Flat()));

        x_init = mean(x_grid) # Initial x
        k_init = mean(k_grid)
        x_init_ind = searchsortedfirst(x_grid, x_init)

        k_sim=zeros(T+1,1)
        k_sim[1]=k_init
        i_sim=zeros(T+1,1)
        d_sim=zeros(T+1,1)
        mc_x = MarkovChain(Π_x)
        x_sim_indices = simulate(mc_x, T+1; init = x_init_ind)
        x_sim_val = zeros(T)
        #Pre return
        r_sim=ones(T+1,1)

        for t in 1:T
            current_x=x_grid[x_sim_indices[t]]
            current_k=k_sim[t]
            k_sim[t+1]=(1-δ)*current_k+I_func(current_k, current_x)
            i_sim[t]=I_func(current_k, current_x)
            d_sim[t]=d(current_x, current_k, i_sim[t])
            r_sim[t+1]=V_func(k_sim[t+1], x_grid[x_sim_indices[t+1]])/(V_func(k_sim[t], x_grid[x_sim_indices[t]])+d_sim[t])
        end

        # Here we plot the outcome of the model

        pygui(true)
        subplot(511)
        PyPlot.plot(1:T, x_grid[x_sim_indices][1:T])
        PyPlot.ylim(x_grid[1],x_grid[n_x])
        xlabel("Time")
        ylabel("X")
        subplot(512)
        PyPlot.plot(1:T, k_sim[1:T])
        PyPlot.ylim(minimum(k_sim),maximum(k_sim))
        xlabel("Time")
        ylabel("K")
        subplot(513)
        PyPlot.plot(1:T, i_sim[1:T])
        PyPlot.ylim(minimum(i_sim),maximum(i_sim))
        xlabel("Time")
        ylabel("I")
        subplot(514)
        PyPlot.plot(1:T, d_sim[1:T])
        PyPlot.ylim(minimum(d_sim),maximum(d_sim))
        xlabel("Time")
        ylabel("D")
        subplot(515)
        PyPlot.plot(1:T, r_sim[1:T])
        PyPlot.ylim(minimum(r_sim),maximum(r_sim))
        xlabel("Time")
        ylabel("R")
        show()

    end

    return (solve = solve_VFI, simulate_model=simulate_model, V=V, i_opt=i_opt)

end

neoclassical=NeoClassical()

solve_model=neoclassical.solve
simulate_model=neoclassical.simulate_model
V=neoclassical.V
i_opt=neoclassical.i_opt

V_opt, i_opt=solve_model(V, i_opt)

simulate_model(V, i_opt)

#end # module


text:
###################################################
## geo/timings/knn.jl
## Data-driven timings estimator, K-nearest-neighbors mean
###################################################

"""
    `KnnTimings`, K-nearest-neighbors timings estimator
"""
mutable struct KnnTimings <: GeoTimings
    trips::GeoData # compulsary attribute

    "number of neighbors to average"
    k::Int
    "kdd tree to search for nearest neighbors"
    tree::KDTree
    "center of projection"
    center::Tuple{Float32, Float32}
    "if we are weighting the mean by the distance"
    weighted::Bool
    "training set indices"
    trainset::Vector{Int}

    function KnnTimings(trips::GeoData, trainSet::AbstractArray{Int}, k::Int; weighted::Bool=false)
        gt = new()
        gt.trips = trips
        gt.k = k
        gt.weighted = weighted
        gt.trainset = trainSet
        # create center for projection
        minLat, maxLat, minLon, maxLon = Inf32, -Inf32, Inf32, -Inf32
        for t in trips
            minLat > t.pLat && (minLat = t.pLat)
            minLat > t.dLat && (minLat = t.dLat)
            minLon > t.pLon && (minLon = t.pLon)
            minLon > t.dLon && (minLon = t.dLon)
            maxLat < t.pLat && (maxLat = t.pLat)
            maxLat < t.dLat && (maxLat = t.dLat)
            maxLon < t.pLon && (maxLon = t.pLon)
            maxLon < t.dLon && (maxLon = t.dLon)
        end
        gt.center = ((maxLon + minLon)/2, (maxLat + minLat)/2)
        println("Creating KDTree...")
        dataPos = Array{Float64}(undef, (4, length(trainSet)))
        for (k,i) in enumerate(trainSet)
            if (k % 100_000) == 0
                @printf("\r%.2f%% trips loaded     ",100*k/length(trainSet))
            end
            px, py = toENU(trips[i].pLon, trips[i].pLat, gt.center)
            dx, dy = toENU(trips[i].dLon, trips[i].dLat, gt.center)
            dataPos[1, k] = px
            dataPos[2, k] = py
            dataPos[3, k] = dx
            dataPos[4, k] = dy
        end
        println("\r100.00% trips loaded      ")
        gt.tree = KDTree(dataPos)

        return gt
    end
end
KnnTimings(s::DataSplit, k::Int; args...) =  KnnTimings(s.geodata, trainSet(s), k; args...)

"""
    `estimateTime`: Implementation of required method for GeoTimings
"""
function estimateTime(gt::KnnTimings, pLon::Float32, pLat::Float32, dLon::Float32, dLat::Float32)
@inbounds begin
    px, py = toENU(pLon, pLat, gt.center)
    dx, dy = toENU(dLon, dLat, gt.center)
    idxs, dists = knn(gt.tree, [px, py, dx, dy], gt.k)
    if gt.weighted
        num = 1.0
        denum = 0.
        for i in eachindex(idxs)
            if dists[i] == 0.
                return gt.trips[gt.trainset[idxs[i]]].time
            else
                num *= gt.trips[gt.trainset[idxs[i]]].time ^ (1/ dists[i])
                denum += 1.0/dists[i]
            end
        end
        return num^(1/denum)
    else
        mean = 1.0
        for id in idxs
            mean *= gt.trips[gt.trainset[id]].time
        end
        return mean^(1.0/length(idxs))
    end
end #inbounds
end

"""
    `estimateTimeWithDistance`: Alternate version of estimateTime that also returns average distance from point to neighbors
"""
function estimateTimeWithDistance(gt::KnnTimings, pLon::Float32, pLat::Float32, dLon::Float32, dLat::Float32)
@inbounds begin
    px, py = toENU(pLon, pLat, gt.center)
    dx, dy = toENU(dLon, dLat, gt.center)
    idxs, dists = knn(gt.tree, [px, py, dx, dy], gt.k)
    if gt.weighted
        num = 1.0
        denum = 0.
        for i in eachindex(idxs)
            if dists[i] == 0.
                return gt.trips[gt.trainset[idxs[i]]].time
            else
                num *= gt.trips[gt.trainset[idxs[i]]].time ^ (1/ dists[i])
                denum += 1.0/dists[i]
            end
        end
        return num^(1/denum), sum(dists)/length(dists)
    else
        mean = 1.0
        for id in idxs
            mean *= gt.trips[gt.trainset[id]].time
        end
        return mean^(1/length(idxs)), sum(dists)/length(dists)
    end
end #inbounds
end

estimateTimeWithDistance(gt::KnnTimings, tripID::Int) = estimateTimeWithDistance(gt, gt.trips[tripID].pLon, gt.trips[tripID].pLat, gt.trips[tripID].dLon, gt.trips[tripID].dLat)


text:
using GRIB
grib_file = "/home/neo/Projekte/Predly/weather/gfs_workflow/download/data/gfs.20210109/06/gfs.t06z.pgrb2.0p25.f002"
function check_messages(grib_file)
    messages = 0
    GribFile(grib_file) do f
        for message in f
            messages+=1
        end
    end
    println(messages)
end

function check_keys(grib_file)
    GribFile(grib_file) do f
        for message in f
            for key in keys(message)
                if key ∉ keylist
                push!(keylist, key)
                end
            end
        end
    end
end

check_messages(grib_file)


text:
function LofOffset(arg0::Frame, arg1::LOFType)
    return LofOffset((Frame, LOFType), arg0, arg1)
end

function LofOffset(arg0::Frame, arg1::LOFType, arg2::RotationOrder, arg3::jdouble, arg4::jdouble, arg5::jdouble)
    return LofOffset((Frame, LOFType, RotationOrder, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function get_attitude(obj::LofOffset, arg0::FieldPVCoordinatesProvider, arg1::FieldAbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", FieldAttitude, (FieldPVCoordinatesProvider, FieldAbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_attitude(obj::LofOffset, arg0::PVCoordinatesProvider, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", Attitude, (PVCoordinatesProvider, AbsoluteDate, Frame), arg0, arg1, arg2)
end


text:
# This file defines the Transition type. This type defines a transition between
#   nodes that a personnel member can make.

export Transition
"""
This type defines a transition between two nodes that a personnel member can perform, along with all the necessary information  about the conditions of the transition.

The type contains the following fields:
* `name::string`: the name of the transition.
* `sourceNode::String`: the name of the node that the personnel member is currently in.
* `targetNode::String`: the name of the node that the personnel member can attain.
* `isOutTransition::Bool`: a flag to indicate if this is a transition out of the organisation. If this is `true`, the end node gets ignored.
* `freq::Float64`: the time between two checks in the transition's schedule. Default = 1.0
* `offset::Float64`: the offset of the transition's schedule with respect to the start of the simulation. Default = 0.0
* `maxAttempts::Int`: the maximum number of tries a personnel member has to undergo the transition. A value of 0 means as many as there are entries in the list of transition probabilites, and a value of -1 means there is no maximum. Default = 1
* `minFlux::Int`: the minimum number of people that must undergo the transition at the same time, if this many people are eligible. Default = 0
* `maxFlux::Int`: the maximum number of people that can undergo the transition at the same time. A value of -1 means there is no maximum. Default = Inf
* `hasPriority::Bool`: a flag indicating that this transition can override the target population of the transition's target node. If the flag is `true`, it means that if the max flux of the node is 15, and only 10 spots are available in the target node, 15 people will undergo the transition nonetheless. If the flag is `false`, 10 persons would.
* `extraConditions::Vector{MPcondition}`: the extra conditions that must be satisfied before the transition can take place.
* `extraChanges::Dict{String,String}`: the extra changes to attributes that happen during the transition.
* `probabilityList::Vector{Float64}`: the list of probabilities for this transition to occur.

The type has two extra field, computed at the start of the simulation:
* `priority::Int`: the priority in the simulation on which the transition gets executed. This priority will be 0 (or > 0) for transitions with the `hasPriority` flag set to `true`, and < 0 otherwise. A priority == 1 means it needs to be determined first.
* `probRNG`: the random stream generating values to determine whhich of the eligible personnel members does the transition.

Constructors:
```
Transition(
    name::String,
    sourceNode::String,
    targetNode::String )
```
This constructor creates a `Transition` object with source node `sourceNode` and target node `targetNode`.

```
Transition(
    name::String,
    sourceNode::String )
```
This constructor creates a `Transition` object with source node `sourceNode` and `isOutTrransition` flag to `true`.
"""
mutable struct Transition

    name::String
    sourceNode::String
    targetNode::String
    isOutTransition::Bool
    freq::Float64
    offset::Float64
    maxAttempts::Int
    minFlux::Int
    maxFlux::Int
    hasPriority::Bool
    extraConditions::Vector{MPcondition}
    extraChanges::Dict{String,String}
    probabilityList::Vector{Float64}

    priority::Int
    probRNG::MersenneTwister

    # Basic constructor.
    function Transition( name::String, sourceNode::String,
        targetNode::String )::Transition

        newTrans = new()
        newTrans.name = name
        newTrans.sourceNode = lowercase( sourceNode ) == "dummy" ?
            "dummy" : sourceNode
        newTrans.targetNode = lowercase( targetNode ) == "dummy" ?
            "dummy" : targetNode
        newTrans.isOutTransition = false
        newTrans.freq = 1.0
        newTrans.offset = 0.0
        newTrans.maxAttempts = 1
        newTrans.minFlux = 0
        newTrans.maxFlux = -1
        newTrans.hasPriority = false
        newTrans.extraConditions = Vector{MPcondition}()
        newTrans.extraChanges = Dict{String,String}()
        newTrans.probabilityList = [1.0]
        newTrans.priority = 1
        newTrans.probRNG = MersenneTwister()
        return newTrans

    end  # Transition( name, sourceNode, targetNode, freq, offset, maxAttempts,
         #   minFlux, maxFlux, hasPriority )

    function Transition( name::String, sourceNode::String )::Transition

        newTrans = Transition( name, sourceNode, "dummy" )
        newTrans.isOutTransition = true
        return newTrans

    end  # Transition( name, sourceNode )

end  # mutable struct Transition


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module PhysQuantities

export PhysQuantity

# Mimic a quantity with a physical unit that is not convertible to a real number
struct PhysQuantity{n,T}   # n is like the exponent of the unit
    val::T
end
PhysQuantity{n}(x::T) where {n,T} = PhysQuantity{n,T}(x)
Base.zero(::Type{PhysQuantity{n,T}}) where {n,T} = PhysQuantity{n,T}(zero(T))
Base.zero(x::PhysQuantity) = zero(typeof(x))
Base.:+(x::PhysQuantity{n}, y::PhysQuantity{n}) where n = PhysQuantity{n}(x.val + y.val)
Base.:-(x::PhysQuantity{n}, y::PhysQuantity{n}) where n = PhysQuantity{n}(x.val - y.val)
Base.:*(x::PhysQuantity{n,T}, y::Int) where {n,T} = PhysQuantity{n}(x.val*y)
Base.:/(x::PhysQuantity{n,T}, y::Int) where {n,T} = PhysQuantity{n}(x.val/y)
Base.:*(x::PhysQuantity{n1,S}, y::PhysQuantity{n2,T}) where {n1,n2,S,T} =
    PhysQuantity{n1+n2}(x.val*y.val)
Base.:/(x::PhysQuantity{n1,S}, y::PhysQuantity{n2,T}) where {n1,n2,S,T} =
    PhysQuantity{n1-n2}(x.val/y.val)
Base.convert(::Type{PhysQuantity{0,T}}, x::Int) where T = PhysQuantity{0}(convert(T, x))
Base.convert(::Type{P}, ::Int) where P<:PhysQuantity =
    error("Int is incommensurate with PhysQuantity")

end #module

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
# %% imports

import Cambrian: run!
export run!

using Cambrian
using JLD

# %% save generation's best individual
function step!(e::AbstractEvolution)
    e.gen += 1
    if e.gen > 1
        populate(e)
    end
    evaluate(e)
    generation(e)
    if ((e.config.log_gen > 0) && mod(e.gen, e.config.log_gen) == 0)
        log_gen(e)
    end
    if ((e.config.save_gen > 0) && mod(e.gen, e.config.save_gen) == 0)
        save_gen(e)
    end
end

#if sNES
# function run!(e::AbstractEvolution, opt::Int64)
#     for i in (e.gen+1):e.config.n_gen
#         step!(e)
#         for idx in eachindex(e.population)
#             if e.population[idx].fitness[1] > 195.
#                 f = jldopen("Julia/Actors/topgensNES.jld", "r+")
#                 write(f, "$(opt)_run", e.gen)
#                 close(f)
#                 return 0
#             end
#         end
#     end
# end

#if xNES
function run!(e::AbstractEvolution, opt::Int64)
    for i in (e.gen+1):e.config.n_gen
        step!(e)
        for idx in eachindex(e.population)
            if e.population[idx].fitness[1] > 195.
                f = jldopen("Julia/Actors/topgenxNES.jld", "r+")
                write(f, "$(opt)_run", e.gen)
                close(f)
                return 0
            end
        end
    end
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
function add_polygon!(m::Model,xs,ys;lc=0.,z=0,transfinite=0,holes=Vector{LineLoop}(0))
    # Add the points
    pts = add_point!.([m],xs,ys,[z],[lc;])
    # Add the lines
    lns = add_line!.([m],pts,circshift(pts,-1),[transfinite;])
    # Add the lineloop
    ll  = add_lineloop!(m,lns)
    # Add and return the surface
    return add_surface!(m,[ll;holes])
end

function add_rectangle!(m::Model,x1,y1,x2,y2;lc=0.,z=0,transfinite=0,holes=Vector{LineLoop}(0))
    return add_polygon!(m,[x1,x2,x2,x1],[y1,y1,y2,y2],lc=lc,z=z,transfinite=transfinite,holes=holes)
end

max_stars_repo_name: jebej/FastCap

max_issues_repo_name: jebej/FastCap

max_forks_repo_name: jebej/FastCap


text:
for op in (:(<=), :(<), :(==), :(!=), :(?), :(>=), :isless, :isequal)
    @eval begin
        ($op){T<:Float65}(a::T, b::Float65) = ($op)(reflect(a.fp), reflect(b.fp))
        ($op){T<:Float65}(a::T, b::Float64) = ($op)(reflect(a.fp), b)
        ($op){T<:Float65}(a::Float64, b::T) = ($op)(a, reflect(b.fp))

        ($op){T<:Float65,I<:Union{Int64,Int32}}(a::T, b::I) = ($op)(reflect(a.fp), b)
        ($op){T<:Float65,I<:Union{Int64,Int32}}(a::I, b::T) = ($op)(a, reflect(b.fp))
        ($op){T<:Float65}(a::T, b::Bool)      = ($op)(reflect(a.fp), convert(Float64,b))
        ($op){T<:Float65}(a::Bool, b::T)      = ($op)(convert(Float64,a), reflect(b.fp))
        ($op){T<:Float65,R<:Real}(a::T, b::R) = ($op)(reflect(a.fp), convert(Float64,b))
        ($op){T<:Float65,R<:Real}(a::R, b::T) = ($op)(convert(Float64,a), reflect(b.fp))
    end
end


text:
module NumberTheoreticTransforms

include("ntt.jl")
include("fnt.jl")

end # module


text:
#### Domain decomposition indexing
export DomainDecomp

import Base
abstract type Domain{N} end
struct GridMean{N} <: Domain{N} end
struct Environment{N} <: Domain{N} end
struct Updraft{N} <: Domain{N} end

"""
    get_param(::Type{Domain})

The number of domains in domain `Domain`
"""
get_param(::Type{GridMean{N}}) where N = N
get_param(::Type{Environment{N}}) where N = N
get_param(::Type{Updraft{N}}) where N = N

"""
    DomainDecomp{GM,EN,UD}

Decomposition of number of domains
 - `gm` number of grid-mean domain
 - `en` number of environment sub-domain
 - `ud` number of updraft sub-domains
"""
struct DomainDecomp{GM,EN,UD}
  function DomainDecomp(;gm::GM=0,en::EN=0,ud::UD=0) where {GM<:Int,EN<:Int,UD<:Int}
    @assert 0<=gm<=1
    @assert 0<=en<=1
    @assert 0<=ud
    return new{GridMean{gm},Environment{en},Updraft{ud}}()
  end
end

Base.sum(   ::DomainDecomp{GM,EN,UD}) where {GM,EN,UD} = get_param(GM)+get_param(EN)+get_param(UD)
get_param(  ::DomainDecomp{GM,EN,UD}) where {GM,EN,UD} = (get_param(GM),get_param(EN),get_param(UD))

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
mutable struct CCMaximinObjectiveCore
    alpha::Float64
    points::Vector{Float64}
    IIF::Array{Float64,3}
    CCMaximinObjectiveCore() =
        new(0.05, Vector{Float64}(), Array{Float64,3}(undef, 0, 0, 0))
end


text:
function _display(plt::Plot{PyPlotBackend})
    plt.o[:show]()
end

_pyplot_mimeformats = Dict(
    "application/eps"         => "eps",
    "image/eps"               => "eps",
    "application/pdf"         => "pdf",
    "image/png"               => "png",
    "application/postscript"  => "ps",
    "image/svg+xml"           => "svg"
)

for (mime, fmt) in _pyplot_mimeformats
    @eval function _show(io::IO, ::MIME{Symbol($mime)}, plt::Plot{PyPlotBackend})
        fig = plt.o
        fig[:canvas][:print_figure](
            io,
            format=$fmt,
            # bbox_inches = "tight",
            # figsize = map(px2inch, plt[:size]),
            facecolor = fig[:get_facecolor](),
            edgecolor = "none",
            dpi = plt[:dpi] * plt[:thickness_scaling]
        )
    end
end

closeall(::PyPlotBackend) = PyPlot.plt[:close]("all")


text:
# This file was generated, do not modify it. # hide
fix_blanks(v) = map(v) do x
    if x == " "
        return "0.0"
    else
        return x
    end
end

df0.TotalCharges = fix_blanks(df0.TotalCharges);

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
@doc raw"""
    device_commitment(canonical_model::CanonicalModel,
                        initial_conditions::Vector{InitialCondition},
                        cons_name::Symbol,
                        var_names::Tuple{Symbol, Symbol, Symbol})

Constructs multi-timestep constraint from initial conditions and binary variable tuple.

# Constraints

``` varstart + varstop <= 1.0 ```

If t = 1:

``` varon[name, 1] == ic.value + varstart[name, 1] - varstop[name, 1] ```

where ic in initial_condtions.

If t > 1:

``` varon[name, t] == varon[name, t-1] + varstart[name, t] - varstop[name, t] ```

# LaTeX

`` x^{on}_t + x^{off}_t \leq 1.0 \forall t ``

`` x^{on}_1 = x^{on}_{init} + x^{start}_1 - x^{stop}_1, \text{ for } t = 1 ``

`` x^{on}_t = x^{on}_{t-1} + x^{start}_t - x^{stop}_t, \forall t \geq 2 ``

# Arguments
* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* initial_conditions::Vector{InitialCondition} : for time zero 'varon'
* cons_name::Symbol : name of the constraint
* var_names::Tuple{Symbol, Symbol, Symbol} : the names of the variables
-  : var_names[1] : varstart
-  : var_names[2] : varstop
-  : var_names[3] : varon
"""
function device_commitment(canonical_model::CanonicalModel,
                        initial_conditions::Vector{InitialCondition},
                        cons_name::Symbol,
                        var_names::Tuple{Symbol, Symbol, Symbol})

    time_steps = model_time_steps(canonical_model)
    varstart = var(canonical_model, var_names[1])
    varstop = var(canonical_model, var_names[2])
    varon = var(canonical_model, var_names[3])
    varstart_names = axes(varstart, 1)
    _add_cons_container!(canonical_model, cons_name, varstart_names, time_steps)
    constraint = con(canonical_model, cons_name)
    aux_cons_name = _middle_rename(cons_name, "_", "aux")
    _add_cons_container!(canonical_model, aux_cons_name, varstart_names, time_steps)
    aux_constraint = con(canonical_model, aux_cons_name)

    for ic in initial_conditions
        name = PSY.get_name(ic.device)
        constraint[name, 1] = JuMP.@constraint(canonical_model.JuMPmodel,
                               varon[name, 1] == ic.value + varstart[name, 1] - varstop[name, 1])
        aux_constraint[name, 1] = JuMP.@constraint(canonical_model.JuMPmodel,
                               varstart[name, 1] + varstop[name, 1] <= 1.0)
    end

    for t in time_steps[2:end], i in initial_conditions
        name = PSY.get_name(i.device)
        constraint[name, t] = JuMP.@constraint(canonical_model.JuMPmodel,
                        varon[name, t] == varon[name, t-1] + varstart[name, t] - varstop[name, t])
        aux_constraint[name, t] = JuMP.@constraint(canonical_model.JuMPmodel,
                                varstart[name, t] + varstop[name, t] <= 1.0)
    end

    return

end


text:
# src is either a string representing the path to the file, or an IO stream
mutable struct SndFileSource{T, S<:Union{String, LengthIO}} <: SampleSource
    src::S
    filePtr::Ptr{Cvoid}
    sfinfo::SF_INFO
    pos::Int64
    readbuf::Array{T, 2}
end

function SndFileSource(src, filePtr, sfinfo, bufsize=4096)
    T = fmt_to_type(sfinfo.format)
    readbuf = zeros(T, sfinfo.channels, bufsize)

    SndFileSource(src, filePtr, sfinfo, Int64(1), readbuf)
end

mutable struct SndFileSink{T, S<:Union{String, LengthIO}} <: SampleSink
    src::S # the stream or filename used to create this
    filePtr::Ptr{Cvoid}
    sfinfo::SF_INFO
    nframes::Int64
    writebuf::Array{T, 2}
end

function SndFileSink(src, filePtr, sfinfo, bufsize=4096)
    T = fmt_to_type(sfinfo.format)
    writebuf = zeros(T, sfinfo.channels, bufsize)
    SndFileSink(src, filePtr, sfinfo, Int64(0), writebuf)
end

const SndFileStream{T, S} = Union{SndFileSource{T, S}, SndFileSink{T, S}}

SampledSignals.nchannels(str::SndFileStream) = Int(str.sfinfo.channels)
SampledSignals.samplerate(str::SndFileStream) = str.sfinfo.samplerate
Base.eltype(str::SndFileStream{T, S}) where {T, S} = T

SampledSignals.nframes(sink::SndFileSink) = sink.nframes
SampledSignals.nframes(source::SndFileSource) = source.sfinfo.frames

function Base.close(s::SndFileStream)
    if s.filePtr != C_NULL
        sf_close(s.filePtr)
        s.filePtr = C_NULL
    else
        @warn "close called more than once on $s"
    end
end

function Base.show(io::IO, s::SndFileStream)
    println(io, typeof(s))
    println(io, "  path: \"$(s.src)\"")
    println(io, "  channels: ", nchannels(s))
    println(io, "  samplerate: ", samplerate(s), "Hz")
    # SndFileSinks don't have a position and we're always at the end
    pos = isa(s, SndFileSource) ? s.pos-1 : nframes(s)
    postime = float((pos)/samplerate(s))
    endtime = float((nframes(s))/samplerate(s))
    println(io, "  position: $(pos) of $(nframes(s)) frames")
    @printf(io, "            %0.2f of %0.2f seconds", postime, endtime)
end

function _seek(source::SndFileSource, offset::Integer, whence::Integer)
    new_offset = sf_seek(source.filePtr, offset, whence)

    if new_offset < 0
        error("Could not seek to $(offset) in file")
    else
        # switch to one-indexing
        source.pos = new_offset + 1
    end

    new_offset
end

# Some convenience methods for easily navigating through a sound file
Base.seek(source::SndFileSource, offset::Integer) = _seek(source, offset-1, SF_SEEK_SET)
Base.skip(source::SndFileSource, offset::Integer) = _seek(source, offset, SF_SEEK_CUR)
Base.position(source::SndFileSource) = source.pos