text:
const use_return  = haskey(ENV, "USE_RETURN" ) ? parse(Bool, ENV["USE_RETURN"] ) : false
const USE_GPU     = haskey(ENV, "USE_GPU"    ) ? parse(Bool, ENV["USE_GPU"]    ) : false
const do_viz      = haskey(ENV, "DO_VIZ"     ) ? parse(Bool, ENV["DO_VIZ"]     ) : true
const do_save     = haskey(ENV, "DO_SAVE"    ) ? parse(Bool, ENV["DO_SAVE"]    ) : false
const do_save_viz = haskey(ENV, "DO_SAVE_VIZ") ? parse(Bool, ENV["DO_SAVE_VIZ"]) : false
const nx          = haskey(ENV, "NX"         ) ? parse(Int , ENV["NX"]         ) : 64
const ny          = haskey(ENV, "NY"         ) ? parse(Int , ENV["NY"]         ) : 64
const nz          = haskey(ENV, "NZ"         ) ? parse(Int , ENV["NZ"]         ) : 64
###
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 3)
else
    @init_parallel_stencil(Threads, Float64, 3)
end
using ImplicitGlobalGrid, Plots, Printf, LinearAlgebra, MAT
import MPI

norm_g(A) = (sum2_l = sum(A.^2); sqrt(MPI.Allreduce(sum2_l, MPI.SUM, MPI.COMM_WORLD)))

@views inn(A) = A[2:end-1,2:end-1,2:end-1]

@parallel function compute_iter_params!(dτ_ρ, D, Re, Vpdτ, max_lxyz)
    @all(dτ_ρ)  = Vpdτ * max_lxyz / @maxloc(D) / Re
    return
end

@parallel function compute_flux!(qHx, qHy, qHz, qHx2, qHy2, qHz2, H, D, θr_dτ, dx, dy, dz)
    @all(qHx)  = (@all(qHx) * θr_dτ - @av_xi(D) * @d_xi(H) / dx) / (1.0 + θr_dτ)
    @all(qHy)  = (@all(qHy) * θr_dτ - @av_yi(D) * @d_yi(H) / dy) / (1.0 + θr_dτ)
    @all(qHz)  = (@all(qHz) * θr_dτ - @av_zi(D) * @d_zi(H) / dz) / (1.0 + θr_dτ)
    @all(qHx2) = -@av_xi(D) * @d_xi(H) / dx
    @all(qHy2) = -@av_yi(D) * @d_yi(H) / dy
    @all(qHz2) = -@av_zi(D) * @d_zi(H) / dz
    return
end

@parallel function compute_update!(H, Hold, qHx, qHy, qHz, dτ_ρ, dt, dx, dy, dz)
    @inn(H) = (@inn(H) + @all(dτ_ρ) * (@inn(Hold) / dt - (@d_xa(qHx) / dx + @d_ya(qHy) / dy + @d_za(qHz)/dz))) / (1.0 + @all(dτ_ρ) / dt)
    return
end

@parallel function check_res!(ResH, H, Hold, qHx2, qHy2, qHz2, dt, dx, dy, dz)
    @all(ResH) = -(@inn(H) - @inn(Hold)) / dt - (@d_xa(qHx2) / dx + @d_ya(qHy2) / dy + @d_za(qHz2) / dz)
    return
end

@parallel_indices (iy,iz) function bc_x!(A)
    A[1  , iy, iz] = A[2    , iy, iz]
    A[end, iy, iz] = A[end-1, iy, iz]
    return
end

@parallel_indices (ix,iz) function bc_y!(A)
    A[ix, 1  , iz] = A[ix, 2    , iz]
    A[ix, end, iz] = A[ix, end-1, iz]
    return
end

@parallel_indices (ix,iy) function bc_z!(A)
    A[ix, iy, 1  ] = A[ix, iy, 2    ]
    A[ix, iy, end] = A[ix, iy, end-1]
    return
end

