text:
# tensor_deriv.jl - tensor derivative utils (using Einstein notation)

const TDIFF_PHS = [:A, :B, :C, :X, :Y, :V, :W, :Z,
                   :i, :j, :k, :m, :n, :p, :q, :r, :s, :t, :l]

const TDIFF_VAR_NAMES = [:V, :W, :X, :Y]

## TensorDerivExpr - an atom of TensorDeriv

"""Helper struct to simplify calculations on tensor derivative expressions"""
mutable struct TensorDeriv
    var::Union{Symbol,Expr}
    wrt::Union{Symbol,Expr}
    ex::Any
    guards::Vector{Expr}
end

"""
Create a TensorDeriv from an expression. E.g. in:

     dy!dx[i,j] = W[i,j]

`x` and `y` will be looked up in `g` to find their sizes and make a derivative, e.g.:

    dy[i] / dx[j] = W[i,j]
"""
function TensorDeriv(g::EinGraph, dex::Expr; guards=nothing)
    full_vname, idxs = split_indexed(dex.args[1])
    vname, wrtname_ = split_deriv_name(full_vname)
    var_idx_len = g[undname(vname)].val |> size |> length
    vidxs, wrtidxs_ = isempty(idxs) ? ([],  []) : (idxs[1:var_idx_len], idxs[var_idx_len+1:end])
    var, wrt = make_indexed(vname, vidxs), make_indexed(wrtname_, wrtidxs_)
    ex = without_guards(dex.args[2])
    guards = guards != nothing ? guards : find_guards(dex.args[2])
    return TensorDeriv(var, wrt, ex, guards)
end

"""
Create a TensorDeriv from a pretty-printed expression. E.g.:

     dy[i] / dx[j] = W[i,j]

`y[i]` will be parsed as var and `z[j]` as w.r.t. var.
"""
function TensorDeriv(dex::Expr; guards=nothing)
    var, wrt = dex.args[1].args[2:3]
    ex = without_guards(dex.args[2])
    guards = guards != nothing ? guards : find_guards(dex.args[2])
    return TensorDeriv(var, wrt, ex, guards)
end

function Base.copy(td::TensorDeriv; var=getvar(td), wrt=getwrtvar(td),
                   ex=getexpr(td), guards=getguards(td))
    return TensorDeriv(var, wrt, ex, guards)
end

Espresso.getvar(td::TensorDeriv) = td.var
Espresso.varname(td::TensorDeriv) = split_indexed(getvar(td))[1]
Espresso.varidxs(td::TensorDeriv) = split_indexed(getvar(td))[2]

getwrtvar(td::TensorDeriv) = td.wrt
wrtname(td::TensorDeriv) = split_indexed(getwrtvar(td))[1]
wrtidxs(td::TensorDeriv) = split_indexed(getwrtvar(td))[2]

lhs_indices(td::TensorDeriv) = vcat(varidxs(td), wrtidxs(td))
rhs_indices(td::TensorDeriv) = convert(Vector{Any}, unique(flatten(find_indices(getexpr(td)))))
free_indices(td::TensorDeriv) = convert(Vector{Any}, setdiff(rhs_indices(td), lhs_indices(td)))

Espresso.getexpr(td::TensorDeriv) = td.ex
Espresso.getguards(td::TensorDeriv) = td.guards

# # var_indices(td::TensorDeriv) = convert(Vector{Any}, td.dvar.args[2:end])
# # wrt_indices(td::TensorDeriv) = convert(Vector{Any}, td.wrt.args[2:end])
# # deriv_indices(td::TensorDeriv) = vcat(var_indices(td), wrt_indices(td))
# # all_indices(td::TensorDeriv) = union(deriv_indices(td),
# #                                     flatten(Any, get_indices(expr(td))))

function all_indices(td::TensorDeriv)
    return union(varidxs(td),
                 wrtidxs(td),
                 flatten(Any, get_indices(getexpr(td))))
end

function single_var(td::TensorDeriv)
    new_name = Symbol("$(varname(td))!$(wrtname(td))")
    new_idxs = vcat(varidxs(td), wrtidxs(td))
    return make_indexed(new_name, new_idxs)
end

function to_expr(td::TensorDeriv)
    lhs = single_var(td)
    lhs_idxs = get_indices(lhs)
    guards = getguards(td)
    rhs = getexpr(td)
    rhs_idxs = get_indices(rhs)
    # TODO: check if it works for blocks too
    return Espresso.apply_guards(:($lhs = $rhs),  guards)
end

function to_expr_pp(td::TensorDeriv)
    lhs = :($(getvar(td)) / $(getwrtvar(td)))
    grds = getguards(td)
    rhs = !isempty(grds) > 0 ? Expr(:call, :*, getexpr(td), grds...) : getexpr(td)
    return :($lhs = $rhs)
end

function Base.show(io::IO, td::TensorDeriv)
    print(io, to_expr_pp(td))
end

# TODO:  moved to Esresso, import from there

"""
Given a set of existing indices and current position of iterator,
find the next index not in the set.
"""
function next_index(existing::Set{T}, pos::Int) where T
    while pos <= length(IDX_NAMES) && in(IDX_NAMES[pos], existing)
        pos += 1
    end
    if pos <= length(IDX_NAMES)
        return IDX_NAMES[pos], pos + 1
    else
        throw(BoundsError("IDX_NAMES"))
    end
end

function next_indices(existing::Set{T}, pos::Int, count::Int) where T
    new_indices = Array{Symbol}(0)
    for i=1:count
        new_idx, pos = next_index(existing, pos)
        push!(new_indices, new_idx)
    end
    return new_indices
end

"""
Given a set of existing indicies and possible duplicates, find for each duplicate
a replacement - index from IDX_NAMES that is not used yet.
"""
function index_replacements(existing::Set{T}, maybedups::Vector{T}) where T
    repls = Dict{Symbol,Symbol}()
    pos = 1
    for idx in maybedups
        # maybedups should also be included in existing index set
        all_existing = union(existing, Set(maybedups), Set(keys(repls)))
        if in(idx, existing) && !in(idx, keys(repls))
            repls[idx], pos = next_index(all_existing, pos)
        end
    end
    return repls
end

function reindex_with_guards(full_ex, guards::Vector{Expr})
    @assert full_ex.head == :(=) && full_ex.args[1].args[1] == :/
    def, ex = full_ex.args
    anchors = Set(get_indices(def) |> flatten)
    pairs = Tuple{Any,Any}[(grd.args[2], grd.args[3]) for grd in guards]
    st, new_pairs = reduce_equalities(pairs, anchors)
    new_guards = [:($i1 == $i2) for (i1, i2) in new_pairs]
    new_ex = subs(ex, st)
    return :($def = $new_ex), new_guards
end

"""
Reindex second tensor derivative (dydx) so that:

    * dydx's var indices match dzdy's w.r.t. indices
    * no other indices in dydx equal any indices in dzdy
"""
function reindex_to_match(dzdy::TensorDeriv, dydx::TensorDeriv)
    common_idxs_st = Dict(zip(varidxs(dydx), wrtidxs(dzdy)))
    other_idxs_st = index_replacements(Set(all_indices(dzdy)), all_indices(dydx))
    st = merge(other_idxs_st, common_idxs_st)
    new_dydx_ex = subs(to_expr_pp(dydx), st) |> sanitize
    new_dydx = TensorDeriv(new_dydx_ex)
    return dzdy, new_dydx
end

"""
Reindex sibling derivatives so that:

    * td1's var indices match td2's w.r.t. indices
    * no other indices in td1 equal any indices in td2
"""
function reindex_siblings_to_match(td1::TensorDeriv, td2::TensorDeriv)
    common_idxs_st = Dict(zip(lhs_indices(td2), lhs_indices(td1)))
    # replace free indices from td2 with anything not in td1
    all_idxs1 = Set{Any}(all_indices(td1))
    free_idxs2_st = index_replacements(all_idxs1, collect(free_indices(td2)))
    st = merge(free_idxs2_st, common_idxs_st)
    new_td2_ex = subs(to_expr_pp(td2), st) |> sanitize
    new_td2 = TensorDeriv(new_td2_ex)
    return td1, new_td2
end

# """
# Elementwise multuplication of tensor derivatives.
# Example:

#     dzdx = dzdy ⊗ dydx

# which may expand to:

#     dz/dy[i] = v[i]
#     dy[i]/dx[j] = w[i,j]
#     dz/dx[j] = v[i] .* w[i,j]
# """
function ⊗(dzdy::TensorDeriv, dydx::TensorDeriv)
    # can only multiply related derivatives, e.g. dz/dy * dy/dx
    @assert wrtname(dzdy) == varname(dydx)
    dzdy, dydx = reindex_to_match(dzdy, dydx)
    guards = vcat(getguards(dzdy), getguards(dydx))
    new_td_lhs = :($(getvar(dzdy)) / $(getwrtvar(dydx)))
    new_td_rhs = :($(single_var(dzdy)) .* $(getexpr(dydx))) |> simplify
    new_td_ex_no_guards = :($new_td_lhs = $new_td_rhs)
    new_td_ex, new_guards = reindex_with_guards(new_td_ex_no_guards, guards)
    new_td = TensorDeriv(new_td_ex; guards=new_guards)
    dzdx = TensorDeriv(new_td)
    return dzdx
end

# presumably, we need to use extend_deriv instead
function ⊕(td1::TensorDeriv, td2::TensorDeriv)
    @assert varname(td1) == varname(td2)
    @assert wrtname(td1) == wrtname(td2)
    dvar_idxs_st = Dict(zip(varidxs(td2), varidxs(td1)))
    wrt_idxs_st = Dict(zip(wrtidxs(td2), wrtidxs(td1)))
    st = merge(dvar_idxs_st, wrt_idxs_st)
    free_idxs = free_indices(td2)
    # TODO: should we also inclue all indicies of expr(td1)?
    all_existing_idxs = Set{Symbol}(vcat(keys(st)..., values(st)..., free_idxs))
    # TODO: use `index_replacements()` instead
    next_idx_pos = 1
    for idx in free_idxs
        if in(idx, values(st))
            st[idx], next_idx_pos = next_index(all_existing_idxs, next_idx_pos)
        end
    end
    wrt2_reindexed = subs(getwrtvar(td2), st)
    ex2_reindexed = subs(getexpr(td2), st)
    guards2_reindexed = Expr[subs(g, st) for g in find_guards(td2)]
    new_ex = simplify(getexpr(td1) ⊕ ex2_reindexed)
    new_guards = vcat(td1.guards, guards2_reindexed)
    new_td = TensorDeriv(getvar(td1), wrt2_reindexed, new_ex, new_guards)
    return reindex_with_guards(new_td)
end

## tensor differentiation rules

struct TensorDiffRule <: AbstractDiffRule
    pat::Expr             # pattern of expression to differentiate
    deriv::TensorDeriv    # pattern of differentiation expression
end

function Base.show(io::IO, rule::TensorDiffRule)
    print(io, "TensorDiffRule($(rule.pat) ==> $(to_expr(rule.deriv)))")
end

