text:
# ========== Transforms On A Group Of Signals ==========
# ----- Discrete Wavelet Transform on a set of signals -----
"""
    dwtall(x, wt[, L])

Computes the discrete wavelet transform (DWT) on each slice of signal. Signals are sliced on
the ``n``-th dimension for an ``n``-dimensional input `x`.

!!! note
    `dwt` is currently available for 1-D, 2-D, and 3-D signals only.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Input signals, where each slice corresponds to
  one signal. For a set of input signals `x` of dimension ``n``, signals are sliced on the
  ``n``-th dimension.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `Wavelets.maxtransformlevels(xᵢ)`) Number of levels of wavelet
  transforms.

# Returns
`::Array{T}`: Slices of transformed signals. Signals are sliced the same way as the input
signal `x`.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# DWT on all signals in x
xw = dwtall(x, wt)
```

**See also:** [`idwtall`](@ref), [`wpdall`](@ref), [`wptall`](@ref)
"""
function dwtall(x::AbstractArray{T}, args...) where T<:Number
    # Sanity check
    @assert ndims(x) > 1

    # Allocate space for transforms
    y = similar(x)
    # Dimension to slice
    dim = ndims(y)
    # Compute transforms
    @inbounds begin
        @views for (yᵢ, xᵢ) in zip(eachslice(y, dims=dim), eachslice(x, dims=dim))
            dwt!(yᵢ, xᵢ, args...)
        end
    end
    return y
end

"""
    idwtall(xw, wt[, L])

Computes the inverse discrete wavelet transform (iDWT) on each slice of signal. Signals are
sliced on the ``n``-th dimension for an ``n``-dimensional input `xw`.

!!! note
    `idwt` is currently available for 1-D, 2-D, and 3-D signals only.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: Input decomposed signals, where each slice
  corresponds to one signal. For a set of input signals `xw` of dimension ``n``, signals are
  sliced on the ``n``-th dimension.
- `wt::OrthoFilter`: Wavelet used.
- `L::Integer`: (Default: `Wavelets.maxtransformlevels(xwᵢ)`) Number of levels of wavelet
  transforms.

# Returns
`::Array{T}`: Slices of reconstructed signals. Signals are sliced the same way as the input
`xw`.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# DWT on all signals in x
xw = dwtall(x, wt)

# iDWT on all signals
x̂ = idwtall(xw, wt)
```

**See also:** [`dwtall`](@ref), [`iwpdall`](@ref), [`iwptall`](@ref)
"""
function idwtall(xw::AbstractArray{T}, args...) where T<:Number
    # Sanity check
    @assert ndims(xw) > 1

    # Setup
    y = similar(xw)      # Allocation for output
    dim = ndims(xw)      # Dimension to slice

    # IDWT
    @inbounds begin
        @views for (yᵢ, xwᵢ) in zip(eachslice(y, dims=dim), eachslice(xw, dims=dim))
            idwt!(yᵢ, xwᵢ, args...)
        end
    end
    return y
end

# ----- Wavelet Packet Transforms on a set of signals -----
"""
    wptall(x, wt[, L; standard])
    wptall(x, wt, tree[; standard])

Computes the wavelet packet transform (WPT) on each slice of signal. Signals are sliced on
the ``n``-th dimension for an ``n``-dimensional input `x`.

!!! note
    `wpt` is currently available for 1-D and 2-D signals only.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Input signals, where each slice corresponds to
  one signal. For a set of input signals `x` of dimension ``n``, signals are sliced on the
  ``n``-th dimension.
- `wt::OrthoFilter`: Wavelet used.
- `L::Integer`: (Default: `Wavelets.maxtransformlevels(xᵢ)`) Number of levels of wavelet
  decomposition.
- `tree::BitVector`: (Default: `Wavelets.maketree(xᵢ, :full)`) Tree to follow for wavelet
  decomposition. Default value is only applicable for 1D signals.

# Keyword Arguments
- `standard::Bool`: (Default: `true`) Whether to compute the standard or non-standard
  wavelet transform. Only applicable for 2D signals.

# Returns
`::Array{T}`: Slices of transformed signals. Signals are sliced the same way as the input
signal `x`.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# WPT on all signals in x
xw = wptall(x, wt)
```

**See also:** [`wpdall`](@ref), [`dwtall`](@ref), [`wpt`](@ref)
"""
function wptall(x::AbstractArray{T}, args...; kwargs...) where T<:Number
    # Sanity check
    @assert ndims(x) > 1

    # Setup
    y = similar(x)      # Allocation for output
    dim = ndims(x)      # Dimension to slice
    # Compute transforms
    @inbounds begin
        @views for (yᵢ, xᵢ) in zip(eachslice(y, dims=dim), eachslice(x, dims=dim))
            wpt!(yᵢ, Array(xᵢ), args..., kwargs...)
        end
    end
    return y
end

"""
    iwptall(xw, wt[, L; standard])
    iwptall(xw, wt, tree[; standard])

Computes the inverse wavelet packet transform (iWPT) on each slice of signal. Signals are
sliced on the ``n``-th dimension for an ``n``-dimensional input `xw`.

!!! note
    `iwpt` is currently available for 1-D and 2-D signals only.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Input signals, where each slice corresponds to one
  signal. For a set of input signals `x` of dimension ``n``, signals are sliced on the
  ``n``-th dimension.
- `wt::OrthoFilter`: Wavelet used.
- `L::Integer`: (Default: `Wavelets.maxtransformlevels(xᵢ)`) Number of levels of wavelet
  decomposition.
- `tree::BitVector`: (Default: `Wavelets.maketree(xᵢ, :full)`) Tree to follow for wavelet
  decomposition. Default value is only applicable for 1D signals.

# Keyword Arguments
- `standard::Bool`: (Default: `true`) Whether to compute the standard or non-standard
  wavelet transform. Only applicable for 2D signals.

# Returns
`::Array{T}`: Slices of transformed signals. Signals are sliced the same way as the input
signal `x`.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# WPT on all signals in x
xw = wptall(x, wt)

# iWPT on all signals on x
x̂ = iwptall(xw, wt)
```

**See also:** [`wpdall`](@ref), [`dwtall`](@ref), [`wpt`](@ref)
"""
function iwptall(xw::AbstractArray{T}, args...; kwargs...) where T<:Number
    # Sanity check
    @assert ndims(xw) > 1

    # Setup
    y = similar(xw)      # Allocation for output
    dim = ndims(xw)      # Dimension to slice

    # IWPT
    @inbounds begin
        @views for (yᵢ, xwᵢ) in zip(eachslice(y, dims=dim), eachslice(xw, dims=dim))
            iwpt!(yᵢ, Array(xwᵢ), args..., kwargs...)
        end
    end
    return y
end

# ----- Wavelet Packet Decomposition on a set of signals -----
"""
    wpdall(x, wt[, L; standard])

Computes the wavelet packet decomposition (WPD) on each slice of signal. Signals are sliced
on the ``n``-th dimension for an ``n``-dimensional input `x`.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Input signals, where each slice corresponds to
  one signal. For a set of input signals `x` of dimension ``n``, signals are sliced on the
  ``n``-th dimension.
- `wt::OrthoFilter`: Wavelet used.
- `L::Integer`: (Default: `Wavelets.maxtransformlevels(xᵢ)`) Number of levels of wavelet
  decomposition.

# Keyword Arguments
- `standard::Bool`: (Default: `true`) Whether to compute the standard or non-standard
  wavelet transform. Only applicable for 2D signals.

# Returns
`::Array{T}`: Array of decomposed signals. Signals are sliced by the final dimension.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# WPT on all signals in x
xw = wpdall(x, wt)
```

**See also:** [`wptall`](@ref), [`dwtall`](@ref), [`wpd`](@ref), [`wpd!`](@ref)
"""
function wpdall(x::AbstractArray{T},
                wt::OrthoFilter,
                L::Integer = maxtransformlevels(x,1);
                kwargs...) where T<:Number
    # Sanity check
    @assert ndims(x) > 1
    @assert 0 ≤ L ≤ maxtransformlevels(selectdim(x, ndims(x), 1))

    # Allocate space for transform
    sz = size(x)[begin:end-1]  # Signal size
    N = size(x)[end]           # Number of signals
    y = Array{T}(undef, (sz...,L+1,N))
    # Dimension to slice
    dim_y = ndims(y)
    dim_x = ndims(x)
    # Compute transforms
    @inbounds begin
        @views for (yᵢ, xᵢ) in zip(eachslice(y, dims=dim_y), eachslice(x, dims=dim_x))
            wpd!(yᵢ, Array(xᵢ), wt, L, kwargs...)
        end
    end
    return y
end

"""
    iwpdall(xw, wt[, L; standard])
    iwpdall(xw, wt, tree[; standard])

Computes the inverse wavelet packet decomposition (IWPD) on each slice of signal. Signals
are sliced on the ``n``-th dimension for an ``n``-dimensional input `xw`.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: Input signals, where each slice corresponds to one
  signal decomposition. For a set of input signals `x` of dimension ``n``, `xw` is sliced on
  the ``n``-th dimension.
- `wt::OrthoFilter`: Wavelet used.
- `L::Integer`: (Default: `Wavelets.maxtransformlevels(xᵢ)`) Number of levels of wavelet
  decomposition.

# Keyword Arguments
- `standard::Bool`: (Default: `true`) Whether to compute the standard or non-standard
  wavelet transform. Only applicable for 2D signals.

# Returns
`::Array{T}`: Array of decomposed signals. Signals are sliced by the final dimension.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# WPD on all signals in x
xw = wpdall(x, wt)

# IWPD on all signals
x̂ = iwpdall(xw, wt)
```

**See also:** [`wptall`](@ref), [`dwtall`](@ref), [`wpd`](@ref), [`wpd!`](@ref)
"""
function iwpdall(xw::AbstractArray{T}, args...; kwargs...) where T<:Number
    # Sanity check
    @assert ndims(xw) > 2

    # Setup
    sz = size(xw)[1:(end-2)]
    N = size(xw)[end]
    x̂ = Array{T}(undef, (sz...,N))
    # Dimension to slice
    dim_xw = ndims(xw)
    dim_x̂ = ndims(x̂)
    # Compute transforms
    @inbounds begin
        @views for (x̂ᵢ, xwᵢ) in zip(eachslice(x̂, dims=dim_x̂), eachslice(xw, dims=dim_xw))
            iwpd!(x̂ᵢ, Array(xwᵢ), args..., kwargs...)
        end
    end
    return x̂
end


text:
using Documenter, Pkg, TwoBit

DocMeta.setdocmeta!(TwoBit, :DocTestSetup, :(using TwoBit); recursive=true)

makedocs(
    format = Documenter.HTML(
        prettyurls = haskey(ENV, "GITHUB_ACTIONS"),
    ),
    modules = [TwoBit],
    sitename = "TwoBit.jl",
    doctest = true,
    strict = false,
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Reading files" => "man/reading.md",
            "Writing files" => "man/writing.md",
            "TwoBit Records" => "man/records.md"
        ],
    ],
    authors = replace(join(Pkg.TOML.parsefile("Project.toml")["authors"], ", "), r" <.*?>" => "" ) * ", The BioJulia Organisation, and other contributors.",
)

deploydocs(
    repo = "github.com/BioJulia/TwoBit.jl.git",
    push_preview = true
)

ext: jl


text:
using OrdinaryDiffEq,DiffEqProblemLibrary, DiffEqDevTools

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear

prob = prob_ode_2Dlinear
sol =solve(prob,Rosenbrock32(),dt=1/2^4)

sol =solve(prob,ExplicitRK(tableau = constructBogakiShampine3()),dt=1/2^4)
val1 = maximum(abs.(sol.u[end] - sol.u_analytic[end]))

sol2 =solve(prob,ExplicitRK(tableau = constructDormandPrince()),dt=1/2^4)
val2 = maximum(abs.(sol2.u[end] - sol2.u_analytic[end]))

sol3 =solve(prob,ExplicitRK(tableau = constructRKF8(Float64)),dt=1/2^4)
val3 = maximum(abs.(sol3.u[end] - sol3.u_analytic[end]))

@test length(sol.t)>length(sol2.t)>=length(sol3.t)
@test max(val1,val2,val3)<2e-3


text:
immutable Spiral
    m::Int
    n::Int
    cmax::Int
    dir::Array{Array{Int,1},1}
    bdelta::Array{Array{Int,1},1}
end

function Spiral(m::Int, n::Int)
    cmax = m*n
    dir = Array{Int,1}[[0,1], [1,0], [0,-1], [-1,0]]
    bdelta = Array{Int,1}[[0,0,0,1], [-1,0,0,0],
                          [0,-1,0,0], [0,0,1,0]]
    Spiral(m, n, cmax, dir, bdelta)
end

function spiral(m::Int, n::Int)
    0<m&&0<n || error("The matrix dimensions must be positive.")
    Spiral(m, n)
end
spiral(n::Int) = spiral(n, n)

type SpState
    cnt::Int
    dirdex::Int
    cell::Array{Int,1}
    bounds::Array{Int,1}
end

Base.length(sp::Spiral) = sp.cmax
Base.start(sp::Spiral) = SpState(1, 1, [1,1], [sp.n,sp.m,1,1])
Base.done(sp::Spiral, sps::SpState) = sps.cnt > sp.cmax

function Base.next(sp::Spiral, sps::SpState)
    s = sub2ind((sp.m, sp.n), sps.cell[1], sps.cell[2])
    if sps.cell[rem1(sps.dirdex+1, 2)] == sps.bounds[sps.dirdex]
        sps.bounds += sp.bdelta[sps.dirdex]
        sps.dirdex = rem1(sps.dirdex+1, 4)
    end
    sps.cell += sp.dir[sps.dirdex]
    sps.cnt += 1
    return (s, sps)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import Base: rand, randn
