text:
# Make our document distances work with the TextAnalysis package

const DefaultDocDistance = WordMoversDistance()

# Find the closest token with an embedding in a DocDistance. If the
# token itself has no embedding we will try to "clean" it until we find
# an embedding which is likely to carry very similar meaning.
function findclosetoken(t::String, dd::AbstractDocumentDistance = DefaultDocDistance;
    skipdot = true, skiphyphen = true, skipcomma = true, makelowercase = true)
    hasembedding(dd, t) && return(t)
    if (skipdot && endswith(t, ".")) || (skiphyphen && endswith(t, "-")) ||
        (skipcomma && endswith(t, ","))
        t = t[1:prevind(t, length(t), 1)]
        hasembedding(dd, t) && return(t)
    end
    if makelowercase
        t = lowercase(t)
    end
    # Return nothing to indicate we couldn't find a close token with an embedding
    return hasembedding(dd, t) ? t : nothing
end

# Filter the tokens of a Document so we keep only the ones that have an embedding.
# We try some simple transformations before we give up on finding an embedding.
function filtertokens(tokens, dd::AbstractDocumentDistance = DefaultDocDistance;
    skipdot = true, skiphyphen = true, skipcomma = true, makelowercase = true)
    res = String[]
    for t in tokens
        t2 = findclosetoken(t, dd; skipdot = skipdot, skiphyphen = skiphyphen, skipcomma = skipcomma, makelowercase = makelowercase)
        if !isnothing(t2)
            push!(res, t2)
        end
    end
    return res
end

# Rather than filter tokens every time we want to calculate the doc distance
# we can pre-filter and cache the ngrams and then more quickly calculate
# the distance.
mutable struct FilteredNGramDocument <: AbstractDocument
    ngramdoc::NGramDocument # NGramDoc where only the filtered/transformed tokens are left
    metadata::TextAnalysis.DocumentMetadata
end

import TextAnalysis.ngrams
ngrams(d::FilteredNGramDocument) = ngrams(d.ngramdoc)

function FilteredNGramDocument(d::AbstractDocument, dd::AbstractDocumentDistance = DefaultDocDistance)
    ftoks = filtertokens(tokens(d), dd)
    ngd = NGramDocument(TextAnalysis.ngramize(d.metadata.language, ftoks, 1), 1)
    FilteredNGramDocument(ngd, d.metadata)
end

FilteredNGramDocument(s::AbstractString, dd::AbstractDocumentDistance = DefaultDocDistance) =
    FilteredNGramDocument(StringDocument(s), dd)

function FilteredNGramDocument(d::NGramDocument, dd::AbstractDocumentDistance = DefaultDocDistance)
    toks = keys(ngrams(d))
    ftoks = filtertokens(toks, dd)
    ngd = NGramDocument(ftoks, 1)
    FilteredNGramDocument(ngd, s.metadata)
end

# Quicker calculation of the marginals and vocabulary since we already have the counts.
function marginals_and_vocabulary(d1::FilteredNGramDocument, d2::FilteredNGramDocument; usesparse = true)
    v = make_vocabulary()
    k1 = keys(d1.ngramdoc.ngrams)
    k2 = keys(d2.ngramdoc.ngrams)
    update_vocabulary!(v, k1)
    update_vocabulary!(v, k2)
    m1 = calc_marginal(d1, k1, v; usesparse = usesparse)
    m2 = calc_marginal(d2, k2, v; usesparse = usesparse)
    return m1, m2, v
end

# Quicker calculation of the marginals since we already have the counts.
function calc_marginal(doc::FilteredNGramDocument, ks::Base.KeySet, vocabulary::AbstractDict{T, V}; usesparse = true) where {T<:AbstractString, V<:Integer}
    N, counts = makecountvector(vocabulary, usesparse)
    for k in ks
        if haskey(vocabulary, k)
            idx = vocabulary[k]
            counts[idx] += doc.ngramdoc[k]
        end
    end
    return frequencies_from_counts(counts)
end

# Find all pdf files (recursively) in a dir, convert them to text files, then read
# in a corpus representing them. Iff onlyifnew is true then convert to text file
# only if the pdf file is newer than the text file.
#function pdf_directory_corpus(dirname::String;
#    recursive = true, onlyifnew = true, wordDistance = WordDistanceCache(),
#    ngramdocs = true)
#    txtfiles = convert_pdf_files_to_txt(dirname; recursive = recursive, onlyifnew = onlyifnew)
#    docs = map(fp -> EmbeddedTokensFileDocument(fp, wordDistance), txtfiles)
#    if ngramdocs
#        docs = map(NGramDocument, docs)
#    end
#    return Corpus(docs)
#end


text:
using Distributions, ForwardDiff
using KernelGoodnessOfFit

# setup the problem
rbf_ = GaussianRBF(1.0)

μ = 0.0
σ² = 1.0
p_ = MultivariateNormal([μ, 0], [σ² 0; 0 σ²])  # true
q_ = MultivariateNormal([μ, 5.0], [σ² 0; 0 σ²])  # false

vs_ = [0.0 ; 0.0]
vs_ = reshape(vs_, 2, 1)

n_ = 100
xs_ = rand(p_, n_)

# simple example
T = KernelGoodnessOfFit.KSDTest(rbf_, 10, 0.05)

@test perform(T, q_, xs_).result == :reject
@test perform(T, p_, xs_).result == :accept

ext: jl


text:
## -------------------------------------------------------------------
function ToyModel2D()
    net = Dict()

    # feistGenomescaleMetabolicReconstruction2007 DOI: 10.1038/msb4100155
    # Using chemostat data for E. coli growing aerobically on
    # glucose (see Supplementary information), we estimated the
    # GAM and NGAM costs (Figure 3D). We found that an NGAM
    # value of 8.39 mmol ATP gDW1 h1 and a GAM value of
    # 59.81 mmol ATP gDW1 best fit the experimental data.
    GAM = 59.81 # EColi GAM 59.81 mmol ATP gDW1
    NGAM = 8.4 # EColi NGAM 8.4 mmolATP/gDW/h

    # https://en.wikipedia.org/wiki/Glycolysis
    # The overall process of glycolysis is:
    # Glucose + 2 NAD+ + 2 ADP + 2 Pi → 2 pyruvate + 2 NADH + 2 H+ + 2 ATP
    # Glucose + 2 NAD+ + 2 ADP + 2 Pi → 2 ac-Coa + 2 NADH + 2 H+ + 2 ATP
    GEY = 2.0 # (glycolisis energy yield) Per glucose molecule

    # https://en.wikipedia.org/wiki/Citric_acid_cycle
    # The theoretical maximum yield of ATP through oxidation of one molecule of glucose in
    # glycolysis, citric acid cycle, and oxidative phosphorylation is 38 (assuming 3 molar
    # equivalents of ATP per equivalent NADH and 2 ATP per UQH2).
    REY = (38.0 - GEY)/2 # (Respiratory energy yield) Per acetate molecule

    # http://ecolistudentportal.org/article_fermentation#_
    # Many fermenting bacteria including E. coli can convert acetyl-CoA to
    # acetate to generate ATP. The two key enzymes are phosphate acetyltransferase and
    # acetate kinase. This mini-pathway generates 1 mole of acetate and
    # 1 mole of ATP (by SLP) per acetyl-CoA consumed.
    # ferm: acetyl-CoA + Pi → acetyl-P + CoA
    # ferm: acetyl-P + ADP → acetate + ATP
    # ferm: AcCoa -> Ac + Atp
    # I will not include the ATP poduction

    # z    ass
    # Data From:
    # Varma, (1993): 2465–73. https://doi.org/10.1128/AEM.59.8.2465-2473.1993.
    # FIG 2
    # GBY = 0.068 # mmol/gDW
    GBY = 1 / 0.068 # gDW/ mmol

    # uptakes limits
    AB = 1000.0 # abs bound
    # This model is bounded by the maximum rates found for EColi.
    # Data From:
    # Varma, (1993): 2465–73. https://doi.org/10.1128/AEM.59.8.2465-2473.1993.
    # Extract max exchages from FIG 3 to form the maximum polytope
    Vg = 20      # (glucose uptake limit) mmol/gDW h
    Vz = 0.5     # D < 0.5 below the acetate switch [23] Nature, 528(7580):99–104, 2015.

    net[:S] =
    # rxns:  glyc  resp  ferm  vatp  atpm    ua   z     ug   # mets
    [       -1.0   0.0   0.0   0.0   0.0    0.0  -GBY  1.0;  #  Glc
             GEY   REY   0.0   0.0  -NGAM   0.0  -GAM  0.0;  #  Atp
             GEY  -1.0  -1.0   0.0   0.0    0.0   0.0  0.0;  #  AcCoa
             0.0   0.0   1.0   0.0   0.0    1.0   0.0  0.0;  #  Ac
             GEY   REY   0.0  -1.0   0.0    0.0   0.0  0.0;  #  AUXvatp
    ]

    net[:mets] = ["Glc", "Atp", "AcCoa", "Ac", "AUXvatp"]
    net[:b] =    [0.0  , 0.0  , 0.0    , 0.0 , 0.0      ] # const exchanges

    net[:rxns] = [ "glyc" , "resp" , "ferm" , "vatp", "atpm", "ua" ,  "z"  , "ug" ]
    net[:lb]   = [ 0.0    , 0.0    , 0.0    , 0.0   , 1.0   , -AB  ,  0.0  , 0.0  ]
    net[:ub]   = [ AB     , AB     , AB     , AB    , 1.0   , 0.0  ,  Vz   , Vg   ]
    net[:c]    = [ 0.0    , 0.0    , 0.0    , 0.0   , 0.0   , 0.0  ,  -1.0 , 0.0  ]

    return reducebox!(MetNet(;net...))
end


