text:
# Script to:
# - load diff bias csv files,
# - average to find most biased document sets,
# - create perturbation cooc files for those document sets

using DataFrames
using CSV
using Statistics
using SparseArrays
using LinearAlgebra
using Random

include("Corpora.jl")
include("GloVe.jl")
include("Bias.jl")
include("utils.jl")

# Command ARGS
const target = get(ARGS, 1, "C0-V15-W8-D75-R0.05-E300")
const diff_bias_dir = get(ARGS, 2, "results/diff_bias")
const embedding_dir = get(ARGS, 3, "embeddings")
const corpus_dir = get(ARGS, 4, "corpora")
const pert_dir = get(ARGS, 5, "results/perturbations")

const corpus_prefix = match(r"C[0-9]+-", target).match == "C0-" ? "simplewiki" : "nyt"
const corpus_path = joinpath(corpus_dir, "$(corpus_prefix)select.txt")
const corpus = Corpora.Corpus(corpus_path)

const vocab_prefix = match(r"C[0-9]+-V[0-9]+", target).match
const vocab_path = joinpath(embedding_dir, "vocab-$vocab_prefix.txt")
const (vocab, ivocab) = GloVe.load_vocab(vocab_path)

const window = parse(Int, match(r"W[0-9]+", target).match[2:end])

const set_sizes = corpus_prefix == "nyt" ? [100, 300, 1000, 3000, 10000] : [10, 30, 100, 300, 1000]

# Load a bunch of diff bias results and average them
function load_df(target, bias_col=:ΔBIF_1)
    df = DataFrame()
    i = 0
    cols = []
    for f in readdir(diff_bias_dir)
        if occursin(Regex(target), f)
            i += 1
            tmp = sort(CSV.read(joinpath(diff_bias_dir, f)), :doc_num)
            if (i == 1)
                df[:doc_num] = tmp[:doc_num]
            end
            col = Symbol("ΔBIF_$(i)")
            df[col] = tmp[bias_col]
            push!(cols, col)
        end
    end
    println("Loaded $i differential bias results")
    df[:ΔBIF_μ] = [mean([row[col] for col in cols]) for row in eachrow(df)]
    df[:ΔBIF_σ] = [std([row[col] for col in cols]) for row in eachrow(df)]
    return df
end

function make_perturbation(name, article_ids, biases, out_dir; verbose=true)
    verbose && println("Making $(name)... ")
    articles = Corpora.get_texts(corpus, article_ids)
    δX = GloVe.docs2cooc(articles, vocab, window)
    bias = sum(biases)
    data = (name, nnz(δX), vecnorm(δX), bias)
    GloVe.save_coocs(joinpath(out_dir, "pert-$name.bin"), δX)
    verbose && println("Done.\n")
end

function make_perturbations(target, bias_col=:ΔBIF_1, dir_name="B1")
    out_dir = joinpath(pert_dir, dir_name)
    mkpath(out_dir)
    df = load_df(target, bias_col)
    sorted = sort(df[[:ΔBIF_μ, :doc_num]])

    for N in set_sizes
        # Correcting Articles
        name = "correct_$(N)_1"
        article_ids, biases = sorted[:doc_num][1:N], sorted[:ΔBIF_μ][1:N]
        make_perturbation(name, article_ids, biases, out_dir)

        # Aggravating Articles
        name = "aggravate_$(N)_1"

        article_ids = sorted[:doc_num][end-(N-1):end]
        biases = sorted[:ΔBIF_μ][end-(N-1):end]
        make_perturbation(name, article_ids, biases, out_dir)

        # Random
        all_article_ids = Array{Int64}(df[:doc_num])
        for run in 1:3
            name = "random_$(N)_$(run)"
            article_ids = shuffle(all_article_ids)[1:N]
            biases = [df[:ΔBIF_μ][df[:doc_num] .== id][1] for id in article_ids]
            make_perturbation(name, article_ids, biases, out_dir)
        end
    end
    return out_dir
end

function make_config_file(target, out_dir)
    params = fileinfo("-" * target * "-")
    open(joinpath(out_dir, "config.txt"), "w") do f
        println(f, "CORPUS_ID=$(params.corpus)")
        println(f, "VOCAB_MIN_COUNT=$(params.min_vocab)")
        println(f, "WINDOW_SIZE=$(params.window)")
        println(f, "VECTOR_SIZE=$(params.dimension)")
        println(f, "ETA=$(params.eta)")
        println(f, "MAX_ITER=$(params.max_iters)")
    end
end

function main()
    for i in 1:length(Bias.WEAT_WORD_SETS)
        out_dir = make_perturbations(target, Symbol("ΔBIF_$i"), "$target-B$i")
        make_config_file(target, out_dir)
    end
end

main()


text:
module VariableTypes

export  vartype,
        VariableType,
            Categorical,
            Quantitative,
                Continuous,
                Discrete

include("vartypes.jl")

end # module

ext: jl


text:
# Implement a hFILE backend using arbitrary julia IO type

struct hFILE_julia
    base::hFILE_base
    ioptr::Ptr{Cvoid}
end

function _restore_hfile_io(fp)
    fp = convert(Ptr{hFILE_julia}, fp)
    unsafe_pointer_to_objref(unsafe_load(fp).ioptr)
end

function _debug_io(msg)
    # @async println(msg)
end

# As per read(2), returning the number of bytes read (possibly 0) or
# negative (and setting errno) on errors.  Front-end code will call this
# repeatedly if necessary to attempt to get the desired byte count.
function hfile_backend_read(fp, buffer, nbytes)
    _hfile_backend_read(_restore_hfile_io(fp), buffer, nbytes)
end
function _hfile_backend_read(io::IO, buffer::Ptr{Cvoid}, nbytes::Csize_t)::Cssize_t
    ## TODO: handle error, set errno and return negative value
    # TranscodingStreams.unsafe_read is different from Base.unsafe_read
    bytes_written = TranscodingStreams.unsafe_read(io, convert(Ptr{UInt8}, buffer), convert(Int, nbytes))
    _debug_io("hfile: read $bytes_written ($nbytes requested)")
    return bytes_written
end
# As per write(2), returning the number of bytes written or negative (and
# setting errno) on errors.  Front-end code will call this repeatedly if
# necessary until the desired block is written or an error occurs.
function hfile_backend_write(fp, buffer, nbytes)
    _hfile_backend_write(_restore_hfile_io(fp), buffer, nbytes)
end
function _hfile_backend_write(io::IO, buffer::Ptr{Cvoid}, nbytes::Csize_t)::Cssize_t
    ## TODO: handle error, set errno and return negative value
    bytes_written = unsafe_write(io, buffer, nbytes)
    _debug_io("hfile: write $bytes_written ($nbytes requested)")
    return bytes_written
end

# As per lseek(2), returning the resulting offset within the stream or
# negative (and setting errno) on errors.
function hfile_backend_seek(fp, offset, whence)
    ans = _hfile_backend_seek(_restore_hfile_io(fp), offset, whence)
    _debug_io("hfile: seek offset $offset with result $ans (whence $whence)")
    return ans
end
function _hfile_backend_seek(io::IO, offset::Coff_t, whence::Cint)::Coff_t
    SEEK_SET = Cint(0)
    SEEK_CUR = Cint(1)
    SEEK_END = Cint(2)
    origin::Coff_t = typemin(Coff_t)
    if whence == SEEK_SET
        origin = 0
    elseif whence == SEEK_CUR
        if !Core.applicable(position, io)
            Base.Libc.errno(Base.Libc.ESPIPE)
            return -1
        end
        origin = position(io)
    elseif whence == SEEK_END
        if !Core.applicable(seekend, io)
            Base.Libc.errno(Base.Libc.ESPIPE)
            return -1
        end
        # trick to get end position
        seekend(io)
        origin = position(io)
    else
        println("invalid whence argument: $whence")
        Base.Libc.errno(Base.Libc.EINVAL)
        return -1
    end
    realoffset = origin + offset

    if !Core.applicable(seek, io, realoffset)
        Base.Libc.errno(Base.Libc.ESPIPE)
        return -1
    end

    seek(io, realoffset)
    return realoffset
end

# Performs low-level flushing, if any, e.g., fsync(2); for writing streams
# only.  Returns 0 for success or negative (and sets errno) on errors.
function hfile_backend_flush(fp)
    _debug_io("hfile: flush")
    _hfile_backend_flush(_restore_hfile_io(fp))
end
function _hfile_backend_flush(io::IO)::Cint
    flush(io)
    return 0
end

# Closes the underlying stream (for output streams, the buffer will
# already have been flushed), returning 0 for success or negative (and
# setting errno) on errors, as per close(2).
function hfile_backend_close(fp)
    _debug_io("hfile: close")
    _hfile_backend_close(_restore_hfile_io(fp))
end
function _hfile_backend_close(io::IO)::Cint
    # if the io object didn't implement close, we shouldn't worry about it
    if !Core.applicable(close, io)
        return 0
    end
    close(io)
    return 0
end

function hfile_julia_backend()
    read_f  = @cfunction(hfile_backend_read,  Cssize_t, (Ptr{hFILE}, Ptr{Cvoid}, Csize_t))
    write_f = @cfunction(hfile_backend_write, Cssize_t, (Ptr{hFILE}, Ptr{Cvoid}, Csize_t))
    seek_f  = @cfunction(hfile_backend_seek,  Coff_t,   (Ptr{hFILE}, Coff_t, Cint))
    flush_f = @cfunction(hfile_backend_flush, Cint,     (Ptr{hFILE}, ))
    close_f = @cfunction(hfile_backend_close, Cint,     (Ptr{hFILE}, ))
    hFILE_backend(read_f, write_f, seek_f, flush_f, close_f)
end


text:
using CGDycore
using LinearAlgebra
n=5
u=zeros(n)
d=zeros(n)
l=zeros(n)
x=zeros(n)
y=zeros(n)
for i=1:n
  u[i] = 2*i
  l[i] = i
  d[i] = 3*i
  x[i] = i*i
end
u .= u .- 0
l .= l .+ 0
d .= d .+ 10

ev = l[1:n-1]
dv = d
biLL = Bidiagonal(dv, ev, :L)

ev = u[2:n]
dv = d
biUU = Bidiagonal(dv, ev, :U)

triSS = biUU * biLL

tri = zeros(3,n)
biU = zeros(2,n)
biL = zeros(2,n)
biL[1,:] = d
biL[2,1:n-1] = l[1:n-1]
biU[1,2:n] = u[2:n]
biU[2,:] = d

CGDycore.mulUL!(tri,biU,biL)

y .= triSS*x
@show tri
@show y
@show x

CGDycore.triSolve!(x,tri,y)

@show x


text:
using KrylovKit
using Rimu
using Test

@testset "Krylov eigsolve with BoseFS{6,6}" begin
    ham = HubbardReal1D(near_uniform(BoseFS{6,6}); u=6.0, t=1.0)

    a_init = near_uniform(ham)
    c_init = DVec(a_init => 1.0)

    all_results = eigsolve(ham, c_init, 1, :SR; issymmetric = true)
    energy = all_results[1][1]

    @test energy ≈ -4.0215 atol=0.0001
end


text:
module Turing

##############
# Dependency #
########################################################################
# NOTE: when using anything from external packages,                    #
#       let's keep the practice of explictly writing Package.something #
#       to indicate that's not implemented inside Turing.jl            #
########################################################################

using Requires
using Reexport
@reexport using Distributions
@reexport using MCMCChain
using ForwardDiff
using StatsFuns

using LinearAlgebra
using ProgressMeter
using Markdown
using Libtask
using MacroTools

@init @require Stan="682df890-35be-576f-97d0-3d8c8b33a550" @eval begin
    using Stan
    import Stan: Adapt, Hmc
end
import Base: ~, convert, promote_rule, rand, getindex, setindex!
import Distributions: sample
import ForwardDiff: gradient
using Flux: Tracker
import MCMCChain: AbstractChains, Chains
@init @require DynamicHMC="bbc10e6e-7c05-544b-b16e-64fede858acb" @eval begin
    using DynamicHMC, LogDensityProblems
    using LogDensityProblems: AbstractLogDensityProblem, ValueGradient

    struct FunctionLogDensity{F} <: AbstractLogDensityProblem
      dimension::Int
      f::F
    end

    LogDensityProblems.dimension(ℓ::FunctionLogDensity) = ℓ.dimension

    LogDensityProblems.logdensity(::Type{ValueGradient}, ℓ::FunctionLogDensity, x) = ℓ.f(x)::ValueGradient
end

##############################
# Global variables/constants #
##############################

const ADBACKEND = Ref(:reverse_diff)
setadbackend(backend_sym) = begin
  @assert backend_sym == :forward_diff || backend_sym == :reverse_diff
  backend_sym == :forward_diff && CHUNKSIZE[] == 0 && setchunksize(40)
  ADBACKEND[] = backend_sym
end

const ADSAFE = Ref(false)
setadsafe(switch::Bool) = begin
  @info("[Turing]: global ADSAFE is set as $switch")
  ADSAFE[] = switch
end

const CHUNKSIZE = Ref(40) # default chunksize used by AD

setchunksize(chunk_size::Int) = begin
  if ~(CHUNKSIZE[] == chunk_size)
    @info("[Turing]: AD chunk size is set as $chunk_size")
    CHUNKSIZE[] = chunk_size
  end
end

const PROGRESS = Ref(true)
turnprogress(switch::Bool) = begin
  @info("[Turing]: global PROGRESS is set as $switch")
  PROGRESS[] = switch
end

# Constants for caching
const CACHERESET  = 0b00
const CACHEIDCS   = 0b10
const CACHERANGES = 0b01

#######################
# Sampler abstraction #
#######################

