text:
module CountingProcesses

using Distributions
using StaticArrays
using Random

export DiscreteCountingProcess, BernoulliScheme, StaggeredBernoulliScheme, NonhomogenousBernoulliScheme

include("poisson_process.jl")
include("bernoulli_process.jl")
include("bernoulli_scheme.jl")
include("staggered_bernoulli_scheme.jl")
include("nonhomgen_bernoulli_scheme.jl")

end


text:
for name in (:real, :imag, :conj)
    @eval begin
        Base.$name(eb::Eigenbrot) = Eigenbrot(Array{ComplexF64}($name(eb.vals)), eb.fft)
    end
end

for name in (:angle, :abs)
    @eval begin
        Base.$name(eb::Eigenbrot) = Eigenbrot(Array{ComplexF64}($name.(eb.vals)), eb.fft)
    end
end

phase(eb::Eigenbrot) = angle(eb)

import Base: (+), (-), (*), (/), (^)

(+)(eb1::Eigenbrot, eb2::Eigenbrot) = Eigenbrot(eb1.vals .+ eb2.vals, eb1.fft)
(+)(eb1::Eigenbrot, n::Number) = Eigenbrot(eb1.vals .+ n, eb1.fft)
(+)(n::Number, eb1::Eigenbrot) = eb1 + n
(-)(eb1::Eigenbrot, eb2::Eigenbrot) = Eigenbrot(eb1.vals .- eb2.vals, eb1.fft)
(-)(eb1::Eigenbrot, n::Number) = Eigenbrot(eb1.vals .- n, eb1.fft)
(-)(n::Number, eb1::Eigenbrot) = Eigenbrot(n .- eb1.vals, eb1.fft)

(*)(eb::Eigenbrot, x::Number) = x * eb
(*)(x::Number, eb::Eigenbrot) = Eigenbrot(x * eb.vals, eb.fft)
(*)(eb1::Eigenbrot, eb2::Eigenbrot) = Eigenbrot(eb1.vals .* eb2.vals, eb1.fft)
(/)(eb::Eigenbrot, x::Number) = Eigenbrot(eb.vals / x, eb.fft)
(/)(eb1::Eigenbrot, eb2::Eigenbrot) = Eigenbrot(eb1.vals ./ eb2.vals, eb1.fft)

#=
Utility method used by pad() methods, not exported.
=#
function _pad(v::Matrix{ComplexF64}, value::Complex,
               left::Integer, top::Integer, right::Integer, bottom::Integer)
    rows, cols = size(v)
    cols_out = cols + left + right
    rows_out = rows + top + bottom
    (cols_out ≤ 0 || rows_out ≤ 0) &&
        error("Invalid resize: $(cols)x$rows -> $(cols_out)x$rows_out")

    col_dest1 = left > 0 ? 1 + left : 1
    col_src1 = left > 0 ? 1 : 1 - left
    col_src2 = right > 0 ? cols : cols + right
    col_dest2 = col_dest1 + col_src2 - col_src1

    row_dest1 = bottom > 0 ? 1 + bottom : 1
    row_src1 = bottom > 0 ? 1 : 1 - bottom
    row_src2 = top > 0 ? rows : rows + top
    row_dest2 = row_dest1 + row_src2 - row_src1

    v2 = fill(complex(value), rows_out, cols_out)
    v2[row_dest1:row_dest2, col_dest1:col_dest2] = v[row_src1:row_src2, col_src1:col_src2]
    return v2
end

"""
    pad(eb::Eigenbrot, margin::Integer)
Pad edges of `eb` by `margin` pixels. If `margin` is negative,
extend by that amount, filling the new pixels with zero.

    pad(eb::Eigenbrot, value, margin::Integer)
Pad edges of `eb` by `margin` pixels. If `margin` is negative,
extend by that amount, filling the new pixels with `value`.

    pad(eb::Eigenbrot, value = 0; left::Integer = 0, top::Integer = 0, right::Integer = 0, bottom::Integer = 0)
Pad the named edges of `eb` by the number of pixels given. Negative
edge values will extend by the relevant amount, filling the new pixels
with `value`.
"""
pad(eb::Eigenbrot, value::Number = 0.0im; left::Integer = 0, top::Integer = 0,
        right::Integer = 0, bottom::Integer = 0) =
    Eigenbrot(_pad(eb.vals, ComplexF64(value), left, top, right, bottom), eb.fft)

pad(eb::Eigenbrot, value::Real, margin::Integer) =
    Eigenbrot(_pad(eb.vals, ComplexF64(value), margin, margin, margin, margin), eb.fft)

pad(eb::Eigenbrot, value::Complex, margin::Integer) =
    Eigenbrot(_pad(eb.vals, value, margin, margin, margin, margin), eb.fft)

pad(eb::Eigenbrot, margin::Integer) =
   pad(eb, 0.0im, margin)

function flipver(eb::Eigenbrot)
    flip = copy(eb)
    h = size(flip, 1)
    for r in 1:h
        reverse!(@view flip.vals[r, :])
    end
    return flip
end

function fliphor(eb::Eigenbrot)
    flip = copy(eb)
    w = size(flip, 2)
    for c in 1:w
        reverse!(@view flip.vals[:, c])
    end
    return flip
end

function swapxy(eb::Eigenbrot, xflip = false)
    h, w = size(eb)
    swap = Eigenbrot(w, h)
    if xflip
        for c in 1:h
            swap.vals[:, c] = @view eb.vals[c, :]
        end
    else
        for c in 1:h
            swap.vals[w:-1:1, h-c+1] = @view eb.vals[c, :]
        end
    end
    return swap
end

ext: jl


text:
# Install/Manage Genome Annotations
# currently downloads annotation files to a local store .GenomeAnnotations.jl/annotations

module GenomeAnnotations
 using Base.Intrinsics
 importall Base

export setup_genome_annotations
export install_ucsc_annotation
export install_gencode_annotation
export current_annotation_files

# Config.jl
export SET_LOCAL_STORE_PATH
export LOCAL_STORE_PATH

include(Pkg.dir("GenomeAnnotations","src","Config.jl"))
include(Pkg.dir("GenomeAnnotations","src","Store.jl"))
include(Pkg.dir("GenomeAnnotations","src","UCSCGB.jl"))
include(Pkg.dir("GenomeAnnotations","src","Gencode.jl"))
end


text:
# This file was generated, do not modify it. # hide
using Pkg  # hideall
Pkg.activate("_literate/D0-loading/Project.toml")
Pkg.update()

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
export something

something = "green"


text:
module MyActor2
export Actor, Message, Context, start, ActorRef, callOn, @startat, stop, _Rep, mergeFuture
import Distributed

abstract type Actor end
abstract type Message end
abstract type Context end
const MessagePair = Tuple{Message, Distributed.RRID}

struct _Rep <: Message end
struct _Stop <: Message end
function handle(a::Actor, m::Message) println("should not be called") end
function handle(a::Actor, m::_Rep) println(a) end
##
function taskLoop(cn::AbstractChannel, target::Actor)
    while true
        msg, rrid = take!(cn)
        fcn = Distributed.lookup_ref(rrid)
        if typeof(msg) <: _Stop
            put!(fcn, nothing)
            break
        end
        try
            res = handle(target, msg)
            put!(fcn, res)
        catch e
            println(e)
            put!(fcn, e)
        end
    end
end
##
function createLoop(actorCons)
    cn = Channel()
    @async taskLoop(cn, actorCons())
    return cn
end

function start(actorCons, node=1)
    cn = Distributed.RemoteChannel(()->createLoop(actorCons), node)
    ActorRef(cn, node)
end

macro startat(id, expr)
    thunk = esc(:(()->($expr)))
    quote
        start($thunk, $(esc(id)))
    end
end

##
struct ActorRef
    cn::Distributed.RemoteChannel
    w::Integer
end
#ActorRef(cn::Distributed.RemoteChannel, id::Integer) =
#    ActorRef(cn, Distributed.worker_from_id(id))

function callOn(ref::ActorRef, mes::Message)
    rr = Distributed.Future(ref.w)
    RRID = Distributed.remoteref_id(rr)
    put!(ref.cn, (mes, RRID))
    return rr
end

function stop(ref::ActorRef)
    callOn(ref, _Stop())
    close(ref.cn)
end

function mergeFuture(futures)
    cn = Channel()
    counter = length(futures)
    for (i, f) in enumerate(futures)
        @async begin
            put!(cn, fetch(f));
            counter-=1;
            if counter == 0
                close(cn)
            end
        end
    end
    cn
end

end ## module

max_stars_repo_name: HidemotoNakada/julia-actor

max_issues_repo_name: HidemotoNakada/julia-actor

max_forks_repo_name: HidemotoNakada/julia-actor


text:
using LinearAlgebra: qr, I, norm
using LowRankApprox: pqrfact
using IterativeSolvers
using PyCall

include("bayesianlinear.jl")

# TODO:
#   - read_dict, write_dict
#   - LSQR
#   - various scikit-learn solvers

@doc raw"""
`struct QR` : linear least squares solver, using standard QR factorisation;
this solver computes
```math
 θ = \arg\min \| A \theta - y \|^2 + \lambda \| P \theta \|^2
```
Constructor
```julia
ACEfit.QR(; λ = 0.0, P = nothing)
```
where
* `λ` : regularisation parameter
* `P` : right-preconditioner / tychonov operator
"""
struct QR
   λ::Number
   P
end

QR(; λ = 0.0, P = I) = QR(λ, P)

function solve_llsq(solver::QR, A, y)
   if solver.λ == 0
      AP = A
      yP = y
   else
      AP = [A; solver.λ * solver.P]
      yP = [y; zeros(eltype(y), size(A, 2))]
   end
   return qr(AP) \ yP
end

@doc raw"""
`struct RRQR` : linear least squares solver, using rank-revealing QR
factorisation, which can sometimes be more robust / faster than the
standard regularised QR factorisation. This solver first transforms the
parameters ``\theta_P = P \theta``, then solves
```math
 θ = \arg\min \| A P^{-1} \theta_P - y \|^2
```
where the truncation tolerance is given by the `rtol` parameter, and
finally reverses the transformation. This uses the `pqrfact` of `LowRankApprox.jl`;
For further details see the documentation of
[`LowRankApprox.jl`](https://github.com/JuliaMatrices/LowRankApprox.jl#qr-decomposition).

Crucially, note that this algorithm is *not deterministic*; the results can change
slightly between applications.

Constructor
```julia
ACEfit.RRQR(; rtol = 1e-15, P = I)
```
where
* `rtol` : truncation tolerance
* `P` : right-preconditioner / tychonov operator
"""
struct RRQR
   rtol::Number
   P
end

RRQR(; rtol = 1e-15, P = I) = RRQR(rtol, P)

function solve_llsq(solver::RRQR, A, y)
   AP = A / solver.P
   θP = pqrfact(AP, rtol = solver.rtol) \ y
   return solver.P \ θP
end

@doc raw"""
LSQR
"""
struct LSQR
   damp::Number
   atol::Number
   P
end

LSQR(damp, atol) = LSQR(damp, atol, I)
LSQR(; damp=5e-3, atol=1e-6) = LSQR(damp, atol)

function solve_llsq(solver::LSQR, A, y)
   println("damp  ", solver.damp)
   println("atol  ", solver.atol)
   # preconditioning
   if !(solver.P == I)
      @info("LSQR: Using preconditioning")
      D_inv = pinv(P)
      mul!(A,A,D_inv)
   end
   c, ch = lsqr(A, y, damp=solver.damp, atol=solver.atol, conlim=1e12, log=true, verbose=false)
   println(ch)
   println("relative RMS error  ", norm(A*c - y) / norm(y))
   return c
end

function solve_llsq(solver::LSQR, A::DArray, y::DArray)
   c = dzeros((size(A,2),), [1])
   lsqr!(c, A, y, damp=solver.damp, atol=solver.atol)
   return convert(Vector, c)
end

@doc raw"""
SKLEARN_BRR
"""
struct SKLEARN_BRR
end

function solve_llsq(solver::SKLEARN_BRR, A, y)
   BRR = pyimport("sklearn.linear_model")["BayesianRidge"]
   clf = BRR()
   clf.fit(A, y)
   c = clf.coef_
   return c
end

@doc raw"""
SKLEARN_ARD
"""
struct SKLEARN_ARD
end

function solve_llsq(solver::SKLEARN_ARD, A, y)
   ARD = pyimport("sklearn.linear_model")["ARDRegression"]
   clf = ARD(n_iter=100000, compute_score=true, verbose=true)
   clf.fit(A, y)
   println("scores length  ", length(clf.scores_))
   c = clf.coef_
   return c
end

@doc raw"""
Bayesian Linear
"""
struct BL
end

function solve_llsq(solver::BL, A, y)
   c, _, _, _ = BayesianRegression.bayesian_fit(y, A; verbose=false)
   return c
end

@doc raw"""
Bayesian ARD
"""
struct BARD
end

function solve_llsq(solver::BARD, A, y)
   c, _, _, _, _ = BayesianRegression.ard_fit(y, A; verbose=false)
   return c
end


text:
__precompile__()
module ConfParser

import Base.merge!
using Compat
import Compat.String

export ConfParse, parse_conf!, erase!,
       save!, retrieve, commit!