"""
Convert scalar diff rule to a tensor diff rule.

 * ew_rule   - elementwise (scalar) rule
 * orig_idxs - indices of full tensor expression,
               e.g. for `z[i] = X[i,j] * y[j]` it's [[:i], [:i, :j], [:j]]
 * idx       - index of input parameter to differentiate w.r.t. it
"""
function to_tensor_rule(ew_rule::DiffRule, orig_idxs::Vector{Vector{T}}, idx::Int) where T
    ew_pat = ew_rule.pat
    op = ew_pat.args[1]
    ew_ex = ew_rule.dpat
    # tensor var names and indices
    tvar_names = TDIFF_VAR_NAMES[1:length(ew_pat.args)-1]
    # tvar_idxs = IDX_NAMES[1:length(orig_idxs[1])]
    tvars = [make_indexed(name, IX) for (name, IX) in zip(tvar_names, orig_idxs[2:end])]
    tvar_idxs = IDX_NAMES[1:length(orig_idxs[1])]
    # tvars = [Expr(:ref, tvar, tvar_idxs...) for tvar in tvar_names]
    # dvar variable
    var_name = :Z
    dvar_name = dname(var_name)
    var = make_indexed(var_name, orig_idxs[1])
    dvar = make_indexed(dvar_name, orig_idxs[1])
    # w.r.t. variable
    wrt_name = tvar_names[idx]
    dwrt_name = dname(wrt_name)
    wrt_idxs = next_indices(Set(flatten(orig_idxs)), 1, length(orig_idxs[idx + 1]))
    dwrt = make_indexed(dwrt_name, wrt_idxs)
    # new pattern
    tpat = Expr(:call, op, tvars...)
    full_tpat = :($var = $tpat)
    # new tensor derivative expression
    tex = rewrite(tpat, ew_pat, ew_ex; phs=DIFF_PHS)
    # constructing tensor derivative
    if length(orig_idxs[idx + 1]) > 0
        # old and new w.r.t. indices must be equal
        tguards = Expr[:($i1 == $i2) for (i1, i2) in zip(orig_idxs[idx+1], wrt_idxs)]
    else
        tguards = Expr[]  # TODO: this should be covered by previous definition too
    end
    tderiv = TensorDeriv(dvar, dwrt, tex, tguards)
    return TensorDiffRule(full_tpat, tderiv)
end

const TENSOR_DIFF_RULES = Dict{Tuple{OpName, Int}, Vector{TensorDiffRule}}()

function push_tensordiff!(op::OpName, deriv_idx::Int, rule::TensorDiffRule)
    if !haskey(TENSOR_DIFF_RULES, (op, deriv_idx))
        TENSOR_DIFF_RULES[(op, deriv_idx)] = TensorDiffRule[]
    end
    push!(TENSOR_DIFF_RULES[(op, deriv_idx)], rule)
end

function _tensordiff(ex, dex)
    ex = sanitize(ex)
    dex = sanitize(dex)
    op = canonical(current_module(), ex.args[2].args[1])
    idxs = get_indices(ex.args[2])
    dvar = dex.args[1].args[2]
    wrt = dex.args[1].args[3]
    deriv_ex = without_guards(sanitize(dex.args[2]))
    grds = find_guards(dex)
    deriv = TensorDeriv(dvar, wrt, deriv_ex, grds)
    diff_var_name = Symbol(string(wrt.args[1])[2:end])
    var_names = [iex.args[1] for iex in ex.args[2].args[2:end]]
    deriv_idx = find(var_names .== diff_var_name)[1]
    rule = TensorDiffRule(ex, deriv)
    push_tensordiff!(op, deriv_idx, rule)
end

macro tensordiff(ex, dex)
    _tensordiff(ex, dex)
    nothing
end

function tfind_rule(fullex::Expr, idx::Int)
    @assert fullex.head == :(=) && fullex.args[2].head == :call
    op = fullex.args[2].args[1]
    haskey(TENSOR_DIFF_RULES, (op, idx)) || return Nullable{TensorDiffRule}()
    rules = TENSOR_DIFF_RULES[(op, idx)]
    matches = pat -> !isnull(matchex(pat, fullex; phs=TDIFF_PHS, allow_ex=false))
    matching = findfirst(matches,
                         [r.pat for r in rules])
    matching != 0 || return Nullable{TensorDiffRule}()
    return Nullable{TensorDiffRule}(rules[matching])
end

dname(var::Symbol) = Symbol("d$var")
undname(dvar::Symbol) = Symbol(string(dvar)[2:end])

"""dZ[i]/dX[j] = ... ==> Z[i]/X[i] = ..."""
function unpack_deriv(ex::Expr)
    @assert ex.head == :(=)
    @assert ex.args[1].head == :call && ex.args[1].args[1] == :/
    dvar, dwrt = [split_indexed(dv)[1] for dv in ex.args[1].args[2:3]]
    var, wrt = undname(dvar), undname(dwrt)
    return subs(ex, Dict(dvar => var, dwrt => wrt))
end

"""Z[i]/X[j] = ... ==> dZ[i]/dX[i] = ..."""
function pack_deriv(ex::Expr)
    @assert ex.head == :(=)
    @assert ex.args[1].head == :call && ex.args[1].args[1] == :/
    var, wrt = [split_indexed(v)[1] for v in ex.args[1].args[2:3]]
    dvar, dwrt = dname(var), dname(wrt)
    lhs = subs(ex.args[1], Dict(var => dvar, wrt => dwrt))
    rhs = ex.args[2]
    return Expr(:(=), lhs, rhs)
end

"""
Derivative of a primitive expression in Einstein notation.
"""
function tderivative(fullex::Expr, idx::Int)
    maybe_rule = tfind_rule(fullex, idx)
    if !isnull(maybe_rule)
        rule = get(maybe_rule)
        unpacked_rule_dex = unpack_deriv(to_expr_pp(rule.deriv))
        unpacked_dex = rewrite(fullex, rule.pat, unpacked_rule_dex; phs=TDIFF_PHS)
        dex = pack_deriv(unpacked_dex)
        return TensorDeriv(dex)
    else
        idxs = get_indices(fullex)
        # elementwise or broadcasting function
        op = opname(current_module(), fullex.args[2].args[1])
        types = [Float64 for i=1:length(fullex.args[2].args)-1]
        ew_maybe_rule = find_rule(op, types, idx)
        ew_rule = (!isnull(ew_maybe_rule) ? get(ew_maybe_rule) :
                   register_rule(op, types, idx))
        trule = to_tensor_rule(ew_rule, idxs, idx)
        push_tensordiff!(op, idx, trule)
        # now rule is registered, recursively call itself
        return tderivative(fullex, idx)
    end
end

function tderivative(fullex::Expr, var::Symbol)
    @assert fullex.head == :(=)
    @assert fullex.args[2].head == :call
    ivars = [var for var in fullex.args[2].args[2:end]]
    vars = [isa(ivar, Expr) ? ivar.args[1] : ivar for ivar in ivars]
    matching = findfirst(vars .== var)
    matching != 0 || error("Variable `$var` isn't present " *
                           "in expression `$fullex`")
    return tderivative(fullex, matching[1])
end


text:
Libdl.dlopen(:libmkl_rt) # necessary to preload - don't know why
const libmkl   = :libmkl_rt
const libvml   = :libmkl_vml_avx

## ?gemv
function mkl_gemv!( layout::Cint, trans::Char,
                    m::Integer, n::Integer,
                    alpha::Float64, A::Union{Ptr{Float64},DenseArray{Float64}}, lda::Integer, X::Union{Ptr{Float64},DenseArray{Float64}}, incx::Integer,
                    beta::(Float64), Y::Union{Ptr{Float64},DenseArray{Float64}}, incy::Integer)

    ccall((:dgemv, libmkl), Void,
            (Ptr{Cint}, Ptr{Int},
            Ptr{Int}, Ptr{Int},
            Ptr{Float64}, Ptr{Float64}, Ptr{Int}, Ptr{Float64}, Ptr{Int},
            Ptr{Float64}, Ptr{Float64}, Ptr{Int}),
            &layout, &trans,
            &m, &n,
            &alpha, A, &lda, X, incx,
            &beta, Y, &incy)
    Y
end

A = [1.0 2.0 3.0 ; 3.0 4.0 5.0]
m,n = size(A,1),size(A,2)

X = [1.0,2.0,3.0]
Y = [1.0,1.0]

# Native Julia
2*A*X+Y

# LinAlg.BLAS.gemv!('N',2.0,A,X,1.0,Y)

convert(Cint,102)

mkl_gemv!(  convert(Cint,102), 'N',
            2, 3,
            2.0, A, 3, X, 1,
            1.0, Y, 1)

ext: jl


text:
using BinaryBuilder
using Pkg

sources = [
	GitSource("https://github.com/jstrube/Minuit2_Julia_Wrapper.git", "10bf344912baa777a940862d6d80631bea34a448")
]

julia_version = v"1.5.3"

name = "Minuit2_Julia_Wrapper"
version = v"0.1"

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/Minuit2_Julia_Wrapper/
mkdir build && cd build
cmake -DJulia_PREFIX=${prefix} -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DCMAKE_FIND_ROOT_PATH=${prefix} ..
VERBOSE=ON cmake --build . --config Release --target install
install_license $WORKSPACE/srcdir/Minuit2_Julia_Wrapper/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")
platforms = libjulia_platforms(julia_version)
# filter!(!Sys.isfreebsd, platforms)
# filter!(!Sys.iswindows, platforms)
# filter!(p -> arch(p) != "armv7l", platforms)

platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libminuit2wrap", :minuit2wrap)
]

# Dependencies that must be installed before this package can be built
dependencies = [
        Dependency(PackageSpec(name="libcxxwrap_julia_jll",version=v"0.8.5")),
        Dependency(PackageSpec(name="Minuit2_jll")),
        BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version))
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"8", julia_compat = "^$(julia_version.major).$(julia_version.minor)")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using DifferentialForms
using DifferentialForms.Forms

using StaticArrays
using Test

@testset "binomial D=$D" for D in 0:10
    for R in 0:D
        @test binomial(Val(D), Val(R)) == binomial(D, R)
    end
end

@testset "bits2uint D=$D" for D in 0:10
    for uint in Unsigned(0):((Unsigned(1) << D) - 1)
        bits = DifferentialForms.Forms.uint2bits(Val(D), uint)
        uint′ = DifferentialForms.Forms.bits2uint(bits)
        @test uint′ == uint
    end
end

@testset "lin2bit D=$D" for D in 0:10
    for R in 0:D
        N = binomial(Val(D), Val(R))
        for lin in 1:N
            bits = DifferentialForms.Forms.lin2bit(Val(D), Val(R), lin)
            @test count(bits) == R
            lin′ = DifferentialForms.Forms.bit2lin(Val(D), Val(R), bits)
            @test lin′ == lin
        end
    end
end

@testset "bit2lst D=$D" for D in 0:10
    for R in 0:D
        N = binomial(Val(D), Val(R))
        for lin in 1:N
            bits = DifferentialForms.Forms.lin2bit(Val(D), Val(R), lin)
            lst = DifferentialForms.Forms.bit2lst(Val(D), Val(R), bits)
            @test length(lst) == R
            @test all(>=(1), lst)
            @test all(<=(D), lst)
            @test issorted(lst)
            @test length(unique(lst)) == R
            bits′ = DifferentialForms.Forms.lst2bit(Val(D), Val(R), lst)
            @test bits′ == bits
        end
    end
end

@testset "lin2lst D=$D" for D in 0:10
    for R in 0:D
        N = binomial(Val(D), Val(R))
        for lin in 1:N
            lst = DifferentialForms.Forms.lin2lst(Val(D), Val(R), lin)
            lin′ = DifferentialForms.Forms.lst2lin(Val(D), Val(R), lst)
            @test lin′ == lin
        end
    end
end


text:
function train!(method::KernelInterpolation)
    X = method.trainingData.Xs
    Fs = method.trainingData.Fs
    k = method.kernel

    N, D = size(X)

    A = ones(N+1, N+1)

    y = zeros(N+1)
    y[1:N] = Fs

    for i = 1:N
        for j = 1:N
            A[i, j] = kernel(k, view(X, i,:), view(X, j,:))
        end

        A[i, i] += method.λ
    end

    A[N+1, N+1] = 0.0
    method.kernelMatrix = A
    method.kernelMatrixInv = inv(A)

    # y = Ab
    # b = A⁻¹ ⋅ y

    b = method.kernelMatrixInv * y

    method.coeffs = b
    method