text:
@testset "setters" begin
    @testset "Mag" begin
        @test π <
              Float64(Arblib.set!(Mag(), π)) ==
              Float64(Arblib.set!(Arblib.radref(Arb()), π)) ==
              Float64(Arblib.set!(Arblib.mag_struct(), π)) <
              3.15
    end

    @testset "Arb" begin
        @test Arblib.set!(Arb(), BigFloat(1.0)) == one(Arb)
        @test Arblib.set!(Arblib.realref(Acb()), BigFloat(1.0)) == one(Arb)
        @test Arblib.equal(Arblib.set!(Arblib.arb_struct(), BigFloat(1.0)), one(Arb))
        @test Arblib.set!(Arb(), one(Mag)) == one(Arb)
        @test Arblib.set!(Arblib.realref(Acb()), one(Mag)) == one(Arb)
        @test Arblib.equal(Arblib.set!(Arblib.arb_struct(), one(Mag).mag), one(Arb))

        @test Arblib.set!(Arb(), 1 // 2) == one(Arb) / 2
        @test Arblib.set!(Arblib.realref(Acb()), 1 // 2) == one(Arb) / 2
        @test Arblib.equal(Arblib.set!(Arblib.arb_struct(), 1 // 2), one(Arb) / 2)

        @test Arb(3) < Arblib.set!(Arb(), π) < Arb(4)
        @test Arblib.equal(Arblib.set!(Arb(), π), Arblib.set!(Arblib.realref(Acb()), π))
        @test Arblib.equal(Arblib.set!(Arb(), π), Arblib.set!(Arblib.arb_struct(), π))

        @test Arb(2) < Arblib.set!(Arb(), ℯ) < Arb(3)
        @test Arblib.equal(Arblib.set!(Arb(), ℯ), Arblib.set!(Arblib.realref(Acb()), ℯ))
        @test Arblib.equal(Arblib.set!(Arb(), ℯ), Arblib.set!(Arblib.arb_struct(), ℯ))

        @test Arb(0) < Arblib.set!(Arb(), MathConstants.γ) < Arb(1)
        @test Arblib.equal(
            Arblib.set!(Arb(), MathConstants.γ),
            Arblib.set!(Arblib.realref(Acb()), MathConstants.γ),
        )
        @test Arblib.equal(
            Arblib.set!(Arb(), MathConstants.γ),
            Arblib.set!(Arblib.arb_struct(), MathConstants.γ),
        )

        @test Arb(0) < Arblib.set!(Arb(), MathConstants.catalan) < Arb(1)
        @test Arblib.equal(
            Arblib.set!(Arb(), MathConstants.catalan),
            Arblib.set!(Arblib.realref(Acb()), MathConstants.catalan),
        )
        @test Arblib.equal(
            Arblib.set!(Arb(), MathConstants.catalan),
            Arblib.set!(Arblib.arb_struct(), MathConstants.catalan),
        )

        @test Arblib.overlaps(
            Arb(BigFloat(MathConstants.φ)),
            Arblib.set!(Arb(), MathConstants.φ),
        )
        @test Arblib.equal(
            Arblib.set!(Arb(), MathConstants.φ),
            Arblib.set!(Arblib.realref(Acb()), MathConstants.φ),
        )
        @test Arblib.equal(
            Arblib.set!(Arb(), MathConstants.φ),
            Arblib.set!(Arblib.arb_struct(), MathConstants.φ),
        )
    end

    @testset "Acb" begin
        @test Arblib.set!(Acb(), BigFloat(1.0)) == one(Acb)
        @test Arblib.equal(Arblib.set!(Arblib.acb_struct(), BigFloat(1.0)), one(Acb))
        x = Arblib.set!(Acb(), ℯ)
        @test isequal(Arblib.realref(x), Arb(ℯ))
        @test Arblib.imagref(x) == zero(Arb)
        @test Arblib.equal(Arblib.set!(Arblib.acb_struct(), BigFloat(1.0)), one(Acb))
        @test Arblib.set!(Acb(), one(Arf).arf) == one(Acb)
        @test Arblib.equal(Arblib.set!(Arblib.acb_struct(), one(Arf).arf), one(Acb))
        @test Arblib.set!(Acb(0, 5), BigFloat(1.0)) == one(Acb)

        @test Arblib.set!(Acb(), BigFloat(1.0), 2) == Acb(1, 2)
        x = Arblib.set!(Acb(), π, Arf(4))
        @test isequal(Arblib.realref(x), Arb(π))
        @test isequal(Arblib.imagref(x), Arb(4))
        @test Arblib.set!(Acb(), Arb(1), Arf(2)) == Acb(1, 2)
        @test Arblib.set!(Acb(), Arb(1).arb, Arf(2).arf) == Acb(1, 2)

        @test Arblib.set!(Acb(), Complex(1, 2)) == Acb(1, 2)
        @test Arblib.set!(Acb(), Complex(BigFloat(1), BigFloat(2))) == Acb(1, 2)
        @test Arblib.equal(Arblib.set!(Arblib.acb_struct(), Complex(1, 2)), Acb(1, 2))

        x = Arblib.set!(Acb(), π)
        @test isequal(Arblib.realref(x), Arb(π))
        @test Arblib.imagref(x) == zero(Arb)
        x = Arblib.set!(Arblib.acb_struct(), π)
        @test Arblib.equal(Arblib.realref(x), Arb(π))
        @test Arblib.equal(Arblib.imagref(x), zero(Arb))
    end

    @testset "Set Rational" begin
        z = Arb(prec = 64)
        z[] = 1 // 2
        @test z == 0.5

        z = Acb(prec = 64)
        z[] = 1 // 2
        @test z == 0.5

        z = Acb(prec = 64)
        z[] = 1 // 2 + 3 // 4 * im
        @test z == 0.5 + 0.75im
    end
end


text:
module GMHPhotoReceptorAnalysis

import GeneralizedMetropolisHastings
import GMHPhotoReceptor
import Distributions
import PyPlot

import Base: show,mean,std,var

using GeneralizedMetropolisHastings
using GMHPhotoReceptor

include("analysis/StreamingStats.jl")
include("analysis/evaluate.jl")
include("analysis/analysis.jl")

end


text:
# Create example data
plot(-2π:0.01:2π, sin)
xlim(-2π, 2π)
ylim((-sqrt(2), sqrt(2)), true)
xticks(π/8, 4)
xticklabels(x -> @sprintf("%0.1f\\pi", float(x/π)))


text:
module test_h3_lib_miscellaneous

using H3.Lib

Lib.degsToRads
Lib.radsToDegs
Lib.hexAreaKm2
Lib.hexAreaM2
Lib.edgeLengthKm
Lib.edgeLengthM
Lib.numHexagons
Lib.getRes0Indexes
Lib.res0IndexCount

end # module test_h3_lib_miscellaneous


text:
using Ripserer
using BenchmarkTools

const SUITE = BenchmarkGroup()
const SKIP = []

for file in readdir(@__DIR__)
    if file in SKIP
        @warn "skipping $file"
    elseif startswith(file, "bench_") && endswith(file, ".jl")
        SUITE[file[(length("bench_") + 1):(end - length(".jl"))]] = include(file)
    end
end


text:
"""
    struct TargetOptions <: Any
Hold a pointer to a `clang::TargetOptions` object.
"""
struct TargetOptions
    ptr::CXTargetOptions
end

Base.unsafe_convert(::Type{CXTargetOptions}, x::TargetOptions) = x.ptr
Base.cconvert(::Type{CXTargetOptions}, x::TargetOptions) = x


text:
#testing the fused dot product

#test constructing a posit from the following tuple:

# a helper function.
# (sign::bool, exponent::Int64, fraction::UInt64)
# relevant formula:  value = ((-1) ^ sign) * (2 ^ exp) + (2^exp) * (fraction / 2^(64))
function float_4_component(sign::Bool, exp::Int64, frac::UInt64)
  ((-1)^sign) * (2.0^exp) + frac * (2.0^(exp - 64 - sign))
end
#test float_construction
@test float_4_component(false, 0, 0x0000_0000_0000_0000) == 1.0
@test float_4_component(false, 0, 0x8000_0000_0000_0000) == 1.5
@test float_4_component(false, 1, 0x0000_0000_0000_0000) == 2.0
@test float_4_component(false, 1, 0x8000_0000_0000_0000) == 3.0
@test float_4_component(false, -1, 0x0000_0000_0000_0000) == 0.5
@test float_4_component(false, -1, 0x8000_0000_0000_0000) == 0.75

@test float_4_component(true, 0, 0x0000_0000_0000_0000) == -1.0
@test float_4_component(true, 0, 0x8000_0000_0000_0000) == -0.75
@test float_4_component(true, 1, 0x0000_0000_0000_0000) == -2.0
@test float_4_component(true, 1, 0x8000_0000_0000_0000) == -1.5
@test float_4_component(true, -1, 0x0000_0000_0000_0000) == -0.5
@test float_4_component(true, -1, 0x8000_0000_0000_0000) == -0.375

#test a fixed point construction.
float_4_fixed(data::UInt8) = reinterpret(Int8, data) / 16.0

@test float_4_fixed(0b0001_0000) == 1.0
@test float_4_fixed(0b0001_1000) == 1.5
@test float_4_fixed(0b0010_0000) == 2.0
@test float_4_fixed(0b0011_0000) == 3.0
@test float_4_fixed(0b0000_1000) == 0.5
@test float_4_fixed(0b0000_1100) == 0.75

@test float_4_fixed(0b1111_0000) == -1.0
@test float_4_fixed(0b1111_0100) == -0.75
@test float_4_fixed(0b1110_0000) == -2.0
@test float_4_fixed(0b1110_1000) == -1.5
@test float_4_fixed(0b1111_1000) == -0.5
@test float_4_fixed(0b1111_1010) == -0.375

#test breaking up a fixed point construction as the descriptive components.
function component_4_fixed(data::UInt8)
  if (data & 0x80) != 0
    #then we're negative.  Meaure exponent by counting leading ones.
    lo = leading_ones(data)
    exponent = 4 - lo
    fraction = UInt64(data << (lo + 1)) << 56
    (true, exponent, fraction)
  else
    #else we're positive
    lz = leading_zeros(data)
    exponent = 3 - lz
    fraction = UInt64(data << (lz + 1)) << 56
    (false, exponent, fraction)
  end
end

#comprehensively test that fixed breakup correctly generates the desired value (except for zero)
for test_val = 0x01:0xff
  @test (float_4_component(component_4_fixed(test_val)...), test_val) == (float_4_fixed(test_val), test_val)
end

#comprehensively test that fixed breakup correctly generates the desired value (except for zero)
for test_val = 0x01:0xff
  @test (Posit{8,0}(component_4_fixed(test_val)...), test_val) == (Posit{8,0}(float_4_fixed(test_val)), test_val)
end

for test_val = 0x01:0xff
  @test (Posit{8,1}(component_4_fixed(test_val)...), test_val) == (Posit{8,1}(float_4_fixed(test_val)), test_val)
end

#test breaking up a posit into descriptive elements.
for test_val = 0x01:0xff
  if (test_val != 0x80)
    @test Posit{8,1}(SigmoidNumbers.posit_components(Posit{8,1}(test_val))...) == Posit{8,1}(test_val)
  end
end

#test loading a posit into a new quire and pulling it back out.
test_q = Quire(Posit{8,1})
for test_val = 0x00:0xff
  set!(test_q, Posit{8,1}(test_val))
  @test Posit{8,1}(test_val) == Posit{8,1}(test_q)
end

#test loading a posit into a new quire and pulling it back out in 32 bits
for idx = 1:10000
  test_val = rand(UInt32)
  set!(test_q, Posit{32,2}(test_val))
  @test Posit{32,2}(test_val) == Posit{32,2}(test_q)
end

################################################################################
# test adding two 32-bit posits.

for test_val1 = 0x00:0xff, test_val2 = 0x00:0xff

  (test_val1 == 0x80 && test_val2 == 0x80) && continue

  p1 = Posit{8,1}(test_val1)
  p2 = Posit{8,1}(test_val2)

  set!(test_q, p1)
  add!(test_q, p2)

  @test (Float64(Posit{32,2}(test_q)), test_val1, test_val2) == (Float64(p1) + Float64(p2), test_val1, test_val2)
end

p1 = Posit{8,1}(0x85)
p2 = Posit{8,1}(0x8b)

SigmoidNumbers.zero!(test_q)
fdp!(test_q, p1, p2)

@test (Float64(Posit{32,2}(test_q))) == Float64(p1) * Float64(p2)

################################################################################
# test multiplying two 32-bit posits.

for test_val1 = 0x00:0xff, test_val2 = 0x00:0xff

  (test_val1 == 0x80 || test_val2 == 0x80) && continue

  p1 = Posit{8,1}(test_val1)
  p2 = Posit{8,1}(test_val2)

  SigmoidNumbers.zero!(test_q)
  fdp!(test_q, p1, p2)

  @test (Float64(Posit{32,2}(test_q)), test_val1, test_val2) == (Float64(p1) * Float64(p2), test_val1, test_val2)
end

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
function strip_escape!(cur_expr::Expr)
  ( cur_expr.head == :escape ) || return cur_expr
  ( length(cur_expr.args) == 1 ) || return cur_expr

  first_expr = first(cur_expr.args)

  cur_expr.head = first_expr.head
  cur_expr.args = first_expr.args

  cur_expr
end


text:
"""
    emptyshim(ContainerType)
    emptyshim(ContainerType, ElementType)

Create an empty "shim" container that is widen to `ContainerType` when
appending elements to it.  Default `ElementType` is `Union{}`.

# Examples
```jldoctest
julia> using MicroCollections, BangBang

julia> @assert append!!(emptyshim(Vector), [0]) == [0]

julia> @assert merge!!(emptyshim(Dict), Dict(:a => 1)) == Dict(:a => 1)

julia> @assert union!!(emptyshim(Set), Set([0])) == Set([0])
```
"""
emptyshim

"""
    singletonshim(ContainerType, x)

Create a "shim" container with one element `x` that is widen to
`ContainerType` when appending elements to it.

# Examples
```jldoctest
julia> using MicroCollections, BangBang

julia> @assert push!!(singletonshim(BitVector, false), true)::BitArray == [false, true]
```
"""
singletonshim

emptyshim(T) = emptyshim(T, Union{})

"""
    vec0(T::Type = Union{})

Create an empty vector shim of element type `T`.
"""
vec0

"""
    vec1(x)

Create a singleton vector with element `x`.
"""
vec1

"""
    dict0(T::Type = Pair{Union{},Union{}})

Create an empty dict shim of element type `T`.
"""
dict0

"""
    dict1(k => v)

Create a singleton dict with key `k` and value `v`.
"""
dict1

"""
    set0(T::Type = Union{})

Create an empty set shim of element type `T`.
"""
set0

"""
    set1(x)

Create a singleton set with element `x`.
"""
set1

vec0() = emptyshim(Vector, Union{})
vec0(::Type{T}) where {T} = emptyshim(Vector, T)
vec1(x) = singletonshim(Vector, x)

dict0() = emptyshim(Dict)
dict0(::Type{T}) where {T} = emptyshim(Dict, T)
dict1(x) = singletonshim(Dict, x)

set0() = emptyshim(Set)
set0(::Type{T}) where {T} = emptyshim(Set, T)
set1(x) = singletonshim(Set, x)

"""
    MicroCollections.upcast(singleton_container::C) -> container::C
"""
upcast

# utils
function validate_empty_iterator(itr, container_name)
    y = iterate(itr)
    y === nothing || throw(ArgumentError(string(
        "Given iterator is not empty. ",
        "`$container_name` expects an empty iterator.",
    )))
end

function validate_singleton_iterator(itr, container_name)
    y = iterate(itr)
    y === nothing && throw(ArgumentError(string(
        "Given iterator is empty. ",
        "`$container_name` requires an iterator of length 1.",
    )))
    x, state = y
    y = iterate(itr, state)
    y === nothing || throw(ArgumentError(string(
        "Given iterator has more than one element. ",
        "`$container_name` requires an iterator of length 1.",
    )))
    return x
end

function safe_length(x)
    if Base.IteratorSize(x) isa Union{Base.HasLength,Base.HasShape}
        return length(x)
    else
        return nothing
    end
end


text:
# lgr_load.jl
#
#   Load Los Gatos Research (LGR) Methane data
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.4.2
# 07.11.2014
# Last Edit: 17.07.2020

# - Programmatically zipped data files have a PGP signature at the end after the last line of data
# - Data files are TXT files withing a ZIP file
# - Sometimes multiple text files and a folder can be contained within a ZIP file
# - Line width is assumed to be fixed

"# lgr_load(source::String,mindate::DateTime,maxdate::DateTime;verbose::Bool=false)

`data = lgr_load(source)` Load a directory of LGR zip files\n
* **data**::DataFrame = Data
* **source**::String = Directory of LGR zip files\n\n

`data = lgr_load(source,minimumdate)` \n
* **minimumdate**::DateTime = Starting date to load data from\n\n

`data = lgr_load(source,minimumdate,maximumdate)` \n
* **maximumdate**::DateTime = Load data before this date\n\n

---\n

#### Keywords:\n
* verbose::Bool = Display information as the function runs, TRUE is default\n\n"
function lgr_load(source::String,mindate::DateTime=DateTime(0),maxdate::DateTime=DateTime(9999);verbose::Bool=false)
	# Check source
	if !isdir(source) & !isfile(source)
		error("source is not a file or directory.")
	end

	# Temporary Directory, Unzipping Destination
	dest = tempdir()

	#############################
	##  Basic Settings Output  ##
	#############################
	if verbose
		println("Loading LGR ZIP Files")
		println("  From: " * string(mindate))
		println("  To  : " * string(maxdate))
		println("  Temp Directory: " * dest)
	end

	##################
	##  List Files  ##
	##################
	(files,folders) = dirlist(source,regex=r"\d{2}[A-Z][a-z][a-z]\d{4}\.zip$")
	times = Array{DateTime}(undef,length(files))
	for i=1:1:length(times)
		times[i] = DateTime(files[i][end-12:end-4],"dduuuyyyy") # Convert file name to a timestamp
	end
	f = sortperm(times)
	times = times[f]
	files = files[f]

	# Remove files out of time range
	f = findall(mindate .<= times .< maxdate)
	times = times[f]
	files = files[f]
	verbose ? println("  ZIP File Count: " * string(length(files)) * "\n") : nothing

	# List and unzip contents
	D = DataFrame()
	for i=1:1:length(files)
		println("-=-==  " * files[i] * "  ==-=-") # Temp
		temp = files[i]
		zipfiles = Vector{String}()
		zipdirectories = Vector{String}()

		# List Zip Contents
		verbose ? println("   # " * string(i) * ": " * basename(temp)) : nothing
		zipfiles = joinpath.(dest,zipList(files[i])) # List the files in the archive
		zipExtractAll(files[i],dest)

		# Load Each Data File
		for j=1:1:length(zipfiles)
			tempD = lgr_read(zipfiles[j],verbose=verbose)
			isempty(D) ? D = tempD : D = [D;tempD]
		end

		# Delete Files and Directories
		for j=1:1:length(zipfiles)
			if isfile(zipfiles[j])
				rm(zipfiles[j])
			end
		end
		for j=1:1:length(zipdirectories)
			if isdir(zipdirectories[j])
				rm(zipdirectories[j])
			end
		end
	end

	verbose ? println("  Data Loading Complete") : nothing

	return D
end

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO


text:
module Discourse

using HTTP, Dates, JSON, CodecZlib, Base64

# HTTP verbs to be used as non string literals
DELETE = "DELETE"
GET = "GET"
POST = "POST"
PUT = "PUT"

mutable struct DiscourseClient
    host::String
    api_username::String
    api_key::String
    timeout::Second
    function DiscourseClient(args...)
        new(host, api_username, api_key, timeout)
    end
end

include("utils/utils.jl")
export
        convert_HTTP_Response_To_JSON
        _get
        _request
include("Categories/Categories.jl")
include("Users/User.jl")
export
        get_user

end # module


text:
function homenowallmind!(r)
    x=0
    y=0
    while isborder(r,Nord)==0
        move!(r,Nord)
        x=x+1
    end
    while isborder(r,West)==0
        move!(r,West)
        y=y+1
    end
    return x,y
end

function homenowall!(r,x,y)
    for i in 1:x
        move!(r,Sud)
    end
    for i in 1:y
        move!(r,Ost)
    end
end
function cornernowall!(r)
    while isborder(r,Nord)==0
        move!(r,Nord)
    end
    while isborder(r,West)==0
        move!(r,West)
    end
end

function counter!(r)
c=1
    while isborder(r,Ost)==0
        move!(r,Ost)
        c=c+1
    end
    return c
end

function movechessW!(r)
    putmarker!(r)
    if isborder(r,West)==0
        move!(r,West)
    end
    if isborder(r,West)==0
        move!(r,West)
    end
end

function chess(r,y,c)

    cc=c
    while isborder(r,Sud)!=1
        if ((c-y)%2==0)
            while cc!=0
                movechessW!(r)
                cc=cc-2
            end
        else
            if isborder(r,West)==0
                move!(r,West)
            end
            while cc!=0
                movechessW!(r)
                cc=cc-2
            end
        end

        if isborder(r,Sud)==0
            move!(r,Sud)
        end
        while isborder(r,Ost)==0
            move!(r,Ost)
        end
        move!(r,Nord)
        cc=c
        if ismarker(r)==0
            move!(r,Sud)
            while cc!=0
                movechessW!(r)
                cc=cc-2
            end
        else
            move!(r,Sud)
            move!(r,West)
            cc=c
            while cc!=0
                movechessW!(r)
                cc=cc-2
            end
        end
    end
end
function task_7(r)
    x,y=homenowallmind!(r)
    c=counter!(r)
    chess(r,y,c)
    cornernowall!(r)
    homenowall!(r,x,y)
end

task_7(r)


text:
module Optimise

export train!, update!,
	SGD, Descent, ADAM, Momentum, Nesterov, RMSProp,
	ADAGrad, AdaMax, ADADelta, AMSGrad, NADAM, ADAMW,RADAM,
	InvDecay, ExpDecay, WeightDecay, stop, Optimiser, Rprop, SparseADAM,
	HeavyBall

include("optimisers.jl")
include("train.jl")

end


text:
## RayMain.jl
# Todd Anderson
# 5 October 2020
#
# Ray-tracing of electromagnetic waves from upper ionosphere to plasmasphere using Haselgrove equations, following method in Rice 1997.  Includes main ODE solver function calls.
#
# Dependencies:
#   Haselgrove.jl:      Ray-tracing equations derived by Haselgrove (1957)
#   RTConstants.jl:     Global constants
#   RTParameters.jl:    Parameter choices for ODE solver
#   RTPlot.jl:          Plotting function calls and options

#= Outline of Python implementation:
    1. Import packages
    2. Read input parameters and initial conditions
        includes parameters in ray_param.dat
        a. load parameters
        b. load constants
        c. bundle initial conditions
        d. select and initialize ionosphere
    3. Calculate RHS of ray-tracing equations
        a. f(t, rtcG, freq)
            rtcG: current values of r, θ, χ, G
        b. calculate μ with call of phase_refractive_index(r,θ,χ,freq,ion[])
        c. calculate dμ/dx, where x = r, θ, χ, freq
            i) note that dμ/dχ = -1*(dμ/dψ)
        d. calculate derivatives w.r.t:
            derivs = [dr/dt, dθ/dt, dχ/dt, dfreq/dt]
    4. ODE solver calls
        a. select solver (e.g. LSODA)
        b. set plot parameters -- see ray_plot.py
        c. loop over initial frequencies to calculate ray paths:
            i.  determine number of evenly-spaced frequencies
            ii. initialize freq, group delay time, nphase arrays
            iii.for loop:
                1. set frequency to new value in range
                2. call ode solver: soln = solve f (see 3a)
                    a. set integration method
                    b. set initial values of radius, latitude, gdt, nphase
                3. while: solver is successful & time < tstop & 1Re < altitude < max alt
                    a.  soln.integrate(soln.t+dt)
                    b.  radius.append(soln.y[0])
                    c.  latitude.apppend(soln.y[1])
                    d.  gdt.append(soln.y[3])
                    e.  nphase_single =
                        phase_refractive_index(radius,latitude,chi,freq,[ion])  [from soln]
                    f.  nphase.append(nphase_single[0]) [= μ; nphase_single[1] is dμ/dψ]
                4. calculate x = radius[:]*cos(latitude[:]); y = radius[:]*sin(latitude[:])
                5. append freq, gdt, nphase to arrays
                    QUESTION: this program uses the initialization of empty (unsized) arrays/lists and subsequent appending to fill them in, resulting in unknown list size and differing size between list elements (i.e. different frequency rays).  How slow is this compared to a pre-allocated implementation?
                    TODO: optimize this ODE solve process.
                6.  plot x,y for this frequency

    5. Show plots

=#
## Read parameters, constants, and initial conditions

## Ray-tracing equations calculation

#= Functions:
	phase_refractive_index!():
		inputs:	μ, dμdψ, r, θ, χ, freq
		outputs: mutates μ, dμdψ

	ddr!():
		inputs:	μ, dμdψ, r, θ, χ, freq, dμdr
		outputs: mutates dμdr
		calls: phase_refractive_index!()

	ddθ!():
		inputs:	μ, dμdψ, r, θ, χ, freq, dμdθ
		outputs: mutates dμdθ
		calls: phase_refractive_index!()

	ddχ!():
		inputs:	μ, dμdψ, r, θ, χ, freq, dμdχ
		outputs: mutates dμdχ
		calls: phase_refractive_index!()

	ddf!():
		inputs:	μ, dμdψ, r, θ, χ, freq, dμdf
		outputs: mutates dμdf
		calls: phase_refractive_index!()

	haselgrove!():
		inputs: μ, dμdψ, r, θ, χ, freq, t, dμ, ddt
				dμdr,dμdθ,dμdχ,dμdf = dμ
				ddt = [drdt, dθdt, dχdt]
		outputs: mutates ddt

=#

#=TODO:
	0. Rewrite particle model for realistic plasmasphere! See Sousa thesis p. 39 for plots of electron density models.
		a. see GCPM at plasmasphere.nasa.gov/models/
		b. the simplified GCPM looks like it could be approximated analytically as an exponential(?) decay multiplied by a dipole field equation
		i.e. n_e = a*exp(-b*r)*c*sqrt(1 + 3cos²θ)
	1. Get functions and argument-passing working with single input (i.e. non-vectorized, non-looping)
	2. Vectorize!
	3. Run and time solver
	4. Switch to in-place functions where possible
	5. Run and time solver, optimize
=#
#const rₑ = 6378100
c = 2.99792458e8
re = 6.3712e6      	# radius of Earth in meters
B0 = 3.0696381e-5   # magnitude of B field at Earth's equator surface in Tesla
e = 1.602e-19 		# electric charge in Coulombs
me = 9.1093e-31		# electron rest mass
mp = 1.6726219e-27	# proton rest mass
eps = 8.854e-12 	# permittivity of free space in mks

using DifferentialEquations

# calculate phase refractive index μ
# DONE! todo: convert this function to mutating form
# 	QUESTION: is this possible, given that phase_refractive_index!() must be called multiple times in a single ddr!() (e.g.) call, before using the result?
# 	ANSWER: yes! see call/assignment in ddr!()

u = function phase_refractive_index(r,θ,χ,freq)
# u = [μ, dμdψ]
    # convert from radial angle to wave normal angle
    dip = atan(2.0*cot(θ)) # dip angle: angle between horizontal and B field;

#	ψ = 2*π - (3π/2 - dip - χ) 		# this should be ≡ to the below, but the plotted ray paths are not all that similar!
	ψ = π/2.0 + dip + χ
#	ψ = χ - (3.0/2.0)*pi + dip		# wave normal angle: angle between wave direction and B; 11/9: fixed error (χ - ϕ) ➡ (ϕ - χ)
	# NOTE 11/24/2020: Rice 1997 uses ψ = χ - 3π/2 + dip, as well as a strange method of calculating χ (see FORTRAN code on p. 112)

    # convert from frequency to angular frequency
    ω = 2.0*pi*freq

    # convert radius to fraction of Earth radius
    rE = r/re

    # find magnetic field at (r,θ) from dipole field model
    # Bmag = B0*((re^3)/(r^3))*sqrt(4.0-3.0*cos((pi/2.)-θ)*cos((pi/2.0)-θ))
	# 	   = B0*((re^3)/(r^3))*sqrt(4 - 3sin²θ)
	# 	   ≡ B0*((re^3)/(r^3))*sqrt(1 + 3cos²θ)
	# dipole field magnitude from Parks p. 61:
	# B(r,λ) = μ₀M/(4πr³)*(1+3sin²λ)^(1/2)
	# 		 = μ₀M/(4πr³)*(1+3sin²(π/2 - θ))^(1/2)
	# 		 = μ₀M/(4πr³)*(1+3cos²θ)^(1/2)
	# here B0*re³ = μ₀M/(4π)
	Bmag = B0*(re^3/(r^3))*sqrt(1+3*cos(θ)*cos(θ))

    # calculate electron and proton density profiles
    n_e = 1.e6*(1.8e5*exp(-4.183119*(rE-1.0471)))
	n_p = 1.e6*(1.8e5*exp(-4.183119*(rE-1.0471)))

	# electron plasma angular frequency squared
    ω_e2 = (n_e*(e^2))/(eps*me)
	# proton plasma angular frequency squared
    ω_p2 = (n_p*(e^2))/(eps*mp)

	# electron angular gyrofrequency
    Ω_e = (e*Bmag)/(me)
	# proton angular gyrofrequency
    Ω_p = (e*Bmag)/(mp)

	# define R, L, P, D, S,
	# R ≡ 1 - Σ(ωₖ²/ω²)(ω/(ω+ϵₖΩₖ))
	R = 1.0 - (ω_e2/ω^2.0)*(ω/(ω - Ω_e)) - (ω_p2/ω^2.0)*(ω/(ω + Ω_p))
	#println("R = ",R)

	# L ≡ 1 - Σ(ωₖ²/ω²)(ω/(ω-ϵₖΩₖ))
	L = 1.0 - (ω_e2/ω^2.0)*(ω/(ω + Ω_e)) - (ω_p2/ω^2.0)*(ω/(ω - Ω_p))
	#println("L = ",L)

	# P ≡ 1 - Σ(ωₖ²/ω²)
	P = 1.0 - (ω_e2/ω^2.0) - (ω_p2/ω^2.0)
	#println("P = ",P)

	# D ≡ ½(R-L); S ≡ ½(R+L)
	D = (R - L)/2.0
	S = (R + L)/2.0

	# define parts of dispersion relation
	# Aμ⁴ - Bμ² + C = 0
	# A = Ssin²ψ + Pcos²ψ
	A = S*sin(ψ)^2.0 + P*cos(ψ)^2.0

	# B = RLsin²ψ + PS(1 + cos²ψ)
	B = R*L*sin(ψ)^2.0 + P*S*(1.0+cos(ψ)^2.0)

	# C  = PRL
	C = P*R*L

	# solve the dispersion relation for μ:
	# μ² = (B +- F)/2A
	# where F² = (RL-PS)²sin⁴ψ + 4P²D²cos²ψ
	F2 = (R*L - P*S)^2.0*sin(ψ)^4.0 + 4.0*(P*D*cos(ψ))^2.0
	F = sqrt(F2)

	# Rice 1997: typical solution to dispersion relation with quadratic formula
	μ2_minus = (B - F)/(2.0*A)
	μ2_plus = (B + F)/(2.0*A)

	# Bortnik 2004:
	# if B > 0
	# 	μ2_minus = (B - F)/(2.0*A)
	# else
	# 	μ2_plus = (2.0*C)/(B + F)
	# end

	μ_minus = try
		sqrt(abs(μ2_minus))	# abs() is not physical! for test only

	catch err
		if isa(err,DomainError)
			println("DomainError in μ_minus!")
			println("B=", B)
			println("F=", F)
			println("A=", A)
		else
			# println(err)
		end
	end

	μ_plus = try
		sqrt(abs(μ2_plus)) # abs() is not physical! for test only

	catch err
		if isa(err,DomainError)
			# println("DomainError in μ_plus!")
			# println("B=", B)
			# println("F=", F)
			# println("A=", A)
		else
			# println(err)
		end
	end

	# Electron whistler case: ψ = 0, μ² = R; this is the μ_plus case
	#μ = μ_minus		# EMIC case
	μ = μ_plus		# electron whistler case

	# Find dA/dψ, dB/dψ, dC/dψ, dμ/dψ
	dAdψ = 2.0*(S-P)*sin(ψ)*cos(ψ)
	dBdψ = 2.0*(R*L-P*S)*sin(ψ)*cos(ψ)
    dCdψ = 0.0
    #dμdψ = ((μ^4.0)*dAdψ-(μ^2.0)*dBdψ+dCdψ)/(4.0*A*(μ^3.0)-2.0*B*μ)

	dFdψ = 1/(2.0*F)*((R*L-P*S)^2 * 4*sin(ψ)^3*cos(ψ) - 8*(P*D)^2*sin(ψ)*cos(ψ))
	#dFdψ = sqrt(abs((R*L-P*S)^2 * 4*sin(ψ)^3*cos(ψ) - 8*(P*D)^2*sin(ψ)*cos(ψ)))
	dμdψ = 1/(2.0*μ)*((dBdψ + dFdψ)/(2*A) - 2*dAdψ*(B + F)/(2*A^2))
	# NOTE 11/24/2020: from Rice 1997; choose '+' for (B +- F) and (dBdψ +- dFdψ)

	#DEBUG check values
	#println("μ = ",μ)
	#println("dμdψ = ", dμdψ)

	# non-mutating output
	u = [μ,dμdψ]

end

# Derivatives of phase refractive index w.r.t. r, θ, χ,
# calculate derivative of μ w.r.t. r
dμdr = function ddr(r,θ,χ,freq)
	#μ,dμdψ = u			# unpack
	dr = 1.0e-11		# r step size

	μ_l = phase_refractive_index(r-dr/2.0,θ,χ,freq)
	#μ_l = μ
	μ_r = phase_refractive_index(r+dr/2.0,θ,χ,freq)
	#μ_r = μ

	dμdr = (μ_r[1] - μ_l[1])/dr
end

dμdθ = function ddθ(r,θ,χ,freq)
	#μ,dμdψ = u			# unpack
	dθ = 1.0e-11		# θ step size

	μ_l = phase_refractive_index(r,θ-dθ/2.0,χ,freq)
	#μ_l = μ
	#println("μ_l = ", μ)
	μ_r = phase_refractive_index(r,θ+dθ/2.0,χ,freq)
	#μ_r = μ
	#println("μ_r = ", μ)

	dμdθ = (μ_r[1] - μ_l[1])/dθ
end

dμdχ = function ddχ(r,θ,χ,freq)
	#μ,dμdψ = u			# unpack
	dχ = 1.0e-11		# χ step size

	μ_l = phase_refractive_index(r,θ,χ-dχ/2.0,freq)
	#μ_l = μ
	μ_r = phase_refractive_index(r,θ,χ+dχ/2.0,freq) # NOTE Nov 3: corrected error 'χ-dχ/2.0'➡'χ+dχ/2.0'
	#μ_r = μ

	dμdχ = (μ_r[1] - μ_l[1])/dχ
end

dμdf = function ddf(r,θ,χ,freq)
	#μ,dμdψ = u			# unpack
	df = 1.0e-5		# f step size

	μ_l = phase_refractive_index(r,θ,χ,freq-df/2.0)
	#μ_l = μ
	μ_r = phase_refractive_index(r,θ,χ,freq+df/2.0)
	#μ_r = μ

	dμdf = (μ_r[1] - μ_l[1])/df
end

# calculate derivatives w.r.t. time

# calculate derivatives w.r.t. time
function haselgrove!(du,u,p,t)
# 	du = drdt, dθdt, dχdt, dfdt
# 	u = r, θ, χ, f
#  	p = freq, μ, dμdψ, dμdr, dμdθ, dμdχ, dμdf
# 		^ this compiles and runs, but parameters
# 	t = time

	r, θ, χ, freq = u 						# unpack u
	#freq = p[1] 	# unpack p

	dμdr = ddr(r,θ,χ,freq)			# mutating form: ddr!([μ,dμdψ],r,θ,χ,freq,dμdr)
	dμdθ = ddθ(r,θ,χ,freq)			# calls phase_refractive_index(r,θ,χ,freq)
	dμdχ = ddχ(r,θ,χ,freq)
	dμdf = ddf(r,θ,χ,freq)

	v = phase_refractive_index(r,θ,χ,freq)
	μ = v[1]
	dμdψ = v[2]

	# 11/20: try dμdχ --> -1*dμdψ; these should be equal but are calcualted differently
    du[1] = 1/(μ^2)*(μ*cos(χ)-dμdψ*sin(χ))
    du[2] = 1/(r*μ^2)*(μ*sin(χ)+dμdψ*cos(χ))
    du[3] = 1/(r*μ^2)*(dμdθ*cos(χ)-(r*dμdr + μ)*sin(χ))
	du[4] = 1/c*(1+(freq/μ)*dμdf)

	#du = [drdt, dθdt, dχdt, dfdt]

	#println("drdt = ", drdt)
	#println("dθdt = ", dθdt)
	#println("dχdt = ", dχdt)

	#DEBUG check values
	#println("drdt = ",drdt)
	#println("dθdt = ",dθdt)
	#println("dχdt = ",dχdt)

    # dμdψ = 1/(2*μ)*((dBdψ + dFdψ)/(2*A) - 2*dAdψ*(B+F)/(2*A^2))
    # dμdω = 1/(2*μ)*((dBdω + dFdω)/(2*A) - 2*dAdω*(B+F)/(2*A^2))
	#
    # dμdr = (μ*(r+δr/2) - μ*(r-δr/2))/δr
    # dμdθ = (μ*(θ+δθ/2) - μ*(θ-δθ/2))/δθ

end

#DEBUG check values
#haselgrove!(1.0,1.0,10000.0,pi/4.0,0.0,5000.0,0.0,[1.0,1.0,1.0,1.0],[1.0,1.0,1.0])

## ODE solver calls

## Plot results
using BenchmarkTools
using LSODA
using Sundials

u0 = [re+1.0e+6, 1.0*pi/4, 0.0, 5000.0]					# r0, θ0, χ0
p = []	# f0, dμdψ, dμdr, dμdθ, dμdχ, dμdf
tspan = (0.0,5.0e+9)

hasel_prob = ODEProblem(haselgrove!,u0,tspan,p)
hasel_soln = solve(hasel_prob, CVODE_BDF(), reltol=1e-7)

using Plots
plotly()
plot(hasel_soln)
plot(hasel_soln, vars=(1,2,3))

t = hasel_soln.t
r = hasel_soln[1,:]
θ = hasel_soln[2,:]
χ = hasel_soln[3,:]
f = hasel_soln[4,:]

x = r.*sin.(θ)
y = r.*cos.(θ)

plot(re.*cos.([0:0.01:2*π;]),re.*sin.([0:0.01:2*π;]), aspect_ratio = 1, legend=:none)
plot!(x,y, aspect_ratio=1)

max_stars_repo_name: andersontodds/ray-trace

max_issues_repo_name: andersontodds/ray-trace

max_forks_repo_name: andersontodds/ray-trace


text:
export
  putstreamfunc,
  putcallbackfunc

function msk_stream_callback_wrapper(userdata::Ptr{Void}, msg :: Ptr{UInt8})
  f = unsafe_pointer_to_objref(userdata) :: Function
  f(unsafe_string(msg))
  convert(Int32,0)::Int32
end

function putstreamfunc(t::MSKtask, whichstream:: Streamtype, f :: Function)
  cbfunc = cfunction(msk_stream_callback_wrapper, Int32, (Ptr{Void},Ptr{UInt8}))
  r = @msk_ccall(linkfunctotaskstream,Int32,(Ptr{Void},Int32, Any, Ptr{Void}),t.task,whichstream.value,f,cbfunc)
  if r != MSK_RES_OK.value
    throw(MosekError(r,getlasterror(t)))
  end
  t.streamcallbackfunc = cbfunc
  t.userstreamcallbackfunc = f
  nothing
end

# MSKcallbackfunc :: ( MSKtask_t, void *, Cint, double *, int32 *, int64 * -> int32 )
# The callback function is called with following parameers:
#  - a 'where' identifying where in the solver we currently are
#  - a 'dinf' array carrying double information on the current state of the solution/solver
#  - a 'iinf' array carrying int32 information on the current state of the solution/solver
#  - a 'linf' array carrying int64 information on the current state of the solution/solver
function msk_info_callback_wrapper(t::Ptr{Void}, userdata::Ptr{Void}, where :: Int32, douinf :: Ptr{Float64}, intinf :: Ptr{Int32}, lintinf :: Ptr{Int64})
    r = if unsafe_load(cglobal(:jl_signal_pending, Cint),1) > Cint(0) # Not so good... jl_signal_pending is actually a sig_atomic_t...
            1
        else
            task = unsafe_pointer_to_objref(userdata) :: MSKtask
            if task.usercallbackfunc == nothing
                0
            else
                dinfa  = unsafe_wrap(Vector{Float64},douinf,(length(Dinfitem),),false)
                iinfa  = unsafe_wrap(Vector{Int32},intinf,(length(Iinfitem),),false)
                liinfa = unsafe_wrap(Vector{Int64},lintinf,(length(Liinfitem),),false)

                task.usercallbackfunc(Callbackcode(where), dinfa, iinfa, liinfa)
            end
        end

    convert(Cint,r)::Cint
end

# f :: where :: Cint, dinf :: Array{Float64,1}, iinf :: Array{Int32,1}, linf :: Array{Int64,1} -> Int32
# NOTE: On Win32 the callback function should be stdcall
function putcallbackfunc(t::MSKtask, f::Function)
    t.usercallbackfunc = f
    nothing
end

function clearcallbackfunc(t::MSKtask)
    t.usercallbackfunc = nothing
    nothing
end


text:
# Marek Lipert 2015

type Point
  x::Float64
  y::Float64
  visited::Int64
  cluster::Int64
end

type Rectangle
  parent::Int64
  x0::Float64
  x1::Float64
  y0::Float64
  y1::Float64
  children::Array{Int64}
  kind::Int64 # 1 - leaf, 0 - non-leaf
end

type Cluster
  points::Array{Int64}
end

distance(x::(Float64,Float64),y::(Float64,Float64)) = sqrt((x[1]-y[1])^2 + (x[2]-y[2])^2)
distance(x,y) = sqrt((x.x-y.x)^2 + (x.y-y.y)^2)
distance(x,y::(Float64,Float64)) = sqrt((x.x-y[1])^2 + (x.y-y[2])^2)

distance(x::(Float64,Float64),y) = sqrt((x[1]-y.x)^2 + (x[2]-y.y)^2)

mean(sequence) = reduce(+,sequence)/length(sequence)
sigma(sequence,mean) = sqrt(reduce(+,map((x) -> (x-mean)^2,sequence))/(length(sequence)-1))

# helpers

# solves quadratic equation for two intersecting circles

function compute_possible_points(r1,r2,p1,p2)
  x1 = p1[1]
  x2 = p2[1]
  y1 = p1[2]
  y2 = p2[2]

  # safety (epsilon cannot be too small for numerical stability)
  epsilon = 0.001
  if abs(x1-x2) < epsilon
    x1 = x1 + epsilon
  end
  if abs(y1-y2) < epsilon
    y1 = y1 + epsilon
  end
  @assert distance(p1,p2) < r1 + r2 "Circles do not cross"

  a = (x2-x1)/(y2-y1)
  A = 1 + a^2

  S = ((r1^2 - r2^2) - ((x1^2-x2^2) + (y1^2-y2^2)))/(2*(y2-y1))
  B = 2*a*(y1 - S) - 2*x1
  C = (S - y1)^2 + x1^2 - r1^2

  delta = B^2 - 4 * A * C
  @assert delta > 0

  x_temp = (-B - sqrt(delta))/(2*A)
  sol_1 = (x_temp , S - a * x_temp)

  x_temp = (-B + sqrt(delta))/(2*A)
  sol_2 = (x_temp , S - a * x_temp)
  (sol_1, sol_2)
end

# checks if a circle of radius eps and center at point intersects the rectangle

function intersects(rectangle,point,eps)
  x0 = point.x-eps
  x1 = point.x+eps
  y0 = point.y-eps
  y1 = point.y+eps
  if distance(point,Point(x0,y0,0,-1)) > eps && distance(point,Point(x0,y1,0,-1)) > eps && distance(point,Point(x1,y0,0,-1)) > eps && distance(point,Point(x1,y1,0,-1)) > eps
    if point.x > x0 && point.x < x1 && point.y > y0 && point.y < y1
      return 1
    else
      return 0
    end
  else
    return 1
  end
end

# searches for all points that are inside circle of 'radius' radius centered on point

function find_neighbours_inside_circle(rectangles,points,point,radius,rectangle)
    ret_val = Int64[]
    if rectangle.kind == 1
      for i in rectangle.children
        if distance(point,points[i]) < radius
          push!(ret_val,i)
        end
      end
      return ret_val
    end

    for i in rectangle.children
      child = rectangles[i]
      if intersects(child,point,eps) == 1
        append!(ret_val,find_eps(rectangles,points,point,radius,child))
      end
    end
    ret_val
end

function DBSCAN(D,eps,MinPts,rectangles)
   clusters = Cluster[]
   for i = 1:length(D)
     point = D[i]
     if point.visited > 0
       continue
     end
     point.visited = 1
     NeighborPoints = find_neighbours_inside_circle(rectangles,D,point,eps,rectangles[end])
     if length(NeighborPoints) < MinPts
       point.visited = 2
     else
       cl = Cluster(Int64[])
       push!(clusters,cl)
       push!(cl.points,i)
       point.cluster = length(clusters)
       while length(NeighborPoints) > 0
         point2 = D[NeighborPoints[1]]
         if point2.visited != 1
           point2.visited = 1
           NeighborPoints2 = find_neighbours_inside_circle(rectangles,D,point2,eps,rectangles[end])
           if length(NeighborPoints2) >= MinPts
             NeighborPoints = append!(NeighborPoints,NeighborPoints2)
           end
         end
         if point2.cluster == -1
           push!(cl.points,NeighborPoints[1])
           point2.cluster = length(clusters)
         end
         shift!(NeighborPoints)
       end
    end
  end
  clusters
end

function build_Rtree(input_points,leaf_rectangle_size)

  # Nearest - X
  dataset = sort(input_points, by = (x) -> x.x)
  rectangles = Rectangle[]

  # building R-tree
  for ci = 1:ceil(length(dataset)/leaf_rectangle_size)
    points = []
    e = (ci == ceil(length(dataset)/leaf_rectangle_size)) ? length(dataset) : leaf_rectangle_size*ci
    points = dataset[leaf_rectangle_size*(ci-1)+1:e]

    x0 = points[1].x
    x1 = points[end].x
    y0 = reduce( (x,y) -> x.y < y.y ? x : y,points).y
    y1 = reduce( (x,y) -> x.y > y.y ? x : y,points).y
    push!(rectangles,Rectangle( -1, x0,y0,x1,y1,[(ci-1)*leaf_rectangle_size + 1:e] , 1))
  end

  c_begin = 1
  c_end   = length(rectangles)
  while (c_end - c_begin + 1) > 1
    no = ceil((c_end - c_begin + 1) / 2)
    for ci = 1:no
      if (c_end - c_begin + 1) % 2 == 0 || ci!=no
        cl = rectangles[2*ci-1 + c_begin - 1]
        cr = rectangles[2*ci + c_begin - 1]
        x0 = min(cl.x0,cr.x0)
        x1 = max(cl.x1,cr.x1)
        y0 = min(cl.y0,cr.y0)
        y1 = max(cl.y1,cr.y1)
        cl.parent = c_end + ci
        cr.parent = c_end + ci
        push!(rectangles,Rectangle(-1,x0,y0,x1,y1,[int(2*ci-1++ c_begin - 1):int(2*ci + c_begin - 1)],0))
      else
        cl = rectangles[2*ci-1 + c_begin - 1]
        x0 = cl.x0
        x1 = cl.x1
        y0 = cl.y0
        y1 = cl.y1
        cl.parent = c_end + ci
        push!(rectangles,Rectangle(-1,x0,y0,x1,y1,[int64((2*ci-1 + c_begin - 1))],0))
      end
    end
    c_begin = c_end + 1
    c_end = c_end + no
  end
  dataset,rectangles
end

# O(N log(N))

function smart_clusterizer(candidates)
  @assert length(candidates) > 1
  #             x        y    visited
  unpacked = Point[]

  for i = 1:length(candidates)
    push!(unpacked,Point(candidates[i][1][1],candidates[i][1][2],0,-1))
    push!(unpacked,Point(candidates[i][2][1],candidates[i][2][2],0,-1))
  end

  dataset , rectangles = build_Rtree(unpacked,2)

  clusters = DBSCAN(dataset,0.3,2,rectangles)

  @assert length(clusters) > 0

  max_cluster = clusters[1]
  for cluster in clusters
    if length(cluster.points) > length(max_cluster.points)
      max_cluster = cluster
    end
  end

  all_x = map( (p) -> dataset[p].x ,max_cluster.points)
  all_y = map( (p) -> dataset[p].y ,max_cluster.points)
  x = mean(all_x)
  y = mean(all_y)
  dx = sigma(all_x,x)
  dy = sigma(all_y,y)
  (x,y, dx, dy)
end

# Naive but simple
#
# Returns tuple (x,y,dx,dy)
#
# We do not control the error in this computation, hence -1

function compute_position_naive(ranges, beacon_positions)
  @assert length(ranges) == length(beacon_positions)
  @assert length(ranges) == 3 "We support only 3 beacons for now"
  bp = beacon_positions
  A = [
        2.0 * (bp[2][1] - bp[1][1])   2.0 * (bp[2][2] - bp[1][2]) ;
        2.0 * (bp[3][1] - bp[1][1])   2.0 * (bp[3][2] - bp[1][2])
      ]
  B = [
        (ranges[1]^2 - ranges[2]^2) - ((bp[1][1]^2-bp[2][1]^2) + (bp[1][2]^2-bp[2][2]^2)) ;
        (ranges[1]^2 - ranges[3]^2) - ((bp[1][1]^2-bp[3][1]^2) + (bp[1][2]^2-bp[3][2]^2))
      ]
  sol = inv(A) * B
  (sol[1],sol[2],-1,-1)
end

# More robust position computation. Many beacons (>= 3) can be used to improve accuracy
#
# Returns tuple (x,y,dx,dy)
#

function compute_position(ranges, beacon_positions)
  @assert length(ranges) == length(beacon_positions)
  @assert length(ranges) >= 3 "At least three points are needed for unambigous position computation"

  candidates = ((Float64,Float64),(Float64,Float64))[]
  bp = beacon_positions
  for i = 1:length(bp)
   for j = 1:length(bp)
     if i >=j
       continue
     end
     push!(candidates,compute_possible_points(ranges[i],ranges[j],bp[i],bp[j]))
   end
  end
  smart_clusterizer(candidates)
end

r = [1.0,sqrt(2.0),1.0]
p = [(1.0,0.0),(0.0,0.0),(0.0,1.0)]

println("Naive position: ")
println(compute_position_naive(r,p))
println("Good position: ")
println(compute_position(r,p))

println("Now with distortion")
eps = 0.1

r = [1+eps,sqrt(2)-eps,1]
p = [(1.0,0.0),(0.0,0.0),(0.0,1.0)]

println("Naive position: ")
println(compute_position_naive(r,p))
println("Good position: ")
println(compute_position(r,p))

max_stars_repo_name: mll/iBeacon-supporting-calculations

max_issues_repo_name: mll/iBeacon-supporting-calculations

max_forks_repo_name: mll/iBeacon-supporting-calculations


text:
module _ExampleData
export ExampleData
using DelimitedFiles
using HTTP

    function ExampleData(SigName)

    Chk = ["uniform","uniform2","randintegers2","randintegers",
        "henon","chirp","gaussian","gaussian2","lorenz",
        "uniform_Mat", "gaussian_Mat", "entropyhub_Mat",
        "mandelbrot_Mat","randintegers_Mat"]

    SigName in Chk ? nothing : error("SigName must be one of the following:\n$Chk")
    url = "https://raw.githubusercontent.com/MattWillFlood/EntropyHub/main/ExampleData/" * SigName * ".txt"
    Temp = HTTP.get(url).body

        if (SigName in ["henon","lorenz"]) || (SigName[end] == '2') || (SigName[end-2:end] == "Mat")
            X = readdlm(Temp, skipstart=2)
        else
            X = readdlm(Temp, skipstart=2);
            size(X)[1] < size(X)[2] ? X = X[1,:] : X = X[:,1]
        end

    return X
    end

end


text:
using Test
include("../src/PortableMinskyMachineNotation.jl")
using .PortableMinskyMachineNotation

@testset "Standard syntax" begin
    standardmachine = pmmn"""
/* Compiles the core grammar and extensions */
inc_by(0, 88);
inc_by(1, 112);
inc_by(2, 120);
inc_by(3, 34);
output(0);
output(1);
output(2);
output(3);
"""
    @test standardmachine isa Function
    @test standardmachine() == "Wow!"
end

@testset "Extended syntax" begin
    julianmachine = pmmn"""
# supports Julia style syntax too
inc(1, 2)
while dec(1)
    input(0)
    output(0)
end
"""
    @test julianmachine(":D") == ":D"

    chunkybaconmachine = pmmn"""
# or pretend it's Ruby if that's your jam
inc 1, 8
while dec 1
  input 0
  dec 0
  output 0
end
"""
    @test chunkybaconmachine("IBM.:111") == "HAL-9000"
end

@testset "Loading files" begin
    @test load("program.txt")() == "Load programs from files!"
end

max_stars_repo_name: stellartux/Esolangs

max_issues_repo_name: stellartux/Esolangs

max_forks_repo_name: stellartux/Esolangs


text:
using MinHash
using DataStructures
using Test

@testset "HashSet" begin
    include("hashset.jl")
end

@testset "MinHashing" begin
    include("minhashing.jl")
end

@testset "Distance" begin
    include("distance.jl")
end


text:
module test_loltools_summonerv4

using Test
using LOLTools.SummonerV4

include("helpers.jl")

struct SummonerController <: ApplicationController
    conn::Conn
end

function summoner_by_name(c::SummonerController)
    # @info :c c.params.summonerName
    render(JSON, (a=1,))
end

routes() do
    get("/lol/summoner/v4/summoners/by-name/:summonerName", SummonerController, summoner_by_name)
end

api_key = ""
platform = "kr"
SummonerV4.summoner_by_name(api_key, platform, "오뚜기3분미트볼"; action=mock_action)

end # module test_loltools_summonerv4


text:
"""
`df2latex(df::DataFrame, outpath::AbstractString)` generates the latex table according to the dataframe.
The result is copied to clipboard; nothing will be returned.
"""
function df2latex(df::DataFrame)
	str = let
		io = IOBuffer();
		show(io, MIME("text/latex"), df);
		str = String(take!(io))
	end
	clipboard(str);
	nothing
end


text:
module TestFunctions

export Benchmark, Unconstrained

"""
    Benchmark

Abstract supertype for all benchmark functions.
"""
abstract type Benchmark <: Function end

"""
    Unconstrained

Abstract supertype for all unconstrained benchmark functions.
"""
abstract type Unconstrained <: Benchmark end

export evaluate, Sphere, Ackley, Rosenbrock, Easom, GoldsteinPrice,
Beale, Levy
include("implementations.jl")

end


text:
using Rasters, Test, Statistics, Dates, Plots
using Rasters.LookupArrays, Rasters.Dimensions
import ArchGDAL, NCDatasets, HDF5, CFTime
using Rasters: layerkeys, SMAPfile, FileArray

testpath = joinpath(dirname(pathof(Rasters)), "../test/")
include(joinpath(testpath, "test_utils.jl"))

# TODO example files without a login requirement
path1 = joinpath(testpath, "data/SMAP_L4_SM_gph_20160101T223000_Vv4011_001.h5")
path2 = joinpath(testpath, "data/SMAP_L4_SM_gph_20160102T223000_Vv4011_001.h5")

smapkeys = (
    :baseflow_flux, :heat_flux_ground, :heat_flux_latent, :heat_flux_sensible,
    :height_lowatmmodlay, :land_evapotranspiration_flux, :land_fraction_saturated,
    :land_fraction_snow_covered, :land_fraction_unsaturated, :land_fraction_wilting,
    :leaf_area_index, :net_downward_longwave_flux, :net_downward_shortwave_flux,
    :overland_runoff_flux, :precipitation_total_surface_flux, :radiation_longwave_absorbed_flux,
    :radiation_shortwave_downward_flux, :sm_profile, :sm_profile_pctl, :sm_profile_wetness,
    :sm_rootzone, :sm_rootzone_pctl, :sm_rootzone_wetness, :sm_surface, :sm_surface_wetness,
    :snow_depth, :snow_mass, :snow_melt_flux, :snowfall_surface_flux, :soil_temp_layer1,
    :soil_temp_layer2, :soil_temp_layer3, :soil_temp_layer4, :soil_temp_layer5, :soil_temp_layer6,
    :soil_water_infiltration_flux, :specific_humidity_lowatmmodlay, :surface_pressure, :surface_temp,
    :temp_lowatmmodlay, :vegetation_greenness_fraction, :windspeed_lowatmmodlay
)

if isfile(path1) && isfile(path2)
    # We need to wrap HDF5 for SMAP, as h5 file may not be SMAP files
    @testset "SMAPhdf5 wrapper" begin
        HDF5.h5open(path1) do f
            ds= Rasters.SMAPhdf5(f)
            @test keys(ds) == layerkeys(ds) == smapkeys
            @test dims(ds) isa Tuple{<:X,<:Y}
        end
    end

    @testset "Raster" begin
        @time smaparray = Raster(path1)

        @testset "lazyness" begin
            @time read(Raster(path1));
            @time lazyarray = Raster(path1; lazy=true);
            @time eagerarray = Raster(path1; lazy=false);
            # Lazy is the default
            @test parent(smaparray) isa FileArray
            @test parent(lazyarray) isa FileArray
            @test parent(eagerarray) isa Array
        end

        @testset "open" begin
            @test all(open(A -> A[Y=1], smaparray) .=== smaparray[:, 1])
        end

        @testset "read" begin
            @time A = read(smaparray);
            @test A isa Raster
            @test parent(A) isa Array
            A2 = zero(A)
            read!(smaparray, A2)
            A3 = zero(A)
            @time read!(path1, A3);
            @test A == A2 == A3
        end

        @testset "array properties" begin
            @test smaparray isa Raster
        end

        @testset "dimensions" begin
            @test ndims(smaparray) == 2
            HDF5.h5open(path1) do ds
                @test size(smaparray) == length.(dims(smaparray)) == size(ds["Geophysical_Data/baseflow_flux"])
            end
            @test dims(smaparray) isa Tuple{<:X,<:Y}
            @test refdims(smaparray) == ()
            @test sampling(smaparray) == (Intervals(Center()), Intervals(Center()))
            @test span(smaparray) == (Irregular((-180.0f0, 180.0f0)), Irregular((-85.04456f0, 85.04456f0)))
            @test bounds(smaparray) == ((-180.0f0, 180.0f0), (-85.04456f0, 85.04456f0))
        end

        @testset "other fields" begin @test missingval(smaparray) == -9999.0
            @test metadata(smaparray) isa Metadata{SMAPfile}
            @test name(smaparray) == :baseflow_flux
        end

        @testset "indexing" begin
            @test smaparray[Ti(1)] isa Raster{<:Any,2}
            @test smaparray[Y(1), Ti(1)] isa Raster{<:Any,1}
            @test smaparray[X(1), Ti(1)] isa Raster{<:Any,1}
            @test smaparray[X(1), Y(1)] == -9999.0
            @test smaparray[X(30), Y(30)] isa Float32
            # Russia
            @test smaparray[X(50), Y(100)] == -9999.0
            # Alaska
            @test smaparray[Y(Near(64.2008)), X(Near(149.4937))] == 0.0f0
        end

        @testset "selectors" begin
            a = smaparray[X(Near(21.0)), Y(Between(50, 52))]
            @test bounds(a) == ((50.08451f0, 51.977905f0),)
            x = smaparray[X(Near(150)), Y(Near(30))]
            @test x isa Float32
            dimz = X(Between(-180.0, 180)), Y(Between(-90, 90))
            @test size(smaparray[dimz...]) == (3856, 1624)
            @test index(smaparray[dimz...]) == index(smaparray)
            nca = smaparray[Y(Between(-80, -25)), X(Between(0, 180))]
        end

        @testset "methods" begin
            @test all(mean(smaparray; dims=Y) .=== mean(parent(smaparray); dims=2))
            @testset "trim, crop, extend" begin
                a = read(smaparray)
                a[X(1:20)] .= missingval(a)
                trimmed = trim(a)
                @test size(trimmed) == (3836, 1512)
                cropped = crop(a; to=trimmed)
                dims(trimmed, X)
                dims(cropped, X)
                @test size(cropped) == (3836, 1512)
                @test all(collect(cropped .=== trimmed))
                dims(cropped)
                extended = extend(cropped; to=a)
                @test all(collect(extended .=== a))
            end
            @testset "slice" begin
                ser = Rasters.slice(smaparray, X)
                @test ser isa RasterSeries
                @test size(ser) == (3856,)
                @test index(ser, X) == index(smaparray, X)
                @test bounds(ser, X) == (-180.0f0, 180.0f0)
                A = ser[1]
                @test bounds(A, Y) == (-85.04456f0, 85.04456f0)
            end
        end

        @testset "save" begin
            @testset "to grd" begin
                # TODO save and load subset
                geoA = read(smaparray)
                @test size(geoA) == size(smaparray)
                filename = tempname() * ".grd"
                write(filename, geoA)
                saved = read(Raster(filename; mappedcrs=EPSG(4326))[Band(1)])
                dims(saved)
                @test size(saved) == size(geoA)
                @test missingval(saved) === missingval(geoA)
                @test map(metadata.(dims(saved)), metadata.(dims(Raster))) do s, g
                    all(s .== g)
                end |> all
                @test Rasters.name(saved) == Rasters.name(geoA)
                @test all(lookup.(dims(saved)) .!= lookup.(dims(geoA)))
                @test all(
                          mappedbounds(saved)[1]
                          .≈
                          bounds(geoA)[1])
                @test all(mappedbounds(saved)[2] .≈ mappedbounds(geoA)[2])
                bounds(saved)
                @test all(parent(saved) .=== parent(geoA))
            end
            @testset "to gdal" begin
                gdalfilename = tempname() * ".tif"
                @time write(gdalfilename, read(smaparray))
                gdalarray = Raster(gdalfilename; mappedcrs=EPSG(4326))
                # These come out with slightly different format
                # @test convert(ProjString, crs(gdalarray)) == crs(smaparray)
                @test all(map((a, b) -> all(a .≈ b),
                              mappedbounds(dims(gdalarray))
                              , (bounds(smaparray)..., (1, 1))))
                # Tiff locus = Start, SMAP locus = Center
                @test mappedindex(DimensionalData.shiftlocus(Center(), dims(gdalarray, Y))) ≈ index(smaparray, Y)
                @test mappedindex(DimensionalData.shiftlocus(Center(), dims(gdalarray, X))) ≈ index(smaparray, X)
                @test all(gdalarray .== read(smaparray))
            end
            @testset "to netcdf" begin
                ncdfilename = tempname() * ".nc"
                @time write(ncdfilename, smaparray)
                reorder(smaparray, ForwardOrdered())
                saved = Raster(ncdfilename)
                @test_broken bounds(saved) == bounds(smaparray)
                @test index(saved, Y) == index(smaparray, Y)
                @test index(saved, X) == index(smaparray, X)
            end
        end

        @testset "show" begin
            sh = sprint(show, MIME("text/plain"), smaparray)
            @test occursin("Raster", sh)
            @test occursin("Y", sh)
            @test occursin("X", sh)
        end

        @testset "plot" begin
            smaparray[Ti(1:3:12)] |> plot # FIXME plot is weird
            smaparray[Ti(1)] |> plot
            smaparray[Y(100), Ti(1)] |> plot
        end

    end

    @testset "stack" begin
        @time smapstack = RasterStack(path1)

        @testset "lazyness" begin
            @time read(RasterStack(path1));
            @time lazystack = RasterStack(path1; lazy=true)
            @time eagerstack = RasterStack(path1; lazy=false);
            # Lazy is the default
            @test parent(smapstack[:snow_mass]) isa FileArray
            @test parent(lazystack[:snow_mass]) isa FileArray
            @test parent(eagerstack[:snow_mass]) isa Array
        end

        @testset "read" begin
            @time st = read(smapstack);
            @test st isa RasterStack
            @test st.data isa NamedTuple
            @test first(st.data) isa Array
            st2 = map(a -> a .* 0, st)
            @time read!(smapstack, st2);
            st3 = map(a -> a .* 0, st)
            @time st = read!(path1, st3);
            @test all(map((a, b, c) -> all(a .== b .== c), st, st2, st3))
        end

        @testset "conversion to Raster" begin
            smaparray = smapstack[:soil_temp_layer1];
            @test smaparray isa Raster{Float32,2}
            @test dims(smaparray) isa Tuple{<:X{<:Mapped{Float32}}, <:Y{<:Mapped{Float32}}}
            @test span(smaparray) isa Tuple{Irregular{Tuple{Float32,Float32}},Irregular{Tuple{Float32,Float32}}}
            @test span(smaparray) == (Irregular((-180.0f0, 180.0f0)), Irregular((-85.04456f0, 85.04456f0)))
            @test bounds(smaparray) == ((-180.0f0, 180.0f0), (-85.04456f0, 85.04456f0))
            @test index(smaparray) isa Tuple{Vector{Float32},Vector{Float32}}
            @test missingval(smaparray) == -9999.0
            @test smaparray[1, 1] == -9999.0
            @test name(smaparray) == :soil_temp_layer1
            dt = DateTime(2016, 1, 1, 22, 30)
            step_ = Hour(3)
            # Currently empty
            @test metadata(smaparray) isa Metadata{SMAPfile}
        end

        @testset "conversion to regular RasterStack" begin
            # Stack Constructors
            geostack = RasterStack(smapstack)
            @test parent(geostack) isa NamedTuple;
            @test Symbol.(Tuple(keys(smapstack))) == keys(geostack)
            geostack = RasterStack(smapstack; keys=(:baseflow_flux, :snow_mass, :soil_temp_layer1))
            @test keys(geostack) == (:baseflow_flux, :snow_mass, :soil_temp_layer1)
        end

        if VERSION > v"1.1-"
            @testset "copy" begin
                geoA = zero(smapstack[:soil_temp_layer1])
                @test geoA isa Raster
                @test geoA != smapstack[:soil_temp_layer1]
                copy!(geoA, smapstack, :soil_temp_layer1)
                @test geoA == smapstack[:soil_temp_layer1]
            end
        end

        @testset "show" begin
            sh1 = sprint(show, MIME("text/plain"), smapstack[:soil_temp_layer1])
            @test occursin("Raster", sh1)
            @test occursin("Y", sh1)
            @test occursin("X", sh1)
        end

    end

    @testset "series" begin
        ser = smapseries([path1, path2])
        val.(dims(ser))
        @test ser[1] isa RasterStack
        @test first(bounds(ser, Ti)) == DateTime(2016, 1, 1, 22, 30)
        @test last(bounds(ser, Ti)) == DateTime(2016, 1, 3, 1, 30)
        @time modified_series = modify(Array, ser)
        stackkeys = keys(modified_series[1])
        @test typeof(modified_series) <: RasterSeries{<:RasterStack{<:NamedTuple{stackkeys,<:Tuple{<:Array{Float32,2,},Vararg}}}}
        @testset "read" begin
            # FIXME: uses too much memory
            # Seems to be a HDF5.jl memory leak?
            # @time geoseries = read(ser)
            # @test geoseries isa RasterSeries{<:RasterStack}
            # @test geoseries.data isa Vector{<:RasterStack}
            # @test first(geoseries.data[1].data) isa Array
        end

        @testset "show" begin
            sh = sprint(show, MIME("text/plain"), ser)
            # Test but don't lock this down too much
            @test occursin("RasterSeries", sh)
        end
    end
end


text:
@testset "Core" begin
    e2 = Edge(1, 3)
    e3 = Edge(1, 4)
    @test @inferred(is_ordered(e2))
    @test @inferred(!is_ordered(reverse(e3)))

    @testset "add_vertices!" begin
        gx = SimpleGraph(10); gdx = SimpleDiGraph(10)
        @testset "$g" for g in testgraphs(gx, gdx)
            gc = copy(g)
            @test add_vertices!(gc, 5) == 5
            @test @inferred(nv(gc)) == 15
        end
    end

    g5w = wheel_graph(5); g5wd = wheel_digraph(5)
    @testset "degree functions" begin
        @testset "$g" for g in testgraphs(g5w)
            @test @inferred(indegree(g, 1)) == @inferred(outdegree(g, 1)) == 4
            @test degree(g, 1) == 4 # explicit codecov
            @test @inferred(indegree(g)) == @inferred(outdegree(g)) == @inferred(degree(g)) == [4, 3, 3, 3, 3]

            @test @inferred(Δout(g)) == @inferred(Δin(g)) == @inferred(Δ(g)) == 4
            @test @inferred(δout(g)) == @inferred(δin(g)) == @inferred(δ(g)) == 3
            z1 = @inferred(degree_histogram(g))
            z2 = @inferred(degree_histogram(g, indegree))
            z3 = @inferred(degree_histogram(g, outdegree))
            @test z1 == z2 == z3 == Dict(4 => 1, 3 => 4)
        end
        @testset "$g" for g in testdigraphs(g5wd)
            @test @inferred(indegree(g, 2)) == 2
            @test @inferred(outdegree(g, 2)) == 1
            @test @inferred(degree(g, 2)) == 3
            @test @inferred(indegree(g)) == [0, 2, 2, 2, 2]
            @test @inferred(outdegree(g)) == [4, 1, 1, 1, 1]
            @test @inferred(degree(g)) == [4, 3, 3, 3, 3]
            @test @inferred(Δout(g)) == @inferred(Δ(g)) == 4
            @test @inferred(Δin(g)) == 2
            @test @inferred(δout(g)) == 1
            @test @inferred(δin(g)) == 0
            @test @inferred(δ(g)) == 3
            z1 = @inferred(degree_histogram(g))
            z2 = @inferred(degree_histogram(g, indegree))
            z3 = @inferred(degree_histogram(g, outdegree))
            @test z1 == Dict(4 => 1, 3 => 4)
            @test z2 == Dict(0 => 1, 2 => 4)
            @test z3 == Dict(4 => 1, 1 => 4)
        end
    end

    @testset "weights" begin
        @testset "$g" for g in testgraphs(g5w, g5wd)
            @test @inferred(weights(g)) == LightGraphs.DefaultDistance(nv(g))
        end
    end

    @testset "neighbor functions" begin
        @testset "$g" for g in testgraphs(g5w)
            @test @inferred(neighbors(g, 2)) == [1, 3, 5]
            @test @inferred(all_neighbors(g, 2)) == [1, 3, 5]
            @test @inferred(common_neighbors(g, 1, 5)) == [2, 4]
        end
        @testset "$g" for g in testgraphs(g5wd)
            @test @inferred(neighbors(g, 2)) == [3]
            @test Set(@inferred(all_neighbors(g, 2))) == Set([1, 3, 5])
            @test @inferred(common_neighbors(g, 1, 5)) == [2]
        end
    end

    @testset "self loops" begin
        @testset "$g" for g in testgraphs(g5w, g5wd)
            gsl = copy(g)
            add_edge!(gsl, 3, 3)
            add_edge!(gsl, 2, 2)
            @test @inferred(!has_self_loops(g))
            @test @inferred(has_self_loops(gsl))
            @test @inferred(num_self_loops(g)) == 0
            @test @inferred(num_self_loops(gsl)) == 2
        end
    end

    @testset "density" begin
        @testset "$g" for g in testgraphs(g5w)
            @test @inferred(density(g)) == 0.8
        end

        @testset "$g" for g in testdigraphs(g5wd)
            @test @inferred(density(g)) == 0.4
        end

    end
    @testset "squash" begin
        @testset "$g" for g in testgraphs(g5w, g5wd)
            @test eltype(squash(g)) == UInt8
        end
    end
end


text:
# Hodgkin Huxley with current pulse setup for testing
using Conductor, ModelingToolkit, OrdinaryDiffEq, Unitful
using IfElse
using Conductor: Na, K
using Unitful: mV, mS, cm, µm, pA, nA, mA, µA, ms

Vₘ = MembranePotential()

nav_kinetics = [
    Gate(AlphaBetaRates,
         αₘ = IfElse.ifelse(Vₘ == -40.0, 1.0, (0.1*(Vₘ + 40.0))/(1.0 - exp(-(Vₘ + 40.0)/10.0))),
         βₘ = 4.0*exp(-(Vₘ + 65.0)/18.0), p = 3)
    Gate(AlphaBetaRates,
         αₕ = 0.07*exp(-(Vₘ+65.0)/20.0),
         βₕ = 1.0/(1.0 + exp(-(Vₘ + 35.0)/10.0)))]

kdr_kinetics = [
    Gate(AlphaBetaRates,
         αₙ = IfElse.ifelse(Vₘ == -55.0, 0.1, (0.01*(Vₘ + 55.0))/(1.0 - exp(-(Vₘ + 55.0)/10.0))),
         βₙ = 0.125 * exp(-(Vₘ + 65.0)/80.0), p = 4)]

@named NaV = IonChannel(Sodium, nav_kinetics, 120mS/cm^2)
@named Kdr = IonChannel(Potassium, kdr_kinetics, 36mS/cm^2)
@named leak = PassiveChannel(Leak, 0.3mS/cm^2)

gradients = Equilibria([Na   =>  50.0mV, K    => -77.0mV, Leak => -54.4mV])
area = 4*pi*(20µm)^2
pulse(t, current) = IfElse.ifelse(t > 100.0, IfElse.ifelse(t < 200.0, 0.0004, 0.0), 0.0)

@named neuron = Soma([NaV,Kdr,leak], gradients, stimulus = pulse, area = ustrip(Float64, cm^2, area));

t = 300.
sim = Simulation(neuron, time = t*ms) # ODESystem -> structural_simplify -> ODAEProblem output

### Hand-written model setup ###
# Steady-state functions
minf(V0) = inv((4.0 + 0.1V0)*inv(1.0 - exp(-4.0 - 0.1V0)) + 4.0*exp(-3.6111111111111107 - 0.05555555555555555V0)) *
           (4.0 + 0.1V0) * inv(1.0 - exp(-4.0 - 0.1V0))
ninf(V0) = (0.55 + 0.01V0) * inv(1.0 - exp(-5.5 - 0.1V0)) * inv(0.125exp(-0.8125 - 0.0125V0) +
                                                                (0.55 + 0.01V0) * inv(1.0 - exp(-5.5 - 0.1V0)))
hinf(V0) = 0.07*exp(-3.25 - 0.05V0) * inv(0.07*exp(-3.25 - 0.05V0) + inv(1.0 + exp(-3.5 - 0.1V0)) )
hand_pulse(t) = ifelse(100. < t < 200. , 0.0004 , 0.0)

# Initial conditions
V0 = -65.0
u0 = zeros(Float64, 4)
p = zeros(Float64, 8)

# Serialized MTK function initial conditions
u0[1] = V0       # Vₘ
u0[2] = minf(V0) # NaV₊m
u0[3] = hinf(V0) # NaV₊h
u0[4] = ninf(V0) # Kdr₊n

# Parameters
p[1]  = 0.001      # cₘ
p[2]  = 5.02655e-5 # aₘ
p[3]  = 50.0       # ENa
p[4]  = -77.0      # EK
p[5]  = -54.4      # El
p[6]  = 120.0      # NaV₊gbar
p[7]  = 36.0       # Kdr₊gbar
p[8]  = 0.3        # leak₊g

# Cleaned up version of what default MTK output looks like
function hodgkin_huxley!(du, u, p, t)
    @inbounds let Vₘ         = u[1],
                  NaV₊m      = u[2],
                  NaV₊h      = u[3],
                  Kdr₊n      = u[4],

                  cₘ          = p[1],
                  aₘ          = p[2],
                  ENa         = p[3],
                  EK          = p[4],
                  El          = p[5],
                  NaV₊gbar    = p[6],
                  Kdr₊gbar    = p[7],
                  leak₊g      = p[8],
                  Iapp        = hand_pulse(t)
    du[1] = (Iapp - aₘ*leak₊g*(Vₘ - El) - Kdr₊gbar*aₘ*Kdr₊n^4*(Vₘ - EK) -
             NaV₊gbar*aₘ*NaV₊h*NaV₊m^3*(Vₘ - ENa))*inv(aₘ)*inv(cₘ)
    du[2] = ifelse(Vₘ == -40.0 , 1.0 , (4.0 + 0.1Vₘ) * inv(1.0 - exp(-4.0 - 0.1Vₘ))) *
            (1 - NaV₊m) - 4.0*NaV₊m*exp(-3.6111111111111107 + -0.05555555555555555Vₘ)

    du[3] = 0.07*exp(-3.25 - 0.05Vₘ)*(1 - NaV₊h) - NaV₊h*inv(1.0 + exp(-3.5 - 0.1Vₘ))
    du[4] = ifelse(Vₘ == -55.0 , 0.1 , (0.55 + 0.01Vₘ) * inv(1.0 - exp(-5.5 - 0.1Vₘ)))*(1 - Kdr₊n) -
            0.125 * Kdr₊n * exp(-0.8125 - 0.0125Vₘ)
    end
    return nothing
end

# Original time
# ~99ns+/-5ns

# Int exponents -- this is a sizeable optimization
# ~87ns+/-5ns

# ifelse
# ~85.5+/-4

# static params
# 84.3 +/- 3


text:
module AutoMerge

using Dates: Dates
# import GitCommand
using GitHub: GitHub
using HTTP: HTTP
using LibGit2: LibGit2
using Pkg: Pkg
using TimeZones: TimeZones
using JSON: JSON
using VisualStringDistances: VisualStringDistances
using StringDistances: StringDistances
using LicenseCheck: LicenseCheck
using TOML: TOML
using Printf: Printf
using RegistryTools: RegistryTools
using ..RegistryCI: RegistryCI
using Tar: Tar

include("types.jl")
include("ciservice.jl")

include("api_rate_limiting.jl")
include("assert.jl")
include("automerge_comment.jl")
include("changed_files.jl")
include("cron.jl")
include("dates.jl")
include("dependency_confusion.jl")
include("github.jl")
include("guidelines.jl")
include("jll.jl")
include("not_automerge_applicable.jl")
include("package_path_in_registry.jl")
include("public.jl")
include("pull_requests.jl")
include("semver.jl")
include("update_status.jl")
include("util.jl")

end # module


text:
struct GeneralWilsonFermiAction{Dim,Dirac,fermion,gauge} <: FermiAction{Dim,Dirac,fermion,gauge}
    hascovnet::Bool
    covneuralnet::Union{Nothing,CovNeuralnet{Dim}}
    diracoperator::Dirac
    _temporary_fermionfields::Vector{fermion}
    _temporary_gaugefields::Vector{gauge}

    function GeneralWilsonFermiAction(D::Dirac_operator{Dim},hascovnet,covneuralnet,parameters_action) where Dim
        num = 6

        x = D._temporary_fermi[1]
        xtype = typeof(x)
        _temporary_fermionfields = Array{xtype,1}(undef,num)
        for i=1:num
            _temporary_fermionfields[i] = similar(x)
        end

        Utemp = D.U[1]
        Utype = typeof(Utemp)
        numU = 2
        _temporary_gaugefields = Array{Utype,1}(undef,numU)
        for i=1:numU
            _temporary_gaugefields[i] = similar(Utemp)
        end

        return new{Dim,typeof(D),xtype,Utype}(hascovnet,
                        covneuralnet,D,
                        _temporary_fermionfields,
                        _temporary_gaugefields)

    end
end

function evaluate_FermiAction(fermi_action::GeneralWilsonFermiAction{Dim,Dirac,fermion,gauge},U,ϕ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge}
    W = fermi_action.diracoperator(U)
    η = fermi_action._temporary_fermionfields[1]
    solve_DinvX!(η,W',ϕ)
    Sf = dot(η,η)
    return real(Sf)
end

function gauss_sampling_in_action!(η::AbstractFermionfields,U,fermi_action::GeneralWilsonFermiAction{Dim,Dirac,fermion,gauge} ) where {Dim,Dirac,fermion,gauge}
    #gauss_distribution_fermion!(η)
    gauss_distribution_fermion!(η,rand)
end

function sample_pseudofermions!(ϕ::AbstractFermionfields,U,fermi_action::GeneralWilsonFermiAction{Dim,Dirac,fermion,gauge} ,ξ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge}
    W = fermi_action.diracoperator(U)
    mul!(ϕ,W',ξ)
    set_wing_fermion!(ϕ)
end


text:
@testset "memory allocator" begin

CuArrays.reclaim(1024)
CuArrays.reclaim()

CuArrays.extalloc(()->())
CuArrays.extalloc(()->(); check=ex->true)
CuArrays.extalloc(()->(); nb=1)

@test_throws OutOfGPUMemoryError CuArray{Int}(undef, 10^20)
@test_throws OutOfGPUMemoryError CuArrays.extalloc() do
    CuArray{Int}(undef, 10^20)
end

end


text:
# A script that compares row storage to column storage for rods
# Also, a good reference for the most highly optimized kinematics
# routines I could come up with for each storage pattern.
# Unsurprisingly, row storage is never slower, and often 2-3x faster
# due to better memory access patterns.

using BenchmarkTools
N = 500

# ------------------------------------------------------------------------------
# Cross Products (Rows Win)
# ------------------------------------------------------------------------------
function cross3_col!(c,a,b)
    @inbounds for i = 1:size(c,2)
        c[1,i] = a[2,i]*b[3,i] - a[3,i]*b[2,i]
        c[2,i] = a[3,i]*b[1,i] - a[1,i]*b[3,i]
        c[3,i] = a[1,i]*b[2,i] - a[2,i]*b[1,i]
    end
end

function cross3_row!(c,a,b)
    a1 = @view a[:,1]
    a2 = @view a[:,2]
    a3 = @view a[:,3]

    b1 = @view b[:,1]
    b2 = @view b[:,2]
    b3 = @view b[:,3]

    c[:,1] .= @. a2*b3 - a3*b2
    c[:,2] .= @. a3*b1 - a1*b3
    c[:,3] .= @. a1*b2 - a2*b1
end

a = rand(3,N)
b = rand(3,N)
c = rand(3,N)
@btime cross3_col!($c,$a,$b)

a = rand(N,3)
b = rand(N,3)
c = rand(N,3)
@btime cross3_row!($c,$a,$b)

# ------------------------------------------------------------------------------
# Dot Products (Rows Win)
# ------------------------------------------------------------------------------
function dot3_col!(c,a,b)
    @inbounds for i = 1:length(c)
        c[i] = a[1,i]*b[1,i] + a[2,i]*b[2,i] + a[3,i]*b[3,i]
    end
end

function dot3_row!(c,a,b)
    @inbounds for i = 1:length(c)
        c[i] = a[i,1]*b[i,1] + a[i,2]*b[i,2] + a[i,3]*b[i,3]
    end
end

a = rand(3,N)
b = rand(3,N)
c = rand(1,N)
@btime dot3_col!($c,$a,$b)

a = rand(N,3)
b = rand(N,3)
c = rand(N,1)
@btime dot3_row!($c,$a,$b)

# ------------------------------------------------------------------------------
# Normalizing In Place (rows win, by less)
# ------------------------------------------------------------------------------
function norm3_col!(a)
    @inbounds for i = 1:size(a,2)
        tmp = sqrt(a[1,i]*a[1,i] + a[2,i]*a[2,i] + a[3,i]*a[3,i])
        a[1,i] /= tmp
        a[2,i] /= tmp
        a[3,i] /= tmp
    end
end

function norm3_row!(a)
    @inbounds for i = 1:size(a,1)
        tmp = sqrt(a[i,1]*a[i,1] + a[i,2]*a[i,2] + a[i,3]*a[i,3])
        a[i,1] /= tmp
        a[i,2] /= tmp
        a[i,3] /= tmp
    end
end

a = rand(3,N)
@btime norm3_col!($a)

a = rand(N,3)
@btime norm3_row!($a)

# ------------------------------------------------------------------------------
# Norms out of Place (rows win, by alot)
# ------------------------------------------------------------------------------
function norm3_col!(out, a)
    @inbounds for i = 1:length(out)
        out[i] = sqrt(a[1,i]*a[1,i] + a[2,i]*a[2,i] + a[3,i]*a[3,i])
    end
end

function norm3_row!(out, a)
    @inbounds for i = 1:length(out)
        out[i] = sqrt(a[i,1]*a[i,1] + a[i,2]*a[i,2] + a[i,3]*a[i,3])
    end
end

a = rand(3,N)
out = rand(1,N)
@btime norm3_col!($out,$a)

a = rand(N,3)
out = rand(N,1)
@btime norm3_row!($out, $a)

# ------------------------------------------------------------------------------
# Edges (rows win, off the charts)
# ------------------------------------------------------------------------------
function edges_col!(e, x)
    @inbounds for i = 1:size(e,2)
        e[1,i] = x[1,i+1] - x[1,i]
        e[2,i] = x[2,i+1] - x[2,i]
        e[3,i] = x[3,i+1] - x[3,i]
    end
end

function edges_row!(e, x)
    nv = size(x, 1)
    e .= view(x, 2:nv, :) .- view(x, 1:nv-1, :)
end

function tangents_col!(t, x)
    edges_col!(t,x)
    norm3_col!(t)
end

function tangents_row!(t, x)
    edges_row!(t,x)
    norm3_row!(t)
end

a = rand(3,N)
e = rand(3,N-1)
@btime edges_col!($e,$a)
@btime tangents_col!($e,$a)

a = rand(N,3)
e = rand(N-1,3)
@btime edges_row!($e,$a)
@btime tangents_row!($e,$a)

# ------------------------------------------------------------------------------
# Parallel Transport (Example, 3x speedup with row storage)
# ------------------------------------------------------------------------------
const PTRANSPORT_TOL2 = 0.00036
function ptransport_inner_coefficient!(out, s2, c)
    out .= 1/2 .+ (1 .- c)./4 .+ (1 .- c).^2 .*(1 .+ s2./4)./8
    @inbounds for i = 1:length(out)
        if s2[i] > PTRANSPORT_TOL2
            out[i] = (1 - c[i])/s2[i]
        end
    end
end

function ptransport_col!(v2, v1, t1, t2, cache)
    X = view(cache, 1:3, :)
    s2 = view(cache, 4:4, :)
    c = view(cache, 5:5, :)
    ptic = view(cache, 6:6, :)

    cross3_col!(X, t1, t2)
    dot3_col!(s2, X, X)
    dot3_col!(c, t1, t2)
    ptransport_inner_coefficient!(ptic,s2,c)

    dXv = view(cache, 4:4, :) # Overwrites s2
    dot3_col!(dXv,X,v1)
    cross3_col!(v2,X,v1)
    v2 .+= c.*v1 .+ ptic.*dXv.*X
end

function ptransport_row!(v2, v1, t1, t2, cache)
    X = view(cache, :, 1:3)
    s2 = view(cache, :, 4)
    c = view(cache, :, 5)
    ptic = view(cache, :, 6)

    cross3_row!(X, t1, t2)
    dot3_row!(s2, X, X)
    dot3_row!(c, t1, t2)
    ptransport_inner_coefficient!(ptic,s2,c)

    dXv = view(cache, :, 4) # Overwrites s2
    dot3_row!(dXv,X,v1)
    cross3_row!(v2,X,v1)
    v2 .+= c.*v1 .+ ptic.*dXv.*X
end

v1 = rand(3,N)
v2 = similar(v1)
t1 = rand(3,N)
t2 = rand(3,N)
cache = zeros(6,N)
@btime ptransport_col!($v2,$v1,$t1,$t2,$cache)

v1 = rand(N,3)
v2 = similar(v1)
t1 = rand(N,3)
t2 = rand(N,3)
cache = zeros(N,6)
@btime ptransport_row!($v2,$v1,$t1,$t2,$cache)

# ------------------------------------------------------------------------------
# Rotation (Example, only minor speedup)
# ------------------------------------------------------------------------------
function rotate_orthogonal_unit_col!(v, t, θ, cache)
    #ne = length(θ)
    #cache = reshape(view(cacheraw,1:5*ne), 5, ne)

    tmp = view(cache, 1:3, :)
    c = view(cache, 4:4, :)
    s = view(cache, 5:5, :)

    c .= cos.(θ)
    s .= sin.(θ)
    cross3_col!(tmp, t, v)
    v .= c.*v .+ s.*tmp
end

function rotate_orthogonal_unit_row!(v, t, θ, cache)
    tmp = view(cache, :, 1:3)
    c = view(cache, :, 4:4)
    s = view(cache, :, 5:5)

    c .= cos.(θ)
    s .= sin.(θ)
    cross3_row!(tmp, t, v)
    v .= c.*v .+ s.*tmp
end

v = rand(3,N)
t = rand(3,N)
norm3_col!(t)
θ = rand(1,N)
cache = zeros(5,N)
@btime rotate_orthogonal_unit_col!($v, $t, $θ, $cache)

v = rand(N,3)
t = rand(N,3)
norm3_row!(t)
θ = rand(N,1)
cache = zeros(N,5)
@btime rotate_orthogonal_unit_row!($v, $t, $θ, $cache)

# ------------------------------------------------------------------------------
# Rod Update (2x Speedup for Rows)
# ------------------------------------------------------------------------------
function rod_update_col!(x, d, Δx, Δθ, cache)
    t0 = view(cache, 1:3, :)
    t = view(cache, 4:6, :)
    d0 = view(cache, 7:9, :)
    cache_pt = view(cache, 10:15, :)
    cache_rt = view(cache, 10:14, :)

    tangents_col!(t0, x)
    x .+= Δx
    tangents_col!(t, x)

    d0 .= d
    ptransport_col!(d, d0, t0, t, cache_pt)
    rotate_orthogonal_unit_col!(d, t, Δθ, cache_pt)
end

function rod_update_row!(x, d, Δx, Δθ, cache)
    t0 = view(cache, :, 1:3)
    t = view(cache, :, 4:6)
    d0 = view(cache, :, 7:9)
    cache_pt = view(cache, :, 10:15)
    cache_rt = view(cache, :, 10:14)

    tangents_row!(t0, x)
    x .+= Δx
    tangents_row!(t, x)

    d0 .= d
    ptransport_row!(d, d0, t0, t, cache_pt)
    rotate_orthogonal_unit_row!(d, t, Δθ, cache_pt)
end

x = rand(3,N)
d = rand(3,N-1)
norm3_col!(d)
Δx = rand(3,N)
Δθ = rand(1,N-1)
cache = zeros(15,N-1)
@btime rod_update_col!($x,$d,$Δx,$Δθ,$cache)

x = rand(N,3)
d = rand(N-1,3)
norm3_row!(d)
Δx = rand(N,3)
Δθ = rand(N-1,1)
cache = zeros(N-1,15)
@btime rod_update_row!($x,$d,$Δx,$Δθ,$cache)

# ------------------------------------------------------------------------------
# Full Kinematics Update (Factor of 2.5 for Row Storage)
# ------------------------------------------------------------------------------
function full_kinematics_col!(l, κ, τ, x, d1, caches)

    # Unpack and divy up caches
    cache1, cache2 = caches
    t = view(cache1, 1:3, :)
    d2 = view(cache1, 4:6, :)
    k = view(cache2, 1:3, :)
    tmp1 = view(cache2, 4:4, :)
    tmp3 = view(cache2, 4:6, :)
    pt_cache = view(cache2, 1:6, :) # overrides k, tmp
    d1trans = view(cache2, 7:9, :)
    cosτ = view(cache2, 1:1, :)
    sinτ = view(cache2, 2:2, :)

    # Lengths
    edges_col!(t, x)
    norm3_col!(l, t)
    t ./= l
    cross3_col!(d2, t, d1)

    # Break tangents and directors into sets
    nv = size(x,2)
    t1 = view(t,:,1:nv-2)
    t2 = view(t,:,2:nv-1)
    d1l = view(d1,:,1:nv-2)
    d1r = view(d1,:,2:nv-1)
    d2l = view(d2,:,1:nv-2)
    d2r = view(d2,:,2:nv-1)

    # Curvature
    cross3_col!(k, t1, t2)
    dot3_col!(tmp1, t1, t2)
    k .*= 2.0./(1.0 .+ tmp1)
    tmp3 .= k.*(d2l .+ d2r)./2
    sum!(view(κ, 1:1, :), tmp3)
    tmp3 .= -1.0.*k.*(d1l .+ d1r)./2
    sum!(view(κ, 2:2, :), tmp3)

    # # Twists (could move in part of ptransport to reduce intermediates)
    ptransport_col!(d1trans, d1l, t1, t2, pt_cache)
    dot3_col!(cosτ, d1r, d1trans)
    dot3_col!(sinτ, d2r, d1trans)
    τ .= -1.0.*atan.(sinτ, cosτ)
end

function full_kinematics_row!(l, κ, τ, x, d1, caches)

    # Unpack and divy up caches
    cache1, cache2 = caches
    t = view(cache1, :, 1:3)
    d2 = view(cache1, :, 4:6)
    k = view(cache2, :, 1:3)
    tmp1 = view(cache2, :, 4:4)
    tmp3 = view(cache2, :, 4:6)
    pt_cache = view(cache2, :, 1:6) # overrides k, tmp
    d1trans = view(cache2, :, 7:9)
    cosτ = view(cache2, :, 1)
    sinτ = view(cache2, :, 2)

    # Lengths
    edges_row!(t, x)
    norm3_row!(l, t)
    t ./= l
    cross3_row!(d2, t, d1)

    # Break tangents and directors into sets
    nv = size(x,1)
    t1 = view(t, 1:nv-2, :)
    t2 = view(t, 2:nv-1, :)
    d1l = view(d1, 1:nv-2, :)
    d1r = view(d1, 2:nv-1, :)
    d2l = view(d2, 1:nv-2, :)
    d2r = view(d2, 2:nv-1, :)

    # Curvature
    cross3_row!(k, t1, t2)
    dot3_row!(tmp1, t1, t2)
    k .*= 2.0./(1.0 .+ tmp1)
    tmp3 .= k.*(d2l .+ d2r)./2
    sum!(view(κ, :, 1:1), tmp3)
    tmp3 .= -1.0.*k.*(d1l .+ d1r)./2
    sum!(view(κ, :, 2:2), tmp3)

    # # Twists (could move in part of ptransport to reduce intermediates)
    ptransport_row!(d1trans, d1l, t1, t2, pt_cache)
    dot3_row!(cosτ, d1r, d1trans)
    dot3_row!(sinτ, d2r, d1trans)
    τ .= -1.0.*atan.(sinτ, cosτ)
end

x = rand(3,N)
d = rand(3,N-1)
norm3_col!(d)
l = zeros(1, N-1)
κ = zeros(2, N-2)
τ = zeros(1, N-2)
caches = (zeros(6, N-1), zeros(9, N-2))
@btime full_kinematics_col!(l, κ, τ, x, d, caches)

x = rand(N,3)
d = rand(N-1,3)
norm3_row!(d)
l = zeros(N-1)
κ = zeros(N-2, 2)
τ = zeros(N-2)
caches = (zeros(N-1,6), zeros(N-2,9))
@btime full_kinematics_row!(l, κ, τ, x, d, caches)

# ------------------------------------------------------------------------------
# Elastic Energy
# ------------------------------------------------------------------------------
# cache is (3 x ne) or larger
function elastic_energy_col(l0,κ0,τ0, l,κ,τ, k,B,β, cache)

    ne = length(l)
    nv = ne - 1
    cache_s = reshape(view(cache, 1:ne), 1, ne)
    cache_s .= k./l0.*(l .- l0).^2
    Es = 1/2*sum(cache_s)

    # Reference voronoi lengths
    vl0 = view(cache, 1:1, 1:nv)
    vl0[1] = l0[1] + l0[2]/2
    vl0[2:nv-1] .= (view(l0, 2:ne-2) .+ view(l0, 3:ne-1))./2
    vl0[nv] = l0[ne-1]/2 + l0[ne]

    cache_t = view(cache, 2:2, 1:nv)
    cache_t .= β./vl0.*(τ .- τ0).^2
    Et = 1/2*sum(cache_t)

    cache_b = view(cache, 2:3, 1:nv)
    cache_b .= B./vl0.*(κ .- κ0).^2
    Eb = 1/2*sum(cache_b)

    Es + Et + Eb
end

# cache is (ne x 3) or larger
function elastic_energy_row(l0,κ0,τ0, l,κ,τ, k,B,β, cache)

    ne = length(l)
    nv = ne - 1
    cache_s = view(cache, 1:ne)
    cache_s .= k./l0.*(l .- l0).^2
    Es = 1/2*sum(cache_s)

    # Reference voronoi lengths
    vl0 = view(cache, 1:nv, 1)
    vl0[1] = l0[1] + l0[2]/2
    vl0[2:nv-1] .= (view(l0, 2:ne-2) .+ view(l0, 3:ne-1))./2
    vl0[nv] = l0[ne-1]/2 + l0[ne]

    cache_t = view(cache, 1:nv, 2)
    cache_t .= β./vl0.*(τ .- τ0).^2
    Et = 1/2*sum(cache_t)

    cache_b = view(cache, 1:nv, 2:3)
    cache_b .= B./vl0.*(κ .- κ0).^2
    Eb = 1/2*sum(cache_b)

    Es + Et + Eb
end

l0 = rand(1, N-1)
κ0 = rand(2, N-2)
τ0 = rand(1, N-2)
l = rand(1, N-1)
κ = rand(2, N-2)
τ = rand(1, N-2)
k = rand()
B = rand(2,1)
β = rand()
cache = zeros(3,N-1)

@btime elastic_energy_col(l, κ, τ, l0, κ0, τ0, k ,B, β, cache)

l0 = rand(N-1)
κ0 = rand(N-2,2)
τ0 = rand(N-2)
l = rand(N-1)
κ = rand(N-2,2)
τ = rand(N-2)
k = rand()
B = rand(1, 2)
β = rand()
cache = zeros(N-1,3)

@btime elastic_energy_row(l, κ, τ, l0, κ0, τ0, k ,B, β, cache)

# ------------------------------------------------------------------------------
# Gravitational Energy
# ------------------------------------------------------------------------------
# cache is at least 1xne
function gravitational_energy_col(x, m, g, cache)
    nv = size(x,2)
    tmp = reshape(view(cache, 1, 1:nv-1), 1, nv-1)
    tmp .= (view(x,3:3,1:nv-1) .+ view(x,3:3,2:nv)).*m
    return sum(tmp)*g/2
end

# cache is at least 1xne
function gravitational_energy_row(x, m, g, cache)
    nv = size(x,1)
    tmp = view(cache, 1:nv-1)
    tmp .= (view(x, 1:nv-1, 3) .+ view(x, 2:nv, 3)).*m
    return sum(tmp)*g/2
end

x = rand(3, N)
m = rand(1, N-1)
g = 10.0
cache = zeros(1, N-1)
@btime gravitational_energy_col(x, m, g, cache)

x = rand(N, 3)
m = rand(N-1)
cache = zeros(N-1)
@btime gravitational_energy_row(x, m, g, cache)


text:
module ReducedimBenchmark
using DataArrays, Benchmark

# seed rng for more consistent timings
Random.seed!(1776)

const TEST_NAMES = [
    "Matrix",
    "DataMatrix No missing skipmissing=false",
    "DataMatrix No missing skipmissing=true",
    "DataMatrix Half missing skipmissing=false",
    "DataMatrix Half missing skipmissing=true"
]

function make_test_types(genfunc, sz)
    mat = genfunc(abs2(sz))
    na = shuffle!([trues(ifloor(abs2(sz)/2)), falses(iceil(abs2(sz)/2))])
    (
        reshape(mat, sz, sz),
        DataArray(reshape(mat, sz, sz)),
        DataArray(reshape(mat, sz, sz), reshape(na, sz, sz))
    )
end

const Data = make_test_types(rand, 10000)

macro perf(fn, dim, replications)
    quote
        println($fn, " (region = ", $dim, ")")
        fns = [()->$fn(Data[1], $dim),
               ()->$fn(Data[2], $dim),
               ()->$fn(Data[2], $dim; skipmissing=true),
               ()->$fn(Data[3], $dim),
               ()->$fn(Data[3], $dim; skipmissing=true)]
        gc_disable()
        df = compare(fns, $replications)
        gc_enable()
        gc()
        df[:Function] = TEST_NAMES
        df[:Relative] = df[:Average]./df[1, :Average]
        println(df)
    end
end

@perf sum 1 10
@perf sum 2 10
@perf maximum 1 10
@perf maximum 2 10
@perf mean 1 10
@perf mean 2 10
@perf var 1 10
@perf var 2 10
end