export srand

function rand(::Type{AFArray{T,N}}, t::NTuple{N,Int}) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_randu,af_lib),af_err,(Ptr{af_array},UInt32,Ptr{dim_t},af_dtype),
                 out,UInt32(N),[t...],af_type(T)))
    AFArray{T,N}(out[])
end
rand(::Type{AFArray{T}}, t::NTuple{N,Int}) where {T,N}= rand(AFArray{T,N}, t)
rand(::Type{AFArray}, t::NTuple{N,Int}) where {N} = rand(AFArray{Float32,N}, t)
rand(::Type{AFArray{T}}, t::Int...) where {T} = rand(AFArray{T}, t)
rand(::Type{AFArray}, t::Int...) = rand(AFArray, t)

function randn(::Type{AFArray{T,N}}, t::NTuple{N,Int}) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_randn,af_lib),af_err,(Ptr{af_array},UInt32,Ptr{dim_t},af_dtype),
                 out,UInt32(N),[t...],af_type(T)))
    AFArray{T,N}(out[])
end
randn(::Type{AFArray{T}}, t::NTuple{N,Int}) where {T,N} = randn(AFArray{T,N}, t)
randn(::Type{AFArray}, t::NTuple{N,Int}) where {N} = randn(AFArray{Float32,N}, t)
randn(::Type{AFArray{T}}, t::Int...) where {T} = randn(AFArray{T}, t)
randn(::Type{AFArray}, t::Int...) = randn(AFArray, t)

srand(::Type{AFArray}, i) = set_seed(UInt(abs(i)))


text:
function ptype_disksize_nodefaults(h::IAEAHeader{Nf, Ni}) where {Nf, Ni}
    1 + # typ
    4 + # energy
    12 + # x,y,z
    8 + # u,v (w is not stored)
    4 + # weight
    4 * Nf +
    4 * Ni
end

function ptype_disksize(h::IAEAHeader{Nf, Ni}) where {Nf, Ni}
    size_reduction_due_to_defaults = sizeof(h.record_contents)
    ptype_disksize_nodefaults(h) - size_reduction_due_to_defaults
end

@generated function read_default(io::IO,
                             ::Val{field},
                             h::IAEAHeader{Nf,Ni,NT}) where {field,Nf,Ni,NT}
    if field in fieldnames(NT)
        :(h.record_contents.$field)
    else
        :(read_(io, Float32))
    end
end

@generated function write_default(io::IO,
                             ::Val{field},
                             p::IAEAParticle,
                             h::IAEAHeader{Nf,Ni,NT}) where {field,Nf,Ni,NT}
    if field in fieldnames(NT)
        quote
            @assert p.$field == h.record_contents.$field
            0
        end
    else
        :(write(io, p.$field))
    end
end
@noinline function read_particle(io::IO, h::IAEAHeader{Nf, Ni}) where {Nf, Ni}

    P = ptype(h)
    typ8 = read_(io, Int8)
    typ = ParticleType(abs(typ8))
    E = read_(io, Float32)
    new_history = E < 0
    E = abs(E)
    x = read_default(io, Val(:x), h)
    y = read_default(io, Val(:y), h)
    z = read_default(io, Val(:z), h)
    u = read_default(io, Val(:u), h)
    v = read_default(io, Val(:v), h)
    weight = read_default(io, Val(:weight), h)

    sign_w = Float32(-1)^(typ8 < 0)
    u,v,w = compute_u_v_w(u,v,sign_w)

    extra_floats = read_(io, NTuple{Nf, Float32})
    extra_ints = read_(io, NTuple{Ni, Int32})
    P(typ,
        E,weight,
        x,y,z,
        u,v,w,
        new_history,
        extra_floats, extra_ints,
    )
end

@noinline function write_particle(io::IO,
                                  p::IAEAParticle{Nf, Ni},
                                  h::IAEAHeader{Nf, Ni}) where {Nf, Ni}
    typ8 = Int8(p.typ)
    sign_typ8 = Int8(-1)^(p.w < 0)
    typ8 = sign_typ8 * typ8
    sign_E = Float32(-1)^(p.new_history)
    E = sign_E * p.E
    ret = 0
    ret += write(io, typ8, E)
    ret += write_default(io, Val(:x), p, h)
    ret += write_default(io, Val(:y), p, h)
    ret += write_default(io, Val(:z), p, h)
    ret += write_default(io, Val(:u), p, h)
    ret += write_default(io, Val(:v), p, h)
    ret += write_default(io, Val(:weight), p, h)
    for f in p.extra_floats
        ret += write(io, f)
    end
    for i in p.extra_ints
        ret += write(io, i)
    end
    ret
end

ptype(h::Type{IAEAHeader{Nf, Ni, NT}}) where {Nf, Ni, NT} = IAEAParticle{Nf, Ni}


text:
using Ramnet.Encoders: Thermometer, CircularThermometer, encode!, encode

@testset "Simple Thermometer" begin
    ## Thermometer instantiation

    @test_throws ArgumentError Thermometer(0.0, -20.0, 10)
    @test_throws DomainError Thermometer(0.0, 1.0, -2)

    ## Encoding scalars

    thermo = Thermometer(0.0, 1.0, 10)
    out_pattern = Vector{Bool}(undef, 10)

    encode!(thermo, 0.5, out_pattern)

    @test out_pattern == Bool[1,1,1,1,1,0,0,0,0,0]
    @test_throws DimensionMismatch encode!(thermo, 0.5, Vector{Bool}(undef, 5))

    @test encode(thermo, 0.5) == Bool[1,1,1,1,1,0,0,0,0,0]
    @test encode(thermo, -1.1) == zeros(Bool, 10)
    @test encode(thermo, 2.0) == ones(Bool, 10)

    ## Encoding vectors

    out_pattern = Matrix{Bool}(undef, 10, 2)
    encode!(thermo, [0.0, 0.5], out_pattern)

    @test out_pattern == hcat(
      zeros(Bool, 10),
      Bool[1,1,1,1,1,0,0,0,0,0]
    )
    @test_throws DimensionMismatch encode!(thermo, [0.0, 0.5], Matrix{Bool}(undef, 5, 2))
    @test_throws DimensionMismatch encode!(thermo, [0.0, 0.5], Matrix{Bool}(undef, 10, 1))

    @test encode(thermo, [0.0, 0.5]; flat=false) == hcat(
      zeros(Bool, 10),
      Bool[1,1,1,1,1,0,0,0,0,0]
    )
    @test encode(thermo, [-1.1, 0.5]; flat=false) == hcat(
      zeros(Bool, 10),
      Bool[1,1,1,1,1,0,0,0,0,0]
    )
    @test encode(thermo, [2.0, 0.5]; flat=false) == hcat(
      ones(Bool, 10),
      Bool[1,1,1,1,1,0,0,0,0,0]
    )

    ### Flat outputs
    out_pattern = Vector{Bool}(undef, 20)
    encode!(thermo, [0.0, 0.5], out_pattern)

    @test out_pattern == vcat(
      zeros(Bool, 10),
      Bool[1,1,1,1,1,0,0,0,0,0]
    )

    @test_throws DimensionMismatch encode!(thermo, [0.0, 0.5], Vector{Bool}(undef, 10))

    @test encode(thermo, [0.0, 0.5]) == vcat(
      zeros(Bool, 10),
      Bool[1,1,1,1,1,0,0,0,0,0]
    )
    @test encode(thermo, [-1.1, 0.5]) == vcat(
      zeros(Bool, 10),
      Bool[1,1,1,1,1,0,0,0,0,0]
    )
    @test encode(thermo, [2.0, 0.5]) == vcat(
      ones(Bool, 10),
      Bool[1,1,1,1,1,0,0,0,0,0]
    )

    ## Encoding matrices (one observation per row)

    out_pattern = Array{Bool,3}(undef, 10, 2, 2)
    encode!(thermo, [0.0 0.5; 0.7 1.0], out_pattern)

    @test out_pattern == cat(
        hcat(
          zeros(Bool, 10),
          Bool[1,1,1,1,1,0,0,0,0,0]
        ),
        hcat(
          Bool[1,1,1,1,1,1,1,0,0,0],
          ones(Bool, 10)
        );
        dims=3
    )
    @test_throws DimensionMismatch encode!(thermo, [0.0 0.5; 0.7 1.0], Array{Bool,3}(undef, 5, 2, 2))
    @test_throws DimensionMismatch encode!(thermo, [0.0 0.5; 0.7 1.0], Array{Bool,3}(undef, 10, 1, 2))
    @test_throws DimensionMismatch encode!(thermo, [0.0 0.5; 0.7 1.0], Array{Bool,3}(undef, 10, 2, 1))

    @test encode(thermo, [-1.0 0.5; 0.7 1.0]; flat=false) == cat(
        hcat(
          zeros(Bool, 10),
          Bool[1,1,1,1,1,0,0,0,0,0]
        ),
        hcat(
          Bool[1,1,1,1,1,1,1,0,0,0],
          ones(Bool, 10)
        );
        dims=3
    )

    ### Flat outputs

    out_pattern = Matrix{Bool}(undef, 2, 20)
    encode!(thermo, [0.0 0.5; 0.7 1.0], out_pattern)

    @test out_pattern ==  vcat(
        reshape(vcat(
          zeros(Bool, 10),
          Bool[1,1,1,1,1,0,0,0,0,0]
        ), (1, :)),
        reshape(vcat(
          Bool[1,1,1,1,1,1,1,0,0,0],
          ones(Bool, 10)
        ), (1, :))
    )
    @test_throws DimensionMismatch encode!(thermo, [0.0 0.5; 0.7 1.0], Matrix{Bool}(undef, 1, 20))
    @test_throws DimensionMismatch encode!(thermo, [0.0 0.5; 0.7 1.0], Matrix{Bool}(undef, 2, 10))

    @test encode(thermo, [-1.0 0.5; 0.7 1.0]) == vcat(
        reshape(vcat(
          zeros(Bool, 10),
          Bool[1,1,1,1,1,0,0,0,0,0]
        ), (1, :)),
        reshape(vcat(
          Bool[1,1,1,1,1,1,1,0,0,0],
          ones(Bool, 10)
        ), (1, :))
    )

    # Specifying different limits for the components of the input
    thermo = Thermometer([0.0, 1.0], [1.0, 10.0], 10)

    # When multiple limits are specified, input is expected to be at least 1D
    @test_throws DomainError encode(thermo, 0.5)

    @test encode(thermo, [0.5, 5]) == vcat(Bool[1,1,1,1,1,0,0,0,0,0], Bool[1,1,1,1,1,0,0,0,0,0])
    @test encode(thermo, [0.5, 5]; flat=false) == hcat(Bool[1,1,1,1,1,0,0,0,0,0], Bool[1,1,1,1,1,0,0,0,0,0])
end

@testset "Circular Thermometer" begin
    @test_throws ArgumentError CircularThermometer(10, 1, 5)
    @test_throws DomainError CircularThermometer(0, 1, 1)

    thermo = CircularThermometer(Float64, 0, 1, 5)

    @test encode(thermo, 0.0) == Bool[1, 1, 0, 0, 0]
    @test encode(thermo, 0.2) == Bool[0, 1, 1, 0, 0]
    @test encode(thermo, 0.4) == Bool[0, 0, 1, 1, 0]
    @test encode(thermo, [0.2, 0.4]) == Bool[0, 1, 1, 0, 0, 0, 0, 1, 1, 0]
    @test encode(thermo, [0.2, 0.4];flat=false) == Bool[
      0 0
      1 0
      1 1
      0 1
      0 0
    ]

    @test encode(thermo, [
      0.0 0.2
      0.4 0.6
    ]) == Bool[
      1 1 0 0 0 0 1 1 0 0
      0 0 1 1 0 0 0 0 1 1
    ]

    @test encode(thermo, [
      0.0 0.2
      0.4 0.6
    ]; flat=false) == cat(
      [
        1 0
        1 1
        0 1
        0 0
        0 0
      ],
      [
        0 0
        0 0
        1 0
        1 1
        0 1
      ];
      dims=3
    )

    # Specifying different limits for the components of the input
    thermo = CircularThermometer([0.0, 1.0], [1.0, 10.0], 5)

    @test encode(thermo, [0.9, 7.0]) == Bool[1, 0, 0, 0, 1, 0, 0, 0, 1, 1]
end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
#: by David Thierry 2021
using JuMP
using Clp
using DataFrames
using StatsPlots
using CSV
tHorz = 30
m = Model()

#vCapCombTurb = 3.
vCapSteamTurb = 2.
vCapTransInter = 5.
vCapPcc = 20.
vCapComp = 1000.

rnum = rand(Float64, tHorz)

#@variable(m, 0 <= vCapCombTurb)
#@variable(m, 0 <= vPowCombTurb[0:tHorz - 1] <= vCapCombTurb)
@variable(m, 0 <= vPowCombTurb[0:tHorz - 1])

#@variable(m, 0 <= ePowSteamTurb[0:tHorz - 1] <= vCapSteamTurb)
@variable(m, 0 <= ePowSteamTurb[0:tHorz - 1])
#@variable(m, 0 <= vCapSteamTurb)

#@variable(m, rnum[i+1] <= ePowOut[i = 0:tHorz - 1] <= vCapTransInter)
@variable(m, ePowOut[i = 0:tHorz - 1] >= rnum[i + 1])
#@variable(m, 0 <= ePowOut[i = 0:tHorz - 1])
@variable(m, 0 <= ePowHp[0:tHorz - 1])
@variable(m, 0 <= ePowIp[0:tHorz - 1])
#@variable(m, 0 <= vCapTransInter)

