text:
export phsp_iterator
export iaea_iterator
export iaea_writer
export egs_iterator
export egs_writer
export egs_write
export iaea_write
export phsp_write

function guess_format_from_path(path::AbstractString)
    stem, ext = splitext(path)
    if startswith(ext, ".IAEA")
        return FormatIAEA()
    elseif startswith(ext, ".egsphsp")
        return FormatEGS()
    else
        msg = "Cannot guess format from file extension for $path."
        throw(ArgumentError(msg))
    end
end

function iaea_writer end

iaea_iterator(path) = iaea_iterator(IAEAPath(path))

function iaea_iterator(path::IAEAPath)
    @argcheck ispath(path.header)
    @argcheck ispath(path.phsp)
    h = load(path.header, IAEAHeader)
    io = FastReadIO(open(path.phsp))
    IAEAPhspIterator(io,h)
end

function egs_iterator(path::AbstractString)
    io = open(path, "r")
    egs_iterator(io)
end

phsp_iterator(path::IAEAPath) = iaea_iterator(path)
function phsp_iterator(path::AbstractString)
    fmt = guess_format_from_path(path)
    phsp_iterator(path, fmt)
end
phsp_iterator(path::AbstractString, ::FormatEGS)  = egs_iterator(path)
phsp_iterator(path::AbstractString, ::FormatIAEA) = iaea_iterator(path)

for xxx_iterator in [:egs_iterator, :iaea_iterator, :phsp_iterator]
    @eval function $(xxx_iterator)(f, path)
        iter = $(xxx_iterator)(path)
        ret = call_fenced(f, iter)
        close(iter)
        ret
    end
end

function header_type(::Type{IAEAParticle{Nf, Ni}}) where {Nf, Ni}
    IAEAHeader{Nf, Ni}
end

function similar_header(ps, args...)
    P = eltype(ps)
    H = header_type(P, args...)
    H()
end

function _iaea_header_like(P::Type{IAEAParticle{Nf, Ni}}) where {Nf, Ni}
    r = IAEAHeader{Nf, Ni}()
end

function iaea_write(path, ps)
    r = _iaea_header_like(eltype(ps))
    iaea_writer(path,r) do w
        for p in ps
            write(w, p)
        end
    end
end

function egs_write(path, ps)
    P = eltype(ps)
    egs_writer(path, P) do w
        for p in ps
            write(w, p)
        end
    end
end

function phsp_write(path, ps)
    P = eltype(ps)
    if P <: IAEAParticle
        iaea_write(path, ps)
    elseif P <: EGSParticle
        egs_write(path, ps)
    else
        msg = "Cannot guess format from eltype(ps) = $P."
        throw(ArgumentError(msg))
    end
end


text:
# neo_import_dir = "/Users/cameronprybol/Library/Application Support/Neo4j Desktop/Application/relate-data/dbmss/dbms-8ab8baac-5dea-4137-bb24-e0b426447940/import"

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Description

```jldoctest
julia> 1 + 1
2
```
"""
function cypher(cmd;
    address="neo4j://localhost:7687",
    username="neo4j",
    password="password",
    format="auto",
    database="system"
    )
    cmd = `cypher-shell --address $(address) --username $(username) --password $(password) --format $(format) --database $(database) $(cmd)`
#    cmd = `/home/jupyter-cjprybol/software/cypher-shell/cypher-shell --address $(address) --username $(username) --password $(password) --format $(format) --database $(database) $(cmd)`
    return cmd
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Description

```jldoctest
julia> 1 + 1
2
```
"""
function load_ncbi_metadata(db)
    if !(db in ["genbank", "refseq"])
        error()
    end
    ncbi_summary_url = "https://ftp.ncbi.nih.gov/genomes/$(db)/assembly_summary_$(db).txt"
    buffer = IOBuffer(HTTP.get(ncbi_summary_url).body)
    ncbi_summary_table = DataFrames.DataFrame(uCSV.read(buffer, comment = "#  ", header=1, delim='\t')...)
    return ncbi_summary_table
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Description

```jldoctest
julia> 1 + 1
2
```
"""
function taxonomic_id_to_children(tax_id)
    # NEO4J_BIN_DIR = "/home/jupyter-cjprybol/software/neo4j-community-4.4.3/bin"
    # if !occursin(NEO4J_BIN_DIR, ENV["PATH"])
    #     ENV["PATH"] = "$(NEO4J_BIN_DIR):" * ENV["PATH"]
    # end
    DOMAIN = "ncbi-taxonomy.cjp.garden"
    # DOMAIN="34.138.169.173"
    USERNAME="neo4j"
    PASSWORD=readline(joinpath(homedir(), ".config", "neo4j", "ncbi-taxonomy.password.txt"));
    DATABASE = "neo4j"
    # ADDRESS="neo4j://$(DOMAIN):7687"
    ADDRESS="bolt://$(DOMAIN):7687"

    # NOTE! *, or 0 distance (e.g. [*0..2]) step range will include source node!!!!
    cmd = "MATCH (n)<-[*]-(n2) WHERE n.tax_id IS NOT NULL AND n.tax_id = \"$(tax_id)\" RETURN DISTINCT n2.tax_id AS tax_id"
    println(cmd)

    cypher = Mycelia.cypher(cmd, address=ADDRESS, username = USERNAME, password = PASSWORD, database = DATABASE)
    tax_ids = readlines(open(cypher))[2:end]
    tax_ids = strip.(tax_ids, '"')
    tax_ids = parse.(Int, tax_ids)
    return unique(tax_ids)
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Description

```jldoctest
julia> 1 + 1
2
```
"""
function load_refseq_metadata()
    return load_ncbi_metadata("refseq")
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Description

```jldoctest
julia> 1 + 1
2
```
"""
function load_genbank_metadata()
    return load_ncbi_metadata("genbank")
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Description

```jldoctest
julia> 1 + 1
2
```
"""
function ncbi_ftp_path_to_url(ftp_path, extension)
    # genomic.fna.gz
    # genomic.gff.gz
    # protein.faa.gz

    # assembly_report.txt
    # assembly_stats.txt
    # cds_from_genomic.fna.gz
    # feature_count.txt.gz
    # feature_table.txt.gz
    # genomic.gbff.gz
    # genomic.gtf.gz
    # protein.gpff.gz
    # translated_cds.faa.gz
    f_name = basename(ftp_path) * "_" * extension
    new_path = joinpath(ftp_path, f_name)
    return new_path
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Description

```jldoctest
julia> 1 + 1
2
```
"""
function countmap_columns(table)
    for n in names(refseq_metadata)
        display(n)
        display(StatsBase.countmap(refseq_metadata[!, n]))
    end
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Get dna (db = "nuccore") or protein (db = "protein") sequences from NCBI
or get fasta directly from FTP site

```jldoctest
julia> 1 + 1
2
```
"""
function get_sequence(;db=""::String, accession=""::String, ftp=""::String)
    if !isempty(db) && !isempty(accession)
        # API will block if we request more than 3 times per second, so set a 1/2 second sleep to set max of 2 requests per second when looping
        sleep(0.5)
        url = "https://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?db=$(db)&report=fasta&id=$(accession)"
        return FASTX.FASTA.Reader(IOBuffer(HTTP.get(url).body))
    elseif !isempty(ftp)
        return FASTX.FASTA.Reader(CodecZlib.GzipDecompressorStream(IOBuffer(HTTP.get(ftp).body)))
    else
        @error "invalid call"
    end
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Get dna (db = "nuccore") or protein (db = "protein") sequences from NCBI
or get fasta directly from FTP site

```jldoctest
julia> 1 + 1
2
```
"""
function get_gff(;db=""::String, accession=""::String, ftp=""::String)
    if !isempty(db) && !isempty(accession)
        # API will block if we request more than 3 times per second, so set a 1/2 second sleep to set max of 2 requests per second when looping
        sleep(0.5)
        url = "https://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?db=$(db)&report=gff3&id=$(accession)"
        return GFF3.Reader(IOBuffer(HTTP.get(url).body))
    elseif !isempty(ftp)
        return GFF3.Reader(CodecZlib.GzipDecompressorStream(IOBuffer(HTTP.get(ftp).body)))
    else
        @error "invalid call"
    end
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Get dna (db = "nuccore") or protein (db = "protein") sequences from NCBI
or get fasta directly from FTP site

```jldoctest
julia> 1 + 1
2
```
"""
function get_genbank(;db=""::String, accession=""::String, ftp=""::String)
    if !isempty(db) && !isempty(accession)
        # API will block if we request more than 3 times per second, so set a 1/2 second sleep to set max of 2 requests per second when looping
        sleep(0.5)
        url = "https://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?db=$(db)&report=genbank&id=$(accession)"
        return GenomicAnnotations.readgbk(IOBuffer(HTTP.get(url).body))
    elseif !isempty(ftp)

        return GenomicAnnotations.readgbk(CodecZlib.GzipDecompressorStream(IOBuffer(HTTP.get(ftp).body)))
    else
        @error "invalid call"
    end
end

ext: jl

max_stars_repo_name: cjprybol/Eisenia

max_issues_repo_name: cjprybol/Eisenia

max_forks_repo_name: cjprybol/Eisenia


text:
ENV["PYTHON"]=""
using Pkg
Pkg.add("Conda")
using Conda
Conda.update()
Conda.add("matplotlib")
Conda.add_channel("conda-forge")
Pkg.add("PyCall")
Pkg.build("PyCall")
Pkg.add("PyPlot")
Pkg.add("SpecialFunctions")
Pkg.add("ForwardDiff")  # 6/30/2020
Pkg.add("DiffResults")  #
Pkg.add(PackageSpec(name="Optim",version="0.22.0"))
Pkg.pin(PackageSpec(name="Optim",version="0.22.0"))
Pkg.add("Documenter")
Pkg.add("GSL")
Pkg.add("DelimitedFiles")
Pkg.add("IterativeSolvers") # 4/6/2020
Pkg.add("JLD2")
Pkg.add("MCMCDiagnostics")
if VERSION <= v"1.0.1"
  Pkg.clone("https://github.com/rodluger/Limbdark.jl.git")
else
  Pkg.add(PackageSpec(url="https://github.com/rodluger/Limbdark.jl.git"))
end


text: using Test, PreSumm


text:
module TestCompress
using Test
using UrlDownload
using DataFrames

@testset "compress" begin
    df = DataFrame(x = [1, 3], y = [2, 4])
    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/test.xz"
    res = urldownload(url) |> DataFrame
    @test res == df
    res = urldownload(url, compress = :xz) |> DataFrame
    @test res == df

    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/test.gz"
    res = urldownload(url) |> DataFrame
    @test res == df
    res = urldownload(url, compress = :gzip) |> DataFrame
    @test res == df

    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/test.bz2"
    res = urldownload(url) |> DataFrame
    @test res == df
    res = urldownload(url, compress = :bzip2) |> DataFrame
    @test res == df

    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/test.lz4"
    res = urldownload(url) |> DataFrame
    @test res == df
    res = urldownload(url, compress = :lz4) |> DataFrame
    @test res == df

    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/test.zst"
    res = urldownload(url) |> DataFrame
    @test res == df
    res = urldownload(url, compress = :zstd) |> DataFrame
    @test res == df
end

