text:
@testset "2-point Green's function" begin
    using LinearAlgebra
    using Atom.Green:GreenN, Gn, Gnc

    # construct model
    # U, μ, β = 1.0, 0.4, 5.0
    U, μ, β = 10.0, 5.0, 100.0
    rtol = 1e-6

    E = [0.0, -μ, -μ, U - 2μ]
    H = zeros(Float64, (4, 4))
    H[diagind(H)] = E

    # |0>=|00>=1, |↑>=|10>=2, |↓>=|01>=3, |↑↓>=|11>=4
    # the first is the forck state for ↑ spin, the second is forck state for ↓

    cpup = zeros(Float64, (4, 4))
    cpdown = zeros(Float64, (4, 4))

    cpup[2, 1], cpup[4, 3] = 1, 1
    cpdown[3, 1], cpdown[4, 2] = 1, 1
    cmup, cmdown = cpup', cpdown'

    @assert abs(tr(cpup * cmdown)) < 1e-16
    @assert abs(tr(cpdown * cmup)) < 1e-16

    m = Green.Model(β, H, [cpup, cpdown], true)

    G(τ) = (exp(μ * τ) + exp(μ * β - (U - μ) * τ)) / (1 + 2 * exp(μ * β) + exp(-(U - 2μ) * β))

    @test Green.density(m, UP) ≈ G(β)
    # @test thermalavg(m.n[UP], m.E, m.β, m.Z) ≈ G(β)

    τi, τo = 0.1, 0.8
    g2 = GreenN(m, [τi, τo], [UP, DOWN])
    @test abs(Gn(m, g2)) < 1e-16
    g2 = GreenN(m, [τi, τo], [UP, UP])
    @test Gn(m, g2) ≈ G(τo - τi)

    # G2(m, g2)
    Gn(m, g2)
    # @time G2(m, g2)
    @time Gn(m, g2)

    τi, τo = 0.8, 0.1
    g2 = GreenN(m, [τi, τo], [UP, DOWN])
    @test abs(Gn(m, g2)) < 1e-16
    g2 = GreenN(m, [τi, τo], [UP, UP])
    @test Gn(m, g2) ≈ -G(τo - τi + β)

    τ = [0.1, 0.8, 0.8 - 1e-10, 0.1 - 1e-10]
    orbital = [UP, UP, UP, UP]
    g4 = GreenN(m, τ, [UP, UP, UP, UP])
    @test Gn(m, g4) ≈ G(β) # density up or down is conserved, and n^2=m
    gc = 0.5 - 0.5^2 - G(τ[2] - τ[1])^2 # <n^2>-<n>^2=<n>-<n>^2
    @test Gnc(m, g4) ≈ gc
    # exit(0)

    Gn(m, g4)
    Gnc(m, g4)
    @time Gn(m, g4)
    @time Gnc(m, g4)

    δ = 1e-6
    # τ = [t for t in LinRange(3δ, β, 16)]
    for (ti, t) in enumerate(LinRange(3δ, β, 16))
    t4 = [0.0, 2δ, δ, t]
    g4 = GreenN(m, t4, [UP, UP, UP, UP])
        # println(Gnc(m, g4, [1,2,3,4]), " vs ", G(t))

        # <Tτ a(t)a(δ)a⁺(2δ)a⁺(0)> ≈ -<Tτ a(t)a(0)>
    @test abs(Gn(m, g4, [1,2,3,4]) + G(t)) < 10δ

        # <Tτ a(t)a(δ)a⁺(2δ)a⁺(0)>_{connected} ≈ 0.0
    @test abs(Gnc(m, g4, [1,2,3,4])) < 10δ
end

    δ = 1e-6
    for (ti, t) in enumerate(LinRange(5δ, β, 16))
    t6 = [0.0, 2δ, 4δ, δ, 3δ, t]
    g6 = GreenN(m, t6, [UP, UP, UP, UP, UP, UP])
    @test abs(Gn(m, g6, [1, 2,3,4,5,6]) + G(t)) < 10δ
    @test abs(Gnc(m, g6, [1, 2,3,4,5,6])) < 10δ
end

    t6 = [0.0, 2δ, 4δ, δ, 3δ, 0.5]
    g6 = GreenN(m, t6, [UP, UP, UP, UP, UP, UP])
    Gn(m, g6, [1,2,3,4,5,6])
    Gnc(m, g6, [1,2,3,4,5,6])
    @time Gn(m, g6, [1,2,3,4,5,6])
    @time Gnc(m, g6, [1,2,3,4,5,6])

    t12 = [0.0, 2δ, 4δ, δ, 3δ, 0.1, 0.15, 0.13, 0.11, 0.12, 0.05, 0.01]
    g12 = GreenN(m, t12, [UP for i in 1:12])
    Gn(m, g12)
    Gnc(m, g12)
    @time Gn(m, g12)
    @time Gnc(m, g12)
    println(Gnc(m, g12))

    # t14 = [0.0, 2δ, 4δ, δ, 3δ, 0.1, 0.15, 0.13, 0.11, 0.12, 0.05, 0.01, 0.06, 0.03]
    # g14 = Green(m, t14, [UP for i in 1:14])
    # Gn(m, g14)
    # Gnc(m, g14)
    # @time Gn(m, g14)
    # @time Gnc(m, g14)
    # println(Gnc(m, g12))

end


text:
function setupexperimentresults(molecule_names,
    ΩS_ppm,
    Δsys_border::Vector{Vector{T}};
    min_dist = 0.15) where T

    metabolite_regions = collect( getmetaboliteintervals(molecule_names[n],
        molecule_names, ΩS_ppm, Δsys_border;
        min_dist = min_dist) for n = 1:length(molecule_names) )

    # get connected graph. A region is defined as a connected subgraph.
    D, sts_set, fins_set, IDs_set,
        G = setupregioninfo(metabolite_regions, Δsys_border)
    # D is the list of graph nodes that form connected components.

    # information from each graph nodes.
    sts = combinevectors(combinevectors(sts_set))
    fins = combinevectors(combinevectors(fins_set))
    IDs = combinevectors(combinevectors(IDs_set))

    regions = Vector{RegionInfoType{T}}(undef, length(D))
    for i = 1:length(D)

        inds = D[i]
        ids = IDs[inds]

        s = collect( ids[k][1:2] for k = 1:length(ids) )
        spin_group_indices = sort(unique(s))

        # find start and fin of region.
        tmp = collect( fins[inds[k]] for k = 1:length(inds) )
        region_fin = maximum(tmp)

        tmp = collect( sts[inds[k]] for k = 1:length(inds) )
        region_st = minimum(tmp)

        regions[i] = RegionInfoType(ids, region_st, region_fin, spin_group_indices)
    end

    return ExperimentInfoType(regions, G, D, sts_set, fins_set, IDs_set,
        sts, fins, IDs)
end

function setupregioninfo(metabolite_regions,
    Δsys_border::Vector{Vector{T}}) where T

    # get connected graph.
    sts_set, fins_set, IDs_set, G = processintervals(metabolite_regions, Δsys_border)
    D0 = LightGraphs.connected_components(G)

    sts = combinevectors(combinevectors(sts_set))
    fins = combinevectors(combinevectors(fins_set))
    IDs = combinevectors(combinevectors(IDs_set))

    # sort D by size.
    region_sizes = collect( length(D0[i]) for i = 1:length(D0))
    inds = sortperm(region_sizes)
    region_sizes = region_sizes[inds]
    D = D0[inds]

    # sanity-check.
    @assert length(unique(combinevectors(D))) == LightGraphs.nv(G)

    # schedule jobs.
    IDs = combinevectors(combinevectors(IDs_set))

    return D, sts_set, fins_set, IDs_set, G
end

# border in ppm, adds that amount to each metabolite interval.
function processintervals(metabolite_regions::Vector{MetaboliteRegionType{T}},
    Δsys_border::Vector{Vector{T}}) where T

    # initialize regions to every interval.
    N_compounds = length(metabolite_regions)
    @assert length(Δsys_border) == N_compounds

    sts_set = Vector{Vector{Vector{T}}}(undef, N_compounds)
    fins_set = Vector{Vector{Vector{T}}}(undef, N_compounds)
    IDs_set = Vector{Vector{Vector{Tuple{Int,Int,Int}}}}(undef, N_compounds)

    for n = 1:N_compounds

        N_groups = length(Δsys_border[n])

        sts_set[n] = Vector{Vector{T}}(undef, N_groups)
        fins_set[n] = Vector{Vector{T}}(undef, N_groups)
        IDs_set[n] = Vector{Vector{Tuple{Int,Int}}}(undef, N_groups)

        for i = 1:N_groups

            sts_set[n][i] = metabolite_regions[n].interval_starts[i] .- Δsys_border[n][i]
            fins_set[n][i] = metabolite_regions[n].interval_fins[i] .+ Δsys_border[n][i]

            N_intervals = length( metabolite_regions[n].interval_starts[i] )
            #IDs_set[n][i] = (n,i,length(metabolite_regions[n].interval_fins[i]))
            IDs_set[n][i] = collect( (n,i,k) for k = 1:N_intervals )
        end
    end

    sts = combinevectors(combinevectors(sts_set))
    fins = combinevectors(combinevectors(fins_set))

    # sanity check.
    M2 = getnumelements(sts_set)
    M3 = getnumelements(IDs_set)
    @assert M2 == M3

    # construct graph.
    G = constructgraph(sts, fins)

    return sts_set, fins_set, IDs_set, G
end

# get total number of elements.
function getnumelements(X::Vector{Vector{Vector{T}}})::Int where T
    return sum( sum( length(X[n][i]) for i = 1:length(X[n]) ) for n = 1:length(X) )
end

#########

function filterfreqpositions(U_cost,
    st_us::Vector{T}, fin_us::Vector{T}) where T

    band_inds = Vector{Int}(undef, 0)
    for m = 1:length(st_us)

        st_u = st_us[m]
        fin_u = fin_us[m]
        inds = findall(xx->(st_u<xx<fin_u), U_cost)
        push!(band_inds, inds...)
    end

    return unique(band_inds)
end

function getcostinds(exp_info::ExperimentInfoType{T}, P_cost0) where T

    band_inds = Vector{Int}(undef, 0)
    band_inds_set = Vector{Vector{Int}}(undef, 0)

    for i = 1:length(exp_info.regions)

        tmp = filterfreqpositions(P_cost0, [exp_info.regions[i].st;], [exp_info.regions[i].fin;])
        push!(band_inds, tmp...)
        push!(band_inds_set, tmp)
    end

    unique!(band_inds)
    return band_inds, band_inds_set
end

ext: jl


text:
"""
    Function SLEPc lib location.
"""
function get_slepc_location()
    SLEPC_DIR   = haskey(ENV,"SLEPC_DIR") ? ENV["SLEPC_DIR"] : "/usr/lib/slepc"
    PETSC_ARCH  = haskey(ENV,"PETSC_ARCH") ? ENV["PETSC_ARCH"] : ""
    SLEPC_LIB = ""

    # Check SLEPC_DIR exists
    if isdir(SLEPC_DIR)

        # Define default paths
        SLEPC_LIB_DIR = joinpath(SLEPC_DIR,PETSC_ARCH,"lib")

        # Check SLEPC_LIB (.../libslepc.so or .../libslepc_real.so file) exists
        if isfile(joinpath(SLEPC_LIB_DIR,"libslepc.so"))
            SLEPC_LIB = joinpath(SLEPC_LIB_DIR,"libslepc.so")
        elseif isfile(joinpath(SLEPC_LIB_DIR,"libslepc_real.so"))
            SLEPC_LIB = joinpath(SLEPC_LIB_DIR,"libslepc_real.so")
        end
    end

    # PETSc lib not found
    if(length(SLEPC_LIB) == 0)
        # Workaround for automerging on RegistryCI or doc deployment
        if(haskey(ENV,"JULIA_REGISTRYCI_AUTOMERGE") || haskey(ENV, "DOC_DEPLOYMENT"))
            @warn "Setting fictive SLEPc path because of RegistryCI or Doc deployment"
            SLEPC_LIB = "JULIA_REGISTRYCI_AUTOMERGE"
        else
            @show ENV
            throw(ErrorException("SLEPc shared library (libslepc.so) not found. Please check that SLEPC_DIR and PETSC_ARCH env. variables are set."))
        end
    end

    return SLEPC_LIB
end

# Absolute path to libslepc.so
const libslepc = get_slepc_location()


text:
function ReadCSV(name; type=Float32)
    pf = open(name)
    line = readline(pf)
    sep = findfirst(" ", line)[1]
    len = parse(Int, line[1:(sep-1)])
    dim = parse(Int, line[(sep+1):end])
    Points = zeros(type, len, dim)

    for i = 1:len
        line = readline(pf)
        for j = 1:(dim-1)
            sep = findfirst(" ", line)[1]
            val = parse(type, line[1:(sep-1)])
            line = line[(sep+1):end]
            Points[i, j] = val
        end
        val = parse(type, line)
        Points[i, dim] = val
    end

    return  Points

end

function AssignPointsToBlock(Points::Array{Float32,2},
    numOfGrids::Int64, dimensions::Int64)
    grids = [numOfGrids for i = 1:dimensions];
    gridMap1D = grids.^((0:(dimensions-1)));
    gridMap1Doffset = [0; -1*(2:dimensions).^0];
    BlockOfPoint = dropdims(
        mapslices(
        x -> Int32(sum( ((1 .+ div.(min.(x, 0.99999999), 1/numOfGrids)) + gridMap1Doffset ) .*gridMap1D )),
            Points, dims=2),
        dims=2);
    return BlockOfPoint
end

function CountPointsPerBlock(Points::Array{Float32,2},
    numOfGrids::Int64, dimensions::Int64)
    grids = [numOfGrids for i = 1:dimensions];
    gridMap1D = grids.^((0:(dimensions-1)));
    gridMap1Doffset = [0; -1*(2:dimensions).^0];
    PointsPerBlock = zeros(Int32, numOfGrids^dimensions)
    num = size(Points)[1]
    for i = 1:num
        x = Points[i,:]
        x = min.(x, 0.99999999)
        gid = Int64(sum(( (div.(x,1/numOfGrids) .+ 1) + gridMap1Doffset ) .*gridMap1D ))
        PointsPerBlock[gid...] += 1
    end
    IntegralPointsPerBlock = similar(PointsPerBlock)
    IntegralPointsPerBlock[1] = 0#PointsPerBlock[0]
    for i = 2:numOfGrids^dimensions
        IntegralPointsPerBlock[i] = PointsPerBlock[i-1] + IntegralPointsPerBlock[i-1]
    end
    return PointsPerBlock, IntegralPointsPerBlock
end

function ReorderPointsByBlock(Points::Array{Float32,2},
    BlockOfPoint::Array{Int32,1})
    BlockPermutationP = sortperm(BlockOfPoint);
    OrderedPoints = Points[BlockPermutationP,:];
    return OrderedPoints
end