@variable(m, 0 <= eFuel[0:tHorz - 1])

@variable(m, 0 <= eSteamIpLp[0:tHorz - 1])
@variable(m, 0 <= eSteamLp1[0:tHorz - 1])

@variable(m, 0 <= ePowLp1[0:tHorz - 1])
@variable(m, 0 <= ePowLp2[0:tHorz - 1])

@variable(m, 0 <= eCo2Fuel[0:tHorz - 1])

# Pcc
#@variable(m, 0 <= eCo2CapPcc[0:tHorz - 1] <= vCapPcc)
@variable(m, 0 <= eCo2CapPcc[0:tHorz - 1])
@variable(m, 0 <= eSteamUsePcc[0:tHorz - 1])
@variable(m, 0 <= ePowUsePcc[0:tHorz - 1])
@variable(m, 0 <= eCo2PccOut[0:tHorz - 1])
#@variable(m, 0 <= vCo2PccVent[0:tHorz - 1] <= 0.1)
@variable(m, 0 <= vCo2PccVent[0:tHorz - 1])
#vCo2PccVent = 0.0
@variable(m, 0 <= eCo2DacFlueIn[0:tHorz - 1])

# Dac
@variable(m, 0 <= a0[0:tHorz - 1])  # We could generalize this for any number of hours.
@variable(m, 0 <= a1[0:tHorz])
#@variable(m, 0 <= a2[0:tHorz - 1])

@variable(m, 0 <= aR0[0:tHorz - 1])
#@variable(m, 0 <= aR1[0:tHorz - 1])

@variable(m, 0 <= sF[0:tHorz])
@variable(m, 0 <= sS[0:tHorz])

@variable(m, 0 <= vAbs[0:tHorz - 1])
@variable(m, 0 <= vReg[0:tHorz - 1])

@variable(m, 0 <= eCo2StoreDac[0:tHorz - 1])
@variable(m, 0 <= eCo2CapDac[0:tHorz - 1])
@variable(m, eSteamUseDac[0:tHorz - 1] >= 0)
@variable(m, 0 <= ePowUseDac[0:tHorz - 1])
# CO2 compression

@variable(m, 0 <= eCo2DacVent[0:tHorz - 1])

@variable(m, 0 <= eCo2Comp[0:tHorz - 1])
#@variable(m, 0 <= ePowUseComp[0:tHorz - 1] <= vCapComp)
@variable(m, 0 <= ePowUseComp[0:tHorz - 1])
#@variable(m, 0 <= vCapComp)
@variable(m, 0 <= eCo2Vent[0:tHorz - 1])

# Parameters
pHeatRateCombTur = 1.
# pSteamRate = 1.
pEmissFactor = 20.1

# pCo2CapRatePcc = 1.

#: Turbine Parameters
pFuelPowHp = 1.
pFuelPowIp = 0.4
pFuelSteamIpLp = 1.
pHeatRateLp1 = 11.11
pFuelPowLp = 0.09
#
#
pSteamUseRatePcc = 0.001
pPowUseRatePcc = 0.0022

pSteamUseRateDacAir = 0.00001
pSteamUseRateDacFlue = 0.00001

pPowUseRateDacAir = 0.0001
pPowUseRateDacFlue = 0.00002

pHeatRateSteamTur = 1.

pCo2PccCapRate = .1
pSorbCo2Cap = 0.1

pPowUseRateDac = 0.0001

pPowUseRateComp = 0.001
#

pCostInvCombTurb = 1e+02
pCostInvSteaTurb = 1e+02
pCostInvTransInter = 1e+02
pCostInvPcc = 1e+02
pCostInvDac = 1e+03
pCostInvComp = 1e+01

# Cost parameters.
pCostFuel = 1e+02
pEmissionPrice = 1e+04
pCo2TranspPrice = 1e+01
pPowBasePrice = 1e+02

# Constraints
# 1
#@constraint(m, capCombTurIn[i = 0:tHorz - 1], vPowCombTurb[i] <= vCapCombTurb)
# 2
# @constraint(m, capSteamTurbIn[i = 0:tHorz - 1], ePowSteamTurb[i] <= vCapSteamTurb)
# 3
# @constraint(m, capTransmIntercIn[i = 0:tHorz - 1], ePowOut[i] <= vCapTransInter)
# 4
@constraint(m, fuelEq[i = 0:tHorz - 1], eFuel[i] == pHeatRateCombTur * vPowCombTurb[i])
# 5
@constraint(m, co2FuelEq[i = 0:tHorz - 1], eCo2Fuel[i] == pEmissFactor * vPowCombTurb[i])
@variable(m, 0 <= ePowGross[0:tHorz - 1])
@constraint(m, powGrossEq[i = 0:tHorz - 1], ePowGross[i] == vPowCombTurb[i] + ePowSteamTurb[i])
# 6
@constraint(m, powOutEq[i = 0:tHorz - 1], ePowOut[i] ==  ePowGross[i] - ePowUsePcc[i] - ePowUseDac[i] - ePowUseComp[i])

# Steam
# 7a
@constraint(m, powHpEq[i = 0:tHorz - 1], ePowHp[i] == pFuelPowHp * eFuel[i])
# 7b
@constraint(m, powIpEq[i = 0:tHorz - 1], ePowIp[i] == pFuelPowIp * eFuel[i])

# 8
@constraint(m, stamIpLpEq[i = 0:tHorz - 1], eSteamIpLp[i] == pFuelSteamIpLp * eFuel[i])
# 9
@constraint(m, steamLp1[i = 0:tHorz - 1], eSteamLp1[i] == eSteamIpLp[i] - eSteamUsePcc[i] - eSteamUseDac[i])
# 10
@constraint(m, powLp1Eq[i = 0:tHorz - 1], ePowLp1[i] == eSteamLp1[i] / pHeatRateLp1)
# 11
@constraint(m, powLp2Eq[i = 0:tHorz - 1], ePowLp2[i] == pFuelPowLp * eFuel[i])
# 12
@constraint(m, powerSteamEq[i = 0:tHorz - 1], ePowSteamTurb[i] == ePowHp[i] + ePowIp[i] + ePowLp1[i] + ePowLp2[i])

# PCC
# 13
#@constraint(m, co2CapPccEq[i = 0:tHorz - 1], eCo2CapPcc[i] == pCo2PccCapRate * eCo2Fuel[i])
@constraint(m, co2CapPccEq[i = 0:tHorz - 1], eCo2CapPcc[i] == 0.95 * eCo2Fuel[i])
# 14
@constraint(m, co2PccOutEq[i = 0:tHorz - 1], eCo2PccOut[i] == eCo2Fuel[i] - eCo2CapPcc[i])
# 15
@constraint(m, co2DacFlueInEq[i = 0:tHorz - 1], eCo2DacFlueIn[i] == eCo2PccOut[i] - vCo2PccVent[i])
# 16
# @constraint(m, co2CapPccIn[i = 0:tHorz - 1], eCo2CapPcc[i] <= vCapPcc)
# 17
@constraint(m, steamUsePccEq[i = 0:tHorz - 1], eSteamUsePcc[i] == pSteamUseRatePcc * eCo2CapPcc[i])
# 18
@constraint(m, powerUsePccEq[i = 0:tHorz - 1], ePowUsePcc[i] == pPowUseRatePcc * eCo2CapPcc[i])

# DAC
@constraint(m, a0Eq[i = 0:tHorz - 1], a0[i] == vAbs[i])
@constraint(m, aR0Eq[i = 0:tHorz - 1], aR0[i] == vReg[i])

@constraint(m, a1dEq[i = 0:tHorz - 1], a1[i + 1] == a0[i])
@constraint(m, storeFeq[i = 0:tHorz - 1], sF[i + 1] == sF[i] - vAbs[i] + aR0[i])
@constraint(m, storeSeq[i = 0:tHorz - 1], sS[i + 1] == sS[i] - vReg[i] + a1[i])
@constraint(m, capDacEq, sF[0] == 100 / pSorbCo2Cap)
# 22
@constraint(m, icA1Eq, a1[0] == 0.)
# @constraint(m, icA2Eq, aR1[0] == 0.)
#@constraint(m, icSfEq, sF[0] == 100.)
@constraint(m, icSsEq, sS[0] == 0.)
# 23
@constraint(m, co2StoreDacEq[i = 0:tHorz - 1], eCo2StoreDac[i] == pSorbCo2Cap * sS[i])
# 24
@constraint(m, co2CapDacEq[i = 0:tHorz - 1], eCo2CapDac[i] == pSorbCo2Cap * aR0[i])
# 25
@constraint(m, steamUseDacEq[i = 0:tHorz - 1], eSteamUseDac[i] == pSteamUseRateDacFlue * eCo2CapDac[i])
# 26
@constraint(m, powUseDacEq[i = 0:tHorz - 1], ePowUseDac[i] == pPowUseRateDacFlue * eCo2CapDac[i])
# equal to the amount vented at least in flue mode.
@constraint(m, co2DacFlueVentEq[i = 0:tHorz - 1], eCo2DacVent[i] == eCo2DacFlueIn[i] - eCo2CapDac[i])
#
# Co2 Compression
# 27
@constraint(m, co2CompEq[i = 0:tHorz - 1], eCo2Comp[i] == eCo2CapPcc[i])
# 28
@constraint(m, powUseCompEq[i = 0:tHorz - 1], ePowUseComp[i] == pPowUseRateComp * eCo2Comp[i])
# 29
# @constraint(m, powUseCompIn[i = 0:tHorz - 1], ePowUseComp[i] <= vCapComp)

@constraint(m, co2VentEq[i = 0:tHorz - 1], eCo2Vent[i] == vCo2PccVent[i] + eCo2DacVent[i])

@expression(m, eObjfExpr, sum(pCostFuel * eFuel[i] +
                              pEmissionPrice * eCo2Vent[i] +
                              pCo2TranspPrice * eCo2Comp[i] -
                              pPowBasePrice * ePowOut[i] for i in 0:tHorz - 1))
#@expression(m, eObjfExpr, sum(-ePowOut[i] for i in 0:tHorz - 1))

@objective(m, Min, eObjfExpr)

set_optimizer(m, Clp.Optimizer)
set_optimizer_attribute(m, "LogLevel", 3)
set_optimizer_attribute(m, "PresolveType", 1)

optimize!(m)
println(termination_status(m))

#f = open("model.lp", "w")
#print(f, m)
#close(f)

write_to_file(m, "lp_mk0.mps")
#format::MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC

# Design decisions.
# vCapCombTurb
# vCapSteamTurb
# vCapTransInter
# vCapPcc
# vCapDac
# vCapComp
#
# Raw materials.
# eFuel
# sF0
# sS0
#
# penalties
# vCo2Vent
#
#
## vCo2PccVent * 2N

# Actual variables
# vPowCombTurb * 2

# Co2 Data Frame
df_co = DataFrame(Symbol("Co2Fuel") => Float64[], # Pairs.
                  Symbol("Co2CapPcc") => Float64[],
                  Symbol("Co2PccOut") => Float64[],
                  Symbol("vCo2PccVent") => Float64[],
                  Symbol("Co2DacFlueIn") => Float64[],
                  Symbol("Co2StoreDac") => Float64[],
                  Symbol("Co2CapDac") => Float64[],
                  Symbol("Co2DacVent") => Float64[],
                  Symbol("Co2Vent") => Float64[],
                 )
for i in 0:tHorz-1
    push!(df_co, (value(eCo2Fuel[i]),
               value(eCo2CapPcc[i]),
               value(eCo2PccOut[i]),
               value(vCo2PccVent[i]),
               value(eCo2DacFlueIn[i]),
               value(eCo2StoreDac[i]),
               value(eCo2CapDac[i]),
               value(eCo2DacVent[i]),
               value(eCo2Vent[i])))
end

print(df_co)

df_pow = DataFrame(Symbol("vPowCombTurb") => Float64[], # Pairs.
                  Symbol("PowSteamTurb") => Float64[],
                  Symbol("PowHp") => Float64[],
                  Symbol("PowIp") => Float64[],
                  Symbol("PowLp1") => Float64[],
                  Symbol("PowLp2") => Float64[],
                  Symbol("PowUsePcc") => Float64[],
                  Symbol("PowUseDac") => Float64[],
                  Symbol("PowUseComp") => Float64[],
                  Symbol("PowGross") => Float64[],
                  Symbol("PowOut") => Float64[],
                  Symbol("Demand") => Float64[]
                 )
for i in 0:tHorz-1
    push!(df_pow, (
                   value(vPowCombTurb[i]),
                   value(ePowSteamTurb[i]),
                   value(ePowHp[i]),
                   value(ePowIp[i]),
                   value(ePowLp1[i]),
                   value(ePowLp2[i]),
                   value(ePowUsePcc[i]),
                   value(ePowUseDac[i]),
                   value(ePowUseComp[i]),
                   value(ePowGross[i]),
                   value(ePowOut[i]),
                   rnum[i + 1]
                  ))
end

df_steam = DataFrame(Symbol("SteamIpLp") => Float64[], # Pairs.
                     Symbol("SteamLp1") => Float64[],
                     Symbol("SteamUsePcc") => Float64[],
                     Symbol("SteamUseDac") => Float64[],
                     Symbol("Fuel") => Float64[]
                    )
for i in 0:tHorz-1
    push!(df_steam, (
                     value(eSteamIpLp[i]),
                     value(eSteamLp1[i]),
                     value(eSteamUsePcc[i]),
                     value(eSteamUseDac[i]),
                     value(eFuel[i])
                    ),
         )
