text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    add_hot_storage_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")

Adds the Storage results to the dictionary passed back from `run_reopt` using the solved model `m` and the `REoptInputs` for node `_n`.
Note: the node number is an empty string if evaluating a single `Site`.

Storage results:
- `size_gal` Optimal TES capacity, by volume [gal]
- `year_one_soc_series_pct` Vector of normalized (0-1) state of charge values over the first year [-]
- `year_one_to_load_series_mmbtu_per_hour` Vector of power used to meet load over the first year [MMBTU/hr]
"""
function add_hot_storage_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict, b::String; _n="")
    delta_T_degF = p.s.storage.attr["HotThermalStorage"].hot_water_temp_degF - p.s.storage.attr["HotThermalStorage"].cool_water_temp_degF
    avg_cp_kj_per_kgK = 998.2
    avg_rho_kg_per_m3 = 4.184 #TODO: add CoolProp reference or perform analogous calculations for water and build lookup tables
    kwh_per_gal = convert_gal_to_kwh(delta_T_degF, avg_rho_kg_per_m3, avg_cp_kj_per_kgK)

    r = Dict{String, Any}()
    size_kwh = round(value(m[Symbol("dvStorageEnergy"*_n)][b]), digits=2)
    r["size_gal"] = size_kwh / kwh_per_gal

    if size_kwh != 0
    	soc = (m[Symbol("dvStoredEnergy"*_n)][b, ts] for ts in p.time_steps)
        r["year_one_soc_series_pct"] = round.(value.(soc) ./ size_kwh, digits=3)

        discharge = (m[Symbol("dvDischargeFromStorage"*_n)][b, ts] for ts in p.time_steps)
        r["year_one_to_load_series_mmbtu_per_hour"] = round.(value.(discharge) / KWH_PER_MMBTU, digits=7)
    else
        r["year_one_soc_series_pct"] = []
        r["year_one_to_load_series_mmbtu_per_hour"] = []
    end

    d[b] = r
    nothing
end

"""
    add_cold_storage_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")

Adds the Storage results to the dictionary passed back from `run_reopt` using the solved model `m` and the `REoptInputs` for node `_n`.
Note: the node number is an empty string if evaluating a single `Site`.

Storage results:
- `size_gal` Optimal TES capacity, by volume [gal]
- `year_one_soc_series_pct` Vector of normalized (0-1) state of charge values over the first year [-]
- `year_one_to_load_series_ton` Vector of power used to meet load over the first year [ton]
"""
function add_cold_storage_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict, b::String; _n="")
    delta_T_degF = p.s.storage.attr["ColdThermalStorage"].hot_water_temp_degF - p.s.storage.attr["ColdThermalStorage"].cool_water_temp_degF
    avg_cp_kj_per_kgK = 998.2
    avg_rho_kg_per_m3 = 4.184 #TODO: add CoolProp reference or perform analogous calculations for water and build lookup tables
    kwh_per_gal = convert_gal_to_kwh(delta_T_degF, avg_rho_kg_per_m3, avg_cp_kj_per_kgK)

    r = Dict{String, Any}()
    size_kwh = round(value(m[Symbol("dvStorageEnergy"*_n)][b]), digits=2)
    r["size_gal"] = size_kwh / kwh_per_gal

    if size_kwh != 0
    	soc = (m[Symbol("dvStoredEnergy"*_n)][b, ts] for ts in p.time_steps)
        r["year_one_soc_series_pct"] = round.(value.(soc) ./ size_kwh, digits=3)

        discharge = (m[Symbol("dvDischargeFromStorage"*_n)][b, ts] for ts in p.time_steps)
        r["year_one_to_load_series_ton"] = round.(value.(discharge) / KWH_THERMAL_PER_TONHOUR, digits=7)
    else
        r["year_one_soc_series_pct"] = []
        r["year_one_to_load_series_ton"] = []
    end

    d[b] = r
    nothing
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
 # Test quadratic module when mindegree of polynomial is strictly bigger than 0

@testset "SOSquartic with $(factory.optimizer_constructor)" for factory in sdp_factories
    @polyvar x y
    m = SOSModel(factory)
    p = x^4 - y^4
    S = @set - x^2 - y^4 >=0
    soscon = @constraint m p >= 0 domain = S

    JuMP.optimize!(m)

    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    Q = gram_matrix(soscon)
    @test mindegree(monomials(Q)) == 2
    @test maxdegree(monomials(Q)) == 4
end

ext: jl


text:
# error codes
const UNKNOWN = -1
const NONE = 0
const OFFSET_OUT_OF_RANGE = 1
const CORRUPT_MESSAGE = 2
const UNKNOWN_TOPIC_OR_PARTITION = 3
const LEADER_NOT_AVAILABLE = 5
const NOT_LEADER_FOR_PARTITION = 6
const REQUEST_TIMED_OUT = 7
const BROKER_NOT_AVAILABLE = 8
const REPLICA_NOT_AVAILABLE = 9
const MESSAGE_TOO_LARGE = 10
const STALE_CONTROLLER_EPOCH = 11
const OFFSET_METADATA_TOO_LARGE = 12
const NETWORK_EXCEPTION = 13
const GROUP_LOAD_IN_PROGRESS = 14
const GROUP_COORDINATOR_NOT_AVAILABLE = 15
const NOT_COORDINATOR_FOR_GROUP = 16
const INVALID_TOPIC_EXCEPTION = 17
const RECORD_LIST_TOO_LARGE = 18
const NOT_ENOUGH_REPLICAS = 19
const NOT_ENOUGH_REPLICAS_AFTER_APPEND = 20
const INVALID_REQUIRED_ACKS = 21
const ILLEGAL_GENERATION = 22
const INCONSISTENT_GROUP_PROTOCOL = 23
const INVALID_GROUP_ID = 24
const UNKNOWN_MEMBER_ID = 25
const INVALID_SESSION_TIMEOUT = 26
const REBALANCE_IN_PROGRESS = 27
const INVALID_COMMIT_OFFSET_SIZE = 28
const TOPIC_AUTHORIZATION_FAILED = 29
const GROUP_AUTHORIZATION_FAILED = 30
const CLUSTER_AUTHORIZATION_FAILED = 31

# API keys (part)
const PRODUCE = 0
const FETCH = 1
const OFFSETS = 2
const METADATA = 3
const API_VERSIONS = 18

# Kafka.jl specific
const DEFAULT_ID = "kafka-jl"


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
Pkg.update()

using JuliaFormatter

main_paths = ["./src", "./test", "./script"]
for main_path in main_paths
    for (root, dir, files) in walkdir(main_path)
        for f in files
        @show file_path = abspath(root, f)
        !occursin(".jl", f) && continue
        occursin("generated", file_path) && continue
        format(file_path;
            whitespace_ops_in_indices = true,
            remove_extra_newlines = true,
            verbose = true,
            always_for_in = true,
            whitespace_typedefs = true,
            # always_use_return = true. # Disabled since it throws a lot of false positives
            )
        end
    end
end

max_stars_repo_name: NREL-SIIP/Examples

max_issues_repo_name: NREL-SIIP/Examples

max_forks_repo_name: NREL-SIIP/Examples


text:
prog = "1 + 1"
ex = Meta.parse(prog)
println(ex)
println(ex.head)
println(ex.args)

ex2 = Expr(:call,:+,1,1)
println(ex2 == ex)
dump(ex2)

max_stars_repo_name: zqsheng/snippet

max_issues_repo_name: zqsheng/snippet

max_forks_repo_name: zqsheng/snippet


text:
# this many digits[bits] (are given by) that many bits[digits]
# sum([abs(i-maximin_digits(maximin_bits(i))) for i=24:4000])

const log10of2 = log10(2)
const log2of10 = log2(10)

@inline digits2bits(x::Int) = floor(Int, log2of10*x)
@inline bits2digits(x::Int) = ceil(Int, log10of2*x)

@inline maximin_digits(nbits) = floor(Int, log10of2*nbits + 0.125)   # match on round trip
@inline maximin_bits(ndigits) = floor(Int, log2of10*ndigits + 2.625) # at least enough, at most 3 more on round trip

@inline function digits4bits(nbits)
    bits = log10of2 * nbits
    bits += 0.125
    floor(Int, bits)
end

@inline function bits4digits(ndigits)
    digs = log2of10 * ndigits
    digs += 2.625
    floor(Int, digs)
end

#     working_precision exceeds evinced_precision

# additional precision (more of the significant bits)
#    absorbs some kinds of numerical jitter
#    before any undesired resonance occurs
const BitsOfStability  =  9

# additional accuracy (more bits of the significand)
#    may compensate for 1,2,3 ulp enclosure widenings
const BitsOfAbsorption = 15

const ExtraBits = Ref(BitsOfStability + BitsOfAbsorption)

@inline workingbits(evincedbits) = evincedbits + ExtraBits.x
@inline evincedbits(workingbits) = workingbits - ExtraBits.x

# preset precisions
const MINIMUM_PRECISION_BASE2 = 24
const MINIMUM_PRECISION_BASE10 = 8
const MINIMUM_PRECISION = MINIMUM_PRECISION_BASE2

const DEFAULT_PRECISION = Ref(workingbits(128 - ExtraBits.x))

# these typed significands have this many signficant bits

workingprecision(::Type{Mag}) = 30 # bits of significand

workingprecision(::Type{ArbFloat}) = DEFAULT_PRECISION.x
workingprecision(::Type{ArbReal}) = DEFAULT_PRECISION.x
workingprecision(::Type{ArbComplex}) = DEFAULT_PRECISION.x

workingprecision(::Type{ArbFloat{P}}) where {P} = P
workingprecision(::Type{ArbReal{P}}) where {P} = P
workingprecision(::Type{ArbComplex{P}}) where {P} = P

workingprecision(x::ArbFloat{P}) where {P} = P
workingprecision(x::ArbReal{P}) where {P} = P
workingprecision(x::ArbComplex{P}) where {P} = P

# these typed significands have this many signficant bits shown

precision(::Type{Mag}) = 30 # bits of significand
precision(::Type{ArbFloat}) = evincedbits(DEFAULT_PRECISION.x)
precision(::Type{ArbReal}) = evincedbits(DEFAULT_PRECISION.x)
precision(::Type{ArbComplex}) = evincedbits(DEFAULT_PRECISION.x)

precision(::Type{ArbFloat{P}}) where {P} = evincedbits(P)
precision(::Type{ArbReal{P}}) where {P} = evincedbits(P)
precision(::Type{ArbComplex{P}}) where {P} = evincedbits(P)

precision(x::ArbFloat{P}) where {P} = evincedbits(P)
precision(x::ArbReal{P}) where {P} = evincedbits(P)
precision(x::ArbComplex{P}) where {P} = evincedbits(P)

function setprecision(::Type{T}; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:ArbNumber}
    if base === 10
        digits = digits > 0 ? bits4digits(digits) : (bits > 0 ? bits : DEFAULT_PRECISION.x)
    elseif base === 2
        digits = bits > 0 ? bits : (digits > 0 ? digits : DEFAULT_PRECISION.x)
    else
        throw(ErrorException("base expects 2 or 10"))
    end
    setprecision(T, digits)
end

function setprecision(::Type{T}, n::Int) where {T<:ArbNumber}
    n < MINIMUM_PRECISION && throw(DomainError("bit precision must be >= $MINIMUM_PRECISION"))
    DEFAULT_PRECISION.x = workingbits(n)
    return n
end

function setworkingprecision(::Type{T}; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:ArbNumber}
    if base === 10
        digits = digits > 0 ? bits4digits(digits) : (bits > 0 ? bits : DEFAULT_PRECISION.x)
    elseif base === 2
        digits = bits > 0 ? bits : (digits > 0 ? digits : DEFAULT_PRECISION.x)
    else
        throw(ErrorException("base expects 2 or 10"))
    end
    setworkingprecision(T, digits)
end

function setworkingprecision(::Type{T}, n::Int) where {T<:ArbNumber}
    n < workingbits(MINIMUM_PRECISION) && throw(DomainError("working bit precision must be >= $(workingbits(MINIMUM_PRECISION))"))
    DEFAULT_PRECISION.x = n
    return n
end

extrabits() = ExtraBits.x

extrabits(::Type{ArbFloat}) = ExtraBits.x
extrabits(::Type{ArbReal}) = ExtraBits.x
extrabits(::Type{ArbComplex}) = ExtraBits.x

extrabits(::Type{ArbFloat{P}}) where {P} = extrabits(ArbFloat)
extrabits(::Type{ArbReal{P}}) where {P} = extrabits(ArbReal)
extrabits(::Type{ArbComplex{P}}) where {P} = extrabits(ArbComplex)

extrabits(x::ArbFloat{P}) where {P} = extrabits(ArbFloat)
extrabits(x::ArbReal{P}) where {P} = extrabits(ArbReal)
extrabits(x::ArbComplex{P}) where {P} = extrabits(ArbComplex)

function setextrabits(n::Int)
    n >=0 || throw(DomainError("extrabits must be >= 0"))
    priorextra = extrabits()
    priorworking = workingprecision(ArbFloat)
    newextra = n
    newworking = priorworking + (newextra - priorextra)
    ExtraBits.x = newextra
    DEFAULT_PRECISION.x = newworking
    return ExtraBits.x
end

setextrabits(::Type{T}, n::Int) where {T<:ArbNumber} = setextrabits(n)


text:
# Linear Quadratic Regulator Objective Functions

using SparseArrays

include("QP_Linear_objectives.jl")

abstract type objectiveLQR_abstract <: objectiveFunc end

"""
    LQR_simple(Q, R)