text:
"""
    BTGFuncs
BTGFuncs contains definitions for various functions and classes of functions used by
the BTG model
"""
module BTGFuncs

using Distributions, Distances

export
    Parametric,
    Correlation,
    IsotropicCorrelation,
    SquaredExponential,
    Transform,
    BoxCox,
    Covariate,
    Identity,
    kernelmatrix,
    prime

"""
    Parametric

Abstract supertype for all parametric functions
"""
abstract type Parametric end

"""
    Correlation

Abstract supertype for all correlation functions
"""
abstract type Correlation <: Parametric end

"""
    IsotropicCorrelation <: Correlation

Abstract supertype for all isotropic correlation functions
"""
abstract type IsotropicCorrelation <: Correlation end

"""
    kernelmatrix(k::IsotropicCorrelation, θ, A, B)

Construct a correlation matrix by applying `k` with
parameters `θ` pairwise to the rows of `A` and `B`
"""
function kernelmatrix(k::IsotropicCorrelation, θ, A, B)
    return k.(θ, pairwise(Euclidean(), A, B, dims=1))
end

"""
    SquaredExponential(γ::Uniform)

The squared exponential kernel (aka the Gaussian kernel or the radial basis
function kernel) with a length scale
parameter `γ` given by a continuous uniform distribution.

```math
k_\\theta(\\tau) = \\exp\\Big\\lbrace -\\Big(\\frac{\\tau}{\\theta}\\Big)^2\\Big\\rbrace
```

```julia
SquaredExponential(Uniform()) # γ is uniform over [0, 1]
SquaredExponential(Uniform(a, b)) # γ is uniform over [a, b]

sampleparam(k) # A random sample of the parameters of k
k(θ, τ) # The kernel evaluated on the distance τ with parameters θ
```

External links
* [RBF Kernel on Wikipedia](https://en.wikipedia.org/wiki/Radial_basis_function_kernel)

"""
struct SquaredExponential <: IsotropicCorrelation
    γ::Union{Uniform, Exponential, Gamma}
end

function (k::SquaredExponential)(θ, τ)
    return exp(-(τ / θ[1])^2)
end

"""
    Transform

Abstract supertype for all transformation functions
"""
abstract type Transform <: Parametric end

"""
    BoxCox(λ::Uniform)

The Box-Cox power transformation function with parameter λ given
by a continuous uniform distribution.

```math
g_\\lambda(x) = \\begin{cases}
    \\frac{x^\\lambda - 1}{\\lambda} & \\lambda \\neq 0\\\\
    \\ln(x) & \\lambda = 0
\\end{cases}
```

```julia
BoxCox(Uniform()) # λ is uniform over [0, 1]
BoxCox(Uniform(a, b)) # λ is uniform over [a, b]

sampleparam(g) # A random sample of the parameters of g
g(λ, x) # The transformation evaluated on x with parameters y
prime(g, λ, x) # The derivative of the transformation evaluated on x with parameters y
```

External links
* [Power Transform on Wikipedia](https://en.wikipedia.org/wiki/Power_transform)

"""
struct BoxCox <: Transform
    λ::Union{Uniform, Normal}
end

function (g::BoxCox)(λ, x)
    return λ[1] == 0 ? log(x) : expm1(log(x) * λ[1]) / λ[1]
end

function prime(g::BoxCox, λ, x)
    return x^(λ[1] - 1)
end

"""
    Covariate

Abstract supertype for all covariate functions
"""
abstract type Covariate end

"""
    Identity

A covariate function which returns its input as-is
"""
struct Identity <: Covariate end

function (f::Identity)(x)
    return x
end

end # module

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module Separators

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, template, register_normal_element

export separator

register_normal_element("q__separator", context = @__MODULE__)

function separator(args...; wrap::Function = StippleUI.DEFAULT_WRAPPER, kwargs...)
  wrap() do
    q__separator(args...; kwargs...)
  end
end

end


text:
using Gadfly
import Cairo, Fontconfig

"""
  clippedStats(reader)
  clippedStats(reader, chrom::String, chromrange::UnitRange)

Returns the number of reads that were clipped from the BAM file, or
from a region as readsClipped, readsTotal, and a DataFrame containing
all of the reads mapped to the region with the read type (soft/hard/both/not)
and mean quality of the read
"""
function clippedStats(reader::XAM.BAM.Reader, chrom::String, chromrange::UnitRange)
  type = String[]
  qual = Float64[]
  readsClipped = 0
  readsTotal = 0
  for record in XAM.BAM.GenomicFeatures.eachoverlap(reader, chrom, chromrange)
    readsTotal = readsTotal + 1
    clippedStatsHelper!(record,qual,type,readsClipped)
  end
  return readsClipped, readsTotal, DataFrames.DataFrame(type=type,qual=qual)
end

function clippedStats(reader::XAM.BAM.Reader)
  type = String[]
  qual = Float64[]
  readsClipped = 0
  readsTotal = 0
  for record in reader
    readsTotal = readsTotal + 1
    clippedStatsHelper!(record,qual,type,readsClipped)
  end
  return readsClipped, readsTotal, DataFrames.DataFrame(type=type,qual=qual)
end

function clippedStatsHelper!(record::XAM.BAM.Record, qual::Vector{Float64},type::Vector{String},
  readsClipped::Int8)
  cigar = XAM.BAM.cigar_rle(record)
  push!(qual,Statistics.mean(XAM.BAM.quality(record)))
  flag = "not"
  if findfirst(x->x==convert(BioAlignments.Operation,'S'),cigar[1]) != nothing
    flag = "soft"
  end
  if findfirst(x->x==convert(BioAlignments.Operation,'H'),cigar[1])!= nothing
    if flag == "soft"
      push!(type,"both")
    else
      push!(type,"hard")
    end
    readsClipped = readsClipped + 1
  else
    if flag == "soft"
      push!(type,"soft")
      readsClipped = readsClipped + 1
    else
      push!(type,"not")
    end
  end
  return type,qual,readsClipped
end

"""
  clippedReadQual(df::DataFrame, outfile::String="clippedReadQual.pdf")

Plots a histogram of mean read quality split by whether the read was soft-clipped,
hard-clipped, both, or not clipped and saves it to outfile.
"""
function clippedReadQual(df::DataFrames.DataFrame, outfile::String="clippedReadQual.pdf")
  p = plot(df,x="qual", color="type", Geom.histogram,
    Guide.title("Quality score distribution for clipped and not in region"),
    Guide.xlabel("Mean quality"),
    Guide.ylabel("Number of reads"))
  p |> PDF(outfile)
end

#using BenchmarkTools
#@benchmark
#reader = open(BAM.Reader, "/Users/aguang/CORE/scratch/Bliss_71991_101.dedup.rg.srtd.bam",
#  index="/Users/aguang/CORE/scratch/Bliss_71991_101.dedup.rg.srtd.bam.bai")
#df = clippedStats("Ca22chr4A_C_albicans_SC5314",726523:727338,reader)
#print(df)
#clippedReadQual(df)
#close(reader)

max_stars_repo_name: compbiocore/BioJulia-cookbook

max_issues_repo_name: compbiocore/BioJulia-cookbook

max_forks_repo_name: compbiocore/BioJulia-cookbook


text:
function compile_romeo(path;
        app_name="romeo",
        filter_stdlibs=true,
        precompile_execution_file=abspath(joinpath(@__DIR__, "..", "test", "romeo_test.jl")),
        clean=true,
        kw...)
    romeopath = pathof("RomeoApp")
    if !isdir(romeopath)
        download_pkg("RomeoApp")
    end
    create_app(romeopath, path; app_name, filter_stdlibs, precompile_execution_file, kw...)
    test_romeo(path, app_name)
    if clean
        clean_app(path, app_name) # remove unneccesary artifacts dir (600MB)
    end
    copy_matlab(path)
    printstyled("Success! Romeo compiled and tested!\n"; color=:green)
    @warn("Relocatability has to be tested manually!")
end

function test_romeo(path, app_name)
    file = tempname()
    phasefile = abspath(joinpath(@__DIR__, "..", "test", "data", "small", "Phase.nii"))
    args = [phasefile, "-o", file, "-t", "1:3", "-k", "nomask"]
    name = app_name * (Sys.iswindows() ? ".exe" : "")
    romeofile = joinpath(path, "bin", name)
    @assert isfile(romeofile)
    cmd = `$romeofile $args`
    @assert success(run(cmd))
end

function update_romeo()
    try
        rm(pathof("RomeoApp"); force=true, recursive=true)
    catch
        @warn "Couldn't remove the old RomeoApp folder! ($(pathof("RomeoApp"))) Maybe it is opened in another App"
    end
    download_pkg("RomeoApp")
end


text:
# Geometry which displays arbitrary shapes at given (x, y) positions.
# note: vertices is a list of shapes
struct ShapeGeometry <: Gadfly.GeometryElement
    vertices::AbstractVector #{Tuple{Float64,Float64}}
    tag::Symbol

    function ShapeGeometry(shape; tag::Symbol=Gadfly.Geom.empty_tag)
        new(shape, tag)
    end
end

# TODO: add for PR
# const shape = ShapeGeometry

function Gadfly.element_aesthetics(::ShapeGeometry)
    [:x, :y, :size, :color]
end

# Generate a form for a shape geometry.
#
# Args:
#   geom: shape geometry.
#   theme: the plot's theme.
#   aes: aesthetics.
#
# Returns:
#   A compose Form.
#
function Gadfly.render(geom::ShapeGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)

    # TODO: add for PR
    # Gadfly.assert_aesthetics_defined("Geom.shape", aes, :x, :y)
    # Gadfly.assert_aesthetics_equal_length("Geom.shape", aes,
    #                                       element_aesthetics(geom)...)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = Gadfly.DataFrames.PooledDataArray(RGBA{Float32}[theme.default_color])
    default_aes.size = Compose.Measure[theme.default_point_size]
    aes = Gadfly.inherit(aes, default_aes)

    lw_hover_scale = 10
    lw_ratio = theme.line_width / aes.size[1]

    aes_x, aes_y = Gadfly.concretize(aes.x, aes.y)

    ctx = Compose.compose!(
        Compose.context(),
        make_polygon(geom, aes.x, aes.y, aes.size),
        Compose.fill(aes.color),
        Compose.linewidth(theme.highlight_width))

    if aes.color_key_continuous != nothing && aes.color_key_continuous
        Compose.compose!(ctx,
            Compose.stroke(map(theme.continuous_highlight_color, aes.color)))
    else
        Compose.compose!(ctx,
            Compose.stroke(map(theme.discrete_highlight_color, aes.color)),
            Compose.svgclass([Gadfly.svg_color_class_from_label(Gadfly.escape_id(aes.color_label([c])[1]))
                      for c in aes.color]))
    end

    return Compose.compose!(Compose.context(order=4), Compose.svgclass("geometry"), ctx)
end

function gadflyshape(sv::Shape)
  ShapeGeometry(Any[vertices(sv)])
end

function gadflyshape(sv::AVec{Shape})
    ShapeGeometry(Any[vertices(s) for s in sv])
end

# create a Compose context given a ShapeGeometry and the xs/ys/sizes
function make_polygon(geom::ShapeGeometry, xs::AbstractArray, ys::AbstractArray, rs::AbstractArray)
  n = max(length(xs), length(ys), length(rs))
  T = Tuple{Compose.Measure, Compose.Measure}
  polys = Array(Vector{T}, n)
  for i in 1:n
    x = Compose.x_measure(xs[mod1(i, length(xs))])
    y = Compose.y_measure(ys[mod1(i, length(ys))])
    r = rs[mod1(i, length(rs))]
    polys[i] = T[(x + r * sx, y + r * sy) for (sx,sy) in _cycle(geom.vertices, i)]
  end
  Gadfly.polygon(polys, geom.tag)
end

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


text:
using GenomicIndexes
using Base.Test

@testset "Tabix" begin
    Tabix = GenomicIndexes.Tabix
    tabixfile = joinpath(dirname(@__FILE__), "TAIR10_GFF3_genes.sorted.gff.bgz.tbi")

    # Read from a file.
    tbi = Tabix(tabixfile)
    @test tbi isa Tabix
    @test GenomicIndexes.seqnames(tbi) == ["Chr1", "Chr2", "Chr3", "Chr4", "Chr5", "ChrC", "ChrM"]

    # Read from a stream.
    tbi = Tabix(open(tabixfile))
    @test tbi isa Tabix
    @test GenomicIndexes.seqnames(tbi) == ["Chr1", "Chr2", "Chr3", "Chr4", "Chr5", "ChrC", "ChrM"]

    # Test overlapchunks.
    @test GenomicIndexes.overlapchunks(tbi, "Chr1", 240057:242608) == [
        (0x0000000000000000, 0x00000006dfe7e426),
        (0x0000000a9781af62, 0x0000000a9781b540)]
    @test GenomicIndexes.overlapchunks(tbi, "Chr1", 30415837:30418323) == [
        (0x0000000000000000, 0x000000000000003a),
        (0x0000000a9781af62, 0x0000000a9781b540),
        (0x0000001052ebb299, 0x000000147216192d)]
    @test GenomicIndexes.overlapchunks(tbi, "Chr4", 677054:679237) == [
        (0x0000002fa1fb670b, 0x0000003376fb1fbc)]
    @test GenomicIndexes.overlapchunks(tbi, "ChrM", 292973:293431) == [
        (0x0000004d38e83e1c, 0x0000004d38e8fa5c)]
    @test_throws ArgumentError GenomicIndexes.overlapchunks(tbi, "ChrZ", 1000:20000)
end

@testset "BAI" begin
    BAI = GenomicIndexes.BAI
    baifile = joinpath(dirname(@__FILE__), "ex1.bam.bai")

    bai = BAI(baifile)
    @test bai isa BAI

    bai = BAI(open(baifile))
    @test bai isa BAI
end

@testset "CSI" begin
    CSI = GenomicIndexes.CSI
    csifile = joinpath(dirname(@__FILE__), "TAIR10_GFF3_genes.sorted.gff.bgz.csi")

    # Load from a file.
    csi = CSI(csifile)
    @test csi isa CSI

    # Load from a stream.
    csi = CSI(open(csifile))
    @test csi isa CSI

    # Test overlapchunks.
    @test GenomicIndexes.overlapchunks(csi, 1, 240057:242608) == [
        (0x0000000000000000, 0x00000006dfe7e426),
        (0x0000000a9781af62, 0x0000000a9781b540)]
    @test GenomicIndexes.overlapchunks(csi, 1, 30415837:30418323) == [
        (0x0000000000000000, 0x000000000000003a),
        (0x0000000a9781af62, 0x0000000a9781b540),
        (0x0000001052ebb299, 0x000000147216192d)]
    @test GenomicIndexes.overlapchunks(csi, 4, 677054:679237) == [
        (0x0000002fa1fb670b, 0x0000003376fb1fbc)]
    @test GenomicIndexes.overlapchunks(csi, 7, 292973:293431) == [
        (0x0000004d38e83e1c, 0x0000004d38e8fa5c)]
    @test_throws ArgumentError GenomicIndexes.overlapchunks(csi, 10, 1000:20000)