end
# Generate Dataframe of the DAC variables.
df_dac = DataFrame(:sF => Float64[],
                     :sS => Float64[],
                     :vAbs => Float64[],
                     :vReg => Float64[],
                     :a0 => Float64[],
                     :a1 => Float64[],
                     :aR0 => Float64[])
for i in 0:tHorz - 1
    push!(df_dac,(value(sF[i]), value(sS[i]),
                  value(vAbs[i]), value(vReg[i]),
                  value(a0[i]), value(a1[i]),
                  value(aR0[i])))
end

print(df_co)
print(df_pow)
print(df_steam)
print(df_dac)

CSV.write("df_co.csv", df_co)
CSV.write("df_pow.csv", df_pow)
CSV.write("df_steam.csv", df_steam)
CSV.write("df_dac.csv", df_dac)

max_stars_repo_name: dthierry/princetonDacLti

max_issues_repo_name: dthierry/princetonDacLti

max_forks_repo_name: dthierry/princetonDacLti


text:
using MLGoodies
using Test

@testset "MLGoodies.jl" begin
    include("utils/utils.jl")
    include("transforms/transforms.jl")
end

max_stars_repo_name: dhonza/MLGoodies


text:
export AbstractLearner

using Flux

"""
    (learner::AbstractLearner)(env)

A learner is usually used to estimate state values, state-action values or distributional values based on experiences.
"""
abstract type AbstractLearner end

function (learner::AbstractLearner)(env) end

"""
    get_priority(p::AbstractLearner, experience)
"""
function RLBase.get_priority(p::AbstractLearner, experience) end

# TODO: deprecate this default function
Flux.testmode!(learner::AbstractLearner, mode = true) =
    Flux.testmode!(learner.approximator, mode)


text:
# These constants were manually copied from the file mpi.h in the Microsoft
# MPI SDK v7

const MPI_BYTE = Int32(0x4c00010d)
const MPI_WCHAR = Int32(0x4c00020e)
const MPI_INT8_T = Int32(0x4c000133)
const MPI_UINT8_T = Int32(0x4c000137)
const MPI_INT16_T = Int32(0x4c000234)
const MPI_UINT16_T = Int32(0x4c000238)
const MPI_INT32_T = Int32(0x4c000435)
const MPI_UINT32_T = Int32(0x4c000439)
const MPI_INT64_T = Int32(0x4c000836)
const MPI_UINT64_T = Int32(0x4c00083a)
const MPI_REAL4 = Int32(0x4c000427)
const MPI_REAL8 = Int32(0x4c000829)
const MPI_COMPLEX8 = Int32(0x4c000828)
const MPI_COMPLEX16 = Int32(0x4c00102a)
const MPI_INTEGER1 = Int32(0x4c00012d)
const MPI_INTEGER2 = Int32(0x4c00022f)
const MPI_INTEGER4 = Int32(0x4c000430)
const MPI_INTEGER8 = Int32(0x4c000831)
const MPI_COMM_NULL = Int32(67108864)
const MPI_COMM_SELF = Int32(1140850689)
const MPI_COMM_WORLD = Int32(1140850688)
const MPI_COMM_TYPE_SHARED = Int32(1)
const MPI_OP_NULL = Int32(402653184)
const MPI_BAND = Int32(1476395014)
const MPI_BOR = Int32(1476395016)
const MPI_BXOR = Int32(1476395018)
const MPI_LAND = Int32(1476395013)
const MPI_LOCK_EXCLUSIVE = Int32(234)
const MPI_LOCK_SHARED = Int32(235)
const MPI_LOR = Int32(1476395015)
const MPI_LXOR = Int32(1476395017)
const MPI_MAX = Int32(1476395009)
const MPI_MAXLOC = Int32(1476395020)
const MPI_MIN = Int32(1476395010)
const MPI_MINLOC = Int32(1476395019)
const MPI_NO_OP = Int32(0x5800000e)
const MPI_PROD = Int32(1476395012)
const MPI_REPLACE = Int32(1476395021)
const MPI_SUM = Int32(1476395011)
const MPI_REQUEST_NULL = Int32(738197504)
const MPI_INFO_NULL = Int32(469762048)
const MPI_STATUS_SIZE = Int32(5)
const MPI_ERROR = Int32(5)
const MPI_SOURCE = Int32(3)
const MPI_TAG = Int32(4)
const MPI_ANY_SOURCE = Int32(-2)
const MPI_ANY_TAG = Int32(-1)
const MPI_TAG_UB = Int32(1681915906)
const MPI_UNDEFINED = Int32(-32766)
const HAVE_MPI_COMM_C2F = false

const libmpi = "msmpi.dll"

const MPI_ABORT = (:MPI_ABORT, libmpi)
const MPI_ACCUMULATE = (:MPI_ACCUMULATE, libmpi)
const MPI_ALLREDUCE = (:MPI_ALLREDUCE, libmpi)
const MPI_INIT = (:MPI_INIT, libmpi)
const MPI_CANCEL = (:MPI_CANCEL, libmpi)
const MPI_COMM_FREE = (:MPI_COMM_FREE, libmpi)
const MPI_COMM_GET_PARENT = (:MPI_COMM_GET_PARENT, libmpi)
const MPI_COMM_RANK = (:MPI_COMM_RANK, libmpi)
const MPI_COMM_SIZE = (:MPI_COMM_SIZE, libmpi)
const MPI_COMM_SPLIT = (:MPI_COMM_SPLIT, libmpi)
const MPI_COMM_SPLIT_TYPE = (:MPI_COMM_SPLIT_TYPE, libmpi)
const MPI_BARRIER = (:MPI_BARRIER, libmpi)
const MPI_FINALIZE = (:MPI_FINALIZE, libmpi)
const MPI_BCAST = (:MPI_BCAST, libmpi)
const MPI_REDUCE = (:MPI_REDUCE, libmpi)
const MPI_IRECV = (:MPI_IRECV, libmpi)
const MPI_RECV = (:MPI_RECV, libmpi)
const MPI_ISEND = (:MPI_ISEND, libmpi)
const MPI_WAITALL = (:MPI_WAITALL, libmpi)
const MPI_ALLGATHER = (:MPI_ALLGATHER, libmpi)
const MPI_ALLGATHERV = (:MPI_ALLGATHERV, libmpi)
const MPI_ALLTOALL = (:MPI_ALLTOALL, libmpi)
const MPI_ALLTOALLV = (:MPI_ALLTOALLV, libmpi)
const MPI_INFO_CREATE = (:MPI_INFO_CREATE, libmpi)
const MPI_INFO_DELETE = (:MPI_INFO_DELETE, libmpi)
const MPI_INFO_FREE = (:MPI_INFO_FREE, libmpi)
const MPI_INFO_GET = (:MPI_INFO_GET, libmpi)
const MPI_INFO_GET_VALUELEN = (:MPI_INFO_GET_VALUELEN, libmpi)
const MPI_INFO_SET = (:MPI_INFO_SET, libmpi)
const MPI_INITIALIZED = (:MPI_INITIALIZED, libmpi)
const MPI_FINALIZED = (:MPI_FINALIZED, libmpi)
const MPI_FETCH_AND_OP = (:MPI_FETCH_AND_OP, libmpi)
const MPI_OP_CREATE = (:MPI_OP_CREATE, libmpi)
const MPI_OP_FREE = (:MPI_OP_FREE, libmpi)
const MPI_SCATTER = (:MPI_SCATTER, libmpi)
const MPI_SCATTERV = (:MPI_SCATTERV, libmpi)
const MPI_SEND = (:MPI_SEND, libmpi)
const MPI_SCAN = (:MPI_SCAN, libmpi)
const MPI_EXSCAN = (:MPI_EXSCAN, libmpi)
const MPI_GET = (:MPI_GET, libmpi)
const MPI_GET_ACCUMULATE = (:MPI_GET_ACCUMULATE, libmpi)
const MPI_GET_ADDRESS = (:MPI_GET_ADDRESS, libmpi)
const MPI_GET_COUNT = (:MPI_GET_COUNT, libmpi)
const MPI_GATHER = (:MPI_GATHER, libmpi)
const MPI_GATHERV = (:MPI_GATHERV, libmpi)
const MPI_COMM_DUP = (:MPI_COMM_DUP, libmpi)
const MPI_INTERCOMM_MERGE = (:MPI_INTERCOMM_MERGE, libmpi)
const MPI_IPROBE = (:MPI_IPROBE, libmpi)
const MPI_PROBE = (:MPI_PROBE, libmpi)
const MPI_PUT = (:MPI_PUT, libmpi)
const MPI_TEST = (:MPI_TEST, libmpi)
const MPI_TESTSOME = (:MPI_TESTSOME, libmpi)
const MPI_TESTANY = (:MPI_TESTANY, libmpi)
const MPI_TESTALL = (:MPI_TESTALL, libmpi)
const MPI_TYPE_CREATE_STRUCT = (:MPI_TYPE_CREATE_STRUCT, libmpi)
const MPI_TYPE_COMMIT = (:MPI_TYPE_COMMIT, libmpi)
const MPI_WAIT = (:MPI_WAIT, libmpi)
const MPI_WAITANY = (:MPI_WAITANY, libmpi)
const MPI_WAITSOME = (:MPI_WAITSOME, libmpi)
const MPI_WIN_ATTACH = (:MPI_WIN_ATTACH, libmpi)
const MPI_WIN_CREATE = (:MPI_WIN_CREATE, libmpi)
const MPI_WIN_CREATE_DYNAMIC = (:MPI_WIN_CREATE, libmpi)
const MPI_WIN_DETACH = (:MPI_WIN_DETACH, libmpi)
const MPI_WIN_FENCE = (:MPI_WIN_FENCE, libmpi)
const MPI_WIN_FLUSH = (:MPI_WIN_FLUSH, libmpi)
const MPI_WIN_FREE = (:MPI_WIN_FREE, libmpi)
const MPI_WIN_LOCK = (:MPI_WIN_LOCK, libmpi)
const MPI_WIN_SYNC = (:MPI_WIN_SYNC, libmpi)
const MPI_WIN_UNLOCK = (:MPI_WIN_UNLOCK, libmpi)
const MPI_WTICK = (:MPI_WTICK, libmpi)
const MPI_WTIME = (:MPI_WTIME, libmpi)

primitive type CComm 32 end
primitive type CInfo 32 end
primitive type CWin 32 end


text:
"""
    mveltsplot(setting; alpha = 0.05, showplot = true)

Generate MVE - LTS plot for visual detecting of regression outliers.

# Arguments
- `setting::RegressionSetting`: A regression setting object.
- `alpha::Float64`: Probability for quantiles of Chi-Squared statistic.
- `showplot::Bool`: Whether a plot is shown or only return statistics.

# References
Van Aelst, Stefan, and Peter Rousseeuw. "Minimum volume ellipsoid." Wiley
Interdisciplinary Reviews: Computational Statistics 1.1 (2009): 71-82.
"""
function mveltsplot(setting::RegressionSetting; alpha=0.05, showplot=true)
    ltsresult = lts(setting)
    mveresult = mve(setting.data)

    n, p = size(setting.data)
    indices = collect(1:n)

    chidist = Chisq(p)
    chicrit = sqrt(quantile(chidist, 1.0 - alpha))

    scaledresiduals = ltsresult["scaled.residuals"]
    robdistances = sqrt.(mveresult["squared.mahalanobis"])

    regularpoints = filter(i -> abs(scaledresiduals[i]) < 2.5 && robdistances[i] < chicrit, indices)
    outliers = filter(i -> scaledresiduals[i] > 2.5 && robdistances[i] < chicrit, indices)
    leverage = filter(i -> abs(scaledresiduals[i]) < 2.5 && robdistances[i] > chicrit, indices)
    outlierandleveragepoints = filter(i -> abs(scaledresiduals[i]) > 2.5 && robdistances[i] > chicrit, indices)

    scplot = nothing
    if showplot
        scplot = scatter(robdistances,
            scaledresiduals,
            legend=false,
            series_annotations=text.(1:n, :bottom),
            tickfont=font(10), guidefont=font(10), labelfont=font(10)
            )
        title!("MVE - LTS Plot")
        xlabel!("Robust distances")
        ylabel!("Scaled residuals")
        hline!([-2.5, 2.5])
        vline!([chicrit])
    end

    result = Dict()
    result["plot"] = scplot
    result["robust.distances"] = robdistances
    result["scaled.residuals"] = scaledresiduals
    result["chi.squared"] = chicrit^2
    result["regular.points"] = regularpoints
    result["outlier.points"] = outliers
    result["leverage.points"] = leverage
    result["outlier.and.leverage.points"] = outlierandleveragepoints
    return result
end

max_stars_repo_name: akadal/LinRegOutliers

max_issues_repo_name: akadal/LinRegOutliers

max_forks_repo_name: akadal/LinRegOutliers


text:
module Matrix2Tabular

export float2sci, matrix2tabular

function float2sci(x::Float64, figs::Int64=2, dollars::Bool=false)
    power = convert(Int,round(log10(x),RoundDown))
    if 0 <= power <= 1
        figs = max(figs-power,2)
        x = string(round(x,digits=figs))
    else
        x = round(x/(10.0^power),digits=figs)
        x = string(x,"\\times10^{",power,"}")
    end

    if dollars
        return string(" \$ ",x," \$ ")
    else
        return x
    end
end

function matrix2tabular(A::Matrix{Float64}, figs::Int64=2, dollars::Bool=false)
    M,N = size(A)

    outstring = " "

    for m = 1:M
        for n = 1:N
            outstring *= float2sci(A[m,n], figs, dollars)
            if n == N
                outstring *= "\\\\ \n "
            end
        end
    end

    return outstring
end
end


text:
#!/usr/bin/env julia

# A trivial example of deferred canvases.