@views function diffusion_3D_()
    # Physics
    lx, ly, lz = 10.0, 10.0, 10.0 # domain size
    D1         = 1.0              # diffusion coefficient
    D2         = 1e-4             # diffusion coefficient
    ttot       = 1.0              # total simulation time
    dt         = 0.2              # physical time step
    # Numerics
    tol        = 1e-8             # tolerance
    itMax      = 1e5              # max number of iterations
    nout       = 10               # tol check
    CFL        = 1/sqrt(3)        # CFL number
    me, dims   = init_global_grid(nx, ny, nz) # MPI initialisation
    @static if USE_GPU select_device() end    # select one GPU per MPI local rank (if >1 GPU per node)
    b_width    = (8, 4, 4)       # boundary width for comm/comp overlap
    # Derived numerics
    dx, dy, dz = lx/nx_g(), ly/ny_g(), lz/nz_g() # cell sizes
    Vpdτ       = CFL * min(dx, dy, dz)
    max_lxyz   = max(lx, ly, lz)
    Re         = π + sqrt(π^2 + (max_lxyz^2 / max(D1,D2)) / dt)
    θr_dτ      = max_lxyz / Vpdτ / Re
    xc, yc, zc = LinRange(dx/2, lx-dx/2, nx), LinRange(dy/2, ly-dy/2, ny), LinRange(dz/2, lz-dz/2, nz)
    # Array allocation
    qHx        = @zeros(nx-1,ny-2,nz-2)
    qHy        = @zeros(nx-2,ny-1,nz-2)
    qHz        = @zeros(nx-2,ny-2,nz-1)
    qHx2       = @zeros(nx-1,ny-2,nz-2)
    qHy2       = @zeros(nx-2,ny-1,nz-2)
    qHz2       = @zeros(nx-2,ny-2,nz-1)
    ResH       = @zeros(nx-2,ny-2,nz-2)
    dτ_ρ       = @zeros(nx-2,ny-2,nz-2)
    # Initial condition
    H0         =   zeros(nx,ny,nz)
    D          = D1*ones(nx,ny,nz)
    Tmp        = [x_g(ix,dx,H0) for ix=1:size(H0,1), iy=1:size(H0,2), iz=1:size(H0,3)]
    D[Tmp.<lx/2.2] .= D2
    Tmp        = [y_g(iy,dy,H0) for ix=1:size(H0,1), iy=1:size(H0,2), iz=1:size(H0,3)]
    D[Tmp.<ly/2.2] .= D2
    D          = Data.Array(D)
    H0         = Data.Array([exp(-(x_g(ix,dx,H0)-0.5*lx+dx/2)*(x_g(ix,dx,H0)-0.5*lx+dx/2) - (y_g(iy,dy,H0)-0.5*ly+dy/2)*(y_g(iy,dy,H0)-0.5*ly+dy/2) - (z_g(iz,dz,H0)-0.5*lz+dz/2)*(z_g(iz,dz,H0)-0.5*lz+dz/2)) for ix=1:size(H0,1), iy=1:size(H0,2), iz=1:size(H0,3)])
    Hold       = @ones(nx,ny,nz) .* H0
    H          = @ones(nx,ny,nz) .* H0
    @parallel compute_iter_params!(dτ_ρ, D, Re, Vpdτ, max_lxyz)
    len_ResH_g = ((nx-2-2)*dims[1]+2)*((ny-2-2)*dims[2]+2)*((nz-2-2)*dims[3]+2)
    if do_viz || do_save_viz
        if (me==0) ENV["GKSwstype"]="nul"; if do_viz !ispath("../../figures") && mkdir("../../figures") end; end
        nx_v, ny_v, nz_v = (nx-2)*dims[1], (ny-2)*dims[2], (nz-2)*dims[3]
        if (nx_v*ny_v*nz_v*sizeof(Data.Number) > 0.8*Sys.free_memory()) error("Not enough memory for visualization.") end
        H_v    = zeros(nx_v, ny_v, nz_v) # global array for visu
        H_inn  = zeros(nx-2, ny-2, nz-2) # no halo local array for visu
        z_sl   = Int(ceil(nz_g()/2))     # Central z-slice
        Xi_g, Yi_g = dx+dx/2:dx:lx-dx-dx/2, dy+dy/2:dy:ly-dy-dy/2 # inner points only
    end
    t = 0.0; it = 0; ittot = 0; nt = Int(ceil(ttot/dt))
    # Physical time loop
    while it<nt
        iter = 0; err = 2 * tol
        # Pseudo-transient iteration
        while err>tol && iter<itMax
            @parallel compute_flux!(qHx, qHy, qHz, qHx2, qHy2, qHz2, H, D, θr_dτ, dx, dy, dz)
            @hide_communication b_width begin # communication/computation overlap
                @parallel compute_update!(H, Hold, qHx, qHy, qHz, dτ_ρ, dt, dx, dy, dz)
                update_halo!(H)
            end
            iter += 1
            if iter % nout == 0
                @parallel check_res!(ResH, H, Hold, qHx2, qHy2, qHz2, dt, dx, dy, dz)
                err = norm_g(ResH) / sqrt(len_ResH_g)
            end
        end
        ittot += iter; it += 1; t += dt
        Hold .= H
        if isnan(err) error("NaN") end
    end
    if (me==0) @printf("Total time = %1.2f, time steps = %d, nx = %d, iterations tot = %d \n", round(ttot, sigdigits=2), it, nx_g(), ittot) end
    # Visualise
    if do_viz || do_save_viz
        H_inn .= inn(H); gather!(H_inn, H_v)
        if me==0 && do_viz
            heatmap(Xi_g, Yi_g, H_v[:,:,z_sl]', dpi=150, aspect_ratio=1, framestyle=:box, xlims=(Xi_g[1],Xi_g[end]), ylims=(Yi_g[1],Yi_g[end]), xlabel="lx", ylabel="ly", c=:viridis, clims=(0,1), title="linear step diffusion (nt=$it, iters=$ittot)")
            savefig("../../figures/diff_3D_linstep_$(nx_g()).png")
        end
    end
    if me==0 && do_save
        !ispath("../../output") && mkdir("../../output")
        open("../../output/out_diff_3D_linstep.txt","a") do io
            println(io, "$(nx_g()) $(ny_g()) $(nz_g()) $(ittot) $(nt)")
        end
    end
    if me==0 && do_save_viz
        !ispath("../../out_visu") && mkdir("../../out_visu")
        matwrite("../../out_visu/diff_3D_linstep.mat", Dict("H_3D"=> Array(H_v), "xc_3D"=> Array(xc), "yc_3D"=> Array(yc), "zc_3D"=> Array(zc)); compress = true)
    end
    finalize_global_grid()
    return xc, yc, zc, H
end

if use_return
    xc, yc, zc, H = diffusion_3D_();
else
    diffusion_3D = begin diffusion_3D_(); return; end
end


text:
using Test
using Suppressor

function countdown(i)
    println(i)
    if i <= 0
        return
    end
    countdown(i-1)
end

result = @capture_out(countdown(10)) # return stdout result
@test result == "10
9
8
7
6
5
4
3
2
1
0
"

ext: jl


text:
"""
Collection of utilities for working with Requests data
"""
module Requests

using Genie, Genie.Router, Genie.Input
using HTTP

export jsonpayload, rawpayload, filespayload, postpayload, getpayload, getrequest, infilespayload, download, filename, payload

"""
    jsonpayload() :: Dict{String,Any}

Processes an `application/json` `POST` request attempting to convert the payload into a `Dict`.
If it fails to successfully parse the `JSON` data this returns `nothing`. The original payload can still be accessed invoking `rawpayload()`
"""
@inline function jsonpayload() :: Dict{String,Any}
  @params(Genie.PARAMS_JSON_PAYLOAD)
end

"""
    rawpayload() :: String

Returns the raw `POST` payload as a `String`.
"""
@inline function rawpayload() :: String
  @params(Genie.PARAMS_RAW_PAYLOAD)
end

"""
    filespayload() :: Dict{String,HttpFile}

Collection of form uploaded files.
"""
@inline function filespayload() :: Dict{String,HttpFile}
  @params(Genie.PARAMS_FILES)
end

"""
    filespayload(filename::Union{String,Symbol}) :: HttpFile

Returns the `HttpFile` uploaded through the `key` input name.
"""
@inline function filespayload(key::Union{String,Symbol}) :: HttpFile
  @params(Genie.PARAMS_FILES)[string(key)]
end

"""
    infilespayload(key::Union{String,Symbol}) :: Bool

Checks if the collection of uploaded files contains a file stored under the `key` name.
"""
@inline function infilespayload(key::Union{String,Symbol}) :: Bool
  haskey(filespayload(), string(key))
end

"""
    Base.download(file::HttpFile, name::String = file.name)

Saves uploaded `HttpFile` `file` to local storage under the `name` filename.
"""
@inline function Base.download(file::HttpFile, name::String = file.name) :: Int
  write(name, IOBuffer(file.data))
end

"""
    filename(file::HttpFile) :: String

Original filename of the uploaded `HttpFile` `file`.
"""
@inline function filename(file::HttpFile) :: String
  file.name
end

"""
    postpayload() :: Dict{Symbol,Any}

A dict representing the POST variables payload of the request (corresponding to a `form-data` request)
"""
@inline function postpayload() :: Dict{Symbol,Any}
  @params(Genie.PARAMS_POST_KEY)
end

"""
    postpayload(key::Symbol) :: Any

Returns the value of the POST variables `key`.
"""
@inline function postpayload(key::Symbol)
  postpayload()[key]
end

"""
    getpayload(key::Symbol, default::Any)

Returns the value of the POST variables `key` or the `default` value if `key` is not defined.
"""
@inline function postpayload(key::Symbol, default::Any)
  haskey(postpayload(), key) ? postpayload(key) : default
end

"""
    getpayload() :: Dict{Symbol,Any}

A dict representing the GET/query variables payload of the request (the part correspoding to `?foo=bar&baz=moo`)
"""
@inline function getpayload() :: Dict{Symbol,Any}
  @params(Genie.PARAMS_GET_KEY)
end

"""
    getpayload(key::Symbol) :: Any

The value of the GET/query variable `key`, as in `?key=value`
"""
@inline function getpayload(key::Symbol) :: Any
  getpayload()[key]
end

"""
    getpayload(key::Symbol, default::Any) :: Any

The value of the GET/query variable `key`, as in `?key=value`. If `key` is not defined, `default` is returned.
"""
@inline function getpayload(key::Symbol, default::Any) :: Any
  haskey(getpayload(), key) ? getpayload(key) : default
end

"""
    getrequest() :: HTTP.Request

Returns the raw HTTP.Request object associated with the request.
"""
@inline function getrequest() :: HTTP.Request
  @params(Genie.PARAMS_REQUEST_KEY)
end

"""
    payload() :: Any

Utility function for accessing the `@params` collection, which holds the request variables.
"""
@inline function payload() :: Any
  @params
end

"""
    payload(key::Symbol) :: Any

Utility function for accessing the `key` value within the `@params` collection of request variables.
"""
@inline function payload(key::Symbol) :: Any
  @params()[key]
end

end


text:
export BruteForce

"""
    bruteforcesearch(data, metric, query, t::SearchType[, skip])

Perform a brute-force search of type `t` against data array `data`
(by calculating the metric for `query` and and every point in `data`).
"""
function bruteforcesearch end

function bruteforcesearch(data, metric, query, t::NeighborNumber, skip=nothing)
    dists = [metric(query, d) for d in data]
    indices = sortperm(dists)
    !isnothing(skip) && filter!(i -> !skip(i), indices)
    length(indices) > t.k && resize!(indices, t.k)
    return indices, dists[indices]
end

function bruteforcesearch(data, metric, query, t::WithinRange, skip=nothing)
    indices = Int[]
    dists = metricreturntype(metric, first(data))[]
    for (i, datum) in enumerate(data)
        !isnothing(skip) && skip(i) && continue
        d = metric(query, datum)
        if d <= t.r
            push!(indices, i)
            push!(dists, d)
        end
    end
    return indices, dists
end

"""
    BruteForce

A "search structure" which simply performs a brute-force search through the
entire data array.
"""
struct BruteForce{T, M}
    data::Vector{T}
    metric::M

    BruteForce(data::Vector, metric) = new{eltype(data), typeof(metric)}(data, metric)
end

searchstructure(::Type{BruteForce}, data, metric) = BruteForce(data, metric)

datatype(::Type{<:BruteForce{T}}) where T = T
getmetric(bf::BruteForce) = bf.metric

function search(bf::BruteForce, query, t::SearchType, skip=alwaysfalse)
    bruteforcesearch(bf.data, bf.metric, query, t, skip)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

#### Specialized matrix types ####

## (complex) symmetric tridiagonal matrices
struct SymTridiagonal{T,V<:AbstractVector{T}} <: AbstractMatrix{T}
    dv::V                        # diagonal
    ev::V                        # subdiagonal
    function SymTridiagonal{T}(dv::V, ev::V) where {T,V<:AbstractVector{T}}
        if !(length(dv) - 1 <= length(ev) <= length(dv))
            throw(DimensionMismatch("subdiagonal has wrong length. Has length $(length(ev)), but should be either $(length(dv) - 1) or $(length(dv))."))
        end
        new{T,V}(dv,ev)
    end
end

"""
    SymTridiagonal(dv::V, ev::V) where V <: AbstractVector

Construct a symmetric tridiagonal matrix from the diagonal (`dv`) and first
sub/super-diagonal (`ev`), respectively. The result is of type `SymTridiagonal`
and provides efficient specialized eigensolvers, but may be converted into a
regular matrix with [`convert(Array, _)`](@ref) (or `Array(_)` for short).

# Examples
```jldoctest
julia> dv = [1, 2, 3, 4]
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> ev = [7, 8, 9]
3-element Array{Int64,1}:
 7
 8
 9

julia> SymTridiagonal(dv, ev)
4×4 SymTridiagonal{Int64,Array{Int64,1}}:
 1  7  ⋅  ⋅
 7  2  8  ⋅
 ⋅  8  3  9
 ⋅  ⋅  9  4
```
"""
SymTridiagonal(dv::V, ev::V) where {T,V<:AbstractVector{T}} = SymTridiagonal{T}(dv, ev)

"""
    SymTridiagonal(A::AbstractMatrix)

Construct a symmetric tridiagonal matrix from the diagonal and
first sub/super-diagonal, of the symmetric matrix `A`.

# Examples
```jldoctest
julia> A = [1 2 3; 2 4 5; 3 5 6]
3×3 Array{Int64,2}:
 1  2  3
 2  4  5
 3  5  6

julia> SymTridiagonal(A)
3×3 SymTridiagonal{Int64,Array{Int64,1}}:
 1  2  ⋅
 2  4  5
 ⋅  5  6
```
"""
function SymTridiagonal(A::AbstractMatrix)
    if diag(A,1) == diag(A,-1)
        SymTridiagonal(diag(A,0), diag(A,1))
    else
        throw(ArgumentError("matrix is not symmetric; cannot convert to SymTridiagonal"))
    end
end

convert(::Type{SymTridiagonal{T}}, S::SymTridiagonal) where {T} =
    SymTridiagonal(convert(AbstractVector{T}, S.dv), convert(AbstractVector{T}, S.ev))
convert(::Type{AbstractMatrix{T}}, S::SymTridiagonal) where {T} =
    SymTridiagonal(convert(AbstractVector{T}, S.dv), convert(AbstractVector{T}, S.ev))
function convert(::Type{Matrix{T}}, M::SymTridiagonal) where T
    n = size(M, 1)
    Mf = zeros(T, n, n)
    @inbounds begin
        @simd for i = 1:n-1
            Mf[i,i] = M.dv[i]
            Mf[i+1,i] = M.ev[i]
            Mf[i,i+1] = M.ev[i]
        end
        Mf[n,n] = M.dv[n]
    end
    return Mf
end
convert(::Type{Matrix}, M::SymTridiagonal{T}) where {T} = convert(Matrix{T}, M)
convert(::Type{Array}, M::SymTridiagonal) = convert(Matrix, M)
full(M::SymTridiagonal) = convert(Array, M)

size(A::SymTridiagonal) = (length(A.dv), length(A.dv))
function size(A::SymTridiagonal, d::Integer)
    if d < 1
        throw(ArgumentError("dimension must be ≥ 1, got $d"))
    elseif d<=2
        return length(A.dv)
    else
        return 1
    end
end

similar(S::SymTridiagonal, ::Type{T}) where {T} = SymTridiagonal{T}(similar(S.dv, T), similar(S.ev, T))

#Elementary operations
broadcast(::typeof(abs), M::SymTridiagonal) = SymTridiagonal(abs.(M.dv), abs.(M.ev))
broadcast(::typeof(round), M::SymTridiagonal) = SymTridiagonal(round.(M.dv), round.(M.ev))
broadcast(::typeof(trunc), M::SymTridiagonal) = SymTridiagonal(trunc.(M.dv), trunc.(M.ev))
broadcast(::typeof(floor), M::SymTridiagonal) = SymTridiagonal(floor.(M.dv), floor.(M.ev))
broadcast(::typeof(ceil), M::SymTridiagonal) = SymTridiagonal(ceil.(M.dv), ceil.(M.ev))
for func in (:conj, :copy, :real, :imag)
    @eval ($func)(M::SymTridiagonal) = SymTridiagonal(($func)(M.dv), ($func)(M.ev))
end
broadcast(::typeof(round), ::Type{T}, M::SymTridiagonal) where {T<:Integer} = SymTridiagonal(round.(T, M.dv), round.(T, M.ev))
broadcast(::typeof(trunc), ::Type{T}, M::SymTridiagonal) where {T<:Integer} = SymTridiagonal(trunc.(T, M.dv), trunc.(T, M.ev))
broadcast(::typeof(floor), ::Type{T}, M::SymTridiagonal) where {T<:Integer} = SymTridiagonal(floor.(T, M.dv), floor.(T, M.ev))
broadcast(::typeof(ceil), ::Type{T}, M::SymTridiagonal) where {T<:Integer} = SymTridiagonal(ceil.(T, M.dv), ceil.(T, M.ev))

transpose(M::SymTridiagonal) = M #Identity operation
adjoint(M::SymTridiagonal) = conj(M)

function diag(M::SymTridiagonal{T}, n::Integer=0) where T
    absn = abs(n)
    if absn == 0
        return M.dv
    elseif absn==1
        return M.ev
    elseif absn<size(M,1)
        return zeros(T,size(M,1)-absn)
    else
        throw(ArgumentError("$n-th diagonal of a $(size(M)) matrix doesn't exist!"))
    end
end

+(A::SymTridiagonal, B::SymTridiagonal) = SymTridiagonal(A.dv+B.dv, A.ev+B.ev)
-(A::SymTridiagonal, B::SymTridiagonal) = SymTridiagonal(A.dv-B.dv, A.ev-B.ev)
*(A::SymTridiagonal, B::Number) = SymTridiagonal(A.dv*B, A.ev*B)
*(B::Number, A::SymTridiagonal) = A*B
/(A::SymTridiagonal, B::Number) = SymTridiagonal(A.dv/B, A.ev/B)
==(A::SymTridiagonal, B::SymTridiagonal) = (A.dv==B.dv) && (A.ev==B.ev)

function A_mul_B!(C::StridedVecOrMat, S::SymTridiagonal, B::StridedVecOrMat)
    m, n = size(B, 1), size(B, 2)
    if !(m == size(S, 1) == size(C, 1))
        throw(DimensionMismatch("A has first dimension $(size(S,1)), B has $(size(B,1)), C has $(size(C,1)) but all must match"))
    end
    if n != size(C, 2)
        throw(DimensionMismatch("second dimension of B, $n, doesn't match second dimension of C, $(size(C,2))"))
    end

    α = S.dv
    β = S.ev
    @inbounds begin
        for j = 1:n
            x₀, x₊ = B[1, j], B[2, j]
            β₀ = β[1]
            C[1, j] = α[1]*x₀ + x₊*β₀
            for i = 2:m - 1
                x₋, x₀, x₊ = x₀, x₊, B[i + 1, j]
                β₋, β₀ = β₀, β[i]
                C[i, j] = β₋*x₋ + α[i]*x₀ + β₀*x₊
            end
            C[m, j] = β₀*x₀ + α[m]*x₊
        end
    end

    return C
end

(\)(T::SymTridiagonal, B::StridedVecOrMat) = ldltfact(T)\B

eigfact!(A::SymTridiagonal{<:BlasReal}) = Eigen(LAPACK.stegr!('V', A.dv, A.ev)...)
function eigfact(A::SymTridiagonal{T}) where T
    S = promote_type(Float32, typeof(zero(T)/norm(one(T))))
    eigfact!(copy_oftype(A, S))
end

eigfact!(A::SymTridiagonal{<:BlasReal}, irange::UnitRange) =
    Eigen(LAPACK.stegr!('V', 'I', A.dv, A.ev, 0.0, 0.0, irange.start, irange.stop)...)
function eigfact(A::SymTridiagonal{T}, irange::UnitRange) where T
    S = promote_type(Float32, typeof(zero(T)/norm(one(T))))
    return eigfact!(copy_oftype(A, S), irange)
end

eigfact!(A::SymTridiagonal{<:BlasReal}, vl::Real, vu::Real) =
    Eigen(LAPACK.stegr!('V', 'V', A.dv, A.ev, vl, vu, 0, 0)...)
function eigfact(A::SymTridiagonal{T}, vl::Real, vu::Real) where T
    S = promote_type(Float32, typeof(zero(T)/norm(one(T))))
    return eigfact!(copy_oftype(A, S), vl, vu)
end

eigvals!(A::SymTridiagonal{<:BlasReal}) = LAPACK.stev!('N', A.dv, A.ev)[1]
function eigvals(A::SymTridiagonal{T}) where T
    S = promote_type(Float32, typeof(zero(T)/norm(one(T))))
    return eigvals!(copy_oftype(A, S))
end

eigvals!(A::SymTridiagonal{<:BlasReal}, irange::UnitRange) =
    LAPACK.stegr!('N', 'I', A.dv, A.ev, 0.0, 0.0, irange.start, irange.stop)[1]
function eigvals(A::SymTridiagonal{T}, irange::UnitRange) where T
    S = promote_type(Float32, typeof(zero(T)/norm(one(T))))
    return eigvals!(copy_oftype(A, S), irange)
end

eigvals!(A::SymTridiagonal{<:BlasReal}, vl::Real, vu::Real) =
    LAPACK.stegr!('N', 'V', A.dv, A.ev, vl, vu, 0, 0)[1]
function eigvals(A::SymTridiagonal{T}, vl::Real, vu::Real) where T
    S = promote_type(Float32, typeof(zero(T)/norm(one(T))))
    return eigvals!(copy_oftype(A, S), vl, vu)
end

#Computes largest and smallest eigenvalue
eigmax(A::SymTridiagonal) = eigvals(A, size(A, 1):size(A, 1))[1]
eigmin(A::SymTridiagonal) = eigvals(A, 1:1)[1]

#Compute selected eigenvectors only corresponding to particular eigenvalues
eigvecs(A::SymTridiagonal) = eigfact(A)[:vectors]

"""
    eigvecs(A::SymTridiagonal[, eigvals]) -> Matrix

Returns a matrix `M` whose columns are the eigenvectors of `A`. (The `k`th eigenvector can
be obtained from the slice `M[:, k]`.)

If the optional vector of eigenvalues `eigvals` is specified, `eigvecs`
returns the specific corresponding eigenvectors.

# Examples
```jldoctest
julia> A = SymTridiagonal([1.; 2.; 1.], [2.; 3.])
3×3 SymTridiagonal{Float64}:
 1.0  2.0   ⋅
 2.0  2.0  3.0
  ⋅   3.0  1.0

julia> eigvals(A)
3-element Array{Float64,1}:
 -2.14005
  1.0
  5.14005

julia> eigvecs(A)
3×3 Array{Float64,2}:
  0.418304  -0.83205      0.364299
 -0.656749  -7.39009e-16  0.754109
  0.627457   0.5547       0.546448

julia> eigvecs(A, [1.])
3×1 Array{Float64,2}:
  0.83205
  4.26351e-17
 -0.5547
```
"""
eigvecs(A::SymTridiagonal{<:BlasFloat}, eigvals::Vector{<:Real}) = LAPACK.stein!(A.dv, A.ev, eigvals)

#tril and triu

istriu(M::SymTridiagonal) = iszero(M.ev)
istril(M::SymTridiagonal) = iszero(M.ev)

function tril!(M::SymTridiagonal, k::Integer=0)
    n = length(M.dv)
    if abs(k) > n
        throw(ArgumentError("requested diagonal, $k, out of bounds in matrix of size ($n,$n)"))
    elseif k < -1
        fill!(M.ev,0)
        fill!(M.dv,0)
        return Tridiagonal(M.ev,M.dv,copy(M.ev))
    elseif k == -1
        fill!(M.dv,0)
        return Tridiagonal(M.ev,M.dv,zeros(M.ev))
    elseif k == 0
        return Tridiagonal(M.ev,M.dv,zeros(M.ev))
    elseif k >= 1
        return Tridiagonal(M.ev,M.dv,copy(M.ev))
    end
end

function triu!(M::SymTridiagonal, k::Integer=0)
    n = length(M.dv)
    if abs(k) > n
        throw(ArgumentError("requested diagonal, $k, out of bounds in matrix of size ($n,$n)"))
    elseif k > 1
        fill!(M.ev,0)
        fill!(M.dv,0)
        return Tridiagonal(M.ev,M.dv,copy(M.ev))
    elseif k == 1
        fill!(M.dv,0)
        return Tridiagonal(zeros(M.ev),M.dv,M.ev)
    elseif k == 0
        return Tridiagonal(zeros(M.ev),M.dv,M.ev)
    elseif k <= -1
        return Tridiagonal(M.ev,M.dv,copy(M.ev))
    end
end

###################
# Generic methods #
###################

#Needed for inv_usmani()
mutable struct ZeroOffsetVector
    data::Vector
end
getindex(a::ZeroOffsetVector, i) = a.data[i+1]
setindex!(a::ZeroOffsetVector, x, i) = a.data[i+1]=x

## structured matrix methods ##
function Base.replace_in_print_matrix(A::SymTridiagonal, i::Integer, j::Integer, s::AbstractString)
    i==j-1||i==j||i==j+1 ? s : Base.replace_with_centered_mark(s)
end

#Implements the inverse using the recurrence relation between principal minors
# a, b, c are assumed to be the subdiagonal, diagonal, and superdiagonal of
# a tridiagonal matrix.
#Reference:
#    R. Usmani, "Inversion of a tridiagonal Jacobi matrix",
#    Linear Algebra and its Applications 212-213 (1994), pp.413-414
#    doi:10.1016/0024-3795(94)90414-6
function inv_usmani(a::V, b::V, c::V) where {T,V<:AbstractVector{T}}
    n = length(b)
    θ = ZeroOffsetVector(zeros(T, n+1)) #principal minors of A
    θ[0] = 1
    n>=1 && (θ[1] = b[1])
    for i=2:n
        θ[i] = b[i]*θ[i-1]-a[i-1]*c[i-1]*θ[i-2]
    end
    φ = zeros(T, n+1)
    φ[n+1] = 1
    n>=1 && (φ[n] = b[n])
    for i=n-1:-1:1
        φ[i] = b[i]*φ[i+1]-a[i]*c[i]*φ[i+2]
    end
    α = Matrix{T}(n, n)
    for i=1:n, j=1:n
        sign = (i+j)%2==0 ? (+) : (-)
        if i<j
            α[i,j]=(sign)(prod(c[i:j-1]))*θ[i-1]*φ[j+1]/θ[n]
        elseif i==j
            α[i,i]=                       θ[i-1]*φ[i+1]/θ[n]
        else #i>j
            α[i,j]=(sign)(prod(a[j:i-1]))*θ[j-1]*φ[i+1]/θ[n]
        end
    end
    α
end

#Implements the determinant using principal minors
#Inputs and reference are as above for inv_usmani()
function det_usmani(a::V, b::V, c::V) where {T,V<:AbstractVector{T}}
    n = length(b)
    θa = one(T)
    if n == 0
        return θa
    end
    θb = b[1]
    for i=2:n
        θb, θa = b[i]*θb-a[i-1]*c[i-1]*θa, θb
    end
    return θb
end

inv(A::SymTridiagonal) = inv_usmani(A.ev, A.dv, A.ev)
det(A::SymTridiagonal) = det_usmani(A.ev, A.dv, A.ev)

function getindex(A::SymTridiagonal{T}, i::Integer, j::Integer) where T
    if !(1 <= i <= size(A,2) && 1 <= j <= size(A,2))
        throw(BoundsError(A, (i,j)))
    end
    if i == j
        return A.dv[i]
    elseif i == j + 1
        return A.ev[j]
    elseif i + 1 == j
        return A.ev[i]
    else
        return zero(T)
    end
end

function setindex!(A::SymTridiagonal, x, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    if i == j
        @inbounds A.dv[i] = x
    else
        throw(ArgumentError("cannot set off-diagonal entry ($i, $j)"))
    end
    return x
end

## Tridiagonal matrices ##
struct Tridiagonal{T,V<:AbstractVector{T}} <: AbstractMatrix{T}
    dl::V    # sub-diagonal
    d::V     # diagonal
    du::V    # sup-diagonal
    du2::V   # supsup-diagonal for pivoting in LU
    function Tridiagonal{T}(dl::V, d::V, du::V) where {T,V<:AbstractVector{T}}
        n = length(d)
        if (length(dl) != n-1 || length(du) != n-1)
            throw(ArgumentError(string("cannot construct Tridiagonal from incompatible ",
                "lengths of subdiagonal, diagonal and superdiagonal: ",
                "($(length(dl)), $(length(d)), $(length(du)))")))
        end
        new{T,V}(dl, d, du)
    end
    # constructor used in lufact!
    function Tridiagonal{T,V}(dl::V, d::V, du::V, du2::V) where {T,V<:AbstractVector{T}}
        new{T,V}(dl, d, du, du2)
    end
end

"""
    Tridiagonal(dl::V, d::V, du::V) where V <: AbstractVector

Construct a tridiagonal matrix from the first subdiagonal, diagonal, and first superdiagonal,
respectively. The result is of type `Tridiagonal` and provides efficient specialized linear
solvers, but may be converted into a regular matrix with
[`convert(Array, _)`](@ref) (or `Array(_)` for short).
The lengths of `dl` and `du` must be one less than the length of `d`.

# Examples
```jldoctest
julia> dl = [1, 2, 3];

julia> du = [4, 5, 6];

julia> d = [7, 8, 9, 0];

julia> Tridiagonal(dl, d, du)
4×4 Tridiagonal{Int64,Array{Int64,1}}:
 7  4  ⋅  ⋅
 1  8  5  ⋅
 ⋅  2  9  6
 ⋅  ⋅  3  0
```
"""
Tridiagonal(dl::V, d::V, du::V) where {T,V<:AbstractVector{T}} = Tridiagonal{T}(dl, d, du)

"""
    Tridiagonal(A)

Construct a tridiagonal matrix from the first sub-diagonal,
diagonal and first super-diagonal of the matrix `A`.

# Examples
```jldoctest
julia> A = [1 2 3 4; 1 2 3 4; 1 2 3 4; 1 2 3 4]
4×4 Array{Int64,2}:
 1  2  3  4
 1  2  3  4
 1  2  3  4
 1  2  3  4

julia> Tridiagonal(A)
4×4 Tridiagonal{Int64,Array{Int64,1}}:
 1  2  ⋅  ⋅
 1  2  3  ⋅
 ⋅  2  3  4
 ⋅  ⋅  3  4
```
"""
Tridiagonal(A::AbstractMatrix) = Tridiagonal(diag(A,-1), diag(A,0), diag(A,1))

size(M::Tridiagonal) = (length(M.d), length(M.d))
function size(M::Tridiagonal, d::Integer)
    if d < 1
        throw(ArgumentError("dimension d must be ≥ 1, got $d"))
    elseif d <= 2
        return length(M.d)
    else
        return 1
    end
end

function convert(::Type{Matrix{T}}, M::Tridiagonal{T}) where T
    A = zeros(T, size(M))
    for i = 1:length(M.d)
        A[i,i] = M.d[i]
    end
    for i = 1:length(M.d)-1
        A[i+1,i] = M.dl[i]
        A[i,i+1] = M.du[i]
    end
    A
end
convert(::Type{Matrix}, M::Tridiagonal{T}) where {T} = convert(Matrix{T}, M)
convert(::Type{Array}, M::Tridiagonal) = convert(Matrix, M)
full(M::Tridiagonal) = convert(Array, M)
function similar(M::Tridiagonal, ::Type{T}) where T
    Tridiagonal{T}(similar(M.dl, T), similar(M.d, T), similar(M.du, T))
end

# Operations on Tridiagonal matrices
copy!(dest::Tridiagonal, src::Tridiagonal) = (copy!(dest.dl, src.dl); copy!(dest.d, src.d); copy!(dest.du, src.du); dest)

#Elementary operations
broadcast(::typeof(abs), M::Tridiagonal) = Tridiagonal(abs.(M.dl), abs.(M.d), abs.(M.du))
broadcast(::typeof(round), M::Tridiagonal) = Tridiagonal(round.(M.dl), round.(M.d), round.(M.du))
broadcast(::typeof(trunc), M::Tridiagonal) = Tridiagonal(trunc.(M.dl), trunc.(M.d), trunc.(M.du))
broadcast(::typeof(floor), M::Tridiagonal) = Tridiagonal(floor.(M.dl), floor.(M.d), floor.(M.du))
broadcast(::typeof(ceil), M::Tridiagonal) = Tridiagonal(ceil.(M.dl), ceil.(M.d), ceil.(M.du))
for func in (:conj, :copy, :real, :imag)
    @eval function ($func)(M::Tridiagonal)
        Tridiagonal(($func)(M.dl), ($func)(M.d), ($func)(M.du))
    end
end
broadcast(::typeof(round), ::Type{T}, M::Tridiagonal) where {T<:Integer} =
    Tridiagonal(round.(T, M.dl), round.(T, M.d), round.(T, M.du))
broadcast(::typeof(trunc), ::Type{T}, M::Tridiagonal) where {T<:Integer} =
    Tridiagonal(trunc.(T, M.dl), trunc.(T, M.d), trunc.(T, M.du))
broadcast(::typeof(floor), ::Type{T}, M::Tridiagonal) where {T<:Integer} =
    Tridiagonal(floor.(T, M.dl), floor.(T, M.d), floor.(T, M.du))
broadcast(::typeof(ceil), ::Type{T}, M::Tridiagonal) where {T<:Integer} =
    Tridiagonal(ceil.(T, M.dl), ceil.(T, M.d), ceil.(T, M.du))

transpose(M::Tridiagonal) = Tridiagonal(M.du, M.d, M.dl)
adjoint(M::Tridiagonal) = conj(transpose(M))

function diag(M::Tridiagonal{T}, n::Integer=0) where T
    if n == 0
        return M.d
    elseif n == -1
        return M.dl
    elseif n == 1
        return M.du
    elseif abs(n) < size(M,1)
        return zeros(T,size(M,1)-abs(n))
    else
        throw(ArgumentError("$n-th diagonal of a $(size(M)) matrix doesn't exist!"))
    end
end

function getindex(A::Tridiagonal{T}, i::Integer, j::Integer) where T
    if !(1 <= i <= size(A,2) && 1 <= j <= size(A,2))
        throw(BoundsError(A, (i,j)))
    end
    if i == j
        return A.d[i]
    elseif i == j + 1
        return A.dl[j]
    elseif i + 1 == j
        return A.du[i]
    else
        return zero(T)
    end
end

function setindex!(A::Tridiagonal, x, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    if i == j
        @inbounds A.d[i] = x
    elseif i - j == 1
        @inbounds A.dl[j] = x
    elseif j - i == 1
        @inbounds A.du[i] = x
    elseif !iszero(x)
        throw(ArgumentError(string("cannot set entry ($i, $j) off ",
            "the tridiagonal band to a nonzero value ($x)")))
    end
    return x
end

## structured matrix methods ##
function Base.replace_in_print_matrix(A::Tridiagonal,i::Integer,j::Integer,s::AbstractString)
    i==j-1||i==j||i==j+1 ? s : Base.replace_with_centered_mark(s)
end

#tril and triu

istriu(M::Tridiagonal) = iszero(M.dl)
istril(M::Tridiagonal) = iszero(M.du)

function tril!(M::Tridiagonal, k::Integer=0)
    n = length(M.d)
    if abs(k) > n
        throw(ArgumentError("requested diagonal, $k, out of bounds in matrix of size ($n,$n)"))
    elseif k < -1
        fill!(M.dl,0)
        fill!(M.d,0)
        fill!(M.du,0)
    elseif k == -1
        fill!(M.d,0)
        fill!(M.du,0)
    elseif k == 0
        fill!(M.du,0)
    end
    return M
end

function triu!(M::Tridiagonal, k::Integer=0)
    n = length(M.d)
    if abs(k) > n
        throw(ArgumentError("requested diagonal, $k, out of bounds in matrix of size ($n,$n)"))
    elseif k > 1
        fill!(M.dl,0)
        fill!(M.d,0)
        fill!(M.du,0)
    elseif k == 1
        fill!(M.dl,0)
        fill!(M.d,0)
    elseif k == 0
        fill!(M.dl,0)
    end
    return M
end

###################
# Generic methods #
###################

+(A::Tridiagonal, B::Tridiagonal) = Tridiagonal(A.dl+B.dl, A.d+B.d, A.du+B.du)
-(A::Tridiagonal, B::Tridiagonal) = Tridiagonal(A.dl-B.dl, A.d-B.d, A.du-B.du)
*(A::Tridiagonal, B::Number) = Tridiagonal(A.dl*B, A.d*B, A.du*B)
*(B::Number, A::Tridiagonal) = A*B
/(A::Tridiagonal, B::Number) = Tridiagonal(A.dl/B, A.d/B, A.du/B)

==(A::Tridiagonal, B::Tridiagonal) = (A.dl==B.dl) && (A.d==B.d) && (A.du==B.du)
==(A::Tridiagonal, B::SymTridiagonal) = (A.dl==A.du==B.ev) && (A.d==B.dv)
==(A::SymTridiagonal, B::Tridiagonal) = (B.dl==B.du==A.ev) && (B.d==A.dv)

inv(A::Tridiagonal) = inv_usmani(A.dl, A.d, A.du)
det(A::Tridiagonal) = det_usmani(A.dl, A.d, A.du)

convert(::Type{Tridiagonal{T}},M::Tridiagonal) where {T} =
    Tridiagonal(convert(AbstractVector{T}, M.dl), convert(AbstractVector{T}, M.d), convert(AbstractVector{T}, M.du))
convert(::Type{AbstractMatrix{T}},M::Tridiagonal) where {T} = convert(Tridiagonal{T}, M)
convert(::Type{Tridiagonal{T}}, M::SymTridiagonal{T}) where {T} = Tridiagonal(M)
function convert(::Type{SymTridiagonal{T}}, M::Tridiagonal) where T
    if M.dl == M.du
        return SymTridiagonal{T}(convert(AbstractVector{T},M.d), convert(AbstractVector{T},M.dl))
    else
        throw(ArgumentError("Tridiagonal is not symmetric, cannot convert to SymTridiagonal"))
    end
end

max_stars_repo_name: visr/julia

max_issues_repo_name: visr/julia

max_forks_repo_name: visr/julia


text:
"""
    run1simulation(input::BranchingMoranInput[, rng::AbstractRNG])

    Simulate a population of cells that grows by a branching process to fixed size then
    switches to a Moran process.

    Take simulation parameters from `input` and return a Simulation.
"""
function run1simulation(input::BranchingMoranInput, rng::AbstractRNG = Random.GLOBAL_RNG)

    #Run branching simulation starting with a single cell.
    #Initially set clonalmutations = 0 and μ = 1. These are expanded later.
    moduletracker =
        branchingprocess(input.b, input.d, input.Nmax, 1, rng, numclones = input.numclones,
            fixedmu = true, clonalmutations = 0, selection = input.selection,
            tevent = input.tevent, maxclonesize = Inf)

    moduletracker =
        moranprocess!(moduletracker, input.bdrate, input.tmax, 1, rng::AbstractRNG;
        numclones = input.numclones, fixedmu = true, selection = input.selection,
        tevent = input.tevent)

    moduletracker =
        processresults!(moduletracker, input.μ, input.clonalmutations, rng, fixedmu=input.fixedmu)

    return Simulation(input, moduletracker)
end

"""
    run1simulation(input::BranchingInput[, rng::AbstractRNG])

    Simulate a population of cells that grows by a branching process to a fixed size.
"""
function run1simulation(input::BranchingInput, rng::AbstractRNG = Random.GLOBAL_RNG)

    #Run branching simulation starting with a single cell.
    #Initially set clonalmutations = 0 and μ = 1. These are expanded later.
    moduletracker =
        branchingprocess(input.b, input.d, input.Nmax, 1, rng, numclones = input.numclones,
            fixedmu = true, clonalmutations = 0, selection = input.selection,
            tevent = input.tevent, maxclonesize = input.maxclonesize)

    #Add mutations and process simulation output to get SimResults.
    #Remove undetectable subclones from moduletracker
    moduletracker =
        processresults!(moduletracker, input.μ, input.clonalmutations, rng, fixedmu=input.fixedmu)

    return Simulation(input, moduletracker)
end

"""
    run1simulation(input::MoranInput[, rng::AbstractRNG])

    Simulate a population of cells according to a Moran process for fixed time.
"""

function run1simulation(input::MoranInput, rng::AbstractRNG = Random.GLOBAL_RNG)

    #Run Moran simulation starting from a population of N identical cells.
    #Initially set clonalmutations = 0 and μ = 1. These are expanded later.
    moduletracker =
        moranprocess(input.N, input.bdrate, input.tmax, 1, rng,
                    numclones = input.numclones, fixedmu = true,
                    clonalmutations = 0, selection = input.selection,
                    tevent = input.tevent)

    #Add mutations and process simulation output to get SimResults.
    #Remove undetectable subclones from moduletracker
    moduletracker =
        processresults!(moduletracker, input.μ, input.clonalmutations, rng, fixedmu=input.fixedmu)
    return Simulation(input,moduletracker)
end

"""
    branchingprocess(input::BranchingInput, rng::AbstractRNG; <keyword arguments>)

Simulate a stochastic branching process with parameters defined by input and return
ModuleTracker.

Simulation is by a rejection-kinetic Monte Carlo algorithm and starts with a single cell.

"""
function branchingprocess(input::BranchingInput, rng::AbstractRNG,
                            fixedmu=input.fixedmu, μ=input.μ,
                            clonalmutations=input.clonalmutations)

    return branchingprocess(input.b, input.d, input.Nmax, μ, rng,
                            numclones = input.numclones, fixedmu = fixedmu,
                            clonalmutations = clonalmutations,
                            selection = input.selection, tevent = input.tevent,
                            maxclonesize = input.maxclonesize)
end

"""
    branchingprocess(b, d, Nmax, μ, rng::AbstractRNG; <keyword arguments>)

Simulate a stochastic branching process, starting with a single cell, with with birth rate
`b`, death rate `d` until population reaches size `Nmax`.

Cells accumulate neutral mutations at division with rate `μ`, until all subclones exceed
`maxclonesize`.

If `numclones` = 0, all cells have the same fitness and there is only one (sub)clone.
Otherwise, `numclones` is the number of fit subclones. The `i`th subclone arises by a single
cell mutating at time `tevent[i]` and has selection coefficient `selection[i]`.

"""
function branchingprocess(b, d, Nmax, μ, rng::AbstractRNG; numclones=0, fixedmu=false,
    clonalmutations=μ, selection=Float64[], tevent=Float64[], maxclonesize=200)

    #initialize arrays and parameters
    moduletracker = initializesim_branching(
        Nmax,
        clonalmutations=clonalmutations
    )

    #run simulation
    moduletracker = branchingprocess!(moduletracker, b, d, Nmax, μ, rng, numclones=numclones,
        fixedmu=fixedmu, selection=selection, tevent=tevent, maxclonesize=maxclonesize)
    return moduletracker
end

"""
    branchingprocess!(moduletracker::ModuleTracker, b, d, Nmax, μ, rng::AbstractRNG;
        <keyword arguments>)

Run branching process simulation, starting in state defined by moduletracker.

See also [`branchingprocess`](@ref)
"""
function branchingprocess!(moduletracker::ModuleTracker, b, d, Nmax, μ, rng::AbstractRNG;
    numclones=0, fixedmu=false, selection=Float64[], tevent=Float64[],
    maxclonesize=200, tmax=Inf)

    t, N = moduletracker.tvec[end], moduletracker.Nvec[end]
    mutID = N == 1 ? 1 : getmutID(moduletracker.cells)

    nclonescurrent = length(moduletracker.subclones) + 1
    executed = false
    changemutrate = BitArray(undef, numclones + 1)
    changemutrate .= 1

    birthrates, deathrates = set_branching_birthdeath_rates(b, d, selection)

    #Rmax starts with b + d and changes once a fitter mutant is introduced, this ensures
    #that b and d have correct units
    Rmax = (maximum(birthrates[1:nclonescurrent])
                                + maximum(deathrates[1:nclonescurrent]))

    while N < Nmax

        #calc next event time and break if it exceeds tmax
        Δt =  1/(Rmax * N) .* exptime(rng)
        t = t + Δt
        if t > tmax
            break
        end

        randcell = rand(rng,1:N) #pick a random cell
        r = rand(rng,Uniform(0,Rmax))
        #get birth and death rates for randcell
        br = birthrates[moduletracker.cells[randcell].clonetype]
        dr = deathrates[moduletracker.cells[randcell].clonetype]

        if r < br
            #cell divides
            moduletracker, mutID = celldivision!(moduletracker, randcell, mutID, μ, rng,
                fixedmu = fixedmu)
            N += 1
            #check if t>=tevent for next fit subclone
            if nclonescurrent < numclones + 1 && t >= tevent[nclonescurrent]
                #if current number clones != final number clones, one of the new cells is
                #mutated to become fitter and form a new clone
                    moduletracker, nclonescurrent = cellmutation!(moduletracker, N, N, t,
                        nclonescurrent)

                    #change Rmax now there is a new fitter mutant
                    Rmax = (maximum(birthrates[1:nclonescurrent])
                                + maximum(deathrates[1:nclonescurrent]))
            end
            moduletracker = update_time_popsize!(moduletracker, t, N)

        elseif r < br + dr
            #cell dies
            moduletracker = celldeath!(moduletracker, randcell)
            N -= 1
            moduletracker = update_time_popsize!(moduletracker, t, N)
            #return empty moduletracker if all cells have died
            if N == 0
                return moduletracker
            end
        end

        #if population of all clones is sufficiently large no new mutations
        #are acquired, can use this approximation as only mutations above 1%
        #frequency can be reliably detected
        if ((maxclonesize !== nothing &&
            executed == false &&
            (getclonesize(moduletracker) .> maxclonesize) == changemutrate))

            μ = 0
            fixedmu = true
            executed = true
        end
    end
    return moduletracker
end

"""
    getclonesize(moduletracker::ModuleTracker)

Return number of cells in each subclone (including wild-type).
"""
function getclonesize(moduletracker::ModuleTracker)
    return getclonesize(moduletracker.Nvec[end], moduletracker.subclones)
end

"""
    getclonesize(N, subclones)
"""
function getclonesize(N, subclones)
   sizevec = [clone.size for clone in subclones]
   prepend!(sizevec, N - sum(sizevec))
end

"""
    update_time_popsize(moduletracker::ModuleTracker, t, N)

Update moduletracker with new time `t` and pop size `N`.
"""
function update_time_popsize!(moduletracker::ModuleTracker, t, N)
    push!(moduletracker.tvec,t)
    push!(moduletracker.Nvec, N)
    return moduletracker
end

"""
    set_branching_birthdeath_rates(b, d, selection)

Return Vectors of birthrates and deathrates for each subclone (including wild-type).
"""
function set_branching_birthdeath_rates(b, d, selection)
    birthrates = [b]
    deathrates = [d]
    #add birth and death rates for each subclone.
    for i in 1:length(selection)
        push!(deathrates, d)
        push!(birthrates, (1 + selection[i]) .* b)
    end
    return birthrates,deathrates
end

"""
    moranprocess(input::MoranInput, rng::AbstractRNG; <keyword arguments>)

Simulate a Moran process with parameters defined by input and return ModuleTracker.

Simulation is by a Gillespie algorithm.

"""
function moranprocess(input::MoranInput, rng::AbstractRNG,
                    fixedmu=input.fixedmu, μ=input.μ,
                    clonalmutations=input.clonalmutations)

    return moranprocess(input.N, input.bdrate, input.tmax, μ, rng,
                        numclones = input.numclones, fixedmu = fixedmu,
                        clonalmutations = clonalmutations,
                        selection = input.selection, tevent = input.tevent)
end

"""
    moranprocess(N, bdrate, tmax, μ, rng::AbstractRNG; <keyword arguments>)

Simulate a Moran process starting with `N` cells until time `tmax`.

Update events comprise of a birth and a death, and occur with rate `bdrate`. Cells
accumulate neutral mutations at division with rate `μ`, until all subclones exceed
`maxclonesize`.

If `numclones` = 0, all cells have the same fitness and there is only one (sub)clone.
Otherwise, `numclones` is the number of fit subclones. The `i`th subclone arises by a single
cell mutating at time `tevent[i]` and has selection coefficient `selection[i]`.

"""
function moranprocess(N, bdrate, tmax, μ, rng::AbstractRNG; numclones = 0, fixedmu = false,
    clonalmutations = μ, selection = Float64[], tevent = Float64[])

    moduletracker = initializesim_moran(N, clonalmutations = clonalmutations)

    #run simulation
    moduletracker = moranprocess!(moduletracker, bdrate, tmax, μ, rng, numclones = numclones,
    fixedmu = fixedmu, selection = selection, tevent = tevent)
    return moduletracker
end
"""
    moranprocess!(moduletracker::ModuleTracker, bdrate, tmax, μ, rng::AbstractRNG;
        <keyword-arguments>)

Run Moran process simulation, starting in state defined by moduletracker.

See also [`moranprocess`](@ref)

"""
function moranprocess!(moduletracker::ModuleTracker, bdrate, tmax, μ, rng::AbstractRNG;
    numclones = 0, fixedmu = false, selection = Float64[], tevent = Float64[])

    t, N = moduletracker.tvec[end], moduletracker.Nvec[end]
    mutID = getmutID(moduletracker.cells)

    nclonescurrent = length(moduletracker.subclones) + 1

    while true

        #calc next event time and break if it exceeds tmax
        Δt =  1/(bdrate*N) .* exptime(rng)
        t = t + Δt
        if t > tmax
            break
        end

        #pick a cell to divide proportional to clonetype
        if nclonescurrent == 1
            randcelldivide = rand(rng, 1:N)
        else
            p = [cell.clonetype==1 ? 1 : 1 + selection[cell.clonetype - 1]
                    for cell in moduletracker.cells]
            p /= sum(p)
            randcelldivide = sample(rng, 1:N, ProbabilityWeights(p))
        end

        #pick a random cell to die
        randcelldie = rand(rng,1:N)

        #cell divides
        moduletracker, mutID = celldivision!(moduletracker, randcelldivide, mutID, μ, rng,
            fixedmu = fixedmu)

        #check if t>=tevent for next fit subclone and more subclones are expected
        if nclonescurrent < numclones + 1 && t >= tevent[nclonescurrent]
            moduletracker, nclonescurrent = cellmutation!(moduletracker, N+1, N, t, nclonescurrent)
        end

        #cell dies
        moduletracker = celldeath!(moduletracker, randcelldie)

        moduletracker = update_time_popsize!(moduletracker, t, N)

    end
    return moduletracker
end

"""
    initializesim(input::BranchingInput, rng::AbstractRNG=Random.GLOBAL_RNG)

Initialise simulation and return a ModuleTracker.
"""
function initializesim(siminput::BranchingInput, rng::AbstractRNG=Random.GLOBAL_RNG)

    return initializesim_branching(
        siminput.Nmax,
        clonalmutations=siminput.clonalmutations,
    )
end

"""
    initializesim(input::MoranInput, rng::AbstractRNG=Random.GLOBAL_RNG)
"""
function initializesim(siminput::MoranInput, rng::AbstractRNG=Random.GLOBAL_RNG)

    return initializesim_moran(
        siminput.N,
        clonalmutations=siminput.clonalmutations,
    )
end

"""
    initializesim_branching(input::BranchingInput, rng::AbstractRNG=Random.GLOBAL_RNG)

Initialise simulation and return a ModuleTracker.
"""
function initializesim_branching(Nmax=nothing; clonalmutations=0)

    #initialize time to zero
    t = 0.0
    tvec = Float64[]
    push!(tvec,t)

    #population starts with one cell
    N = 1
    Nvec = Int64[]
    push!(Nvec,N)

    #Initialize array of cell type that stores mutations for each cell and their clone type
    #clone type of 1 is the host population with selection=0
    cells = Cell[]
    if Nmax !== nothing
        sizehint!(cells, Nmax)
    end
    push!(cells, Cell([],1))

    #need to keep track of mutations, assuming infinite sites, new mutations will be unique,
    #we assign each new muation with a unique integer by simply counting up from one
    mutID = 1
    cells[1],mutID = newmutations!(cells[1], clonalmutations, mutID)

    subclones = CloneTracker[]

    moduletracker = ModuleTracker(
        Nvec,
        tvec,
        cells,
        subclones,
        1,
        0
    )
    return moduletracker
end

function initializesim_moran(N; clonalmutations=0)

    #initialize time to zero
    t = 0.0
    tvec = Float64[]
    push!(tvec,t)

    #Initialize array of cell type that stores mutations for each cell and their clone type
    #clone type of 1 is the host population with selection=0
    cells = [Cell(Int64[],1) for _ in 1:N]

    #need to keep track of mutations, assuming infinite sites, new mutations will be unique,
    #we assign each new muation with a unique integer by simply counting up from one
    for cell in cells
        cell.mutations = collect(1:clonalmutations)
    end

    subclones = CloneTracker[]

    moduletracker = ModuleTracker(
        [N],
        tvec,
        cells,
        subclones,
        1,
        0
    )
    return moduletracker
end

function initializesim_from_cells(cells::Array{Cell,1}, subclones::Array{CloneTracker, 1},
    id, parentid;
    inittime=0.0)

    #initialize time to zero
    t = inittime
    tvec = Float64[]
    push!(tvec,t)

    #population starts from list of cells
    N = length(cells)
    Nvec = Int64[]
    push!(Nvec, N)

    for (i, subclone) in enumerate(subclones)
        subclone.size = sum(cell.clonetype .== 1)
    end

    moduletracker = ModuleTracker(
        Nvec,
        tvec,
        cells,
        subclones,
        id,
        parentid
    )
    return moduletracker
end

function newmutations!(cell, μ, mutID, rng::AbstractRNG; fixedmu = false)
    #function to add new mutations to cells based on μ
    numbermutations = fixedmu ? μ : rand(rng,Poisson(μ))
    return newmutations!(cell, numbermutations, mutID)
end

function newmutations!(cell, μ, mutID)
    #function to add new mutations to cells based on μ
    newmutations = mutID:mutID + μ - 1
    append!(cell.mutations, newmutations)
    mutID = mutID + μ
    return cell, mutID
end

function celldivision!(moduletracker::ModuleTracker, parentcell, mutID, μ,
    rng::AbstractRNG; fixedmu=true)

    push!(moduletracker.cells, copycell(moduletracker.cells[parentcell])) #add new copy of parent cell to cells
    #add new mutations to both new cells
    if μ > 0.0
        moduletracker.cells[parentcell],mutID = newmutations!(moduletracker.cells[parentcell], μ,
            mutID, rng, fixedmu = fixedmu)
        moduletracker.cells[end],mutID = newmutations!(moduletracker.cells[end], μ, mutID, rng,
            fixedmu=fixedmu)
    end
    clonetype = moduletracker.cells[parentcell].clonetype
    if clonetype > 1
        moduletracker.subclones[clonetype - 1].size += 1
    end

    return moduletracker, mutID
end

function cellmutation!(moduletracker::ModuleTracker, mutatingcell, N, t, nclonescurrent)

    #add new clone
    parenttype = moduletracker.cells[mutatingcell].clonetype
    mutations = deepcopy(moduletracker.cells[mutatingcell].mutations)
    Ndivisions = length(moduletracker.cells[mutatingcell].mutations)
    avdivisions = mean(map(x -> length(x.mutations), moduletracker.cells))
    clone = CloneTracker(parenttype, moduletracker.id, t, mutations, N, Ndivisions,
        avdivisions, 1)
    push!(moduletracker.subclones, clone)

    #change clone type of new cell and update clone sizes
    nclonescurrent += 1
    moduletracker.cells[mutatingcell].clonetype = nclonescurrent

    if parenttype > 1
        moduletracker.subclones[parenttype - 1].size -= 1
    end

    return moduletracker, nclonescurrent
end

function celldeath!(moduletracker::ModuleTracker, deadcell::Int64)
    #frequency of cell type decreases
    clonetype = moduletracker.cells[deadcell].clonetype
    if clonetype > 1
        moduletracker.subclones[clonetype - 1].size -= 1
    end
    #remove deleted cell
    deleteat!(moduletracker.cells, deadcell)

    return moduletracker
end

function celldeath!(moduletracker::ModuleTracker, deadcells::Array{Int64, 1})
    for deadcell in deadcells
        clonetype = moduletracker.cells[deadcell].clonetype
        if clonetype > 1
            moduletracker.subclones[clonetype - 1].size -= 1
        end
    end
    deleteat!(moduletracker.cells, sort(deadcells))

    return moduletracker
end

function getmutID(cells::Vector{Cell})
    if all(no_mutations.(cells))
        return 1
    else
        allmutations = reduce(vcat, [cell.mutations for cell in cells])
        return maximum(allmutations)
    end
end

function copycell(cellold::Cell)
    return Cell(copy(cellold.mutations), cellold.clonetype)
  end

function exptime(rng::AbstractRNG)
    rand(rng, Exponential(1))
end

function exptime(rng::AbstractRNG, λ)
    rand(rng, Exponential(1/λ))
end

function no_mutations(cell)
    return length(cell.mutations) == 0
end


text:
#***************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
# Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#***************************************************************************

function readMDP(filename, mode)

# natom, lattice, pbc, energy, stress, id, type, group, move, atomic number, atomic mass, charge, position, velocity, force

if mode==0 # binary
    frame = reinterpret(Float64,read(filename));
else # text
    tmp = Main.DelimitedFiles.readdlm(filename, ',');
    tmp = String.(reduce(vcat, split.(tmp))); # convert to 1D array of strings
    frame = parse.(Float64, tmp);  # convert to 1D array of Float64
end

config = initializeconfig(0);

natomall = 0;
nframe = 0;
n2 = 0;
n = length(frame);
while (n2 < n)
    nframe = nframe + 1;

    # natom, lattice, pbc, energy, stress, id, type, group, move, atomic number, atomic mass, charge, position, velocity, force
    natom = Int64(frame[n2+1]);  # number of atoms
    ncl = Int64(frame[n2+2]);    # lattice
    ncp = Int64(frame[n2+3]);    # periodic boundary conditions
    nce = Int64(frame[n2+4]);    # energy
    ncs = Int64(frame[n2+5]);    # stress
    nci = Int64(frame[n2+6]);    # atom id
    nct = Int64(frame[n2+7]);    # atom type
    ncg = Int64(frame[n2+8]);    # atom groups
    nco = Int64(frame[n2+9]);    # atom moved/fixed flags
    ncz = Int64(frame[n2+10]);   # atomic numbers
    ncm = Int64(frame[n2+11]);   # atomic masses
    ncq = Int64(frame[n2+12]);   # atom charges
    ncx = Int64(frame[n2+13]);   # positions
    ncv = Int64(frame[n2+14]);   # atom velocities
    ncf = Int64(frame[n2+15]);   # atom forces
    n2 = n2+15;

    if nframe==1
        config.ncx = ncx;
        config.ncv = ncv;
        config.ncq = ncq;
        config.nce = nce;
        config.ncf = ncf;
        config.ncs = ncs;
        config.nci = nci;
        config.nct = nct;
        config.ncg = ncg;
        config.ncz = ncz;
        config.ncm = ncm;
        config.nco = nco;
        config.ncl = ncl;
        config.ncp = ncp;
        config.nconfigs = 1;
        config.dim = ncx;
    end

    config.natom = hcat(config.natom, reshape([natom],1,1))

    # lattice, pbc, energy, stress
    n1 = n2 + 1; n2 = n2 + ncl;
    if nframe == 1 & ncl>0
        config.lattice = reshape(frame[n1:n2], (ncl, 1));
    elseif ncl>0
        config.lattice = hcat(config.lattice, reshape(frame[n1:n2], (ncl, 1)));
    end

    n1 = n2 + 1; n2 = n2 + ncp;
    if nframe == 1 & ncp>0
        config.pbc = reshape(frame[n1:n2], (ncp, 1));
    elseif ncp>0
        config.pbc = hcat(config.pbc, reshape(frame[n1:n2], (ncp, 1)));
    end

    n1 = n2 + 1; n2 = n2 + nce;
    if nframe == 1 & nce>0
        config.e = reshape(frame[n1:n2], (nce, 1));
    elseif nce>0
        config.e = hcat(config.e, reshape(frame[n1:n2], (nce, 1)));
    end

    n1 = n2 + 1; n2 = n2 + ncs;
    if nframe == 1 & ncs>0
        config.stress = reshape(frame[n1:n2], (ncs, 1));
    elseif ncs>0
        config.stress = hcat(config.stress, reshape(frame[n1:n2], (ncs, 1)));
    end

    k = nci + nct + ncg + nco + ncz + ncm + ncq + ncx + ncv + ncf;
    n1 = n2 + 1; n2 = n2 + natom*k;
    tmp = reshape(frame[n1:n2], (k, natom));

    if nframe == 1 & nci>0
        config.tags = tmp[1:nci,:];
    elseif nci>0
        config.tags = hcat(config.tags, tmp[1:nci,:]);
    end

    if nframe == 1 & nct>0
        config.t = tmp[(1+nci):(nci+nct),:];
    elseif nct>0
        config.t = hcat(config.t, tmp[(1+nci):(nci+nct),:]);
    end

    if nframe == 1 & ncg>0
        config.group = tmp[(1+nci+nct):(nci+nct+ncg),:];
    elseif ncg>0
        config.group = hcat(config.group, tmp[(1+nci+nct):(nci+nct+ncg),:]);
    end

    if nframe == 1 & nco>0
        config.move = tmp[(1+nci+nct+ncg):(nci+nct+ncg+nco),:];
    elseif nco>0
        config.move = hcat(config.move, tmp[(1+nci+nct+ncg):(nci+nct+ncg+nco),:]);
    end

    if nframe == 1 & ncz>0
        config.Z = tmp[(1+nci+nct+ncg+nco):(nci+nct+ncg+nco+ncz),:];
    elseif ncz>0
        config.Z = hcat(config.Z, tmp[(1+nci+nct+ncg+nco):(nci+nct+ncg+nco+ncz),:]);
    end

    if nframe == 1 & ncm>0
        config.mass = tmp[(1+nci+nct+ncg+nco+ncz):(nci+nct+ncg+nco+ncz+ncm),:];
    elseif ncm>0
        config.mass = hcat(config.mass, tmp[(1+nci+nct+ncg+nco+ncz):(nci+nct+ncg+nco+ncz+ncm),:]);
    end

    if nframe == 1 & ncq>0
        config.q = tmp[(1+nci+nct+ncg+nco+ncz+ncm):(nci+nct+ncg+nco+ncz+ncm+ncq),:];
    elseif ncq>0
        config.q = hcat(config.q, tmp[(1+nci+nct+ncg+nco+ncz+ncm):(nci+nct+ncg+nco+ncz+ncm+ncq),:]);
    end

    if nframe == 1 & ncx>0
        config.x = tmp[(1+nci+nct+ncg+nco+ncz+ncm+ncq):(nci+nct+ncg+nco+ncz+ncm+ncq+ncx),:];
    elseif ncx>0
        config.x = hcat(config.x, tmp[(1+nci+nct+ncg+nco+ncz+ncm+ncq):(nci+nct+ncg+nco+ncz+ncm+ncq+ncx),:]);
    end

    if nframe == 1 & ncv>0
        config.v = tmp[(1+nci+nct+ncg+nco+ncz+ncm+ncq+ncx):(nci+nct+ncg+nco+ncz+ncm+ncq+ncx+ncv),:];
    elseif ncv>0
        config.v = hcat(config.v, tmp[(1+nci+nct+ncg+nco+ncz+ncm+ncq+ncx):(nci+nct+ncg+nco+ncz+ncm+ncq+ncx+ncv),:]);
    end

    if nframe == 1 & ncf>0
        config.f = tmp[(1+nci+nct+ncg+nco+ncz+ncm+ncq+ncx+ncv):(nci+nct+ncg+nco+ncz+ncm+ncq+ncx+ncv+ncf),:];
    elseif ncf>0
        config.f = hcat(config.f, tmp[(1+nci+nct+ncg+nco+ncz+ncm+ncq+ncx+ncv):(nci+nct+ncg+nco+ncz+ncm+ncq+ncx+ncv+ncf),:]);
    end

    natomall = natomall + natom;
end
config.natomall = natomall;

return config

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# code to compile the Stock Market Execss and the State of the Union Address texts
# into a matching DataFrame and Document-Term matrix

using Lasso, CSV, DataDeps, DataFrames, FamaFrenchData, TextAnalysis, SparseArrays

register(DataDep("sotu",
        """
        Dataset: State of the Union Adresses
        Website: http://stateoftheunion.onetwothree.net
        Author: Brad Borevitz
        """,
        "http://stateoftheunion.onetwothree.net/texts/stateoftheunion1790-2020.txt.zip",
        "f61dfaf9b487872173e664488e6daad1783fa42691ed8009b2c5245866939b97",
        post_fetch_method = unpack
    ))

# fast way to remove sparse terms using the dtm, returns a new dtm with only terms that occur at least minoccurances times in the corpus
function remove_sparse_terms(D::TextAnalysis.DocumentTermMatrix, minoccurances = 3)
  m,n=size(D.dtm)
  rows = rowvals(D.dtm)
  vals = nonzeros(D.dtm)
  newrows = Int[]
  newcols = Int[]
  newvals = Int[]
  newterms = String[]
  new_column_indices = Dict{String, Int}()
  colsums = sum(D.dtm, dims=1)
  col = 1
  for i = 1:n
    if colsums[i]>=minoccurances
      for j in nzrange(D.dtm, i)
        row = rows[j]
        val = vals[j]
        push!(newrows, row)
        push!(newcols, col)
        push!(newvals, val)
      end
      term = D.terms[i]
      push!(newterms, term)
      new_column_indices[term] = col
      col += 1
    end
  end
  newterms
  if length(newrows) > 0
      newdtm = sparse(newrows, newcols, newvals)
  else
      newdtm = spzeros(Int,m, 0)
  end
  DocumentTermMatrix(newdtm, newterms, new_column_indices)
end

# get the SOTU text
function loadsotu()
    sotutxt = datadep"sotu/stateoftheunion1790-2020.txt"
    sotustr = read(sotutxt, String)
    addresses = strip.(split(sotustr, "***", keepempty=false)[2:end-1])

    # typical structure:
    #
    # ***
    #
    # State of the Union Address
    # Harry S. Truman
    # January 6, 1947
    #
    # Mr. President, Mr. Speaker, Members of the Congress of the United States:
    #
    # It looks like a good many of you have moved over to the left since I was
    # here last! ...

    # vectors to save parsed documents and dates
    docs = NGramDocument[]
    years = Int[]

    for a in addresses[1:end]
        # an IOBuffer provides an efficient way of reading lines
        strbuf = IOBuffer(a)

        # read header
        title = readline(strbuf)
        speaker = readline(strbuf)
        datestr = readline(strbuf)
        year = parse(Int, datestr[end-3:end])
        if 1927 <= year <= 2019

            # skip "Dear Congress" at start
            readline(strbuf)
            readline(strbuf)
            readline(strbuf)

            speechtext = read(strbuf, String)
            # remove everything after "NOTE:"
            ix = findfirst("NOTE:", speechtext)
            if !isnothing(ix)
                speechtext = speechtext[1:ix[1]-1]
            end

            # initial preprocessing
            sdoc = StringDocument(speechtext)
            prepare!(sdoc, strip_non_letters | strip_case)
            prepare!(sdoc, strip_stopwords)
            prepare!(sdoc, strip_whitespace | strip_corrupt_utf8)

            # build TextAnalysis NGramDocument of bigrams (complexity=2)
            doc = NGramDocument(text(sdoc), 2)
            author!(doc, speaker)
            timestamp!(doc, datestr)

            # save parsed documents into vectors
            if length(years) > 0 && year == years[end]
                # overwrite last one added becasue sometimes (e.g. 1953) the outgoing president addresses too
                docs[end] = doc
            else
                push!(docs, doc)
                push!(years, year)
            end
        end
    end

    corpus = Corpus(docs)
    @info "built SOTU corpus with $(length(corpus)) speeches"

    update_lexicon!(corpus)
    D = DocumentTermMatrix(corpus)
    D = remove_sparse_terms(D, 20)
    @info "built Document Term Matrix with bigrams occuring at least 20 times in the corpus"

    sotu = DataFrame(
            Date=[parse(Int, ts[end-3:end]) for ts in timestamps(corpus)],
            President=author.(corpus)
        )

    sotu, D
end

# load fama french annual market excess Returns
function loadreturns()
    tables, tablenotes, filenotes = readFamaFrench("F-F_Research_Data_Factors")
    FF3_annual = tables[2]
    rename!(FF3_annual, "Mkt-RF"=>"Rem")
    FF3_annual
end

# load merged sotu + market excess returns dataset
function loadsoturet(;cache=true)
    csvfile = joinpath(datadep"sotu", "sotu.ret.csv")
    dtmfile = joinpath(datadep"sotu", "sotu.dtm.csv")
    termsfile = joinpath(datadep"sotu", "sotu.terms.csv")

    if cache && isfile(csvfile) && isfile(dtmfile) && isfile(termsfile)
        df = CSV.read(csvfile, DataFrame)

        countsIJV = CSV.read(dtmfile, DataFrame)
        counts = sparse(countsIJV.I, countsIJV.J, countsIJV.V)

        terms = CSV.read(termsfile, DataFrame).term
        @info "SOTU + Rem loaded from cached $(datadep"sotu")"
    else
        sotu, D = loadsotu()
        ret = loadreturns()

        # merged dataset
        ixret = findall(in(sotu.Date), ret.Date)
        @assert sotu.Date == ret.Date[ixret] "sotu and return years do not lineup as expected"

        df = innerjoin(ret[ixret, ["Date", "Rem"]], sotu, on="Date")

        # save to csv file
        df |> CSV.write(csvfile)

        # save dtm in IJV sparse format
        counts = dtm(D)
        I, J, V = findnz(counts)
        DataFrame(I=I, J=J, V=V) |> CSV.write(dtmfile)

        # save terms corresponding to the columns of the counts matrix
        terms = D.terms
        DataFrame(index=1:length(terms), term=terms) |> CSV.write(termsfile)

        @info "SOTU + Rem loaded and saved to $(datadep"sotu")"
    end

    df, counts, terms
end

"""
converts a python/scipy sprase matrix to julia
Code from https://github.com/JuliaPy/PyCall.jl/issues/735#issuecomment-577336332
"""
function scipyCSC_to_julia(A)
    m, n = A.shape
    colPtr = Int[i+1 for i in PyArray(A."indptr")]
    rowVal = Int[i+1 for i in PyArray(A."indices")]
    nzVal = Vector{Float64}(PyArray(A."data"))
    B = SparseMatrixCSC{Float64,Int}(m, n, colPtr, rowVal, nzVal)
    return PyCall.pyjlwrap_new(B)
end

"""
converts a python/pandas dataframe to julia
Code from https://discourse.julialang.org/t/converting-pandas-dataframe-returned-from-pycall-to-julia-dataframe/43001/3
"""
function pd_to_df(df_pd)
    colnames = map(Symbol, df_pd.columns)
    df = DataFrame(Any[Array(df_pd[c].values) for c in colnames], colnames)
end

covarsdf, counts, terms = loadsoturet(cache=false)


text:
module Ising

export Ensamble, medicionesIndependientes,magnetizaciónensamble,Energíaensamble_Fperiódica, correSweeps

type Ensamble
    configuración::Array{Int}
    H::Float64
    β::Float64
    M::Float64
    δM::Float64
    E::Float64
    δE::Float64
    J::Float64
    function Ensamble(a::Tuple,β::Float64)
	    configuración = rand(-1:2:1,a)
	    new(configuración,0,β,0,0,0,0,1)
    end
end

medicionesIndependientes(a,τ) = a[2*τ:2*τ:(length(a)-mod(length(a),2*τ))]

function eligeEspín(ensamble)
    configuración = ensamble.configuración
    tamaño = size(configuración)
    m = tamaño[1]
    n = tamaño[2]
    ñ = tamaño[3]
    espín = (rand(1:m), rand(1:n), rand(1:ñ))
end

function Erenglónensamble_Flibre(renglón)
    E=0
    l = length(renglón)
    for i in 1:l-1
        E += - renglón[i] * renglón[i+1]
    end
    E
end

function Energíaensamble_Flibre(ensamble)
    configuración = ensamble.configuración
    E = 0
    tamaño = size(configuración)
    m = tamaño[1]
    n = tamaño[2]
    ñ = tamaño[3]
    for i in 1:m, j in 1:n
        E += ensamble.J*Erenglónensamble_Flibre(configuración[i,j,1:end])
    end
    for i in 1:n, j in 1:ñ
        E += ensamble.J*Erenglónensamble_Flibre(configuración[1:end,i,j])
    end
    for i in 1:m, j in 1:ñ
        E += ensamble.J*Erenglónensamble_Flibre(configuración[i,1:end,j])
    end
    E
end

function magnetizaciónensamble(ensamble::Ensamble)
    configuración = ensamble.configuración
    sum(configuración)
end
function magnetizaciónensamble(configuración::Array)
    sum(configuración)
end

function δEnergía(ensamble,espín)
	configuración = ensamble.configuración
	tamaño = size(configuración)
	m = tamaño[1]
	n = tamaño[2]
	ñ = tamaño[3]
	i=espín[1]
	j=espín[2]
	k=espín[3]
	δE = (ensamble.J*configuración[i,j,mod1(k+1,ñ)]+
	ensamble.J*configuración[i,j,mod1(k-1,ñ)]+
	ensamble.J*configuración[i,mod1(j+1,n),k]+
	ensamble.J*configuración[i,mod1(j-1,n),k]+
	ensamble.J*configuración[mod1(i+1,m),j,k]+
	ensamble.J*configuración[mod1(i-1,m),j,k]+
	ensamble.H*ensamble.J/abs(ensamble.J)
	)*2*configuración[i,j,k]
end

function δMagnetización(ensamble,espín)
	configuración = ensamble.configuración
	δM = -2*configuración[espín...]
end

function cambioestado(ensamble,espín)
       ensamble.δE += Ising.δEnergía(ensamble,espín)
       ensamble.δM += -2*ensamble.configuración[espín...]
       end

function paso(ensamble)
	espín = eligeEspín(ensamble)
	δE = δEnergía(ensamble,espín)
	if exp(-ensamble.β*δE)>rand()
		cambioestado(ensamble,espín)
		flipensamble!(ensamble,espín)
	end
end

function sweep(ensamble)
	for j in 1:length(ensamble.configuración)
		paso(ensamble)
	end
end

function correSweeps(ensamble,sweeps)
	Ms = zeros(sweeps)
	Es = zeros(sweeps)
	ensamble.E = Energíaensamble_Fperiódica(ensamble)
	ensamble.M = magnetizaciónensamble(ensamble)
	ensamble.δE = 0
	ensamble.δM = 0
	for i in 1:sweeps
		sweep(ensamble)
		Ms[i] = ensamble.M + ensamble.δM
		Es[i] = ensamble.E + ensamble.δE
	end
	{"M"=>Ms,"E"=>Es}
end

function flipensamble!(ensamble,espín)
    configuración = ensamble.configuración
    configuración[espín...]*= -1
    return espín
end
function Erenglónensamble_Fperiódica(renglón)
    l = length(renglón)
    E = - renglón[1] * renglón[l]
    for i in 1:l-1
        E += - renglón[i] * renglón[i+1]
    end
    E
end

function Energíaensamble_Fperiódica(ensamble::Ensamble)
	configuración=ensamble.configuración
    E = 0
    tamaño = size(configuración)
    m = tamaño[1]
    n = tamaño[2]
    ñ = tamaño[3]
    for i in 1:m, j in 1:n
        E += ensamble.J*Erenglónensamble_Fperiódica(configuración[i,j,1:end])
    end
    for i in 1:n, j in 1:ñ
        E += ensamble.J*Erenglónensamble_Fperiódica(configuración[1:end,i,j])
    end
    for i in 1:m, j in 1:ñ
        E += ensamble.J*Erenglónensamble_Fperiódica(configuración[i,1:end,j])
    end
    E+magnetizaciónensamble(ensamble)*ensamble.H*ensamble.J/abs(ensamble.J)
end

end


text:
using PackageCompiler, Libdl

PackageCompiler.create_library(".", "bundle";
                                lib_name="diffusion",
                                # project=joinpath(@__DIR__, ".."),
                                precompile_execution_file=[joinpath(@__DIR__, "generate_precompile.jl")],
                                precompile_statements_file=[joinpath(@__DIR__, "additional_precompile.jl")],
                                incremental=false,
                                filter_stdlibs=true,
                                include_lazy_artifacts=false,
                                include_transitive_dependencies=false)

max_stars_repo_name: JuliaLabs/CallingJuliaFromCWithMPI

max_issues_repo_name: JuliaLabs/CallingJuliaFromCWithMPI

max_forks_repo_name: JuliaLabs/CallingJuliaFromCWithMPI


text:
# make things more self-consistent later

const kboltz_cgs = 1.380649e-16 # erg/K
const hplanck_cgs = 6.62607015e-27 # erg*s
const c_cgs = 2.99792458e10 #cm/s

const electron_mass_cgs = 9.1093897e-28 #g
const electron_charge_cgs = 4.80320425e-10 # statcoulomb or cm^3/2 * g^1/2 / s
const amu_cgs = 1.6605402e-24 #g

# 2018 CODATA recommended value:
const bohr_radius_cgs =  5.29177210903e-9 # cm

const kboltz_eV = 8.617333262145e-5 #eV/K
const hplanck_eV = 4.135667696e-15 #eV*s
const RydbergH_eV = 13.595 #eV - this could be updated
const Rydberg_eV = 13.605693122994 #eV 2018 CODATA via wikipedia


text:
module LFkDTree

abstract type Node end

mutable struct INode <: Node
    i::Float64
    c::Float64
    lt::Ref{Node}
    rt::Ref{Node}
    pr::Ref{Node}
end

mutable struct LNode{K} <: Node where {K}
    k::K
end

end # module


text:
# v0.11 deprecations
@deprecate poisson poisson_loss
@deprecate hinge hinge_loss
@deprecate squared_hinge squared_hinge_loss
@deprecate binarycrossentropy(ŷ, y) bce_loss(ŷ, y, agg=identity)
@deprecate logitbinarycrossentropy(ŷ, y) logitbce_loss(ŷ, y, agg=identity)
@deprecate normalise(x) normalise(x, dims=1)


text:
#################### Heidelberger and Welch Diagnostic ####################

function heideldiag(x::Vector{T}; alpha::Real=0.05, eps::Real=0.1,
                             etype=:imse, start::Integer=1, args...) where {T<:Real}
  n = length(x)
  delta = trunc(Int, 0.10 * n)
  y = x[trunc(Int, n / 2):end]
  S0 = length(y) * mcse(y, etype; args...)^2
  i, pvalue, converged, ybar = 1, 1.0, false, NaN
  while i < n / 2
    y = x[i:end]
    m = length(y)
    ybar = mean(y)
    B = cumsum(y) - ybar * collect(1:m)
    Bsq = (B .* B) ./ (m * S0)
    I = sum(Bsq) / m
    pvalue = 1.0 - pcramer(I)
    converged = pvalue > alpha
    if converged
      break
    end
    i += delta
  end
  halfwidth = sqrt(2.0) * erfinv(1.0 - alpha) * mcse(y, etype; args...)
  passed = halfwidth / abs(ybar) <= eps
  [i + start - 2, converged, round(pvalue, digits = 4), ybar, halfwidth, passed]
end

function heideldiag(chn::AbstractChains;
                    alpha = 0.05,
                    eps = 0.1,
                    etype = :imse,
                    sections::Vector{Symbol}=[:parameters],
                    showall=false,
                    args...
                   )
    c = showall ?
       sort(chn) :
       Chains(chn, _clean_sections(chn, sections); sorted=true)

    # Preallocate.
    _, p, m = size(c.value)
    vals = [Array{Float64}(undef, p, 6) for i in 1:m]

    # Perform tests.
    for j in 1:p, k in 1:m
        vals[k][j, :] = heideldiag(
                            collect(skipmissing(c.value[:, j, k])),
                            alpha=alpha,
                            eps=eps,
                            etype=etype,
                            start=first(c);
                            args...
                           )
    end

    colnames = Symbol.(["parameters", "Burn-in", "Stationarity", "p-value", "Mean",
        "Halfwidth", "Test"])

    # Round values.
    pnames = Symbol.(names(c))
    vals = map(x -> round.(x, digits=4), vals)
    columns = [vcat([pnames], [vals[k][:,i] for i in 1:6]) for k in 1:m]

    dfs = [DataFrame(columns[k], colnames) for k in 1:m]
    dfs_wrapped = [ChainDataFrame("Heidelberger and Welch Diagnostic - Chain $k",
                   dfs[k]) for k in 1:m]
    return dfs_wrapped
end


text:
using Celeste: DeterministicVI, SensitiveFloats

using Base.Test
using Distributions
using DerivativeTestUtils
using StaticArrays

import SampleData: gen_two_body_dataset, true_star_init

function test_set_hess()
    sf = SensitiveFloat{Float64}(length(ids), 1, true, true)
    set_hess!(sf, 2, 3, 5.0)
    @test_approx_eq sf.h[2, 3] 5.0
    @test_approx_eq sf.h[3, 2] 5.0

    set_hess!(sf, 4, 4, 6.0)
    @test_approx_eq sf.h[4, 4] 6.0
end

function test_bvn_cov()
    e_axis = .7
    e_angle = pi/5
    e_scale = 2.

    manual_11 = e_scale^2 * (1 + (e_axis^2 - 1) * (sin(e_angle))^2)
    util_11 = DeterministicVI.get_bvn_cov(e_axis, e_angle, e_scale)[1,1]
    @test_approx_eq util_11 manual_11

    manual_12 = e_scale^2 * (1 - e_axis^2) * (cos(e_angle)sin(e_angle))
    util_12 = DeterministicVI.get_bvn_cov(e_axis, e_angle, e_scale)[1,2]
    @test_approx_eq util_12 manual_12

    manual_22 = e_scale^2 * (1 + (e_axis^2 - 1) * (cos(e_angle))^2)
    util_22 = DeterministicVI.get_bvn_cov(e_axis, e_angle, e_scale)[2,2]
    @test_approx_eq util_22 manual_22
end

function test_active_sources()
    # Test that the derivatives of the expected brightnesses partition in
    # active_sources.

    images, ea, bodies = gen_two_body_dataset()

    s = 1
    n = 5
    p = ea.patches[s, n]

    # for subsequent tests, ensure that the second light source
    # has a radius at least as large as both the height and the
    # width of image 3, and that its center is within the image
    @test p.radius_pix >= 23
    @test 0.5 <= p.center[1] <= 20.5
    @test 0.5 <= p.center[2] <= 23.5

    # this patch is huge, the bottom left corner should be the
    # bottom left of the image
    @test p.bitmap_offset == [0, 0]

    # this patch is huge, all the pixels should be active pixels
    @test size(p.active_pixel_bitmap) == size(images[n].pixels)
    @test all(p.active_pixel_bitmap)

    # lets make the second source have no active pixels
    s2 = 2
    for i in 1:ea.N
        fill!(ea.patches[s2, i].active_pixel_bitmap, false)
    end

    P = length(CanonicalParams)
    ea_no2 = ElboArgs(ea.images, ea.vp, ea.patches, [1, 2])
    elbo_no2 = DeterministicVI.elbo_likelihood(ea_no2)
    @test_approx_eq elbo_no2.d[:, 2] zeros(P, 1)

    # lets make the second source have only a few active pixels
    # in one band, so it and source 1 don't overlap completely
    p2 = ea.patches[s2, n]
    p2.active_pixel_bitmap[10:11,10:11] = true

    # now on to a main test--active source order shouldn't matter
    ea12 = ElboArgs(ea.images, ea.vp, ea.patches, [1, 2])
    elbo_lik_12 = DeterministicVI.elbo_likelihood(ea12)

    ea21 = ElboArgs(ea.images, ea.vp, ea.patches, [2, 1])
    elbo_lik_21 = DeterministicVI.elbo_likelihood(ea21)
    @test_approx_eq elbo_lik_12.v[] elbo_lik_21.v[]
    @test_approx_eq elbo_lik_12.d[:,1] elbo_lik_21.d[:,2]
    @test_approx_eq elbo_lik_12.d[:,2] elbo_lik_21.d[:,1]

    # next main test--deriviatives for active sources don't
    ea1 = ElboArgs(ea.images, ea.vp, ea.patches, [1,])
    elbo_lik_1 = DeterministicVI.elbo_likelihood(ea1)
    # source 1 includes all of source 2's active pixels
    @test_approx_eq elbo_lik_1.v[] elbo_lik_12.v[]

    ea2 = ElboArgs(ea.images, ea.vp, ea.patches, [2,])
    elbo_lik_2 = DeterministicVI.elbo_likelihood(ea2)

    @test_approx_eq elbo_lik_12.d[:, 1] elbo_lik_1.d[:, 1]
    @test_approx_eq elbo_lik_12.d[:, 2] elbo_lik_2.d[:, 1]

    @test_approx_eq elbo_lik_12.h[1:P, 1:P] elbo_lik_1.h
    @test_approx_eq elbo_lik_12.h[(1:P) + P, (1:P) + P] elbo_lik_2.h
end

function test_that_variance_is_low()
    # very peaked variational distribution---variance for F(m) should be low
    images, ea, body = true_star_init()
    n = 1

    star_mcs, gal_mcs = Model.load_bvn_mixtures(ea.S, ea.patches,
                                ea.vp, ea.active_sources,
                                ea.psf_K, n)
    sbs = DeterministicVI.SourceBrightness{Float64}[
        DeterministicVI.SourceBrightness(ea.vp[s]) for s in 1:ea.S]

    elbo_vars = DeterministicVI.ElboIntermediateVariables(
      Float64, ea.S, length(ea.active_sources))

    clear!(elbo_vars.E_G)
    clear!(elbo_vars.var_G)

    h, w = 10, 12
    DeterministicVI.get_expected_pixel_brightness!(
      elbo_vars, n, h, w, sbs, star_mcs, gal_mcs, ea)

    @test 0 < elbo_vars.var_G.v[] < 1e-2 * elbo_vars.E_G.v[]^2
end

function test_that_star_truth_is_most_likely()
    images, ea, body = true_star_init()
    best = DeterministicVI.elbo_likelihood(ea)

    for bad_a in [.3, .5, .9]
        ea_a = deepcopy(ea)
        ea_a.vp[1][ids.a[:, 1]] = [ 1.0 - bad_a, bad_a ]
        bad_a_lik = DeterministicVI.elbo_likelihood(ea_a)
        @test best.v[] > bad_a_lik.v[]
    end

    for h2 in -2:2
        for w2 in -2:2
            if !(h2 == 0 && w2 == 0)
                ea_mu = deepcopy(ea)
                ea_mu.vp[1][ids.u] += [h2 * .5, w2 * .5]
                bad_mu = DeterministicVI.elbo_likelihood(ea_mu)
                @test best.v[] > bad_mu.v[]
            end
        end
    end

    for delta in [.7, .9, 1.1, 1.3]
        ea_r1 = deepcopy(ea)
        ea_r1.vp[1][ids.r1] += log(delta)
        bad_r1 = DeterministicVI.elbo_likelihood(ea_r1)
        @test best.v[] > bad_r1.v[]
    end

    for b in 1:4
        for delta in [-.3, .3]
            ea_c1 = deepcopy(ea)
            ea_c1.vp[1][ids.c1[b, 1]] += delta
            bad_c1 = DeterministicVI.elbo_likelihood(ea_c1)
            @test best.v[] > bad_c1.v[]
        end
    end

end

function test_that_galaxy_truth_is_most_likely()
    images, ea, body = gen_sample_galaxy_dataset(perturb=false)
    ea.vp[1][ids.a[:, 1]] = [ 0.01, .99 ]
    best = DeterministicVI.elbo_likelihood(ea)

    for bad_a in [.3, .5, .9]
        ea_a = deepcopy(ea)
        ea_a.vp[1][ids.a[:, 1]] = [ 1.0 - bad_a, bad_a ]
        bad_a = DeterministicVI.elbo_likelihood(ea_a)
        @test best.v[] > bad_a.v[];
    end

    for h2 in -2:2
        for w2 in -2:2
            if !(h2 == 0 && w2 == 0)
                ea_mu = deepcopy(ea)
                ea_mu.vp[1][ids.u] += [h2 * .5, w2 * .5]
                bad_mu = DeterministicVI.elbo_likelihood(ea_mu)
                @test best.v[] > bad_mu.v[]
            end
        end
    end

    for bad_scale in [.8, 1.2]
        ea_r1 = deepcopy(ea)
        ea_r1.vp[1][ids.r1] += 2 * log(bad_scale)
        bad_r1 = DeterministicVI.elbo_likelihood(ea_r1)
        @test best.v[] > bad_r1.v[]
    end

    for n in [:e_axis, :e_angle, :e_scale]
        for bad_scale in [.8, 1.2]
            ea_bad = deepcopy(ea)
            ea_bad.vp[1][getfield(ids, n)] *= bad_scale
            bad_elbo = DeterministicVI.elbo_likelihood(ea_bad)
            @test best.v[] > bad_elbo.v[]
        end
    end

    for b in 1:4
        for delta in [-.3, .3]
            ea_c1 = deepcopy(ea)
            ea_c1.vp[1][ids.c1[b, 2]] += delta
            bad_c1 = DeterministicVI.elbo_likelihood(ea_c1)
            @test best.v[] > bad_c1.v[]
        end
    end
end

function test_coadd_cat_init_is_most_likely()  # on a real stamp
    stamp_id = "5.0073-0.0739_2kpsf"
    images = SampleData.load_stamp_blob(datadir, stamp_id)

    cat_entries = SampleData.load_stamp_catalog(datadir, "s82-$stamp_id", images)
    bright(ce) = sum(ce.star_fluxes) > 3 || sum(ce.gal_fluxes) > 3
    cat_entries = filter(bright, cat_entries)

    ce_pix_locs =
      [ [ WCS.world_to_pix(images[b].wcs, ce.pos) for b=1:5 ]
        for ce in cat_entries ]

    function ce_inbounds(ce)
        pix_locs = [ WCS.world_to_pix(images[b].wcs, ce.pos) for b=1:5 ]
        inbounds(pos) = pos[1] > -10. && pos[2] > -10 &&
                        pos[1] < 61 && pos[2] < 61
        reduce(&, [inbounds(pos) for pos in pix_locs])
    end
    cat_entries = filter(ce_inbounds, cat_entries)

    ea = SampleData.make_elbo_args(images, cat_entries)
    for s in 1:length(cat_entries)
        ea.vp[s][ids.a[2, 1]] = cat_entries[s].is_star ? 0.01 : 0.99
        ea.vp[s][ids.a[1, 1]] = 1.0 - ea.vp[s][ids.a[2, 1]]
    end
    best = DeterministicVI.elbo_likelihood(ea)

    # s is the brightest source.
    s = 1

    for bad_scale in [.7, 1.3]
        ea_r1 = deepcopy(ea)
        ea_r1.vp[s][ids.r1] += 2 * log(bad_scale)
        bad_r1 = DeterministicVI.elbo_likelihood(ea_r1)
        @test best.v[] > bad_r1.v[]
    end

    for n in [:e_axis, :e_angle, :e_scale]
        for bad_scale in [.6, 1.8]
            ea_bad = deepcopy(ea)
            ea_bad.vp[s][getfield(ids, n)] *= bad_scale
            bad_elbo = DeterministicVI.elbo_likelihood(ea_bad)
            @test best.v[] > bad_elbo.v[]
        end
    end

    for bad_a in [.3, .7]
        ea_a = deepcopy(ea)
        ea_a.vp[s][ids.a[:, 1]] = [ 1.0 - bad_a, bad_a ]

        bad_a = DeterministicVI.elbo_likelihood(ea_a)
        @test best.v[] > bad_a.v[]
    end

    for h2 in -2:2
        for w2 in -2:2
            if !(h2 == 0 && w2 == 0)
                ea_mu = deepcopy(ea)
                ea_mu.vp[s][ids.u] += [0.5h2, 0.5w2]
                bad_mu = DeterministicVI.elbo_likelihood(ea_mu)
                @test best.v[] > bad_mu.v[]
            end
        end
    end

    for b in 1:4
        for delta in [-2., 2.]
            ea_c1 = deepcopy(ea)
            ea_c1.vp[s][ids.c1[b, :]] += delta
            bad_c1 = DeterministicVI.elbo_likelihood(ea_c1)
            info("$(best.v[])  >  $(bad_c1.v[])")
            @test best.v[] > bad_c1.v[]
        end
    end
end

function test_num_allowed_sd()
    images, ea, body = gen_two_body_dataset()

    ea.num_allowed_sd = Inf
    elbo_inf = DeterministicVI.elbo(ea)

    ea.num_allowed_sd = 3
    elbo_4sd = DeterministicVI.elbo(ea)

    @test_approx_eq elbo_inf.v[] elbo_4sd.v[]
    @test_approx_eq elbo_inf.d elbo_4sd.d
    @test_approx_eq elbo_inf.h elbo_4sd.h
end

function test_populate_fsm!()
    images, ea, body = gen_two_body_dataset()

    n = 3
    h = w = 5
    s = 2
    ea.active_sources = [s]

    star_mcs, gal_mcs = Model.load_bvn_mixtures(ea.S, ea.patches,
                                ea.vp, ea.active_sources,
                                ea.psf_K, n)
    Model.populate_fsm_vecs!(ea.elbo_vars.bvn_derivs,
                             ea.elbo_vars.fs0m_vec,
                             ea.elbo_vars.fs1m_vec,
                             ea.patches,
                             ea.active_sources,
                             ea.num_allowed_sd,
                             n, h, w,
                             gal_mcs, star_mcs)

    fs0m = SensitiveFloat{Float64}(length(StarPosParams), 1, true, true)
    fs1m = SensitiveFloat{Float64}(length(GalaxyPosParams), 1, true, true)

    x = @SVector Float64[h, w]
    elbo_vars = ea.elbo_vars
    Model.populate_fsm!(elbo_vars.bvn_derivs,
                        fs0m, fs1m,
                        s, x, true,
                        ea.num_allowed_sd,
                        ea.patches[s, n].wcs_jacobian,
                        gal_mcs, star_mcs)

    @test_approx_eq fs0m.v[] elbo_vars.fs0m_vec[s].v[]
    @test_approx_eq fs0m.d elbo_vars.fs0m_vec[s].d
    @test_approx_eq fs0m.h elbo_vars.fs0m_vec[s].h

    @test_approx_eq fs1m.v[] elbo_vars.fs1m_vec[s].v[]
    @test_approx_eq fs1m.d elbo_vars.fs1m_vec[s].d
    @test_approx_eq fs1m.h elbo_vars.fs1m_vec[s].h
end

test_active_sources()
test_set_hess()
test_bvn_cov()
test_num_allowed_sd()
#test_that_variance_is_low()
test_that_star_truth_is_most_likely()
test_that_galaxy_truth_is_most_likely()
test_coadd_cat_init_is_most_likely()
test_populate_fsm!()


text:
using Test
using MCPhyloTree
@testset "MCPhyloTree" begin
    include("basics.jl")
    println("done basics")

    include("search.jl")
    println("done search")

    include("traversal.jl")
    println("done traversal")

    include("consensus.jl")
    println("done consensus")

    include("ladderize.jl")
    println("done ladderize")

    include("pruning.jl")
    println("done pruning")

    include("spr.jl")
    println("done spr")

    include("moves.jl")
    println("done moves")

    include("building.jl")
    println("done building")

    include("distances.jl")
    println("done distances")
end


text:
@enum Shape CIRCLE RECTANGLE ROUNDED_RECTANGLE DISTANCEFIELD TRIANGLE
@enum RaymarchAlgorithm IsoValue Absorption MaximumIntensityProjection
@enum CubeSides TOP BOTTOM FRONT BACK RIGHT LEFT

immutable Grid{N, T <: Range}
    dims::NTuple{N, T}
end
Base.ndims{N,T}(::Grid{N,T}) = N

Grid(ranges::Range...) = Grid(ranges)
function Grid{N, T}(a::Array{T, N})
	s = Vec{N, Float32}(size(a))
	smax = maximum(s)
	s = s./smax
	Grid(ntuple(Val{N}) do i
		linspace(0, s[i], size(a, i))
	end)
end

Grid(a::AbstractArray, ranges...) = Grid(a, ranges)

"""
This constructor constructs a grid from ranges given as a tuple.
Due to the approach, the tuple `ranges` can consist of NTuple(2, T)
and all kind of range types. The constructor will make sure that all ranges match
the size of the dimension of the array `a`.
"""
function Grid{T, N}(a::AbstractArray{T, N}, ranges::Tuple)
    length(ranges) =! N && throw(ArgumentError(
        "You need to supply a range for every dimension of the array. Given: $ranges
        given Array: $(typeof(a))"
    ))
	Grid(ntuple(Val{N}) do i
		linspace(first(ranges[i]), last(ranges[i]), size(a, i))
	end)
end

Base.length(p::Grid) = prod(size(p))
Base.size(p::Grid) = map(length, p.dims)
function Base.getindex{N,T}(p::Grid{N,T}, i)
    inds = ind2sub(size(p), i)
    Point{N, eltype(T)}(ntuple(Val{N}) do i
        p.dims[i][inds[i]]
    end)
end

Base.start(g::Grid) = 1
Base.done(g::Grid, i) = i > length(g)
Base.next(g::Grid, i) = g[i], i+1

GLAbstraction.isa_gl_struct(x::Grid) = true
GLAbstraction.toglsltype_string{N,T}(t::Grid{N,T}) = "uniform Grid$(N)D"
function GLAbstraction.gl_convert_struct{N,T}(g::Grid{N,T}, uniform_name::Symbol)
    return Dict{Symbol, Any}(
        Symbol("$uniform_name.minimum") => Vec{N,Float32}(map(first, g.dims)),
        Symbol("$uniform_name.maximum") => Vec{N,Float32}(map(last, g.dims)),
        Symbol("$uniform_name.dims")    => Vec{N,Cint}(map(length, g.dims)),
        Symbol("$uniform_name.multiplicator") => Vec{N,Float32}(map(x->1/x.divisor, g.dims)),
    )
end
function GLAbstraction.gl_convert_struct{T}(g::Grid{1,T}, uniform_name::Symbol)
    return Dict{Symbol, Any}(
        Symbol("$uniform_name.minimum") => Float32(first(g.dims[1])),
        Symbol("$uniform_name.maximum") => Float32(last(g.dims[1])),
        Symbol("$uniform_name.dims")    => Cint(length(g.dims[1])),
        Symbol("$uniform_name.multiplicator") => Float32(1/g.dims[1].divisor),

    )
end
import Base: getindex, length, next, start, done

to_cpu_mem(x) = x
to_cpu_mem(x::GPUArray) = gpu_data(x)

typealias ScaleTypes Union{Vector, Vec, AbstractFloat, Void, Grid}
typealias PositionTypes Union{Vector, Point, AbstractFloat, Void, Grid}

type ScalarRepeat{T}
    scalar::T
end
Base.ndims(::ScalarRepeat) = 1
Base.getindex(s::ScalarRepeat, i...) = s.scalar
#should setindex! really be allowed? It will set the index for the whole row...
Base.setindex!{T}(s::ScalarRepeat{T}, value, i...) = (s.scalar = T(value))
Base.eltype{T}(::ScalarRepeat{T}) = T

Base.start(::ScalarRepeat) = 1
Base.next(sr::ScalarRepeat, i) = sr.scalar, i+1
Base.done(sr::ScalarRepeat, i) = false

immutable Instances{P,T,S,R}
    primitive::P
    translation::T
    scale::S
    rotation::R
end

function _Instances(position,px,py,pz, scale,sx,sy,sz, rotation, primitive)
    args = (position,px,py,pz, scale,sx,sy,sz, rotation, primitive)
    args = map(to_cpu_mem, args)
    p = const_lift(ArrayOrStructOfArray, Point3f0, args[1:4]...)
    s = const_lift(ArrayOrStructOfArray, Vec3f0, args[5:8]...)
    r = const_lift(ArrayOrStructOfArray, Vec3f0, args[9])
    const_lift(Instances, args[10], p, s, r)
end
function _Instances(position, scale, rotation, primitive)
    p = const_lift(ArrayOrStructOfArray, Point3f0, position)
    s = const_lift(ArrayOrStructOfArray, Vec3f0, scale)
    r = const_lift(ArrayOrStructOfArray, Vec3f0, rotation)
    const_lift(Instances, primitive, p, s, r)
end

immutable GridZRepeat{G,T,N} <: AbstractArray{Point{3,T}, N}
    grid::G
    z::Array{T, N}
end
Base.size(g::GridZRepeat) = size(g.z)
Base.size(g::GridZRepeat, i) = size(g.z, i)
Base.linearindexing{T<:GridZRepeat}(::Type{T}) = Base.LinearFast()
Base.getindex{G,T}(g::GridZRepeat{G,T}, i) = Point{3, T}(g.grid[i], g.z[i])

function ArrayOrStructOfArray{T}(::Type{T}, array::Void, a, elements...)
    StructOfArrays(T, a, elements...)
end
function ArrayOrStructOfArray{T}(::Type{T}, array::FixedVector, a, elements...)
    StructOfArrays(T, a, elements...)
end
function ArrayOrStructOfArray{T}(::Type{T}, scalar::FixedVector, a::Void, elements::Void...)
    ScalarRepeat(transformation_convert(T, scalar))
end
function ArrayOrStructOfArray{T1,T2}(::Type{T1}, array::Array{T2}, a::Void, elements::Void...)
    array
end
function ArrayOrStructOfArray{T1<:Point}(::Type{T1}, grid::Grid, x::Void, y::Void, z::Array)
    GridZRepeat(grid, z)
end
function ArrayOrStructOfArray{T1<:Point}(::Type{T1}, array::Grid, a::Void, elements::Void...)
    array
end
function ArrayOrStructOfArray{T}(::Type{T}, scalar::T)
    ScalarRepeat(scalar)
end
function ArrayOrStructOfArray{T}(::Type{T}, array::Array)
    array
end

transformation_convert{T}(::Type{T}, scalar) = convert(T, scalar)
function transformation_convert{T1<:FixedVector,T2<:FixedVector}(
        ::Type{T1}, scalar::T2
    )
    T1(scalar)
end

function transformation_convert{T1,T2,N1,N2}(
        PT::Type{Point{N1, T1}}, scalar::FixedVector{N2,T2}
    )
    PT(scalar, ntuple(FixedSizeArrays.ConstFunctor(T1(0)), Val{N1-N2})...)
end
function transformation_convert{T1,T2,N1,N2}(
        VT::Type{Vec{N1, T1}}, scalar::FixedVector{N2, T2}
    )
    VT(scalar, ntuple(FixedSizeArrays.ConstFunctor(T1(1)), Val{N1-N2})...)
end

immutable TransformationIterator{T,S,R}
    translation::T
    scale::S
    rotation::R
end
function TransformationIterator(instances::Instances)
    TransformationIterator(
        instances.translation,
        instances.scale,
        instances.rotation
    )
end
function start(t::TransformationIterator)
    start(t.translation),start(t.scale),start(t.rotation)
end
function done(t::TransformationIterator, state)
    done(t.translation, state[1]) ||
    done(t.scale, state[2]) ||
    done(t.rotation, state[3])
end
function next(t::TransformationIterator, state)
    _translation, st = next(t.translation, state[1])
    _scale, ss = next(t.scale, state[2])
    _rotation, sr = next(t.rotation, state[3])

    translation = transformation_convert(Point3f0, _translation)
    scale = transformation_convert(Vec3f0, _scale)
    rotation = transformation_convert(Vec3f0, _rotation)
    v,u = FixedSizeArrays.normalize(rotation), Vec3f0(0,0,1)
    # Unfortunately, we have to check for when u == -v, as u + v
    # in this case will be (0, 0, 0), which cannot be normalized.
    T = Float32
    if (u == -v)
        # 180 degree rotation around any orthogonal vector
        other = (abs(dot(u, Vec{3, T}(1,0,0))) < 1.0) ? Vec{3, T}(1,0,0) : Vec{3, T}(0,1,0)
        q = Quaternions.qrotation(FixedSizeArrays.normalize(cross(u, other)), T(180))
    else
        half = FixedSizeArrays.normalize(u+v)
        q = Quaternions.Quaternion(dot(u, half), cross(u, half)...)
    end
    (translation, scale, Mat{4,4,T}(q)), (st, ss, sr)
end

immutable Intensity{N, T} <: FixedVector{N, T}
	_::NTuple{N, T}
end
typealias GLIntensity Intensity{1, Float32}
export Intensity,GLIntensity

NOT(x) = !x
immutable GLVisualizeShader <: AbstractLazyShader
    paths  ::Tuple
    kw_args::Vector
    function GLVisualizeShader(paths...; kw_args...)
        view = filter(kv->kv[1]==:view, kw_args)
        if isempty(view)
            view = Dict{Compat.UTF8String, Compat.UTF8String}()
        else
            view = view[1][2]
        end

        # TODO properly check what extensions are available
        if !is_apple()
            merge!(view, Dict{Compat.String, Compat.String}(
                "GLSL_EXTENSIONS" => "#extension GL_ARB_conservative_depth: enable",
                "SUPPORTED_EXTENSIONS" => "#define DETPH_LAYOUT"
            ))
        end

        paths = map(paths) do shader
            isa(shader, String) ? loadasset("shader", shader) : shader
        end
        new(paths, vcat(kw_args, [
        	(:fragdatalocation, [(0, "fragment_color"), (1, "fragment_groupid")]),
    		(:updatewhile, current_screen().inputs[:window_open]),
    		(:update_interval, 1.0),
            (:view, view)
        ]))
    end
end


text:
# calculating error between two poses

using Plots
pyplot()

function calc_lon_lat_error(pose::Array, pose_ref::Array)
  error_x = pose[1] - pose_ref[1]
  error_y = pose[2] - pose_ref[2]
  error_theta = pose[3] - pose_ref[3]
  error_lon = cos(pose_ref[3])*error_x + sin(pose_ref[3])*error_y
  error_lat = -sin(pose_ref[3])*error_x + cos(pose_ref[3])*error_y
  return [round(error_lon, digits=2), round(error_lat, digits=2),
          round(error_theta, digits=2)]
end

function draw_lon_lat_error!(pose::Array, pose_ref::Array)
  error = calc_lon_lat_error(pose, pose_ref)
  annotate!(pose_ref[1], pose_ref[2]+0.1,
            text("($(error[1]), $(error[2]), $(error[3]))",
            :left, 8))
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using IntervalArithmetic
using DifferentialEquations
using Plots

function SHO(du,u,p,t)
du[1]=u[2]+@interval(pi)
du[2]=-u[1]
end

T_end=2*pi
u0=[1.,0.]
tspan=(0.,T_end)
prob = ODEProblem(SHO,u0,tspan)
t,u=solve(prob, RK4(),reltol=1e-8,abstol=1e-10,maxiters=1e15)
# plot(u[1,:],u[2,:])
diff=u[1,end]-u[1,1]
println(diff)


text:
using GLPK

@testset "Variables" begin
    @testset "@variable" begin
        sp = Model()
        data = Dict("hydro" => Dict("nHyd" => 4,
        "Hydrogenerators" => [Dict("max_turn" => 200) for r =1:4]))

        # inflow
        HydroPowerModels.variable_inflow(sp, data)
        @test size(sp[:inflow],1) == 4
        # spillage
        HydroPowerModels.variable_spillage(sp, data)
        @test size(sp[:spill],1) == 4
        # outflow
        HydroPowerModels.variable_outflow(sp, data)
        @test size(sp[:outflow],1) == 4
    end

    @testset "@state" begin
        sp = Model()
        data = Dict("hydro" => Dict("nHyd" => 1,
        "Hydrogenerators" => [Dict("min_volume" => 0, "max_volume" => 200, "initial_volume" => 100)]))

        m = SDDP.LinearPolicyGraph(
            sense   = :Min,
            stages  = 1,
            optimizer  = with_optimizer(GLPK.Optimizer),
            lower_bound = 0.0,
                    direct_mode=false
                                            ) do sp,t
            HydroPowerModels.variable_volume(sp, data)
            @stageobjective(sp,0)
        end
        SDDP.train(m;iteration_limit = 60);
        # state variable volume
        @test JuMP.lower_bound(m[1].subproblem[:reservoir][1].out) == 0.0
        @test JuMP.upper_bound(m[1].subproblem[:reservoir][1].out) == 200

        # initial value of state variable volume
        JuMP.optimize!(m[1].subproblem)
        @test JuMP.value(m[1].subproblem[:reservoir][1].in) == 100
    end

end


text:
@testset "#853" begin

"Model parameter"
struct yObj
  p
end

"Univariable Functor"
function (s::yObj)(x)
  return x*s.p
end

"Multivariable Functor with array input"
function (s::yObj)(x_v)
  return x_v[1]*s.p
end

model = yObj(1.0)

Optim.optimize(model, -1, 2, [1.]) # already worked
Optim.optimize(model, -1, 2) # didn't work
Optim.optimize(model, -1., 2.) # didn't work
end