text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

"""
    Plans

Module part of StateSpaceEcon. This module implements the `Plan` data structure,
which is used in simulations. The plan object contains information about the
range of the simulation and which variables and shocks are exogenous or
endogenous at each period of the range.

### Constructors
  * [`Plan(model, range)`](@ref Plan)

### Modify the plan
  * [`exogenize!`](@ref), [`endogenize!`](@ref) - make variables exogenous or
    endogenous
  * [`exog_endo!`](@ref), [`endo_exog!`](@ref) - swap exogenous and endogenous
    variables
  * [`autoexogenize!`](@ref) - exogenize and endogenize variables according to
    the list in the model

### Prepare data for simulation
  * [`zeroarray`](@ref), [`zerodict`](@ref), [`zerodata`](@ref) - prepare a
    matrix or a dictionary or a [`SimData`](@ref) of data for the simulation
    containing zeros.
  * [`steadystatearray`](@ref), [`steadystatedict`](@ref),
    [`steadystatedata`](@ref) - prepare a matrix or a dictionary or a [`SimData`](@ref) of data for the
    simulation containing the steady state.

"""
module Plans

using TimeSeriesEcon
using ModelBaseEcon

# exported from this module that are also exported from StateSpaceEcon
export Plan,
    exogenize!, endogenize!, exog_endo!, endo_exog!, autoexogenize!,
    zeroarray, zerodict, steadystatearray, steadystatedict, zerodata, steadystatedata

# exported from this module, but meant for internal use in StateSpaceEcon
export plansum

"""
    Plan{T <: MIT}

A data structure representing the simulation plan. It holds information about
the time range of the simulation and which variables/shocks are exogenous at
each period.

"""
struct Plan{T <: MIT} <: AbstractVector{Vector{Symbol}}
    range::AbstractUnitRange{T}
    varshks::NamedTuple
    exogenous::BitArray{2}
end

Plan(model::Model, r::Union{MIT,Int}) = Plan(model, r:r)

"""
    Plan(model, range)

Create a default simulation plan for the given model over the given range. The
range of the plan is augmented to include periods before and after the given
range, over which initial and final conditions will be applied.

Instead of a range, one could also pass in a single moment in time
([`MIT`](@ref TimeSeriesEcon.MIT)) instance, in which case it is interpreted as
a range of length 1.

"""
function Plan(model::Model, range::AbstractUnitRange)
    if !(eltype(range) <: MIT)
        range = UnitRange{MIT{Unit}}(range)
    end
    range = (first(range) - model.maxlag):(last(range) + model.maxlead)
    local varshks = model.varshks
    local N = length(varshks)
    local names = tuple(Symbol.(varshks)...) # force conversion of Vector{ModelSymbol} to NTuple{N,Symbol}
    p = Plan(range, NamedTuple{names}(1:N), falses(length(range), length(varshks)))
    for (ind, var) = enumerate(varshks)
        if isshock(var) || isexog(var)
            p.exogenous[:, ind] .= [true]
        end
    end
    return p
end

#######################################
# AbstractVector interface

Base.size(p::Plan) = size(p.range)
Base.axes(p::Plan) = (p.range,)
Base.length(p::Plan) = length(p.range)
Base.IndexStyle(::Plan) = IndexLinear()
Base.similar(p::Plan) = Plan(p.range, p.varshks, similar(p.exogenous))
Base.copy(p::Plan) = Plan(p.range, p.varshks, copy(p.exogenous))

@inline _offset(p::Plan{T}, idx::T) where T = idx - first(p.range) + 1
@inline _offset(p::Plan{T}, idx::AbstractUnitRange{T}) where T = idx .- first(p.range) .+ 1

## Integer index is taken as is
Base.getindex(p::Plan, idx::Int) = p[p.range[idx]]
Base.getindex(p::Plan, idx::AbstractUnitRange{Int}) = p[p.range[idx]]
Base.getindex(p::Plan, idx::Int, ::Val{:inds}) = findall(p.exogenous[idx,:])

# Index of the frequency type returns the list of exogenous symbols
Base.getindex(p::Plan{T}, idx::T) where T <: MIT = [keys(p.varshks)[p[_offset(p, idx),Val(:inds)]]...,]
function Base.getindex(p::Plan{T}, idx::T, ::Val{:inds}) where T <: MIT
    first(p.range) ≤ idx ≤ last(p.range) || throw(BoundsError(p, idx))
    return p[_offset(p, idx), Val(true)]
end

# A range returns the plan trimmed over that exact range.
Base.getindex(p::Plan{MIT{Unit}}, rng::AbstractUnitRange{Int}) = p[UnitRange{MIT{Unit}}(rng)]
function Base.getindex(p::Plan{T}, rng::AbstractUnitRange{T}) where T <: MIT
    rng.start < p.range.start && throw(BoundsError(p, rng.start))
    rng.stop > p.range.stop && throw(BoundsError(p, rng.stop))
    return Plan{T}(rng, p.varshks, p.exogenous[_offset(p, rng), :])
end

# A range with a model returns a plan trimmed over that range and extended for initial and final conditions.
Base.getindex(p::Plan{MIT{Unit}}, rng::AbstractUnitRange{Int}, m::Model) = p[UnitRange{MIT{Unit}}(rng), m]
function Base.getindex(p::Plan{T}, rng::AbstractUnitRange{T}, m::Model) where T <: MIT
    rng = (rng.start - m.maxlag):(rng.stop + m.maxlead)
    return p[rng]
end

Base.setindex!(p::Plan, x, i...) = error("Cannot assign directly. Use `exogenize` and `endogenize` to alter plan.")

#######################################
# Pretty printing

Base.summary(io::IO, p::Plan) = print(io, typeof(p), " with range ", p.range)

# Used in the show() implementation below
function collapsed_range(p::Plan{T}) where T <: MIT
    ret = Pair{Union{T,UnitRange{T}},Vector{Symbol}}[]
    i1 = first(p.range)
    i2 = i1
    val = p[i1]
    make_key() = i1 == i2 ? i1 : i1:i2
    for i in p.range[2:end]
        val_i = p[i]
        if val_i == val
            i2 = i
        else
            push!(ret, make_key() => val)
            i1 = i2 = i
            val = val_i
        end
    end
    push!(ret, make_key() => val)
end

Base.show(io::IO, ::MIME"text/plain", p::Plan) = show(io, p)
function Base.show(io::IO, p::Plan)
    # 0) show summary before setting :compact
    summary(io, p)
    isempty(p) && return
    # print(io, ":")
    nrow, ncol = displaysize(io)
    limit = get(io, :limit, true)
    cp = collapsed_range(p)
    # find the longest string left of "=>" for padding
    maxl = maximum(length("$k") for (k, v) in cp)
    if limit
        dcol = ncol - maxl - 6
    else
        dcol = typemax(Int)
    end
    function print_exog(names)
        if isempty(names)
            print(io, "∅")
            return
        end
        lens = cumsum(map(x -> length("$x, "), names))
        show = lens .< dcol
        show[1] = true
        print(io, join(names[show], ", "))
        if !all(show)
            print(io, ", …")
        end
    end
    if !limit || length(cp) <= nrow - 5
        for (r, v) in cp
            print(io, "\n  ", lpad("$r", maxl, " "), " → ")
            print_exog(v)
        end
    else
        top = div(nrow - 5, 2)
        bot = length(cp) - nrow + 6 + top
        for (r, v) in cp[1:top]
            print(io, "\n  ", lpad("$r", maxl, " "), " → ")
            print_exog(v)
        end
        print(io, "\n   ⋮")
        for (r, v) in cp[bot:end]
            print(io, "\n  ", lpad("$r", maxl, " "), " → ")
            print_exog(v)
        end
    end
end

#######################################
# The user interface to modify plan instances.

"""
    setplanvalue!(plan, value, vars, range)
    setplanvalue!(plan, value, vars, date)
    setplanvalue!(plan, value, vars, dates)

Modify the status of the given variable(s) on the given date(s). If `value` is
`true` then variables become exogenous, otherwise they become endogenous.

"""
function setplanvalue!(p::Plan{T}, val::Bool, vars::Array{Symbol,1}, date::AbstractUnitRange{T}) where T <: MIT
    firstindex(p) ≤ first(date) && last(date) ≤ lastindex(p) || throw(BoundsError(p, date))
    idx1 = date .- firstindex(p) .+ 1
    for v in vars
        idx2 = get(p.varshks, v, 0)
        if idx2 == 0
            throw(ArgumentError("Unknown variable or shock name: $v"))
        end
        p.exogenous[idx1, idx2] .= val
    end
    return p
end
setplanvalue!(p::Plan{MIT{Unit}}, val::Bool, vars::AbstractArray{Symbol,1}, date::AbstractUnitRange{Int}) = setplanvalue!(p, val, vars, UnitRange{MIT{Unit}}(date))
setplanvalue!(p::Plan{MIT{Unit}}, val::Bool, vars::AbstractArray{Symbol,1}, date::Integer) = setplanvalue!(p, val, vars, ii(date):ii(date))
setplanvalue!(p::Plan{MIT{Unit}}, val::Bool, vars::AbstractArray{Symbol,1}, date::MIT{Unit}) = setplanvalue!(p, val, vars, date:date)
setplanvalue!(p::Plan{T}, val::Bool, vars::AbstractArray{Symbol,1}, date::T) where T <: MIT = setplanvalue!(p, val, vars, date:date)
setplanvalue!(p::Plan, val::Bool, vars::AbstractArray{Symbol,1}, date) = (foreach(d -> setplanvalue!(p, val, vars, d), date); p)

"""
    exogenize!(plan, vars, date)

Modify the given plan so that the given variables will be exogenous on the given
dates. `vars` can be a `Symbol` or a `String` or a `Vector` of such. `date` can
be a moment in time (same type as the plan), or a range or an iterable or a
container.

"""
exogenize!(p::Plan, var, date) = setplanvalue!(p, true, Symbol[var,], date)
exogenize!(p::Plan, vars::AbstractVector, date) = setplanvalue!(p, true, Symbol[vars...], date)
exogenize!(p::Plan, vars::Tuple, date) = setplanvalue!(p, true, Symbol[vars...], date)

"""
    endogenize!(plan, vars, date)

Modify the given plan so that the given variables will be endogenous on the
given dates. `vars` can be a `Symbol` or a `String` or a `Vector` of such.
`date` can be a moment in time (same type as the plan), or a range or an
iterable or a container.

"""
endogenize!(p::Plan, var, date) = setplanvalue!(p, false, Symbol[var,], date)
endogenize!(p::Plan, vars::AbstractVector, date) = setplanvalue!(p, false, Symbol[vars...], date)
endogenize!(p::Plan, vars::Tuple, date) = setplanvalue!(p, false, Symbol[vars...], date)

"""
    exog_endo!(plan, exog_vars, endo_vars, date)

Modify the given plan so that the given variables listed in `exog_vars` will be
exogenous and the variables listed in `endo_vars` will be endogenous on the
given dates. `exog_vars` and `endo_vars` can each be a `Symbol` or a `String` or
a `Vector` of such. `date` can be a moment in time (same type as the plan), or a
range or an iterable or a container.

"""
function exog_endo!(p::Plan, exog, endo, date)
    setplanvalue!(p, true, Symbol[exog...], date)
    setplanvalue!(p, false, Symbol[endo...], date)
end

"""
    endo_exog!(plan, endo_vars, exog_vars, date)

Modify the given plan so that the given variables listed in `exog_vars` will be
exogenous and the variables listed in `endo_vars` will be endogenous on the
given dates. `exog_vars` and `endo_vars` can each be a `Symbol` or a `String` or
a `Vector` of such. `date` can be a moment in time (same type as the plan), or a
range or an iterable or a container.

"""
@inline endo_exog!(p::Plan, endo, exog, date) = exog_endo!(p, exog, endo, date)

"""
autoexogenize!(plan, model, date)

Modify the given plan according to the "autoexogenize" protocol defined in the
given model. All variables in the autoexogenization list become endogenous and
their corresponding shocks become exogenous over the given date or range. `date`
can be a moment in time (same frequency as the given plan), a range, an
iterable, or a container.

"""
function autoexogenize!(p::Plan, m::Model, date)
    auto_vars = keys(m.autoexogenize)
    auto_shks = values(m.autoexogenize)
    exog_endo!(p, auto_vars, auto_shks, date)
end

#######################################
# The user interface to prepare data for simulation.

@inline TimeSeriesEcon.frequencyof(p::Plan) = frequencyof(p.range)
@inline TimeSeriesEcon.firstdate(p::Plan) = first(p.range)
@inline TimeSeriesEcon.lastdate(p::Plan) = last(p.range)

#######################################
# The internal interface to simulations code.

"""
    plansum(model, plan)

Return the total number of exogenous variables in the simulation plan. Periods
over which initial and final conditions are imposed are not counted in this sum.

"""
plansum(m::Model, p::Plan) = sum(p.exogenous[(1 + m.maxlag):(end - m.maxlead), :])

export setexog!
"""
    setexog!(plan, t, vinds)

Modify the plan at time t such that `vinds` are exogenous and the rest are
endogenous.

"""
function setexog!(p::Plan, tt::Int, vinds)
    p.exogenous[tt, :] .= false
    p.exogenous[tt, vinds] .= true
end
setexog!(p::Plan{T}, tt::T, vinds) where T <: MIT = setexog!(p, _offset(p, tt), vinds)