A simple LQR struct where `Q` and `R` are constant for all `x` and `u`.
"""
struct LQR_simple <: objectiveLQR_abstract
    Q
    R
end

"""
    LQR_QP(QR_Full)

An LQR what acts like a QP objective function.
"""
struct LQR_QP <: objectiveLQR_abstract
    QR_Full
end

"""
    LQR_QP_Referenced(lqr_qp::LQR_QP, xRef)

An LQR what acts like a QP objective function, but with a reference
trajectory.
"""
struct LQR_QP_Referenced <: objectiveLQR_abstract
    lqr_qp::LQR_QP
    xRef
end

"""
    LQR_QP_Ref_Terminal(lqr_qp::LQR_QP, xRef)

An LQR what acts like a QP objective function, but with a reference
trajectory.
"""
struct LQR_QP_Ref_Terminal <: objectiveLQR_abstract
    lqr_qp::LQR_QP
    xRef
end

#####################################
#  LQR -> ObjectiveQP Constructors  #
#####################################
@doc raw"""
    makeLQR_TrajSimple(Q, R, NSteps::Int64)

Pseudo-constructor for an LQR objective function (which is really an
[`objectiveQ`](@ref) objective function).

An LQR struct for a trajectory is one that holds the sparse matrix matching:

```math
B_{QR} = [Q 0 0 0 0; 0 R 0 0 0; 0 0 Q 0 0; 0 0 0 R 0; 0 0 0 0 Q]
```

Where size corresponds to the size of the trajectory (or trajectory horizon).

```math
f(y) = \frac{1}{2} \ y^{\top} * B_{QR} * y
```

See also: [`LQR_simple`](@ref) and [`LQR_QP`](@ref)
"""
function makeLQR_TrajSimple(Q, R, NSteps::Int64)
    QSize = size(Q, 1)
    RSize = size(R, 1)
    totSize = (QSize + RSize) * NSteps + QSize

    # Initialize and empty sparse array
    QRFull = spzeros(totSize, totSize)

    rStart = 1

    for k in 1:NSteps
        rEnd = rStart + QSize - 1
        QRFull[rStart:rEnd, rStart:rEnd] = Q

        rStart += QSize
        rEnd = rStart + RSize - 1

        QRFull[rStart:rEnd, rStart:rEnd] = R

        rStart += RSize
    end
    rEnd = rStart + QSize - 1
    QRFull[rStart:rEnd, rStart:rEnd] = Q

    return LQR_QP(QRFull)
end

function makeLQR_TrajSimple(lqr::LQR_simple, NSteps::Int64)
    return makeLQR_TrajSimple(lqr.Q, lqr.R, NSteps)
end

function makeLQR_TrajReferenced(Q, R, NSteps::Int64, xRef)
    baseLQR = makeLQR_TrajSimple(Q, R, NSteps)
    return LQR_QP_Referenced(baseLQR, xRef)
end

function makeLQR_TrajReferenced(lqr::LQR_simple, NSteps::Int64, xRef)
    baseLQR = makeLQR_TrajSimple(lqr.Q, lqr.R, NSteps)
    return LQR_QP_Referenced(baseLQR, xRef)
end

@doc raw"""
    makeLQR_TrajTerminal(Q, R, NSteps::Int64)

Pseudo-constructor for an LQR objective function (which is really an
[`objectiveQ`](@ref) objective function).

An LQR struct for a trajectory is one that holds the sparse matrix matching:

```math
B_{QR} = [0 0 0 0 0; 0 R 0 0 0; 0 0 0 0 0; 0 0 0 R 0; 0 0 0 0 Q]
```

Where size corresponds to the size of the trajectory (or trajectory horizon).

```math
f(y) = \frac{1}{2} \ y^{\top} * B_{QR} * y
```

See also: [`LQR_simple`](@ref) and [`LQR_QP`](@ref)
"""
function makeLQR_TrajTerminal(Q, R, NSteps::Int64)
    QSize = size(Q, 1)
    RSize = size(R, 1)
    totSize = (QSize + RSize) * NSteps + QSize

    # Initialize and empty sparse array
    QRFull = spzeros(totSize, totSize)

    rStart = 1 + QSize

    for k in 1:NSteps
        rEnd = rStart + RSize - 1

        QRFull[rStart:rEnd, rStart:rEnd] = R

        rStart += RSize + QSize
    end
    rStart -= QSize
    rEnd = rStart + QSize - 1
    QRFull[rStart:rEnd, rStart:rEnd] = Q

    return LQR_QP(QRFull)
end

function makeLQR_TrajRefTerminal(Q, R, NSteps::Int64, xRef)
    baseLQR = makeLQR_TrajTerminal(Q, R, NSteps)
    return LQR_QP_Referenced(baseLQR, xRef)
end

function makeLQR_TrajRefTerminal(lqr::LQR_simple, NSteps::Int64, xRef)
    baseLQR = makeLQR_TrajTerminal(lqr.Q, lqr.R, NSteps)
    return LQR_QP_Referenced(baseLQR, xRef)
end

#######################################
###       Evaluate An LQR QP        ###
#######################################

"""
    fObjQP(lqr::LQR_QP, x)

Evaluates a full QP-like LQR. See [`LQR_QP`](@ref).
"""
function fObjQP(lqr::LQR_QP, x)
    return (1/2) * x' * (lqr.QR_Full) * x
end

"""
    fObjQP(lqr::LQR_QP, x, xRef)

Evaluates a full QP-like LQR relative to a reference. (i.e. x -> (x - xRef))

See [`LQR_QP`](@ref).
"""
function fObjQP(lqr::LQR_QP, x, xRef)
    xDiff = x - xRef
    return fObjQP(lqr, xDiff)
end

"""
    fObjQP(lqrRef::LQR_QP_Referenced, x)

Evaluates a full QP-like LQR relative to a reference. (i.e. x -> (x - xRef))

See [`LQR_QP`](@ref), [`LQR_QP_Referenced`](@ref)
"""
function fObjQP(lqrRef::LQR_QP_Referenced, x)
    # println("x = $(size(x))")
    # println("xRef = $(size(lqrRef.xRef))")
    xDiff = x - lqrRef.xRef
    return fObjQP(lqrRef.lqr_qp, xDiff)
end

#######################################
###      Gradients of LQR QP        ###
#######################################
"""
    dfdxQP(lqr::LQR_QP, x)

Evaluates the derivative of a full QP-like LQR at input `x`.

See [`LQR_QP`](@ref), [`dfdxQP`](@ref)
"""
function dfdxQP(lqr::LQR_QP, x)
    return (lqr.QR_Full) * x
end

"""
    dfdxQP(lqr::LQR_QP, x, xRef)

Evaluates the derivative of a full QP-like LQR at input `x` with respect to
a reference `xRef`.

See [`LQR_QP`](@ref), [`dfdxQP`](@ref)
"""
function dfdxQP(lqr::LQR_QP, x, xRef)
    xDiff = x - xRef
    return dfdxQP(lqr, xDiff)
end

"""
    dfdxQP(lqrRef::LQR_QP_Referenced, x)

Evaluates the derivative of a full QP-like LQR at input `x` with respect to
a reference `xRef`.

See [`LQR_QP_Referenced`](@ref), [`LQR_QP`](@ref), [`dfdxQP`](@ref)
"""
function dfdxQP(lqrRef::LQR_QP_Referenced, x)
    return dfdxQP(lqrRef.lqr_qp, x, lqrRef.xRef)
end

#######################################
###      Hessians of LQR QP         ###
#######################################
"""
    hessQP(lqr::LQR_QP)

Evaluates the hessian of a full QP-like LQR, which is simple the base matrix.

See [`LQR_QP`](@ref)
"""
function hessQP(lqr::LQR_QP)
    return lqr.QR_Full
end

"""
    hessQP(lqr::LQR_QP)

Evaluates the hessian of a full QP-like LQR, which is simple the base matrix.