abstract type InferenceAlgorithm end
abstract type Hamiltonian <: InferenceAlgorithm end
abstract type AbstractSampler end
"""
    Sampler{T}

Generic interface for implementing inference algorithms.
An implementation of an algorithm should include the following:

1. A type specifying the algorithm and its parameters, derived from InferenceAlgorithm
2. A method of `sample` function that produces results of inference, which is where actual inference happens.

Turing translates models to chunks that call the modelling functions at specified points. The dispatch is based on the value of a `sampler` variable. To include a new inference algorithm implements the requirements mentioned above in a separate file,
then include that file at the end of this one.
"""
mutable struct Sampler{T<:InferenceAlgorithm} <: AbstractSampler
  alg   ::  T
  info  ::  Dict{Symbol, Any}         # sampler infomation
end

"""
Robust initialization method for model parameters in Hamiltonian samplers.
"""
struct HamiltonianRobustInit <: AbstractSampler end
struct SampleFromPrior <: AbstractSampler end

# This can be removed when all `spl=nothing` is replaced with
#   `spl=SampleFromPrior`
const AnySampler = Union{Nothing, AbstractSampler}

include("utilities/resample.jl")
@init @require Stan="682df890-35be-576f-97d0-3d8c8b33a550" @eval begin
    include("support/stan-interface.jl")
end
include("utilities/helper.jl")
include("utilities/transform.jl")
include("utilities/robustinit.jl")
include("utilities/util.jl")         # utility functions
include("utilities/io.jl")           # I/O
include("models/distributions.jl")
include("core/varinfo.jl")  # core internal variable container
include("core/trace.jl")   # to run probabilistic programs as tasks

using Turing.VarReplay

###########
# Exports #
###########

# Turing essentials - modelling macros and inference algorithms
export @model, @~, @VarName                   # modelling
export MH, Gibbs                              # classic sampling
export HMC, SGLD, SGHMC, HMCDA, NUTS          # Hamiltonian-like sampling
export DynamicNUTS
export IS, SMC, CSMC, PG, PIMH, PMMH, IPMCMC  # particle-based sampling
export sample, setchunksize, resume           # inference
export auto_tune_chunk_size!, setadbackend, setadsafe # helper
export turnprogress  # debugging
export consume, produce

# Turing-safe data structures and associated functions
export TArray, tzeros, localcopy, IArray

export @sym_str

export Flat, FlatPos

##################
# Inference code #
##################

include("core/compiler.jl")     # compiler
include("core/container.jl")    # particle container
include("samplers/sampler.jl")  # samplers
include("core/ad.jl")           # Automatic Differentiation

end


text:
export A2CLearner

using Flux

"""
    A2CLearner(;kwargs...)

# Keyword arguments

- `approximator`::[`ActorCritic`](@ref)
- `γ::Float32`, reward discount rate.
- `actor_loss_weight::Float32`
- `critic_loss_weight::Float32`
- `entropy_loss_weight::Float32`
"""
Base.@kwdef mutable struct A2CLearner{A<:ActorCritic} <: AbstractLearner
    approximator::A
    γ::Float32
    max_grad_norm::Union{Nothing,Float32} = nothing
    norm::Float32 = 0.0f0
    actor_loss_weight::Float32
    critic_loss_weight::Float32
    entropy_loss_weight::Float32
    actor_loss::Float32 = 0.0f0
    critic_loss::Float32 = 0.0f0
    entropy_loss::Float32 = 0.0f0
    loss::Float32 = 0.0f0
end

function (learner::A2CLearner)(env::MultiThreadEnv)
    learner.approximator.actor(send_to_device(
        device(learner.approximator),
        get_state(env),
    )) |> send_to_host
end

function (learner::A2CLearner)(env)
    s = get_state(env)
    s = Flux.unsqueeze(s, ndims(s) + 1)
    s = send_to_device(device(learner.approximator), s)
    learner.approximator.actor(s) |> vec |> send_to_host
end

function RLBase.update!(learner::A2CLearner, t::AbstractTrajectory)
    isfull(t) || return

    states = t[:state]
    actions = t[:action]
    rewards = t[:reward]
    terminals = t[:terminal]
    next_state = select_last_frame(t[:next_state])

    AC = learner.approximator
    γ = learner.γ
    w₁ = learner.actor_loss_weight
    w₂ = learner.critic_loss_weight
    w₃ = learner.entropy_loss_weight
    D = device(AC)
    states = send_to_device(D, states)
    next_state = send_to_device(D, next_state)

    states_flattened = flatten_batch(states) # (state_size..., n_thread * update_step)
    actions = flatten_batch(actions)
    actions = CartesianIndex.(actions, 1:length(actions))

    next_state_values = AC.critic(next_state)
    gains = discount_rewards(
        rewards,
        γ;
        dims = 2,
        init = send_to_host(next_state_values),
        terminal = terminals,
    )
    gains = send_to_device(D, gains)

    ps = Flux.params(AC)
    gs = gradient(ps) do
        logits = AC.actor(states_flattened)
        probs = softmax(logits)
        log_probs = logsoftmax(logits)
        log_probs_select = log_probs[actions]
        values = AC.critic(states_flattened)
        advantage = vec(gains) .- vec(values)
        actor_loss = -mean(log_probs_select .* Zygote.dropgrad(advantage))
        critic_loss = mean(advantage .^ 2)
        entropy_loss = -sum(probs .* log_probs) * 1 // size(probs, 2)
        loss = w₁ * actor_loss + w₂ * critic_loss - w₃ * entropy_loss
        ignore() do
            learner.actor_loss = actor_loss
            learner.critic_loss = critic_loss
            learner.entropy_loss = entropy_loss
            learner.loss = loss
        end
        loss
    end
    if !isnothing(learner.max_grad_norm)
        learner.norm = clip_by_global_norm!(gs, ps, learner.max_grad_norm)
    end
    update!(AC, gs)
end

function (agent::Agent{<:QBasedPolicy{<:A2CLearner},<:CircularCompactSARTSATrajectory})(
    ::Training{PreActStage},
    env,
)
    action = agent.policy(env)
    state = get_state(env)
    push!(agent.trajectory; state = state, action = action)
    update!(agent.policy, agent.trajectory)

    # the main difference is we'd like to flush the buffer after each update!
    if isfull(agent.trajectory)
        empty!(agent.trajectory)
        push!(agent.trajectory; state = state, action = action)
    end

    action
end


text:
#custom extenstion of CudaArray in CUDArt for sparse vectors/matrices
#using CSC format for interop with Julia's native sparse functionality

import Base: length, size, ndims, eltype, similar, pointer, stride,
    copy, convert, reinterpret, show, summary, copy!, get!, fill!, issym,
    ishermitian, isupper, islower
import Base.LinAlg: BlasFloat, Hermitian, HermOrSym
import CUDArt: device, to_host, free

abstract AbstractCudaSparseArray{Tv,N} <: AbstractSparseArray{Tv,Cint,N}
typealias AbstractCudaSparseVector{Tv} AbstractCudaSparseArray{Tv,1}
typealias AbstractCudaSparseMatrix{Tv} AbstractCudaSparseArray{Tv,2}

type CudaSparseVector{Tv} <: AbstractCudaSparseVector{Tv}
    iPtr::CudaArray{Cint,1}
    nzVal::CudaArray{Tv,1}
    dims::NTuple{2,Int}
    nnz::Cint
    dev::Int
    function CudaSparseVector{Tv}(iPtr::CudaVector{Cint}, nzVal::CudaVector{Tv}, dims::Int, nnz::Cint, dev::Int)
        new(iPtr,nzVal,(dims,1),nnz,dev)
    end
end

type CudaSparseMatrixCSC{Tv} <: AbstractCudaSparseMatrix{Tv}
    colPtr::CudaArray{Cint,1}
    rowVal::CudaArray{Cint,1}
    nzVal::CudaArray{Tv,1}
    dims::NTuple{2,Int}
    nnz::Cint
    dev::Int

    function CudaSparseMatrixCSC{Tv}(colPtr::CudaVector{Cint}, rowVal::CudaVector{Cint}, nzVal::CudaVector{Tv}, dims::NTuple{2,Int}, nnz::Cint, dev::Int)
        new(colPtr,rowVal,nzVal,dims,nnz,dev)
    end
end

type CudaSparseMatrixCSR{Tv} <: AbstractCudaSparseMatrix{Tv}
    rowPtr::CudaArray{Cint,1}
    colVal::CudaArray{Cint,1}
    nzVal::CudaArray{Tv,1}
    dims::NTuple{2,Int}
    nnz::Cint
    dev::Int

    function CudaSparseMatrixCSR{Tv}(rowPtr::CudaVector{Cint}, colVal::CudaVector{Cint}, nzVal::CudaVector{Tv}, dims::NTuple{2,Int}, nnz::Cint, dev::Int)
        new(rowPtr,colVal,nzVal,dims,nnz,dev)
    end
end

type CudaSparseMatrixBSR{Tv} <: AbstractCudaSparseMatrix{Tv}
    rowPtr::CudaArray{Cint,1}
    colVal::CudaArray{Cint,1}
    nzVal::CudaArray{Tv,1}
    dims::NTuple{2,Int}
    blockDim::Cint
    dir::SparseChar
    nnz::Cint
    dev::Int

    function CudaSparseMatrixBSR(rowPtr::CudaVector{Cint}, colVal::CudaVector{Cint}, nzVal::CudaVector{Tv}, dims::NTuple{2,Int},blockDim::Cint, dir::SparseChar, nnz::Cint, dev::Int)
        new(rowPtr,colVal,nzVal,dims,blockDim,dir,nnz,dev)
    end
end

typealias cusparseHybMat_t Ptr{Void}
type CudaSparseMatrixHYB{Tv} <: AbstractCudaSparseMatrix{Tv}
    Mat::cusparseHybMat_t
    dims::NTuple{2,Int}
    nnz::Cint
    dev::Int

    function CudaSparseMatrixHYB(Mat::cusparseHybMat_t, dims::NTuple{2,Int}, nnz::Cint, dev::Int)
        new(Mat,dims,nnz,dev)
    end
end

typealias CompressedSparse{T} Union{CudaSparseMatrixCSC{T},CudaSparseMatrixCSR{T},HermOrSym{T,CudaSparseMatrixCSC{T}},HermOrSym{T,CudaSparseMatrixCSR{T}}}
typealias CudaSparseMatrix{T} Union{CudaSparseMatrixCSC{T},CudaSparseMatrixCSR{T}, CudaSparseMatrixBSR{T}, CudaSparseMatrixHYB{T}}

Hermitian{T}(Mat::CudaSparseMatrix{T}) = Hermitian{T,typeof(Mat)}(Mat,'U')

length(g::CudaSparseVector) = prod(g.dims)
size(g::CudaSparseVector) = g.dims
ndims(g::CudaSparseVector) = 1
length(g::CudaSparseMatrix) = prod(g.dims)
size(g::CudaSparseMatrix) = g.dims
ndims(g::CudaSparseMatrix) = 2

function size{T}(g::CudaSparseVector{T}, d::Integer)
    d >= 1 ? (d <= 1 ? g.dims[d] : 1) : throw(ArgumentError("dimension must be 1, got $d"))
end
function size{T}(g::CudaSparseMatrix{T}, d::Integer)
    d >= 1 ? (d <= 2 ? g.dims[d] : 1) : throw(ArgumentError("dimension must be ≥ 1, got $d"))
end

issym{T}(M::Union{CudaSparseMatrixCSC{T},CudaSparseMatrixCSR{T}})       = false
ishermitian{T}(M::Union{CudaSparseMatrixCSC{T},CudaSparseMatrixCSR{T}}) = false
issym{T}(M::Symmetric{T,CudaSparseMatrixCSC{T}})       = true
ishermitian{T}(M::Hermitian{T,CudaSparseMatrixCSC{T}}) = true

isupper{T}(M::UpperTriangular{T,CudaSparseMatrixCSC{T}}) = true
islower{T}(M::UpperTriangular{T,CudaSparseMatrixCSC{T}}) = false
isupper{T}(M::UpperTriangular{T,CudaSparseMatrixCSR{T}}) = true
islower{T}(M::UpperTriangular{T,CudaSparseMatrixCSR{T}}) = false
isupper{T}(M::UpperTriangular{T,CudaSparseMatrixHYB{T}}) = true
islower{T}(M::UpperTriangular{T,CudaSparseMatrixHYB{T}}) = false
isupper{T}(M::UpperTriangular{T,CudaSparseMatrixBSR{T}}) = true
islower{T}(M::UpperTriangular{T,CudaSparseMatrixBSR{T}}) = false
isupper{T}(M::LowerTriangular{T,CudaSparseMatrixCSC{T}}) = false
islower{T}(M::LowerTriangular{T,CudaSparseMatrixCSC{T}}) = true
isupper{T}(M::LowerTriangular{T,CudaSparseMatrixCSR{T}}) = false
islower{T}(M::LowerTriangular{T,CudaSparseMatrixCSR{T}}) = true
isupper{T}(M::LowerTriangular{T,CudaSparseMatrixHYB{T}}) = false
islower{T}(M::LowerTriangular{T,CudaSparseMatrixHYB{T}}) = true
isupper{T}(M::LowerTriangular{T,CudaSparseMatrixBSR{T}}) = false
islower{T}(M::LowerTriangular{T,CudaSparseMatrixBSR{T}}) = true

eltype{T}(g::CudaSparseMatrix{T}) = T
device(A::CudaSparseMatrix)       = A.dev
device(A::SparseMatrixCSC)        = -1  # for host

if VERSION >= v"0.5.0-dev+742"
    function to_host{T}(Vec::CudaSparseVector{T})
        SparseVector(Vec.dims[1], to_host(Vec.iPtr), to_host(Vec.nzVal))
    end
else
    function to_host{T}(Vec::CudaSparseVector{T})
        sparsevec(to_host(Vec.iPtr), to_host(Vec.nzVal), Vec.dims[1])
    end
end