end


text:
@testset "Mesh Basics" begin
    @testset "DG Mesh basics" begin
        @test typeof(dgmesh)<:DGMesh
        @test n_intervals(dgmesh)==length(nodes)-1
        @test Δ(dgmesh) == nodes[2:end]-nodes[1:end-1]
        @test Δ(dgmesh,1) == nodes[2]-nodes[1]
        @test n_bases_per_phase(dgmesh) == (length(nodes)-1)*nbases
        @test n_bases_per_cell(dgmesh) == 3
        @test cell_nodes(dgmesh)≈
            [nodes[1:end-1]';(nodes[1:end-1]'+nodes[2:end]')/2;nodes[2:end]'] atol=1e-5
        @test DiscretisedFluidQueues.basis(dgmesh) == "lagrange"
        # @test local_dg_operators(dgmesh) == ???
        # test QBDidx
    end

    @testset "FV Mesh basics" begin
        @test typeof(fvmesh)<:FVMesh
        @test n_intervals(fvmesh)==length(nodes)-1
        @test Δ(fvmesh) == nodes[2:end]-nodes[1:end-1]
        @test Δ(fvmesh,1) == nodes[2]-nodes[1]
        @test n_bases_per_phase(fvmesh) == (length(nodes)-1)
        @test n_bases_per_cell(fvmesh) == 1
        @test DiscretisedFluidQueues._order(fvmesh) == fv_order
        @test cell_nodes(fvmesh)≈Array(((fvmesh.nodes[1:end-1] + fvmesh.nodes[2:end]) / 2 )') atol=1e-5
        @test DiscretisedFluidQueues.basis(fvmesh) == ""
    end

    @testset "FRAP Mesh basics" begin
        @test typeof(frapmesh)<:FRAPMesh
        @test n_intervals(frapmesh)==length(nodes)-1
        @test Δ(frapmesh) == nodes[2:end]-nodes[1:end-1]
        @test Δ(frapmesh,1) == nodes[2]-nodes[1]
        @test n_bases_per_phase(frapmesh) == (length(nodes)-1)*order
        @test n_bases_per_cell(frapmesh) == order
        @test cell_nodes(frapmesh)≈Array(((frapmesh.nodes[1:end-1] + frapmesh.nodes[2:end]) / 2 )') atol=1e-5
        @test DiscretisedFluidQueues.basis(frapmesh) == ""
    end
end

max_stars_repo_name: angus-lewis/StochasticFluidQueues

max_issues_repo_name: angus-lewis/StochasticFluidQueues

max_forks_repo_name: angus-lewis/StochasticFluidQueues


text:
module TestFoldLC
include("preamble.jl")

using Bifurcations: resolved_points, LimitCycleProblem
using Bifurcations.Codim1LimitCycle
using Bifurcations.Codim2LimitCycle: FoldLimitCycleProblem
using Bifurcations.Examples.Bautin

prob = Bautin.make_prob(Bautin.BautinParam(β₂ = 1.0))
solver0 = init(prob)
solve!(solver0)

hopf_points = resolved_points(solver0)
@test length(hopf_points) == 1
@test hopf_points[1].point_type === Codim1.PointTypes.hopf

prob_lc = LimitCycleProblem(
    hopf_points[1], solver0;
    num_mesh = 20,
    degree = 3,
)
solver_lc = init(
    prob_lc;
    start_from_nearest_root = true,
    bidirectional_first_sweep = false,
    max_branches = 0,
    nominal_angle_rad = 0.01,
)
@info "Solving LimitCycleProblem..."
@time solve!(solver_lc)

@info "Testing LimitCycleProblem solution..."
@test length(solver_lc.sol.sweeps[1].super.simple_bifurcation) == 0

@time flc_points = resolved_points(solver_lc)
@test length(flc_points) == 1
@test flc_points[1].point_type === Codim1LimitCycle.PointTypes.saddle_node

@info "Preparing FoldLimitCycleProblem..."
prob_flc = BifurcationProblem(
    flc_points[1],
    solver_lc,
    (@lens _.β₂),
    (-2.0, 2.0);
)
solver_flc = init(
    prob_flc;
    start_from_nearest_root = true,
    bidirectional_first_sweep = false,
    max_branches = 0,
    nominal_angle_rad = 0.01,
    verbose = true,
)
@info "Solving FoldLimitCycleProblem..."
@time solve!(solver_flc)

@info "Done!"

end  # module


text:
function update_automerge_comment!(data::GitHubAutoMergeData, body::AbstractString)::Nothing
    if data.read_only
        @info "`read_only` mode; skipping updating automerge comment."
        return nothing
    end
    api = data.api
    repo = data.registry
    pr = data.pr
    auth = data.auth
    whoami = data.whoami

    my_comments = my_retry(
        () -> get_all_my_pull_request_comments(api, repo, pr; auth=auth, whoami=whoami)
    )
    num_comments = length(my_comments)
    _body = string(
        strip(
            string(
                body, "\n", "<!---\n", "this_is_the_single_automerge_comment\n", "--->\n"
            ),
        ),
    )
    if num_comments > 1
        for i in 2:num_comments
            comment_to_delete = my_comments[i]
            try
                my_retry(() -> delete_comment!(api, repo, pr, comment_to_delete; auth=auth))
            catch ex
                @error("Ignoring error: ", exception = (ex, catch_backtrace()))
            end
        end
        comment_to_update = my_comments[1]
        if strip(comment_to_update.body) != _body
            my_retry(
                () -> edit_comment!(api, repo, pr, comment_to_update, _body; auth=auth)
            )
        end
    elseif num_comments == 1
        comment_to_update = my_comments[1]
        if strip(comment_to_update.body) != _body
            my_retry(
                () -> edit_comment!(api, repo, pr, comment_to_update, _body; auth=auth)
            )
        end
    else
        always_assert(num_comments < 1)
        my_retry(() -> post_comment!(api, repo, pr, _body; auth=auth))
    end
    return nothing
end


text:
using CuNESSie: _project2d, _regularyukawapot_double, _regularyukawapot_single

@testset "_project2d" begin
    function  _project2d_kernel!(
        dst     ::CuDeviceVector{T},
        elements::CuTriangleVector{T},
        refs    ::CuDeviceVector{T},
        eidx    ::Int,
        numref  ::Int
    ) where T
        for ridx in 1:numref
            pr = _project2d(elements[eidx], refs[(ridx-1) * 2 + 1], refs[(ridx-1) * 2 + 2])
            dst[(ridx-1) * 3 + 1] = pr[1]
            dst[(ridx-1) * 3 + 2] = pr[2]
            dst[(ridx-1) * 3 + 3] = pr[3]
        end
        nothing
    end

    for T in [Float32, Float64]
        elements = TriangleVector([
            Triangle(T[0, 0, 0], T[0, 1, 0], T[0, 0, 1]),
            Triangle(T[0, 0, 0], T[0, 2, 0], T[0, 0, 2]),
            Triangle(T[1, 0, 0], T[1, 2, 0], T[1, 0, 2])
        ])
        refs = CuArray(T[0, 0, 0, .5, .5, 0, .5, .5, 0, 1, 1, 0])
        dst = CuArray{T}(undef, 3 * 6)

        @cuda _project2d_kernel!(dst, elements, refs, 1, 6)
        res = Array(dst)
        @test res[1:3] ≈ T[0, 0, 0]
        @test res[4:6] ≈ T[0, 0, .5]
        @test res[7:9] ≈ T[0, .5, 0]
        @test res[10:12] ≈ T[0, .5, .5]
        @test res[13:15] ≈ T[0, 0, 1]
        @test res[16:18] ≈ T[0, 1, 0]

        @cuda _project2d_kernel!(dst, elements, refs, 2, 6)
        res = Array(dst)
        @test res[1:3] ≈ T[0, 0, 0]
        @test res[4:6] ≈ T[0, 0, 1]
        @test res[7:9] ≈ T[0, 1, 0]
        @test res[10:12] ≈ T[0, 1, 1]
        @test res[13:15] ≈ T[0, 0, 2]
        @test res[16:18] ≈ T[0, 2, 0]

        @cuda _project2d_kernel!(dst, elements, refs, 3, 6)
        res = Array(dst)
        @test res[1:3] ≈ T[1, 0, 0]
        @test res[4:6] ≈ T[1, 0, 1]
        @test res[7:9] ≈ T[1, 1, 0]
        @test res[10:12] ≈ T[1, 1, 1]
        @test res[13:15] ≈ T[1, 0, 2]
        @test res[16:18] ≈ T[1, 2, 0]
    end
end

@testset "_regularyukawapot_single" begin
    function _regularyukawapot_single_kernel!(
        dst  ::CuDeviceVector{T},
        vals ::CuDeviceVector{T},
        nvals::Int,
        yuk  ::T
    ) where T
        for i in 1:nvals
            dst[i] = _regularyukawapot_single(vals[i], yuk)
        end
        nothing
    end

    for T in [Float32, Float64]
        x    = ones(T, 3)
        o    = zeros(T, 3)
        yuk  = T(7)
        vals = CuArray(T[0, √12, √3, √(3 * .001^2), √(3 * .0001^2), √(3 * .00001^2)])
        dst  = CuArray{T}(undef, 6)

        @cuda _regularyukawapot_single_kernel!(dst, vals, 6, yuk)
        res = Array(dst)
        @test res[1] ≈ Radon.regularyukawapot(x, x, yuk)
        @test res[2] ≈ Radon.regularyukawapot(x, .-x, yuk)
        @test res[3] ≈ Radon.regularyukawapot(x, o, yuk)
        @test res[4] ≈ Radon.regularyukawapot(T(.001) .* x, o, yuk)
        @test res[5] ≈ Radon.regularyukawapot(T(.0001) .* x, o, yuk)
        @test res[6] ≈ Radon.regularyukawapot(T(.00001) .* x, o, yuk)
    end
end

@testset "_regularyukawapot_double" begin
    function _regularyukawapot_double_kernel!(
        dst  ::CuDeviceVector{T},
        vals ::CuDeviceVector{T},
        coss ::CuDeviceVector{T},
        nvals::Int,
        yuk  ::T
    ) where T
        for i in 1:nvals
            dst[i] = _regularyukawapot_double(vals[i], coss[i], yuk)
        end
        nothing
    end

    for T in [Float32, Float64]
        x    = ones(T, 3)
        o    = zeros(T, 3)
        n    = T[1, 0, 0]
        yuk  = T(7)
        vals = CuArray(T[0, √12, √3, √(3 * .001^2), √(3 * .0001^2), √(3 * .00001^2)])
        coss = CuArray(T[0, 2, 1, .001, .0001, .00001])
        dst  = CuArray{T}(undef, 6)

        @cuda _regularyukawapot_double_kernel!(dst, vals, coss, 6, yuk)
        res = Array(dst)
        @test res[1] ≈ Radon.∂ₙregularyukawapot(x, x, yuk, n)
        @test res[2] ≈ Radon.∂ₙregularyukawapot(x, .-x, yuk, n)
        @test res[3] ≈ Radon.∂ₙregularyukawapot(x, o, yuk, n)
        @test res[4] ≈ Radon.∂ₙregularyukawapot(T(.001) .* x, o, yuk, n)
        @test res[5] ≈ Radon.∂ₙregularyukawapot(T(.0001) .* x, o, yuk, n)
        @test res[6] ≈ Radon.∂ₙregularyukawapot(T(.00001) .* x, o, yuk, n)
    end
end

@testset "regularyukawapot_single" begin
    function _regularyukawapot_single_kernel!(
        dst      ::CuDeviceVector{T},
        Ξ        ::CuPositionVector{T},
        elements ::CuTriangleVector{T},
        eidx     ::Int,
        yuk      ::T
    ) where T
        elem = elements[eidx]
        for ξidx in 1:length(Ξ)
            dst[ξidx] = regularyukawapot_single(Ξ[ξidx], elem, yuk)
        end
        nothing
    end

    for T in [Float32, Float64]
        yuk = T(7)
        Ξ = PositionVector([
            T[0, 0, 0],  T[0, 1, 0],   T[0, 0, 1],
            T[0, -1, 0], T[0, -1, -1],
            T[1, 0, 0],  T[1, 1, 0],   T[1, 0, 1],
            T[1, -1, 0], T[1, -1, -1]
        ])
        dst  = CuArray{T}(undef, length(Ξ))

        for elem in [
            Triangle(T[0, 0, 0], T[0, 1, 0], T[0, 0, 1]),
            Triangle(T[0, 0, 0], T[0, 2, 0], T[0, 0, 2]),
            Triangle(T[1, 0, 0], T[1, 2, 0], T[1, 0, 2])
        ]
            @cuda _regularyukawapot_single_kernel!(
                dst, Ξ, TriangleVector([elem]), 1, yuk)
            res = Array(dst)
            @test res[1] ≈ Radon.regularyukawacoll(SingleLayer, T[0, 0, 0], elem, yuk)
            @test res[2] ≈ Radon.regularyukawacoll(SingleLayer, T[0, 1, 0], elem, yuk)
            @test res[3] ≈ Radon.regularyukawacoll(SingleLayer, T[0, 0, 1], elem, yuk)
            @test res[4] ≈ Radon.regularyukawacoll(SingleLayer, T[0, -1, 0], elem, yuk)
            @test res[5] ≈ Radon.regularyukawacoll(SingleLayer, T[0, -1, -1], elem, yuk)
            @test res[6] ≈ Radon.regularyukawacoll(SingleLayer, T[1, 0, 0], elem, yuk)
            @test res[7] ≈ Radon.regularyukawacoll(SingleLayer, T[1, 1, 0], elem, yuk)
            @test res[8] ≈ Radon.regularyukawacoll(SingleLayer, T[1, 0, 1], elem, yuk)
            @test res[9] ≈ Radon.regularyukawacoll(SingleLayer, T[1, -1, 0], elem, yuk)
            @test res[10] ≈ Radon.regularyukawacoll(SingleLayer, T[1, -1, -1], elem, yuk)
        end
    end
end

@testset "regularyukawapot_double" begin
    function _regularyukawapot_double_kernel!(
        dst      ::CuDeviceVector{T},
        Ξ        ::CuPositionVector{T},
        elements ::CuTriangleVector{T},
        eidx     ::Int,
        yuk      ::T
    ) where T
        elem = elements[eidx]
        for ξidx in 1:length(Ξ)
            dst[ξidx] = regularyukawapot_double(Ξ[ξidx], elem, yuk)
        end
        nothing
    end

    for T in [Float32, Float64]
        yuk = T(7)
        Ξ = PositionVector([
            T[0, 0, 0],  T[0, 1, 0],   T[0, 0, 1],
            T[0, -1, 0], T[0, -1, -1],
            T[1, 0, 0],  T[1, 1, 0],   T[1, 0, 1],
            T[1, -1, 0], T[1, -1, -1]
        ])
        dst  = CuArray{T}(undef, length(Ξ))

        for elem in [
            Triangle(T[0, 0, 0], T[0, 1, 0], T[0, 0, 1]),
            Triangle(T[0, 0, 0], T[0, 2, 0], T[0, 0, 2]),
            Triangle(T[1, 0, 0], T[1, 2, 0], T[1, 0, 2])
        ]
            @cuda _regularyukawapot_double_kernel!(
                dst, Ξ, TriangleVector([elem]), 1, yuk)
            res = Array(dst)
            @test res[1] ≈ Radon.regularyukawacoll(DoubleLayer, T[0, 0, 0], elem, yuk)
            @test res[2] ≈ Radon.regularyukawacoll(DoubleLayer, T[0, 1, 0], elem, yuk)
            @test res[3] ≈ Radon.regularyukawacoll(DoubleLayer, T[0, 0, 1], elem, yuk)
            @test res[4] ≈ Radon.regularyukawacoll(DoubleLayer, T[0, -1, 0], elem, yuk)
            @test res[5] ≈ Radon.regularyukawacoll(DoubleLayer, T[0, -1, -1], elem, yuk)
            @test res[6] ≈ Radon.regularyukawacoll(DoubleLayer, T[1, 0, 0], elem, yuk)
            @test res[7] ≈ Radon.regularyukawacoll(DoubleLayer, T[1, 1, 0], elem, yuk)
            @test res[8] ≈ Radon.regularyukawacoll(DoubleLayer, T[1, 0, 1], elem, yuk)
            @test res[9] ≈ Radon.regularyukawacoll(DoubleLayer, T[1, -1, 0], elem, yuk)
            @test res[10] ≈ Radon.regularyukawacoll(DoubleLayer, T[1, -1, -1], elem, yuk)
        end
    end
end


text:
function is_horizontal_or_vertical(startx, starty, endx, endy)
    return (startx==endx) || (starty == endy)
end

function test_is_horizontal_or_vertical()
    @assert is_horizontal_or_vertical(3,4,3,7)
    @assert is_horizontal_or_vertical(2,4,3,4)
    @assert !is_horizontal_or_vertical(1,2,3,4)
end

test_is_horizontal_or_vertical()

function draw_grid(startxs, startys, endxs, endys)
    xmax = maximum([maximum(startxs), maximum(endxs)])
    ymax = maximum([maximum(startys), maximum(endys)])
    grid =  zeros(Int, xmax, ymax)
    # assume all arrays are same length
    for i in 1:length(startxs)
        if ! is_horizontal_or_vertical(startxs[i],
                                       startys[i],
                                       endxs[i],
                                       endys[i])
            continue
        end

        for ix in minimum([startxs[i], endxs[i]]):maximum([startxs[i], endxs[i]])
            for iy in minimum([startys[i], endys[i]]):maximum([startys[i], endys[i]])
                grid[ix,iy] += 1
            end
        end
    end
    return grid
end

inputs = readlines("input.txt")

starts = map(x->split(x)[1], inputs)
ends = map(x->split(x)[3], inputs)

startxs = map(x->parse(Int,split(x,",")[1]), starts)
startys = map(x->parse(Int, split(x,",")[2]), starts)
endxs = map(x->parse(Int, split(x,",")[1]), ends)
endys = map(x->parse(Int, split(x,",")[2]), ends)

grid = draw_grid(startxs, startys, endxs, endys)

overlaps = count(x->(x>1), grid)
print("Count of overlaps ", overlaps)


text:
"""
    VAREstimate

# Fields
- `Y::Union{DataFrame,TimeArray}`
- `X::Union{DataFrame,TimeArray,Nothing}`
- `ynames::Array{String}`
- `xnames::Array{String}`
- `lags::Int`
- `constant::Bool`
- `trend::Bool`
- `obs::Int`
- `Z::Matrix`
- `B::Matrix`
- `seB::Matrix`
- `U::Matrix`
- `ΣU::Matrix`
- `Yhat::Matrix`
"""
struct VAREstimate
    Y::Union{DataFrame,TimeArray}
    X::Union{DataFrame,TimeArray,Nothing}
    ynames::Array{String}
    xnames::Array{String}
    lags::Int
    constant::Bool
    trend::Bool
    obs::Int
    Z::Matrix
    B::Matrix
    seB::Matrix
    U::Matrix
    ΣU::Matrix
    Yhat::Matrix
end

"""
    VAR(endog, lags; constant = true, trend = false, exog = nothing) -> VAREstimate

Estimate an unrestricted vector autoregression (VAR) using OLS.

# Arguments
- `endog` : `DataFrame` or `TimeArray` of observations on endogenous variables
- `lags::Int` : the number of lags
- `constant::Bool = true` : include an intercept term
- `trend::Bool = false` : include a linear trend
- `exog` : `DataFrame` or `TimeArray` of observations on exogenous variables
"""
function VAR(
    endog::DataFrame,
    lags;
    constant::Bool = true,
    trend::Bool = false,
    exog::Union{DataFrame,Nothing} = nothing,
)
    if lags > size(endog, 1)
        error("Number of lags is greater than number of observations.")
    end
    if !isa(exog, Nothing) && size(exog, 1) !== size(endog, 1)
        error("Unequal number of observations for exogenous and endogenous variables.")
    end
    ynames = String.(names(endog))
    xnames = !isa(exog, Nothing) ? String.(names(exog)) : [""]
    datamat = Matrix(endog)
    xmat = !isa(exog, Nothing) ? Matrix(exog) : nothing
    obs = size(datamat, 1) - lags
    Z, B, U, seB, ΣU, Yhat = varols(datamat, lags, constant, trend, xmat)
    VAREstimate(
        endog,
        exog,
        ynames,
        xnames,
        lags,
        constant,
        trend,
        obs,
        Z,
        B,
        seB,
        U,
        ΣU,
        Yhat,
    )
end

function VAR(
    endog::TimeArray,
    lags;
    constant::Bool = true,
    trend::Bool = false,
    exog = nothing,
)
    endog = DataFrame(endog)[:, 2:end]
    if !isa(exog, Nothing)
        exog = DataFrame(exog)[:, 2:end]
    end
    VAR(endog, lags; constant = constant, trend = trend, exog = exog)
end

function varols(y, ylag, constant, trend, x)
    # Set up RHS Matrix
    Z = rhs_matrix(y, ylag, constant, trend, x)
    # Set up LHS Matrix
    Y = y[ylag+1:end, :]
    # Coefficient estimates: each column of B corresponds to an equation
    B = (Z' * Z) \ (Z' * Y)
    Yhat = Z * B
    U = Y - Yhat
    ΣU = (U' * U) / (size(Y, 1) - size(B, 1))
    ΣB = kron(ΣU, inv(Z' * Z))  # see Lutkepohl p.80
    seB = sqrt.(reshape(diag(ΣB), size(B, 1), size(B, 2)))
    return Z, B, U, seB, ΣU, Yhat
end

function show(io::IO, v::VAREstimate)
    ct = coeftable(v)
    K = length(ct)
    println(io, typeof(v))
    println(io)
    print(io, "Endogenous variables: ")
    for k = 1:K-1
        print(io, "$(v.ynames[k]), ")
    end
    println(io, "$(v.ynames[K])")
    print(io, "Deterministic variables: ")
    if !isa(v.X, Nothing)
        for xn = 1:length(v.xnames)
            print(io, "$(v.xnames[xn]), ")
        end
    end
    v.constant && v.trend ? println(io, "intercept, trend") :
    (v.constant ? println(io, "intercept") : println(io))
    println(io, "Lags: $(v.lags)")
    println(io, "Sample size: $(v.obs)")
    for k = 1:K
        println(io)
        println(io, "Estimates for equation $(v.ynames[k]):")
        show(io, ct[k])
        println(io)
    end
end


text:
import Base: ==
==(x,y,z,ws...) = x==y && ==(y,z,ws...)

"""
Return the type's fields as a tuple
"""
@generated fieldvalues(x) = Expr(:tuple, (:(x.$f) for f=fieldnames(x))...)
@generated fields(x) = Expr(:tuple, (:($f=x.$f) for f=fieldnames(x))...)
firstfield(x) = first(fieldvalues(x))

"""
Rewrites `@! x = f(args...)` to `x = f!(x,args...)`

Special cases for `*` and `\\` forward to `mul!` and `ldiv!`, respectively.
"""
macro !(ex)
    if @capture(ex, x_ = f_(args__; kwargs_...))
        esc(:($(Symbol(string(f,"!")))($x,$(args...); $kwargs...)))
    elseif @capture(ex, x_ = f_(args__))
        if f == :*
            f = :mul
        elseif f==:\
            f = :ldiv
        end
        esc(:($x = $(Symbol(string(f,"!")))($x,$(args...))::typeof($x))) # ::typeof part helps inference sometimes
    else
        error("Usage: @! x = f(...)")
    end
end

nan2zero(x::T) where {T} = !isfinite(x) ? zero(T) : x
nan2zero(x::Diagonal{T}) where {T} = Diagonal{T}(nan2zero.(x.diag))
nan2inf(x::T) where {T} = !isfinite(x) ? T(Inf) : x

""" Return a tuple with the expression repeated n times """
macro repeated(ex,n)
    :(tuple($(repeated(esc(ex),n)...)))
end

"""
Pack some variables in a dictionary

```julia
> x = 3
> y = 4
> @dictpack x y z=>5
Dict(:x=>3,:y=>4,:z=>5)
```
"""
macro dictpack(exs...)
    kv(ex::Symbol) = :($(QuoteNode(ex))=>$(esc(ex)))
    kv(ex) = isexpr(ex,:call) && ex.args[1]==:(=>) ? :($(QuoteNode(ex.args[2]))=>$(esc(ex.args[3]))) : error()
    :(Dict($((kv(ex) for ex=exs)...)))
end

"""
Pack some variables into a NamedTuple. E.g.:

```julia
> x = 3
> y = 4
> @namedtuple(x, y, z=5)
(x=3,y=4,z=5)
```
"""
macro namedtuple(exs...)
    if length(exs)==1 && isexpr(exs[1],:tuple)
        exs = exs[1].args
    end
    kv(ex::Symbol) = :($(esc(ex))=$(esc(ex)))
    kv(ex) = isexpr(ex,:(=)) ? :($(esc(ex.args[1]))=$(esc(ex.args[2]))) : error()
    Expr(:tuple, (kv(ex) for ex=exs)...)
end

# these allow pinv and sqrt of SMatrices of Diagonals to work correctly, which
# we use for the T-E block of the covariance. hopefully some of this can be cut
# down on in the futue with some PRs into StaticArrays.
permutedims(A::SMatrix{2,2}) = @SMatrix[A[1] A[3]; A[2] A[4]]
function sqrt(A::SMatrix{2,2,<:Diagonal})
    a,b,c,d = A
    s = @. sqrt(a*d-b*c)
    t = pinv(@. sqrt(a+d+2s))
    @SMatrix[t*(a+s) t*b; t*c t*(d+s)]
end
function pinv(A::SMatrix{2,2,<:Diagonal})
    a,b,c,d = A
    idet = pinv(@. a*d-b*c)
    @SMatrix[d*idet -(b*idet); -(c*idet) a*idet]
end

# some usefule tuple manipulation functions:

# see: https://discourse.julialang.org/t/efficient-tuple-concatenation/5398/10
# and https://github.com/JuliaLang/julia/issues/27988
@inline tuplejoin(x) = x
@inline tuplejoin(x, y) = (x..., y...)
@inline tuplejoin(x, y, z...) = (x..., tuplejoin(y, z...)...)

# see https://discourse.julialang.org/t/any-way-to-make-this-one-liner-type-stable/10636/2
using Base: tuple_type_cons, tuple_type_head, tuple_type_tail, first, tail
map_tupleargs(f,::Type{T}) where {T<:Tuple} =
    (f(tuple_type_head(T)), map_tupleargs(f,tuple_type_tail(T))...)
map_tupleargs(f,::Type{T},::Type{S}) where {T<:Tuple,S<:Tuple} =
    (f(tuple_type_head(T),tuple_type_head(S)), map_tupleargs(f,tuple_type_tail(T),tuple_type_tail(S))...)
map_tupleargs(f,::Type{T},s::Tuple) where {T<:Tuple} =
    (f(tuple_type_head(T),first(s)), map_tupleargs(f,tuple_type_tail(T),tail(s))...)
map_tupleargs(f,::Type{<:Tuple{}}...) = ()
map_tupleargs(f,::Type{<:Tuple{}},::Tuple) = ()

# returns the base parametric type with all type parameters stripped out
basetype(::Type{T}) where {T} = T.name.wrapper
@generated function basetype(t::UnionAll)
    unwrap_expr(s::UnionAll, t=:t) = unwrap_expr(s.body, :($t.body))
    unwrap_expr(::DataType, t) = t
    :($(unwrap_expr(t.parameters[1])).name.wrapper)
end

function ensuresame(args...)
    @assert all(args .== Ref(args[1]))
    args[1]
end

tuple_type_len(::Type{<:NTuple{N,Any}}) where {N} = N

ensure1d(x::Union{Tuple,AbstractArray}) = x
ensure1d(x) = (x,)

# see https://discourse.julialang.org/t/dispatching-on-the-result-of-unwrap-unionall-seems-weird/25677
# for why we need this
# to use, just decorate the custom show_datatype with it, and make sure the args
# are named `io` and `t`.
macro show_datatype(ex)
    def = splitdef(ex)
    def[:body] = quote
        isconcretetype(t) ? $(def[:body]) : invoke(Base.show_datatype, Tuple{IO,DataType}, io, t)
    end
    esc(combinedef(def))
end

"""
    # symmetric in any of its final arguments except for bar:
    @sym_memo foo(bar, @sym(args...)) = <body>
    # symmetric in (i,j), but not baz
    @sym_memo foo(baz, @sym(i, j)) = <body>

The `@sym_memo` macro should be applied to a definition of a function
which is symmetric in some of its arguments. The arguments in which its
symmetric are specified by being wrapping them in @sym, and they must come at
the very end. The resulting function will be memoized and permutations of the
arguments which are equal due to symmetry will only be computed once.
"""
macro sym_memo(funcdef)

    sfuncdef = splitdef(funcdef)

    asymargs = sfuncdef[:args][1:end-1]
    symargs = collect(@match sfuncdef[:args][end] begin
        Expr(:macrocall, [head, _, ex...]), if head==Symbol("@sym") end => ex
        _ => error("final argument(s) should be marked @sym")
    end)
    sfuncdef[:args] = [asymargs..., symargs...]

    sfuncdef[:body] = quote
        symargs = [$(symargs...)]
        sorted_symargs = sort(symargs)
        if symargs==sorted_symargs
            $((sfuncdef[:body]))
        else
            $(sfuncdef[:name])($(asymargs...), sorted_symargs...)
        end
    end

    esc(:(@memoize $(combinedef(sfuncdef))))

end

@doc doc"""
```
@subst sum(x*$(y+1) for x=1:2)
```

becomes

```
let tmp=(y+1)
    sum(x*tmp for x=1:2)
end
```

to aid in writing clear/succinct code that doesn't recompute things
unnecessarily.
"""
macro subst(ex)

    subs = []
    ex = postwalk(ex) do x
        if isexpr(x, Symbol(raw"$"))
            var = gensym()
            push!(subs, :($(esc(var))=$(esc(x.args[1]))))
            var
        else
            x
        end
    end

    quote
        let $(subs...)
            $(esc(ex))
        end
    end

end

"""
    @invokelatest expr...

Rewrites all non-broadcasted function calls anywhere within an expression to use
`Base.invokelatest`. This means functions can be called that have a newer world
age, at the price of making things non-inferrable.
"""
macro invokelatest(ex)
    function walk(x)
        if isdef(x)
            x.args[2:end] .= map(walk, x.args[2:end])
            x
        elseif @capture(x, f_(args__; kwargs__)) && !startswith(string(f),'.')
            :(Base.invokelatest($f, $(map(walk,args)...); $(map(walk,kwargs)...)))
        elseif @capture(x, f_(args__)) && !startswith(string(f),'.')
            :(Base.invokelatest($f, $(map(walk,args)...)))
        elseif isexpr(x)
            x.args .= map(walk, x.args)
            x
        else
            x
        end
    end
    esc(walk(ex))
end

"""
    @ondemand(Package.function)(args...; kwargs...)
    @ondemand(Package.Submodule.function)(args...; kwargs...)

Just like calling `Package.function` or `Package.Submodule.function`, but
`Package` will be loaded on-demand if it is not already loaded. The call is no
longer inferrable.
"""
macro ondemand(ex)
    get_root_package(x) = @capture(x, a_.b_) ? get_root_package(a) : x
    quote
        @eval import $(get_root_package(ex))
        (args...; kwargs...) -> Base.invokelatest($(esc(ex)), args...; kwargs...)
    end
end

function tmap(f, args...)
    @static if nthreads()==1 || VERSION<v"1.3"
        map(f, args...)
    else
        map(fetch,map((args...)->Threads.@spawn(f(args...)),args...))
    end
end

get_kwarg_names(func::Function) = Vector{Symbol}(kwarg_decl(first(methods(func)), typeof(methods(func).mt.kwsorter)))
kwarg_decl(m::Method,kw::DataType) = VERSION<=v"1.3.999" ? Base.kwarg_decl(m,kw) : Base.kwarg_decl(m)

# maps a function recursively across all arguments of a Broadcasted expression,
# using the function `broadcasted` to reconstruct the `Broadcasted` object at
# each point.
map_bc_args(f, bc::Broadcasted) = broadcasted(bc.f, map(arg->map_bc_args(f, arg), bc.args)...)
map_bc_args(f, arg) = f(arg)

# adapting a closure adapts the captured variables
# this could probably be a PR into Adapt.jl
@generated function adapt_structure(to, f::F) where {F<:Function}
    if fieldcount(F) == 0
        :f
    else
        quote
            captured_vars = $(Expr(:tuple, (:(adapt(to, f.$x)) for x=fieldnames(F))...))
            $(Expr(:new, :($(F.name.wrapper){map(typeof,captured_vars)...}), (:(captured_vars[$i]) for i=1:fieldcount(F))...))
        end
    end
end

adapt_structure(to, d::Dict) = Dict(k => adapt(to, v) for (k,v) in d)


text:
@testset "Random" begin
    @testset "rand $T" for T in (Arf, Arb, Acb)
        @test rand(T) isa T
        @test precision(rand(T)) == precision(T)
        @test precision(rand(T(prec = 128))) == 128
    end
    a = rand(Arb(prec = 128))
    @test rand(Arblib.midref(a)) isa Arf
    @test precision(rand(Arblib.midref(a))) == precision(a)

    z = Acb(rand(Arb(prec = 128), 2)...)
    @test rand(Arblib.realref(z)) isa Arb
    @test precision(rand(Arblib.realref(z))) == precision(z)

    V = AcbVector(rand(Acb(prec = 128), 3))
    @test rand(Arblib.ref(V, 1)) isa Acb
    @test precision(rand(Arblib.ref(V, 1))) == precision(V)
end


text:
# If want to use an empty function for normalizing/testing validity of parameters
noop(x::AbstractArray) = return true
noop(x::Real) = return true

# Copy population
import Base.copy
copy(x::abc_population_type) = abc_population_type(copy(x.theta),copy(x.weights),copy(x.dist), x.epsilon, copy(x.repeats))

# Compute weighted covariance of sample (from origin)
function cov_weighted(x::AbstractArray{Float64,2}, w::AbstractArray{Float64,1} )
  @assert(size(x,1)==length(w) )
  sumw = sum(w)
  @assert( sumw > 0. )
  if(sumw!= 1.0)
     w /= sum(w)
     sumw = 1.0
  end
  sumw2 = sum(w.*w)
  xbar = [ sum(x[:,i].*w) for i in 1:size(x,2) ]
  covar = zeros(size(x,2),size(x,2))
  for k in 1:size(x,2)
    for j in 1:size(x,2)
        for i in 1:size(x,1)
            @inbounds covar[j,k] += (x[i,j]-xbar[j])*(x[i,k]-xbar[k]) *w[i]
        end
    @inbounds covar[j,k] *= sumw/(sumw*sumw-sumw2)
    end
  end
  covar
end

# Compute weighted variances of sample (from origin)
#function var_weighted(x::AbstractArray{Float64,2}, w::AbstractArray{Float64,1} )
function var_weighted(x::Array{Float64,2}, w::Array{Float64,1} )
  @assert(size(x,1)==length(w) )
  sumw = sum(w)
  @assert( sumw > 0. )
  if(sumw!= 1.0)
     w /= sum(w)
     sumw = 1.0
  end
  sumw2 = sum(w.*w)
  xbar = [ sum(x[:,i].*w) for i in 1:size(x,2) ]
  covar = zeros(size(x,2))
    for j in 1:size(x,2)
        for i in 1:size(x,1)
            @inbounds covar[j] += (x[i,j]-xbar[j])*(x[i,j]-xbar[j]) *w[i]
        end
    @inbounds covar[j] *= sumw/(sumw*sumw-sumw2)
    end
  covar
end

function make_matrix_pd(A::AbstractArray{Float64,2}; epsabs::Float64 = 0.0, epsfac::Float64 = 1.0e-6)
  @assert(size(A,1)==size(A,2))
  println("make_matrix_pd input: ", A)
  flush(Base.stdout)
  B = 0.5*(A+A')
  itt = 1
  while !isposdef(B)
	eigvalB,eigvecB = eig(B)
        pos_eigval_idx = eigvalB.>0.0
	neweigval = (epsabs == 0.0) ? epsfac*minimum(abs.(eigvalB[pos_eigval_idx])) : epsabs
	eigvalB[.!pos_eigval_idx] = neweigval
	B = eigvecB *diagm(eigvalB)*eigvecB'
	println(itt,": ",B)
        flush(Base.stdout)
        #cholB = chol(B)
	itt +=1
	if itt>size(A,1)
	  error("There's a problem in make_matrix_pd.\n")
  	  break
	end
  end
  return B
end

# Compute Effective Sample Size for array of weights
function ess(w::AbstractArray{Float64,1})
  sumw = sum(w)
  sumw2 = sum(w.*w)
  return sumw*sumw/sumw2
end

# Compute Effective Sample Size for array of weights, ignorring any repeated elements
function ess(w::AbstractArray{Float64,1}, repeat::AbstractArray{Int64,1} )
  @assert(length(w)==length(repeat))
  sumw = sum(w[findall(x->x==0,repeat)])
  sumw2 = sum(w[findall(x->x==0,repeat)].^2)
  return sumw*sumw/sumw2
end

# Common summary stats and distances
function calc_summary_stats_mean_var(x::AbstractArray{T,2} )  where T <:Real
 @inbounds m = mean(x,dims=2);
 @inbounds v = var(x,dims=2,mean=vec(m));
 return vcat(m, v)
end

function calc_summary_stats_mean_var(x::AbstractArray{T,1} )  where T <:Real
 @inbounds m = mean(x);
 @inbounds v = var(x,mean=m);
 return vcat(m, v)
end

calc_dist_max(x::AbstractArray{Float64},y::AbstractArray{Float64}) = maximum(abs.(x.-y))
dist_scale = Array{Float64}(undef,0)
calc_scaled_dist_max(x::AbstractArray{Float64,1},y::AbstractArray{Float64,1}, scale::AbstractArray{Float64,1} = dist_scale) = maximum(abs.(x.-y)./scale)

function set_distance_scale(ds::Array{Float64,1})
  global dist_scale
  dist_scale = copy(ds)
end
function set_distance_scale( plan::abc_pmc_plan_type, theta::AbstractArray{Float64,1}; num_draw::Integer = 40,
                  ss::AbstractArray{T,1} = plan.calc_summary_stats(plan.gen_data(theta)) ) where T
  dist = Array{Float64}(undef,length(ss),num_draw)
  for i in 1:num_draw
    dist[:,i] = abs.(ss.-plan.calc_summary_stats(plan.gen_data(theta)))
  end
  set_distance_scale(vec(median(dist,dims=2)))
end

# Summary Stats & Distances based on Empirical CDF
# using Distributions
calc_summary_stats_ecdf(x::AbstractArray) = x # EmpiricalUnivariateDistribution(x)

#function calc_dist_ks(x::EmpiricalUnivariateDistribution, y::EmpiricalUnivariateDistribution)
function calc_dist_ks(x::AbstractArray{T1,1}, y::AbstractArray{T2,1}) where {T1<:Real, T2<:Real}
   maxd = 0.0
   cdf_x = ecdf(x)
   cdf_y = ecdf(y)
   #for v in x.values
   for v in x
     #d = abs(cdf(x,v)-cdf(y,v))
     d = abs(cdf_x(v)-cdf_y(v))
     if d>maxd maxd = d end
   end
   #for v in y.values
   for v in y
     #d = abs(cdf(x,v)-cdf(y,v))
     d = abs(cdf_x(v)-cdf_y(v))
     if d>maxd maxd = d end
   end
   maxd
end

function normalize_vector_of_log_values(x::Vector{Float64})
      logsum = logsumexp(x)
      param_vector = exp(param_vector-logsum)
end


text:
import Base: *, sign
using Dates

const Price = Float64    # quote bid/ask, traded price
const Volume = Float64   # trade volume / number of shares
const Return = Price     # same as price

@enum TradeMode::Int16 LongShort=0 LongOnly=1 ShortOnly=-1
@enum TradeDir::Int16 Undefined=0 Long=1 Short=-1

@inline sign(x::TradeDir) = Int64(x)

@inline *(x::Volume, dir::TradeDir) = Volume(x * sign(dir))
@inline *(dir::TradeDir, x::Volume) = Volume(x * sign(dir))

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

struct Instrument
    symbol  ::String
    data    ::Any       # user-defined data field for instrument
    __hash  ::UInt64    # precomputed/cached hash
    Instrument(symbol; data::Any=nothing) = new(symbol, data, hash(symbol))
end

# custom hash (cached hash)
Base.hash(inst::Instrument) = inst.__hash

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

struct BidAsk
    dt      ::DateTime
    bid     ::Price
    ask     ::Price
    BidAsk() = new(DateTime(0), Price(0.0), Price(0.0))
    BidAsk(dt::DateTime, bid::Price, ask::Price) = new(dt, bid, ask)
end

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

@enum CloseReason::Int64 Unspecified=0 StopLoss=1 TakeProfit=2

mutable struct Position
    inst            ::Instrument
    size            ::Volume        # negative = short selling
    dir             ::TradeDir
    open_quote      ::BidAsk
    open_dt         ::DateTime
    open_price      ::Price
    last_quote      ::BidAsk
    last_dt         ::DateTime
    last_price      ::Price
    stop_loss       ::Price
    take_profit     ::Price
    close_reason    ::CloseReason
    pnl             ::Price
    data            ::Any           # user-defined data field for position instance
end

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

abstract type Order
end

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

struct OpenOrder <: Order
    inst            ::Instrument
    size            ::Volume
    dir             ::TradeDir
    stop_loss       ::Price
    take_profit     ::Price
    data            ::Any           # user-defined data field for position instance
    OpenOrder(
        inst::Instrument,
        size::Volume,
        dir::TradeDir;
        stop_loss::Price=NaN,
        take_profit::Price=NaN,
        data::Any=nothing) = new(inst, size, dir, stop_loss, take_profit, data)
end

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

struct CloseOrder <: Order
    pos             ::Position
    close_reason    ::CloseReason
    CloseOrder(pos::Position) = new(pos, Unspecified::CloseReason)
    CloseOrder(pos::Position, close_reason::CloseReason) = new(pos, close_reason)
end

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

struct CloseAllOrder <: Order
end

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

mutable struct Account
    open_positions      ::Vector{Position}
    closed_positions    ::Vector{Position}
    initial_balance     ::Price
    balance             ::Price
    equity              ::Price
    data                ::Any           # user-defined data field for account instance

    Account(initial_balance::Price; data::Any=missing) = new(
        Vector{Position}(),
        Vector{Position}(),
        initial_balance,
        initial_balance,
        initial_balance,
        data)
end

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


text:
#=
mdEmbed.jl:

a Julia based alternative for running multi-dimensional recurrence analysis.

this function is used for multidimensional embedding. It takes a n-dimensional matrix and outputs an n*D dimensional Dataset. Important that Dataset format is state of art for DynamicalSystems functions.

example:

emb_series = mdEmbed(data, D, τ)

using DynamicalSystems
RecurrenceMatrix(emb_series, .05; fixedrate = true)

requires:
Pkg.add("ElasticArrays"); # efficent matrices
using ElasticArrays;
=#

function mdEmbed(data, D, τ)

    if typeof(data) <: Dataset
        data = data |> Matrix
    end

    nrows = size(data,1)
    embed_series = ElasticArray{Float64}(undef, nrows-(D-1)*τ, 0);

    for i = 1:D
        append!(embed_series,data[1+(i-1)*τ:nrows-(D-i)*τ,:])
    end
    return(Dataset(embed_series))
end


text:
module SaddleSearch

using Parameters

export run!

# logging, preconditioner transformation, weidth norms and dots,
include("misc.jl")

# some line-search related methods
include("linesearch.jl")

# numerical ode dolvers
include("ode.jl")

# =========== Walker-type saddle search methods ============

include("dimer_types.jl")

# dimer, bb, ode
include("dimer.jl")

include("superlineardimer.jl")

# ============ Sting type methods ===================
# string, neb type definitions
include("string_types.jl")

# path manipulation: linear algebra in the path's framework
include("pathmisc.jl")

include("momentum_descent_schemes.jl")

# string method
include("string.jl")

# NEB method
include("neb.jl")

end # module


text:
"""
    SineAudioGenerator

An audio generator for a simple sine wave.  The waveform will be continuous even
if the rate (frequency) is modified during playback.
"""
mutable struct SineAudioGenerator <: AbstractAudioGenerator
    phase::Float64
    amp::Float64
    rate::Float64
end

function reset!(gen::SineAudioGenerator)
    gen.phase = 0
    nothing
end
function has_samples(gen::SineAudioGenerator, sample_rate::Float64)
    true
end
function next_samples!(gen::SineAudioGenerator, max_n::Int,
                       sample_rate::Float64)
    dphase = gen.rate / sample_rate
    pi2 = 2*pi
    samples = zeros(Float64, max_n, 2)
    p = gen.phase
    for i in 1:max_n
        samples[i, 1] = samples[i, 2] = gen.amp * sin(pi2*p)
        p += dphase
    end
    gen.phase = p
    samples
end


text:
using MD5
function restricted_md5(fname)
    if (!isfile(fname))
        return randn(10);
    end
    f=open(fname,"r");
    lines=readlines(fname);
    for (i,line) in enumerate(lines)
        if occursin("Created: ", line)
            lines[i]=""
        end
    end
    md5(join(lines,"\n"));
end

include("exp_set_list_all.jl");

println("Saving simulation");

rho=mono.params[:rho];
m=count(values(state_list[1].graph.operations) .== :mult)

println("Domain: $(rho), m=$m",);
for (i,state)=enumerate(state_list)
    if (!isnothing(state))
        name=state.params[:graphname]
        @show name
        if contains(name,"_opt") || contains(name,"ps") || contains(name,"mono")
            graph=state.graph;
            rho_str=replace(string(rho),"."=>"_");
            fname=joinpath("..","..","data","exp","exp_$(name)_rho$(rho_str).cgr");
            err=Float64(showerr(state,output=false,n=state.params[:target_n]*2))

            isreal=false;
            if (norm(imag.(get_coeffs(graph)))==0)
                isreal=true;
                g=Compgraph(real(eltype(graph)),graph);
            end
            println("Saving $name $fname isreal=$isreal");
            optimized=""
            if (contains(name,"_opt"))
                optimized=" optimized starting from "*replace(mono.params[:graphname], "_opt" => "");
            end

            tmpname = "/tmp/tmpfile.jl"
            export_compgraph(g, tmpname;
                             fun="exp", dom="$(rho)D", err="$err",
                             order=get_topo_order(g)[1],
                             descr="Matrix exponential approximation with degree-optimal form $(optimized). This graph contains redundancies. Run compress_graph! before generating code.",
                             user="Elias Jarlebring");

            # Only save it if the file is actually different (beside time-stamp)
            if (restricted_md5(tmpname) != restricted_md5(fname))
                cp(tmpname,fname,force=true);
            else
                println("Skip saving $name (since file unchanged)");
            end

        else
            println("Skip saving $name (since it is independent of the disk)");
        end

    else
        println("Skip saving $i");
    end
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
infile = "./src/Day17/input.txt"
input = readlines(infile)

using DataStructures

cube = Dict{CartesianIndex, Bool}()
for x = 1:length(input[1]), y in 1:length(input)
    z = 0
    if input[y][x] == '#'
        cube[CartesianIndex(x,y,z)] = true
    end
end

neighbourgrid = CartesianIndices((-1:1, -1:1, -1:1))

function countactiveneighbours(cube, loc, neighbourgrid)
    grid = loc .+ neighbourgrid
    neighbours = sum(cube[i] for i in grid if haskey(cube, i))
    if haskey(cube, loc)
        count = cube[loc] ? neighbours - 1 : neighbours
    else
        count = neighbours
    end
    return count
end

test = countactiveneighbours(cube, CartesianIndex(8,6,0), neighbourgrid)

function getnext(cube, neighbourgrid)
    nextcube = Dict{CartesianIndex, Bool}()
    for centre in keys(cube)
        # println("Centre is $centre")
        # ps = centre .+ neighbourgrid
        for p in (centre .+ neighbourgrid)
            # println("Current point: $p")
            if haskey(nextcube, p)
                # already pre-calculated
                # println("Pre-done")
                continue
            end
            mystate = haskey(cube, p) ? cube[p] : false
            neighbourcount = countactiveneighbours( cube, p, neighbourgrid)
            if mystate
                if (neighbourcount == 2) || ( neighbourcount == 3)
                    active = true
                else
                    active = false
                end
            else
                if neighbourcount == 3
                    active = true
                else
                    active = false
                end
            end
            nextcube[p] = active
        end
    end
    return nextcube
end

t = getnext(cube, neighbourgrid)

cycles = 6
function part1(cube, neighbourgrid, cycles)
    for i in 1:cycles
        cube = getnext(cube, neighbourgrid)
    end
    return cube
end

cube
sixthcube = part1(cube, neighbourgrid, cycles)
active6 = sum(v for v in values(sixthcube))

neighbourgrid2 = CartesianIndices((-1:1, -1:1, -1:1, -1:1))

cube2 = Dict{CartesianIndex, Bool}()
for x = 1:length(input[1]), y in 1:length(input)
    z = 0
    w = 0
    if input[y][x] == '#'
        cube2[CartesianIndex(x,y,z,w)] = true
    end
end

using BenchmarkTools
@profview part2 = part1(cube2, neighbourgrid2, cycles)
activepart2 = sum(v for v in values(part2))
@benchmark part1(cube2, neighbourgrid2, 3)


text:
NI(x...) = error("This function is not implemented.")

# filemap is filled in the format specific source files
const filemap = Dict{Symbol, Tuple{Function, Function, Function, Function}}()
        # :gml        => (readgml, writegml, readnetgml, writenetgml)

"""
    readgraph(filename, G=Graph)
    readgraph(filename, t, G=Graph; compressed=false)

Reads a graph from  `filename` in the format `t`. Returns a graph of type `G`
or the corresponding digraph/graph type.
Compressed files can eventually be read.

Supported formats are `:gml, :dot, :graphml, :gexf, :net, :gt`.

If no format is provided, it will be inferred from `filename`.

    readgraph(s::Symbol, G=Graph)

Read a graph identified by `s` from Erdos datasets collection (e.g. `s=:karate`).
They are stored in the `gt` binary format in the `datasets` directory of the package.
For a list of available graph refer to the documentation.
"""
function readgraph(fn::String, t::Symbol, ::Type{G}=Graph; compressed=false) where G<:AGraphOrDiGraph
    t in keys(filemap) || error("Please select a supported graph format: one of $(keys(filemap))")
    io = compressed ? GZip.open(fn,"r") : open(fn,"r")
    g = readgraph(io, t, G)
    close(io)
    return g
end

function readgraph(io::IO, t::Symbol, ::Type{G}=Graph) where G<:AGraphOrDiGraph
    return filemap[t][1](io, G)
end

function readgraph(fn::String, ::Type{G}=Graph) where G<:AGraphOrDiGraph
    compressed = false
    ft = split(fn,'.')[end]
    if ft == "gz"
        compressed = true
        ft = split(fn,'.')[end-1]
    end
    Symbol(ft) in keys(filemap) || error("Could not infer file format.")
    return readgraph(fn, Symbol(ft), G; compressed=compressed)
end

const DATASETS_DIR = joinpath(Base.source_dir(), "..", "..", "datasets")
function readgraph(s::Symbol, ::Type{G}=Graph) where G<:AGraphOrDiGraph
    readgraph(joinpath(DATASETS_DIR, string(s) * ".gt.gz"), G)
end

"""
    readnetwork(filename, G=Network)
    readnetwork(filename, t, G=Network; compressed=false)

Read a network from  `filename` in the format `t`.
Returns a network of type `G`
(or the corresponding directed/undirected type if needed).
Compressed files can eventually be read.

Supported formats are `:gml, :dot, :graphml, :gexf, :net, :gt`.
When possible, graph, edge, and vertex properties will be read as well.

If no format is provided, it will be inferred from the `filename`.

    readnetwork(s::Symbol, G=Network)

Read a network identified by `s` from Erdos' datasets collection (e.g. `s=:karate`).
They are stored in the `gt` binary format in the `datasets` directory of the package.
For a list of available graph refer to the documentation.
"""
function readnetwork(fn::String, t::Symbol, ::Type{G}=Network; compressed=false) where G<:ANetOrDiNet
    t in keys(filemap) || error("Please select a supported graph format: one of $(keys(filemap))")
    io = compressed ? GZip.open(fn,"r") : open(fn,"r")
    g = readnetwork(io, t, G)
    close(io)
    return g
end

function readnetwork(io::IO, t::Symbol, ::Type{G}=Network) where G<:ANetOrDiNet
    return filemap[t][3](io, G)
end

function readnetwork(fn::String, ::Type{G}=Network) where G<:ANetOrDiNet
    compressed = false
    ft = split(fn,'.')[end]
    if ft == "gz"
        compressed = true
        ft = split(fn,'.')[end-1]
    end
    Symbol(ft) in keys(filemap) || error("Could not infer file format.")
    return readnetwork(fn, Symbol(ft), G; compressed=compressed)
end

function readnetwork(s::Symbol, ::Type{G}=Network) where G<:AGraphOrDiGraph
    return readnetwork(joinpath(DATASETS_DIR, string(s) * ".gt.gz"), G)
end

"""
    writegraph(file, g)
    writegraph(file, g, t; compress=false)

Save a graph `g` to `file` in the format `t`.

Eventually the resulting file can be compressed in the gzip format.

Currently supported formats are `:gml, :graphml, :gexf, :dot, :net, :gt`.

If no format is provided, it will be inferred from `file` along with compression.
"""
function writegraph(io::IO, g::AGraphOrDiGraph, t::Symbol)
    filemap[t][2](io, g)
    return nothing
end

function writegraph(fn::String, g::AGraphOrDiGraph, t::Symbol; compress::Bool=false)
    t in keys(filemap) || error("Please select a supported graph format: one of $(keys(filemap))")
    io = compress ? GZip.open(fn,"w") : open(fn,"w")
    writegraph(io, g, t)
    close(io)
    return nothing
end

function writegraph(fn::String, g::AGraphOrDiGraph)
    compress = false
    ft = split(fn,'.')[end]
    if ft == "gz"
        compress = true
        ft = split(fn,'.')[end-1]
    end
    Symbol(ft) in keys(filemap) || error("Could not infer file format.")
    writegraph(fn, g, Symbol(ft), compress=compress)
    return nothing
end

"""
    writenetwork(file, g)
    writenetwork(file, g, t; compress=false)

Save a network `g` to `file` in the format `t`.

Eventually the resulting file can be compressed in the gzip format.

Currently supported formats are `:gml, :graphml, :gexf, :dot, :net, :gt`.
When possible, graph, edge, and vertex properties will be written as well.

If no format is provided, it will be inferred from `file` along with compression.
"""
function writenetwork(io::IO, g::ANetOrDiNet, t::Symbol)
    filemap[t][4](io, g)
    return nothing
end

function writenetwork(fn::String, g::ANetOrDiNet, t::Symbol; compress::Bool=false)
    t in keys(filemap) || error("Please select a supported graph format: one of $(keys(filemap))")
    io = compress ? GZip.open(fn,"w") : open(fn,"w")
    writenetwork(io, g, t)
    close(io)
    return nothing
end

function writenetwork(fn::String, g::ANetOrDiNet)
    compress = false
    ft = split(fn,'.')[end]
    if ft == "gz"
        compress = true
        ft = split(fn,'.')[end-1]
    end
    Symbol(ft) in keys(filemap) || error("Could not infer file format.")
    writenetwork(fn, g, Symbol(ft), compress=compress)
    return nothing
end

function getchild(el::EzXML.Node, s::String)
    childs = elements(el)
    i = findfirst(x->nodename(x)==s, childs)
    i === nothing && error("no child $s")
    return childs[i]
end

"adds a child `s` to `el` if it doesn't exist"
function getchild!(el::EzXML.Node, s::String)
    childs = elements(el)
    i = findfirst(x->nodename(x)==s, childs)
    return i !== nothing ? childs[i] : addelement!(el, s)
end


text:
using Test
using LinearAlgebra
using OrdinaryDiffEq: ODEProblem, solve, SSPRK33, Tsit5

import ClimaCore:
    Fields,
    Domains,
    Topologies,
    Meshes,
    DataLayouts,
    Operators,
    Geometry,
    Spaces

convergence_rate(err, Δh) =
    [log(err[i] / err[i - 1]) / log(Δh[i] / Δh[i - 1]) for i in 2:length(Δh)]

@testset "1D ∂T/t = α ∇²T ODE solve" begin

    FT = Float64
    n_elems_seq = 2 .^ (5, 6, 7)
    err, Δh = zeros(length(n_elems_seq)), zeros(length(n_elems_seq))

    for (k, n) in enumerate(n_elems_seq)
        z₀ = FT(0)
        z₁ = FT(10)
        t₀ = FT(0)
        t₁ = FT(10)
        μ = FT(-1 / 2)
        ν = FT(5)
        𝓌 = FT(1)
        δ = FT(1)

        domain = Domains.IntervalDomain(z₀, z₁, x3boundary = (:bottom, :top))
        zp = (z₀ + z₁ / n / 2):(z₁ / n):(z₁ - z₁ / n / 2)

        function gaussian(z, t; μ = -1 // 2, ν = 1, 𝓌 = 1, δ = 1)
            return exp(-(z - μ - 𝓌 * t)^2 / (4 * ν * (t + δ))) / sqrt(1 + t / δ)
        end

        function ∇gaussian(z, t; μ = -1 // 2, ν = 1, 𝓌 = 1, δ = 1)
            return -2 * (z - μ - 𝓌 * t) / (4 * ν * (δ + t)) *
                   exp(-(z - μ - 𝓌 * t)^2 / (4 * ν * (δ + t))) / sqrt(1 + t / δ)
        end

        z₀ = FT(0)
        z₁ = FT(10)

        mesh = Meshes.IntervalMesh(domain, nelems = n)

        cs = Spaces.CenterFiniteDifferenceSpace(mesh)
        fs = Spaces.FaceFiniteDifferenceSpace(cs)
        zc = Fields.coordinate_field(cs)

        T = gaussian.(zc, -0; μ = μ, δ = δ, ν = ν, 𝓌 = 𝓌)
        V = ones(FT, fs)

        function ∑tendencies!(dT, T, z, t)
            bc_vb = Operators.SetValue(
                FT(gaussian(z₀, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ)),
            )
            bc_gt = Operators.SetGradient(
                FT(∇gaussian(z₁, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ)),
            )
            A = Operators.AdvectionC2C(
                bottom = bc_vb,
                top = Operators.Extrapolate(),
            )
            gradc2f = Operators.GradientC2F(bottom = bc_vb, top = bc_gt)
            gradf2c = Operators.GradientF2C()
            return @. dT = gradf2c(ν * gradc2f(T)) - A(V, T)
        end

        # Solve the ODE operator
        Δt = 0.001
        prob = ODEProblem(∑tendencies!, T, (t₀, t₁))
        sol = solve(
            prob,
            Tsit5(),
            reltol = 1e-8,
            abstol = 1e-8,
            #dt = Δt,
        )
        computed_result = sol.u[end]
        analytical_result = gaussian.(zp, t₁; μ = μ, δ = δ, ν = ν, 𝓌 = 𝓌)
        Δh[k] = cs.Δh_c2c[1]
        err[k] =
            norm(parent(computed_result) .- analytical_result) /
            length(analytical_result)
    end
    conv = convergence_rate(err, Δh)
    # conv should be approximately 2 for second order-accurate stencil.
    @test 1.4 ≤ conv[1] ≤ 2.1
    @test 1.4 ≤ conv[2] ≤ 2.1
    @test conv[1] ≤ conv[2]
    @test err[3] ≤ err[2] ≤ err[1] ≤ 1e-2
end


text:
module TestEnumerationArray

using FactCheck
using DataCubes
import DataCubes: wrap_array

facts("EnumerationArray tests") do
  @fact EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"])).elems --> Int[1 0 2 2;0 2 1 3]
  @fact EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"]).elems --> Int[2 0 3 3;0 3 2 1]
  @fact EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]))[2,2].value --> "b"
  @fact EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"])[2,2].value --> "b"
  @fact dcube.wrap_array(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]))[1:3]) --> dcube.wrap_array(EnumerationArray(@nalift(["a",NA,NA])))
  @fact dcube.wrap_array(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"])[1:3]) --> dcube.wrap_array(EnumerationArray(@nalift(["a",NA,NA]), ["c","a","b"]))
  @fact size(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"])) --> (2, 4)
  @fact dcube.wrap_array(reshape(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"]), 4, 2)) --> dcube.wrap_array(permutedims(EnumerationArray(@nalift(["a" NA  NA "b";"b" "a" "b" "c"]), ["c","a","b"]), [2,1]))
  @fact dcube.wrap_array(permutedims(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"]), [2,1])) --> dcube.wrap_array(EnumerationArray(@nalift(["a" NA;NA "b";"b" "a";"b" "c"])))
  @fact dcube.wrap_array(permutedims(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"]), [2,1])) --> dcube.wrap_array(permutedims(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"]), (2,1)))
  @fact dcube.wrap_array(reverse(enumeration(["hello", "hello", "hi"]))) --> dcube.wrap_array(enumeration(["hi", "hello", "hello"]))
  @fact dcube.wrap_array((a=enumeration([:a,:b,:c]);a[2]=Nullable{Symbol}();a)) --> dcube.wrap_array(@enumeration([:a,NA,:c],[:a,:b,:c]))
  @fact dcube.wrap_array((a=enumeration([:x,:y]);copy!(a, enumeration([:m,:n]));a)) --> dcube.wrap_array(enumeration([:m,:n]))
  @fact_throws enumeration([1 2.0 NA])
  @fact_throws enumeration([1 2.0 NA],[2.0,1])
  @fact dcube.wrap_array(sort(enumeration([6,5,3,1,2,1]))) --> nalift([6,5,3,1,1,2])
  @fact dcube.wrap_array((a=enumeration([6,5,3,1,2,1]);sort!(a);a)) --> nalift([6,5,3,1,1,2])
  @fact dcube.wrap_array((a=enumeration([6,5,3,1,2,1],1:6);sort!(a);a)) --> nalift([1,1,2,3,5,6])
  @fact peel(@enumeration([NA :x :y;:x :z NA])) --> [0 1 3;1 2 0]

  context("enumerations in labeled array tests") do
    arr = larr(
                :a=>enumeration(reshape(repeat(collect(1:4),inner=[1],outer=[5]),5,4)),
                :b=>2.0*nalift(reshape(1:20,5,4)),
                2=>enumeration(map(i->Symbol(:sym,i), reshape(repeat(collect(1:4),inner=[1],outer=[5]),5,4))),
                :third=>enumeration(map(i->string(:str,i), reshape(repeat(collect(1:4),inner=[1],outer=[5]),5,4))),
          axis1=DictArray(k1=nalift(["a","a","b","b","b"]), k2=nalift(collect(101:105))),
          axis2=DictArray(r1=nalift([:alpha,:beta,:gamma,:delta])))
    @fact arr[3,2] --> LDict(:a=>Nullable(4), :b=>Nullable(16.0), 2=>Nullable(:sym4), :third=>Nullable("str4"))
    @fact dcube.getindexpair(arr,3,2).second.second --> LDict(:a=>Nullable(4), :b=>Nullable(16.0), 2=>Nullable(:sym4), :third=>Nullable("str4"))
    @fact arr[1:3,2] --> LabeledArray(
                            DictArray(:a=>nalift([2,3,4]),
                                      :b=>nalift([12.0,14.0,16.0]),
                                      2=>nalift([:sym2,:sym3,:sym4]),
                                      :third=>nalift(["str2","str3","str4"])),
                            axis1=DictArray(k1=nalift(["a","a","b"]), k2=nalift(collect(101:103))))
    @fact typeof(arr[1:3,2].data.data[:a]) --> DataCubes.AbstractArrayWrapper{Nullable{Int64},1,DataCubes.EnumerationArray{Int64,1,DataCubes.AbstractArrayWrapper{Int64,1,Array{Int64,1}},Int64}}
    @fact arr[2:2,2:4] --> LabeledArray(
                            DictArray(:a=>nalift([3 4 1]),
                                      :b=>nalift([14.0 24.0 34.0]),
                                      2=>nalift([:sym3 :sym4 :sym1]),
                                      :third=>nalift(["str3" "str4" "str1"])),
                            axis1=DictArray(k1=nalift(["a"]), k2=nalift([102])),
                            axis2=DictArray(r1=nalift([:beta,:gamma,:delta])))
    @fact size(arr) --> (5,4)
    #@fact eltype(arr) --> LDict{Any}
    @fact pick(arr,:a)[1,1].value --> 1
    @fact pick(arr,:b)[1,1].value --> 2.0
    @fact pick(arr,(:a,))[1][1,1].value --> 1
    @fact pick(arr,(:b,))[1][1,1].value --> 2.0
    @fact endof(arr) --> length(arr)
    @fact permutedims(arr, [2,1]) --> LabeledArray(DictArray(mapvalues(x->permutedims(x,[2,1]), arr.data.data)),
                                           axis1=arr.axes[2],
                                           axis2=arr.axes[1])
    @fact permutedims(arr, (2,1)) --> permutedims(arr, [2,1])
    @fact copy(arr) !== arr --> true
    @fact size(cat(1, arr, arr)) --> ntuple(n->n==1 ? size(arr,1)*2 : size(arr,n), ndims(arr))
    @fact size(cat(2, arr, arr)) --> ntuple(n->n==2 ? size(arr,2)*2 : size(arr,n), ndims(arr))
    @fact size(repeat(arr, inner=[5,2], outer=[3,4])) --> (size(arr,1)*5*3, size(arr,2)*2*4)

    d1 = larr(:a=>enumeration(reshape(1:20, 5, 4)),
              2=>enumeration(fill(:sym, 5, 4)),
              :third=>enumeration(reshape(map(i->string("str_",i), 1:20), 5, 4)))
    d2 = larr(:x=>enumeration(reshape(1:20, 5, 4)),
              :third=>enumeration(fill(:sym, 5, 4)),
              :z=>enumeration(reshape(map(i->string("str_",i), 1:20), 5, 4)))
    @fact pick(merge(d1, d2), [:a,2]) --> pick(d1, [:a,2])
    @fact pick(merge(d1, d2), [:x,:third,:z]) --> pick(d2, [:x,:third,:z])
    @fact keys(peel(merge(d1, d2))) --> Any[:a,2,:third,:x,:z]
    @fact peel(delete(merge(d1, d2), :a,2,:third)) --> pick(d2, [:x, :z])
    @fact pick(d1, [:a]) --> DictArray(a=pick(d1, :a))
    @fact mapslices(x->LDict(:c=>Nullable(length(x))),@larr(a=enumeration([1 2 3;4 5 6]),b=["a" "b" "c";"d" "e" "f"],axis1[k=[:x,:y]],axis2[r=[:m,:n,:p]]),[1]) --> @larr(c=[2,2,2], axis1[r=[:m,:n,:p]])
    @fact mapslices(x->LDict(:c=>Nullable(length(x))),@larr(a=enumeration([1 2 3;4 5 6]),b=["a" "b" "c";"d" "e" "f"],axis1[k=[:x,:y]],axis2[r=[:m,:n,:p]]),[2]) --> @larr(c=[3,3], axis1[k=[:x,:y]])
    @fact mapslices(x->Nullable(length(x)),@larr(a=enumeration([1 2 3;4 5 6]),b=["a" "b" "c";"d" "e" "f"]),[1]) --> LabeledArray(@nalift([2,2,2]), axes1=@nalift([1,2,3]))
    @fact size(mapslices(identity, @larr(a=[1 2 3;4 5 6]), [])) --> (2,3)
    @fact mapslices(x->LDict(:c1=>DataCubes.naop_plus(x[:a],x[:b]),:c2=>Nullable(10)), @larr(a=enumeration([1 2 3;4 5 6]),b=[1.0 2.0 3.0;4.0 5.0 6.0]), []) --> @larr(c1=[2.0 4.0 6.0;8.0 10.0 12.0], c2=@rap reshape(_,(2,3)) fill (10,6)...)
    @fact mapslices(x -> x, @larr(a=enumeration([1 2 3;4 5 6])), []) --> @larr(a=[1 2 3;4 5 6])
    @fact map(x->x[:a], @larr(a=enumeration([1,2,3]),b=[4,5,6],axis1[[:x,:y,:z]])) --> @larr([1,2,3],axis1[[:x,:y,:z]])
    @fact dcube.create_dict(@larr(a=@enumeration([1 NA NA;4 5 6]),b=[NA NA 6; 7 8 9],axis2[r=[:x,:y,:z]]))[Nullable(2)][LDict(:r=>Nullable(:z))][:a].value --> 6
    @fact Set(keys(dcube.create_dict(@larr(a=@enumeration([1 NA NA;4 5 6]),b=[NA NA 6; 7 8 9],axis2[r=[:x,:y,:z]]))[Nullable(1)])) --> Set([LDict(:r=>Nullable(:x)), LDict(:r=>Nullable(:z))])
    @fact dcube.create_dict(@larr(a=[1 NA NA;4 5 6],b=[NA NA 6; 7 8 9],axis2[r=[:x,:y,:z]]))[Nullable(1)][LDict(:r=>Nullable(:x))] --> LDict(:a=>Nullable(1), :b=>Nullable{Int}())
    @fact reverse(@larr(a=enumeration([1 2 3;4 5 6]),axis2[r=enumeration([:x,:y,:z])]),[1]) --> @larr(a=[4 5 6;1 2 3],axis2[r=[:x,:y,:z]])
    @fact reverse(@larr(a=enumeration([1 2 3;4 5 6]),axis2[r=[:x,:y,:z]]),1:2) --> @larr(a=[6 5 4;3 2 1],axis2[r=[:z,:y,:x]])
    @fact flipdim(@larr(a=[1 2 3;4 5 6],axis2[r=enumeration([:x,:y,:z])]),2) --> @larr(a=[3 2 1;6 5 4],axis2[r=[:z,:y,:x]])
    @fact flipdim(@larr(a=enumeration([1 2 3;4 5 6]),axis2[r=[:x,:y,:z]]),1,2) --> @larr(a=[6 5 4;3 2 1],axis2[r=[:z,:y,:x]])
    @fact reshape(larr(a=enumeration([1 2 3;4 5 6]),b=[10 11 12;13 14 15],axis1=darr(k1=[:a,:b],k2=[100,101]),axis2=[:m,:n,:p]),6) --> larr(a=[1,4,2,5,3,6],b=[10,13,11,14,12,15], axis1=darr(k1=repmat([:a,:b],3),k2=repmat([100,101],3),x1=[:m,:m,:n,:n,:p,:p]))
    @fact reshape(larr(a=[1 2 3;4 5 6],b=enumeration([10 11 12;13 14 15]),axis1=darr(k1=[:a,:b],k2=[100,101]),axis2=[:m,:n,:p]),1,6) --> @rap permutedims(_,[2,1]) reshape(_,6,1) larr(a=[1,4,2,5,3,6],b=[10,13,11,14,12,15], axis1=darr(k1=repmat(enumeration([:a,:b]),3),k2=repmat([100,101],3),x1=[:m,:m,:n,:n,:p,:p]))
    @fact reshape(larr([1 2 3;4 5 6],axis1=[10,11]),1,6) --> larr([1 4 2 5 3 6], axis2=repmat([10,11],3))
    @fact reducedim((x,y)->x+y[:a].value, larr(a=enumeration(reshape(1:24,2,3,4))),[1],0) --> @larr([3 15 27 39;7 19 31 43;11 23 35 47])
    @fact reducedim((x,y)->x+y[:a].value, larr(a=enumeration(reshape(24:-1:1,2,3,4))),[1,2,3],0).value --> 300
    @fact reducedim((x,y)->x+y[:a].value, larr(a=enumeration(reshape(1:24,2,3,4)),axis1=darr(k1=[:a,:b]),axis3=[:x,:y,:z,:w]),[1,2,3],0).value --> 300
    # runnability test when the result is empty.
    @fact size(@rap permutedims(_,[2,1]) @select((@rap permutedims(_,[2,1]) larr(r=enumeration(rand(8,20)), axis1=darr(a=rand(8),b=101:108))), :r=broadcast(*, _r, 100), where[broadcast(>, _r, 1)])) --> (0,0)
    @fact cat(1, larr(a=enumeration([1 2 3;4 5 6]), b=['a' 'b' 'c';'d' 'e' 'f'], axis1=[:u,:v]), larr(b=['x' 'y' 'z'], d=[:m :n :p], axis1=[3])) --> larr(reshape(@larr(a=[1,4,NA,2,5,NA,3,6,NA],b=['a','d','x','b','e','y','c','f','z'],d=[NA,NA,:m,NA,NA,:n,NA,NA,:p]), 3, 3), axis1=[:u,:v,3])
    @fact cat(2, larr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f'], axis1=[:u,:v]), larr(b=['x','y'], d=[:m,:n], axis1=[:u,:v])) --> larr(reshape(@larr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f','x','y'], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4), axis1=[:u,:v])
    @fact cat(2, larr(a=[1 2 3;4 5 6], b=enumeration(['a' 'b' 'c';'d' 'e' 'f'])), larr(b=[10,11], d=[:m,:n])) --> reshape(@larr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f',10,11], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4)
    @fact cat(2, larr(a=enumeration([1 2 3;4 5 6]), b=['a' 'b' 'c';'d' 'e' 'f'], axis2=darr(r=[:x,:y,:z])), larr(b=[10,11], d=[:m,:n])) --> larr(reshape(@darr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f',10,11], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4), axis2=@darr(r=[:x,:y,:z,NA]))
    #@fact sub(@larr(a=enumeration([1 2 3;4 5 6]),axis2[r=[:x,:y,:z]]),[2,1],2) --> getindex(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),[2,1],2)
    #@fact sub(@larr(a=enumeration([1 2 3;4 5 6]),axis2[r=[:x,:y,:z]]),[2,1],2:3) --> getindex(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),[2,1],2:3)
    @fact view(@larr(a=enumeration([1 2 3;4 5 6]),axis2[r=[:x,:y,:z]]),[2,1],2) --> larr(a=[5,2])
    @fact view(@larr(a=enumeration([1 2 3;4 5 6]),axis2[r=[:x,:y,:z]]),1, 2:3) --> larr(a=[2,3], axis1=darr(r=[:y,:z]))
    @fact view(@larr(a=enumeration([1 2 3;4 5 6]),axis2[r=enumeration([:x,:y,:z])]),1, 2:3) --> larr(a=[2,3], axis1=darr(r=[:y,:z]))
    @fact endof(@enumeration([1,2,3,NA,5])) --> 5
    @fact DataCubes.wrap_array(reshape(@enumeration(6:-1:1),(2,3))) --> nalift([6 4 2;5 3 1])
    @fact nalift((a=enumeration([10,20,30,10,20]);a[1]=3;a)) --> nalift([30,20,30,10,20])
    @fact DataCubes.wrap_array((a=enumeration([1,2,3]);copy!(a, nalift([3,2,1]));a)) --> DataCubes.wrap_array(enumeration([3,2,1],[1,2,3]))
    @fact DataCubes.wrap_array((a=enumeration([1,2,3]);copy!(a, enumeration([3,2,1]));a)) --> DataCubes.wrap_array(enumeration([3,2,1]))
    @fact typeof(similar(@enumeration([1.0 2.0 NA]))) --> EnumerationArray{Float64,2,DataCubes.AbstractArrayWrapper{Int64,2,Array{Int64,2}},Int64}
    @fact typeof(similar(enumeration([:a :b]))) --> EnumerationArray{Symbol,2,DataCubes.AbstractArrayWrapper{Int64,2,Array{Int64,2}},Int64}
    @fact typeof(similar(enumeration([:a :b]),2,2,1)) --> EnumerationArray{Symbol,3,DataCubes.AbstractArrayWrapper{Int64,3,Array{Int64,3}},Int64}
    @fact typeof(similar(enumeration([:a :b]),(2,2,1))) --> EnumerationArray{Symbol,3,DataCubes.AbstractArrayWrapper{Int64,3,Array{Int64,3}},Int64}
    @fact_throws cat(1, enumeration([:a,:b]), enumeration([:c,:d]))
    @fact wrap_array(vcat(enumeration([:a :b :c],[:c,:b,:a,:d]), enumeration([:b :c :d],[:c,:b,:a,:d]), enumeration([:c :d :a],[:c,:b,:a,:d]))) --> wrap_array(enumeration([:a :b :c;:b :c :d;:c :d :a],[:c,:b,:a,:d]))
    @fact wrap_array(hcat(enumeration([:a :b :c],[:c,:b,:a,:d]), enumeration([:b :c :d],[:c,:b,:a,:d]), enumeration([:c :d :a],[:c,:b,:a,:d]))) --> wrap_array(enumeration([:a :b :c :b :c :d :c :d :a],[:c,:b,:a,:d]))
    # seems there is a bug in FactCheck. Comment it for now.
    #temp1 = enumeration([:a,:c,:b])
    #@show typeof(temp1)
    #@show temp1
    #temp2 = enumeration(temp1, [:a,:b])
    #@fact temp2.elems --> [1,3,2]
    #@fact temp2.pool --> [:a,:b,:c]
    #@fact enumeration(enumeration([:a,:c,:b]),[:a,:b]).elems --> [1,3,2]
    #@fact enumeration(enumeration([:a,:c,:b]),[:a,:b]).pool --> [:a,:b,:c]
  end
end

end


text:
# Geodesy

# WGS84
# etc

function getPosition()

end


text:
using PtFEM

data = Dict(
  # GenericSolid(ndim, nst, nels, nn, nip, finite_element(nod, nodof), axisymmentrix)
  :struc_el => GenericSolid(3, 6, 7, 40, 27, Hexahedron(14, 3), false),
  :properties => [1.0 0.49 0.0;],
  :g_coord => [
    0.0000  0.0000  0.0000;     0.0000  1.0000  0.0000;
    0.0000  0.0000 -1.0000;     0.0000  1.0000 -1.0000;
    1.0000  0.0000  0.0000;     1.0000  1.0000  0.0000;
    1.0000  0.0000 -1.0000;     1.0000  1.0000 -1.0000;
    0.1650  0.2980 -0.2550;     0.2720  0.7700 -0.2500;
    0.3200  0.3570 -0.8140;     0.2490  0.8080 -0.6580;
    0.7880  0.3560 -0.3070;     0.8500  0.7370 -0.3510;
    0.6770  0.3170 -0.6950;     0.8260  0.7120 -0.7120;
    0.0000  0.5000 -0.5000;     0.5000  0.5000  0.0000;
    0.5000  0.0000 -0.5000;     0.5000  1.0000 -0.5000;
    0.5000  0.5000 -1.0000;     1.0000  0.5000 -0.5000;
    0.2515  0.5582 -0.4942;     0.5188  0.5403 -0.2908;
    0.4875  0.3320 -0.5177;     0.5493  0.7568 -0.4928;
    0.5180  0.5485 -0.7197;     0.7853  0.5305 -0.5162;
    0.1093  0.5170 -0.1263;     0.5305  0.8767 -0.1503;
    0.9095  0.5232 -0.1645;     0.4882  0.1635 -0.1405;
    0.1213  0.1637 -0.5173;     0.1303  0.8945 -0.4770;
    0.9190  0.8623 -0.5158;     0.8663  0.1682 -0.5005;
    0.1423  0.5412 -0.8680;     0.5187  0.8800 -0.8425;
    0.8758  0.5072 -0.8518;     0.4992  0.1685 -0.8773
  ],
  :g_num => [
    11,  9, 13, 15, 25, 23, 24, 28, 27, 12, 10, 14, 16, 26,
     3,  1,  9, 11, 33, 17, 29, 23, 37,  4,  2, 10, 12, 34,
     9,  1,  5, 13, 32, 29, 18, 31, 24, 10,  2,  6, 14, 30,
    15, 13,  5,  7, 36, 28, 31, 22, 39, 16, 14,  6,  8, 35,
     3, 11, 15,  7, 40, 37, 27, 39, 21,  4, 12, 16,  8, 38,
    12, 10, 14, 16, 26, 34, 30, 35, 38,  4,  2,  6,  8, 20,
     3,  1,  5,  7, 19, 33, 32, 36, 40, 11,  9, 13, 15, 25
  ],
  :support => [
    (1, [0 1 1]), (2, [0 0 1]), ( 3, [0 1 0]), ( 4, [0 0 0]), ( 6, [1 0 1]),
    (7, [1 1 0]), (8, [1 0 0]), (17, [0 1 1]), (20, [1 0 1]), (21, [1 1 0])
  ],
  :fixed_freedoms => [
    ( 1, 2, 0.01),
    ( 3, 2, 0.01),
    ( 5, 2, 0.01),
    ( 7, 2, 0.01),
    (19, 2, 0.01)
  ]
)

data |> display
println()

@time m, dis_df, fm_df = p54(data)
println()

dis_df |> display
println()

println("\nThe integration point (nip = $(data[:struc_el].nip)) stresses are:\n")

fm_df |> display
println()


text:
function _firewall(input, delay::Int=0)
    n = size(input, 1)
    totalseverity = 0
    caught = 0
    for i in 1:n
        depth = input[i, 1]
        scanrange = input[i, 2]

        x = vcat(collect(1:scanrange), collect((scanrange-1):-1:2))
        x = circshift(x, depth + delay)

        if first(x)==1
            caught += 1
            totalseverity += depth*scanrange
        end

    end
    return (totalseverity, caught)
end

function puzzle13(path::String=joinpath(@__DIR__, "..", "data/13.txt"))
    input = readdlm(path, ':', Int)
    delay = 0
    notfound = true
    partone = nothing
    parttwo = nothing

    while notfound
        totalseverity, caught = _firewall(input, delay)

        if caught == 0
            notfound = false
            parttwo = delay
        end

        if delay == 0
            partone = totalseverity[1]
        end

        delay += 1
    end

    return [partone, parttwo]
end


text:
export NormalOGIP, logOGIP,
    PriceProcess, TestPriceProcess,
    SmallLeasedProblem,
    FullLeasedProblem,
    RoyaltyRates,
    RoyaltyRatesSmall,
    RoyaltyRatesFull,
    RoyaltyExogVars,
    RealDiscountRate,
    PsiSpace,
    TestPsiSpace,
    TestDynamicDrillModel,
    TestDataRoyalty,
    TestDataProduce,
    CoefLinks,
    ThetaRho, AlphaPsi, AlphaG, AlphaT, Alpha0,
    Theta,
    DefaultDrillReward,
    BaseYear

# ----------------------------
# Time
# ----------------------------

BaseYear() = ShaleDrillingLikelihood.TIME_TREND_BASE

# ----------------------------
# Geology
# ----------------------------

NormalOGIP(;mu=4.68, sigma=0.31) = Normal(mu,sigma)

logOGIP(num_i, d=NormalOGIP()) = rand(d, num_i)

# ----------------------------
# Price process
# ----------------------------

struct PriceProcess{TS<:TimeSeriesProcess, SRL<:Tuple, ETV<:ExogTimeVars, M<:AbstractMatrix}
    process::TS
    grid::SRL
    series::ETV
    transition::M

    function PriceProcess(process::TS, grid::SRL, series::ETV, transition::M) where {TS,SRL,ETV,M}
        all(minimum.(grid) .<= minimum(series)) || throw(error())
        all(maximum.(grid) .>= maximum(series)) || throw(error())
        length(product(grid...)) == checksquare(transition) || throw(DimensionMismatch())
        new{TS,SRL,ETV,M}(process,grid,series,transition)
    end
end

process(   x::PriceProcess) = x.process
grid(      x::PriceProcess) = x.grid
series(    x::PriceProcess) = x.series
transition(x::PriceProcess) = x.transition
zchars(    x::PriceProcess) = series(x)

function TestAR1Process(;rho=0.8, mean=1.33*(1-rho), var=0.265^2*(1-rho^2))
    return AR1process(mean, rho, var)
end

function TestPriceDates(n=150;start=Date(2003,10), step=Dates.Quarter(1))
    return range(start; step=step, length=n)
end

function TestPriceGrid(zvec, length=15; process=TestPriceProcess())
    zmin = min(minimum(zvec), lrmean(process)-3*lrstd(process))
    zmax = max(maximum(zvec), lrmean(process)+3*lrstd(process))
    zrng = range(zmin, zmax; length=length)
    return tuple(zrng)
end

function TestPriceTransition(grid::StepRangeLen; process=TestPriceProcess(), minp=1e-4)
    Pdense = tauchen_1d(process, grid)
    Pdense_zero = zero_out_small_probs(Pdense, minp)
    return sparse(Pdense_zero)
end

function TestPriceProcess(; process=TestAR1Process(), nsim=150, ngrid=13)
    dates = TestPriceDates(nsim)
    series = simulate(process, nsim)
    etv = ExogTimeVars(tuple.(series), dates)
    grid = TestPriceGrid(series, ngrid; process=process)
    transition = TestPriceTransition(grid...; process=process)
    return PriceProcess(process, grid, etv, transition)
end

# ----------------------------
# Lease problem
# ----------------------------

SmallLeasedProblem() = LeasedProblem(3, 4, 5, 3, 2)
FullLeasedProblem() = LeasedProblem(8, 8, 12, 12, 8)

# ----------------------------
# Royalty rates
# ----------------------------

RoyaltyRatesSmall() = [1/8, 3/16, 1/4,]
RoyaltyRatesFull() = [1/8, 1/6, 3/16, 1/5, 9/40, 1/4]

function RoyaltyRates(n)
    if n ∈ (3, :small)
        rr = RoyaltyRatesSmall()
    elseif n ∈ (6, :full)
        rr = RoyaltyRatesFull()
    else
        throw(error("don't know this n=$n"))
    end
    return rr::Vector{Float64}
end

function RoyaltyExogVars(num_i::Integer, rates, theta, log_ogip=zeros(num_i,0), model=RoyaltyModel())
    L = length(rates)
    k = length(theta) - (L-1) - extra_parm(model)
    @assert k > 0
    kminus = k - size(log_ogip',1)
    @assert kminus > 0
    Xminus = randn(kminus, num_i)
    X = vcat(log_ogip', Xminus)
    return X
end

# ----------------------------
# Test Data
# ----------------------------

function TestDataRoyalty(u, v, rates, theta, args...)
    num_i = length(u)
    num_i == length(v) || throw(DimensionMismatch())
    Xroyalty = RoyaltyExogVars(num_i, rates, theta, args...)
    data_roy = DataRoyalty(u, v, Xroyalty, theta, rates)
    return data_roy
end

function TestDataProduce(u, datadrill, obs_per_well, theta, log_ogip)
    nwells = total_wells_drilled(datadrill)
    return DataProduce(u, nwells, obs_per_well, theta, log_ogip)
end

# ----------------------------
# Model
# ----------------------------

RealDiscountRate() = ((0x1.006b55c832502p+0)^12 / 1.125) ^ (1/4)
PsiSpace(n; start=-4.5, stop=4.5) = range(start; stop=stop, length=n)
TestPsiSpace() = PsiSpace(13)

ichars(log_ogip, royalty_rates) = collect(zip(log_ogip, royalty_rates))

DefaultDrillReward() = DrillReward(
    DrillingRevenue(Unconstrained(), NoTrend(), NoTaxes() ),
    DrillingCost_constant(),
    ExtensionCost_Constant()
)

function TestDynamicDrillModel(z::PriceProcess, anticipate_t1ev::Bool = true;
    reward = DefaultDrillReward(),
    discount = RealDiscountRate(),
    statespace = FullLeasedProblem(),
    psispace = TestPsiSpace()
)

    zs = grid(z)
    ztrans = transition(z)
    ddm = DynamicDrillModel(
        reward, discount, statespace, zs, ztrans, psispace, anticipate_t1ev
    )
    return ddm
end


text:
using Gadfly, DataFrames

set_default_plot_size(6inch, 3inch)

df = DataFrame()
df[:pb]=rand(1:50, 50);
df[:t]=rand(1:50, 50);
df[:h]=[0.0 for i in 1:50];

plot(df, x = :pb, y = :t, color = :h, Geom.rectbin)