end

train(method::KernelInterpolation) = train!(method)

function evaluate(x::Real, method::KernelInterpolation)
    if length(method.coeffs) == 0
        @info("Training method...")
        train!(method)
    end

    evaluate(Float64[x], method)
end

function evaluate(x, method::KernelInterpolation)
    if length(method.coeffs) == 0
        @info("Training method...")
        train!(method)
    end

    X = method.trainingData.Xs
    k = method.kernel

    s = 0.0
    for i = 1:size(X,1)
        s += kernel(k,x, view(X, i,:)) * method.coeffs[i]
    end

    s + method.coeffs[end]
end

function evaluate(X::Array{Float64, 2}, method::KernelInterpolation)
    if length(method.coeffs) == 0
        @info("Training method...")
        train!(method)
    end

    Float64[ evaluate(view(X, i,:), method) for i=1:size(X,1) ]
end


text:
function contract(T::BlockSparseTensor,
                  labelsT,
                  C::CombinerTensor,
                  labelsC)
  # Get the label marking the combined index
  # By convention the combined index is the first one
  # TODO: consider storing the location of the combined
  # index in preperation for multiplce combined indices
  cpos_in_labelsC = 1
  clabel = labelsC[cpos_in_labelsC]
  c = combinedindex(C)
  labels_uc = deleteat(labelsC,cpos_in_labelsC)
  if labelsC[1] ∉ labelsT
    # Combine
    labelsRc = contract_labels(labelsC,labelsT)
    cpos_in_labelsRc = findfirst(==(clabel),labelsRc)
    labelsRuc = insertat(labelsRc,labels_uc,cpos_in_labelsRc)
    indsRc = contract_inds(inds(C),labelsC,inds(T),labelsT,labelsRc)
    perm = getperm(labelsRuc,labelsT)
    ucpos_in_labelsT = Tuple(findall(x->x in labels_uc,labelsT))
    Rc = permutedims_combine(T,indsRc,perm,ucpos_in_labelsT,blockperm(C),blockcomb(C))
    return Rc
  else
    # Uncombine
    labelsRc = labelsT
    cpos_in_labelsRc = findfirst(==(clabel),labelsRc)
    # Move combined index to first position
    #if cpos_in_labelsRc != 1
    #  labelsRc_orig = labelsRc
    #  labelsRc = deleteat(labelsRc,cpos_in_labelsRc)
    #  labelsRc = insertafter(labelsRc,clabel,0)
    #  cpos_in_labelsRc = 1
    #  perm = getperm(labelsRc,labelsRc_orig)
    #  T = permutedims(T,perm)
    #  labelsT = permute(labelsT,perm)
    #end
    labelsRuc = insertat(labelsRc,labels_uc,cpos_in_labelsRc)
    indsRuc = contract_inds(inds(C),labelsC,inds(T),labelsT,labelsRuc)
    Ruc = uncombine(T,indsRuc,cpos_in_labelsRc,blockperm(C),blockcomb(C))
    return Ruc
  end
end

contract(C::CombinerTensor,
         labelsC,
         T::BlockSparseTensor,
         labelsT) = contract(T,labelsT,C,labelsC)

# Special case when no indices are combined
contract(T::BlockSparseTensor,
         labelsT,
         C::CombinerTensor{<:Any,0},
         labelsC) = copy(T)


text:
module SparseMatricesCSR

using SparseArrays
using LinearAlgebra
using SuiteSparse

import Base: convert, copy, size, getindex, setindex!, show, count, *, IndexStyle
import LinearAlgebra: mul!, lu, lu!
import SparseArrays: nnz, getnzval, nonzeros, nzrange
import SparseArrays: findnz, rowvals, getnzval, issparse

export SparseMatrixCSR
export SymSparseMatrixCSR
export sparsecsr, symsparsecsr
export colvals, getBi, getoffset

include("SparseMatrixCSR.jl")

include("SymSparseMatrixCSR.jl")

end # module

# DONE sparse
# DONE issparse
# DONE nnz
# DONE nonzeros
# DONE colvals
# DONE nzrange
# DONE findnz
# DONE getrowptr
# DONE getnzval
# DONE getcolval
# TODO spzeros
# TODO spdiagm
# TODO blockdiag
# TODO sprand
# TODO sprandn
# TODO droptol!
# TODO dropzeros
# TODO dropzeros!
# TODO permute
# TODO permute!


text:
using Revise, BenchmarkTools
using DifferentialEquations, Optim, DiffEqFlux, DiffEqSensitivity, Plots
using Zygote, ForwardDiff, ReverseDiff
using LinearAlgebra, RecursiveArrayTools, Measurements, Statistics
using DataFrames
using BlackBoxOptim
# -------------------------------------------------------------------------

function run_model(model,p,data,densities,constants,t) # version of run function with multiple models
    # f = (du,u,p,t) -> model(du,u,10 .^ p,t,densities,constants)
    f = (du,u,p,t) -> model(du,u,p,t,densities,constants)
    tmp_prob = ODEProblem(f,u₀,tspan,p)
    tmp_sol = solve(tmp_prob,AutoVern7(KenCarp4()),saveat=t, abstol=1e-8,reltol=1e-8,verbose=false)
    # tmp_sol = solve(tmp_prob,Vern7(),saveat=t, abstol=1e-8,reltol=1e-8,verbose=false)
    tmp_sol
end

function loss(p,data,model,dens,cons,t,unique_t)
    Δy_all = [] # keep track of errors for all curves (change this to array of arrays!)
    # Δy_all = Zygote.Buffer([[0.,0.],[0.,0.]],length(data))
    for i in 1:length(data)
      # run model with densities and parameter/dens/cons combinations
      tmp_sol = run_model(model,p,data[i],dens[i],cons[i],unique_t[i])

      if tmp_sol.retcode != :Success #if there is an error return infitite loss
        Δy_sol = fill(Inf, length(data[i]))
      else
        Δy_sol = find_diffs(tmp_sol,t[i],unique_t[i],data[i])
      end
      Δy_all = vcat(Δy_all,Δy_sol)
      # Δy_all[i] = Δy_sol
    end
    # Δy = copy(Δy_all
    # Losses = [sum(Δ .^ 2) for Δ in Δy]
    # Losses = [sum( (Δy[i] .^ 2) .* t[i] .^2 ) for i in 1:length(Δy)]
    t_all = reduce(vcat,t)

    # Losses = Δy_all .^2 .* t_all .^2
    Losses = Δy_all .^2 .* t_all
    L = sum(Losses) # calculate error across all differences from all curves

    return L
end

function find_diffs(tmp_sol,t,unique_t,data) # solutions use unique time points,
# need to evaluate the solution at each experimental time point as well as
# calculating the difference between fit and experiment

    Σ_sol = sum(Array(tmp_sol),dims=1)
    # Σ_sol = sum(tmp_sol,dims=1)
    # Δy = Zygote.Buffer(t,length(t))
    Δy = []
    # Δy = zeros(length(t))
    for (i, tᵢ) in enumerate(t)
        ind = isequal.(unique_t,tᵢ)
        yᵢ_exp = Σ_sol[ind][1]
        yᵢ_obs = data[i]
        Δyᵢ = yᵢ_obs .- yᵢ_exp
        # Δy[i] = Δyᵢ
        # push!(Δy,Δyᵢ)
        Δy=vcat(Δy,Δyᵢ)
    end
    return copy(Δy) #Δy[1:end] # copy(Δy)
end

# function callback(p,l,tmp_sol,data,model,dens,cons,t)
function callback(p,l,data,model,dens,cons,t,unique_t)
  # @show l, p
  @show l
  # Σ_sol = run_model(model,p,data,dens,cons,unique_t)

  # fig = plot(unique_t,Σ_sol')
  # scatter!(fig,t,data,legend=:bottomright)
  # display(fig)

  false
end

function scan_param_space(p,data,model,dens,cons,t,unique_t,lb,ub)
      #create rand mesh
      n = 5000
      Δbounds = ub - lb
      # ps = rand(n,size(p,1)) .* 10 .- 8
      ps = rand(n,size(p,1)) .* Δbounds  .+ lb
      # ps = rand(n,size(p,1)) .* 100
      errors = zeros(n)
      for i in 1:n
          # error, a = loss(ps[i,:],data,model,dens,cons,t)
          error = loss(ps[i,:],data,model,dens,cons,t,unique_t)
          errors[i] = error
          # @show error
          if mod(i,100) == 0
            @show i, error
          end
      end
      min_ind = argmin(errors)
      min_error = minimum(errors)
      @show findmin(errors)
      p_best = ps[min_ind,:]
      @show p_best
      if size(p,1) == 2
          fig = scatter(ps[:,1],ps[:,2],log10.(errors),
          zcolor = log10.(errors),camera=[0,90])
          scatter!([p_best[1]],[p_best[2]],[log10(min_error)],markercolor=:green,markersize=5)
          display(fig)
      end
      return p_best
end

function optimization(model,data,dens,cons,options,t)
# fit function - given model, data, densities, options get best parameters
  unique_t = [unique(times) for times in t]
  loss_in = (p) -> loss(p,data,model,dens,cons,t,unique_t)
  callback_in = (p,l) -> callback(p,l,data,model,dens,cons,t,unique_t)

  # @unpack lb,ub = options
  lb,ub = options

  # rates = [0.41751687626687767, -0.23866855433293688, 0.06689410353314687, -0.5510289313898404, -0.09606981254493874, -0.1155499030032302, -3.9999999999999964, -1.7908513448187287]
  # hes_zyg = Zygote.hessian(loss_in,rates)
  @show loss_in(pinit)
  # start = scan_param_space(pinit,data,model,dens,cons,t,unique_t,lb,ub) # !! change behavior so that it doesnt need pinit
  # start = DiffEqFlux.sciml_train(loss_in,start,ADAM(0.001),cb=callback_in,maxiters=50)
  # start = start.minimizer
  # @show start
  # res = optimize(loss_in,start,BFGS(),Optim.Options(f_tol = 1e-8))
  # res = optimize(loss_in,start,NelderMead(),Optim.Options(f_tol = 1e-8))
  res_bbo = bboptimize(loss_in; SearchRange = (lb, ub),
                                NumDimensions = length(pinit),
                                Method = :adaptive_de_rand_1_bin_radiuslimited,
                                # Method = :separable_nes,
                                NThreads=Threads.nthreads()-1,
                                MaxSteps = 10000)#,

                                # # MaxFuncEvals = 10000)#,
                                # TargetFitness=0.05)#,
                                # NThreads=Threads.nthreads()-1)
  p₀_bbo = best_candidate(res_bbo)
  @show p₀_bbo
  # p₀_bbo = pinit
  # bbo_minimizer = [0.829897, -3.39453, -3.83109, 0.215222, 0.848485, -3.08371, 0.606198, -0.906295]
  # res = optimize(loss_in,lb,ub,p₀_bbo, Fminbox(NelderMead()),
  res = optimize(loss_in,lb,ub,p₀_bbo, Fminbox(BFGS()), autodiff = :forward,
                                Optim.Options(show_trace=true,
                                              f_tol = 1e-5,
                                              outer_iterations = 10))

  # res = optimize(loss_in,p₀_bbo, LBFGS(),
  #               Optim.Options(show_trace=true,
  #                             f_tol = 1e-4))
  # res = optimize(loss_in,lb,ub,bbo_minimizer, Fminbox(NelderMead()) , Optim.Options(show_trace=true))

  # SCIML_TRAIN otpimizers
  # res = DiffEqFlux.sciml_train(loss_in,p₀_bbo,ADAM(),cb=callback_in,maxiters=100)
  # res = DiffEqFlux.sciml_train(loss_in,start,BBO(),cb=callback_in,maxiters=100)

  rates = res.minimizer
  l = res.minimum

  # bootstrap_rates = bootstrapping(rates,data,model,dens,cons,t)
  # bootstrap_rates = []

  @show model
  # @show start
  @show l
  @show rates
  # @show best_fitness(res_bbo), best_candidate(res_bbo)
  # hes_for = ReverseDiff.hessian(loss_in,rates)
  # # hes_for = ForwardDiff.hessian(loss_in,rates)
  # hes_zyg = Zygote.hessian(loss_in,rates)
  # @show hes_for ≈ hes_zyg
  hes_zyg = []
  return rates , hes_zyg, l