function to_host{T}(Mat::CudaSparseMatrixCSC{T})
    SparseMatrixCSC(Mat.dims[1], Mat.dims[2], to_host(Mat.colPtr), to_host(Mat.rowVal), to_host(Mat.nzVal))

end
function to_host{T}(Mat::CudaSparseMatrixCSR{T})
    rowPtr = to_host(Mat.rowPtr)
    colVal = to_host(Mat.colVal)
    nzVal = to_host(Mat.nzVal)
    #construct Is
    I = similar(colVal)
    counter = 1
    for row = 1 : size(Mat)[1], k = rowPtr[row] : (rowPtr[row+1]-1)
        I[counter] = row
        counter += 1
    end
    return sparse(I,colVal,nzVal,Mat.dims[1],Mat.dims[2])
end

summary(g::CudaSparseMatrix) = string(g)
summary(g::CudaSparseVector) = string(g)

CudaSparseVector{T<:BlasFloat,Ti<:Integer}(iPtr::Vector{Ti}, nzVal::Vector{T}, dims::Int) = CudaSparseVector{T}(CudaArray(convert(Vector{Cint},iPtr)), CudaArray(nzVal), dims, convert(Cint,length(nzVal)), device())
CudaSparseVector{T<:BlasFloat,Ti<:Integer}(iPtr::CudaArray{Ti}, nzVal::CudaArray{T}, dims::Int) = CudaSparseVector{T}(iPtr, nzVal, dims, convert(Cint,length(nzVal)), device())

CudaSparseMatrixCSC{T<:BlasFloat,Ti<:Integer}(colPtr::Vector{Ti}, rowVal::Vector{Ti}, nzVal::Vector{T}, dims::NTuple{2,Int}) = CudaSparseMatrixCSC{T}(CudaArray(convert(Vector{Cint},colPtr)), CudaArray(convert(Vector{Cint},rowVal)), CudaArray(nzVal), dims, convert(Cint,length(nzVal)), device())
CudaSparseMatrixCSC{T<:BlasFloat,Ti<:Integer}(colPtr::CudaArray{Ti}, rowVal::CudaArray{Ti}, nzVal::CudaArray{T}, dims::NTuple{2,Int}) = CudaSparseMatrixCSC{T}(colPtr, rowVal, nzVal, dims, convert(Cint,length(nzVal)), device())
CudaSparseMatrixCSC{T<:BlasFloat,Ti<:Integer}(colPtr::CudaArray{Ti}, rowVal::CudaArray{Ti}, nzVal::CudaArray{T}, nnz, dims::NTuple{2,Int}) = CudaSparseMatrixCSC{T}(colPtr, rowVal, nzVal, dims, nnz, device())

CudaSparseMatrixCSR{T}(rowPtr::CudaArray, colVal::CudaArray, nzVal::CudaArray{T}, dims::NTuple{2,Int}) = CudaSparseMatrixCSR{T}(rowPtr, colVal, nzVal, dims, convert(Cint,length(nzVal)), device())
CudaSparseMatrixCSR{T}(rowPtr::CudaArray, colVal::CudaArray, nzVal::CudaArray{T}, nnz, dims::NTuple{2,Int}) = CudaSparseMatrixCSR{T}(rowPtr, colVal, nzVal, dims, nnz, device())

CudaSparseMatrixBSR{T}(rowPtr::CudaArray, colVal::CudaArray, nzVal::CudaArray{T}, blockDim, dir, nnz, dims::NTuple{2,Int}) = CudaSparseMatrixBSR{T}(rowPtr, colVal, nzVal, dims, blockDim, dir, nnz, device())

if VERSION >= v"0.5.0-dev+742"
    CudaSparseVector(Vec::SparseVector)    = CudaSparseVector(Vec.nzind, Vec.nzval, size(Vec)[1])
    CudaSparseMatrixCSC(Vec::SparseVector)    = CudaSparseMatrixCSC([1], Vec.nzind, Vec.nzval, size(Vec))
end
CudaSparseVector(Mat::SparseMatrixCSC) = size(Mat,2) == 1 ? CudaSparseVector(Mat.rowval, Mat.nzval, size(Mat)[1]) : throw(ArgumentError())
CudaSparseMatrixCSC(Mat::SparseMatrixCSC) = CudaSparseMatrixCSC(Mat.colptr, Mat.rowval, Mat.nzval, size(Mat))
CudaSparseMatrixCSR(Mat::SparseMatrixCSC) = switch2csr(CudaSparseMatrixCSC(Mat))

similar(Vec::CudaSparseVector) = CudaSparseVector(copy(Vec.iPtr), similar(Vec.nzVal), Vec.dims[1])
similar(Mat::CudaSparseMatrixCSC) = CudaSparseMatrixCSC(copy(Mat.colPtr), copy(Mat.rowVal), similar(Mat.nzVal), Mat.nnz, Mat.dims)
similar(Mat::CudaSparseMatrixCSR) = CudaSparseMatrixCSR(copy(Mat.rowPtr), copy(Mat.colVal), similar(Mat.nzVal), Mat.nnz, Mat.dims)
similar(Mat::CudaSparseMatrixBSR) = CudaSparseMatrixBSR(copy(Mat.rowPtr), copy(Mat.colVal), similar(Mat.nzVal), Mat.blockDim, Mat.dir, Mat.nnz, Mat.dims)

function copy!(dst::CudaSparseVector, src::CudaSparseVector; stream=null_stream)
    if dst.dims != src.dims
        throw(ArgumentError("Inconsistent Sparse Vector size"))
    end
    copy!( dst.iPtr, src.iPtr )
    copy!( dst.nzVal, src.nzVal )
    dst.nnz = src.nnz
    dst
end

function copy!(dst::CudaSparseMatrixCSC, src::CudaSparseMatrixCSC; stream=null_stream)
    if dst.dims != src.dims
        throw(ArgumentError("Inconsistent Sparse Matrix size"))
    end
    copy!( dst.colPtr, src.colPtr )
    copy!( dst.rowVal, src.rowVal )
    copy!( dst.nzVal, src.nzVal )
    dst.nnz = src.nnz
    dst
end

function copy!(dst::CudaSparseMatrixCSR, src::CudaSparseMatrixCSR; stream=null_stream)
    if dst.dims != src.dims
        throw(ArgumentError("Inconsistent Sparse Matrix size"))
    end
    copy!( dst.rowPtr, src.rowPtr )
    copy!( dst.colVal, src.colVal )
    copy!( dst.nzVal, src.nzVal )
    dst.nnz = src.nnz
    dst
end

function copy!(dst::CudaSparseMatrixBSR, src::CudaSparseMatrixBSR; stream=null_stream)
    if dst.dims != src.dims
        throw(ArgumentError("Inconsistent Sparse Matrix size"))
    end
    copy!( dst.rowPtr, src.rowPtr )
    copy!( dst.colVal, src.colVal )
    copy!( dst.nzVal, src.nzVal )
    dst.dir = src.dir
    dst.nnz = src.nnz
    dst
end

function copy!(dst::CudaSparseMatrixHYB, src::CudaSparseMatrixHYB; stream=null_stream)
    if dst.dims != src.dims
        throw(ArgumentError("Inconsistent Sparse Matrix size"))
    end
    dst.Mat = src.Mat
    dst.nnz = src.nnz
    dst
end

copy(Vec::CudaSparseVector; stream=null_stream) = copy!(similar(Vec),Vec;stream=null_stream)
copy(Mat::CudaSparseMatrixCSC; stream=null_stream) = copy!(similar(Mat),Mat;stream=null_stream)
copy(Mat::CudaSparseMatrixCSR; stream=null_stream) = copy!(similar(Mat),Mat;stream=null_stream)
copy(Mat::CudaSparseMatrixBSR; stream=null_stream) = copy!(similar(Mat),Mat;stream=null_stream)


text:
function random_samples(rho::Float64, N::Int)
  if !(0.0 <= rho < 1.0)
    throw(DomainError())
  end
  rhobar = sqrt(1.0-rho*rho)
  ret = zeros(N)
  ret[1] = randn()
  for i in 2:N
    ret[i] = ret[i-1]*rho + rhobar * randn()
  end
  return ret
end

function testBinningObservable(numtest, MCS, rho = 0.5)
  inform("Start testing BinningObservable")

  nsigma = 3
  counts = zeros(Int, nsigma)

  progress = 0.05
  for i in 1:numtest
    sobs = SimpleObservable()
    bobs = BinningObservable()
    for x in random_samples(rho, MCS)
      sobs << x
      bobs << x
    end
    m = mean(bobs)
    @test_approx_eq_eps mean(sobs) m 1e-10
    m = abs(m)
    er = stderror(bobs)
    for isigma in 1:nsigma
      if m < isigma * er
        counts[isigma] += 1
      end
    end

    if i/numtest >= progress
      inform(round(Int, progress * 100), "% done.")
      progress += 0.05
    end
  end

  sigma = counts ./ numtest

  @show sigma

  ## 1sigma, 2sigma, 3sigma
  ## 0.6827, 0.9545, 0.9973
  @test 0.67 < sigma[1] < 0.69
  @test 0.95 < sigma[2] < 0.96
  @test 0.99 < sigma[3]
  inform("Succeed testing BinningObservable")
end

testBinningObservable(10_000, 16384)


text:
# ensure stdlib is in load path
function ensure_stdlib()
    separator = Sys.iswindows() ? ';' : ':'
    if "JULIA_LOAD_PATH" in keys(ENV)
        _add_path = false
        comps = split(ENV["JULIA_LOAD_PATH"], separator)
        for needed in split("@:@v#.#:@stdlib", separator)
            if !(needed in comps)
                push!(comps, needed)
                _add_path = true
            end
        end
        if _add_path
            ENV["JULIA_LOAD_PATH"] = join(comps, separator)
        end
    end
    nothing
end

const PROTOC_TESTCASES = [
    (
        files=["t1.proto"],
        check="include(\"t1_pb.jl\")",
        envs=[],
    ),
    (
        files=["t2.proto"],
        check="include(\"t2_pb.jl\")",
        envs=[],
    ),
    (
        files=["recursive.proto"],
        check="include(\"recursive_pb.jl\")",
        envs=[],
    ),
    (
        files=["recursive_julia_keyword.proto"],
        check="include(\"recursive_julia_keyword_pb.jl\")",
        envs=[],
    ),
    (
        files=["a.proto", "b.proto"],
        check="include(\"AB.jl\"); using .AB; using .AB.A, .AB.B",
        envs=[],
    ),
    (
        files=["p1proto.proto", "p2proto.proto"],
        check="include(\"P1.jl\"); include(\"P2.jl\"); using .P1; using .P2",
        envs=[],
    ),
    (
        files=["module_type_name_collision.proto"],
        check="include(\"Foo_pb.jl\"); using .Foo_pb",
        envs=["JULIA_PROTOBUF_MODULE_POSTFIX"=>"1"],
    ),
    (
        files=["packed2.proto"],
        check="include(\"packed2_pb.jl\")",
        envs=[],
    ),
    (
        files=["inf_nan.proto"],
        check="include(\"inf_nan_pb.jl\")",
        envs=[],
    ),
    (
        files=["map3.proto"],
        check="include(\"map3_pb.jl\"); @test meta(MapTest).ordered[3].jtyp <: Dict; mv(\"map3_pb.jl\", \"map3_dict_pb.jl\"; force=true);",
        envs=[],
    ),
    (
        files=["map3.proto"],
        check="include(\"map3_pb.jl\"); @test meta(MapTest).ordered[3].jtyp <: Array; mv(\"map3_pb.jl\", \"map3_array_pb.jl\"; force=true);",
        envs=["JULIA_PROTOBUF_MAP_AS_ARRAY"=>"1"],
    ),
    (
        files=["oneof3.proto"],
        check="include(\"oneof3_pb.jl\")",
        envs=[],
    ),
    (
        files=["packed3.proto"],
        check="include(\"packed3_pb.jl\")",
        envs=[],
    ),
    (
        files=["any_test.proto"],
        check="include(\"any_test_pb.jl\")",
        envs=[],
    ),
    (
        files=["svc3.proto"],
        check="include(\"svc3_pb.jl\")",
        envs=[],
    ),
]

function protoc_test(files, check, envs, outdir)
    srcdir = joinpath(@__DIR__, "proto")
    well_known_proto_srcdir = abspath(joinpath(@__DIR__, "..", "gen"))
    srcpaths = join([joinpath(srcdir, file) for file in files], ' ')
    testscript = joinpath(outdir, "testprotoc_run.jl")

    is_ci = "CI" in keys(ENV)

    open(testscript, "w") do os
        if "JULIA_LOAD_PATH" in keys(ENV)
            loadpath = ENV["JULIA_LOAD_PATH"]
            if Sys.iswindows()
                # escape separator for printing
                loadpath = replace(loadpath, "\\"=>"\\\\")
                outdir = replace(outdir, "\\"=>"\\\\")
                srcdir = replace(srcdir, "\\"=>"\\\\")
                well_known_proto_srcdir = replace(well_known_proto_srcdir, "\\"=>"\\\\")
                srcpaths = replace(srcpaths, "\\"=>"\\\\")
            end
            println(os, "ENV[\"JULIA_LOAD_PATH\"] = \"$loadpath\"")
        end
        println(os, "using ProtoBuf, Test")
        for (env_name, env_val) in envs
            println(os, "ENV[\"$env_name\"]=\"$env_val\"")
        end
        if is_ci
            println(os, "ENV[\"COVERAGE\"]=\"--code-coverage=user --inline=no\"")
        end
        println(os, "cd(\"$outdir\")")

        # escape ~ character in Windows temp folder names
        escaped_outdir = replace(outdir, "~"=>"\\~")

        println(os, "ProtoBuf.protoc(`-I=$srcdir -I=$well_known_proto_srcdir --julia_out=$escaped_outdir $srcpaths`)")
        println(os, check)
    end
    julia_fullpath = joinpath(Sys.BINDIR, Base.julia_exename())
    julia = is_ci ? `$julia_fullpath` : `$julia_fullpath --code-coverage=user --inline=no`
    proc = run(`$julia $testscript`)
    proc.exitcode