@testset "zip compress" begin
    df = DataFrame(x = [1, 3], y = [2, 4])
    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/ext.csv.zip"
    res = urldownload(url) |> DataFrame
    @test res == df

    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/test2.zip"
    @test_logs (:warn, "More than one file in zip archive, returning first.") urldownload(url)

    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/test2.zip"
    @test_logs (:warn, "Data format unknown is not supported.") urldownload(url, multifiles = true)
end

@testset "compress overrides" begin
    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/test.zst"
    urldownload(url, parser = identity)

    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/test.zst"
    urldownload(url, compress = :none, parser = identity)

    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/test2.zip"
    @test_logs (:warn, "More than one file in zip archive, returning first.") urldownload(url, parser = identity)

    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/test2.zip"
    urldownload(url, compress = :none, parser = identity)
end

end # module


text:
# # Dispatch Analysis
#
# ```@meta
# CurrentModule = JET
# ```

# ## Motivation
#
# When Julia compiles your code and type inference on it was not so successful, the compiler
# is likely to be unable to determine which method should be called at each generic function
# callsite, and then it will be looked up at runtime.
# That is called "runtime dispatch", which is known as a common source of performance problem
# – the compiler can't do various kinds of optimizations including inlining when it can't
# determine a single matching method, and method lookup itself can also be a bottleneck if
# the call happens many times.
#
# In order to avoid this problem, we usually use `code_typed`, inspect its output, and check
# if there is anywhere type is not well inferred (i.e. where is "type-instable") and
# optimization was not successful.
# But the problem is that `code_typed` can only present the "final" output of inference or
# optimization, and we can't inspect an entire call graph and _may_ not be able to find where
# a problem happened and how the "type instability" has been propagated.
#
# There is a nice package called [Cthulhu.jl](https://github.com/JuliaDebug/Cthulhu.jl),
# which allows us to inspect the output of `code_typed` by _descending_ into a call tree,
# recursively and interactively.
# The workflow with Cthulhu is much more powerful, but still, it's tedious.
#
# So, why not automate it ?
# We can use JET's pluggable analysis framework and create such an analyzer that automatically
# analyzes your code and alarms you when it detects anywhere Julia can't determine matching
# method statically and thus runtime dispatch will happen at runtime.

# ## Implementation
#
# In this analysis, the analyzer will be designed to detect:
# 1. where Julia compiler gives up optimization
# 2. where a runtime dispatch will happen
#
# The case 1. will happen when there are (mutually) recursive calls and Julia compiler
# decided not to do inference in order to make sure the inference's termination.
# In such a case, optimization won't happen and method dispatches aren't resolved statically,
# so we will just report it (as `OptimizationFailureReport`). In order to detect the case 2.,
# we will inspect the optimized IR and look for `:call` expressions. `:call` expressions are
# such calls that were not resolved statically and will be dispatched at runtime (as opposed
# to [`:invoke`](https://docs.julialang.org/en/v1/devdocs/ast/#Expr-types) expressions, that
# represent staticall resolved generic function calls).
#
# We will define `DispatchAnalyzer <: AbstractAnalyzer`, and overload some of `Core.Compiler` methods with it:
# - `Core.Compiler.finish(frame::CC.InferenceState, analyzer::DispatchAnalyzer)` to check if optimization will happen or not (the case 1.)
# - `Core.Compiler.finish!(analyzer::DispatchAnalyzer, caller::CC.InferenceResult)` to inspect an optimized IR (the case 2.)

using JET.JETInterface
const CC = Core.Compiler
import JET:
    JET,
    @invoke,
    get_source,
    @isexpr

struct DispatchAnalyzer{T} <: AbstractAnalyzer
    state::AnalyzerState
    opts::BitVector
    frame_filter::T
    __cache_key::UInt
end
function DispatchAnalyzer(;
    ## a predicate, which takes `CC.InfernceState` and returns whether we want to analyze the call or not
    frame_filter = x::CC.InferenceState->true,
    jetconfigs...)
    state = AnalyzerState(; jetconfigs...)
    ## we want to run different analysis with a different filter, so include its hash into the cache key
    cache_key = state.param_key
    cache_key = hash(frame_filter, cache_key)
    return DispatchAnalyzer(state, BitVector(), frame_filter, cache_key)
end

## maybe filtering by module would be most useful
function module_filter(m)
    return function (frame::CC.InferenceState)
        frame.mod === m
    end
end

## AbstractAnalyzer API requirements
JETInterface.AnalyzerState(analyzer::DispatchAnalyzer)                          = analyzer.state
JETInterface.AbstractAnalyzer(analyzer::DispatchAnalyzer, state::AnalyzerState) = DispatchAnalyzer(state, analyzer.opts, analyzer.frame_filter, analyzer.__cache_key)
JETInterface.ReportPass(analyzer::DispatchAnalyzer)                             = DispatchAnalysisPass()
JETInterface.get_cache_key(analyzer::DispatchAnalyzer)                          = analyzer.__cache_key

struct DispatchAnalysisPass <: ReportPass end

## ignore all reports defined by JET, since we'll just define our own reports
(::DispatchAnalysisPass)(T::Type{<:InferenceErrorReport}, @nospecialize(_...)) = return

function CC._typeinf(analyzer::DispatchAnalyzer, frame::CC.InferenceState)
    @assert isempty(analyzer.opts)
    ret = @invoke CC._typeinf(analyzer::AbstractAnalyzer, frame::CC.InferenceState)
    @assert isempty(analyzer.opts)
    return ret
end

function CC.finish(frame::CC.InferenceState, analyzer::DispatchAnalyzer)
    ret = @invoke CC.finish(frame::CC.InferenceState, analyzer::AbstractAnalyzer)

    if !analyzer.frame_filter(frame)
        push!(analyzer.opts, false)
    else
        if isa(get_source(frame.result), CC.OptimizationState)
            push!(analyzer.opts, true)
        else # means, compiler decides not to do optimization
            ReportPass(analyzer)(OptimizationFailureReport, analyzer, frame.result)
            push!(analyzer.opts, false)
        end
    end

    return ret
end

@reportdef struct OptimizationFailureReport <: InferenceErrorReport end
JETInterface.get_msg(::Type{OptimizationFailureReport}, args...) =
    return "failed to optimize" #: signature of this MethodInstance

function (::DispatchAnalysisPass)(::Type{OptimizationFailureReport}, analyzer::DispatchAnalyzer, result::CC.InferenceResult)
    add_new_report!(result, OptimizationFailureReport(result.linfo))
end

function CC.finish!(analyzer::DispatchAnalyzer, frame::CC.InferenceState)
    caller = frame.result

    ## get the source before running `finish!` to keep the reference to `OptimizationState`
    src = get_source(caller)

    ## run `finish!(::AbstractAnalyzer, ::CC.InferenceState)` first to convert the optimized `IRCode` into optimized `CodeInfo`
    ret = @invoke CC.finish!(analyzer::AbstractAnalyzer, frame::CC.InferenceState)

    if popfirst!(analyzer.opts) # optimization happened
        if isa(src, Core.Const) # the optimization was very successful, nothing to report
        elseif isa(src, CC.OptimizationState) # the compiler cached the optimized IR, just analyze it
            ReportPass(analyzer)(RuntimeDispatchReport, analyzer, caller, src)
        else
            ## we should already report `OptimizationFailureReport` for this case,
            ## and thus this pass should never happen
            throw("got $src, unexpected source found")
        end
    end

    return ret
end

@reportdef struct RuntimeDispatchReport <: InferenceErrorReport end
JETInterface.get_msg(::Type{RuntimeDispatchReport}, _) =
    return "runtime dispatch detected" #: call signature

function (::DispatchAnalysisPass)(::Type{RuntimeDispatchReport}, analyzer::DispatchAnalyzer, caller::CC.InferenceResult, opt::CC.OptimizationState)
    (; sptypes, slottypes) = opt
    for (pc, x) in enumerate(opt.src.code)
        if @isexpr(x, :call)
            ft = CC.widenconst(CC.argextype(first(x.args), opt.src, sptypes, slottypes))
            ft <: Core.Builtin && continue # ignore `:call`s of language intrinsics
            add_new_report!(caller, RuntimeDispatchReport((opt, pc)))
        end
    end
end

# ## Usages
#
# So we defined our analyzer.
# Let's setup utility analysis entries first:

using JET, InteractiveUtils # to use analysis entry points

function report_dispatch(@nospecialize(f), @nospecialize(types = Tuple{});
                         analyzer = DispatchAnalyzer,
                         jetconfigs...)
    @assert analyzer === DispatchAnalyzer "analyzer is fixed to $DispatchAnalyzer"
    report_call(f, types; analyzer, jetconfigs...)
end
macro report_dispatch(ex0...)
    return InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :report_dispatch, ex0)
end

# Now we can just call `@report_dispatch f(args...)` and check if there are any problematic
# part within the entire call tree of `f(args...)`.

# ### Simple cases
#
# First, let's play with simple and factitious examples and check if `DispatchAnalyzer`
# works as expected.

f(a) = a
f(a::Number) = a

# `f(::Int)` a concrete call and just type stable and anything shouldn't be reported:
@report_dispatch f(10) # should be ok

# But if the argument type isn't well typed, compiler can't determine which method to call,
# and it will lead to runtime dispatch:
report_dispatch((Any,)) do a
    f(a) # runtime dispatch !
end

# Note that even if a call is not "well-typed", i.e. it's not a concrete call, runtime
# dispatch won't happen as far as a single method can be resovled statically:
report_dispatch((Integer,)) do a
    f(a) # this isn't so good, but ok
end

# Ok, working nicely so far. Let's move on to a bit more complicated examples.
# If we annotate `@noinline` to a function, then its call won't be inlined and will be
# dispatched runtime. We will confirm this:

@inline   g1(a) = return a # WARNING use `@inline` just for demonstrative purpose, usually we really don't need to use it for a very trivial function like this
@noinline g2(a) = return a
report_dispatch((Any,)) do a
    r1 = g1(a) # this call should be statically resolved and inlined
    r2 = g2(a) # this call should be statically resolved but not inlined, and will be dispatched
    return (r1, r2)
end

# We can assert this report by looking at the output of `code_typed`, where `g1(a)` has been
# just disappeared by inlining and optimizations, but `g2(a)` still remains as a `:call` expression:
code_typed((Any,)) do a
    r1 = g1(a) # this call should be statically resolved and inlined
    r2 = g2(a) # this call should be statically resolved but not inlined, and will be dispatched
    return (r1, r2)
end |> first

# ### Real-world targets
#
# Let's run `DispatchAnalyzer` on real-world code and check how it works.
# Here we will test with Julia's `Base` module.

# Numerical computations are usually written to be very type-stable, and so we expects e.g.
# `sin(10)` to be dispatch-free and run fast. Let's check it.
@report_dispatch sin(10)

# Oh no, so runtime dispatch happens there even in `Base`. Well, actually, this specific dispatch
# is expected. Especially, <https://github.com/JuliaLang/julia/pull/35982> implements an
# heuristic to intentionally disable inference (and so succeeding optimizations too) in
# order to ease [the latency problem, a.k.a. "first-time-to-plot"](https://julialang.org/blog/2020/08/invalidations/).
# The report trace certainly suggests a dispatch was detected where `DomainError` can be thrown.
# We can turn off the heuristic by turning off [the `unoptimize_throw_blocks::Bool` configuration](@ref abstractinterpret-config),
# and this time any runtime dispatch won't be reported:
@report_dispatch unoptimize_throw_blocks=false sin(10)