See [`LQR_QP`](@ref)
"""
function hessQP(lqrRef::LQR_QP_Referenced)
    return hessQP(lqrRef.lqr_qp)
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
export FourRooms

mutable struct FourRooms <: AbstractGridWorld
    world::GridWorldBase{Tuple{Empty,Wall,Goal}}
    agent_pos::CartesianIndex{2}
    agent::Agent
    goal_reward::Float64
    reward::Float64
end

function FourRooms(;n=9, agent_start_pos=CartesianIndex(2,2), agent_start_dir=RIGHT, goal_pos=CartesianIndex(n-1, n-1))
    objects = (EMPTY, WALL, GOAL)
    world = GridWorldBase(objects, n, n)

    world[WALL, [1,n], 1:n] .= true
    world[WALL, 1:n, [1,n]] .= true
    world[WALL, ceil(Int,n/2), vcat(2:ceil(Int,n/4)-1,ceil(Int,n/4)+1:ceil(Int,n/2)-1,ceil(Int,n/2):ceil(Int,3*n/4)-1,ceil(Int,3*n/4)+1:n)] .= true
    world[WALL, vcat(2:ceil(Int,n/4)-1,ceil(Int,n/4)+1:ceil(Int,n/2)-1,ceil(Int,n/2):ceil(Int,3*n/4)-1,ceil(Int,3*n/4)+1:n), ceil(Int,n/2)] .= true
    world[EMPTY, :, :] .= .!world[WALL, :, :]
    world[GOAL, goal_pos] = true
    world[EMPTY, goal_pos] = false

    goal_reward = 1.0
    reward = 0.0

    env = FourRooms(world,agent_start_pos,Agent(dir=RIGHT), goal_reward, reward)

    reset!(env, agent_start_pos = agent_start_pos, agent_start_dir = agent_start_dir, goal_pos = goal_pos)

    return env
end

function (env::FourRooms)(::MoveForward)
    dir = get_dir(env.agent)
    dest = dir(env.agent_pos)
    env.reward = 0.0
    if !env.world[WALL, dest]
        env.agent_pos = dest
        if env.world[GOAL, env.agent_pos]
            env.reward = env.goal_reward
        end
    end
    env
end

RLBase.get_terminal(env::FourRooms) = env.world[GOAL, env.agent_pos]

function RLBase.reset!(env::FourRooms; agent_start_pos = CartesianIndex(2, 2), agent_start_dir = RIGHT, goal_pos = CartesianIndex(size(env.world[end]) - 1, size(env.world)[end] - 1))
    n = size(env.world)[end]
    env.reward = 0.0
    env.agent_pos = agent_start_pos
    agent = get_agent(env)
    set_dir!(agent, agent_start_dir)

    env.world[GOAL, :, :] .= false
    env.world[GOAL, goal_pos] = true
    env.world[EMPTY, :, :] .= .!env.world[WALL, :, :]
    env.world[EMPTY, goal_pos] = false
    return env
end


text:
function ExperienceReplay(;π,
                          S,
                          𝒫=(;),
                          A=action_space(π),
                          N_experience_replay,
                          solver=TD3,
                          required_columns=Symbol[],
                          c_opt=(;),
                          a_opt=(;),
                          replay_store_weight = (D)->1f0,
                          kwargs...)

    buffer_er = ExperienceBuffer(S, A, N_experience_replay, [:episode_end, :weight, :value, required_columns...])
    𝒫 = (;buffer_er, 𝒫...)
    function cb(D; 𝒮, info=Dict())
        D[:value] = mean(value(𝒮.agent.π, D[:s], D[:a]))
        D[:weight] .= replay_store_weight(D)
        push_reservoir!(buffer_er, D, weighted=true)
    end

    solver(;  π=π,
              S=S,
              𝒫=𝒫,
              required_columns=unique([:episode_end, :weight, required_columns...]),
              post_sample_callback=cb,
              c_opt=(;regularizer=BatchRegularizer(buffers=[buffer_er], batch_size=128, λ=0.5f0, loss=action_value_regularization, c_opt...)),
              a_opt=(;regularizer=BatchRegularizer(buffers=[buffer_er], batch_size=128, λ=0.5f0, loss=action_regularization, a_opt...)),
              extra_buffers=[buffer_er],
              buffer_fractions=[0.5, 0.5],
              kwargs...
              )
end


text:
#=
Problem 101
=#

function u_gen(k::Array, n)
    sum = 0
    for i = length(k):-1:1
        sum = k[i] + sum * n
    end
    return sum
end

N = 11

# first 11 terms
u_n = Vector{Float64}(N)
for i=1:N
    u_n[i] = u_gen([1 -1 1 -1 1 -1 1 -1 1 -1 1], i)
end

# Matrix of coefficients
A = Array{Float64}(N,N)
for i=1:N, j=1:N
    A[i,j] = i^(j-1)
end

sum = 0
for i=1:N-1
    coeff = Array{Int64}(round(inv(A[1:i, 1:i]) * u_n[1:i]))
    FIT = u_gen(coeff, i+1)
    @show coeff
    @show FIT
    sum += FIT
end

println("Sum of FITs is $(sum)")

# Alternatively u_n = A*[1 -1 1 -1 1 -1 1 -1 1 -1 1]'

max_stars_repo_name: mebrunet/project-euler

max_issues_repo_name: mebrunet/project-euler

max_forks_repo_name: mebrunet/project-euler


text:
#= Julia (v1) function to retrieve frequencies as well as Stefan Th. Gries' DP (Deviation of Proportions) word dispersion algorithm for words in a directory of TXT files.
Earl K. Brown, ekbrown byu edu (add appropriate characters to create email)
=#

using DataFrames

function get_dispersion(input_dir, min_freq = 1, sort_by = "dp")
    #= The function takes as input a string with the pathway to the input directory with .txt files and returns as output a Julia DataFrame with columns: word, freq, log freq, range, DP, and DP normalized.

    input_dir: A string with the pathway (whether absolute or relative) to the directory with the .txt files; other file types are ignored.
    min_freq: An integer indicating the minimum frequency for words to include (default is 1).
    sort_by: A string indicating how the results are ordered, from among "freq", "range" or "dp" (default is "dp").
    return value: A Julia DataFrame.

    For information about the DP dispersion algorithm, see:
    Gries, Stefan Th. 2008. Dispersions and adjusted frequencies in corpora. International Journal of Corpus Linguistics 13(4), 403–437.
    Lijffijt, Jeffrey & Stefan Th. Gries. 2012. Correction to "Dispersions and adjusted frequencies in corpora". International Journal of Corpus Linguistics 17(1), 147-149.
    =#

    # get filenames
    cd(input_dir)
    filenames = filter(f -> occursin(r"\.txt$"i, f), readdir())

    # create empty array to collect number of words in each file
    s = Array{Int64, 1}()

    # create dictionary to collect frequency of words in each file
    freq_by_file = Dict{Tuple, Int64}()

    # create dictionary to collect frequency of words across all files
    all_freq = Dict{String, Int64}()

    # counter for number of files
    counter = 0

    # loop over files in directory
    for f in filenames

        # increment file counter
        counter += 1

        # open connection to current file
        open(f) do infile

            # read whole file as string and make it uppercase
            whole_file = uppercase(read(infile, String))

            # get words in current file
            wds = split(whole_file, r"[^-'a-záéíóúüñ0-9]+"i, keepempty = false)

            # push number of words in current file to collector
            push!(s, length(wds))

            # loop over words in current file, pushing them to the collector dictionaries
            for wd in wds
                freq_by_file[(wd, counter)] = get(freq_by_file, (wd, counter), 0) + 1
                all_freq[wd] = get(all_freq, wd, 0) + 1
            end  # next word
        end  # close connection to current file
    end  # next file

    # get number of files in directory
    num_files = length(s)

    # convert number of words in each file to ratios
    s /= sum(s)

    # create a data frame collector to format output
    output = DataFrame(wd = String[], freq = Int64[], freq_log10 = Float64[], range = Int64[], dp = Float64[], dp_norm = Float64[])

    # loop over words and pushes DP of each word to collector
    for (k, v) in all_freq
        observed = Array{Int64, 1}()
        for j in 1:num_files
            push!(observed, get(freq_by_file, (k, j), 0))
        end
        cur_range = length(filter(x -> x > 0, observed))
        cur_dp = sum(abs.((observed / all_freq[k]) - s)) / 2
        cur_dp_norm = cur_dp / (1 - minimum(s))
        push!(output, [k, v, log10(v), cur_range, cur_dp, cur_dp_norm])
    end

    # filter by minimum frequency desired by user
    output = output[(output[!, :freq].>= min_freq),:]

    # sort data frame
    if sort_by == "freq"
        sort!(output, [order(:freq, rev = true), :wd])
    elseif sort_by == "range"
        sort!(output, [order(:range, rev = true), :wd])
    elseif sort_by == "dp"
        sort!(output, [order(:dp, rev = false), :wd])
    end

    return output

end  # end function definition

#########################
### test the function ###
#########################

input_dir = "/Users/ekb5/Corpora/USA/California/Salinas/transcripts"
min_freq = 1
sort_by = "freq"  # one of "freq", "range" or "dp"
@time results = get_dispersion(input_dir, min_freq, sort_by)  # runs function
println(results[1:10,:])


text:
abstract type AbstractBroadcastedNormal{T} <: ContinuousBatchDistribution end

function rand(rng::AbstractRNG, bd::AbstractBroadcastedNormal, dims::Int...)
    return bd.m .+ std(bd) .* randnsimilar(rng, bd.m, dims...)
end

function logpdf(bd::AbstractBroadcastedNormal{T}, x) where {T}
    diff = x .- bd.m
    return -(log(2 * T(pi)) .+ logvar(bd) .+ diff .* diff ./ var(bd)) / 2
end

mean(bd::AbstractBroadcastedNormal) = bd.m
mode(bd::AbstractBroadcastedNormal) = bd.m

_vlogv(bd::AbstractBroadcastedNormal) = (var(bd), logvar(bd))

function kldiv(bd1::AbstractBroadcastedNormal{T}, bd2::AbstractBroadcastedNormal{T}) where {T}
    diff = bd2.m .- bd1.m
    v1, logv1 = _vlogv(bd1)
    v2, logv2 = _vlogv(bd2)
    return (logv2 .- logv1 .- 1 .+ v1 ./ v2 .+ diff .* diff ./ v2) / 2
end

"""
Broadcasted Normal distribution with standard deviation.
"""
struct BroadcastedNormalStd{T<:Number,Tm<:AbstractArray{T},Ts<:AbstractArray{T}} <: AbstractBroadcastedNormal{T}
    m::Tm
    s::Ts
end

Broadcast.broadcastable(bd::BroadcastedNormalStd) = Ref(bd)

Normal(m::AbstractArray, s::AbstractArray) = BroadcastedNormalStd(m, s)
NormalStd(m::AbstractArray, s::AbstractArray) = BroadcastedNormalStd(m, s)

   std(bd::BroadcastedNormalStd) = bd.s
   var(bd::BroadcastedNormalStd) = bd.s.^2
logvar(bd::BroadcastedNormalStd) = 2log.(bd.s)

"""
Broadcasted Normal distribution with variance.
"""
struct BroadcastedNormalVar{T<:Number,Tm<:AbstractArray{T},Tv<:AbstractArray{T}} <: AbstractBroadcastedNormal{T}
    m::Tm
    v::Tv
end

Broadcast.broadcastable(bd::BroadcastedNormalVar) = Ref(bd)

NormalVar(m::AbstractArray, v::AbstractArray) = BroadcastedNormalVar(m, v)

   std(bd::BroadcastedNormalVar) = sqrt.(bd.v)
   var(bd::BroadcastedNormalVar) = bd.v
logvar(bd::BroadcastedNormalVar) = log.(bd.v)

"""
Broadcasted Normal distribution with log standard deviation.
"""
struct BroadcastedNormalLogStd{T<:Number,Tm<:AbstractArray{T},Ts<:AbstractArray{T}} <: AbstractBroadcastedNormal{T}
       m::Tm    # mean
    logs::Ts    # log std
end

Broadcast.broadcastable(bd::BroadcastedNormalLogStd) = Ref(bd)

NormalLogStd(m::AbstractArray, logs::AbstractArray) = BroadcastedNormalLogStd(m, logs)

   std(bd::BroadcastedNormalLogStd) = exp.(bd.logs)
   var(bd::BroadcastedNormalLogStd) = exp.(2bd.logs)
logvar(bd::BroadcastedNormalLogStd) = 2bd.logs

function _vlogv(bd::BroadcastedNormalLogStd)
    logv = logvar(bd)
    v = exp.(logv)
    return v, logv
end

"""
Broadcasted Normal distribution with log variance.
"""
struct BroadcastedNormalLogVar{T<:Number,Tm<:AbstractArray{T},Tv<:AbstractArray{T}} <: AbstractBroadcastedNormal{T}
       m::Tm    # mean
    logv::Tv    # log std
end

Broadcast.broadcastable(bd::BroadcastedNormalLogVar) = Ref(bd)

NormalLogVar(m::AbstractArray, logv::AbstractArray) = BroadcastedNormalLogVar(m, logv)

   std(bd::BroadcastedNormalLogVar) = exp.(bd.logv ./ 2)
   var(bd::BroadcastedNormalLogVar) = exp.(bd.logv)
logvar(bd::BroadcastedNormalLogVar) = bd.logv

function _vlogv(bd::BroadcastedNormalLogVar)
    logv = logvar(bd)
    v = exp.(logv)
    return v, logv
end

###

# using LinearAlgebra: det, tr, inv
# using Distributions: MvNormal

# """
#     kldiv(mvn1::MvNormal, mvn2::MvNormal)

# Compute ``KL(MvNormal_1||MvNormal_2)``.
# """
# function kldiv(mvn1::MvNormal, mvn2::MvNormal)
#     d = length(mvn1.μ)
#     diff = mvn2.μ .- mvn1.μ
#     Σ1 = Matrix(mvn1.Σ)
#     Σ2 = Matrix(mvn2.Σ)
#     return (log(det(Σ2)) - log(det(Σ1)) - d + tr(inv(Σ2) * Σ1) + diff' * inv(Σ2) * diff) ./ 2
# end


text:
#! /usr/bin/julia

# Rosetta Code, Twelve statements

function showflaggedbits{T<:BitArray{1}}(a::T, f::T)
    tf = map(x->x ? "T" : "F", a)
    flg = map(x->x ? "*" : " ", f)
    join(tf .* flg, " ")
end

const props = [s -> length(s) == 12,
               s -> sum(s[7:12]) == 3,
               s -> sum(s[2:2:end]) == 2,
               s -> !s[5] || (s[6] & s[7]),
               s -> !any(s[2:4]),
               s -> sum(s[1:2:end]) == 4,
               s -> s[2] $ s[3],
               s -> !s[7] || (s[5] & s[6]),
               s -> sum(s[1:6]) == 3,
               s -> s[11] & s[12],
               s -> sum(s[7:9]) == 1,
               s -> sum(s[1:end-1]) == 4]

const NDIG = length(props)
NDIG < WORD_SIZE || println("WARNING, too many propositions!")

mhist = zeros(Int, NDIG+1)

println("Checking the ", NDIG, " statements against all possibilities.\n")
print(" "^15)
for i in 1:NDIG
    print(@sprintf "%3d" i)
end
println()

for i in 0:(2^NDIG-1)
    s = bitpack(digits(i, 2, NDIG))
    t = bitpack([p(s) for p in props])
    misses = s$t
    mcnt = sum(misses)
    mhist[NDIG-mcnt+1] += 1
    mcnt < 2 || mcnt == NDIG || continue
    if mcnt == 0
        print("    Exact Match: ")
    elseif mcnt == NDIG
        print("     Total Miss: ")
    else
        print("      Near Miss: ")
    end
    println(showflaggedbits(t, misses))
end

println()
println("Distribution of matches")
println(" Matches  Cases")
for i in (NDIG+1):-1:1
    println(@sprintf "    %2d => %4d" i-1 mhist[i])
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
export FixedNeighborhood

struct FixedNeighborhood <: NeighborhoodAlgorithm
    k::Int
end

function FixedNeighborhood()
    return FixedNeighborhood(8)
end

function optimize_neighborhood!(algo::FixedNeighborhood, index::SearchGraph{T}, dist::Function, perf, recall) where T
end

function neighborhood(algo::FixedNeighborhood, index::SearchGraph{T}, dist::Function, item::T) where T
    nbuffer::Vector{Int32} = Vector{Int}(undef, 0)
    knn = search(index, dist, item, KnnResult(algo.k))
    visible = Set{Int32}()

    for p in knn
        push!(nbuffer, p.objID)
    end

    return knn, nbuffer
end


text:
#! /usr/bin/env julia

# ---------------------------------------------------------------------------- #
# Parsing of command line arguments

using ArgParse

function ArgParse.parse_item{R <: Real}(::Type{Tuple{R, R}}, x::AbstractString)
  try
    a = parse(R, x)
    return (a, a)
  catch ArgumentError
    return ((parse(R, y) for y in split(x, '-'))...) :: Tuple{R, R}
  end
end

s = ArgParseSettings()

@add_arg_table s begin
  "train"
    help = "train a model to fit synthesized images that imitate microscopy data"
    action = :command
  "test"
    help = "show an example image and test the model performance on it"
    action = :command
  "generate"
    help = "generate a synthetic image with given parameters"
    action = :command
end

@add_arg_table s["train"] begin
  "--model", "-m"
    help = "model to be used - one of 'unetlike', 'multiscale3', fcrna'"
    arg_type = String
    default = "multiscale3"
  "--res", "-r"
    help = "resolution of the quadratic image patches used for training / testing"
    arg_type = Int
    default = 256
  "--sample-size", "-s"
    help = "number of patches used for training"
    arg_type = Int
    default = 10
  "--test-size", "--testset-size", "-t"
    help = "number of patches used for testing"
    arg_type = Int
    default = 5
  "--batch-size", "-b"
    help = "number of patches per mini-batch during training"
    arg_type = Int
    default = 1
  "--batch-normalization", "-B"
    help = "whether to activate batch normalization"
    action = :store_true
  "--no-gpu"
    help = "prevent usage of gpu acceleration even if gpu support is detected"
    action = :store_true
  "--epochs", "-e"
    help = "number of training epochs"
    arg_type = Int
    default = 10
  "--learning-rate", "-l"
    help = "initial learning rate for the chosen optimizer [-o]"
    arg_type = Float64
    default = 5e-5
  "--shuffle", "-p"
    help = "shuffle patches in the training set for each epoch"
    action = :store_true
  "--model-path", "-M"
    help = "save the model in this file after training"
    arg_type = String
    default = ""
  "--jitter", "-j"
    help = "randomness of the annotated positions in the label"
    arg_type = Int
    default = 0
  "--cell-number", "-n"
    help = "range of allowed cell-numbers per patch"
    arg_type = Tuple{Int, Int}
    default = (50, 80)
  "--cell-size", "-a"
    help = "range of the size of the two main axis of a cell in pixels"
    arg_type = Tuple{Float64, Float64}
    default = (3., 15.)
  "--cell-intensity", "-c"
    help = "range of the intensity a single cell may have"
    arg_type = Tuple{Float64, Float64}
    default = (0.2, 4.0)
  "--soften", "-S"
    help = "filtersize of the softening post-processing on the patches"
    arg_type = Int
    default = 2
  "--noise", "-N"
    help = "magnitude of the applied random noise on the patches"
    arg_type = Float64
    default = 0.2
  "--proximity-kernel-size", "-k"
    help = "size of the peaks in the proximity map"
    arg_type = Int
    default = 7
  "--proximity-kernel-height", "-K"
    help = "height of the peaks in the proximity map"
    arg_type = Float64
    default = 100.
  "--optimizer", "-o"
    help = "optimizer used - one of 'adam', 'rmsprop', 'nesterov'"
    arg_type = String
    default = "adam"
