text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# weak key dictionaries

"""
    WeakKeyDict([itr])

`WeakKeyDict()` constructs a hash table where the keys are weak
references to objects, and thus may be garbage collected even when
referenced in a hash table.

See [`Dict`](@ref) for further help.  Note, unlike [`Dict`](@ref),
`WeakKeyDict` does not convert keys on insertion.
"""
mutable struct WeakKeyDict{K,V} <: AbstractDict{K,V}
    ht::Dict{WeakRef,V}
    lock::ReentrantLock
    finalizer::Function

    # Constructors mirror Dict's
    function WeakKeyDict{K,V}() where V where K
        t = new(Dict{Any,V}(), ReentrantLock(), identity)
        t.finalizer = function (k)
            # when a weak key is finalized, remove from dictionary if it is still there
            if islocked(t)
                finalizer(t.finalizer, k)
                return nothing
            end
            delete!(t, k)
        end
        return t
    end
end
function WeakKeyDict{K,V}(kv) where V where K
    h = WeakKeyDict{K,V}()
    for (k,v) in kv
        h[k] = v
    end
    return h
end
WeakKeyDict{K,V}(p::Pair) where V where K = setindex!(WeakKeyDict{K,V}(), p.second, p.first)
function WeakKeyDict{K,V}(ps::Pair...) where V where K
    h = WeakKeyDict{K,V}()
    sizehint!(h, length(ps))
    for p in ps
        h[p.first] = p.second
    end
    return h
end
WeakKeyDict() = WeakKeyDict{Any,Any}()

WeakKeyDict(kv::Tuple{}) = WeakKeyDict()
copy(d::WeakKeyDict) = WeakKeyDict(d)

WeakKeyDict(ps::Pair{K,V}...)           where {K,V} = WeakKeyDict{K,V}(ps)
WeakKeyDict(ps::Pair{K}...)             where {K}   = WeakKeyDict{K,Any}(ps)
WeakKeyDict(ps::(Pair{K,V} where K)...) where {V}   = WeakKeyDict{Any,V}(ps)
WeakKeyDict(ps::Pair...)                            = WeakKeyDict{Any,Any}(ps)

function WeakKeyDict(kv)
    try
        Base.dict_with_eltype((K, V) -> WeakKeyDict{K, V}, kv, eltype(kv))
    catch
        if !isiterable(typeof(kv)) || !all(x->isa(x,Union{Tuple,Pair}),kv)
            throw(ArgumentError("WeakKeyDict(kv): kv needs to be an iterator of tuples or pairs"))
        else
            rethrow()
        end
    end
end

empty(d::WeakKeyDict, ::Type{K}, ::Type{V}) where {K, V} = WeakKeyDict{K, V}()

islocked(wkh::WeakKeyDict) = islocked(wkh.lock)
lock(f, wkh::WeakKeyDict) = lock(f, wkh.lock)
trylock(f, wkh::WeakKeyDict) = trylock(f, wkh.lock)

function setindex!(wkh::WeakKeyDict{K}, v, key) where K
    !isa(key, K) && throw(ArgumentError("$(limitrepr(key)) is not a valid key for type $K"))
    finalizer(wkh.finalizer, key)
    lock(wkh) do
        wkh.ht[WeakRef(key)] = v
    end
    return wkh
end

function getkey(wkh::WeakKeyDict{K}, kk, default) where K
    return lock(wkh) do
        k = getkey(wkh.ht, kk, secret_table_token)
        k === secret_table_token && return default
        return k.value::K
    end
end

map!(f,iter::ValueIterator{<:WeakKeyDict})= map!(f, values(iter.dict.ht))
get(wkh::WeakKeyDict{K}, key, default) where {K} = lock(() -> get(wkh.ht, key, default), wkh)
get(default::Callable, wkh::WeakKeyDict{K}, key) where {K} = lock(() -> get(default, wkh.ht, key), wkh)
function get!(wkh::WeakKeyDict{K}, key, default) where {K}
    !isa(key, K) && throw(ArgumentError("$(limitrepr(key)) is not a valid key for type $K"))
    lock(() -> get!(wkh.ht, WeakRef(key), default), wkh)
end
function get!(default::Callable, wkh::WeakKeyDict{K}, key) where {K}
    !isa(key, K) && throw(ArgumentError("$(limitrepr(key)) is not a valid key for type $K"))
    lock(() -> get!(default, wkh.ht, WeakRef(key)), wkh)
end
pop!(wkh::WeakKeyDict{K}, key) where {K} = lock(() -> pop!(wkh.ht, key), wkh)
pop!(wkh::WeakKeyDict{K}, key, default) where {K} = lock(() -> pop!(wkh.ht, key, default), wkh)
delete!(wkh::WeakKeyDict, key) = lock(() -> delete!(wkh.ht, key), wkh)
empty!(wkh::WeakKeyDict) = (lock(() -> empty!(wkh.ht), wkh); wkh)
haskey(wkh::WeakKeyDict{K}, key) where {K} = lock(() -> haskey(wkh.ht, key), wkh)
getindex(wkh::WeakKeyDict{K}, key) where {K} = lock(() -> getindex(wkh.ht, key), wkh)
isempty(wkh::WeakKeyDict) = isempty(wkh.ht)
length(t::WeakKeyDict) = length(t.ht)

function iterate(t::WeakKeyDict{K,V}) where V where K
    gc_token = Ref{Bool}(false) # no keys will be deleted via finalizers until this token is gc'd
    finalizer(gc_token) do r
        if r[]
            r[] = false
            unlock(t.lock)
        end
    end
    s = lock(t.lock)
    iterate(t, (gc_token,))
end
function iterate(t::WeakKeyDict{K,V}, state) where V where K
    gc_token = first(state)
    y = iterate(t.ht, tail(state)...)
    y === nothing && return nothing
    wkv, i = y
    kv = Pair{K,V}(wkv[1].value::K, wkv[2])
    return (kv, (gc_token, i))
end

filter!(f, d::WeakKeyDict) = filter_in_one_pass!(f, d)

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

using StatsBase

@testset "histogram addition" begin
	# null histogram
	h1 = fit(Histogram, [0], 0:1)
	h2 = NullHist()
	h3 = h1+h2
	@test h3 == h1
	@test h3 == h2+h1 # commutative

	# histograms with UnitRange edges
	h1 = fit(Histogram, [0], 0:2)
	h2 = fit(Histogram, [0,2], 0:3)
	h3 = h1+h2
	@test h3 == h2+h1 # commutative
	@test h3.edges == h2.edges
	@test h3.weights == [2,0,1] # [1,0] + [1,0,1]

	# histograms with StepRange edges
	h1 = fit(Histogram, [0], 0:2:4)
	h2 = fit(Histogram, [0,5], 0:2:6)
	h3 = h1+h2
	@test h3 == h2+h1 # commutative
	@test h3.edges == h2.edges
	@test h3.weights == [2,0,1] # [1,0] + [1,0,1]

	# histograms with StepRangeLen edges
	h1 = fit(Histogram, [0], 0.0:2.0)
	h2 = fit(Histogram, [0,2], 0.0:3.0)
	h3 = h1+h2
	@test h3 == h2+h1 # commutative
	@test h3.edges == h2.edges
	@test h3.weights == [2,0,1] # [1,0] + [1,0,1]
end

ext: jl


text:
using Documenter
using MCMCChains

DocMeta.setdocmeta!(
    MCMCChains,
    :DocTestSetup,
    :(using MCMCChains);
    recursive=true
)

makedocs(
    sitename = "MCMCChains.jl",
    pages = [
        "MCMCChains" => "index.md",
        "Getting started" => "getting-started.md",
        "Plotting" => [
            "StatsPlots.jl" => "statsplots.md",
            "Gadfly.jl" => "gadfly.md"
        ],
        "API" => [
            "Chains" => "chains.md",
            "Diagnostics" => "diagnostics.md",
            "Posterior statistics" => "stats.md",
            "Model selection" => "modelstats.md",
            "Summarize" => "summarize.md"
        ]
    ],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [MCMCChains],
    strict = true,
    checkdocs = :exports,
    # https://github.com/JuliaLang/julia/pull/37085#issuecomment-683356098
    doctestfilters = [
        r"{([a-zA-Z0-9]+,\s?)+[a-zA-Z0-9]+}",
        r"(Array{[a-zA-Z0-9]+,\s?1}|Vector{[a-zA-Z0-9]+})",
        r"(Array{[a-zA-Z0-9]+,\s?2}|Matrix{[a-zA-Z0-9]+})",
    ]
)

deploydocs(repo = "github.com/TuringLang/MCMCChains.jl.git", push_preview=true)


text:
"""
$(SIGNATURES)

Create and/or update `__site/tag` folders/files. It takes the set of tags to
refresh (in which case only the pages associated to those tags will be
refreshed) or an empty set in which case all tags will be (re)-generated.
Note that if a tag is removed from a page, it may be that `refresh_tags`
contains a tag which should effectively be removed because no page have it.
(This only matters at the end in the call to `write_tag_page`).
Note: the `clean_tags` cleans up orphan tags (tags that wouldn't be pointing
to any page anymore).
"""
function generate_tag_pages(refresh_tags=Set{String}())::Nothing
    # if there are no page tags, cleanup and finish
    PAGE_TAGS = globvar("fd_page_tags")
    isnothing(PAGE_TAGS) && return clean_tags()
    # if there are page tags, eliminiate the pages that don't exist anymore
    for rpath in keys(PAGE_TAGS)
        isfile(rpath * ".md") || delete!(PAGE_TAGS, rpath)
    end
    # if there's nothing left, clean up and finish
    isempty(PAGE_TAGS) && return clean_tags()

    # Here we have a non-empty PAGE_TAGS dictionary where the keys are
    # paths of existing files, and values are the tags associated with
    # these files.
    # Get the inverse dictionary {tag -> [rp1, rp2...]}
    TAG_PAGES = invert_dict(PAGE_TAGS)
    # store it in globvar
    set_var!(GLOBAL_VARS, "fd_tag_pages", TAG_PAGES)
    all_tags = collect(keys(TAG_PAGES))

    # some tags may have been given to refresh which don't have any
    # pages linking to them anymore, these tags will have to be cleaned up
    rm_tags = filter(t -> t ∉ all_tags, refresh_tags)

    # check which tags should be refreshed
    update_tags = isempty(refresh_tags) ? all_tags :
                    filter(t -> t ∈ all_tags, refresh_tags)

    # Generate the tag folder if it doesn't exist already
    isdir(path(:tag)) || mkpath(path(:tag))
    # Generate the tag layout page if it doesn't exist (it should...)
    isfile(joinpath(path(:layout), "tag.html")) || write_default_tag_layout()
    # write each tag, note that they are necessarily in TAG_PAGES
    for tag in update_tags
        write_tag_page(tag)
    end
    return clean_tags(rm_tags)
end

"""
    clean_tags(rm_tags=nothing)

Check the content of the tag folder and remove orphan pages.
"""
function clean_tags(rm_tags=nothing)::Nothing
    isdir(path(:tag)) || return nothing
    PAGE_TAGS = globvar("fd_page_tags")
    if isnothing(PAGE_TAGS) || isempty(PAGE_TAGS)
        all_tags = []
    else
        all_tags = union(values(PAGE_TAGS)...)
    end
    # check what folders are present
    all_tag_folders = readdir(path(:tag))
    # check what folders shouldn't be present
    orphans = setdiff(all_tags, all_tag_folders)
    if !isnothing(rm_tags)
        orphans = union(orphans, rm_tags)
    end
    # clean up
    for dirname in orphans
        dirpath = joinpath(path(:tag), dirname)
        isdir(dirpath) && rm(dirpath, recursive=true)
    end
    return nothing
end

"""
$SIGNATURES

Internal function to (re)write the tag pages corresponding to `update_tags`.
"""
function write_tag_page(tag)::Nothing
    FD_ENV[:SOURCE] = "generated - tag"
    # make `fd_tag` available to that page generation
    set_var!(LOCAL_VARS, "fd_tag", tag)

    layout  = path(:layout)
    content = read(joinpath(layout, "tag.html"), String)

    dir = joinpath(path(:tag), tag)
    isdir(dir) || mkdir(dir)

    h = convert_html(content) |> postprocess_page

    write(joinpath(dir, "index.html"), h)

    # reset `fd_tag`
    set_var!(LOCAL_VARS, "fd_tag", "")
    return nothing
end

"""
    write_default_tag_layout()

If `_layout/tag.html` is not defined, input a basic default one indicating
that the user has to modify it. This will help users transition when they may
have used a template that did not define `_layout/tag.html`.
"""
function write_default_tag_layout()::Nothing
    dc = globvar("div_content")
    dc = ifelse(isempty(dc), globvar("content_class"), dc)
    html = """
        <!doctype html>
        <html lang="en">
        <head>
          <meta charset="UTF-8">
          <meta name="viewport" content="width=device-width, initial-scale=1">
          <title>Tag: {{fill fd_tag}}</title>
        </head>
        <body>
          <div class="$dc tagpage">
            {{taglist}}
          </div>
        </body>
        </html>
        """
    write(joinpath(path(:layout), "tag.html"), html)
    return nothing
end


text:
import Distributions: Poisson

bsfactor = 0.0

struct SingleBinnedSpectra
  f
  p::Parameter
  o::Observable
  df::Symbol
end