end # module Plans

using .Plans

export Plan,
    exogenize!, endogenize!, exog_endo!, endo_exog!, autoexogenize!


text:
struct MeshPointNM{U,D,C,N,T}
    patch::Simplex{U,D,C,N,T}
    bary::SVector{D,T}
    cart::SVector{U,T}
end

Base.length(m::MeshPointNM) = length(m.cart)
Base.getindex(p::MeshPointNM, i::Int) = p.cart[i]

cartesian(mp::MeshPointNM) = mp.cart
parametric(mp::MeshPointNM) = mp.bary
chart(nbd::MeshPointNM) = mp.patch

"Return the barycentric coordinates of `mp`"
barycentric(mp::MeshPointNM) = SVector(mp.bary[1], mp.bary[2], 1-mp.bary[1]-mp.bary[2])

"""
A number defines a neighborhood in euclidian space
"""
jacobian(x::Number) = one(x)

jacobian(mp::MeshPointNM) = volume(mp.patch) * factorial(dimension(mp.patch))
tangents(mp::MeshPointNM, i) = mp.patch.tangents[i]
normal(mp::MeshPointNM) = mp.patch.normals[1]

function neighborhood(p::Simplex, bary)
  D = dimension(p)
  T = coordtype(p)
  P = SVector{D,T}
  cart = barytocart(p, bary)
  MeshPointNM(p, P(bary), cart)
end

@generated function center(p::Simplex{U,D,C,N,T}) where {U,D,C,N,T}
    uv = ones(T,D)/(D+1)
    :(neighborhood(p, $uv))
end

ext: jl


text:
using TravelingSalesmanExact, GLPK
using TravelingSalesmanHeuristics
using TravelingSalesmanExact: ATT, euclidean_distance
using TravelingSalesmanBenchmarks
using Plots
using Printf, Random
gr(fmt=:svg)

function get_plt_coords(cities)
    n = length(cities)
    inc(a) = a == n ? one(a) : a + 1
    return [cities[inc(j)][1] for j = 0:n], [cities[inc(j)][2] for j = 0:n]
end
TravelingSalesmanExact.plot_cities(cities; kwargs...) = plot(get_plt_coords(cities)...; kwargs...)
plot_cities!(cities; kwargs...) = plot!(get_plt_coords(cities)...; kwargs...)

function plot_tours(cities, pairs; kwargs...)
    plts = []
    colors = sequential_palette(0, length(pairs)+1)[2:end]
    for (index, (tour, label)) in enumerate(pairs)
        plt = plot_cities(cities[tour], label = label, linewidth = 2,
                            color = colors[index])
        plot_cities!(cities; title="Comparison of tours", st=:scatter, label="City locations", markersize = 5, kwargs...)
        push!(plts, plt)
    end
    return plts
end

Random.seed!(155);
N = 50
cities = [100*rand(2) for _ = 1:N]
cost = [ euclidean_distance(cities[i], cities[j]) for i=1:N, j=1:N ]

t_exact, c_exact = get_optimal_tour(cost, with_optimizer(GLPK.Optimizer))

c(q) = solve_tsp(cost; quality_factor = q)[2]

qs = range(10, stop = 100, step = 10)

plot(qs, c, xlabel="quality", ylabel="Cost", label="solve_tsp 1", title="random cities")
for j = 2:5
    plot!(qs, c, label="solve_tsp $j")
end
hline!([c_exact], label="Exact cost")

t_heuristic_100, c_heuristic_100 = solve_tsp(cost; quality_factor = 100)
t_heuristic_40, c_heuristic_40 = solve_tsp(cost; quality_factor = 40)

plts = plot_tours(cities, [
    (t_heuristic_40, @sprintf("Heuristic; quality 40, cost=%.2f",   c_heuristic_40)),
    (t_heuristic_100, @sprintf("Heuristic; quality 100, cost=%.2f", c_heuristic_100)),
    (t_exact, @sprintf("Optimal tour; cost=%.2f", c_exact))
    ])
for plt in plts
    display(plt)
end

function compare_cities(N)
    cities = [100*rand(2) for _ = 1:N]
    cost = [ euclidean_distance(cities[i], cities[j]) for i=1:N, j=1:N ]
    t_exact, c_exact = get_optimal_tour(cost, with_optimizer(GLPK.Optimizer))

    c(q) = solve_tsp(cost; quality_factor = q)[2]

    qs = range(10, stop = 100, step = 10)

    line_plt = plot(qs, c, xlabel="quality", ylabel="Cost", label="solve_tsp 1", title="random cities")
    for j = 2:5
        plot!(qs, c, label="solve_tsp $j")
    end
    hline!([c_exact], label="Exact cost")

    t_heuristic_100, c_heuristic_100 = solve_tsp(cost; quality_factor = 100)
    t_heuristic_40, c_heuristic_40 = solve_tsp(cost; quality_factor = 40)
    city_plts = plot_tours(cities, [
    (t_heuristic_40, @sprintf("Heuristic; quality 40, cost=%.2f",   c_heuristic_40)),
    (t_heuristic_100, @sprintf("Heuristic; quality 100, cost=%.2f", c_heuristic_100)),
    (t_exact, @sprintf("Optimal tour; cost=%.2f", c_exact))
    ])
    return line_plt, city_plts...
end

for j = 1:5
    plts = compare_cities(50)
    for plt in plts
        display(plt)
    end
end

file = @isdefined(WEAVE_ARGS) ? WEAVE_ARGS[:file] : nothing
TravelingSalesmanBenchmarks.bench_footer(file)

max_stars_repo_name: ericphanson/TravelingSalesmanBenchmarks


text:
module MultinomialBayes
  using ...Statistics
  using StatsBase

  function getCounts(dataset)
    return Statistics.sumColumn(dataset)
  end

  function getDictionarySize(dataset)
    return size(dataset)[2]
  end

  function getClassifierCount(labels)
    return countmap(labels)
  end

  function getDictionary(data, labels)
    classifier_counts = getClassifierCount(labels)
    dictionarySize = getDictionarySize(data)
    dictionary = Dict{Integer, Array{Float64,2}}();
    dictionaryProbabilities = Dict{Integer, Array{Float64,2}}()
    for key in keys(classifier_counts)
      dictionary[key] = getCounts(data[find(x-> x==key,labels),:])
    #  println((sum(dictionary[key])))
      dictionaryProbabilities[key] = (dictionary[key].+1)./(sum(dictionary[key])+dictionarySize)
    end

    return dictionary, dictionaryProbabilities
  end

  function classifyDocuments(documents, dictionaryProb)
    predictedLabels = zeros(Float64,size(documents)[1],1)
    for i=1:size(documents)[1]
      probabilities = Dict()
      for key in keys(dictionaryProb)
        println(filter!(x->x!=0,dictionaryProb[key].^documents[i,:]))
        probabilities[key] = prod(dictionaryProb[key].^documents[i,:])
      end
      #println(probabilities)
      predictedLabels[i]=collect(keys(probabilities))[indmax(collect(values(probabilities)))]
    end

    return predictedLabels
  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
using DrWatson
@quickactivate "StatReth"

# %%
using CSV
using DataFrames
using Turing
using Plots

include(srcdir("quap.jl"))
include(srcdir("tools.jl"))

# %%
d = DataFrame(CSV.File(datadir("exp_raw/Howell_1.csv")))
sort!(d, :weight)  # to make the plots look nice
d2 = d[d.age .< 18, :]

z_weight, unz_weight = zscore_transform(d2.weight)
d2.weight_s = z_weight(d2.weight)

@model function height_linear(weight_s, height)
    α ~ Normal(100, 20)
    β ~ LogNormal(2, 1)
    σ ~ Exponential(5)
    μ = f_linear(weight_s, α, β)
    height ~ MvNormal(μ, σ)
end

f_linear(weight_s, α, β) = α .+ β .* weight_s

q = quap(height_linear(d2.weight_s, d2.height))
post = DataFrame(rand(q.distr, 1000)', ["α"; "β"; "σ"])

# %%
estimmean, estimlower, estimupper = estimparam(post)
"""
All values heavily rounded.
α = 108 ± 1 cm
β = 23 ± 1 cm/stddev(kg)
σ = 8.5 ± 0.5 cm
For β we need to convert this result into cm/kg
β = 23 / std(d2.weight) = 2.5

Average height of the children is 108 cm. On average they are 2.5 cm taller for every 1 kg
they are heavier. Children with the same weight have an expected standard deviation of their
height of 8.5 cm
"""

# %%
mu = f_linear(d2.weight_s, post.α', post.β')
mu_m, mu_lower, mu_upper = meanlowerupper(mu)

scatter(d2.weight, d2.height, alpha = 0.3, legend = false)
plot!(d2.weight, mu_m, ribbon = (mu_m .- mu_lower, mu_upper .- mu_m))

pred = rand.(Normal.(mu, post.σ'))
_, pred_lower, pred_upper = meanlowerupper(pred)

plot!(d2.weight, pred_lower, fillrange = pred_upper, alpha = 0.3, linealpha = 0.0, c = 2)

# %%
"""
Well, I'm concerned about the fact that I'm using a linear model even though the data
doesn't seem to be particularly linear. I mean, it's predicting there's going to be a 60 cm
tall child that doesn't weigh anything. Negative weights seen to be a-ok, too.
I suspect/hope that doing something like sqrt/cbrt or the likes would be better.
"""


text:
using Distributions, Plots, LaTeXStrings; pyplot()

mu, sig = 2, 3
eta = sqrt(3)*sig/pi
n, N = 15, 10^4
dNormal   = Normal(mu, sig)
dLogistic = Logistic(mu, eta)
alphaUsed = 0.001:0.001:0.1

function alphaSimulator(dist, n, alpha)
    popVar        = var(dist)
    coverageCount = 0
    for _ in 1:N
        sVar = var(rand(dist, n))
        L = (n - 1) * sVar / quantile(Chisq(n-1),1-alpha/2)
        U = (n - 1) * sVar / quantile(Chisq(n-1),alpha/2)
        coverageCount +=  L < popVar && popVar < U
    end
    1 - coverageCount/N
end

scatter(alphaUsed, alphaSimulator.(dNormal,n,alphaUsed),
	c=:blue, msw=0, label="Normal")
scatter!(alphaUsed, alphaSimulator.(dLogistic, n, alphaUsed),
	c=:red, msw=0, label="Logistic")
plot!([0,0.1],[0,0.1],c=:black, label="1:1 slope",
	xlabel=L"\alpha"*" used", ylabel=L"\alpha"*" actual",
	legend=:topleft, xlim=(0,0.1), ylims=(0,0.2))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
VERSION < v"0.7.0-beta2.199" && __precompile__()

"""
The `VersionParsing` package implements flexible parsing of
version-number strings into Julia's built-in `VersionNumber` type, via
the `vparse(string)` function.

Unlike the `VersionNumber(string)` constructor, `vparse(string)` can
handle version-number strings in a much wider range of formats
without throwing an exception.
"""
module VersionParsing
export vparse

"""
    vparse(s::AbstractString)

Returns a `VersionNumber` representing as much as possible of the version
information in the string `s`.  A much wider range of formats is supported
than the semver styles recognized by `VersionNumber(s)`.

For example,

* Non-numeric prefixes are stripped along with any invalid version characters.
  Commas are treated as decimal points.
* Text following whitespace after the version number is ignored.
* `major.minor.patch.x.y.z` is supported, with `x.y.z` prepended to the
  semver build identifier, i.e. it is parsed like `major.minor.patch+x.y.z`.
* Multiple `+x+y` build identifiers are concatenated as if they were `+x.y`.
* A leading `0` is prepended if needed, e.g. `.x` is treated as `0.x`.
* When all else fails, everything except the first `major.minor.patch`
  digits found are ignored.
"""
vparse(s::AbstractString) = vparse(String(s))

digits2num(s::AbstractString) = all(isdigit, s) ? parse(Int, s) : s
splitparts(s::AbstractString) = map(digits2num, filter!(!isempty, split(s, '.')))

function vparse(s_::String)
    s = replace(s_, r"^[^\d]*[^.\d](\.?\d)"=>s"\1") # strip non-numeric prefix
    s = replace(s, ','=>".") # treat , as . (e.g MS resource-file syntax)
    isempty(s) && throw(ArgumentError("non-numeric version string $s_"))
    occursin(r"^\.\d", s) && (s = "0" * s) # treat .x as 0.x
    if occursin(r"^\d:\d+", s) # debian-style version number
        s = replace(s, r"^\d:"=>"") # strip epoch
    end
    i = findfirst(isspace, s)
    if i !== nothing
        s = s[1:prevind(s,i)]
    end
    s = replace(s, r"[^A-Za-z0-9.+-]*"=>"") # strip invalid version chars
    try
        return VersionNumber(s) # first look for semver-compliant version
    catch
        m = match(r"^(\d+)(\.\d+)?(\.\d+)?", s)
        m === nothing && throw(ArgumentError("unrecognized version string $s_"))
        major = parse(Int, m[1])
        minor = m[2] === nothing ? 0 : parse(Int, m[2][2:end])
        patch = m[3] === nothing ? 0 : parse(Int, m[3][2:end])
        pre = build = ()
        s = s[m.offset + lastindex(m.match):end] # remaining string after the match
        @assert !isempty(s) # otherwise VersionNumber(s) would have succeeded
        if s[1] == '.' # treat following x.y.z as build part
            m = match(r"^(\.[0-9a-zA-Z]*)+", s)
            build = (build..., splitparts(m.match)...)
            s = s[m.offset + lastindex(m.match):end]
        end
        splitplus = split(s, '+')
        if !isempty(splitplus[1])
            pre = (pre..., splitparts(splitplus[1][1]=='-' ? splitplus[1][2:end] : splitplus[1])...)
        end
        for j = 2:length(splitplus) # build version(s)
            build = (build..., splitparts(splitplus[j])...)
        end
        return VersionNumber(major, minor, patch, pre, build)
    end
end

end # module


text:
include("benchmarks.jl")
RESULTS = run(SUITE; verbose = true)
display(RESULTS)


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities

mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
config = MOIT.TestConfig()

@testset "QP" begin
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [4 / 7, 3 / 7, 6 / 7],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [5 / 7, 6 / 7],
        ),
    )
    MOIT.qp1test(mock, config)
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [4 / 7, 3 / 7, 6 / 7],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [5 / 7, 6 / 7],
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [4 / 7, 3 / 7, 6 / 7],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [10 / 7, 12 / 7],
        ),
    )
    MOIT.qp2test(mock, config)
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1 / 4, 3 / 4],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [11 / 4],
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1.0, 0.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [-2.0],
        ),
    )
    MOIT.qp3test(mock, config)