end

@add_arg_table s["test"] begin
  "model"
    help = "path to the model that is to be used for testing"
  "--res", "-r"
    help = "resolution of the image patches to be tested"
    arg_type = Int
    default = 256
  "--cell-number", "-n"
    help = "range of allowed cell-numbers per patch"
    arg_type = Tuple{Int, Int}
    default = (50, 80)
  "--cell-size", "-a"
    help = "range of the size of the two main axis of a cell in pixels"
    arg_type = Tuple{Float64, Float64}
    default = (3., 15.)
  "--cell-intensity", "-c"
    help = "range of the intensity a single cell may have"
    arg_type = Tuple{Float64, Float64}
    default = (0.2, 4.0)
  "--soften", "-S"
    help = "filtersize of the softening post-processing on the patches"
    arg_type = Int
    default = 2
  "--noise", "-N"
    help = "magnitude of the applied random noise on the patches"
    arg_type = Float64
    default = 0.2
end

@add_arg_table s["generate"] begin
  "name"
    help = "filename of the generated image"
  "number"
    help = "number of image files to create"
    arg_type = Int
    default = 1
  "--res", "-r"
    help = "resolution of the image patches to be tested"
    arg_type = Int
    default = 256
  "--jitter", "-j"
    help = "randomness of the annotated positions in the label"
    arg_type = Int
    default = 0
  "--cell-number", "-n"
    help = "range of allowed cell-numbers per patch"
    arg_type = Tuple{Int, Int}
    default = (50, 80)
  "--cell-size", "-a"
    help = "range of the size of the two main axis of a cell in pixels"
    arg_type = Tuple{Float64, Float64}
    default = (3., 15.)
  "--cell-intensity", "-c"
    help = "range of the intensity a single cell may have"
    arg_type = Tuple{Float64, Float64}
    default = (0.2, 4.0)
  "--soften", "-S"
    help = "filtersize of the softening post-processing on the patches"
    arg_type = Int
    default = 2
  "--noise", "-N"
    help = "magnitude of the applied random noise on the patches"
    arg_type = Float64
    default = 0.2
end

args = parse_args(s)

# ---------------------------------------------------------------------------- #
# Parsing done. Now load DCellC and execute the program

using Knet

include("../src/DCellC.jl")
using DCellC

# ---------------------------------------------------------------------------- #
# First mode: Train a network, print the performance, and save it

if args["%COMMAND%"] == "train"

  args = args["train"]

  if args["model"] == "unetlike"
    model = UNetLike(GreyscaleImage, bn=args["batch-normalization"])
  elseif args["model"] == "multiscale3"
    model = Multiscale3(GreyscaleImage, bn=args["batch-normalization"])
  elseif args["model"] == "fcrna"
    model = FCRNA(GreyscaleImage, bn=args["batch-normalization"])
  else
    println(STDERR, "Invalid argument for option '--model' given. See --help")
  end

  res = args["res"]
  gen() = synthesize(args["res"], args["res"], args["cell-number"],
                     jitter = args["jitter"],
                     cell = SharpEllipticCell(args["cell-size"], args["cell-intensity"]),
                     pp = [PP.soften(args["soften"]), PP.noise(args["noise"])])

  train_data = [ gen() for i in 1:args["sample-size"] ]
  test_data = [ gen() for i in 1:args["test-size"] ]

  path = args["model-path"]
  if path == ""
    path = nothing
  end

  if args["optimizer"] == "adam"
    opt = Adam
  elseif args["optimizer"] == "rmsprop"
    opt = Rmsprop
  elseif args["optimizer"] == "nesterov"
    opt = Nesterov
  else
    println(STDERR, "Optimizer $(args["optimizer"]) not recognized")
  end

  if args["no-gpu"] || gpu() < 0
    println("# CPU mode")
    at = Array{Float32}
  else
    println("# GPU mode ($(gpu()))")
    at = KnetArray{Float32}
  end

  train!(model, train_data,
         epochs = args["epochs"],
         lr = args["learning-rate"],
         modelpath = path,
         batchsize = args["batch-size"],
         shuffle = args["shuffle"],
         testset = test_data,
         record = false,
         opt = opt,
         kernelsize = args["proximity-kernel-size"],
         peakheight = args["proximity-kernel-height"],
         at = at)

# ---------------------------------------------------------------------------- #
# Second mode: Create test image, test a given network on it,
# and display the result

elseif args["%COMMAND%"] == "test"
  import ImageView
  using Gtk.ShortNames

  args = args["test"]

  res = args["res"]
  gen() = synthesize(args["res"], args["res"], args["cell-number"],
                     cell = SharpEllipticCell(args["cell-size"], args["cell-intensity"]),
                     pp = [PP.soften(args["soften"]), PP.noise(args["noise"])])

  img, lbl = gen()

  lbls = [ lbl ]

  if args["model"] != nothing
    @printf "load model from %s...\n" args["model"]
    model = load(args["model"])
    dens  = densitymap(model, img)
    dlbl, cds = label(dens, candidates = true)
    adj  = adjacency(dlbl, lbl)

    @printf "calculate metrics...\n"
    @printf "loss:    %.3f\n" loss(model, img, lbl)
    @printf "count:   %d / %d (%.3f)\n" length(dlbl) length(lbl) (length(dlbl) / length(lbl))
    @printf "meanadj: %.3f\n" mean(adj)
    @printf "maxadj:  %.3f\n" maximum(adj)

    @printf "display test image\n"

    push!(lbls, dlbl)
    push!(lbls, cds)
    guidict = imshow(img, lbls)
    imshow(GreyscaleImage(proximitymap(args["res"], args["res"], lbl)))
    imshow(GreyscaleImage(dens))

  else
    @printf "no model provided - only display test image\n"
    guidict = imshow(img, lbls)
  end

  if (!isinteractive())
    c = Condition()
    win = guidict["gui"]["window"]
    signal_connect(win, :destroy) do widget
        notify(c)
    end
    wait(c)
  end

# ---------------------------------------------------------------------------- #
# Third mode: Generate synthetic images and save them

elseif args["%COMMAND%"] == "generate"

  args = args["generate"]

  res = args["res"]
  gen() = synthesize(args["res"], args["res"], args["cell-number"],
                     jitter = args["jitter"],
                     cell = SharpEllipticCell(args["cell-size"], args["cell-intensity"]),
                     pp = [PP.soften(args["soften"]), PP.noise(args["noise"])])

  # Create images
  number = args["number"]
  images = [ gen() for i in 1:number ]

  # Create files
  fname = args["name"]
  if length(images) == 1
    lmgsave(fname, images[1])
  elseif contains(fname, "%n")
    for i in 1:number
      name = replace(fname, "%n", (@sprintf "%03d" i))
      lmgsave(name, images[i])
    end
  else
    bname, ext = splitext(fname)
    for i in 1:number
      name = bname * (@sprintf "%03d" i) * "." * ext
      lmgsave(name, images[i])
    end
  end
end

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC


text:
"""
    build_logger(logging_stream, log_level)

Builds a logger using the stream `logging_stream`and `log_level` provided.

# Arguments
  * `logging_stream::String` is the output stream and can take the values:
  `"null"` logs to `/dev/null`, `"stdout"` (default) logs to standard output,
  `"/path/to/existing/file"` logs to an existing file and
  `"/path/to/non-existing/file"` creates the log file. If no valid option
  is provided, the default stream is the standard output.
  * `log_level::String` is the log level can take the values `"debug"`,
  `"info"`, `"error"` and defaults to `"info"` if no valid option is provided.
"""
function build_logger(logging_stream="stdout", log_level="info")
    # Quickly return null logger if the case
    logging_stream == "null" && return NullLogger()
    # Build stream logger
    logger_type = ConsoleLogger
    # Handle logging stream option
    if logging_stream == "stdout"
        _stream = stdout
    else
        try
            _stream = open(logging_stream, "a+")
        catch e
            @warn "Failed to create logfile $logging_stream.\n$e"
            _stream = stdout
        end
    end
    # Logging
    if log_level == "debug"
        _level = Logging.Debug
    elseif log_level == "info"
        _level = Logging.Info
    elseif log_level == "error"
        _level = Logging.Error
    else
        @warn "Wrong logging level, defaulting to $DEFAULT_LOG_LEVEL."
        _level = DEFAULT_LOG_LEVEL
    end
    # Output logger
    return logger_type(_stream,
                       _level,
                       meta_formatter=garamond_log_formatter,
                       show_limited=true,
                       right_justify=0)
end

"""
    garamond_log_formatter(level, _module, group, id, file, line)

Garamond -specific log message formatter. Takes a fixed set of input arguments
and returns the color, prefix and suffix for the log message.
"""
function garamond_log_formatter(level, _module, group, id, file, line)
    color = :normal
    _file = split(file, "/")[end]  # take only the file name, not full path
    prefix = "[$(Dates.format(Dates.now(), "yyyy-mm-yy HH:MM:SS"))]" *
             "[$(uppercase(string(level)))][$_file:$line]"
    suffix=""
    return color, prefix, suffix  # fixed signature
end


text:
"Generate k-nearest neighborhood graph edges with distances"
function find_nn(X::AbstractMatrix{T}, k::Int=12) where T<:Real
    m, n = size(X)
    r = Array{T}(undef, (n, n))
    d = Array{T}(undef, k, n)
    e = Array{Int}(undef, k, n)

    mul!(r, transpose(X), X)
    sa2 = sum(X.^2, dims=1)

    @inbounds for j = 1 : n
        @inbounds for i = 1 : j-1
             r[i,j] = r[j,i]
        end
        r[j,j] = 0
        @inbounds for i = j+1 : n
            v = sa2[i] + sa2[j] - 2 * r[i,j]
            r[i,j] = isnan(v) ? NaN : sqrt(max(v, 0.))
        end
        e[:, j] = sortperm(r[:,j])[2:k+1]
        d[:, j] = r[e[:, j],j]
    end

    return (d, e)
end

function adjmat(D::AbstractMatrix{T}, E::AbstractMatrix{<:Integer}) where {T<:Real}
    @assert size(D) == size(E) "Distance and edge matrix must be of the same size"
    n = size(D, 2)
    W = spzeros(T, n, n)
    @inbounds for II in CartesianIndices(D)
        d  = D[II]
        ii = E[II]
        j = II[2]
        W[ii, j] = d
        W[j, ii] = d
    end
    return W
end

"Crate a graph with largest connected component of adjacency matrix `W`"
function largest_component(G)
    CC = connected_components(G)
    if length(CC) > 1
        @warn "Found $(length(CC)) connected components. Largest component is selected."
        C = CC[argmax(map(length, CC))]
        G = first(induced_subgraph(G, C))
    else
        C = first(CC)
    end
    return G, C
end

"Generate a swiss roll dataset"
function swiss_roll(n::Int = 1000, noise::Float64=0.05)
    t = (3 * pi / 2) * (1 .+ 2 * rand(n, 1))
    height = 30 * rand(n, 1)
    X = [t .* cos.(t) height t .* sin.(t)] + noise * randn(n, 3)
    labels = vec(rem.(sum([round.(Int, t / 2) round.(Int, height / 12)], dims=2), 2))
    return collect(transpose(X)), labels
end

"Perform spectral decomposition for Ax=λI"
function decompose(M::AbstractMatrix{<:Real}, d::Int)
    W = isa(M, AbstractSparseArray) ? Symmetric(Matrix(M)) : Symmetric(M)
    F = eigen!(W)
    idx = sortperm(F.values)[2:d+1]
    return F.values[idx], F.vectors[:,idx]
end

"Perform spectral decomposition for Ax=λB"
function decompose(A::AbstractMatrix{<:Real}, B::AbstractMatrix{<:Real}, d::Int)
    AA = isa(A, AbstractSparseArray) ? Symmetric(full(A)) : Symmetric(A)
    BB = isa(B, AbstractSparseArray) ? Symmetric(full(B)) : Symmetric(B)
    F = eigen(AA, BB)
    idx = sortperm(F.values)[2:d+1]
    return F.values[idx], F.vectors[:,idx]
end


text:
using JLD2, NNlib, BenchmarkTools

# TODO organize and compare benchmarks using BenchmarkGroups

# We need things to go quickly here
BenchmarkTools.DEFAULT_PARAMETERS.samples = 20
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 2.5

results = Dict()

function add_result(val, keys...)
    r = results
    for k in keys[1:end-1]
        if !haskey(r, k)
            r[k] = Dict()
        end
        r = r[k]
    end
    r[keys[end]] = val
    return r
end