mutable struct ConfParse
    _fh::IO
    _filename::String
    _syntax::String
    _data::Dict
    _is_modified::Bool

    function ConfParse(filename::String, syntax::String = "")
        if (isempty(filename))
            error("no file name specified")
        end

        _filename = filename
        _fh = open_fh(filename, "r") # _fh was defined inside if-statement: so if syntax!="", then _fh is not defined
        if (isempty(syntax))
            _syntax = guess_syntax(_fh)
        else
            if ((syntax != "ini")  &&
                (syntax != "http") &&
                (syntax != "simple"))
                error("unknown configuration syntax: $(syntax)")
            end
            _syntax = syntax
        end
        _data        = Dict()
        _is_modified = false
        new(_fh, _filename, _syntax,
            _data, _is_modified)
    end # function ConfParse

end # type ConfParse

#----------
# open file handler for IO
#----------
function open_fh(filename::String, mode::String)
    try
        fh = open(filename, mode)
        return fh
    catch
        error("configuration file could not be opened")
    end
end

#------------
# attempts to guess the configuration file syntax using
# regular expressions
#----------
function guess_syntax(fh::IO)
    syntax = ""
    for line in eachline(fh)

        # is a commented line
        if (occursin(r"^\s*(?:#|$)", line))
            continue
        end

        # is not alphanumeric
        if (!occursin(r"\w", line))
            continue
        end

        # remove \n
        line = chomp(line)

        # contains a [block]; ini
        if (occursin(r"^\s*\[\s*[^\]]+\s*\]\s*$", line))
            syntax = "ini"
            break
        end

        # key/value pairs are seperated by a '='; ini
        if (occursin(r"^\s*[\w-]+\s*=\s*.*\s*$", line))
            syntax = "ini"
            break
        end

        # key/value pairs are seperated by a ':'; http
        if (occursin(r"^\s*[\w-]+\s*:\s*.*\s*$", line))
            syntax = "http"
            break
        end

        # key/value pairs are seperated by whitespace; simple
        if (occursin(r"^\s*[\w-]+\s+.*$", line))
            syntax = "simple"
            break
        end
    end

    if (syntax != "")
        return syntax
    end

    error("unable to identify the configuration file syntax")
end # function guess_syntax

#----------
# tasks appropriate parser function based on configuration
# syntax
#----------
function parse_conf!(s::ConfParse)
    if (s._syntax == "ini")
        parse_ini(s)
    elseif (s._syntax == "http")
        parse_http(s)
    elseif (s._syntax == "simple")
        parse_simple(s)
    else
        error("unknown configuration syntax: $(s._syntax)")
    end

    close(s._fh)                # no need to keep file opened
end # function parse_conf

#----------
# Sperates by commas, removes newlines and such
#----------
function parse_line(line::String)
    parsed   = (AbstractString)[]
    splitted = split(line, ",")
    for raw = splitted
        if (occursin(r"\S+", raw))
            clean = match(r"\S+", raw)
            push!(parsed, clean.match)
        end
    end

    parsed
end # function parse_line

#----------
# parses configuration files utilizing ini sytnax.
# Populate the ConfParser.data dictionary
#----------
function parse_ini(s::ConfParse)
    blockname = "default"
    seekstart(s._fh)
    for line in eachline(s._fh)
        # skip comments and newlines
        if (occursin(r"^\s*(\n|\#|;)", line))
            continue
        end

        if (!occursin(r"\w", line))
            continue
        end

        line = chomp(line)

        # parse blockname
        m = match(r"^\s*\[\s*([^\]]+)\s*\]$", line)
        if (m != nothing)
            blockname = lowercase(m.captures[1])
            continue
        end

        # parse key/value
        m = match(r"^\s*([^=]*\w)\s*=\s*(.*)\s*$", line)
        if (m != nothing)
            key::String, values::String = m.captures
            if (!haskey(s._data, blockname))
                s._data[blockname] = Dict(key => parse_line(values))
            else
                merge!(s._data[blockname], Dict(key => parse_line(values)))
            end
            continue
        end
        error("invalid syntax on line: $(line)")
    end
end # function parse_ini

#----------
# parses configuration files utilizing http sytnax.
# Populate the ConfParser.data dictionary
#----------
function parse_http(s::ConfParse)
    seekstart(s._fh)
    for line in eachline(s._fh)
        # skip comments and newlines
        if (occursin(r"^\s*(\n|\#|;)", line))
            continue
        end

        if (!occursin(r"\w", line))
            continue
        end

        line = chomp(line)

        m = match(r"^\s*([\w-]+)\s*:\s*(.*)$", line)
        if (m != nothing)
            key::String, values::String = m.captures
            s._data[key] = parse_line(values)
            continue
        end

        error("invalid syntax on line: $(line)")
    end
end # function parse_http

#----------
# parses configuration files utilizing simple syntax.
# Populates the ConfParser.data dictionary
#----------
function parse_simple(s::ConfParse)
    seekstart(s._fh)
    for line in eachline(s._fh)
        # skip comments and newlines
        if (occursin(r"^\s*(\n|\#|;)", line))
            continue
        end

        if (!occursin(r"\w", line))
            continue
        end

        line = chomp(line)

        m = match(r"^\s*([\w-]+)\s+(.*)\s*$", line)
        if (m != nothing)
            key::String, values::String = m.captures
            s._data[key] = parse_line(values)
            continue
        end

        error("invalid syntax on line: $(line)")
    end
end # function parse_simple

#----------
# craft content strings from data array for saved config
#----------
function craft_content(s::ConfParse)
   content = ""
   if (s._syntax == "ini")
        for (block, key_values) = s._data
            content *= "[$block]\n"
            for (key, values) = key_values
                if (typeof(values) == Array{AbstractString, 1})
                    content *= "$key=$(join(values, ","))\n"
                else
                    content *= "$key=$values\n"
                end
            end
            content *= "\n"
        end

    elseif (s._syntax == "http")
        for (key, values) = s._data
            if (typeof(values) == Array{AbstractString, 1})
                content *= "$key: $(join(values, ","))\n"
            else
                content *= "$key: $values\n"
            end
        end

    elseif (s._syntax == "simple")
        for (key, values) = s._data
            if (typeof(values) == Array{AbstractString, 1})
                content *= "$key $(join(values, ","))\n"
            else
                content *= "$key $values\n"
            end
        end

    else
        error("unknown syntax type: $(s._syntax)")
    end

    content
end # function craft_content

#----------
# remove entry from inside ini block
#----------
function erase!(s::ConfParse, block::String, key::String)
    block_key = getkey(s._data, block, nothing)
    if (block_key != nothing)
        if (haskey(s._data[block_key], key))
            delete!(s._data[block_key], key)
        end
    end

    s._is_modified = true
end # function erase!

#----------
# remove entry from config (outside of block if ini)
#----------
function erase!(s::ConfParse, key::String)
    if (haskey(s._data, key))
        delete!(s._data, key)
    end

    s._is_modified = true
end # function erase!

#----------
# for writing out new or modified configuration files
#----------
function save!(s::ConfParse, filename::Any = nothing)
    # if data has not been modified and a new file has not
    # been specified, don't write out
    if (s._is_modified == false) && (filename == nothing)
        return
    end

    # if there is no content to write out, don't create an
    # empty file
    content = craft_content(s)
    if (content == nothing)
        return
    end

    if (filename == nothing)
        s._fh = open_fh(s._filename, "w")
    else
        s._fh = open_fh(filename, "w")
    end

    write(s._fh, content)
    close(s._fh)                # if not closed, content is written when julia-session finishes
end # function save

#----------
# for retrieving data outside of a block
#----------
function retrieve(s::ConfParse, key::String)
    if (length(s._data[key]) == 1)
        return s._data[key][1]
    end

    s._data[key]
end # function retrieve

#----------
# for retrieving data from an ini config file block
#----------
function retrieve(s::ConfParse, block::String, key::String)
    if (length(s._data[block][key]) == 1)
        return s._data[block][key][1]
    end

    s._data[block][key]
end # function retrieve

#----------
# for retrieving data outside of a block and converting to type
#----------
function retrieve(s::ConfParse, key::String, t::Type)
    if (length(s._data[key]) == 1)
        return parse(t, s._data[key][1])
    end

    parse(t, s._data[key])
end # function retrieve

#----------
# for retrieving data from an ini config file block and converting to type
#----------
function retrieve(s::ConfParse, block::String, key::String, t::Type)
    if (length(s._data[block][key]) == 1)
        return parse(t, s._data[block][key][1])
    end

    parse(t, s._data[block][key])
end # function retrieve

#----------
# for inserting data in a config file
#----------
function commit!(s::ConfParse, key::String, value::Any)
    s._data[key]   = value
    s._is_modified = true
end # function commit!

#----------
# for inserting data inside an ini file block
#----------
function commit!(s::ConfParse, block::String, key::String, values::String)
    if (s._syntax != "ini")
        error("invalid setter function called for syntax type: $(s._syntax)")
    end

    s._data[block][key] = [values]
    s._is_modified      = true
end # function commit!

#----------
# for merging data of two configuration files
#----------
function merge!(s::ConfParse, t::ConfParse)
    for (block, key_values) = t._data
        if (!haskey(s._data, block))
            s._data[block] = key_values
        else
            merge!(s._data[block], key_values)
        end
    end
    s._is_modified = true
end # function merge!

end # module ConfParser


text:
for F in (:asin, :acos, :atan, :acsc, :asec, :acot)
    @eval begin
      $F(x::Double64) = Double64Float128($F, x)
      $F(x::Double32) = Double32Float128($F, x)
      $F(x::Double16) = Double16Float128($F, x)
    end
end

atan(y::Double64, x::Double64) = Double64(atan(Quadmath.Float128(y), Quadmath.Float128(x)))
atan(y::Double32, x::Double32) = Double32(atan(Quadmath.Float128(y), Quadmath.Float128(x)))
atan(y::Double16, x::Double16) = Double16(atan(Quadmath.Float128(y), Quadmath.Float128(x)))

#=
const k_pio2 = Double64(1.5707963267948966, 6.123233995736766e-17)

function atan(x::DoubleFloat{T}) where {T<:IEEEFloat}
   signbit(x) && return -atan(abs(x))
   isinf(x) && return copysign(DoubleFloat{T}(pi)/2, x)
   iszero(x) && return x
   y = DoubleFloat{T}(atan(x.hi))
   s, c = sin(y), cos(y)
   c2 = square(c)
   t = s/c
   z = y - (t - x)*c2 # z - (tan(x)-x)*(cos(x)^2)
   s, c = sin(z), cos(z)
   t = s/c
   z = z - (t - x)*c2 # z - (tan(x)-x)*(cos(x)^2)
   return z
end

function atan(y::DoubleFloat{T}, x::DoubleFloat{T}) where {T<:IEEEFloat}
   iszero(y) && iszero(x) && return x
   iszero(x) && return copysign(pi1o2(DoubleFloat{T}), y)
   iszero(y) && return signbit(x) ? pi1o1(DoubleFloat{T}) : zero(DoubleFloat{T})
   atanyx = atan(y/x)
   if x > 0
      atanyx
   elseif x < 0
      if y >= 0
         atanyx + T(pi)
      else
         atanyx - T(pi)
      end
   elseif y > 0
       T(pi)/2
   elseif y < 0
      -T(pi)/2
   else
      zero(T)
   end
end

function asin(x::DoubleFloat{T}) where {T<:IEEEFloat}
   (abs(x) > 1.0 || isinf(x)) && throw(DomainError("$x"))
   signbit(x) && return -asin(abs(x))
   y = x
   y = y / (1.0 + sqrt(1.0 - square(y)))
   z = atan(y)
   return DoubleFloat{T}(z.hi+z.hi, z.lo+z.lo)
end

function acos(x::DoubleFloat{T}) where {T<:IEEEFloat}
   (abs(x) > 1.0 || isinf(x)) && throw(DomainError("$x"))
   signbit(x) && return DoubleFloat{T}(onepi - acos(abs(x)))
   y = x
   y = sqrt(1.0 - square(y)) / (1.0 + y)
   z = atan(y)
   return DoubleFloat{T}(z.hi+z.hi, z.lo+z.lo)
end