end

@info("testing protoc")
@testset "protoc" begin
    ensure_stdlib()
    mktempdir() do outdir
        for testcase in PROTOC_TESTCASES
            if Sys.iswindows() && !isempty(testcase.envs)
                @info("skipping protoc tests that need environment variables on windows")
            else
                @test protoc_test(testcase.files, testcase.check, testcase.envs, outdir) == 0
            end
        end
    end
end


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# ╔═╡ 4de0ee11-f8af-4865-8737-ba0dc5c3404e
begin
	using Pkg
	cd(joinpath(dirname(@__FILE__),".."))
    Pkg.activate(pwd())
	using Logging
	# for simulation
	using Random
	using SimJulia, ResumableFunctions
end

# ╔═╡ 9cda8f1c-2394-42bf-883c-4dbe5df8ae56
# Make cells wider
html"""<style>
/*              screen size more than:                     and  less than:                     */
@media screen and (max-width: 699px) { /* Tablet */
  /* Nest everything into here */
    main { /* Same as before */
        max-width: 1200px !important; /* Same as before */
        margin-right: 100px !important; /* Same as before */
    } /* Same as before*/

}

@media screen and (min-width: 700px) and (max-width: 1199px) { /* Laptop*/
  /* Nest everything into here */
    main { /* Same as before */
        max-width: 1200px !important; /* Same as before */
        margin-right: 100px !important; /* Same as before */
    } /* Same as before*/
}

@media screen and (min-width:1000px) and (max-width: 1920px) { /* Desktop */
  /* Nest everything into here */
    main { /* Same as before */
        max-width: 1000px !important; /* Same as before */
        margin-right: 100px !important; /* Same as before */
    } /* Same as before*/
}

@media screen and (min-width:1921px) { /* Stadium */
  /* Nest everything into here */
    main { /* Same as before */
        max-width: 1200px !important; /* Same as before */
        margin-right: 100px !important; /* Same as before */
    } /* Same as before*/
}
</style>
"""

# ╔═╡ c395df98-145a-11eb-1716-2de187df1a1a
md"""
# Logging
The [Logging](https://docs.julialang.org/en/v1/stdlib/Logging/index.html) module will be used for efficient debugging and testing during development.

A logger has its own lower bound on the `LogLevel` that it can show. In addition to this, there is a global setting that determines the lowest level that will be registered.

Below you find some practical examples of using this module.
"""

# ╔═╡ dafa45ae-1462-11eb-3338-037167917f4d
disable_logging(LogLevel(-5001))

# ╔═╡ 0560f32a-1462-11eb-0685-09e4f341ddf5
begin
	"""
		logdemo1()

	A small demo where everything is run on a single logger. Keep in mind when using the global logger, that its lowest level is `Info` (`LogLevel(0)`), so you won't see anything below.

	Also keep in mind that you need to modify the global settings before you can see anything below `Debug`. You can do this with `disable_logging(LogLevel(N))`.

	### keywords
	* logger: the logger you want to use. Defaults to the global logger.
	"""
	function logdemo1(args...; kwargs...)
		# direct all the following messages to my logger
		logger = get(kwargs, :logger, Logging.global_logger())
		with_logger(logger) do
			# print some information about the function
			@info "logdemo1 was invoked with:\n\t- args: $(args)\n\t- kwargs: $(kwargs)\n\t- using logger <$(logger)> (current logger's lowest level: $(Logging.min_enabled_level(logger)))"
			@debug "logdemo1 lowest level message that allows the debug level"
			@logmsg LogLevel(-2000) "not visible by default"
		end
	end

	nothing;
end

# ╔═╡ 8ddac37c-1465-11eb-17d4-fdce80dc78fe
begin
	println("DEMO 1a - USING LOGGING DEFAULT SETTINGS\n$("-"^70)\n\n")
	logdemo1(1,2, goedemorgen="bonjour")
	println("-"^70)
end

# ╔═╡ b05f57c2-1466-11eb-272f-5fc93c04c744
begin
	println("DEMO 1b - USING A SPECIFIC LOGGER THAT WILL SHOW THE LAST MESSAGE\n$("-"^70)\n\n")
	customlogger = Logging.SimpleLogger(stdout, LogLevel(-2000))
	logdemo1(3,4,logger=customlogger, goedenavond="bonsoir")
	println("-"^70)
end

# ╔═╡ 19b12fc6-1464-11eb-2fc1-cbb3f82479c5
begin
	"""
		logdemo2()

	A small demo where it is possible to direct the logs generated by a specific function to a file. The can be very handy for debugging purposes or analysis after a simulation. By default everything is run on a single logger.

	### keywords
	* logger: the logger you want to use. Defaults to the global logger.
	* myspecialfunlogfilename: if you want to log `myspecialfun` to a file specify its name. When not specified the global logger is used
	* myspecialfunlogfilemode: the mode you want to use to write to a file. Defaults to "w" (cf. [write modes](https://docs.julialang.org/en/v1/base/io-network/#Base.open))
	"""
	function logdemo2(args...; kwargs...)
		# direct all the following messages to my logger
		logger = get(kwargs, :logger, Logging.global_logger())
		with_logger(logger) do
			# log message from the
			@info "logdemo2 was invoked with:\n\t- args: $(args)\n\t- kwargs: $(kwargs)\n\t- using logger <$(logger)> (current logger's lowest level: $(Logging.min_enabled_level(logger)))"
			# verify if a special logger should be used
			if haskey(kwargs, :myspecialfunlogfilename)
				logname = kwargs[:myspecialfunlogfilename]
				logmode = get(kwargs, :myspecialfunlogfilemode, "w")
				io = open(logname, logmode)
				speciallogger = SimpleLogger(io)
			else
				speciallogger = logger
			end

			# run the function with the appropriate logger
			with_logger(speciallogger) do
				for i in 1:get(kwargs,:maxrep, 20)
					myspecialfun(i; speciallogger=speciallogger)
				end
			end

			# close io if required
			if haskey(kwargs, :myspecialfunlogfilename)
				close(io)
			end
		end
	end

	"""
		myspecialfun(args...; kwargs...)

	a function that generates log messages

	### keywords
	"""
	function myspecialfun(args...; kwargs...)
		@info "myspecialfun was invoked with:\n\t- args: $(args)\n\t- kwargs: $(kwargs)\n\t- using logger <$(kwargs[:speciallogger])> "
	end

	nothing
end

# ╔═╡ e3eae0ce-1462-11eb-2e02-fd2f746569d1
begin
	println("DEMO 2a - USING LOGGING DEFAULT SETTINGS\n$("-"^70)\n\n\n")
	logdemo2(1,2, goedemorgen="bonjour")
	println("-"^70)
end

# ╔═╡ 12710aa4-146b-11eb-034a-97b515563abc
begin
	println("DEMO 2b - USING LOGGING TO FILE SETTINGS\n$("-"^70)\n\n")
	logdemo2(1,2, myspecialfunlogfilename="demo2.log", maxrep=10)
	println("-"^70)
end

# ╔═╡ dd63ff16-146d-11eb-059c-0586e1f972a5
md"""
## Working with resumable functions
Using a resumable function, try to implement:
1. the pascal triangle: each iteration should return a line of the Pascal triangle.
2. a root finding method (e.g. square root of a number)

*Note*: a resumable function returns an iterator you need to call.
"""

# ╔═╡ 010af16a-1474-11eb-10dc-292a149988f1
begin
	@resumable function pascal()
		# two initial values
		a = [1]
		@yield a
		a = [1,1]
		@yield a
		# all the following values
		while true
			a = vcat([1],a[1:end-1] .+ a[2:end],[1])
			@yield a
		end
	end

	p = pascal()
	println("Pascal's triangle:")
	for i in 1:10
		println(p())
	end
end

# ╔═╡ 363eba24-1474-11eb-26b8-718eed4e5f21
begin
	@resumable function findroot(x0::Number)
		# note: based on Newton's method
		res = x0
		@yield res
		while true
			res = res - (res^2 - x0) / (2*res)
			@yield res
		end
	end

	r = findroot(5)
	println("evolution of absolute error")
	for i in 0:10
		println("Iteration $(i): $(abs(sqrt(5) - r()) )")
	end
end

# ╔═╡ 0368ea70-145b-11eb-0b5c-fb3a33bf2027
md"""
# SimJulia
Before starting a larger project, we will look into some SimJulia tricks.

There are some compatibility issues between Pluto and more complex SimJulia constructions, which is why you will find specific examples in a seperate file.

You can execute these file by running them from the REPL by using
```julia
include("path/to/file.jl")
```

"""

# ╔═╡ 9da38750-146d-11eb-0757-2318eb6520ca
md"""
## Working with `containers`

Containers represent a level of something (e.g. liquid level, energy ...). If you want to store a specific type of object, you will be better of using a `store`.

Experiment a bit with containers (::Container). Discover their attributes (environment, capacity, level, get\_queue, put\_queue, seid) and find out how to use them. Generate a simple setting with:
1. a fill process that waits for a random time $t < 10 \in \mathbb{N}$ and then adds 1 unit to a container. This process repeats forever.
2. an empty process that waits for a random time $t < 10 \in \mathbb{N}$ and then requires a random amount from the container. This process repeats forever.
3. a monitor proces that periodically prints an info message detailing the current level of the container. This process repeats forever.
"""

# ╔═╡ 73f24a8e-146f-11eb-2978-cfef033adae1
let
	@resumable function fill(sim::Simulation, c::Container)
		while true
			@yield timeout(sim,rand(1:10))
			@yield put(c,1)
			@info "item added to the container on time $(now(sim))"
		end
	end

	@resumable function empty(sim::Simulation, c::Container)
		while true
			@yield timeout(sim,rand(1:10))
			n = rand(1:3)
			@info "Filed my request for $(n) items on time $(now(sim))"
			@yield get(c,n)
			@info "Got my $(n) items on time $(now(sim))"
		end
	end

	@resumable function monitor(sim::Simulation, c::Container)
		while true
			@info "$(now(sim)) - current container level: $(c.level)/$(c.capacity)"
			@yield timeout(sim,1)
		end
	end

	# setup the simulation
	# fix random seed for reproduction

	Random.seed!(174)

	@info "\n$("-"^70)\nWorking with containers\n$("-"^70)\n"
	sim = Simulation()
	c = Container(sim,10)
	@process fill(sim,c)
	@process monitor(sim,c)
	@process empty(sim,c)
	run(sim,30)
end

# ╔═╡ 71911828-1470-11eb-3519-bb52522ed2c9
md"""
## Working with `Stores`
A store can hold objects (struct) that can be used by other processes. Let's reconsider the same small scale application we did with the containers, i.e. generate a simple setting and verify everything works as intended (e.g. a fill, empty and monitor process).
"""

# ╔═╡ fbfed3c4-1470-11eb-0f28-231b891d0d9b
let
	# our own type of object
	struct Object
		id::Int
	end

	@resumable function fill(sim::Simulation, s::Store)
		i = 0
		while true
			i += 1
			item = Object(i)
			@yield timeout(sim,rand(1:10))
			@yield put(s,item)
			@info "item $(item) added to the store on time $(now(sim))"
		end
	end

	@resumable function empty(sim::Simulation, s::Store)
		while true
			@yield timeout(sim,rand(1:10))
			n = rand(1:3)
			@info "Filed my request for $(n) items on time $(now(sim))"
			for _ in 1:n
				@yield get(s)
			end
			@info "Got my $(n) items on time $(now(sim))"
		end
	end

	@resumable function monitor(sim::Simulation, s::Store)
		while true
			@info "$(now(sim)) - current store level: $(length(s.items))/$(s.capacity)"
			@yield timeout(sim,1)
		end
	end

	# fix random seed for reproduction
	Random.seed!(174)

	# setup the simulation
	@info "\n$("-"^70)\nWorking with stores\n$("-"^70)\n"
	sim = Simulation()
	s = Store{Object}(sim, capacity=UInt(10))

	@process fill(sim, s)
	@process empty(sim, s)
	@process monitor(sim, s)
	run(sim,30)

end

# ╔═╡ 118eecdc-146d-11eb-1d4b-71b301d4d5e6
md"""
## Process dependencies
Below you have an illustration of a process waiting for another one to terminate before continuing.
"""

# ╔═╡ 363592fc-146d-11eb-2dde-d56c18095702
let
	@resumable function basic(sim::Simulation)
    	@info "Basic goes to work on time $(now(sim))"
    	p = @process bottleneck(sim)
    	@yield p
    	@info "Basic continues after bottleneck completion on time $(now(sim))"
	end

	@resumable function bottleneck(sim::Simulation)
   		@yield timeout(sim, 10)
	end

	@info "\n$("-"^70)\nProcess dependencies\n$("-"^70)\n"
	sim = Simulation()
	@process basic(sim)
	run(sim)
end

# ╔═╡ ddd43e36-1473-11eb-2544-8f9e1ac0f59c
md"""
## Linking a process to a type
We want to simulate the life of a puppy. It's life consists of three activities:
eating, sleeping and playing. This simple life can be disturbed (i.e. interrupted)
when it gets picked up by a human. If a puppy likes you, it might lick your face.
After being picked up, a puppy continues its life as before.

We create our own type `Puppy` that has an associated process that models its life.
We also create a type `Human` that has an associated process to pick up a random
puppy from time to time.

This application illustrates how you can interrupt an ongoing process and even do
something with the cause of the interruption (in this case keeping track of who
got liked by a puppy).

```julia
include("path/to/puppies.jl")
```

"""