# Modify these as needed
for rank in (2,),
    N in (20, 40, 80),
    C_in in (1,),
    C_out in (1,),
    K in (3,),
    stride in (1,),
    dilation in (1,),
    padding in (0, 2)

    benchmark_items = [
            (NNlib.conv_direct!, NNlib.∇conv_data_direct!, NNlib.∇conv_filter_direct!, DenseConvDims, "direct"),
            (NNlib.conv_im2col!, NNlib.∇conv_data_im2col!, NNlib.∇conv_filter_im2col!, DenseConvDims, "im2col"),
            (NNlib.depthwiseconv_direct!, NNlib.∇depthwiseconv_data_direct!, NNlib.∇depthwiseconv_filter_direct!, DepthwiseConvDims, "direct"),
            (NNlib.depthwiseconv_im2col!, NNlib.∇depthwiseconv_data_im2col!, NNlib.∇depthwiseconv_filter_im2col!, DepthwiseConvDims, "im2col"),
    ]

    if NNlib.is_nnpack_available()
        push!(benchmark_items, (NNlib.conv_nnpack!, NNlib.∇conv_data_nnpack!, NNlib.∇conv_filter_nnpack!, DenseConvDims, "nnpack"))
    end

    for (conv!, ∇conv_data!, ∇conv_filter!, cT, backend) in benchmark_items

        x = zeros(Float32, repeat([N], rank)..., C_in, 1)
        if cT == DenseConvDims
            w = zeros(Float32, repeat([K], rank)..., C_in, C_out)
        else
            w = zeros(Float32, repeat([K], rank)..., C_out, C_in)
        end
        cdims = try
            cT(x, w; stride=stride, dilation=dilation, padding=padding)
        catch
            continue
        end

        if cT == DenseConvDims
            y = zeros(Float32, NNlib.output_size(cdims)..., C_out, 1)
        else
            y = zeros(Float32, NNlib.output_size(cdims)..., C_out*C_in, 1)
        end

        dx = similar(x)
        dw = similar(w)
        dy = similar(y)

        t_fwd = @benchmark $(conv!)($y, $x, $w, $cdims)
        t_dx = @benchmark $(∇conv_data!)($dx, $y, $w, $cdims)
        t_dw = @benchmark $(∇conv_filter!)($dw, $x, $y, $cdims)

        add_result(t_fwd, "conv$(rank)d", backend, cdims)
        add_result(t_dx, "conv$(rank)d_data", backend, cdims)
        add_result(t_dw, "conv$(rank)d_filter", backend, cdims)

        @show(cdims)
        @save "results.jld2" results
    end
end

# Modify these as needed
for rank in (2,),
    N in (20,),
    K in (2, 4),
    stride in (1, 2, 4)

    x = zeros(Float32, repeat([N], rank)..., 1, 1)
    pdims = PoolDims(x, K; stride=stride)
    y = zeros(Float32, NNlib.output_size(pdims)..., 1, 1)
    dx = similar(x)

    for (pool, ∇pool, name) in (
            (NNlib.maxpool!, NNlib.∇maxpool!, "maxpool"),
            (NNlib.meanpool!, NNlib.∇meanpool!, "meanpool"),
        )

        t_fwd  = @benchmark $(pool)( $y, $x, $pdims)
        t_data = @benchmark $(∇pool)($dx, $y, $y, $x, $pdims)

        add_result(t_fwd, "$(name)$(rank)d", "direct", pdims)
        add_result(t_data, "$(name)$(rank)d_data", "direct", pdims)

        @show(pdims)
        @save "results.jld2" results
    end

    if NNlib.is_nnpack_available()
        if NNlib.nnpack_supported_operation(pdims)
            t_fwd  = @benchmark NNlib.maxpool_nnpack!($y, $x, $pdims)

            add_result(t_fwd, "maxpool2d", "nnpack", pdims)

            @show(pdims)
            @save "results.jld2" results
        end
    end
end


text:
using ThreeBodyTB
using Plots

#setup chosen backend
#pyplot()
#gr()

# load tb_crys from file

TESTDIR = ThreeBodyTB.TESTDIR

filname = "$TESTDIR/data_forces/znse.in_vnscf_vol_2/projham.xml.gz"
tbc = read_tb_crys(filname)

println(tbc)

#these are all legal plotting commands

#plot_bandstr(tbc) #default plotting

#plot_bandstr(tbc, efermi=-0.3, align="vbm", proj_types=["Mg"]) #projection align
#plot_bandstr(tbc, efermi=-0.3, align="min") #align
#plot_bandstr(tbc, efermi=-0.3, align="fermi") #align

kpath=[0.0 0.0 0.0;0.5 0 0; 0.5 0.25 0.75; 0.0 0.0 0.0]
knames=["Γ", "X", "W", "Γ"]

plot_bandstr(tbc, kpath=kpath, names=knames, efermi=-0.3, align="vbm", proj_types=["Mg"], proj_orbs=[:s])

savefig("mgs_proj.pdf")

#compare two tb_crys
plot(reuse=false)
plot_compare_tb(tbc,tbc, efermi=-0.3, align="fermi") # compare to self

savefig("mgs_compare_self.pdf")


text:
"""
    from_toml(::Type{T}, filename::String; kw...) where T

Convert a given TOML file `filename` to an option type `T`. Valid fields
can be override by keyword arguments. See also [`from_dict`](@ref).
"""
function from_toml(::Type{T}, filename::String; kw...) where T
    is_option(T) || error("not an option type")
    d = TOML.parsefile(filename)
    d["#filename#"] = filename
    return from_dict(T, d; kw...)
end

"""
    from_toml_if_exists(::Type{T}, filename::String; kw...) where T

Similar to [`from_toml`](@ref) but will create the option instance
via `from_kwargs(T;kw...)` instead of error if the file does not exist.
"""
function from_toml_if_exists(::Type{T}, filename::String; kw...) where T
    if isfile(filename)
        return from_toml(T, filename; kw...)
    else
        return from_kwargs(T; kw...)
    end
end


text:
using BumpStdlibs
using Documenter

makedocs(;
    modules=[BumpStdlibs],
    authors="Dilum Aluthge and contributors",
    repo="https://github.com/JuliaLang/BumpStdlibs.jl/blob/{commit}{path}#L{line}",
    sitename="BumpStdlibs.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaLang.github.io/BumpStdlibs.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Usage" => "usage.md",
        "API" => "api.md",
    ],
    strict=true,
    linkcheck=true,
)

deploydocs(;
    repo="github.com/JuliaLang/BumpStdlibs.jl",
)


text:
#=
--------------------------------------------------------------------------------------------------------------------------------
Base and generic math
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    verb_message(verb::Bool, message::String)

@info `message` if `verb` is true.
"""
verb_message(verb::Bool, message::String) = verb ? @info(message) : nothing;

"""
    check_bounds(X::Real, LB::Real, UB::Real)

Check whether `X` is larger or equal than `LB` and lower or equal than `UB`
"""
check_bounds(X::Real, LB::Real, UB::Real) = X < LB || X > UB ? throw(DomainError) : nothing

"""
    check_bounds(X::Real, LB::Real)

Check whether `X` is larger or equal than `LB`
"""
check_bounds(X::Real, LB::Real) = X < LB ? throw(DomainError) : nothing

"""
    error_info(err::Exception)
    error_info(err::RemoteException)

Return error main information
"""
error_info(err::Exception) = (err, err.msg, stacktrace(catch_backtrace()));
error_info(err::RemoteException) = (err.captured.ex, err.captured.ex.msg, [err.captured.processed_bt[i][1] for i=1:length(err.captured.processed_bt)]);

"""
    nan_to_missing!(X::JArray{Float64})

Replace NaN with missing in `X`.
"""
function nan_to_missing!(X::JArray{Float64})
    X[isnan.(X) .=== true] .= missing;
end

"""
    trimmed_mean(X::AbstractArray{Float64,1}, α::Float64)

Compute the trimmed mean of `X` (i.e., the sample average of `X` having removed its `α` smallest and largest values).
"""
function trimmed_mean(X::AbstractArray{Float64,1}, α::Float64)

    # Check bounds for α
    check_bounds(α, 0, 0.5);

    # Compute trimmed mean
    trimmed_sample = @view X[quantile(X, α) .<= X .<= quantile(X, 1-α)];
    return mean(trimmed_sample);
end

"""
    sum_skipmissing(X::AbstractVector{Float64})
    sum_skipmissing(X::AbstractVector{Union{Missing, Float64}})

Compute the sum of the observed values in `X`.

# Examples
```jldoctest
julia> sum_skipmissing([1.0; missing; 3.0])
4.0
```
"""
sum_skipmissing(X::AbstractVector{Float64}) = sum(X);
sum_skipmissing(X::AbstractVector{Union{Missing, Float64}}) = sum(skipmissing(X));

"""
    sum_skipmissing(X::AbstractMatrix{Float64})
    sum_skipmissing(X::AbstractMatrix{Union{Missing, Float64}})

Compute the sum of the observed values in `X` column wise.

# Examples
```jldoctest
julia> sum_skipmissing([1.0 2.0; missing 3.0; 3.0 5.0])
3-element Array{Float64,1}:
 3.0
 3.0
 8.0
```
"""
sum_skipmissing(X::AbstractMatrix{Float64}) = sum(X, dims=2);

function sum_skipmissing(X::AbstractMatrix{Union{Missing, Float64}})
    n = size(X,1);
    output = zeros(n,1);
    for i=1:n
        Xi = @view X[i,:];
        output[i] = sum(skipmissing(Xi));
    end

    return output;
end

"""
    mean_skipmissing(X::AbstractVector{Float64})
    mean_skipmissing(X::AbstractVector{Union{Missing, Float64}})

Compute the mean of the observed values in `X`.

# Examples
```jldoctest
julia> mean_skipmissing([1.0; missing; 3.0])
2.0
```
"""
mean_skipmissing(X::AbstractVector{Float64}) = mean(X);
mean_skipmissing(X::AbstractVector{Union{Missing, Float64}}) = mean(skipmissing(X));

"""
    mean_skipmissing(X::AbstractMatrix{Float64})
    mean_skipmissing(X::AbstractMatrix{Union{Missing, Float64}})

Compute the mean of the observed values in `X` column wise.

# Examples
```jldoctest
julia> mean_skipmissing([1.0 2.0; missing 3.0; 3.0 5.0])
3-element Array{Float64,1}:
 1.5
 3.0
 4.0
```
"""
mean_skipmissing(X::AbstractMatrix{Float64}) = mean(X, dims=2);

function mean_skipmissing(X::AbstractMatrix{Union{Missing, Float64}})
    n = size(X,1);
    output = zeros(n,1);
    for i=1:n
        Xi = @view X[i,:];
        output[i] = mean(skipmissing(Xi));
    end

    return output;
end

"""
    std_skipmissing(X::AbstractVector{Float64})
    std_skipmissing(X::AbstractVector{Union{Missing, Float64}})

Compute the standard deviation of the observed values in `X`.

# Examples
```jldoctest
julia> std_skipmissing([1.0; missing; 3.0])
1.4142135623730951
```
"""
std_skipmissing(X::AbstractVector{Float64}) = std(X);
std_skipmissing(X::AbstractVector{Union{Missing, Float64}}) = std(skipmissing(X));

"""
    std_skipmissing(X::AbstractMatrix{Float64})
    std_skipmissing(X::AbstractMatrix{Union{Missing, Float64}})

Compute the standard deviation of the observed values in `X` column wise.

# Examples
```jldoctest
julia> std_skipmissing([1.0 2.0; missing 3.0; 3.0 5.0])
3-element Array{Float64,1}:
   0.7071067811865476
 NaN
   1.4142135623730951
```
"""
std_skipmissing(X::AbstractMatrix{Float64}) = std(X, dims=2);

function std_skipmissing(X::AbstractMatrix{Union{Missing, Float64}})
    n = size(X,1);
    output = zeros(n,1);
    for i=1:n
        Xi = @view X[i,:];
        output[i] = std(skipmissing(Xi));
    end

    return output;
end

"""
    is_vector_in_matrix(vect::IntVector, matr::IntMatrix)
    is_vector_in_matrix(vect::FloatVector, matr::FloatMatrix)

Check whether the vector `vect` is included in the matrix `matr`.

# Examples
julia> is_vector_in_matrix([1; 2], [1 2; 2 3])
true
"""
is_vector_in_matrix(vect::IntVector, matr::IntMatrix) = sum(sum(vect .== matr, dims=1) .== length(vect)) > 0;
is_vector_in_matrix(vect::FloatVector, matr::FloatMatrix) = sum(sum(vect .== matr, dims=1) .== length(vect)) > 0;

"""
    isconverged(new::Float64, old::Float64, tol::Float64, ε::Float64, increasing::Bool)

Check whether `new` is close enough to `old`.

# Arguments
- `new`: new objective or loss
- `old`: old objective or loss
- `tol`: tolerance
- `ε`: small Float64
- `increasing`: true if `new` increases, at each iteration, with respect to `old`
"""
isconverged(new::Float64, old::Float64, tol::Float64, ε::Float64, increasing::Bool) = increasing ? (new-old)./(abs(old)+ε) <= tol : -(new-old)./(abs(old)+ε) <= tol;

"""
    soft_thresholding(z::Float64, ζ::Float64)

Soft thresholding operator.
"""
soft_thresholding(z::Float64, ζ::Float64) = sign(z)*max(abs(z)-ζ, 0.0);

"""
    square_vandermonde_matrix(λ::FloatVector)

Construct square vandermonde matrix on the basis of a vector of eigenvalues λ.
"""
square_vandermonde_matrix(λ::FloatVector) = λ'.^collect(length(λ)-1:-1:0);

"""
    solve_discrete_lyapunov(A::AbstractArray{Float64,2}, Q::SymMatrix)

Use a bilinear transformation to convert the discrete Lyapunov equation to a continuous Lyapunov equation, which is then solved using BLAS.

The notation used for representing the discrete Lyapunov equation is

``P - APA' = Q``,

where ``P`` and ``Q`` are symmetric. This equation is transformed into

``B'P + PB = -C``