require("Compose")
using Compose

function sierpinski(n)
    if n == 0
        canvas() << polygon((1,1), (0,1), (1/2, 0))
    else
        t = sierpinski(n - 1)
        compose(canvas(),
                (canvas(1/4,   0, 1/2, 1/2), t),
                (canvas(  0, 1/2, 1/2, 1/2), t),
                (canvas(1/2, 1/2, 1/2, 1/2), t))
    end
end

img = SVG("sierpinski.svg", 4inch, 4(sqrt(3)/2)inch)
draw(img, compose(canvas(), fill(nothing), linewidth(0.1mm),
                  deferredcanvas((box, unit) -> sierpinski(8))))
finish(img)


text:
#+t           L:foo
# foo
# [[foo][]]
# [[foo2][]]
# [[faa][org]]
# "[[foo][]]"
foo()=0

#+t           L:foo2
foo2()=0


text:
using Test
using Pkg;
Pkg.activate(".")

using ArgParse
using Profile
using dlgo
using dlgo.Agent

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--board_size", "-b"
            help = "The game will have dimensions BOARD_SIZExBOARD_SIZE"
            arg_type = Int
            default = 9
        "--num_games", "-n"
            help = "Number of games to play"
            arg_type = Int
            default = 1
        "--komi", "-k"
            help = "White player handicap for scoring"
            arg_type = Float64
            default = 7.5
        "--verbose", "-v"
            help = "print results of game to command line"
            arg_type = Bool
            default = true
        "--profile", "-p"
            help = "profile shows time spent per call"
            arg_type = Bool
            default = false
    end
    return parse_args(s, as_symbols=true)
end

function main(board_size::Int, num_games::Int, komi::Float64, verbose::Bool)
    verbose && println("Playing ", num_games, " game(s) with board size ", board_size, "x", board_size)
    for it in 1:num_games
        game = new_game(board_size)
        bot_black = RandomAgent()
        bot_white = RandomAgent()
        nmoves = 0
        while !is_over(game)
            bot = (game.next_player == black) ? bot_black : bot_white
            move = select_move(bot, game)
            game = apply_move(game, move)
            nmoves += 1
        end
        result = compute_game_result(game, komi)
        verbose && println("Finished game ", it, " in " , nmoves, " moves with result ", result)
        verbose && print_board(game.board)
    end
end

args = parse_commandline()
@time begin
    if args[:profile]
        @profile main(args[:board_size], args[:num_games], args[:komi], args[:verbose])
        Profile.print(format=:flat)
    else
        main(args[:board_size], args[:num_games], args[:komi], args[:verbose])
    end
end


text:
using Spectra
using Random
using Test

x = collect(50:1.0:500)

# gaussian peaks
p1 = 20.0 .* exp.(-log(2) .* ((x.-150.0)./15.0).^2)
p2 = 100.0 .* exp.(-log(2) .* ((x.-250.0)./5.0).^2)
p3 = 50.0 .* exp.(-log(2) .* ((x.-450.0)./1.0).^2)
p4 = 20.0 .* exp.(-log(2) .* ((x.-350.0)./30.0).^2)
p5 = 30.0 .* exp.(-log(2) .* ((x.-460.0)./5.0).^2)

# background: a large gaussian + linear
bkg = 60.0 .* exp.(-log(2) .* ((x.-250.0)./200.0).^2) + 0.1.*x

#noise
noise = 2.0 * randn!(ones(size(x,1)))

#observation
y = p1 + p2 + p3 + p4 + p5 + noise +bkg

roi = [0. 100.] # provide whatever matrix

# for the ALS algorithm, 10^2-10^5 lambda and 0.001-0.1 p values are recommended
y_als, bas_als = baseline(x,y,roi,"als",p=0.04,lam=10^6,niter=10)

rmse = sqrt.(sum((y - y_als).^2))./sum(y)

@test rmse <0.05 # error inferior to 5%, should actually be arround 0.026.


text:
include("GeneralMethods.jl")
include("SuperOperator.jl")

#using SymPy
using FileIO
#
# function buildKernels(maxLen, Hcommuting)
#     TrRight=buildTraceMatrixRight(maxLen-1);
#     TrLeft=buildTraceMatrixLeft(maxLen-1);
#     ###For the left link
#     leftCommuter=buildSuperOperatorLeft(kron(speye(2^(maxLen-1)), Hcommuting))
#     MatJoin=sparse(vcat(TrRight, TrLeft, TrLeft*leftCommuter));
#     kernelLeft=nullspace(full(MatJoin))
#     ###For the right link
#     rightCommuter=buildSuperOperatorLeft(kron(Hcommuting, speye(2^(maxLen-1))))
#     MatJoin=sparse(vcat(TrLeft, TrRight, TrRight*rightCommuter));
#     kernelRight=nullspace(full(MatJoin))
#     return kernelLeft, kernelRight
# end