"""
    BinnedSpectralMonofit()

Multidimensional fit to one-dimensional distributions (neglecting
correlations between observables).

# Example
```julia-repl
m = BinnedSpectralMonofit()
```
"""
mutable struct BinnedSpectralMonofit <: Model
  params
  pdflist::Array{NLogPDF}
  nll::NLogLikelihood
  # Additional parameters
  spectra::Array{SingleBinnedSpectra}
  datasets::Array{Symbol}
  observables::Array{Observable}
  function BinnedSpectralMonofit(observables::Symbol...)
    new( [], [], NLogLikelihood(),
        Array{SingleBinnedSpectra}(undef, 0), Array{DataFrame}(undef, 0), Array{Observable}(undef, 0) )
  end
end

struct BinnedSpectralPDF
  p::Parameter
  f::String
end

"""
    generate_mock_dataset(m::BinnedSpectralMonofit; kwargs...)

words

# Arguments
- `m::BinnedSpectralMonofit`: Specific model type
- `time::Float64`: Length of the dataset
"""
function generate_mock_dataset(m::BinnedSpectralMonofit; kwargs...)
  kwargs = Dict(kwargs)
  time = get(kwargs, :time, 1.0)
  # Clean out the old data
  for d in m.datasets
    df = eval(d)
    deleterows!(df, 1:size(df, 1))
  end
  # Poisson random number of events
  for p in m.params
    p.fit = rand(Poisson(p.init))
  end
  for sp in m.spectra
    append!(sp.df[!, sp.o.sname],
            randomFromSpectrum(sp.f; count=sp.p.fit, min=sp.o.min, max=sp.o.max))
  end
  # Should global the data
  #for d in m.datasets
  #  add_dataset(d, d)
  #end
end

mutable struct EventCounter
  name::String
  count::Int64
  spectra::Array{SingleBinnedSpectra}
  obs::Dict{Symbol,Array{Float64}}
  function EventCounter(name::String, count::Int64)
    new( name, count, Array{SingleBinnedSpectra}(undef, 0), Dict() )
  end
end

function build_dataframe(ev::EventCounter)
  df = DataFrame()
  for (k,v) in ev.obs
    df[!, k] = v
  end
  df[!, :meta] = repeat([ev.name], ev.count)
  df
end

function compare_and_add!(df1::DataFrame, df2)
  names1 = names(df1)
  names2 = names(df2)
  for n1 in names1
    if !(n1 in names2)
      df2[!, n1] = typeof(df1[!, n1])(undef, size(df2,1))
    end
  end
  for n2 in names2
    if !(n2 in names1)
      df1[!, n2] = typeof(df2[!, n2])(undef, size(df1,1))
    end
  end
  append!(df1, df2)
end

function genmo(m::BinnedSpectralMonofit; kwargs...)
  kwargs = Dict(kwargs)
  time = get(kwargs, :time, 1.0)
  verbose   = get(kwargs, :verbose, false)
  countdict = get(kwargs, :counts, Dict())

  for d in m.datasets
    df = eval(d)
    deleterows!(df, 1:size(df, 1))
    @debug "Deleting from" d
  end
  dictDPCount = Dict{Symbol, Dict{String, EventCounter}}(sp.df=>Dict() for sp in m.spectra)
  for sp in m.spectra
    param = sp.p
    name = param.name
    count_init = get(countdict, name, param.init)
    count::Int64 = rand(Poisson(param.efficiencies[sp.df] * time * (count_init)))
    ### OVERRIDE
    ### count = floor(param.efficiencies[sp.df] * time * param.init)
    ###
    if verbose
      println("Name: ", param.name, " Eff: ", param.efficiencies, " Init: ", param.init, " time ", time, " count ", count)
    end
    if !haskey( dictDPCount[sp.df], name )
      dictDPCount[sp.df][name] = EventCounter(name, count)
    end
    push!( dictDPCount[sp.df][name].spectra, sp )
  end

  for (dset, evdict ) in dictDPCount
    for (pname, ev) in evdict
      count = ev.count
      for sp in ev.spectra
        newdata = randomFromSpectrum(sp.f; count=count, min=sp.o.min, max=sp.o.max)
        ev.obs[sp.o.sname] = newdata
      end
      dfnew = build_dataframe(ev)
      compare_and_add!( getfield(Batman, dset), dfnew )
    end
  end
  updateDataFunctions!(m)
  dictDPCount
end

function updateDataFunctions!(m::BinnedSpectralMonofit)
  for sp in m.spectra
    pname = sp.p.name
    obs = sp.o.name
    fSymbol = Symbol(pname*"_"*obs*"_"*String(sp.df))
    add_array(fSymbol, sp.f( eval(sp.df)[!, sp.o.sname] ))
  end
end

function constructPDF!(m::BinnedSpectralMonofit, p::Parameter, shapes, observables, ds::Symbol)
  if !(ds in m.datasets)
    push!(m.datasets, ds)
  end
  seval = ""
  for (shp, obs) in zip(shapes, observables)
    fSymbol = Symbol(p.name*"_"*String(obs)*"_"* String(ds))
    add_array(fSymbol, shp(eval(ds)[!, obs]))
    #push!(m.spectra, SingleBinnedSpectra(p, shp, obs, df))
    push!(m.spectra, SingleBinnedSpectra(shp, p, observableFromSymbol(m,obs), ds))
    # Lets add an array
    seval *= "$fSymbol.*"
  end
  seval = seval[1:end-2]
  @debug seval
  unique_name = "spectralPDF_"*p.name*"_"*String(ds)
  add_function(unique_name, eval(Meta.parse("()->"*seval)))
  BinnedSpectralPDF(p, unique_name)
end

function combinePDFs!(m::BinnedSpectralMonofit, spectralpdfs::Array{BinnedSpectralPDF},
                      ds::Symbol; kwargs... )
  kwargs = Dict(kwargs)
  livetime = get(kwargs, :livetime, nothing)
  eff = get(kwargs, :efficiency, Array{Nothing}(undef, size(spectralpdfs)))

  seval = "-sum(batlog.("
  fcall = "("
  extended = "+sum("
  ncon = "+( sum("
  plist = []
  flist = []
  if livetime != nothing
    push!(plist, livetime)
    tname = livetime.name
    fcall *= "$tname,"
  end
  for (spdf,e) in zip(spectralpdfs, eff)
    p = spdf.p
    f = spdf.f
    push!(plist, p)
    push!(flist, f)
    vname = p.name
    if e != nothing
      ename = e.name
      seval *= "$ename*"
      fcall *= "$ename,"
      extended *= "$ename*"
      ncon *= "$ename*"
      push!(plist, e)
    end
    if livetime != nothing
      tname = livetime.name
      seval *= "$tname*" ## Is this wrong???!?!?!
      extended *= "$tname*"
      ncon *= "$tname*"
    end
    vinit = p.init
    seval *= "($vname-$bsfactor)*$vinit*$f().+"
    fcall *= "$vname,"
    extended *= "($vname-$bsfactor)*$vinit+"
    ncon *= "($vname-$bsfactor)*$vinit+"
  end
  seval=seval[1:end-2]
  seval*="))"
  fcall=fcall[1:end-1]
  fcall*=")->"
  extended=extended[1:end-1]
  ncon=ncon[1:end-1]
  n = size(eval(ds), 1)
  #extended*=")"*"+$n*log($n)-$n"
  extended*=")*1.0"*"+size($ds,1)*log(size($ds,1))-size($ds,1)"
  ncon *=") - size($ds,1) )^2"
  unique_name = "fullBinnedSpectralPDF_"*String(ds)
  @debug fcall*seval*extended*ncon
  add_function(unique_name, eval(Meta.parse(fcall*seval*extended)))
  push!(m.pdflist, NLogPDF(unique_name, plist...))
end

function add_parameter!(m::BinnedSpectralMonofit, name::String, init; σ=Inf, kwargs... )
  kwargs = Dict(kwargs)
  dsEfficiency = get(kwargs, :dsEfficiency, (Symbol(name)=>1.0))
  p = getparam(m, name)
  if p != nothing
    @warn name*" already exists in model."
    push!( p.efficiencies, dsEfficiency )
    return p
  end
  @debug "Init" init+bsfactor
  p = Parameter(name; init=(init + bsfactor))
  push!( p.efficiencies, dsEfficiency )
  p.samplewidth = σ
  push!(m.params, p)
  if σ != Inf
    #a = Constant(name*"_mu", init+bsfactor)
    a = Constant(name*"_mu", 1.0)
    b = Constant(name*"_sig", σ/init)
    nlp = NLogPDF("lognormal", p, a, b)
    push!(m.pdflist, nlp)
  end
  return p
end

function add_constant!(m::BinnedSpectralMonofit, name::String, init)
  k = getparam(m, name)
  if k != nothing
    @warn name*" already exists in model."
    return k
  end
  Constant(name, init)
end

function add_observable!(m::BinnedSpectralMonofit, name::Symbol, min::Float64, max::Float64)
  push!(m.observables, Observable(String(name), name, min, max))
end

function observableFromSymbol(m::BinnedSpectralMonofit, name::Symbol)
  for obs in m.observables
    if obs.sname == name
      return obs
    end
  end
  return nothing
end

function minimize!( m::BinnedSpectralMonofit; options=Dict() )
  likelihood = NLogLikelihood( m.pdflist )
  add_likelihood!( m, likelihood )
  optimize_model!( m, likelihood; options=options )
end

"""
    randomFromSpectrum(f; min=0.0, max=1.0, top=1.0, count=1)

Return an `Array{Float64}` set of random elements from a 1D pdf given by `f`
between the values of `min` and `max` with a maximum probability of `top`.
"""
function randomFromSpectrum(f; min=0.0, max=1.0, top=10.0, count=1)
  xArr = Array{Float64}(undef, 0)
  for c in collect(1:count)
    while true
      x = rand()*(max - min) + min
      y = rand()*top
      if y < f(x)
        push!(xArr, x)
        break
      end
    end
  end
  return xArr
end


text:
# Unit Tests
using Test
using MLKernels

using MLKernels:
    BaseFunction,
        InnerProduct,
            ScalarProduct,
        PreMetric,
            ChiSquared,
            SineSquared,
            Metric,
                SquaredEuclidean

using SpecialFunctions:
    besselk, gamma

import LinearAlgebra
import Statistics

FloatingPointTypes = (Float32, Float64)
IntegerTypes = (Int32, Int64)
MLK = MLKernels
HP = MLKernels.HyperParameters

include("reference.jl")

include("HyperParameters/hyperparameter.jl")

include("utils.jl")
include("basefunctions.jl")
include("basematrix.jl")

include("kernelfunctions.jl")
include("kernelmatrix.jl")
include("nystrom.jl")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module VirtualNetworkGatewayClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_LocalNetworkGatewaysApi.jl")
include("api_VirtualNetworkGatewayConnectionsApi.jl")
include("api_VirtualNetworkGatewaysApi.jl")

# export models
export convert, AddressSpace
export convert, AuthenticationMethod
export convert, BgpPeerStatus
export convert, BgpPeerStatusListResult
export convert, BgpSettings
export convert, CloudError
export convert, CloudErrorBody
export convert, ConnectionProtocol
export convert, ConnectionResetSharedKey
export convert, ConnectionSharedKey
export convert, DhGroup
export convert, Error
export convert, ErrorDetails
export convert, GatewayRoute
export convert, GatewayRouteListResult
export convert, IPAllocationMethod
export convert, IPConfigurationBgpPeeringAddress
export convert, IkeEncryption
export convert, IkeIntegrity
export convert, IpsecEncryption
export convert, IpsecIntegrity
export convert, IpsecPolicy
export convert, LocalNetworkGateway
export convert, LocalNetworkGatewayListResult
export convert, LocalNetworkGatewayPropertiesFormat
export convert, P2SVpnConnectionRequest
export convert, PfsGroup
export convert, ProvisioningState
export convert, RadiusServer
export convert, Resource
export convert, SubResource
export convert, TagsObject
export convert, TrafficSelectorPolicy
export convert, TunnelConnectionHealth
export convert, VirtualNetworkConnectionGatewayReference
export convert, VirtualNetworkGateway
export convert, VirtualNetworkGatewayConnection
export convert, VirtualNetworkGatewayConnectionListEntity
export convert, VirtualNetworkGatewayConnectionListEntityPropertiesFormat
export convert, VirtualNetworkGatewayConnectionListResult
export convert, VirtualNetworkGatewayConnectionPropertiesFormat
export convert, VirtualNetworkGatewayConnectionStatus
export convert, VirtualNetworkGatewayConnectionType
export convert, VirtualNetworkGatewayIPConfiguration
export convert, VirtualNetworkGatewayIPConfigurationPropertiesFormat
export convert, VirtualNetworkGatewayListConnectionsResult
export convert, VirtualNetworkGatewayListResult
export convert, VirtualNetworkGatewayPropertiesFormat
export convert, VirtualNetworkGatewaySku
export convert, VpnClientConfiguration
export convert, VpnClientConnectionHealthDetail
export convert, VpnClientConnectionHealthDetailListResult
export convert, VpnClientIPsecParameters
export convert, VpnClientParameters
export convert, VpnClientRevokedCertificate
export convert, VpnClientRevokedCertificatePropertiesFormat
export convert, VpnClientRootCertificate
export convert, VpnClientRootCertificatePropertiesFormat
export convert, VpnDeviceScriptParameters
export convert, VpnPacketCaptureStartParameters
export convert, VpnPacketCaptureStopParameters

# export operations
export convert, LocalNetworkGatewaysApi, VirtualNetworkGatewayConnectionsApi, VirtualNetworkGatewaysApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
using Pkg
path = @__DIR__
cd(path)
Pkg.activate(path)

using Mitosis
using MitosisStochasticDiffEq
using LinearAlgebra, Statistics, Random, StatsBase
using DelimitedFiles
using StaticArrays
using NewickTree
using DifferentialEquations
using DiffEqNoiseProcess
using Test