end

# function bootstrapping(rates,data,model,dens,cons,t)
#   #create sample
#   data = n
#   B_size = 10
#   # rand_mat = rand(1:length(data),(length(data),B_size))
#
#   rand_mat = [rand(1:length(data),length(data)) for i in 1:B_size]
#
#   bootstrap_rates = zeros(size(rates))
#   # bootstrap_rates = []
#   for i in 1:B_size
#     resampled_data = data[rand_mat[i]]
#     resampled_t = t[rand_mat[i]]
#     sort_ind = sortperm(resampled_t)
#     resampled_data = resampled_data[sort_ind]
#     resampled_t = resampled_t[sort_ind]
#     unique_t = unique(resampled_t)
#
#     # loss_in = (p) -> loss_2(p,resampled_data,model,dens,cons,resampled_t,unique_t)
#     # callback_in = (p,l) -> callback(p,l,resampled_data,model,dens,cons,resampled_t,unique_t)
#     loss_in = (p) -> loss(p,data,model,dens,cons,t,unique(t))
#     callback_in = (p,l) -> callback(p,l,data,model,dens,cons,t,unique(t))
#
#     # @show loss_in(rates)
#     res = optimize(loss_in,-9, 2, rates,Fminbox(NelderMead()))
#     # @time res = DiffEqFlux.sciml_train(loss_in,rates,BFGS(),
#     #                                 lower_bounds = -9,
#     #                                 upper_bounds = 2,
#     #                                 cb = callback_in,
#     #                                 maxiters=50,
#     #                                 f_tol = 1e-4)
#
#     # @show res.minimizer
#     bootstrap_rates= hcat(bootstrap_rates,res.minimizer)
#     # push!(bootstrap_rates,res.minimizer)
#     # Σ_sol = run_model_2(model,p,resampled_data,dens,cons,unique_t)
#   end
#
#   return bootstrap_rates[:,2:end]
# end
#


text:
module ConstrainedFESpacesTests

using Test
using Gridap
using Gridap.ConstrainedFESpaces: VectorOfTwoParts

D = 2
dom = fill(4,D)
model = CartesianDiscreteModel(partition=tuple(dom...))

order = 1
diritag = "boundary"
_fespace = CLagrangianFESpace(Float64,model,order,diritag)

fixeddofs = [2,5]
fespace = ConstrainedFESpace(_fespace,fixeddofs)

ufun(x) = x[1] + x[2] + 2.0

trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2)

bh = FEBasis(fespace)
uh = zero(fespace)

cellmat = integrate(inner(bh,bh),trian,quad)
cellvec = integrate(inner(bh,uh),trian,quad)

nfree = 7

ndiri = 16

test_fe_space(fespace, nfree, ndiri, cellmat, cellvec, ufun)

g(x) = 2.0

U = TrialFESpace(fespace,g)

u(x) = x[1] + x[2] + 2.0

uh = interpolate(fespace,u)

#writevtk(trian,"trian",cellfields=["uh"=>uh])

nfree = 10
ndiri = 4
fixeddofs = [3,1,5,7]
nfixed = length(fixeddofs)
nfree_new = nfree - nfixed
dof_to_new_dof = zeros(Int,nfree)
is_fixed = fill(false,nfree)
is_fixed[fixeddofs] .= true
is_free = is_fixed.==false
dof_to_new_dof[is_free] .= 1:nfree_new
dof_to_new_dof[is_fixed] .= (-(ndiri+1)):-1:(-(ndiri+nfixed))

freevals=[20,40,60,80,90,100]
fixedvals=[30,10,50,70]

v = VectorOfTwoParts(dof_to_new_dof,freevals,fixedvals,ndiri)
@test all(v[is_free] .== freevals)
@test all(v[is_fixed] .== fixedvals)

end # module


text:
function two_species_approx_wavenumber(ω::Number, medium, species)
  sp1 = species[1]
  sp2 = species[2]

  if sp1.r > sp2.r
    @warn("method two_species_approx was designed for species[1] to be the smallest.")
  end
  vol = pi*sp1.num_density*sp1.r^2
  rhoeff =  medium.ρ*(medium.ρ + sp1.ρ - vol*(medium.ρ - sp1.ρ))/
  (medium.ρ + sp1.ρ + vol*(medium.ρ - sp1.ρ))

  kTS = wavenumber_low_volfrac(ω, medium, [sp1])
  mS = Medium(ρ=rhoeff, c= ω/kTS)
  wavenumber_low_volfrac(ω, mS, [sp2])
end

function two_species_approx_wavenumber(ωs::AbstractArray, medium, species)
  sp1 = species[1]
  sp2 = species[2]

  if sp1.r > sp2.r
    @warn("method two_species_approx was designed for species[1] to be the smallest.")
  end
  vol = pi*sp1.num_density*sp1.r^2
  rhoeff =  medium.ρ*(medium.ρ + sp1.ρ - vol*(medium.ρ - sp1.ρ))/
  (medium.ρ + sp1.ρ + vol*(medium.ρ - sp1.ρ))

  kTSs = wavenumber_low_volfrac(ωs, medium, [sp1])
  kTLs = [
  begin
    mS = Medium(ρ=rhoeff, c= ωs[i]/kTSs[i])
    wavenumber_low_volfrac(ωs[i], mS, [sp2])
  end
  for i in eachindex(ωs)]
end


text:
#=== elastic ===#
#=== write MSB ===#
function wMSB(
    model::Union{elmod2d},
    sou::Array{MTsource},
    sourcenumber::Int64,
    ParbinPath::Union{String, Void})

    nwf = model.nwf
    pkf = model.medium.pkf
    iflag = model.medium.iflag
    medium = model.medium
    ext = model.medium.ext

    if ParbinPath != nothing
        if ParbinPath[end-2:end] != "bin"
            error("Check output file extensions. It has to be in 'bin' format.")
        else
            fid = open(ParbinPath,"a+")
            write(fid,iflag)
            write(fid,ext)
            write(fid,nwf.BDntpp)
            write(fid,medium.pvel)
            write(fid,medium.svel)
            write(fid,medium.rho)
            write(fid,Float64(medium.dx))
            write(fid,Float64(medium.dz))
            write(fid,Float64(medium.dt))
            write(fid,medium.nT)
            write(fid,Float64(pkf))

            write(fid,sourcenumber)
            for isn = 1:sourcenumber
                write(fid,sou[i].nloc)
            end
            for isn = 1:sourcenumber
                write(fid,sou[i].waveform)
            end
            for isn = 1:sourcenumber
                write(fid,sou[i].not)
            end
            for isn = 1:sourcenumber
                write(fid,sou[i].mt)
            end
            close(fid)
        end
    end
end

#=== acoustic ===#
#=== write MSB ===#
function wMSB(
    model::Union{acmod2d},
    sou::Array{source},
    sourcenumber::Int64,
    ParbinPath::Union{String, Void})

    nwf = model.nwf
    pkf = model.medium.pkf
    iflag = model.medium.iflag
    medium = model.medium
    ext = model.medium.ext

    if ParbinPath != nothing
        if ParbinPath[end-2:end] != "bin"
            error("Check output file extensions. It has to be in 'bin' format.")
        else
            fid = open(ParbinPath,"a+")
            write(fid,iflag)
            write(fid,ext)
            write(fid,nwf.BDntpp)
            write(fid,medium.pvel)
            write(fid,medium.rho)
            write(fid,Float64(medium.dx))
            write(fid,Float64(medium.dz))
            write(fid,Float64(medium.dt))
            write(fid,medium.nT)
            write(fid,Float64(pkf))

            write(fid,sourcenumber)
            for i = 1:sourcenumber
                write(fid,sou[i].nloc)
            end
            for i = 1:sourcenumber
                write(fid,sou[i].waveform)
            end
            for i = 1:sourcenumber
                write(fid,sou[i].not)
            end
            for i = 1:sourcenumber
                if sou[i].tp == "expl"
                    write(fid,1)
                elseif sou[i].tp == "sfx"
                    write(fid,2)
                elseif sou[i].tp == "sfz"
                    write(fid,4)
                elseif sou[i].tp == "scx"
                    write(fid,5)
                elseif sou[i].tp == "scz"
                    write(fid,7)
                end
            end
            close(fid)
        end
    end
end


text:
"""
    SplitRightHandSideFunction <: AbstractRightHandSideFunction

A composite type for the split right-hand side of an [`AbstractInitialValueProblem`](@ref).

# Constructors
```julia
SplitRightHandSideFunction(L, rhs)
SRHS(args...; kwargs...)
```

# Arguments
- `L :: Union{Number, Matrix}` : linear part of right-hand side derivative.
- `rhs :: RightHandSideFunction` : nonlinear part of right-hand side derivative.
"""
mutable struct SplitRightHandSideFunction{L_T, rhs_T} <: AbstractRightHandSideFunction
    L::L_T
    rhs::rhs_T
end

@doc (@doc SplitRightHandSideFunction) SRHS(args...; kwargs...) = SplitRightHandSideFunction(args...; kwargs...)


text:
import Base: show

abstract type AbstractWizardDependency; end
struct InlineBuildDependency <: AbstractWizardDependency
    script::String
end
struct RemoteBuildDependency <: AbstractWizardDependency
    url::String
    script::Union{String, Nothing}
end
struct TarballDependency <: AbstractWizardDependency
    url::String
    hash::String
end

"""
    WizardState

Building large dependencies can take a lot of time. This state object captures
all relevant state of this function. It can be passed back to the function to
resume where we left off. This can aid debugging when code changes are
necessary.  It also holds all necessary metadata such as input/output streams.
"""
@Base.kwdef mutable struct WizardState
    step::Symbol = :step1
    ins::IO = stdin
    outs::IO = stdout

    # Filled in by step 1
    platforms::Union{Nothing, Vector{Platform}} = nothing

    # Filled in by step 2
    workspace::Union{Nothing, String} = nothing
    source_urls::Union{Nothing, Vector{String}} = nothing
    source_files::Union{Nothing, Vector{SetupSource}} = nothing
    dependencies::Union{Nothing, Vector{Dependency}} = nothing
    compilers::Union{Nothing, Vector{Symbol}} = nothing
    preferred_gcc_version::Union{Nothing, VersionNumber} = nothing
    preferred_llvm_version::Union{Nothing, VersionNumber} = nothing

    # Filled in by step 3
    history::Union{Nothing, String} = nothing
    files::Union{Nothing, Vector{String}} = nothing
    file_kinds::Union{Nothing, Vector{Symbol}} = nothing
    file_varnames::Union{Nothing, Vector{Symbol}} = nothing

    # Filled in by step 5c
    failed_platforms::Set{Platform} = Set{Platform}()
    # Used to keep track of which platforms we already visited
    visited_platforms::Set{Platform} = Set{Platform}()
    # Used to keep track of which platforms we have shown to work
    # with the current script. This gets reset if the script is edited.
    validated_platforms::Set{Platform} = Set{Platform}()

    # Filled in by step 7
    name::Union{Nothing, String} = nothing
    version::Union{Nothing, VersionNumber} = nothing
    github_api::GitHub.GitHubAPI = GitHub.DEFAULT_API
    travis_endpoint::String = "https://api.travis-ci.org/"
end

function serializeable_fields(::WizardState)
    # We can't serialize TTY's, in general.
    bad_fields = [:ins, :outs, :github_api]
    return [f for f in fieldnames(WizardState) if !(f in bad_fields)]
end

# Serialize a WizardState out into a JLD2 dictionary-like object
function serialize(io, x::WizardState)
    for field in serializeable_fields(x)
        io[string(field)] = getproperty(x, field)
    end

    # For unnecessarily complicated fields (such as `x.github_api`) store the internal data raw:
    io["github_api"] = string(x.github_api.endpoint)

    # For non-serializable fields (such as `x.ins` and `x.outs`) we just recreate them in unserialize().
end

function unserialize(io)
    x = WizardState()

    for field in serializeable_fields(x)
        setproperty!(x, field, io[string(field)])
    end

    # Manually recreate `ins` and `outs`.  Note that this just sets them to their default values
    x.ins = stdin
    x.outs = stdout
    x.github_api = GitHub.GitHubWebAPI(HTTP.URI(io["github_api"]))

    return x
end

function show(io::IO, x::WizardState)
    print(io, "WizardState [$(x.step)]")
end


text:
using OceanTurbulenceParameterEstimation.InverseProblems: vectorize, forward_run!, transpose_model_output
using CairoMakie

include("visualize_profile_predictions_utils.jl")

"""
    visualize!(ip::InverseProblem, parameters;
                    field_names = [:u, :v, :b, :e],
                    directory = pwd(),
                    filename = "realizations.png"
                    )

    For visualizing 1-dimensional time series predictions.
"""
function visualize!(ip::InverseProblem, parameters;
                    field_names = [:u, :v, :b, :e],
                    directory = pwd(),
                    filename = "realizations.png"
                    )

    isdir(directory) || makedir(directory)

    model = ip.simulation.model

    n_fields = length(field_names)

    observations = vectorize(ip.observations)

    forward_run!(ip, parameters)

    # Vector of SyntheticObservations objects, one for each observation
    predictions = transpose_model_output(ip.time_series_collector, observations)

    fig = Figure(resolution = (200*(length(field_names)+1), 200*(length(ip.observations)+1)), font = "CMU Serif")
    colors = [:black, :red, :blue]

    function empty_plot!(fig_position)
        ax = fig_position = Axis(fig_position)
        hidedecorations!(ax)
        hidespines!(ax, :t, :b, :l, :r)
    end

    calibration.observations[1].field_time_serieses.b[1].grid

    for (oi, observation) in enumerate(observations)

        i = oi + 1
        prediction = predictions[oi]

        targets = observation.times
        snapshots = round.(Int, range(1, length(targets), length=3))

        Qᵇ = model.tracers.b.boundary_conditions.top.condition[1,oi]
        Qᵘ = model.velocities.u.boundary_conditions.top.condition[1,oi]
        fv = model.coriolis[1,oi].f

        empty_plot!(fig[i,1])
        text!(fig[i,1], "Qᵇ = $(tostring(Qᵇ)) m⁻¹s⁻³\nQᵘ = $(tostring(Qᵘ)) m⁻¹s⁻²\nf = $(tostring(fv)) s⁻¹",
                    position = (0, 0),
                    align = (:center, :center),
                    textsize = 15,
                    justification = :left)

        for (j, field_name) in enumerate(field_names)

            middle = j > 1 && j < n_fields
            remove_spines = j == 1 ? (:t, :r) : j == n_fields ? (:t, :l) : (:t, :l, :r)
            axis_position = j == n_fields ? (ylabelposition=:right, yaxisposition=:right) : NamedTuple()

            j += 1 # reserve the first column for row labels

            info = field_guide[field_name]

            grid = observation.grid

            z = field_name ∈ [:u, :v] ? grid.zᵃᵃᶠ[1:grid.Nz] : grid.zᵃᵃᶜ[1:grid.Nz]

            to_plot = field_name ∈ keys(observation.field_time_serieses)

            if to_plot

                # field data for each time step
                obsn = getproperty(observation.field_time_serieses, field_name)
                pred = getproperty(prediction.field_time_serieses, field_name)

                ax = Axis(fig[i,j]; xlabelpadding=0, xtickalign=1, ytickalign=1,
                                            merge(axis_position, info.axis_args)...)

                hidespines!(ax, remove_spines...)

                middle && hideydecorations!(ax, grid=false)

                lins = []
                for (color_index, t) in enumerate(snapshots)
                    l = lines!([obsn[1,1,1:grid.Nz,t] .* info.scaling ...], z; color = (colors[color_index], 0.4))
                    push!(lins, l)
                    l = lines!([pred[1,1,1:grid.Nz,t] .* info.scaling ...], z; color = (colors[color_index], 1.0), linestyle = :dash)
                    push!(lins, l)
                end

                times = @. round((targets[snapshots] - targets[snapshots[1]]) / 86400, sigdigits=2)

                legendlabel(time) = ["LES, t = $time days", "Model, t = $time days"]
                Legend(fig[1,2:3], lins, vcat([legendlabel(time) for time in times]...), nbanks=2)
                lins = []
            else

                empty_plot!(fig[i,j])
            end
        end
    end

    save(joinpath(directory, filename), fig, px_per_unit = 2.0)
    return nothing
end

# visualize!(ip::InverseProblem, parameters; kwargs...) = visualize!(model_time_series(ip, parameters); kwargs...)

# function visualize_and_save!(calibration, validation, parameters, directory; fields=[:u, :v, :b, :e])
#         isdir(directory) || makedir(directory)

#         path = joinpath(directory, "results.txt")
#         o = open_output_file(path)
#         write(o, "Training relative weights: $(calibration.relative_weights) \n")
#         write(o, "Validation relative weights: $(validation.relative_weights) \n")
#         write(o, "Training default parameters: $(validation.default_parameters) \n")
#         write(o, "Validation default parameters: $(validation.default_parameters) \n")

#         write(o, "------------ \n \n")
#         default_parameters = calibration.default_parameters
#         train_loss_default = calibration(default_parameters)
#         valid_loss_default = validation(default_parameters)
#         write(o, "Default parameters: $(default_parameters) \nLoss on training: $(train_loss_default) \nLoss on validation: $(valid_loss_default) \n------------ \n \n")

#         train_loss = calibration(parameters)
#         valid_loss = validation(parameters)
#         write(o, "Parameters: $(parameters) \nLoss on training: $(train_loss) \nLoss on validation: $(valid_loss) \n------------ \n \n")

#         write(o, "Training loss reduction: $(train_loss/train_loss_default) \n")
#         write(o, "Validation loss reduction: $(valid_loss/valid_loss_default) \n")
#         close(o)

#         parameters = calibration.parameters.ParametersToOptimize(parameters)

#         for inverse_problem in [calibration, validation]

#             all_data = inverse_problem.observations
#             simulation = inverse_problem.simulation
#             set!(simulation.model, parameters)

#             for data_length in Set(length.(getproperty.(all_data, :t)))

#                 observations = [d for d in all_data if length(d.t) == data_length]
#                 days = observations[1].t[end]/86400

#                 new_ip = InverseProblem()

#                 visualize!(simulation, observations, parameters;
#                             fields = fields,
#                             filename = joinpath(directory, "$(days)_day_simulations.png"))
#             end
#         end

#     return nothing
# end


text:
module Twilert

using HTTPClient.HTTPC
using JSON

export alert

# Constants

const TWILIO_BASE = "https://api.twilio.com/2010-04-01"

# Structures

type Credentials
    sid::String
    token::String
end

type Message
    body::String
    to::String
    from::String
    creds::Credentials
end

# Helper functions

function authheader(creds::Credentials)
    sid = creds.sid
    token = creds.token
    encoded = base64("$sid:$token")

    value = "Basic $encoded"

    return ("Authorization", value)
end

function sendmsg(msg::Message)
    sid = msg.creds.sid

    body = msg.body
    to = msg.to
    from = msg.from

    url = "$TWILIO_BASE/Accounts/$sid/Messages"
    data = [("Body", body), ("To", to), ("From", from)]
    headers = [authheader(msg.creds)]
    req = HTTPC.RequestOptions(headers=headers)

    resp = HTTPC.post(url, data, req)

    if resp.http_code == 201
        return :ok
    else
        return :error
    end
end

function loadcfg()
    cfgfile = Nothing
    try
        cfgfile = open(".twilertrc", "r")
    catch
        try
            home = homedir()
            cfgfile = open("$home/.twilertrc", "r")
        catch
            error("Failed to load config")
        end
    end

    cfg = JSON.parse(cfgfile)
    close(cfgfile)

    # Test for minimal set of keys
    if ! ("sid" in keys(cfg) && "token" in keys(cfg))
        error("Missing keys in config")
    end

    return cfg
end

# Exported functions

function alert(body::String, to::String, from::String, creds::Credentials)
    # Basic validation
    if first(to) != '+' || first(from) != '+'
        error("Phone number must begin with '+'")
    end
    if length(body) > 1600
        error("Message body too long")
    end

    msg = Message(body, to, from, creds)
    result = sendmsg(msg)

    return result
end

function alert(body::String, to::String, from::String)
    cfg = loadcfg()
    creds = Credentials(cfg["sid"], cfg["token"])
    return alert(body, to, from, creds)
end

function alert(body::String, to::String)
    cfg = loadcfg()
    if ! ("from" in keys(cfg))
        error("No 'from' key in config")
    end
    from = cfg["from"]
    return alert(body, to, from)
end

function alert(body::String)
    cfg = loadcfg()
    if ! ("to" in keys(cfg))
        error("No 'to' key in config")
    end
    to = cfg["to"]
    return alert(body, to)
end

end


text:
module SimpleFields
using StaticArrays
import Base: eltype, zero, (+)

abstract type Field{U<:AbstractFloat} end
eltype(::Type{Field{U}}) where {U<:AbstractFloat} = U

struct CompositeField{U<:AbstractFloat} <: Field{U}
    components::Vector{Field}
    λ::U
    T::U
    ω::U
    tmax::U
    zero
end

function (field::CompositeField{U})(t::U) where {U<:AbstractFloat}
    f = field.zero
    for c in field.components
        f += c(t)
    end
    f
end

function (+)(a::F1, b::F2) where {F1<:Field,F2<:Field}
    zero(F1) == zero(F2) || error("Can only add fields of same type")
    a.λ == b.λ || error("Can only add fields of the same fundamental wavelength")
    CompositeField{eltype(a)}([a, b], a.λ, a.T, a.ω, max(a.tmax, b.tmax), zero(F1))
end

function (+)(a::CompositeField, b::F) where {F<:Field}
    a.zero == zero(typeof(b)) || error("Can only add fields of same type")
    a.λ == b.λ || error("Can only add fields of the same fundamental wavelength")
    CompositeField{eltype(a)}([a.components..., b], a.λ, a.T, a.ω, max(a.tmax, b.tmax), a.zero)
end

struct DelayedField{U<:AbstractFloat,F<:Field} <: Field{U}
    f::F
    tau::U
    λ::U
    T::U
    ω::U
    tmax::U
end
eltype(::Type{DelayedField{U,F}}) where {U<:AbstractFloat,F<:Field} = eltype(F)
(field::DelayedField{U})(t::U) where {U<:AbstractFloat} = field.f(t-field.tau)
zero(::Type{DelayedField{U,F}}) where {U<:AbstractFloat,F<:Field} = zero(F)

delay(f::F, tau::U) where {U<:AbstractFloat, F<:Field} = DelayedField{U,F}(f, tau, f.λ, f.T, f.ω, f.tmax+tau)

