text:
module LBV

using StaticArrays

export AbstractLBV, LBVLeaf, LBVNode
export @define_node

abstract type AbstractLBV{T, D<:AbstractVector{T}} <: AbstractVector{T} end

#note: it is still unclear whether there is any benefit in allowing zero-length
#blocks. it may allow treating systems more homogeneously, regardless of whether
#they have inputs/states or not. when this becomes clear, remove support for
#zero length blocks if it does not afford anything

########################### LBVLeaf ############################

struct LBVLeaf{L, T, D <: AbstractVector{T}} <: AbstractLBV{T, D}
    data::D
    function LBVLeaf{L,T,D}(data) where {L,T,D} #for some reason, you can't restrict type parameters here!
        if !(L > 0)
            throw(ArgumentError("LBVLeaf length must be positive"))
        end
        if length(data) != length(LBVLeaf{L,T,D})
            throw(ArgumentError("Got input length $(length(data)), expected $(length(LBVLeaf{L,T,D}))"))
        end
        new{L,T,D}(data)
    end
end
LBVLeaf{L,T}(data::D) where {L,T,D} = LBVLeaf{L,T,D}(data) #T/D inconsistencies caught by the inner const
LBVLeaf{L}(data::D) where {L,D} = LBVLeaf{L,eltype(D),D}(data)
LBVLeaf(data::D) where {D} = LBVLeaf{length(data),eltype(data),typeof(data)}(data) #avoid for efficiency

LBVLeaf{L,T}(::UndefInitializer) where {L,T} = LBVLeaf{L,T}(MVector{L,T}(undef))
LBVLeaf{L}(::UndefInitializer) where {L} = LBVLeaf{L,Float64}(undef)

LBVLeaf{L,T}() where {L,T} = LBVLeaf{L,T}(undef)
LBVLeaf{L}() where {L} = LBVLeaf{L,Float64}()

Base.similar(::Type{<:LBVLeaf{L,T}}) where {L,T} = LBVLeaf{L,T}(MVector{L,T}(undef))
Base.similar(::LBVLeaf{L,T}) where {L,T} = similar(LBVLeaf{L,T})
Base.copy(x::LBVLeaf{L}) where {L} = LBVLeaf{L}(copy(getfield(x, :data)))
# Base.convert(::Type{<:LBVLeaf{L}}, v::AbstractVector) where {L} = LBVLeaf{L}(v)

###### Abstract Array #######

Base.length(::Type{<:LBVLeaf{L}}) where {L} = L
Base.size(::LBVLeaf{L}) where {L} = (L,)
Base.getindex(x::LBVLeaf, i) = getindex(getfield(x,:data), i)
Base.setindex!(x::LBVLeaf, v, i) = setindex!(getfield(x,:data), v, i)

####### Custom Broadcasting #######

struct LBVLeafStyle{L,T} <: Broadcast.AbstractArrayStyle{1} end

LBVLeafStyle{L,T}(::Val{1}) where {L,T} = LBVLeafStyle{L,T}()
Base.BroadcastStyle(::Type{LBVLeaf{L,T,D}}) where {L,T,D} = LBVLeafStyle{L,T}()

function Base.similar(::Broadcast.Broadcasted{LBVLeafStyle{L,T}}, ::Type{ElType}) where {L,T,ElType}
    similar(LBVLeaf{L,T})
end

function Base.BroadcastStyle(::LBVLeafStyle{L,T1}, ::LBVLeafStyle{L,T2}) where {L,T1,T2}
    LBVLeafStyle{L,promote_type(T1, T2)}()
end

########################### LBVNode ############################

struct LBVNode{S, T, D <: AbstractVector{T}} <: AbstractLBV{T, D} #S: identifier
    data::D
    function LBVNode{S,T,D}(data::D) where {S,T,D} #for some reason, you can't restrict type parameters here
        assert_symbol(S)
        if length(data) != length(LBVNode{S,T,D})
            throw(ArgumentError("Got input length $(length(data)), expected $(length(LBVNode{S,T,D}))"))
        end
        new{S,T,D}(data)
    end
end
@generated assert_symbol(x) = x<:Symbol ? nothing : :(throw(TypeError(:LBVNode, "inner constructor", Symbol, $x)))

LBVNode{S,T}(data::D) where {S,T,D} = LBVNode{S,T,D}(data) #T/D inconsistencies caught by the inner const
LBVNode{S}(data::D) where {S,D} = LBVNode{S,eltype(D),D}(data)

LBVNode{S,T}(::UndefInitializer) where {S,T} = LBVNode{S,T}(MVector{length(LBVNode{S,T}) ,T}(undef))
LBVNode{S}(::UndefInitializer) where {S} = LBVNode{S,Float64}(undef)

LBVNode{S,T}() where {S,T} = LBVNode{S,T}(undef)
LBVNode{S}() where {S} = LBVNode{S,Float64}()

####### Abstract Array #############

Base.size(x::LBVNode) = size(getfield(x,:data))
Base.getindex(x::LBVNode, i) = getindex(getfield(x,:data), i)
Base.setindex!(x::LBVNode, v, i) = setindex!(getfield(x,:data), v, i)

Base.similar(::Type{<:LBVNode{S,T}}) where {S,T} = LBVNode{S,T}(MVector{length(LBVNode{S,T}),T}(undef))
Base.similar(::LBVNode{S,T}) where {S,T} = similar(LBVNode{S,T})
Base.copy(x::LBVNode{S}) where {S} = LBVNode{S}(copy(getfield(x, :data)))
# Base.convert(::Type{<:LBVNode{S}}, v::AbstractVector) where {S} = LBVNode{S}(v)

####### Custom Broadcasting #######

struct LBVNodeStyle{S,T} <: Broadcast.AbstractArrayStyle{1} end

LBVNodeStyle{S,T}(::Val{1}) where {S,T} = LBVNodeStyle{S,T}()
Base.BroadcastStyle(::Type{LBVNode{S,T,D}}) where {S,T,D} = LBVNodeStyle{S,T}()

function Base.similar(::Broadcast.Broadcasted{LBVNodeStyle{S,T}}, ::Type{ElType}) where {S,T,ElType}
    similar(LBVNode{S,T})
end

function Base.BroadcastStyle(::LBVNodeStyle{S,T1}, ::LBVNodeStyle{S,T2}) where {S,T1,T2}
    LBVNodeStyle{S,promote_type(T1, T2)}()
end

######### Code Generation #########

Base.getproperty(x::LBVNode, s::Symbol) = getproperty(x, Val(s))
Base.setproperty!(x::LBVNode, s::Symbol, v) = setproperty!(x, Val(s), v)

#define and register a LBVNode
macro define_node(name, children)

    ex = quote

        println("Generating code for $($(QuoteNode(name))) = ",
                "LBVNode{:$($(QuoteNode(name)))}...")

        if isdefined($__module__, $(QuoteNode(name)))
            println("WARNING: Type $($(QuoteNode(name))) already defined ",
            "in module $($(QuoteNode(__module__)))")
        end

        #use a let block to create a local scope so that:

        #1) variables captured in methods do not inhabit the global module
        #   scope. this avoids the performance hit from capturing global
        #   variables in a closure. the global keyword is required to make
        #   methods and variables defined within the let block visible in the
        #   global module scope.

        #2) macro hygiene is respected (we have escaped the whole quote)

        let children = $(children)

            global const $(name) = LBVNode{$(QuoteNode(name))}

            println(children)
            child_labels = keys(children)
            child_types = values(children)
            @assert all([t <: Union{LBVLeaf, LBVNode} for t in child_types]) "Invalid child type"
            node_length = sum(length.(child_types))

            global Base.propertynames(::$(name), private::Bool = false) = child_label
            global Base.length(::Type{<:$(name)}) = node_length
            global get_children(::Type{<:$(name)}) = children
            global get_children(::$(name)) = children

            offset = 0
            for (child_label, child_type) in zip(child_labels, child_types)

                child_length = length(child_type)
                child_range = (1 + offset):(child_length + offset)
                offset += child_length
                println("Generating helper methods for child $child_label, range $child_range")

                global function Base.getproperty(x::$(name), ::Val{child_label})
                    return (child_type)(view(getfield(x,:data), child_range))
                end
                global function Base.setproperty!(x::$(name), ::Val{child_label}, v)
                    return setindex!(getfield(x, :data), v, child_range)
                end

            end

            global function Base.show(io::IO, ::MIME"text/plain", x::$(name))

                println(io, "$node_length-element ", typeof(x), " with blocks:")
                for (child_label, child_type) in zip(child_labels, child_types)
                    print(io, "\t", child_label, ": ", getproperty(x, child_label)[:])
                    child_label == child_labels[end] ? break : println()
                end
            end

        end

    end
    return Base.remove_linenums!(esc(ex)) #esc: escapes the whole expression
    # return esc(ex)
end

end #module


text:
const DATE_FORMAT = dateformat"yyyy-mm-dd H\hM\mS\s ZZZZ"
const DATE_FORMAT2 = dateformat"yyyy-mm-dd HH:MM:SS ZZZZ"

DIR = @__DIR__

PREAMBLE = """<meta charset="utf-8">

"""

PART1 = mt"""# Analysis Results for Sample *{{{:samplename}}}*

Sample *{{:samplename}}* was found to contain the following.

| **Compound** | **Amount (w%)** |
|:------------ | -----------:|
{{#:matches}}| {{name}}          | {{percent}}      |
{{/:matches}}

## Sample Information
{{#:info}}
- **{{{key}}}:** {{{value}}}
{{/:info}}
"""

PLOT_HEADER = mt"""### {{{:component}}} component
- Library source: {{{:source}}}
- Fit Score: {{{:overall_fit}}}
"""

PART2 = mt"""## Signal Reconstruction and Residue

Portion of signal recovered by reconstruction: {{:recovery}}.

Solvent, internal standard, and chemical shift reference omitted for clarity.
"""