# ╔═╡ 4f9fee58-1536-11eb-1a4d-bd94541cb767
md"""
### Waiting on more than one event

We want to simulate a number of machines that each produce a specific product and
put it in a warehouse (a `Store`). There is a combination process that needs one
of each generated products and combines it into another. The simulation ends when
the fictive container of combined goods is full.

We create our own type `Product` with two fields that allow to identify the kind of
product and to identify its "creator" by means of a serial number.
We also create a type `Machine` that has an associated process. This machine has
an ID an make one type of products. All generated product are put into the same
store.

The combiner process works as follows:
* generate the events matching the availability of one of each product
* wait until all of these events are realised. *Note*: if you only have two events, you can simply use `ev1` & `ev2`.
* generate the event of putting a fictive combined product in a container
* @yield the event (i.e. time-out until done)
* verify if the container is full and if so stop the simulation on this event.

This application illustrates how you can wait on multiple other events before
continuing the simulation. Keep in mind that this requires ALL events to be
processed.

```julia
include("path/to/machines.jl")
```
"""

# ╔═╡ 63db05d2-1546-11eb-169c-7b23fe0009c4
md"""
### What if only one event needs to be realised?
Suppose an agent requests a resource but only has a limited amount of patience before no longer wanting/needing the resource.

For the example, a simulation is made with a `::Resource` with a capacity of $0$. So the agent can never obtain the requested resource. In the `agent` function the following happens:
1. A request for `r::Resource` is made. The type of `req` is `SimJulia.Put`. This event will be triggered by an `@yield`
2. the variable `res` is a dictionary with the events as key and the `::StateValue` as value. The first event to have been processed will have its `::StateValue` equal to `SimJulia.processed`
3. the `if` conditions test whether the `::StateValue` of our request is equal to `SimJulia.processed`.
  1. If this is the case, the agent obtains the `::Resource`, uses it for 1 time unit and releases it back for further use.
  2. If this is NOT the case, the other event will have taken place (in this case the timeout) and we remove the request from the `::Resource` queue with `cancel`.
4. the simulation terminates since no more processes are active on time 4.0.

"""

# ╔═╡ 9abf31ea-1546-11eb-3ff2-41ad2484a04b
let
	@resumable function agent(env::Environment,r::Resource)
		req = request(r)
		res = @yield req | timeout(env, 4)
		if res[req].state == SimJulia.processed
			@info "$(env.time) - Agent is using the resource..."
			@yield timeout(env,1)
			release(r)
			@info "$(env.time) - Agent released the resource."
		else
			@info "$(env.time) - Patience ran out..."
			cancel(r, req)
		end
	end

	function runsim()
		@info "\n$("-"^70)\nOne of both events\n$("-"^70)\n"
		sim = Simulation()
		r = Resource(sim,0)
		@process agent(sim, r)
		run(sim,10)
	end

	runsim()
end

# ╔═╡ a955fd5c-1536-11eb-0405-9bc433188115
md"""
### Using the first available resource

We want to simulate a number of warehouses that store the same product. At a regular
interval, a product is required. But the origin of the product does not matter.

We create our own type `Warehouse` with two field that allow to identify the
warehouse and that allow to track its stock by by means of a `Store`.

The production process works adds a random quantity to a random warehouse
(the first available) and works as follows:
* generate a product
* generate the requests for all resources
* yield the requests the will occur first. *note*: if two events occur at the same time, they will both happen. We deal with this later.
* cancel all the other requests that have not occured yet. For those that have occured,
we decrement the store with the value it was increased by.

The simulation stops when either all warehouses are full (or cannot handle the
produced quantity).

A similar approach can be followed when dealing with `get` requests instead of
`put` requests.

This application illustrates how you can deal with resource concurrency,
i.e. taking whatever resource(s) come(s) available first without blocking the
other ones or introducing unwanted artifacts.

```julia
include("path/to/warehouse.jl")
```

"""

# ╔═╡ 414b4ee2-1473-11eb-3d7b-ef4f49e7efa0
md"""
## Exercises

"""

# ╔═╡ 4bf1fce0-1470-11eb-1290-63d06c8246a2
md"""
#### Application 1
Consider a candy machines that is continuously being monitored by a supervisor.  If the level is below a given treshold, the supervisor fills the machine up.
* Client arrival follows an exponential distribution with parameter $\theta = 1$ and each client takes two candies at a time.
* Look at the mean time between refills. Is this what you would expect?
* What happens when the amount of candy varies?  Is this still what you would expect? E.g. a clients takes one, two or three candies.
"""

# ╔═╡ f3c8a4ba-1474-11eb-06b0-7f0e5ba47670

# ╔═╡ c7b95ece-150e-11eb-0058-c15fde632ea0
md"""
### Application 2
When modeling physical things such as cables, RF propagation, etc. encapsulation of this process is better in order to keep the propagation mechanism outside of the sending and receiving processes.

Consider the following:
* a sender sends messages on a regular interval (e.g. every 5 minutes)
* a receiver is listening on a permanent basis for new messages
* the transfer between both of them is not instantaneous, but takes some time. To model this, you can store (hint: use a `::Store`) the messages on the cable for later reception.

"""

# ╔═╡ 25a78e3a-1511-11eb-2239-1b938a2129fa

# ╔═╡ Cell order:
# ╟─9cda8f1c-2394-42bf-883c-4dbe5df8ae56
# ╠═4de0ee11-f8af-4865-8737-ba0dc5c3404e
# ╟─c395df98-145a-11eb-1716-2de187df1a1a
# ╠═dafa45ae-1462-11eb-3338-037167917f4d
# ╠═0560f32a-1462-11eb-0685-09e4f341ddf5
# ╠═8ddac37c-1465-11eb-17d4-fdce80dc78fe
# ╠═b05f57c2-1466-11eb-272f-5fc93c04c744
# ╠═19b12fc6-1464-11eb-2fc1-cbb3f82479c5
# ╠═e3eae0ce-1462-11eb-2e02-fd2f746569d1
# ╠═12710aa4-146b-11eb-034a-97b515563abc
# ╟─dd63ff16-146d-11eb-059c-0586e1f972a5
# ╠═010af16a-1474-11eb-10dc-292a149988f1
# ╠═363eba24-1474-11eb-26b8-718eed4e5f21
# ╟─0368ea70-145b-11eb-0b5c-fb3a33bf2027
# ╟─9da38750-146d-11eb-0757-2318eb6520ca
# ╠═73f24a8e-146f-11eb-2978-cfef033adae1
# ╟─71911828-1470-11eb-3519-bb52522ed2c9
# ╠═fbfed3c4-1470-11eb-0f28-231b891d0d9b
# ╟─118eecdc-146d-11eb-1d4b-71b301d4d5e6
# ╠═363592fc-146d-11eb-2dde-d56c18095702
# ╟─ddd43e36-1473-11eb-2544-8f9e1ac0f59c
# ╟─4f9fee58-1536-11eb-1a4d-bd94541cb767
# ╟─63db05d2-1546-11eb-169c-7b23fe0009c4
# ╠═9abf31ea-1546-11eb-3ff2-41ad2484a04b
# ╟─a955fd5c-1536-11eb-0405-9bc433188115
# ╟─414b4ee2-1473-11eb-3d7b-ef4f49e7efa0
# ╟─4bf1fce0-1470-11eb-1290-63d06c8246a2
# ╠═f3c8a4ba-1474-11eb-06b0-7f0e5ba47670
# ╟─c7b95ece-150e-11eb-0058-c15fde632ea0
# ╠═25a78e3a-1511-11eb-2239-1b938a2129fa

max_stars_repo_path: Exercises/PS07 - SimJulia.jl

max_issues_repo_path: Exercises/PS07 - SimJulia.jl

max_forks_repo_path: Exercises/PS07 - SimJulia.jl


text:
# use conventional beam forming

using Cairo, Fontconfig
using Gadfly
using JLD2, FileIO

using DSP
using Caesar

AMP.setForceEvalDirect!(true)

include(joinpath(dirname(@__FILE__),"ps3utils.jl"))

## Load the data
dataparent = joinpath(dirname(@__FILE__),"..","..","..","..","..","..","data","ps3")
datadir = joinpath(dataparent,"equilateral_22inches") #
# datadir = joinpath(dataparent,"demo")

# reload data

leftwavdata = loadPs3WavEnvelopes(joinpath(datadir,"left_wavdata_corrected.jld2"))
rightwavdata = loadPs3WavEnvelopes(joinpath(datadir,"right_wavdata_corrected.jld2"))

## Load the templates
templdir = joinpath(dataparent, "templateX200ms")
lefttemplate = loadPs3WavEnvelopes(joinpath(templdir, "template_left.jld2"), variable="templ", cutmax=false)
righttemplate = loadPs3WavEnvelopes(joinpath(templdir, "template_right.jld2"), variable="templ", cutmax=false)

leftwav = [leftwavdata[0] leftwavdata[1] leftwavdata[2] leftwavdata[3]]
rightwav = [rightwavdata[0] rightwavdata[1] rightwavdata[2] rightwavdata[3]]

## Spot check

Gadfly.plot(y=leftwav[:,1], Geom.line)
Gadfly.plot(y=rightwav[:,1], Geom.line)

## generate nonparametric bearing intensities from waveform data

intensityLeft = beamformBasic(leftwav, lefttemplate, soundSpeed=340)
intensityRight = beamformBasic(rightwav, righttemplate, soundSpeed=340)
# Gadfly.plot(x=azimuths, y=intensity, Geom.line)

pcLeft = intensityToCircularKDE(intensityLeft, SNRfloor=0.80, N=100)
pcRight = intensityToCircularKDE(intensityRight, SNRfloor=0.80, N=100)

pl = AMP.plotKDECircular([pcLeft; pcRight], scale=0.4, c=["blue";"green"])
pl.coord = Coord.Cartesian(aspect_ratio=1.0, xmin=-2, xmax=2, ymin=-2,ymax=2)
pl

##

pl |> PDF("/tmp/test.pdf",30cm,20cm)
pl |> SVG("/tmp/test.svg",10cm,10cm)

##

@async run(`evince /tmp/test.pdf`)

##

pl = AMP.plotKDECircular([pcLeft], scale=0.4, c=["red";"green"])

#


text:
using DifferentialEquations
using PyCall
using PyPlot

function aginglaw(v, theta, dc)
    return 1 - theta * v / dc
end

function sliplaw(v, theta, dc)
    return -v * theta / dc * log(v * θ / dc)
end

function calcdvθclassic!(du, u, p, t)
    dc, η, σn, a, b, μ, Vp, L, ρ, statelaw = p
    θ = u[1]
    v = u[2]
    du[1] = statelaw(v, θ, dc)
    du[2] = 1 / (η / σn + a / v) * (μ * (Vp - v) / (L * σn) - b * du[1] / θ)
    return nothing
end

function plottimeseries(sol, siay, titlelabel)
    t1 = [x / siay for x in sol.t]
    θ1 = [x[1] for x in sol.u]
    v1 = [x[2] for x in sol.u]

    figure(figsize = (12, 6))
    subplot(2, 2, 1)
    plot(t1, θ1, "-b", linewidth = 0.5)
    yscale("log")
    ylabel(L"\theta")

    subplot(2, 2, 2)
    plot(1:1:length(t1), θ1, "-b", linewidth = 0.5)
    yscale("log")
    ylabel(L"\theta")

    subplot(2, 2, 3)
    plot(t1, v1, "-b", linewidth = 0.5)
    yscale("log")
    xlabel("time (years)")
    ylabel("v (m/s)")

    subplot(2, 2, 4)
    plot(1:1:length(t1), v1, "-b", linewidth = 0.5)
    yscale("log")
    xlabel("time (step)")
    ylabel("v (m/s)")
    suptitle(titlelabel)

    return nothing
end

function sliding()
    # Model parameters
    siay = 365.25 * 24 * 60 * 60
    tspan = (0.0, siay * 1000.0) # This will run
    # tspan = (0.0, siay * 302.0) # This will run but fail for the a, b case
    μ = 3e10
    ν = 0.25
    ρ = 2700.0
    η = μ / (2.0 * sqrt(μ / ρ))
    L = 60 * 1e3
    a = 0.015
    b = 0.02
    vp = 1e-9
    σn = 30e6
    dc = 0.1
    abstol = 1e-4
    reltol = 1e-4

    # Parameters need for d(a, b σn)/dt
    fstar = 0.6 # ???
    vstar = vp # ???
    θstar = 1e9 # ???

    # Show Tushar scalings for a and b
    # tusharabmodels()

    # Time integrate - classic
    icsclassic = [1e8; vp / 1000]
    pclassic = (dc, η, σn, a, b, μ, vp, L, ρ, aginglaw)
    probclassic = ODEProblem(calcdvθclassic!, icsclassic, tspan, pclassic)
    solclassic = solve(probclassic, RK4(), abstol = abstol, reltol = reltol)
    plottimeseries(solclassic, siay, "RSF (classic)")

    # Time integrate - a, b can evolve
    # icsab = [1e8 ; vp / 1000 ; a ; b]
    # pab = (dc, η, σn, μ, vp, L, ρ, aginglaw, fstar, vstar, θstar)
    # probab = ODEProblem(calcdvθab!, icsab, tspan, pab)
    # solab = solve(probab, RK4(), abstol = abstol, reltol = reltol)
    # plottimeseries(solab, siay, "RSF (ab)")

    return nothing
end
sliding()


text:
__precompile__(true)
module CLBLAS

using OpenCL: cl

export axpy!, scal!, gemm!, gemv!, gbmv!

# why there is a type assertion at context.jl line 38
depsfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("CLBLAS not properly installed. Please run Pkg.build(\"CLBLAS\") then restart Julia.") # now that this is decoupled from images, should this be an error?
end

include("constants.jl")
include("macros.jl")

include("L1/L1.jl")
include("L2/L2.jl")
include("L3/L3.jl")

include("highlevel.jl")

function LocalMem(::Type{T}, len::Integer) where T
    @assert len > 0
    nbytes = sizeof(T) * len
    return LocalMem{T}(convert(Csize_t, nbytes))