seed = 10
Random.seed!(seed)

const d = 2
const 𝕏 = SVector{d,Float64}
# needed for StochasticDiffEq package where we simulate d+1 states in forwardguiding
# to compute the likelihood.
#const 𝕏_ = SVector{d+1,Float64}
const MS = Mitosis
const MSDE = MitosisStochasticDiffEq

# using MeasureTheory
# import MeasureTheory.logdensity
MS.dim(p::WGaussian{(:μ,:Σ,:c)}) = length(p.μ)
MS._logdet(p::WGaussian{(:μ,:Σ,:c)}) = MS._logdet(p.Σ, MS.dim(p))
MS.whiten(p::WGaussian{(:μ,:Σ,:c)}, x) = MS.lchol(p.Σ)\(x - p.μ)
MS.sqmahal(p::WGaussian, x) = MS.norm_sqr(MS.whiten(p, x))
MS.logdensity(p::WGaussian{(:μ,:Σ,:c)}, x) = p.c - (MS.sqmahal(p,x) + MS._logdet(p) + MS.dim(p)*log(2pi))/2

include("tree.jl")
include("sdetree.jl")

## Read tree
tree = [Tree(S20coal), Tree(S)][1]
#print_tree(tree.newick)

## define model (two traits, as in presentation)
bθ = @SVector [.5, 1.2]
σ0 = @SVector [log(0.2), log(1.05)]
θ0 = (bθ, σ0)

#const M = [-1.0 1.0; 1.0 -1.0]
const M = SMatrix{d,d}(-1.0, 1.0, 1.0, -1.0)
fns(u,θ,t) = (tanh.(Diagonal(θ[1]) * M * u) )
f(u,θ,t) = 𝕏(tanh.(Diagonal(θ[1]) * M * u) )  # f(u,θ,t) = Diagonal(θ[1]) * M * u
gns(u,θ,t) = Diagonal((exp.(θ[2])))
dt0 = 0.001
u0ns = zeros(2)
u0 = zero(𝕏)  # value at root node

## forward sample no tree non-static arrays
ts = 0.0:dt0:1.0
Wsns = cumsum([[zero(u0ns)];[sqrt(ts[i+1]-ti)*randn(size(u0ns))
  for (i,ti) in enumerate(ts[1:end-1])]])
NGns = NoiseGrid(ts,Wsns)
κns = MSDE.SDEKernel(fns, gns, ts, θ0, zeros(d,d))
xT1ns, _ = MSDE.sample(κns, u0ns, MSDE.EulerMaruyama!(), NGns)
xT2ns, _ = MSDE.sample(κns, u0ns, EM(false), NGns)
@test xT1ns ≈ xT2ns

## forward sample no tree static arrays
ts = 0.0:dt0:1.0
Ws = 𝕏.(Wsns)
NG = NoiseGrid(ts,Ws)
κ = MSDE.SDEKernel(f, g, ts, θ0, zeros(d,d))
xT1, _ = MSDE.sample(κ, u0, MSDE.EulerMaruyama!(), NG)
xT2, _ = MSDE.sample(κ, u0, EM(false), NGns) # needs oop
@test xT1 ≈ xT1ns
@test xT1 ≈ xT2

NG = myinnov(ts, 𝕏)
xT1, _ = MSDE.sample(κ, u0, MSDE.EulerMaruyama!(), NG)
xT2, _ = MSDE.sample(κ, u0, EM(false), NG)
@test xT1 ≈ xT2
# dpesn't have same randomness?
Random.seed!(seed)
xT1, (t1,x1,n1) = MSDE.sample(κ, u0, MSDE.EulerMaruyama!())
Random.seed!(seed)
xT2, (t2,x2,n2) = MSDE.sample(κ, u0, EM(false); save_noise=true)
@test xT1 ≈ xT1ns
@test t1 ≈ t2
@test_broken n1 ≈ n2
@test_broken xT1 ≈ xT2

## forward sample on the tree to simulate the observations at the tips

Random.seed!(10)
Xd1, segs1 = forwardsample(tree, u0, θ0, dt0, f, g, MSDE.EulerMaruyama!())
Random.seed!(10)
Xd2, segs2 = forwardsample(tree, u0, θ0, dt0, f, g, EM(false))

@test_broken Xd1 ≈ Xd2
@test segs1[3][1] ≈ segs2[3][1]
@test_broken segs1[3][2][258:260] ≈ segs2[3][2][258:260]
@show segs1[3][2][258:260] - segs2[3][2][258:260]
@show segs1[3][3][258:260] - segs2[3][3][258:260]

# backward filtering
leavescov = inv(10e5) * SA_F64[1 0; 0 1]
Q = [i in tree.lids ? WGaussian{(:μ,:Σ,:c)}(Xd[i], leavescov, 0.0) : missing for i in tree.ids]
θlin = (B(θinit), zeros(d), σ̃(θinit))
dt = 0.01
Q, messages = bwfiltertree!(Q, tree, θlin, dt)

# forward guiding
guidedsegs = Vector{Any}(undef, tree.n) # save all guided segments
X = zeros(𝕏, tree.n)  # values at nodes
for id in tree.lids
    X[id] = Xd[id]
end

Z = [myinnov(messages[i].ts, 𝕏) for i ∈ 2:tree.n]
fwguidtree!(X, guidedsegs, Q, messages, tree, f, g, θ0, Z, EM(false))


text:
# Autogenerated wrapper script for wolfSSL_jll for x86_64-apple-darwin
export libwolfssl