# We can also confirm that the same thing would happen for `rand(1:1000)`:
@report_dispatch rand(1:1000) # an runtiem dispatch will be detected within a `throw` block
#
@report_dispatch unoptimize_throw_blocks=false rand(1:1000) # nothing should be reported

# Finally, let's see an example of very "type-instable" code maintained within `Base`.
# Typically, anything involved with I/O is written in a very dynamic way for good reasons,
# and certainly, e.g. `println(QuoteNode(nothing))` will yield bunch of runtime dispatches:
# ```julia
# @report_dispatch println(QuoteNode(nothing))
# ═════ 585 possible errors found ═════
# ┌ @ coreio.jl:4 Base.println(Core.tuple(Core.typeassert(Base.stdout, Base.IO)), xs...)
# │┌ @ strings/io.jl:73 Base.print(Core.tuple(io), xs, Core.tuple("\n")...)
# ││┌ @ strings/io.jl:43 Base.lock(io)
# │││┌ @ show.jl:334 Base.lock(Base.getproperty(io, :io))
# # so many reports follow ...
# ```

# But what if your code contains a single `println` call, which you're absolutely okay with
# the type instabilities involved with it (e.g. it's only called once or only in debug mode,
# or such), but still you want to assert that any other part of code is type-stable and dispatch-free ?
# `DispatchAnalyzer`'s `frame_filter` option can be useful for this, by allowing us to
# specificy where it should and shouldn't run analysis.
# For example, we can check type-stabilities of anything in the current module like this:

## problem: when ∑1/n exceeds 30 ?
function compute(x)
    r = 1
    s = 0.0
    n = 1
    @time while r < x
        s += 1/n
        if s ≥ r
            println("round $r/$x has been finished") # we're not interested type-instabilities within this call
            r += 1
        end
        n += 1
    end
    return n, s
end

## NOTE:
## `compute(30)` will take more than hours in actual execution, according to https://twitter.com/genkuroki/status/1401332946707963909,
## but `@report_dispatch` will just do abstract interpretation of the call, so will finish instantly
@report_dispatch frame_filter=module_filter(@__MODULE__) compute(30)


text:
@withkw immutable HHParameter
    Cm::Float = 1uF * cm^(-2) * 20000um^2
    gl::Float = 5e-5siemens * cm^(-2) * 20000um^2
    El::Float = -65mV
    Ek::Float = -90mV
    En::Float = 50mV
    gn::Float = 100msiemens * cm^(-2) * 20000um^2
    gk::Float = 30msiemens * cm^(-2) * 20000um^2
    Vt::Float = -63mV
    τe::Float = 5ms
    τi::Float = 10ms
    Ee::Float = 0mV
    Ei::Float = -80mV
end

@withkw type HH
    param::HHParameter = HHParameter()
    N::Int = 100
    v::Vector{Float} = param.El + 5(randn(N) - 1)
    m::Vector{Float} = zeros(N)
    n::Vector{Float} = zeros(N)
    h::Vector{Float} = ones(N)
    ge::Vector{Float}  = (1.5randn(N) + 4) * 10nS
    gi::Vector{Float}  = (12randn(N) + 20) * 10nS
    fire::Vector{Bool} = zeros(Bool, N)
    I::Vector{Float} = zeros(N)
    records::Dict = Dict()
end

@replace function integrate!(p::HH, param::HHParameter, dt::Float)
    @inbounds for i = 1:N
        m[i] += dt * (0.32f0 * (13f0 - v[i] + Vt) / (exp((13f0 - v[i] + Vt) / 4f0) - 1f0) * (1f0 - m[i]) -
        0.28f0 * (v[i] - Vt - 40f0) / (exp((v[i] - Vt - 40f0) / 5f0) - 1f0) * m[i])
        n[i] += dt * (0.032f0 * (15f0 - v[i] + Vt) / (exp((15f0 - v[i] + Vt) / 5f0) - 1f0) * (1f0 - n[i]) -
        0.5f0 * exp((10f0 - v[i] + Vt) / 40f0) * n[i])
        h[i] += dt * (0.128f0 * exp((17f0 - v[i] + Vt) / 18f0) * (1f0 - h[i]) -
        4f0 / (1f0 + exp((40f0 - v[i] + Vt) / 5f0)) * h[i])
        v[i] += dt / Cm * ( I[i] + gl * (El - v[i]) + ge[i] * (Ee - v[i]) + gi[i] * (Ei - v[i]) +
        gn * m[i]^3 * h[i] * (En - v[i]) + gk * n[i]^4 * (Ek - v[i]) )
        ge[i] += dt * -ge[i] / τe
        gi[i] += dt * -gi[i] / τi
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > -20f0
    end
end


text:
include("diffusion_mnist.jl")

function test()
    @info "Begin tests for diffusion_mnist.jl"

    gaussfourierproj_test = GaussianFourierProjection(32, 20.0f0)
    # GaussianFourierProjection(embed_dim, ⋅)(batch) => (embed_dim, batch)
    @assert gaussfourierproj_test(randn(Float32, 32)) |> size == (32, 32)
    # W is fixed wrt. repeated calls
    @assert gaussfourierproj_test(
        ones(Float32, 32)) ==
            gaussfourierproj_test(ones(Float32, 32)
    )
    # W is not trainable
    @assert params(gaussfourierproj_test) == Flux.Params([])

    @assert expand_dims(ones(Float32, 32), 3) |> size == (1, 1, 1, 32)

    unet_test = UNet()
    x_test = randn(Float32, (28, 28, 1, 32))
    t_test = rand(Float32, 32)
    score_test = unet_test(x_test, t_test)
    @assert score_test |> size == (28, 28, 1, 32)
    @assert typeof(score_test) == Array{Float32,4}

    # Test gradient computation
    grad_test = gradient(
        () -> model_loss(unet_test, x_test), params(unet_test)
    )
    @assert grad_test.params == params(unet_test)

    train(save_path="test", epochs=1, batch_size=4096, tblogger=false)

    @info "Tests complete for diffusion_mnist.jl"
end

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


text:
using DataFrames, Statistics, DelimitedFiles, CSV, LinearAlgebra, Random, Distributions
using Pkg
Pkg.update()
#pseudo data generation 5x5 matrix
b = randn(5)
X0 = zeros(5)
A = [1 2 0 1 3; 0 1 2 1 0; 1 2 1 2 2; 1 4 2 2 1; 2 1 2 0 1]