# λ_SI in meters
function fundamental(λ_SI::U) where {U<:AbstractFloat}
    λ = λ_SI/U(5.2917721067e-11)
    T = λ/U(137.035999)
    ω = U(2π/T)
    λ,T,ω
end

# I_SI in W/cm²
intensity(I_SI::U) where {U<:AbstractFloat} = I_SI/U(3.5094452e16)

# Gaussian envelope for the amplitude, fwhm of the intensity envelope
# in cycles.
gaussian(t::U, fwhm::U) where {U<:AbstractFloat} = exp(-t.^2/(2(fwhm/U(2*√(log(2))))^2))

function top_hat(t, ramp, tmax)
    ((t.>=ramp) .* (t.<tmax-ramp)) +
        ((t.>=0) .* (t.<ramp)).*(t/ramp) +
        (((t.>=tmax-ramp) .* (t.<tmax))).*((tmax-t)/ramp)
end

function top_hat_trunc(t, ramp, tmax)
    ((t.>=ramp) .* (t.<tmax)) +
        ((t.>=0) .* (t.<ramp)).*(t/ramp)
end

box(t::U, tmax::U, c::Bool) where {U<:AbstractFloat} = !c || (t >= 0 && t <= tmax) ? one(U) : zero(U)

struct LinearField{U<:AbstractFloat} <: Field{U}
    λ::U
    T::U
    ω::U
    tmax::U
    Ez::Function
    vanish::Bool
end

(field::LinearField{U})(t::U) where {U<:AbstractFloat} = field.Ez(t)*box(t,field.tmax,field.vanish)
zero(::Type{LinearField{U}}) where {U<:AbstractFloat} = zero(U)

struct TransverseField{U<:AbstractFloat} <: Field{U}
    λ::U
    T::U
    ω::U
    tmax::U
    Ez::Function
    Ex::Function
    vanish::Bool
end

(field::TransverseField{U})(t::U) where {U<:AbstractFloat} = SVector{2,U}(field.Ez(t),field.Ex(t)).*box(t,field.tmax,field.vanish)
zero(::Type{TransverseField{U}}) where {U<:AbstractFloat} = SVector{2,U}(zero(U), zero(U))

function gdd_params(λ_SI::U, τ₀::U, η::U = zero(U);
                    gdd_phase = false) where {U<:AbstractFloat}
    γ = τ₀^2/(8log(2))
    γ² = γ^2
    η² = η^2

    A = √(γ/√(γ² + η²))
    gdd_phase && (A *= exp(im*atan2(-η,γ)/2))

    a = one(U)./4*(γ/(γ² + η²))
    b = one(U)./2*(η/(γ² + η²))
    A,a,b
end

"""`λ_SI` in meters, `I_SI` in W/cm², `tmax` in cycles,
(Fourier-limited) `fwhm` in seconds, `cep` in units of π, `gdd` in
seconds squared. The pulse will be centered at `tmax/2`. If `vanish`
is `true`, the pulse will identically vanish outside the interval
[0,`tmax`]. `ξ` is the ratio of the minor to major axes of the ellipse
(if non-zero, will be 3d pulse automatically, if zero, will be 2d by
default, but can be overridden)."""
function pulse(λ_SI::U, I_SI::U,
               tmax::U, fwhm::U,
               q::U = one(U),
               cep::U = zero(U), gdd::U = zero(U);
               ξ::U = zero(U),
               threeD::Bool = ξ != zero(U),
               gdd_phase = false,
               vanish = true) where {U<:AbstractFloat}
    ξ != 0 && !threeD && error("Can't specify non-linearly polarized pulses in 2d")

    λ,T,ω = fundamental(λ_SI)
    T2 = U(2.41888430e-17T)^2
    I = intensity(I_SI)
    E₀ = √(I)

    A,a,b = gdd_params(λ_SI, fwhm, gdd; gdd_phase = gdd_phase)

    # Transfer to cycles as time base
    a *= T2
    b *= T2

    φ₀ = cep*π

    E = t -> E₀/√(1+ξ^2)*exp(im*2π*q*(t-tmax/2) + im*q*φ₀)*A*exp(-(a-im*b)*(t-tmax/2)^2)

    threeD ?
        TransverseField(λ, T, ω, tmax, t -> U(real(E(t))), t -> U(real(im*ξ*E(t))), vanish) :
        LinearField(λ, T, ω, tmax, t -> U(real(E(t))), vanish)
end

function strong_field_params(λ_SI::U, I_SI::U, Ip::U) where {U<:AbstractFloat}
    λ,T,ω = fundamental(λ_SI)
    I = intensity(I_SI)
    E₀ = √(I)
    Up = I/(4*ω^2)
    cutoff = 3.17Up - Ip
    keldysh = ω*√(2abs(Ip))/E₀
    Dict(:lambda => λ,
         :lambda_SI => λ_SI,
         :T => T,
         :omega => ω,
         :I => I,
         :I_SI => I_SI,
         :E0 => E₀,
         :A0 => E₀/ω,
         :alpha0 => E₀/ω^2,
         :Ip => Ip,
         :Up => Up,
         :keldysh => keldysh,
         :cutoff => cutoff,
         :cutoff_HO => cutoff/ω)
end

export Field, CompositeField, fundamental, delay, gdd_params, pulse, eltype, (+), top_hat, top_hat_trunc, strong_field_params
end # module


text:
function tree_from_leaves(leaf_nodes::Vector{String}, final_length::Int64)::Tuple{Vector{GeneralNode}, Int}
    my_node_list::Array{GeneralNode,1} = []

    # first create a list of leaf nodes
    for node_name in leaf_nodes
        nn = Node(node_name)

        push!(my_node_list,nn)
    end # for

    # Internal nodes are created using integers as names.
    temp_name::Int = length(my_node_list)+1

    # shuffle the node list to get a random tree
    Random.shuffle!(my_node_list)

    while length(my_node_list) > final_length
        # get two nodes
        # create a new mother node to which the two first nodes are added as children
        # add the new mother node to the list and reshuffle
        first_child::GeneralNode = pop!(my_node_list)
        first_child.inc_length = rand()#*0.1
        second_child::GeneralNode = pop!(my_node_list)
        second_child.inc_length = rand()
        curr_node::GeneralNode = Node(string(temp_name))

        add_child!(curr_node, first_child)
        add_child!(curr_node, second_child)
        push!(my_node_list, curr_node)
        temp_name += 1
        Random.shuffle!(my_node_list)
    end # while

    return my_node_list, temp_name
end

"""
    function create_tree_from_leaves(leaf_nodes::Vector{String}, rooted::Bool=false<:AbstractNode

Build a random tree from a list of leaf names. The tree is unrooted by default.

Returns the root node of the new tree.

* `leaf_nodes` : A list of strings which are used as the names of the leaves.

* `rooted` : Boolean indicating if the tree should be rooted
"""
function create_tree_from_leaves(leaf_nodes::Vector{String}, rooted::Bool=false)::GeneralNode{Float64, Int64}

    my_node_list, temp_name = rooted ? tree_from_leaves(leaf_nodes, 2) : tree_from_leaves(leaf_nodes, 3)
    root = Node(string(temp_name))
    if rooted
        lchild = pop!(my_node_list)
        lchild.inc_length = rand()
        rchild = pop!(my_node_list)
        rchild.inc_length = rand()
        add_child!(root, lchild)
        add_child!(root, rchild)
    else
        lchild = pop!(my_node_list)
        lchild.inc_length = rand()
        mchild = pop!(my_node_list)
        mchild.inc_length = rand()
        rchild = pop!(my_node_list)
        rchild.inc_length = rand()
        add_child!(root, lchild)
        add_child!(root, rchild)
        add_child!(root, mchild)
    end

    initialize_tree!(root)

    return root
end # function create_tree_from_leaves

"""
    function from_df(df::Array{T,2}, name_list::Vector{String})::GeneralNode{T, Int64} where T<:Real

This function takes an adjacency matrix and a vector of names
and turns it into a tree. No checks are performed.

Returns the root node of the tree.

* `df` : matrix with edge weights
* `name_list` : a list of names such that they match the column indices of the matrix
"""
function from_df(df::Array{T,2}, name_list::Vector{String})::GeneralNode{T, Int64} where T<:Real

    node_list::Vector{GeneralNode} = [Node(String(i)) for i in name_list]

    for (col_index, col) in enumerate(eachcol(df))
        for (row_index, entry) in enumerate(col)
            if entry != 0
                # there is a branch connecting these nodes
                node_list[col_index].inc_length = entry
                add_child!(node_list[row_index], node_list[col_index])
            end # end if
        end # end for
    end # end for
    i::Int = 0
    # find the root
    for (ind, n) in enumerate(node_list)
        if n.root == true
            i = ind
            break
        end # end if
    end # for
    node = node_list[i]

    # do some bookeeping here and set the binary representation of the nodes
    initialize_tree!(node)
    return node
end # function from_df

struct division
    arr::Array{<:Real, 2}
    names::Vector{<:AbstractString}
    number::Vector{Int}
end

Base.length(divi::division) = length(divi.names)

function single_inds(x::Array{Tuple{R, R}})::Vector{R} where R<:Real
    y = R[]
    for i in 1:length(x)
        if x[i][1] ∉ y && x[i][2] ∉ y
            push!(y, x[i][1])
        end
    end
    y
end

function create_division(covmat, names, numbers, tol)::Tuple{division, division}
    true_inds = single_inds(getfield.(findall(covmat .<= tol), :I))
    false_inds = [i for i in 1:size(covmat, 1) if i ∉ true_inds]
    @assert length(true_inds) > 0
    @assert length(false_inds) > 0
    d1 = division(covmat[true_inds, true_inds], names[true_inds], numbers[true_inds])
    d2 = division(covmat[false_inds, false_inds], names[false_inds], numbers[false_inds])
    d1, d2
end

function div2node(d1::division)
    n1 = Node(d1.names[1], d1.arr[1])
    n1.num = d1.number[1]
    n1
end

"""
    function cov2tree(covmat::Array{<:T, 2}, names::Vector{<:AbstractString}, numbers::Vector{Int64}; tol::Real=1e-7)::GeneralNode{T, Int64} where T<:Real

This function reconstructs a tree from a covariance matrix. It takes a covariance matrix,
a vector of leaf names and a vector of node numbers as mandatory arguments. The order of the two vectors must
correspond to the order of rows and columns in the covariance matrix. Optionally, the `tol` paramter indicates
the boundary below which all values are treated as zero.

Returns the root node of the tree corresponding to the supplied covariance matrix.

* `covmat` : covariance matrix
* `names` : a list of names such that they match the column/row indices of the matrix
* `numbers` : a list of Integers such that they match the column/row indices of the matrix
* `tol` : cut off value below which all values are treated as zero
"""
function cov2tree(covmat::Array{<:T, 2}, names::Vector{<:AbstractString}, numbers::Vector{Int64}; tol::Real=1e-7)::GeneralNode{T, Int64} where T<:Real

    if all(covmat .> tol)
        covmat = covmat .- minimum(covmat)
    elseif any(covmat .< 0.0)
        covmat[covmat .< 0] .= 0.0
    end
    try
        cov2tree_int(covmat, names, numbers, tol=tol)
    catch
        println(covmat)
    end

end