JLLWrappers.@generate_wrapper_header("wolfSSL")
JLLWrappers.@declare_library_product(libwolfssl, "@rpath/libwolfssl.24.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libwolfssl,
        "lib/libwolfssl.24.3.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
####################
### Rule classes ###
####################

struct ContextFreeRule{C1,C2}
    lhs :: C1                # left hand side
    rhs :: Tuple{Vararg{C2}} # right hand side
end

lhs(r::ContextFreeRule) = r.lhs
rhs(r::ContextFreeRule) = r.rhs

function can_create_as_prefix(r::ContextFreeRule, cats)
    @assert length(cats) <= length(r.rhs)
    for (i,cat) in enumerate(cats)
        if cat != r.rhs[i]
            return false
        end
    end
    return true
end

function can_create(r::ContextFreeRule, cats)
    length(r.rhs) == length(cats) && can_create_as_prefix(r,cats)
end
can_create(r::ContextFreeRule, cats...) = can_create(r, cats)

function completion(r::ContextFreeRule, cats)
    @assert cats == r.rhs
    r.lhs
end

isapplicable(r::ContextFreeRule, cat) = lhs(r) == cat
(r::ContextFreeRule)(cat) = rhs(r)

domain(r::ContextFreeRule) = r.lhs

# R is the rule Type
struct MetaRule{R}
    name  :: String
    rules :: Vector{R}
end
name(r::MetaRule) = r.name
show(io::IO, r::MetaRule) = print(io, name(r))

start(mr::MetaRule) = start(mr.rules)
next(mr::MetaRule, state) = next(mr.rules, state)
done(mr::MetaRule, state) = done(mr.rules, state)
length(mr::MetaRule) = length(mr.rules)
eltype(mr::MetaRule) = eltype(mr.rules)

function isapplicable(mr::MetaRule, cat)
    for r in mr.rules
        if isapplicable(r, cat)
            return true
        end
    end
    return false
end

function (mr::MetaRule)(cat)
    for r in mr.rules
        if isapplicable(r, cat) # only one cf-rule can be isapplicable at a time
            return r(cat)
        end
    end
end

function can_create_as_prefix(mr::MetaRule, cats)
    for r in mr.rules
        if can_create_as_prefix(r, cats)
            return true
        end
    end
    return false
end

function can_create(mr::MetaRule, cats)
    for r in mr.rules
        if can_create(r, cats)
            return true
        end
    end
    return false
end

can_create(r::MetaRule, cats...) = can_create(r, cats)

function completion(mr::MetaRule, cats)
    for r in mr.rules
        if can_create(r,cats)
            return completion(r, cats)
        end
    end
    error("No completion possible.")
end

domain(mr::MetaRule) = [domain(r) for r in mr]

###################
### State class ###
###################

# States are hashed by their object ID
mutable struct State{C, CR} # category and rule
    comp    :: Vector{Tuple{C, CR}}  # completions: tuples of categories and rules
    trans   :: Dict{C, State{C, CR}} # possible transitions
    isfinal :: Bool
end

function show(io::IO, state::State)
    print(io ,"State(")
    for comp in state.comp print(io, comp, ", ") end
    for t in state.trans print(io, t, ", ") end
    print(io, state.isfinal, ")")
end

State(C, CR) = State(Vector{Tuple{C, CR}}(), Dict{C, State{C, CR}}(), true)

isfinal(state::State) = state.isfinal
completions(state::State) = state.comp
is_possible_transition(state::State, cat) = !isfinal(state) && haskey(state.trans, cat)
is_possible_transition(grammar, state, cat) = is_possible_transition(state, cat)
transition(state::State, cat) = state.trans[cat]
transition(grammar, state, cat) = transition(state, cat)

function add_rule!(state::State{C,CR}, rule, head, cats) where {C,CR}
    s = state
    for c in cats
        if is_possible_transition(s, c)
            s = transition(s, c)
        else
            s.isfinal = false
            s = s.trans[c] = State(C, CR)
        end
    end
    push!(s.comp, (head, rule))
end

#####################
### Grammar Class ###
#####################

mutable struct Grammar{C, T, CR, TR, S, Cond}
    startstate    :: State{C, CR}
    startsymbols  :: Vector{C}
    terminal_dict :: Dict{T, Vector{Tuple{C, TR}}}
    rule_cond     :: Cond
    prior_cond    :: Cond # needed for GibbsSampling
end

Grammar(startstate    :: State{C, CR},
                            startsymbols  :: Vector{C},
                            terminal_dict :: Dict{T, Vector{Tuple{C, TR}}},
                            rule_cond     :: Cond,
                            prior_cond    :: Cond,
                            Score         :: DataType) where {C, T, CR, TR, Cond} =
    Grammar{C, T, CR, TR, Score, Cond}(startstate,
                                       startsymbols,
                                       terminal_dict,
                                       rule_cond,
                                       prior_cond)

category_type(grammar::Grammar{C, T, CR, TR, S, Cond}) where {C, T, CR, TR, S, Cond} = C
terminal_type(grammar::Grammar{C, T, CR, TR, S, Cond}) where {C, T, CR, TR, S, Cond} = T
category_rule_type(grammar::Grammar{C, T, CR, TR, S, Cond}) where {C, T, CR, TR, S, Cond} = CR
terminal_rule_type(grammar::Grammar{C, T, CR, TR, S, Cond}) where {C, T, CR, TR, S, Cond} = TR
score_type(grammar::Grammar{C, T, CR, TR, S, Cond}) where {C, T, CR, TR, S, Cond} = S
state_type(grammar) = typeof(startstate(grammar))

show(io::IO, grammar::Grammar{C, T, CR, TR, Cond}) where {C, T, CR, TR, Cond} =
    print("Grammar{$C, $T, $(CR), $(TR), $(Cond)}()")

# terminal rules must be unary
function Grammar(category_rules::Vector{MetaRule{ContextFreeRule{C,C}}},
                      terminal_rules::Vector{MetaRule{ContextFreeRule{C,T}}},
                      startsymbols,
                      Score,
                      dependent_components::Function) where {C,T}
    CR = MetaRule{ContextFreeRule{C,C}}
    TR = MetaRule{ContextFreeRule{C,T}}

    # rule_cond  = UrnModelCond(MetaRule[category_rules;terminal_rules], Tuple{ScaleDeg})
    # prior_cond = UrnModelCond(MetaRule[category_rules;terminal_rules], Tuple{ScaleDeg})
    meta_rules = MetaRule[category_rules;terminal_rules]
    cats = union(map(domain, meta_rules)...)
    rule_cond = SimpleCond(
        Dict(
            dependent_components(cat) => DirCat(
                Dict(
                    mr => Float64(isapplicable(mr, cat))
                    for mr in meta_rules
                )
            )
            for cat in cats
        )
    )
    prior_cond = SimpleCond(
        Dict(
            dependent_components(cat) => DirCat(
                Dict(
                    mr => Float64(isapplicable(mr, cat))
                    for mr in meta_rules
                )
            )
            for cat in cats
        )
    )

    terminal_dict = Dict{T, Vector{Tuple{C, TR}}}()
    for mr in terminal_rules
        for r in mr
            t = rhs(r)[1]
            if haskey(terminal_dict, t)
                push!(terminal_dict[t], (lhs(r), mr))
            else
                terminal_dict[t] = [(lhs(r), mr)]
            end
        end
    end

    startstate = State(C, CR)
    for mr in category_rules
        for r in mr
            add_rule!(startstate, mr, lhs(r), rhs(r))
        end
    end

    Grammar(startstate, startsymbols, terminal_dict, rule_cond, prior_cond, Score)
end

function Grammar(category_rules::Vector{ContextFreeRule{C,C}},
                      terminal_rules::Vector{ContextFreeRule{C,T}},
                      startsymbols,
                      Score) where {C,T}
    CR = ContextFreeRule{C,C}
    TR = ContextFreeRule{C,T}

    # rule_cond  = PYPBackoffConditional(0.0, 1.0, ContextFreeRule[category_rules;terminal_rules])
    # prior_cond = PYPBackoffConditional(0.0, 1.0, ContextFreeRule[category_rules;terminal_rules])

    rule_cond =  UrnModelCond(ContextFreeRule[category_rules;terminal_rules], Tuple{ScaleDeg})
    prior_cond = UrnModelCond(ContextFreeRule[category_rules;terminal_rules], Tuple{ScaleDeg})

    terminal_dict = Dict{T, Vector{Tuple{C, TR}}}()
    for r in terminal_rules
        t = rhs(r)[1]
        if haskey(terminal_dict, t)
            push!(terminal_dict[t], (lhs(r), r))
        else
            terminal_dict[t] = [(lhs(r), r)]
        end
    end

    startstate = State(C, CR)
    for r in category_rules
        add_rule!(startstate, r, lhs(r), rhs(r))
    end

    Grammar(startstate, startsymbols, terminal_dict, rule_cond, prior_cond, Score)
end

function Grammar(
        category_rules_string::AbstractString,
        terminal_rules_string::AbstractString,
        startsymbols,
        Score
    )
    Grammar(
        map(x->String.(x), map(split, split(category_rules_string, "\n"))),
        map(x->String.(x), map(split, split(terminal_rules_string, "\n"))),
        startsymbols,
        Score
    )
end

function split_category_from_terminal_rules(rules_string::AbstractString)
    rule_string_lists = map(split, split(rules_string, "\n"))
    terminal_rule_stringlists = Vector{String}[]
    category_rule_stringlists = Vector{String}[]
    for lst in rule_string_lists
        if !isempty(lst)
            if lst[end][1] == '_'
                corrected_lst = [s[1] == '_' ? s[2:end] : s for s in lst]
                push!(terminal_rule_stringlists, corrected_lst)
            else
                push!(category_rule_stringlists, lst)
            end
        end
    end
    category_rule_stringlists, terminal_rule_stringlists
end

function Grammar(rules_string::AbstractString, startsymbols, Score=LogProb)
    category_rules, terminal_rules = split_category_from_terminal_rules(rules_string)
    Grammar(category_rules, terminal_rules, startsymbols, Score)
end

function Grammar(
        category_rules_stringlists::Vector{Vector{S}},
        terminal_rules_stringlists::Vector{Vector{T}},
        startsymbols,
        Score
    ) where {S<:AbstractString, T<:AbstractString}

    if Meta.parse(category_rules_stringlists[1][1]) isa Symbol
        category_rules = [
            ContextFreeRule{String, String}(s[1], Tuple{String, Vararg{String}}(s[2:end]))
            for s in category_rules_stringlists
        ]
        terminal_rules = [
            ContextFreeRule{String, String}(s[1], Tuple{String, Vararg{String}}(s[2:end]))
            for s in terminal_rules_stringlists
        ]
        all_rules = [category_rules; terminal_rules]
        nonterminal_cats = unique(lhs.(all_rules))
        rule_cond = SimpleCond(
            Dict(
                cat => DirCat(
                    Dict(
                        rule => 1.0
                        for rule in all_rules
                        if lhs(rule) == cat
                    )
                )
            for cat in nonterminal_cats
            )
        )
        prior_cond = SimpleCond(
            Dict(
                cat => DirCat(
                    Dict(
                        rule => 1.0
                        for rule in all_rules
                        if lhs(rule) == cat
                    )
                )
            for cat in nonterminal_cats
            )
        )
        terminal_dict = Dict{String, Vector{Tuple{String, ContextFreeRule{String,String}}}}()
        for r in terminal_rules
            t = rhs(r)[1]
            if haskey(terminal_dict, t)
                push!(terminal_dict[t], (lhs(r), r))
            else
                terminal_dict[t] = [(lhs(r), r)]
            end
        end

        startstate = State(String, ContextFreeRule{String,String})
        for r in category_rules
            add_rule!(startstate, r, lhs(r), rhs(r))
        end
    else
        category_rules_with_probs = [
            (ContextFreeRule{String,String}(s[2],Tuple{String, Vararg{String}}(s[3:end])), LogProb(eval(Meta.parse(s[1]))))
            for s in category_rules_stringlists
        ]
        terminal_rules_with_probs = [
            (ContextFreeRule{String,String}(s[2],Tuple{String, Vararg{String}}(s[3:end])), LogProb(eval(Meta.parse(s[1]))))
            for s in terminal_rules_stringlists
        ]
        all_rules_with_probs = [category_rules_with_probs;terminal_rules_with_probs]
        rule_cond = SimpleCond(Dict(
            cat => CatDist(Dict(
                r => p for (r,p) in all_rules_with_probs if lhs(r) == cat
            ))
            for cat in unique(map(x->lhs(x[1]), all_rules_with_probs))
        ))
        prior_cond = SimpleCond(Dict(
            cat => CatDist(Dict(
                r => p for (r,p) in all_rules_with_probs if lhs(r) == cat
            ))
            for cat in unique(map(x->lhs(x[1]), all_rules_with_probs))
        ))

        terminal_dict = Dict{String, Vector{Tuple{String, ContextFreeRule{String,String}}}}()
        for (r,p) in terminal_rules_with_probs
            t = rhs(r)[1]
            if haskey(terminal_dict, t)
                push!(terminal_dict[t], (lhs(r), r))
            else
                terminal_dict[t] = [(lhs(r), r)]
            end
        end

        startstate = State(String, ContextFreeRule{String,String})
        for (r,p) in category_rules_with_probs
            add_rule!(startstate, r, lhs(r), rhs(r))
        end
    end

    Grammar(startstate, startsymbols, terminal_dict, rule_cond, prior_cond, Score)
end

dependent_components(str::String) = str

# Score == LogProb
completions(grammar::Grammar{C, T, CR, TR, LogProb, Cond}, state::State) where {C, T, CR, TR, Cond} =
    ((cat, mr, prob(grammar, cat, mr)) for (cat, mr) in completions(state))
completions(grammar::Grammar{C, T, CR, TR, LogProb, Cond}, t::T) where {C, T, CR, TR, Cond} =
    ((cat, rule, prob(grammar, cat, rule)) for (cat, rule) in grammar.terminal_dict[t])

# Score is a Integer
completions(grammar::Grammar{C, T, CR, TR, S, Cond}, state::State) where {C, T, CR, TR, S <: Integer, Cond} =
    ((cat, mr, one(S)) for (cat, mr) in completions(state))
completions(grammar::Grammar{C, T, CR, TR, S, Cond}, t::T) where {C, T, CR, TR, S <: Integer, Cond} =
    ((cat, rule, one(S)) for (cat, rule) in grammar.terminal_dict[t])

# Score == Product{LogProb,BigInt}
completions(grammar::Grammar{C, T, CR, TR, Product{LogProb,BigInt}, Cond}, state::State) where {C, T, CR, TR, Cond} =
    ((cat, mr, Product(prob(grammar, cat, mr),BigInt(1))) for (cat, mr) in completions(state))
completions(grammar::Grammar{C, T, CR, TR, Product{LogProb,BigInt}, Cond}, t::T) where {C, T, CR, TR, Cond} =
    ((cat, rule, Product(prob(grammar, cat, rule),BigInt(1))) for (cat, rule) in grammar.terminal_dict[t])

startstate(grammar::Grammar) = grammar.startstate
startsymbols(grammar::Grammar) = grammar.startsymbols
prob(grammar::Grammar, cat, rule) =
    logscore(grammar.rule_cond, rule, dependent_components(cat))
priorprob(grammar::Grammar, cat, rule) =
    logscore(grammar.prior_cond, rule, dependent_components(cat))
add_obs!(grammar::Grammar, rule, cat) =
    add_obs!(grammar.rule_cond, rule, dependent_components(cat))
remove_obs!(grammar::Grammar, obs) =
    remove_obs!(grammar.rule_cond, obs)
add_obs!(grammar::Grammar, parse_tree) =
    # Customer[add_obs!(grammar, rule, cat) for (cat, rule) in cats_and_rules(parse_tree)]
    [add_obs!(grammar, rule, cat) for (cat, rule) in cats_and_rules(parse_tree)]
remove_obs!(grammar::Grammar, tree_obs::Vector) =
    for obs in tree_obs remove_obs!(grammar, obs) end
logscore(grammar::Grammar, parse_tree) =
    sum(prob(grammar, cat, rule) for (cat,rule) in cats_and_rules(parse_tree))
priorlogscore(grammar::Grammar, parse_tree) =
    sum(priorprob(grammar, cat, rule) for (cat,rule) in cats_and_rules(parse_tree))


text:
# This file was automatically generated by the ILK-compiler on Wed Apr 10 14:11:52 2019

module ur5

import LinearAlgebra
import ILKBackend

struct ConstPosesContainer
    elbow__upperarm::ILKBackend.MatrixType
    forearm__wr1::ILKBackend.MatrixType
    shoulder__shoulder_lift::ILKBackend.MatrixType
    shoulder_lift__shoulder::ILKBackend.MatrixType
    shoulder_pan__base::ILKBackend.MatrixType
    upperarm__elbow::ILKBackend.MatrixType
    wr1__forearm::ILKBackend.MatrixType
    wr2__wrist_1::ILKBackend.MatrixType
    wr3__wrist_2::ILKBackend.MatrixType
    wrist_1__wr2::ILKBackend.MatrixType
    wrist_2__wr3::ILKBackend.MatrixType
end

function ConstPosesInit()
    elbow__upperarm = ILKBackend.pose()
    ILKBackend.setPosition(elbow__upperarm,0.425,0,-0.1197)
    ILKBackend.setRotation(elbow__upperarm,1,-0,0,0,1,-0,0,0,1)
    forearm__wr1 = ILKBackend.pose()
    ILKBackend.setPosition(forearm__wr1,-0.39225,0,-0.09315)
    ILKBackend.setRotation(forearm__wr1,1,-0,0,0,1,-0,0,0,1)
    shoulder__shoulder_lift = ILKBackend.pose()
    ILKBackend.setPosition(shoulder__shoulder_lift,0,0,-0.13585)
    ILKBackend.setRotation(shoulder__shoulder_lift,0,-0,1,-0,-1,0,1,-0,-0)
    shoulder_lift__shoulder = ILKBackend.pose()
    ILKBackend.setPosition(shoulder_lift__shoulder,0.13585,0,0)
    ILKBackend.setRotation(shoulder_lift__shoulder,0,-0,1,-0,-1,0,1,-0,-0)
    shoulder_pan__base = ILKBackend.pose()
    ILKBackend.setPosition(shoulder_pan__base,0,0,0.089159)
    ILKBackend.setRotation(shoulder_pan__base,1,-0,0,0,1,-0,0,0,1)
    upperarm__elbow = ILKBackend.pose()
    ILKBackend.setPosition(upperarm__elbow,-0.425,0,0.1197)
    ILKBackend.setRotation(upperarm__elbow,1,-0,0,0,1,-0,0,0,1)
    wr1__forearm = ILKBackend.pose()
    ILKBackend.setPosition(wr1__forearm,0.39225,0,0.09315)
    ILKBackend.setRotation(wr1__forearm,1,-0,0,0,1,-0,0,0,1)
    wr2__wrist_1 = ILKBackend.pose()
    ILKBackend.setPosition(wr2__wrist_1,0.09475,0,0)
    ILKBackend.setRotation(wr2__wrist_1,0,-0,1,-0,-1,0,1,-0,-0)
    wr3__wrist_2 = ILKBackend.pose()
    ILKBackend.setPosition(wr3__wrist_2,0.0825,0,0)
    ILKBackend.setRotation(wr3__wrist_2,-0,0,1,0,1,0,-1,0,-0)
    wrist_1__wr2 = ILKBackend.pose()
    ILKBackend.setPosition(wrist_1__wr2,0,0,-0.09475)
    ILKBackend.setRotation(wrist_1__wr2,0,0,1,0,-1,0,1,0,0)
    wrist_2__wr3 = ILKBackend.pose()
    ILKBackend.setPosition(wrist_2__wr3,0,0,-0.0825)
    ILKBackend.setRotation(wrist_2__wr3,0,0,-1,0,1,0,1,0,-0)

    ret = ConstPosesContainer(
        elbow__upperarm,
        forearm__wr1,
        shoulder__shoulder_lift,
        shoulder_lift__shoulder,
        shoulder_pan__base,
        upperarm__elbow,
        wr1__forearm,
        wr2__wrist_1,
        wr3__wrist_2,
        wrist_1__wr2,
        wrist_2__wr3,
    )
    return ret
end

const ModelConstants = ConstPosesInit()

function solver1(mc, q, qd)

    elbow__forearm = ILKBackend.rot_z__b_x_a(q[3])
    forearm__elbow = ILKBackend.rot_z__a_x_b(q[3])
    shoulder__shoulder_pan = ILKBackend.rot_z__a_x_b(q[1])
    shoulder_lift__upperarm = ILKBackend.rot_z__b_x_a(q[2])
    upperarm__shoulder_lift = ILKBackend.rot_z__a_x_b(q[2])
    wr1__wrist_1 = ILKBackend.rot_z__b_x_a(q[4])
    wr2__wrist_2 = ILKBackend.rot_z__b_x_a(q[5])
    wr3__wrist_3 = ILKBackend.rot_z__b_x_a(q[6])
    wrist_1__wr1 = ILKBackend.rot_z__a_x_b(q[4])
    wrist_2__wr2 = ILKBackend.rot_z__a_x_b(q[5])
    wrist_3__wr3 = ILKBackend.rot_z__a_x_b(q[6])

    shoulder_lift__shoulder_pan = shoulder__shoulder_pan * mc.shoulder_lift__shoulder
    shoulder_lift__base = mc.shoulder_pan__base * shoulder_lift__shoulder_pan
    elbow__shoulder_lift = upperarm__shoulder_lift * mc.elbow__upperarm
    elbow__base = shoulder_lift__base * elbow__shoulder_lift
    wr1__elbow = forearm__elbow * mc.wr1__forearm
    wr1__base = elbow__base * wr1__elbow
    wr2__wr1 = wrist_1__wr1 * mc.wr2__wrist_1
    wr2__base = wr1__base * wr2__wr1
    wr3__wr2 = wrist_2__wr2 * mc.wr3__wrist_2
    wr3__base = wr2__base * wr3__wr2
    forearm__wrist_1 = wr1__wrist_1 * mc.forearm__wr1
    upperarm__forearm = elbow__forearm * mc.upperarm__elbow
    wrist_2__wrist_3 = wr3__wrist_3 * mc.wrist_2__wr3
    shoulder__upperarm = shoulder_lift__upperarm * mc.shoulder__shoulder_lift
    wrist_1__wrist_2 = wr2__wrist_2 * mc.wrist_1__wr2
    wrist_3__base = wr3__base * wrist_3__wr3

    v__shoulder__base = ILKBackend.matrix(6,1);
    v__shoulder__base[3] = qd[1];

    v__upperarm__base = ILKBackend.matrix(6,1);
    ILKBackend.ct_twist( shoulder__upperarm, v__shoulder__base, v__upperarm__base );
    v__upperarm__base[3] += qd[1];
    v__forearm__base = ILKBackend.matrix(6,1);
    ILKBackend.ct_twist( upperarm__forearm, v__upperarm__base, v__forearm__base );
    v__forearm__base[3] += qd[2];
    v__wrist_1__base = ILKBackend.matrix(6,1);
    ILKBackend.ct_twist( forearm__wrist_1, v__forearm__base, v__wrist_1__base );
    v__wrist_1__base[3] += qd[3];
    v__wrist_2__base = ILKBackend.matrix(6,1);
    ILKBackend.ct_twist( wrist_1__wrist_2, v__wrist_1__base, v__wrist_2__base );
    v__wrist_2__base[3] += qd[4];
    v__wrist_3__base = ILKBackend.matrix(6,1);
    ILKBackend.ct_twist( wrist_2__wrist_3, v__wrist_2__base, v__wrist_3__base );
    v__wrist_3__base[3] += qd[5];

    J_wrist_3_base = ILKBackend.matrix(6,6)
    poi_J_wrist_3_base = ILKBackend.positionView(wrist_3__base)

    ILKBackend.geometricJacobianColumn_revolute(
      poi_J_wrist_3_base,
      ILKBackend.positionView(mc.shoulder_pan__base),
      ILKBackend.zaxisView(mc.shoulder_pan__base),
      view(J_wrist_3_base, :,1) )
    ILKBackend.geometricJacobianColumn_revolute(
      poi_J_wrist_3_base,
      ILKBackend.positionView(shoulder_lift__base),
      ILKBackend.zaxisView(shoulder_lift__base),
      view(J_wrist_3_base, :,2) )
    ILKBackend.geometricJacobianColumn_revolute(
      poi_J_wrist_3_base,
      ILKBackend.positionView(elbow__base),
      ILKBackend.zaxisView(elbow__base),
      view(J_wrist_3_base, :,3) )
    ILKBackend.geometricJacobianColumn_revolute(
      poi_J_wrist_3_base,
      ILKBackend.positionView(wr1__base),
      ILKBackend.zaxisView(wr1__base),
      view(J_wrist_3_base, :,4) )
    ILKBackend.geometricJacobianColumn_revolute(
      poi_J_wrist_3_base,
      ILKBackend.positionView(wr2__base),
      ILKBackend.zaxisView(wr2__base),
      view(J_wrist_3_base, :,5) )
    ILKBackend.geometricJacobianColumn_revolute(
      poi_J_wrist_3_base,
      ILKBackend.positionView(wr3__base),
      ILKBackend.zaxisView(wr3__base),
      view(J_wrist_3_base, :,6) )

    return wrist_3__base, v__wrist_3__base, J_wrist_3_base
end

end

max_stars_repo_name: ESROCOS/ilk-compiler

max_issues_repo_name: ESROCOS/ilk-compiler

max_forks_repo_name: ESROCOS/ilk-compiler


text:
export Exemplar, SumAgg, LastValueAgg, HistogramAgg, DROP

const N_MAX_POINTS_PER_METRIC = 2_000

#####

"""
    Exemplar(;kw...)

Exemplars are example data points for aggregated data. Read [the
specification](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/metrics/sdk.md#exemplar)
to understand its relation to trace and metric.

## Keyword arguments:

  - `value`
  - `time_unix_nano`
  - `filtered_attributes`::[`StaticAttrs`](@ref), extra attributes of a
    [`Measurement`](@ref) that are not included in a [`Metric`](@ref)'s
    `:attribute_keys` field.
  - `trace_id`, the `trace_id` in the span context when the measurement happens.
  - `span_id`, the `span_id` in the span context when the measurement happens.
"""
Base.@kwdef struct Exemplar{T}
    value::T
    time_unix_nano::UInt
    filtered_attributes::StaticAttrs
    trace_id::TraceIdType
    span_id::SpanIdType
end

abstract type AbstractExemplarReservoir end

# TODO
struct SimpleFixedSizeExemplarReservoir <: AbstractExemplarReservoir end

# TODO
struct AlignedHistogramBucketExemplarReservoir <: AbstractExemplarReservoir end

#####

if VERSION >= v"1.7.0"
    include("datapoint_atomic.jl")
else
    include("datapoint_lock.jl")
end

struct AggregationStore{D<:DataPoint}
    points::Dict{StaticAttrs,D}
    unique_points::Dict{StaticAttrs,D}
    n_max_points::UInt
    n_max_attrs::UInt
    lock::ReentrantLock
end

"""
    AggregationStore{D}(;kw...) where D<:DataPoint

The `AggregationStore` holds all the aggregated datapoints in a
[`Metric`](@ref).

## Keyword arguments

  - `n_max_points = N_MAX_POINTS_PER_METRIC`, the maximum number of data points.
  - `n_max_attrs = 2 * N_MAX_POINTS_PER_METRIC`, the maximum number of allowed
    attributes. Note that each datapoint may have several attributes pointing to
    them. (Those attributes have the same key-value pair but with different
    order)
"""
function AggregationStore{D}(;
    n_max_points = N_MAX_POINTS_PER_METRIC,
    n_max_attrs = 2 * N_MAX_POINTS_PER_METRIC,
) where {D<:DataPoint}
    AggregationStore{D}(
        Dict{StaticAttrs,D}(),
        Dict{StaticAttrs,D}(),
        UInt(n_max_points),
        UInt(n_max_attrs),
        ReentrantLock(),
    )
end

Base.iterate(a::AggregationStore, args...) = iterate(a.unique_points, args...)
Base.length(m::AggregationStore) = length(m.unique_points)

function Base.getindex(m::AggregationStore, k)
    v = get(m.points, k, nothing)
    if isnothing(v)
        k_sorted = sort(k)
        get(m.points, k_sorted)
    else
        v
    end
end

function Base.get!(f, agg::AggregationStore, attrs)
    point = get(agg.points, attrs, nothing)
    if isnothing(point)
        sorted_attrs = sort(attrs)
        point = get(agg.points, sorted_attrs, nothing)
        if isnothing(point)
            if length(agg.unique_points) < agg.n_max_points
                lock(agg.lock) do
                    if haskey(agg.points, attrs)
                        agg.points[attrs]
                    elseif haskey(agg.points, sorted_attrs)
                        agg.points[sorted_attrs]
                    else
                        p = f()
                        agg.points[attrs] = p
                        agg.points[sorted_attrs] = p
                        agg.unique_points[sorted_attrs] = p
                        p
                    end
                end
            else
                @warn "maximum number of attributes reached, dropped."
                nothing
            end
        else
            if length(agg.points) < agg.n_max_attrs
                agg.points[attrs] = point
            else
                @warn "maximum cached keys in agg store reached, please consider increase `n_max_points`"
            end
            point
        end
    else
        point
    end
end

#####

abstract type AbstractAggregation end

Base.iterate(a::AbstractAggregation, args...) = iterate(a.agg_store, args...)
Base.getindex(m::AbstractAggregation, k) = getindex(m.agg_store, k)
Base.length(m::AbstractAggregation) = length(m.agg_store)

"""
    SumAgg(agg_store::AggregationStore, exemplar_reservoir_factory)

When `exemplar_reservoir_factory` set to `nothing`, no exemplar will be stored.
See more details in [the specification](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/metrics/sdk.md#sum-aggregation).
"""
struct SumAgg{T,E,F} <: AbstractAggregation
    agg_store::AggregationStore{DataPoint{T,E}}
    exemplar_reservoir_factory::F
end

"""
    SumAgg{T}()
"""
SumAgg{T}() where {T} = SumAgg(AggregationStore{DataPoint{T,Nothing}}(), () -> nothing)

function (agg::SumAgg{T,E})(e::Exemplar{<:Measurement}) where {T,E}
    point = get!(agg.agg_store, e.value.attributes) do
        DataPoint{T}(agg.exemplar_reservoir_factory())
    end
    if !isnothing(point)
        _add_to_datapoint!(point, e.value.value, e.time_unix_nano)
    end
end

#####

"""
    LastValueAgg(agg_store::AggregationStore, exemplar_reservoir_factory)

When `exemplar_reservoir_factory` set to `nothing`, no exemplar will be stored.

See more details in [the specification](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/metrics/sdk.md#last-value-aggregation).
"""
struct LastValueAgg{T,E,F} <: AbstractAggregation
    agg_store::AggregationStore{DataPoint{T,E}}
    exemplar_reservoir_factory::F
end

"""
    LastValueAgg{T}()
"""
LastValueAgg{T}() where {T} =
    LastValueAgg(AggregationStore{DataPoint{T,Nothing}}(), () -> nothing)

function (agg::LastValueAgg{T,E})(e::Exemplar{<:Measurement}) where {T,E}
    point = get!(agg.agg_store, e.value.attributes) do
        DataPoint{T}(agg.exemplar_reservoir_factory())
    end
    if !isnothing(point)
        _update_to_datapoint!(point, e.value.value, e.time_unix_nano)
    end
end

#####

const DEFAULT_HISTOGRAM_BOUNDARIES =
    (0.0, 5.0, 10.0, 25.0, 50.0, 75.0, 100.0, 250.0, 500.0, 1000.0)

struct HistogramValue{T,M,N}
    boundaries::NTuple{M,Float64}
    counts::NTuple{N,UInt}
    sum::Union{T,Nothing}
    min::Union{T,Nothing}
    max::Union{T,Nothing}
end

function HistogramValue{T}(;
    boundaries = DEFAULT_HISTOGRAM_BOUNDARIES,
    is_record_sum = true,
    is_record_min = true,
    is_record_max = true,
) where {T}
    M = length(boundaries)
    N = M + 1
    HistogramValue{T,M,N}(
        boundaries,
        ntuple(i -> UInt(0), N),
        is_record_sum ? zero(T) : nothing,
        is_record_min ? typemax(T) : nothing,
        is_record_max ? typemin(T) : nothing,
    )
end

function Base.:(+)(v::HistogramValue, x)
    n = 0
    for (i, b) in enumerate(v.boundaries)
        if x <= b
            n = i
            break
        end
    end
    if n == 0
        n = length(v.boundaries) + 1
    end

    HistogramValue(
        v.boundaries,
        (v.counts[1:n-1]..., v.counts[n] + UInt(1), v.counts[n+1:end]...),
        v.sum + x,
        isnothing(v.min) ? nothing : min(v.min, x),
        isnothing(v.max) ? nothing : max(v.max, x),
    )
end

"""
    HistogramAgg(args...)

## Arguments:

  - `boundaries::NTuple{M, Float64} where M`, the boundaries to calculate
    histogram buckets. Note that `-Inf` and `Inf` shouldn't be included.
  - `is_record_min`
  - `is_record_max`
  - `agg_store`
  - `exemplar_reservoir_factory`, when set to `nothing`, no exemplar will be stored.
"""
struct HistogramAgg{T,E,F,M,N} <: AbstractAggregation
    boundaries::NTuple{M,Float64}
    is_record_min::Bool
    is_record_max::Bool
    agg_store::AggregationStore{DataPoint{HistogramValue{T,M,N},E}}
    exemplar_reservoir_factory::F
end

"""
    HistogramAgg{T}(;boundaries = DEFAULT_HISTOGRAM_BOUNDARIES, is_record_min = true, is_record_max = true)
"""
HistogramAgg{T}(;
    boundaries = DEFAULT_HISTOGRAM_BOUNDARIES,
    is_record_min = true,
    is_record_max = true,
) where {T} = HistogramAgg(
    boundaries,
    is_record_min,
    is_record_max,
    AggregationStore{
        DataPoint{HistogramValue{T,length(boundaries),length(boundaries) + 1},Nothing},
    }(),
    () -> nothing,
)

function (agg::HistogramAgg{T,E})(e::Exemplar{<:Measurement}) where {T,E}
    point = get!(agg.agg_store, e.value.attributes) do
        t = UInt(time() * 10^9)
        DataPoint(
            HistogramValue{T}(;
                boundaries = agg.boundaries,
                is_record_min = agg.is_record_min,
                is_record_max = agg.is_record_max,
            ),
            t,
            t,
            agg.exemplar_reservoir_factory(),
        )
    end
    if !isnothing(point)
        _add_to_datapoint!(point, e.value.value, e.time_unix_nano)
    end
end

#####

struct Drop <: AbstractAggregation end

"""
All measurement will be dropped.
"""
const DROP = Drop()

#####

default_aggregation(ins::Counter{T}) where {T} = SumAgg{T}()
default_aggregation(ins::ObservableCounter{T}) where {T} = SumAgg{T}()
default_aggregation(ins::UpDownCounter{T}) where {T} = SumAgg{T}()
default_aggregation(ins::ObservableUpDownCounter{T}) where {T} = SumAgg{T}()
default_aggregation(ins::ObservableGauge{T}) where {T} = LastValueAgg{T}()
default_aggregation(ins::Histogram{T}) where {T} = HistogramAgg{T}()


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);