function acsc(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isinf(x) && return copysign(zero(DoubleFloat{T}), x)
    return asin(inv(x))
end

function asec(x::DoubleFloat{T}) where {T<:IEEEFloat}
   isinf(x) && return DoubleFloat{T}(halfpi)
   return acos(inv(x))
end

function acot(x::DoubleFloat{T}) where {T<:IEEEFloat}
   isinf(x) && return copysign(zero(DoubleFloat{T}), x)
   signbit(x) && return DoubleFloat{T}(onepi - acot(abs(x)))
   iszero(x) && return k_pio2
   z = k_pio2 - atan(abs(x))
   if signbit(x.hi)
      z = -z
   end
   return z
end
=#


text:
# a. prints in singe line
print("This ")
print("is ")
print("sinle ")
print("line .\n\n")
#: This is sinle line .

# b. prints in multiple line
println("This ")
println("is ")
println("multiple ")
println("lines .\n\n")
#: This
#: is
#: multiple
#: lines .

# c. args
for iarg in ARGS
    println(iarg, " " ,typeof(iarg))
end
print("\n\n")
#: > julia file.jl 1 2 3 a
#: 1 String
#: 2 String
#: 3 String
#: a String


text:
mutable struct FactorGraph <: ANetwork
    nvars::Int
    nfacts::Int
    ne::Int
    edge_index_range::Int
    out_edges::Vector{Vector{Pair{Int,Int}}}  # Unordered adjlist
    epos::Vector{Pair{Int,Int}}      # Position of the edge in out_edges
    free_indexes::Vector{Int}       # Indexes of deleted edges
    props::PropertyStore
end

function FactorGraph(nvars::Int=0, nfacts::Int=0)
    n = nvars + nfacts
    out_edges = [Vector{Pair{Int,Int}}() for _=1:n]
    epos = Vector{Pair{Int,Int}}()
    free_indexes = Vector{Int}()
    return FactorGraph(nvars, nfacts, 0, 0, out_edges, epos, free_indexes, PropertyStore())
end

FactorGraph(nvars::Int, nfacts::Int, nedges::Int) = bipartite_erdos_renyi(nvars, nfacts, nedges)

edgetype(::Type{G}) where {G<:FactorGraph} = IndexedEdge
vertextype(::Type{G}) where {G<:FactorGraph} = Int

nv(g::FactorGraph) = length(g.out_edges)
ne(g::FactorGraph) = g.ne

variables(g::FactorGraph) = 1:g.nvars
factors(g::FactorGraph) = g.nvars+1:nv(g)

# TODO can be improved (see graph/digraph)
function edge(g::FactorGraph, i::Integer, j::Integer)
    (i > nv(g) || j > nv(g)) && return IndexedEdge(i, j, -1)
    oes = g.out_edges[i]
    pos = findfirst(e->e.first==j, oes)
    if pos !== nothing
        return IndexedEdge(i, j, oes[pos].second)
    else
        return IndexedEdge(i, j, -1)
    end
end

function out_edges(g::FactorGraph, i::Integer)
    oes = g.out_edges[i]
    return (IndexedEdge(i, j, idx) for (j, idx) in oes)
end

function out_neighbors(g::FactorGraph, i::Integer)
    oes = g.out_edges[i]
    return (j for (j, idx) in oes)
end

pop_vertex!(g::FactorGraph) = (clean_vertex!(g, nv(g)); pop!(g.out_edges); nv(g)+1)

# function add_vertex!(g::FactorGraph)
#     push!(g.out_edges, Vector{Pair{Int,Int}}())
#     return nv(g)
# end

function add_factor!(g::FactorGraph)
    push!(g.out_edges, Vector{Pair{Int,Int}}())
    g.nfacts += 1
    return nv(g)
end

function add_variable!(g::FactorGraph)
    push!(g.out_edges, Vector{Pair{Int,Int}}())
    g.nvars += 1
    swap_vertices!(g, g.nvars, nv(g))
    return g.nvars
end

function add_edge!(g::FactorGraph, u::Integer, v::Integer)
    (u ∈ variables(g) && v ∈ factors(g)) ||
    (v ∈ variables(g) && u ∈ factors(g)) ||
        throw(ArgumentError("Can add edge only between variable and factor."))
    u, v = u <= v ? (u, v) : (v, u)
    (u in vertices(g) && v in vertices(g)) || return (false, IndexedEdge(u,v,-1))
    has_edge(g, u, v) && return (false, IndexedEdge(u,v,-1)) # could be removed for multigraphs

    if isempty(g.free_indexes)
        g.edge_index_range += 1
        idx = g.edge_index_range
    else
        idx = pop!(g.free_indexes)
    end
    oes = g.out_edges[u]
    ies = g.out_edges[v]
    push!(oes, Pair(v, idx))
    if u != v
        push!(ies, Pair(u, idx))
    end
    g.ne += 1

    length(g.epos) < idx && resize!(g.epos, idx)
    g.epos[idx] = Pair(length(oes), length(ies))

    return (true, IndexedEdge(u,v,idx))
end

rem_edge!(g::FactorGraph, s::Integer, t::Integer) = rem_edge!(g, edge(g, s, t))

function rem_edge!(g::FactorGraph, e::IndexedEdge)
    s = e.src
    t = e.dst
    if s > t
        s,t = t,s
    end
    idx = e.idx
    idx <= 0 && return false
    oes = g.out_edges[s]
    ies = g.out_edges[t]
    idx > length(g.epos) && return false
    length(oes) == 0 && return false
    p1 = g.epos[idx].first
    p1 < 0 && return false

    back = last(oes)
    if back.first > s
        p2 = g.epos[back.second].second
        g.epos[back.second] = Pair(p1 , p2)
    elseif back.first == s #fix self-edges
        g.epos[back.second] = Pair(p1, p1)
    else
        p2 = g.epos[back.second].first
        g.epos[back.second] = Pair(p2 , p1)
    end
    oes[p1] = back
    pop!(oes)

    if s != t
        back = last(ies)
        p1 = g.epos[idx].second
        if back.first > t
            p2 = g.epos[back.second].second
            g.epos[back.second] = Pair(p1 , p2)
        elseif back.first == t #fix self-edges
            g.epos[back.second] = Pair(p1 , p1)
        else
            p2 = g.epos[back.second].first
            g.epos[back.second] = Pair(p2 , p1)
        end
        ies[p1] = back
        pop!(ies)
    end

    g.epos[idx] = Pair(-1,-1)

    g.ne -= 1
    push!(g.free_indexes, idx)
    return true
end

function in_edges(g::FactorGraph, i::Integer)
    ies = g.out_edges[i]
    return (IndexedEdge(j, i, idx) for (j, idx) in ies)
end

function swap_vertices!(g::FactorGraph, u::Integer, v::Integer)
    if u != v
        #TODO copying to avoid problems with self edges
        # maybe can copy only one of the two
        neigu = deepcopy(g.out_edges[u])
        neigv = deepcopy(g.out_edges[v])

        for (k,p) in enumerate(neigu)
            j, idx = p
            kj = j <= u ? g.epos[idx].first : g.epos[idx].second
            g.out_edges[j][kj] = Pair(v, idx)
            g.epos[idx] = j <= v ? Pair(kj, k) : Pair(k, kj)
        end

        for (k,p) in enumerate(neigv)
            j, idx = p
            kj = j <= v ? g.epos[idx].first : g.epos[idx].second
            g.out_edges[j][kj] = Pair(u, idx)
            g.epos[idx] = j <= u ? Pair(kj, k) : Pair(k, kj)
        end

        g.out_edges[u], g.out_edges[v] = g.out_edges[v], g.out_edges[u]

        swap_vertices!(g.props, u, v)
    end
end

function bipartite_erdos_renyi(n1::Int, n2::Int, p::Real; seed::Int=-1)
    m = n1 * n2
    rng = Erdos.getRNG(seed)
    nedg = rand(rng, Binomial(m, p))
    seed = seed >= 0 ? seed + 1 : seed
    return erdos_renyi(n1, n2, nedg; seed)
end

function bipartite_erdos_renyi(n1::Int, n2::Int, m::Int; seed::Int = -1)
    m > n1 * n2 && throw(ArgumentError("Cannot insert more than n1*n2 edges."))
    rng = Erdos.getRNG(seed)
    g = FactorGraph(n1, n2)
    while ne(g) < m
        source = rand(rng, 1:n1)
        dest = rand(rng, n1+1:n1+n2)
        res = add_edge!(g, source, dest)
    end
    return g
end

function random_bipartite_regular_graph(n1::Int, n2::Int, k::Int;
                                    seed::Int = -1, first_regular=true)

    if (first_regular && k > n2) || (!first_regular && k > n1)
        throw(ArgumentError("Not enough vertices for the required degree."))
    end

    rng = Erdos.getRNG(seed)
    g = FactorGraph(n1, n2)
    if first_regular
        for i in 1:n1
            while degree(g, i) < k
                j = rand(rng, n1+1:n1+n2)
                add_edge!(g, i, j)
            end
        end
    else
        for j in n1+1:n1+n2
            while degree(g, j) < k
                i = rand(rng, 1:n1)
                add_edge!(g, i, j)
            end
        end
    end
    return g
end


text:
using DataFrames
using CSV, DelimitedFiles

struct Generator
    Unit::Array{Float64}
    FixedCost::Array{Float64}
    VariableCost::Array{Float64}
    StartUpCost::Array{Float64}
    ShutdownCost::Array{Float64}
    MinPowerOut::Array{Float64}
    MaxPowerOut::Array{Float64}
    RampDownLimit::Array{Float64}
    RampShutDownLimit::Array{Float64}
    RampUpLimit::Array{Float64}
    RampStartUpLimit::Array{Float64}
    Uini::Array{Float64}
end

struct XGenerator
    Unit::Float64
    FixedCost::Float64
    VariableCost::Float64
    StartUpCost::Float64
    ShutdownCost::Float64
    MinPowerOut::Float64
    MaxPowerOut::Float64
    RampDownLimit::Float64
    RampShutDownLimit::Float64
    RampUpLimit::Float64
    RampStartUpLimit::Float64
    Uini::Float64
end

struct YGenerator
    Unit::Float64
    FixedCost::Float64
end

@timev dataGenerator = CSV.read("data_generators.csv", DataFrame)
@timev dataDemand = CSV.read("data_generators.csv", DataFrame)

println(df);

@timev data_generators = DelimitedFiles.readdlm("data_generators.csv", ',', header=true, Float64);
@timev data_demand = DelimitedFiles.readdlm("demand_reserves.csv", ',', header=true, Float64);

#println(data_generators)
#println(data_demand)

generators = data_generators[1]

#Generator 1
println("Size of generators[1,:]", size(generators) )

println(generators[1,:])
mydata =[1.0 5.0 0.1 20.0 0.5 50.0 350.0 300.0 300.0 200.0 200.0 0.0];
mydata2 =[1.0 5.0 0.1 20.0 0.5]
println("My data ",mydata[:, 1])

S = view(generators, 1:5)
println("S: ", S)

generator1 = YGenerator(1.0, 5.0)

println("Size of my data", size(mydata))
println("Size of my data2", size(mydata2))

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
using Oceananigans.Architectures: architecture, device_event, arch_array
using Oceananigans.AbstractOperations: KernelFunctionOperation
using Oceananigans.BuoyancyModels: ∂z_b
using Oceananigans.Operators: ℑzᵃᵃᶜ

struct ConvectiveAdjustmentVerticalDiffusivity{TD, CK, CN, BK, BN} <: AbstractScalarDiffusivity{TD, VerticalFormulation}
    convective_κz :: CK
    convective_νz :: CN
    background_κz :: BK
    background_νz :: BN

    function ConvectiveAdjustmentVerticalDiffusivity{TD}(convective_κz::CK,
                                                         convective_νz::CN,
                                                         background_κz::BK,
                                                         background_νz::BN) where {TD, CK, CN, BK, BN}

        return new{TD, CK, CN, BK, BN}(convective_κz, convective_νz, background_κz, background_νz)
    end
end

"""
    ConvectiveAdjustmentVerticalDiffusivity([time_discretization = VerticallyImplicitTimeDiscretization(), FT=Float64;]
                                            convective_κz = 0,
                                            convective_νz = 0,
                                            background_κz = 0,
                                            background_νz = 0)

Return a convective adjustment vertical diffusivity closure that applies different values of diffusivity and/or viscosity depending
whether the region is statically stable (positive or zero buoyancy gradient) or statically unstable (negative buoyancy gradient).

Arguments
=========

* `time_discretization`: Either `ExplicitTimeDiscretization()` or `VerticallyImplicitTimeDiscretization()`;
                         default `VerticallyImplicitTimeDiscretization()`.

* `FT`: Float type; default `Float64`.

Keyword arguments
=================

* `convective_κz`: Vertical tracer diffusivity in regions with negative (unstable) buoyancy gradients. Either
                   a single number, function, array, field, or tuple of diffusivities for each tracer.

* `background_κz`: Vertical tracer diffusivity in regions with zero or positive (stable) buoyancy gradients.

* `convective_νz`: Vertical viscosity in regions with negative (unstable) buoyancy gradients. Either
                  a number, function, array, or field.

* `background_κz`: Vertical viscosity in regions with zero or positive (stable) buoyancy gradients.

Example
=======

```jldoctest
julia> using Oceananigans

julia> cavd = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1)
ConvectiveAdjustmentVerticalDiffusivity{VerticallyImplicitTimeDiscretization}(background_κz=0.0 convective_κz=1 background_νz=0.0 convective_νz=0.0)
```
"""
function ConvectiveAdjustmentVerticalDiffusivity(time_discretization = VerticallyImplicitTimeDiscretization(), FT = Float64;
                                                 convective_κz = zero(FT),
                                                 convective_νz = zero(FT),
                                                 background_κz = zero(FT),
                                                 background_νz = zero(FT))

    return ConvectiveAdjustmentVerticalDiffusivity{typeof(time_discretization)}(convective_κz, convective_νz,
                                                                                background_κz, background_νz)
end

ConvectiveAdjustmentVerticalDiffusivity(FT::DataType; kwargs...) = ConvectiveAdjustmentVerticalDiffusivity(VerticallyImplicitTimeDiscretization(), FT; kwargs...)

const CAVD = ConvectiveAdjustmentVerticalDiffusivity

#####
##### Diffusivity field utilities
#####

# Support for "ManyIndependentColumnMode"
const CAVDArray = AbstractArray{<:CAVD}
const FlavorOfCAVD = Union{CAVD, CAVDArray}

with_tracers(tracers, closure::FlavorOfCAVD) = closure
DiffusivityFields(grid, tracer_names, bcs, closure::FlavorOfCAVD) = (; κ = ZFaceField(grid), ν = ZFaceField(grid))
@inline viscosity_location(::FlavorOfCAVD) = (Center(), Center(), Face())
@inline diffusivity_location(::FlavorOfCAVD) = (Center(), Center(), Face())
@inline viscosity(::FlavorOfCAVD, diffusivities) = diffusivities.ν
@inline diffusivity(::FlavorOfCAVD, diffusivities, id) = diffusivities.κ

function calculate_diffusivities!(diffusivities, closure::FlavorOfCAVD, model)

    arch = model.architecture
    grid = model.grid
    tracers = model.tracers
    buoyancy = model.buoyancy

    event = launch!(arch, grid, :xyz,
                    ## If we can figure out how to only precompute the "stability" of a cell:
                    # compute_stability!, diffusivities, grid, closure, tracers, buoyancy,
                    compute_convective_adjustment_diffusivities!, diffusivities, grid, closure, tracers, buoyancy,
                    dependencies = device_event(arch))

    wait(device(arch), event)

    return nothing
end

@inline is_stableᶜᶜᶠ(i, j, k, grid, tracers, buoyancy) = ∂z_b(i, j, k, grid, buoyancy, tracers) >= 0

@kernel function compute_convective_adjustment_diffusivities!(diffusivities, grid, closure, tracers, buoyancy)
    i, j, k, = @index(Global, NTuple)

    # Ensure this works with "ensembles" of closures, in addition to ordinary single closures
    closure_ij = getclosure(i, j, closure)

    stable_cell = is_stableᶜᶜᶠ(i, j, k, grid, tracers, buoyancy)

    @inbounds diffusivities.κ[i, j, k] = ifelse(stable_cell,
                                                closure_ij.background_κz,
                                                closure_ij.convective_κz)

    @inbounds diffusivities.ν[i, j, k] = ifelse(stable_cell,
                                                closure_ij.background_νz,
                                                closure_ij.convective_νz)
end

#=
## If we can figure out how to only precompute the "stability" of a cell:
@kernel function compute_stability!(diffusivities, grid, closure, tracers, buoyancy)
    i, j, k, = @index(Global, NTuple)
    @inbounds diffusivities.unstable_buoyancy_gradient[i, j, k] = is_unstableᶜᶜᶠ(i, j, k, grid, tracers, buoyancy)
end
=#

#####
##### Show
#####

function Base.summary(closure::ConvectiveAdjustmentVerticalDiffusivity{TD}) where TD
    return string("ConvectiveAdjustmentVerticalDiffusivity{$TD}" *
        "(background_κz=", prettysummary(closure.background_κz), " convective_κz=", prettysummary(closure.convective_κz),
        " background_νz=", prettysummary(closure.background_νz), " convective_νz=", prettysummary(closure.convective_νz), ")")
end

Base.show(io::IO, closure::ConvectiveAdjustmentVerticalDiffusivity) = print(io, summary(closure))


text:
using .Plots

import .Plots.plot
"""
    plot(anova::AnovaData)

Creates factor plots for ANOVA data.

Enables inspection of ANOVA for the effect size of interactions and relationships between factor levels.

Import Plots package to use this.
"""
function plot(anova::AnovaData)
    nfactors = ndims(anova.crossedcellmeans)
    if nfactors == 1
        error("Not useful to plot just one factor.")
    end
    nfactorlevels = size(anova.crossedcellmeans)
    plots = Array{Plots.Plot, 2}(undef, nfactors, nfactors)
    for i in 1:nfactors
        otherfactorindexes = (1:nfactors)[Not(i)]
        factormidpoint = (1 + nfactorlevels[i]) / 2
        factormeans = [mean(anova.crossedcellmeans, dims = (1:nfactors)[Not(i)]) |> vec for i = 1:nfactors]
        for jindex in 1:(nfactors - 1)
            j = otherfactorindexes[jindex]
            pairfactormeans = mean(anova.crossedcellmeans, dims = otherfactorindexes[Not(jindex)])
            plots[i,j] = plot([selectdim(pairfactormeans, j, k) |> vec for k = 1:nfactorlevels[j]], legend = false)
            scatter!(factormeans[i], markershape = :+, markercolor = :gray)
            scatter!(repeat([factormidpoint], nfactorlevels[j]), factormeans[j], markershape = :x, markercolor = :gray)
        end
        plots[i,i] = plot()

    end
    plot(plots..., layout = (nfactors, nfactors))
end


text:
# ----------------------------------------------------
# --- Utils.jl
# ----------------------------------------------------
# We put here common usefull stuff for SDR managment

# --- Container for Radio use
# We will have functions from different origin and different supported keywords
# This function parse the input keywords and returns the ones supported by the function, listed in iteration
function parseKeyword(kwargs,iteration)
    # We populate a new dictionnary based on the input keywords and the supported ones
    # In order not to create keywords that are not supported (i.e leaving the default value)
    # we only evaluate the keywords defined both in the 2 dictionnaries
    # This means that the default fallback should never happen
    kwargs = Dict(key=>get(kwargs,key,0) for key in intersect(iteration,keys(kwargs)))
    return kwargs
end

# --- Conversion
# Adalm Pluto structure is based on Int parameters, and AbstractSDRs use massively Float. We need to convert just before dispatching. As some parameter may be float (as gain) we should round before conversion. The following function does that.
_toInt(x) = Int(round(x));


text:
using EconPDEs, Distributions

struct AchdouHanLasryLionsMollModel
    # income process parameters
    κy::Float64
    ybar::Float64
    σy::Float64

    r::Float64

    # utility parameters
    ρ::Float64
    γ::Float64

    amin::Float64
    amax::Float64
end

function AchdouHanLasryLionsMollModel(;κy = 0.1, ybar = 1.0, σy = 0.07, r = 0.03, ρ = 0.05, γ = 2.0, amin = 0.0, amax = 500.0)
    AchdouHanLasryLionsMollModel(κy, ybar, σy, r, ρ, γ, amin, amax)
end

function initialize_stategrid(m::AchdouHanLasryLionsMollModel; yn = 5, an = 50)
    κy = m.κy ; ybar = m.ybar ; σy = m.σy  ; ρ = m.ρ ; γ = m.γ ; amin = m.amin ; amax = m.amax

    distribution = Gamma(2 * κy * ybar / σy^2, σy^2 / (2 * κy))
    ymin = quantile(distribution, 0.001)
    ymax = quantile(distribution, 0.999)
    ys = range(ymin, stop = ymax, length = yn)
    as = range(amin, stop = amax, length = an)
    OrderedDict(:y => ys, :a => as)
end

function initialize_y(m::AchdouHanLasryLionsMollModel, stategrid)
    OrderedDict(:v => [log(y + a) for y in stategrid[:y], a in stategrid[:a]])
end

function (m::AchdouHanLasryLionsMollModel)(state::NamedTuple, value::NamedTuple)
    κy = m.κy ; σy = m.σy ; ybar = m.ybar ; r = m.r ; ρ = m.ρ ; γ = m.γ ; amin = m.amin ; amax = m.amax
    y, a = state.y, state.a
    v, vy, va, vyy, vya, vaa = value.v, value.vy, value.va, value.vyy, value.vya, value.vaa
    μy = κy * (ybar - y)
    va = max(va, eps())

    # There is no second derivative at 0 so just specify first order derivative
    c = va^(-1 / γ)
    μa = y + r * a - c
    if (a ≈ amin) && (μa <= 0.0)
        va = (y + r * amin)^(-γ)
        c = y + r * amin
        μa = 0.0
    end

    vt = c^(1 - γ) / (1 - γ) + μa * va + μy * vy + 0.5 * vyy * σy^2 - ρ * v
    return (vt,), (μy, μa), (v = v, c = c, va = va, vy = vy, y = y, a = a, μa = μa, vaa = vaa)
end

m = AchdouHanLasryLionsMollModel()
stategrid = initialize_stategrid(m)
y0 = initialize_y(m, stategrid)
y, result, distance = pdesolve(m, stategrid, y0)

# Check marginal value of wealth converges to 1.0 at infinity
#b = ((m.r + (m.ρ - m.r)/m.γ))^(1/(1 - 1/m.γ))
#pw = (result[:v] * (1-m.γ)).^(1/(1-m.γ)-1) .* result[:va] ./ b


text:
export
    DebugEnvDeterministicSingleStep,
    reset,
    step,
    observation_space_spec,
    action_space_spec

#=
Description:
    Takes a single step from initial state 0.
    Actions are 1 or 2.
    Action 1 yields 1 reward.
    Action 2 yields 2 reward.
    Terminates after the single action.
=#
mutable struct DeterministicSingleStepDebugEnv <: Env
    function DeterministicSingleStepDebugEnv(params::Dict=Dict())
        return new()
    end
end
Base.reset(env::DeterministicSingleStepDebugEnv) = 0
function Base.step(env::DeterministicSingleStepDebugEnv, action::Int)
    @assert action == 1 || action == 2
    return (0, action, true, Dict())
end
observation_space_spec(env::DeterministicSingleStepDebugEnv) = (1,), "Discrete"
action_space_spec(env::DeterministicSingleStepDebugEnv) = (2,), "Discrete"
render(env::DeterministicSingleStepDebugEnv) = nothing


text:
using DtPlot

#using Winston
using Winston
using PosixCalendar
using Predictors
using DeepThought.DTPredictors
using DeepThought.DTPDFS
using Verification
using DataFrames
using Datetime
import DeepThought.DTPDFS.PDF

REPOSITORY = joinpath(Pkg.dir("DeepThoughtTestData"), "data-root")
STATION = "93439"
FVAR = "TTTTT"
CANDIDATE = "vbagging_TTTTTnz8km_TTTTT"

OBS = Observation(REPOSITORY,STATION,"TTTTT","target","observation")
F = PDFForecast(REPOSITORY,STATION,FVAR, CANDIDATE)

hourinsecs = 60*60
dayinsecs = 24*3600
basestep = 12*hourinsecs
range = (0, 48 * hourinsecs)
# Example 1 --------------------------------------------------------------------
period = Request((Time(2013,7,1,12),
                  Time(2013,7,1,12),
                  basestep),
                 range)
#Monday 1st July 2013 12:00:00 AM

hint(OBS, period)
obs = getdata(OBS, period)

vts = None
pps = None
vals = None
bts = sort(collect(keys(obs)))

# Pdf forecast
df = prepareDataFramePdf(F, period)
fp = leaf(df, CANDIDATE)
display(fp)

# Obs
for bt = bts
    pps = sort(collect(keys(obs[bt])))
    vts = sort([bt + pp for pp = pps])
    vals = [obs[bt][pp] for pp = pps]
    pps = [pp.seconds for pp in pps]
    vts = [vt.seconds/3600 for vt in vts]
    #plot(vts, vals,"k")
end

obsp = Curve(vts, vals, color="grey")
add(fp, obsp)
#-------------------------------------------------------------------------------

# Example 2 --------------------------------------------------------------------
period = Request((Time(2013,7,1,12),
                  Time(2013,8,5,12),
                  basestep),
                 range)

df = prepareDataFramePdf(F, period)
figure()
fp2 = FramedPlot()
#figure()
colors = ["red", "blue", "green", "magenta", "orange"]
cnt = 0
for fc in groupby(df, :basetime)
    cnt += 1
    x = convert(Array{Float64,1}, fc[:x3])
    y = convert(Array{Float64,1}, fc[:median])

    color = colors[1+cnt%length(colors)]
    p = Points(x, y)
    style(p, kind="filled circle", size=0.3, color=color)

    c = Curve(x, y, color=color)
    add(fp2, p)
    add(fp2, c)
    #oplot(x,y)
    ## errorbar(fc[:x3], fc[:median], yerr=[fc[:q25], fc[:q75]],
    ##          fmt="o", ecolor='g', capthick=2)
end

# Obs
obs = getdata(OBS, period)
vts = None
pps = None
vals = None
bts = sort(collect(keys(obs)))
for bt = bts
    pps = sort(collect(keys(obs[bt])))
    vts = sort([bt + pp for pp = pps])
    vals = [obs[bt][pp] for pp = pps]
    pps = [pp.seconds for pp in pps]
    vts = [vt.seconds/3600 for vt in vts]
    #plot(vts, vals,"k")
end

obsp = Points(vts, vals, color="black")
style(obsp, kind="circle", size=0.3)

obsc = Curve(vts, vals, color="grey")
#fp2 = FramedPlot()
add(fp2, obsp)
add(fp2, obsc)

display(fp2)

## # adjust xbounds
## setattr(fp2, "xrange",(vts[1]-20, vts[end]))
## display(fp2)

# Example 2b -------------------------------------------------------------------
fp2a = FramedPlot(fp2)
# FramedPlot()

# adjust xbounds
setattr(fp2a, "xrange",(vts[1]-20, vts[end]))
display(fp2a)

# Example 3 --------------------------------------------------------------------
period = Request((Time(2013,7,1,12),
                  Time(2013,7,5,12),
                  basestep*2),
                 range)
figure()
df = prepareDataFramePdf(F, period)
fp3 = FramedPlot(title = "Basehour = 12Z")
#figure()
colors = ["red", "blue", "green", "magenta", "orange"]
cnt = 0
for fc in groupby(df, :basetime)
    cnt += 1
    x = convert(Array{Float64,1}, fc[:x2])
    y = convert(Array{Float64,1}, fc[:median])
    #println(x, "\n",y)
    println(cnt%length(colors))
    color = colors[1+cnt%length(colors)]
    p = Points(x, y)
    setattr(p, label=fc[:basetime][1])
    style(p, kind="filled circle", size=0.3, color=color)

    #l = Legend(.1, .9, {a,b,s})

    c = Curve(x, y, color=color)
    add(fp3, p)
    add(fp3, c)
    #oplot(x,y)
    ## errorbar(fc[:x3], fc[:median], yerr=[fc[:q25], fc[:q75]],
    ##          fmt="o", ecolor='g', capthick=2)
end

display(fp3)

# Example 4 --------------------------------------------------------------------
period = Request((Time(2013,7,1,12),
                  Time(2013,7,5,12),
                  basestep),
                 (0,basestep))
figure()
df = prepareDataFramePdf(F, period)
fp4 = FramedPlot(title = "Overlapping regions prognoses omitted")
#figure()
colors = ["red", "blue", "green", "magenta", "orange"]
cnt = 0
for fc in groupby(df, :basetime)
    cnt += 1
    x = convert(Array{Float64,1}, fc[:x3])
    y = convert(Array{Float64,1}, fc[:median])
    #println(x, "\n",y)
    println(cnt%length(colors))
    color = colors[1+cnt%length(colors)]
    p = Points(x, y)
    setattr(p, label=fc[:basetime][1])
    style(p, kind="filled circle", size=0.3, color=color)

    #l = Legend(.1, .9, {a,b,s})

    c = Curve(x, y, color=color)
    add(fp4, p)
    add(fp4, c)
    #oplot(x,y)
    ## errorbar(fc[:x3], fc[:median], yerr=[fc[:q25], fc[:q75]],
    ##          fmt="o", ecolor='g', capthick=2)
end

display(fp4)


text:
import .Utils: capacity
import .Utils: update!

include("action_selectors/action_selectors.jl")
include("approximators/approximators.jl")
include("buffers/buffers.jl")
include("environment_models/environment_models.jl")
include("learners/learners.jl")
include("policies/policies.jl")
include("agents/agents.jl")
include("preprocessors.jl")


text:
using WebIO
using WebSockets
using Test

@testset "HTTP provider" begin
    output = sprint(io-> show(io, WebIO.WEBIO_APPLICATION_MIME(), node(:div, "hello, world")))
    @test occursin("_webIOWebSocketURL = ", output)
    @test occursin("ws://127.0.0.1:8081/webio_websocket/", output)
    @test occursin("""hello, world""", output)
    @test WebIO.webio_server_config[] == (url = "127.0.0.1", http_port = 8081, ws_url = "ws://127.0.0.1:8081/webio_websocket/")
    @test isassigned(WebIO.singleton_instance)
end

using WebSockets, WebIO

app = Ref{Any}(node(:div, "hi"))
function serve_app(req)
    if req.target == "/"
        return sprint() do io
            print(io,
                """
                <!doctype html>
                <html>
                <head>
                <meta charset="UTF-8">
                </head>
                <body>
                """
            )
            show(io, WebIO.WEBIO_APPLICATION_MIME(), app[])
            print(io,
                """
                </body>
                </html>
                """
            )
        end
    else
        return missing
    end
end
server = WebIO.WebIOServer(serve_app, logger = stdout, verbose = true)
server.serve_task

w = Scope()

obs = Observable(w, "rand-value", 0.0)

on(obs) do x
   println("JS sent $x")
end

app[] = w(
 dom"button"(
   "generate random",
   events=Dict("click" => js"""function() { _webIOScope.setObservableValue('rand-value', Math.random()); }"""),
 ),
);


text:
"""creates outflow variables specified in data"""
function variable_inflow(sp, data::Dict)
    @variables(sp, begin
        inflow[r=1:data["hydro"]["nHyd"]]
    end)
end

# TODO: add data["hydro"]["Hydrogenerators"][r]["min_turn"] as penalized constraint
"""creates outflow variables specified in data"""
function variable_outflow(sp, data::Dict)
    @variables(sp, begin
        0 <= outflow[r=1:data["hydro"]["nHyd"]] <= data["hydro"]["Hydrogenerators"][r]["max_turn"]
    end)
end

"""creates spillage variables specified in data"""
function variable_spillage(sp, data::Dict)
    @variables(sp, begin
        spill[r=1:data["hydro"]["nHyd"]]      >= 0
    end)
end

# TODO: add data["hydro"]["Hydrogenerators"][r]["min_volume"] as penalized constraint
"""creates volume variables specified in data"""
function variable_volume(sp, data::Dict)
    @variable(sp, 0 <= reservoir[r = 1:data["hydro"]["nHyd"]] <= data["hydro"]["Hydrogenerators"][r]["max_volume"], SDDP.State, initial_value = data["hydro"]["Hydrogenerators"][r]["initial_volume"])
end

"""creates deficit variables"""
function variable_deficit(sp, data::Dict, pm::AbstractPowerModel)
    PowerModels.var(pm, 0, 1)[:deficit] = @variable(sp, deficit[i in collect(1:length(PowerModels.con(pm, 0, 1, :kcl_p)))] >= 0)
end

"""creates dict of cost variables"""
function variable_cost(sp, data::Dict)
    cost = Dict()
    sp.ext[:cost] = cost
end


text:
include("utils.jl")
@info "Manifolds.jl Test settings:\n\n" *
      "Testing Float32:  $(TEST_FLOAT32)\n" *
      "Testing Double64: $(TEST_DOUBLE64)\n" *
      "Testing Static:   $(TEST_STATIC_SIZED)\n\n" *
      "Check test/utils.jl if you wish to change these settings."

@testset "Manifolds.jl" begin
    include_test("differentiation.jl")

    @testset "Ambiguities" begin
        # TODO: reduce the number of ambiguities for Julia 1.6
        if VERSION.prerelease == () && !Sys.iswindows() && VERSION < v"1.6.0"
            @test length(Test.detect_ambiguities(ManifoldsBase)) <= 58
            @test length(Test.detect_ambiguities(Manifolds)) == 0
            @test length(our_base_ambiguities()) <= 24
        else
            @info "Skipping Ambiguity tests for pre-release versions"
        end
    end

    @testset "utils test" begin
        Random.seed!(42)
        @testset "usinc_from_cos" begin
            @test Manifolds.usinc_from_cos(-1) == 0
            @test Manifolds.usinc_from_cos(-1.0) == 0.0
        end
        @testset "log_safe!" begin
            n = 8
            Q = qr(randn(n, n)).Q
            A1 = Matrix(Hermitian(Q * Diagonal(rand(n)) * Q'))
            @test exp(Manifolds.log_safe!(similar(A1), A1)) ≈ A1 atol = 1e-6
            A1_fail = Matrix(Hermitian(Q * Diagonal([-1; rand(n - 1)]) * Q'))
            @test_throws DomainError Manifolds.log_safe!(similar(A1_fail), A1_fail)

            T = triu!(randn(n, n))
            T[diagind(T)] .= rand.()
            @test exp(Manifolds.log_safe!(similar(T), T)) ≈ T atol = 1e-6
            T_fail = copy(T)
            T_fail[1] = -1
            @test_throws DomainError Manifolds.log_safe!(similar(T_fail), T_fail)

            A2 = Q * T * Q'
            @test exp(Manifolds.log_safe!(similar(A2), A2)) ≈ A2 atol = 1e-6
            A2_fail = Q * T_fail * Q'
            @test_throws DomainError Manifolds.log_safe!(similar(A2_fail), A2_fail)

            A3 = exp(randn(n, n))
            @test exp(Manifolds.log_safe!(similar(A3), A3)) ≈ A3 atol = 1e-6

            A3_fail = Float64[1 2; 3 1]
            @test_throws DomainError Manifolds.log_safe!(similar(A3_fail), A3_fail)

            A4 = randn(ComplexF64, n, n)
            @test exp(Manifolds.log_safe!(similar(A4), A4)) ≈ A4 atol = 1e-6
        end
        @testset "isnormal" begin
            @test !Manifolds.isnormal([1.0 2.0; 3.0 4.0])
            @test !Manifolds.isnormal(complex.(reshape(1:4, 2, 2), reshape(5:8, 2, 2)))

            # diagonal
            @test Manifolds.isnormal(diagm(randn(5)))
            @test Manifolds.isnormal(diagm(randn(ComplexF64, 5)))
            @test Manifolds.isnormal(Diagonal(randn(5)))
            @test Manifolds.isnormal(Diagonal(randn(ComplexF64, 5)))

            # symmetric/hermitian
            @test Manifolds.isnormal(Symmetric(randn(3, 3)))
            @test Manifolds.isnormal(Hermitian(randn(3, 3)))
            @test Manifolds.isnormal(Hermitian(randn(ComplexF64, 3, 3)))
            x = Matrix(Symmetric(randn(3, 3)))
            x[3, 1] += eps()
            @test !Manifolds.isnormal(x)
            @test Manifolds.isnormal(x; atol=sqrt(eps()))

            # skew-symmetric/skew-hermitian
            skew(x) = x - x'
            @test Manifolds.isnormal(skew(randn(3, 3)))
            @test Manifolds.isnormal(skew(randn(ComplexF64, 3, 3)))

            # orthogonal/unitary
            @test Manifolds.isnormal(Matrix(qr(randn(3, 3)).Q); atol=sqrt(eps()))
            @test Manifolds.isnormal(
                Matrix(qr(randn(ComplexF64, 3, 3)).Q);
                atol=sqrt(eps()),
            )
        end
        @testset "realify/unrealify!" begin
            # round trip real
            x = randn(3, 3)
            @test Manifolds.realify(x, ℝ) === x
            @test Manifolds.unrealify!(similar(x), x, ℝ) == x

            # round trip complex
            x2 = randn(ComplexF64, 3, 3)
            x2r = Manifolds.realify(x2, ℂ)
            @test eltype(x2r) <: Real
            @test size(x2r) == (6, 6)
            x2c = Manifolds.unrealify!(similar(x2), x2r, ℂ)
            @test x2c ≈ x2

            # matrix multiplication is preserved
            x3 = randn(ComplexF64, 3, 3)
            x3r = Manifolds.realify(x3, ℂ)
            @test x2 * x3 ≈ Manifolds.unrealify!(similar(x2), x2r * x3r, ℂ)
        end
    end

    include_test("groups/group_utils.jl")
    include_test("notation.jl")
    # starting with tests of simple manifolds
    include_test("centered_matrices.jl")
    include_test("circle.jl")
    include_test("cholesky_space.jl")
    include_test("elliptope.jl")
    include_test("euclidean.jl")
    include_test("fixed_rank.jl")
    include_test("generalized_grassmann.jl")
    include_test("generalized_stiefel.jl")
    include_test("grassmann.jl")
    include_test("hyperbolic.jl")
    include_test("lorentz.jl")
    include_test("multinomial_doubly_stochastic.jl")
    include_test("multinomial_symmetric.jl")
    include_test("positive_numbers.jl")
    include_test("probability_simplex.jl")
    include_test("projective_space.jl")
    include_test("rotations.jl")
    include_test("skewhermitian.jl")
    include_test("spectrahedron.jl")
    include_test("sphere.jl")
    include_test("sphere_symmetric_matrices.jl")
    include_test("stiefel.jl")
    include_test("symmetric.jl")
    include_test("symmetric_positive_definite.jl")
    include_test("symmetric_positive_semidefinite_fixed_rank.jl")

    include_test("essential_manifold.jl")
    include_test("multinomial_matrices.jl")
    include_test("oblique.jl")
    include_test("torus.jl")

    #meta manifolds
    include_test("product_manifold.jl")
    include_test("power_manifold.jl")
    include_test("vector_bundle.jl")
    include_test("graph.jl")

    include_test("metric.jl")
    include_test("statistics.jl")
    include_test("approx_inverse_retraction.jl")

    # Lie groups and actions
    include_test("groups/groups_general.jl")
    include_test("groups/array_manifold.jl")
    include_test("groups/circle_group.jl")
    include_test("groups/translation_group.jl")
    include_test("groups/general_linear.jl")
    include_test("groups/special_linear.jl")
    include_test("groups/special_orthogonal.jl")
    include_test("groups/product_group.jl")
    include_test("groups/semidirect_product_group.jl")
    include_test("groups/special_euclidean.jl")
    include_test("groups/group_operation_action.jl")
    include_test("groups/rotation_action.jl")
    include_test("groups/translation_action.jl")
    include_test("groups/metric.jl")

    include_test("recipes.jl")
end


text:
using Test
using Distributed
using KronLinInv

# get all the functions
include("test_suite.jl")

nwor = nworkers()

@testset "Tests " begin
    println("\n Number of workers available: $nwor")
    println()

    printstyled("Testing 2D example \n", bold=true,color=:cyan)
    @test test2D()

    printstyled("Testing 3D example \n", bold=true,color=:cyan)
    @test test3D()

    println()
end


text:
using Dates, Test
using MultiStartLH

mdl = MultiStartLH;

function history_test()
    @testset "History" begin
        fPath = mdl.test_history_path();
        nPoints = 5;
        nParams = 4;
        h = mdl.make_test_history(nPoints, nParams);
        @test h isa History;

        solnV = collect(1 : nParams) ./ 2;
        startV = solnV ./ 2;
        fVal = 0.8;
        p1 = Point(startV, solnV, fVal + 1.0, fVal, :a, Second(210), 35);
        mdl.add_to_history!(h, 1, p1);
        @test n_solved(h) == 1;

        save_history(h);
        h2 = load_history(fPath);
        @test n_solved(h) == 1;
        p11 = h2.pointV[1];
        @test all(p11.solnV .≈ solnV);

        p2 = Point(startV .+ 0.2, solnV .+ 0.2,
            fVal + 1.1, fVal - 0.1, :b, Second(120), 45);
        mdl.add_to_history!(h, 2, p2);
        p3 = Point(startV .+ 0.1, solnV .+ 0.1,
            fVal + 2.1, fVal + 0.1, :c, Second(90), 66);
        mdl.add_to_history!(h, 3, p3);
        @test n_solved(h) == 3;
        pBest = best_point(h);
        @test all(isapprox.(pBest.solnV, solnV .+ 0.2));
        @test isapprox(pBest.fVal, fVal - 0.1);
        @test pBest.exitFlag == :b;
    end
end

@testset "History All" begin
    history_test();
end

# ------------------


text:
module DGFEOperatorsTests

using Test
using Gridap
import Gridap: ∇

u(x) = x[1] + x[2]
∇u(x) = VectorValue(1.0,1.0)
∇(::typeof(u)) = ∇u
f(x) = 0.0
g(x) = u(x)

L = 1.0
limits = (0.0, L, 0.0, L)
ncellx = 4
model = CartesianDiscreteModel(domain=limits, partition=(ncellx,ncellx))

h = L / ncellx

γ = 10

order = 3
fespace = DLagrangianFESpace(Float64,model,order)

V = TestFESpace(fespace)
U = TrialFESpace(fespace)

trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2*order)

btrian = BoundaryTriangulation(model)
bquad = CellQuadrature(btrian,degree=2*order)
nb = NormalVector(btrian)

strian = SkeletonTriangulation(model)
squad = CellQuadrature(strian,degree=2*order)
ns = NormalVector(strian)

a_Ω(v,u) = inner(∇(v), ∇(u))
b_Ω(v) = inner(v,f)
t_Ω = AffineFETerm(a_Ω,b_Ω,trian,quad)

a_∂Ω(v,u) = (γ/h) * inner(v,u) - inner(v, ∇(u)*nb ) - inner(∇(v)*nb, u)
b_∂Ω(v) = (γ/h) * inner(v,g) - inner(∇(v)*nb, g)
t_∂Ω = AffineFETerm(a_∂Ω,b_∂Ω,btrian,bquad)

a_Γ(v,u) = (γ/h) * inner( jump(v*ns), jump(u*ns)) - inner( jump(v*ns), mean(∇(u)) ) - inner( mean(∇(v)), jump(u*ns) )
t_Γ = LinearFETerm(a_Γ,strian,squad)

op = LinearFEOperator(V,U,t_Ω,t_∂Ω,t_Γ)

uh = solve(op)

e = u - uh

#writevtk(trian,"trian",cellfields=["uh"=>uh,"e"=>e])

l2(u) = inner(u,u)
h1(u) = a_Ω(u,u) + l2(u)

el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))

@test el2 < 1.e-8
@test eh1 < 1.e-8

end # module


text:
GreetString = "Hello     World"
grtMat = fill('A',3, 5)#creating a char matrix
k=1
#filling up
for i in 1:3, j in 1:5
    grtMat[i,j] = GreetString[k]
    k+=1
end
display(grtMat)

max_stars_repo_path: Julia/Matrix Hello World.jl

max_stars_repo_name: saurabhcommand/Hello-world

max_issues_repo_path: Julia/Matrix Hello World.jl

max_issues_repo_name: saurabhcommand/Hello-world

max_forks_repo_path: Julia/Matrix Hello World.jl

max_forks_repo_name: saurabhcommand/Hello-world


text:
using Distributed, Test, CUDA

# GPUArrays has a testsuite that isn't part of the main package.
# Include it directly.
import GPUArrays
gpuarrays = pathof(GPUArrays)
gpuarrays_root = dirname(dirname(gpuarrays))
include(joinpath(gpuarrays_root, "test", "testsuite.jl"))
testf(f, xs...; kwargs...) = TestSuite.compare(f, CuArray, xs...; kwargs...)

using Random

# detect cuda-memcheck
const memcheck = haskey(ENV, "CUDA_MEMCHECK")

# precompile the runtime library
CUDA.precompile_runtime()

# for when we include tests directly
CUDA.allowscalar(false)

## entry point

function runtests(f, name, time_source=:cuda, snoop=nothing)
    old_print_setting = Test.TESTSET_PRINT_ENABLE[]
    Test.TESTSET_PRINT_ENABLE[] = false

    if snoop !== nothing
        io = open(snoop, "w")
        ccall(:jl_dump_compiles, Nothing, (Ptr{Nothing},), io.handle)
    end

    try
        # generate a temporary module to execute the tests in
        mod_name = Symbol("Test", rand(1:100), "Main_", replace(name, '/' => '_'))
        mod = @eval(Main, module $mod_name end)
        @eval(mod, using Test, Random, CUDA)

        let id = myid()
            wait(@spawnat 1 print_testworker_started(name, id))
        end

        ex = quote
            Random.seed!(1)
            CUDA.allowscalar(false)

            if $(QuoteNode(time_source)) == :cuda
                CUDA.@timed @testset $"$name" begin
                    $f()
                end
            elseif $(QuoteNode(time_source)) == :julia
                res = @timed @testset $"$name" begin
                    $f()
                end
                res..., 0, 0, 0
            else
                error($"Unknown time source: ",$(QuoteNode(time_source)))
            end
        end
        data = Core.eval(mod, ex)
        #data[1] is the testset

        # process results
        cpu_rss = Sys.maxrss()
        gpu_rss = if has_nvml()
            cuda_dev = device()
            nvml_dev = NVML.Device(uuid(cuda_dev))
            try
                gpu_processes = NVML.compute_processes(nvml_dev)
                if haskey(gpu_processes, getpid())
                    gpu_processes[getpid()].used_gpu_memory
                else
                    # happens when we didn't do compute, or when using containers:
                    # https://github.com/NVIDIA/gpu-monitoring-tools/issues/63
                    missing
                end
            catch err
                (isa(err, NVML.NVMLError) && err.code == NVML.ERROR_NOT_SUPPORTED) || rethrow()
                missing
            end
        else
            missing
        end
        passes,fails,error,broken,c_passes,c_fails,c_errors,c_broken =
            Test.get_test_counts(data[1])
        if data[1].anynonpass == false
            data = ((passes+c_passes,broken+c_broken),
                    data[2],
                    data[3],
                    data[4],
                    data[5],
                    data[6],
                    data[7],
                    data[8])
        end
        res = vcat(collect(data), cpu_rss, gpu_rss)

        device_reset!()
        res
    finally
        if snoop !== nothing
            ccall(:jl_dump_compiles, Nothing, (Ptr{Nothing},), C_NULL)
            close(io)
        end

        Test.TESTSET_PRINT_ENABLE[] = old_print_setting
    end
end

## auxiliary stuff

# NOTE: based on test/pkg.jl::capture_stdout, but doesn't discard exceptions
macro grab_output(ex)
    quote
        mktemp() do fname, fout
            ret = nothing
            open(fname, "w") do fout
                redirect_stdout(fout) do
                    ret = $(esc(ex))
                end
            end
            ret, read(fname, String)
        end
    end
end

# variant on @test_throws that checks the CuError error code
macro test_throws_cuerror(code, ex)
    # generate a test only returning CuError if it is the correct one
    test = quote
        try
            $(esc(ex))
        catch err
            isa(err, CuError) || rethrow()
            err.code == $code || error("Wrong CuError: ", err.code, " instead of ", $code)
            rethrow()
        end
    end

    # now re-use @test_throws (which ties into @testset, etc)
    quote
        @test_throws CuError $test
    end
end

# @test_throw, with additional testing for the exception message
macro test_throws_message(f, typ, ex...)
    quote
        msg = ""
        @test_throws $(esc(typ)) try
            $(esc(ex...))
        catch err
            msg = sprint(showerror, err)
            rethrow()
        end

        if !$(esc(f))(msg)
            # @test should return its result, but doesn't
            @error "Failed to validate error message\n$msg"
        end
        @test $(esc(f))(msg)
    end
end

# @test_throw, peeking into the load error for testing macro errors
macro test_throws_macro(ty, ex)
    return quote
        Test.@test_throws $(esc(ty)) try
            $(esc(ex))
        catch err
            @test err isa LoadError
            @test err.file === $(string(__source__.file))
            @test err.line === $(__source__.line + 1)
            rethrow(err.error)
        end
    end
end

# Run some code on-device, returning captured standard output
macro on_device(ex)
    @gensym kernel
    esc(quote
        let
            function $kernel()
                $ex
                return
            end

            @cuda $kernel()
            synchronize()
        end
    end)
end

# helper function for sinking a value to prevent the callee from getting optimized away
@inline sink(i::Int32) =
    Base.llvmcall("""%slot = alloca i32
                     store volatile i32 %0, i32* %slot
                     %value = load volatile i32, i32* %slot
                     ret i32 %value""", Int32, Tuple{Int32}, i)
@inline sink(i::Int64) =
    Base.llvmcall("""%slot = alloca i64
                     store volatile i64 %0, i64* %slot
                     %value = load volatile i64, i64* %slot
                     ret i64 %value""", Int64, Tuple{Int64}, i)

function julia_script(code, args=``)
    # FIXME: this doesn't work when the compute mode is set to exclusive
    script = """using CUDA
                device!($(device()))

                $code"""
    cmd = Base.julia_cmd()
    if Base.JLOptions().project != C_NULL
        cmd = `$cmd --project=$(unsafe_string(Base.JLOptions().project))`
    end
    cmd = `$cmd --color=no --eval $script $args`

    out = Pipe()
    err = Pipe()
    proc = run(pipeline(cmd, stdout=out, stderr=err), wait=false)
    close(out.in)
    close(err.in)
    wait(proc)
    proc.exitcode, read(out, String), read(err, String)
end

# tests that are conditionall broken
macro test_broken_if(cond, ex...)
    quote
        if $(esc(cond))
            @test_broken $(map(esc, ex)...)
        else
            @test $(map(esc, ex)...)
        end
    end
end

# some tests are mysteriously broken with certain hardware/software.
# use a horrible macro to mark those tests as "potentially broken"
@eval Test begin
    export @test_maybe_broken

    macro test_maybe_broken(ex, kws...)
        test_expr!("@test_maybe_broken", ex, kws...)
        orig_ex = Expr(:inert, ex)
        result = get_test_result(ex, __source__)
        quote
            x = $result
            if x.value
                do_test(x, $orig_ex)
            else
                do_broken_test(x, $orig_ex)
            end
        end
    end
end

nothing # File is loaded via a remotecall to "include". Ensure it returns "nothing".


text:
using LowDimNearestNeighbors
using Base.Test

# Define multidimensional vector types for testing
typealias Vec2{T <: Unsigned} Tuple{T, T}
Base.call{T}(::Type{Vec2{T}}, x::T, y::T) = (unsigned(x), unsigned(y))
Base.rand{T}(::Type{Vec2{T}}) = Vec2(rand(T), rand(T))

#

lessmsb = LowDimNearestNeighbors.lessmsb
@test  lessmsb(1, 2) # msb(001) <  msb(010)
@test !lessmsb(2, 1) # msb(010) >  msb(001)

@test !lessmsb(2, 3) # msb(010) == msb(011)
@test !lessmsb(3, 2) # msb(011) == msb(010)

@test  lessmsb(1, 5) # msb(001) <  msb(101)
@test !lessmsb(5, 1) # msb(101) >  msb(001)

#

shuffdim = LowDimNearestNeighbors.shuffdim
# The shuffdim should be 1 for identical points
@test shuffdim([1, 2, 3], [1, 2, 3]) == 1
@test shuffdim([123, 123], [123, 123]) == 1

# The first dimension is greater in p
@test shuffdim([1, 0], [0, 0]) == 1
@test shuffdim([1, 0], [0, 1]) == 1

# The first dimension is greater in q
@test shuffdim([0, 0], [1, 0]) == 1
@test shuffdim([0, 1], [1, 0]) == 1

# The second dimension is greater in p
@test shuffdim([1, 2], [1, 1]) == 2
# The second dimension is greater in q
@test shuffdim([1, 1], [1, 2]) == 2

#

@test  shuffless([1, 1, 3], [1, 2, 3])
@test !shuffless([1, 2, 3], [1, 2, 3])
@test !shuffless([1, 3, 3], [1, 2, 3])

@test !shuffmore([1, 1, 3], [1, 2, 3])
@test !shuffmore([1, 2, 3], [1, 2, 3])
@test  shuffmore([1, 3, 3], [1, 2, 3])

@test !shuffeq([1, 3, 3], [1, 2, 3])
@test shuffeq([1, 2, 3], [1, 2, 3])

# Test preprocess!
let
	arr = [Vec2(1, 1), Vec2(0, 1), Vec2(1, 0), Vec2(0, 0)]
	preprocess!(arr)
	@test arr[1] == Vec2(0, 0)
	@test arr[2] == Vec2(0, 1)
	@test arr[3] == Vec2(1, 0)
	@test arr[4] == Vec2(1, 1)

	# Error on negative coordinates
	@test_throws ErrorException preprocess!([1, 2, -3])

	# Error on noninteger coordinates
	@test_throws ErrorException preprocess!([1, 2.2, 3])
	@test_throws ErrorException preprocess!([1, 2.0, 3])
end

# Test Shifted indexing, clamping, and length
let
	el = Vec2{UInt64}(1, 2)
	shifted = LowDimNearestNeighbors.ShiftedPos{Vec2{UInt64}}(el, 5)
	@test shifted[1] == 6
	@test shifted[2] == 7
	@test length(shifted) == 2

	el = Vec2{UInt64}(4, 6)
	shifted = LowDimNearestNeighbors.ShiftedNeg{Vec2{UInt64}}(el, 5)
	@test shifted[1] == 0
	@test shifted[2] == 1
	@test length(shifted) == 2
end

# Test sqdist
let
	sqdist = LowDimNearestNeighbors.sqdist
	@test sqdist(2, 2) == 0
	@test sqdist(2, 3) == 1
	@test sqdist(2, 4) == 4
	@test sqdist(Vec2(3, 0), Vec2(0, 4)) == 5*5

	# Test saturation cases
	@test sqdist(0, typemax(UInt)) == typemax(UInt)
	@test sqdist(typemax(UInt), 0) == typemax(UInt)
	@test sqdist(Vec2(3, 3), Vec2(typemax(UInt), typemax(UInt))) == typemax(UInt)
end

# Test sqdist_to_quadtree_box
let
	sqdist_to_quadtree_box = LowDimNearestNeighbors.sqdist_to_quadtree_box
	@test sqdist_to_quadtree_box(Vec2(0, 0), Vec2(3, 0), Vec2(0, 3)) == 0
	@test sqdist_to_quadtree_box(Vec2(2, 2), Vec2(3, 0), Vec2(0, 3)) == 0

	@test sqdist_to_quadtree_box(Vec2(0, 0), Vec2(0, 0), Vec2(3, 3)) == 0
	@test sqdist_to_quadtree_box(Vec2(2, 2), Vec2(0, 0), Vec2(3, 3)) == 0

	@test sqdist_to_quadtree_box(Vec2(0, 0), Vec2(2, 2), Vec2(3, 3)) == 8
	@test sqdist_to_quadtree_box(Vec2(4, 4), Vec2(0, 0), Vec2(1, 1)) == 8
end

# Nearest: Deterministic test
let
	arr = preprocess!([Vec2(3, 3), Vec2(10, 0), Vec2(5, 8)])
	@test nearest(arr, Vec2(0, 0)) == Vec2(3, 3)
	@test nearest(arr, Vec2(8, 1)) == Vec2(10, 0)
	@test nearest(arr, Vec2(1000, 1000)) == Vec2(5, 8)
	@test nearest(arr, Vec2(5, 8)) == Vec2(5, 8)
end

# Nearest: Randomized test
let
	sqdist = LowDimNearestNeighbors.sqdist
	function linear_nearest(arr, q)
		local best
		best_sqdist = Inf
		for point in arr
			if sqdist(point, q) < best_sqdist
				best, best_sqdist = point, sqdist(point, q)
			end
		end
		best
	end

	function test(arr, queries; verbose=false)
		for q in queries
			result = nearest(arr, q)
			result_sqdist = sqdist(q, result)

			correct_result = linear_nearest(arr, q)
			correct_sqdist = sqdist(q, correct_result)

			if verbose && result_sqdist != correct_sqdist
				result_dist = sqrt(result_sqdist)
				correct_dist = sqrt(correct_sqdist)
				println("Mismatch when searching for ", q, ":")
				println("\t Result: ", result, "\t", result_dist)
				println("\tCorrect: ", correct_result, "\t", correct_dist)
				println("\t% error: ", 100 * (1 - correct_dist / result_dist), "%")
				println()
			end

			@test result_sqdist == correct_sqdist
		end
	end

	srand(0)
	numelements = numqueries = 1000

	test(
		preprocess!([rand(Vec2{UInt64}) for i in 1:numelements]),
		[rand(Vec2{UInt64}) for i in 1:numqueries]
	)

	test(
		preprocess!([rand(Vec2{UInt32}) for i in 1:numelements]),
		[rand(Vec2{UInt32}) for i in 1:numqueries]
	)

	test(
		preprocess!([rand(Vec2{UInt16}) for i in 1:numelements]),
		[rand(Vec2{UInt16}) for i in 1:numqueries]
	)

	test(
		preprocess!([rand(Vec2{UInt8}) for i in 1:numelements]),
		[rand(Vec2{UInt8}) for i in 1:numqueries]
	)

end

# Benchmarks
let
	function benchmark(numelements, numqueries)
		for i in 1:10
			arr = preprocess!([rand(Vec2{UInt8}) for i in 1:numelements])
			queries = [rand(Vec2{UInt8}) for i in 1:numqueries]
			@time for q in queries
				nearest(arr, q)
			end
		end
	end

	# benchmark(100000, 100000)
end


text:
abstract type Style end

mutable struct TextStyle <: Style
    font ::Option{String}
    hei  ::Option{F64}
    color::Option{Color}
end
TextStyle() = default(TextStyle)

mutable struct LineStyle <: Style
    lstyle::Option{Int}
    lwidth::Option{F64}
    smooth::Option{Bool}
    color ::Option{Color}
end
LineStyle() = default(LineStyle)

mutable struct MarkerStyle <: Style
    marker::Option{String}
    msize ::Option{F64}
    color ::Option{Color}
end
MarkerStyle() = default(MarkerStyle)

mutable struct BarStyle <: Style
    color::Option{Color}
    fill ::Colorant
#   pattern::Option{String}   =  .... see page 39 of manual, test first
end
BarStyle() = default(BarStyle)

mutable struct FillStyle <: Style
    fill::Colorant
end
FillStyle() = default(FillStyle)

mutable struct BoxplotStyle <: Style
    bwidth  ::F64          # width of the box
    wwidth  ::F64          # width of the whiskers
    wrlength::F64          # whisker length is wrlength * IQR, if INF will be min-max
    blstyle ::LineStyle    # box line style
    mlstyle ::LineStyle    # median line style
    mmstyle ::MarkerStyle  # mean marker style
    omstyle ::MarkerStyle  # outlier marker style
    mshow   ::Bool         # show the mean
    oshow   ::Bool         # show outliers
end
BoxplotStyle() = BoxplotStyle(
    GLE_DEFAULTS[:boxplotstyle_bwidth],
    GLE_DEFAULTS[:boxplotstyle_wwidth],
    GLE_DEFAULTS[:boxplotstyle_wrlength],
    LineStyle(
        GLE_DEFAULTS[:boxplotstyle_blstyle_lstyle],
        GLE_DEFAULTS[:boxplotstyle_blstyle_lwidth],
        false,
        GLE_DEFAULTS[:boxplotstyle_blstyle_color]
    ),
    LineStyle(
        GLE_DEFAULTS[:boxplotstyle_mlstyle_lstyle],
        GLE_DEFAULTS[:boxplotstyle_mlstyle_lwidth],
        false,
        GLE_DEFAULTS[:boxplotstyle_mlstyle_color]
    ),
    MarkerStyle(
        GLE_DEFAULTS[:boxplotstyle_mmstyle_marker],
        GLE_DEFAULTS[:boxplotstyle_mmstyle_msize],
        GLE_DEFAULTS[:boxplotstyle_mmstyle_color]
    ),
    MarkerStyle(
        GLE_DEFAULTS[:boxplotstyle_omstyle_marker],
        GLE_DEFAULTS[:boxplotstyle_omstyle_msize],
        GLE_DEFAULTS[:boxplotstyle_omstyle_color]
        ),
    GLE_DEFAULTS[:boxplotstyle_mshow],
    GLE_DEFAULTS[:boxplotstyle_oshow]
)


text:
import Base.isempty

export AbstractCentrallySymmetricPolytope,
       center,
       an_element,
       vertices_list

"""
    AbstractCentrallySymmetricPolytope{N} <: AbstractPolytope{N}

Abstract type for centrally symmetric, polytopic sets.
It combines the `AbstractCentrallySymmetric` and `AbstractPolytope` interfaces.
Such a type combination is necessary as long as Julia does not support
[multiple inheritance](https://github.com/JuliaLang/julia/issues/5).

### Notes

Every concrete `AbstractCentrallySymmetricPolytope` must define the following
functions:
- from `AbstractCentrallySymmetric`:
  - `center(::AbstractCentrallySymmetricPolytope)` -- return the center point
- from `AbstractPolytope`:
  - `vertices_list(::AbstractCentrallySymmetricPolytope)`
     -- return a list of all vertices

```jldoctest; setup = :(using LazySets: subtypes)
julia> subtypes(AbstractCentrallySymmetricPolytope)
2-element Array{Any,1}:
 AbstractZonotope
 Ball1
```
"""
abstract type AbstractCentrallySymmetricPolytope{N} <: AbstractPolytope{N} end

isconvextype(::Type{<:AbstractCentrallySymmetricPolytope}) = true

# --- common AbstractCentrallySymmetric functions (copy-pasted) ---

"""
    dim(P::AbstractCentrallySymmetricPolytope)

Return the ambient dimension of a centrally symmetric, polytopic set.

### Input

- `P` -- centrally symmetric, polytopic set

### Output

The ambient dimension of the polytopic set.
"""
@inline function dim(P::AbstractCentrallySymmetricPolytope)
    return length(center(P))
end

"""
    an_element(P::AbstractCentrallySymmetricPolytope)

Return some element of a centrally symmetric polytope.

### Input

- `P` -- centrally symmetric polytope

### Output

The center of the centrally symmetric polytope.
"""
function an_element(P::AbstractCentrallySymmetricPolytope)
    return center(P)
end

"""
    isempty(P::AbstractCentrallySymmetricPolytope)

Return if a centrally symmetric, polytopic set is empty or not.

### Input

- `P` -- centrally symmetric, polytopic set

### Output

`false`.
"""
function isempty(::AbstractCentrallySymmetricPolytope)
    return false
end

"""
    isuniversal(S::AbstractCentrallySymmetricPolytope{N},
                [witness]::Bool=false) where {N}

Check whether a centrally symmetric polytope is universal.

### Input

- `S`       -- centrally symmetric polytope
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `false`
* If `witness` option is activated: `(false, v)` where ``v ∉ S``

### Algorithm

A witness is obtained by computing the support vector in direction
`d = [1, 0, …, 0]` and adding `d` on top.
"""
function isuniversal(S::AbstractCentrallySymmetricPolytope{N},
                     witness::Bool=false) where {N}
    if witness
        d = SingleEntryVector{N}(1, dim(S))
        w = σ(d, S) + d
        return (false, w)
    else
        return false
    end
end

"""
    center(S::AbstractCentrallySymmetricPolytope, i::Int)

Return the center along a given dimension of a centrally symmetric polytope.

### Input

- `S` -- centrally symmetric polytope
- `i` -- dimension of interest

### Output

The center along a given dimension of the centrally symmetric polytope.
"""
@inline function center(S::AbstractCentrallySymmetricPolytope, i::Int)
    return center(S)[i]
end


text:
using ModelingToolkit
using DifferentialEquations
using Distributions
using Tables
using DataFrames
using StatsPlots
using BenchmarkTools

@inline function rate_to_proportion(r,t)
    1-exp(-r*t)
end;

@parameters β=0.05 c=10.0 γ=0.25 N=1000.0 δt=0.1

@variables t SI(t)=0.0 IR(t)=0.0 S(t)=990.0 I(t)=10.0 R(t)=0.0

D = DiscreteUpdate(t; dt=δt)
eqs = [D(SI) ~ rate_to_proportion(β*c*I/N,δt)*S,
       D(IR) ~ rate_to_proportion(γ,δt)*I,
       D(S) ~ S-SI,
       D(I) ~ I+SI-IR,
       D(R) ~ R+IR]

@named sys = DiscreteSystem(eqs)

@named sys = DiscreteSystem(eqs, t, [SI, IR, S, I, R], [β,c,γ,N,δt])

prob = DiscreteProblem(sys, [], (0,400), [])

u0 = [S => 990.0,
      I => 10.0,
      R => 0.0];
p = [β=>0.05,
     c=>10.0,
     γ=>0.25,
     N=>1000.0,
     δt=>0.1];
prob = DiscreteProblem(sys, u0, (0,400), p)

sol = solve(prob,solver=FunctionMap);

df = DataFrame(Tables.table(sol'))
rename!(df,["SI", "IR", "S", "I", "R"])
df[!,:t] = 0:0.1:40.0;

@df df plot(:t,
    [:S :I :R],
    xlabel="Time",
    ylabel="Number")

C = cumsum(df[!,:SI])
cases = vcat(C[1:9],C[10:end] .- C[1:(end-9)])
df[!,"cases"] = cases

@df df plot(:t,
    [:cases],
    xlabel="Time",
    ylabel="Cases per day")

@benchmark solve(prob,solver=FunctionMap)

include(joinpath(@__DIR__,"tutorials","appendix.jl"))
appendix()


text:
include("ModuleA.jl")

using MyModule: foo, MyType
#euivalent to
#using MyModule.foo, MyModule.MyType

mytype = MyType(3)
@show foo(mytype)
@show MyType
#@show bar(2)

@show MyModule.foo(mytype)
@show MyModule.MyType
@show MyModule.bar(2)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
mutable struct TimeGradientWrapper{fType,uType,P} <: Function
  f::fType
  uprev::uType
  p::P
end
(ff::TimeGradientWrapper)(t) = (du2 = similar(ff.uprev); ff.f(du2,ff.uprev,ff.p,t); du2)
(ff::TimeGradientWrapper)(du2,t) = ff.f(du2,ff.uprev,ff.p,t)

mutable struct UJacobianWrapper{fType,tType,P} <: Function
  f::fType
  t::tType
  p::P
end

(ff::UJacobianWrapper)(du1,uprev) = ff.f(du1,uprev,ff.p,ff.t)
(ff::UJacobianWrapper)(uprev) = (du1 = similar(uprev); ff.f(du1,uprev,ff.p,ff.t); du1)

mutable struct TimeDerivativeWrapper{F,uType,P} <: Function
  f::F
  u::uType
  p::P
end
(ff::TimeDerivativeWrapper)(t) = ff.f(ff.u,ff.p,t)

mutable struct UDerivativeWrapper{F,tType,P} <: Function
  f::F
  t::tType
  p::P
end
(ff::UDerivativeWrapper)(u) = ff.f(u,ff.p,ff.t)

mutable struct ParamJacobianWrapper{fType,tType,uType} <: Function
  f::fType
  t::tType
  u::uType
end

function (ff::ParamJacobianWrapper)(du1,p)
  ff.f(du1,ff.u,p,ff.t)
end

function (ff::ParamJacobianWrapper)(p)
  du1 = similar(p, size(ff.u))
  ff.f(du1,ff.u,p,ff.t)
  return du1
end

mutable struct JacobianWrapper{fType,pType}
  f::fType
  p::pType
end

(uf::JacobianWrapper)(u) = uf.f(u, uf.p)
(uf::JacobianWrapper)(res, u) = uf.f(res, u, uf.p)


text:
macro gcinfo(x); end
# Uncomment to print gc information:
#macro gcinfo(x); esc(:(print(($x)))); end

# KnetPtr type holds a gpu (dev>=0) or cpu (dev=-1) allocated pointer.
# We try to minimize the number of actual allocations, which are slow,
# by reusing preallocated but garbage collected pointers.

type KnetPtr
    ptr::Cptr                   # actual pointer
    len::Int                    # size in bytes
    dev::Int                    # id of the device the pointer belongs to
    parent                      # used to implement shared memory pointers
end

# We use the KnetPtrs type to keep track of allocated and garbage
# collected pointers: We keep one KnetPtrs struct per size per device.

type KnetPtrs
    used::Int                   # number of allocated pointers
    free::Array{Cptr,1}         # pointers available for reuse
    KnetPtrs()=new(0,Array(Cptr,0))
end

# KnetFree[dev+2] will hold a dictionary from sizes to KnetPtrs for
# device dev.  KnetFree[1] reserved for the cpu in case we decide to
# support it in the future.  It is initialized in KnetPtr(nbytes).

KnetFree = nothing
initKnetFree()=(global KnetFree=[ Dict{Int,KnetPtrs}() for i=1:gpuCount()+1 ])

# When Julia gc reclaims a KnetPtr object, the following special
# finalizer does not actually release the memory, but inserts it in
# the KnetFree[dev+2] dict keyed by length in bytes so it can be
# reused.

function freeKnetPtr(p::KnetPtr)
    ptrs = KnetFree[p.dev+2][p.len]
    push!(ptrs.free,p.ptr)
end

# This is the low level KnetPtr constructor, it adds the finalizer and
# sets parent to `nothing` which is only needed for shared pointers.

function KnetPtr(ptr::Cptr,len::Int,dev::Int)
    kp = KnetPtr(ptr,len,dev,nothing)
    finalizer(kp, freeKnetPtr)
    return kp
end

# This constructor is used to create a shared pointer.  We need to
# keep the parent field to prevent premature gc of the parent.  The
# child does not need a special finalizer.

function KnetPtr(parent::KnetPtr, offs::Integer, len::Integer)
    if len < 0 || offs < 1 || offs+len-1 > parent.len
        throw(BoundsError())
    end
    KnetPtr(parent.ptr+offs-1, len, parent.dev, parent)
end

# This the main KnetPtr constructor.  It tries to avoid actual
# allocation which is slow.  First it tries to find a previously
# allocated and garbage collected pointer in KnetFree[dev+2].  If not
# available it tries to allocate a new one (about 10 μs).  Otherwise
# it tries running gc() and see if we get a pointer back (about 75
# ms).  Finally if all else fails, it calls knetgc which cleans up all
# allocated and garbage collected KnetPtrs on the current device and
# tries allocation one last time.

function KnetPtr(nbytes::Integer)
    KnetFree==nothing && initKnetFree()
    dev = gpu()
    if dev < 0
        error("KnetPtr: bad device id $dev.")
    end
    ptrs = get!(KnetPtrs,KnetFree[dev+2],nbytes)
    if !isempty(ptrs.free)
        return KnetPtr(pop!(ptrs.free),nbytes,dev)
    end
    ptr = knetMalloc(nbytes)
    if ptr != nothing
        ptrs.used += 1
        return KnetPtr(ptr,nbytes,dev)
    end
    gc(); @gcinfo(".")
    if !isempty(ptrs.free)
        return KnetPtr(pop!(ptrs.free),nbytes,dev)
    end
    knetgc(); @gcinfo("+")
    ptr = knetMalloc(nbytes)
    if ptr != nothing
        ptrs.used += 1
        return KnetPtr(ptr,nbytes,dev)
    end
    error("Out of gpu memory")
end

# This does the actual allocation, returns `nothing` in case of error
function knetMalloc(nbytes::Int)
    # we no longer support cpu pointers, all overloaded ops rely on KnetPtr being on a GPU
    # gpu() >= 0 || return(convert(Cptr, pointer(Array(UInt8,nbytes))))
    ptr = Cptr[0]
    ret = @cuda1(cudart,cudaMalloc,(Ptr{Cptr},Csize_t),ptr,nbytes)
    if ret == 0
        return ptr[1]
    else # error("cudaMalloc($nbytes) error $ret")
        return nothing
    end
end

# If you really want to clean up memory you need to call knetgc: Note
# that this only cleans the current device.  It frees all pointers
# previously allocated and garbage collected.
function knetgc(dev=gpu())
    gc_enable(false)
    for v in values(KnetFree[dev+2])
        if dev >= 0
            for p in v.free
                @cuda(cudart,cudaFree,(Cptr,),p)
            end
        end
        v.used -= length(v.free)
        empty!(v.free)
    end
    gc_enable(true)
end

# Some utilities
meminfo(i=gpu())=[(k,v.used,length(v.free)) for (k,v) in KnetFree[i+2]]
gpufree(i=gpu())=nvmlDeviceGetMemoryInfo(i)[2]

function gpuinfo(msg="",dev=gpu())
    msg != "" && print("$msg ")
    g = nvmlDeviceGetMemoryInfo(dev)
    println((:dev,dev,:total,g[1],:free,g[2],:used,g[3]))
    for (s,u,f) in sort(meminfo(dev), by=(x->x[1]), rev=true)
        println((:size,s,:alloc,u,:avail,f))
    end
end

function memdbg(msg="")
    m = nvmlDeviceGetMemoryInfo()
    c = cudaGetMemInfo()
    x = [0,0]
    for (s,u,f) in meminfo()
        x[1] += s*u
        x[2] += s*f
    end
    println("$msg: ntotal: $(m[1]) nfree: $(m[2]) nused: $(m[3]) ktotal: $(x[1]) kfree: $(x[2])\n nfree+used: $(m[2]+m[3]) ctotal: $(c[2]) cfree: $(c[1]) ktotal-kfree: $(x[1]-x[2]) nused-ktotal: $(m[3]-x[1])")
end


text:
"""
Gridap, grid-based approximation of PDEs in the Julia programming language

This module provides rich set of tools for the numerical solution of PDE, mainly based
on finite element methods.

The module is structured in the following sub-modules:

- [`InplaceArrays.Helpers`](@ref)
- [`InplaceArrays.Inference`](@ref)
- [`InplaceArrays.TensorValues`](@ref)
- [`InplaceArrays.Arrays`](@ref)
- [`InplaceArrays.Fields`](@ref)
- [`InplaceArrays.Polynomials`](@ref)
- [`InplaceArrays.ReferenceFEs`](@ref)

The exported names are:
$(EXPORTS)
"""
module InplaceArrays

using DocStringExtensions

include("Helpers/Helpers.jl")

include("Inference/Inference.jl")

include("TensorValues/TensorValues.jl")

include("Arrays/Arrays.jl")

include("Fields/Fields.jl")

include("Polynomials/Polynomials.jl")

include("ReferenceFEs/ReferenceFEs.jl")

end # module


text:
# @testset "Indices" begin
#     s = ContactImplicitMPC.get_simulation("hopper_2D", "flat_2D_lc", "flat")
#
#     model = s.model
#     env = s.env
#     for quat ∈ (false, true)
#         terms, vars = ContactImplicitMPC.get_bilinear_indices(model, env, quat = quat)
#         idyn, irst, ibil, ialt = ContactImplicitMPC.linearization_term_index(model, env, quat = quat)
#         ix, iy1, iy2 = ContactImplicitMPC.linearization_var_index(model, env, quat = quat)
#         nquat = 0
#         @test iy1 == vcat(Vector.([vars[1][1], vars[2][1], vars[3][1]])...)
#         @test iy2 == vcat(Vector.([vars[1][2], vars[2][2], vars[3][2]])...)
#         @test ix == setdiff(1:ContactImplicitMPC.num_var(model, env) - nquat, vcat(iy1, iy2))
#         @test ibil == vcat(Vector.([terms[1], terms[2], terms[3]])...)
#     end
# end


text:
module FluPredictibility

include("BioTools/src/BioTools.jl")
include("TreeTools/src/TreeTools.jl")
include("EarthMoversDistance/src/EarthMoversDistance.jl")

using FluPredictibility.BioTools
using FluPredictibility.TreeTools
using BioSequences
using Dates, DelimitedFiles, StatsBase, JSON, DataStructures, DataFrames
using Polynomials

import StatsBase.entropy
import Base.length

export FluPop, AAFluPop
export PosEvo
export bin_by_date!

# Should I implement a flu strain that has segment / flustrain fields and a ntseq field?
# For now, I'll do this for AAs only.
# One way might also be to do codon alphabets in BioSequence
mutable struct FluPop{T<:AbstractStrain}
	strains::Dict{String, T}
	datebin::Dict{Tuple{Date,Date},Array{T,1}}
end
function FluPop(; strains = Dict{String,Strain}(),
				datebin = Dict{Tuple{Date,Date},Array{eltype(values(strains)),1}}())
	return FluPop(strains, datebin)
end
function FluPop(S::Array{<:AbstractStrain,1})
	return FluPop(strains = Dict(String(x[:strain])=>x for x in S))
end
"""
	add_strain_field!(fp::FluPop, field, default=0.)
"""
function add_strain_field!(fp::FluPop, field, default=0.)
	for s in values(fp.strains)
		s.data[field] = default
	end
	nothing
end
"""
	consensus(fp::FluPop)
"""
consensus(fp::FluPop) = consensus([x for x in values(fp.strains)])

mutable struct PosEvo{A}
	i::Int64
	alphabet::Array{A,1}
	data::Dict{Tuple{Date,Date}, SiteFrequency{A}}
end
function PosEvo(A::DataType, i::Int64)
	return PosEvo(i, Array{A,1}(undef,0), Dict{Tuple{Date,Date}, SiteFrequency{A}}())
end

"""
	mutable struct FrequencyTraj{A}
"""
mutable struct FrequencyTraj{A}
	# All fields except `.data ` are initialized when calling the `all_trajectories` function.
	# Info on what trajectory represents
	i::Int64
	val::A # Symbol that trajectory represents
	# Data
	date::Date # Date at t=0
	t::Array{Day,1} # Days
	freq::Array{Float64,1}
	pop::Array{Int64,1} # Size of (sampled) population at that time.
	strains::Array{<:Array{<:AbstractString, 1}, 1}
	# Fixation
	index::Dict{Symbol, Union{Int64, Missing}} # Storing indices of start / end / activity (ie t=0)
	fixation::Symbol # :fixed, :lost, :poly
	# Misc. data
	data::Dict # Expected fields --
end
length(t::FrequencyTraj) = length(t.t)
# getproperty(t, field) = hasfield(t, field) ? getfield(t, field) : t.data[field]

include("IO.jl")
include("global.jl")
include("filtering.jl")
include("frequencies.jl")
include("trajectories.jl")
include("lbi.jl")
include("fitness.jl")
include("plot_functions.jl")
include("prediction_functions.jl")
include("misc.jl")

end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility