text:
using WignerSymbols: wigner3j
using .PhysicalConstants: e, ca40_qubit_transition_frequency, m_ca40, ħ, α, μB

export mass, level_structure, selected_level_structure, stark_shift,
       selected_matrix_elements, matrix_element, get_basis, ion_number, ion_position,
       gJ, zeeman_shift, matrix_elements, zero_stark_shift, Ion, Ca40

#############################################################################################
# Ion - the physical parameters defining an ion's structure
#############################################################################################

"""
    Ion
The physical parameters defining an isolated ion's internal structure.
"""
abstract type Ion <: IonSimBasis end

# required fields
mass(I::Ion)::Real = I.mass
level_structure(I::Ion)::OrderedDict{String,NamedTuple} = I.level_structure
selected_level_structure(I::Ion)::OrderedDict{String,NamedTuple} = I.selected_level_structure
matrix_elements(I::Ion)::OrderedDict = I.matrix_elements
selected_matrix_elements(I::Ion)::OrderedDict = I.selected_matrix_elements
ion_number(I::Ion)::Union{Int,Missing} = I.number
ion_position(I::Ion)::Union{Real,Missing} = I.position
stark_shift(I::Ion)::OrderedDict{String,Real} = I.stark_shift

#############################################################################################
# Ca40 Ion
#############################################################################################

"""
    Ca40(selected_level_structure::Vector{String}[, stark_shift])

**user-defined fields**
* `selected_level_structure`:
    keys ⊂ `["S-1/2", "S+1/2", "D-5/2", "D-3/2", "D-1/2", "D+1/2", "D+3/2", "D+5/2"]`.
    Values are a `NamedTuple` with:
    * `l`: orbital angular momentum
    * `j`: total angular momentum
    * `mⱼ`: projection of total angular momentum along quantization axis
    * `E`: relative energies
    Note: indexing the instantiated structure with one of these strings will return
    the corresponding `Ket`.
* `stark_shift`: A dictionary with keys denoting the selected levels and values, a real
    number for describing a shift of the level's energy. This is just a convenient way to add
    Stark shifts to the simulation without additional resources.
**fixed fields**
* `mass::Real`: The ion's mass in kg.
* `level_structure`: A full description of the ion's electronic structure.
* `matrix_elements::OrderedDict{Tuple,Function}`: Same as `selected_matrix_elements` but for
    all of the ion's allowed transitions.
**derived fields**
* `selected_matrix_elements`: Functions for the allowed transitions (contained in the
    selected levels) that return the corresponding coupling strengths. These functions take
    as arguments:
    * `Efield`: magnitude of the electric field at the position of the ion [V/m]
    * `γ`: ``ϵ̂⋅B̂`` (angle between laser polarization and B-field)
    * `ϕ`: ``k̂⋅B̂`` (angle between laser k-vector and B-field)
* `shape::Vector{Int}`: Indicates the dimension of the used Hilbert space.
* `number`: When the ion is added to an `IonConfiguration`, this value keeps track of its
    order in the chain.
* `position`: When the ion is added to an `IonConfiguration`, this value keeps track of its
    physical position in meters.
"""
mutable struct Ca40 <: Ion
    mass::Real
    level_structure::OrderedDict{String,NamedTuple}
    selected_level_structure::OrderedDict{String,NamedTuple}
    shape::Vector{Int}
    matrix_elements::OrderedDict{Tuple,Function}
    selected_matrix_elements::OrderedDict{Tuple,Function}
    stark_shift::OrderedDict{String,Real}
    number::Union{Int,Missing}
    position::Union{Real,Missing}
    function Ca40(selected_level_structure; ss=Dict())
        fls, sls_dict, me, me_dict=_structure(selected_level_structure)
        shape = [length(sls_dict)]
        ss_full = OrderedDict{String,Float64}()
        for level in keys(sls_dict)
            haskey(ss, level) ? ss_full[level] = ss[level] : ss_full[level] = 0.
        end
        new(m_ca40, fls, sls_dict, shape, me, me_dict, ss_full, missing, missing)
    end
    Ca40(;ss=Dict()) = Ca40("default", ss=ss)
    # for copying
    function Ca40(
            mass, level_structure, selected_level_structure, shape, matrix_elements,
            selected_matrix_elements, stark_shift, number, position
        )
        selected_level_structure = deepcopy(selected_level_structure)
        shape = copy(shape)
        selected_matrix_elements = deepcopy(selected_matrix_elements)
        stark_shift = deepcopy(stark_shift)
        new(mass, level_structure, selected_level_structure, shape, matrix_elements,
            selected_matrix_elements, stark_shift, number, position)
    end
end

function _structure(selected_level_structure)
    ED = ca40_qubit_transition_frequency
    fls = OrderedDict(
            "S-1/2" => (l=0, j=1//2, mⱼ=-1//2, E=0),
            "S+1/2" => (l=0, j=1//2, mⱼ=1//2, E=0),
            "D-5/2" => (l=2, j=5//2, mⱼ=-5//2, E=ED),
            "D-3/2" => (l=2, j=5//2, mⱼ=-3//2, E=ED),
            "D-1/2" => (l=2, j=5//2, mⱼ=-1//2, E=ED),
            "D+1/2" => (l=2, j=5//2, mⱼ=1//2, E=ED),
            "D+3/2" => (l=2, j=5//2, mⱼ=3//2, E=ED),
            "D+5/2" => (l=2, j=5//2, mⱼ=5//2, E=ED)
        )
    me = OrderedDict{Tuple,Function}()
    k = fls.keys
    for i in eachindex(k), j in i+1:length(k)
        t1, t2 = fls[k[i]], fls[k[j]]
        if ! (typeof(_ca40_matrix_elements((t1, t2), 0, 0, 0)) <: Nothing)
            f(Efield, γ, ϕ) = _ca40_matrix_elements((t1, t2), Efield, γ, ϕ)
            push!(me, (k[i], k[j]) => f)
        end
    end
    if selected_level_structure == "default"
        sls = collect(keys(fls))
    else
        sls = selected_level_structure
    end
    sls_dict = OrderedDict{String,NamedTuple}()
    for k in sls
        @assert k in fls.keys "invalid level $k"
        push!(sls_dict, k => fls[k])
    end
    me_dict = OrderedDict{Tuple,Function}()
    for (k, v) in me
        if k[1] in sls && k[2] in sls
            push!(me_dict, k => v)
        end
    end
    return fls, sls_dict, me, me_dict
end

# geometric part of the matrix element for 40Ca S1/2 <-> D5/2 transitions,
# assuming linearly polarized light
_ca40_geo = [
        (γ, ϕ) -> begin
                γ = deg2rad(γ)
                ϕ = deg2rad(ϕ)
                (1/2)abs(cos(γ)sin(2ϕ))
            end,
        (γ, ϕ) -> begin
                γ = deg2rad(γ)
                ϕ = deg2rad(ϕ)
                sqrt(1/6)abs(cos(γ)cos(2ϕ) + im*sin(γ)cos(ϕ))
            end,
        (γ, ϕ) -> begin
                γ = deg2rad(γ)
                ϕ = deg2rad(ϕ)
                sqrt(1/6)abs((1/2)cos(γ)sin(2ϕ) + im*sin(γ)sin(ϕ))
            end
    ]

function _ca40_matrix_elements(
        transition::Tuple{NamedTuple,NamedTuple}, Efield::Real, γ::Real, ϕ::Real
    )
    t1 = transition[1]
    t2 = transition[2]
    Δl = t2.l - t1.l
    Δm = Int(abs(t2.mⱼ - t1.mⱼ))
    if Δl ≡ 0 || abs(Δm) > 2
        return nothing
    end
    λ = c / ca40_qubit_transition_frequency
    Ω = (e * Efield / ħ) * √(5λ^3 * (1/1.17) / (2 * π^3 * c * α)) / (2π)
    wig = abs(wigner3j(t1.j, t2.j - t1.j, t2.j, -t1.mⱼ, t1.mⱼ - t2.mⱼ, t2.mⱼ))
    Ω * _ca40_geo[Δm+1](γ, ϕ) * wig
end

"""
    matrix_element(transition::Vector{String}, Efield::Real, γ::Real, ϕ::Real)

Computes the coupling strengths of the various S ⟷ D transitions in ⁴⁰Ca.
See e.g. page 30 of
[Roos's thesis](https://quantumoptics.at/images/publications/dissertation/roos_diss.pdf).
Only considers linearly polarized light fields.

**args**
* `C`: Ca40 ion
* `transition`: i.e. ["S-1/2", "D-1/2"]
* `Efield`: magnitude of the electric field at the position of the ion [V/m]
* `γ`: ``ϵ̂⋅B̂`` (angle between laser polarization and B-field)
* `ϕ`: ``k̂⋅B̂`` (angle between laser k-vector and B-field)
"""
function matrix_element(C::Ca40, transition::Vector{String}, Efield::Real, γ::Real, ϕ::Real)
    t1 = C.level_structure[transition[1]]
    t2 = C.level_structure[transition[2]]
    _ca40_matrix_elements((t1, t2), Efield, γ, ϕ)
end

function Base.print(I::Ca40)
    println("⁴⁰Ca\n")
    for (k, v) in I.selected_level_structure
        println(k, ": ", v)
    end
end

Base.show(io::IO, I::Ca40) = println(io, "⁴⁰Ca")  # suppress long output

#############################################################################################
# general functions
#############################################################################################

"""
    gJ(l::real, j::real; s::Real=1/2)
Landé g-factor

**args**
* `l`: orbital angular momentum quantum number
* `j`: total angular momentum quantum number
* `s`: spin angular momentum quantum number (defaults to 1/2)
"""
gJ(l::Real, j::Real; s::Real=1/2) = 3/2 + (s * (s + 1) - l * (l + 1)) / (2j * (j + 1))

"""
    zeeman_shift(B::Real, l::Real, j::Real, mⱼ::Real)
``ΔE = (μ_B/ħ) ⋅ g_J(l, j) ⋅ B ⋅ mⱼ / 2π``
**args**
* `B`: magnitude of B-field at ion
* `l`: orbital angular momentum quantum number
* `j`: total angular momentum quantum number
* `mⱼ`: projection of total angular momentum along quantization axis
"""
zeeman_shift(B::Real, l::Real, j::Real, mⱼ::Real) = zeeman_shift(B, (l=l, j=j, mⱼ=mⱼ))
zeeman_shift(B::Real, p::NamedTuple) = (μB/ħ) * gJ(p.l, p.j) * B * p.mⱼ / 2π
zeeman_shift(B::Real, p::Tuple) = zeeman_shift(B, (l=p[1], j=p[2], mⱼ=p[3]))
zeeman_shift(;B::Real, l::Real, j::Real, mⱼ::Real) = zeeman_shift(B, (l=l, j=j, mⱼ=mⱼ))

Base.getindex(I::Ion, state::String) = ionstate(I, state)
Base.getindex(I::Ion, state::Int) = ionstate(I, state)

function Base.getproperty(I::Ion, s::Symbol)
    if s == :number || s == :position
        if typeof(getfield(I, s)) <: Missing
            @warn "ion has not been added to a configuration"
        return missing
        end
    end
    getfield(I, s)
end

function zero_stark_shift(I::Ion)
    for k in keys(I.stark_shift)
        I.stark_shift[k] = 0.0
    end
end

function Base.setproperty!(I::Ion, s::Symbol, v::Tv) where{Tv}
    if (s == :mass ||
        s == :level_structure ||
        s == :shape ||
        s == :matrix_elements ||
        s == :selected_matrix_elements ||
        s == :number ||
        s == :position)
        return
    elseif s == :selected_level_structure
        @assert Tv == Vector{String} "type must be Vector{String}"
        _, sls_dict, _, me_dict = _structure(v)
        Core.setproperty!(I, :selected_level_structure, sls_dict)
        Core.setproperty!(I, :selected_matrix_elements, me_dict)
        Core.setproperty!(I, :shape, [length(sls_dict)])
        I.stark_shift = OrderedDict{String,Real}()
        for key in v
            I.stark_shift[key] = 0.0
        end
        return
    end
    Core.setproperty!(I, s, v)
end

function Base.:(==)(b1::T, b2::T) where {T<:Ion}
    (
        b1.mass == b2.mass &&
        b1.selected_level_structure == b2.selected_level_structure &&
        b1.shape == b2.shape &&
        b1.stark_shift == b2.stark_shift
    )
end


text:
module FQR

using LinearAlgebra, Printf
using StaticArrays
# using GenericLinearAlgebra

const Float = Float64

### faster, a couple of less digits
using DoubleFloats
# const Float = Double64
const Double = Double64

# similar speed as DoubleFloats
# using MultiFloats
# const Float = Float64x2
# const Double = Float64x2

### a couple of more digits, but slower
# using Quadmath
# const Float = Float128

### 64 digits by default, but a lot more slower
# const Float = BigFloat

###################### traits to the functional QR  ############################
abstract type Grid end
abstract type FineMesh end

dot(mesh, g1, g2) = error("QR.dot is not implemented!")
mirror(g) = error("QR.mirror for $(typeof(g)) is not implemented!")
irreducible(g) = error("QR.irreducible for $(typeof(g)) is not implemented!")
#################################################################################

mutable struct Basis{D,Grid,Mesh}
    ############    fundamental parameters  ##################
    Λ::Float  # UV energy cutoff * inverse temperature
    rtol::Float # error tolerance

    ###############     DLR grids    ###############################
    N::Int # number of basis
    grid::Vector{Grid} # grid for the basis
    error::Vector{Float}  # the relative error achieved by adding the current grid point

    ###############  linear coefficients for orthognalization #######
    Q::Matrix{Double} # , Q = R^{-1}, Q*R'= I
    R::Matrix{Double}

    ############ fine mesh #################
    mesh::Mesh

    function Basis{d,Grid}(Λ, rtol, mesh::Mesh) where {d,Grid,Mesh}
        _Q = Matrix{Float}(undef, (0, 0))
        _R = similar(_Q)
        return new{d,Grid,Mesh}(Λ, rtol, 0, [], [], _Q, _R, mesh)
    end
end

function addBasis!(basis::Basis{D,G,M}, grid, verbose) where {D,G,M}
    basis.N += 1
    push!(basis.grid, grid)

    basis.Q, basis.R = GramSchmidt(basis)

    # println(maximum(basis.mesh.residual))
    # update the residual on the fine mesh
    updateResidual!(basis)

    # println(maximum(basis.mesh.residual))
    # the new rtol achieved by adding the new grid point
    push!(basis.error, sqrt(maximum(basis.mesh.residual)))

    (verbose > 0) && @printf("%3i %s -> error=%16.8g, Rmin=%16.8g\n", basis.N, "$(grid)", basis.error[end], basis.R[end, end])
end

function addBasisBlock!(basis::Basis{D,G,M}, idx, verbose) where {D,G,M}
    _norm = sqrt(basis.mesh.residual[idx]) # the norm derived from the delta update in updateResidual
    addBasis!(basis, basis.mesh.candidates[idx], verbose)
    _R = basis.R[end, end] # the norm derived from the GramSchmidt

    @assert abs(_norm - _R) < basis.rtol * 100 "inconsistent norm on the grid $(basis.grid[end]) $_norm - $_R = $(_norm-_R)"
    if abs(_norm - _R) > basis.rtol * 10
        @warn("inconsistent norm on the grid $(basis.grid[end]) $_norm - $_R = $(_norm-_R)")
    end

    ## set the residual of the selected grid point to be zero
    basis.mesh.selected[idx] = true
    basis.mesh.residual[idx] = 0 # the selected mesh grid has zero residual

    # println(mirror(basis.mesh, idx))
    for grid in mirror(basis.mesh, idx)
        addBasis!(basis, grid, verbose)
    end
end

function updateResidual!(basis::Basis{D}) where {D}
    mesh = basis.mesh

    # q = Float.(basis.Q[end, :])
    q = Double.(basis.Q[:, end])

    Threads.@threads for idx in 1:length(mesh.candidates)
        if mesh.selected[idx] == false
            candidate = mesh.candidates[idx]
            pp = sum(q[j] * dot(mesh, basis.grid[j], candidate) for j in 1:basis.N)
            _residual = mesh.residual[idx] - pp * pp
            # @assert isnan(_residual) == false "$pp and $([q[j] for j in 1:basis.N]) => $([dot(mesh, basis.grid[j], candidate) for j in 1:basis.N])"
            # println("working on $candidate : $_residual")
            if _residual < 0
                if _residual < -basis.rtol
                    @warn("warning: residual smaller than 0 at $candidate got $(mesh.residual[idx]) - $(pp)^2 = $_residual")
                end
                mesh.residual[idx] = 0
            else
                mesh.residual[idx] = _residual
            end
        end
    end
end

"""
Gram-Schmidt process to the last grid point in basis.grid
"""
function GramSchmidt(basis::Basis{D,G,M}) where {D,G,M}
    _Q = zeros(Double, (basis.N, basis.N))
    _Q[1:end-1, 1:end-1] = basis.Q

    _R = zeros(Double, (basis.N, basis.N))
    _R[1:end-1, 1:end-1] = basis.R
    _Q[end, end] = 1

    newgrid = basis.grid[end]

    overlap = [dot(basis.mesh, basis.grid[j], newgrid) for j in 1:basis.N-1]

    for qi in 1:basis.N-1
        _R[qi, end] = basis.Q[:, qi]' * overlap
        _Q[:, end] -= _R[qi, end] * _Q[:, qi]  # <q, qnew> q
    end

    _norm = dot(basis.mesh, newgrid, newgrid) - _R[:, end]' * _R[:, end]
    _norm = sqrt(abs(_norm))

    @assert _norm > eps(Double(1)) * 100 "$_norm is too small as a denominator!\nnewgrid = $newgrid\nexisting grid = $(basis.grid)\noverlap=$overlap\nR=$_R\nQ=$_Q"

    _R[end, end] = _norm
    _Q[:, end] /= _norm

    return _Q, _R
end

function test(basis::Basis{D}) where {D}
    println("testing orthognalization...")
    KK = zeros(Double, (basis.N, basis.N))
    Threads.@threads for i in 1:basis.N
        g1 = basis.grid[i]
        for (j, g2) in enumerate(basis.grid)
            KK[i, j] = dot(basis.mesh, g1, g2)
        end
    end
    maxerr = maximum(abs.(KK - basis.R' * basis.R))
    println("Max overlap matrix R'*R Error: ", maxerr)

    maxerr = maximum(abs.(basis.R * basis.Q - I))
    println("Max R*R^{-1} Error: ", maxerr)

    II = basis.Q' * KK * basis.Q
    maxerr = maximum(abs.(II - I))
    println("Max Orthognalization Error: ", maxerr)

    # KK = zeros(Float, (basis.N, basis.N))
    # Threads.@threads for i in 1:basis.N
    #     g1 = basis.grid[i]
    #     for (j, g2) in enumerate(basis.grid)
    #         KK[i, j] = dot(basis.mesh, g1, g2)
    #     end
    # end
    # println(maximum(abs.(KK' - KK)))
    # A = cholesky(KK, Val{true}())
    # println(maximum(abs.(A.L * A.U - KK)))
    # println(maximum(abs.(A.L' - A.U)))
end

# function testResidual(basis, proj)
#     # residual = [Residual(basis, proj, basis.grid[i, :]) for i in 1:basis.N]
#     # println("Max deviation from zero residual: ", maximum(abs.(residual)))
#     println("Max deviation from zero residual on the DLR grids: ", maximum(abs.(basis.residualFineGrid[basis.gridIdx])))
# end

function qr!(basis::Basis{dim,G,M}; initial = [], N = 10000, verbose = 0) where {dim,G,M}
    #### add the grid in the idx vector first

    for i in initial
        addBasisBlock!(basis, i, verbose)
    end

    ####### add grids that has the maximum residual
    maxResidual, idx = findmax(basis.mesh.residual)
    while sqrt(maxResidual) > basis.rtol && basis.N < N
        addBasisBlock!(basis, idx, verbose)
        maxResidual, idx = findmax(basis.mesh.residual)
    end
    @printf("rtol = %.16e\n", sqrt(maxResidual))
    return basis
end

end

ext: jl


text:
include("new_SimulationFunctions.jl")

data_dictionary=Dict()

TSTART = 0.0;
Ts = 1.0;
TSTOP = 100.0;

(TSIM,X,ICA) =  SolveModelAddTGFB(TSTART,TSTOP,Ts,data_dictionary)

writedlm("/Users/dbassen/work/varner_gould_emt_model/emt_grn_bassen_model/Results/TGFbSimulation/TSIM.dat", TSIM)
writedlm("/Users/dbassen/work/varner_gould_emt_model/emt_grn_bassen_model/Results/TGFbSimulation/X.dat", X)


text:
using IntervalRootFinding
using Test

#include("bisect.jl")
#include("findroots.jl")
include("roots.jl")
include("test_smiley.jl")
include("newton1d.jl")
include("quadratic.jl")
include("linear_eq.jl")
include("slopes.jl")


text:
using Printf: @printf

function debug(verbosity, string...)
    if verbosity > 0
        println(string...)
    end
end

function debug_inline(verbosity, string...)
    if verbosity > 0
        print(string...)
    end
end

function getTime()::Int
    return round(Int, 1e3*(time()-1.6e9))
end

function check_numeric(n)
    return tryparse(Float64, n) !== nothing
end

function is_anonymous_function(op)
	op_string = string(nameof(op))
	return length(op_string) > 1 && op_string[1] == '#' && check_numeric(op_string[2:2])
end


text:
struct CSV_Foo
end

@testset "CSV.File basics" begin

#test on non-existent file
@test_throws ArgumentError CSV.File("");

#test where datarow > headerrow
@test_throws ArgumentError CSV.File(joinpath(dir, "test_no_header.csv"); datarow=1, header=2);

#test bad types
@test_throws CSV.Error CSV.File(joinpath(dir, "test_float_in_int_column.csv"); types=[Int, Int, Int], strict=true)

# Integer overflow; #100
@test_throws CSV.Error CSV.File(joinpath(dir, "int64_overflow.csv"); types=[Int8], strict=true)

# #172
@test_throws ArgumentError CSV.File(joinpath(dir, "test_newline_line_endings.csv"), types=Dict(1=>Integer))

# #289
tmp = CSV.File(IOBuffer(" \"a, b\", \"c\" "), datarow=1) |> DataFrame
@test size(tmp) == (1, 2)
@test tmp.Column1[1] == "a, b"
@test tmp.Column2[1] == "c"

tmp = CSV.File(IOBuffer(" \"2018-01-01\", \"1\" ,1,2,3"), datarow=1) |> DataFrame
@test size(tmp) == (1, 5)
@test tmp.Column1[1] == Date(2018, 1, 1)
@test tmp.Column2[1] == 1
@test tmp.Column3[1] == 1
@test tmp.Column4[1] == 2
@test tmp.Column5[1] == 3

# #329
df = CSV.read(joinpath(dir, "test_types.csv"), types=Dict(:string=>Union{Missing,DateTime}), silencewarnings=true)
@test df.string[1] === missing

# #352
@test_throws KeyError first(CSV.File(joinpath(dir, "test_types.csv"))).a

@test_throws ArgumentError CSV.File(IOBuffer("a\0b\n1\02\n"); delim='\0')
@test_throws ArgumentError CSV.File(IOBuffer("a\0b\n1\02\n"); delim="\0")

f = CSV.File(IOBuffer("a,b\n1,2\n"); header=3)
@test isempty(f.types)
@test f.rows == 0

f = CSV.File(IOBuffer("a,b\n1,2\n"); datarow=3)
@test f.names == [:a, :b]
@test f.rows == 0

f = CSV.File(IOBuffer("a,b\n1,2\n"); skipto=3)
@test f.names == [:a, :b]
@test f.rows == 0

f = CSV.File(IOBuffer("a,b\n1,2\n"); limit=0)
@test f.names == [:a, :b]
@test f.rows == 0

# delim same as quotechar
f = CSV.File(IOBuffer("a\"b\n1\"2\n"); delim='"')
@test f.rows == 0
@test f.cols == 2
@test f.types == [Missing, Missing]

# delim same as quotechar w/ quoted field
f = CSV.File(IOBuffer("a\"b\n1\"\"2\"\n"); delim='"')
@test f.rows == 0
@test f.cols == 2
@test f.types == [Missing, Missing]

# 387
f = CSV.File(IOBuffer("a,b\n1,1\n1,2\n1,3\n"))

for rows in Iterators.partition(f, 2)
    for row in rows
        # get last property
        getproperty(row, :b)
        @test row.a == 1
    end
end

# 391
rows = collect(f)
@test rows[1].a == 1
@test rows[1].a == 1

# 459
rows = collect(CSV.File(joinpath(dir, "time.csv"); dateformat="H:M:S"))
@test rows[1].time == Time(0)
@test rows[2].time == Time(0, 10)

# 388
df = CSV.read(joinpath(dir, "GSM2230757_human1_umifm_counts.csv"))
@test size(df) == (3, 20128)

# 386
f = CSV.File(IOBuffer("fullVisitorId,PredictedLogRevenue\n18966949534117,0\n39738481224681,0\n"), limit=3)
@test f.rows == 2
rows = collect(f)
@test rows[1].fullVisitorId == 18966949534117
@test rows[2].fullVisitorId == 39738481224681

df = CSV.read(IOBuffer("col1,col2\n1.0,hi"), limit=3)
@test size(df) == (1, 2)

# 288
df = CSV.read(IOBuffer("x\n\",\"\n\",\""))
@test size(df) == (2, 1)

# delimiter auto-detection
df = CSV.read(IOBuffer("x\n1\n2\n"))
@test size(df) == (2, 1)

df = CSV.read(IOBuffer("x,y\n1,\n2,\n"))
@test size(df) == (2, 2)

df = CSV.read(IOBuffer("x y\n1 \n2 \n"))
@test size(df) == (2, 2)

df = CSV.read(IOBuffer("x\ty\n1\t\n2\t\n"))
@test size(df) == (2, 2)

df = CSV.read(IOBuffer("x|y\n1|\n2|\n"))
@test size(df) == (2, 2)

# type promotion
# int => float
df = CSV.read(IOBuffer("x\n1\n3.14"))
@test size(df) == (2, 1)
@test df.x[1] === 1.0
@test df.x[2] === 3.14

# int => missing
df = CSV.read(IOBuffer("x\n1\n\n"))
@test size(df) == (2, 1)
@test df.x[1] == 1
@test df.x[2] === missing

# missing => int
df = CSV.read(IOBuffer("x\n\n1\n"))
@test size(df) == (2, 1)
@test df.x[1] === missing
@test df.x[2] == 1

# missing => int => float
df = CSV.read(IOBuffer("x\n\n1\n3.14\n"))
@test size(df) == (3, 1)
@test df.x[1] === missing
@test df.x[2] === 1.0
@test df.x[3] === 3.14

# int => missing => float
df = CSV.read(IOBuffer("x\n1\n\n3.14\n"))
@test size(df) == (3, 1)
@test df.x[1] === 1.0
@test df.x[2] === missing
@test df.x[3] === 3.14

# int => float => missing
df = CSV.read(IOBuffer("x\n1\n3.14\n\n"))
@test size(df) == (3, 1)
@test df.x[1] === 1.0
@test df.x[2] === 3.14
@test df.x[3] === missing

# int => string
df = CSV.read(IOBuffer("x\n1\nabc"))
@test size(df) == (2, 1)
@test df.x[1] == "1"
@test df.x[2] == "abc"

# float => string
df = CSV.read(IOBuffer("x\n3.14\nabc"))
@test size(df) == (2, 1)
@test df.x[1] == "3.14"
@test df.x[2] == "abc"

# int => float => string
df = CSV.read(IOBuffer("x\n1\n3.14\nabc"))
@test size(df) == (3, 1)
@test df.x[1] == "1"
@test df.x[2] == "3.14"
@test df.x[3] == "abc"

# missing => int => float => string
df = CSV.read(IOBuffer("x\n\n1\n3.14\nabc"))
@test size(df) == (4, 1)
@test df.x[1] === missing
@test df.x[2] == "1"
@test df.x[3] == "3.14"
@test df.x[4] == "abc"

# missing => catg
df = CSV.read(IOBuffer("x\n\na\n"), pool=true)
@test size(df) == (2, 1)
@test df.x[1] === missing
@test df.x[2] == "a"

# catg => missing
df = CSV.read(IOBuffer("x\na\n\n"), pool=true)
@test size(df) == (2, 1)
@test df.x[1] == "a"
@test df.x[2] === missing

# catg => string
df = CSV.read(IOBuffer("x\na\nb\na\nb\na\nb\na\nc\nd\ne\nf\ng\nh\ni\nj\nk\nl\nn\nm\no\np\nq\nr\n"), pool=0.5)
@test typeof(df.x) == CSV.Column{String, String}

# a few corner cases for escape strings
df = CSV.read(IOBuffer("\"column name with \"\" escape character inside\"\n1\n"))
@test names(df)[1] == Symbol("column name with \" escape character inside")

df = CSV.read(IOBuffer("\"column name with \"\" escape character inside\",1\n,2"), transpose=true)
@test names(df)[1] == Symbol("column name with \" escape character inside")
@test names(df)[2] == :Column2

df = CSV.read(IOBuffer("x\na\nb\n\"quoted field with \"\" escape character inside\"\n"), pool=true)
@test df.x[1] == "a"
@test df.x[3] == "quoted field with \" escape character inside"

df = CSV.read(IOBuffer("x\na\nb\n\"quoted field with \"\" escape character inside\"\n"), pool=true)
@test df.x[1] == "a"
@test df.x[3] == "quoted field with \" escape character inside"

# invalid quoted field is fatal error
@test_throws CSV.Error CSV.read(IOBuffer("x\n\"quoted field that never ends"))
@test_throws CSV.Error CSV.read(IOBuffer("x\nhey\n\"quoted field that never ends"))
@test_throws CSV.Error CSV.read(IOBuffer("x\n\n\"quoted field that never ends"))
@test_throws CSV.Error CSV.read(IOBuffer("x\n1\n\"quoted field that never ends"))
@test_throws CSV.Error CSV.read(IOBuffer("x\n1.0\n\"quoted field that never ends"))
@test_throws CSV.Error CSV.read(IOBuffer("x\na\n\"quoted field that never ends"), pool=true)

# invalid integer
df = CSV.read(IOBuffer("x\nabc\n"), type=Int)
@test size(df) == (1, 1)
@test df.x[1] === missing

@test_throws CSV.Error CSV.read(IOBuffer("x\nabc\n"), type=Int, strict=true)

# transpose corner cases
df = CSV.read(IOBuffer("x,y,1\nx2,y2,2\n"), transpose=true, header=2)
@test names(df) == [:y, :y2]
@test size(df) == (1, 2)
@test df.y[1] == 1
@test df.y2[1] == 2

df = CSV.read(IOBuffer("x,y,1\nx2,y2,2\n"), transpose=true, header=1, datarow=3)
@test names(df) == [:x, :x2]
@test size(df) == (1, 2)
@test df.x[1] == 1
@test df.x2[1] == 2

df = CSV.read(IOBuffer("x,y,1\nx2,y2,2\n"), transpose=true, header=false, datarow=3)
@test names(df) == [:Column1, :Column2]
@test size(df) == (1, 2)
@test df.Column1[1] == 1
@test df.Column2[1] == 2

df = CSV.read(IOBuffer(""), transpose=true, header=false)
@test size(df) == (0, 0)

df = CSV.read(IOBuffer(""), transpose=true, header=Symbol[])
@test size(df) == (0, 0)

# providing empty header vector
df = CSV.read(IOBuffer("x\nabc\n"), header=Symbol[])
@test names(df) == [:Column1]

# Union{Bool, Missing}
df = CSV.read(IOBuffer("x\ntrue\n\n"))
@test size(df) == (2, 1)
@test df.x[1] === true
@test df.x[2] === missing

# Union{Date, Missing}
df = CSV.read(IOBuffer("x\n2019-01-01\n\n"))
@test size(df) == (2, 1)
@test df.x[1] === Date(2019, 1, 1)
@test df.x[2] === missing

# use_mmap=false
df = CSV.read(IOBuffer("x\n2019-01-01\n\n"), use_mmap=false)
@test size(df) == (2, 1)
@test df.x[1] === Date(2019, 1, 1)
@test df.x[2] === missing

# types is Dict{String, Type}
df = CSV.read(IOBuffer("x\n2019-01-01\n\n"), types=Dict("x"=>Date))
@test size(df) == (2, 1)
@test df.x[1] === Date(2019, 1, 1)
@test df.x[2] === missing

# various CSV.File/CSV.Row properties
f = CSV.File(IOBuffer("int,float,date,datetime,bool,null,str,catg,int_float\n1,3.14,2019-01-01,2019-01-01T01:02:03,true,,hey,abc,2\n2,NaN,2019-01-02,2019-01-03T01:02:03,false,,there,abc,3.14\n"), pool=0.3)
@test Tables.istable(f)
@test Tables.rowaccess(typeof(f))
@test Tables.columnaccess(typeof(f))
@test Tables.schema(f) == Tables.Schema([:int, :float, :date, :datetime, :bool, :null, :str, :catg, :int_float], [Int64, Float64, Date, DateTime, Bool, Missing, String, String, Float64])
@test Tables.rows(f) === f
@test eltype(f) <: CSV.Row
row = first(f)
@test propertynames(row) == [:int, :float, :date, :datetime, :bool, :null, :str, :catg, :int_float]
@test row.int == 1
@test row.float == 3.14
@test row.date == Date(2019, 1, 1)
@test row.datetime == DateTime(2019, 1, 1, 1, 2, 3)
@test row.bool == true
@test row.null === missing
@test row.str == "hey"
@test row.catg == "abc"
@test typeof(row.catg) == String
@test row.int_float === 2.0
row = iterate(f, 2)[1]
@test row.int_float === 3.14

# 448
@test_throws ArgumentError CSV.File(IOBuffer("x\n1\n2\n3\n#4"), ignorerepeated=true)

# reported by oxinabox on slack; issue w/ manually specified pool column type and 0 rows
df = CSV.read(IOBuffer("x\n"), types=[CSV.PooledString], copycols=true)
@test size(df) == (0, 1)

df = CSV.read(IOBuffer("x\n"), types=[Union{CSV.PooledString, Missing}], copycols=true)
@test size(df) == (0, 1)

f = CSV.File(IOBuffer("x\n1\n2\n3\n#4"), comment="#")
@test length(f.x) == 3
@test f.x[end] == 3

# 453
@test_throws ArgumentError CSV.File(IOBuffer("x\n1\n2\n3\n#4"), types=[CSV_Foo])
@test_throws ArgumentError CSV.File(IOBuffer("x\n1\n2\n3\n#4"), types=Dict(:x=>CSV_Foo))

# 447
df = CSV.read(IOBuffer("a,b,c\n1,2,3\n\n"), ignoreemptylines=true)
@test size(df) == (1, 3)

df = CSV.read(IOBuffer("zip\n11111-1111\n"), dateformat = "y-m-dTH:M:S.s")
@test size(df) == (1, 1)

# Supporting commands across multiple platforms cribbed from julia/test/spawn.jl
catcmd = `cat`
havebb = false
if Sys.iswindows()
    busybox = download("https://frippery.org/files/busybox/busybox.exe", joinpath(tempdir(), "busybox.exe"))
    havebb = try # use busybox-w32 on windows, if available
        success(`$busybox`)
        true
    catch
        false
    end
    if havebb
        catcmd = `$busybox cat`
    end
end
@test CSV.read(`$(catcmd) $(joinpath(dir, "test_basic.csv"))`) == CSV.read(joinpath(dir, "test_basic.csv"))

#476
df = CSV.read(GzipDecompressorStream(open(joinpath(dir, "randoms.csv.gz"))))
@test size(df) == (70000, 7)

# int sentinel cycling
import Random; Random.seed!(1234)
df = CSV.read(joinpath(dir, "test_int_sentinel.csv"))
@test df.id.sentinel == 0x6690d6bed2a7da46
@test df.id[end] == CSV.INT_SENTINEL
@test df.id[end-1] === missing

Random.seed!(1234)
df = CSV.read(joinpath(dir, "test_int_sentinel.csv"), threaded=true)
@test df.id[end] == CSV.INT_SENTINEL
@test df.id[end-1] === missing
@test df.id[50] == 4495412479959714370

df = CSV.read(IOBuffer("thistime\n10:00:00.0\n12:00:00.0"))
@test typeof(df.thistime[1]) <: Dates.Time
@test df.thistime[1] === Time(10)

end


text:
module GlobalDiffEq

using Reexport
@reexport using DiffEqBase

import Distributions, LinearAlgebra, ModelingToolkit, OrdinaryDiffEq, Richardson

include("vector_expectation.jl")
include("adjoint_sol_stepsize.jl")

abstract type GlobalDiffEqAlgorithm <: DiffEqBase.AbstractODEAlgorithm end

struct GlobalRichardson{A} <: GlobalDiffEqAlgorithm
    alg::A
end

function DiffEqBase.__solve(prob::Union{DiffEqBase.AbstractODEProblem,DiffEqBase.AbstractDAEProblem},
                            alg::GlobalRichardson, args...;
                            dt,kwargs...)
    opt = Dict(kwargs)
    otheropts = delete!(copy(opt), :dt)
    tstops = get(opt, :tstops, range(prob.tspan[1], stop=prob.tspan[2], step=dt))
    local sol
    val, err = Richardson.extrapolate(dt, rtol=get(opt, :reltol, 1e-3), atol=get(opt, :abstol, 1e-6), contract=0.5) do _dt
        sol = solve(prob, alg.alg, args...; dt=_dt, adaptive=false, otheropts...)
        sol.(tstops)
    end
    return sol
end

struct GlobalAdjoint{A} <: GlobalDiffEqAlgorithm
    alg::A
end

function DiffEqBase.__solve(prob::Union{DiffEqBase.AbstractODEProblem,DiffEqBase.AbstractDAEProblem},
                            alg::GlobalAdjoint, args...;
                            dt,kwargs...)

    opt = Dict(kwargs)
    otheropts = delete!(copy(opt), :dt)
    tstops = get(opt, :tstops, range(prob.tspan[1], stop=prob.tspan[2], step=dt))
    local sol

    sol_init = solve(prob, alg.alg, args...; dt=dt, otheropts...)
    n = length(prob.u0)
    C = 11/6 # for k = 2, 2 random vectors
    E₂ = 2/π
    Eₙ = vec_expt_full(n)
    rtol = get(opt, :reltol, 1e-3)
    atol = get(opt, :abstol, 1e-6)
    gtol = atol # target gtol is abstol? abstol cannot be zero?
    z₁, z₂ = orth_vec(n)

    t_ctrl = tstops[1]
    tstops_ctrl = [t_ctrl]
    λ₁ = adjoint_solve(prob, sol_init, z₁, tstops, prob.p, rtol, atol)
    λ₂ = adjoint_solve(prob, sol_init, z₂, tstops, prob.p, rtol, atol)
    while t_ctrl < tstops[end]
        h = adjoint_step_ctrl(λ₁(t_ctrl), λ₂(t_ctrl), E₂, Eₙ, C, gtol)
        t_ctrl += h
        push!(tstops_ctrl, t_ctrl)
    end
    push!(tstops_ctrl, tstops[end])

    sol = solve(prob, alg.alg, args...; dt=dt, tstops=tstops_ctrl, adaptive=false, otheropts...) # control stepsize or LTOL?
    return sol
end

export GlobalRichardson, GlobalAdjoint

end


text:
using DataStructures: CircularBuffer
using SHA

mutable struct TierrianOrganism
    key::UInt64

    a::UInt16
    b::UInt16
    c::UInt16
    d::UInt16

    stack::CircularBuffer{UInt16}

    ip::UInt16

    error_flag::Bool

    start_address::UInt16
    length::UInt16

    has_daughter::Bool
    daughter_address::UInt16
    daughter_length::UInt16

    hash::String
    parent_hash::String

    daughters::Vector{TierrianOrganism}

    time_birth::UInt64
    age::UInt16

    function TierrianOrganism(start_address::UInt16, length::UInt16)
        new(0, 0, 0, 0, 0, CircularBuffer{UInt16}(10), start_address, false, start_address, length, false, 0, 0, "", "", [], 0, 0)
    end
end

function _address_within(organism, address)
    mod(address - organism.start_address, organism.length) + organism.start_address
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
"""
```
compute_meansbands(m, input_type, cond_type, output_vars; forecast_string = "",
    verbose = :low, kwargs...)

compute_meansbands(m, input_type, cond_type, output_var, df; forecast_string = "",
    population_data = DataFrame(), population_forecast = DataFrame(),
    verbose = :none, kwargs...)

compute_meansbands(m, input_type, cond_type, output_var, var_name, df;
    forecast_string = "", population_data = DataFrame(),
    population_forecast = DataFrame(), verbose = :low, kwargs...)
```

Compute means and bands for pseudo-observables, observables, and shocks, and write
the results to a file. Other methods are for one `output_var` and one `var_name` respectively.

### Keyword Arguments

- `forecast_string::String`: forecast identifier string (the value
  \"fcid=value\" in the forecast output filename). Required when
  `input_type == :subset`

- `density_bands::Vector{Float64}`: a vector of percent values (between 0 and 1) for
  which to compute density bands

- `minimize::Bool`: if `true`, choose shortest interval, otherwise just chop off
  lowest and highest (percent/2)

- `verbose`: level of error messages to be printed to screen. One of `:none`,
  `:low`, `:high`
"""
function compute_meansbands(m::AbstractDSGEModel, input_type::Symbol,
                            cond_type::Symbol, output_vars::Vector{Symbol};
                            forecast_string::String = "",
                            verbose::Symbol = :low, df::DataFrame = DataFrame(),
                            check_empty_columns::Bool = true,
                            kwargs...)

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        output_dir = workpath(m, "forecast")
        println()
        @Base.info "Computing means and bands for input_type = $input_type, cond_type = $cond_type..."
        println("Start time: $(now())")
        println("Means and bands will be saved in $output_dir")
    end
    elapsed_time = @elapsed let
        # Determine full set of output_vars necessary for plotting desired result
        output_vars = add_requisite_output_vars(output_vars)
        if input_type == :prior
            output_vars = setdiff(output_vars, [:bddforecastobs])
        end
        # Load population data and main dataset (required for some transformations)
        if all(var -> get_product(var) == :irf, output_vars)
            population_data, population_forecast = DataFrame(), DataFrame()
        else
            population_data, population_forecast = load_population_growth(m, verbose = verbose)
            isempty(df) && (df = load_data(m, check_empty_columns = check_empty_columns, verbose = :none))
        end
        for output_var in output_vars
            prod = get_product(output_var)
            if VERBOSITY[verbose] >= VERBOSITY[:high]
                if prod in [:shockdec, :irf]
                    println("Computing " * string(output_var) * " for shocks:")
                else
                    print("Computing " * string(output_var) * "... ")
                end
            end

            # Compute means and bands
            mb = compute_meansbands(m, input_type, cond_type, output_var, df;
                                    forecast_string = forecast_string,
                                    population_data = population_data,
                                    population_forecast = population_forecast,
                                    verbose = verbose,
                                    kwargs...)
            GC.gc()
        end
    end
    if VERBOSITY[verbose] >= VERBOSITY[:low]
        total_mb_time     = elapsed_time
        total_mb_time_min = total_mb_time/60

        println("\nTotal time to compute means and bands: " * string(total_mb_time_min) * " minutes")
        println("Computation of means and bands complete: " * string(now()))
    end
end

function compute_meansbands(m::AbstractDSGEModel, input_type::Symbol, cond_type::Symbol,
                            output_var::Symbol, df::DataFrame;
                            forecast_string::String = "",
                            population_data::DataFrame = DataFrame(),
                            population_forecast::DataFrame = DataFrame(),
                            verbose::Symbol = :none,
                            kwargs...)

    # Determine class and product
    class   = get_class(output_var)
    product = get_product(output_var)

    # Read in forecast metadata
    metadata = get_mb_metadata(m, input_type, cond_type, output_var; forecast_string = forecast_string)

    date_list      = product == :irf ? Date[] : collect(keys(metadata[:date_inds]))
    variable_names = collect(keys(metadata[:indices]))
    pop_growth     = get_mb_population_series(product, population_data, population_forecast, date_list)

    # Compute means and bands
    if product in [:hist, :histut, :hist4q, :forecast, :forecastut, :forecast4q,
                   :bddforecast, :bddforecastut, :bddforecast4q, :dettrend, :trend]
        # Get to work!
        # pmap produces an error for trendobs sometimes, so just doing this iteratively
        mb_vec = Vector{Any}(undef,length(variable_names))
        for i in 1:length(mb_vec)
            mb_vec[i] = compute_meansbands(m, input_type, cond_type, output_var,
                                           variable_names[i], df; pop_growth = pop_growth,
                                           forecast_string = forecast_string, kwargs...)
        end
        # mb_vec = pmap(var_name -> compute_meansbands(m, input_type, cond_type, output_var, var_name, df;
        #                               pop_growth = pop_growth, forecast_string = forecast_string, kwargs...),
        #               variable_names)

        # Re-assemble pmap outputs
        means = DataFrame(date = date_list)
        bands = Dict{Symbol,DataFrame}()

        for (var_name, (var_means, var_bands)) in zip(variable_names, mb_vec)
            means[!,var_name] = var_means
            bands[var_name] = var_bands
            bands[var_name][!,:date] = date_list
        end

    elseif product in [:shockdec, :irf]
        means = product == :irf ? DataFrame() : DataFrame(date = date_list)
        bands = Dict{Symbol, DataFrame}()

        # Get to work!
        for shock_name in keys(metadata[:shock_indices])
            println(verbose, :high, "  * " * string(shock_name))

            mb_vec = pmap(var_name -> compute_meansbands(m, input_type, cond_type, output_var, var_name, df;
                                          pop_growth = pop_growth, shock_name = Nullables.Nullable(shock_name),
                                          forecast_string = forecast_string, kwargs...),
                          variable_names)

            # Re-assemble pmap outputs
            for (var_name, (var_means, var_bands)) in zip(variable_names, mb_vec)
                means[!, Symbol(var_name, DSGE_SHOCKDEC_DELIM, shock_name)] = var_means
                bands[Symbol(var_name, DSGE_SHOCKDEC_DELIM, shock_name)] = var_bands
                if product != :irf
                    bands[Symbol(var_name, DSGE_SHOCKDEC_DELIM, shock_name)][!, :date] = date_list
                end
            end
        end

    else
        error("Invalid product: $product")
    end # of if product

    mb = MeansBands(metadata, means, bands)

    # Write to file
    filepath = get_meansbands_output_file(m, input_type, cond_type, output_var,
                                          forecast_string = forecast_string)
    dirpath = dirname(filepath)
    isdir(dirpath) || mkpath(dirpath)
    JLD2.jldopen(filepath, true, true, true, IOStream) do file
        write(file, "mb", mb)
    end

    sep = prod in [:shockdec, :irf] ? "  " : ""
    println(verbose, :high, sep * "wrote " * basename(filepath))

    return mb
end

function compute_meansbands(m::AbstractDSGEModel, input_type::Symbol, cond_type::Symbol,
                            output_var::Symbol, var_name::Symbol, df::DataFrame;
                            forecast_string::String = "",
                            pop_growth::AbstractVector{Float64} = Float64[],
                            shock_name::Nullable{Symbol} = Nullables.Nullable{Symbol}(),
                            density_bands::Vector{Float64} = [0.5,0.6,0.7,0.8,0.9],
                            minimize::Bool = false,
                            compute_shockdec_bands::Bool = false)

    # Return only one set of bands if we read in only one draw
    if input_type in [:init, :mode, :mean]
        density_bands = [.5]
    end

    # Determine class and product
    class = get_class(output_var)
    product = get_product(output_var)

    # Read in forecast draws
    fcast_series, transform = read_forecast_output(m, input_type, cond_type,
                                                   output_var, var_name, shock_name,
                                                   forecast_string = forecast_string)

    # Reverse transform
    y0_index = get_y0_index(m, product)
    yt_index = get_yt_index(m, product)
    data = class == :obs && product != :irf ? Float64.(collect(Missings.replace(Vector{Union{Missing, Float64}}(df[!,var_name]), NaN))) : fill(NaN, size(df, 1))
    # data = class == :obs && product != :irf ? Float64.(collect(Missings.replace(df[:,var_name], NaN))) : fill(NaN, size(df, 1))
    transformed_series = mb_reverse_transform(fcast_series, transform, product, class,
                                              y0_index = y0_index, yt_index = yt_index,
                                              data = data,
                                              pop_growth = pop_growth)

    # Compute means and bands
    means = vec(mean(transformed_series, dims= 1))
    bands = if product in [:shockdec, :dettrend, :trend] && !compute_shockdec_bands
        Dict{Symbol,DataFrame}()
    else
        find_density_bands(transformed_series, density_bands, minimize = minimize)
    end
    return means, bands
end

function mb_reverse_transform(fcast_series::AbstractArray, transform::Function,
                              product::Symbol, class::Symbol;
                              y0_index::Int = -1, yt_index::Int = -1,
                              data::AbstractVector{Float64} = Float64[],
                              pop_growth::AbstractVector{Float64} = Float64[])
    # No transformation
    if product in [:histut, :forecastut, :bddforecastut]
        return fcast_series
    end

    if product in [:hist4q, :forecast4q, :bddforecast4q]
        transform4q = get_transform4q(transform)
        use_data = class == :obs && product != :hist4q

        y0s = if use_data && transform4q in [loggrowthtopct_4q_percapita, loggrowthtopct_4q]
            # Sum growth rates y_{t-3}, y_{t-2}, y_{t-1}, and y_t
            data[y0_index+1:yt_index] #end]
        elseif use_data && transform4q in [logleveltopct_4q_percapita, logleveltopct_4q]
            # Divide log levels y_t by y_{t-4}
            data[y0_index:yt_index] #end]
        else
            Float64[]
        end

        reverse_transform(fcast_series, transform4q;
                          fourquarter = true, y0s = y0s,
                          pop_growth = pop_growth)
    else
        # Use transformation that doesn't add back population growth for
        # products which are given in deviations
        if product in [:shockdec, :dettrend, :irf]
            transform = get_nopop_transform(transform)
            y0 = 0.0
        else
            y0 = class == :obs ? data[y0_index] : NaN
        end

        reverse_transform(fcast_series, transform;
                          y0 = y0, pop_growth = pop_growth)
    end
end


text:
# Fusion trees:
#==============================================================================#
struct FusionTree{G<:Sector,N,M,L,T}
    uncoupled::NTuple{N,G}
    coupled::G
    innerlines::NTuple{M,G} # M = N-2
    vertices::NTuple{L,T} # L = N-1
end
FusionTree{G}(uncoupled::NTuple{N,Any},
                coupled,
                innerlines,
                vertices = ntuple(n->nothing, StaticLength(N)-StaticLength(1))
                ) where {G<:Sector,N} =
    fusiontreetype(G, StaticLength(N))(map(s->convert(G,s),uncoupled),
        convert(G,coupled), map(s->convert(G,s), innerlines), vertices)
FusionTree(uncoupled::NTuple{N,G},
            coupled::G,
            innerlines,
            vertices = ntuple(n->nothing, StaticLength(N)-StaticLength(1))
            ) where {G<:Sector,N} =
    fusiontreetype(G, StaticLength(N))(uncoupled, coupled, innerlines, vertices)

function FusionTree{G}(uncoupled::NTuple{N}, coupled = one(G)) where {G<:Sector, N}
    FusionStyle(G) isa Abelian ||
        error("fusion tree requires inner lines if `FusionStyle(G) <: NonAbelian`")
    FusionTree{G}(map(s->convert(G,s), uncoupled), convert(G, coupled),
                    _abelianinner(map(s->convert(G,s),(uncoupled..., dual(coupled)))))
end
function FusionTree(uncoupled::NTuple{N,G}, coupled::G = one(G)) where {G<:Sector, N}
    FusionStyle(G) isa Abelian ||
        error("fusion tree requires inner lines if `FusionStyle(G) <: NonAbelian`")
    FusionTree{G}(uncoupled, coupled, _abelianinner((uncoupled..., dual(coupled))))
end

# Properties
sectortype(::Type{<:FusionTree{G}}) where {G<:Sector} = G
FusionStyle(::Type{<:FusionTree{G}}) where {G<:Sector} = FusionStyle(G)
Base.length(::Type{<:FusionTree{<:Sector,N}}) where {N} = N

sectortype(t::FusionTree) = sectortype(typeof(t))
FusionStyle(t::FusionTree) = FusionStyle(typeof(t))
Base.length(t::FusionTree) = length(typeof(t))

# Hashing, important for using fusion trees as key in a dictionary
function Base.hash(f::FusionTree{G}, h::UInt) where {G}
    if FusionStyle(G) isa Abelian
        hash(f.uncoupled, hash(f.coupled, h))
    elseif FusionStyle(G) isa SimpleNonAbelian
        hash(f.innerlines, hash(f.uncoupled, hash(f.coupled, h)))
    else
        hash(f.vertices, hash(f.innerlines, hash(f.uncoupled, hash(f.coupled, h))))
    end
end
function Base.isequal(f1::FusionTree{G,N}, f2::FusionTree{G,N}) where {G,N}
    f1.coupled == f2.coupled || return false
    @inbounds for i = 1:N
        f1.uncoupled[i] == f2.uncoupled[i] || return false
    end
    if FusionStyle(G) isa SimpleNonAbelian
        @inbounds for i=1:N-2
            f1.innerlines[i] == f2.innerlines[i] || return false
        end
    end
    if FusionStyle(G) isa DegenerateNonAbelian
        @inbounds for i=1:N-1
            f1.vertices[i] == f2.vertices[i] || return false
        end
    end
    return true
end
Base.isequal(f1::FusionTree{G1,N1}, f2::FusionTree{G2,N2}) where {G1,G2,N1,N2} = false

# Facilitate getting correct fusion tree types
Base.@pure fusiontreetype(::Type{G}, ::StaticLength{0}) where {G<:Sector} =
    FusionTree{G, 0, 0, 0, vertex_labeltype(G)}
Base.@pure fusiontreetype(::Type{G}, ::StaticLength{1}) where {G<:Sector} =
    FusionTree{G, 1, 0, 0, vertex_labeltype(G)}
Base.@pure fusiontreetype(::Type{G}, ::StaticLength{2}) where {G<:Sector} =
    FusionTree{G, 2, 0, 1, vertex_labeltype(G)}
Base.@pure fusiontreetype(::Type{G}, ::StaticLength{N}) where {G<:Sector, N} =
    _fusiontreetype(G, StaticLength(N),
        StaticLength(N) - StaticLength(2), StaticLength(N) - StaticLength(1))
Base.@pure _fusiontreetype(::Type{G}, ::StaticLength{N}, ::StaticLength{M},
                            ::StaticLength{L}) where {G<:Sector, N, M, L} =
    FusionTree{G,N,M,L,vertex_labeltype(G)}

# converting to actual array
function Base.convert(::Type{Array}, f::FusionTree{G, 0}) where {G}
    T = eltype(fusiontensor(one(G), one(G), one(G)))
    return fill(one(T), 1)
end
function Base.convert(::Type{Array}, f::FusionTree{G, 1}) where {G}
    T = eltype(fusiontensor(one(G), one(G), one(G)))
    return copyto!(Matrix{T}(undef, dim(f.coupled), dim(f.coupled)), I)
end
Base.convert(::Type{Array}, f::FusionTree{G,2}) where {G} =
    fusiontensor(f.uncoupled[1], f.uncoupled[2], f.coupled, f.vertices[1])

function Base.convert(::Type{Array}, f::FusionTree{G}) where {G}
    tailout = (f.innerlines[1], TupleTools.tail2(f.uncoupled)...)
    ftail = FusionTree(tailout, f.coupled, Base.tail(f.innerlines), Base.tail(f.vertices))
    Ctail = convert(Array, ftail)
    C1 = fusiontensor(f.uncoupled[1], f.uncoupled[2], f.innerlines[1], f.vertices[1])
    dtail = size(Ctail)
    d1 = size(C1)
    C = reshape(C1, d1[1]*d1[2], d1[3]) *
            reshape(Ctail, dtail[1], prod(Base.tail(dtail)))
    return reshape(C, (d1[1], d1[2], Base.tail(dtail)...))
end

# Show methods
function Base.show(io::IO, t::FusionTree{G,N,M,K,Nothing}) where {G<:Sector,N,M,K}
    print(IOContext(io, :typeinfo => G), "FusionTree{", G, "}(",
        t.uncoupled, ", ", t.coupled, ", ", t.innerlines, ")")
end
function Base.show(io::IO, t::FusionTree{G}) where {G<:Sector}
    print(IOContext(io, :typeinfo => G), "FusionTree{", G, "}(",
        t.uncoupled, ", ", t.coupled, ", ", t.innerlines, ", ", t.vertices, ")")
end

# Manipulate fusion trees
include("manipulations.jl")

# Fusion tree iterators
include("iterator.jl")

# auxiliary routines
# _abelianinner: generate the inner indices for given outer indices in the abelian case
_abelianinner(outer::Tuple{}) = ()
_abelianinner(outer::Tuple{G}) where {G<:Sector} =
    outer[1] == one(G) ? () : throw(SectorMismatch())
_abelianinner(outer::Tuple{G,G}) where {G<:Sector} =
    outer[1] == dual(outer[2]) ? () : throw(SectorMismatch())
_abelianinner(outer::Tuple{G,G,G}) where {G<:Sector} =
    first(⊗(outer...)) == one(G) ? () : throw(SectorMismatch())
function _abelianinner(outer::NTuple{N,G}) where {G<:Sector,N}
    c = first(outer[1] ⊗ outer[2])
    return (c, _abelianinner((c, TupleTools.tail2(outer)...))...)
end


text:
import HSARuntime: hsa_memory_allocate, hsa_memory_free, check, DEFAULT_AGENT

mutable struct ROCArray{T,N} <: AbstractArray{T,N}
    size::Dims{N}
    handle::Ptr{T}
end

const ROCVector{T} = ROCArray{T,1}
const ROCMatrix{T} = ROCArray{T,2}
const ROCVecOrMat{T} = Union{ROCVector{T},ROCMatrix{T}}

# TODO: Support non-isbitstype allocations
function ROCArray(agent::HSAAgent, ::Type{T}, size::NTuple{N,Int}) where {T,N}
    @assert isbitstype(T) "$T is not a primitive type"
    @assert all(x->x>0, size) "Invalid array size: $size"
    region = get_region(agent, :finegrained)
    nbytes = sizeof(T) * prod(size)
    handle = Ref{Ptr{T}}()
    hsa_memory_allocate(region[], nbytes, handle) |> check
    arr = ROCArray{T,N}(size, handle[])
    finalizer(arr) do arr
        hsa_memory_free(arr.handle) |> check
    end
    return arr
end
ROCArray(::Type{T}, size::NTuple{N,Int}) where {T,N} =
    ROCArray(DEFAULT_AGENT[], T, size)

function ROCArray(agent::HSAAgent, arr::Array{T,N}) where {T,N}
    rarr = ROCArray(agent, T, size(arr))
    for idx in eachindex(arr)
        rarr[idx] = arr[idx]
    end
    return rarr
end
ROCArray(arr::Array{T,N}) where {T,N} =
    ROCArray(DEFAULT_AGENT[], arr)

function Array(rarr::ROCArray{T,N}) where {T,N}
    arr = Array{T}(undef, size(rarr))
    # FIXME: Use Mem
    ref_arr = Ref(arr)
    GC.@preserve ref_arr begin
        ccall(:memcpy, Cvoid,
            (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
            ref_arr, rarr.handle, sizeof(arr))
    end
    return rarr
end

Base.pointer(arr::ROCArray) = arr.handle
Base.cconvert(::Type{Ptr{T}}, x::ROCArray{T}) where T = x.handle
Base.cconvert(::Type{Ptr{Nothing}}, x::ROCArray) = x.handle
Base.IndexStyle(::Type{<:ROCArray}) = Base.IndexLinear()
Base.IndexStyle(::ROCArray) = Base.IndexLinear()
function Base.iterate(A::ROCArray, i=1) # copy-pasta from Base
    Base.@_inline_meta
    (i % UInt) - 1 < length(A) ? (@inbounds A[i], i + 1) : nothing
end
Base.similar(arr::ROCArray{T,N}) where {T,N} =
    ROCArray(T, size(arr))
Base.similar(agent::HSAAgent, arr::ROCArray{T,N}) where {T,N} =
    ROCArray(agent, T, size(arr))
Base.similar(arr::ROCArray{T1,N}, ::Type{T2}, dims::Dims) where {T1,N,T2} =
    similar(DEFAULT_AGENT[], arr, T2, dims)
function Base.similar(agent::HSAAgent, arr::ROCArray{T1,N}, ::Type{T2}, dims::Dims) where {T1,N,T2}
    ROCArray(agent, T2, dims)
end

# copy-pasta from Base
function Base.stride(arr::ROCArray, i::Int)
    if i > ndims(arr)
        return length(ar)
    end
    s = 1
    for n = 1:(i-1)
        s *= size(arr, n)
    end
    return s
end

Base.size(arr::ROCArray) = arr.size
Base.length(arr::ROCArray) = prod(size(arr))

function Base.fill!(arr::ROCArray{T,N}, value::T) where {T,N}
    for idx in 1:length(arr)
        arr[idx] = value
    end
end
@inline function Base.getindex(arr::ROCArray{T,N}, idx) where {T,N}
    @boundscheck checkbounds(arr, idx)
    Base.unsafe_load(pointer(arr), idx)::T
end
@inline function Base.setindex!(arr::ROCArray{T,N}, value, idx) where {T,N}
    @boundscheck checkbounds(arr, idx)
    Base.unsafe_store!(pointer(arr), value, idx)
end

#=

mutable struct ROCArray{T,N} <: GPUArray{T,N}
  buf::Mem.Buffer
  own::Bool

  dims::Dims{N}
  offset::Int

  function ROCArray{T,N}(buf::Mem.Buffer, dims::Dims{N}; offset::Integer=0, own::Bool=true) where {T,N}
    xs = new{T,N}(buf, own, dims, offset)
    if own
      Mem.retain(buf)
      finalizer(unsafe_free!, xs)
    end
    return xs
  end
end

function unsafe_free!(xs::ROCArray)
  Mem.release(xs.buf) && dealloc(xs.buf, prod(xs.dims)*sizeof(eltype(xs)))
  return
end

## construction

# type and dimensionality specified, accepting dims as tuples of Ints
ROCArray{T,N}(::UndefInitializer, dims::Dims{N}) where {T,N} =
  ROCArray{T,N}(alloc(prod(dims)*sizeof(T)), dims)

# type and dimensionality specified, accepting dims as series of Ints
ROCArray{T,N}(::UndefInitializer, dims::Integer...) where {T,N} = ROCArray{T,N}(undef, dims)

# type but not dimensionality specified
ROCArray{T}(::UndefInitializer, dims::Dims{N}) where {T,N} = ROCArray{T,N}(undef, dims)
ROCArray{T}(::UndefInitializer, dims::Integer...) where {T} =
  ROCArray{T}(undef, convert(Tuple{Vararg{Int}}, dims))

# empty vector constructor
ROCArray{T,1}() where {T} = ROCArray{T,1}(undef, 0)

Base.similar(a::ROCArray{T,N}) where {T,N} = ROCArray{T,N}(undef, size(a))
Base.similar(a::ROCArray{T}, dims::Base.Dims{N}) where {T,N} = ROCArray{T,N}(undef, dims)
Base.similar(a::ROCArray, ::Type{T}, dims::Base.Dims{N}) where {T,N} = ROCArray{T,N}(undef, dims)

"""
  unsafe_wrap(::ROCArray, pointer{T}, dims; own=false, agent::HSAAgent)

Wrap a `ROCArray` object around the data at the address given by `pointer`. The pointer
element type `T` determines the array element type. `dims` is either an integer (for a 1d
array) or a tuple of the array dimensions. `own` optionally specified whether Julia should
take ownership of the memory, calling `free` when the array is no longer referenced. The
`agent` argument determines the ROC agent where the data is allocated in.
"""
function Base.unsafe_wrap(::Union{Type{ROCArray},Type{ROCArray{T}},Type{ROCArray{T,N}}},
                          p::Ptr{T}, dims::NTuple{N,Int};
                          own::Bool=false, agent::HSAAgent=HSARuntime.get_default_agent()) where {T,N}
  buf = Mem.Buffer(convert(Ptr{Cvoid}, p), prod(dims) * sizeof(T), agent)
  return ROCArray{T, length(dims)}(buf, dims; own=own)
end
function Base.unsafe_wrap(Atype::Union{Type{ROCArray},Type{ROCArray{T}},Type{ROCArray{T,1}}},
                          p::Ptr{T}, dim::Integer;
                          own::Bool=false, agent::HSAAgent=HSARuntime.get_default_agent()) where {T}
  unsafe_wrap(Atype, p, (dim,); own=own, agent=agent)
end

## array interface

Base.elsize(::Type{<:ROCArray{T}}) where {T} = sizeof(T)

Base.size(x::ROCArray) = x.dims
Base.sizeof(x::ROCArray) = Base.elsize(x) * length(x)

## interop with other arrays

ROCArray{T,N}(xs::AbstractArray{T,N}) where {T,N} =
  isbits(xs) ?
    (ROCArray{T,N}(undef, size(xs)) .= xs) :
    copyto!(ROCArray{T,N}(undef, size(xs)), collect(xs))

ROCArray{T,N}(xs::AbstractArray{S,N}) where {T,N,S} = ROCArray{T,N}((x -> T(x)).(xs))

# underspecified constructors
ROCArray{T}(xs::AbstractArray{S,N}) where {T,N,S} = ROCArray{T,N}(xs)
(::Type{ROCArray{T,N} where T})(x::AbstractArray{S,N}) where {S,N} = ROCArray{S,N}(x)
ROCArray(A::AbstractArray{T,N}) where {T,N} = ROCArray{T,N}(A)

# idempotency
ROCArray{T,N}(xs::ROCArray{T,N}) where {T,N} = xs

## conversions

Base.convert(::Type{T}, x::T) where T <: ROCArray = x

function Base._reshape(parent::ROCArray, dims::Dims)
  n = length(parent)
  prod(dims) == n || throw(DimensionMismatch("parent has $n elements, which is incompatible with size $dims"))
  return ROCArray{eltype(parent),length(dims)}(parent.buf, dims;
                                              offset=parent.offset, own=parent.own)
end

## interop with C libraries

"""
  buffer(array::ROCArray [, index])

Get the native address of a ROCArray, optionally at a given location `index`.
Equivalent of `Base.pointer` on `Array`s.
"""
function buffer(xs::ROCArray, index=1)
  extra_offset = (index-1) * Base.elsize(xs)
  Mem.Buffer(xs.buf.ptr + xs.offset + extra_offset,
             sizeof(xs) - extra_offset,
             xs.buf.agent)
end

Base.cconvert(::Type{Ptr{T}}, x::ROCArray{T}) where T = buffer(x)
Base.cconvert(::Type{Ptr{Nothing}}, x::ROCArray) = buffer(x)

## interop with AMDGPUnative

function Base.convert(::Type{ROCDeviceArray{T,N,AS.Global}}, a::ROCArray{T,N}) where {T,N}
  ptr = Base.unsafe_convert(Ptr{T}, a.buf)
  ROCDeviceArray{T,N,AS.Global}(a.dims, DevicePtr{T,AS.Global}(ptr+a.offset))
end

Adapt.adapt_storage(::AMDGPUnative.Adaptor, xs::ROCArray{T,N}) where {T,N} =
  convert(ROCDeviceArray{T,N,AS.Global}, xs)

## interop with CPU array

# We don't convert isbits types in `adapt`, since they are already
# considered GPU-compatible.

Adapt.adapt_storage(::Type{<:ROCArray}, xs::AbstractArray) =
  isbits(xs) ? xs : convert(ROCArray, xs)

Adapt.adapt_storage(::Type{<:ROCArray{T}}, xs::AbstractArray{<:Real}) where T <: AbstractFloat =
  isbits(xs) ? xs : convert(ROCArray{T}, xs)

Adapt.adapt_storage(::Type{<:Array}, xs::ROCArray) = convert(Array, xs)

Base.collect(x::ROCArray{T,N}) where {T,N} = copyto!(Array{T,N}(undef, size(x)), x)

function Base.unsafe_copyto!(dest::ROCArray{T}, doffs, src::Array{T}, soffs, n) where T
  Mem.upload!(buffer(dest, doffs), pointer(src, soffs), n*sizeof(T))
  return dest
end

function Base.unsafe_copyto!(dest::Array{T}, doffs, src::ROCArray{T}, soffs, n) where T
  Mem.download!(pointer(dest, doffs), buffer(src, soffs), n*sizeof(T))
  return dest
end

function Base.unsafe_copyto!(dest::ROCArray{T}, doffs, src::ROCArray{T}, soffs, n) where T
  Mem.transfer!(buffer(dest, doffs), buffer(src, soffs), n*sizeof(T))
  return dest
end

function Base.deepcopy_internal(x::ROCArray, dict::IdDict)
  haskey(dict, x) && return dict[x]::typeof(x)
  return dict[x] = copy(x)
end

## utilities

roc(xs) = adapt(ROCArray{Float32}, xs)
Base.getindex(::typeof(roc), xs...) = ROCArray([xs...])

roczeros(T::Type, dims...) = fill!(ROCArray{T}(undef, dims...), 0)
rocones(T::Type, dims...) = fill!(ROCArray{T}(undef, dims...), 1)
roczeros(dims...) = roczeros(Float32, dims...)
rocones(dims...) = rocones(Float32, dims...)
rocfill(v, dims...) = fill!(ROCArray{typeof(v)}(undef, dims...), v)
rocfill(v, dims::Dims) = fill!(ROCArray{typeof(v)}(undef, dims...), v)

# optimized implementation of `fill!` for types that are directly supported by memset
const MemsetTypes = Dict(1=>UInt8, 2=>UInt16, 4=>UInt32)
const MemsetCompatTypes = Union{UInt8, Int8,
                                UInt16, Int16, Float16,
                                UInt32, Int32, Float32}
function Base.fill!(A::ROCArray{T}, x) where T <: MemsetCompatTypes
  y = reinterpret(MemsetTypes[sizeof(T)], convert(T, x))
  Mem.set!(buffer(A), y, length(A))
  A
end

## generic linear algebra routines

function LinearAlgebra.tril!(A::ROCMatrix{T}, d::Integer = 0) where T
  function kernel!(_A, _d)
    li = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    m, n = size(_A)
    if 0 < li <= m*n
      i, j = Tuple(CartesianIndices(_A)[li])
      if i < j - _d
        _A[i, j] = 0
      end
    end
    return nothing
  end

  grid, group = rocdims(A)
  @roc gridsize=grid groupsize=group kernel!(A, d)
  return A
end

function LinearAlgebra.triu!(A::ROCMatrix{T}, d::Integer = 0) where T
  function kernel!(_A, _d)
    li = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    m, n = size(_A)
    if 0 < li <= m*n
      i, j = Tuple(CartesianIndices(_A)[li])
      if j < i + _d
        _A[i, j] = 0
      end
    end
    return nothing
  end

  grid, group = rocdims(A)
  @roc gridsize=grid groupsize=group kernel!(A, d)
  return A
end
=#


text:
using TimeZones.TZData: MIN_OFFSET, MAX_OFFSET

# TimeZone concepts used to disambiguate context of DateTimes
# abstract type UTC <: TimeZone end # Defined in Dates
abstract type Local <: TimeZone end

function transition_range(local_dt::DateTime, tz::VariableTimeZone, ::Type{Local})
    transitions = tz.transitions

    # Determine the earliest and latest possible UTC DateTime
    # that this local DateTime could be.
    # TODO: Alternatively we should only look at the range of offsets available within
    # this TimeZone.
    earliest = local_dt + MIN_OFFSET
    latest = local_dt + MAX_OFFSET

    # Determine the earliest transition the local DateTime could
    # occur within.
    start = searchsortedlast(
        transitions, earliest,
        by=el -> isa(el, Transition) ? el.utc_datetime : el,
    )
    start = max(start, 1)
    finish = length(transitions)
    for i in start:finish
        if transitions[i].utc_datetime > latest
            finish = i - 1
            break
        end
    end

    return start:finish
end

function transition_range(utc_dt::DateTime, tz::VariableTimeZone, ::Type{UTC})
    index = searchsortedlast(
        tz.transitions, utc_dt,
        by=el -> isa(el, Transition) ? el.utc_datetime : el,
    )
    index = max(index, 1)
    return index:index
end

"""
    transition_range(dt::DateTime, tz::VariableTimeZone, context::Type{Union{Local,UTC}}) -> UnitRange

Finds the indexes of the `tz` transitions which may be applicable for the `dt`. The given
DateTime is expected to be local to the time zone or in UTC as specified by `context`. Note
that UTC context will always return a range of length one.
"""
transition_range(::DateTime, ::VariableTimeZone, ::Type{Union{Local,UTC}})

function interpret(local_dt::DateTime, tz::VariableTimeZone, ::Type{Local})
    interpretations = ZonedDateTime[]
    t = tz.transitions
    n = length(t)
    for i in transition_range(local_dt, tz, Local)
        # Convert the local DateTime into UTC
        utc_dt = local_dt - t[i].zone.offset

        if utc_dt >= t[i].utc_datetime && (i == n || utc_dt < t[i + 1].utc_datetime)
            push!(interpretations, ZonedDateTime(utc_dt, tz, t[i].zone))
        end
    end

    return interpretations
end

function interpret(utc_dt::DateTime, tz::VariableTimeZone, ::Type{UTC})
    range = transition_range(utc_dt, tz, UTC)
    length(range) == 1 || error("Internal TimeZones error: A UTC DateTime should only have a single interpretation")
    i = first(range)
    return [ZonedDateTime(utc_dt, tz, tz.transitions[i].zone)]
end

"""
    interpret(dt::DateTime, tz::VariableTimeZone, context::Type{Union{Local,UTC}}) -> Array{ZonedDateTime}

Produces a list of possible `ZonedDateTime`s given a `DateTime` and `VariableTimeZone`.
The result will be returned in chronological order. Note that `DateTime`s in the local
context typically return 0-2 results while the UTC context will always return 1 result.
"""
interpret(::DateTime, ::VariableTimeZone, ::Type{Union{Local,UTC}})

"""
    shift_gap(local_dt::DateTime, tz::VariableTimeZone) -> Array{ZonedDateTime}

Given a non-existent local `DateTime` in a `TimeZone` produces two valid `ZonedDateTime`s
that span the gap. Providing a valid local `DateTime` returns an empty array. Note that this
function does not support passing in a UTC `DateTime` since there are no non-existent UTC
`DateTime`s.

Aside: the function name refers to a period of invalid local time (gap) caused by daylight
saving time or offset changes (shift).
"""
function shift_gap(local_dt::DateTime, tz::VariableTimeZone)
    boundaries = ZonedDateTime[]
    t = tz.transitions
    n = length(t)
    delta = eps(local_dt)
    for i in transition_range(local_dt, tz, Local)
        # Convert the local DateTime into UTC
        utc_dt = local_dt - t[i].zone.offset

        # Essentially: t[i].utc_datetime <= utc_dt < t[i + 1].utc_datetime
        starts_after = utc_dt >= t[i].utc_datetime
        ends_before = i == n || utc_dt < t[i + 1].utc_datetime

        # No boundaries should be produced when the given UTC DateTime exists
        if starts_after && ends_before
            empty!(boundaries)
            break

        # UTC DateTime proceeds the end of the transition range
        elseif !ends_before
            push!(boundaries, ZonedDateTime(t[i + 1].utc_datetime - delta, tz, t[i].zone))

        # UTC DateTime preceeds the start of the transition range
        elseif !starts_after
            push!(boundaries, ZonedDateTime(t[i].utc_datetime, tz, t[i].zone))
        end

        # A slower but much easier to understand version of the above code:
        #
        # if starts_after && ends_before
        #     empty!(boundaries)
        #     break
        # elseif !starts_after
        #     push!(
        #         boundaries,
        #         ZonedDateTime(t[i].utc_datetime - eps(t[i].utc_datetime), tz, from_utc=true),
        #         ZonedDateTime(t[i].utc_datetime, tz, from_utc=true),
        #     )
        # end
    end

    # In time zones with hidden transitions we could end up with more than two "bounds".
    # Note this is more of a theoretical issue and would probably only ever occur with hand-
    # crafted VariableTimeZones.
    if length(boundaries) > 2
        boundaries = [first(boundaries), last(boundaries)]
    end

    # Although we are using an array the only valid output from this function should be an
    # empty array or a 2-element array.
    return boundaries
end

"""
    first_valid(local_dt::DateTime, tz::VariableTimeZone, step::Period)

Construct a valid `ZonedDateTime` by adjusting the local `DateTime`. If the local `DateTime`
is non-existent then it will be adjusted using the `step` to be *after* the gap. When the
local `DateTime` is ambiguous the *first* ambiguous `DateTime` will be returned.
"""
function first_valid(local_dt::DateTime, tz::VariableTimeZone, step::Period)
    possible = interpret(local_dt, tz, Local)

    # Skip all non-existent local datetimes.
    while isempty(possible)
        local_dt += step
        possible = interpret(local_dt, tz, Local)
    end

    return first(possible)
end

"""
    last_valid(local_dt::DateTime, tz::VariableTimeZone, step::Period)

Construct a valid `ZonedDateTime` by adjusting the local `DateTime`. If the local `DateTime`
is non-existent then it will be adjusted using the `step` to be *before* the gap. When the
local `DateTime` is ambiguous the *last* ambiguous `DateTime` will be returned.
"""
function last_valid(local_dt::DateTime, tz::VariableTimeZone, step::Period)
    possible = interpret(local_dt, tz, Local)

    # Skip all non-existent local datetimes.
    while isempty(possible)
        local_dt -= step
        possible = interpret(local_dt, tz, Local)
    end

    return last(possible)
end

function first_valid(local_dt::DateTime, tz::VariableTimeZone)
    possible = interpret(local_dt, tz, Local)
    return isempty(possible) ? last(shift_gap(local_dt, tz)) : first(possible)
end

function last_valid(local_dt::DateTime, tz::VariableTimeZone)
    possible = interpret(local_dt, tz, Local)
    return isempty(possible) ? first(shift_gap(local_dt, tz)) : last(possible)
end


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot(x=[0.0533077, -0.0631488, -0.0132793, -0.0316036, 0.00812427, -0.00512818],
     xmin=[-0.566499, -0.373052, -0.168231, -0.109079, -0.0306136, -0.0244971],
     xmax=[0.673114, 0.246754, 0.141672, 0.0458722, 0.0468622, 0.0142408],
     y=[0.0533077, -0.0631488, -0.0132793, -0.0316036, 0.00812427, -0.00512818],
     ymin=[-0.566499, -0.373052, -0.168231, -0.109079, -0.0306136, -0.0244971],
     ymax=[0.673114, 0.246754, 0.141672, 0.0458722, 0.0468622, 0.0142408],
     color=String["1", "2", "3", "4", "5", "6"],
     Geom.point, Geom.errorbar)


text:
# workspace()
ENV["JuLearnBackend"]="ArrayFire"
cd("/home/astupidbear/.julia/v0.4/JuLearn/src/ArrayFire/examples/rnn/")
using JuLearn,ArrayFire,JLD
@Eval lr=1e-2 batchsize=50

include("predict.jl")
include("dataProvider.jl")

inputBase,outputBase=dataProvider();
T=10
net=rnn([:d,:h,:h,:l],[1 1;1 10;10 1;1 1],T,:regression);

# backprop
@time err=backprop(net,AFArray(inputBase),AFArray(outputBase),nIter=100)
using PyPlot;plot(err)

# prediction
input=inputBase[:,100+(0:T÷2-1),1]
target=outputBase[:,100+(0:T-1),1]'
prediction=predict(net,AFArray(input))
prediction=map(i->Array(prediction[i])[1],1:T)
plot([target prediction],label=["target" "prediction"]);


text:
module transformations

export transform

import ..particle: transform

using ..polynomials, ..particle, ..operators_dense, ..fock

"""
    transform(b1::PositionBasis, b2::FockBasis; x0=1)
    transform(b1::FockBasis, b2::PositionBasis; x0=1)

Transformation operator between position basis and fock basis.

The coefficients are connected via the relation
```math
ψ(x_i) = \\sum_{n=0}^N ⟨x_i|n⟩ ψ_n
```
where ``⟨x_i|n⟩`` is the value of the n-th eigenstate of a particle in
a harmonic trap potential at position ``x``, i.e.:
```math
⟨x_i|n⟩ = π^{-\\frac{1}{4}} \\frac{e^{-\\frac{1}{2}\\left(\\frac{x}{x_0}\\right)^2}}{\\sqrt{x_0}}
            \\frac{1}{\\sqrt{2^n n!}} H_n\\left(\\frac{x}{x_0}\\right)
```
"""
function transform(b1::PositionBasis, b2::FockBasis; x0::Real=1)
    T = Matrix{ComplexF64}(undef, length(b1), length(b2))
    xvec = samplepoints(b1)
    A = hermite.A(b2.N)
    delta_x = particle.spacing(b1)
    c = 1.0/sqrt(x0*sqrt(pi))*sqrt(delta_x)
    for i in 1:length(b1)
        u = xvec[i]/x0
        C = c*exp(-u^2/2)
        for n=0:b2.N
            T[i,n+1] = C*horner(A[n+1], u)
        end
    end
    DenseOperator(b1, b2, T)
end

function transform(b1::FockBasis, b2::PositionBasis; x0::Real=1)
    T = Matrix{ComplexF64}(undef, length(b1), length(b2))
    xvec = samplepoints(b2)
    A = hermite.A(b1.N)
    delta_x = particle.spacing(b2)
    c = 1.0/sqrt(x0*sqrt(pi))*sqrt(delta_x)
    for i in 1:length(b2)
        u = xvec[i]/x0
        C = c*exp(-u^2/2)
        for n=0:b1.N
            T[n+1,i] = C*horner(A[n+1], u)
        end
    end
    DenseOperator(b1, b2, T)
end

end # module


text:
#! /usr/bin/env julia

using PyPlot

x = linspace(0,2*pi,1000); y = sin(3*x + 4*cos(2*x));
plot(x, y, color="red", linewidth=2.0, linestyle="--")
title("A sinusoidally modulated sinusoid")

max_stars_repo_name: sjkelly/JuliaPlayground

max_issues_repo_name: sjkelly/JuliaPlayground

max_forks_repo_name: sjkelly/JuliaPlayground


text:
import LinearAlgebra: LU, checknonsingular

export QuotientSpace

struct QuotientSpace{S,O,D,R} <: Space{D,R}
    space::S
    bcs::O
    F::LU{R,Matrix{R}}
    x::Vector{R}
    function QuotientSpace{S,O,D,R}(space::S, bcs::O) where {S,O,D,R}
        n = size(bcs, 1)
        @assert isfinite(n)
        new{S,O,D,R}(space, bcs, lu!(Matrix{R}(I, n, n)), zeros(R, n))
    end
end

QuotientSpace(sp::Space{D,R}, bcs::Operator{T}) where {D,R,T} = QuotientSpace{typeof(sp), typeof(bcs), D, promote_type(R, T)}(sp, bcs)
QuotientSpace(bcs::Operator) = QuotientSpace(domainspace(bcs), bcs)

domain(QS::QuotientSpace) = domain(QS.space)
canonicalspace(QS::QuotientSpace) = QS.space

spacescompatible(a::QuotientSpace,b::QuotientSpace) = spacescompatible(a.space,b.space)
hasconversion(a::Space,b::QuotientSpace) = hasconversion(a,b.space)
hasconversion(a::QuotientSpace,b::Space) = hasconversion(a.space,b)
hasconversion(a::QuotientSpace,b::QuotientSpace) = hasconversion(a.space,b)

Conversion(Q::QuotientSpace{S}, sp::S) where S<:Space = ConcreteConversion(Q, sp)
@inline bandwidths(C::ConcreteConversion{QuotientSpace{S,O,D,R},S}) where {S,O,D,R} = (size(C.domainspace.F, 1), 0)

function getindex(C::ConcreteConversion{QuotientSpace{S,O,D,R},S}, i::Integer, j::Integer) where {S,O,D,R}
    sp = domainspace(C)
    B = sp.bcs
    F = sp.F
    A = F.factors
    n = size(A, 1)
    x = sp.x
    if i == j
        return one(R)
    elseif j < i ≤ j+n
        @inbounds for jj = 1:n, ii = 1:n
            A[ii,jj] = B[ii,j+jj]
        end
        @inbounds for ii = 1:n
            x[ii] = -B[ii,j]
        end
        if norm(x) > 8*norm(A)*eps(R)
            mutable_lu!(F)
            ldiv!(F, x)
        end
        return @inbounds x[i-j]
    else
        return zero(R)
    end
end

# This function is modified from Julia's lu.jl in stdlib to work on a single `LU`.
# License is MIT: https://julialang.org/license
function mutable_lu!(F::LU{T}, ::Val{Pivot} = Val(true);
                         check::Bool = true) where {T,Pivot}
    A = F.factors
    m, n = size(A)
    minmn = min(m,n)
    info = 0
    ipiv = F.ipiv
    @inbounds begin
        for k = 1:minmn
            # find index max
            kp = k
            if Pivot
                amax = abs(zero(T))
                for i = k:m
                    absi = abs(A[i,k])
                    if absi > amax
                        kp = i
                        amax = absi
                    end
                end
            end
            ipiv[k] = kp
            if !iszero(A[kp,k])
                if k != kp
                    # Interchange
                    for i = 1:n
                        tmp = A[k,i]
                        A[k,i] = A[kp,i]
                        A[kp,i] = tmp
                    end
                end
                # Scale first column
                Akkinv = inv(A[k,k])
                for i = k+1:m
                    A[i,k] *= Akkinv
                end
            elseif info == 0
                info = k
            end
            # Update the rest
            for j = k+1:n
                for i = k+1:m
                    A[i,j] -= A[i,k]*A[k,j]
                end
            end
        end
    end
    check && checknonsingular(info)
    return F
end

import LinearAlgebra.BLAS.@blasfunc
import LinearAlgebra: chkstride1, BlasInt
import LinearAlgebra.LAPACK.chklapackerror

export PathologicalQuotientSpace

struct PathologicalQuotientSpace{S,O<:Operator,DD,T,RT} <: Space{DD,T}
    space::S
    bcs::O
    A::Matrix{T}
    x::Vector{T}
    b::Vector{T}
    c::Vector{T}
    U::Matrix{T}
    Σ::Diagonal{T,Vector{T}}
    VT::Matrix{T}
    work::Vector{T}
    iwork::Vector{BlasInt}
    rwork::Vector{RT}
    info::Ref{BlasInt}
    function PathologicalQuotientSpace{S,O,DD,T,RT}(space::S, bcs::O) where {S,O,DD,T,RT}
        n = size(bcs, 1)
        @assert isfinite(n)
        A = zeros(T, n, 2n)
        x = zeros(T, 2n)
        b = zeros(T, n)
        c = zeros(T, n)
        U = zeros(T, n, n)
        Σ = Diagonal(zeros(T, n))
        VT = zeros(T, n, 2n)
        work = Vector{T}(undef,max((6n+7)*2n,68))
        iwork = Vector{BlasInt}(undef,16n)
        rwork = Vector{RT}(undef,(10n+7)*2n)
        info = Ref{BlasInt}()
        new{S,O,DD,T,RT}(space, bcs, A, x, b, c, U, Σ, VT, work, iwork, rwork, info)
    end
end

PathologicalQuotientSpace(sp::Space, bcs::Operator) = PathologicalQuotientSpace{typeof(sp), typeof(bcs), domaintype(sp), rangetype(sp), real(rangetype(sp))}(sp, bcs)

PathologicalQuotientSpace(bcs::Operator) = PathologicalQuotientSpace(domainspace(bcs), bcs)

domain(QS::PathologicalQuotientSpace) = domain(QS.space)
dimension(QS::PathologicalQuotientSpace) = dimension(QS.space)

Conversion(Q::PathologicalQuotientSpace{SP}, S::SP) where {SP<:Space} = ConcreteConversion(Q, S)

bandwidths(C::ConcreteConversion{PathologicalQuotientSpace{SP,O,DD,T,RT},SP}) where {SP,O,DD,T,RT} = (2size(C.domainspace.A, 1), 0)

function getindex(C::ConcreteConversion{PathologicalQuotientSpace{SP,O,DD,T,RT},SP}, i::Integer, j::Integer) where {SP,O,DD,T,RT}
    sp = domainspace(C)
    n = size(sp.A, 1)
    A = sp.A
    x = sp.x
    b = sp.b
    c = sp.c
    B = sp.bcs
    U = sp.U
    Σ = sp.Σ
    VT = sp.VT
    work = sp.work
    iwork = sp.iwork
    rwork = sp.rwork
    info = sp.info

    if i == j
        one(T)
    elseif j < i ≤ j+2n
        for jj = 1:2n, ii = 1:n
            A[ii,jj] = B[ii,j+jj]
        end
        for ii = 1:n
            b[ii] = -B[ii,j]
        end
        in_place_gesdd!(A, U, Σ.diag, VT, work, iwork, rwork, info)
        mul!(c, transpose(U), b)
        mul!(b, pinv!(Σ), c)
        mul!(x, transpose(VT), b)
        x[i-j]
    else
        zero(T)
    end
end

function pinv!(D::Diagonal{T}) where T
    d = D.diag
    @inbounds @simd for i = 1:length(d)
        isfinite(inv(d[i])) ? d[i]=inv(d[i]) : d[i]=zero(T)
    end
    D
end

for (gesdd, elty, relty) in ((:dgesdd_,:Float64,:Float64),
                      (:sgesdd_,:Float32,:Float32),
                      (:zgesdd_,:ComplexF64,:Float64),
                      (:cgesdd_,:ComplexF32,:Float32))
    @eval begin
        #    SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
        #                   LWORK, IWORK, INFO )
        #*     .. Scalar Arguments ..
        #      CHARACTER          JOBZ
        #      INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
        #*     ..
        #*     .. Array Arguments ..
        #      INTEGER            IWORK( * )
        #      DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
        #                        VT( LDVT, * ), WORK( * )
        function in_place_gesdd!(A::Matrix{$elty}, U::Matrix{$elty}, S::Vector{$elty}, VT::Matrix{$elty}, work::Vector{$elty}, iwork::Vector{BlasInt}, rwork::Vector{$relty}, info::Ref{BlasInt})
            chkstride1(A)
            m, n   = size(A)
            minmn  = min(m, n)
            job = 'S'
            lwork  = BlasInt(-1)
            cmplx  = $elty != $relty
            for i = 1:2
                if cmplx
                    ccall((@blasfunc($gesdd), liblapack), Nothing,
                          (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty},
                           Ref{BlasInt}, Ptr{$relty}, Ptr{$elty}, Ref{BlasInt},
                           Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                           Ptr{$relty}, Ptr{BlasInt}, Ptr{BlasInt}),
                          job, m, n, A,
                          max(1,stride(A,2)), S, U, max(1,stride(U,2)),
                          VT, max(1,stride(VT,2)), work, lwork,
                          rwork, iwork, info)
                else
                    ccall((@blasfunc($gesdd), liblapack), Nothing,
                          (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty},
                           Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt},
                           Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                           Ptr{BlasInt}, Ptr{BlasInt}),
                          job, m, n, A,
                          max(1,stride(A,2)), S, U, max(1,stride(U,2)),
                          VT, max(1,stride(VT,2)), work, lwork,
                          iwork, info)
                end
                chklapackerror(info[])
                if i == 1
                    # Work around issue with truncated Float32 representation of lwork in
                    # sgesdd by using nextfloat. See
                    # http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4587&p=11036&hilit=sgesdd#p11036
                    # and
                    # https://github.com/scipy/scipy/issues/5401
                    lwork = round(BlasInt, nextfloat(real(work[1])))
                    #work = Vector{$elty}(lwork)
                end
            end
        end
    end
end

function BandedMatrix(S::SubOperator{T,ConcreteConversion{QuotientSpace{SP,O,D,R},SP,T},Tuple{UnitRange{Int},UnitRange{Int}}}) where {SP,O,D,R,T}
    kr,jr = parentindices(S)
    C = parent(S)
    #ret = BandedMatrix{eltype(S)}(undef, size(S), bandwidths(S))
    ret = BandedMatrix{eltype(S)}(undef, (last(kr),last(jr)), bandwidths(C))
    #@assert first(kr) == first(jr) == 1

    sp = domainspace(C)
    F = sp.F
    A = F.factors
    n = size(A, 1)
    B = sp.bcs[1:n,1:last(jr)+n]
    x = sp.x
    @inbounds for j in 1:last(jr)
        kk = colrange(ret, j)
        if j in kk
            ret[j,j] = one(R)
        end
        for jj = 1:n, ii = 1:n
            A[ii,jj] = B[ii,j+jj]
        end
        for ii = 1:n
            x[ii] = -B[ii,j]
        end
        if norm(x) > 8*norm(A)*eps(R)
            mutable_lu!(F)
            ldiv!(F, x)
        end
        for k = first(kk)+1:last(kk)
            ret[k,j] = x[k-j]
        end
    end

    ret[kr,jr]
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function textbox_tests(fname)
    Drawing(1000, 1000, fname)
    origin()
    setopacity(0.8)
    sethue("black")
    fontface("Georgia-Bold")
    tiles = Tiler(1000, 1000, 4, 3, margin=50)
    for (pos, n) in tiles
        fsize = 10
        fontsize(fsize)
        @layer begin
            translate(pos)
            sethue("ivory")
            box(O, tiles.tilewidth-2, tiles.tileheight-2, :fill)
            sethue("black")
            # even-numbered boxes split and count
            if iseven(n)
                opacity = 1.0
                # split on ., (, and )
                textbox(split("Luxor is a city in Upper (southern) Egypt. " ^ 6, r"\.|\(|\)"),
                O - (tiles.tilewidth/2 - 2, tiles.tileheight/2 - 2),
                # l is linenumber, t is the text of a line, p is the position of the start of the line, h is the interline height
                linefunc = (l, t, p, h) ->
                    begin
                        # add length of line at end
                        @layer begin
                            fontsize(8)
                            fontface("Menlo")
                            text(string(length(t)), p + (tiles.tilewidth - 20, 0))
                        end
                    end,
                    leading=10)
            else
                # odd-numbered boxes number and fade
                # split on the letter "y"
               textbox(split("Luxor is a city in Upper (southern) Egypt. " ^ 6, "y"),
                O - (tiles.tilewidth/2 - 2, tiles.tileheight/2 - 2),
                linefunc = (l, t, p, h) ->
                    begin
                        text(string(l), p + (tiles.tilewidth - 20, 0))
                        setopacity(rescale(l, 1, 10, 1.0, 0.3))
                    end,
                    leading=12)
            end
        end
    end

    fontsize(80)
    textbox("LUXOR", alignment=:center)

    @test finish() == true
    println("...finished test: output in $(fname)")
end

textbox_tests("textbox-tests.png")


text:
module Potentials

export PowerLawPotential, mass, density, potential

include("constants.jl")
using .Constants
c = Constants

abstract Potential
abstract SphericalPotential <: Potential

function density(phi::SphericalPotential, r::Real)
    error("Must use subtype of SphericalPotential with defined density!")
end

function mass(phi::SphericalPotential, r::Real)
    error("Numeric integration of density profile not implemented yet.")
end

function potential(phi::SphericalPotential, r::Real)
    M = mass(phi, r)
    return -c.G * M ./ r
end

immutable PowerLawPotential <: SphericalPotential
    #=
    r0 : scale radius
    rho0 : scale density
    alpha : power law index of density
    =#
    r0::Real
    rho0::Real
    alpha::Real
end

function mass(phi::PowerLawPotential, r::Real)
    r0 = phi.r0
    rho0 = phi.rho0
    alpha = phi.alpha
    4pi * rho0 * r0 ^ 3 / (3 + alpha) * (r / r0) .^ (3 + alpha)
end

function density(phi::PowerLawPotential, r::Real)
    r0 = phi.r0
    rho0 = phi.rho0
    alpha = phi.alpha
    rho0 * (r / r0) .^ alpha
end

end


text:
# Julia wrapper for header: /usr/include/scip/scip_prop.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeProp(scip, name, desc, priority, freq, delay, timingmask, presolpriority, presolmaxrounds, presoltiming, propcopy, propfree, propinit, propexit, propinitpre, propexitpre, propinitsol, propexitsol, proppresol, propexec, propresprop, propdata)
    ccall((:SCIPincludeProp, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Cint, Cint, UInt32, SCIP_PROPTIMING, Cint, Cint, SCIP_PRESOLTIMING, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_PROPDATA}), scip, name, desc, priority, freq, delay, timingmask, presolpriority, presolmaxrounds, presoltiming, propcopy, propfree, propinit, propexit, propinitpre, propexitpre, propinitsol, propexitsol, proppresol, propexec, propresprop, propdata)
end

function SCIPincludePropBasic(scip, propptr, name, desc, priority, freq, delay, timingmask, propexec, propdata)
    ccall((:SCIPincludePropBasic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_PROP}}, Cstring, Cstring, Cint, Cint, UInt32, SCIP_PROPTIMING, Ptr{Cvoid}, Ptr{SCIP_PROPDATA}), scip, propptr, name, desc, priority, freq, delay, timingmask, propexec, propdata)
end

function SCIPsetPropCopy(scip, prop, propcopy)
    ccall((:SCIPsetPropCopy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PROP}, Ptr{Cvoid}), scip, prop, propcopy)
end

function SCIPsetPropFree(scip, prop, propfree)
    ccall((:SCIPsetPropFree, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PROP}, Ptr{Cvoid}), scip, prop, propfree)
end

function SCIPsetPropInit(scip, prop, propinit)
    ccall((:SCIPsetPropInit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PROP}, Ptr{Cvoid}), scip, prop, propinit)
end

function SCIPsetPropExit(scip, prop, propexit)
    ccall((:SCIPsetPropExit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PROP}, Ptr{Cvoid}), scip, prop, propexit)
end

function SCIPsetPropInitsol(scip, prop, propinitsol)
    ccall((:SCIPsetPropInitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PROP}, Ptr{Cvoid}), scip, prop, propinitsol)