## original example, following here

pat = pattern_create_linear(0.0, 0.0,  0.0, 256.0);
pattern_add_color_stop_rgba(pat, 1, 0, 0, 0, 1);
pattern_add_color_stop_rgba(pat, 0, 1, 1, 1, 1);
rectangle(cr, 0, 0, 256, 256);
set_source(cr, pat);
fill(cr);
destroy(pat);

pat = pattern_create_radial(115.2, 102.4, 25.6,
                            102.4,  102.4, 128.0);
pattern_add_color_stop_rgba(pat, 0, 1, 1, 1, 1);
pattern_add_color_stop_rgba(pat, 1, 0, 0, 0, 1);
set_source(cr, pat);
arc(cr, 128.0, 128.0, 76.8, 0, 2 * pi);
fill(cr);
destroy(pat);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_gradient.png");


text:
"""
    expon_rng(λ, shape=1; seed=nothing)

Generate a `shape` element array of random variables from a Exponential(`λ`) distribution. Optionally you can set a specific seed.

# Notes

The pdf of an Exponential(λ) distribution is given as:

``f(x, λ) = λe^{-λx} \\quad x ≥ 0``

!!! note
    Some Exponential distributions are parameterized by ``β = \\frac{1}{λ}``, where λ is the number of events in an interval.
    In such cases, β represents the mean interarrival time; here we use λ to represent mean arrival rate per unit of time.

# Examples

```julia-repl
julia> expon_rng(3)
1-element Vector{Float64}:
 0.07033135663980515

julia> expon_rng(1.2, seed=42)
1-element Vector{Float64}:
 0.3296112244200808

julia> expon_rng(1.2, (2, 2))
2×2 Matrix{Float64}:
 1.9327    0.134739
 0.746861  0.155614
```

# References

D.P. Kroese, T. Taimre, Z.I. Botev. Handbook of Monte Carlo Methods.
  Wiley Series in Probability and Statistics, John Wiley & Sons, New York, 2011.
"""
function expon_rng(λ::Real, shape::Union{Int, Tuple{Vararg{Int}}}=1; seed::Union{Int, Nothing}=nothing)
    U = get_std_uniform(shape, seed=seed)
    X = zeros(shape)
    X = (-1/λ) .* log.(1 .- U)  # could also use just U
    return X
end


text:
using DataFrames

A = readdlm("u1.base",'\t';header=false)
file = "u1.base"

function pearsonCorrelation(x,y)
        @assert(length(x) == length(y))
        n = length(x)
        @assert(n>0)
        avg_x = mean(x)
        avg_y = mean(y)
        prod_diff = 0
        xdiff2 = 0
        ydiff2 = 0
        for id = 1:n
                xdiff = x[id] - avg_x
                ydiff = y[id] - avg_y
                prod_diff += xdiff*ydiff
                xdiff2 += xdiff^2
                ydiff2 += ydiff^2
        end
        return prod_diff/sqrt(xdiff2*ydiff2)
end

function PearsonCorrelationHelper(file)
     df = readtable(file,separator = '\t',header = false);
     rename!(df,:x1,:userId)
     rename!(df,:x2,:itemId)
     rename!(df,:x3,:rating)
     #rename!(df,:x4,:timestamp)
     M1 =Set(df[df[:userId] .== 1, :itemId])
     CorrelationArray = Float64[]
     (N,ind) = findmax(df[:,1])
     for(i = 2:N)
            M_i = Set(df[df[:userId] .== i, :itemId])
            if(length(intersect(M1,M_i)) == 0) push!( CorrelationArray , NA)
            else
                  X = Float64[]
                  Y = Float64[]
                  for(movie in intersect(M1,M_i))
                      push!(X , (df[(df[:userId]   .== 1)  & (df[:itemId] .== movie), :rating])[1])
                      push!(Y , (df[(df[:userId]   .== i)  & (df[:itemId] .== movie), :rating])[1])
                  end
                  push!( CorrelationArray , pearsonCorrelation(X,Y))
            end
     end
end


text:
"""
aasol(GFix!, x0, m, Vstore; maxit=20,
      rtol=1.e-10, atol=1.e-10, beta=1.0, pdata=nothing, keepsolhist = false)

C. T. Kelley, 2021

Julia code for Anderson acceleration. Nothing fancy.

Solvers fixed point problems x = G(x).

You must allocate storage for the function and fixed point map
history --> in the calling program <-- in the array Vstore.

For an n dimensional problem with Anderson(m), Vstore must have
at least 2m + 4 columns and 3m + 3 is better.  If m=0 (Picard) then
V must have at least 4 columns.

Inputs:\n

- GFix!: fixed-point map, the ! indicates that GFix! overwrites xout, your
    preallocated storage for the function value xout=G(xin).\n
    So xout=GFix!(xout,xin) or xout=GFix!(xout,xin,pdata) returns
    xout=G(xin)

- sol: Initial iterate

- m: depth for Anderson acceleration. m=0 is Picard iteration

- Vstore: Working storage array. For an n dimensional problem Vstore
  should have at least 3m+3 columns unless you are storage bound. If storage
  is a problem, then you can allocate a minimum of 2m+4 columns. The smaller
  allocation exacts a performance penalty, especially for small problems
  and small values of m. So for Anderson(3), Vstore should be no smaller
  than zeros(n,8) with zeros(n,11) a better choice. Vstore needs to
  allocate for the history of differences of the residuals and fixed
  point maps. The extra m-1 columns are for storing intermediate results
  in the downdating phase of the QR factorization for the coefficeint
  matrix of the optimization problem. See the notebook or the print book
  for the details of this mess.

  If m=0, then Vstore needs 4 columns.

Keyword Arguments (kwargs):\n

maxit: default = 20\n
limit on nonlinear iterations\n

rtol and atol: default = 1.e-10\n
relative and absolute error tolerances\n

beta:\n
Anderson mixing parameter. Changes G(x) to (1-beta)x + beta G(x).
Equivalent to accelerating damped Picard iteration. The history
vector is the one for the damped fixed point map, not the original
one. Keep this in mind when comparing results.

pdata:\n
precomputed data for the fixed point map.
Things will go better if you use this rather than hide the data
in global variables within the module for your function.

keepsolhist: default = false\n
Set this to true to get the history of the iteration in the output
tuple. This is on by default for scalar equations and off for systems.
Only turn it on if you have use for the data, which can get REALLY LARGE.

Output:\n
- A named tuple (solution, functionval, history, stats, idid,
               errcode, solhist)
where

   -- solution = converged result

   -- functionval = G(solution)
      You might want to use functionval as your solution since it's
      a Picard iteration applied to the converged Anderson result. If G
      is a contraction it will be better than the solution.

   -- history = the vector of residual norms (||x-G(x)||) for the iteration

   -- stats = named tuple (condhist, alphanorm) of the history of the
              condition numbers of the optimization problem
              and l1 norm of the coefficients.
This is only for diagosing
problems and research. Condihist[k] and alphanorm[k] are
the condition number and coefficient norm for the optimization
problem that computes iteration k+1 from iteration k.

I record this for iterations k=1, ... until the final iteration
K. So I do not record the stats for k=0 or the final iteration.
We did record the data for the final iteration in Toth/Kelley
2015 at the cost of an extra optimiztion problem solve.
Since we've already terminated, there's not any point in
collecting that data.\n
Bottom line: if history has length K+1 for iterations
0 ... K, then condhist and alphanorm have length K-1.

   -- idid=true if the iteration succeeded and false if not.

   -- errcode = 0 if if the iteration succeeded

        = -1 if the initial iterate satisfies the termination criteria

        = 10 if no convergence after maxit iterations

   -- solhist:\n
This is the entire history of the iteration if you've set
keepsolhist=true\n

solhist is an N x K array where N is the length of x and K is the number of
iterations + 1.

### Examples for aasol

#### Duplicate Table 1 from Toth-Kelley 2015.

The final entries in the condition number and coefficient norm statistics
are never used in the computation and we don't compute them in Julia.
See the docstrings, notebook, and the print book for the story on this.

```jldoctest
julia> function tothk!(G, u)
       G[1]=cos(.5*(u[1]+u[2]))
       G[2]=G[1]+ 1.e-8 * sin(u[1]*u[1])
       return G
       end
tothk! (generic function with 1 method)

julia> u0=ones(2,); m=2; vdim=3*m+3; Vstore = zeros(2, vdim);
julia> aout = aasol(tothk!, u0, m, Vstore; rtol = 1.e-10);
julia> aout.history
8-element Vector{Float64}:
 6.50111e-01
 4.48661e-01
 2.61480e-02
 7.25389e-02
 1.53107e-04
 1.18513e-05
 1.82466e-08
 1.04725e-13

julia> [aout.stats.condhist aout.stats.alphanorm]
6×2 Matrix{Float64}:
 1.00000e+00  1.00000e+00
 2.01556e+10  4.61720e+00
 1.37776e+09  2.15749e+00
 3.61348e+10  1.18377e+00
 2.54948e+11  1.00000e+00
 3.67694e+10  1.00171e+00
```

Now we put a mixing or damping paramter in there with beta = .5. This
example is nasty enough to make mixing do ok. Keep in mind
that the history is for the damped residual, not the original one.

```
julia> bout=aasol(tothk!, u0, m, Vstore; rtol = 1.e-10, beta=.5);

julia> bout.history
7-element Vector{Float64}:
 3.25055e-01
 3.70140e-02
 1.81111e-03
 9.55308e-04
 1.25936e-05
 1.40854e-09
 2.18196e-12
```

#### H-equation example with m=2. This takes more iterations than
Newton, which should surprise no one.

```jldoctest
julia> n=16; x0=ones(n,); Vstore=zeros(n,20); m=2;
julia> hdata=heqinit(x0,.99);
julia> hout=aasol(HeqFix!, x0, m, Vstore; pdata=hdata);
julia> hout.history
12-element Vector{Float64}:
 1.47613e+00
 7.47800e-01
 2.16609e-01
 4.32017e-02
 2.66867e-02
 6.82965e-03
 2.70779e-04
 6.51027e-05
 7.35581e-07
 1.85649e-09
 4.94803e-10
 5.18866e-12

```
"""
function aasol(
    GFix!,
    x0,
    m,
    Vstore;
    maxit = 20,
    rtol = 1.e-10,
    atol = 1.e-10,
    beta = 1.0,
    pdata = nothing,
    keepsolhist = false,
)
    #
    # Startup
    #
    # Set up the storage
    #
    (sol, gx, df, dg, res, DG, QP, Qd, solhist) =
           Anderson_Init(x0, Vstore, m, maxit, beta, keepsolhist)
    res = copy(sol)
    #
    #   Iteration 1
    #
    k = 0
    ~keepsolhist || (@views solhist[:, k+1] .= sol)
    gx = EvalF!(GFix!, gx, sol, pdata)
    (beta == 1.0) || (gx=betafix!(gx, sol, beta))
    res .= gx - sol
    resnorm = norm(res)
    tol = rtol * resnorm + atol
    ItData = ItStatsA(resnorm)
    toosoon = (resnorm <= tol)