# References
Kailath (1980, page 180)
"""
function solve_discrete_lyapunov(A::AbstractArray{Float64,2}, Q::SymMatrix)

    # Compute tranformed parameters
    inv_A_plus_I = inv(A+I);
    B_tr = inv_A_plus_I*(A-I); # alias for B'
    C = 2*inv_A_plus_I*Q*inv_A_plus_I'; # the scalar `2` is correct

    # Return solution
    return Symmetric(lyap(B_tr, C))::SymMatrix;
end

#=
--------------------------------------------------------------------------------------------------------------------------------
Parameter transformations
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    get_bounded_log(Θ_unbound::Float64, MIN::Float64)

Compute parameters with bounded support using a generalised log transformation.
"""
get_bounded_log(Θ_unbound::Float64, MIN::Float64) = exp(Θ_unbound) + MIN;

"""
    get_unbounded_log(Θ_bound::Float64, MIN::Float64)

Compute parameters with unbounded support using a generalised log transformation.
"""
get_unbounded_log(Θ_bound::Float64, MIN::Float64) = log(Θ_bound - MIN);

"""
    get_bounded_logit(Θ_unbound::Float64, MIN::Float64, MAX::Float64)

Compute parameters with bounded support using a generalised logit transformation.
"""
get_bounded_logit(Θ_unbound::Float64, MIN::Float64, MAX::Float64) = (MIN + (MAX * exp(Θ_unbound))) / (1 + exp(Θ_unbound));

"""
    get_unbounded_logit(Θ_bound::Float64, MIN::Float64, MAX::Float64)

Compute parameters with unbounded support using a generalised logit transformation.
"""
get_unbounded_logit(Θ_bound::Float64, MIN::Float64, MAX::Float64) = log((Θ_bound - MIN) / (MAX - Θ_bound));

#=
--------------------------------------------------------------------------------------------------------------------------------
Time series
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    demean(X::FloatVector)
    demean(X::FloatMatrix)
    demean(X::JVector{Float64})
    demean(X::JMatrix{Float64})

Demean `X`.

# Examples
```jldoctest
julia> demean([1.0; 1.5; 2.0; 2.5; 3.0])
5-element Array{Float64,1}:
 -1.0
 -0.5
  0.0
  0.5
  1.0

julia> demean([1.0 3.5 1.5 4.0 2.0; 4.5 2.5 5.0 3.0 5.5])
2×5 Array{Float64,2}:
 -1.4   1.1  -0.9   1.6  -0.4
  0.4  -1.6   0.9  -1.1   1.4
```
"""
demean(X::FloatVector) = X .- mean(X);
demean(X::FloatMatrix) = X .- mean(X,dims=2);
demean(X::JVector{Float64}) = X .- mean_skipmissing(X);
demean(X::JMatrix{Float64}) = X .- mean_skipmissing(X);

"""
    standardise(X::FloatVector)
    standardise(X::FloatMatrix)
    standardise(X::JVector{Float64})
    standardise(X::JMatrix{Float64})

Standardise `X`.

# Examples
```jldoctest
julia> standardise([1.0; 1.5; 2.0; 2.5; 3.0])
5-element Array{Float64,1}:
 -1.2649110640673518
 -0.6324555320336759
  0.0
  0.6324555320336759
  1.2649110640673518

julia> standardise([1.0 3.5 1.5 4.0 2.0; 4.5 2.5 5.0 3.0 5.5])
2×5 Array{Float64,2}:
 -1.08173    0.849934  -0.695401   1.23627   -0.309067
  0.309067  -1.23627    0.695401  -0.849934   1.08173
```
"""
standardise(X::FloatVector) = (X .- mean(X))./std(X);
standardise(X::FloatMatrix) = (X .- mean(X,dims=2))./std(X,dims=2);
standardise(X::JVector{Float64}) = (X .- mean_skipmissing(X))./std_skipmissing(X);
standardise(X::JMatrix{Float64}) = (X .- mean_skipmissing(X))./std_skipmissing(X);

"""
    interpolate_series(X::JMatrix{Float64}, n::Int64, T::Int64)

Interpolate each series in `X`, in turn, by replacing missing observations with the sample average of the non-missing values.

# Arguments
- `X`: observed measurements (`nxT`)
- `n` and `T` are the number of series and observations
"""
function interpolate_series(X::JMatrix{Float64}, n::Int64, T::Int64)
    data = copy(X);
    for i=1:n
        data[i, ismissing.(X[i, :])] .= mean_skipmissing(X[i, :]);
    end
    output = data |> FloatMatrix;
    return output;
end

"""
    interpolate_series(X::FloatMatrix, n::Int64, T::Int64)

Return `X`.

# Arguments
- `X`: observed measurements (`nxT`)
- `n` and `T` are the number of series and observations
"""
interpolate_series(X::FloatMatrix, n::Int64, T::Int64) = X;

"""
    interpolate_series(X::Union{FloatMatrix, JMatrix{Float64}})

Interpolate each series in `X`, in turn, by replacing missing observations with the sample average of the non-missing values.

# Arguments
- `X`: observed measurements (`nxT`)
"""
interpolate_series(X::Union{FloatMatrix, JMatrix{Float64}}) = interpolate_series(X, size(X)...);

"""
    forward_backwards_rw_interpolation(X::JMatrix{Float64}, n::Int64, T::Int64)

Interpolate each non-stationary series in `X`, in turn, using a random walk logic both forward and backwards in time.

# Arguments
- `X`: observed measurements (`nxT`)
- `n` and `T` are the number of series and observations
"""
function forward_backwards_rw_interpolation(X::JMatrix{Float64}, n::Int64, T::Int64)

    data = copy(X);

    # Loop over n and T
    for i=1:n

        # Forward interpolation
        drift_forward = mean_skipmissing(diff(data, dims=2))[:];
        for t=2:T
            if ismissing(data[i,t]) && ~ismissing(data[i,t-1])
                data[i,t] = data[i,t-1] + drift_forward[i];
            end
        end

        # Backwards interpolation
        drift_backwards = mean_skipmissing(diff(data, dims=2))[:];
        for t=T:-1:1
            if ismissing(data[i,t]) && ~ismissing(data[i,t+1])
                data[i,t] = data[i,t+1] - drift_backwards[i];
            end
        end
    end

    # Return interpolated data
    output = data |> FloatMatrix;
    return output;
end

"""
    forward_backwards_rw_interpolation(X::FloatMatrix, n::Int64, T::Int64)

Return `X`.

# Arguments
- `X`: observed measurements (`nxT`)
- `n` and `T` are the number of series and observations
"""
forward_backwards_rw_interpolation(X::FloatMatrix, n::Int64, T::Int64) = X;

"""
    forward_backwards_rw_interpolation(X::JMatrix{Float64})

Interpolate each non-stationary series in `X`, in turn, using a random walk logic both forward and backwards in time.

# Arguments
- `X`: observed measurements (`nxT`)
"""
forward_backwards_rw_interpolation(X::JMatrix{Float64}) = forward_backwards_rw_interpolation(X, size(X)...);

"""
    centred_moving_average(X::Union{FloatMatrix, JMatrix{Float64}}, n::Int64, T::Int64, window::Int64)

Compute the centred moving average of `X`.

# Arguments
- `X`: observed measurements (`nxT`)
- `n` and `T` are the number of series and observations
- `window` is the total number of observations (lagging, current and leading) included in the average
"""
function centred_moving_average(X::Union{FloatMatrix, JMatrix{Float64}}, n::Int64, T::Int64, window::Int64)

    # Checks on the moving average window
    check_bounds(T, 3);
    check_bounds(window, 3, T);

    # window_adj must be odd
    window_adj = copy(window);
    if mod(window-1,2) != 0
        window_adj = window-1;
    end
    one_sided_window = (window_adj-1)/2 |> Int64;

    # Compute centred moving average
    data = missing .* zeros(n,T) |> JMatrix{Float64};
    for t=one_sided_window+1:T-one_sided_window
        for i=1:n
            # Convenient view
            X_window = @view X[i, t-one_sided_window:t+one_sided_window];
            if sum(ismissing.(X_window)) < window_adj
                data[i,t] = mean_skipmissing(X[i, t-one_sided_window:t+one_sided_window]);
            end
        end
    end

    return data;
end

"""
    centred_moving_average(X::Union{FloatMatrix, JMatrix{Float64}}, window::Int64)

Compute the centred moving average of `X`.

# Arguments
- `X`: observed measurements (`nxT`)
- `window` is the total number of observations (lagging, current and leading) included in the average
"""
function centred_moving_average(X::Union{FloatMatrix, JMatrix{Float64}}, window::Int64)
    return centred_moving_average(X, size(X)..., window);
end

"""
    lag(X::FloatArray, p::Int64)

Construct the data required to run a standard vector autoregression.

# Arguments
- `X`: observed measurements (`nxT`), where `n` and `T` are the number of series and observations.
- `p`: number of lags in the vector autoregression

# Output
- `X_{t}`
- `X_{t-1}`
"""
function lag(X::FloatArray, p::Int64)

    # VAR(p) data
    X_t = X[:, 1+p:end];
    X_lagged = vcat([X[:, p-j+1:end-j] for j=1:p]...);

    # Return output
    return X_t, X_lagged;
end

"""
    companion_form(Ψ::AbstractArray{Float64,2}; extended::Bool=false)

Construct the companion form matrix from the generic coefficients Ψ.

If extended is true, it increases the typical dimension of the companion matrix by n rows.
"""
function companion_form(Ψ::AbstractArray{Float64,2}; extended::Bool=false)

    # Dimensions
    n = size(Ψ,1);
    p = Int64(size(Ψ,2)/n);

    # Number of rows of the block of identities in the companion form matrix
    standard_rows = n*(p-1);
    extra_rows = n*(extended == true);
    total_rows = standard_rows + extra_rows;

    # Return companion matrix
    companion = [Ψ zeros(n, extra_rows); Matrix(1.0I, total_rows, total_rows) zeros(total_rows, n)]::FloatMatrix;
    return companion;
end

#=
-------------------------------------------------------------------------------------------------------------------------------
Combinatorics and probability
-------------------------------------------------------------------------------------------------------------------------------
=#

"""
    no_combinations(n::Int64, k::Int64)

Compute the binomial coefficient of `n` observations and `k` groups, for big integers.

# Examples
```jldoctest
julia> no_combinations(1000000,100000)
7.333191945934207610471288280331309569215030711272858517142085449641265002716664e+141178
```
"""
no_combinations(n::Int64, k::Int64) = factorial(big(n))/(factorial(big(k))*factorial(big(n-k)));

"""
    rand_without_replacement(rng::StableRNGs.LehmerRNG, nT::Int64, d::Int64)

Draw `length(P)-d` elements from the positional vector `P` without replacement.

`P` is permanently changed in the process.

# Examples
```jldoctest
julia> rand_without_replacement(StableRNG(1), 20, 5)
5-element Array{Int64,1}:
  1
  9
 14
 19
 20
```
"""
function rand_without_replacement(rng::StableRNGs.LehmerRNG, nT::Int64, d::Int64)

    # Positional vector
    P = collect(1:nT);

    # Draw without replacement d times
    for i=1:d
        deleteat!(P, findall(P.==rand(rng, P)));
    end

    # Return output
    return setdiff(1:nT, P);
end

"""
    rand_without_replacement(rng::StableRNGs.LehmerRNG, n::Int64, T::Int64, d::Int64)

Draw `length(P)-d` elements from the positional vector `P` without replacement.