end

function SCIPsetPropExitsol(scip, prop, propexitsol)
    ccall((:SCIPsetPropExitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PROP}, Ptr{Cvoid}), scip, prop, propexitsol)
end

function SCIPsetPropInitpre(scip, prop, propinitpre)
    ccall((:SCIPsetPropInitpre, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PROP}, Ptr{Cvoid}), scip, prop, propinitpre)
end

function SCIPsetPropExitpre(scip, prop, propexitpre)
    ccall((:SCIPsetPropExitpre, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PROP}, Ptr{Cvoid}), scip, prop, propexitpre)
end

function SCIPsetPropPresol(scip, prop, proppresol, presolpriority, presolmaxrounds, presoltiming)
    ccall((:SCIPsetPropPresol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PROP}, Ptr{Cvoid}, Cint, Cint, SCIP_PRESOLTIMING), scip, prop, proppresol, presolpriority, presolmaxrounds, presoltiming)
end

function SCIPsetPropResprop(scip, prop, propresprop)
    ccall((:SCIPsetPropResprop, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PROP}, Ptr{Cvoid}), scip, prop, propresprop)
end

function SCIPfindProp(scip, name)
    ccall((:SCIPfindProp, libscip), Ptr{SCIP_PROP}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetProps(scip)
    ccall((:SCIPgetProps, libscip), Ptr{Ptr{SCIP_PROP}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNProps(scip)
    ccall((:SCIPgetNProps, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPsetPropPriority(scip, prop, priority)
    ccall((:SCIPsetPropPriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PROP}, Cint), scip, prop, priority)
end

function SCIPsetPropPresolPriority(scip, prop, presolpriority)
    ccall((:SCIPsetPropPresolPriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PROP}, Cint), scip, prop, presolpriority)
end


text:
# Declarations mirroring the C declarations in xsum.h

struct xsum_small_accumulator
    chunk::NTuple{67,Int64}
    inf::Int64
    nan::Int64
    adds_until_propagate::Cint
end

struct xsum_large_accumulator
    chunk::NTuple{4096,UInt64}
    count::NTuple{4096,Int16}
    chunks_used::NTuple{64,UInt64}
    used_used::UInt64
    sacc::xsum_small_accumulator
end

# The following low-level functions do little checking of their arguments, e.g. they don't check
# that the array strides == 1 or that the dot product is of vectors with equal lengths.
# They do check that length(array) ≤ typemax(Cint), however.
for kind in ("xsum_small_", "xsum_large_")
    T = Symbol(kind, "accumulator")
    @eval begin
        $(Symbol(kind,"init"))(acc) =
            ccall(($(QuoteNode(Symbol(kind,"init"))),libxsum), Cvoid, (Ref{$T},), acc)
        $(Symbol(kind,"addv"))(acc, a::StridedVector{Float64}) =
            ccall(($(QuoteNode(Symbol(kind,"addv"))),libxsum), Cvoid, (Ref{$T},Ptr{Float64},Cint), acc, a, length(a))
        $(Symbol(kind,"add_sqnorm"))(acc, a::StridedVector{Float64}) =
            ccall(($(QuoteNode(Symbol(kind,"add_sqnorm"))),libxsum), Cvoid, (Ref{$T},Ptr{Float64},Cint), acc, a, length(a))
        $(Symbol(kind,"add_dot"))(acc, a::StridedVector{Float64}, b::StridedVector{Float64}) =
            ccall(($(QuoteNode(Symbol(kind,"add_dot"))),libxsum), Cvoid, (Ref{$T},Ptr{Float64},Ptr{Float64},Cint), acc, a, b, length(a))
        $(Symbol(kind,"round"))(acc) =
            ccall(($(QuoteNode(Symbol(kind,"round"))),libxsum), Cdouble, (Ref{$T},), acc)
        $(Symbol(kind,"chunks_used"))(acc) =
            ccall(($(QuoteNode(Symbol(kind,"chunks_used"))),libxsum), Cint, (Ref{$T},), acc)
    end
end

xsum_small_add1(acc, x::Real) =
    ccall((:xsum_small_add1,libxsum), Cvoid, (Ref{xsum_small_accumulator}, Cdouble), acc, x)


text:
using CommonLH
using Test

test_dir() = joinpath(@__DIR__, "test_files");

@testset "CommonLH" begin
	include("kwargs_test.jl")
	include("check_test.jl")
	include("display_test.jl")
	include("grid_test.jl")
	include("vector_test.jl")
	include("probabilities_test.jl")
	include("arrays_test.jl");
	# include("struct_test.jl");
	include("statistics/discretize_test.jl");
	# include("check_test.jl")
	# # include("files_test.jl")
	# include("matrix_test.jl")
	# include("statistics_test.jl")
end

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


text:
using Test
using DataWrangler

@testset "outlie" begin

    x = 1:100
    y = Array{Union{Missing,Float64}}(undef,100);
    y[:] = collect(x);
    mid = vcat(10,30,60,90);
    y[mid] .= y[mid] .+ [5,10,-5,-10];

    @test outlie(x,y) == [10, 30, 60, 81, 82, 85, 86, 87, 89, 90, 93]

    ϵ = Array{Union{Missing,Float64}}(undef,100);
    ϵ[:] = rand(length(y))/100;

    yϵ = Array{Union{Missing,Float64}}(undef,100);
    yϵ[:] = y.+ϵ

    @test outlie(x,yϵ) == [10, 30, 60, 90]

end


text:
module HonlyPettableCatModule

using ..Ahorn, Maple

@mapdef Entity "HonlyHelper/PettableCat" HonlyPettableCat(x::Integer, y::Integer, catFlag::String = "CatHasBeenPet")

const placements = Ahorn.PlacementDict(
    "Pettable Cat (HonlyHelper)" => Ahorn.EntityPlacement(
        HonlyPettableCat
    )
)

sprite = "characters/HonlyHelper/pettableCat/spoons_idle00.png"

function Ahorn.selection(entity::HonlyPettableCat)
    x, y = Ahorn.position(entity)
    return Ahorn.Rectangle(x-8, y-8, 16, 16)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::HonlyPettableCat, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 4)

end

max_stars_repo_name: HollyMagala/HonlyHelper

max_issues_repo_name: HollyMagala/HonlyHelper

max_forks_repo_name: HollyMagala/HonlyHelper


text:
using VLKeggSDK
using DataFrames
using ProgressMeter

# get all the compounds in the database -
master_compound_dictionary = build_metabolite_keggid_matching_table() |> check

# so we need to make sure that *+ compounds have the * synonum in the collection -
for (key, value) in master_compound_dictionary

    if (contains(key, "+") == true)

        # build a new key -
        new_key = replace(key, "+" => "")

        # add the new key to the dict -
        master_compound_dictionary[new_key] = value
    end
end


text:
# Data is from  ``Data analysis recipes: Fitting a model to data'' by Hogg, Bovy & Lang
# available on https://arxiv.org/abs/1008.4686.

if !isdef(:versions); versions = [1] end
if !isdef(:penalty_names); penalty_names = ["huber"] end
if !isdef(:ovf_formulations); ovf_formulations = ["equilibrium"] end

using DelimitedFiles

N = length(y)

@testset "loss test: penalty = $penalty_name; version = $version; ovf_formulation = $ovf_formulation" for penalty_name in penalty_names, version in versions, ovf_formulation in ovf_formulations

m = direct_model(ReSHOP.Optimizer(;ovf_formulation=ovf_formulation))

fit_model = EMP.Model(m)

if version == 1
    # primary unknonwn
    @variable(m, c1[j=1:1])
    @variable(m, d1)

    # OVF variable
    @variable(m, loss_var1)

    # OVF arguments
    @variable(m, fit1[i=1:N]);

    @constraint(m, fiteqn[i=1:N], (y[i] - sum(c1[j]*x[i, j] for j=1:M) - d1) / sigma_y[i] == fit1[i]);

    addovf!(fit_model, loss_var1, fit1, penalty_name, all_params[penalty_name])

    @objective(m, Min, loss_var1)

    c = c1
    d = d1
    fit = fit1
    loss_var = loss_var1
elseif version == 2
    # primary unknonwn
    @variable(m, c2[j=1:1])
    @variable(m, d2)

    # OVF variable
    @variable(m, loss_var2)

    # OVF arguments
    @variable(m, fit2[i=1:N]);

    @NLconstraint(m, fiteqn[i=1:N], log(exp((y[i] - sum(c2[j]*x[i, j] for j=1:M) - d2) / sigma_y[i])) == fit2[i]);

    addovf!(fit_model, loss_var2, fit2, penalty_name, all_params[penalty_name])

    @objective(m, Min, loss_var2)

    c = c2
    d = d2
    fit = fit2
    loss_var = loss_var2
elseif version == 3
    # primary unknonwn
    @variable(m, c3[j=1:1])
    @variable(m, d3)

    # OVF variable
    @variable(m, loss_var3)

    # OVF arguments
    @variable(m, fit3[i=1:N]);

    @constraint(m, fiteqn[i=1:N], fit3[i] == (y[i] - sum(c3[j]*x[i, j] for j=1:M) - d3) / sigma_y[i]);

    addovf!(fit_model, loss_var3, fit3, penalty_name, all_params[penalty_name])

    @objective(m, Min, loss_var3)

    c = c3
    d = d3
    fit = fit3
    loss_var = loss_var3
elseif version == 4
    # primary unknonwn
    @variable(m, c4[j=1:1])
    @variable(m, d4)

    # OVF variable
    @variable(m, loss_var4)

    # OVF arguments
    @variable(m, fit4[i=1:N]);

    @NLconstraint(m, fiteqn[i=1:N], fit4[i] == log(exp((y[i] - sum(c4[j]*x[i, j] for j=1:M) - d4) / sigma_y[i])));

    addovf!(fit_model, loss_var4, fit4, penalty_name, all_params[penalty_name])

    @objective(m, Min, loss_var4)

    c = c4
    d = d4
    fit = fit4
    loss_var = loss_var4
elseif version == 5
    # primary unknonwn
    @variable(m, c5[j=1:1])
    @variable(m, d5)

    # OVF variable
    @variable(m, loss_var5)

    # OVF arguments
    @variable(m, fit5[i=1:N]);

    @constraint(m, fiteqn[i=1:N], pi*(y[i] - sum(c5[j]*x[i, j] for j=1:M) - d5) / sigma_y[i] == pi*fit5[i]);

    addovf!(fit_model, loss_var5, fit5, penalty_name, all_params[penalty_name])

    @objective(m, Min, loss_var5)

    c = c5
    d = d5
    fit = fit5
    loss_var = loss_var5
end

cmp_dir = joinpath(dirname(@__FILE__), "test_res")
unknown_ref = readdlm(joinpath(cmp_dir, "mcp_" * penalty_name * "_v1.out"))
file_ref = readdlm(joinpath(cmp_dir, "mcp_" * penalty_name * "_v1_fit.out"))

solveEMP(fit_model)

occursin("hinge", penalty_name) && continue

@test isapprox(value.(c), [unknown_ref[1]], rtol=1e-4)
@test isapprox(value(d), unknown_ref[2], rtol=1e-4)
@test isapprox(value.(fit), file_ref[1:N], rtol=1e-4)

@test isapprox(value(loss_var), unknown_ref[3], rtol=1e-4)
#@test isapprox(getobjectivevalue(m), unknown_ref[3], rtol=1e-4)

end


text:
using .JuMP

struct DiscretizedFunctionVariable{S<:AbstractRealLineDiscretizer,V,M}
    discretizer::S
    finite_param::V
    model::M
end

function Base.show(io::IO,  discvar::DiscretizedFunctionVariable)
    println(io, "DiscretizedFunctionVariable with JuMP param:")
    Base.show(io, discvar.finite_param)
end

function add_discretized_function!(model, discretizer)
    n = length(discretizer)
    tmp_vars = @variable(model, [i = 1:n])
    DiscretizedFunctionVariable(discretizer, tmp_vars, model)
end

function (f::DiscretizedFunctionVariable)(z)
    idx = findfirst(f.discretizer, z)
    f.finite_param[idx]
end

function JuMP.value(f::DiscretizedFunctionVariable)
    DiscretizedFunction(f.discretizer, JuMP.value.(f.finite_param))
end

function JuMP.value(f::DiscretizedFunctionVariable, flift)
    DiscretizedFunction(discretizer = f.discretizer,
        dictionary = Dictionary(flift.(f.discretizer), JuMP.value.(f.finite_param)))
end