#    if toosoon
#        println("aasol terminates on entry")
#    else
    if ~toosoon
    #
    #   If we need more iterations, get organized.
    #
    sol .= gx
    alpha = zeros(m + 1)
    k = k + 1
    ~keepsolhist || (@views solhist[:, k+1] .= sol)
    (gx, dg, df, res, resnorm) =
        aa_point!(gx, GFix!, sol, res, dg, df, beta, pdata)
    updateHist!(ItData, resnorm)
    end
    n=length(x0)
    RF=zeros(m,m)
    RP=zeros(m,m)
    while (k < maxit) && resnorm > tol && ~toosoon
        if m == 0
            alphanrm = 1.0
            condit = 1.0
            sol .= gx
        else
            BuildDG!(DG,m,k+1,dg)
            (QP, RP) = aa_qr_update!(QP, RP, df, m, k-1, Qd)
            mk = min(m, k)
            @views QA=QP[:,1:mk]
            @views RA=RP[1:mk,1:mk]
            theta = RA\ (QA'*res)
            condit = cond(RA)
            alphanrm = falpha(alpha, theta, min(m, k))
            copy!(sol, gx)
            @views sol .-= DG[:, 1:mk] * theta
        end
        updateStats!(ItData, condit, alphanrm)
        k += 1
        ~keepsolhist || (@views solhist[:, k+1] .= sol)
        (gx, dg, df, res, resnorm) =
            aa_point!(gx, GFix!, sol, res, dg, df, beta, pdata)
        updateHist!(ItData, resnorm)
    end
    (idid, errcode) = AndersonOK(resnorm, tol, k, toosoon)
    aaout=CloseIteration(sol, gx, ItData, idid, errcode, keepsolhist, solhist)
    return aaout
end

"""
BuildDG!(DG,m,k,dg)

Keeps the history of the fixed point map differences
"""
function BuildDG!(DG,m,k,dg)
    if m == 1
        @views copy!(DG[:, 1], dg)
    elseif k > m + 1
        for ic = 1:m-1
            @views DG[:, ic] .= DG[:, ic+1]
        end
        @views copy!(DG[:, m], dg)
    else
        @views copy!(DG[:, k-1], dg)
    end
end

"""
aa_point!(gx, gfix, sol, res, dg, df, pdata)

Evaluate the fixed point map at the new point.
Keep the books to get ready to update the coefficient matrix
for the optimization problem.
"""
function aa_point!(gx, gfix, sol, res, dg, df, beta, pdata)
    dg .= -gx
    gx = EvalF!(gfix, gx, sol, pdata)
    (beta == 1.0) || (gx=betafix!(gx, sol, beta))
    dg .+= gx
    df .= -res
    res .= gx - sol
    df .+= res
    resnorm = norm(res)
    return (gx, dg, df, res, resnorm)
end

"""
betafix(gx, sol, dg, beta)

Put the mixing parameter beta in the right place.
"""
function betafix!(gx, sol, beta)
gx=axpby!((1.0-beta),sol,beta,gx)
return gx
end


text:
"""
`module  JAC.MultiPhotonDoubleIon`
... a submodel of JAC that contains all methods for computing multi-photon two-electron (double) ionization cross sections.
"""
module MultiPhotonDoubleIon

    using ..Basics, ..ManyElectron, ..Radial

    """
    `struct  MultiPhotonDoubleIon.Settings`
        ... defines a type for the settings in estimating multi-photon two-electron (double) ionization cross sections.

        + multipoles              ::Array{EmMultipole}           ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge}              ... Specifies the gauges to be included into the computations.
        + NoPhotons               ::Int64                        ... Number of photons in the multi-photon ionization
        + photonEnergies          ::Array{Float64,1}             ... List of photon energies.
        + printBeforeComputation  ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + selectLines             ::Bool                         ... True, if lines are selected individually for the computations.
        + selectedLines           ::Array{Tuple{Int64,Int64},1}  ... List of lines, given by tupels (inital-level, final-level).

    """
    struct Settings
        multipoles                ::Array{EmMultipole}
        gauges                    ::Array{UseGauge}
        NoPhotons                 ::Int64
        photonEnergies            ::Array{Float64,1}
        printBeforeComputation    ::Bool
        selectLines               ::Bool
        selectedLines             ::Array{Tuple{Int64,Int64},1}
    end

    """
    `MultiPhotonDoubleIon.Settings()`  ... constructor for the default values of multi-photon two-electron (double) ionization estimates.
    """
    function Settings()
        Settings(EmMultipole[], UseGauge[], 0, Float64[], false, false, Tuple{Int64,Int64}[])
    end

    # `Base.show(io::IO, settings::MultiPhotonDoubleIon.Settings)`
    #		... prepares a proper printout of the variable settings::MultiPhotonDoubleIon.Settings.
    function Base.show(io::IO, settings::MultiPhotonDoubleIon.Settings)
        println(io, "multipoles:               $(settings.multipoles)  ")
        println(io, "gauges:                   $(settings.gauges)  ")
        println(io, "NoPhotons:                $(settings.NoPhotons)  ")
        println(io, "photonEnergies:           $(settings.photonEnergies)  ")
        println(io, "printBeforeComputation:   $(settings.printBeforeComputation)  ")
        println(io, "selectLines:              $(settings.selectLines)  ")
        println(io, "selectedLines:            $(settings.selectedLines)  ")
    end

    """
    `struct  MultiPhotonDoubleIon.Channel`
        ... defines a type for a multi-photon channel to help characterize multi-photon (single-electron) ionization with
            well-defined multipolarities.

        + NoPhotons      ::EmMultipole            ... Number of photons in the multi-photon ionization
        + multipoles     ::Array{EmMultipole,1}   ... Multipoles of all N incoming photons.
        + gauge          ::EmGauge                ... Gauge for dealing with the (coupled) radiation field.
        + kappa1         ::Int64                  ... partial-wave of the first free electron
        + kappa2         ::Int64                  ... partial-wave of the second free electron
        + phase1         ::Float64                ... phase of the first partial wave
        + phase2         ::Float64                ... phase of the second partial wave
        + symmetry1      ::LevelSymmetry          ... total angular momentum and parity of scattering state after coupling of first electron
        + symmetry2      ::LevelSymmetry          ... total angular momentum and parity of scattering state after coupling of second electron
        + amplitude      ::Complex{Float64}       ... multi-photon ionization amplitude associated with the given channel.
   """
    struct  Channel
        NoPhotons        ::EmMultipole
        multipoles       ::Array{EmMultipole,1}
        gauge            ::EmGauge
        kappa1           ::Int64
        kappa2           ::Int64
        phase1           ::Float64
        phase2           ::Float64
        symmetry1        ::LevelSymmetry
        symmetry2        ::LevelSymmetry
        amplitude        ::Complex{Float64}
    end

    # MultiPhotonDoubleIon line between initial and final (bound-state) levels
    """
    `struct  MultiPhotonDoubleIon.Line`  ... defines a type for a multi-photon ionization line that may include the definition of channels.

        + initialLevel     ::Level           ... initial-(state) level
        + finalLevel       ::Level           ... final-(state) level
        + NoPhotons        ::Int64           ... Number of photons in the multi-photon process
        + photonEnergy     ::Float64         ... Energy of the incoming photons; all photons are assumed to have equal energy.
        + crossSection     ::EmProperty      ... Cross section for this multi-photon ionization.
        + hasChannels      ::Bool            ... Determines whether the individual (sub-) channels are defined in terms of
                                                 their multipolarities, etc., or not.
        + channels         ::Array{MultiPhotonDoubleIon.Channel,1}  ... List of MultiPhotonDoubleIon.Channels of this line.
    """
    struct  Line
        initialLevel       ::Level
        finalLevel         ::Level
        NoPhotons          ::Int64
        photonEnergy       ::Float64
        crossSection       ::EmProperty
        hasChannels        ::Bool
        channels           ::Array{MultiPhotonDoubleIon.Channel,1}
    end

    # `Base.show(io::IO, line::MultiPhotonDoubleIon.Line)`  ... prepares a proper printout of the variable line::MultiPhotonDoubleIon.Line.
    function Base.show(io::IO, line::MultiPhotonDoubleIon.Line)
        println(io, "initialLevel:      $(line.initialLevel)  ")
        println(io, "finalLevel:        $(line.finalLevel)  ")
        println(io, "NoPhotons:         $(line.NoPhotons)  ")
        println(io, "photonEnergy:      $(line.photonEnergy)  ")
        println(io, "crossSection:      $(line.crossSection)  ")
        println(io, "hasChannels:       $(line.hasChannels)  ")
        println(io, "channels:          $(line.channels)  ")
    end

    """
    `MultiPhotonDoubleIon.computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                                       settings::MultiPhotonDoubleIon.Settings; output=true)`
        ... to compute the multiphoton transition amplitudes and all properties as requested by the given settings.
            A list of lines::Array{MultiPhotonDoubleIon.Lines} is returned.
    """
    function  computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                           settings::MultiPhotonDoubleIon.Settings; output=true)
        println("")
        printstyled("MultiPhotonDoubleIon.computeLines(): The computation of multi-photon double ionization amplitudes starts now ... \n", color=:light_green)
        printstyled("----------------------------------------------------- ---------------------------------------------------------- \n", color=:light_green)
        println("")
        #
        #
        pathways = "Not yet implemented !"
        #
        if    output    return( pathways )
        else            return( nothing )
        end
    end

end # module


text:
"""
check https://github.com/thautwarm/MLStyle.jl/blob/master/matrix_benchmark.jl
"""


text:
"""
Load segment size list
"""
function load_seg_list(fn)
    return readdlm(fn, ';', Int64)
end

function load_semmap(fname)
    table = readdlm(fname, ';', Int)
    Dict( table[i,1] => table[i,2] for i in 1:size(table,1) )
end

"""
Load Synaptor edge csv & parse columns appropriately
"""
function load_edges(fn)
    tbl = readdlm(fn, ';')[:,1:4]
    tbl[:,2] = [map(parse, split(tbl[i,2][2:end-1],",")) for i in 1:size(tbl,1)]
    tbl[:,3] = [map(parse, split(tbl[i,3][2:end-1],",")) for i in 1:size(tbl,1)]
    return tbl
end

"""
Load Omni working/valid/uncertain
"""
function load_valid_list(fn)
	return readdlm(fn, ',', Int64)
end

function read_h5( fname, read_whole_dset=true, h5_dset_name="/main" )

  if read_whole_dset
    d = h5read( fname, h5_dset_name );
  else
    f = h5open( fname );
    d = f[h5_dset_name];
  end

  d
end

function write_h5( dset, fname, dset_name="/main" )
  if isfile(fname) rm(fname) end
  HDF5.h5write(fname, dset_name, dset)
end

function write_map_file( output_fname, dicts... )

  open(output_fname, "w+") do f
    if length(dicts) > 0
    for k in keys(dicts[1])

      vals = ["$(d[k]);" for d in dicts ];
      write(f, "$k;$(vals...)\n")

    end #for
    end #if
  end#open

end

"""
Load sparse matrix from file with format i, j, v
"""
function read_sparse(fn)
    d = readdlm(fn)
    n = Int64(maximum(d[:,1]))
    m = Int64(maximum(d[:,2]))
    S = spzeros(n,m)
    for i in 1:size(d,1)
        S[Int64(d[i,1]), Int64(d[i,2])] = d[i,3]
    end
    return S
end

"""
Write sparse matrix to file with format i, j, v
"""
function write_sparse(fn, arr)
    r = []
    c = []
    v = []
    rows = rowvals(arr)
    vals = nonzeros(arr)
    n = size(arr,1)
    for j in 1:n
        for i in nzrange(arr, j)
            push!(r, rows[i])
            push!(c, j)
            push!(v, vals[i])
        end
    end
    if length(nzrange(arr, n)) == 0
        push!(r, n)
        push!(c, n)
        push!(v, 0)
    end
    if !(n in rows)
        push!(r, n)
        push!(c, n)
        push!(v, 0)
    end

    writedlm(fn, hcat(r,c,v))
end


text:
# indexing
export setindex, addindex

@generated function setindex(x::NamedTuple,y,v::Val)
    k = first(v.parameters)
    k ∉ x.names ? :x : :( (x..., $k=y) )
end

@generated function addindex(x::NamedTuple,y,v::Val)
    k = first(v.parameters)
    :( (x..., $k=y) )
end

# git
function get_hash(dir::AbstractString)
    read(`git -C $dir rev-parse --short HEAD`, String) |> strip
end

# Conversions
mmtoin(x::Real) = x/25.399986284
intomm(x::Real) = x*25.399986284

# Rounding
roundup(x,n)   = round(x,digits=n) >= x ? round(x,digits=n) : round(x+1/(10^n),digits=n)
rounddown(x,n) = round(x,digits=n) <= x ? round(x,digits=n) : round(x+1/(10^n),digits=n)

# DSP
using DSP

function twodresample(mat, rate, fs=1250)
    out = zeros(floor(Int, size(mat,1)*rate), size(mat,2))
    for i in 1:size(mat,2)
        out[:,i] = resample(mat[:,i],rate)
    end
    floor(Int, fs*rate), out
end

# Wavelets
using Wavelets
function wavedec2d(clips::Array{T,2}, wavtype = WT.haar) where T<:AbstractFloat
    # level of wavelet transform
    wt = wavelet(wavtype)
    N = floor(Int,log2(size(clips,2)))
    # intialise array
    components_all = zeros(clips)
    for i in 1:size(components_all,1)
        components_all[i,:] = dwt(clips[i,:], wt, N)
    end
    return components_all::Array{T,2}
end

# Math
from_choose(n,k) = factorial(n)/(factorial(k)*factorial(n-k))

# NaN Math
using NaNMath; const nm = NaNMath
nanmaximum(x) = nm.maximum(x)
nanminimum(x) = nm.minimum(x)

nan2zero!(arr) = arr[isnan.(arr)] .= zero(eltype(arr))

nans(x) = (z = similar(x); fill!(z,NaN))

nanmap(f,x,args...;kwargs...) = f(filter(!isnan,x),args...,kwargs...)

function nanmapslices(f,x,args...;dims=1, kwargs...)
    function f2(x_slice)
        x_slice_nonan = filter(!isnan, x_slice)
        f(x_slice_nonan, args..., kwargs...)
    end
    mapslices(f2, x, dims=dims)
end


text:
#Heatmap plots

function PlotHeatmap(x::AbstractArray, rows, cols, scale_min = 0.0, scale_max = 1.0;
                 label = "", label_fmt="%.1f", bounds_min = LibCImPlot.ImPlotPoint(0.0,0.0),
                 bounds_max = LibCImPlot.ImPlotPoint(1.0,1.0))

    if eltype(x) == Float64
        LibCImPlot.PlotHeatmapdoublePtr(label, x, Cint(rows), Cint(cols), Cdouble(scale_min),
                                    Cdouble(scale_max), label_fmt, bounds_min, bounds_max)
    elseif eltype(x) == Float32
        LibCImPlot.PlotHeatmapFloatPtr(label, x, Cint(rows), Cint(cols), Cfloat(scale_min),
                                   Cfloat(scale_max), label_fmt, bounds_min, bounds_max)
    else
        throw("Only 32/64 bit floating point values supported")
    end
end


text:
#
# Define tests
#

function fill_ghost_cells!(c, κtop, κbottom, model, fieldbcs)
    fill_bottom_ghost_cell!(fieldbcs.bottom, c, κbottom, model)
    fill_top_ghost_cell!(fieldbcs.top, c, κtop, model, c.grid.N)
    return nothing
end

function test_cell_field_construction(T, N, L)
    grid = UniformGrid(T, N, L)
    c = CellField(grid)
    typeof(c) <: CellField