end

# setup function
@blasfun clblasSetup()

#destroy
clblasTeardown() = ccall((:clblasTeardown, libCLBLAS), Void, ())

global compute_context_holder = []
global next_compute_context = -1
global clblas_initialized = false

function get_next_compute_context()
    global next_compute_context = next_compute_context + 1
    return get_next_compute_context(next_compute_context)
end

function get_next_compute_context(next_compute_context)
    return compute_context_holder[( next_compute_context % length(compute_context_holder)) + 1]
end

function supress_context_error(error_info, private_info)
    warn(error_info)
    return C_NULL
end

function setup(profile = false)
    if isempty(cl.platforms())
        throw(cl.OpenCLException("No OpenCL.Platform available"))
    end
    #TODO: osx only due to https://github.com/clMathLibraries/clBLAS/issues/25
    available_devices = is_apple() ? cl.devices(:gpu) : cl.devices()
    if isempty(available_devices)
        throw(cl.OpenCLException("Unable to create any OpenCL.Context, no available devices"))
    end
    if !isempty(available_devices)
        for dev in available_devices
            local ctx::cl.Context
            local queue::cl.CmdQueue
            try
                ctx = cl.Context(dev, properties=nothing, callback=supress_context_error)
                queue  = profile ? cl.CmdQueue(ctx, :profile) : cl.CmdQueue(ctx);
                push!(compute_context_holder, (dev, ctx, queue))
            catch err
                error(err)
            end
        end
    end
    if isempty(compute_context_holder)
        throw(cl.OpenCLException("Unable to create any OpenCL.Context, not worked..."))
    end
    return
end

function teardown()
    if clblas_initialized
        clblasTeardown()
    end
    global compute_context_holder = []
    global next_compute_context = -1
    global clblas_initialized = false
    return
end

function __init__()
    clblasSetup()
    global clblas_initialized = true
    atexit(teardown)
end

end # module


text:
module Model

export C, V, param_values, get_derivatives, get_steady_state

using SteadyStateDiffEq;
using Sundials;

include("name2idx/parameters.jl")
include("name2idx/species.jl")
include("param_const.jl");
include("param_derivative.jl");
include("steady_state.jl");

end # module

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


text:
# This file is mmodified from Julia. License is MIT: https://julialang.org/license

const PosInfinity = Union{Infinity, InfiniteCardinal{0}}

(:)(start::T, stop::PosInfinity) where {T<:Integer} = InfUnitRange{T}(start)
(:)(start::PosInfinity, stop::Integer) = start+1:start
function (:)(start::T, step::T, stop::RealInfinity) where {T<:Real}
    signbit(step) == signbit(stop) || throw(ArgumentError("InfStepRange must have infinite length"))
    InfStepRange(start, step)
end
(:)(start::T, step::Real, stop::RealInfinity) where {T<:Real} = (:)(promote(start, step)..., stop)
(:)(start::Real, step, stop::PosInfinity)= (:)(start, step, RealInfinity(stop))
(:)(::PosInfinity, _, ::Real) = throw(ArgumentError("Cannot create range starting at infinity"))

# AbstractFloat specializations
(:)(a::T, b::Union{PosInfinity,RealInfinity}) where {T<:Real} = (:)(a, T(1), b)

function (:)(start::T, step::T, stop::PosInfinity) where {T<:Real}
    sign(step) == sign(stop) || throw(ArgumentError("InfStepRange must have infinite length"))
    InfStepRange(start,step)
end

# this is needed for showarray
(:)(::PosInfinity, ::PosInfinity) = 1:0

# Range of a given length: range(a, [step=s,] length=l), no stop
_range(a::Real,          ::Nothing,         ::Nothing, len::InfiniteCardinal{0}) = InfUnitRange{typeof(a)}(a)
_range(a::AbstractFloat, ::Nothing,         ::Nothing, len::InfiniteCardinal{0}) = _range(a, oftype(a, 1),   nothing, len)
_rangestyle(::Ordered, ::ArithmeticWraps, a::T, step::S, len::InfiniteCardinal{0}) where {T,S} =
    InfStepRange{T,S}(a, step)
_rangestyle(::Ordered, ::ArithmeticUnknown, a::T, step::S, len::InfiniteCardinal{0}) where {T,S} =
    InfStepRange{T,S}(a, step)
_rangestyle(::Any, ::Any, a::T, step::S, len::InfiniteCardinal{0}) where {T,S} =
    InfStepRange{T,S}(a, step)
_range(a::T, st::T, ::Nothing, ::InfiniteCardinal{0}) where T<:Union{Float16,Float32,Float64} =
    InfStepRange{T,T}(a, st)
_range(a::T, st::T, ::Nothing, ::InfiniteCardinal{0}) where T<:AbstractFloat =
    InfStepRange{T,T}(a, st)

# Construct range for rational start=start_n/den, step=step_n/den
floatrange(::Type{T}, start_n::Integer, step_n::Integer, ::PosInfinity, den::Integer) where T =
    InfStepRange(T(start_n)/den,T(step_n)/den)

floatrange(a::AbstractFloat, st::AbstractFloat, ::PosInfinity, divisor::AbstractFloat) =
    InfStepRange(a/divisor,st/divisor)

## 1-dimensional ranges ##

struct InfStepRange{T,S} <: OrdinalRange{T,S}
    start::T
    step::S
    function InfStepRange{T,S}(start::T, step::S) where {T,S}
        new{T,S}(start, step)
    end
end

InfStepRange(start::T, step::S) where {T,S} = InfStepRange{T,S}(start,step)
InfStepRange{T,S}(start, step) where {T,S} = InfStepRange{T,S}(convert(T,start),convert(S,step))

abstract type AbstractInfUnitRange{T<:Real} <: AbstractUnitRange{T} end

done(r::AbstractInfUnitRange{T}, i) where {T} = false
unitrange_last(start, stop::PosInfinity) = ∞

struct InfUnitRange{T<:Real} <: AbstractInfUnitRange{T}
    start::T
end

InfUnitRange(a::InfUnitRange) = a
InfUnitRange{T}(a::AbstractInfUnitRange) where T<:Real = InfUnitRange{T}(first(a))
InfUnitRange(a::AbstractInfUnitRange{T}) where T<:Real = InfUnitRange{T}(first(a))
unitrange(a::AbstractInfUnitRange) = InfUnitRange(a)

AbstractArray{T}(a::InfUnitRange) where T<:Real = InfUnitRange{T}(a.start)
AbstractVector{T}(a::InfUnitRange) where T<:Real = InfUnitRange{T}(a.start)
AbstractArray{T}(a::InfStepRange) where T<:Real = InfStepRange(convert(T,a.start), convert(T,a.step))
AbstractVector{T}(a::InfStepRange) where T<:Real = InfStepRange(convert(T,a.start), convert(T,a.step))

const InfRanges{T} = Union{InfStepRange{T},AbstractInfUnitRange{T}}
const InfAxes = Union{InfRanges{<:Integer},Slice{<:AbstractInfUnitRange{<:Integer}},IdentityUnitRange{<:AbstractInfUnitRange{<:Integer}}}

AbstractArray{T}(ac::Adjoint{<:Any,<:InfRanges}) where T<:Real = AbstractArray{T}(parent(ac))'
AbstractMatrix{T}(ac::Adjoint{<:Any,<:InfRanges}) where T<:Real = AbstractVector{T}(parent(ac))'
AbstractArray{T}(ac::Transpose{<:Any,<:InfRanges}) where T<:Real = transpose(AbstractArray{T}(parent(ac)))
AbstractMatrix{T}(ac::Transpose{<:Any,<:InfRanges}) where T<:Real = transpose(AbstractVector{T}(parent(ac)))

copy(ac::AdjOrTrans{<:Any,<:InfRanges}) = ac
reverse(a::InfRanges) = throw(ArgumentError("Cannot reverse infinite range"))

if VERSION ≥ v"1.7-"
    Base.range_start_step_length(a::T, st::T, len::PosInfinity) where T<:Union{Float16,Float32,Float64} = InfStepRange(a, st)
    for op in (:+, :-)
        @eval $op(a::InfRanges, b::InfRanges) = InfStepRange($op(first(a),first(b)), $op(step(a),step(b)))
    end
    broadcasted(::DefaultArrayStyle{1}, ::typeof(*), x::Number, r::InfRanges) = InfStepRange(x*first(r), x*step(r))
    broadcasted(::DefaultArrayStyle{1}, ::typeof(*), r::InfRanges, x::Number) = InfStepRange(first(r)*x, step(r)*x)
    broadcasted(::DefaultArrayStyle{1}, ::typeof(*), x::AbstractFloat, r::InfRanges) = InfStepRange(x*first(r), x*step(r))
    broadcasted(::DefaultArrayStyle{1}, ::typeof(*), r::InfRanges, x::AbstractFloat) = InfStepRange(first(r)*x, step(r)*x)
end

"""
    OneToInf(n)

Define an `AbstractInfUnitRange` that behaves like `1:∞`, with the added
distinction that the limits are guaranteed (by the type system) to
be 1 and ∞.
"""
struct OneToInf{T<:Integer} <: AbstractInfUnitRange{T} end

OneToInf() = OneToInf{Int}()
oneto(::PosInfinity) = OneToInf()
function oneto(x::ComplexInfinity)
    iszero(angle(x)) && return oneto(∞)
    throw(ArgumentError("Cannot create infinite range with negative length"))
 end
 function oneto(x::RealInfinity)
    signbit(x) || return oneto(∞)
    throw(ArgumentError("Cannot create infinite range with negative length"))
 end

AbstractArray{T}(a::OneToInf) where T<:Integer = OneToInf{T}()
AbstractVector{T}(a::OneToInf) where T<:Integer = OneToInf{T}()
AbstractArray{T}(a::OneToInf) where T<:Real = InfUnitRange{T}(a)
AbstractVector{T}(a::OneToInf) where T<:Real = InfUnitRange{T}(a)

(==)(::OneToInf, ::OneToInf) = true

## interface implementations

size(r::InfRanges) = (ℵ₀,)

isempty(r::InfRanges) = false

step(r::InfStepRange) = r.step

length(r::InfRanges) = ℵ₀
unsafe_length(r::InfRanges) = ℵ₀

first(r::OneToInf{T}) where {T} = oneunit(T)

last(r::AbstractInfUnitRange) = ℵ₀
last(r::InfStepRange) = RealInfinity(signbit(step(r)))

minimum(r::InfUnitRange) = first(r)
maximum(r::InfUnitRange) = ℵ₀

## iteration
start(r::InfStepRange) = oftype(r.start + r.step, r.start)
next(r::InfStepRange{T}, i) where {T} = (convert(T,i), i+r.step)

start(r::InfUnitRange{T}) where {T} = oftype(r.start + oneunit(T), r.start)
start(r::OneToInf{T}) where {T} = oneunit(T)

done(r::InfStepRange{T}, i) where {T} = false

## iteration with zip + finite iterator
#  allows axes(zip(...)) and size(zip(...))
Base.Iterators._promote_tuple_shape((a,)::Tuple{OneToInf}, (b,)::Tuple{OneTo}) =
    (intersect(a, b),)
Base.Iterators._promote_tuple_shape((a,)::Tuple{OneTo}, (b,)::Tuple{OneToInf}) =
    (intersect(a, b),)

## indexing

unsafe_indices(S::Slice{<:OneToInf}) = (S.indices,)
axes(S::Slice{<:OneToInf}) = (S.indices,)
_sub2ind(inds::Tuple{OneToInf}, i::Integer)    = i

to_shape(::OneToInf) = ℵ₀

# used for linear indexing
function _ind2sub_recurse(inds::Tuple{OneToInf{Int},Vararg{Any}}, ind::Integer)
    @_inline_meta
    (ind+1, _ind2sub_recurse(tail(inds), 0)...)
end

function getindex(v::InfUnitRange{T}, i::Integer) where T
    @boundscheck i > 0 || Base.throw_boundserror(v, i)
    convert(T, first(v) + i - 1)
end
function getindex(v::OneToInf{T}, i::Integer) where T
    @boundscheck i > 0 || Base.throw_boundserror(v, i)
    convert(T, i)
end
function getindex(v::InfStepRange{T}, i::Integer) where T
    @boundscheck i > 0 || Base.throw_boundserror(v, i)
    convert(T, first(v) + (i - 1)*step(v))
end

function getindex(x::AbstractUnitRange, y::PosInfinity)
    isinf(length(x)) || throw(BoundsError(x,y))
    ℵ₀
end

function getindex(x::OneToInf{T}, y::PosInfinity) where T
    isinf(length(x)) || throw(BoundsError(x,y))
    ℵ₀
end

function getindex(x::InfStepRange{T}, y::PosInfinity) where T
    isinf(length(x)) || throw(BoundsError(x,y))
    ℵ₀
end
function getindex(x::InfUnitRange{T}, y::PosInfinity) where T
    isinf(length(x)) || throw(BoundsError(x,y))
    ℵ₀
end

getindex(::AbstractInfUnitRange, ::Infinity) = ℵ₀
getindex(::OneToInf, ::Infinity) = ℵ₀
getindex(v::InfUnitRange, i::Infinity) = ℵ₀
getindex(v::InfStepRange, i::Infinity) = ℵ₀

function getindex(r::AbstractInfUnitRange, s::AbstractInfUnitRange{<:Integer})
    f = first(r)
    @boundscheck first(s) ≥ 1 || Base.throw_boundserror(r, first(s))
    st = oftype(f, f + first(s)-1)
    st:ℵ₀
end

function getindex(r::AbstractInfUnitRange, s::AbstractUnitRange{<:Integer})
    f = first(r)
    @boundscheck first(s) ≥ 1 || Base.throw_boundserror(r, first(s))
    st = oftype(f, f + first(s)-1)
    range(st; length=length(s))
end