function hermitianed(vect)
    len=getLength(length(vect));
    return vect+view(reshape(vect, (2^len, 2^len))', :)
end

# function makeHermitian(kernel)
#     herm=zeros(Complex128, size(kernel));
#     herm[:, 1]=hermitianed(kernel[:, 1])
#     count=2
#     for j=2:(size(kernel)[2])
#         temp=symmetrize(kernel[:, j])
#         herm[:, count]=copy(temp)
#         if rank(herm)==count
#             count+=1
#         end
#     end
#     println(count)
#     herm[:, count]=hermitianed(1im*kernel[:, 1])
#      count+=1
#      for j=2:(size(kernel)[2])
#          temp=symmetrize(1im*kernel[:, j])
#          herm[:, count]=copy(temp)
#          if rank(herm)==count
#              count+=1
#              if count>size(kernel)[2]
#                  break
#              end
#          end
#      end
#     return herm
# end
#
#
# function generateHermitianKernel(maxLen)
#     sx=[0 1; 1 0]; sz=-[1  0 ; 0 -1]; sy=[0 -1im; 1im 0]
#     kernelLeft, kernelRight=buildKernels(maxLen, sparse(sz))
#     hermitianKernelLeft=makeHermitian(kernelLeft)
#     hermitianKernelRight=makeHermitian(kernelRight)
#     return hermitianKernelLeft, hermitianKernelRight
# end

function makeHermitian(kernel)
    herm=zeros(Complex{Int64}, length(kernel[1]), length(kernel))
    herm[:, 1]=N(hermitianed(kernel[1]))
    count=2
    for j=2:length(kernel)
        println(count)
        herm[:, count]=N(hermitianed(kernel[j]))
        if rank(herm)==count
            count+=1
        end
    end
    herm[:, count]= convert(Array{Complex{Int64},1},  N(hermitianed(1im*kernel[1])))
     count+=1
     for j=2:length(kernel)
         herm[:, count]= convert(Array{Complex{Int64},1},  N(hermitianed(1im*kernel[j])))
         println(count)
         if rank(herm)==count
             count+=1
             if count>length(kernel)
                 break
             end
         end
     end
    return herm
end

function buildKernels(maxLen, Hcommuting)
    TrRight=buildTraceMatrixRight(maxLen-1, typed=SymPy.Sym);
    TrLeft=buildTraceMatrixLeft(maxLen-1, typed=SymPy.Sym);
    ###For the left link
    leftCommuter=buildSuperOperatorLeft(kron(speye(SymPy.Sym, 2^(maxLen-1)), Hcommuting))
    MatJoin=sparse(vcat(TrRight, TrLeft, TrLeft*leftCommuter));
    kernelLeft=sparse.(MatJoin[:nullspace]())
    ###For the right link
    rightCommuter=buildSuperOperatorLeft(kron(Hcommuting, speye(SymPy.Sym,2^(maxLen-1))))
    MatJoin=sparse(vcat(TrLeft, TrRight, TrRight*rightCommuter));
    kernelRight=sparse.(MatJoin[:nullspace]())
    return kernelLeft, kernelRight
end

function generateHermitianKernel(maxLen)
    sz=-Sym[1  0 ; 0 -1];
    kernelLeft, kernelRight=buildKernels(maxLen, sparse(sz))
    hermitianKernelLeft=sparse(makeHermitian(kernelLeft))
    hermitianKernelRight=sparse(makeHermitian(kernelRight))
    return hermitianKernelLeft, hermitianKernelRight
end

function mainGeneration(maxLen)
    hermLeft, hermRight=generateHermitianKernel(maxLen)
    FileIO.save(string("kernelIsing-len-", maxLen, ".jld2"), "len", maxLen, "kernelLeft", hermLeft, "kernelRight", hermRight)
end

function extractLowerEigenvalue(vect)
    len=getLength(length(vect));
    a, _=eigs(reshape(view(vect, :), (2^len, 2^len))-10*I, nev=1, which=:LM)
    return real(a[1])+10.0
end

function improveRho_algo!(rho, kernel, gradient, temprho; de=0.0001, step=100, dx=0.1)
    for x=1:step
        target=extractLowerEigenvalue(rho)
        for j=1:length(gradient)
            temprho=rho+de*kernel[:, j]
            gradient[j]=(extractLowerEigenvalue(temprho)-target)/de*dx
        end
        for j=1:length(gradient)
            @. rho=rho+gradient[j]*kernel[:, j]
        end
    end
end

function improveRho_training!(rho, kernel; de=0.0001, step=100, dx=0.1)
    gradient=zeros(size(kernel)[2])
    temprho=similar(rho)
    totalmodif=zeros(size(kernel)[2])
    for x=1:step
        target=extractLowerEigenvalue(rho)
        for j=1:length(gradient)
            temprho=rho+de*kernel[:, j]
            gradient[j]=(extractLowerEigenvalue(temprho)-target)/de*dx
        end
        for j=1:length(gradient)
            @. rho=rho+gradient[j]*kernel[:, j]
            @. totalmodif=totalmodif+gradient
        end
    end
end

matrho=rand(Complex128, 16, 16); matrho=matrho+matrho';
rho=reshape(copy(matrho), length(matrho))
target=extractLowerEigenvalue(rho)
 improveRho!(rho, kernelLeft)
finished=extractLowerEigenvalue(rho)

max_stars_repo_name: LoicHerviou/MLQuantumMarginal

max_issues_repo_name: LoicHerviou/MLQuantumMarginal

max_forks_repo_name: LoicHerviou/MLQuantumMarginal


text:
using Distributions

m = constructdice()
run(m)

# Uncertain parameters defined in "Projections and Uncertainties about Climate
# Change in an Era of Minimal Climate Policies" (Nordhaus, 2018) EViews file
# provided by the AEJ : https://www.aeaweb.org/articles?id=10.1257/pol.20170046

mcs = @defsim begin
    t2xco2 = LogNormal(1.106, 0.2646)
	ga0 = Normal(0.076, 0.056)
	gsigma1 = Normal(-0.0152, 0.0032)
	mueq = LogNormal(5.851, 0.2649)
	a2 = Normal(0.00227, 0.001135)

    save(damages.DAMAGES)
	save(grosseconomy.YGROSS)
	save(emissions.E)
	save(co2cycle.MU)
	save(damages.DAMFRAC)
end

run(mcs, m, 1000; trials_output_filename = "/tmp/dice-2016/trialdata.csv", results_output_dir="/tmp/dice-2016")


text:
function convert_common_kwargs(opt::NonconvexJuniper.JuniperIpoptAlg, opt_kwargs;
    cb=nothing,
    maxiters=nothing,
    maxtime=nothing,
    abstol=nothing,
    reltol=nothing)

    conv_opt_kwargs = (; opt_kwargs...)

    if !isnothing(cb)
        @warn "common callback argument is currently not used by $(opt)"
    end

    if !isnothing(maxiters)
        @warn "common maxiters argument is currently not used by $(opt)"
    end

    if !isnothing(maxtime)
        conv_opt_kwargs = (; conv_opt_kwargs..., time_limit = maxtime)
    end

    if !isnothing(abstol)
        conv_opt_kwargs = (; conv_opt_kwargs..., atol = abstol)
    end

    if !isnothing(reltol)
        @warn "common reltol argument is currently not used by $(opt)"
    end

    return conv_opt_kwargs
end

function _create_options(opt::NonconvexJuniper.JuniperIpoptAlg;
    opt_kwargs=nothing,
    sub_options=nothing,
    convergence_criteria=nothing)

    if !isnothing(sub_options)
        options = (; options = !isnothing(opt_kwargs) ? NonconvexJuniper.JuniperIpoptOptions(;subsolver_options = IpoptOptions(sub_options...), opt_kwargs...) : NonconvexJuniper.JuniperIpoptOptions(subsolver_options = IpoptOptions(sub_options...)))
    else
        options = (; options = !isnothing(opt_kwargs) ? NonconvexJuniper.JuniperIpoptOptions(;subsolver_options = IpoptOptions(),opt_kwargs...) : NonconvexJuniper.JuniperIpoptOptions(subsolver_options=IpoptOptions()))
    end
    return options
end

check_optimizer_backend(opt::NonconvexJuniper.JuniperIpoptAlg) = false

include("nonconvex.jl")


text:
@diffusion_process LotkaVolterraAux{K,R} begin
    :dimensions
    process --> 2
    wiener --> 2

    :parameters
    (α, β, γ, δ, σ1, σ2) --> K
    (t, T) --> Float64
    (u, v) --> R

    :additional
    constdiff --> true
    linear --> true
end

DiffusionDefinition.B(t, P::LotkaVolterraAux) = @SMatrix [-0.0 -P.β*P.γ/P.δ; P.α*P.δ/P.β -0.0]
DiffusionDefinition.β(t, P::LotkaVolterraAux) = ℝ{2}(P.γ/P.δ*P.α, -P.α/P.β*P.γ)
DiffusionDefinition.σ(t, P::LotkaVolterraAux) = SDiagonal(P.σ1, P.σ2)


text:
# This file is part of InvertedFiles.jl

import SparseArrays: sparse

"""
    sparse(idx::BinaryInvertedFile, one::Type{RealType}=1f0)

Creates an sparse matrix (from SparseArrays) from `idx` using `one` as value.

```
   I
   ↓    1 2 3 4 5 … n  ← J
 L[1] = 0 1 0 0 1 … 0
 L[2] = 1 0 0 1 0 … 1
 L[3] = 1 0 1 0 0 … 1
 ⋮
 L[m] = 0 0 1 1 0 … 0
```
"""
function sparse(idx::BinaryInvertedFile, one::Type{RealType}=1f0) where {RealType<:Real}
    n = length(idx)
    I = eltype(idx.lists[1])[]
    J = eltype(idx.lists[1])[]
    F = RealType[]
    sizehint!(I, n)
    sizehint!(J, n)
    sizehint!(F, n)

    for i in eachindex(idx.lists)
        L = idx.lists[i]
        for j in L
            push!(I, i)
            push!(J, j)
            push!(F, one)
        end
    end

    sparse(I, J, F, length(idx.lists), n)
end

"""
    sparse(idx::WeightedInvertedFile)

Creates an sparse matrix (from SparseArrays) from `idx`
"""
function sparse(idx::WeightedInvertedFile)
    n = length(idx)
    I = eltype(idx.lists[1])[]
    J = eltype(idx.lists[1])[]
    F = eltype(idx.weights[1])[]
    sizehint!(I, n)
    sizehint!(J, n)
    sizehint!(F, n)

    for i in eachindex(idx.lists)
        L = idx.lists[i]
        W = idx.weights[i]
        for j in L
            push!(I, i)
            push!(J, j)
            push!(F, W[j])
        end
    end

    sparse(I, J, F, length(idx.lists), n)
end


text:
GRCh38_highld = [
    (chr = "1", range1 = 47761740, range2 = 51761740),
    (chr = "1", range1 = 125169943, range2 = 125170022),
    (chr = "1", range1 = 144106678, range2 = 144106709),
    (chr = "1", range1 = 181955019, range2 = 181955047),
    (chr = "2", range1 = 85919365, range2 = 100517106),
    (chr = "2", range1 = 135275091, range2 = 135275210),
    (chr = "2", range1 = 182427027, range2 = 189427029),
    (chr = "3", range1 = 47483505, range2 = 49987563),
    (chr = "3", range1 = 83368158, range2 = 86868160),
    (chr = "5", range1 = 44464140, range2 = 51168409),
    (chr = "5", range1 = 129636407, range2 = 132636409),
    (chr = "6", range1 = 25391792, range2 = 33424245),
    (chr = "6", range1 = 57788603, range2 = 58453888),
    (chr = "6", range1 = 61109122, range2 = 61424451),
    (chr = "6", range1 = 139637169, range2 = 142137170),
    (chr = "7", range1 = 54964812, range2 = 66897578),
    (chr = "8", range1 = 8105067, range2 = 12105082),
    (chr = "8", range1 = 43025699, range2 = 48924888),
    (chr = "8", range1 = 110918594, range2 = 113918595),
    (chr = "9", range1 = 64198500, range2 = 64200392),
    (chr = "10", range1 = 36671065, range2 = 43184546),
    (chr = "11", range1 = 88127183, range2 = 91127184),
    (chr = "12", range1 = 32955798, range2 = 41319931),
    (chr = "20", range1 = 33948532, range2 = 36438183),
]

GRCh37_highld = [
    (chr = "1", range1 = 48000000, range2 = 52000000),
    (chr = "2", range1 = 86000000, range2 = 100500000),
    (chr = "2", range1 = 134500000, range2 = 138000000),
    (chr = "2", range1 = 183000000, range2 = 190000000),
    (chr = "3", range1 = 47500000, range2 = 50000000),
    (chr = "3", range1 = 83500000, range2 = 87000000),
    (chr = "3", range1 = 89000000, range2 = 97500000),
    (chr = "5", range1 = 44500000, range2 = 50500000),
    (chr = "5", range1 = 98000000, range2 = 100500000),
    (chr = "5", range1 = 129000000, range2 = 132000000),
    (chr = "5", range1 = 135500000, range2 = 138500000),
    (chr = "6", range1 = 24000000, range2 = 36000000),
    (chr = "6", range1 = 57000000, range2 = 64000000),
    (chr = "6", range1 = 140000000, range2 = 142500000),
    (chr = "7", range1 = 55000000, range2 = 66000000),
    (chr = "8", range1 = 7000000, range2 = 13000000),
    (chr = "8", range1 = 43000000, range2 = 50000000),
    (chr = "8", range1 = 112000000, range2 = 115000000),
    (chr = "10", range1 = 37000000, range2 = 43000000),
    (chr = "11", range1 = 46000000, range2 = 57000000),
    (chr = "11", range1 = 87500000, range2 = 90500000),
    (chr = "12", range1 = 33000000, range2 = 40000000),
    (chr = "12", range1 = 109500000, range2 = 112000000),
    (chr = "20", range1 = 32000000, range2 = 34500000),
    (chr = "17", range1 = 42900000, range2 = 45100000)
]
# https://genome.sph.umich.edu/wiki/Regions_of_high_linkage_disequilibrium_(LD)

GRCh38_totlength = Dict(
    "1" => 248956422,
    "2" => 242193529,
    "3" => 198295559,
    "4" => 190214555,
    "5" => 181538259,
    "6" => 170805979,
    "7" => 159345973,
    "8" => 145138636,
    "9" => 138394717,
    "10" => 133797422,
    "11" => 135086622,
    "12" => 133275309,
    "13" => 114364328,
    "14" => 107043718,
    "15" => 101991189,
    "16" => 90338345,
    "17" => 83257441,
    "18" => 80373285,
    "19" => 58617616,
    "20" => 64444167,
    "21" => 46709983,
    "22" => 50818468,
    "X" => 156040895,
    "Y" => 57227415,
)

GRCh37_totlength = Dict(
    "1" => 249250621,
    "2" => 243199373,
    "3" => 198022430,
    "4" => 191154276,
    "5" => 180915260,
    "6" => 171115067,
    "7" => 159138663,
    "8" => 146364022,
    "9" => 141213431,
    "10" => 135534747,
    "11" => 135006516,
    "12" => 133851895,
    "13" => 115169878,
    "14" => 107349540,
    "15" => 102531392,
    "16" => 90354753,
    "17" => 81195210,
    "18" => 78077248,
    "19" => 59128983,
    "20" => 63025520,
    "21" => 48129895,
    "22" => 51304566,
    "X" => 155270560,
    "Y" => 59373566,
)
# https://www.ncbi.nlm.nih.gov/grc/human/data?asm=GRCh37


text:
import GitHub
using Nanosoldier, Test, BenchmarkTools
using Nanosoldier: BuildRef, JobSubmission, Config, BenchmarkJob, PkgEvalJob, AbstractJob
using BenchmarkTools: TrialEstimate, Parameters
using DataFrames

#########
# setup #
#########

vinfo = """
Julia Version 0.4.3-pre+6
Commit adffe19 (2015-12-11 00:38 UTC)
Platform Info:
  System: Darwin (x86_64-apple-darwin13.4.0)
  CPU: Intel(R) Core(TM) i5-4288U CPU @ 2.60GHz
  WORD_SIZE: 64
  uname: Darwin 13.4.0 Darwin Kernel Version 13.4.0: Wed Dec 17 19:05:52 PST 2014; root:xnu-2422.115.10~1/RELEASE_X86_64 x86_64 i386
Memory: 8.0 GB (1025.07421875 MB free)
Uptime: 646189.0 sec
Load Avg:  1.38427734375  1.416015625  1.41455078125
Intel(R) Core(TM) i5-4288U CPU @ 2.60GHz:
       speed         user         nice          sys         idle          irq
#1  2600 MHz     186848 s          0 s     114241 s    1462955 s          0 s
#2  2600 MHz     114716 s          0 s      43882 s    1605408 s          0 s
#3  2600 MHz     189864 s          0 s      77629 s    1496513 s          0 s
#4  2600 MHz     114652 s          0 s      42497 s    1606856 s          0 s

  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.3
"""
auth = if haskey(ENV, "GITHUB_AUTH")
    GitHub.authenticate(ENV["GITHUB_AUTH"])
else
    GitHub.AnonymousAuth()
end
latest_commit = GitHub.commits("JuliaLang/julia"; auth, page_limit=1)[1][1].sha
primary = BuildRef("JuliaLang/julia", latest_commit, vinfo)
against = BuildRef("JuliaLang/julia", "bb73f3489d837e3339fce2c1aab283d3b2e97a4c", vinfo*"_against")
config = Config("user", Dict(Any => [getpid()]), auth, "test", trackrepo=primary.repo);
tagpred = "ALL && !(\"tag1\" || \"tag2\")"
pkgsel = "[\"Example\"]"

#####################################
# submission parsing and validation #
#####################################

function build_test_submission(jobtyp, submission_string)
    func, args, kwargs = Nanosoldier.parse_submission_string(submission_string)
    submission = JobSubmission(config, primary, primary.sha, "https://www.test.com", :commit, nothing, func, args, kwargs)
    @test Nanosoldier.isvalid(submission, jobtyp)
    return jobtyp(submission)
end

job = build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL)`")
@test job.against === nothing
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(\"tag\")`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks($tagpred)`")

job = build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL, vs = \"JuliaLang/julia:master\")`")
@test job.against !== nothing
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(\"tag\", vs = \"JuliaLang/julia:master\")`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks($tagpred, vs = \"JuliaLang/julia:master\")`")

job = build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL, vs = \"%self\")`")
@test job.against == primary
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(\"tag\", vs = \"%self\")`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks($tagpred, vs = \"%self\")`")

@test_throws ErrorException("invalid commit to run isdaily") build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL, isdaily = true, vs = \"JuliaLang/julia:master\")`")
@test_throws ErrorException("invalid commit to run isdaily") build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(\"tag\", isdaily = true, vs = \"JuliaLang/julia:master\")`")
@test_throws ErrorException("invalid commit to run isdaily") build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks($tagpred, isdaily = true, vs = \"JuliaLang/julia:master\")`")

@test_throws ErrorException("invalid commit to run isdaily") build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL; isdaily = true, vs = \"JuliaLang/julia:master\")`")
@test_throws ErrorException("invalid commit to run isdaily") build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(\"tag\"; isdaily = true, vs = \"JuliaLang/julia:master\")`")
@test_throws ErrorException("invalid commit to run isdaily") build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks($tagpred; isdaily = true, vs = \"JuliaLang/julia:master\")`")

job = build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL, vs = \"JuliaLang/julia#v1.0.0\")`")
@test job.against !== nothing
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(\"tag\", vs = \"JuliaLang/julia#v1.0.0\")`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks($tagpred, vs = \"JuliaLang/julia#v1.0.0\")`")

job = build_test_submission(PkgEvalJob, "@nanosoldier `runtests(ALL)`")
@test job.against === nothing
build_test_submission(PkgEvalJob, "@nanosoldier `runtests(\"pkg\")`")
build_test_submission(PkgEvalJob, "@nanosoldier `runtests($pkgsel)`")

job = build_test_submission(PkgEvalJob, "@nanosoldier `runtests($pkgsel)`")
@test job.compiled == :none
job = build_test_submission(PkgEvalJob, "@nanosoldier `runtests($pkgsel, compiled=:primary)`")
@test job.compiled == :primary
job = build_test_submission(PkgEvalJob, "@nanosoldier `runtests($pkgsel, compiled=:against)`")
@test job.compiled == :against
job = build_test_submission(PkgEvalJob, "@nanosoldier `runtests($pkgsel, compiled=:both)`")
@test job.compiled == :both

#############################
# retrieval from job queue  #
#############################

non_daily_job = build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL)`")
daily_job = build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL, isdaily = true)`")

queue = [daily_job, daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, true)
@test job !== nothing && job.isdaily
@test length(queue) == 1

queue = [non_daily_job, daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, true)
@test job !== nothing && !job.isdaily
@test length(queue) == 1

queue = [daily_job, non_daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, true)
@test job !== nothing && job.isdaily
@test length(queue) == 1

queue = [daily_job, daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, false)
@test job === nothing
@test length(queue) == 2

queue = [non_daily_job, daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, false)
@test job !== nothing && !job.isdaily
@test length(queue) == 1

queue = [daily_job, non_daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, false)
@test job !== nothing && !job.isdaily
@test length(queue) == 1

queue = [daily_job, non_daily_job]
job = Nanosoldier.retrieve_job!(queue, PkgEvalJob, true)
@test job === nothing
@test length(queue) == 2

queue = [daily_job, non_daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, true)
@test job !== nothing && job.isdaily
@test length(queue) == 1

queue = [daily_job, non_daily_job]
job = Nanosoldier.retrieve_job!(queue, Any, false)
@test job !== nothing && !job.isdaily
@test length(queue) == 1

#########################
# job report generation #
#########################