#jacobi solver
# """function for jacobi solver
# b: is the response vector
# A: is the design matrix for intercept + slope
# iter: is the number of the iterations in the jacobi
# Tol: Tolerance
# X0: Initial guess for k variables
# """
function jacobi(A,b,N,x,TOL)
    """Solves the equation Ax=b via the Jacobi iterative method."""
    # Create an initial guess if needed
    # if x
    #     x = zeros(size(A, 2))
    # end
    # Create a vector of the diagonal elements of A
    # and subtract them from A
    D = diag(A)
    R = A .- D
    println(size(R))

    # Iterate for N times
    for i in 1:N
        X0 = x
        x = inv(Diagonal(D))*(b .- R*x)
        if abs.((x .- X0)'*(x .- X0)) .< TOL
            return round.(x, digits=2)
        end
    end
    return round.(x, digits=2)
end

#small genotype matrix 14000 x 100
geno = readdlm("../SSBR-JWAS-Implementation\\SSBR-vs-BR-Example\\data\\Input-Data\\Geno_Small.txt")
println(size(geno))
ped = readdlm("../SSBR-JWAS-Implementation\\SSBR-vs-BR-Example\\data\\Input-Data\\Pedigree_XSim_Raw1.txt")

#get A matrix
s=ped[:,2]
println(s)
d=ped[:,3]
n = length(s)
s=(s .== 0)*n .+s
d=(d .== 0)*n .+d;
A = zeros(n,n);
println(s[2])
for i in 1:n
    A[i,i] = 1 + A[Int64(s[i]), Int64(d[i])]/2
    for j in (i+1):n
        A[i,j] = ( A[i, Int64(s[j])] + A[i, Int64(d[j])] ) / 2
        A[j,i] = A[i,j]
    end
end
println(A[1, 1])

#Define genetic and residual variances
vara = Diagonal(repeat([2.0], size(geno, 2)))
vare = Diagonal(repeat([1.0], size(geno, 1)))

# Define random phenotype of animals
y = randn(size(geno, 1))

#Construct MME
lhs = (geno'*inv(vare)*geno + inv(vara))
rhs = geno'*inv(vare)*y

#Solve MME using jacobi solver
X0 = randn(size(geno, 2))
alpha = jacobi(lhs, rhs, 50, X0, 0.00002)
println(alpha[1])
#Get y pred and compare to y observed
y_pred = geno*alpha .+ diag(vare)

#Condition of MME matrix
cond_mme = cond(lhs)
#calculate correlation between y_pred and y_obs
accuracy = cor(y, y_pred)


text:
@testset "extreme_filter" begin
    @testset "interface" begin
        img = rand(1:10, 7, 7)

        out = extreme_filter(max, img)
        @test out == extreme_filter(max, img; dims=(1, 2), r=1)

        se = strel_box(img, (1, 2); r=(1, 1))
        @test out == extreme_filter(max, img, se)

        se = strel_box(img, (1, 2); r=(1, 3))
        @test extreme_filter(max, img; dims=(1, 2), r=(1, 3)) == extreme_filter(max, img, se)

        out = similar(img)
        rst = extreme_filter!(max, out, img)
        @test rst === out
        @test out == extreme_filter!(max, out, img)

        out = similar(img)
        rst = extreme_filter!(max, out, img; dims=(1, 2), r=1)
        @test rst === out
        se = strel_box(img, (1, 2); r=(1, 1))
        ref_out = similar(img)
        extreme_filter!(max, ref_out, img, se)
        @test ref_out == out

        out = similar(img)
        rst = extreme_filter!(max, out, img; dims=(1, 2), r=(1, 3))
        @test rst === out
        se = strel_box(img, (1, 2); r=(1, 3))
        ref_out = similar(img)
        extreme_filter!(max, ref_out, img, se)
        @test ref_out == out
    end

    @testset "numerical" begin
        # Bool
        img = fill(false, 5, 5)
        img[3, 3] = 1
        ref = collect(strel_box((5, 5); r=1))

        img_d = @inferred extreme_filter(max, img)
        @test img_d == ref
        @test ref == .!extreme_filter(min, .!img) # dual property

        # Int
        A = Int[4 6 5 3 4; 8 6 9 4 8; 7 8 4 9 6; 6 2 2 1 7; 1 6 5 2 6]
        ref = [8 6 9 5 8; 8 9 9 9 8; 8 8 9 9 9; 7 8 5 9 7; 6 6 6 6 7]
        Ad = @inferred extreme_filter(max, A, strel_diamond(A))
        @test Ad == ref
        ref = [8 9 9 9 8; 8 9 9 9 9; 8 9 9 9 9; 8 8 9 9 9; 6 6 6 7 7]
        Ad = @inferred extreme_filter(max, A)
        @test Ad == ref

        # Gray{Float32}
        img = Gray{Float32}.(A ./ 9)
        img_d = @inferred extreme_filter(max, img)
        @test ref == gray.(img_d) .* 9
    end

    @testset "strel" begin
        # ensure it supports both strel representations
        A = rand(32, 32)
        se_mask = centered(collect(strel_diamond(A)))
        out = extreme_filter(max, A, se_mask)
        se_offsets = strel(CartesianIndex, se_mask)
        @test out == extreme_filter(max, A, se_offsets)
    end

    @testset "offset arrays" begin
        img = centered(rand(1:10, 32, 32))
        out = extreme_filter(max, img)
        @test out == centered(extreme_filter(max, OffsetArrays.no_offset_view(img)))
    end

    @testset "multi channel image" begin
        # supporting color image requires some reduced ordering function
        img = rand(RGB, 32, 32)
        msg = "function `max` is not a well-defined select function on type `RGB{Float64}`: does `f(x::T, y::T)` work as expected?"
        @test_throws ArgumentError(msg) extreme_filter(max, img)

        # the test functions below don't make much sense in practice, but it's good to test them
        _to_number(c::AbstractRGB) = red(c)
        _to_number(x::Number) = x
        _select(x, y) = _to_number(x) > _to_number(y) ? x : y
        out = extreme_filter(_select, img)
        @test issubset(Set(unique(out)), Set(unique(img)))

        _select_new(x, y) = max(_to_number(x), _to_number(y))
        out = extreme_filter(_select_new, img)
        @test eltype(out) == eltype(img)
        @test sum(abs, channelview(RGB.(Gray.(out))) - channelview(out)) < 1e-4
    end

    @testset "optimization: diamond" begin
        # ensure the optimized implementation work equivalently to the generic fallback implementation
        for T in Any[Bool, Int, N0f8, Gray{N0f8}, Gray{Float64}, Float64]
            for N in (1, 2, 3)
                sz = ntuple(_ -> 32, N)
                img = T == Int ? rand(1:10, sz...) : rand(T, sz...)
                for r in (1, 3)
                    dims_list = ntuple(i -> ntuple(identity, i), N)
                    for dims in dims_list
                        se = strel_diamond(ntuple(_ -> 2r + 1, N), dims)
                        ref = ImageMorphology._extreme_filter_generic!(max, similar(img), img, se)
                        out = extreme_filter(max, img, se)
                        @test out == ref

                        imgc = centered(img)
                        outc = extreme_filter(max, imgc, se)
                        @test outc == centered(out)
                    end
                end
            end
        end
    end

    @testset "optimization: bool" begin
        # ensure the optimized implementation work equivalently to the generic fallback implementation
        for N in (1, 2, 3)
            sz = ntuple(_ -> 32, N)
            img = rand(Bool, sz...)
            for r in (1, 3)
                dims_list = ntuple(i -> ntuple(identity, i), N)
                for dims in dims_list, select in (max, min)
                    se = strel_diamond(ntuple(_ -> 2r + 1, N), dims)
                    ref = ImageMorphology._extreme_filter_generic!(select, similar(img), img, se)
                    out = extreme_filter(select, img, se) # SEDiamond method
                    @test out == ref

                    se = centered(collect(se))
                    out = extreme_filter(select, img, se) # MorphologySE method
                    @test out == ref
                end
            end
        end
    end
end


text:
"""
    UltimateMortality(vector; start_age=0)

Given a vector of rates, returns an `OffsetArray` that is indexed by attained age.

Give the optional keyword argument to start the indexing at an age other than zero.

# Examples
```julia-repl
julia> m = UltimateMortality([0.1,0.3,0.6,1]);

julia> m[0]
0.1

julia> m = UltimateMortality([0.1,0.3,0.6,1], start_age = 18);

julia> m[18]
0.1

```
"""
function UltimateMortality(v::Array{<:Real,1}; start_age=0)
    return OffsetArray(v, start_age - 1)
end

"""
    SelectMortality(select, ultimate; start_age=0)

Given a matrix rates, where the first row represents the select rates for a risk, will create a an `OffsetArray` that is indexed by issue age, containing a vector of rate indexed by attained age. The ultimate mortality vector is used for rates in the post-select period.

Give the optional keyword argument to start the indexing at an age other than zero.

# Examples
```
julia> ult = UltimateMortality([x / 100 for x in 0:100]);

julia> matrix = rand(50,10); # represents random(!) mortality rates with a select period of 10 years

julia> sel = SelectMortality(matrix,ult,start_age=0);

julia> sel[0] # the mortality vector for a select life with issue age 0
 0.12858960119349439
 0.1172480189376135
 0.8237661916705163
 ⋮
 0.98
 0.99
 1.0

julia> sel[0][95] # the mortality rate for a life age 95, that was issued at age 0
0.95
```
"""
function SelectMortality(select, ultimate; start_age=0)

    # iterate down the rows (issue ages)
    vs = map(enumerate(eachrow(select))) do (row_num, row)
        end_age = start_age + (row_num - 1) + (length(row) - 1)
        OffsetArray([row ; ultimate[end_age + 1:end]], (start_age - 1) + (row_num - 1))
    end

    return OffsetArray(vs, start_age - 1)
end

"""
    MortalityTable(ultimate)
    MortalityTable(select, ultimate)
    MortalityTable(select, ultimate; metadata::MetaData)

Constructs a container object which can hold either:
    - ultimate-only rates (an `UltimateTable`)
    - select and ultimate rates (a `SelectUltimateTable`)

Also pass a keyword argument `metadata=MetaData(...)` to store relevant information (source, notes, etc) about the table itself.

# Examples
```julia
# first construct the underlying data
ult = UltimateMortality([x / 100 for x in 0:100]); # first ma
matrix = rand(10,50); # represents random mortality rates with a select period of 10 years
sel = SelectMortality(matrix,ult,start_age=0);

table = MortalityTable(sel,ult)

# can now get rates, indexed by attained age:

table.select[10] # the vector of rates for a risk issued select at age 10

table.ultimate[99] # 0.99

```
"""
abstract type MortalityTable end

struct SelectUltimateTable{S,U} <: MortalityTable
    select::S
    ultimate::U
    metadata::TableMetaData
end

struct UltimateTable{U} <: MortalityTable
    ultimate::U
    metadata::TableMetaData
end

Base.getindex(u::UltimateTable,x) = u.ultimate[x]
Base.lastindex(u::UltimateTable) = lastindex(u.ultimate)

function MortalityTable(select, ultimate; metadata=TableMetaData())
    return SelectUltimateTable(select, ultimate, metadata)
end

function MortalityTable(ultimate; metadata=TableMetaData())
    return UltimateTable(ultimate, metadata)
end

Base.show(io::IO, ::MIME"text/plain", mt::MortalityTable) = print(
    io,
    """
    MortalityTable ($(mt.metadata.content_type)):
       Name:
           $(mt.metadata.name)
       Fields:
           $(fieldnames(typeof(mt)))
       Provider:
           $(mt.metadata.provider)
       mort.SOA.org ID:
           $(mt.metadata.id)
       mort.SOA.org link:
           https://mort.soa.org/ViewTable.aspx?&TableIdentity=$(mt.metadata.id)
       Description:
           $(mt.metadata.description)
    """,
)

"""
    survival(mortality_vector,to_age)
    survival(mortality_vector,from_age,to_age)

Returns the survival through attained age `to_age`. The start of the calculation is either the start of the vector, or attained_age `from_age`. `from_age` and `to_age` need to be Integers. Add a DeathDistribution as the last argument to handle floating point and non-whole ages:

    survival(mortality_vector,to_age,::DeathDistribution)
    survival(mortality_vector,from_age,to_age,::DeathDistribution)

If given a negative `to_age`, it will return `1.0`. Aside from simplifying the code, this makes sense as for something to exist in order to decrement in the first place, it must have existed and surived to the point of  being able to be decremented.

# Examples
```julia-repl
julia> qs = UltimateMortality([0.1,0.3,0.6,1]);

julia> survival(qs,0)
1.0
julia> survival(qs,1)
0.9

julia> survival(qs,1,1)
1.0
julia> survival(qs,1,2)
0.7

julia> survival(qs,0.5,Uniform())
0.95
```
"""
function survival(v, to_age)
    return survival(v, firstindex(v), to_age)
end
function survival(v, to_age, dd::DeathDistribution)
    return survival(v, firstindex(v), to_age, dd)
end

function survival(v::T, from_age::Int, to_age::Int) where {T <: AbstractArray}
    if from_age == to_age
        return 1.0
    else
        return @views reduce(*,
            1 .- v[from_age:(to_age - 1)],
            init=1.0
            )
    end
end

function survival(v::T, from_age, to_age, dd::DeathDistribution) where {T <: AbstractArray}
    # calculate the survival for the rounded ages, and then the high and low high_residual
    age_low = ceil(Int, from_age)
    age_high = floor(Int, to_age)

    # if from_age and to_age are fractional parts of the same attained age, then age_high will round down to
    # be below the rounded-up age_low. This line will short circuit the rest and just return the fractional year survival
    age_high < age_low && return 1 - decrement_partial_year(v, from_age, to_age, dd)

    if age_low == from_age
        low_residual = 1.0
    else
        low_residual = 1 - decrement_partial_year(v, from_age, age_low, dd)
    end

    if age_high == to_age
        high_residual = 1.0
    else
        high_residual = 1 - decrement_partial_year(v, age_high, to_age, dd)
    end

    if from_age == to_age
        return 1.0
    else

        whole = @views reduce(*,
            1 .- v[age_low:(age_high - 1)],
            init=1.0
            )

        return whole * low_residual * high_residual
    end
end

# Reference: Experience Study Calculations, 2016, Society of Actuaries
# https://www.soa.org/globalassets/assets/Files/Research/2016-10-experience-study-calculations.pdf
function decrement_partial_year(v, from_age, to_age, dd::Uniform)
    return v[floor(Int, from_age)] * (to_age - from_age)
end

function decrement_partial_year(v, from_age, to_age, dd::Constant)
    return 1 - (1 - v[floor(Int, from_age)])^(to_age - from_age)
end

function decrement_partial_year(v, from_age, to_age, dd::Balducci)
    q′ = v[floor(Int, from_age)]
    frac = (to_age - from_age)
    return 1 - (1 - q′) / (1 - (1 - frac) * q′)
end

"""
    decrement(mortality_vector,to_age)
    decrement(mortality_vector,from_age,to_age)

Returns the cumulative decrement through attained age `to_age`. The start of the calculation is either the start of the vector, or attained_age `from_age`. `from_age` and `to_age` need to be Integers. Add a DeathDistribution as the last argument to handle floating point and non-whole ages:

    decrement(mortality_vector,to_age,::DeathDistribution)
    decrement(mortality_vector,from_age,to_age,::DeathDistribution)

# Examples
```julia-repl
julia> qs = UltimateMortality([0.1,0.3,0.6,1]);

julia> decrement(qs,0)
0.0
julia> decrement(qs,1)
0.1

julia> decrement(qs,1,1)
0.0
julia> decrement(qs,1,2)
0.3

julia> decrement(qs,0.5,Uniform())
0.05
```
"""
decrement(v,to_age) = 1 .- survival(v, to_age)
decrement(v,to_age,dd::DeathDistribution) = 1 .- survival(v, to_age, dd)
decrement(v,from_age,to_age) = 1 .- survival(v, from_age, to_age)
decrement(v,from_age,to_age,dd::DeathDistribution) = 1 .- survival(v, from_age, to_age, dd)

"""
    omega(x)
    ω(x)

Returns the last index of the given vector. For mortality vectors this means the last attained age for which a rate is defined.

Note that `omega` can vary depending on the issue age for a select table, and that a select `omega` may differ from the table's ultimate `omega`.

ω is aliased to omega, but unexported. To use, do `using MortalityTables: ω` when importing or call `MortalityTables.ω()`

# Examples

```julia-repl
julia> qs = UltimateMortality([0.1,0.3,0.6,1]);
julia> omega(qs)
3

julia> qs = UltimateMortality([0.1,0.3,0.6,1],start_age=10);
julia> omega(qs)
13

```
"""
function omega(x)
    return lastindex(x)
end

ω = omega

"""
    mortality_vector(vec; start_age=0)

A convienience constructor to create an OffsetArray which has the array indexed by attained age rather than always starting from `1`. The package and JuliaActuary ecosystem assume that the rates are indexed by attained age and this allows transformation of tables without a direct dependency on **OffsetArrays.jl**.

Equivalent to doing:
```
using OffsetArrays
OffsetArray(vec,start_age-1)
```
"""
mortality_vector(vec; start_age=0) = return OffsetArray(vec, start_age - 1)


text:
include("./layer.jl")
include("./spiking.jl")
using LoopVectorization, Polyester

mutable struct Network
    layers::Array

    runtime::Int64
    spiking_trains::Array{Float32}

    default_input_size::Int64

    add_layer
    activate
    train

    function Network(;runtime::Int64, input_size::Int64)
        spiking_trains = zeros(Float32, runtime, input_size)
        new(Any[], runtime, spiking_trains, input_size, add_layer, activate, train)
    end

    function add_layer(self::Network; args...)
        push!(self.layers, Layer(;input_size=self.default_input_size, args...))
        self.default_input_size = self.layers[end].layer_size
    end

    function activate(self::Network, input_data::Array{Float32}, output_data::Array{Float32}, table::Array{Int64}; η, A₊, A₋, γᵥ, γₚ, γₘ, γₜ, w₊, w₋, r, scales, guide, record)
        for layer in self.layers
            layer.reset(layer, r)
        end
        @avxt self.spiking_trains .= 0.0f0

        θ = Spiking!(input_data, self.spiking_trains, scales, r, self.runtime)*γₜ
        γᵥ = 1-γᵥ/θ

        self.layers[guide].guide = true
        self.layers[guide].sample_output = output_data
        self.layers[record].record = true

        first_spike = false
        for t in 1:self.runtime
            first_spike = self.layers[1].activate(self.layers[1], self.spiking_trains[t,:], η, A₊, A₋, γᵥ, γₚ, γₘ, w₊, w₋, θ, t, r, first_spike, table)
            for l in 2:length(self.layers)
                first_spike = self.layers[l].activate(self.layers[l], self.layers[l-1].output, η, A₊, A₋, γᵥ, γₚ, γₘ, w₊, w₋, θ, t, r, first_spike, table)
            end
        end
    end

    function train(self::Network; input_data::Array{Float32}, output_data::Array{Float32}, epochs::Int64=20, batch_size::Int64=500, parameters...)
        current_input_data = zeros(Float32, size(input_data,1))
        current_output_data = zeros(Float32, size(output_data,1))
        table = zeros(Int64, self.layers[end].layer_size,self.layers[end].layer_size)
        accuracy = 0.0

        for e in 1:epochs
            print("Epoch ", e, "\n[")

            @time begin
                for t in 1:batch_size
                    index = rand(1:size(input_data)[end])
                    @avxt current_input_data .= input_data[:,index]
                    @avxt current_output_data .= output_data[:,index]

                    if t%ceil(batch_size/50)==0
                        print("=")
                    end
                    self.activate(self, current_input_data, current_output_data, table; parameters...)
                end
                @avxt for i in 1:10
                    accuracy += table[i,i]
                end
                accuracy = accuracy/batch_size*100
                print("] with accuracy $accuracy%\nTime usage: ")
                accuracy = 0
            end
            display(table)
            println("\n")
            @avxt table .= 0
        end
    end
end


text:
export reconstruction

################################################################################
include("DirectReconstruction.jl")
include("IterativeReconstruction.jl")
include("RecoParameters.jl")

"""
    reconstruction(acqData::AcquisitionData, recoParams::Dict)

takes an AcquisitionData object and a parameter
dictionary and calculates an image from the given raw data.

Reconstruction types are specified by the symbol `:reco`.
Valid reconstruction names are:
* :direct - direct Fourier reconstruction
* :standard           - iterative reconstruction for all contrasts, coils & slices independently
* :multiEcho          - iterative joint reconstruction of all echo images
* :multiCoil          - SENSE-type iterative reconstruction
* :multiCoilMultiEcho - SENSE-type iterative reconstruction of all echo images
"""
function reconstruction(acqData::AcquisitionData, recoParams::Dict)
  encodingDims = dims(trajectory(acqData))
  if encodingDims == 2 #"2D"
    return reconstruction_2d(acqData,recoParams)
  elseif encodingDims == 3 #"3D"
    return reconstruction_3d(acqData,recoParams)
  else
    error("$(encodingDims)d encoding not yet supported")
  end
  return reconstruction_2d(acqData,recoParams)
end

"""
    reconstruction_2d(acqData::AcquisitionData, recoParams::Dict)

Performs image reconstruction of a 2d encoded AcquisitionData object.
Parameters are specified in a dictionary.

Reconstruction types are specified by the symbol `:reco`.
Valid reconstruction names are:
* :direct - direct Fourier reconstruction
* :standard           - iterative reconstruction for all contrasts, coils & slices independently
* :multiEcho          - iterative joint reconstruction of all echo images
* :multiCoil          - SENSE-type iterative reconstruction
* :multiCoilMultiEcho - SENSE-type iterative reconstruction of all echo images
"""
function reconstruction_2d(acqData::AcquisitionData, recoParams::Dict)
  recoParams = merge(defaultRecoParams(), recoParams)

  # direct reco
  if recoParams[:reco] == "direct"
    reconSize, weights, cmap = setupDirectReco(acqData, recoParams)
    return reconstruction_direct_2d(acqData, reconSize, weights, cmap)
  end

  # iterative reco
  par = setupIterativeReco(acqData, recoParams)
  if recoParams[:reco] == "standard"
    return reconstruction_simple(acqData, par.reconSize, par.reg, par.sparseTrafo, par.weights, par.solvername, par.correctionMap, par.method, par.normalize, recoParams)
  elseif recoParams[:reco] == "multiEcho"
    return reconstruction_multiEcho(acqData, par.reconSize, par.reg, par.sparseTrafo, par.weights, par.solvername, par.correctionMap, par.method, par.normalize, recoParams)
  elseif recoParams[:reco] == "multiCoil"
    return reconstruction_multiCoil(acqData, par.reconSize, par.reg, par.sparseTrafo, par.weights, par.solvername, par.senseMaps, par.correctionMap, par.method, par.normalize, recoParams)
  elseif recoParams[:reco] == "multiCoilMultiEcho"
    return reconstruction_multiCoilMultiEcho(acqData, par.reconSize, par.reg, par.sparseTrafo, par.weights, par.solvername, par.senseMaps, par.correctionMap, par.method, par.normalize, recoParams)
  else
    @error "RecoModel $(recoParams[:reco]) not found."
  end

  return reconstruction_direct_2d(acqData, recoParams)
end

"""
    reconstruction_3d(acqData::AcquisitionData, recoParams::Dict)

Performs image reconstruction of a 3d encoded AcquisitionData object.
Parameters are specified in a dictionary.

Reconstruction types are specified by the symbol `:reco`.
Valid reconstruction names are:
* :direct - direct Fourier reconstruction
* :standard           - iterative reconstruction for all contrasts, coils & slices independently
* :multiEcho          - iterative joint reconstruction of all echo images
* :multiCoil          - SENSE-type iterative reconstruction
* :multiCoilMultiEcho - SENSE-type iterative reconstruction of all echo images
"""
function reconstruction_3d(acqData::AcquisitionData, recoParams::Dict)
  recoParams = merge(defaultRecoParams(), recoParams)
  if recoParams[:reco] == "direct"
    reconSize, weights, cmap = setupDirectReco(acqData, recoParams)
    return reconstruction_direct_3d(acqData, reconSize, weights, cmap)
  end

  acqData2d = convert3dTo2d(acqData)
  par = setupIterativeReco(acqData2d, recoParams)
  if recoParams[:reco] == "standard"
    Ireco = reconstruction_simple(acqData2d, par.reconSize, par.reg, par.sparseTrafo, par.weights, par.solvername, par.correctionMap, par.method, par.normalize, recoParams)
  elseif recoParams[:reco] == "multiEcho"
    Ireco = reconstruction_multiEcho(acqData2d, par.reconSize, par.reg, par.sparseTrafo, par.weights, par.solvername, par.correctionMap, par.method, par.normalize, recoParams)
  elseif recoParams[:reco] == "multiCoil"
    Ireco = reconstruction_multiCoil(acqData2d, par.reconSize, par.reg, par.sparseTrafo, par.weights, par.solvername, par.senseMaps, par.correctionMap, par.method, par.normalize, recoParams)
  elseif recoParams[:reco] == "multiCoilMultiEcho"
    Ireco = reconstruction_multiCoilMultiEcho(acqData2d, par.reconSize, par.reg, par.sparseTrafo, par.weights, par.solvername, par.senseMaps, par.correctionMap, par.method, par.normalize, recoParams)
  else
    @error "RecoModel $(recoParams[:reco]) not found."
  end
  Ireco = permutedims(Ireco,[3,1,2,4,5])

  return Ireco
end

"""
    reconstruction(acqData::AcquisitionData, recoParams::Dict,filename::String; force=false)

performs the same image reconstrucion as `reconstruction(acqData::AcquisitionData, recoParams::Dict)`
and saves the image in a file with name `filename`.
If `force=false`, the reconstructed image is loaded from the the file `filename` if the latter is
present.
"""
function reconstruction(acqData::AcquisitionData, recoParams::Dict, filename::String;
                        force=false)
  if !force && isfile(filename)
    return recoImage( RecoFileIBI(filename) )
  else
    I = reconstruction(acqData, recoParams)
    saveasRecoFile(filename, I, recoParams)
    return I
  end
end


text:
function mform(n1::Int64,l1::Int64,m1::Int64,n2::Int64,l2::Int64,ρ::Float64,q::Float64)

    # The function wigner3j_f produces a vector of wigner symbols where each index corresponds to a value of l between
    # abs(l1-l2) and l1+l2 (and assumes m = -m1-m2) meaning that:
    #    collect(wigner3j_f(l1, l2, m1, m2)) ≈ [wigner3j(l1, l2, l3, m1, m2, -m1-m2) for l ∈ abs(l1-l2):l1+l2]

    wigner3j_l1_l2__0_0 = wigner3j_f(l1, l2, 0, 0)
    wigner3j_l1_l2__m1_negm2 = wigner3j_f(l1, l2, m1, -m1)

    δ = (ρ*n2+ n1)/(n1*ρ*n2)
    β = 2/n1
    α = 2*l1+1
    dd = -q^2/δ^2
    N1 = (4π*4/(ρ^(3/2) * (n1 * n2)^2)
          * √(gamma(n1-l1) * gamma(n2-l2)
              / (gamma(n1+l1+1) * gamma(n2+l2+1)))
          * (2/(n1))^l1 * (2/(ρ * n2))^l2 * √(π/(2q)))

    if n1 > l1 >= m1 && n2 > l2 >= m1
        # @floop from the package Floops.jl is able to parallelize a loop for us. If the loop contains a 'reduction'
        # (in our case, the sum), we just tell it which variable is being reduced with the @reduce macro and it takes
        # care of initializing the reduction variables and will make sure the multi-threading does not introduce a 'race
        # condition' which could cause the sum to be incorrect.
        @floop for l in abs(l1-l2):2:l1+l2
            for k in 0:n2-l2-1
                for k1 in 0:n1-l1-1
                    if (l1+l2) % 2 == 0
                        @reduce Sum += mfun(n1,l1,m1,n2,l2,ρ,q,k,k1,l, δ, β, α, dd, N1,
                                            wigner3j_l1_l2__0_0, wigner3j_l1_l2__m1_negm2)
                    else
                        @reduce Sum += mfun(n1,l1,m1,n2,l2,ρ,q,k,k1,l,  δ, β, α, dd, N1,
                                                 wigner3j_l1_l2__0_0, wigner3j_l1_l2__m1_negm2) #without the i factor
                    end
                end
            end
        end
        return N1 * Sum
    else
        "Incorrect Combination"
    end
end

function mfun(n1::Int64,l1::Int64,m1::Int64,n2::Int64,l2::Int64,ρ::Float64,q::Float64,k::Int64,k1::Int64,l::Int64,
              δ, β, α, dd, N1,
              wigner3j_l1_l2__0_0, wigner3j_l1_l2__m1_negm2)
    m=0
    n = n1-l1-1
    γ = l1+l2+k+3/2
    ν = l+1/2
    aa = (ν+γ+k1+1)/2
    bb = (ν+γ+k1+2)/2
    cc = 1+ν

    I1 = gamma(n2+l2+1) * gamma(n+α+1) * ( (-2)^k * (-β)^k1 * q^ν *gamma(ν+γ+k1+1)
          / (gamma(k+1) * gamma(n2-l2-k) * gamma(2*l2+k+2) * (ρ*n2)^k * gamma(k1+1)
             * gamma(n-k1+1) * gamma(α+k1+1) * 2^ν * gamma(ν+1) * δ^(ν+γ+k1+1)) * _₂F₁(aa,bb,cc,dd))

    if (l1+l2) % 2 == 0
        A1 = (-1)^(l/2+m1) * √((2*l+1) * (2*l1+1) * (2*l2+1) * (2*l+1) / (4π*4π))
    else
        A1 = (-1)^((l-1)/2+m1) * √((2*l+1) * (2*l1+1) * (2*l2+1) * (2*l+1) / (4π*4π))
    end
    return I1 * A1 * wigner3j_l1_l2__0_0[l] * wigner3j_l1_l2__m1_negm2[l]
end


text:
let
    n = 50
    p = 3
    t = (1:n)/n
    Ut,Vt = egrss.generators(t,p);
    Wt,c = egrss.potrf(Ut,Vt,1e-2*ones(n))
    Yt,Zt = egrss.trtri(Ut,Wt,c)

    Lref = cholesky(egrss.full(Ut,Vt) + Diagonal(1e-2*ones(n))).L
    Lref_inv = LowerTriangular(tril(inv(Lref)));

    L_inv = egrss.full_tril(Yt,Zt,1.0./c)
    @test ( isapprox(L_inv,Lref_inv) )

    L_inv = egrss.trtri2(Ut,Wt,c)
    @test ( isapprox(L_inv,Lref_inv) )
end

max_stars_repo_name: martinandersen/EGRSS

max_issues_repo_name: martinandersen/EGRSS

max_forks_repo_name: martinandersen/EGRSS


text:
 # This code tests the function kepler_driftij2

@testset "kepler_driftij" begin
for drift_first in [true,false]
# Next, try computing two-body Keplerian Jacobian:

n = 3
H = [3,1,1]
t0 = 7257.93115525
#t0 = -300.0
h  = 0.0000005
hbig  = big(h)
tmax = 600.0
#dlnq = 1e-8
dlnq = big(1e-20)

elements = readdlm("elements.txt",',')
#elements[2,1] = 0.75
elements[2,1] = 1.0
elements[3,1] = 1.0

m =zeros(n)
x0=zeros(NDIM,n)
v0=zeros(NDIM,n)

for k=1:n
  m[k] = elements[k,1]
end
for iter = 1:2

init = ElementsIC(elements,H,t0)
x0,v0,_ = init_nbody(init)
 if iter == 2
   # Reduce masses to trigger hyperbolic routine:
    m[1:n] *= 1e-1
    h = 0.0000005
    hbig = big(h)
 end
# Tilt the orbits a bit:
x0[2,1] = 5e-1*sqrt(x0[1,1]^2+x0[3,1]^2)
x0[2,2] = -5e-1*sqrt(x0[1,2]^2+x0[3,2]^2)
v0[2,1] = 5e-1*sqrt(v0[1,1]^2+v0[3,1]^2)
v0[2,2] = -5e-1*sqrt(v0[1,2]^2+v0[3,2]^2)

jac_ij = zeros(Float64,14,14)
dqdt_ij = zeros(Float64,14)
dqdt_num = zeros(BigFloat,14)

println("Initial values: ",x0,v0)
println("masses: ",m)
i=1 ; j=2
x = copy(x0) ; v=copy(v0)
xerror = zeros(NDIM,n); verror = zeros(NDIM,n)
d = Derivatives(Float64)
# Predict values of s:
kepler_driftij!(m,x,v,xerror,verror,i,j,h,jac_ij,dqdt_ij,drift_first,d)
x0 = copy(x) ; v0 = copy(v)
xerror = zeros(NDIM,n); verror = zeros(NDIM,n)
xbig = big.(x) ; vbig=big.(v); mbig = big.(m)
xerr_big = zeros(BigFloat,NDIM,n); verr_big = zeros(BigFloat,NDIM,n)
kepler_driftij!(m,x,v,xerror,verror,i,j,h,jac_ij,dqdt_ij,drift_first,d)
# Now compute Jacobian with BigFloat precision:
jac_ij_big = zeros(BigFloat,14,14)
dqdt_ij_big = zeros(BigFloat,14)
KEPLER_TOL = sqrt(eps(big(1.0)))
dBig = Derivatives(BigFloat)
kepler_driftij!(mbig,xbig,vbig,xerr_big,verr_big,i,j,hbig,jac_ij_big,dqdt_ij_big,drift_first,dBig)
#println("jac_ij: ",convert(Array{Float64,2},jac_ij_big))
#println("jac_ij - jac_ij_big: ",convert(Array{Float64,2},jac_ij_big)-jac_ij)
println("max(jac_ij - jac_ij_big): ",maxabs(convert(Array{Float64,2},jac_ij_big)-jac_ij))

# Now, compute the derivatives numerically:
jac_ij_num = zeros(BigFloat,14,14)
xsave = big.(x)
vsave = big.(v)
msave = big.(m)

# Compute the time derivatives:
# Initial positions, velocities & masses:
xm = big.(x0)
vm = big.(v0)
mm = big.(msave)
dq = dlnq * hbig
hbig -= dq
kepler_driftij!(mm,xm,vm,i,j,hbig,drift_first)
xp = big.(x0)
vp = big.(v0)
hbig += 2dq
kepler_driftij!(mm,xp,vp,i,j,hbig,drift_first)
# Now x & v are final positions & velocities after time step
for k=1:3
  dqdt_num[   k] = .5*(xp[k,i]-xm[k,i])/dq
  dqdt_num[ 3+k] = .5*(vp[k,i]-vm[k,i])/dq
  dqdt_num[ 7+k] = .5*(xp[k,j]-xm[k,j])/dq
  dqdt_num[10+k] = .5*(vp[k,j]-vm[k,j])/dq
end
hbig = big(h)

# Compute position, velocity & mass derivatives:
for jj=1:3
  # Initial positions, velocities & masses:
  xm = big.(x0)
  vm = big.(v0)
  mm = big.(msave)
  dq = dlnq * xm[jj,i]
  if xm[jj,i] != 0.0
    xm[jj,i] -=  dq
  else
    dq = dlnq
    xm[jj,i] = -dq
  end
  kepler_driftij!(mm,xm,vm,i,j,hbig,drift_first)
  xp = big.(x0)
  vp = big.(v0)
  if xm[jj,i] != 0.0
    xp[jj,i] +=  dq
  else
    dq = dlnq
    xp[jj,i] = dq
  end
  kepler_driftij!(mm,xp,vp,i,j,hbig,drift_first)
  # Now x & v are final positions & velocities after time step
  for k=1:3
    jac_ij_num[   k,  jj] = .5*(xp[k,i]-xm[k,i])/dq
    jac_ij_num[ 3+k,  jj] = .5*(vp[k,i]-vm[k,i])/dq
    jac_ij_num[ 7+k,  jj] = .5*(xp[k,j]-xm[k,j])/dq
    jac_ij_num[10+k,  jj] = .5*(vp[k,j]-vm[k,j])/dq
  end
  xm = big.(x0)
  vm = big.(v0)
  mm  = big.(msave)
  dq = dlnq * vm[jj,i]
  if vm[jj,i] != 0.0
    vm[jj,i] -=  dq
  else
    dq = dlnq
    vm[jj,i] = -dq
  end
  kepler_driftij!(mm,xm,vm,i,j,hbig,drift_first)
  xp = big.(x0)
  vp = big.(v0)
  mm  = big.(msave)
  if vp[jj,i] != 0.0
    vp[jj,i] +=  dq
  else
    dq = dlnq
    vp[jj,i] = dq
  end
  kepler_driftij!(mm,xp,vp,i,j,hbig,drift_first)
  for k=1:3
    jac_ij_num[   k,3+jj] = .5*(xp[k,i]-xm[k,i])/dq
    jac_ij_num[ 3+k,3+jj] = .5*(vp[k,i]-vm[k,i])/dq
    jac_ij_num[ 7+k,3+jj] = .5*(xp[k,j]-xm[k,j])/dq
    jac_ij_num[10+k,3+jj] = .5*(vp[k,j]-vm[k,j])/dq
  end
end

# Now vary mass of inner planet:
xm= big.(x0)
vm= big.(v0)
mm= big.(msave)
dq = mm[i]*dlnq
mm[i] -= dq
kepler_driftij!(mm,xm,vm,i,j,hbig,drift_first)
xp= big.(x0)
vp= big.(v0)
mp= big.(msave)
dq = mp[i]*dlnq
mp[i] += dq
kepler_driftij!(mp,xp,vp,i,j,hbig,drift_first)
for k=1:3
  jac_ij_num[   k,7] = .5*(xp[k,i]-xm[k,i])/dq
  jac_ij_num[ 3+k,7] = .5*(vp[k,i]-vm[k,i])/dq
  jac_ij_num[ 7+k,7] = .5*(xp[k,j]-xm[k,j])/dq
  jac_ij_num[10+k,7] = .5*(vp[k,j]-vm[k,j])/dq
end
# The mass doesn't change:
jac_ij_num[7,7] =  1.0

for jj=1:3
  # Now vary parameters of outer planet:
  xm = big.(x0)
  vm = big.(v0)
  mm = big.(msave)
  dq = dlnq * xm[jj,j]
  if xm[jj,j] != 0.0
    xm[jj,j] -=  dq
  else
    dq = dlnq
    xm[jj,j] = -dq
  end
  kepler_driftij!(mm,xm,vm,i,j,hbig,drift_first)
  xp = big.(x0)
  vp = big.(v0)
  if xp[jj,j] != 0.0
    xp[jj,j] +=  dq
  else
    dq = dlnq
    xp[jj,j] = dq
  end
  kepler_driftij!(mm,xp,vp,i,j,hbig,drift_first)
  for k=1:3
    jac_ij_num[   k,7+jj] = .5*(xp[k,i]-xm[k,i])/dq
    jac_ij_num[ 3+k,7+jj] = .5*(vp[k,i]-vm[k,i])/dq
    jac_ij_num[ 7+k,7+jj] = .5*(xp[k,j]-xm[k,j])/dq
    jac_ij_num[10+k,7+jj] = .5*(vp[k,j]-vm[k,j])/dq
  end
  xm= big.(x0)
  vm= big.(v0)
  mm = big.(msave)
  dq = dlnq * vm[jj,j]
  if vm[jj,j] != 0.0
    vm[jj,j] -=  dq
  else
    dq = dlnq
    vm[jj,j] = -dq
  end
  kepler_driftij!(mm,xm,vm,i,j,hbig,drift_first)
  xp= big.(x0)
  vp= big.(v0)
  if vp[jj,j] != 0.0
    vp[jj,j] +=  dq
  else
    dq = dlnq
    vp[jj,j] = dq
  end
  kepler_driftij!(mm,xp,vp,i,j,hbig,drift_first)
  for k=1:3
    jac_ij_num[   k,10+jj] = .5*(xp[k,i]-xm[k,i])/dq
    jac_ij_num[ 3+k,10+jj] = .5*(vp[k,i]-vm[k,i])/dq
    jac_ij_num[ 7+k,10+jj] = .5*(xp[k,j]-xm[k,j])/dq
    jac_ij_num[10+k,10+jj] = .5*(vp[k,j]-vm[k,j])/dq
  end
end

# Now vary mass of outer planet:
xm = big.(x0)
vm = big.(v0)
mm = big.(msave)
dq = mm[j]*dlnq
mm[j] -= dq
kepler_driftij!(mm,xm,vm,i,j,hbig,drift_first)
xp = big.(x0)
vp = big.(v0)
mp = big.(msave)
dq = mp[j]*dlnq
mp[j] += dq
kepler_driftij!(mp,xp,vp,i,j,hbig,drift_first)
for k=1:3
  jac_ij_num[   k,14] = .5*(xp[k,i]-xm[k,i])/dq
  jac_ij_num[ 3+k,14] = .5*(vp[k,i]-vm[k,i])/dq
  jac_ij_num[ 7+k,14] = .5*(xp[k,j]-xm[k,j])/dq
  jac_ij_num[10+k,14] = .5*(vp[k,j]-vm[k,j])/dq
end
# The mass doesn't change:
jac_ij_num[14,14] =  1.0

#println(jac_ij)
#println(jac_ij_num)
#println(jac_ij./jac_ij_num)
emax = 0.0; imax = 0; jmax = 0
emax_big = big(0.0); imax_big = 0; jmax_big = 0
for i=1:14
  jac_ij[i,i] += 1
  jac_ij_big[i,i] += 1
end
for i=1:14, j=1:14
  if jac_ij[i,j] != 0.0
    diff = abs(convert(Float64,jac_ij_num[i,j])/jac_ij[i,j]-1.0)
    if  diff > emax
      emax = diff; imax = i; jmax = j
    end
    diff_big = abs(convert(Float64,jac_ij_num[i,j])/jac_ij_big[i,j]-1.0)
    if  diff_big > emax_big
      emax_big = diff; imax_big = i; jmax_big = j
    end
  end
end
println("Maximum fractional error: ",emax," ",imax," ",jmax)
println("Maximum fractional error big: ",emax_big," ",imax_big," ",jmax_big)
#println(jac_ij)
#println(convert(Array{Float64,2},jac_ij_num))
println("Maximum jac_ij error:   ",maxabs(convert(Array{Float64,2},asinh.(jac_ij_num))-asinh.(jac_ij)))
println("Maximum jac_ij_big-jac_ij_num:   ",maxabs(convert(Array{Float64,2},asinh.(jac_ij_num)-asinh.(jac_ij_big))))
println("Max dqdt error: ",maxabs(dqdt_ij-convert(Array{Float64,1},dqdt_num)))

@test isapprox(jac_ij_num,jac_ij;norm=maxabs)
#@test isapprox(dqdt_ij,convert(Array{Float64,1},dqdt_num);norm=maxabs)

end
end
end


text:
# Completed: 8/13

import Statistics: mean, var, std, varm, stdm

# TODO:
# cor
# cov
# mean done
# mean! overwrites.
# median
# median! overwrites.
# middle
# quantile
# quantile! overwrites.
# std done.
# stdm done.
# var done.
# varm done.

@primitive mean(x;d...),dy  (dy.*one.(x)./(length(x)÷length(dy)))
@primitive mean(f::typeof(abs),x;d...),dy   nothing  (dy.*sign.(x)./(length(x)÷length(dy)))
@primitive mean(f::typeof(abs2),x;d...),dy  nothing  (dy.*(2x)./(length(x)÷length(dy)))

std(x::Value, args...; kws...) = sqrt.(var(x, args...; kws...))
stdm(x::Value, args...; kws...) = sqrt.(varm(x, args...; kws...))

var(x::Value; corrected::Bool=true, mean=nothing, dims=:)=_varm(x, mean; corrected=corrected, dims=dims)
varm(x::Value, m; corrected::Bool=true, dims=:)=_varm(x, m; corrected=corrected, dims=dims)

function _varm(x, m; corrected::Bool=true, dims=:)
    if m === nothing; m = mean(x,dims=dims); end
    s = sum(abs2, x .- m; dims=dims)
    r = length(x) ÷ length(s) - Int(corrected)
    s ./ r
end

# I don't think making _varm a primitive improves performance, replaces a single sum op with a single _varm op.
# @primitive _varm(x, m; corrected::Bool=true, dims=:),dy  (dy.*(x.-m).*(2/(length(x)÷length(dy)-Int(corrected))))  nothing
# @primitive sum(x;d...),dy  (dy.*one.(x))
# @primitive sum(f::typeof(abs),x;d...),dy   nothing  (dy.*sign.(x))
# @primitive sum(f::typeof(abs2),x;d...),dy  nothing  (dy.*(2x))


text:
# Account

import HTTP: request
import JSON: parse

const valid_types =
    ["messages", "rules", "rule-records", "alerts", "alert-records"]

struct UsageTypeException <: Exception
    msg::AbstractString

    function UsageTypeException()
        msg = "quota_type must be one of " * join(valid_types, ", ")
        new(msg)
    end
end

"""
    usage(get_iex; quota_type="messages")

Get usage for account associated with IEX connection conn
"""
function usage(conn::Connection; quota_type::String = "messages")
    if !(quota_type in valid_types)
        throw(UsageTypeException())
    end

    url = conn.url * conn.version * "/account/usage/" * quota_type
    r = request("GET", url; query = Dict("token" => conn.token))
    return parse(String(r.body))
end

"""
    metadata(conn::Connection)

Get metadata for account associated wth IEX connection conn
"""
function metadata(conn::Connection)
    url = conn.url * conn.version * "/account/metadata"
    r = request("GET", url; query = Dict("token" => conn.token))
    return parse(String(r.body))
end

max_stars_repo_name: GodolphinCapitalManagement/Iex

max_issues_repo_name: GodolphinCapitalManagement/Iex

max_forks_repo_name: GodolphinCapitalManagement/Iex


text:
using Oceananigans.Architectures
using Oceananigans.Architectures: device
import Oceananigans.Architectures: architecture
using CUDA, CUDA.CUSPARSE
using KernelAbstractions: @kernel, @index

using LinearAlgebra, SparseArrays, IncompleteLU
using SparseArrays: nnz

import LinearAlgebra.ldiv!

"""
`ILUFactorization` (`preconditioner_method = :ILUFactorization`)
    stores two sparse lower and upper trianguilar matrices `L` and `U` such that `LU ≈ A`
    is applied to `r` with `forward_substitution!(L, r)` followed by `backward_substitution!(U, r)`
    constructed with `ilu(A, τ = drop_tolerance)`

`SparseInversePreconditioner` (`preconditioner_method = :SparseInverse` and `:AsymptoticInverse`)
    stores a sparse matrix `M` such that `M ≈ A⁻¹`
    is applied to `r` with a matrix multiplication `M * r`
    constructed with
    `asymptotic_diagonal_inverse_preconditioner(A)`
        -> is an asymptotic expansion of the inverse of A assuming that A is diagonally dominant
        -> it is possible to choose order 0 (Jacobi), 1 or 2
    `sparse_approximate_preconditioner(A, ε = tolerance, nzrel = relative_maximum_number_of_elements)`
        -> same formulation as Grote M. J. & Huckle T, "Parallel Preconditioning with sparse approximate inverses"
        -> starts constructing the sparse inverse of A from identity matrix until, either a tolerance (ε) is met or nnz(M) = nzrel * nnz(A)

The suggested preconditioners are

on the `CPU`
`ilu()` (superior to everything always and in every situation!)

on the `GPU`
`aymptotic_diagonal_inverse_preconditioner()` (if `Δt` is variable or large problem_sizes)
`sparse_inverse_preconditioner()` (if `Δt` is constant and problem_size is not too large)

as a rule of thumb, for poisson solvers:
`sparse_inverse_preconditioner` is better performing than `asymptotic_diagonal_inverse_preconditioner` only if `nzrel >= 2.0`
As such, we urge to use `sparse_inverse_preconditioner` only when
- Δt is constant (we don't have to recalculate the preconditioner during the simulation)
- it is feasible to choose `nzrel = 2.0` (for not too large problem sizes)

Note that `asymptotic_diagonal_inverse_preconditioner` assumes the matrix to be diagonally dominant, for this reason it could
be detrimental when used on non-diagonally dominant system (cases where Δt is very large). In this case it is better
to use `sparse_inverse_preconditioner`

`ilu()` cannot be used on the GPU because preconditioning the solver with a direct LU (or Choleski) type
of preconditioner would require too much computation for the `ldiv!(P, r)` step completely hindering the performances
"""

validate_settings(T, arch, settings)                                  = settings
validate_settings(::Val{:Default}, arch, settings)                    = arch isa CPU ? (τ = 0.001, ) : (order = 1, )
validate_settings(::Val{:SparseInverse}, arch, settings::Nothing)     = (ε = 0.1, nzrel = 2.0)
validate_settings(::Val{:ILUFactorization}, arch, settings::Nothing)  = (τ = 0.001, )
validate_settings(::Val{:AsymptoticInverse}, arch, settings::Nothing) = (order = 1, )

validate_settings(::Val{:ILUFactorization}, arch, settings)  = haskey(settings, :τ) ?
                                                                      settings :
                                                                      throw(ArgumentError("τ has to be specified for ILUFactorization"))
validate_settings(::Val{:SparseInverse}, arch, settings)     = haskey(settings, :ε) && haskey(settings, :nzrel) ?
                                                                      settings :
                                                                      throw(ArgumentError("both ε and nzrel have to be specified for SparseInverse"))
validate_settings(::Val{:AsymptoticInverse}, arch, settings) = haskey(settings, :order) ?
                                                                      settings :
                                                                      throw(ArgumentError("and order ∈ [0, 1, 2] has to be specified for AsymptoticInverse"))

function build_preconditioner(::Val{:Default}, matrix, settings)
    default_method = architecture(matrix) isa CPU ? :ILUFactorization : :AsymptoticInverse
    return build_preconditioner(Val(default_method), matrix, settings)
end

build_preconditioner(::Val{nothing},            A, settings)  = Identity()
build_preconditioner(::Val{:SparseInverse},     A, settings)  = sparse_inverse_preconditioner(A, ε = settings.ε, nzrel = settings.nzrel)
build_preconditioner(::Val{:AsymptoticInverse}, A, settings)  = asymptotic_diagonal_inverse_preconditioner(A, asymptotic_order = settings.order)

function build_preconditioner(::Val{:ILUFactorization},  A, settings)
    if architecture(A) isa GPU
        throw(ArgumentError("the ILU factorization is not available on the GPU! choose another method"))
    else
        return ilu(A, τ = settings.τ)
    end
end

@inline architecture(::CuSparseMatrixCSC) = CUDAGPU()
@inline architecture(::SparseMatrixCSC)   = CPU()

abstract type AbstractInversePreconditioner{M} end

function  LinearAlgebra.ldiv!(u, precon::AbstractInversePreconditioner, v)
    mul!(u, matrix(precon), v)
end

function  LinearAlgebra.ldiv!(precon::AbstractInversePreconditioner, v)
    mul!(v, matrix(precon), v)
end

struct SparseInversePreconditioner{M} <: AbstractInversePreconditioner{M}
    Minv :: M
end

@inline matrix(p::SparseInversePreconditioner)  = p.Minv

"""
The diagonally dominant inverse preconditioner is constructed with an asymptotic expansion of `A⁻¹` up to the second order
If `I` is the Identity matrix and `D` is the matrix containing the diagonal of `A`, then

the 0th order expansion is the Jacobi preconditioner `M = D⁻¹ ≈ A⁻¹`

the 1st order expansion corresponds to `M = D⁻¹(I - (A - D)D⁻¹) ≈ A⁻¹`

the 2nd order expansion corresponds to `M = D⁻¹(I - (A - D)D⁻¹ + (A - D)D⁻¹(A - D)D⁻¹) ≈ A⁻¹`

all preconditioners are calculated on CPU and then moved to the GPU.
Additionally the first order expansion has a method to calculate the preconditioner directly on the GPU
`asymptotic_diagonal_inverse_preconditioner_first_order(A)` in case of variable time step where the preconditioner
has to be recalculated often
"""

function asymptotic_diagonal_inverse_preconditioner(A::AbstractMatrix; asymptotic_order)

    arch                  = architecture(A)
    constr                = deepcopy(constructors(arch, A))
    colptr, rowval, nzval = copy_unpack_constructors(arch, constr)
    dev                   = device(arch)

    M = size(A, 1)

    invdiag = arch_array(arch, zeros(eltype(nzval), M))

    loop! = _get_inv_diag!(dev, 256, M)
    event = loop!(invdiag, colptr, rowval, nzval; dependencies=Event(dev))
    wait(dev, event)

    if asymptotic_order == 0
        Minv_cpu = spdiagm(0=>arch_array(CPU(), invdiag))
        Minv     = arch_sparse_matrix(arch, Minv_cpu)
    elseif asymptotic_order == 1
        loop! = _initialize_asymptotic_diagonal_inverse_preconditioner_first_order!(dev, 256, M)
        event = loop!(nzval, colptr, rowval, invdiag; dependencies=Event(dev))
        wait(dev, event)

        constr_new = (colptr, rowval, nzval)
        Minv = arch_sparse_matrix(arch, constructors(arch, M, constr_new))
    else
        D   = spdiagm(0=>diag(arch_sparse_matrix(CPU(), A)))
        D⁻¹ = spdiagm(0=>arch_array(CPU(), invdiag))
        Minv_cpu = D⁻¹ * (I - (A - D) * D⁻¹ + (A - D) * D⁻¹ * (A - D) * D⁻¹)
        Minv = arch_sparse_matrix(arch, Minv_cpu)
    end

    return SparseInversePreconditioner(Minv)
end

@kernel function _initialize_asymptotic_diagonal_inverse_preconditioner_first_order!(nzval, colptr, rowval, invdiag)
    col = @index(Global, Linear)

    for idx = colptr[col] : colptr[col+1] - 1
        if rowval[idx] == col
            nzval[idx] = invdiag[col]
        else
            nzval[idx] = - nzval[idx] * invdiag[rowval[idx]] * invdiag[col]
        end
    end
end

function sparse_inverse_preconditioner(A::AbstractMatrix; ε, nzrel)

   # let's choose an initial sparsity => diagonal
   A_cpu    = arch_sparse_matrix(CPU(), A)
   Minv_cpu = sparse_approximate_inverse(A_cpu, ε = ε, nzrel = nzrel)

   Minv = arch_sparse_matrix(architecture(A), Minv_cpu)
   return SparseInversePreconditioner(Minv)
end


text:
function toCartesian(lon, lat; r = 1, cxyz = (0,0,0) )
    lat, lon = lat*π/180, lon*π/180
    cxyz[1] + r * cos(lat) * cos(lon), cxyz[2] + r * cos(lat) * sin(lon), cxyz[3] + r *sin(lat)
end
function lonlat3D(lon, lat, data; cxyz = (0,0,0))
    xyzw = zeros(size(data)..., 3)
    for (i,lon) in enumerate(lon), (j,lat) in enumerate(lat)
        x, y, z = toCartesian(lon, lat; cxyz = cxyz)
        xyzw[i,j,1] = x
        xyzw[i,j,2] = y
        xyzw[i,j,3] = z
    end
    xyzw[:,:,1], xyzw[:,:,2], xyzw[:,:,3]
end

function SphereBall(; r = 0.995, cxyz = (0,0,0))
    Θ = LinRange(0, 2π, 500) # 50
    Φ = LinRange(0, π, 500)
    x0 = [cxyz[1] + r * cos(θ) * sin(ϕ)      for θ in Θ, ϕ in Φ]
    y0 = [cxyz[2] + r * sin(θ) * sin(ϕ)      for θ in Θ, ϕ in Φ]
    z0 = [cxyz[3] + r * cos(ϕ) for θ in Θ, ϕ in Φ]
    x0, y0, z0
end


text:
include("typographyaux.jl")

export text, textFont, textSize
export textWidth, textHeight, textAscent, textDescent

## Loading & Displaying

#loadFont

function text(str::AbstractString, x, y)
    xo = x
    yo = y
	x = ((x+1)/2)*state.width
	y = ((y+1)/2)*state.height

	shader("fontDrawing")
	glActiveTexture(GL_TEXTURE1)
    glEnable(GL_CULL_FACE)

    posData = zeros(GLfloat, 2*6*length(str))
    texData = zeros(GLfloat, 2*6*length(str))

    n = 0
    totyadv = 0
    nforline = 0
	@inbounds for c in str
        # if we hit a newline character or we reach the max line
        # length, then move to next line, where the next line is
        # determined by the average y advance from the characters
        # of the previous line.
        if c == '\n' || nforline == fontState.maxLineLength
            cht = fontState.characters[str[1]]
            x = xo
            y -= totyadv/nforline
            totyadv = 0
            nforline = 0
            continue
        end

		ch = fontState.characters[c]

		xpos = x + ch.bearing[1] * state.textSize
        ypos = y - (ch.size[2] - ch.bearing[2]) * state.textSize

        w = ch.size[1] * state.textSize
        h = ch.size[2] * state.textSize

        totyadv += ch.advance[2] * state.textSize

        if w == 0 || h == 0
            continue
        end

        posData[n+1] = xpos
        posData[n+2] = ypos + h

        posData[n+3] = xpos
        posData[n+4] = ypos

        posData[n+5] = xpos + w
        posData[n+6] = ypos

        posData[n+7] = xpos
        posData[n+8] = ypos + h

        posData[n+9] = xpos + w
        posData[n+10] = ypos

        posData[n+11] = xpos + w
        posData[n+12] = ypos + h

        texData[n+1] = ch.atlasOffset
        texData[n+2] = 0.0

        texData[n+3] = ch.atlasOffset
        texData[n+4] = ch.size[2] / fontState.atlasHeight

        texData[n+5] = ch.atlasOffset + ch.size[1] / fontState.atlasWidth
        texData[n+6] = ch.size[2] / fontState.atlasHeight

        texData[n+7] = ch.atlasOffset
        texData[n+8] = 0.0

        texData[n+9] = ch.atlasOffset + ch.size[1] / fontState.atlasWidth
        texData[n+10] = ch.size[2] / fontState.atlasHeight

        texData[n+11] = ch.atlasOffset + ch.size[1] / fontState.atlasWidth
        texData[n+12] = 0.0

		x += ch.advance[1] * state.textSize
        n += 2*6
        nforline += 1
	end
    glBindBuffer(GL_ARRAY_BUFFER, globjs.posvbos[globjs.posind])
    glBufferData(GL_ARRAY_BUFFER, sizeof(posData), posData, GL_STATIC_DRAW)
    glBindBuffer(GL_ARRAY_BUFFER, globjs.texvbos[globjs.texind])
    glBufferData(GL_ARRAY_BUFFER, sizeof(texData), texData, GL_STATIC_DRAW)
    glDrawArrays(GL_TRIANGLES, 0, 6*length(str))

    glDisable(GL_CULL_FACE)
	shader("basicShapes")

    sfx = state.fbSize[1]/state.width; sfy = state.fbSize[2]/state.height
    x = (x./state.width - 0.5)*sfx
    y = ((y - totyadv/nforline)./state.height - 0.5)*sfy
    return [xo, yo, x, y]
end

function textFont(fontname::AbstractString)
    # for the time being, we only allow system fonts or those that have been
    # installed directly into the main system storage
    @windows_only state.fontFace = "C:/Windows/Fonts/"*fontname
    @linux_only state.fontFace = "/usr/share/fonts/"*fontname
    @osx_only state.fontFace = "/System/Library/Fonts/"*fontname
    fontState = fontStruct(newface(state.fontFace), Dict(' ' => blankChar), 0, 12)
	setpixelsize(fontState.face, fontState.fontWidth, fontState.fontHeight)
	setupFontCharacters()
end

## Attributes

#function textAlign()
#
#end

#textLeading
#textMode

function textSize(size)
	state.textSize = size
end

function textWidth(str::AbstractString)
    # extents = Cairo.text_extents(cr, str)
    # return extents[1]
end

function textHeight(str::AbstractString)
    # extents = Cairo.text_extents(cr, str)
    # return extents[2]
end

## Metrics

function textAscent(str::AbstractString)
   # extents = Cairo.scaled_font_extents(cr, str)
   # return extents[1]
end

function textDescent(str::AbstractString)
   # extents = Cairo.scaled_font_extents(cr, str)
   # return extents[2]
end

function textLineLength(l)
    fontState.maxLineLength = l
end