end
@testset "QCP" begin
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1 / 2, 7 / 4],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                [zeros(2)],
            (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1.0],
        ),
    )
    MOIT.qcp1test(mock, config)
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [√2],
            (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1 / (2 * √2)],
        ),
    )
    MOIT.qcp2test(mock, config)
    MOIT.qcp3test(mock, config)
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1.0, 1.0],
            (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1 / 3],
        ),
    )
    MOIT.qcp4test(mock, config)
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1.0, 1.0],
            (MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}) => [1 / 3],
        ),
    )
    MOIT.qcp5test(mock, config)
end
@testset "Non-convex QCP" begin
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, [4.0, 1.0], MOI.FEASIBLE_POINT),
    )
    MOIT.ncqcp1test(mock, config)
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, [2.0, 2.0], MOI.FEASIBLE_POINT),
    )
    MOIT.ncqcp2test(mock, config)
end
@testset "SOCP" begin
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1 / 2, 1 / 2, 1 / √2],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [1 / √2],
            (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1 / √2],
        ),
    )
    MOIT.socp1test(mock, config)
end


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
import StatsBase: fit!, predict
using StatsBase
using BSON
using Flux
using MLDataPattern
using IPMeasures
using GenerativeModels
using Distributions
using LinearAlgebra
using ValueHistories
using FileIO
using DistributionsAD

"""
	save_results(parameters, training_info, results, savepath, data, ac, modelname, seed, dataset,
		contamination)

this computes and saves score and model files
"""
function save_results(parameters, training_info, results, savepath, data, ac, modelname, seed, dataset,
	contamination)
	# save the model separately
	tagsave(joinpath(savepath, savename("model", parameters, "bson", digits=5)),
		Dict("model"=>training_info.model,
			 "tr_encodings"=>get(training_info, :tr_encodings, nothing),
			 "val_encodings"=>get(training_info, :val_encodings, nothing),
			 "tst_encodings"=>get(training_info, :tst_encodings, nothing),
			 "fit_t"=>training_info.fit_t,
			 "history"=>training_info.history,
			 "parameters"=>parameters
			 ),
		safe = true)
	training_info = merge(training_info,
		(model=nothing,tr_encodings=nothing,val_encodings=nothing,tst_encodings=nothing))

	# here define what additional info should be saved together with parameters, scores, labels and predict times
	save_entries = merge(training_info, (modelname = modelname, seed = seed, dataset = dataset,
		anomaly_class = ac,
		contamination=contamination))

	# now loop over all anomaly score funs
	for result in results
		GenerativeAD.experiment(result..., data, savepath; save_entries...)
	end
end


text:
@testset "booleans" begin
    line1 = LineString([[9.170356, 45.477985], [9.164434, 45.482551], [9.166644, 45.484003]])
    line2 = LineString([[9.169356, 45.477985], [9.163434, 45.482551], [9.165644, 45.484003]])

    line3 = LineString([[
						-111.544189453125,
						24.186847428521244
					],
					[
						-110.687255859375,
						24.966140159912975
					],
					[
						-110.4510498046875,
						24.467150664739002
					],
					[
						-109.9951171875,
						25.180087808990645
					]])

	line4 = LineString([[
						-111.4617919921875,
						24.05148034322011
					],
					[
						-110.8795166015625,
						24.681961205014595
					],
					[
						-110.841064453125,
						24.14174098050432
					],
					[
						-109.97863769531249,
						24.617057340809524
					]])

    @test parallel(line1, line2) == true
	@test parallel(line3, line4) == false

	poly1 = Polygon([[[0, 0], [1, 0], [1, 1], [0.5, 0.5], [0, 1], [0, 0]]])
	poly2 = Polygon([[[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]]])

	@test concave(poly1) == true
	@test concave(poly2) == false

	l1 = LineString([[0, 0], [1, 1], [1, 0], [0, 0]])
	l2 = LineString([[0, 0], [1, 0], [1, 1], [0, 0]])

	@test clockwise(l1) == true
	@test clockwise(l2) == false

	l3 = LineString([[0, 0], [3, 3], [4, 4]])
	p1 = Point([1,1])

	l4 = LineString([[0, 0], [3, 3]])
	p2 = Point([0, 0])

	p3 = Point([20, 20])
	l5 = LineString([[0, 0], [3, 3], [38.32, 5.96]])

	@test point_on_line(p2, l4, true) == false
	@test point_on_line(p3, l5, true) == false
	@test point_on_line(p1, l3) == true

	pt = Point([-77, 44])
	poly = Polygon([[ [-81, 41], [-81, 47], [-72, 47], [-72, 41], [-81, 41]]])

	@test point_in_polygon(pt, poly) == true

	poly3 = Polygon([[[1, 1], [1, 10], [10, 10], [10, 1], [1, 1]]])
	poly4 = Polygon([[[1, 1], [2, 2], [3, 2], [1, 1]]])
	line5 = LineString([[1, 1], [2, 3], [2, 3.5]])

	line6 = LineString([[1, 1], [1, 2], [1, 3], [1, 4]])
	poly5 = Polygon( [[[1, 1], [1, 20], [1, 3], [1, 4], [1, 1]]])
	line7 = LineString([[1, 2], [1, 3], [1, 3.5]])

	@test contains(poly3, poly4) == true
	@test contains(poly3, line5) == true
	@test contains(line6, Point([1, 2])) == true
	@test contains(poly3, poly5) == false
	@test contains(poly3 , line7) == false

	@test within(poly4, poly3) == true
	@test within(line5, poly3) == true
	@test within(poly5, poly3) == false
	@test within(Point([1, 2]), line6) == true
	@test within(line7, poly3) == false

	poly6 = Polygon([[[-11, -12], [-13, -12], [-13, -13], [-11, -13], [-11, -12]]])
	poly7 = Polygon([[[-1, 2], [3, 2], [3, 3], [-1, 3], [-1, 2]]])

	#@test disjoint(poly7, poly6) == true # <- fails
	@test disjoint(poly7, Point([1, 1])) == true
	@test disjoint(poly7, LineString([[0, 0], [12, 2], [12, 3], [12, 4]])) == true

	poly8 = Polygon([[[-1, 2], [-13, -12], [-13, -13], [-11, -13], [-1, 2]]])

	@test disjoint(poly8, poly7) == false

	line8 = LineString([[124.584961,-12.768946],[126.738281,-17.224758]])
	line9 = LineString([[123.354492,-15.961329],[127.22168,-14.008696]])

	@test line_intersects(line8, line9).coordinates ≈ [125.583754,-14.835723]

	line10 = LineString([
					[
						142.03125,
						-11.695273
					],
					[
						138.691406,
						-16.804541
					],
					[
						136.40625,
						-14.604847
					],
					[
						135.966797,
						-12.039321
					],
					[
						131.308594,
						-11.436955
					],
					[
						128.232422,
						-15.36895
					],
					[
						125.947266,
						-13.581921
					],
					[
						121.816406,
						-18.729502
					],
					[
						117.421875,
						-20.632784
					],
					[
						113.378906,
						-23.402765
					],
					[
						114.169922,
						-26.667096
					]])

	line11 = LineString([[
						117.861328,
						-15.029686
					],
					[
						122.124023,
						-24.886436
					],
					[
						132.583008,
						-22.309426
					],
					[
						132.890625,
						-7.754537]])

	points = line_intersects(line10, line11)
	@test points[1] ≈ [132.808697,-11.630938]
	@test points[2] ≈ [119.832884,-19.58857]

	@test crosses(LineString([[-2, 2], [4, 2]]), line6) == true
	@test crosses(LineString([[0.5, 2.5], [1, 1]]), poly7) == true
	@test crosses(MultiPoint([[1, 2], [12, 12]]), LineString([[1, 1], [1, 2], [1, 3], [1, 4]])) == true
	@test crosses(MultiPoint([[1, 0], [12, 12]]), LineString([[1, 1], [1, 2], [1, 3], [1, 4]])) == false
	@test crosses(LineString([[-2, 2], [-4, 2]]), poly7) == false

    pl1 = Polygon([[[0,0],[0,5],[5,5],[5,0],[0,0]]])
	pl2 = Polygon([[[1,1],[1,6],[6,6],[6,1],[1,1]]])

	@test overlap(pl1, pl2) == true
	@test_throws ErrorException overlap(pl1, Point([1, 1]))
	@test_throws ErrorException overlap(Point([1, 1]), pl2)

	pl3 = pl4 = Polygon([
					[
						[
							-53.57208251953125,
							28.287451910503744
						],
						[
							-53.33038330078125,
							28.29228897739706
						],
						[
							-53.34136962890625,
							28.430052892335723
						],
						[
							-53.57208251953125,
							28.287451910503744
						]
					]
				])
	@test overlap(pl3, pl4) == false

	mp1 = MultiPoint([
					[
						-36.05712890625,
						26.480407161007275
					],
					[
						-35.7220458984375,
						27.137368359795584
					],
					[
						-35.13427734375,
						26.83387451505858
					],
					[
						-35.4638671875,
						27.254629577800063
					],
					[
						-35.5462646484375,
						26.86328062676624
					],
					[
						-35.3924560546875,
						26.504988828743404
					]
				])
	mp2 = MultiPoint([
					[
						-35.4638671875,
						27.254629577800063
					],
					[
						-35.5462646484375,
						26.86328062676624
					],
					[
						-35.3924560546875,
						26.504988828743404
					],
					[
						-35.2001953125,
						26.12091815959972
					],
					[
						-34.9969482421875,
						26.455820238459893
					]
				])

	@test overlap(mp1, mp2) == true
	@test overlap(mp1, mp2) == overlap(mp2, mp1)
end


text:
# Packages used in this chapter
# To check for installed packages in V1.0 we now need to use the PKG API

using Pkg
ipks = Pkg.installed()
pkgs = ["BenchmarkTools","MacroTools", "Lazy"]
for p in pkgs
    !haskey(ipks,p) && Pkg.add(p)
end

# We can define these here or as required
using BenchmarkTools, MacroTools, Lazy, PyPlot

# We need these as some functions have moved from Base to Stdlib
#
using Printf, Statistics, SpecialFunctions

# We can define these here or as required
using BenchmarkTools, MacroTools, PyPlot

# We need these as some functions have moved from Base to Stdlib
#
using Printf, Statistics

recip(x::Number) =
(x == zero(typeof(x))) ? error("Invalid reciprocal") : one(typeof(x)) / x

recip(2)

recip(recip(2))