EMPTY = mt"""# Analysis Results for Sample *{{{:samplename}}}*

Sample *{{:samplename}}* did not match any of the spectra in the reference library.

## Sample Information
{{#:info}}
- **{{{key}}}:** {{{value}}}
{{/:info}}
"""

"""
    report(s::NMR.Spectrum, lib::Array{NMR.Spectrum,1}, d::NMR.DecompositionResult, quants, darks, is, path, short = true)

Generate an HTML report from decomposition results `d` and quantitation results `quants`, obtained from the analysis of NMR spectrum `s` analyzed against library group `lib`.

The report will be written to the directory `path`. `darks`, an array of NMR spectra containing 'dark' resonances, and `is`, the internal standard spectrum, are passed for reference.
If `short` is **true** a short report is generated.
"""
function report(s::NMR.Spectrum, lib::Array{NMR.Spectrum,1}, d::NMR.DecompositionResult, quants, darks, is,
                path, short = true)
    sampinf = sampleinfo(s)
    if haskey(sampinf, "Long report")
        short = !(lowercase(sampinf["Long report"]) ∈ ["true", "yes"])
    end
    mp = match_plots(s, lib, d, quants, short)
    plots = Dict( refno => IOBuffer() for refno=keys(mp) ) # Buffers to hold plot output
    for (k,v) in mp
        show(plots[k], MIME("text/html"), v)
    end
    matches = [ Dict( "name" => name(lib[k]),"percent" => statistics(v))
                for (k,v) in quants ]

    #augment sample info
    if haskey(s.acqu, "AUTOPOS")
        sampinf["Autosampler position"] = string(s["AUTOPOS"])
    end
    sampinf["NMRquant algorithm version"] = findversion("NMRquant")
    sampinf["Library version"] = findlibversion(lib)
    sampinf["Experiment number"] = string(s.expno)
    sampinf["Date of data collection"] = string(Dates.format(astimezone(ZonedDateTime(unix2datetime(s.acqu["DATE"]), tz"UTC"), localzone()), DATE_FORMAT2))
    sampinf["Analysis method"] = s.acqu["EXP"]

    info = [ Dict( "key" => k, "value" => v) for (k,v) in sorted_kv(sampinf) ]
    overall_fit = Dict{Int64, Float64}()
    for refno in keys(mp)
        refindex = findfirst(x->x==refno,d.refnums)
        overall_fit[refno] = *((fit_score for fit_score in d.fit_scores[refindex])...)^(1/length(d.fit_scores[refindex]))
    end
    _, recon, residue = NMR.decompose(d)
    recon_plt = reconplot(s, recon, residue, darks)
    recovery = @sprintf "%.1f" 100.0 * norm(recon,1) / norm(recon .+ residue,1)
    md1 = render(PART1, samplename=s.name, matches=matches, info=info)
    md2 = render(PART2, recovery=recovery)
    #testing
    #v3.3.2:
    #outfile = joinpath(path, "report " * zonedtime(now(), DATE_FORMAT) * ".html")
    #v3.4.0: Sample name in report file name
    outfile = joinpath(path, "$(get(sampinf, "Sample", "report")) " * zonedtime(now(), DATE_FORMAT) * ".html")
    open(outfile, "w") do f
        write(f, PREAMBLE)
        if !isempty(quants)
            write(f, Markdown.html(Markdown.parse(md1)))
        else
            write(f, Markdown.html(Markdown.parse(render(EMPTY, samplename=s.name, info=info))))
        end
        for (refno,p) in plots
            write(f, Markdown.html(Markdown.parse(render(PLOT_HEADER, component=name(lib[refno]), source=get(sampleinfo(lib[refno]),"Source",""), overall_fit=score_format(overall_fit[refno])))))
            write(f, p.data)
        end
        if !short
            write(f, Markdown.html(Markdown.parse(md2)))
            show(f, MIME("text/html"), recon_plt)
        end
    end
    outfile
end

LIBRARY_PREAMBLE = """<meta charset="utf-8">

"""

LIBRARY_TITLE = mt"""# DAS NMRquant {{:name}} Library v{{:version}} Integrals
"""

COLUMN_WIDTHS = """
    <col width='16%'>
    <col width='16%'>
    <col width='12%'>
    <col width='56%'>
"""

TABLE_TITLE = """

        Compound
        Number of Integral(s)
        Number 1H(s) represented in each Integral Window
        Integral Windows

"""

LIBRARY_TABLE="""

        {{:name}}
        {{:numint}}
        {{:intprotons}}
"""

function library_report(lib::Array{NMR.Spectrum,1}, path)
    lp = library_plots(lib)
    plots = Dict( refno => IOBuffer() for refno=keys(lp) ) # Buffers to hold plot output
    for (k,v) in lp
        show(plots[k], MIME("text/html"), v)
    end

    outfile = joinpath(path, "DAS NMRquant $(lib[1].name) Library v$(findlibversion(lib)) " * zonedtime(now(), DATE_FORMAT) * ".html")
    title = render(LIBRARY_TITLE, name=lib[1].name, version=findlibversion(lib))
    open(outfile, "w") do f
        write(f, LIBRARY_PREAMBLE)
        write(f, Markdown.html(Markdown.parse(title)))
        write(f, "
", COLUMN_WIDTHS, TABLE_TITLE)
        write(f, "
", COLUMN_WIDTHS)

        libnum = Dict(name(lib[x]) => x for x = keys(plots))
        for n in sort(collect(keys(libnum)))
            ref = lib[libnum[n]]

            intprotons = ""
            intprotons = "[" * string(round(Int, protons(ref)[1]))
            if length(protons(ref)) > 1
                for i = 2:length(protons(ref))
                    intprotons *= ","*string(round(Int, protons(ref)[i]))
                end
            end
            intprotons *= "]"

            write(f, render(LIBRARY_TABLE, name=n, numint=length(protons(ref)), intprotons=intprotons), "", plots[libnum[n]].data, "
")
        end
        write(f, "
")
    end
    outfile
end

VALIDATION_PREAMBLE = """<meta charset="utf-8">

"""

VRPART1 = """# Validation Report

- Date: {{{:reportdate}}}
- Program version: {{{:version}}}
- {{{:libname}}} Library version: {{{:libversion}}}

## Summary:

| **Test** | **Result** |
|:------------ | :-----------:|
{{#:r}}| {{test}}          | {{{result}}}      |
{{/:r}}
"""

VRPART2 = """### Quantitation Results:

| **Sample** | **Expected** | **Result** |
| :--------- | :---------- | :-------- |
{{#:r}}|  {{{sample}}} | {{#expected}} {{.}}
 {{/expected}} |  {{#result}} {{.}}
 {{/result}}  |
{{/:r}}

## Test log:
"""

function validation_report(lib::Array{NMR.Spectrum,1}, results::Dict{String, Bool}, quantlist::Dict, quantresults, testlogs::Dict{String,IOBuffer}, outpath)
    version = findversion("NMRquant")
    libname = lib[1].name
    libversion = findlibversion(lib)
    reportdate = now()
    outfile = joinpath(outpath, "Validation Report NMRquant v$version $libname v$libversion $(zonedtime(reportdate, DATE_FORMAT)).html")
    open(outfile, "w") do f
        write(f, VALIDATION_PREAMBLE)

        r = [ Dict( "test" => k,"result" => v ? "Passed" : "Failed")
            for (k,v) in sorted_kv(results) ]
        vrp1 = Markdown.html(Markdown.parse(render(VRPART1, reportdate=zonedtime(reportdate, DATE_FORMAT2),version=version, libname=libname, libversion=libversion,r=r)))
        write(f, restorehtml(vrp1))

        r = [ Dict( "sample" => k,
                    "expected" => haskey(quantlist[k],0) ? ["None"] : ["$(ev)% $(sampleinfo(lib[en])["Name"])" for (en,ev) in quantlist[k]],
                    "result" => haskey(quantresults[k][1],0) ? ["None"] : ["$(statistics(qv))% $(sampleinfo(lib[qn])["Name"])" for (qn, qv) in quantresults[k][1]],
                    "passed" => quantresults[k][2] ? "true" : "")
            for (k,_) in sorted_kv(quantlist) ]
        temp = Markdown.html(Markdown.parse(render(VRPART2, r=r)))
        write(f, restorehtml(temp))

        write(f, "
")
        for k in sort(collect(keys(testlogs)))
            write(f, Markdown.html(Markdown.parse("### $(k)\n")))

            #Writes with ANSI codes converted to HTML
            write(f, ansitohtml(String(take!(copy(testlogs[k])))))
        end
        write(f, "
")
        close(f)
    end
    outfile
end

ext: jl


text:
module UnivariateWrapper

using JavaCall

const AbstractConvergenceChecker = @jimport org.hipparchus.optim.AbstractConvergenceChecker
const BaseOptimizer = @jimport org.hipparchus.optim.BaseOptimizer
const BracketFinder = @jimport org.hipparchus.optim.univariate.BracketFinder
const BrentOptimizer = @jimport org.hipparchus.optim.univariate.BrentOptimizer
const Class = @jimport java.lang.Class
const ConvergenceChecker = @jimport org.hipparchus.optim.ConvergenceChecker
const GoalType = @jimport org.hipparchus.optim.nonlinear.scalar.GoalType
const JString = @jimport java.lang.String
const MultiStartUnivariateOptimizer = @jimport org.hipparchus.optim.univariate.MultiStartUnivariateOptimizer
const Object = @jimport java.lang.Object
const OptimizationData = @jimport org.hipparchus.optim.OptimizationData
const RandomGenerator = @jimport org.hipparchus.random.RandomGenerator
const SearchInterval = @jimport org.hipparchus.optim.univariate.SearchInterval
const SimpleUnivariateValueChecker = @jimport org.hipparchus.optim.univariate.SimpleUnivariateValueChecker
const UnivariateFunction = @jimport org.hipparchus.analysis.UnivariateFunction
const UnivariateObjectiveFunction = @jimport org.hipparchus.optim.univariate.UnivariateObjectiveFunction
const UnivariateOptimizer = @jimport org.hipparchus.optim.univariate.UnivariateOptimizer
const UnivariatePointValuePair = @jimport org.hipparchus.optim.univariate.UnivariatePointValuePair

export BracketFinder
export BrentOptimizer
export MultiStartUnivariateOptimizer
export SearchInterval
export SimpleUnivariateValueChecker
export UnivariateObjectiveFunction
export UnivariateOptimizer
export UnivariatePointValuePair
export converged
export equals
export get_absolute_threshold
export get_class
export get_convergence_checker
export get_evaluations
export get_f_hi
export get_f_lo
export get_f_mid
export get_goal_type
export get_hi
export get_iterations
export get_lo
export get_max
export get_max_evaluations
export get_max_iterations
export get_mid
export get_min
export get_objective_function
export get_optima
export get_point
export get_relative_threshold
export get_start_value
export get_value
export hash_code
export notify
export notify_all
export optimize
export search
export to_string
export wait

include("bracket_finder.jl")
include("brent_optimizer.jl")
include("multi_start_univariate_optimizer.jl")
include("search_interval.jl")
include("simple_univariate_value_checker.jl")
include("univariate_objective_function.jl")
include("univariate_optimizer.jl")
include("univariate_point_value_pair.jl")

end


text:
module PseudospectraView

using Plots
using Pseudospectra
using QML
if VERSION >= v"1.6"
    using Pkg: dependencies
    function _getver(pname)
        pd = dependencies()
        for p in values(pd)
            if p.name == pname
                return p.version
            end
        end
        return nothing
    end
    if _getver("QML") >= v"0.7"
        using Qt5QuickControls_jll
    end
end

export psagui

# CHECKME: why do QML examples use dirname(@__FILE__) instead?
const pvsrcdir = realpath(@__DIR__)

# this is for tracking GUI-related activity, distinct from PSA option.
# 0: quiet, 1: reassure us that something is going on, 2: help the developer
const _verbosity = Ref(0)

# backend for Plots
const _pbackend = Ref(:unset)
const _supported_backends = [:pyplot, :gr]

include("psadata.jl")

function psagui(; backend::Symbol=:pyplot)
    (backend in _supported_backends) || throw(ArgumentError("unsupported Plots backend"))
    _pbackend[] = backend
    include(joinpath(pvsrcdir,"psagui.jl"))
    results = Base.invokelatest(PSApp.runme)
    return (length(keys(results)) == 0) ? nothing : results
end

"""
    psagui()
    psagui(A[,opts])

Start the QML GUI for Pseudospectra, optionally loading matrix `A` first.
`opts` is a `Dict{Symbol,Any}` of Pseudospectra options (many of which may
be reset by GUI objects).
Returns a `Dict` with variables saved by the GUI, or `nothing`.
"""
function psagui(A::AbstractMatrix,
                opts::Dict{Symbol,Any}=Dict{Symbol,Any}();
                backend::Symbol=:pyplot)
    (backend in _supported_backends) || throw(ArgumentError("unsupported Plots backend"))
    _pbackend[] = backend
    Base.invokelatest(PSAData.setdefaultmatrix,A)
    isempty(opts) || Base.invokelatest(PSAData.addopts,opts)
    include(joinpath(pvsrcdir,"psagui.jl"))
    results = Base.invokelatest(PSApp.runme)
    return (length(keys(results)) == 0) ? nothing : results
end

# this is for automated testing
const _doomed = Ref(false)

end # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "GSL"
version = v"1.16"

# Collection of sources required to build GSL
sources = [
    "http://ftp.gnu.org/gnu/gsl/gsl-$(version.major).$(version.minor).tar.gz" =>
    "73bc2f51b90d2a780e6d266d43e487b3dbd78945dd0b04b14ca5980fe28d2f53",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gsl-*/
update_configure_scripts

# We need to massage configure script to convince it to build the shared library
# for PowerPC.
if [[ "${target}" == powerpc64le-* ]]; then
    autoreconf -vi
fi

./configure --prefix=$prefix --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built.
# Note that the products we are declaring here should be kept in-sync with those of `GSL@2`,
# so that users that don't care about versions can simply use (e.g. `libgsl`) without having
# to worry about whether `libgsl` is called something else in `GSL@2` versus `GSL@1`.
products = [
    LibraryProduct("libgsl", :libgsl)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
module TwoDTurb

export
  Problem,
  set_zeta!,
  updatevars!,

  energy,
  enstrophy,
  dissipation,
  work,
  drag

using
  FFTW,
  Reexport

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!
using FourierFlows: getfieldspecs, structvarsexpr, parsevalsum, parsevalsum2

abstract type TwoDTurbVars <: AbstractVars end

nothingfunction(args...) = nothing

"""
    Problem(; parameters...)

Construct a 2D turbulence problem.
"""
function Problem(;
    # Numerical parameters
          nx = 256,
          Lx = 2π,
          ny = nx,
          Ly = Lx,
          dt = 0.01,
    # Drag and/or hyper-/hypo-viscosity
           ν = 0,
          nν = 1,
           μ = 0,
          nμ = 0,
    # Timestepper and eqn options
     stepper = "RK4",
       calcF = nothingfunction,
  stochastic = false,
           T = Float64,
         dev = CPU()
)

  gr = TwoDGrid(dev, nx, Lx, ny, Ly; T=T)
  pr = Params{T}(ν, nν, μ, nμ, calcF)
  vs = calcF == nothingfunction ? Vars(dev, gr) : (stochastic ? StochasticForcedVars(dev, gr) : ForcedVars(dev, gr))
  eq = Equation(pr, gr)
  FourierFlows.Problem(eq, stepper, dt, gr, vs, pr, dev)
end

# ----------
# Parameters
# ----------

"""
    Params(ν, nν, μ, nμ, calcF!)

Returns the params for two-dimensional turbulence.
"""
struct Params{T} <: AbstractParams
       ν :: T         # Vorticity viscosity
      nν :: Int       # Vorticity hyperviscous order
       μ :: T         # Bottom drag or hypoviscosity
      nμ :: Int       # Order of hypodrag
  calcF! :: Function  # Function that calculates the forcing F
end
Params(ν, nν) = Params(ν, nν, typeof(ν)(0), 0, nothingfunction)

# ---------
# Equations
# ---------

"""
    Equation(p, g)

Returns the equation for two-dimensional turbulence with params p and grid g.
"""
function Equation(p::Params, g::AbstractGrid{T}) where T
  L = @. - p.ν*g.Krsq^p.nν - p.μ*g.Krsq^p.nμ
  L[1, 1] = 0
  FourierFlows.Equation(L, calcN!, g)
end

# ----
# Vars
# ----

struct Vars{Aphys, Atrans, F, P} <: TwoDTurbVars
     zeta :: Aphys
        u :: Aphys
        v :: Aphys
    zetah :: Atrans
       uh :: Atrans
       vh :: Atrans
       Fh :: F
  prevsol :: P
end

const ForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, Nothing}
const StochasticForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, <:AbstractArray}

"""
    Vars(dev, g)

Returns the vars for unforced two-dimensional turbulence on device dev and with
grid g.
"""
function Vars(::Dev, g::AbstractGrid{T}) where {Dev, T}
  @devzeros Dev T (g.nx, g.ny) zeta u v
  @devzeros Dev Complex{T} (g.nkr, g.nl) zetah uh vh
  Vars(zeta, u, v, zetah, uh, vh, nothing, nothing)
end

"""
    ForcedVars(dev, g)

Returns the vars for forced two-dimensional turbulence on device dev and with
grid g.
"""
function ForcedVars(dev::Dev, g::AbstractGrid{T}) where {Dev, T}
  @devzeros Dev T (g.nx, g.ny) zeta u v
  @devzeros Dev Complex{T} (g.nkr, g.nl) zetah uh vh Fh
  Vars(zeta, u, v, zetah, uh, vh, Fh, nothing)
end

"""
    StochasticForcedVars(dev, g)

Returns the vars for stochastically forced two-dimensional turbulence on device
dev and with grid g.
"""
function StochasticForcedVars(dev::Dev, g::AbstractGrid{T}) where {Dev, T}
  @devzeros Dev T (g.nx, g.ny) zeta u v
  @devzeros Dev Complex{T} (g.nkr, g.nl) zetah uh vh Fh prevsol
  Vars(zeta, u, v, zetah, uh, vh, Fh, prevsol)
end

# -------
# Solvers
# -------

"""
    calcN_advection(N, sol, t, cl, v, p, g)

Calculates the advection term.
"""
function calcN_advection!(N, sol, t, cl, v, p, g)
  @. v.uh =   im * g.l  * g.invKrsq * sol
  @. v.vh = - im * g.kr * g.invKrsq * sol
  @. v.zetah = sol

  ldiv!(v.u, g.rfftplan, v.uh)
  ldiv!(v.v, g.rfftplan, v.vh)
  ldiv!(v.zeta, g.rfftplan, v.zetah)

  @. v.u *= v.zeta # u*zeta
  @. v.v *= v.zeta # v*zeta

  mul!(v.uh, g.rfftplan, v.u) # \hat{u*zeta}
  mul!(v.vh, g.rfftplan, v.v) # \hat{v*zeta}

  @. N = - im*g.kr*v.uh - im*g.l*v.vh
  nothing
end

function calcN!(N, sol, t, cl, v, p, g)
  calcN_advection!(N, sol, t, cl, v, p, g)
  addforcing!(N, sol, t, cl, v, p, g)
  nothing
end

addforcing!(N, sol, t, cl, v::Vars, p, g) = nothing

function addforcing!(N, sol, t, cl, v::ForcedVars, p, g)
  p.calcF!(v.Fh, sol, t, cl, v, p, g)
  @. N += v.Fh
  nothing
end

function addforcing!(N, sol, t, cl, v::StochasticForcedVars, p, g)
  if t == cl.t # not a substep
    @. v.prevsol = sol # sol at previous time-step is needed to compute budgets for stochastic forcing
    p.calcF!(v.Fh, sol, t, cl, v, p, g)
  end
  @. N += v.Fh
  nothing
end

# ----------------
# Helper functions
# ----------------

"""
    updatevars!(prob)

Update the vars in v on the grid g with the solution in sol.
"""
function updatevars!(prob)
  v, g, sol = prob.vars, prob.grid, prob.sol
  @. v.zetah = sol
  @. v.uh =   im * g.l  * g.invKrsq * sol
  @. v.vh = - im * g.kr * g.invKrsq * sol
  ldiv!(v.zeta, g.rfftplan, deepcopy(v.zetah))
  ldiv!(v.u, g.rfftplan, deepcopy(v.uh))
  ldiv!(v.v, g.rfftplan, deepcopy(v.vh))
  nothing
end

"""
    set_zeta!(prob, zeta)

Set the solution sol as the transform of zeta and update variables v
on the grid g.
"""
function set_zeta!(prob, zeta)
  p, v, g, sol = prob.params, prob.vars, prob.grid, prob.sol
  mul!(sol, g.rfftplan, zeta)
  sol[1, 1] = 0 # zero domain average
  updatevars!(prob)
  nothing
end

"""
    energy(prob)

Returns the domain-averaged kinetic energy in the Fourier-transformed vorticity
solution `sol`.
"""
@inline function energy(prob)
  sol, v, g = prob.sol, prob.vars, prob.grid
  @. v.uh = g.invKrsq * abs2(sol)
  1/(2*g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

"""
    enstrophy(prob)

Returns the domain-averaged enstrophy in the Fourier-transformed vorticity
solution `sol`.
"""
@inline function enstrophy(prob)
  sol, g = prob.sol, prob.grid
  1/(2*g.Lx*g.Ly)*parsevalsum2(sol, g)
end

"""
    dissipation(prob)

Returns the domain-averaged dissipation rate. nν must be >= 1.
"""
@inline function dissipation(prob)
  sol, v, p, g = prob.sol, prob.vars, prob.params, prob.grid
  @. v.uh = g.Krsq^(p.nν-1) * abs2(sol)
  v.uh[1, 1] = 0
  p.ν/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

"""
    work(prob)
    work(sol, v, g)

Returns the domain-averaged rate of work of energy by the forcing Fh.
"""
@inline function work(sol, v::ForcedVars, g)
  @. v.uh = g.invKrsq * sol * conj(v.Fh)
  1/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

@inline function work(sol, v::StochasticForcedVars, g)
  @. v.uh = g.invKrsq * (v.prevsol + sol)/2 * conj(v.Fh) # Stratonovich
  # @. v.uh = g.invKrsq * v.prevsol * conj(v.Fh)           # Ito
  1/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

@inline work(prob) = work(prob.sol, prob.vars, prob.grid)

"""
    drag(prob)

Returns the extraction of domain-averaged energy by drag/hypodrag μ.
"""
@inline function drag(prob)
  sol, v, p, g = prob.sol, prob.vars, prob.params, prob.grid
  @. v.uh = g.Krsq^(p.nμ-1) * abs2(sol)
  v.uh[1, 1] = 0
  p.μ/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

end # module


text:
using DrWatson
@quickactivate "Medical segmentation evaluation"
"""
module based on https://github.com/Visceral-Project/EvaluateSegmentation/blob/master/source/ContingencyTable.h
"""
module ContTable
using ..BasicStructs, Parameters

"""
giving set of constant associated with  the image that will be usefull for calculating
    I - image the we are trying to segment - 3 dimensional array of supplied type imageNumb
    return ImageConstants - some constants related to image iself (nor masks)
    """
function getImageConstants(::Type{imageNumb}
                      I::Array{imageNumb, 3})
                      ::ImageConstants
                      where{imageNumb}

  return ImageConstants(
            mvspx = # Voxelspacing x
            mvspy= # Voxelspacing y
            mvspz= #mean Voxelspacing z
            isZConst = # true if slices thickness is the same in all image
            numberOfVox= # number of voxels in image

  )

end #getImageConstants

"""
calculating tn, fn, fp and tp WITHOUT voxel volume correction
Type{maskNumb}  - type of the numbers hold in mask
G - 3 dimensional array holding ground truth segmentation
T - 3 dimensional array holding segmentation that we want to compare to ground truth
isVariedSlice- true if slices can have varying thickness
imageConst - image dependent constants
return array of doubles with entry 1)TN- true negative 2)  TP - true positive 3) FN - false negative  4) FP - false positive
"""
function getTnTpFpFn(::Type{maskNumb}
                ,G::Array{maskNumb, 3}
                 ,T::Array{maskNumb, 3}
                  ,isVariedSlice::Bool
                  ,imageConst::ImageConstants
    )::Vector{Float64}
     where{maskNumb}

    x1::Float64 =((double)values_f[i])/((double)PIXEL_VALUE_RANGE_MAX)
    y1::Float64  =1-x1
    x2::Float64 =1-x2
    tn += min(y1,y2)
    tp += min(x1,x2)
    fn += x1>x2 ? x1-x2 : 0
    fp += x2>x1 ? x2-x1 : 0

	look inot https://github.com/JuliaArrays/StructArrays.jl

return [tn,tp,fp,fn ]
end   #TnTpFpFn

"""
calculating  a,b,c,d constants - needed for some particulary pairwise comparison metrics WITHOUT voxel volume correction
isVariedSlice- true if slices can have varying thickness
imageConst - image dependent constants
TnTpFpFns - calculated Tn Tp Fp and Fn in that order
return array of doubles with entry 1) a  2)  b 3) c  4) FP d
"""

function getAbcdConsts(isVariedSlice::Bool
                  ,imageConst::ImageConstants
                  ,TnTpFnFps::Vector{Bool}
    )::Vector{Float64}
	n = std::min(numberElements_f, numberElements_m);
	coltot1 = tn + fp;
	coltot2 = fn + tp;
	rowtot1 = tn + fn;
	rowtot2 = fp + tp;
    nis = rowtot1*rowtot1 + rowtot2*rowtot2;
    njs = coltot1*coltot1 + coltot2*coltot2;
    s = tp*tp + tn*tn + fp*fp + fn*fn ;

    a = ( binomial(tn) + binomial(fn) + binomial(fp) + binomial(tp) )/2.0;

    b = (njs - s)/2.0;
    c = (nis - s)/2.0;
    d = ( n*n + s - nis - njs )/2.0;

return [a,b,c,d]

end   #getAbcdConsts

"""
helper function calculating binomial
"""
function binomial(val::Float64)::Float64
    return  val*(val-1)
end #binomial

end #ContTable


text:
"""
Write out to a file a visualization for a WaterModels network dictionary parsed from an
EPANET file. `basefilename` should not include an extension and will be appended with
`_w_cb.pdf` in the final output file, which is a multi-page PDF. The `layout` option equates
to the layout functions of graphviz (dot, neato, etc.). Use `del_files=false` to keep the
intermediate files. If a `solution` dict is provided, it should be for the same time as that
of the `data` dict.
"""
function write_visualization(data::Dict{String,Any}, basefilename::String,
                             solution::Union{Nothing, Dict{String,Any}}=nothing;
                             layout::String="dot", args::String="",
                             sep_page::Bool=false, del_files::Bool=true)

    gpdffile = basefilename*"_graph.pdf"
    cbfile = basefilename*"_cbar.pdf"
    outfile = basefilename*"_graph_w_cb.pdf"

    G = build_graph(data, solution, layout=layout, args=args)
    write_graph(G, gpdffile)
    colorbar(G, cbfile)

    # note that `stack_bar` is a python function from `wntr_vis.py`
    stack_cbar(gpdffile, cbfile, outfile, sep_page)

    # delete the intermediate files
    if del_files
        run(`rm $gpdffile $cbfile`)
    end
end

"""
Create a visualizations for each time. `wmdata` and `solution` are multi-time WM data
objects.
"""
function write_multi_time_viz(wmdata::Dict{String,Any}, solution::Dict{String,Any},
                              basefilename::String; layout::String="dot", args::String="",
                              del_files::Bool=true)
    cbfile = basefilename*"_cbar.pdf"
    outfile = basefilename*"_graph_w_cb.pdf"

    duration = wmdata["duration"]/3600
    step = wmdata["time_step"]/3600
    tarr = collect(range(1, stop=duration, step=step))
    tarr = string.(Int.(tarr))
    N = length(tarr)

    k = 0
    filenames = Array{String,1}(undef, N+1)
    filenames[k+=1] = cbfile

    # build graph at t = 1
    t = tarr[1]
    G = build_graph(wmdata["nw"][t], solution["nw"][t], layout=layout, args=args)
    filename = "$(basefilename)_$(t)_graph.pdf"
    write_graph(G, filename)
    filenames[k+=1] = filename

    # loop over remaining timepoints
    for t in tarr[2:end]
        update_graph!(G, wmdata["nw"][t], solution["nw"][t])
        filename = "$(basefilename)_$(t)_graph.pdf"
        write_graph(G, filename)
        filenames[k+=1] = filename
    end

    colorbar(G, cbfile)

    # collate the pages
    collate_viz(filenames, outfile)

    # delete files
    if del_files
        for filename in filenames
            run(`rm $filename`)
        end
    end
end

"""
Use graphviz (via pygraphviz) to output a visualization to a file for a graph.

The layout may be redone here but will add computational cost. Note that absolute coordinate
layouts (`neato -n` will likely not work here as the coordinates were overwritten during the
call to `build_graph`.
"""
function write_graph(G::PyCall.PyObject, filename::String,
                     layout::Union{Nothing, Dict{String, <:Any}} = nothing,
                     args::Union{Nothing, Dict{String, <:Any}} = nothing)

    if layout == nothing
        G.draw(filename) # uses layout that was already created in `build_graph`
        if args != nothing
            @warn "args are ignored if a layout is not given"
        end
    else
        try
            G.draw(filename, prog=layout)
        catch
            G.draw(filename, prog="dot")
            # by putting this warning after, it should display only if `G.draw()` was
            # actually successful using "dot"; if it fails for another reason, that error
            # should be shown
            @warn "$layout is not a supported layout; dot was used instead"
        end
    end
end

"""
create and save a colorbar that represents the node elevations
"""
function colorbar(G::PyCall.PyObject, filename::String)
    elmin = parse(Float64, get(G.graph_attr, "elmin", 0.0))
    elmax = parse(Float64, get(G.graph_attr, "elmax", 1.0))
    elmid = elmin + 0.5 * (elmax - elmin)

    x = reshape(collect(range(0.0, stop=1.0, length=100)), (1,:))
    Plots.heatmap(x, c=:viridis, size=(500,100), legend=:none, yaxis=false)
    Plots.plot!(xticks=(0:50:100, [elmin, elmid, elmax]))
    Plots.plot!(yticks=false) # a regression requires this for GR, JJS 1/4/21,
                              # https://github.com/JuliaPlots/Plots.jl/issues/3019
    Plots.title!("Elevation")
    Plots.savefig(filename)
end


text:
using LinearAlgebra, DelimitedFiles, StatsBase

crabs = vec(readdlm("data_7.dat", ',', Int))

max_pos = maximum(crabs)
min_pos = minimum(crabs)
dists = Int[]
fuels = Int[]
for p in min_pos:max_pos
    dist = broadcast(abs, crabs .- p)
    fuel = 0
    for d in dist
        fuel += (d * (d + 1)) / 2
    end
    append!(fuels, fuel)
end

println(minimum(fuels))


text:
export ∅, 𝔻, NaturalNumber, @natural

@symbols ϕ
∅ = ExtensionalSet(Set([]))
𝔻 = IntensionalSet(ϕ, ¬(ϕ ∈ ∅))

# von Neumann ordinals
function NaturalNumber(x::Int)
    if x == 0
        return ∅
    end

    ExtensionalSet([NaturalNumber(n) for n ∈ 0:(x-1)])
end

macro natural(assignment)
    @assert assignment.head == :(=) "Operator must be the assignment operator (=)"
    quote
        $(esc(assignment.args[1])) = NaturalNumber($(assignment.args[2]))
    end
end

# special definitions
# natural numbers:
#  ℕ = ...


text:
# iterate over vals array to find the most generic non-NA type
function typeloop(vals)
    toptype = None
    for i = 1:length(vals)
        if vals[i] != :NA
            toptype = promote_type(toptype, typeof(vals[i]))
        end
    end
    return toptype
end

# Enumerate a BitArray. This is faster than using enumerate()
macro bitenumerate(ba, i, x, code)
    esc(quote
        $i = 1
        for $x in $ba
            $code
            $i += 1
        end
    end)
end


text:
struct ReidIdealParam <: EoSParam
    a::SingleParam{Float64}
    b::SingleParam{Float64}
    c::SingleParam{Float64}
    d::SingleParam{Float64}
end

abstract type ReidIdealModel <: IdealModel end
@newmodelsimple ReidIdeal ReidIdealModel ReidIdealParam

export ReidIdeal
function ReidIdeal(components::Array{String,1}; userlocations::Array{String,1}=String[], verbose=false)
    params = getparams(components, ["ideal/ReidIdeal.csv"]; userlocations=userlocations, verbose=verbose)
    a = params["a"]
    b = params["b"]
    c = params["c"]
    d = params["d"]
    packagedparams = ReidIdealParam(a, b, c, d)
    references = String[] #  Fill this up.
    return ReidIdeal(packagedparams; references=references)
end

function a_ideal(model::ReidIdealModel, V, T, z)
    x = z/sum(z)
    a = model.params.a.values
    b = model.params.b.values
    c = model.params.c.values
    d = model.params.d.values
    polycoeff = [a, b, c, d]
    return 1/sum(x[i]*(log(z[i]/V) + 1/(R̄*T)*(sum(polycoeff[k][i]/k*(T^k-298^k) for k in 1:4)) -
        1/R̄*((a[i]-R̄)*log(T/298)+sum(polycoeff[k][i]/(k-1)*(T^(k-1)-298^(k-1)) for k in 2:4))) for i in @comps)
end


text:
module Config

end


text:
using Test
import LanguageDetect: detect

@test detect("最近怎样呢？")[1].language == "zh-cn"
@test detect("こんにちは")[1].language == "ja"
@test detect("안녕하세요")[1].language == "ko"


text:
#-------------fft_slow.jl------------------------------------------------------#
#
# Purpose: Implement a slow version of an fft, to see how it's done.
#          For this, we will be using the Cooley-Tukey algorithm
#
#------------------------------------------------------------------------------#

using PyPlot

# Implementing the Cooley-Tukey Algorithm
function cooley_tukey(x)
    N = length(x)

    #println(N)

    if(N%2 !=0)
        #println("Must be a power of 2!")
        exit(0)
    end
    if(N <= 2)
        #println("DFT_slow")
        return DFT_slow(x)
    else
        x_even = cooley_tukey(x[1:2:N])
        x_odd = cooley_tukey(x[2:2:N])
        n = 0:N-1
        #n = n'
        half = div(N,2)
        factor = exp(-2im*pi*n/N)
        #println(length(factor))
        #println(length(x_even))
        #println(length(x_odd))
        return vcat(x_even + factor[1:half] .* x_odd,
                     x_even + factor[half+1:N] .* x_odd)
    end

end

# A super slow DFT to start
function DFT_slow(x)
    N = length(x)

    # range is automatically a column vector, and we want a row vector
    n = 0:N-1
    n = n'
    k = n'
    M = exp(-2im * pi *k *n / N)
    return M * x
end

# function to greyscale image
function greyscale(img)
    a = zeros(width(img), height(img))
    for i = 1:width(img)
        for j = 1:height(img)
            a[i,j] = (img[i,j].r + img[i,j].b + img[i,j].g) / 3
        end
    end

    return a
end

# function to determine whether a point is inside of a decided norm
function in_norm(x, y, power)
    if (abs(x)^power + abs(y)^power)^(1/power) < 1
        return true
    else
        return false
    end
end

# function to create array of points to work with
function monte_carlo(res, power)
    a = []
    for i = 1:res
        x = rand() * 2 - 1
        y = rand() * 2 - 1
        if (in_norm(x, y, power))
            push!(a, [x,y])
        end
    end

    x = map(x -> x[1], a)
    y = map(y -> y[2], a)
    scatter(x,y)
    #return x,y
end

max_issues_repo_name: shiffman/simuleios

max_forks_repo_name: shiffman/simuleios


text:
# This file was generated, do not modify it. # hide
using StatsPlots

times_d = [dμσ, dμ1, d0σ, d01]
times_m = [mμσ, mμ1, m0σ, m01]

times = [times_d ;times_m] / 1e3 * 1e9

pkg = repeat(["Distributions.jl", "MeasureTheory.jl"], inner=4)
dist = repeat(["Normal(μ,σ)", "Normal(μ,1)", "Normal(0,σ)", "Normal(0,1)"], outer=2)
groupedbar(dist, times, group=pkg, legend=:topleft)
ylabel!("Time per element (ns)")
savefig(joinpath(@OUTPUT, "dists-measuretheory-times.svg")) # hide

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
# ------------ Dynamic map ------------ #

struct DynamicMap{K} <: AddressMap{K}
    tree::Dict{Any, AddressMap{<:K}}
    DynamicMap{K}() where K = new{K}(Dict{Any, AddressMap{K}}())
    DynamicMap{K}(tree::Dict{Any, AddressMap{K}}) where K = new{K}(tree)
end
DynamicMap(tree::Dict{Any, AddressMap{K}}) where K = DynamicMap{K}(tree)
Zygote.@adjoint DynamicMap(tree) = DynamicMap(tree), ret_grad -> (nothing, )

@inline shallow_iterator(dm::DynamicMap) = dm.tree

@inline function get(dm::DynamicMap{Value}, addr, fallback)
    haskey(dm, addr) || return fallback
    return getindex(dm, addr)
end

@inline function get_sub(dm::DynamicMap, addr::A) where A <: Address
    haskey(dm.tree, addr) && return getindex(dm.tree, addr)
    Empty()
end
@inline get_sub(dm::DynamicMap, addr::Tuple{}) = Empty()

@inline getindex(dm::DynamicMap, addrs...) = get_value(get_sub(dm, addrs))

@inline Base.isempty(dm::DynamicMap) = isempty(dm.tree)

# This is a fallback for subtypes of DynamicMap.
@inline haskey(dm::DynamicMap, addr::A) where A <: Address = haskey(dm.tree, addr) && has_value(get_sub(dm, addr))

# This is a fallback for subtypes of DynamicMap.
@inline has_sub(dm::DynamicMap, addr::A) where A <: Address = haskey(dm.tree, addr)

function set_sub!(dm::DynamicMap{K}, addr::A, v::AddressMap{<: K}) where {K, A <: Address}
    delete!(dm.tree, addr)
    if !isempty(v)
        dm.tree[addr] = v
    end
end
function set_sub!(dm::DynamicMap{K}, addr::A, v::V) where {K, V <: Union{Value, Choice}, A <: Address}
    delete!(dm.tree, addr)
    if !isempty(v)
        dm.tree[addr] = convert(K, v)
    end
end
@inline set_sub!(dm::DynamicMap{K}, addr::Tuple{A}, v::AddressMap{<: K}) where {A <: Address, K} = set_sub!(dm, addr[1], v)
function set_sub!(dm::DynamicMap{K}, addr::Tuple, v::AddressMap{<:K}) where K
    hd, tl = addr[1], addr[2 : end]
    if !haskey(dm.tree, hd)
        dm.tree[hd] = DynamicMap{K}()
    end
    set_sub!(dm.tree[hd], tl, v)
end

function merge(sel1::DynamicMap{K},
               sel2::DynamicMap{K}) where K
    tree = Dict{Any, AddressMap{K}}()
    inter = intersect(keys(sel1.tree), keys(sel2.tree))
    check = false
    for k in setdiff(keys(sel2.tree), keys(sel1.tree))
        tree[k] = sel2.tree[k]
    end
    for k in setdiff(keys(sel1.tree), keys(sel2.tree))
        tree[k] = sel1.tree[k]
    end
    for k in inter
        tree[k], b = merge(sel1.tree[k], sel2.tree[k])
        check = check || b
    end
    return DynamicMap(tree), !isempty(inter) || check
end
function merge(sel1::DynamicMap{T},
               sel2::DynamicMap{K}) where {T, K}
    tree = Dict{Any, AddressMap{K}}()
    inter = intersect(keys(sel1.tree), keys(sel2.tree))
    check = false
    for k in setdiff(keys(sel2.tree), keys(sel1.tree))
        tree[k] = sel2.tree[k]
    end
    for k in setdiff(keys(sel1.tree), keys(sel2.tree))
        tree[k] = convert(K, sel1.tree[k])
    end
    for k in inter
        tree[k], b = merge(sel1.tree[k], sel2.tree[k])
        check = check || b
    end
    return DynamicMap(tree), !isempty(inter) || check
end
merge(dm::DynamicMap, ::Empty) = deepcopy(dm), false
merge(::Empty, dm::DynamicMap) = deepcopy(dm), false

function merge!(sel1::DynamicMap{K},
                sel2::DynamicMap{K}) where K
    inter = intersect(keys(sel1.tree), keys(sel2.tree))
    check = false
    for k in setdiff(keys(sel2.tree), keys(sel1.tree))
        set_sub!(sel1, k, get_sub(sel2, k))
    end
    for k in inter
        b = merge!(get_sub(sel1, k), get_sub(sel2, k))
        check = check || b
    end
    !isempty(inter) || check
end
function merge!(sel1::DynamicMap{T},
                sel2::DynamicMap{K}) where {T, K}
    inter = intersect(keys(sel1.tree), keys(sel2.tree))
    check = false
    for k in setdiff(keys(sel2.tree), keys(sel1.tree))
        set_sub!(sel1, k, get_sub(sel2, k))
    end
    for k in inter
        b = merge!(get_sub(sel1, k), get_sub(sel2, k))
        check = check || b
    end
    !isempty(inter) || check
end
merge!(dm::DynamicMap, ::Empty) = dm, false
merge!(::Empty, dm::DynamicMap) = dm, false
Zygote.@adjoint merge!(a, b) = merge!(a, b), s_grad -> (nothing, nothing)

# Used for pretty printing.
function collect!(par::T, addrs::Vector{Any}, chd::Dict{Any, Any}, tr::D, meta) where {T <: Tuple, D <: DynamicMap}
    iterate(tr) do (k, v)
        collect!((par..., k), addrs, chd, v, meta)
    end
end
function collect!(addrs::Vector{Any}, chd::Dict{Any, Any}, tr::D, meta) where D <: DynamicMap
    iterate(tr) do (k, v)
        collect!((k, ), addrs, chd, v, meta)
    end
end

# Utility constructor.
@inline target() = DynamicMap{Value}()
function target(v::Vector{Pair{T, K}}) where {T <: Tuple, K}
    #check = false
    #for (k, _) in v
    #    if length(k) == 1 && k[1] isa Symbol
    #        check = true
    #    end
    #end
    #check && begin
    #    @info "Converting to static choice map."
    #    return static(v)
    #end
    tg = DynamicMap{Value}()
    for (k, v) in v
        set_sub!(tg, k, Value(v))
    end
    tg
end
function target(v::Pair{T, K}) where {T <: Tuple, K}
    k, _ = v
    #length(k) == 1 && k[1] isa Symbol && begin
    #    @info "Converting to static choice map."
    #    return static(v)
    #end
    tg = DynamicMap{Value}()
    set_sub!(tg, v[1], Value(v[2]))
    tg
end

# Filter.
function filter(fn, dm::DynamicMap{K}) where K
    new = DynamicMap{K}()
    for (k, v) in shallow_iterator(dm)
        if fn((k, ))
            set_sub!(new, k, v)
        else
            ret = filter(fn, (k, ), v)
            !isempty(ret) && set_sub!(new, k, filter(fn, (k, ), v))
        end
    end
    new
end
function filter(fn, par, dm::DynamicMap{K}) where K
    new = DynamicMap{K}()
    for (k, v) in shallow_iterator(dm)
        if fn((par..., k))
            set_sub!(new, k, v)
        else
            ret = filter(fn, (par..., k), v)
            !isempty(ret) && set_sub!(new, k, filter(fn, (par..., k), v))
        end
    end
    new
end

# ------------ Conversion to Selection ------------ #

function select(dm::DynamicMap{K}) where K
    new = DynamicTarget()
    for (k, v) in shallow_iterator(dm)
        set_sub!(new, k, select(v))
    end
    new
end

# ------------ Convert from Dict ------------ #

@inline function convert(::Type{DynamicMap{Value}}, d::Dict{Symbol, V}) where V
    target(map(collect(d)) do (k, v)
               (k, ) => v
           end)
end


text:
export MultiVector
export localLength, globalLength, numVectors
export scale, scale!
export getVectorView, getVectorCopy, getLocalArray
export commReduce

using LinearAlgebra

# otherwise, the use is forced to import LinearAlgebra
export dot, norm, copyto!

# need to add methods to these functions
import LinearAlgebra: dot, norm

"""
`MultiVector`s represent a group of vectors to be processed together.
They are a subtype of [`AbstractArray{Data, 2}`] and support the [`DistObject`], and [`SrcDistObject`] for transfering between any two `MultiVectors`.
Required methods:

    getMap(::MultiVector)
    numVectors(::MultiVector)
    getLocalArray(::MultiVector{Data})::AbstractMatrix{Data}
    similar(::MultiVector{Data})

`commReduce(::MultiVector)` may need to be overridden if `getLocallArray(multiVector)` doesn't return a type useable by `sumAll`.

See [`DenseMultiVector`] for a concrete implementation.
"""
abstract type MultiVector{Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer} <: AbstractArray{Data, 2}
end

"""
    globalLength(::MultiVector{Data, GID, PID, LID})::GID

Returns the global length of the vectors in the mutlivector
"""
globalLength(mVect::MultiVector) = numGlobalElements(getMap(mVect))

"""
    localLength(::MultiVector{Data, GID, PID, LID})::LID

Returns the local length of the vectors in the MultiVector
"""
localLength(mVect::MultiVector) = numMyElements(getMap(mVect))

function Base.fill!(mVect::MultiVector, values)
    fill!(getLocalArray(mVect), values)
    mVect
end

function scale!(mVect::MultiVector, alpha::Number)
    rmul!(getLocalArray(mVect), alpha)
    mVect
end

function scale!(mVect::MultiVector{Data, GID, PID, LID}, alpha::AbstractArray{<:Number, 1}) where {Data, GID, PID, LID}
    for v in LID(1):numVectors(mVect)
        @inbounds getVectorView(mVect, v)[:] *= alpha[v]
    end
    mVect
end

function dot(vect1::MultiVector{Data, GID, PID, LID}, vect2::MultiVector{Data, GID, PID, LID}
        )::AbstractArray{Data, 2} where {Data, GID, PID, LID}
    numVects = numVectors(vect1)
    length = localLength(vect1)
    @boundscheck if numVects != numVectors(vect2)
        throw(InvalidArgumentError("MultiVectors must have the same number of vectors to take the dot product of them"))
    end
    @boundscheck if length != localLength(vect2)
        throw(InvalidArgumentError("Vectors must have the same length to take the dot product of them"))
    end
    dotProducts = Array{Data, 2}(undef, 1, numVects)

    data1 = getLocalArray(vect1)
    data2 = getLocalArray(vect2)

    @inbounds for vect in LID(1):numVects
        sum = Data(0)
        for i = LID(1):length
            sum += data1[i, vect]*data2[i, vect]
        end
        dotProducts[vect] = sum
    end

    sumAll(getComm(vect1), dotProducts)::Array{Data, 2}
end

function norm(mVect::MultiVector{Data, GID, PID, LID}, n::Real) where {Data, GID, PID, LID}
    numVects = numVectors(mVect)
    localVectLength = localLength(mVect)
    norms = Array{Data, 2}(undef, 1, numVects)

    data = getLocalArray(mVect)

    if n == 2
        @inbounds for vect in LID(1):numVects
            sum = Data(0)
            for i = LID(1):localVectLength
                val = data[i, vect]
                sum += val*val
            end
            norms[vect] = sum
        end

        norms = sumAll(getComm(getMap(mVect)), norms)::Matrix{Data}
        @. norms = sqrt(norms)
    else
        @inbounds for vect in LID(1):numVects
            sum = Data(0)
            for i = LID(1):localVectLength
                sum += data[i, vect]^n
            end
            norms[vect] = sum
        end

        norms = sumAll(getComm(getMap(mVect)), norms)::Matrix{Data}
        @. norms = norms^(1/n)
    end
end

"""
    commReduce(::MultiVector)

Elementwise reduces the content of the MultiVector across all processes.
Note that the MultiVector cannot be distributed globally.
"""
function commReduce(mVect::MultiVector{Data}) where Data
    #can only reduce locally replicated mutlivectors
    if distributedGlobal(mVect)
        throw(InvalidArgumentError("Cannot reduce distributed MultiVector"))
    end
    view = getLocalArray(mVect)::AbstractMatrix{Data}
    view .= sumAll(getComm(mVect), view)
end

"""
    getVectorView(::DenseMultiVector{Data}, columns)::AbstractArray{Data}

Gets a view of the requested column vector(s) in this DenseMultiVector
"""
getVectorView(mVect::MultiVector, column) = view(getLocalArray(mVect), :, column)

"""
    getVectorCopy(::MultiVector{Data}, columns)::Array{Data}

Gets a copy of the requested column vector(s) in this MultiVector
"""
function getVectorCopy(mVect::MultiVector{Data}, column)::Array{Data} where {Data}
    view = getVectorView(mVect, column)
    copyto!(Array{Data}(undef, size(view)), view)
end

#### DistObject Interface ####

function checkSizes(source::MultiVector{Data, GID, PID, LID},
        target::MultiVector{Data, GID, PID, LID})::Bool where {
            Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer}
    (numVectors(source) == numVectors(target)
        && globalLength(source) == globalLength(target))
end

function copyAndPermute(source::MultiVector{Data, GID, PID, LID},
        target::MultiVector{Data, GID, PID, LID}, numSameIDs::LID,
        permuteToLIDs::AbstractArray{LID, 1}, permuteFromLIDs::AbstractArray{LID, 1}
        ) where {Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer}
    numPermuteIDs = length(permuteToLIDs)
    sourceData = getLocalArray(source)
    targetData = getLocalArray(target)
    @inbounds for vect in LID(1):numVectors(source)
        for i in LID(1):numSameIDs
            targetData[i, vect] = sourceData[i, vect]
        end

        #don't need to sort permute[To/From]LIDs, since the orders match
        for i in LID(1):numPermuteIDs
            targetData[permuteToLIDs[i], vect] = sourceData[permuteFromLIDs[i], vect]
        end
    end
end

function packAndPrepare(source::MultiVector{Data, GID, PID, LID},
        target::MultiVector{Data, GID, PID, LID}, exportLIDs::AbstractArray{LID, 1},
        distor::Distributor{GID, PID, LID})::Array where {
            Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer}
    numVects = Int(numVectors(source))
    packAndPrepare_helper(source, target, exportLIDs, distor, Val{numVects})
end
#Use helper function to get numVects as a compile time constant
function packAndPrepare_helper(source::MultiVector{Data, GID, PID, LID},
            target::MultiVector{Data, GID, PID, LID}, exportLIDs::AbstractArray{LID, 1},
            distor::Distributor{GID, PID, LID}, ::Type{Val{numVects}})::Array where {
                Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer, numVects}
    exports = Vector{NTuple{numVects, Data}}(undef, length(exportLIDs))
    sourceData = getLocalArray(source)
    bytes_ptr = convert(Ptr{Data}, Base.unsafe_convert(Ptr{NTuple{numVects, Data}}, exports))
    for i in 1:length(exportLIDs)
        i_base = (i-1)*numVects
        for j in 1:numVects
            @inbounds unsafe_store!(bytes_ptr, sourceData[exportLIDs[i], j], i_base+j)
        end
    end
    exports
end

function unpackAndCombine(target::MultiVector{Data, GID, PID, LID},
        importLIDs::AbstractArray{LID, 1}, imports::AbstractArray,
        distor::Distributor{GID, PID, LID},cm::CombineMode) where {
            Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer}
    numVects = Int(numVectors(target))
    unpackAndCombine_helper(target, importLIDs, imports, distor, cm, Val{numVects})
end
#Use helper function to get numVects as a compile time constant
function unpackAndCombine_helper(target::MultiVector{Data, GID, PID, LID},
        importLIDs::AbstractArray{LID, 1}, imports::AbstractArray,
        distor::Distributor{GID, PID, LID},cm::CombineMode, ::Type{Val{numVects}}) where {
            Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer, numVects}
    targetData = getLocalArray(target)
    bytes_ptr = convert(Ptr{Data}, Base.unsafe_convert(Ptr{NTuple{numVects, Data}}, imports))
    for i in 1:length(importLIDs)
        i_base = (i-1)*numVects
        for j in 1:numVects
            @inbounds targetData[importLIDs[i], j] = unsafe_load(bytes_ptr, i_base+j)
        end
    end
end

### Julia Array API ###

Base.eltype(::MultiVector{Data}) where Data = Data

Base.size(A::MultiVector) = (Int(globalLength(A)), Int(numVectors(A)))

#TODO this might break for funky maps, however indices needs to return a unit range
Base.axes(A::MultiVector) = (minMyGID(getMap(A)):maxMyGID(getMap(A)), 1:numVectors(A))

function Base.getindex(A::MultiVector, row::Integer, col::Integer)
    @boundscheck begin
        if !(1<=col<=numVectors(A))
            throw(BoundsError(A, (row, col)))
        end
    end

    lRow = lid(getMap(A), row)

    @boundscheck begin
        if lRow < 1
            throw(BoundsError(A, (row, col)))
        end
    end

    @inbounds value = getLocalArray(A)[lRow, col]
    value
end

function Base.getindex(A::MultiVector, i::Integer)
    if numVectors(A) != 1
        throw(ArgumentError("Can only use single index if there is just 1 vector"))
    end

    lRow = lid(getMap(A), i)

    @boundscheck begin
        if lRow < 1
            throw(BoundsError(A, I))
        end
    end

    @inbounds value = getLocalArray(A)[lRow, 1]
    value
end

function Base.setindex!(A::MultiVector, v, row::Integer, col::Integer)
    @boundscheck begin
        if !(1<=col<=numVectors(A))
            throw(BoundsError(A, (row, col)))
        end
    end

    lRow = lid(getMap(A), row)

    @boundscheck begin
        if lRow < 1
            throw(BoundsError(A, (row, col)))
        end
    end

    @inbounds getLocalArray(A)[lRow, 1] = v
    v
end

function Base.setindex!(A::MultiVector, v, i::Integer)
    if numVectors(A) != 0
        throw(ArgumentError("Can only use single index if there is just 1 vector"))
    end

    lRow = lid(getMap(A), i)

    @boundscheck begin
        if lRow < 1
            throw(BoundsError(A, I))
        end
    end

    @inbounds getLocalArray(A)[lRow, 1] = v
    v
end

import Base: ==

function ==(A::MultiVector, B::MultiVector)
    localEquality = localLength(A) == localLength(B) &&
                    numVectors(A) == numVectors(B) &&
                    getLocalArray(A) == getLocalArray(B) &&
                    sameAs(getMap(A), getMap(B))
    minAll(getComm(A), localEquality)
end

struct MultiVectorBroadcastStyle <: Broadcast.AbstractArrayStyle{2} end
Base.BroadcastStyle(::Type{<:MultiVector}) = MultiVectorBroadcastStyle()
Base.BroadcastStyle(::MultiVectorBroadcastStyle, ::Broadcast.AbstractArrayStyle) = MultiVectorBroadcastStyle()
Base.Broadcast.BroadcastStyle(::MultiVectorBroadcastStyle, ::Broadcast.DefaultArrayStyle) = MultiVectorBroadcastStyle()

"`A = find_mv(As)` returns the first MultiVector among the arguments."
find_mv(bc::Base.Broadcast.Broadcasted) = find_mv(bc.args)
find_mv(args::Tuple) = find_mv(find_mv(args[1]), Base.tail(args))
find_mv(x) = x
find_mv(a::MultiVector, rest) = a
find_mv(::Any, rest) = find_mv(rest)

@inline function Broadcast.instantiate(bc::Broadcast.Broadcasted{MultiVectorBroadcastStyle})
    bc
end

@inline function Base.copy(bc::Broadcast.Broadcasted{MultiVectorBroadcastStyle})
    flattened = Broadcast.flatten(bc)
    args = map(mv->if isa(mv, MultiVector) getLocalArray(mv) else mv end, flattened.args)
    result = broadcast(flattened.f, args...)
    mv = find_mv(flattened.args)
    DenseMultiVector(getMap(mv), result)
end

@inline function Base.copyto!(dest::MultiVector, bc::Broadcast.Broadcasted{MultiVectorBroadcastStyle})
    flattened = Broadcast.flatten(bc)
    args = map(mv->if isa(mv, MultiVector) getLocalArray(mv) else mv end, flattened.args)
    broadcast!(flattened.f, getLocalArray(dest), args...)
    flattened.args[1]
end

#### Required Method documentation stubs ####

"""
    numVectors(::MultiVector{Data, GID, PID, LID})::LID

Returns the number of vectors in this `MultiVector`
"""
function numVectors end

"""
    getLocalArray(::MultiVector{Data})::AbstractMatrix{Data}

Returns the array holding the `MultiVector`'s local elements.
Changes to the array content are be reflected in the `MultiVector`
"""
function getLocalArray end


text:
using Pkg

Pkg.develop(PackageSpec(path=joinpath(dirname(@__FILE__), "..")))
# when first running instantiate
Pkg.instantiate()

using Documenter, ExaPF

makedocs(
    sitename = "ExaPF.jl",
    format = Documenter.HTML(
        prettyurls = Base.get(ENV, "CI", nothing) == "true",
        mathengine = Documenter.KaTeX()
    ),
    modules = [ExaPF],
    repo = "https://github.com/exanauts/ExaPF.jl/blob/{commit}{path}#{line}",
    strict = true,
    checkdocs = :exports,
    pages = [
        "Home" => "index.md",
        "Quick start" => "quickstart.md",
        "Manual" => [
            "AutoDiff" => "man/autodiff.md",
            "Benchmark" => "man/benchmark.md",
            "Linear Solver" => "man/linearsolver.md",
            "PowerSystem" => "man/powersystem.md",
            "Formulations" => "man/formulations.md",
            "Evaluators" => "man/evaluators.md",
        ],
        "Library" => [
            "AutoDiff" => "lib/autodiff.md",
            "Linear Solver" => "lib/linearsolver.md",
            "PowerSystem" => "lib/powersystem.md",
            "Formulations" => "lib/formulations.md",
            "Evaluators" => "lib/evaluators.md",
        ]
    ]
)


text:
using Base.Iterators: drop, take
using IterTools: iterated, imap

"""
    fixedpoint(f, guess, stopfun)

return the result of applying `guess = f(guess)`
until convergence. Convergence is decided by applying
`stopfun(guess)` which returns a Boolean.
"""
function fixedpoint(f, guess, stopfun)
    for state in iterated(f, guess)
        stopfun(state) && return state
    end
end

mutable struct StopFunction{T,S}
    oldvals::T
    counter::Int
    tol::S
    maxit::Int
end

"""
    (st::StopFunction)(state)

stopfunction for ctmrg, returning true if singular values are converged or the maximum
number of iterations is reached. Implemented as a closure since it needs to remember
the last singular values it saw for comparison.
"""
function (st::StopFunction)(state)
    st.counter += 1
    st.counter > st.maxit && return true

    vals = state[2]
    diff = norm(vals - st.oldvals)
    diff <= st.tol && return true
    st.oldvals = vals

    return false
end


text:
"""
    `individual_summary(df,xvar,yvar; summary = mean, err = :MOUSE)`
    computes a summary statistic (default = mean) for each individual
    (default = :MOUSE) used to later calculate the error across individual
"""

function individual_summary(df,xvar,yvar; summary = mean, err = :MOUSE)
    #evaluate if it's going over multiple grouping columns
    multiple = typeof(xvar) <: AbstractVector && sizeof(xvar) > 1
    if multiple
        gdc = groupby(df,vcat(xvar,err))
    else
        gdc = groupby(df,[xvar,err])
    end
    df1 = combine(gdc,yvar => summary => yvar)
    sort!(df1,xvar)
    if multiple
        label_df = unique(df1[:,xvar])
    else
        label_df = unique(df1[:,[xvar]])
    end
    label_df[!,:xpos] = 1:nrow(label_df)
    leftjoin(df1,label_df; on = xvar)
    try
        firstval = union(df1[:,xvar])[1]
        df1[!,:xpos] = [v == firstval ? 1 : 2  for v in df1[:,xvar]]
    catch
        println("Can't define x position")
    end
    return df1
end

"""
    `group_summary(df1,xvar,yvar; normality = true)`
    using bootstrap computes central measure plus 95% CI
    if normality is true it uses mean as central measure
    else it uses median
"""

function group_summary(df1,xvar,yvar; normality = true)
    if normality
        central = mean
    else
        central = median
    end
    df2 = combine(groupby(df1,xvar)) do dd
        # using bootstrap to calculate the central point and the 95% CI
        b = bootstrap(central, dd[:,yvar], BasicSampling(100000))
        #bootstrap confint returns a tuple containing a tuple with the statistic
        #estimates and the CI lower and uppper bounds
        m, lower, upper = confint(b,PercentileConfInt(0.95))[1]
        ci1 = m - lower
        ci2 = upper - m
        (Central = m, ERR = (ci1,ci2))
    end
    return df2
end

function survivalrate_algorythm(var; step = 0.5, xaxis = nothing)
    isnothing(xaxis) && (xaxis = range(extrema(var)..., step = step))
    survival = 1 .- ecdf(var).(xaxis)
    return (Xaxis = collect(xaxis), fy = survival)
end

function cumulative_algorythm(var; step = 0.5, xaxis = nothing)
    isnothing(xaxis) && (xaxis = range(extrema(var)..., step = step))
    cum = ecdf(var).(xaxis)
    return (Xaxis = collect(xaxis), fy = cum)
end

function hazardrate_algorythm(var; step = 0.5, xaxis = nothing)
    isnothing(xaxis) && (xaxis = range(extrema(var)..., step = step))
    survival = 1 .- ecdf(var).(xaxis)
    hazard = -pushfirst!(diff(survival),0)./survival
    return (Xaxis = collect(xaxis), fy = hazard)
end

"""
    `function_analysis(df,variable, f; grouping = nothing, step =0.05, calc = :basic,
            color = [:auto], linestyle = [:auto])`
    Apply the function f over the vaariable var per each value of grouping and
    plots the result over the variable var
"""

function function_analysis(df,var, f; grouping = nothing, step = 0.5, calc = :basic,
        xaxis = nothing,
        color = [:auto], linestyle = [:auto])
    subgroups = isnothing(grouping) ? [:MOUSE] : vcat(:MOUSE,grouping)
    isnothing(xaxis) && (xaxis = range(extrema(df[:,var])..., step = step))
    dd1 = combine(groupby(df,subgroups), var => (t-> f(t,xaxis = xaxis, step = step)) => AsTable)
    rename!(dd1, Dict(:Xaxis => var))
    sort!(dd1,[:MOUSE,var])
    if calc == :bootstrapping
        dd2 = combine(groupby(dd1,grouping)) do dd3
            group_summary(dd3,var,:fy; normality = false)
        end
        dd2[!,:low] = [x[1] for x in dd2.ERR]
        dd2[!,:up] = [x[2] for x in dd2.ERR]
    elseif calc == :quantiles
        dd2 = combine(groupby(dd1,[grouping,var]), :fy =>(t-> (Central = mean(t),
        low= abs(mean(t) - quantile(t,0.25),
        up = abs(quantile(t,0.975)-mean(t))),
        # ERR = (abs(mean(t) - quantile(t,0.25)) + abs(quantile(t,0.975)-mean(t)))/2,
        SEM = sem(t))) => AsTable)
    elseif calc == :basic
        dd2 = combine(groupby(dd1,[grouping,var]), :fy =>(t-> (Central = mean(t),up = sem(t), low = sem(t))) => AsTable)
    end
    sort!(dd2,var)

    plt = @df dd2 plot(cols(var),:Central, ribbon = (:low, :up), group = cols(grouping), linecolor = :auto, color = color, linestyle = linestyle)
    return plt, dd2
end


text:
function test_disc_log_picard(P, mP, O::NfOrd)
  # principal ideals should always be invertible
  i = 1
  while i <= 5
    I = O(i)*O
    if !iszero(mP\I)
      return false
    end
    i += 1
  end
  for i = 1:5
    a = rand(O, 10)
    while iszero(a)
      a = rand(O, 10)
    end
    I = a*O
    if !iszero(mP\I)
      return false
    end
  end

  if ngens(P) == 0
    return true
  end

  for i = 1:3
    c = rand(1:10, ngens(P))
    p = P(c)
    I1 = mP(p)
    I2 = mP(P[1])^c[1]
    for j = 2:ngens(P)
      I2 *= mP(P[j])^c[j]
    end
    if mP\I1 != mP\I2 || mP\I1 != p
      return false
    end
  end
  return true
end

function test_disc_log_units(U, mU, O::NfOrd)
  if !iszero(mU\O(1))
    return false
  end

  for i = 1:3
    c = rand(-1:1, ngens(U))
    u = U(c)
    a1 = mU(u)
    a2 = mU(U[1])^c[1]
    for j = 2:ngens(U)
      a2 *= mU(U[j])^c[j]
    end
    if mU\a1 != mU\a2 || mU\a1 != u
      return false
    end
  end
  return true
end

@testset "Picard group and unit group of non maximal orders" begin
  Qx, x = FlintQQ["x"]
  AF = ArbField(20)

  f = x^3 - 2
  K, a = number_field(f, "a", cached = false)
  O = Order(K, [ K(1), 10*a, 100*a^2 ])
  P, mP = @inferred picard_group(O)
  @test is_snf(P)
  @test P.snf == fmpz[ 24 ]
  @test test_disc_log_picard(P, mP, O)

  U, mU = Hecke.unit_group_non_maximal(O)
  @test is_snf(U)
  @test U.snf == fmpz[ 2, 0 ]
  @test contains(AF(53.89509393317), Hecke.regulator([ K(mU(U[2])) ], 1))
  @test test_disc_log_units(U, mU, O)

  f = x^3 - 12*x^2 - 6324*x + 459510
  K, a = number_field(f, "a", cached = false)
  O = equation_order(K)
  P, mP = picard_group(O)
  @test is_snf(P)
  @test P.snf == fmpz[ 3, 6, 6, 18 ]
  @test test_disc_log_picard(P, mP, O)

  U, mU = Hecke.unit_group_non_maximal(O)
  @test is_snf(U)
  @test U.snf == fmpz[ 2, 0 ]
  @test contains(AF(169.7695458895), Hecke.regulator([ K(mU(U[2])) ], 1))
  @test test_disc_log_units(U, mU, O)

  f = x^3-9270*x^2-6226*x-2617
  K, a = number_field(f, "a", cached = false)
  O = equation_order(K)
  P, mP = picard_group(O)
  @test is_snf(P)
  @test P.snf == fmpz[ 2, 6, 24 ]
  @test test_disc_log_picard(P, mP, O)

  U, mU = Hecke.unit_group_non_maximal(O)
  @test is_snf(U)
  @test U.snf == fmpz[ 2, 0 ]
  @test contains(AF(31293.8558289993733), Hecke.regulator([ K(mU(U[2])) ], 1))
  @test test_disc_log_units(U, mU, O)

  f = x^4-3072*x^3+7926*x^2-3920*x-9063
  K, a = number_field(f, "a", cached = false)
  O = equation_order(K)
  P, mP = picard_group(O)
  @test is_snf(P)
  @test P.snf == fmpz[ 2, 2, 2, 2, 4 ]
  @test test_disc_log_picard(P, mP, O)

  U, mU = @inferred Hecke.unit_group_non_maximal(O)
  @test is_snf(U)
  @test U.snf == fmpz[ 2, 0, 0 ]
  @test contains(AF(455982050.1598537651), Hecke.regulator(map( x -> K(mU(x)), [ U[2], U[3] ]), 1))
  @test test_disc_log_units(U, mU, O)

  f = x^3+4064*x^2-1608*x-2816
  K, a = number_field(f, "a", cached = false)
  O = equation_order(K)
  P, mP = @inferred picard_group(O)
  @test is_snf(P)
  @test P.snf == fmpz[ 3, 12 ]
  @test test_disc_log_picard(P, mP, O)

  U, mU = Hecke.unit_group_non_maximal(O)
  @test is_snf(U)
  @test U.snf == fmpz[ 2, 0, 0 ]
  @test contains(AF(124666.2260696), Hecke.regulator(map( x -> K(mU(x)), [ U[2], U[3] ]), 1))
  @test test_disc_log_units(U, mU, O)
end


text:
using AssociativeArray
using Test

include("issetequal.jl")
include("Dic.jl")
include("Arr.jl")
include("indexing.jl")


text:
include("test_IO.jl")
include("test_bin.jl")
include("test_patch.jl")