function cov2tree_int(covmat::Array{R, 2}, names::Vector{<:AbstractString}, numbers; tol::Real=1e-7)::GeneralNode{R, Int64} where R<:Real

    d1, d2 = create_division(covmat, names, numbers, tol)
    if length(d1) == 1 && length(d2) == 1
        # two singlular nodes left
        n1 = div2node(d1)
        n2 = div2node(d2)

    elseif length(d1) == 1 && length(d2) > 1
        # d1 is a leaf node, d2 is a cluster of nodes
        n1 = div2node(d1)
        min = minimum(d2.arr)
        n2 = cov2tree_int(d2.arr .- min, d2.names, d2.number, tol=tol)
        n2.inc_length = min

    elseif length(d1) > 1 && length(d2) == 1
        # d1 is a cluster of node2, d2 is a leaf node
        min = minimum(d1.arr)
        n1 = cov2tree_int(d1.arr .- min , d1.names, d1.number, tol=tol)
        n1.inc_length = min
        n2 = div2node(d2)

    else
        # d1 and d2 are a cluster of nodes
        min = minimum(d1.arr)
        n1 = cov2tree_int(d1.arr .- min, d1.names, d1.number, tol=tol)
        n1.inc_length = min
        min = minimum(d2.arr)
        n2 = cov2tree_int(d2.arr .- min, d2.names, d2.number, tol=tol)
        n2.inc_length = min
    end
    n = Node("root", one(R))
    add_child!(n, n1)
    add_child!(n, n2)
    return n
end

"""
      parsimony(tree::N, char::Dict{String,String}; gap::String="-")::Float64 where N<: GeneralNode

Do parsimony reconstruction for a tree and a set of characters. The characters are
supplied as the char dictionary with the key being the name of a leave and the value
string representation of the character.

Returns cost of the most parsimonius reconstruction.

* `tree` : root node of the tree
* `char` : characters for the reconstruction
* `gap` : optional identifier for gaps
"""
function parsimony(tree::N, char::Dict{String,String}; gap::String="-")::Float64 where N<: GeneralNode
      po = post_order(tree)
      states = filter(x -> x != "-", unique(collect(values(char))))
      nStates = length(states)
      if nStates == 0
            return 0.0
      end
      nNodes = length(po)
      costMatrix = zeros(nNodes, nStates) .- 1
      for node in po
            if node.nchild == 0
                  # i am a leave node
                  s = char[node.name]
                  costMatrix[node.num, :] .= s == gap ? 0 : abs.(log.(states .== s))
            else
                  # i am an internal node
                  daughters = node.children
                  localCosts = zeros(nStates, nStates, length(daughters))
                  for (j, s1) in enumerate(states),
                        (k, s2) in enumerate(states),
                        (l, d) in enumerate(daughters)

                        localCosts[j, k, l] = costMatrix[d.num, k] + Int(s1 != s2)
                  end
                  costMatrix[i, :] = vec(
                        mapslices(
                              sum,
                              mapslices(minimum, localCosts, dims = 2),
                              dims = 3,
                        ),
                  )
            end
      end
      minimum(costMatrix[tree.num, :])
end


text:
# Implementation of Dongarra et al., "Improving the accuracy...," SINUM 1983

# find peak index
function _normalizeInf!(x)
    n = length(x)
    s=1
    xm = abs(x[1])
    for j=1:n
        t = abs(x[j])
        if t > xm
            xm = t
            s = j
        end
    end
    x ./= xm
    s
end

# simple version for an isolated eigenvalue
"""
    rfeigen(A,x,λ,DT) => λnew, xnew, status

Improve the precision of a computed eigenpair `(x,λ)` for matrix `A`
via multi-precision iterative refinement, using more-precise real type `DT`.

The higher precision `DT` is only used for residual computation
(i.e. matrix-vector products), so this can be much faster than a full
eigensystem solution with precise eltype.  This method works on a
single eigenpair, and can fail spectacularly if there is another
eigenvalue nearby.
"""
function rfeigen(A::AbstractMatrix{T},
                 x::AbstractVector{Tx},
                 λ::Tλ,
                 DT::Type{<:AbstractFloat} = widen(real(T));
                 maxiter=5,
                 tol=1,
                 factor = lu,
                 scale = true,
                 verbose = false
                 ) where {T,Tλ,Tx}
    Tr = promote_type(promote_type(Tx,DT),Tλ)
    res = _rfeigen(A, x, λ, Tr, factor, maxiter, tol, scale, verbose)
    return res
end

"""
    rfeigen(A,λ,DT) => λnew, xnew, status

Like `rfeigen(A,x,λ,DT)`, but initialize `x` via one step of inverse
iteration.
"""
function rfeigen(A::AbstractMatrix{T},
                 λ::Tλ,
                 DT::Type{<:AbstractFloat} = widen(real(T));
                 maxiter=5,
                 tol=1,
                 factor = lu,
                 scale = true,
                 verbose = false
                 ) where {T,Tλ}
    # CHECKME: is this condition adequate?
    if issymmetric(A) && (Tλ <: Real)
        Tx = Tλ
    else
        Tx = complex(Tλ)
    end
    # There may not be a sampler for types of interest (hello, Quadmath)
    # so let's promote.
    x = normalize!((A - λ * I) \ Tx.(rand(size(A,1))))
    Tr = promote_type(promote_type(Tx,DT))
    res = _rfeigen(A, x, λ, Tr, factor, maxiter, tol, scale, verbose)
    return res
end

function _rfeigen(A::AbstractMatrix{T},
                      x::AbstractVector{Tx},
                      λ::Tλ,
                      ::Type{DT},
                      factor,
                      maxiter,
                      tol,
                      scale,
                      verbose
                      ) where {T,Tx,Tλ,DT}

    status = :unfinished
    λd = convert(DT,λ)
    n = LinearAlgebra.checksquare(A)
    tol1 = tol * eps(real(DT)) # CHECKME: factor of n?
    Ad = DT.(A)
    B = Ad - λd * I

    s = _normalizeInf!(x)
    xd = DT.(x)
    B[:,s] .= -xd
    Btmp = A - λ * I # do it over to get the type right
    Btmp[:,s] .= -x
    FB = factor(Btmp)

    # initial residual
    r::Vector{DT} = λd * xd - Ad * xd
    y = zeros(Tx,n)
    ys = zero(Tx)
    δ = similar(y)
    δp = similar(y)
    yp = similar(y)
    rt = similar(y)
    prevnorm = convert(real(DT),Inf)
    for p = 1:maxiter
        verbose && println("iter $p resnorm: ",norm(r))
        δ = FB \ Tx.(r) # ldiv!(δ,FB,Tx.(r))
        y .= y .+ δ
        δnorm = norm(δ)
        ynorm = norm(y)
        yp .= y
        ys = y[s]
        yp[s] = zero(T)
        if δnorm > prevnorm
            if δnorm > 10.0 * prevnorm
                status = :diverging
            else
                status = :stalled
            end
            verbose && println("$status at iter $p; early exit")
            break
        end
        prevnorm = δnorm
        if δnorm / ynorm < tol1
            status = :converged
            verbose && println("converged")
            # println("iter $p ratio ",δnorm / ynorm)
            break
        end
        δp .= δ
        δs = δ[s]

        r .= r .- B * DT.(δ)
        δp[s] = zero(T)
        r .= r .+ ys * δp .+ δs * y
    end
    xnew = x + yp
    λnew = λ + ys
    return λnew, xnew, status
end

"""
    rfeigen(A, S::Schur, idxλ, DT, maxiter=5) -> vals, vecs, status

Improves the precision of a cluster of eigenvalues of matrix `A`
via multi-precision iterative refinement, using more-precise real type `DT`.
Returns improved estimates of eigenvalues and vectors generating
the corresponding invariant subspace.

This method works on the set of eigenvalues in `S.values` indexed by `idxλ`.
It is designed to handle (nearly) defective cases, but will fail if
the matrix is extremely non-normal or the initial estimates are poor.
Note that `S` must be a true Schur decomposition, not a "real Schur".
"""
function rfeigen(A::AbstractMatrix{T}, S::Schur{TS}, idxλ,
                 DT = widen(real(T)), maxiter=5;
                 tol = 1, verbose = false) where {T, TS <: Complex}
    n = size(A,1)
    m = length(idxλ)
    λ = [S.values[idxλ[i]] for i in 1:m]
    Tw = promote_type(T,eltype(λ))
    DTw = promote_type(DT,Tw)
    tol1 = tol * eps(real(DT)) # CHECKME: factor of n?
    status = :unfinished

    # compute X, M
    # X is an adequately-conditioned set spanning the invariant subspace
    # M is an upper-triangular matrix of mixing coefficients

    # Most of the work is in the space of Schur vectors
    Z = zeros(Tw, n, m)
    z = zeros(Tw, n)
    idxz = Vector{Int}(undef, m)

    k = idxλ[1]
    if k==1
        z[1] = one(Tw)
    else
        x0 = (S.T[1:k-1,1:k-1] - λ[1] * I) \ S.T[1:k-1,k]
        z[1:k-1] .= -x0
        z[k] = one(Tw)
    end
    zm,zi = findmax(abs.(z))
    z ./= zm
    idxz[1] = zi
    Z[:,1] .= z
    M = zeros(Tw, m, m)
    M[1,1] = λ[1]
    for l=2:m
        kp = k
        k = idxλ[l]
        @assert k > kp
        x0 = (S.T[1:k-1,1:k-1] - λ[l]*I) \ S.T[1:k-1,k]
        X1 = (S.T[1:k-1,1:k-1] - λ[l]*I) \ Z[1:k-1,1:l-1]
        z[1:k-1] .= -x0
        z[k] = one(Tw)
        # pick mixing coeffts so that each vector has a good dominant index
        rhs = [-z[idxz[i]] for i=1:l-1]
        mtmp = [X1[idxz[i],j] for i=1:l-1,j=1:l-1]
        mv = mtmp \ rhs
        M[l,l] = λ[l]
        M[1:l-1,l] .= mv
        z[1:k-1] .= z[1:k-1] .+ X1 * mv
        zm, zi = findmax(abs.(z))
        idxz[l] = zi
        Z[:,l] .= z
    end
    X = S.Z * Z

    s = Int[]
    for j=1:m
        xm = zero(real(Tw))
        xi = 0
        for i=1:n
            xt = abs(X[i,j])
            if xt > xm && i ∉ s
                xm = xt
                xi = i
            end
        end
        push!(s, xi)
    end

    λd = DTw.(λ)
    Ad = DTw.(A)
    Xd = DTw.(X)
    Md = DTw.(M)
    # TODO: if cluster is tight enough, only need a singleton B
    # How tight is tight enough?
    B = Vector{Matrix{DTw}}(undef, m)
    for j=1:m
        B[j] = Ad - λd[j] * I
    end
    for j=1:m
        for i=1:m
            B[j][:,s[i]] .= -Xd[:,i]
        end
    end

    r = zeros(DTw, n, m)
    for j=1:m
        r[:,j] = λd[j] * Xd[:,j] - Ad * Xd[:,j]
        for i=1:j-1
            r[:,j] .= r[:,j] .+ M[i,j] * Xd[:,i]
        end
    end
    verbose && println("at iter 0 res norm = ", norm(r))
    FB0 = lu(Tw.(B[1]))
    FB = Vector{typeof(FB0)}(undef, m)
    FB[1] = FB0
    for j=1:m
        FB[j] = lu(Tw.(B[j]))
    end

    y = zeros(Tw, n, m)
    yp = zeros(Tw, n, m)
    ys = zeros(Tw, m, m)
    δp = zeros(Tw, n, m)
    δs = zeros(Tw, m, m)
    prevnorm = convert(real(DT),Inf)
    for p=1:maxiter
        δnorm = zero(real(DT))
        for j=1:m
            rhs = Tw.(r[:,j])
            # for jj=1:j
            #    rhs -= M[jj,j] * yp[:,jj]
            # end
            δ = FB[j] \ rhs
            δnorm += norm(δ)
            y[:,j] .+= δ
            yp[:,j] .= y[:,j]
            δp[:,j] .= δ
            for i=1:m
                δs[i,j] = δ[s[i]]
                δp[s[i],j] = zero(Tw)
                yp[s[i],j] = zero(Tw)
            end
            r[:,j] .= (r[:,j] - B[j] * DTw.(δ))
        end
        # there are occasional strange transients
        if (p > 3) && (δnorm > prevnorm)
            if δnorm > 10.0 * prevnorm
                status = :diverging
            else
                status = :stalled
            end
            verbose && println("$status at iter $p; early exit")
            break
        end
        prevnorm = δnorm
        for j=1:m
            for jj=1:m
                r[:,j] .= r[:,j] + (DTw(ys[jj,j]) * DTw.(δp[:,jj])
                                    + DTw(δs[jj,j]) * DTw.(yp[:,jj]))
            end
            for jj=1:j-1
                r[:,j] .= r[:,j] + Md[jj,j] * DTw.(δp[:,jj])
            end
        end
        # this update is done late to avoid doubly adding δ δ terms
        for j=1:m
            for i=1:m
                ys[i,j] = y[s[i],j]
            end
        end
        verbose && println("at iter $p res norm = ", norm(r))
        # println("ew: ",eigvals(Md + DTw.(ys)))
        verbose && println("DP subspace error: ",
                norm((Xd + DTw.(yp))*(Md+DTw.(ys)) - Ad * (Xd + DTw.(yp))))
        ynorm = norm(y)
        if δnorm / ynorm < tol1
            status = :converged
            verbose && println("converged")
            # println("iter $p ratio ",δnorm / ynorm)
            break
        end
    end
    Xnew = X + yp
    verbose && println("final subspace error norm: ", norm(Xnew*(M+ys) - A * Xnew))
    λbar = (1/m)*sum(λ)
    Mnew = Tw.(Md + DTw.(ys) - DTw(λbar) * I)
    dλ = eigvals(Mnew)
    λnew = λbar .+ dλ
    λnew, Xnew, status