recip(11//17)

recip(11 + 17im)

aa = rand(3)
recip(aa)

## recip(a::Array) = [a[i] = recip(a[i]) for i = 1:length(a)]
## OR

recip(a::Array) = map(recip,a)
recip(aa)

map(sin,recip(aa)) # and can map functions to the array

bb = [2.1 3.2 4.3; 9.8 8.7 7.6]

recip(bb) # our definition works but returns a 2-D array

cc = recip(aa)'.*recip(bb)  # however we can do still do this.

incr(x) = x + 1

@code_native incr(2)

@code_native incr(2.7)

@code_native incr(2//7)

@code_native incr(2.0 + 7.0im)

# Look at some intermediate stages
dump(:(incr(2.7)))

@code_lowered incr(2.7)

@code_typed incr(2.7)

@code_llvm(incr(2.7))

ex1 = :((x^2 + y^2 - 2*x*y)^0.5)

ex2 = quote
    (x^2 + y^2 - 2*x*y)^0.5
end

# If we instantiate the variables x,y,z then we can evaluate the expression
# This is NOT a function call, so values need to be known before hand
#
x = 1.1; y = 2.5;
eval(ex1)

# Note that the @eval macro does NOT behave as the function call
@eval ex1

# Unless the expression is prefixed with '$'
@eval $ex1

# i.e.  @eval $(ex) === eval(ex)

# Although ex1 and ex2 are equivalent the AST of the expression
# is slightly different
# Here it is for ex1, we leave ex2 to the reader
#
dump(ex1)

# Dump is very verbose, so using show_expr creates an S-expression version
#
Meta.show_sexpr(ex1)

# Traversing a tree
#
function traverse!(ex, symbols) end

function traverse!(ex::Symbol, symbols)
    push!(symbols, ex)
end

function traverse!(ex::Expr, symbols)
    if ex.head == :call  # function call
        for arg in ex.args[2:end]
            traverse!(arg, symbols)  # recursive
        end
    else
        for arg in ex.args
            traverse!(arg, symbols)  # recursive
        end
    end
end

function traverse(ex::Expr)
    symbols = Symbol[]
    traverse!(ex, symbols)
    return unique(symbols)  # Don't output duplicate
end

traverse(ex1)

macro pout(n)
    if typeof(n) == Expr
       println(n.args)
    end
    return n
end

@pout x

@pout (x^2 + y^2 - 2*x*y)^0.5

macro dotimes(n, body)
    quote
        for i = 1:$(esc(n))
            $(esc(body))
        end
    end
end

@dotimes 3 print("Hi")

i = 0; @dotimes 3 [global i += 1; println(i*i)]

# Expand the @assert function
macroexpand(Main,:(@assert n > 0))

# Expand our @dotimes function, whic is somewhat simpler
#
macroexpand(Main,:(@dotimes 3 [global i += 1; println(i*i)]))

# Not all macro expansions produce short boiler plate
#
using Printf, Statistics
aa = [rand() for i = 1:100000]
@printf "The average value is %f.3 over %d trials"  mean(aa) length(aa)

# A more useful version of dotimes macro is until
# Creates a loop and breaks out when a condition fails
#
macro until(condition, block)
    quote
        while true
            $(esc(block))
            if $(esc(condition))
                break
            end
        end
    end
end

i = 0; @until (i >= 3) [global i += 1; println(i*i)]

# The until macro can be used to implement a simple if-then-else
#
macro iif(cond, body1, body2)
    :(if !$cond
        $(esc(body1))
    else
        $(esc(body2))
    end)
end

fac(n::Integer) = (n == 1) ? 1 : n*fac(n-1)

n = 10
@iif (n < 1) fac(n) ArgumentError("$n not positive")

n = -1; @iif (n < 1) fac(n) ArgumentError("$n not positive")

function fib(n::Integer)
  @assert n > 0
  (n == 1 || n == 2) ? 1 : fib(n-1) + fib(n-2)
end

macroexpand(Main,:(@time(fac(big(402)))))

# A timing macro

# First defined a modified Kempner function
# This converts very slowly

function kempner(n::Integer)
  @assert n > 0
  s = 0.0
  r9  = r"9"
  r9x = r"9{2}"
  for i in 1:n
    if (match(r9,string(i)) == nothing) || (match(r9x,string(i)) != nothing)
      s += 1.0/float(i)
    end
  end
  return s
end

[kempner(10^i) for i in 1:7]

macro bmk(fex, n::Integer)
    quote
      let s = 0.0
      if $(esc(n)) > 0
        val = $(esc(fex))
        for i = 1:$(esc(n))
          local t0 = Base.time_ns()
          local val = $(esc(fex))
          s += Base.time_ns() - t0
        end
        return s/($(esc(n)) * 10e9)
      else
        Base.error("Number of trials must be positive")
      end
    end
  end
end

@bmk kempner(10^7) 10

# Evaluate a polynomial with a function
# This is NOT a macro version

function poly_native(x, a...)
  p=zero(x)
  for i = 1:length(a)
    p = p + a[i] *  x^(i-1)
  end
  return p
end

f_native(x) = poly_native(x,1,2,3,4,5)
f_native(2.1)

# Neither IS this!

function poly_horner(x, a...)
  b=zero(x)
  for i = length(a):-1:1
    b = a[i] + b * x
  end
  return b
end

# f -> (((5*x + 4)*x + 3)*x + 2)*x + 1
#
f_horner(x) = poly_horner(x,1,2,3,4,5)
f_horner(2.1)

# for the macro version we need a help function
# Define mad(x,a,b)
# [Julia has this function too :- muladd(x,a,b)]

mad(x,a,b) = a*x + b
mad(2.1,5,4)

# And NOW use Horner's method in a macro
# [ p... is a variable list of arguments passed as an array]

macro horner(x, p...)
    ex = esc(p[end])
    for i = length(p)-1:-1:1
        ex = :(mad(t, $ex, $(esc(p[i]))))
    end
    Expr(:block, :(t = $(esc(x))), ex)
end

@horner 2.1 1 2 3 4 5

# The saving of eliminating the loop for a 5th order polynomial is not great
# However for larger arrays and complex calculations this can be substantial
#
macroexpand(Main,:(@horner 2.1 1 2 3 4 5))

macro reverse(ex)
    if isa(ex, Expr) && ex.head == :call
        return Expr(:call, ex.args[1], reverse(ex.args[2:end])...)
    else
        return ex
    end
end

@reverse  8//11 - 5//11 * 3//11

# i.e evaluates as :
3//11 * 5//11 - 8//11

macro q(s)
  s0 = eval(s)
  try
    if length(s0) > 0
      ss = split(ss,'\n')
      for i = 1:length(ss)
         println(reverse(ss[i]))
      end
    end
  catch
    return
  end
end

# The postwalk function splits an expression into symbols and then
# reconstructs it, so we can apply different operations to each symbol

using MacroTools:postwalk

ex = :(1 + (2 + 3) + 4)
p = postwalk(ex) do x
    x isa Integer ? fac(x) : x
end

# Evaluate the expression
eval(p)

map(x -> @show(x), [1,2,3,4]);

postwalk(ex) do x
         @show x
end

@capture(ex, a_ + b_ + c_)

b

a*eval(b) + c    # => 1*5 + 4

reduce(+, 1:10)

plus(a, b) = :($a + $b)
p = reduce(plus, 1:10)

eval(p)

# Do something useful
# This is the series expansion for the SINE function
# We need the factorial function and will use the version from the STDLIB
#
using SpecialFunctions
k = 2
pp = [:($((-1)^k) * x^$(1+2k) / $(factorial(1+2k))) for k = 0:5]

# We can reduce this to a single expression
reduce(plus,pp)

# ... and evaluate it for a specific value of x
x = 2.1; eval(reduce(plus,pp))

using Lazy
import Lazy: cycle, range, drop, take

# Create a list of Fibonacci numbers and pick off first 20, using the @lazy macro.
#
fibs = @lazy 0:1:(fibs + drop(1, fibs));
take(20, fibs)

# Function style, pass the argument to a function
#
@> π/6 sin exp   # ==> exp(sin(π/6))

# Currying: The @> can also have functional arguments
#
f(x,μ) = -(x - μ)^2
@> π/6 f(1.6) exp

# The @>> macro reverse the order of the arguments
# Use this toutput the first 15 EVEN squares
#
esquares = @>> range() map(x -> x^2) filter(iseven);
take(15, esquares)

# Next create a list of primes
# The takewhile function is defined in Lazy
#
isprime(n) =
  @>> primes begin
    takewhile(x -> x<=sqrt(n))
    map(x -> n % x == 0)
    any; !
  end;

# We need to initialise the primes list
#
primes = filter(isprime, range(2));
take(20, primes)

isprime(113)

# A simple generated function to execute a multiply + add
# The result returned is a symbol
#
@generated function mad(a,b,c)
    Core.println("Calculating: a*b + c")
    return :(a * b + c)
end

# Call the function
mad(2.3,1.7,1.1)

# ... and again
# ... this time it is not re-evaluated
mad(2.3,1.7,1.1)

# With different types of arguments - the function is reevaluated
mad(2.3,1.7,1)

# To illustrate the use consider the following function ...
# ... which multiplies the size of dimensions of a n-D array
#
# First a 'normal' function definition

function pdims(x::Array{T,N}) where {T,N}
  s = 1
  for i = 1:N
    s = s * size(x, i)
  end
  return s
end

# And then a 'generated function' version
@generated function gpdims(x::Array{T,N}) where {T,N}
  ex = :(1)
  for i = 1:N
     ex = :(size(x, $i) * $ex)
  end
  return ex
end

# We need an array to test the function
aa = reshape([rand() for i = 1:1000],10,5,5,4); size(aa)

# Both functions provide the same result
#
pdims(aa) == gpdims(aa)

# But if we look at the lowered code
@code_lowered pdims(aa)

# But if we look at the lowered code
@code_lowered gpdims(aa)

# Which therefore results in very different generated native code
@code_native pdims(aa)

# This version does not have the looping
@code_native gpdims(aa)

# Previous versions of Julia had an example section, absent in version v1.0
# An interesting one was modular integer, which required a little re-editing
#
module ModInts
export ModInt

import Base: +, -, *, /, inv

struct ModInt{n} <: Integer
  k::Int
  ModInt{n}(k) where n = new(mod(k,n))
end

Base.show(io::IO, k::ModInt{n}) where n =
    print(io, get(io, :compact, false) ? k.k : "$(k.k) mod $n")

(a::ModInt{n} + b::ModInt{n}) where n = ModInt{n}(a.k+b.k)
(a::ModInt{n} - b::ModInt{n}) where n = ModInt{n}(a.k-b.k)
(a::ModInt{n} * b::ModInt{n}) where n = ModInt{n}(a.k*b.k)
-(a::ModInt{n}) where n = ModInt{n}(-a.k)

inv(a::ModInt{n}) where n = ModInt{n}(invmod(a.k, n))
(a::ModInt{n} / b::ModInt{n}) where n = a*inv(b)

Base.promote_rule(::Type{ModInt{n}}, ::Type{Int}) where n = ModInt{n}
Base.convert(::Type{ModInt{n}}, i::Int) where n = ModInt{n}(i)

end

# Test the module
#
using Main.ModInts
m1 = ModInt{11}(2)
m2 = ModInt{11}(7)
m3 = 3*m1 + m2  # => mod(13,11) => 2

# Because of multiple dispatch we can do the following
#
mm = reshape([ModInt{11}(rand(0:10)) for i = 1:100],10,10)
ma = [ModInt{11}(rand(0:10)) for i = 1:10]

mm.*ma'

## An Ordered Pair Module
#
module OrdPairs

import Base: +,-,*,/,==,!=,>,<,>=,<=
import Base: abs,conj,inv,zero,one,show
import LinearAlgebra: transpose,adjoint,norm

export OrdPair

struct OrdPair{T<:Number}
    a::T
    b::T
end

OrdPair(x::Number) = OrdPair(x, zero(T))

value(u::OrdPair)   = u.a
epsilon(u::OrdPair) = u.b

zero(::Type{OrdPairs.OrdPair}) = OrdPair(zero(T),zero(T))
one(::Type{OrdPairs.OrdPair})  = OrdPair(one(T),zero(T))

abs(u::OrdPair)  = abs(value(u))
norm(u::OrdPair) = norm(value(u))

+(u::OrdPair, v::OrdPair) = OrdPair(value(u) + value(v), epsilon(u) + epsilon(v))
-(u::OrdPair, v::OrdPair) = OrdPair(value(u) - value(v), epsilon(u) - epsilon(v))
*(u::OrdPair, v::OrdPair) = OrdPair(value(u)*value(v), epsilon(u)*value(v) + value(u)*epsilon(v))
/(u::OrdPair, v::OrdPair) = OrdPair(value(u)/value(v),(epsilon(u)*value(v) - value(u)*epsilon(v))/(value(v)*value(v)))

==(u::OrdPair, v::OrdPair) = norm(u) == norm(v)
!=(u::OrdPair, v::OrdPair) = norm(u) != norm(v)
>(u::OrdPair, v::OrdPair)  = norm(u) > norm(v)
>=(u::OrdPair, v::OrdPair) = norm(u) >= norm(v)
<(u::OrdPair, v::OrdPair)  = norm(u) < norm(v)
<=(u::OrdPair, v::OrdPair) = norm(u) <= norm(v)

+(x::Number, u::OrdPair) = OrdPair(value(u) + x, epsilon(u))
+(u::OrdPair, x::Number) = x + u

-(x::Number, u::OrdPair) = OrdPair(x - value(u), epsilon(u))
-(u::OrdPair, x::Number) = OrdPair(value(u) - x, epsilon(u))

*(x::Number, u::OrdPair) = OrdPair(x*value(u), x*epsilon(u))
*(u::OrdPair, x::Number) = x*u

/(u::OrdPair, x::Number) = (1.0/x)*u

conj(u::OrdPair)  = OrdPair(value(u),-epsilon(u))
inv(u::OrdPair)   = one(OrdPair)/u

transpose(u::OrdPair) = u
transpose(uu::Array{OrdPairs.OrdPair,2}) = [uu[j,i] for i=1:size(uu)[1],j=1:size(uu)[2]]
adjoint(u::OrdPair) = u

convert(::Type{OrdPair}, x::Number) = OrdPair(x,zero(x))
promote_rule(::Type{OrdPair}, ::Type{<:Number}) = OrdPair

## show(io::IO,u::OrdPair) = print(io,value(u)," + (",epsilon(u),")ϵ")

function show(io::IO,u::OrdPair)
 op::String = (epsilon(u) < 0.0) ? " - " : " + ";
 print(io,value(u),op,abs(epsilon(u)),"ϵ")
end

end

# Exercise the module
# Notice how the show() routine provides pretty-print output

using Main.OrdPairs
p1 = OrdPair(2.3,-1.7)
p2 = OrdPair(4.4,0.9)
p1 * p2

p2/p1

# Again we can operate on array of Ordered Pairs

using Statistics
pp = [OrdPair(rand(),rand()) for i in 1:100]
mean(pp)

# And promote a rational in a mixed OP
p3 = OrdPair(2.3, 11/7)

# The module is not a full implementation
# So this FAILS
p4 = OrdPair(2.3, 11.0 + 7.2im)

# Also a number of other arithmetic functions need to be imported from Base
std(p1)


text:
using JUGRNModelGenerator

# setup paths -
path_to_input_file = "/Users/jeffreyvarner/Desktop/julia_work/JUGRNModelGenerator.jl/test/data/Model.json"
path_to_output_dir = "/Users/jeffreyvarner/Desktop/julia_work/JUGRNModelGenerator.jl/test/generated_code"

# build the model object -
julia_model_object = build_julia_model_object(path_to_input_file, path_to_output_dir) |> check

# generate the ir dictionary -
ir_dictionary = read_model_document(julia_model_object) |> check


text:
module TestSources

using MLJBase
using Test

X = 7

Xs = source(X)
@test Xs() == X
@test Xs(8) == 8
@test elscitype(Xs) == Count
@test scitype(Xs) == MLJBase.CallableReturning{Count}
rebind!(Xs, nothing)
@test isempty(Xs)
@test Xs.scitype == Nothing

end
true


text:
export PDPage,
    pdPageGetContents,
    pdPageIsEmpty,
    pdPageGetCosObject,
    pdPageGetContentObjects,
    pdPageGetFonts,
    pdPageGetMediaBox,
    pdPageGetCropBox,
    pdPageExtractText,
    pdPageGetPageNumber,
    pdPageEvalContent

using ..Cos

"""
```
    pdPageGetCosObject(page::PDPage) -> CosObject
```
PDF document format is developed in two layers. A logical PDF document
information is represented over a physical file structure called COS. This method
provides the internal COS object associated with the page object.
"""
pdPageGetCosObject(page::PDPage) = page.cospage

"""
```
    pdPageGetContents(page::PDPage) -> CosObject
```
Page rendering objects are normally stored in a `CosStream` object in a PDF file.
This method provides access to the stream object.

Please refer to the PDF specification for further details.

# Example
```
julia> pdPageGetContents(page)

448 0 obj
<<
	/Length	437
	/FFilter	/FlateDecode
	/F	(/tmp/tmpZnGGFn/tmp5J60vr)
>>
stream
...
endstream
endobj
```
"""
function pdPageGetContents(page::PDPage)
    if (page.contents === CosNull)
        ref = get_page_content_ref(page)
        page.contents = get_page_contents(page, ref)
    end
    return page.contents
end

"""
```
    pdPageGetMediaBox(page::PDPage) -> CDRect{Float32}
    pdPageGetCropBox(page::PDPage) -> CDRect{Float32}
```
    Returns the media box associated with the page. See 14.11.2 PDF 1.7 Spec.
It's typically, the designated size of the paper for the page. When a crop box
is not defined, it defaults to the media box.

# Example
```
julia> pdPageGetMediaBox(page)
Rect:[0.0 0.0 595.0 792.0]

julia> pdPageGetCropBox(page)
Rect:[0.0 0.0 595.0 792.0]
```
"""
function pdPageGetMediaBox(page::PDPage)
    arr = page_find_attribute(page, cn"MediaBox")::CosArray
    return CDRect{Float32}(CDRect(arr))::CDRect{Float32}
end

function pdPageGetCropBox(page::PDPage)
    box = page_find_attribute(page, cn"CropBox")
    box === CosNull && return pdPageGetMediaBox(page)
    return CDRect{Float32}(CDRect(box))::CDRect{Float32}
end

"""
```
    pdPageIsEmpty(page::PDPage) -> Bool
```
Returns `true` when the page has no associated content object.

# Example
```
julia> pdPageIsEmpty(page)
false
```
"""
function pdPageIsEmpty(page::PDPage)
    return page.contents === CosNull && get_page_content_ref(page) === CosNull
end

"""
```
    pdPageGetContentObjects(page::PDPage) -> CosObject
```
Page rendering objects are normally stored in a `CosStream` object in a PDF file.
This method provides access to the stream object.
"""
function pdPageGetContentObjects(page::PDPage)
    page.content_objects === nothing && load_page_objects(page)
    return page.content_objects
end

"""
```
    pdPageGetFonts(page::PDPage) -> Dict{CosName, PDFont}()
```
Returns a dictionary of fonts in the page.

#Example

```
julia> pdPageGetFonts(page)
Dict{CosName,PDFIO.PD.PDFont} with 4 entries:
  /F0 => PDFont(…
  /F4 => PDFont(…
  /F8 => PDFont(…
  /F9 => PDFont(…

```
"""
function pdPageGetFonts(page::PDPage)
    cosfonts = find_resource(page, cn"Font", CosNull)
    dres = Dict{CosName, PDFont}()
    for (name, val) in cosfonts.val
        dres[name] = PDFont(page.doc, val)
    end
    return dres
end

function pdPageEvalContent(page::PDPage, state::GState=GState{:PDFIO}())
    state[:source] = page
    evalContent!(pdPageGetContentObjects(page), state)
    return state
end

"""
```
    pdPageExtractText(io::IO, page::PDPage) -> IO
```
Extracts the text from the `page`. This extraction works best for tagged PDF
files.
For PDFs not tagged, some line and word breaks will not be extracted properly.

# Example

Following code will extract the text from a full PDF file.

```
function getPDFText(src, out)
    doc = pdDocOpen(src)
    docinfo = pdDocGetInfo(doc)
    open(out, "w") do io
		npage = pdDocGetPageCount(doc)
        for i=1:npage
            page = pdDocGetPage(doc, i)
            pdPageExtractText(io, page)
        end
    end
    pdDocClose(doc)
    return docinfo
end
```
"""
function pdPageExtractText(io::IO, page::PDPage)
    state = pdPageEvalContent(page)
    show_text_layout!(io, state)
    return io
end

"""
```
    pdPageGetPageNumber(page::PDPage)
```
Returns the page number of the document page.

# Example

```
julia> pdPageGetPageNumber(page)
1
```
"""
pdPageGetPageNumber(page::PDPage) =
    pd_doc_get_pagenum(page.doc, CosIndirectObjectRef(page.cospage))

mutable struct PDPageImpl <: PDPage
    doc::PDDocImpl
    cospage::ID{CosDict}
    contents::CosObject
    content_objects::Union{Nothing, PDPageObjectGroup}
    fonts::Dict{CosName, PDFont}
    xobjs::Dict{CosName, PDXObject}
    PDPageImpl(doc, cospage, contents) =
        new(doc, cospage, contents,
            nothing,
            Dict{CosName,PDFont}(),
            Dict{CosName,PDXObject}())
end

PDPageImpl(doc::PDDocImpl, cospage::ID{CosDict}) =
    PDPageImpl(doc, cospage, CosNull)

#=This function is added as non-exported type. PDPage may need other attributes
which will make the constructor complex. This is the default with all default
values.
=#
create_pdpage(doc::PDDocImpl, cospage::ID{CosDict}) =
    PDPageImpl(doc, cospage)
create_pdpage(doc::PDDocImpl, cospage::CosNullType) =
    throw(ErorException(E_INVALID_OBJECT))
#=
This will return a CosArray of ref or ref to a stream. This needs to be
converted to an actual stream object
=#
get_page_content_ref(page::PDPageImpl) = get(page.cospage, cn"Contents")

function get_page_contents(page::PDPageImpl, contents::CosArray)
    len = length(contents)
    arr = get(contents)
    for i = 1:len
        ref = splice!(arr, 1)
        cosstm = get_page_contents(page, ref)
        cosstm !== CosNull && push!(arr, cosstm)
    end
    stm = merge_streams(page.doc.cosDoc, contents)
    return stm
end

get_page_contents(page::PDPageImpl, contents::CosIndirectObjectRef) =
    cosDocGetObject(page.doc.cosDoc, contents)

get_page_contents(page::PDPage, obj::IDD{CosStream}) = obj

@inline function load_page_objects(page::PDPageImpl)
    contents = pdPageGetContents(page)
    page.content_objects === nothing &&
        (page.content_objects = PDPageObjectGroup())
    return load_page_objects(page, contents)
end

load_page_objects(page::PDPageImpl, stm::CosNullType) = nothing

@inline function load_page_objects(page::PDPageImpl, stm::IDD{CosStream})
    bufstm = decode(stm)
    try
        load_objects(page.content_objects, bufstm)
    finally
        util_close(bufstm)
    end
    return nothing
end

@inline function load_page_objects(page::PDPageImpl, stms::IDD{CosArray})
    stm = merge_streams(page.doc.cosDoc, stms)
    page.contents = stm
    return load_page_objects(page, stm)
end

function find_resource(page::PDPageImpl,
                       restype::CosName,
                       fontname::Union{CosName, CosNullType})
    res = CosNull
    cosdoc = page.doc.cosDoc
    pgnode = page.cospage

    while ((fontname !== CosNull && res === CosNull) ||
           (fontname === CosNull)) && pgnode !== CosNull

        resref = get(pgnode, cn"Resources")
        if resref === CosNull
            pgnode = cosDocGetObject(cosdoc, pgnode, cn"Parent")
            continue
        end
        resources = cosDocGetObject(cosdoc, resref)
        if resources === CosNull
            pgnode = cosDocGetObject(cosdoc, pgnode, cn"Parent")
            continue
        end
        ress = cosDocGetObject(cosdoc, resources, restype)
        if ress === CosNull
            pgnode = cosDocGetObject(cosdoc, pgnode, cn"Parent")
            continue
        end
        if fontname !== CosNull
            res = cosDocGetObject(cosdoc, ress, fontname)
        else
            resdict = cosDocGetObject(cosdoc, ress, fontname)
            res === CosNull && (res = CosDict())
            for (k, v) in resdict.val
                set!(res, k, v)
            end
        end
        pgnode = cosDocGetObject(cosdoc, pgnode, cn"Parent")
    end
    return res
end

get_font(page::PDPageImpl, fontname::CosName) =
    get!(page.fonts, fontname,
         get_pd_font!(page.doc, find_resource(page, cn"Font", fontname)))

get_xobject(page::PDPageImpl, xobjname::CosName) =
    get!(page.xobjs, xobjname,
         get_pd_xobject!(page.doc,
                         find_resource(page, cn"XObject", xobjname)))

function page_find_attribute(page::PDPageImpl, resname::CosName)
    res = CosNull
    cosdoc = page.doc.cosDoc
    pgnode = page.cospage

    while pgnode !== CosNull
        res = cosDocGetObject(cosdoc, pgnode, resname)
        res !== CosNull && break
        pgnode = cosDocGetObject(cosdoc, pgnode, cn"Parent")
    end
    return res
end

get_encoded_string(s::CosString, fontname::CosNullType, page::PDPage) =
    CDTextString(s)

get_encoded_string(s::CosString, fontname::CosName, page::PDPage) =
    get_encoded_string(s, get(page.fonts, fontname, nothing))


text:
using OceanGrids, Test, Unitful, SparseArrays

@testset "Small grid tests" begin
    nlat, nlon, ndepth = 2, 3, 4
    grd = OceanGrid(nlat, nlon, ndepth)
    @test grd isa OceanRectilinearGrid
    @test grd isa OceanGrid
    @test size(grd) == (nlat, nlon, ndepth)
    @test length(grd) == *(nlat, nlon, ndepth)
    @testset "OceanGridBox" for box in grd
        @test box isa OceanGridBox
        show(box)
        @test OceanGrids.area(box) isa Quantity
        @test unit(OceanGrids.area(box)) == u"km^2"
        @test unit(OceanGrids.volume(box)) == u"m^3"
    end
end

@testset "edges to grid" begin
    elat = -90:90:90
    elats = [-90:90:90,
             -90.0:90.0:90.0,
             [-90,0,90],
             [-90.0,0.0,90.0],
             (-90.0:90.0:90.0) * u"°",
             [-90.0,0.0,90.0] * u"°",
             (-90:90:90) * u"°",
             [-90,0,90] * u"°",]
    elons = [-180:90:180,
             [-180,-90,0,90,180],
             (-180:90:180) * u"°",
             [-180,-90,0,90,180] * u"°"]
    edepths = [[0, 100, 200, 500, 1000, 3000] * u"m",
               [0.0, 100.0, 200.0, 500.0, 1000.0, 3000.0] * u"m",
               [0.0, 100.0, 200.0, 500.0, 1000.0, 3000.0],
               [0, 100, 200, 500, 1000, 3000],
               [0, 0.1, 0.2, 0.5, 1, 3] * u"km"]
    base_grd = OceanGrid(elats[1], elons[1], edepths[1])
    @testset "Base constructors" for elat in elats, elon in elons, edepth in edepths
        grd = OceanGrid(elat,elon,edepth)
        @test grd isa OceanRectilinearGrid
        @test grd isa OceanGrid
        @test grd == base_grd
        show(grd)
        @testset "OceanGridBox" for box in grd
            @test box isa OceanGridBox
        end
    end
    @testset "Functions" begin
        x = ones(count(iswet(base_grd)))
        @testset "base functions" begin
            @test iswet(base_grd) isa BitArray{3}
            @test latvec(base_grd) isa Vector
            @test lonvec(base_grd) isa Vector
            @test depthvec(base_grd) isa Vector
            @test topdepthvec(base_grd) isa Vector
            @test bottomdepthvec(base_grd) isa Vector
            @test volumevec(base_grd) isa Vector
            @test isseafloorvec(base_grd) isa BitArray{1}
            @test latlondepthvecs(base_grd) isa Tuple{T,T,T} where {T<:Vector}
            x3D = rearrange_into_3Darray(x, base_grd)
            @test x3D isa Array{Float64,3}
            x2 = vectorize(x3D, base_grd)
            @test x2 == x
            @test_throws ErrorException vectorize(rand(10,10,10), base_grd)
        end
        @testset "interpolations" begin
        end
        @testset "utils" begin
            @test buildIabove(base_grd) isa SparseMatrixCSC
            @test buildIbelow(base_grd) isa SparseMatrixCSC
        end
        @testset "Slices and views" begin
            @test horizontalslice(x, base_grd, depth=100) isa SubArray{Float64,2}
            @test meridionalslice(x, base_grd, lon=330) isa PermutedDimsArray{Float64,2}
            @test zonalslice(x, base_grd, lat=15) isa PermutedDimsArray{Float64,2}
            @test depthprofile(x, base_grd, lonlat=(330,15)) isa SubArray{Float64,1}
        end
        @testset "Integrals, averages, and stds" begin
            @test ∫dxdy(x, base_grd) isa Vector{T} where {T<:Quantity}
            @test ∫dz(x, base_grd) isa Array{T,2} where {T<:Quantity}
            @test ∫dy(x, base_grd) isa Array{T,2} where {T<:Quantity}
            @test ∫dx(x, base_grd) isa Array{T,2} where {T<:Quantity}
            @test ∫dxdydz(x, base_grd) isa Quantity
            @test horizontalmean(x, base_grd) isa Vector{Float64}
            @test verticalmean(x, base_grd) isa Array{Float64,2}
            @test meridionalmean(x, base_grd) isa Array{Float64,2}
            @test zonalmean(x, base_grd) isa Array{Float64,2}
            @test totalmean(x, base_grd) isa Float64
            @test horizontalstd(x, base_grd) isa Vector{Float64}
            @test verticalstd(x, base_grd) isa Array{Float64,2}
            @test meridionalstd(x, base_grd) isa Array{Float64,2}
            @test zonalstd(x, base_grd) isa Array{Float64,2}
            @test totalstd(x, base_grd) isa Float64
            @test zcumsum(x, base_grd) isa Vector{Float64}
        end
        @testset "Regridding arrays" begin
            lat, lon = collect(-80:20:80), collect(15:30:360)
            y2D = cosd.(lat) * sind.(lon)'
            @test regrid(y2D, lat, lon, base_grd) isa Array{Float64,2}
            @test_broken regridandpaintsurface(y2d, lat, lon, base_grd) isa Array{Float64,3}
            depth = [10, 20, 30, 50, 100]
            y3D = y2D .* reshape(sqrt.(depth), (1,1,length(depth)))
            @test regrid(y3D, lat, lon, depth, base_grd) isa Array{Float64,3}
            y3D_with_nans = copy(y3D)
            y3D_with_nans[1:end÷2] .= NaN
            y_regridded_with_nans = regrid(y3D_with_nans, lat, lon, depth, base_grd; interpolate_nans=true)
            @test y_regridded_with_nans isa Array{Float64,3}
            @test any(isnan, y_regridded_with_nans)
            y_regridded_no_nans = regrid(y3D_with_nans, lat, lon, depth, base_grd; interpolate_nans=false)
            @test y_regridded_no_nans isa Array{Float64,3}
            @test all(!isnan, y_regridded_no_nans)
        end
        @testset "Regridding vectors" begin
            lats, lons, depths = [-10, 10, 80], [-100, -50, 120], [10, 100, 1000]
            vs = [-10, -5, -4]
            foo = regrid(vs, lats, lons, depths, base_grd)
            @test foo isa Vector{eltype(vs)}
        end
    end

end


text:
#"Live-Slice" the results of a parametric sinusoidal "simulation
#-------------------------------------------------------------------------------
module CMDimData_SampleUsage

using CMDimData
using CMDimData.EasyPlot
using MDDatasets
import InspectDR
CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
LBL_AXIS_TIME = "Time (s)"
LBL_AXIS_POSITION = "Position (m)"
LBL_AXIS_NORMPOSITION = "Normalized Pos (m/m)"
LBL_AXIS_SPEED = "Speed (m/s)"
lstylesweep = cons(:a, line = set(style=:solid, width=2)) #Default is a bit thin
dfltglyph = cons(:a, glyph = set(shape=:o, size=1.5))

#==Read in results of parametric "simulation"
===============================================================================#
signal = evalfile("../analysis_fmtfiles/parametric_sin_data.jl")
signal = convert(DataHR, signal) #Need DataHR for Live-Slice

#==Parameter extraction
===============================================================================#

#Generate DataRS with value of amplitude (A) in each parameter combination:
ampvalue = parameter(signal, "A")

#Normalize signal amplitudes for ALL signals SIMULTANEOUSLY:
signal_norm = signal / ampvalue

#Compute continuous-time signal rate for ALL signals SIMULTANEOUSLY:
#   signal(𝜙, A, 𝑓, t) -> rate(𝜙, A, 𝑓):
rate = deriv(signal)

#Compute maximum signal rate for ALL signals SIMULTANEOUSLY:
maxrate = maximum(abs(rate))

#Compute first fall-crossing point for ALL signals SIMULTANEOUSLY:
#   signal(𝜙, A, 𝑓, t) -> fallx(𝜙, A, 𝑓):
fallx = xcross1(signal, xstart=0, allow=CrossType(:fall))

#==Generate plots
===============================================================================#
xred1 = paramlist(fallx)[end] #x-value of 1st reduction operation (finding 1st-xing time)
#Fixed yaxis values to place results in better context:
yext_pos = (min=-5e-3, max=5e-3)
yext_speed = (min=-400, max=400)
yext_firstx = (min=0, max=800e-6)
yext_maxrate = (min=0, max=500)

#plot #1: Complete time-domain signals
#-------------------------------------------------------------------------------
p1 = cons(:plot, nstrips = 3,
	ystrip1 = set(axislabel=LBL_AXIS_POSITION, striplabel="Sinusoidal response"; yext_pos...),
	ystrip2 = set(axislabel=LBL_AXIS_NORMPOSITION, striplabel="Normalized response"),
	ystrip3 = set(axislabel=LBL_AXIS_SPEED, striplabel="Rate of change"; yext_speed...),
	xaxis = set(label=LBL_AXIS_TIME)
)
push!(p1,
	cons(:wfrm, signal, lstylesweep, label="", strip=1),
	cons(:wfrm, signal_norm, lstylesweep, label="", strip=2),
	cons(:wfrm, rate, lstylesweep, label="", strip=3),
)

#plot #2: Extracted parameter values: fallx, maxrate
#-------------------------------------------------------------------------------
p2 = cons(:plot, nstrips = 2, legend=false,
	ystrip1 = set(axislabel=LBL_AXIS_TIME, striplabel="Time to 1st fall-crossing"; yext_firstx...),
	ystrip2 = set(axislabel=LBL_AXIS_SPEED, striplabel="Maximum signal rate"; yext_maxrate...),
	xaxis = set(label=xred1),
)
push!(p2,
	cons(:wfrm, fallx, lstylesweep, dfltglyph, label="", strip=1),
	cons(:wfrm, maxrate, lstylesweep, dfltglyph, label="", strip=2),
)

println()
@info("Displaying plot...")
pcoll = cons(:plotcoll, title="Parametric sin() - Live-Slice Results", ncolumns=2)
	push!(pcoll, p1, p2)
gplot = EasyPlot.displaygui(:InspectDR, pcoll)

#Activate Live-Slice GUI:
println()
@info("Activating Live-Slice GUI...")
include("tools/liveslice_blink.jl")
slicelist = sweeps(fallx)[1:end-1] #Use 2 first dimensions for slice
wfrmlist = append!(append!(EasyPlot.Waveform[], p1.wfrmlist), p2.wfrmlist)
LiveSlice.autoslice(slicelist, wfrmlist) do
	EasyPlot.build(gplot.src, pcoll)
	InspectDR.refresh(gplot)
end

end #module
:SampleCode_Executed


text:
using Printf
using DifferentialEquations
using LSODA
using Statistics
using Random
import Distributions: Uniform, LogNormal, fit
using PyCall
plt = pyimport("matplotlib.pyplot")
vis = PyCall.pyimport("visual") # Import python module for visualization. Module file needs to be in python lib folder.

μ_for_mean(m, σ) = log(m) - σ^2/2 # log-normal distribution

function prob_func_single(prob,i,repeat) # problem function for ODE
    prob.p["x"] = xss[i]
    prob
end

function set_pars(prob, pars) # set parameter values of a problem
    for (k, v) in pars
        prob.p[k] = v
    end
end

function prob_func(prob,i,repeat) # problem function for ensemble
    if rem(i, ncols) == 0
        prob.p["x"] = xss[ncols]
    else
        prob.p["x"] = xss[rem(i, ncols)]
    end
    prob.p["y"] = i ÷ ncols + 1
    prob
end

function transi_width(data) # computer transition width
    x = data[1,:,:,end]
    y = data[2,:,:,end]
    p5_pos = x .> minimum([0.05*maximum(x), 1.0])
    p8_pos = y .> minimum([0.05*maximum(y), 1.0])
    if maximum(x) > 300
        p5_pos .= 1
    end
    if maximum(y) > 300
        p8_pos .= 1
    end
    if all(x->x==1, p5_pos) || all(x->x==0, p5_pos) || all(x->x==1, p8_pos) || all(x->x==0, p8_pos)
        tw5, tw8 = size(data[2]), size(data[2])
    end
    p5_pos_sum = sum(p5_pos, dims=2);
    p8_pos_sum = sum(p8_pos, dims=2);
    tw5 = sum( (p5_pos_sum .< nrows) .* (p5_pos_sum .> 0) );
    tw8 = sum( (p8_pos_sum .< nrows) .* (p8_pos_sum .> 0) );
    twdp = sum( sum(p5_pos .* p8_pos, dims=2) .> 0 )
    twdn = sum( sum((p5_pos .+ p8_pos) .== 0 , dims=2) .> 0 )
    if tw5 == 0 && tw8 == 0
        xm, ym = mean(x, dims=2), mean(y, dims=2)
        if maximum(abs.(xm[2:end] - xm[1:end-1])) < (maximum(xm) - minimum(xm))*0.35
            if maximum(abs.(ym[2:end] - ym[1:end-1])) < (maximum(ym) - minimum(ym))*0.35
            tw5, tw8 = 40, 40
            end
        end
    end
    tw = max(tw5, tw8, twdp, twdn)
    return tw
end

function get_fmt_pars(pd, name_only=false) # print formated parameter values
    if name_only == true
        return rstrip(join([rpad(k,8, " ") for (k,v) in sort(collect(pd)) ], "\t"))
    else
        return rstrip(join([rpad(@sprintf("%.2f", v), 8, " ") for (k,v) in sort(collect(pd)) ], "\t"))
    end
end

mxstep(x) = maximum(abs.(mean(x, dims=2)[2:end] ./ mean(x, dims=2)[1:end-1]))
function get_mpdif(data) # compute the segregation index
    p5, p8 = data[1,:,:,end], data[2,:,:,end]
    m5, m8 = data[3,:,:,end], data[4,:,:,end]
    mpdif = (mxstep(m5) + mxstep(m8)) / (mxstep(p8) + mxstep(p5))
    return mpdif
end

include("./models/model_mmiS.jl"); # Select a model: TCR, TUR, TmiUR, TmiFB, mmiS
p_cp = copy(p); # in case we need to rerun the simulation later
prob = ODEProblem(model, u0, tspan, p) # define ODE problem
xs = LinRange(0, 100, 100) # x coordinates
ncols = 40 # number of columns in the space
xss = LinRange(xs[11], xs[81], ncols) # x coorindates for all columns
nrows = 10 # number of rows
ensemble_prob = EnsembleProblem(prob,prob_func=prob_func) # define ensemble problem
set_pars(prob, p_cp) # in case the parameters were changed in simulations
@time sim = solve(ensemble_prob, lsoda(), reltol=1e-6, abstol=1e-6, saveat=1.0
            ,trajectories=size(xss)[1]*nrows, dt=0.01); # run simulation
num_tps = 100 # number of time points to export
d = reshape(hcat(cat(EnsembleAnalysis.componentwise_vectors_timepoint(sim,LinRange(0,1000, num_tps))..., dims=(2))...), size(u0)[1], ncols, nrows, num_tps);
println("Segregation index: ", get_mpdif(d)) # Segregation index
println("Transition width: ", transi_width(d)) # Transition width
vis.anim(d, svlist) # Visualize simulation result

# Perform random sampling of parameter values
su = include("SamplingUtils.jl") # file storing the ranges of parameter values
p_all = Vector{Dict}()
Random.seed!(2020)
n = 10 # number of samples 10-10000
for i in 1:n
    p_i = Dict()
    for (k, v) in p
        if startswith(k, "K") && v > 10
            p_i[k] = v
            continue
        end
        if haskey(su.p_rng, k)
            prng = su.p_rng[k]
            if isa(prng, Dict)
                m, σ = prng["m"], prng["σ"]
                p_i[k] = rand(LogNormal(μ_for_mean(m, σ), σ))
            elseif isa(prng, Array)
                p_i[k] = prng[rand(1:end)]
            else
                p_i[k] = v
            end
        else
            println("Parameter $k not found in Sampling file.")
        end
    end
    push!(p_all, p_i)
end

num_tps = 5 # output time points

d_all = Array{Float64,5}(undef, size(p_all)[1], size(u0)[1], ncols, nrows, num_tps);

open("perf.txt", "w") do f # output file storing the performance
println(f, " "^20*get_fmt_pars(p, true))

for i in 1:size(p_all)[1]
    p_id = lpad(i,5,"0")
    set_pars(prob, p_all[i])
    sim = solve(ensemble_prob, lsoda(), reltol=1e-6, abstol=1e-6, saveat=1.0
                ,trajectories=size(xss)[1]*nrows, dt=0.01);
    d = reshape(hcat(cat(EnsembleAnalysis.componentwise_vectors_timepoint(sim,LinRange(0,tspan[2], num_tps))..., dims=(2))...), size(u0)[1], ncols, nrows, num_tps);
    d_all[i,:,:,:,:] = d
    #save("./sim_data/sol_data_" * p_id * ".jld", "data", d[:,:,:,:], "parvs",p_all[i])
    tw = transi_width(d)
    mpdif = get_mpdif(d)
    println(f, p_id * "\t" * string(tw) * "\t" * @sprintf("%.2f", mpdif) *"\t"* get_fmt_pars(p_all[i]))
end
end


text:
# This file includes code that was formerly a part of Julia. License is MIT: https://julialang.org/license

## reversed strings without data movement ##

struct RevString{T<:AbstractString} <: AbstractString
    string::T
end

lastindex(s::RevString) = lastindex(s.string)
length(s::RevString) = length(s.string)
sizeof(s::RevString) = sizeof(s.string)

function next(s::RevString, i::Int)
    n = lastindex(s); j = n-i+1
    (s.string[j], n-prevind(s.string,j)+1)
end

codeunit(s::RevString) = codeunit(s.string)
ncodeunits(s::RevString) = ncodeunits(s.string)

if isdefined(Base, :iterate)
    function iterate(s::RevString, i::Int = firstindex(s))
        i > lastindex(s) && return nothing
        return next(s, i)
    end
end

function isvalid(s::RevString, i::Int)
    1 ≤ i ≤ ncodeunits(s) || return false
    j = 1
    while j < i
        _, j = iterate(s, j)
    end
    return j == i
end

reverse(s::RevString) = s.string
reverseind(s::RevString, i::Integer) = lastindex(s) - i + 1


text:
"""
A wrapper around a parameterized function of the form f(θ, x...)
    - θ are parameters to be learned
    - x is per-iteration inputs... there may not be any (for
        example the Rosenbrock function is only f(θ))

The type parameter DIFFNUM is the number of times it can be differentiated.
"""
type Differentiable{T,DIFFNUM,P<:Params} <: Minimizable
    f::Function
    df::NTuple{DIFFNUM,Function}
    input::SumNode{T,1}
    output::OutputNode{T,1}
    params::P
end
Differentiable(args...) = Differentiable(Float64, args...)
function Differentiable{T}(::Type{T}, f::Function, nθ::Int, nx::Int, df::Function...)
    input = InputNode(T, nx)
    output = OutputNode(T, 1)
    params = Params(T, nθ)
    Differentiable(f, df, input, output, params)
end

# a convenience constructor for Differentiable types
function tfunc(f::Function, nΘ::Int, nx::Int, df::Function...; T = Float64)
    Differentiable(T, f, nΘ, nx, df...)
end
tfunc(f::Function, nΘ::Int, df::Function...; kw...) = tfunc(f, nΘ, 0, df...; kw...)

typealias NonDifferentiable{T,P} Differentiable{T,0,P}
typealias OnceDifferentiable{T,P} Differentiable{T,1,P}
typealias TwiceDifferentiable{T,P} Differentiable{T,2,P}

function apply_function(t::Differentiable, f::Function)
    if input_length(t) == 0
        # not a function of x
        f(params(t))
    else
        f(params(t), input_value(t))
    end
end

function transform!(t::Differentiable)
    t.output.val[1] = apply_function(t, t.f)
end

totalcost(t::Differentiable) = t.output.val[1]

grad!(t::NonDifferentiable) = error()

function grad!(t::Differentiable)
    grad(t)[:] = apply_function(t, t.df[1])
    return
end

# -----------------------------------------------------------------------
# some common test functions

module TestTransforms

using ..Transformations

export
    rosenbrock,
    rosenbrock_gradient
    # rosenbrock_transform

function rosenbrock(θ::AbstractVector)
    sum(100(θ[i+1] - θ[i]^2)^2 + (θ[i] - 1)^2 for i=1:length(θ)-1)
end

function rosenbrock(θs::Number...)
    rosenbrock(collect(θs))
end

function rosenbrock_gradient(θ::AbstractVector)
    ∇ = zeros(θ)
    n = length(θ)
    for i=1:n
        if i > 1
            ∇[i] += 200*(θ[i] - θ[i-1]^2)
        end
        if i < n
            ∇[i] += 2*(θ[i]-1) - 400θ[i] * (θ[i+1] - θ[i]^2)
        end
    end
    ∇
end

function rosenbrock_gradient(θs::Number...)
    rosenbrock_gradient(collect(θs))
end

# function rosenbrock_transform(nθ::Int = 2)
#     @assert nθ > 1
#     Differentiable(rosenbrock, 0, nθ, rosenbrock_gradient)
# end

end #TestTransforms


text:
# PREMPolyModels

"""
    PREMPolyModel <: SeisModel1D

Type describing the Earth as a set of layers within which properties vary
according to a set of polynomials.

Physical parameters are represented by arrays of size `(order+1, n)`, where `n`
is the number of layers, and `order` is the order of polynomial which is
used to represent the parameter.  Hence a constant layer has size `(1,n)`,
and to compute the value of `x` in layer `i`, for an Earth radius of `a` km,
at a radius of `r` km, the expression is:

    val_x = x[i,1] + (r/a)*x[i,2] + (r/a)^2*x[i,3] ... (r/a)^order*x[i,order+1]

The reference frequency in Hz for the velocities in case of attenuation is
given by the `fref` field.  If `fref` is `NaN`, then no reference frequency
is defined for the model.
"""
struct PREMPolyModel <: SeisModel1D
    "Earth radius in km"
    a :: Float64
    "Number of distinct layers"
    n :: Int
    "Radii in km of top of each layer"
    r :: Vector{Float64}
    "Physical parameters"
    vp :: Array{Float64,2}
    vs :: Array{Float64,2}
    density :: Array{Float64,2}
    "Optional anisotropic parameters"
    aniso :: Bool
    vph :: Array{Float64,2}
    vpv :: Array{Float64,2}
    vsh :: Array{Float64,2}
    vsv :: Array{Float64,2}
    eta :: Array{Float64,2}
    "Optional attenuation parameters"
    attenuation :: Bool
    Qμ :: Array{Float64,2}
    Qκ :: Array{Float64,2}
    "Reference frequency in Hz"
    fref :: Float64
end

"""
    PREMPolyModel(; kwargs...) -> m

Construct a `PREMPolyModel` by specifying keyword arguments.  The length of
arrays defines the number of layers and the presence of arguments determines
whether the model has density, anisotropy and attenuation.

Note that radii are given in km, velocities in km/s, densities in g/cm³
and frequency in Hz.

### Keyword arguments
#### Radii
- `r`: The upper radii of each layer in km.  This is the only mandatory argument.
#### Isotropic properties
These must be given each as an `(norder + 1)`×`nlayers` `Matrix` of polynomial
coefficients.  Each parameter (`vp`, `density`, `Qμ`, etc.) may have a
different polynomial order, but all must have the correct number of layers
- `vp`: P-wave velocity
- `vs`: S-wave velocity
- `density`: Material density
#### Anisotropic parameters
If the model is anisotropic and the radial anisotropy parameters are all
provided, then the Voigt average velocities are constructed automatically
unless `vp` and/or `vs` are also given separately.
- `vpv`: Vertical P-wave velocity
- `vph`: Horizontal P-wave velocity
- `vsv`: Vertical S-wave velocity
- `vsh`: Horizontally-polarised S-wave velocity
- `eta`: The ratio C₁₃/(Vph - 2Vsv), where cᵢⱼ is the Voigt elasticity
  matrix, the 1-direction is horizontal and the 3-direction is radial
#### Attenuation parameters
ASCII versions (`Qmu`, `Qkappa`) of these exist for ease of use as well as
the Unicode versions (`Qμ`, `Qκ`).  If both are supplied, an error is thrown.
- `Qμ` or `Qmu`: Shear modulus quality factor
- `Qκ` or `Qkappa`: Bulk modulus quality factor
- `fref`: Reference Q frequency
"""
function PREMPolyModel(; r, vp=nothing, vs=nothing, density=nothing,
                       vph=nothing, vpv=nothing, vsh=nothing, vsv=nothing,
                       eta=nothing, Qμ=nothing, Qmu=nothing,
                       Qκ=nothing, Qkappa=nothing, fref=nothing)
    n = length(r)
    all(x->x>=0, diff(r)) || throw(ArgumentError("radii must increase monotonically"))
    R = last(r)
    aniso = attenuation = false
    # Anisotropy
    if all(!isnothing, (vph, vpv, vsh, vsv, eta))
        aniso = true
        # Set isotropic velocities automatically if not provided
        if any(isnothing, (vp, vs))
            voigt_vels = voigt_velocities.(vpv, vsv, vph, vsh, eta)
            vp === nothing && (vp = first.(voigt_vels))
            vs === nothing && (vs = last.(voigt_vels))
        end
    end
    # Attenuation
    # Default to ASCII versions if given
    Qmu !== nothing && (Qμ = Qmu)
    Qkappa !== nothing && (Qκ = Qkappa)
    if all(!isnothing, (Qμ, Qκ))
        attenuation = true
    end
    # Convert to correct shape, and check for length
    vp = _set_coeffs_size(PREMPolyModel, n, vp)
    vs = _set_coeffs_size(PREMPolyModel, n, vs)
    density = _set_coeffs_size(PREMPolyModel, n, density)
    vph = _set_coeffs_size(PREMPolyModel, n, vph)
    vsh = _set_coeffs_size(PREMPolyModel, n, vsh)
    vpv = _set_coeffs_size(PREMPolyModel, n, vpv)
    vsv = _set_coeffs_size(PREMPolyModel, n, vsv)
    eta = _set_coeffs_size(PREMPolyModel, n, eta)
    Qμ = _set_coeffs_size(PREMPolyModel, n, Qμ)
    Qκ = _set_coeffs_size(PREMPolyModel, n, Qκ)
    fref === nothing && (fref = NaN)
    PREMPolyModel(R, n, r, vp, vs, density, aniso, vph, vpv, vsh, vsv, eta,
        attenuation, Qμ, Qκ, fref)
end

# Check dimensions of input and resize if needed
_set_coeffs_size(::Type{PREMPolyModel}, n, v::Nothing) = Array{Float64}(undef, 0, 0)
function _set_coeffs_size(::Type{PREMPolyModel}, n, v::AbstractVector)
    length(v) == n && return permutedims(v)
    throw(ArgumentError("vector of coefficients must have length $n"))
end
function _set_coeffs_size(::Type{PREMPolyModel}, n, v::AbstractArray{A,2} where A)
    size(v, 2) == n ||
        throw(ArgumentError("second dimension of coefficients array must be $n"))
    v
end
_set_coeffs_size(::Type{PREMPolyModel}, n, v::A) where A =
    throw(ArgumentError("unexpected type of coefficients: $A"))

"""
    _evalpoly(x, coeffs, i)

Evaluate a polynomial whose coefficients are given in `coeffs[:,i]`
at `x`, using Horner's method.
Coefficients should be orderd `a₀, a₁, ..., aₙ`.
"""
function _evalpoly(x, coeffs::AbstractArray{<:Real,2}, i)
    N = size(coeffs, 1)
    val = coeffs[N,i]
    for k in N-1:-1:1
        @inbounds val = muladd(x, val, coeffs[k,i])
    end
    val
end

# Evaluation routines--all documented here
for (sym, name, unit) in zip(model_variables_SeisModel1D, model_names_SeisModel1D, model_units_SeisModel1D)
    # Velocities, which can be corrected for period if we have attenuation
    # and a reference frequency
    if sym in (:vp, :vs, :vph, :vpv, :vsh, :vsv)
        VP, VS = if sym in (:vp, :vs)
            :vp, :vs
        elseif sym in (:vph, :vsh)
            :vph, :vsh
        else
            :vpv, :vsv
        end
        _correct = Symbol(:_correct_attenuation_, String(sym)[2])
        @eval begin
            """
                $(split(string($sym), ".")[end])(m::SeisModel1D, r; depth=false, freq=reffrequency(m)) -> $($name)

            Return the value of $($name)$($unit) for model `m` at radius `r` km.

            If `depth` is `true`, then `r` is given as a depth in km instead.

            If `freq` is given, then the velocity is corrected for attenuation.
            This requires the model has attenuation and a reference frequency.
            """
            function ($sym)(m::PREMPolyModel, r::Real; depth::Bool=false, freq=nothing)
                length(m.$sym) > 0 || throw(ArgumentError("$($name) not defined for model"))
                depth && (r = radius(m, r))
                ir = findlayer(m, r)
                x = r/m.a
                val = _evalpoly(x, m.$(sym), ir)
                freq === nothing && return val
                hasattenuation(m) ||
                    throw(ArgumentError("cannot correct a nonattenuating model " *
                                        "for attenuation"))
                hasreffrequency(m) ||
                    throw(ArgumentError("no reference frequency defined for model"))
                freq == reffrequency(m) && return val
                _vp = $(sym == VP ? :val : :(_evalpoly(x, m.$(VP), ir)))
                _vs = $(sym == VS ? :val : :(_evalpoly(x, m.$(VS), ir)))
                E = 4/3*(_vs/_vp)^2
                qκ = 1/_evalpoly(x, m.Qκ, ir)
                qμ = 1/_evalpoly(x, m.Qμ, ir)
                $(_correct)(val, freq, reffrequency(m), E, qμ, qκ)
            end
        end
    else
        @eval begin
            """
                $(split(string($sym), ".")[end])(m::SeisModel1D, r; depth=false) -> $($name)

            Return the value of $($name)$($unit) for model `m` at radius `r` km.

            If `depth` is `true`, then `r` is given as a depth in km instead.
            """
            function ($sym)(m::PREMPolyModel, r::Real; depth::Bool=false)
                length(m.$sym) > 0 || throw(ArgumentError("$($name) not defined for model"))
                depth && (r = radius(m, r))
                ir = findlayer(m, r)
                x = r/m.a
                val = _evalpoly(x, m.$(sym), ir)
                val
            end
        end
    end
end
# Alternative names
const ρ = density
const Qmu = Qμ
const Qkappa = Qκ

"""
    evaluate(m::SeisModel1D, field::Symbol, r; depth=false) -> vals

Evaluate the model `m` at radius `r` km for the different property/ies in `field`,
returning a scalar for scalar input, and an array for array input.

If `depth` is `true`, `r` is treated as a depth in km instead.
"""
function evaluate(m::PREMPolyModel, field::Symbol, r; depth::Bool=false)
    y = getfield(m, field)
    length(y) > 0 || throw(ArgumentError("'$field' not defined for model"))
    depth && (r = radius(m, r))
    ir = findlayer(m, r)
    x = r/m.a
    val = _evalpoly(x, y, ir)
    val
end

_correct_attenuation_s(vs, freq, reffreq, E, qμ, qκ) = vs*(1 - log(reffreq/freq)/π*qμ)
_correct_attenuation_p(vp, freq, reffreq, E, qμ, qκ) = vp*(1 - log(reffreq/freq)/π*((1 - E)*qκ + E*qμ))

function mass(m::PREMPolyModel, r; depth::Bool=false)
    hasdensity(m) || throw(ArgumentError("model does not contain density"))
    depth && (r = radius(m, r))
    l = findlayer(m, r)
    r *= 1.e3 # SI
    M = 0.
    for i = 1:l
        if i == 1
            R0 = 0.
        else
            R0 = m.r[i-1]*1.e3
        end
        if i < l
            R = m.r[i]*1.e3
        else
            R = r
        end
        for k in 1:size(m.density, 1)
            rho = m.density[k,i]/m.a^(k-1)*1.e3^(2-k)/(k+2)
            M += rho*(R^(k+2) - R0^(k+2))
        end
    end
    4*π*M
end


text:
using RealSense

include("rs-utils.jl")

err = Ref{Ptr{rs2_error}}(C_NULL)

ctx = rs2_create_context(RS2_API_VERSION, err)
checkerror(err)

device = get_a_realsense_device(ctx)
print_device_information(device)

sensor = get_a_sensor_from_a_device(device)

profile = choose_a_streaming_profile(sensor)

is_extendable = rs2_stream_profile_is(profile, RS2_EXTENSION_VIDEO_PROFILE, err)
checkerror(err)
if is_extendable > 0
    intrinsicsRef = Ref{rs2_intrinsics}(rs2_intrinsics(0,0,0,0,0,0,0,(-1,-1,-1,-1,-1)))
    rs2_get_video_stream_intrinsics(profile, intrinsicsRef, err)
    checkerror(err)
    intrinsics = intrinsicsRef[]
    @info "Principal Point: $(intrinsics.ppx), $(intrinsics.ppy)"
    @info "Focal Length: $(intrinsics.fx), $(intrinsics.fy)"
    @info "Distortion Model: $(intrinsics.model)"
    @info "Distortion Coefficients: " intrinsics.coeffs
end

is_extendable = rs2_stream_profile_is(profile, RS2_EXTENSION_MOTION_PROFILE, err)
checkerror(err)
if is_extendable > 0
    intrinsicsRef = Ref{rs2_motion_device_intrinsic}(rs2_motion_device_intrinsic(ntuple(x->ntuple(x->0,Val{4}),Val{3}), ntuple(x->0,Val{3}), ntuple(x->0,Val{3})))
    rs2_get_motion_intrinsics(profile, intrinsicsRef, err)
    checkerror(err)
    intrinsics = intrinsicsRef[]
    println(" Scale X      cross axis      cross axis  Bias X")
    println(" cross axis    Scale Y        cross axis  Bias Y")
    println(" cross axis    cross axis     Scale Z     Bias Z")
    for i in intrinsics.data
        for j in i
            print("  $j    ")
        end
        print("\n")
    end
    println("Variance of noise for X, Y, Z axis:")
    @info intrinsics.noise_variances
    println("Variance of bias for X, Y, Z axis:")
    @info intrinsics.bias_variances
end

rs2_delete_sensor(sensor)
rs2_delete_device(device)
rs2_delete_context(ctx)


text:
# Reads the conllu formatted file w/o xpostag
function load_conllu(file,v::Vocab)
    corpus = Any[]
    s = Sentence(v)
    for line in eachline(file)
        if line == ""
            push!(corpus, s)
            s = Sentence(v)
        elseif (m = match(r"^\d+\t(.+?)\t.+?\t(.+?)\t.+?\t.+?\t(.+?)\t(.+?)(:.+)?\t", line)) != nothing # modify that to use different columns
            #                id   word   lem  upos   xpos feat head   deprel

            word = m.captures[1]
            push!(s.word, word)

            postag = get(v.postags, m.captures[2], 0)
            if postag==0
                Base.warn_once("Unknown postags")
            end
            push!(s.postag, postag)

            head = tryparse(Position, m.captures[3])
            head = isnull(head) ? -1 : head.value
            if head==-1
                Base.warn_once("Unknown heads")
            end
            push!(s.head, head)

            deprel = get(v.deprels, m.captures[4], 0)
            if deprel==0
                Base.warn_once("Unknown deprels")
            end
            push!(s.deprel, deprel)
        end
    end
    return corpus
end

# add xpos-tagged version
function load_conllu2(file,v::Vocab)
    corpus = Any[]; xpos = Dict{String, Int}();
    s = Sentence2(v)
    for line in eachline(file)
        if line == ""
            push!(corpus, s)
            s = Sentence2(v)
        elseif (m = match(r"^\d+\t(.+?)\t.+?\t(.+?)\t(.+?)\t.+?\t(.+?)\t(.+?)(:.+)?\t", line)) != nothing # modify that to use different columns
            #                id   word   lem  upos   xpos feat head   deprel

            word = m.captures[1]
            push!(s.word, word)

            postag = get(v.postags, m.captures[2], 0)
            if postag==0
                Base.warn_once("Unknown postags")
            end
            push!(s.postag, postag)

            xpostag = get!(xpos, m.captures[3], 1+length(xpos))
            push!(s.xpostag, xpostag)

            head = tryparse(Position, m.captures[4])
            head = isnull(head) ? -1 : head.value
            if head==-1
                Base.warn_once("Unknown heads")
            end
            push!(s.head, head)

            deprel = get(v.deprels, m.captures[5], 0)
            if deprel==0
                Base.warn_once("Unknown deprels")
            end
            push!(s.deprel, deprel)
        end
    end
    (length(xpos) == 1) && Base.warn("No xpostag $file")
    return corpus, xpos
end

# Takes the features of a single word, and parses them into 1-by-1 features
# Modifies fdict returns the (features, val) pair, i.e., ("Case", 2), ("Person", 1)
function parse_feats!(fdict::Dict{String, Vector{String}}, feats; testmode=false)
    getfun = (testmode ? get : get!) # modify fdict in train

    res = []
    if  feats == "_" # underscore for no feats
        return res
    end

    farray = split(feats, "|")
    for fs in farray
        fkey, fval = split(fs, "=")
        falls = getfun(fdict, fkey, Vector{String}())
        idx = findfirst(x->x==fval, falls)
        if idx == 0 && !testmode
            push!(falls, fval)
            idx = length(falls)
        end
        if testmode && idx == 0 # Do not count test features
            Base.warn_once("Unkown feats in test set")
            continue
        end
        push!(res, [fkey, idx])
    end
    return res
end

# first version is buggy - (k, v) in Assoc. is not safe
# Shift morphological features to right columns for embedding mtrx implementation
function shift_cfeats2!(corpus, fnav)
    f(x) = (x[2] += fnav[x[1]]-1) # helper to shift features
    for sentence in corpus
        for fgivens in sentence.feats; map(x->f(x), fgivens);end;
        modify_feats!(sentence) # test them
    end
end

# Caches all the feats as an int array
function modify_feats!(s3::Sentence3)
    fs=[]
    for feat in s3.feats
        fset = Int[]
        for f in feat
            push!(fset, f[2])
        end
        push!(fs, fset)
    end
    s3.feats = fs
end

# numof columns in feature embedding matrix
featdim(f::Dict{String,Array{String,1}}) = mapreduce(length, +, 0, values(f))

# xpos-tag and feats version
function load_conllu3(file, v::Vocab,
                      fdict=nothing,
                      xpos=nothing, testmode=true
                      )

    corpus = Any[];
    if fdict == nothing # create train feats dictionary
        fdict = Dict{String, Vector{String}}()
        testmode = false
    end

    getfun = get # not to modify xpos dict in devset
    if xpos == nothing # use train-file's xpos dict
        xpos = Dict{String, Int}();
        getfun = get!
    end
    s = Sentence3(v)
    for line in eachline(file)
        if line == ""
            push!(corpus, s)
            s = Sentence3(v)
        elseif (m = match(r"^\d+\t(.+?)\t.+?\t(.+?)\t(.+?)\t(.+?)\t(.+?)\t(.+?)(:.+)?\t", line)) != nothing # modify that to use different columns
            #                id   word   lem  upos   xpos  feat   head   deprel

            word = m.captures[1]
            push!(s.word, word)

            postag = get(v.postags, m.captures[2], 0)
            if postag==0
                Base.warn_once("Unknown postags")
            end
            push!(s.postag, postag)

            xpostag = getfun(xpos, m.captures[3], 1+length(xpos))
            push!(s.xpostag, xpostag)

            wordfeats = m.captures[4]
            feats = parse_feats!(fdict, wordfeats, testmode=testmode)
            push!(s.feats, feats)

            head = tryparse(Position, m.captures[5])
            head = isnull(head) ? -1 : head.value
            if head==-1
                Base.warn_once("Unknown heads")
            end
            push!(s.head, head)

            deprel = get(v.deprels, m.captures[6], 0)
            if deprel==0
                Base.warn_once("Unknown deprels")
            end
            push!(s.deprel, deprel)
        end
    end
    (length(xpos) == 1) && Base.warn("No xpostag $file")
    return corpus, xpos, fdict
end

# To create vocabulary from pre-trained lstm model, modify that to use different cols
function create_vocab(d)
    Vocab(d["char_vocab"],
          Dict{String, Int}(),
          d["word_vocab"],
          d["sosword"],
          d["eosword"],
          d["unkword"],
          d["sowchar"],
          d["eowchar"],
          d["unkchar"],
          get(d, "postags", UPOSTAG),
          get(d, "deprels", UDEPREL)
          )
end

# After filling word, forward and backward vectors cache the concated version
function cachelmvec!(corpus)
    for sent in corpus
        for i in 1:length(sent)
            push!(sent.cavec, vcat(sent.wvec[i], sent.fvec[i], sent.bvec[i]))
        end
    end
end

function minibatch(corpus, batchsize; maxlen=typemax(Int), minlen=1, shuf=false)
    data = Any[]
    sorted = sort(corpus, by=length)
    i1 = findfirst(x->(length(x) >= minlen), sorted)
    if i1==0; error("No sentences >= $minlen"); end
    i2 = findlast(x->(length(x) <= maxlen), sorted)
    if i2==0; error("No sentences <= $maxlen"); end
    for i in i1:batchsize:i2
        j = min(i2, i+batchsize-1)
        push!(data, sorted[i:j])
    end
    if shuf
        data=shuffle(data)
    end
    return data
end

createfnav(fdict)=(fnav = Dict{String, Int}();counter = 1;for (k,v) in fdict; fnav[k]=counter; counter+=length(v);end; fnav)

# Dead code
# # Shift morphological features to right columns for embedding mtrx implementation
# function shift_cfeats!(corpus, fdict)
#     fnav = Dict{String, Int}() # Shifting values
#     counter = 1;for (k,v) in fdict; fnav[k]=counter; counter+=length(v);end;
#     f(x) = (x[2] += fnav[x[1]]-1) # helper to shift features
#     for sentence in corpus
#         for fgivens in sentence.feats; map(x->f(x), fgivens);end;
#         modify_feats!(sentence) # test them
#     end
# end