job = build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks($tagpred)`")
@test job.tagpred == tagpred
@test job.against == nothing
job.against = against

results = Dict(
    "primary" => BenchmarkGroup([],
        "g" => BenchmarkGroup([],
            "h" => BenchmarkGroup([],
                "x"      => TrialEstimate(Parameters(), 1.0, 3.5, 1.0, 1.0),  # invariant
                ("y", 1) => TrialEstimate(Parameters(memory_tolerance = 0.03), 2.0, 1.0, 0.0, 1.0),  # regression/improvement
                ("y", 2) => TrialEstimate(Parameters(time_tolerance = 0.04), 0.5, 1.0, 1.0, 1.0),  # improvement
                "z"      => TrialEstimate(Parameters(memory_tolerance = 0.27, time_tolerance = 0.6), 1.0, 1.0, 5.0, 1.0),  # regression
                "∅"      => TrialEstimate(Parameters(), 1.0, 1.0, 1.0, 1.0) # not in "against" group
            )
        )
    ),
    "against" => BenchmarkGroup([],
        "g" => BenchmarkGroup([],
            "h" => BenchmarkGroup([],
                "x"      => TrialEstimate(Parameters(), 1.0, 1.0, 1.0, 1.0),
                ("y", 1) => TrialEstimate(Parameters(), 1.0, 1.0, 1.0, 1.0),
                ("y", 2) => TrialEstimate(Parameters(), 1.0, 1.0, 1.0, 1.0),
                "z"      => TrialEstimate(Parameters(), 1.0, 1.0, 1.0, 1.0)
            )
        )
    )
)

results["judged"] = BenchmarkTools.judge(results["primary"], results["against"])

@test begin
    mdpath = joinpath(@__DIR__, "report.md")
    md = replace(read(mdpath, String), "PRIMARY" => latest_commit)
    md2 = sprint(io->Nanosoldier.printreport(io, job, results))
    chomp.(eachline(IOBuffer(md))) == chomp.(eachline(IOBuffer(md2)))
end

@testset "Markdown" begin
    @test Nanosoldier.markdown_escaped("abc") == "abc"
    @test Nanosoldier.markdown_escaped(raw"a\`*_#+-.!{}[]()<>|b") == raw"a\\\`\*\_\#\+\-\.\!\{\}\[\]\(\)\<\>\|b"

    @test Nanosoldier.markdown_escaped_code("abc") == "`abc`"
    @test Nanosoldier.markdown_escaped_code("a`b`c") == "``a`b`c``"
    @test Nanosoldier.markdown_escaped_code("``ab`c") == "``` ``ab`c```"
    @test Nanosoldier.markdown_escaped_code("a`bc```") == "````a`bc``` ````"
end

@testset "PkgEvalJob" begin
    job = build_test_submission(PkgEvalJob, "@nanosoldier `runtests($pkgsel)`")
    @test job.against == nothing

    results = Dict{String,Any}(
        "primary" => DataFrame(
            julia=v"1.6",
            name="Example",
            uuid=Base.UUID("12345678-1234-1234-1234-123456789abc"),
            version=v"0.1.0",
            status=:ok,
            reason=missing,
            duration=0.0,
            log="everything is fine",
        )
    )

    report = sprint(io -> Nanosoldier.printreport(io, job, results))

    job.against = against
    results["against"] = DataFrame(
            julia=v"1.7",
            name="Example",
            uuid=Base.UUID("12345678-1234-1234-1234-123456789abc"),
            version=v"0.1.0",
            status=:fail,
            reason=:test_failures,
            duration=0.0,
            log="this one failed",
        )

    report = sprint(io -> Nanosoldier.printreport(io, job, results))

    # XXX: actually test contents?
end

nothing


text:
# ---
# title: 1217. Minimum Cost to Move Chips to The Same Position
# id: problem1217
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Math, Greedy
# link: <https://leetcode.com/problems/minimum-cost-to-move-chips-to-the-same-position/description/>
# hidden: true
# ---
#
# We have `n` chips, where the position of the `ith` chip is `position[i]`.
#
# We need to move all the chips to **the same position**. In one step, we can
# change the position of the `ith` chip from `position[i]` to:
#
#   * `position[i] + 2` or `position[i] - 2` with `cost = 0`.
#   * `position[i] + 1` or `position[i] - 1` with `cost = 1`.
#
# Return _the minimum cost_ needed to move all the chips to the same position.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/15/chips_e1.jpg)
#
#
#
#     Input: position = [1,2,3]
#     Output: 1
#     Explanation: First step: Move the chip at position 3 to position 1 with cost = 0.
#     Second step: Move the chip at position 2 to position 1 with cost = 1.
#     Total cost is 1.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/15/chip_e2.jpg)
#
#
#
#     Input: position = [2,2,2,3,3]
#     Output: 2
#     Explanation: We can move the two chips at position  3 to position 2. Each move has cost = 1. The total cost = 2.
#
#
# **Example 3:**
#
#
#
#     Input: position = [1,1000000000]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= position.length <= 100`
#   * `1 <= position[i] <= 10^9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#=

 From http:  reasoning.eas.asu.edu/lpmln/examples/lpmln2al/taxi.lpmln
  """
  Conside a program where you know that you are going to be Late
  if you don't get a Taxi.

  Since we have put a hard constraint that we cannot be late the probability of being
  Late is going to be 0. The probability of there being a NoTaxi comes out to be 71 %.
  {NoTaxi = T, Late = F} satisfies the second soft rule and the hard constraint. Hence
  the probability for this model comes out to be e\textsuperscript{2}/(
  e\textsuperscript{2} + e\textsuperscript{1}) which comes out to be 73\%
  \begin{lstlisting}
  Late 4.9995e-05
  NoTaxi 0.711979
  \end{lstlisting}
  """

  Note: Since lpmln use weights instead of probability, I don't know how to get the no_taxi as 0.71
  Hmm in the lpmln there's an example where the weights are corresponding to percentages:
     0.5 -> -1.6094
     0.8 -> -0.2231
  which is ln(0.5) and ln(0.8)

  But this don't help us here, since the weights are positive, 1 and 2....

  Cf ~/cplint/taxi.pl
     ~/blog/taxi.blog
     ~/psi/taxi.psi
     ~/webpppl/taxi.wppl

=#

using Turing, StatsPlots, DataFrames

@model function taxi()
    # There is a chance for a taxi to be not available
    no_taxi ~ Bernoulli(0.3)

    # You are going to be late if you don't get a taxi
    late ~ no_taxi ? Bernoulli(0.65) : Bernoulli(0);

    # But you certainly cannot be late no matter what situation.
    # Hence this constraint is encoded as a hard rule
    late == false || begin Turing.@addlogprob! -Inf; return end
end

model = taxi()

num_chains = 4
# chains = sample(model, Prior(), MCMCThreads(), 1000, num_chains)

chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains) # 1.47s
# chains = sample(model, IS(), MCMCThreads(), 10_000, num_chains) # 1.48s

# chains = sample(model, PG(20), MCMCThreads(), 10_000, num_chains) # Too slow!
# chains = sample(model, SMC(), MCMCThreads(), 10_000, num_chains) # 4.1s, don't handle observes!

# chains = sample(model, HMC(0.05,10), MCMCThreads(), 10_000, num_chains) # Nope!
# chains = sample(model, NUTS(), MCMCThreads(), 10_000, num_chains) # Nope!

display(chains)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Prjtrellis"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/YosysHQ/prjtrellis.git", "dff1cbcb1bd30de7e96f8a059f2e19be1bb2e44d")
]

dependencies = [
    Dependency("Python_jll"; compat="~3.8.1"),
    Dependency("boost_jll"; compat="=1.76.0"), # max gcc7
]