In the sampling process, no more than n-1 elements are removed for each point in time. `P` is permanently changed in the process.
"""
function rand_without_replacement(rng::StableRNGs.LehmerRNG, n::Int64, T::Int64, d::Int64)

    # Positional vector
    P = collect(1:n*T);

    # Full set of coordinates
    coord = [repeat(1:n, T) kron(1:T, convert(Array{Int64}, ones(n)))];

    # Counter
    coord_counter = convert(Array{Int64}, zeros(T));

    # Loop over d
    for i=1:d

        while true

            # New candidate draw
            draw = rand(rng, P);
            coord_draw = @view coord[draw, :];

            # Accept the draw if all observations are not missing for time t = coord[draw, :][2]
            if coord_counter[coord_draw[2]] < n-1
                coord_counter[coord_draw[2]] += 1;

                # Draw without replacement
                deleteat!(P, findall(P.==draw));
                break;
            end
        end
    end

    # Return output
    return setdiff(1:n*T, P);
end


text:
# Used "Conjugate Bayesian analysis of the Gaussian distribution" by Murphy as
# a reference.  Note that there were some typos in that document so the code
# here may not correspond exactly.

immutable NormalGamma <: Distribution
    mu::Float64
    nu::Float64     # scales precision of Normal
    shape::Float64
    rate::Float64

    function NormalGamma(mu::Real, nu::Real, sh::Real, r::Real)
    	nu > zero(nu) && sh > zero(sh) && r > zero(r) || error("Both shape and scale must be positive")
    	new(mu, nu, sh, r)
    end
end

mu(d::NormalGamma) = d.mu
nu(d::NormalGamma) = d.nu
shape(d::NormalGamma) = d.shape
scale(d::NormalGamma) = 1. / d.rate
rate(d::NormalGamma) = d.rate

insupport(::Type{NormalGamma}, x::Real, tau2::Real) =
    isfinite(x) && zero(tau2) <= tau2 < Inf

# Probably should guard agains dividing by and taking the log of 0.
function pdf(d::NormalGamma, x::Real, tau2::Real)
    Zinv = d.rate.^d.shape / gamma(d.shape) * sqrt(d.nu / (2.*pi))
    return Zinv * tau2.^(d.shape-0.5) * exp(-0.5*tau2*(d.nu*(x-d.mu).^2 + 2.*d.rate))
end
function logpdf(d::NormalGamma, x::Real, tau2::Real)
    lZinv = d.shape*log(d.rate) - lgamma(d.shape) + 0.5*(log(d.nu) - log(2.*pi))
    return lZinv + (d.shape-0.5)*log(tau2) - 0.5*tau2*(d.nu*(x-d.mu).^2 + 2*d.rate)
end

function rand(d::NormalGamma)
    # Guard against invalid precisions
    tau2 = rand(Gamma(d.shape, d.rate))
    if tau2 <= zero(Float64)
        tau2 = eps(Float64)
    end
    mu = rand(Normal(d.mu, 1./(tau2*d.nu)))
    return mu, tau2
end


text:
"""
Discrete pairwise Markov random field
"""
mutable struct DiscretePairwiseMRF{T<:AbstractFloat}
    varlist::SVarList
    gr::UGraph
    cards::Vector{Int}
    epdims::Vector{IntPair}

    # vertex potentials: unary potentials associated with vertices
    vpots::Vector{Vector{T}}

    # edge potentials: binary potentials associated with edge types
    epots::Vector{Matrix{T}}

    function DiscretePairwiseMRF{T}(vs::SVarList, gr::UGraph) where T<:AbstractFloat
        nv = length(vs)
        nvertices(gr) == nv || throw(
            ArgumentError("nvertices(gr) does not match #variables."))
        nt = nedgetypes(gr)

        cards = Vector{Int}(nv)
        for (i, v) in enumerate(vs)
            cards[i] = cardinality(v)
        end

        epdims = fill((0, 0), nt)
        for i = 1:nt
            es = edges(gr, i)
            if isempty(es)
                continue
            end
            (u1, v1) = first(es)
            cu, cv = cards[u1], cards[v1]
            for (u, v) in es
                (cards[u] == cu && cards[v] == cv) || throw(
                    DimensionMismatch("Inconsistent edge cardinalities."))
            end
            epdims[i] = (cu, cv)
        end
        new(vs, gr, cards, epdims, Vector{T}[], Matrix{T}[])
    end
end

nvars(mrf::DiscretePairwiseMRF) = length(mrf.varlist)
nvertices(mrf::DiscretePairwiseMRF) = nvertices(mrf.gr)
nedges(mrf::DiscretePairwiseMRF) = nedges(mrf.gr)
nedgetypes(mrf::DiscretePairwiseMRF) = nedgetypes(mrf.gr)
edges(mrf::DiscretePairwiseMRF, t::Int) = edges(mrf.gr, t)
vars(mrf::DiscretePairwiseMRF) = mrf.varlist
var_cardinalities(mrf::DiscretePairwiseMRF, v::Int) = mrf.cards[v]
edge_pdims(mrf::DiscretePairwiseMRF, t::Int) = mrf.epdims[t]

has_vertex_potentials(mrf::DiscretePairwiseMRF) = !isempty(mrf.vpots)
has_edge_potentials(mrf::DiscretePairwiseMRF) = !isempty(mrf.epots)
vertex_potentials(mrf::DiscretePairwiseMRF) = mrf.vpots
edge_potentials(mrf::DiscretePairwiseMRF) = mrf.epots

function set_vertex_potentials!{T<:AbstractFloat}(mrf::DiscretePairwiseMRF{T},
                                                  vpots::Vector{Vector{T}})
    nv = nvertices(mrf)
    length(vpots) == nv || throw(DimensionMismatch(
        "length(vpots) does not match the number of vertices."))

    vs = vars(mrf)
    for i = 1:nv
        if !isempty(vpots[i])
            length(vpots[i]) == mrf.cards[i] || throw(
                DimensionMismatch("Mismatched unary potential size."))
        end
    end
    mrf.vpots = vpots
end

function set_edge_potentials!{T<:AbstractFloat}(mrf::DiscretePairwiseMRF{T},
                                                epots::Vector{Matrix{T}})
    nt = nedgetypes(mrf)
    length(epots) == nt || throw(DimensionMismatch(
        "length(epots) does not match the number of edge types."))

    for i = 1:nt
        if !isempty(epots[i])
            size(epots[i]) == mrf.epdims[i] || throw(
                DimensionMismatch("Mismatched binary potential size."))
        end
    end
    mrf.epots = epots
end

function pwmrf{T<:AbstractFloat}(vs::SVarList, gr::UGraph,
                                 vpots::Vector{Vector{T}},
                                 epots::Vector{Matrix{T}})
    mrf = DiscretePairwiseMRF{T}(vs, gr)
    set_vertex_potentials!(mrf, vpots)
    set_edge_potentials!(mrf, epots)
    mrf
end

function tpotential{T<:AbstractFloat}(mrf::DiscretePairwiseMRF{T},
                                      x::Vector{Int})
    nv = nvertices(mrf)
    nt = nedgetypes(mrf)
    length(x) == nv || throw(DimensionMismatch("Invalid sample dimension."))

    # accumulate vertex potentials
    r1 = zero(T)
    if has_vertex_potentials(mrf)
        vpots = vertex_potentials(mrf)
        for i = 1:nv
            vp = vpots[i]
            isempty(vp) && continue
            r1 += vp[x[i]]
        end
    end

    # accumulate edge potentials
    r2 = zero(T)
    if has_edge_potentials(mrf)
        epots = edge_potentials(mrf)
        for t = 1:nt
            ep = epots[t]
            isempty(ep) && continue
            r2t = zero(T)
            es = edges(mrf, t)
            for (u, v) in es
                r2t += ep[x[u], x[v]]
            end
            r2 += r2t
        end
    end

    # total potential
    r1 + r2
end

function tpotentials{T<:AbstractFloat}(mrf::DiscretePairwiseMRF{T})
    nv = nvertices(mrf)
    nt = nedgetypes(mrf)

    # prepare storage
    r = zeros(T, mrf.cards...)

    # accumulate vertex potentials
    if has_vertex_potentials(mrf)
        vpots = vertex_potentials(mrf)
        for i = 1:nv
            vp = vpots[i]
            if !isempty(vp)
                shp = ones(Int, nv)
                shp[i] = mrf.cards[i]
                r .+= reshape(vp, shp...)
            end
        end
    end

    # accumulate edge potentials
    if has_edge_potentials(mrf)
        epots = edge_potentials(mrf)
        for t = 1:nt
            ep = epots[t]
            if !isempty(ep)
                ep_t = transpose(ep)
                es = edges(mrf, t)
                for (u, v) in es
                    @assert u != v
                    shp = ones(Int, nv)
                    shp[u] = size(ep, 1)
                    shp[v] = size(ep, 2)
                    if u < v
                        r .+= reshape(ep, shp...)
                    else
                        r .+= reshape(ep_t, shp...)
                    end
                end
            end
        end
    end
    r
end


text:
using Test
include("../src/Simulate.jl")

#########Test we can correctly find overlap#########

#Make 3D detected array
NumSimulations = 1
numGenes = 1
numCells = 1

detectedArr = Array{Float64}(undef, NumSimulations, numGenes, numCells)

#Check we can correctly find mean for 1 x 1 x 1 array
detectedArr[1,1,1] = 1.0
@test FindOverlap(detectedArr, 1, 1, 1) == [1.0]

detectedArr[1,1,1] = 0
@test FindOverlap(detectedArr, 1, 1, 1) == [0]

#Check we can correctly find mean for 2 x 2 x 2 array
NumSimulations = 2
numGenes = 2
numCells = 2

detectedArr = Array{Float64}(undef, NumSimulations, numGenes, numCells)
detectedArr[1,1,1] = 1.0
detectedArr[1,1,2] = 0.0
detectedArr[1,2,1] = 1.0
detectedArr[1,2,2] = 0.0
detectedArr[2,1,1] = 1.0
detectedArr[2,1,2] = 0.0
detectedArr[2,2,1] = 1.0
detectedArr[2,2,2] = 0.0

@test FindOverlap(detectedArr, 1, NumSimulations, numGenes) == [0.5, 0.5, 0.5, 0.5]

detectedArr = Array{Bool}(undef, 1, 2, 2, 2)
rankArr = Array{Int64}(undef, 1 ,2, 2, 2)

detectedArr[1,1,1,1] = true
detectedArr[1,1,1,2] = false
detectedArr[1,1,2,1] = true
detectedArr[1,1,2,2] = false
detectedArr[1,2,1,1] = true
detectedArr[1,2,1,2] = false
detectedArr[1,2,2,1] = true
detectedArr[1,2,2,2] = false

rankArr[1,1,1,1] = 1
rankArr[1,1,1,2] = 2
rankArr[1,1,2,1] = 1
rankArr[1,1,2,2] = 2
rankArr[1,2,1,1] = 2
rankArr[1,2,1,2] = 1
rankArr[1,2,2,1] = 2
rankArr[1,2,2,2] = 1

@test countRanks(detectedArr, rankArr, 2,2,2) == [-0.7745966692414833, -0.25819888974716115]

#Check we can correctly find P(dropout)
test_dict = Dict("I1" => [1,1,1], "I2" => [5,0,0], "I3" => [1,1,0])
@test FindPDropout(test_dict, 1.0) == [0.5,(1-15/24), (1-6/15)]

#Check we can correctly rank isoforms
@test FindRank(test_dict) == [2,3,1]

#Check we can correctly find p(Obs)
@test FindPObs(test_dict) == [1, 1/3, 2/3]

#Check we can correctly find the harmonic mean
@test FindHarmonicMean(2) == exp(-2.25) + (0.5 * exp(-4))

#Check we can correctly find the median frequencies
@test FindMedianFrequencies(2, 1) == [exp(-(1.5 * 1.5)), (exp(-4))/2]

#Check we can correctly scale frequencies
@test ScaleFrequencies([0.15, 0.1]) == [0.6, 0.4]
@test ScaleFrequencies([1.2, 0.8]) == [0.6, 0.4]

#Check we can correctly find cumulative frequencies
@test FindCumulativeFrequencies([0.5, 0.1, 0.2, 0.2]) ≈ [0.5, 0.7,0.9,1.0]

#Check can return rank of picked isoform correctly
@test PickRank([0.5,0.75,1.0], 0.8) == 3
@test PickRank([0.5,0.75,1.0], 0.1) == 1

@test RescaleFreqs([0.1,0.4,1.0], 3) == [0.8, 1.0]

@test RescaleRankings([4,1,3,2], 1, 2) == [3,2,1]
@test RescaleRankings([4,1,3,2], 4, 1) == [1,3,2]

Random.seed!(1234)
out = RandomIsoformChoiceProbabilities(2,2)

@test out[1,1] ≈ 0.5106336456707476 atol=0.001
@test out[1,2] ≈ 0.4893663543292524 atol=0.001
@test out[2,1] ≈ 0.624996535934181 atol=0.001
@test out[2,2] ≈ 0.3750034640658189 atol=0.001

#Tests for QuantErrorsForExpressed
pFP = 0.05
pFN = 0.1
@test QuantErrorsForExpressed(true, 0.05) == false
@test QuantErrorsForExpressed(true, 0.1) == true
@test QuantErrorsForExpressed(false, 0.01) == true
@test QuantErrorsForExpressed(false, 0.1) == false

Random.seed!(1234)
@test RandomQuantErrorsForExpressed(true, 0.5) == true
@test RandomQuantErrorsForExpressed(true, 0.2) == false
@test RandomQuantErrorsForExpressed(false, 0.2) == true
@test RandomQuantErrorsForExpressed(false, 0.3) == false

Random.seed!(1234)
@test RandomQuantErrorsForUnexpressed(0.3) == false
@test RandomQuantErrorsForUnexpressed(0.3) == true

Random.seed!(1234)
unexprRandArr = rand(1,1,1,2)
NumSimulations = 1
genesList = ["G1"]
numCells = 1
unexpr_dict = Dict("G1" => 2)
filteringThreshold = 4
NumIsoformsToSimulate =4
out = rand(1,1,1)
out[1,1,1] = 0

@test UnexpressedErrorsIsoformDependence(unexprRandArr, NumSimulations,
genesList, numCells, unexpr_dict, NumIsoformsToSimulate, filteringThreshold) == out

unexprRandArr[1,1,1,1] = 0.05
out[1,1,1] = 1

@test UnexpressedErrorsIsoformDependence(unexprRandArr, NumSimulations,
genesList, numCells, unexpr_dict, NumIsoformsToSimulate, filteringThreshold) == out

readCapturedIsoforms = Array{Bool}(undef, 1,1,1,4)
readCapturedIsoforms[1,1,1,1] = true
readCapturedIsoforms[1,1,1,2] = true
readCapturedIsoforms[1,1,1,3] = true
readCapturedIsoforms[1,1,1,4] = false

Random.seed!(1234)
randIsoChoiceArr = rand(1,1,1,4)

@test ExpressedErrorsIsoformDependence(readCapturedIsoforms, randIsoChoiceArr,
NumSimulations, genesList, numCells, unexpr_dict, NumIsoformsToSimulate,
filteringThreshold) == readCapturedIsoforms

randIsoChoiceArr[1,1,1,1] = 0.03
output = Array{Bool}(undef, 1,1,1,4)
output[1,1,1,1] = false
output[1,1,1,2] = true
output[1,1,1,3] = true
output[1,1,1,4] = false

@test ExpressedErrorsIsoformDependence(readCapturedIsoforms, randIsoChoiceArr,
NumSimulations, genesList, numCells, unexpr_dict, NumIsoformsToSimulate,
filteringThreshold) == output

randIsoChoiceArr[1,1,1,4] = 0.05
output[1,1,1,4] = true

@test ExpressedErrorsIsoformDependence(readCapturedIsoforms, randIsoChoiceArr,
NumSimulations, genesList, numCells, unexpr_dict, NumIsoformsToSimulate,
filteringThreshold) == output

###############################################################################
#Test that we can correctly pick isoforms

Random.seed!(1234)
randIsoChoiceArr = rand(1,1,1,1) #0.59
IsoChoiceProbabilities = [0.1, 0.2]
pObsDict = Dict("G1" => IsoChoiceProbabilities)
NumSimulations = 1
genesList = ["G1"]
rankingDict = Dict("G1" => [2,1])
pDropoutDict = Dict("G1" => [0.1, 0.05])
numCells = 1
NumIsoformsToSimulate = 1

Random.seed!(1234) #Seed so sampling is reproducible

testArr = pickIsoforms(IsoChoiceProbabilities, NumSimulations,
    genesList, rankingDict, pDropoutDict, numCells,
    NumIsoformsToSimulate)

@test testArr[1,1,1,1] == 0.1

IsoChoiceProbabilities = rand(1, 2)
IsoChoiceProbabilities[1,1] = 0.1
IsoChoiceProbabilities[1,2] = 0.2

Random.seed!(1234) #Seed so sampling is reproducible

testArr = pickRandIsoforms(IsoChoiceProbabilities, NumSimulations,
    genesList, pDropoutDict, numCells,
    NumIsoformsToSimulate)

@test testArr[1,1,1,1] == 0.05

Random.seed!(1234)
testArr[1,1,1,1] = 0.0

testArr = PickUnifObsIsoforms(NumIsoformsToSimulate, genesList, pObsDict,
pDropoutDict, numCells, NumSimulations)

@test testArr[1,1,1,1] == 0.05

Random.seed!(1234)
testArr[1,1,1,1] = 0.0

testArr = PickCellVarIsoforms(NumIsoformsToSimulate, genesList, pObsDict,
pDropoutDict, numCells, NumSimulations)

@test testArr[1,1,1,1] == 0.05

################################################################################
#Test can correctly simulate software FPs
NumIsoformsToSimulate = 4
filteringThreshold =4
Random.seed!(1234)
MaxUnexpressedDetectedIsoforms = map(x -> x > 0.6, rand(1,1,1,4)) #[0.59, 0.76, 0.57, 0.46]
unexpr_dict = Dict("G1" => 1)

#Test can do negatives
testArr = pickUnexpressedErrors(MaxUnexpressedDetectedIsoforms, NumSimulations,
genesList, numCells, unexpr_dict, NumIsoformsToSimulate, filteringThreshold)
@test testArr[1,1,1,1] == 0

#Test can do negatives
testArr = pickUnexpressedErrors(MaxUnexpressedDetectedIsoforms, NumSimulations,
genesList, numCells, unexpr_dict, 3, filteringThreshold)
@test testArr[1,1,1,1] == 1

#Test can do postives
unexpr_dict = Dict("G1" => 2)
testArr = pickUnexpressedErrors(MaxUnexpressedDetectedIsoforms, NumSimulations,
genesList, numCells, unexpr_dict, NumIsoformsToSimulate, filteringThreshold)
@test testArr[1,1,1,1] == 1

#Test can do random isoform choice correctly

###############################################################################
#Test globalArraySimulation function

Random.seed!(1234)
#randIsoChoiceArr = 0.59
#IsoChoiceArr = 0.1
#randIsoChoiceArr =0.85
#unexprRandArr = all false
unexpr_dict = Dict("G1" => 1)
filteringThreshold = 2
MaxNumIsos = 6

testArr = globalArraySimulation(NumSimulations, genesList, rankingDict,
    pDropoutDict, numCells, NumIsoformsToSimulate, filteringThreshold,
    unexpr_dict, MaxNumIsos)

@test testArr == [1]

################################################################################
#Test big simulation function
cumulativeFreqs = [0.2, 1.0]
rankings = [1,2]
pDropouts = [0.9, 0.1]
numCells = 1
IsoformsPerGenePerCell = 1
test = true
pFP = 0.17
pFN = 0.007

#Test the simplest case
@test Simulate(cumulativeFreqs,rankings, pDropouts, numCells,
IsoformsPerGenePerCell, 2, test) == 1

#Test when two isoforms are expressed
@test Simulate(cumulativeFreqs,rankings, pDropouts, numCells,
2, 2, test) == 1

#Test when we swap cumulativeFreqs
@test Simulate([0.8, 1.0], [2,1], pDropouts, numCells,
IsoformsPerGenePerCell, 2, test) == 1

#Test when we swap cumulativeFreqs
@test Simulate([0.8, 1.0], [2,1], pDropouts, numCells,
2, 2, test) == 1

#Test the simplest case when there is more than 1 cell
@test Simulate(cumulativeFreqs,rankings, [0.9, 0.8], 2,
IsoformsPerGenePerCell, 2, test) == 0.5

#Test when we have 3 isoforms
cumulativeFreqs = [0.4, 0.8, 1.0]
rankings = [1,3,2]
pDropouts = [0.9, 0.1, 0.5]

#Test for 1 isoform out of 3
@test Simulate(cumulativeFreqs,rankings, pDropouts, numCells,
IsoformsPerGenePerCell, 3, test) == 1

#Test for 2 isoforms out of 3
@test Simulate(cumulativeFreqs,rankings, pDropouts, numCells,
2, 3, test) == 1

#Test for 3 isoforms out of 3
@test Simulate(cumulativeFreqs,rankings, pDropouts, numCells,
3, 3, test) == 2

#Test that FP and FN logic works
cumulativeFreqs = [0.2, 1.0]
rankings = [1,2]
pDropouts = [0.9, 0.1]
numCells = 1
IsoformsPerGenePerCell = 1
test = true
pFP = 0.6
pFN = 0.007

@test Simulate(cumulativeFreqs,rankings, pDropouts, numCells,
IsoformsPerGenePerCell, 2, test) == 2

pFP = 0.0
pFN = 0.8

@test Simulate(cumulativeFreqs,rankings, pDropouts, numCells,
IsoformsPerGenePerCell, 2, test) == 0

#Test for when ObsExIndex != index
cumulativeFreqs = [0.25,0.5,0.75, 1.0]
rankings = [1,2,3,4]
pDropouts = [0.0,0.0,0.0,0.0]
IsoformsPerGenePerCell = 3
pFP = 0.0
pFN = 0.0

@test Simulate(cumulativeFreqs,rankings, pDropouts, numCells,
IsoformsPerGenePerCell, 4, test) == 3

#Test that number of isoforms logic works
pFP = 1.0
pFN = 0.0

@test Simulate(cumulativeFreqs,rankings, pDropouts, numCells,
IsoformsPerGenePerCell, 10, test) == 10

pFP = 0.0
pFN = 1.0

@test Simulate(cumulativeFreqs,rankings, pDropouts, numCells,
IsoformsPerGenePerCell, 10, test) == 0

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
using Trilinos
using Test

L = 10
teuchosarray = Teuchos.Array(L, "test")
@test size(teuchosarray) == (L,)
@test length(teuchosarray) == L
@test axes(teuchosarray) == (1:L,)

for (i,str) in enumerate(teuchosarray)
  teuchosarray[i] *= string(i)
end

println(teuchosarray)

for (i,str) in enumerate(teuchosarray)
  @test teuchosarray[i] == "test$i"
end

@test convert(Int32,Teuchos.null) == 0


text:
"""
	`mutable struct GLShader`