getindex(r::AbstractInfUnitRange, s::Slice{<:AbstractInfUnitRange{<:Integer}}) = r

getindex(r::OneToInf{T}, s::OneTo) where T = OneTo(T(s.stop))

function getindex(r::AbstractInfUnitRange, s::InfStepRange{<:Integer})
    @_inline_meta
    @boundscheck (step(s) > 0 && first(s) ≥ 1) || throw(BoundsError(r, minimum(s)))
    st = oftype(first(r), first(r) + s.start-1)
    new_step = step(s)
    InfStepRange(st,new_step)
end

function getindex(r::AbstractInfUnitRange, s::StepRange{<:Integer})
    @_inline_meta
    @boundscheck minimum(s) ≥ 1 || throw(BoundsError(r, minimum(s)))
    st = oftype(first(r), first(r) + s.start-1)
    range(st; step=step(s), length=length(s))
end

function getindex(r::InfStepRange, s::InfAxes)
    @_inline_meta
    @boundscheck (step(s) > 0 && first(s) ≥ 1) || throw(BoundsError(r, minimum(s)))
    st = oftype(r.start, r.start + (first(s)-1)*step(r))
    new_step = step(r)*step(s)
    InfStepRange(st,new_step)
end

function getindex(r::InfStepRange, s::AbstractRange{<:Integer})
    @_inline_meta
    @boundscheck isempty(s) || minimum(s) ≥ 1 || throw(BoundsError(r, minimum(s)))
    st = oftype(r.start, r.start + (first(s)-1)*step(r))
    range(st; step=step(r)*step(s), length=length(s))
end

function getindex(Ac::AdjOrTrans{<:Any,<:InfRanges}, k::Integer, j)
    @boundscheck k == 1 || throw(BoundsError(Ac, k))
    parent(Ac)[j]
end

function getindex(Ac::AdjOrTrans{<:Any,<:InfRanges}, k::Integer, j::InfAxes)
    @boundscheck k == 1 || throw(BoundsError(Ac, k))
    parent(Ac)[j]
end

show(io::IO, r::InfUnitRange) = print(io, repr(first(r)), ':', repr(∞))
show(io::IO, r::OneToInf{Int}) = print(io, "OneToInf()")
show(io::IO, r::OneToInf{T}) where T = print(io, "OneToInf{$T}()")

intersect(r::OneToInf{T}, s::OneToInf{V}) where {T,V} = OneToInf{promote_type(T,V)}()
intersect(r::OneToInf{T}, s::OneTo{T}) where T = s
intersect(r::OneTo{T}, s::OneToInf{T}) where T = r

intersect(r::AbstractInfUnitRange{<:Integer}, s::AbstractInfUnitRange{<:Integer}) =
    InfUnitRange(max(first(r),first(s)))
intersect(r::AbstractInfUnitRange{<:Integer}, s::AbstractUnitRange{<:Integer}) =
    max(first(r),first(s)):last(s)
intersect(r::AbstractUnitRange{<:Integer}, s::AbstractInfUnitRange{<:Integer}) =
    intersect(s, r)

intersect(i::Integer, r::AbstractInfUnitRange{<:Integer}) =
    i < first(r) ? (first(r):i) : (i:i)

intersect(r::AbstractInfUnitRange{<:Integer}, i::Integer) = intersect(i, r)

function intersect(r::AbstractInfUnitRange{<:Integer}, s::StepRange{<:Integer})
    if isempty(s)
        range(first(r); length=0)
    elseif step(s) == 0
        intersect(first(s), r)
    elseif step(s) < 0
        intersect(r, reverse(s))
    else
        sta = first(s)
        ste = step(s)
        sto = last(s)
        lo = first(r)
        i0 = max(sta, lo + mod(sta - lo, ste))
        i1 = sto
        i0:ste:i1
    end
end

intersect(r::StepRange{<:Integer}, s::AbstractInfUnitRange{<:Integer}) =
    intersect(s, r)

function intersect(r::InfStepRange, s::InfStepRange)
    sign(step(r)) == sign(step(s)) || throw(ArgumentError("Can only take intersection of infinite ranges"))
    start1 = first(r)
    step1 = step(r)
    start2 = first(s)
    step2 = step(s)
    a = lcm(step1, step2)

    g, x, y = gcdx(step1, step2)

    if rem(start1 - start2, g) != 0
        # Unaligned, no overlap possible.
        throw(ArgumentError("Cannot take intersection of InfStepRange that has no elements"))
    end

    z = div(start1 - start2, g)
    b = start1 - x * z * step1
    # Possible points of the intersection of r and s are
    # ..., b-2a, b-a, b, b+a, b+2a, ...
    # Determine where in the sequence to start and stop.
    m = max(start1 + mod(b - start1, a), start2 + mod(b - start2, a))
    m:a:∞
end

intersect(r::InfStepRange, s::AbstractInfUnitRange) = intersect(r, InfStepRange(s))
intersect(r::AbstractInfUnitRange, s::InfStepRange) = intersect(InfStepRange(r), s)

function intersect(r::InfStepRange, s::StepRange)
    if isempty(s)
        return range(first(r); step=step(r), length=0)
    elseif step(s) < 0
        return intersect(r, reverse(s))
    end

    start1 = first(r)
    step1 = step(r)
    start2 = first(s)
    step2 = step(s)
    stop2 = last(s)
    a = lcm(step1, step2)

    g, x, y = gcdx(step1, step2)

    if rem(start1 - start2, g) != 0
        # Unaligned, no overlap possible.
        return range(start1; step=a, length=0)
    end

    z = div(start1 - start2, g)
    b = start1 - x * z * step1
    # Possible points of the intersection of r and s are
    # ..., b-2a, b-a, b, b+a, b+2a, ...
    # Determine where in the sequence to start and stop.
    m = max(start1 + mod(b - start1, a), start2 + mod(b - start2, a))
    n = stop2 - mod(stop2 - b, a)
    m:a:n
end

intersect(s::StepRange, r::InfStepRange) = intersect(r, s)
intersect(s::AbstractRange, r::InfStepRange) = intersect(StepRange(s), r)
intersect(s::InfStepRange, r::AbstractRange) = intersect(s, StepRange(r))

function union(ain::InfRanges, bin::InfRanges)
    a,b = promote(ain, bin)
    first(b) ∈ a && iszero(mod(step(b), step(a))) && return a
    first(a) ∈ b && iszero(mod(step(a), step(b))) && return b
    throw(ArgumentError("Cannot take union of $a and $b"))
end

promote_rule(a::Type{InfUnitRange{T1}}, b::Type{InfUnitRange{T2}}) where {T1,T2} =
    InfUnitRange{promote_type(T1,T2)}
InfUnitRange{T}(r::InfUnitRange{T}) where {T<:Real} = r
InfUnitRange{T}(r::UnitRange) where {T<:Real} = InfUnitRange{T}(r.start)

promote_rule(a::Type{OneToInf{T1}}, b::Type{OneToInf{T2}}) where {T1,T2} =
    OneToInf{promote_type(T1,T2)}
OneToInf{T}(r::OneToInf{T}) where {T<:Integer} = r
OneToInf{T}(r::OneToInf) where {T<:Integer} = OneToInf{T}()

promote_rule(::Type{InfStepRange{T1a,T1b}}, ::Type{InfStepRange{T2a,T2b}}) where {T1a,T1b,T2a,T2b} =
    InfStepRange{promote_type(T1a,T2a), promote_type(T1b,T2b)}
promote_rule(a::Type{InfStepRange{T1a,T1b}}, ::Type{UR}) where {T1a,T1b,UR<:AbstractInfUnitRange} =
    promote_rule(a, InfStepRange{eltype(UR), eltype(UR)})

InfStepRange{T1,T2}(r::InfStepRange{T1,T2}) where {T1,T2} = r

InfStepRange{T1,T2}(r::AbstractRange) where {T1,T2} =
    InfStepRange{T1,T2}(convert(T1, first(r)), convert(T2, step(r)))
InfStepRange(r::InfUnitRange{T}) where {T} =
    InfStepRange{T,T}(first(r), step(r))
(::Type{InfStepRange{T1,T2} where T1})(r::AbstractRange) where {T2} =
    InfStepRange{eltype(r),T2}(r)

## sorting ##

sum(r::InfRanges{<:Real}) = last(r)
mean(r::InfRanges{<:Real}) = last(r)
median(r::InfRanges{<:Real}) = last(r)

in(x::Union{Infinity,RealInfinity}, r::InfRanges) = false # never reach it...
in(x::Infinity, r::InfRanges{<:Integer}) = false # never reach it...
in(x::Real, r::InfRanges{<:Real}) = _in_range(x, r)
# This method needs to be defined separately since -(::T, ::T) can be implemented
# even if -(::T, ::Real) is not
in(x::T, r::InfRanges{T}) where {T} = _in_range(x, r)

in(x::Integer, r::AbstractInfUnitRange{<:Integer}) = (first(r) <= x)

in(x::Real, r::InfRanges{T}) where {T<:Integer} =
    isinteger(x) && !isempty(r) && ifelse(step(r) > 0, x ≥ first(r), x ≤ first(r)) &&
            (mod(convert(T,x),step(r))-mod(first(r),step(r)) == 0)

# Addition/subtraction of ranges
-(r1::OneToInf{T}, r2::OneToInf{V}) where {T,V} = Zeros{promote_type(T,V)}(∞)
-(r1::AbstractInfUnitRange, r2::AbstractInfUnitRange) = Fill(first(r1)-first(r2), ∞)

function sort!(a::InfStepRange)
    step(a) > 0 || throw(ArgumentError("Cannot sort $a"))
    a
end

### Lazy broadcasting

BroadcastStyle(::Type{<:InfRanges}) = LazyArrayStyle{1}()
BroadcastStyle(::Type{<:Adjoint{<:Any,<:InfRanges}}) = LazyArrayStyle{2}()
BroadcastStyle(::Type{<:Transpose{<:Any,<:InfRanges}}) = LazyArrayStyle{2}()
const InfSubVector = SubArray{<:Any,1,<:Any,<:Tuple{InfRanges}}
BroadcastStyle(::Type{<:InfSubVector}) = LazyArrayStyle{1}()
BroadcastStyle(::Type{<:SubArray{<:Any,1,<:Any,<:Tuple{InfSubVector}}}) = LazyArrayStyle{1}()

BroadcastStyle(::Type{<:Base.Slice{<:InfRanges}}) = LazyArrayStyle{1}()

const InfIndexRanges{T<:Integer} = Union{InfStepRange{T},AbstractInfUnitRange{T},Slice{OneToInf{T}}}

BroadcastStyle(::Type{<:SubArray{<:Any,1,<:Any,Tuple{<:InfIndexRanges}}})= LazyArrayStyle{1}()
BroadcastStyle(::Type{<:SubArray{<:Any,1,<:Any,<:Tuple{<:InfIndexRanges,<:Any}}})= LazyArrayStyle{1}()
BroadcastStyle(::Type{<:SubArray{<:Any,1,<:Any,<:Tuple{<:Any,<:InfIndexRanges}}})= LazyArrayStyle{1}()
BroadcastStyle(::Type{<:SubArray{<:Any,1,<:LazyMatrix,<:Tuple{<:InfIndexRanges,<:Any}}})= LazyArrayStyle{1}()
BroadcastStyle(::Type{<:SubArray{<:Any,1,<:LazyMatrix,<:Tuple{<:Any,<:InfIndexRanges}}})= LazyArrayStyle{1}()
BroadcastStyle(::Type{<:SubArray{<:Any,2,<:Any,<:Tuple{<:InfIndexRanges,<:InfIndexRanges}}})= LazyArrayStyle{2}()
BroadcastStyle(::Type{<:SubArray{<:Any,2,<:Any,<:Tuple{<:InfIndexRanges,<:Any}}})= LazyArrayStyle{2}()
BroadcastStyle(::Type{<:SubArray{<:Any,2,<:Any,<:Tuple{<:Any,<:InfIndexRanges}}})= LazyArrayStyle{2}()

broadcasted(::BroadcastStyle, f, r::Adjoint{<:Any,<:InfRanges}) = broadcast(f,parent(r))'
broadcasted(::BroadcastStyle, f, r::Transpose{<:Any,<:InfRanges}) = transpose(broadcast(f,parent(r)))
broadcasted(::BroadcastStyle, f, a::Number, r::Adjoint{<:Any,<:InfRanges}) = broadcast(f,a,parent(r))'
broadcasted(::BroadcastStyle, f, a::Number, r::Transpose{<:Any,<:InfRanges}) = transpose(broadcast(f,a,parent(r)))
broadcasted(::BroadcastStyle, f, r::Adjoint{<:Any,<:InfRanges}, a::Number) = broadcast(f,parent(r),a)'
broadcasted(::BroadcastStyle, f, r::Transpose{<:Any,<:InfRanges}, a::Number) = transpose(broadcast(f,parent(r),a))

broadcast(f, r::Adjoint{<:Any,<:InfRanges}) = broadcast(f,parent(r))'
broadcast(f, r::Transpose{<:Any,<:InfRanges}) = transpose(broadcast(f,parent(r)))
broadcast(f, a::Number, r::Adjoint{<:Any,<:InfRanges}) = broadcast(f,a,parent(r))'
broadcast(f, a::Number, r::Transpose{<:Any,<:InfRanges}) = transpose(broadcast(f,a,parent(r)))
broadcast(f, r::Adjoint{<:Any,<:InfRanges}, a::Number) = broadcast(f,parent(r),a)'
broadcast(f, r::Transpose{<:Any,<:InfRanges}, a::Number) = transpose(broadcast(f,parent(r),a))

# cumsum(r::InfRanges) = OneToInf() .* (first(r) .+ r) .÷ 2
cumsum(r::InfRanges) = RangeCumsum(r)
diff(r::InfRanges) = Fill(step(r),∞)
diff(r::OneToInf{T}) where T = Ones{T}(∞)
Base.@propagate_inbounds getindex(c::RangeCumsum, kr::OneToInf) = RangeCumsum(c.range[kr])
getindex(c::RangeCumsum{<:Any,<:OneToInf}, k::Integer) = k * (k+1) ÷ 2