end


text:
type Menu
    choice::Int
end

function Menu()

    ShowMenu()
    choice=ReadChoice()
    #println(typeof(choice))
    @match choice begin
        #49 => {LinearAlgebra()}
        #49 => {require("home/abhijith/.julia/JuliaTutorial/src/LinearAlgebra.jl")}
        #50 => {include("MatrixFactorizations.jl")}
        #51 => {include("TextMining.jl")}
        #52 => {include("RecSys.jl")}
        49 => println("You can include the module by typing: using JuliaTutorial.LinearAlgebra")
        50 => println("You can include the module by typing: using JuliaTutorial.MatrixFactoriations")
        51 => println("You can include the module by typing: using JuliaTutorial.TextMining")
        52 => println("You can include the module by typing: using JuliaTutorial.RecSys")
    end
end

function ShowMenu()

    println("Welcome to interactive Julia Tutorial series.")
    println("1.Linear Algebra.")
    println("2.Matrix Applications")
    println("3.Text Mining")
    println("4.Recommender Systems")
    println("Please select valid choice.")

end

function ReadChoice()

    cc=readbytes(STDIN,1) #choice is ascii of the value enetered

    #Check is the value of choice is between 1 and 4, i.e., 48 to 52
    cc=parseint("$(cc[1])")
    if cc < 49 || cc > 52
        Base.warn_once("Enter valid choice: 1,2,3,4")
        #println(choice)
        #ReadChoice()
    end
    return cc

end

#Menu()
println("Please use the method JuliaTutorial.Menu().")


text:
using Clustering

struct Filter
    w::Array{Float64, 3}
    n::Int
    threshold::Float64
    Filter(n, t) = new(weights(n), n, t)
end
Filter(n) = Filter(n, 0.5)

function cluster(w::AbstractArray{T, 3}, radius=3; min_neighbors=20, min_cluster_size=30) where {T}
    (n, m, t) = size(w)
    fg_clusters = zeros(w)
    fg_points = hcat(map(x -> collect(convert.(Float64, ind2sub(w, x))), find(w))...)
    C = dbscan(fg_points, radius; min_neighbors=min_neighbors, min_cluster_size=min_cluster_size)
    for x in C
        for y in x.core_indices
            index = convert.(Int, fg_points[:, y])
            fg_clusters[index...] = 1
        end
        # for y in x.boundary_indices
        #     index = convert.(Int, fg_points[:, y])
        #     video[index...] = colorant"red"
        # end
    end
    return fg_clusters
end

function weights(n)
    w = zeros(Float64, 2n+1, 2n+1, 2n+1)
    for i in 1:2n+1, j in 1:2n+1, k in 1:2n+1
        if n * sqrt(3) > sqrt((i-n-1)^2 + (j-n-1)^2 + (k-n-1)^2)
            w[i, j, k] = 1.0
        end
        # w[i, j, k] = n * sqrt(3) - sqrt((i-n-1)^2 + (j-n-1)^2 + (k-n-1)^2)
    end
    w .= w ./ sum(w)
    return w
end

function morphological_close(w::AbstractArray{T, 2}, n=1) where {T}
    out = zeros(w)
    temp = zeros(w)
    morphological_close!(out, temp, w, n)
end

function morphological_close!(out::AbstractArray{T, 2},
                              temp::AbstractArray{S, 2},
                              w::AbstractArray{R, 2}, n=1) where {R, S, T}
    @assert size(out) == size(w) == size(temp)
    (a, b) = size(w)
    n = 3
    m = 6
    for j in 1:b, i in 1:a
        if w[i, j] == zero(R)
            continue
        end
        h_start, h_end = calc_bounds(i, a, n)
        w_start, w_end = calc_bounds(j, b, m)
        for q in w_start:w_end, p in h_start:h_end
            temp[i+p, j+q] = one(S)
        end
    end
    for j in 1:b, i in 1:a
        if temp[i, j] == zero(S)
            continue
        end
        out[i, j] = one(T)
        h_start, h_end = calc_bounds(i, a, n)
        w_start, w_end = calc_bounds(j, b, m)
        for q in w_start:w_end, p in h_start:h_end
            if temp[i+p, j+q] == zero(S)
                out[i, j] = zero(T)
                break
            end
        end
    end
    return out
end

function apply(f::Filter, w::Array{T, 3}) where {T}
    result = zeros(Float64, size(w))
    (a, b, c) = size(w)
    n = f.n
    for k in 1:c, j in 1:b, i in 1:a
        if w[i, j, k] == zero(T)
            continue
        end
        h_start, h_end = calc_bounds(i, a, n)
        w_start, w_end = calc_bounds(j, b, n)
        d_start, d_end = calc_bounds(k, c, n)
        for r in d_start:d_end, q in w_start:w_end, p in h_start:h_end
            result[i+p, j+q, k+r] += f.w[n+1+p, n+1+q, n+1+r]
        end
    end
    for (idx, x) in enumerate(result)
        if x < f.threshold
            result[idx] = 0.0
        else
            result[idx] = 1.0
        end
    end
    return result
end

function scale(w::Array{UInt8, 3}, s::Array{Float64, 3})
    @assert size(w) == size(s)
    result = zeros(w)
    for (idx, x) in enumerate(w)
        result[idx] = convert(UInt8, round(float(x) * s[idx]))
    end
    return result
end

function calc_bounds(idx, upper_limit, filter_size)
    return (max(-filter_size, 1 - idx), min(filter_size, upper_limit - idx))
end


text:
# Test for gausslegendre().

    tol = 1e-14
    n = 42
    x = gausslegendre(n)
    @test length(x[1]) == n && length(x[2]) == n
    @test abs(x[1][37] - 0.910959724904127) < tol
    @test abs(x[2][37] - 0.030479240699603) < tol

    # Test a larger n (using ASY)
    n = 251
    x = gausslegendre(n)
    w = x[2]; x = x[1]
    @test all(length(x) == n) && all(length(w) == n)
    @test abs(x[37] + 0.896467746955729) < tol
    @test abs(w[37] - 0.005535005742012) < tol

    x = gausslegendre(1013);
    @test norm( x[1][2] - -0.999985167586110, Inf ) < tol
    @test norm( x[1][13] - -0.999218995240887, Inf ) < tol
    @test norm( x[2][2] - 1.681691163200592e-05, Inf ) < tol
    @test norm( x[2][13] - 1.224755309137936e-04, Inf ) < tol

    x = gausslegendre(10013);
    @test norm( x[1][2] - -0.999999848054223, Inf ) < tol
    @test norm( x[1][13] - -0.999991998242661, Inf ) < tol
    @test norm( x[2][2] - 1.722757320118474e-07, Inf ) < tol
    @test norm( x[2][13] - 1.254980540032470e-06, Inf ) < tol


text:
@testset "ContinuousView" begin
    for node in (Node{Node{LeafNode{Float64, 2, 2}, 2, 3}, 2, 3}(),
                 HashNode{Node{LeafNode{Float64, 2, 2}, 2, 3}, 2, 3}(),
                 DynamicNode{Node{LeafNode{Float64, 2, 2}, 2, 3}, 2}(8,8),
                 DynamicHashNode{Node{LeafNode{Float64, 2, 2}, 2, 3}, 2}(8,8))
        A = TreeView(node)

        v = continuousview(A, 11:30, 41:70)
        v .= OffsetArray(reshape(1:600, 20, 30), 11:30, 41:70)
        @test size(v) == (20, 30)
        @test A[11:30, 41:70] == reshape(1:600, 20, 30)

        v = spotview(A, 11, 41)
        @test size(v) == (5, 5)
        @test v == OffsetArray(reshape(1:600, 20, 30)[1:5, 1:5], 11:15, 41:45)
    end
end


text:
using MATLAB
using Base.Test

@test mstatement(:abc) == "abc"

@test mstatement(:(-a)) == "-(a)"

@test mstatement(:(a + b)) == "(a) + (b)"
@test mstatement(:(a - b)) == "(a) - (b)"
@test mstatement(:(a * b)) == "(a) * (b)"
@test mstatement(:(a / b)) == "(a) / (b)"
@test mstatement(:(a \ b)) == "(a) \\ (b)"
@test mstatement(:(a ^ b)) == "(a) ^ (b)"

@test mstatement(:(a .+ b)) == "(a) .+ (b)"
@test mstatement(:(a .- b)) == "(a) .- (b)"
@test mstatement(:(a .* b)) == "(a) .* (b)"
@test mstatement(:(a ./ b)) == "(a) ./ (b)"
@test mstatement(:(a .^ b)) == "(a) .^ (b)"

@test mstatement(:(sin(x))) == "sin(x)"
@test mstatement(:(hypot(x, y))) == "hypot(x, y)"
@test mstatement(:(plot3(x, y, z))) == "plot3(x, y, z)"

@test mstatement(:(x * 2)) == "(x) * (2)"

@test mstatement(:(y = x(1) + x(2))) == "y = (x(1)) + (x(2))"

@test mstatement(:([1 2 3])) == "[1, 2, 3]"
@test mstatement(:([x; y; z])) == "[x; y; z]"
@test mstatement(:([1 2 3; 4 5 6])) == "[[1, 2, 3]; [4, 5, 6]]"

@test mstatement(:(a(1))) == "a(1)"
@test mstatement(:(a(1, 1))) == "a(1, 1)"
@test mstatement(:(a[1, 1])) == "a(1, 1)"
@test mstatement(:(a[1:end])) == "a(1:end)"
@test mstatement(:(a{1})) == "a{1}"
@test mstatement(:(a{1, 1})) == "a{1, 1}"

@test mstatement(:(a.b)) == "a.b"
@test mstatement(:(1:2)) == "1:2"

@test mstatement(:(x')) == "(x)'"
@test mstatement(:(x.')) == "(x).'"

@test mstatement(:(plot(x, y, 'r'))) == "plot(x, y, 'r')"
@test mstatement(:(plot(x, y, "r"))) == "plot(x, y, 'r')"

@test mstatement(:((x, y) = fun(z))) == "[x, y] = fun(z)"