Container for vertex and fragment shaders, their ids, and program id.
"""
mutable struct GLShader

	vertex_source
	frag_source

	program_id::Int32
	vertex_shader_id::Int32
	frag_shader_id::Int32

	# constructor
	function GLShader(vertex, fragment)
		ret=new(vertex,fragment,-1,-1,-1)
		finalizer(releaseGpuResources, ret)
		return ret
	end

end

"""
	releaseGpuResources(shader::GLShader)

Management of shader program memory on the GPU.
"""
function releaseGpuResources(shader::GLShader)

	global __release_gpu_resources__

	if shader.vertex_shader_id>=0
	 	id=shader.vertex_shader_id
		shader.vertex_shader_id=-1
		glDeleteLater(function()  glDeleteShader(id) end)
	end

	if shader.frag_shader_id>=0
		id=shader.frag_shader_id
		shader.frag_shader_id=-1
		glDeleteLater(function()  glDeleteShader(id) end)

	end

	if shader.program_id>=0
		id=shader.program_id
		shader.program_id=-1
		glDeleteLater(function()  glDeleteProgram(id) end)
	end

end

"""
	createShader(type,source)

Shaders creation and compilation by `type`. Return `shader_id`
"""
function createShader(type,source)
	shader_id = glCreateShader(type)::GLuint
	glCheckError()
	glShaderSource(shader_id, 1, convert(Ptr{UInt8}, pointer([convert(Ptr{GLchar}, pointer(source))])) , C_NULL)
	glCompileShader(shader_id)
	status = GLint[0]
	glGetShaderiv(shader_id, GL_COMPILE_STATUS, status)
	if status[1] == GL_FALSE
		maxlength = 8192
		buffer = zeros(GLchar, maxlength)
		sizei = GLsizei[0]
		glGetShaderInfoLog(shader_id, maxlength, sizei, buffer)
		len = sizei[]
		error_msg=unsafe_string(pointer(buffer), len)
		error("shader compilation failed\n",error_msg,"\nsource\n",source)
	end
	return shader_id
end

"""
	enableProgram(shader)

Compiling and linking program shaders. Return status.

A program object will contain an executable that will run on the vertex processor if it contains one or more shader objects of type `GL_VERTEX_SHADER` that have been successfully compiled and linked.
Similarly, a program object will contain an executable that will run on the fragment processor if it contains one or more shader objects of type `GL_FRAGMENT_SHADER` that have been successfully compiled and linked.

# Example
```
shader = GL.GLPhongShader(true,true)
GL.enableProgram(shader)
```

"""
function enableProgram(shader)

	if (shader.program_id<0)

		shader.program_id = glCreateProgram()
		glCheckError()

		shader.vertex_shader_id=createShader(GL_VERTEX_SHADER,shader.vertex_source)
		glAttachShader(shader.program_id, shader.vertex_shader_id)
		glCheckError()

		shader.frag_shader_id=createShader(GL_FRAGMENT_SHADER,shader.frag_source)
		glAttachShader(shader.program_id, shader.frag_shader_id)
		glCheckError()

		glLinkProgram(shader.program_id)
		glCheckError()
		status = GLint[0]
		glGetProgramiv(shader.program_id, GL_LINK_STATUS, status)
		if status[1] == GL_FALSE
		 	maxlength = 8192
			buffer = zeros(GLchar, maxlength)
			sizei = GLsizei[0]
			glGetProgramInfoLog(shader.program_id, maxlength, sizei, buffer)
			len = sizei[]
			error_msg = unsafe_string(pointer(buffer), len)
			error("Error linking program\n",error_msg)
		end
		glCheckError()
	end

	glUseProgram(shader.program_id)
end

"""
	disableProgram(shader)

Make sure that the last used shader is not active anymore.

`glUseProgram` is available only if the GL version is 2.0 or greater. It specifies the handle of the program object whose executables are to be used as part of current rendering state.
"""
function disableProgram(shader)
	glUseProgram(0)
end


text:
@testset "heapify!" begin
    @testset for order = (Forward,Reverse)
        @testset for n = 0
            v = []
            heapify!(v, order)
            @test v == []
        end
        @testset for n = 1
            v = [0]
            heapify!(v, order)
            @test v == [0]
        end
        @testset for n = 2:10
            Random.seed!(42)
            for _ = 1:10
                v = rand(1:100,n)
                heapify!(v,order)
                @test isheap(v,order)
            end
        end
    end
end

@testset "bubble!" begin
    @testset for order = (Forward,Reverse)
        @testset for n = 0
            v = []
            i = bubble!((v,0),1,0,order) == []
            @test v == []
        end
        @testset for n = 1
            v = [0]
            i = bubble!((v,1),1,1,order)
            @test i == 1
            @test v == [1]
        end
        @testset for n = 2:10
            Random.seed!(42)
            for _ = 1:10
                v = rand(1:100,n)
                heapify!(v,order)
                vi = rand(1:100)
                i = bubble!((v,n),vi,rand(1:n),order)
                @test v[i] == vi
                @test isheap(v,order)
            end
        end
    end
end


text:
"""
Add variables to the OptimizationContainer for any component.
"""
function add_variables!(
    container::OptimizationContainer,
    ::Type{T},
    devices::Union{Vector{U}, IS.FlattenIteratorWrapper{U}},
    formulation::Union{AbstractDeviceFormulation, AbstractServiceFormulation},
) where {T <: AuxVariableType, U <: PSY.Component}
    add_variable!(container, T(), devices, formulation)
    return
end

@doc raw"""
Default implementation of adding auxiliary variable to the model.
"""
function add_variable!(
    container::OptimizationContainer,
    var_type::AuxVariableType,
    devices::U,
    formulation,
) where {U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}} where {D <: PSY.Component}
    @assert !isempty(devices)
    time_steps = get_time_steps(container)
    add_aux_variable_container!(
        container,
        var_type,
        D,
        [PSY.get_name(d) for d in devices],
        time_steps,
    )
    return
end


text:
# This is part of Quiz.jl. License is MIT:
# https://raw.githubusercontent.com/cnaak/Quiz.jl/main/LICENSE

module Quiz

include("interface.jl")

using Reexport
@reexport using Unitful

include("types.jl")

include("engines.jl")

end # module