###
# MemoryLayout
####

MemoryLayout(::Type{<:AbstractInfUnitRange}) = LazyLayout()
@inline MemoryLayout(::Type{<:InfStepRange}) = LazyLayout()

##
# findfirst
##

# from array.jl
function _step_findfirst(p, r::InfStepRange{T,S}) where {T,S}
    first(r) <= p.x || return nothing
    d = convert(S, p.x - first(r))
    iszero(d % step(r)) || return nothing
    return d ÷ step(r) + 1
end

for op in (:isequal, :(==))
    @eval begin
        findfirst(p::Fix2{typeof($op),T}, r::InfStepRange{T,S}) where {T,S} =
            _step_findfirst(p, r)

        findfirst(p::Fix2{typeof($op),T}, r::InfStepRange{T,<:Integer}) where {T<:Integer} =
            _step_findfirst(p, r)

        findfirst(p::Fix2{typeof($op),<:Integer}, r::InfStepRange{<:Integer,<:Integer}) =
            _step_findfirst(p, r)

        function findfirst(p::Fix2{typeof($op),<:Integer}, r::AbstractInfUnitRange{<:Integer})
            first(r) <= p.x || return nothing
            p.x - first(r) + 1
        end

        findfirst(p::Fix2{typeof($op),<:Number}, r::AbstractInfUnitRange{<:Integer}) =
            isinteger(p.x) ? findfirst($op(convert(Integer, p.x)), r) : nothing
        findfirst(p::Fix2{typeof($op),<:Number}, r::InfStepRange{<:Integer,<:Integer}) =
            isinteger(p.x) ? findfirst($op(convert(V, p.x)), r) : nothing
    end
end

FillArrays._range_convert(::Type{AbstractVector{T}}, r::InfRanges) where T = convert(AbstractVector{T}, r)

permutedims(D::Diagonal{<:Any,<:InfRanges}) = D
copy(D::Diagonal{<:Any,<:InfRanges}) = D
broadcasted(::LazyArrayStyle{2}, ::typeof(*), a::Number, D::Diagonal{<:Any,<:InfRanges}) = a*D
broadcasted(::LazyArrayStyle{2}, ::typeof(*), D::Diagonal{<:Any,<:InfRanges}, a::Number) = D*a


text:
#minor utility code that does belong anywhere else

#used in implementation of CSRGraph and CSRMatrix
function computeOffsets(rowPtrs::AbstractArray{<: Integer, 1}, numEnts::Integer)
    numOffsets = length(rowPtrs)
    @simd for i = 1:numOffsets
        @inbounds rowPtrs[i] = numEnts*(i-1)+1
    end
    rowPtrs
end
function computeOffsets(rowPtrs::AbstractArray{<: Integer, 1}, numEnts::Array{<: Integer, 1})
    numOffsets = length(rowPtrs)
    numCounts = length(numEnts)
    if numCounts >= numOffsets
        throw(InvalidArgumentError("length(numEnts) = $numCounts "
                * ">= length(rowPtrs) = $numOffsets"))
    end
    sum = 1
    for i = 1:numCounts
        @inbounds rowPtrs[i] = sum
        @inbounds sum += numEnts[i]
    end
    @inbounds rowPtrs[numCounts+1:numOffsets] .= sum
    sum-1
end

export InvalidArgumentError
"""
    InvalidArgumentError(msg)

The values passed as arguments are not valid.  Argument `msg`
is a descriptive error string.
"""
struct InvalidArgumentError <: Exception
    msg::AbstractString
end

export InvalidStateError
"""
    InvalidStateError(msg)

An object is not in a valid state for this method.  Argument `msg`
is a descriptive error string.
"""
struct InvalidStateError <: Exception
    msg::AbstractString
end


text:
using SimpleUpdate
using Test
using TensorOperations, TensorNetworkTensors, LinearAlgebra

@testset "SimpleUpdate.jl" begin
    @testset "Ising" begin
        u = tfisingpropagator(0,1)
        @test toarray(fuselegs(u,((1,2),(3,4)))[1]) ≈ I

        u = tfisingpropagator(1e-3,1)
        @tensor u2[1,2,3,4] := u[1,2,-1,-2] * u[-1,-2,3,4]
        @test u2 ≈ tfisingpropagator(2e-3,1)
        @test ishermitian(toarray(fuselegs(u,((1,2),(3,4)))[1]))
    end

    @testset "Weights" begin
        χ = 4
        a, b = DTensor(rand(χ,χ,χ,χ,2)), DTensor(rand(χ,χ,χ,χ,2))
        #rotate back and forth
        @tensor begin
            arot[1,2,3,4,5] := a[4,1,2,3,5]
            brot[1,2,3,4,5] := b[4,1,2,3,5]
        end
        @tensor begin
            a2[1,2,3,4,5] := arot[2,3,4,1,5]
            b2[1,2,3,4,5] := brot[2,3,4,1,5]
        end
        @test a2 ≈ a
        @test b2 ≈ b

        #[a]--[b] → [a]-c-[b]
        @tensor n1 = a[-1,lab,-2,-3,aa] * conj(a)[-1,lab2,-2,-3,aa] *
                     b[-4,-5,-6,lab,bb] * conj(b)[-4,-5,-6,lab2,bb]
        c, = extract_weight!(a,b)
        @tensor a[1,2,3,4,5] = a[1,-2,3,4,5] * c[-2,2]
        @tensor n2 = a[-1,lab,-2,-3,aa] * conj(a)[-1,lab2,-2,-3,aa] *
                     b[-4,-5,-6,lab,bb] * conj(b)[-4,-5,-6,lab2,bb]
        @test n1 ≈ n2

        #[b]--[a] → [b]-c-[a]
        b, a = DTensor(rand(χ,χ,χ,χ,2)), DTensor(rand(χ,χ,χ,χ,2))
        @tensor n1 = b[-1,lba,-2,-3,bb] * conj(b)[-1,lba2,-2,-3,bb] *
                     a[-4,-5,-6,lba,aa] * conj(a)[-4,-5,-6,lba2,aa]
        c, = extract_weight!(b,a)
        @tensor b[1,2,3,4,5] = b[1,-2,3,4,5] * c[-2,2]
        @tensor n2 = b[-1,lba,-2,-3,bb] * conj(b)[-1,lba2,-2,-3,bb] *
                     a[-4,-5,-6,lba,aa] * conj(a)[-4,-5,-6,lba2,aa]
        @test n1 ≈ n2

        #w/ rotation
        a, b = DTensor(rand(χ,χ,χ,χ,2)), DTensor(rand(χ,χ,χ,χ,2))
        @tensor n1 = a[-1,-2,lab,-3,aa] * conj(a)[-1,-2,lab2,-3,aa] *
                     b[lab,-5,-6,-4,bb] * conj(b)[lab2,-5,-6,-4,bb]
        @tensor begin
            arot[1,2,3,4,5] := a[4,1,2,3,5]
            brot[1,2,3,4,5] := b[4,1,2,3,5]
        end

        c, = extract_weight!(arot, brot)
        @tensor arot[1,2,3,4,5] = arot[1,-2,3,4,5] * c[-2,2]
        @tensor n2 = arot[-1,lab,-2,-3,aa] * conj(arot)[-1,lab2,-2,-3,aa] *
                     brot[-4,-5,-6,lab,bb] * conj(brot)[-4,-5,-6,lab2,bb]
        @test n1 ≈ n2
    end
end


text:
module Curiosity

using Reexport
using GVFHordes

# abstract type Learner end

@reexport using MinimalRLCore

range_check(v, min, max) = v >= min && v <= max

using Flux
import Flux.Optimise: update!

# Abstract type of FeatureProjector used in the learner utils
abstract type FeatureCreator end
# export ValueFeatureProjector, ActionValueFeatureProjector
include("utils/replay.jl")
include("utils/learners.jl")
include("utils/SRHorde.jl")

# Code utils
include("utils/dpi.jl")

export Auto
include("optimizers/Auto.jl")

export QLearner, LinearQLearner, VLearner, SRLearner, GPI, predict, predict_SF
include("learner.jl")

abstract type IntrinsicReward end
include("agent/intrinsic_rewards.jl")

export TabularRoundRobin, update!
include("updates/TabularRoundRobin.jl")

abstract type ExplorationStrategy end
export EpsilonGreedy, ϵGreedyDecay
include("agent/exploration.jl")

export Agent, AgentER, step!, get_demon_prediction, μ_π, proc_input
include("agent/agent.jl")
include("agent/agent_er.jl")

export TileCoder, create_features
include("./agent/tile_coder.jl")

export TabularTMaze, MountainCar, OneDTMaze, valid_state_mask
include("environments.jl")

# logger
export Logger, logger_start!, logger_step!, logger_episode_end!, LoggerKey, LoggerInitKey
include("logger/logger.jl")

#utils
export get_active_action_state_vector, ValueFeatureProjector
include("utils/SRCreation.jl")
include("utils/tmaze.jl")
include("utils/1d-tmaze.jl")
include("utils/2d-gridworld.jl")
include("utils/mountain_car.jl")
include("utils/experiment.jl")
include("utils/agent.jl")
include("utils/features.jl")
include("utils/baseline_learners.jl")

using GVFHordes
export monte_carlo_returns,monte_carlo_return
include("monte_carlo.jl")

end # module

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
module MaximumLikelihood

using Optim
using PDMats
using ValidatedNumerics
using StatsFuns
using AutoAligns
using ForwardDiff

export estimate_ML, estimate_ML_parameters

"""
Maximum likelihood estimator.

θ  ML estiamte
Σ  covariance matrix (estimated at the mode, using Fisher information)
varnames  variable names
optimization_result  contains convergence information from Optim.optimize
"""
immutable ML_estimator{T, ΣT <: AbstractPDMat, oT}
    θ::Vector{T}
    Σ::ΣT
    varnames::Vector{String}
    optimization_result::oT
end

"""
Default variable names, when not provided.
"""
default_varnames(n) = [string("θ", i) for i in 1:n]

"Calculate a default scaling for the log likelihood."
function calculate_log_ℓ_scale(log_ℓ_scale, log_ℓ, initial_θ)
    if log_ℓ_scale == :default
        Float64(1 / (1+abs(log_ℓ(initial_θ))))
    elseif log_ℓ_scale == :none
        1.0
    elseif isa(log_ℓ_scale, Real)
        convert(Float64, log_ℓ_scale)
    else
        error("Invalid log_ℓ_scale $(log_ℓ_scale)")
    end
end

"""
Internal function for estimating the maximum likelihood parameters of
the log likelihood function `log_ℓ`, with initial value
`initial_θ`. Return the result of `optimize`.

The log likelihood is scaled with `log_ℓ_scale`, for which `:default`
uses information estimates from the function at `initial_θ`, while
`:none` uses a unit scaling.

`method` and `options` are passed to `Optim.optimize`.
"""
function _estimate_ML(log_ℓ, initial_θ;
                     log_ℓ_scale = :default,
                     method = BFGS(),
                     options = Optim.Options(autodiff = true))
    γ = -calculate_log_ℓ_scale(log_ℓ_scale, log_ℓ, initial_θ)
    o = optimize(θ -> γ * log_ℓ(θ), initial_θ,
                 method, Optim.Options(autodiff = true))
    if !Optim.converged(o)
        error("Maximum likelihood did not converge. Check concavity and initial value.")
    end
    o
end

"""
Maximum likelihood estimator. Returns just the parameters as a vector.
"""
function estimate_ML_parameters(log_ℓ, initial_θ,
                                log_ℓ_scale = :default,
                                method = BFGS(),
                                options = Optim.Options(autodiff = true))
    Optim.minimizer(_estimate_ML(log_ℓ, initial_θ,
                                 log_ℓ_scale = log_ℓ_scale,
                                 method = method,
                                 options = options))
end

"""
Estimate parameters using maximum likelihood, with the log likelihood
function `log_ℓ` and initial guess `initial_θ`.

Return an `ML_estimator`, with `varnames` optimally
specified. `method` and `options` are passed to
`Optim.optimize`. `log_ℓ` needs to work with ForwardDiff.
"""
function estimate_ML(log_ℓ, initial_θ;
                     log_ℓ_scale = :default,
                     θ_scale = :default,
                     method = BFGS(),
                     options = Optim.Options(autodiff = true),
                     varnames = default_varnames(length(initial_θ)))
    o = _estimate_ML(log_ℓ, initial_θ;
                     log_ℓ_scale = log_ℓ_scale,
                     method = method,
                     options = options)
    θ = Optim.minimizer(o)
    Σ = inv(PDMat(Symmetric(-ForwardDiff.hessian(log_ℓ, θ))))
    ML_estimator(θ, Σ, varnames, o)
end

"""
Univariate confidence intervals.
"""
function confidence_intervals(ml::ML_estimator, p=0.025)
    @assert p < 0.5
    s = norminvcdf(1-p)
    [θ ± σ*s for (θ, σ) in zip(ml.θ, diag(ml.Σ))]
end

function Base.show(io::IO, ml::ML_estimator; p = 0.025)
    @assert p < 0.5
    aa = AutoAlign(align = Dict(1 => left, :default => right))
    p_string(p) = string(round(100*p, 1), "%")
    v_string(v) = @sprintf("%.5g", v) # significant digits hardcoded
    println(aa, "name", "  ", "est", "  ", p_string(p), "  ", p_string(1-p))
    for (varname, θ, ci) in zip(ml.varnames, ml.θ, confidence_intervals(ml, 0.025))
        println(aa, varname, "", v_string(θ), "", v_string(ci.lo), "", v_string(ci.hi))
    end
    print(io, aa)
end

end # module