text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

#############
# test_gthr #
#############

function test_gthr!(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    d_y = CUSPARSE.gthr!(d_x,d_y,'O')
    h_x = to_host(d_x)
    if VERSION >= v"0.5.0-dev+742"
        @test h_x ≈ SparseVector(m,x.nzind,y[x.nzind])
    else
        @test h_x ≈ sparsevec(x.rowval,y[x.rowval],m)
    end
end

function test_gthr(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    d_z = CUSPARSE.gthr(d_x,d_y,'O')
    h_z = to_host(d_z)
    if VERSION >= v"0.5.0-dev+742"
        @test h_z ≈ SparseVector(m,x.nzind,y[x.nzind])
    else
        @test h_z ≈ sparsevec(x.rowval,y[x.rowval],m)
    end
end

##############
# test_gthrz #
##############

function test_gthrz!(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    d_x,d_y = CUSPARSE.gthrz!(d_x,d_y,'O')
    h_x = to_host(d_x)
    h_y = to_host(d_y)
    if VERSION >= v"0.5.0-dev+742"
        @test h_x ≈ SparseVector(m,x.nzind,y[x.nzind])
        y[x.nzind] = zero(elty)
        @test h_y ≈ y
    else
        @test h_x ≈ sparsevec(x.rowval,y[x.rowval],m)
        y[x.rowval] = zero(elty)
        @test h_y ≈ y
    end
end

function test_gthrz(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    d_z,d_w = CUSPARSE.gthrz(d_x,d_y,'O')
    h_w = to_host(d_w)
    h_z = to_host(d_z)
    if VERSION >= v"0.5.0-dev+742"
        @test h_z ≈ SparseVector(m,x.nzind,y[x.nzind])
        y[x.nzind] = zero(elty)
        @test h_w ≈ y
    else
        @test h_z ≈ sparsevec(x.rowval,y[x.rowval],m)
        y[x.rowval] = zero(elty)
        @test h_w ≈ y
    end
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_gthr!(elty)
    println("gthr! took ", toq(), " for ", elty)
    tic()
    test_gthr(elty)
    println("gthr took ", toq(), " for ", elty)
    tic()
    test_gthrz!(elty)
    println("gthrz! took ", toq(), " for ", elty)
    tic()
    test_gthrz(elty)
    println("gthrz took ", toq(), " for ", elty)
end


text:
module Strided

import Base: parent, size, strides, tail, setindex
using Base: @propagate_inbounds, RangeIndex, Dims
const SliceIndex = Union{RangeIndex,Colon}

using LinearAlgebra

using TupleTools
using TupleTools: StaticLength

export StridedView, @strided, @unsafe_strided, sreshape, sview

# function __init__()
#     LinearAlgebra.BLAS.set_num_threads(1)
#     Threads.nthreads() == 1 && warn("Strided disables BLAS multithreading, enable Julia threading (`export JULIA_NUM_THREADS = N`) to benefit from multithreaded matrix multiplication and more")
# end

const _use_threads = Ref(true)
use_threads() = _use_threads[]

"""
    disable_threads()

Disable the use of threading in Strided.
Also see [`enable_threads()`](@ref)
"""
function disable_threads()
    _use_threads[] = false
    return
end

"""
    enable_threads()

(Re)-enable threading in Strided.
"""
function enable_threads()
    _use_threads[] = true
    return
end

# used to factor the number of threads
function simpleprimefactorization(n::Int)
    factors = Vector{Int}()
    k = 2
    while k <= n
        d, r = divrem(n, k)
        if r == 0
            push!(factors, k)
            n = d
        else
            k += 1
        end
    end
    return factors
end

const factors = Vector{Int}(undef, 0)

function __init__()
    f = reverse!(simpleprimefactorization(Threads.nthreads()))
    resize!(factors, length(f))
    copyto!(factors, f)
end

include("abstractstridedview.jl")
include("stridedview.jl")
include("unsafestridedview.jl")
include("mapreduce.jl")
include("broadcast.jl")
include("macros.jl")

# function _precompile_()
#     ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
#     for T in (Float64, ComplexF64)
#         USMatrix = UnsafeStridedView{T, 2, T, typeof(Base.identity)}
#         SMatrix = StridedView{T, 2, Vector{T}, typeof(Base.identity)}
#
#         precompile(Tuple{typeof(mul!), USMatrix, USMatrix, USMatrix})
#         precompile(Tuple{typeof(mul!), USMatrix, USMatrix, USMatrix, T, T})
#
#         precompile(Tuple{typeof(mul!), SMatrix, SMatrix, SMatrix})
#         precompile(Tuple{typeof(mul!), SMatrix, SMatrix, SMatrix, T, T})
#
#         if T <: Complex
#             CUSMatrix = UnsafeStridedView{T, 2, T, typeof(Base.conj)}
#             CSMatrix = StridedView{T, 2, Vector{T}, typeof(Base.conj)}
#
#             precompile(Tuple{typeof(mul!), USMatrix, USMatrix, CUSMatrix})
#             precompile(Tuple{typeof(mul!), USMatrix, CUSMatrix, USMatrix})
#             precompile(Tuple{typeof(mul!), USMatrix, CUSMatrix, CUSMatrix})
#
#             precompile(Tuple{typeof(mul!), SMatrix, SMatrix, CSMatrix})
#             precompile(Tuple{typeof(mul!), SMatrix, CSMatrix, SMatrix})
#             precompile(Tuple{typeof(mul!), SMatrix, CSMatrix, CSMatrix})
#         end
#
#         for N = 1:4
#             USArray = UnsafeStridedView{T, N, T, typeof(Base.identity)}
#             SArray = StridedView{T, N, Vector{T}, typeof(Base.identity)}
#
#             for ASArray in (USArray, SArray)
#                 precompile(Tuple{typeof(copyto!), ASArray, ASArray})
#                 precompile(Tuple{typeof(mul!), ASArray, T, ASArray})
#                 precompile(Tuple{typeof(mul!), ASArray, ASArray, T})
#                 precompile(Tuple{typeof(rmul!), ASArray, T})
#                 precompile(Tuple{typeof(lmul!), T, ASArray})
#                 precompile(Tuple{typeof(axpy!), T, ASArray, ASArray})
#                 precompile(Tuple{typeof(axpby!), T, ASArray, T, ASArray})
#             end
#
#             if T <: Complex
#                 CUSArray = UnsafeStridedView{T, N, T, typeof(Base.conj)}
#                 CSArray = StridedView{T, N, Vector{T}, typeof(Base.conj)}
#
#                 for (ASArray, CASArray) in ((USArray, CUSArray), (SArray, CSArray))
#                     precompile(Tuple{typeof(copyto!), ASArray, CASArray})
#                     precompile(Tuple{typeof(mul!), ASArray, T, CASArray})
#                     precompile(Tuple{typeof(mul!), ASArray, CASArray, T})
#                     precompile(Tuple{typeof(rmul!), CASArray, T})
#                     precompile(Tuple{typeof(lmul!), T, CASArray})
#                     precompile(Tuple{typeof(axpy!), T, CASArray, ASArray})
#                     precompile(Tuple{typeof(axpby!), T, CASArray, T, ASArray})
#                 end
#             end
#         end
#     end
# end
# _precompile_()

end

ext: jl


text:
###################################################################################################
##### Input #####
###################################################################################################
"""
  loadFile(filename::String)
Load a file of any format libAntimony knows about (potentially Antimony, SBML, or CellML).
"""
function loadFile(filename::String)
  status = ccall(dlsym(antlib, :loadFile), cdecl, Int64, (Ptr{UInt8},), filename)
  if status == -1
    error(getLastError())
  end
end

"""
  loadString(model::String)
Load a string of any format libAntimony knows about (potentially Antimony, SBML, or CellML).
"""
function loadString(model::String)
  status = ccall(dlsym(antlib, :loadString), cdecl, Int64, (Ptr{UInt8},), model)
  if status == -1
    error(getLastError())
  end
end

"""
  loadAntimonyFile(filename::String)
Loads a file and parses it as an Antimony file.
"""
function loadAntimonyFile(filename::String)
  status = ccall(dlsym(antlib, :loadAntimonyFile), cdecl, Int64, (Ptr{UInt8},), filename)
  if status == -1
    error(getLastError())
  end
end

"""
  loadAntimonyString(antModel::String)
Loads a string and parses it as an Antimony set of modules.
"""
function loadAntimonyString(antModel::String)
  status = ccall(dlsym(antlib, :loadAntimonyString), cdecl, Int64, (Ptr{UInt8},), antModel)
  if status == -1
    error(getLastError())
  end
  return status
end

"""
  loadSBMLFile(filename::String)
Load a file known to be SBML.
"""
function loadSBMLFile(filename::String)
  status = ccall(dlsym(antlib, :loadSBMLFile), cdecl, Int64, (Ptr{UInt8},), filename)
  if status == -1
    error(getLastError())
  end
  return status
end

"""
  loadSBMLString(model::String)
Load a string known to be SBML.
"""
function loadSBMLString(model::String)
  status = ccall(dlsym(antlib, :loadSBMLString), cdecl, Int64, (Ptr{UInt8},), model)
  if status == -1
    error(getLastError())
  end
end

"""
  loadSBMLStringWithLocation(model::String, location::String)
Load a string known to be SBML with its file location.
"""
function loadSBMLStringWithLocation(model::String, location::String)
  status = ccall(dlsym(antlib, :loadSBMLStringWithLocation), cdecl, Int64, (Ptr{UInt8}, Ptr{UInt8}), model, location)
  if status == -1
    error(getLastError())
  end
end

"""
  loadCellMLFile(filename::String)
Load a file known to be CellML.
"""
function loadCellMLFile(filename::String)
  status = ccall(dlsym(antlib, :loadCellMLFile), cdecl, Int64, (Ptr{UInt8},), filename)
  if status == -1
    error(getLastError())
  end
end

"""
  loadCellMLString(model::String)
Load a string known to be CellML.
"""
function loadCellMLString(model::String)
  status = ccall(dlsym(antlib, :loadCellMLString), cdecl, Int64, (Ptr{UInt8},), model)
  if status == -1
    error(getLastError())
  end
end

"""
  getNumFiles()
Returns the number of files loaded into memory so far.
"""
function getNumFiles()
  return ccall(dlsym(antlib, :getNumFiles), cdecl, UInt64, ())
end

"""
  revertTo(Index)
Change the 'active' set of modules to the ones from the given index (as received from 'load<file/string>').
"""
function revertTo(Index)
  status = ccall(dlsym(antlib, :revertTo), cdecl, Bool, (Int64,), Index)
  if status == false
    error(getLastError())
  end
end

"""
  checkModule(moduleName::String)
Returns 'true' if the submitted module name exists in the current active set, 'false' if not.
"""
function checkModule(moduleName::String)
  status = ccall(dlsym(antlib, :checkModule), cdecl, Bool, (Ptr{UInt8},), moduleName)
  if status == false
    error(getLastError())
  end
end

"""
  clearPreviousLoads()
Clears memory of all files loaded.
"""
function clearPreviousLoads()
  ccall(dlsym(antlib, :clearPreviousLoads), cdecl, Cvoid, ())
end

"""
  addDirectory(directory::String)
Add a directory in which imported files may be found,
and in which to look for a '.antimony' file
(which contains rules about where to look locally for imported antimony and sbml files).
"""
function addDirectory(directory::String)
  status = ccall(dlsym(antlib, :addDirectory), cdecl, Int64, (Ptr{UInt8},), directory)
  if status == -1
    error(getLastError())
  end
end

"""
  clearDirectories()
Clears the list of directories added with the 'addDirectory' function.
"""
function clearDirectories()
  status = ccall(dlsym(antlib, :clearDirectories), cdecl, Int64, ())
  if status == -1
    error(getLastError())
  end
end

###################################################################################################
##### Output #####
###################################################################################################
"""
  writeAntimonyFile(filename::String, moduleName::String)
Writes out an antimony-formatted file containing the given module.
"""
function writeAntimonyFile(filename::String, moduleName::String)
  status = ccall(dlsym(antlib, :writeAntimonyFile), cdecl, Int64, (Ptr{UInt8}, Ptr{UInt8}), filename, moduleName)
  if status == 0
    error(getLastError())
  end
  return status
end

"""
  getAntimonyString(moduleName::String)
Returns the same output as writeAntimonyFile, but to a char* array instead of to a file.
"""
function getAntimonyString(moduleName::String)
  return unsafe_string(ccall(dlsym(antlib, :getAntimonyString), cdecl, Ptr{UInt8}, (Ptr{UInt8},), moduleName))
end

"""
  writeSBMLFile(filename::String, moduleName::String)
Writes out a SBML-formatted XML file to the file indicated.
"""
function writeSBMLFile(filename::String, moduleName::String)
  status = ccall(dlsym(antlib, :writeSBMLFile), cdecl, Int64, (Ptr{UInt8}, Ptr{UInt8}), filename, moduleName)
  if status == 0
    error(getLastError())
  end
  return status
end

"""
  getSBMLString(moduleName::String)
Returns the same output as writeSBMLFile, but to a char* array instead of to a file.
"""
function getSBMLString(moduleName::String)
  return unsafe_string(ccall(dlsym(antlib, :getSBMLString), cdecl, Ptr{UInt8}, (Ptr{UInt8},), moduleName))
end

"""
  writeCompSBMLFile(filename::String, moduleName::String)
Writes out a SBML-formatted XML file to the file indicated,
using the 'Hierarchichal Model Composition' package
"""
function writeCompSBMLFile(filename::String, moduleName::String)
  status = ccall(dlsym(antlib, :writeCompSBMLFile), cdecl, Int64, (Ptr{UInt8}, Ptr{UInt8}), filename, moduleName)
  if status == 0
    error(getLastError())
  end
  return status
end

"""
  getCompSBMLString(moduleName::String)
Returns the same output as writeSBMLFile, but to a char* array instead of to a file,
using the 'Hierarchichal Model Composition' package.
"""
function getCompSBMLString(moduleName::String)
  return unsafe_string(ccall(dlsym(antlib, :getCompSBMLString), cdecl, Ptr{UInt8}, (Ptr{UInt8},), moduleName))
end

"""
  writeCellMLFile(filename::String, moduleName::String)
Writes out a CellML-formatted XML file to the file indicated,
retaining the same Antimony hierarchy using the CellML 'component' hieararchy.
"""
function writeCellMLFile(filename::String, moduleName::String)
  status = ccall(dlsym(antlib, :writeCellMLFile), cdecl, Int64, (Ptr{UInt8}, Ptr{UInt8}), filename, moduleName)
  if status == 0
    error(getLastError())
  end
  return status
end

"""
  getCellMLString(moduleName::String)
Writes out a CellML-formatted XML file to the file indicated,
retaining the same Antimony hierarchy using the CellML 'component' hieararchy.
"""
function getCellMLString(moduleName::String)
  return unsafe_string(ccall(dlsym(antlib, :getCellMLString), cdecl, Ptr{UInt8}, (Ptr{UInt8},), moduleName))
end

"""
  printAllDataFor(moduleName::String)
An example function that will print to stdout all the information in the given module.
"""
function printAllDataFor(moduleName::String)
  status = ccall(dlsym(antlib, :printAllDataFor), cdecl, Cvoid, (Ptr{UInt8},), moduleName)
  if status == 0
    error(getLastError())
  end
  return status
end

###################################################################################################
##### Errors and Warnings #####
###################################################################################################

"""
  antimony_getLastError()
When any function returns an error condition, a longer description of the problem is
stored in memory, and is obtainable with this function.
"""
function getLastError()
  return unsafe_string(ccall(dlsym(antlib, :getLastError), cdecl, Ptr{UInt8}, ()))
end

"""
  getWarnings()
When translating some other format to Antimony, elements that are unable to be
translated are saved as warnings, retrievable with this function
(returns NULL if no warnings present).
"""
function getWarnings()
  return unsafe_string(ccall(dlsym(antlib, :getWarnings), cdecl, Ptr{UInt8}, ()))
end

"""
  getSBMLInfoMessages(moduleName::String)
Returns the 'info' messages from libSBML.
"""
function getSBMLInfoMessages(moduleName::String)
  return unsafe_string(ccall(dlsym(antlib, :getSBMLInfoMessages), cdecl, Ptr{UInt8}, (Ptr{UInt8},), moduleName))
end

"""
  getSBMLWarnings(moduleName::String)
Returns the 'warning' messages from libSBML.
"""
function getSBMLWarnings(moduleName::String)
  return unsafe_string(ccall(dlsym(antlib, :getSBMLWarnings), cdecl, Ptr{UInt8}, (Ptr{UInt8},), moduleName))
end

###################################################################################################
##### Modules #####
###################################################################################################
"""
  getNumModules()
Returns the number of modules in the current active set
(the last file successfully loaded, or whichever file was returned to with 'revertTo').
"""
function getNumModules()
  return ccall(dlsym(antlib, :getNumModules), cdecl, Int64, ())
end

"""
  getModuleNames()
Returns an array of all the current module names.
"""
function getModuleNames()
  numModules = getNumModules()
  moduleNames = Array{String}(undef, numModules)
  for i = 1:numModules
    moduleNames[i] = unsafe_string(unsafe_load(ccall(dlsym(antlib, :getModuleNames), cdecl, Ptr{Ptr{UInt8}}, (UInt64,), i)))
  end
  return moduleNames
end

"""
  getNthModuleName(n::Int64)
Returns the nth module name.
"""
function getNthModuleName(n::Int64)
  return unsafe_string(ccall(dlsym(antlib, :getNthModuleName), cdecl, Ptr{UInt8}, (Int64,), n))
end

"""
  getMainModuleName()
Returns the 'main' module name.
"""
function getMainModuleName()
  return unsafe_string(ccall(dlsym(antlib, :getMainModuleName), cdecl, Ptr{UInt8}, ()))
end

###################################################################################################
##### Module Interface #####
###################################################################################################

###################################################################################################
##### Replacements #####
###################################################################################################

###################################################################################################
##### Symbols and symbol information #####
###################################################################################################

###################################################################################################
##### Reaction #####
###################################################################################################

###################################################################################################
##### Interactions #####
###################################################################################################

###################################################################################################
##### Stoichiometry matrix information #####
###################################################################################################

###################################################################################################
##### Events #####
###################################################################################################

###################################################################################################
##### DNA strands #####
###################################################################################################

###################################################################################################
##### Memory management #####
###################################################################################################

"""
  freeAll()
Frees all pointers handed to you by libAntimony.
"""
function freeAll()
  ccall(dlsym(antlib, :freeAll), cdecl, Cvoid, ())
end

###################################################################################################
##### Defaults #####
###################################################################################################


text:
struct SnpBitMatrix{T} <: AbstractMatrix{T}
    B1::BitMatrix
    B2::BitMatrix
    model::Union{Val{1}, Val{2}, Val{3}}
    center::Bool
    scale::Bool
    μ::Vector{T}
    σinv::Vector{T}
    storagev1::Vector{T}
    storagev2::Vector{T}
end

"""
    SnpBitMatrix{T}(s; model=ADDITIVE_MODEL, center=false, scale=false)

Store an `AbstractSnpArray` `s` as two `BitMatrix`es to simplify linear algebra.

# Arguments
- s: an `AbstractSnpArray`.
- model: one of `ADDITIVE_MODEL`(default), `DOMINANT_MODEL`, `RECESSIVE_MODEL`.
- center: whether to center (default: false).
- scale: whether to scale to standard deviation 1 (default: false).
"""
function SnpBitMatrix{T}(
    s::AbstractSnpArray;
    model = ADDITIVE_MODEL,
    center::Bool = false,
    scale::Bool = false) where T <: AbstractFloat
    if model == ADDITIVE_MODEL
        B1 = s .≥ 0x02
        B2 = s .≥ 0x03
        if center || scale
            μ = Vector{T}(undef, size(s, 2))
            μ[:] = mean(s; dims=1, model=ADDITIVE_MODEL)
        else
            μ = T[]
        end
        if scale
            σinv = Vector{T}(undef, size(s, 2))
            @inbounds @simd for j in 1:size(s, 2)
                σinv[j] = sqrt(μ[j] * (1 - μ[j] / 2))
                σinv[j] = σinv[j] > 0 ? inv(σinv[j]) : one(T)
            end
        else
            σinv = T[]
        end
    elseif model == DOMINANT_MODEL
        B1 = s .≥ 0x02
        B2 = falses(0, 0)
        if center || scale
            μ = Vector{T}(undef, size(s, 2))
            μ[:] = mean(s; dims=1, model=DOMINANT_MODEL)
        else
            μ = T[]
        end
        if scale
            σinv = Vector{T}(undef, size(s, 2))
            @inbounds @simd for j in 1:size(s, 2)
                σinv[j] = sqrt(μ[j] * (1 - μ[j]))
                σinv[j] = σinv[j] > 0 ? inv(σinv[j]) : one(T)
            end
        else
            σinv = T[]
        end
    elseif model == RECESSIVE_MODEL
        B1 = s .== 0x03
        B2 = falses(0, 0)
        if center || scale
            μ = Vector{T}(undef, size(s, 2))
            μ[:] = mean(s; dims=1, model=RECESSIVE_MODEL)
        else
            μ = T[]
        end
        if scale
            σinv = Vector{T}(undef, size(s, 2))
            @inbounds @simd for j in 1:size(s, 2)
                σinv[j] = sqrt(μ[j] * (1 - μ[j]))
                σinv[j] = σinv[j] > 0 ? inv(σinv[j]) : one(T)
            end
        else
            σinv = T[]
        end
    else
        throw(ArgumentError("unrecognized model $model"))
    end
    storagev1 = Vector{T}(undef, size(s, 1))
    storagev2 = Vector{T}(undef, size(s, 2))
    SnpBitMatrix{T}(B1, B2, model, center, scale, μ, σinv, storagev1, storagev2)
end

Base.size(bm::SnpBitMatrix) = size(bm.B1)
Base.size(bm::SnpBitMatrix, k::Integer) = size(bm.B1, k)

eltype(bm::SnpBitMatrix) = eltype(bm.μ)
issymmetric(bm::SnpBitMatrix) = issymmetric(bm.B2) && issymmetric(bm.B1)

function _gemv_avx_sized!(c, A, b, rows, cols)
    @avx for j in 1:cols, i in 1:rows
        c[i] += A[i, j] * b[j]
    end
end

function _gemv_tile!(c, A, b)
    vstep = 512
    hstep = 512
    vstep_log2 = 9
    hstep_log2 = 9
    M, N = size(A)
    Miter = M >>> vstep_log2 # fast div(M, 512)
    Mrem = M & (vstep-1) # fast rem(M, 512)
    Niter = N >>> hstep_log2
    Nrem = N & (hstep-1)
    GC.@preserve c A b for n in 0:Niter-1
        for m in 0:Miter-1
            _gemv_avx_sized!(
                gesp(stridedpointer(c), (vstep*m,)),
                gesp(stridedpointer(A), (vstep*m, hstep*n)),
                gesp(stridedpointer(b), (hstep*n,)),
                vstep,
                hstep
            )
        end
        if Mrem != 0
            _gemv_avx_sized!(
                gesp(stridedpointer(c), (vstep*Miter,)),
                gesp(stridedpointer(A), (vstep*Miter, hstep*n)),
                gesp(stridedpointer(b), (hstep*n,)),
                Mrem, hstep
            )
        end
    end
    if Nrem != 0
        _gemv_avx_sized!(
            gesp(stridedpointer(c), (0,)),
            gesp(stridedpointer(A), (0, hstep*Niter)),
            gesp(stridedpointer(b), (hstep*Niter,)),
            length(c),
            Nrem
        )
    end
end

function mul!(c::Vector{T}, A::BitMatrix, b::Vector{T}) where T
    fill!(c, zero(eltype(c)))
    _gemv_tile!(c, A, b)
end

function mul!(c::Vector{T}, A::Union{Transpose{Bool,BitArray{2}}, Adjoint{Bool, BitArray{2}}}, b::Vector{T}) where T
    tA = transpose(A)
    fill!(c, zero(eltype(c)))
    @avx for i in 1:size(tA, 2)
        for j in 1:size(tA, 1)
            c[i] += tA[j, i] * b[j]
        end
    end
end

"""
    LinearAlgebra.mul!(out, s::SnpBitMatrix, v)

In-place matrix-vector multiplication.
"""
function mul!(
    out::AbstractVector{T},
    s::SnpBitMatrix{T},
    v::AbstractVector{T}) where T <: AbstractFloat
    @assert length(out) == size(s, 1) && length(v) == size(s, 2)
    if s.scale
        s.storagev2 .= s.σinv .* v
        w = s.storagev2
    else
        w = v
    end
    if s.model == ADDITIVE_MODEL
        mul!(out, s.B1, w)
        mul!(s.storagev1, s.B2, w)
        out .+= s.storagev1
    else
        mul!(out, s.B1, w)
    end
    if s.center
        return out .-= dot(s.μ, w)
    else
        return out
    end
end

"""
    LinearAlgebra.mul!(out, s::Union{Transpose{T, SnpBitMatrix{T}}, Adjoint{T, SnpBitMatrix{T}}}, v)

In-place matrix-vector multiplication, with transposed `SnpBitMatrix`.
"""
function mul!(
    out::AbstractVector{T},
    st::Union{Transpose{T, SnpBitMatrix{T}}, Adjoint{T, SnpBitMatrix{T}}},
    v::AbstractVector{T}) where T <: AbstractFloat
    s = st.parent
    @assert length(out) == size(st, 1) && length(v) == size(st, 2)
    if s.model == ADDITIVE_MODEL
        mul!(out, transpose(s.B1), v)
        mul!(s.storagev2, transpose(s.B2), v)
        out .+= s.storagev2
    else
        mul!(out, transpose(s.B1), v)
    end
    if s.center
        out .-= sum(v) .* s.μ
    end
    if s.scale
        return out .*= s.σinv
    else
        return out
    end
end


text:
macro callhook4(hook, r, req, id)
    quote
        for f! in $r.hooks[$hook]
            result = f!($r, $(esc(req)), $id)
            if isa(result, Dict{Symbol, Any})
                $(esc(req)) = result
            elseif isa(result, Response)
                return result
            end
        end
    end
end

macro callhook5(hook, r, req, id, res)
    quote
        for f! in $r.hooks[$hook]
            result = f!($r, $req, $id, $(esc(res)))
            if isa(result, Response)
                $(esc(res)) = result
            end
        end
        $(esc(res))
    end
end

(r::Resource)(req::Request) = req |> parserequest |> r

function (r::Resource)(req::Dict{Symbol, Any}, id::AbstractString="/")
    @callhook4(:onroute, r, req, id)
    res = if isempty(req[:path]) # leaf node
        @callhook4(:onhandle, r, req, id)
        raw = callmethod(r, req, id)
        @callhook5(:onresponse, r, req, id, raw)
    else
        req[Symbol(r.name * "id")] = id
        route(r.subresources, req, shift!(req[:path]))
    end
    @callhook5(:onreturn, r, req, id, res)
end

function route(candidates::Vector{Resource}, req, id)
    ismatch(p::AbstractString) = p == "*" || p == id
    ismatch(p::Function) = p(id)
    ismatch(p::Regex) = Base.ismatch(p, id)

    i = findfirst(x->ismatch(x.route), candidates)
    if i == 0
        Response(404)
    else
        candidates[i](req, id)
    end
end

function parserequest(req::Request)
    uri = URI(req.resource)
    Dict{Symbol, Any}(
        :method  => req.method |> Symbol,
        :headers => req.headers,
        :body    => req.data,
        :path    => uri.path   |> splitpath,
        :query   => uri.query  |> parsequerystring # this should be moved out of Rest.jl. Maybe add a mixin system?
    )
end

splitpath(p::AbstractString) = split(p, '/', keep=false)

function callmethod(r::Resource, req::Dict{Symbol, Any}, id::AbstractString)
    let m = r.methods, v = req[:method],
        h = [r.hooks[v]; (cb, r, req, id) -> haskey(m, v) ? m[v](req, id) : Response(405)]
        callone(i) = (req::Dict{Symbol, Any}, id::AbstractString) ->
            h[i](callone(i+1), r, req, id)
        callone(1)(req, id)
    end
end


text:
# Join two leaf nodes into one.
function merge!(left::LeafNode{K, V, B}, right::LeafNode{K, V, B}) where {K, V, B}
    leftlen, rightlen = length(left), length(right)
    @assert leftlen + rightlen <= B
    resize!(left, leftlen + rightlen)
    copyto!(left.entries, leftlen+1, right.entries, 1, length(right))
    copyto!(left.keys, leftlen+1, right.keys, 1, length(right))
    left.maxend = nodemaxend(left)
    return
end

# Join two internal nodes into one
function merge!(left::InternalNode{K, V, B}, right::InternalNode{K, V, B}) where {K, V, B}
    leftlen, rightlen = length(left), length(right)
    @assert length(left) + length(right) <= B
    resize!(left.keys, leftlen + rightlen - 1)
    resize!(left.children, leftlen + rightlen)
    resize!(left.maxends, leftlen + rightlen)
    copyto!(left.children, leftlen+1, right.children, 1, length(right.children))
    copyto!(left.maxends, leftlen+1, right.maxends, 1, length(right.maxends))
    for i in leftlen+1:length(left.children)
        left.children[i].parent = left
    end
    left.keys[leftlen] = minkey(left.children[leftlen+1])
    copyto!(left.keys, leftlen+1, right.keys, 1, length(right.keys))
    left.maxend = nodemaxend(left)
    return
end


text:
#=
This module contains functions that allow operations in extension fields
GF(2^m). It also contains a dictionary with primitive polynomials for
small extension fields.

Author: Ivan A. Moreno Soto
Last updated: 2019/July/02
=#

module ExtensionField

export primitivePolynomials,
    add,
    multiply,
    divide,
    modulo,
    power,
    inverse,
    squareroot,
    polygcd

primitivePolynomials = Dict(
    2 => UInt32(7),
    3 => UInt32(11),
    4 => UInt32(19),
    5 => UInt32(37),
    6 => UInt32(103),
    7 => UInt32(131),
    8 => UInt32(487),
    9 => UInt32(901),
    10 => UInt32(1033),
    11 => UInt32(2053),
    12 => UInt32(4179),
    13 => UInt32(8219),
    14 => UInt32(16707),
    15 => UInt32(32785),
)

"""
    add(a::UInt32, b::UInt32)

XORs two polynomials in GF(2).
"""
function add(a::UInt32, b::UInt32)
    return a ⊻ b
end

"""
    multiply(a::UInt32, b::UInt32)

Multiplies two polynomials in GF(2). It doesn't reduce modulo
an irreducible polynomial.
"""
function multiply(a::UInt32, b::UInt32)
    partial = Array{UInt32, 1}()
    for (i, bit) in enumerate(digits(b, base=2))
        if bit == 1
            push!(partial, a << (i - 1))
        end
    end

    c = UInt32(0)
    for i in partial
        c ⊻= i
    end

    return c
end

"""
    division(a::UInt32, b::UInt32)

Computes a/b. Returns the quotient and the remainder.
"""
function division(a::UInt32, b::UInt32)
    if a == b
        return UInt32(1), UInt32(0)
    elseif length(digits(a, base=2)) < length(digits(b, base=2)) #a < b
        return UInt32(0), a
    end

    quotient = UInt32(0)
    remainder = a
    for i in (length(digits(a, base=2)) - length(digits(b, base=2))):-1:0
        if length(digits(remainder, base=2)) >= length(digits(b << i, base=2))
            remainder ⊻= b << i
            quotient ⊻= UInt32(2^i)
        end
    end

    return quotient, remainder
end

"""
    divide(a::UInt32, b::UInt32)

Computes a/b but only returns the quotient.
"""
function divide(a::UInt32, b::UInt32)
    quotient, _ = division(a, b)
    return quotient
end

"""
    modulo(a::UInt32, b::UInt32)

Computes a/b but only returns the remainder.
"""
function modulo(a::UInt32, b::UInt32)
    _, remainder = division(a, b)
    return remainder
end

"""
    power(x::UInt32, n::Int, g::UInt32)

Computes and returns x^n modulo g in GF(2^m).
"""
function power(x::UInt32, n::Int, g::UInt32)
    if n == 0
        return UInt32(1)
    end

    y = UInt32(1)

    while n > 1
        if n % 2 == 0
            n = n / 2
        else
            y = modulo(multiply(x, y), g)
            n = (n - 1)/2
        end

        x = modulo(multiply(x, x), g)
    end

    return modulo(multiply(x, y), g)
end

"""
    inverse(a::UInt32, g::UInt32)

Computes and returns the inverse of a mod g(x).
"""
function inverse(a::UInt32, g::UInt32)
    # Extended euclidean algorithm.
    t = UInt32(0)
    T = UInt32(1)
    r = g
    R = a

    while R != 0
        q = divide(r, R)
        x = T; T = add(t, multiply(q, T)); t = x
        x = R; R = add(r, multiply(q, R)); r = x
    end

    if r > 1
        error("The element $(a) doesn't have an inverse mod $(g).")
    end

    t = modulo(t, g)

    return t
end

"""
    squareroot(p::UInt32, m)

Computes the square root of p in GF(2^m).
"""
function squareroot(p::UInt32, m)
    if p == 0
        return UInt32(0)
    elseif p == 1
        return UInt32(1)
    end

    return power(p, 2^(m - 1), primitivePolynomials[m])
end

"""
    polygcd(a::UInt32, b::UInt32)

Computes the greatest common divisor of polynomials a and b.
"""
function polygcd(a::UInt32, b::UInt32)
    while b != 0
        t = b
        b = modulo(a, b)
        a = t
    end

    return a
end

end # ExtensionField

max_stars_repo_name: rexemin/McEliece

max_issues_repo_name: rexemin/McEliece

max_forks_repo_name: rexemin/McEliece


text:
using ApproxFun, SingularIntegralEquations, Test

@testset "HierarchicalVector" begin
    a = rand(5);b = rand(5);v = HierarchicalVector((a,b)); V = HierarchicalVector((v,v));r = rand(10);

    @test isa(v,HierarchicalVector)

    @test isa(V,HierarchicalVector)

    @test norm(v-[a;b]) < 10eps()

    @test full(1+2v) == 1+2full(v)

    @test v.^2 == v.*v

    @test norm((v+r)⊖(r⊕v)) < 10eps()

    @test 2(v⋅v) == V⋅V

    @test 1+nlevels(v) == nlevels(V)

    @test abs(cumsum(v)[end]-sum(v)) < 10eps()

    @test partition(V) == (v,v)

    w = v+im*v
    u = copy(w)

    @test conj!(w) == conj(u)

    @test 1+zero(v) == ones(v)
end


text:
using SafeTestsets
using Pkg

function activate_gpu_env()
    Pkg.activate("gpu")
    Pkg.develop(PackageSpec(path=dirname(@__DIR__)))
    Pkg.instantiate()
end

const LONGER_TESTS = false

const GROUP = get(ENV, "GROUP", "All")

const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")

@time begin
  if GROUP == "All" || GROUP == "Interface"
    @time @safetestset "First Rand Tests" begin include("first_rand_test.jl") end
    @time @safetestset "Inference Tests" begin include("inference_test.jl") end
    @time @safetestset "Linear RODE Tests" begin include("rode_linear_tests.jl") end
    @time @safetestset "Complex Number Tests" begin include("complex_tests.jl") end
    @time @safetestset "Static Array Tests" begin include("static_array_tests.jl") end
    @time @safetestset "Noise Type Tests" begin include("noise_type_test.jl") end
    @time @safetestset "Mass matrix tests" begin include("mass_matrix_tests.jl") end
    #@time @safetestset "Sparse Jacobian tests" begin include("sparsediff_tests.jl") end
    @time @safetestset "Outofplace Arrays Tests" begin include("outofplace_arrays.jl") end
    @time @safetestset "tdir Tests" begin include("tdir_tests.jl") end
    @time @safetestset "tstops Tests" begin include("tstops_tests.jl") end
    @time @safetestset "saveat Tests" begin include("saveat_tests.jl") end
    @time @safetestset "Oval2" begin include("oval2_test.jl") end
  end

  if GROUP == "All" || GROUP == "Interface2"
    @time @safetestset "Basic Tau Leaping Tests" begin include("tau_leaping.jl") end
    @time @safetestset "Linear SDE Tests" begin include("sde/sde_linear_tests.jl") end
    @time @safetestset "Two-dimensional Linear SDE Tests" begin include("sde/sde_twodimlinear_tests.jl") end
    @time @safetestset "Element-wise Tolerances Tests" begin include("tolerances_tests.jl") end
    @time @safetestset "Zero'd Noise Tests" begin include("zerod_noise_test.jl") end
    @time @safetestset "Scalar Tests" begin include("scalar_noise.jl") end
    @time @safetestset "Stiffness Detection Test" begin include("stiffness_detection_test.jl") end
    @time @safetestset "Adaptive SDE Linear Tests" begin include("adaptive/sde_linearadaptive_tests.jl") end
  end

  if GROUP == "All" || GROUP == "Interface3"
    @time @safetestset "Composite Tests" begin include("composite_algorithm_test.jl") end
    @time @safetestset "Events Tests" begin include("events_test.jl") end
    @time @safetestset "Cache Tests" begin include("cache_test.jl") end
    @time @safetestset "Adaptive Complex Mean Test" begin include("adaptive/sde_complex_adaptive_mean_test.jl") end
    @time @safetestset "Utility Tests" begin include("utility_tests.jl") end
    @time @safetestset "Non-diagonal SDE Tests" begin include("nondiagonal_tests.jl") end
    @time @safetestset "No Index Tests" begin include("noindex_tests.jl") end
    @time @safetestset "Multiple Dimension Linear Adaptive Test" begin include("adaptive/sde_twodimlinearadaptive_tests.jl") end
    @time @safetestset "Autostepsize Test" begin include("adaptive/sde_autostepsize_test.jl") end
    @time @safetestset "Additive Lorenz Attractor Test" begin include("adaptive/sde_lorenzattractor_tests.jl") end
    @time @safetestset "Stochastic iterated integrals" begin include("levy_areas.jl") end
  end

  if !is_APPVEYOR && (GROUP == "All" || GROUP == "AlgConvergence")
    @time @safetestset "Convergence Tests" begin include("sde/sde_convergence_tests.jl") end
    @time @safetestset "Dynamical SDE Tests" begin include("sde/sde_dynamical.jl") end
  end

  if !is_APPVEYOR && GROUP == "AlgConvergence2"
    @time @safetestset "IIF Convergence Tests" begin include("iif_methods.jl") end
    @time @safetestset "Commutative Noise Methods Tests" begin include("commutative_tests.jl") end
    @time @safetestset "Multivariate Geometric Tests" begin include("multivariate_geometric.jl") end
  end

  if !is_APPVEYOR && GROUP == "AlgConvergence3"
    @time @safetestset "Rossler Order Tests" begin include("sde/sde_rosslerorder_tests.jl") end
    @time @safetestset "ODE Convergence Regression Tests" begin include("ode_convergence_regression.jl") end
    @time @safetestset "Additive SDE Tests" begin include("sde/sde_additive_tests.jl") end
    @time @safetestset "Split Tests" begin include("split_tests.jl") end
    @time @safetestset "Stratonovich Convergence Tests" begin include("stratonovich_convergence_tests.jl") end
  end

  if !is_APPVEYOR && GROUP == "WeakConvergence1"
    @time @safetestset "Multidimensional IIP Weak Convergence Tests" begin include("weak_convergence/multidim_iip_weak.jl") end
    @time @safetestset "Platen's PL1WM weak second order" begin include("weak_convergence/PL1WM.jl") end
  end

  if !is_APPVEYOR && GROUP == "WeakConvergence2"
    @time @safetestset "Roessler weak SRK Tests" begin include("weak_convergence/srk_weak_final.jl") end
  end

  if !is_APPVEYOR && GROUP == "WeakConvergence3"
    @time @safetestset "Roessler weak SRK (non-diagonal) Tests" begin include("weak_convergence/srk_weak_final_non_diagonal.jl") end
    @time @safetestset "Weak Stratonovich Tests" begin include("weak_convergence/weak_strat.jl") end
  end

  if !is_APPVEYOR && GROUP == "WeakConvergence4"
    @time @safetestset "Weak Stratonovich (non-diagonal) Tests" begin include("weak_convergence/weak_strat_non_diagonal.jl") end
    @time @safetestset "SIE SME weak Tests" begin include("weak_convergence/SIE_SME.jl") end
  end

  if !is_APPVEYOR && GROUP == "OOPWeakConvergence"
    @time @safetestset "OOP Weak Convergence Tests" begin include("weak_convergence/oop_weak.jl") end
    @time @safetestset "Additive Weak Convergence Tests" begin include("weak_convergence/additive_weak.jl") end
  end

  if !is_APPVEYOR && GROUP == "IIPWeakConvergence"
    #activate_gpu_env()
    @time @safetestset "IIP Weak Convergence Tests" begin include("weak_convergence/iip_weak.jl") end
  end

  if !is_APPVEYOR && GROUP == "WeakAdaptiveCPU"
    @time @safetestset "CPU Weak adaptive step size Brusselator " begin include("adaptive/sde_weak_brusselator_adaptive.jl") end
    @time @safetestset "CPU Weak adaptive" begin include("adaptive/sde_weak_adaptive.jl") end
  end

  if !is_APPVEYOR && GROUP == "WeakAdaptiveGPU"
    activate_gpu_env()
    @time @safetestset "GPU Weak adaptive step size scalar noise SDE" begin include("gpu/sde_weak_scalar_adaptive_gpu.jl") end
    @time @safetestset "GPU Weak adaptive" begin include("gpu/sde_weak_adaptive_gpu.jl") end
  end

end


text:
# this generates data according to classical model
# and shows the OLS fit
using Plots
n = 20
x = 1:n

X = [ones(n,1) x]
beta = [10, -1]
e = 6.0*randn(n)
PopRegLine = X*beta
y = PopRegLine + e

# the OLS coefficients
b = inv(X'*X)*X'*y
fit = X*b

# Plot the fitted line
plot(x,PopRegLine, label = "population regression line")
plot!(x, fit, label = "OLS fitted line")
xlabel= "X"
scatter!(x, y, label = "data")
#savefig("OlsFit.png")
#gui()

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
# Mean/variance calculations #
# ========================== #
function outcome_mean(subproblems::Vector{JuMP.Model}, probabilities::AbstractVector, sense::Union{MOI.OptimizationSense, Nothing} = nothing)
    # Sanity check
    length(subproblems) == length(probabilities) || error("Inconsistent number of subproblems and probabilities")
    N = length(subproblems)
    N == 0 && return 0.0
    if sense === nothing || sense == MOI.FEASIBILITY_SENSE
        sense = objective_sense(subproblems[1])
    end
    return mapreduce(+, 1:N) do k
        outcome = subproblems[k]
        val = try
            optimize!(outcome)
            status = termination_status(outcome)
            if !(status ∈ AcceptableTermination)
                if status == MOI.INFEASIBLE
                    objective_sense(outcome) == MOI.MAX_SENSE ? -Inf : Inf
                elseif status == MOI.DUAL_INFEASIBLE
                    objective_sense(outcome) == MOI.MAX_SENSE ? Inf : -Inf
                else
                    error("Outcome model could not be solved, returned status: $status")
                end
            else
                probabilities[k]*objective_value(outcome)
            end
        catch error
            if isa(error, NoOptimizer)
                @warn "No optimizer set, cannot solve outcome model."
                rethrow(NoOptimizer())
            else
                @warn "Outcome model could not be solved."
                rethrow(error)
            end
        end
        val = sense == objective_sense(outcome) ? val : -val
    end
end
function welford(subproblems::Vector{JuMP.Model}, probabilities::AbstractVector, sense::Union{MOI.OptimizationSense, Nothing} = nothing)
    # Sanity check
    length(subproblems) == length(probabilities) || error("Inconsistent number of subproblems and probabilities")
    N = length(subproblems)
    N == 0 && return 0.0, 0.0, 0.0, N
    if sense === nothing || sense == MOI.FEASIBILITY_SENSE
        sense = objective_sense(subproblems[1])
    end
    Q̄ₖ = Sₖ = wₖ = w²ₖ = 0
    for k = 1:N
        π = probabilities[k]
        wₖ = wₖ + π
        Q̄ₖ₋₁ = Q̄ₖ
        problem = subproblems[k]
        try
            optimize!(problem)
            status = termination_status(problem)
            Q = if !(status ∈ AcceptableTermination)
                Q = if status == MOI.INFEASIBLE
                    objective_sense(problem) == MOI.MAX_SENSE ? -Inf : Inf
                elseif status == MOI.DUAL_INFEASIBLE
                    objective_sense(problem) == MOI.MAX_SENSE ? Inf : -Inf
                else
                    error("Subproblem could not be solved during welford algorithm, returned status: $status")
                end
            else
                Q = objective_value(problem)
            end
            Q = sense == objective_sense(problem) ? Q : -Q
            Q̄ₖ = Q̄ₖ + (π / wₖ) * (Q - Q̄ₖ)
            Sₖ = Sₖ + π * (Q - Q̄ₖ) * (Q - Q̄ₖ₋₁)
        catch error
            if isa(error, NoOptimizer)
                @warn "No optimizer set, cannot solve subproblem during welford algorithm."
                rethrow(NoOptimizer())
            else
                @warn "Subproblem could not be solved during welford algorithm."
                rethrow(error)
            end
        end
    end
    correction = N / ((N - 1) * wₖ)
    if isinf(Q̄ₖ)
        correction = 0.0
    end
    return Q̄ₖ, correction * Sₖ, wₖ, N
end
function aggregate_welford(left::Tuple, right::Tuple)
    x̄ₗ, σₗ², wₗ, nₗ = left
    x̄ᵣ, σᵣ², wᵣ, nᵣ = right
    δ = x̄ᵣ - x̄ₗ
    w = wₗ + wᵣ
    n = nₗ + nᵣ
    x̄ = (wₗ * x̄ₗ + wᵣ * x̄ᵣ) / w
    Sₗ = σₗ² * ((nₗ - 1) * wₗ) / nₗ
    Sᵣ = σᵣ² * ((nᵣ - 1) * wᵣ) / nᵣ
    S = Sₗ + Sᵣ + (wₗ * wᵣ / w) * δ * δ
    correction = n / ((n - 1) * w)
    return x̄, correction * S, w, n
end


text:
@testset "floattype" begin
    function _is_fixed_type(x::Symbol)
        try
            @eval $(x) isa Type && $(x) <: FixedPoint && return true
        catch
            return false
        end
    end

    fixed_types = setdiff(filter(_is_fixed_type, names(FixedPointNumbers)), [:Fixed, :Normed, :FixedPoint])
    fixed_types = [@eval $(x) for x in fixed_types]

    exact_types = vcat([UInt8, UInt16, UInt32, UInt64, UInt128, Bool,
                       Int8, Int16, Int32, Int64, Int128],
                       fixed_types)
    for T in exact_types
        @test typemax(T) <= maxintfloat(floattype(T))
    end
end


text:
abstract type Combinatroics{T} <: AbstractArray{T,1} end

"""
    Combination{T} <: Combinatroics{T}

Structure to generated all Combinations of the given moduels with given length
low_level structure that a user don't need to use normaly
mod = Array with modulels or where combination can be generated from
len = lenght of the made constructs
"""
struct Combination{T} <: Combinatroics{T}
    mod::Array{T}
    len::Int
end

# base types of
Base.length(c::Combination) = binomial(length(c.mod),c.len);
Base.size(c::Combination) = (length(c),1);
#Base.eltype(K::Combination{T}) where {T} = Unordered_Construct{N} where N <: eltype(T) ;
Base.eltype(K::Combination{T}) where T = Unordered_Construct{T} ;

"""
    Base.iterate(c::Combination{T} where T, state = [i for i in c.len:-1:1])

Iterater for Combinations, printed in lexicografical order.
"""
function Base.iterate(c::Combination{T} where T, state = [i for i in c.len:-1:1])
    max = length(c.mod)

    if c.len == 1 # special case len = 1 ( needed?)
        if state[1] <= max
            construct = [c.mod[state[1]]]
            state[1] +=1
            return (Unordered_Construct(construct),state)
        else
            return
        end
    end

    # check when done, last constructed is printed
    if state[1] > max
        return nothing
    end

    # make the current construct based on given state
    construct = sum(c.mod[state])

    # update the state
    # i start at the last index, length of the construct
    i = c.len
    #evaluated if the difference between the last and the next is higher than 1 (than their is still a construct that can be made)
    # if not than lower the index with one and recheck if the difference is 1 . i>1 prevents error
    while i > 1  && state[i-1] - state[i] == 1
        i -=1
    end
    # upate the state at position i and reset all higher index to lowest values
    if i < c.len
        state = _restate(state,i)
    else
        state[i] +=1
    end
    return (construct,state)
end;

"""
    _restate(state,i)

Update state at index i and reset all values behind to lowsted allow values.
function used in the Base.iterate(c::Combination{T} where T, state = [i for i in c.len:-1:1])
"""
function _restate(state,i)
    state[i] +=1;
    state[(i+1):end] = collect(length(state)-i:-1:1);
    return state
end

"""
    function getindex(c::Combination,pos::Int)

Get the index of `c` at position `pos`.

Index based on rank Combination.
https://en.wikipedia.org/wiki/Combinatorial_number_system
"""
function getindex(c::Combination,pos::Int)

    # spacial case of lenght is 1
    if c.len == 1
        return Unordered_Construct([c.mod[pos]])
    end
    # if not 1
    com = Array{Int}(undef, c.len, 1)
    for i in c.len:-1:1
        n = i-1
        b_val = binomial(n,i)
        while b_val < pos
            n += 1
            b_val = binomial(n,i)
        end
        pos -= binomial(n-1,i)
        com[c.len-i+1] = n-1
    end
    index =  com .+=1
    return (sum(c.mod[index]))
end;


text:
# cartpole_play.jl
using Gym
using Knet
using ArgParse

include(joinpath(@__DIR__, "dqn", "dqn.jl"))

main() = main(ARGS)
main(args::String) = main(split(args))
function main(args::Vector{<:AbstractString})
    s = ArgParseSettings()
    # s.description = "(c) Ozan Arkan Can, 2018. An implementation of the deep q network."
    @add_arg_table s begin
        ("model"; arg_type=String; help="model name"; required=true)
        ("--frames"; arg_type=Int; default=2000; help="number of frames")
        # ("--lr"; arg_type=Float64; default=0.001; help="learning rate")
        # ("--gamma"; arg_type=Float64; default=0.99; help="discount factor")
        # ("--hiddens"; arg_type=Int; nargs='+'; default=[32]; help="number of units in the hiddens for the mlp")
        ("--env_id"; default="CartPole-v0")
        ("--render"; action=:store_true)
        # ("--memory"; arg_type=Int; default=10000; help="memory size")
        # ("--bs"; arg_type=Int; default=32; help="batch size")
        ("--stack"; arg_type=Int; default=4; help="length of the frame history")
        # ("--save"; default=""; help="model name")
        # ("--load"; default=""; help="model name")
        # ("--play"; action=:store_true; help="only play")
        ("--printinfo"; action=:store_true; help="print the playing messages")
        # ("--tupdate"; arg_type=Int; default=500; help="update frequency for the target network")
    end
    if "--help" ∈ args || "-h" ∈ args
        ArgParse.show_help(s; exit_when_done=false)
        return
    end
    o = parse_args(args, s)
    o["atype"] = Array{Float32}
    # o["load"] = o["model"]
    o["play"] = true

    env = GymEnv(o["env_id"])
    cartpoleenv = env.gymenv
    cartpoleenv[:_max_episode_steps] = typemax(Int)
    cartpoleenv[:_max_episode_seconds] = typemax(Int)

    w = DQN.load_model(o["model"], o["atype"])
    o["n_hiddens"] = length(w) ÷ 2 - 1
    o["hiddens"] = [0]  # DUMMY
    opts = Dict{String, Array{Float32}}()   # DUMMY
    buffer = DQN.ReplayBuffer(1)    # DUMMY
    exploration = DQN.PiecewiseSchedule([(0, 1.0)]) # DUMMY

    _rewards, _frames = DQN.dqn_learn(w, opts, env, buffer, exploration, o)

    close!(env)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

max_stars_repo_name: JuliaTagBot/JuliaGymDemo

max_issues_repo_name: JuliaTagBot/JuliaGymDemo

max_forks_repo_name: JuliaTagBot/JuliaGymDemo


text:
using Turing: libtask

# Utility function for self-copying mechanism
n_copies() = n_copies(current_task())
n_copies(t::Task) = begin
  isa(t.storage, Nothing) && (t.storage = IdDict())
  if haskey(t.storage, :n_copies)
    t.storage[:n_copies]
  else
    t.storage[:n_copies] = 0
  end
end

function Base.copy(t::Task)
  t.state != :runnable && t.state != :done &&
    error("Only runnable or finished tasks can be copied.")
  newt = ccall((:jl_clone_task, libtask), Any, (Any,), t)::Task
  if t.storage != nothing
    n = n_copies(t)
    t.storage[:n_copies]  = 1 + n
    newt.storage = copy(t.storage)
  else
    newt.storage = nothing
  end
  # copy fields not accessible in task.c
  newt.code = t.code
  newt.state = t.state
  newt.result = t.result
  newt.parent = t.parent
  if :last in fieldnames(typeof(t))
    newt.last = nothing
  end
  newt
end

produce(v) = begin
    ct = current_task()

    if ct.storage == nothing
        ct.storage = IdDict()
    end

    haskey(ct.storage, :consumers) || (ct.storage[:consumers] = nothing)
    local empty, t, q
    while true
        q = ct.storage[:consumers]
        if isa(q,Task)
            t = q
            ct.storage[:consumers] = nothing
            empty = true
            break
        elseif isa(q,Condition) && !isempty(q.waitq)
            t = popfirst!(q.waitq)
            empty = isempty(q.waitq)
            break
        end
        wait()
    end

    t.state == :runnable || throw(AssertionError("producer.consumer.state == :runnable"))
    if empty
        Base.schedule_and_wait(t, v)
        ct = current_task() # When a task is copied, ct should be updated to new task ID.
        while true
            # wait until there are more consumers
            q = ct.storage[:consumers]
            if isa(q,Task)
                return q.result
            elseif isa(q,Condition) && !isempty(q.waitq)
                return q.waitq[1].result
            end
            wait()
        end
    else
        schedule(t, v)
        # make sure `t` runs before us. otherwise, the producer might
        # finish before `t` runs again, causing it to see the producer
        # as done, causing done(::Task, _) to miss the value `v`.
        # see issue #7727
        yield()
        return q.waitq[1].result
    end
end

consume(p::Task, values...) = begin

    if p.storage == nothing
        p.storage = IdDict()
    end
    haskey(p.storage, :consumers) || (p.storage[:consumers] = nothing)

    if istaskdone(p)
        return wait(p)
    end

    ct = current_task()
    ct.result = length(values)==1 ? values[1] : values

    #### un-optimized version
    #if P.consumers === nothing
    #    P.consumers = Condition()
    #end
    #push!(P.consumers.waitq, ct)
    # optimized version that avoids the queue for 1 consumer
    consumers = p.storage[:consumers]
    if consumers === nothing || (isa(consumers,Condition)&&isempty(consumers.waitq))
        p.storage[:consumers] = ct
    else
        if isa(consumers, Task)
            t = consumers
            p.storage[:consumers] = Condition()
            push!(p.storage[:consumers].waitq, t)
        end
        push!(p.storage[:consumers].waitq, ct)
    end

    p.state == :runnable ? Base.schedule_and_wait(p) : wait() # don't attempt to queue it twice
end


text:
using SymPyCall
using Test

@testset "Logical" begin
    ## More logical expressions
    @syms x y
    (x ≪ 0) ∧ Le(x*y,1) ∨ (x ⩵ y) ∧ (¬(x ≫ 3))
end

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
export SWAP, FLIP
export ROT, IROT
export INC, DEC, NEG, INV, AddConst, SubConst
export HADAMARD
export PUSH!, POP!, COPYPOP!, COPYPUSH!

"""
    NoGrad{T} <: IWrapper{T}
    NoGrad(x)

A `NoGrad(x)` is equivalent to `GVar^{-1}(x)`, which cancels the `GVar` wrapper.
"""
struct NoGrad{T} <: IWrapper{T}
    x::T
end
NoGrad(x::NoGrad{T}) where T = x # to avoid ambiguity error
NoGrad{T}(x::NoGrad{T}) where T = x # to avoid ambiguity error
(_::Type{Inv{NoGrad}})(x) = x.x
@fieldview value(x::NoGrad) = x.x

const NullType{T} = Union{NoGrad{T}, Partial{T}}

NEG(a!) = -(a!)
@selfdual NEG
@selfdual -

INV(a!) = inv(a!)
@selfdual INV

@inline FLIP(b::Bool) = !b
@selfdual FLIP

"""
    INC(a!) -> a! + 1
"""
@inline function INC(a!::Number)
    a! + one(a!)
end

"""
    DEC(a!) -> a! - 1
"""
@inline function DEC(a!::Number)
    a! - one(a!)
end
@dual INC DEC

"""
    SWAP(a!, b!) -> b!, a!
"""
@inline function SWAP(a!::T, b!::T) where T
    b!, a!
end
@selfdual SWAP

"""
    ROT(a!, b!, θ) -> a!', b!', θ

```math
\\begin{align}
    {\\rm ROT}(a!, b!, \\theta)  = \\begin{bmatrix}
        \\cos(\\theta) & - \\sin(\\theta)\\\\
        \\sin(\\theta)  & \\cos(\\theta)
    \\end{bmatrix}
    \\begin{bmatrix}
        a!\\\\
        b!
    \\end{bmatrix},
\\end{align}
```
"""
@inline function ROT(i::Real, j::Real, θ::Real)
    a, b = rot(i, j, θ)
    a, b, θ
end

"""
    IROT(a!, b!, θ) -> ROT(a!, b!, -θ)
"""
@inline function IROT(i::Real, j::Real, θ::Real)
    i, j, _ = ROT(i, j, -θ)
    i, j, θ
end
@dual ROT IROT

"""
    HADAMARD(x::Real, y::Real)

Hadamard transformation that returns `(x + y)/√2, (x - y)/√2`
"""
function HADAMARD(x::Real, y::Real)
    sqrt(0.5) * (x + y), sqrt(0.5) * (x - y)
end

@selfdual HADAMARD

# more data views
for (DT, OP, NOP) in [(:AddConst, :+, :-), (:SubConst, :-, :+)]
    @eval struct $DT{T}
        x::T
    end

    @eval function (f::$DT)(y::Real)
        $OP(y, f.x)
    end

    @eval NiLangCore.chfield(x::T, ac::$DT, xval::T) where T<:Real = $NOP(xval, ac.x)
end

for F1 in [:(Base.:-), :NEG, :(ac::AddConst), :(sc::SubConst)]
    @eval @inline function $F1(a!::NullType)
        @instr $F1(a! |> value)
        a!
    end
end

for (OP, F, f) in [(:(PlusEq{typeof(identity)}), :(PlusEq(identity)), :+), (:(MinusEq{typeof(identity)}), :(MinusEq(identity)), :-)]
    @eval @inline @generated function (::$OP)(x::T, y::T) where T
        if isprimitivetype(T)
            Expr(:tuple, Expr(:call, $f, :x, :y), :y)
        else
            res = gensym("results")
            computes = Any[:($($F)(x.$field, y.$field)) for field in fieldnames(T)]
            comp = Expr(:(=), res, Expr(:tuple, computes...))
            res1 = Expr(:new, T, [:($res[$i][1]) for i=1:length(computes)]...)
            res2 = Expr(:new, T, [:($res[$i][2]) for i=1:length(computes)]...)
            quote
                $comp
                ($res1, $res2)
            end
        end
    end
    @eval (f::$OP)(x::T, y::T) where T<:Tuple = invoke(f, Tuple{T,T} where T, x, y)
    @eval (f::$OP)(x::T, y::T) where T<:Real = $f(x, y), y
end

for F2 in [:SWAP, :HADAMARD, :((inf::PlusEq)), :((inf::MinusEq)), :((inf::XorEq))]
    @eval @inline function $F2(a::NullType, b::Real)
        @instr $(NiLangCore.get_argname(F2))(a |> value, b)
        a, b
    end
    @eval @inline function $F2(a::NullType, b::NullType)
        @instr $(NiLangCore.get_argname(F2))(a |> value, b |> value)
        a, b
    end
    @eval @inline function $F2(a::Real, b::NullType)
        @instr $(NiLangCore.get_argname(F2))(a, b |> value)
        a, b
    end
end

function type_except(::Type{TT}, ::Type{T2}) where {TT, T2}
    N = length(TT.parameters)
    setdiff(Base.Iterators.product(zip(TT.parameters, repeat([T2], N))...), [ntuple(x->T2, N)])
end

for F3 in [:ROT, :IROT, :((inf::PlusEq)), :((inf::MinusEq)), :((inf::XorEq))]
    PS = (:a, :b, :c)
    for PTS in type_except(Tuple{NullType, NullType, NullType}, Real)
        params = map((P,PT)->PT <: NullType ? :($P |> value) : P, PS, PTS)
        params_ts = map((P,PT)->:($P::$PT), PS, PTS)
        @eval @inline function $F3($(params_ts...))
            @instr $F3($(params...))
            ($(PS...),)
        end
    end
end

# patch for fixed point numbers
function (f::PlusEq{typeof(/)})(out!::T, x::Integer, y::Integer) where T<:Fixed
    out!+T(x)/y, x, y
end

function (f::MinusEq{typeof(/)})(out!::T, x::Integer, y::Integer) where T<:Fixed
    out!-T(x)/y, x, y
end

for F in [:exp, :log, :sin, :sinh, :asin, :cos, :cosh, :acos, :tan, :tanh, :atan]
    @eval Base.$F(x::Fixed43) = Fixed43($F(Float64(x)))
    @eval (f::PlusEq{typeof($F)})(out!::Fixed43, x::Real) = out! + Fixed43($F(x)), x
    @eval (f::MinusEq{typeof($F)})(out!::Fixed43, x::Real) = out! - Fixed43($F(x)), x
end

Base.:^(x::Integer, y::Fixed43) = Fixed43(x^(Float64(y)))
Base.:^(x::Fixed43, y::Fixed43) = Fixed43(x^(Float64(y)))
Base.:^(x::T, y::Fixed43) where T<:AbstractFloat = x^(T(y))

function (::PlusEq{typeof(convert)})(out!::T, y) where T<:Real
    out! + convert(T, y), y
end

function (::MinusEq{typeof(convert)})(out!::T, y) where T<:Real
    out! - convert(T, y), y
end

Base.:~(ac::AddConst) = SubConst(ac.x)
Base.:~(ac::SubConst) = AddConst(ac.x)
@dualtype AddConst SubConst

for F in [:INV, :NEG, :FLIP, :INC, :DEC]
    @eval NiLangCore.chfield(x::T, ::typeof($F), xval::T) where T<:Real = (~$F)(xval)
end

#### The following functions are not safe!
@i @inline function PUSH!(x::T) where T
    PUSH!((@skip! GLOBAL_STACK), x)
end

@i @inline function POP!(x::T) where T
    POP!((@skip! GLOBAL_STACK), x)
end

@i @inline function COPYPUSH!(x)
    COPYPUSH!((@skip! GLOBAL_STACK), x)
end

@i @inline function COPYPOP!(x)
    COPYPOP!((@skip! GLOBAL_STACK), x)
end

# reversibility turned off, in principle, we can not deallocate `GVar{T}` to `T`
@i @inline function PUSH!(st, x::T) where T
    @invcheckoff st[end+1] ↔ x
    @invcheckoff x ← _zero(T)
end

@i @inline function POP!(st, x::T) where T
    @invcheckoff x → _zero(T)
    @invcheckoff st[end] ↔ (x::T)::∅
end

@i @inline function COPYPUSH!(st, x)
    @invcheckoff st[end+1] ← x
end

@i @inline function COPYPOP!(st, x)
    @invcheckoff st[end] → x
end

# accumulation on arrays: initially for Bennett algorithm
# TODO: also define it for composite types. or maybe a macro for it.
@i function :(+=)(identity)(target::AbstractArray, source::AbstractArray)
    @safe @assert length(target) == length(source)
    @inbounds for i=1:length(target)
        target[i] += source[i]
    end
end


text:
# Relays driver

mutable struct Relays
    handle
	initialized
    instr_dict
end
Relays(handle) = Relays(
                handle,
				false,
                Dict(
                    "on"=>"1",
                    "off"=>"0",
					"1"=>"on",
					"0"=>"off",
					"C1"=>"1",
					"C2"=>"2",
					"C3"=>"3",
					"C4"=>"4",
					"C5"=>"5",
					"C6"=>"6",
					"C7"=>"7",
					"C8"=>"8",
					"C9"=>"9",
                    "C10"=>"10",
                    "C11"=>"11",
                    "C12"=>"12",
                    "C13"=>"13",
                    "C14"=>"14",
                    "C15"=>"15",
                    "C16"=>"16"
                    ))

# MISC
function connect!(obj::Relays, handle)
	obj.handle = handle
    obj.connected = 1
	return nothing
end
function disconnect!(obj::Relays)
	obj.handle = 0
    obj.connected = 0
	return nothing
end
#reset_instr(obj::Relays) = write(obj.handle, "*RST")
#clear_instr(obj::Relays) = write(obj.handle, "*CLS")
get_idn(obj::Relays) = query(obj.handle, "*IDN?")

function set_state(obj::Relays, ch, st)
	ch = obj.instr_dict[ch]
	st = obj.instr_dict[st]
	cmd = ":SET:PIN$ch $st"
	write(obj.handle, cmd)
end

function get_state(obj::Relays, ch)
	ch = obj.instr_dict[ch]
	cmd = ":SET:PIN$ch?"
	st = query(obj.handle, cmd)
	st = obj.instr_dict[st]
end

function all_off(obj::Relays)
	for i=1:16
		set_state(obj, "C$i", "off")
		sleep(0.25)
	end
end

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
# Main idea of Knet8: let the user write and run an ordinary Julia
# program for the forward computation.  Extra steps: (1) certain
# variables designated as parameters, (2) certain results are paired
# with gold values and a loss function, (3) automatic gradient and
# updating.

# Motivation: be able to use all control operations of Julia without
# worrying about compiling them.

using Knet
importall Base
isdefined(:MNIST) || include("mnist.jl")

### 1. Here is the simple linear model:

w1 = convert(Array{Float32},0.1*randn(10,784)) # TODO: randn should take Float32 like rand
b1 = zeros(Float32,10)
x1 = reshape(MNIST.xtst, (784,10000)) # TODO: '*' should accept 4D arrays
y1 = w1*x1.+b1
display(y1)
@show quadloss(y1, MNIST.ytst)

### 2. Now try with special types:

# abstract Arr8                   # internal array
# type Par8 <: Arr8; a; end       # internal parameter array
# type Dat8 <: Arr8; a; end       # internal data array

# # If any argument is Arr8, the output is a Dat8.

# for op in (:*, :.+)
#     @eval begin
#         $op(a::Arr8,b::Arr8)=Dat8($op(a.a,b.a))
#         $op(a::Arr8,x)=Dat8($op(a.a,x))
#         $op(x,a::Arr8)=Dat8($op(x,a.a))
#     end
# end

# w2 = Par8(w1)
# b2 = Par8(b1)
# x2 = x1
# y2 = w2*x2.+b2
# display(y2)

# TODO: do a full register, define, minibatch, train, etc. from the
# users perspective to make sure everything works.  Then worry about
# memory efficiency.

### 3. operations with parameter arrays and anything they touch (call
# # these internal arrays) need to be recorded

# type Rec8; op; out; args; end
# Rec8(op,out,args...)=Rec8(op,out,args)

# ST8 = Rec8[]

# function call8(f, args...)
#     out=Dat8(f(args...))
#     push!(ST8, Rec8(f,out,args))
#     return out
# end

# for op in (:*, :.+)
#     @eval begin
#         $op(a::Arr8,b::Arr8)=call8($op,a.a,b.a) # TODO: lost the identity of a & b, recording arrays instead? where are the gradients going to be?
#         $op(a::Arr8,x)=call8($op,a.a,x)         # TODO: ignoring possible overwrites.
#         $op(x,a::Arr8)=call8($op,x,a.a)
#     end
# end

# w3 = Par8(w1)
# b3 = Par8(b1)
# x3 = x1
# y3 = w3*x3.+b3
# display(y3)

### 4. gold output needs to be recorded
# # during function call?
# # after function call, loss(ypred,ygold)?

# function want8(ypred::Arr8, ygold; loss=softloss)
#     # do we return loss here?  loss+grad leads to inefficiency...
#     push!(ST8, Rec8(loss, ypred, ygold))
# end

# ygold = MNIST.ytst
# want8(y3, ygold)                # TODO: wrong loss function without softmax

# # gradients:
# # should we combine gradients and back?
# # need a gradient field in Arr8
# # where do we store the back functions?  gradient registration?  function or method based?
# # In general think carefully about extending the system (by registering more gradients and optimization methods?)

### 5. new design:
# a. combine Par, Dat, and Rec into Reg.
# b. solve registration of ops and loss functions.

# """
# Reg8: DataType for marking parameters and descendents and recording operations.
# For now we'll assume Reg8.out and arrays in Reg8.args do not get overwritten.
# op: operation that created the variable.  This is typically a function.  For parameters and constants it may have a dummy value.
# args: inputs to the operation.  Only pointers are kept (not copies).
# out: the output of the operation.
# dif: gradient array, same size/type with out.
# """
# type Reg8; out; op; args; dif; end

# function writemime(io::IO, ::MIME"text/plain", r::Reg8)
#     @printf("%s,%s,%s",r.op,summary(r.out),join(map(summary,r.args),","))
# end

# # St8: global stack
# isdefined(:St8) || (St8 = Reg8[])

# "Par8: constructor for parameters."
# Par8(out)=push!(St8, Reg8(out, Par8, [], nothing))[end] # TODO: do we need to record parameters?

# # function registration example: TODO: turn this into a macro.

# for op in (:*, :.+)
#     @eval function $op(args...)
#         regs = false
#         getout(x)=(isa(x,Reg8) ? (regs=true; x.out) : x)
#         outs = map(getout, args)
#         regs || die("No registers found.")
#         out = $op(outs...)
#         reg = Reg8(out, $op, args, nothing)
#         push!(St8, reg)
#         reg
#     end
# end

# w5 = Par8(w1)
# b5 = Par8(b1)
# x5 = x1
# y5 = w5*x5.+b5
# display(y5)

# # TF implements loss calculation as part of the TF graph
# # Then it calls minimize instead of loss
# # Can we do the same?
# # We allow for loops etc but not within arrays!  i.e. kernels that modify arrays should be defined outside of Knet8.
# # Still, a loss function can be like any other operation?
# # But what if there are multiple losses like in a sequence model, need to add the losses and minimize the result.
# # Can we have scalar Knet variables, add them, and minimize the result?
# # loss forw can compute the loss value, back can compute the loss gradient.
# # loss calls also mark the places in comp graph the values we are trying to minimize. (we always minimize the sum)

# # Recording problems:
# # - arrays not changed need not be copied
# # - arrays not needed for gradient need not be copied
# # - support for both inplace and allocating ops
# # - need identity (for cumulative gradient) and current contents of the arrays
# # - cumulative vs one-shot gradients?
# # - arrays that get overwritten?
# # - can assume parameters not overwritten during a session.
# # - simple copy works for sparse, cuda etc. no need to deepcopy a.a
# # - need op, out, ins, both id & contents, also gradient arrays
# # - some of out, ins could be regular arrays, Par, or Dat.
# # - Par never change, no need to copy, gradient always cumulative.
# # - Dat and Array can be overwritten, or be used multiple times.
# # - We just need to treat overwriting as just a new array or ban it.
# # - That means id of the same array may be different after overwriting?
# # - There is also reshape we need to think about.
# # - Will we ever need multiple sessions?
# # - 3-place vs 2-place operators, memory management optimization?
# # - Function registration: a .= b*c calls A_mul_B!, do we have to register both? how do we register gradients?
# # TODO:
# # - Record
# # - Loss
# # - Backprop
# # - InplaceOps: julia 0.5 has built-in ones!
# # - Optimization

### 6. Autograd style: params=args, inputs=global, return=loss

function model6(w)
    x = reshape(MNIST.xtst, (784,10000))
    ygold = MNIST.ytst
    ypred = w[1]*x .+ w[2]
    quadloss(ypred, ygold)
end

@show model6((w1,b1))

# OK, now we have to make this work:
# grad_model6 = grad(model6)
# grad_weight = grad_model6((w1,b1))


text:
using Registrator.RegEdit: RegEdit,
    DEFAULT_REGISTRY_URL,
    parse_registry,
    showsafe,
    registration_branch,
    get_registry
using LibGit2
using Pkg.TOML
using Pkg.Types: Project

using Test

const TEST_GITCONFIG = Dict(
    "user.name" => "RegistratorTests",
    "user.email" => "ci@juliacomputing.com",
)
const TEST_SIGNATURE = LibGit2.Signature(
    TEST_GITCONFIG["user.name"],
    TEST_GITCONFIG["user.email"],
)

@testset "RegEdit" begin

@testset "Utilities" begin
    @testset "showsafe" begin
        @test string(showsafe(3)) == "3"
        @test string(showsafe(nothing)) == "nothing"
    end

    @testset "registration_branch" begin
        example = Project(Dict(
            "name" => "Example", "version" => "1.10.2"
        ))

        @test registration_branch(example) == "register/Example/v1.10.2"

        example = Project(Dict(
            "name" => "Example", "version" => nothing
        ))

        @test_throws ArgumentError registration_branch(example)
    end
end

@testset "RegistryCache" begin
    @testset "get_registry" begin
        mktempdir(@__DIR__) do temp_cache_dir
            # test when registry is not in the cache and not downloaded
            cache = RegEdit.RegistryCache(temp_cache_dir)
            repo = get_registry(DEFAULT_REGISTRY_URL, cache=cache, gitconfig=TEST_GITCONFIG)
            @test LibGit2.path(repo) == RegEdit.path(cache, DEFAULT_REGISTRY_URL)
            @test LibGit2.branch(repo) == "master"
            @test !LibGit2.isdirty(repo)
            @test LibGit2.url(LibGit2.lookup_remote(repo, "origin")) == DEFAULT_REGISTRY_URL

            # test when registry is in the cache but not downloaded
            registry_path = RegEdit.path(cache, DEFAULT_REGISTRY_URL)
            rm(registry_path, recursive=true, force=true)
            @test !ispath(registry_path)
            repo = get_registry(DEFAULT_REGISTRY_URL, cache=cache, gitconfig=TEST_GITCONFIG)
            @test LibGit2.path(repo) == RegEdit.path(cache, DEFAULT_REGISTRY_URL)
            @test LibGit2.branch(repo) == "master"
            @test !LibGit2.isdirty(repo)
            @test LibGit2.url(LibGit2.lookup_remote(repo, "origin")) == DEFAULT_REGISTRY_URL

            # test when registry is in the cache, downloaded, but mutated
            orig_hash = LibGit2.GitHash()
            LibGit2.branch!(repo, "newbranch", force=true)
            LibGit2.remove!(repo, "Registry.toml")
            LibGit2.commit(
                repo,
                "Removing Registry.toml in Registrator tests";
                author=TEST_SIGNATURE,
                committer=TEST_SIGNATURE,
            )
            @test LibGit2.GitObject(repo, "HEAD") != LibGit2.GitObject(repo, "master")
            @test ispath(registry_path)
            repo = get_registry(DEFAULT_REGISTRY_URL, cache=cache, gitconfig=TEST_GITCONFIG)
            @test LibGit2.path(repo) == RegEdit.path(cache, DEFAULT_REGISTRY_URL)
            @test LibGit2.branch(repo) == "master"
            @test !LibGit2.isdirty(repo)
            @test LibGit2.url(LibGit2.lookup_remote(repo, "origin")) == DEFAULT_REGISTRY_URL
        end
    end
end

@testset "RegistryData" begin
    blank = RegEdit.RegistryData("BlankRegistry", "d4e2f5cd-0f48-4704-9988-f1754e755b45")

    example = Project(Dict(
        "name" => "Example", "uuid" => "7876af07-990d-54b4-ab0e-23690620f79a"
    ))

    @testset "I/O" begin
        registry = """
            name = "General"
            uuid = "23338594-aafe-5451-b93e-139f81909106"
            repo = "https://github.com/JuliaRegistries/General.git"

            description = \"\"\"
            Official general Julia package registry where people can
            register any package they want without too much debate about
            naming and without enforced standards on documentation or
            testing. We nevertheless encourage documentation, testing and
            some amount of consideration when choosing package names.
            \"\"\"

            [packages]
            00701ae9-d1dc-5365-b64a-a3a3ebf5695e = { name = "BioAlignments", path = "B/BioAlignments" }
            00718b61-6157-5045-8849-3d4c4093d022 = { name = "Convertible", path = "C/Convertible" }
            0087ddc6-3964-5e57-817f-9937aefb0357 = { name = "MathOptInterfaceMosek", path = "M/MathOptInterfaceMosek" }
            """

        registry_data = parse_registry(IOBuffer(registry))
        @test registry_data isa RegEdit.RegistryData
        written_registry = sprint(TOML.print, registry_data)
        written_registry_data = parse_registry(IOBuffer(written_registry))

        @test written_registry_data == registry_data
        @test written_registry == registry
    end

    @testset "Package Operations" begin
        registry_data = copy(blank)

        @test isempty(registry_data.packages)
        @test push!(registry_data, example) == registry_data
        @test length(registry_data.packages) == 1
        @test haskey(registry_data.packages, string(example.uuid))
        @test registry_data.packages[string(example.uuid)]["name"] == "Example"
        @test registry_data.packages[string(example.uuid)]["path"] == joinpath("E", "Example")
    end
end

end


text:
# This file defines the required data structures and methods for vector-valued
# Taylor models.
# For Differentiation the `diff_wrappers.jl` are needed.
# Also, abstract types SurrogateConfig, SurrogateModel, SurrogateMeta must be defined.

@with_kw mutable struct TaylorConfig <: SurrogateConfig
    degree :: Int64 = 1;

    gradients :: Union{Symbol, Nothing, Vector{<:Function}, Function } = :fdm
    hessians ::  Union{Symbol, Nothing, Vector{<:Function}, Function } = gradients

    # alternative to specifying individual gradients
    jacobian :: Union{Symbol, Nothing, Function} = nothing

    max_evals :: Int64 = typemax(Int64);
    @assert !( ( isa(gradients, Vector) && isempty( gradients ) ) && isnothing(jacobian) ) "Provide either `gradients` or `jacobian`."
    @assert !( isnothing(gradients) && isnothing(jacobian) ) "Provide either `gradients` or `jacobian`."
    @assert 1 <= degree <= 2 "Can only construct linear and quadratic polynomial Taylor models."
end

@with_kw mutable struct TaylorModel <: SurrogateModel
    degree :: Int64 = 2;
    # reference to mop to have unscaling availabe;
    mop :: AbstractMOP
    # reference to objective(s) to evaluate
    objf :: AbstractObjective
    diff_fn :: Union{DiffFn,Nothing} = nothing
    hess_fn :: Union{HessFromGrads, DiffFn, HessWrapper, Nothing} = nothing

    # expansion point and value
    x0 :: RVec = Real[];
    fx0 :: RVec = Real[];
    # gradient(s) at x0
    g :: RVecArr = RVec[];
    H :: Vector{<:RMat} = RMat[];
end

struct TaylorMeta <: SurrogateMeta end   # no construction meta data needed

max_evals( cfg :: TaylorConfig ) = cfg.max_evals;

fully_linear( tm :: TaylorModel ) = true;
combinable( :: TaylorConfig ) = false;      # TODO think about this

# Same method as for ExactModel; duplicated for tidyness...
"Modify/initialize thec exact model `mod` so that we can differentiate it later."
function set_gradients!( mod :: TaylorModel, objf :: AbstractObjective, mop :: AbstractMOP ) :: Nothing
    cfg = model_cfg(objf);
    if isa( cfg.gradients, Symbol )
        if cfg.gradients == :autodiff
            mod.diff_fn = AutoDiffWrapper( objf )
        elseif cfg.gradients == :fdm
            mod.diff_fn = FiniteDiffWrapper( objf );
        end
    else
        if isa(cfg.gradients, Vector)
            @assert length(cfg.gradients) == num_outputs(objf) "Provide as many gradient functions as the objective has outputs."
        elseif isa(cfg.gradients, Function)
            @assert num_outputs(objf) == 1 "Only one gradient provided for $(num_outputs(objf)) outputs."
        end
        mod.diff_fn = GradWrapper(mop, cfg.gradients, cfg.jacobian )
    end
    nothing
end

function set_hessians!( mod :: TaylorModel, objf :: AbstractObjective, mop :: AbstractMOP) :: Nothing
    cfg = model_cfg(objf);
    if isa( cfg.hessians, Symbol )
        if isa( mod.diff_fn, GradWrapper )
            mod.hess_fn = HessFromGrads( mod.diff_fn, cfg.hessians );
        else
            if cfg.hessians == :autodiff
                if isa( mod.diff_fn, AutoDiffWrapper )
                    mod.hess_fn = mod.diff_fn
                else
                    mod.hess_fn = AutoDiffWrapper(objf);
                end
            elseif cfg.hessians == :fdm
                if isa( mod.diff_fn, FiniteDiffWrapper)
                    mod.hess_fn = mod.diff_fn
                else
                    mod.hess_fn = FiniteDiffWrapper( objf );
                end
            end
        end
    else
        if isa(cfg.hessians, Vector)
            @assert length(cfg.hessians) == num_outputs(objf) "Provide as many hessian functions as the objective has outputs."
        elseif isa(cfg.hessians, Function)
            @assert num_outputs(objf) == 1 "Only one hessian function provided for $(num_outputs(objf)) outputs."
        end
        mod.hess_fn = HessWrapper(mop, cfg.hessians )
    end
    nothing
end

@doc "Return a TaylorModel build from a VectorObjectiveFunction `objf`."
function _init_model( cfg ::TaylorConfig, objf :: AbstractObjective,
    mop :: AbstractMOP, id :: AbstractIterData, ac :: AbstractConfig ) :: Tuple{TaylorModel, TaylorMeta}
    tm0 = TaylorModel(; mop = mop, objf = objf );
    set_gradients!( tm0, objf, mop );
    if cfg.degree >= 2
        set_hessians!( tm0, objf, mop );
    end
    tmeta0 = TaylorMeta()
    return update_model( tm0, objf, tmeta0, mop, id, ac);
end

function update_model( tm :: TaylorModel, objf :: AbstractObjective, tmeta :: TaylorMeta,
    mop :: AbstractMOP, id :: AbstractIterData, :: AbstractConfig; ensure_fully_linear :: Bool = false ) :: Tuple{TaylorModel,TaylorMeta}
    @info "Building Taylor model(s)."
    tm.x0 = xᵗ(id);
    tm.fx0 = fxᵗ(id)[output_indices(objf,mop)];

    # set gradients
    empty!(tm.g)
    for ℓ = 1 : num_outputs(objf)
        push!(tm.g, get_gradient(tm.diff_fn, tm.x0, ℓ))
    end

    # and hessians if needed
    if !isnothing(tm.hess_fn)
        empty!(tm.H)
        for ℓ = 1 : num_outputs(objf)
            hess_mat = Matrix(get_hessian(tm.hess_fn, tm.x0, ℓ));
            push!(tm.H, hess_mat);
        end
    end
    @info "Done building Taylor model(s)."
    return tm, tmeta
end

function improve_model(tm::TaylorModel, ::AbstractObjective, tmeta :: TaylorMeta,
    ::AbstractMOP, id :: AbstractIterData, :: AbstractConfig;
    ensure_fully_linear :: Bool = false ) :: Tuple{TaylorModel, TaylorMeta}
    tm, tmeta
end

function _eval_models( tm :: TaylorModel, h :: RVec, ℓ :: Int ) :: Real
    ret_val = tm.fx0[ℓ] + tm.g[ℓ]'h
    if !isempty(tm.H)
        ret_val += .5 * h'tm.H[ℓ]*h
    end
    return ret_val
end

function eval_models( tm :: TaylorModel, x̂ :: RVec, ℓ :: Int ) :: Real
    h = x̂ .- tm.x0;
    return _eval_models( tm, h, ℓ);
 end

function eval_models( tm :: TaylorModel, x̂ :: RVec ) :: RVec
    h = x̂ .- tm.x0;
    return vcat( [_eval_models(tm, h, ℓ) for ℓ=1:num_outputs(tm.objf)]... )
end

function get_gradient( tm :: TaylorModel, x̂ :: RVec, ℓ :: Int) :: RVec
    if isempty(tm.H)
        return tm.g[ℓ]
    else
        h = x̂ .- tm.x0;
        return tm.g[ℓ] .+ .5 * ( tm.H[ℓ]' + tm.H[ℓ] ) * h
    end
end

function get_jacobian( tm :: TaylorModel, x̂ :: RVec )
    grad_list = [get_gradient(tm, x̂, ℓ) for ℓ=1:num_outputs( tm.objf ) ]
    return transpose( hcat( grad_list... ) )
end


text:
export SthenoKriging

mutable struct SthenoKriging{X, Y, GP, TΣy, GP_P} <: AbstractSurrogate
    x::X
    y::Y
    gp::GP
    Σy::TΣy
    gp_posterior::GP_P
 end

 """
     SthenoKriging(x::X, y::Y, GP::Stheno.GP=Stheno.GP(Stheno.EQ(), Stheno.GPC()), σ²=1e-18)

Returns a Kriging (or Gaussian process) surrogate conditioned on the data points `x` and `y`.
The `GP` is the base Gaussian process defined with Stheno.

# Arguments
- `x::X`: Vector containing the observation points
- `y::Y`: Vector containing the observed points
- `GP::Stheno.GP`: The base Gaussian process used to condition over. A simple prior is
provided as default. If there are multiple observation dimensions, a `Tuple` of Gaussian
processes can be passed, otherwise the same process is used across all dimensions.
- `σ²`=1e-18: Variance of the observation noise, default is equivalent to no noise
"""
function SthenoKriging(x, y, gp::Stheno.AbstractGP=Stheno.GP(Stheno.EQ(), Stheno.GPC()), σ²=1e-18)
    gpc = gp.gpc
    gps = ntuple(i -> Stheno.GP(Stheno.EQ(), gpc), length(y[1]))

    return SthenoKriging(x, y, gps, σ²)
end
function SthenoKriging(x, y, gps, σ²=1e-18)
    gp_posteriors = _prepare_gps(x, y, gps, σ²)
    return SthenoKriging(x, y, gps, σ², gp_posteriors)
end

"""
    function (k::SthenoKriging)(val)

Gives the mean predicted value for the `SthenoKriging` object.
"""
function (k::SthenoKriging)(x)
    X = _query_to_colvec(x)
    nobs = length(k.y[1])

    ŷ = [first(Stheno.mean_vector(k.gp_posterior[i], X)) for i=1:nobs]

    return _match_container(ŷ, first(k.y))
end

function std_error_at_point(k::SthenoKriging, x)
    X = _query_to_colvec(x)
    nobs = length(k.y[1])

    std_err = [sqrt(abs(first(Stheno.cov(k.gp_posterior[i], X)))) for i=1:nobs]

    return _match_container(std_err, first(k.y))
end
_query_to_colvec(x::Number) = Stheno.ColVecs(reshape([x], length(x), 1))
_query_to_colvec(x::Tuple) = Stheno.ColVecs(reshape([x...], length(x), 1))
_query_to_colvec(x) = Stheno.ColVecs(reshape(x, length(x), 1))

"""
    add_point!(k::SthenoKriging, x_new, y_new)

Adds the new point(s) and its respective value(s) to the sample points, re-conditioning
the surrogate on the updated values.
"""
function add_point!(k::SthenoKriging, x_new, y_new)
    eltype(x_new) == eltype(k.x) ? append!(k.x, x_new) : push!(k.x, x_new)
    eltype(y_new) == eltype(k.y) ? append!(k.y, y_new) : push!(k.y, y_new)

    k.gp_posterior = _prepare_gps(k.x, k.y, k.gp, k.Σy)
    return nothing
end

function _prepare_gps(x, y, gps, σ²)
    X = Stheno.ColVecs([x[j][i] for i = 1:length(x[1]), j = 1:length(x)])
    Y = [[y[i][j] for i in 1:length(y)] for j in 1:length(first(y))]
    gp_posteriors = _condition_gps(X, Y, gps, σ²)

    return gp_posteriors
end
function _condition_gps(X, Y, gp::Stheno.AbstractGP, σ²)
    nobs = length(Y)
    gp_posts = ntuple(i -> gp, nobs) | ntuple(i -> Obs(gp(X, σ²), Y[i]), nobs)

    return gp_posts
end
function _condition_gps(X, Y, gps, σ²)
    nobs = length(Y)
    @assert length(gps) == nobs
    gp_posts = gps | ntuple(i -> Obs(gps[i](X, σ²), Y[i]), nobs)

    return gp_posts
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

using MolecularGraph.Graph: merge!

@testset "graph.planarity" begin

@testset "merge" begin
    cotree = Dict(1 => 1, 2 => 2, 3 => 3, 4 => 4)
    # Any branches can be merged to empty trunk
    ds1 = Vector{Vector{Int}}[]
    ds2 = Vector{Vector{Int}}[[[1], []], [[2, 4], [3]]]
    @test merge!(ds1, ds2, cotree)
    @test ds1 == [[[1], []], [[2, 4], [3]]]

    # ds1 edge has a smaller low(e) than ds2 edges, so ds2 edges are fused
    ds1 = Vector{Vector{Int}}[[[1], []]]
    ds2 = Vector{Vector{Int}}[[[2], []], [[3], []], [[4], []]]
    @test merge!(ds1, ds2, cotree)
    @test ds1 == [[[1], []], [[2, 3, 4], []]]

    # ds1 edges have smaller low(e) than ds2 edge, so ds1 edges are not fused
    ds1 = Vector{Vector{Int}}[[[1], []], [[2], []], [[3], []]]
    ds2 = Vector{Vector{Int}}[[[4], []]]
    @test merge!(ds1, ds2, cotree)
    @test ds1 == [[[1], []], [[2], []], [[3], []], [[4], []]]

    # ds1 edges with greater low(e) than ds2 are fused. d2 is set to the opposite side
    ds1 = Vector{Vector{Int}}[[[1], []], [[3], []], [[4], []]]
    ds2 = Vector{Vector{Int}}[[[2], []]]
    @test merge!(ds1, ds2, cotree)
    @test ds1 == [[[1], []], [[3, 4], [2]]]

    # ds2 edges are fused and set to the opposite side of the ds1 cell that have greater low(e)
    ds1 = Vector{Vector{Int}}[[[1], []], [[3], []]]
    ds2 = Vector{Vector{Int}}[[[2], []], [[4], []]]
    @test merge!(ds1, ds2, cotree)
    @test ds1 == [[[1], []], [[3], [2, 4]]]

    # dichromatic branches can not be merged
    ds1 = Vector{Vector{Int}}[[[1], []]]
    ds2 = Vector{Vector{Int}}[[[2], []], [[4], [3]]]
    @test !merge!(ds1, ds2, cotree)

    # ds2 edge with smaller low(e) can not be merged to the dichromatic ds1 cell
    ds1 = Vector{Vector{Int}}[[[1], []], [[4], [3]]]
    ds2 = Vector{Vector{Int}}[[[2], []]]
    @test !merge!(ds1, ds2, cotree)

    # if low(e) of ds2 bottom is the same as that of ds1 bottom, they are merged into the same cell
    cotree = Dict(11 => 1, 12 => 1, 2 => 2, 3 => 3)
    ds1 = Vector{Vector{Int}}[[[11], []]]
    ds2 = Vector{Vector{Int}}[[[12], []], [[2], []], [[3], []]]
    @test merge!(ds1, ds2, cotree)
    @test ds1 == [[[11, 12], []], [[2, 3], []]]

    # only the top cell in ds1 has greater low(e) than ds2. d2 is set to the opposite side
    cotree = Dict(1 => 1, 21 => 2, 22 => 2, 3 => 3)
    ds1 = Vector{Vector{Int}}[[[1], []], [[21], []], [[3], []]]
    ds2 = Vector{Vector{Int}}[[[22], []]]
    @test merge!(ds1, ds2, cotree)
    @test ds1 == [[[1], []], [[21], []], [[3], [22]]]

    # bottom cells that have the same low(e) are merged into the same cell
    cotree = Dict(11 => 1, 12 => 1, 13 => 1, 14 => 1)
    ds1 = Vector{Vector{Int}}[[[11], []]]
    ds2 = Vector{Vector{Int}}[[[12, 13], []], [[14], []]]
    @test merge!(ds1, ds2, cotree)
    @test ds1 == [[[11, 12, 13], []], [[14], []]]

    # bottom cells that have the same low(e) are merged. rest ds2 edges are fused
    cotree = Dict(11 => 1, 12 => 1, 21 => 2, 22 => 2, 3 => 3)
    ds1 = Vector{Vector{Int}}[[[11], []]]
    ds2 = Vector{Vector{Int}}[[[12], []], [[21], []], [[22], []], [[3], []]]
    @test merge!(ds1, ds2, cotree)
    @test ds1 == [[[11, 12], []], [[21, 22, 3], []]]

    # ds2 edges are fused and set to the opposite side of the ds1 cell that have greater low(e)
    cotree = Dict(1 => 1, 21 => 2, 22 => 2, 23 => 2, 3 => 3)
    ds1 = Vector{Vector{Int}}[[[1], []], [[21], []], [[3], [22]]]
    ds2 = Vector{Vector{Int}}[[[23], []]]
    @test merge!(ds1, ds2, cotree)
    @test ds1 == [[[1], []], [[21], []], [[3], [22, 23]]]

    # if both side of ds1 dichromatic cell have greater low(e) than ds2 bottom low(e), ds2 can not be merged
    cotree = Dict(1 => 1, 21 => 2, 22 => 2, 3 => 3, 4 => 4)
    ds1 = Vector{Vector{Int}}[[[1], []], [[21], []], [[3], [4]]]
    ds2 = Vector{Vector{Int}}[[[22], []]]
    @test !merge!(ds1, ds2, cotree)

    # ds2 edges are fused and merged into top of ds1 cell if ds2 bottom has greater low(e) than that
    cotree = Dict(1 => 1, 2 => 2, 3 => 3, 4 => 4, 5 => 5, 6 => 6)
    ds1 = Vector{Vector{Int}}[[[1], []], [[3], [2]]]
    ds2 = Vector{Vector{Int}}[[[4, 5], []], [[6], []]]
    @test merge!(ds1, ds2, cotree)
    @test ds1 == [[[1], []], [[3], [2]], [[4, 5, 6], []]]
end

@testset "isplanar" begin
    @test planaritytest(completebipartite(2,3))
    @test !planaritytest(completebipartite(3,3))
    @test planaritytest(completegraph(4))
    @test !planaritytest(completegraph(5))
    @test planaritytest(laddergraph(20))
    @test planaritytest(circularladder(20))
    @test !planaritytest(moebiusladder(20))
    @test planaritytest(pathgraph(20))
    disconn = plaingraph(11, [
        (1, 2), (2, 3), (3, 4), (6, 7), (6, 8),
        (6, 9), (10, 7), (10, 8), (10, 9), (11, 7), (11, 8), (11, 9)
    ])
    @test !planaritytest(disconn)
end

@testset "isouterplanar" begin
    @test outerplanaritytest(completebipartite(2,2))
    @test !outerplanaritytest(completebipartite(2,3))
    @test outerplanaritytest(completegraph(3))
    @test !outerplanaritytest(completegraph(4))
    @test outerplanaritytest(laddergraph(20))
    @test !outerplanaritytest(circularladder(20))
    @test outerplanaritytest(pathgraph(20))
    disconn = plaingraph(10, [
        (1, 2), (2, 3), (3, 4), (6, 7), (6, 8),
        (6, 9), (10, 7), (10, 8), (10, 9)
    ])
    @test !outerplanaritytest(disconn)
end

end # graph.planarity


text:
"""
simulate.jl
Combine agent and stock level functions into a simulation framework
"""

function simulate(carrying_capacity::Vector, effort::Vector, bump::Vector, s_db::StockDB, s_a::StockAssumptions, a_a::AgentAssumptions, e_a::EnvironmentAssumptions, reduced=true::Bool, progress=true::Bool, limit=250000::Int64)
  """
  Brings together all of the functions necessary for a life cycle simulation
  """
  @assert(all(carrying_capacity .> 0.), "There is at least one negative carrying capacity")
  @assert(length(effort)<=length(carrying_capacity), "The effort vector must be equal or less than the length of the simulation")
  @assert(length(bump)<=length(carrying_capacity), "The bump vector must be equal or less than the length of the simulation")
  years = length(carrying_capacity)
  a_db = AgentDB(years, a_a, reduced)
  bumpvec = fill(0, years)
  bumpvec[1:length(bump)] = bump
  harvest_effort = fill(0., years)
  harvest_effort[1:length(effort)] = effort
  if reduced
    c=fill(0, 104)
    for i = 1:104
      c[i]=find([0, a_a.growth[1:end-1]] .<= i)[end]
    end
  else
    c = 1:104
  end
  for y = 1:1
    spawn!(a_db, s_db, s_a, e_a, y, carrying_capacity[y])
    totalagents=size(a_db[y,1],1)
    if progress
      progressbar = Progress(52, 5, "Year $y simulation progress ($totalagents agents)", 30)
    end
    for w = 1:52
      if w > 1 && c[w] - c[w-1] == 1
        a_db[y,c[w]] = deepcopy(a_db[y,c[w-1]])
      end
      kill!(a_db, e_a, a_a, y, c[w])
      move!(a_db, a_a, e_a, y, c[w])
      if w==52
        harvest!(harvest_effort[y], s_db, s_a)
        ageadults!(s_db, s_a, carrying_capacity[y], bumpvec[y])
      end
      graduate!(a_db, s_db, a_a, y, w, c[w])
      if progress
        next!(progressbar)
      end
    end
  end
  for y = 2:years
    spawn!(a_db, s_db, s_a, e_a, y, carrying_capacity[y])
    totalagents=size(a_db[y,1],1)
    @assert(totalagents < limit, "> $limit agents in current simulation, stopping here.")
    if progress
      progressbar = Progress(52, 5, "Year $y simulation progress ($totalagents agents)", 30)
    end
    for w = 1:52
      if w > 1 && c[w] - c[w-1] == 1
        a_db[y,c[w]] = deepcopy(a_db[y,c[w-1]])
      end
      if c[w+52] - c[w+51] == 1
        a_db[y-1,c[w+52]] = deepcopy(a_db[y-1,c[w+51]])
      end
      kill!(a_db, e_a, a_a, y, c[w])
      kill!(a_db, e_a, a_a, y-1, c[w+52])
      move!(a_db, a_a, e_a, y, c[w])
      move!(a_db, a_a, e_a, y-1, c[w+52])
      if w==52
        harvest!(harvest_effort[y], s_db, s_a)
        ageadults!(s_db, s_a, carrying_capacity[y], bumpvec[y])
      end
      graduate!(a_db, s_db, a_a, y, w, c[w])
      graduate!(a_db, s_db, a_a, y-1, w+52, c[w+52])
      if progress
        next!(progressbar)
      end
    end
  end
  return a_db
end


text:
## ------------------------------------------------------------------
mutable struct Ref
    title::String
    bibkey::String
    year::String
    doi::String
    authors::String
    keywords::String
    abstract::String
    refdict::Dict

    Ref() = new("", "", "", "", "", "", "", Dict())
end

## ------------------------------------------------------------------
# Meta

get_refdict(r::Ref) = getproperty(r, :refdict)
set_refdict!(r::Ref, dict::AbstractDict) = setproperty!(r, :refdict, dict)

get_bibkey(r::Ref) = getproperty(r, :bibkey)
set_bibkey!(r::Ref, key::AbstractString) = setproperty!(r, :bibkey, strip(key))

get_year(r::Ref) = getproperty(r, :year)
set_year!(r::Ref, year::AbstractString) = setproperty!(r, :year, strip(year))

get_title(r::Ref) = getproperty(r, :title)
set_title!(r::Ref, title::AbstractString) = setproperty!(r, :title, strip(title))

get_doi(r::Ref) = getproperty(r, :doi)
set_doi!(r::Ref, doi::AbstractString) = setproperty!(r, :doi, _doi_to_url(doi))

get_keywords(r::Ref) = getproperty(r, :keywords)
set_keywords!(r::Ref, keywords::AbstractString) = setproperty!(r, :keywords, strip(keywords))

get_authors(r::Ref) = getproperty(r, :authors)
set_authors!(r::Ref, authors::AbstractString) = setproperty!(r, :authors, strip(authors))

get_abstract(r::Ref) = getproperty(r, :abstract)
set_abstract!(r::Ref, abstract::AbstractString) = setproperty!(r, :abstract, strip(abstract))

## ------------------------------------------------------------------
# Base
function Base.show(io::IO, r::Ref)
    print(io, "Ref")
    for meta in [:title, :bibkey, :authors, :year, :doi]
        val = getproperty(r, meta)
        isempty(val) && continue
        print(io, "\n - ", meta, ": \"", val, "\"")
    end
end

## ------------------------------------------------------------------
# Dict API
Base.getindex(r::Ref, keys...) = getindex(get_refdict(r), keys...)
Base.get(r::Ref, key, default) = get(get_refdict(r), key, default)


text:
using Gadfly

Φs = [:Φground, :Φexcited]
ss = [0.001, 0.01, 0.1]                    # wall thickness

for Φ0 in Φs
    for s in ss

    # Plot the log of escape probability
    decay = readdlm("trans_pocket_$(string(Φ0))\_$s.dat")
    tstart = 5
    ts = 50*(tstart-1+[tstart:size(decay)[1]])*1e-6
    lin = linreg(ts, log(decay[tstart:end,1]))
    bestfit = x->lin[2]*x + lin[1]          # Linear Regression
    p = plot(
        layer(
            x = ts,
            y = log(decay[tstart:end,1]),
            Geom.line,
            color = ["log(T)"]
        ),
#       layer(
#           x = ts,
#           y = log(decay[tstart:end,2]),
#           Geom.line,
#           color = ["log(R)"]
#       ),
        layer(
            x = ts,
            y = bestfit(ts),
            Geom.line,
            color = ["Linear Regression of log(T), a = $(lin[2])"]
        ),
        Guide.xlabel("Number of Steps"),
        Guide.ylabel("log(Probability)"),
        Guide.title("Log of Leakage Probability vs. Time for $(string(Φ0)) initial state, s = $s")
    )
    draw(SVG("decay_$(string(Φ0))\_$s.svg", 20cm, 12cm), p)

    # Plot final vs initial configuration
    phis = readdlm("phis_pocket_re_$(string(Φ0))\_$s.dat")
    xs = 1e-3*[0:2000]-1.0
    t = 50size(phis)[2]*1e-6
    q = plot(
        layer(
            x = xs,
            y = phis[:,1],
            Geom.line,
            color = ["Wavefunction at t=0"]
        ),
        layer(
            x = xs,
            y = phis[:,end],
            Geom.line,
            color = ["Wavefunction at t=$t"]
        ),
        Guide.xlabel("x"),
        Guide.ylabel("Re(ϕ)"),
        Guide.title("ϕ initial and final starting from $(string(Φ0)) with wall thickness s = $s")
    )
    draw(SVG("wave_$(string(Φ0))\_$s.svg", 20cm, 12cm), q)

    end
end