# Bash recipe for building across all platforms
script = raw"""
cd prjtrellis/libtrellis
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> Sys.islinux(p), supported_platforms())
platforms = expand_cxxstring_abis(platforms)
# For some reason, building for CXX03 string ABI doesn't actually work, skip it
filter!(x -> cxxstring_abi(x) != "cxx03", platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("ecpunpack", :ecpunpack),
    ExecutableProduct("ecppll", :ecppll),
    ExecutableProduct("ecpbram", :ecpbram),
    ExecutableProduct("ecppack", :ecppack),
    ExecutableProduct("ecpmulti", :ecpmulti),
    ExecutableProduct("ecpmulti", :ecpmulti),
    LibraryProduct("libtrellis", :libtrellis, ["lib/trellis", "lib64/trellis"]),
    LibraryProduct("pytrellis", :pytrellis, ["lib/trellis", "lib64/trellis"])
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# ---
# title: 73. Set Matrix Zeroes
# id: problem73
# author: Qling
# date: 2021-03-25
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/set-matrix-zeroes/description/>
# hidden: true
# ---
#
# Given an ` _m_ x _n_` matrix. If an element is **0** , set its entire row and
# column to **0**. Do it [**in-place**](https://en.wikipedia.org/wiki/In-
# place_algorithm).
#
# **Follow up:**
#
#   * A straight forward solution using O( _m_ _n_ ) space is probably a bad idea.
#   * A simple improvement uses O( _m_ \+ _n_ ) space, but still not the best solution.
#   * Could you devise a constant space solution?
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/17/mat1.jpg)
#
#
#
#     Input: matrix = [[1,1,1],[1,0,1],[1,1,1]]
#     Output: [[1,0,1],[0,0,0],[1,0,1]]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/17/mat2.jpg)
#
#
#
#     Input: matrix = [[0,1,2,0],[3,4,5,2],[1,3,1,5]]
#     Output: [[0,0,0,0],[0,4,5,0],[0,3,1,0]]
#
#
#
#
# **Constraints:**
#
#   * `m == matrix.length`
#   * `n == matrix[0].length`
#   * `1 <= m, n <= 200`
#   * `-231 <= matrix[i][j] <= 231 - 1`
#
#
## @lc code=start
using LeetCode

function set_zeroes(matrix::Vector{Vector{Int}})::Vector{Vector{Int}}
    is_col = false
    m, n = length(matrix), length(matrix[1])

    for i = 1:m
        matrix[i][1] == 0 && (is_col = true)

        for j = 2:n
            matrix[i][j] == 0 && (matrix[i][1] = matrix[1][j] = 0)
        end
    end

    for i = 2:m, j = 2:n
        (matrix[i][1] == 0 || matrix[1][j] == 0) && (matrix[i][j] = 0)
    end

    ## See if the first row needs to be set to zero as well
    if matrix[1][1] == 0
        for j = 1:n
            matrix[1][j] = 0
        end
    end

    ## See if the first column needs to be set to zero as well
    if is_col
        for i = 1:m
            matrix[i][1] = 0
        end
    end

    return matrix
end
## @lc code=end


text:
using MCTS
include("roadnet_pursuer_generator_MDP.jl")
using ProgressMeter
#  using PyPlot
using JLD
using MicroLogging
using Base.Markdown
using DataFrames, CSV
include("road_net_visualize.jl")
include("utilities.jl")
include("self_confidence.jl")

rwd_global = []

####### NEED TO RE-VISIT HOW THESE FILES WORK, THERE ARE A LOT OF MULTIPLE IMPORTS
####### I PROBABLY NEED TO PACKAGE THIS STUFF TO AVOID ISSUES WHEN I TRY TO IMPORT @EVERYWHERE
####### THIS IS CAUSING ISSUES WHEN TRYING TO RUN IN PARALLEL, I GOT A "method too new to be called from this world context" error

function get_reward_distribution(g::MetaGraph, mdp::roadnet_with_pursuer; max_steps::Int64=25,
                        its_vals::Array{Int64}=[5],d_vals::Array{Int64}=[2],exp_vals::Array{Float64}=[5.],
                        num_repeats::Int64=10,discounted_rwd::Bool=true,
                        initial_state::roadnet_pursuer_state=roadnet_pursuer_state(1,4))
    rnet = mdp.road_net.gprops[:net_stats]
    net_diam = rnet.diam

    if discounted_rwd
        dis_str = "Discounted Rwd"
    else
        dis_str = "Rwd"
    end

    @info "calculating distribution"
    @debug "Using:" max_steps its_vals d_vals num_repeats rnet mdp.tprob mdp.discount mdp.road_net.gprops[:reward_dict] mdp.road_net.gprops[:exit_nodes]

    # TODO: investigate what the exploration constant does, i.e. does it scale with N?
    # This instantiates empty MCTS policies because MCTS runs "online"
    policy_tilde = Vector{MCTS.MCTSPlanner}(0)
    hist_tilde = Vector{HistoryRecorder}(0)
    its_axis = Vector{Int64}(0)
    d_axis = Vector{Int64}(0)

    for i in 1:length(its_vals)
        for d in d_vals
            its = its_vals[i]
            s = MCTSSolver(n_iterations=its,depth=d,exploration_constant=exp_vals[i],enable_tree_vis=true)
            push!(its_axis,its)
            push!(d_axis,d)
            push!(policy_tilde,solve(s,mdp))
            push!(hist_tilde,HistoryRecorder(max_steps=max_steps))
        end
    end
    @debug "Total number of solvers being tested:" length(policy_tilde)

    function reward_grab(mdp::roadnet_with_pursuer,policy::MCTS.MCTSPlanner,hist::HistoryRecorder,s::roadnet_pursuer_state;discounted::Bool=true,return_hist::Bool=false)
        hist = simulate(hist,mdp,policy,s)
        if return_hist
            return hist
        end

        # discounted reward takes time to solution into account
        if discounted
            return discounted_reward(hist)
        else
            return undiscounted_reward(hist)
        end
    end

    PT = repmat(policy_tilde,1,num_repeats)
    HT = repmat(hist_tilde,1,num_repeats)
    IA = repmat(its_axis,1,num_repeats)
    DA = repmat(d_axis,1,num_repeats)
    @debug "PT $(size(PT)),HT $(size(HT)),IA $(size(IA)),DA $(size(DA))"
    rewards = Array{Float64}(size(PT))
    sim_time = Array{Int64}(size(PT))
    idx = 1

    @info md"# Simulations"
    tic()

    if nprocs() > 1
        q = []
        for pol in PT
            value = Sim(mdp,pol,initial_state,max_steps=max_steps)
            push!(q,value)
        end
        sim_results = run_parallel(q) do sim,hist
            display(discounted_reward(hist))
            return [:steps=>n_steps(hist), :reward=>discounted_reward(hist)]
        end

        rewards = [float(x) for x in sim_results[:reward]] #need to convert to float, because sim_results[:reward] is a union of floats and missing values, even though there aren't any missing values...

    else
        pm = Progress(length(PT[:]),1)
        for (p,h) in zip(PT,HT)
            @info "Running Simulations" progress=idx/length(PT)
            rewards[idx] = reward_grab(mdp,p,h,initial_state,discounted=discounted_rwd)
            hist = reward_grab(mdp,p,h,initial_state,discounted=discounted_rwd,return_hist=true)
            sim_time[idx] = length(hist)
            if false
              @debug print_hist(mdp,hist)
            end
            next!(pm)
            idx += 1
        end
    end
    @info "Completed $(length(PT)) Simulations in $(toq()) seconds"

    R = reshape(rewards,length(its_axis),:)
    ST = reshape(sim_time,length(its_axis),:)

    # Save off data
    jld_fname = "data/$(num_repeats)_$(DateTime(now())).jld"
    @info "saving data to `$jld_fname`"
    JLD.save(jld_fname,"its_axis",its_axis,"d_axis",d_axis,"rewards",rewards,
         "u_vals",mean(rewards,2),"max_steps",max_steps,"R",R,"ST",ST,"IA",IA,"DA",DA)

    return vec(R) #convert to vector
end

function netprops2array(g::roadnet_with_pursuer)
    discount = g.discount
    exits = g.exit_nodes
    rwd_states = g.reward_states # not used right now
    rwd_vals = g.reward_vals # not used right now
    tprob = g.tprob

    num_exit_nodes = length(exits)

    rnet = g.road_net
    rnet_props = rnet.gprops
    exit_rwd = rnet_props[:reward_dict][:exit]
    caught_rwd = rnet_props[:reward_dict][:caught]
    sensor_rwd = rnet_props[:reward_dict][:sensor]
    rnet_stats = rnet_props[:net_stats]
    avg_degree = rnet_stats.average_degree
    deg_variance = rnet_stats.degree_variance
    diam = rnet_stats.diam
    max_degree = rnet_stats.maximum_degree
    N = rnet_stats.nv
    E = rnet_stats.ne

    netprops_ary = [discount tprob num_exit_nodes exit_rwd caught_rwd sensor_rwd avg_degree deg_variance diam max_degree N E]

    return netprops_ary

end

function make_training_data(net::Dict)
    net_dict = net[:sims]
    reps = net_dict[:repeats]
    log_to_file = net_dict[:log_to_file]
    log_fname = net_dict[:log_fname]
    log_lvl = net_dict[:log_lvl]
    sim_steps = net_dict[:steps]
    dis_rwd = net_dict[:dis_rwd]
    fname = net[:fname]

    fn = joinpath(fldr,fname)

    make_training_data(data_fname=fn,repeats=reps,sim_steps=sim_steps,dis_rwd=dis_rwd)
end

function make_training_data(;data_fname::String="nets",logtofile::Bool=false, logfname::String="logs/$(now()).log",loglvl::Symbol=:debug,repeats::Int64=25,sim_steps::Int64=-1,dis_rwd::Bool=false)

    # initialize the DataFrame
    training_data = DataFrame(graphID=Int64[],discount=Float64[],tprob=Float64[],num_exit_nodes=Float64[],exit_rwd=Float64[],
                              caught_rwd=Float64[],sensor_rwd=Float64[],avg_degree=Float64[],deg_variance=Float64[],
                              diam=Float64[],max_degree=Float64[],N=Float64[],E=Float64[],its=Float64[],
                              e_mcts=Float64[],d_mcts=Float64[],steps=Float64[],repeats=Float64[],
                              exit_distance=Float64[],pursuer_distance=Float64[],X3_1=Float64[],X3_2=Float64[],X4=Float64[],
                              mean=Float64[],median=Float64[],moment_2=Float64[],moment_3=Float64[],
                              moment_4=Float64[],moment_5=Float64[],moment_6=Float64[],moment_7=Float64[],
                              moment_8=Float64[],moment_9=Float64[],moment_10=Float64[])
    #### Logging
    if logtofile
        buffer = IOBuffer()
        logger = SimpleLogger(buffer)
        configure_logging(min_level=loglvl)
    else
        logger = global_logger()
        configure_logging(min_level=loglvl)
    end

    # get list of graphs
    problem_dict = jldopen("$data_fname.jld", "r") do file
        @info "Grabbing Pre-Generated Networks"
        read(file,"problem_dict")
    end

    global_rwd_range = Array{Float64}(2,1) # idx 1 is minimun, idx 2 is maximum

    i = 0
    for problem in problem_dict
        i += 1
        display(problem_dict)
        println()
        display(problem)
        println()
        net_num = problem[1]
        problem = problem[2]
        if :error in keys(problem)
            println(problem)
            continue
        else
            @info md"# Processing Network $i of $(length(problem_dict))"
            g = problem[:graph]
            t = problem[:trans_prob]
            d = problem[:discount]

            # do stuff
            # i'm thinking using the same network multiple times using different discount factors, and transition probabilities, this will do two things: 1) make training data better, 2) increase total amount of training data
            mdp = roadnet_with_pursuer(g,tp=t,d=d)

            its = [problem[:mcts_its]] # perhaps base this on transition probability
            d_mcts = [problem[:mcts_depth]]
            e_mcts = [problem[:mcts_e]]

            if sim_steps == -1
                # -1 signifies auto calculation
                # throw the book at it!
                steps = 20*mdp.road_net.gprops[:net_stats].diam
            else
                steps = sim_steps
            end

            start = problem[:evader_start]
            escape = problem[:exit_loc]
            pursuer = problem[:pursuer_start]
            start_state = roadnet_pursuer_state(start,pursuer)
            exit_distance = a_star(mdp.road_net,start,escape)
            pursuer_distance = a_star(mdp.road_net,start,pursuer)

            r_dist = []
            with_logger(logger) do
                r_dist = get_reward_distribution(g,mdp,its_vals=its,d_vals=d_mcts,exp_vals=e_mcts,max_steps=steps,num_repeats=repeats,discounted_rwd=dis_rwd,initial_state=start_state)
            end

            # update global reward range
            minr = minimum(r_dist)
            maxr = maximum(r_dist)
            if i==1
                global_rwd_range[1] = minr
                global_rwd_range[2] = maxr
            else
                minr < global_rwd_range[1] && (global_rwd_range[1] = minr)
                maxr > global_rwd_range[2] && (global_rwd_range[2] = maxr)
            end

            @info "calculating X3 data"
            SQ_data = X3(r_dist)
            @debug SQ_data

            @info "calculating X4 data"
            OA_data = X4(r_dist,threshold=mean(r_dist))
            @debug OA_data

            @info "calculating statistical moments"
            r_dist_moments = [mean(r_dist) median(r_dist) moment(r_dist,2) moment(r_dist,3) moment(r_dist,4) moment(r_dist,5) moment(r_dist,6) moment(r_dist,7) moment(r_dist,8) moment(r_dist,9) moment(r_dist,10)]
            @debug r_dist_moments

            @info "collecting net props"
            mdp_props_ary = netprops2array(mdp)
            @debug mdp_props_ary

            @info "collecting solver props"
            solver_props_ary = [its[1] e_mcts[1] d_mcts[1] steps repeats length(exit_distance) length(pursuer_distance)]
            @debug solver_props_ary

            @info "combining to a row"
            data_entry = [net_num mdp_props_ary solver_props_ary SQ_data OA_data r_dist_moments]
            @debug length(net_num) length(mdp_props_ary) length(solver_props_ary) length(SQ_data) length(r_dist_moments)
            @debug data_entry
            @debug length(data_entry)
            @debug training_data
            push!(training_data,data_entry)

            # store training data with individual network for future use
            training_data_dict = Dict()
            n_cntr = 1
            for n in names(training_data)
                training_data_dict[n] =  data_entry[n_cntr]
                n_cntr += 1
            end
            problem[:training_data] = training_data_dict
            problem[:training_data][:r_dist] = r_dist
        end
    end
    CSV.write("$data_fname.csv",training_data)

    # re-write original file with "r_dist" object added to each network
    jldopen("$data_fname.jld", "w") do file
        write(file, "problem_dict", problem_dict)
        write(file, "global_rwd_range", global_rwd_range)
    end

    if logtofile
        # using the "open(...) do f ... syntax printed the length of the file at the end, this ways doesn't do that
        f = open(logfname,"w")
        logtxt = strip(String(take!(buffer)));
        logtxt = replace(logtxt, r"^\e\[1m\e\[..m(.- )\e\[39m\e\[22m", s"\1")
        write(f,logtxt)
        close(f)
    end
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
type Sequence
    chain::Task
end

function (seq::Sequence)(f::Function, g::Function)
    seq.chain = init_task(bridge, seq.chain, init_task(f, g))
    return seq
end

doc"""
As the constructor of Sequence, turning Vector of something into a iterator Task
wrapped in Sequence.
"""
function from(s::Vector)
    function start()
        produce()
        for x in s
            produce(Subject(x))
        end
        return :done
    end
    return Sequence(init_task(start))
end

function create(func::Function)

end

function from_iter()

end

function to_list(s::Sequence; debug=false)
    result = Vector()
    while true
        x = consume(s.chain)
        if debug
            println(x)
        end
        if x == :done
            break
        else
            append!(result, get(x))
        end
    end
    return result
end

function to_set(s::Sequence; debug=false)
    result = Set()
    while true
        x = consume(s.chain)
        if debug
            println(x)
        end
        if x == :done
            break
        else
            push!(result, get(x))
        end
    end
    return result
end

function to_dict(s::Sequence; debug=false)
    result = Dict()
    while true
        x = consume(s.chain)
        if debug
            println(x)
        end
        if x == :done
            break
        else
            k, v = get(x)
            result[k] = v
        end
    end
    return result
end


text:
import Base.+, Base.-, Base./, Base.*, Base.^
import SchumakerSpline.evaluate
function evaluate(f::Undefined_Function, point::Union{Real,Date,DateTime,DatePeriod})
    missing
end
function (s::Undefined_Function)(x::Union{Real,Date,DateTime,DatePeriod})
    return missing
end

"""
    derivative(f::UnivariateFunction)
This calculates the derivative of the function and returns it as a `UnivariateFunction`.

### Inputs
* `f` - A `UnivariateFunction`.
### Returns
* A `UnivariateFunction`.
"""
function derivative(f::Undefined_Function)
    return f
end

"""
    indefinite_integral(f::UnivariateFunction)
This calculates the indefinite integral of a `UnivariateFunction`.

### Inputs
* `f` - A `UnivariateFunction`.
### Returns
* A `UnivariateFunction`.
"""
function indefinite_integral(f::Undefined_Function)
    return f
end

function +(f::Undefined_Function, number::Real)
    return f
end
function -(f::Undefined_Function, number::Real)
    return f
end
function *(f::Undefined_Function, number::Real)
    return f
end
function /(f::Undefined_Function, number::Real)
    return f
end
function ^(f::Undefined_Function, number::Integer)
    return f
end

function +(f1::Undefined_Function, f2::Undefined_Function)
    return f1
end
function +(f1::Undefined_Function, f2::PE_Function)
    return f1
end
function +(f1::Undefined_Function, f2::Sum_Of_Functions)
    return f1
end
function +(f1::Undefined_Function, f2::Piecewise_Function)
    return f1
end

function +(f1::PE_Function, f2::Undefined_Function)
    return +(f2,f1)
end
function +(f1::Sum_Of_Functions, f2::Undefined_Function)
    return +(f2,f1)
end
function +(f1::Piecewise_Function, f2::Undefined_Function)
    return +(f2,f1)
end

function -(f1::Undefined_Function, f2::Undefined_Function)
    return f1
end
function -(f1::Undefined_Function, f2::PE_Function)
    return f1
end
function -(f1::Undefined_Function, f2::Sum_Of_Functions)
    return f1
end
function -(f1::Undefined_Function, f2::Piecewise_Function)
    return f1
end

function -(f1::PE_Function, f2::Undefined_Function)
    return -(f2,f1)
end
function -(f1::Sum_Of_Functions, f2::Undefined_Function)
    return -(f2,f1)
end
function -(f1::Piecewise_Function, f2::Undefined_Function)
    return -(f2,f1)
end

function *(f1::Undefined_Function, f2::Undefined_Function)
    return f1
end
function *(f1::Undefined_Function, f2::PE_Function)
    return f1
end
function *(f1::Undefined_Function, f2::Sum_Of_Functions)
    return f1
end
function *(f1::Undefined_Function, f2::Piecewise_Function)
    return f1
end

function *(f1::PE_Function, f2::Undefined_Function)
    return *(f2,f1)
end
function *(f1::Sum_Of_Functions, f2::Undefined_Function)
    return *(f2,f1)
end
function *(f1::Piecewise_Function, f2::Undefined_Function)
    return *(f2,f1)
end

function /(f1::Undefined_Function, f2::Undefined_Function)
    return f1
end
function /(f1::Undefined_Function, f2::PE_Function)
    return f1
end
function /(f1::Undefined_Function, f2::Sum_Of_Functions)
    return f1
end
function /(f1::Undefined_Function, f2::Piecewise_Function)
    return f1
end

function /(f1::PE_Function, f2::Undefined_Function)
    return /(f2,f1)
end
function /(f1::Sum_Of_Functions, f2::Undefined_Function)
    return /(f2,f1)
end
function /(f1::Piecewise_Function, f2::Undefined_Function)
    return /(f2,f1)
end