end

function test_face_field_construction(T, N, L)
    grid = UniformGrid(T, N, L)
    f = FaceField(grid)
    typeof(f) <: FaceField
end

function test_field_indexing()
    grid = UniformGrid(3, 4.2)
    c = CellField(grid)
    f = FaceField(grid)
    val = 2.1
    c[2] = val
    f[2] = val
    c[2] == val && f[2] == val
end

function test_cell_∂z(T)
    grid = UniformGrid(T, 2, 2.0)
    c = CellField([2, 4], grid)
    cz = ∂z(c)
    cz_answer = FaceField([0, 2, 0], grid)
    cz.data == cz_answer.data
end

function test_face_∂z(T)
    grid = UniformGrid(T, 2, 2.0)
    f = FaceField([2, 4, 6], grid)
    fz = ∂z(f)
    fz_answer = CellField([2, 2], grid)
    fz.data[1:2] == fz_answer.data[1:2]
end

function test_cell_plus(T)
    grid = UniformGrid(T, 4, 2.0)
    c1 = CellField([0, 1, 2, 0], grid)
    c2 = CellField([0, 3, 4, 0], grid)
    c3 = CellField([0, 4, 6, 0], grid)
    c1_plus_c2 = c1 + c2
    c1_plus_c2.data == c3.data
end

function test_cell_times(T)
    grid = UniformGrid(T, 4, 2.0)
    c1 = CellField([0, 1, 2, 0], grid)
    c2 = CellField([0, 3, 4, 0], grid)
    c3 = CellField([0, 3, 8, 0], grid)
    c1_times_c2 = c1 * c2
    c1_times_c2.data == c3.data
end

function test_face_plus(T)
    grid = UniformGrid(T, 2, 2.0)
    c1 = FaceField([1, 2, 2], grid)
    c2 = FaceField([3, 4, 3], grid)
    c3 = FaceField([4, 6, 5], grid)
    c1_plus_c2 = c1 + c2
    c1_plus_c2.data == c3.data
end

function test_face_times(T)
    grid = UniformGrid(T, 2, 2.0)
    f1 = FaceField([1, 2, 2], grid)
    f2 = FaceField([3, 4, 3], grid)
    f3 = FaceField([3, 8, 6], grid)
    f1_times_f2 = f1 * f2
    f1_times_f2.data == f3.data
end

function test_set_scalar_field(loc, T)
    g = UniformGrid(T, 2, 2.0)
    f = Field(loc, g)
    a = 7
    set!(f, a)
    !any(@. !(f.data == a))
end

function test_set_array_field(loc, T)
    g = UniformGrid(T, 2, 2.0)
    f = Field(loc, g)
    data = rand(1:10, length(f))
    set!(f, data)
    f.data[1:length(f)] == data
end

function test_set_function_field(loc, T)
    g = UniformGrid(T, 2, 2.0)
    f = Field(loc, g)
    fcn(z) = z^2
    z = nodes(f)
    set!(f, fcn)
    i = 2
    f[i] == fcn(z[i])
end

function test_integral(T)
    grid = UniformGrid(T, 3, 3.0)
    c = CellField([1, 2, 3], grid)
    integral(c) / grid.L == 2
end

function test_integral_range(T, N, L, z₋, z₊)
    grid = UniformGrid(T, N, L)
    c = CellField(ones(N), grid)
    integral(c, z₋, z₊) == z₊ - z₋
end

struct DummyModel end

function test_ghost_cell_value(T)
    model = DummyModel()
    grid = UniformGrid(T, 3, 3.0)
    c = CellField([1, 2, 3], grid)
    c_top = T(0.2)
    c_bottom = T(0.11)
    top_bc = ValueBoundaryCondition(c_top)
    bottom_bc = ValueBoundaryCondition(c_bottom)
    fieldbcs = FieldBoundaryConditions(bottom_bc, top_bc)
    fill_ghost_cells!(c, 0, 0, model, fieldbcs)

    onface(c, 1) ≈ c_bottom && onface(c, grid.N+1) ≈ c_top
end

function test_ghost_cell_gradient(T)
    model = DummyModel()
    grid = UniformGrid(T, 3, 3.0)
    c = CellField([1, 2, 3], grid)
    cz_top = T(0.2)
    cz_bottom = T(0.11)
    top_bc = GradientBoundaryCondition(cz_top)
    bottom_bc = GradientBoundaryCondition(cz_bottom)
    fieldbcs = FieldBoundaryConditions(bottom_bc, top_bc)

    fill_ghost_cells!(c, 0, 0, model, fieldbcs)

    ∂z(c, 1) ≈ cz_bottom && ∂z(c, grid.N+1) ≈ cz_top
end

function test_ghost_cell_flux(T)
    model = DummyModel()
    grid = UniformGrid(T, 3, 3.0)
    c = CellField([1, 2, 3], grid)
    Fc_top = T(0.2)
    Fc_bottom = T(0.11)
    top_bc = FluxBoundaryCondition(Fc_top)
    bottom_bc = FluxBoundaryCondition(Fc_bottom)
    fieldbcs = FieldBoundaryConditions(bottom_bc, top_bc)

    κ = T(0.47)
    fill_ghost_cells!(c, κ, κ, model, fieldbcs)

    (OceanTurb.flux(0, κ, c, 1) ≈ Fc_bottom
        && OceanTurb.flux(0, κ, c, grid.N+1) ≈ Fc_top)
end

function test_absolute_error(T)
    grid = UniformGrid(T, 3, 3)
    c = CellField([0, 0, 0], grid)
    d = CellField([2, 2, 2], grid)
    absolute_error(c, d) == 2
end

function test_relative_error(T)
    grid = UniformGrid(T, 3, 1)
    c = CellField([0, 0, 0], grid)
    d = CellField([2, 2, 2], grid)
    relative_error(c, d) == 1
end

@testset "Fields" begin
    nz, Lz = 3, 4.2
    for T in (Float64,)
        @test test_cell_field_construction(T, nz, Lz)
        @test test_face_field_construction(T, nz, Lz)
        @test test_cell_∂z(T)
        @test test_face_∂z(T)
        @test test_cell_plus(T)
        @test test_cell_times(T)
        for loc in (Face, Cell)
            @test test_set_scalar_field(loc, T)
            @test test_set_array_field(loc, T)
            @test test_set_function_field(loc, T)
        end
        @test test_integral(T)
        @test test_integral_range(T, 10, 10, -5, -3)
        @test test_integral_range(T, 10, 10, -5.1, -4.5)
        @test test_integral_range(T, 17, 3, -2.5, -1.3)
        @test test_integral_range(T, 3, 30, -2.5, -1.3)
        @test test_ghost_cell_value(T)
        @test test_ghost_cell_gradient(T)
        @test test_ghost_cell_flux(T)
        @test test_absolute_error(T)
        @test test_relative_error(T)
    end
    @test test_field_indexing()
end


text:
code = [1,12,2,3,1,1,2,3,1,3,4,3,1,5,0,3,2,9,1,19,1,5,19,23,2,9,23,27,1,27,5,31,2,31,13,35,1,35,9,39,1,39,10,43,2,43,9,47,1,47,5,51,2,13,51,55,1,9,55,59,1,5,59,63,2,6,63,67,1,5,67,71,1,6,71,75,2,9,75,79,1,79,13,83,1,83,13,87,1,87,5,91,1,6,91,95,2,95,13,99,2,13,99,103,1,5,103,107,1,107,10,111,1,111,13,115,1,10,115,119,1,9,119,123,2,6,123,127,1,5,127,131,2,6,131,135,1,135,2,139,1,139,9,0,99,2,14,0,0]

function perform(opcode, input1, input2)
    if opcode == 1
        return input1 + input2
    elseif opcode == 2
        return input1 * input2
    end
end

for i in 1:4:length(code)
    if code[i] != 99
        code[code[i + 3] + 1] = perform(code[i], code[code[i + 1] + 1], code[code[i + 2] + 1])
    else
        println(code[1])
        exit()
    end
end

max_stars_repo_path: Day 2/Day 2 Part 1.jl

max_issues_repo_path: Day 2/Day 2 Part 1.jl

max_forks_repo_path: Day 2/Day 2 Part 1.jl


text:
using DataStructures
using VectorLogging

export TraceRecord, analyze_traces, show_trace

mutable struct TraceRecord
    enter
    exit
    parent::Union{Nothing, TraceRecord}
    children # ::Tuple{TraceRecord}
end

istrace(le::VectorLogging.LogEntry) =
    le._module == NahaJuliaLib &&
    le.group == :trace

istraceenter(le::VectorLogging.LogEntry) =
    le.message == TRACE_ENTER

istraceexit(le::VectorLogging.LogEntry) =
    le.message == TRACE_EXIT

function sametrace(le1::VectorLogging.LogEntry,
                   le2::VectorLogging.LogEntry)
    le1.keys[:id] == le2.keys[:id]
end

tracethread(le::VectorLogging.LogEntry) = le.keys[:threadid]

struct StackEntry
    enter
    children

    StackEntry(enter) =
        new(enter, Vector{TraceRecord}())
end

"""
    analyze_traces(log::VectorLogger)
Given a log containinig log records produced by `@trace`, return a
vector of `TraceRecord`s.  Each of those `TraceRecord`s is the root of
a tree of `TraceRecord`s that represent the call tree.

Use `show_trace` to print the call hierarchy in a human readable form.
"""
function analyze_traces(log::VectorLogger)::Vector{TraceRecord}
    result = Vector{TraceRecord}()
    # stacks is indexed by the threadid of a log entry:
    stacks = Dict{Int, Stack{StackEntry}}()
    function getstack(threadid)
        if haskey(stacks, threadid)
            stacks[threadid]
        else
            stack = Stack{StackEntry}()
            stacks[threadid] = stack
            stack
        end
    end
    for le in log.log
        if !istrace(le)
            continue
        end
        stack = getstack(tracethread(le))
        if istraceenter(le)
            push!(stack, StackEntry(le))
        elseif istraceexit(le)
            # If a traced function is thrown from intentionally or due
            # to an error, that function's exit won't be logged.
            #
            # Trace IDs are monotonically increasing, so the closer a
            # frame is to the top of the stack, the greater its ID
            # when compared to lower frames in the same dynamic
            # context.
            #
            # This means that when we see an "exit" entry, we must
            # construct a TraceRecord for each "enter" entry on the
            # stack until and including the one with the same traceid
            # as the "exit".
            while true
                top = pop!(stack)
                match = le.keys[:id] == top.enter.keys[:id]
                new_tr = TraceRecord(top.enter,
                                     match ? le : :nothing,
                                     nothing,
                                     Tuple(top.children))
                for child in new_tr.children
                    child.parent = new_tr
                end
                if isempty(stack)
                    push!(result, new_tr)
                else
                    push!(first(stack).children, new_tr)
                end
                if match
                    break
                end
            end
        end
    end
    result
end

"""
    show_trace(trace::TraceRecord)
Print the specified `TraceRecord` hierarchy in a human readable form.
"""
show_trace(trace::TraceRecord) = show_trace(Base.stdout, trace)

function show_trace(io::IO, trace::TraceRecord)
    function st(trace, level)
        println(io, repeat("  ", level),
                trace.enter.keys[:fcall],
                if trace.exit == nothing
                    ""
                else
                    string(" => ", trace.exit.keys[:result])
                end
                )
        for c in trace.children
            st(c, level + 1)
        end
    end
    st(trace, 0)
end


text:
const TEST_NAME = "Xorshift1024Star"

include("common.jl")

using RNG.Xorshifts

r = Xorshift1024Star(123)

test_all(r, 100_000_000)


text:
using BinDeps
using Compat

@BinDeps.setup

zmq = library_dependency("zmq", aliases = ["libzmq"])

provides(Sources,URI("https://archive.org/download/zeromq_3.2.4/zeromq-3.2.4.tar.gz"),zmq)
provides(BuildProcess,Autotools(libtarget = "src/.libs/libzmq."*BinDeps.shlib_ext),zmq)

if is_windows()
    using WinRPM
    provides(WinRPM.RPM, "zeromq", [zmq], os = :Windows)
end

if is_apple()
    if Pkg.installed("Homebrew") === nothing
        error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")  end
    using Homebrew
    provides( Homebrew.HB, "zeromq32", zmq, os = :Darwin )
end

@BinDeps.install @compat Dict(:zmq => :zmq)


text:
"""
	struct ExtractKeyAsField{S,V} <: AbstractExtractor
		key::S
		item::V
	end

extracts all items in `vec` and in `other` and return them as a ProductNode.
"""
struct ExtractKeyAsField{S,V} <: AbstractExtractor
	key::S
	item::V
end

function (e::ExtractKeyAsField)(v::V) where {V<:Union{Missing,Nothing}}
	BagNode(ProductNode((key = e.key(missing), item = e.item(missing)))[1:0], [0:-1])
end

function (e::ExtractKeyAsField)(v::ExtractEmpty)
	BagNode(ProductNode((key = e.key(v), item = e.item(v))), Mill.AlignedBags(Array{UnitRange{Int64},1}()))
end

function (e::ExtractKeyAsField)(vs::Dict)
	isempty(vs) && return(e(missing))
	items = map(collect(vs)) do (k,v)
		ProductNode((key = e.key(k), item = e.item(v)))
	end
	BagNode(reduce(catobs, items), [1:length(vs)])
end

Base.hash(e::ExtractKeyAsField, h::UInt) = hash((e.key, e.item), h)
Base.:(==)(e1::ExtractKeyAsField, e2::ExtractKeyAsField) = e1.key == e2.key && e1.item == e2.item