text:
using Knet
using ArgParse
using JLD
using MAT
using Images

include("lib/vocab.jl")
include("lib/convnet.jl")
include("lib/init.jl")
include("lib/base.jl")
include("lib/model.jl")
include("lib/imgproc.jl")
include("lib/util.jl")

function main(args)
    s = ArgParseSettings()
    s.description = "Caption generation script for the model (raw images)"

    @add_arg_table s begin
        ("--image"; help="input image")
        ("--vocabfile"; help="vocab file")
        ("--modelfile"; help="trained model file")
        ("--cnnfile"; help="convnet file for non-finetuned model")
        ("--beamsize"; arg_type=Int; default=1)
        ("--maxlen"; arg_type=Int; default=30; help="max sentence length")
        ("--nogpu"; action=:store_true)
        ("--lastlayer"; default="relu7")
        ("--imsize"; arg_type=Int; nargs=2; default=[224,224];
         help="new image sizes")
        ("--rgbmean"; arg_type=Float32; nargs=3;
         default=map(Float32, [123.68, 116.779, 103.939]))
    end

    # parse args
    println("Datetime: ", now()); flush(STDOUT)
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); println(o); flush(STDOUT)
    o[:loadfile] = o[:modelfile]

    newsize = tuple(o[:imsize]...)
    rgbmean = reshape(o[:rgbmean], (1,1,3))

    # load model
    o[:atype] = !o[:nogpu] ? KnetArray{Float32} : Array{Float32}
    w = load(o[:modelfile], "w")
    w = convert_weight(o[:atype], w)
    vocab = load(o[:vocabfile], "vocab")
    lossval = load(o[:modelfile], "lossval")
    s = initstate(o[:atype], size(w["wsoft"], 1), 1)

    # load convnet
    if o[:cnnfile] != nothing
        vggmat = matread(abspath(o[:cnnfile]))
        wcnn = get_vgg_weights(vggmat; last_layer=o[:lastlayer])
        w["wcnn"] = wcnn
    end

    @printf("Data and model loaded [%s]\n", now()); flush(STDOUT)

    # generate caption
    ti = now()
    @printf("Generation started (loss=%g,date=%s)\n", lossval, ti)
    flush(STDOUT)
    img = load(abspath(o[:image]))
    img = process_image(img, newsize; rgbmean=rgbmean)
    generated = generate(
        w, copy(s), img, vocab; maxlen=o[:maxlen], beamsize=o[:beamsize])
    report_generation(o[:image], generated, o[:beamsize])
    tf = now()
    @printf("\nTime elapsed: %s [%s]\n", tf-ti, tf)
end

function report_generation(filename, generated, beamsize)
    @printf("\nFilename: %s\n", filename)
    @printf("Generated: %s\n", generated)
    @printf("Beamsize: %d\n", beamsize)
    flush(STDOUT)
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
import MechGlueInterpolations
using Test
using MechGlueInterpolations: m
using MechGlueInterpolations: matrix_to_function
@testset "MechGlueInterpolations.jl" begin

    compm = [(complex(x-1, y-1.5)m) for x in 0.0:1:3, y in 0.0:1.5:3]
    xmin = 0.0m
    xmax = 3.0m
    ymin = 0.0m
    ymax = 3.0m
    nx= 3
    ny = 4
    matrix_to_function(compm, xmin, xmax, ymin, ymax, nx, ny)

end

ext: jl


text:
function operator_pass(x, state)
    if x.typ === CSTParser.BinaryOpCall
        if CSTParser.precedence(x.args[2]) in (8, 13, 14, 16) || x.args[2].fullspan == 0
            ensure_no_space_after(x.args[1], state, state.offset)
            ensure_no_space_after(x.args[2], state, state.offset + x.args[1].fullspan)
        else
            ensure_single_space_after(x.args[1], state, state.offset)
            ensure_single_space_after(x.args[2], state, state.offset + x.args[1].fullspan)
        end
    elseif x.typ === CSTParser.WhereOpCall
        ensure_single_space_after(x.args[2], state, state.offset + x.args[1].fullspan)
        n = length(x.args)
        offset = state.offset + x.args[1].fullspan + x.args[2].fullspan
        for i = 3:n
            a = x.args[i]
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    elseif x.typ === CSTParser.ColonOpCall
        offset = state.offset
        n = length(x.args)
        for (i, a) in enumerate(x.args)
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    elseif x.typ === CSTParser.ChainOpCall || x.typ == CSTParser.Comparison
        offset = state.offset
        n = length(x.args)
        for (i, a) in enumerate(x.args)
            if i != n
                ensure_single_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function tuple_pass(x, state)
    if x.typ === CSTParser.TupleH
        offset = state.offset
        n = length(x)
        for (i, a) in enumerate(x)
            i == n && continue
            if a.typ === CSTParser.PUNCTUATION && a.kind === Tokens.COMMA && !(x.args[i + 1].typ === CSTParser.PUNCTUATION)
                ensure_single_space_after(a, state, offset)
            elseif !(x.args[i + 1].typ === CSTParser.Parameters)
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function curly_pass(x, state)
    if x.typ === CSTParser.Curly
        offset = state.offset
        n = length(x)
        for (i, a) in enumerate(x)
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function call_pass(x, state)
    if x.typ === CSTParser.Call
        offset = state.offset + x.args[1].fullspan
        n = length(x)
        for (i, a) in enumerate(x)
            i == 1 && continue
            if a.typ === CSTParser.PUNCTUATION && a.kind === Tokens.COMMA
                ensure_single_space_after(a, state, offset)
            elseif i != n && !(x.args[i + 1].typ === CSTParser.Parameters)
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    elseif x.typ === CSTParser.Kw
        ensure_single_space_after(x.args[1], state, state.offset)
        ensure_single_space_after(x.args[2], state, state.offset + x.args[1].fullspan)
    end
end

function forloop_pass(x, state)
    if x.typ === CSTParser.For
        offset = state.offset + x.args[1].fullspan
        for a in x.args[2]
            # convert iter = I into iter in I
            if a.typ === CSTParser.BinaryOpCall && CSTParser.is_eq(a.args[2])
                offset += a.args[1].fullspan
                push!(state.edits, Edit(offset + 1:offset + 2, "in "))
                offset += a.args[2].fullspan
                offset += a.args[3].fullspan
            else
                offset += a.fullspan
            end
        end
    end
end

# TODO: move this to CSTParser?
function str_value(x)
    if x.typ === CSTParser.PUNCTUATION
        x.kind == Tokens.LPAREN && return "("
        x.kind == Tokens.LBRACE && return "{"
        x.kind == Tokens.LSQUARE && return "["
        x.kind == Tokens.RPAREN && return ")"
        x.kind == Tokens.RBRACE && return "}"
        x.kind == Tokens.RSQUARE && return "]"
        x.kind == Tokens.COMMA && return ","
        x.kind == Tokens.SEMICOLON && return ";"
        x.kind == Tokens.AT_SIGN && return "@"
        x.kind == Tokens.DOT && return "."
        return ""
    elseif x.typ === CSTParser.IDENTIFIER || x.typ === CSTParser.LITERAL || x.typ === CSTParser.OPERATOR || x.typ === CSTParser.KEYWORD
        return CSTParser.str_value(x)
    else
        s = ""
        for a in x
            s *= str_value(a)
        end
        return s
    end
end

function doc_pass(x, state)
    return
    if x.typ === CSTParser.MacroCall && x.args[1].typ === CSTParser.GlobalRefDoc
        # Align global docstring to:
        #
        # """
        # doc
        # """
        #
        # If the doc is single quoted i.e. "doc", they will be replaced with triple quotes.
        offset = state.offset + x.args[1].fullspan
        doc = x.args[2]

        val = str_value(doc)

        # s = escape_string(strip(val, ['\n']), "\$")
        s = strip(val, ['\n'])
        ds = string("\"\"\"\n", s, "\n", "\"\"\"\n")

        # Check if docstring needs to be edited
        if length(ds) != doc.fullspan || s != val
            # Remove previous docstring
            push!(state.edits, Edit(offset + 1:offset + doc.fullspan, ""))
            # Append newly formatted docstring
            push!(state.edits, Edit(offset, ds))
        end
    end
end

function comments_pass(text, state)
    ts = tokenize(text)
    while !Tokenize.Lexers.eof(ts)
        t = Tokenize.Lexers.next_token(ts)
        if Tokens.kind(t) == Tokens.COMMENT
            val = Vector{UInt8}(t.val)
            if length(val) > 1 && val[2] == 0x3d # multiline
                if !(val[3] in (0x20, 0x09)) # ensure single space at start
                    push!(state.edits, Edit(t.startbyte + 2, " "))
                end
                if length(val) > 5 # ensure single space at
                    for i = length(val) - 2:-1:3
                        if val[i] in (0x20, 0x09, 0x0a, 0x0d)
                            continue
                        else
                            push!(state.edits, Edit(t.startbyte .+ (i + 1:length(val) - 2), " "))
                            break
                        end
                    end
                end

            elseif length(val) > 1
                if !(val[2] in (0x20, 0x09)) # ensure single space at start
                    push!(state.edits, Edit(t.startbyte + 1, " "))
                end
            end
        end
    end
end

function lineends_pass(text, state)
    n = sizeof(text)
    io = IOBuffer(reverse(text))
    while !eof(io)
        c = read(io, Char)
        if c === '\n' && !eof(io)
            Base.peek(io) == 0x0d && read(io, Char) # crlf
            i1 = i2 = position(io)
            pc = Base.peek(io)
            while !eof(io) && pc in (0x20, 0x09)
                i2 = position(io)
                pc = read(io, UInt8)
            end
            if i1 != i2
                push!(state.edits, Edit((n - i2) + 1:(n - i1), ""))
            end
        end
    end
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/optimize/nonconvex_branch/stack_management.jl
# Contains the subroutines used for stack management. Namely, node_selection!,
# single_storage!, branch_node!, and fathom!.
#############################################################################

function _variable_infeasibility(m::GlobalOptimizer, i::Int)
    tsum = zero(Float64); tmin = typemax(Float64); tmax = typemin(Float64)
    d = m._branch_cost
    for j in _sparsity(BranchVar(), m, i)
        v = m._constraint_infeasiblity[j]
        tsum += v
        (v > tmax) && (tmax = v)
        (v < tmin) && (tmin = v)
    end
    return d.μ1*tsum + d.μ2*tmin + d.μ3*tmax
end

function _store_pseudocosts!(m::GlobalOptimizer, n::NodeBB)
    k = n.last_branch
    d = m._branch_cost
    if n.branch_direction == BD_POS
        d.ηp[k] += 1
        d.𝛹p[k] = n.lower_bound - m._lower_objective_value
        d.δp[k] = n.branch_extent
    elseif n.branch_direction == BD_NEG
        d.ηn[k] += 1
        d.𝛹n[k] = n.lower_bound - m._lower_objective_value
        d.δn[k] = n.branch_extent
    end
    return
end

function _lo_extent(m::GlobalOptimizer, xb::Float64, k::Int)
    c = _branch_cost(m)
    (c == BC_INFEASIBLE)   && return _variable_infeasibility(m, k)

    l = _lower_bound(BranchVar(), m, k)
    u = _upper_bound(BranchVar(), m, k)
    (c == BC_INTERVAL)     && return xb - l
    (c == BC_INTERVAL_REV) && return u - xb

    xlp = _lower_solution(BranchVar(), m, k)
    ρ = _cost_offset_β(m)*(u - l)
    y = max(min(xlp, u - ρ), l + ρ)
    return (c == BC_INTERVAL_LP) ? (y - l) : (u - y)
end

function _hi_extent(m::GlobalOptimizer, xb::Float64, k::Int)
    c = _branch_cost(m)
    (c == BC_INFEASIBLE)   && return _variable_infeasibility(m, k)

    l = _lower_bound(BranchVar(), m, k)
    u = _upper_bound(BranchVar(), m, k)
    (c == BC_INTERVAL)     && return u - xb
    (c == BC_INTERVAL_REV) && return xb - l

    xlp = _lower_solution(BranchVar(), m, k)
    ρ = _cost_offset_β(m)*(u - l)
    y = max(min(xlp, u - ρ), l + ρ)
    return (c == BC_INTERVAL_LP) ? (u - y) : (y - l)
end

@inline _score(x::T, y::T, μ::T) where T<:Real = (one(T) - μ)*min(x, y) + max(x, y)
@inline function _score(d::BranchCostStorage{T}, i::Int) where T<:Real
    _score(d.𝛹n[i]*d.δn[i], d.𝛹p[i]*d.δp[i], d.μ_score)
end

function _select_branch_variable_cost(m::GlobalOptimizer)
    return map_argmax(i -> score(m.branch_cost, i), 1:_variable_num(BranchVar(),m))
end

# In the case of diam(X) = Inf variables rel_diam may equal NaN, this should be set to
# a zero relative diameter to prevent "branching" on these variables for a functionally infinite
# amount of time
function rel_diam(m::GlobalOptimizer, i::Int)
    rd = _diam(BranchVar(), m, i)/diam(_working_variable_info(m, _bvi(m, i)))
    return isnan(rd) ? 0.0 : rd
end
function _select_branch_variable_width(m::GlobalOptimizer)
    map_argmax(i -> rel_diam(m,i), 1:_variable_num(BranchVar(), m))
end

"""
$(SIGNATURES)

Selects the variable to branch on psuedocost branching is used if
(parameter: `branch_pseudocost_on` = true).
"""
function select_branch_variable(t::ExtensionType, m::GlobalOptimizer)
    _branch_pseudocost_on(m) && return _select_branch_variable_cost(m)
    return _select_branch_variable_width(m)
end

"""
$(SIGNATURES)

Selects a point `xb` which is a convex combination (parameter:
`branch_cvx_factor`) of the solution to the relaxation and the midpoint of the
node. If this solution lies within (parameter: `branch_offset`) of a bound then
the branch point is moved to a distance of `branch_offset` from the bound.
"""
function select_branch_point(t::ExtensionType, m::GlobalOptimizer, i)
    l = _lower_bound(BranchVar(), m, i)
    u = _upper_bound(BranchVar(), m, i)
    s = _lower_solution(BranchVar(), m, i)
    α = _branch_cvx_α(m)
    b = _branch_offset_β(m)*(u - l)
    return max(l + b, min(u - b, α*s + (one(Float64) - α)*_mid(BranchVar(), m, i)))
end

"""
$(SIGNATURES)

Creates two nodes from `current_node` and stores them to the stack. Calls
`select_branch_variable(t, m)` and `select_branch_point(t, m, k)`.
"""
function branch_node!(t::ExtensionType, m::GlobalOptimizer)

    k = select_branch_variable(t, m)
    x = select_branch_point(t, m, k)
    n = m._current_node

    isfinite(n.last_branch) && _store_pseudocosts!(m, n)

    l_bound = max(n.lower_bound, m._lower_objective_value)
    u_bound = min(n.upper_bound, m._upper_objective_value)

    l_lbd = copy(n.lower_variable_bounds);  u_lbd = copy(n.lower_variable_bounds)
    l_ubd = copy(n.upper_variable_bounds);  u_ubd = copy(n.upper_variable_bounds)
    l_int = copy(n.is_integer);             u_int = copy(n.is_integer)

    flag = is_integer(BranchVar(), m, k)
    if flag
        l_int[k] = floor(x) != n.lower_variable_bounds[k]
        u_int[k] = ceil(x)  != n.upper_variable_bounds[k]
    end
    l_cont = flag ? !any(l_int) : true
    u_cont = flag ? !any(u_int) : true
    lx = flag ? floor(x) : x
    ux = flag ? ceil(x)  : x
    l_ubd[k] = lx
    u_lbd[k] = ux

    psuedo_cost_flag = _branch_pseudocost_on(m)
    l_ext = psuedo_cost_flag ? _lo_extent(m, lx, k) : zero(Float64)
    u_ext = psuedo_cost_flag ? _hi_extent(m, ux, k) : zero(Float64)

    l_idepth = (flag && l_cont) ? n.depth + 1 : n.cont_depth
    u_idepth = (flag && u_cont) ? n.depth + 1 : n.cont_depth

    n1 = NodeBB(l_lbd, l_ubd, l_int, l_cont, l_bound, u_bound, n.depth + 1, l_idepth, n.id + 1, BD_NEG, k, l_ext)
    n2 = NodeBB(u_lbd, u_ubd, u_int, u_cont, l_bound, u_bound, n.depth + 1, u_idepth, n.id + 2, BD_POS, k, u_ext)
    push!(m._stack, n1)
    push!(m._stack, n2)

    m._node_repetitions = 1
    m._maximum_node_id += 2
    m._node_count += 2
    return
end
branch_node!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = branch_node!(_ext_typ(m), m)

"""
$(SIGNATURES)

Stores the current node to the stack after updating lower/upper bounds.
"""
function single_storage!(t::ExtensionType, m::GlobalOptimizer)
    y = m._current_node
    m._node_repetitions += 1
    m._node_count += 1
    lower_bound = max(y.lower_bound, m._lower_objective_value)
    upper_bound = min(y.upper_bound, m._upper_objective_value)
    push!(m._stack, NodeBB(y.lower_variable_bounds, y.upper_variable_bounds, y.is_integer, y.continuous,
                           lower_bound, upper_bound, y.depth, y.cont_depth, y.id, y.branch_direction, y.last_branch, y.branch_extent))
    return
end
single_storage!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = single_storage!(_ext_typ(m), m)

"""
$(SIGNATURES)

Selects node with the lowest lower bound in stack.
"""
function node_selection!(t::ExtensionType, m::GlobalOptimizer)
    m._node_count -= 1
    m._current_node = popmin!(m._stack)
    return
end
node_selection!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = node_selection!(_ext_typ(m), m)

"""
$(SIGNATURES)

Selects and deletes nodes from stack with lower bounds greater than global
upper bound.
"""
function fathom!(t::ExtensionType, m::GlobalOptimizer)
    u = m._global_upper_bound
    continue_flag = !isempty(m._stack)
    while continue_flag
        n = maximum(m._stack)
        max_check = n.lower_bound > u
        if max_check
            popmax!(m._stack)
            m._node_count -= 1
        end
        continue_flag = !isempty(m._stack) && max_check
    end
    return
end
fathom!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = fathom!(_ext_typ(m), m)

"""
$(TYPEDSIGNATURES)

Creates an initial node with initial box constraints and adds it to the stack.
"""
function initialize_stack!(m::GlobalOptimizer)
    d = _working_variable_info.(m, m._branch_to_sol_map)
    push!(m._stack, NodeBB(lower_bound.(d), upper_bound.(d), is_integer.(d)))
    m._node_count = 1
    m._maximum_node_id += 1
    return
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Flux
using WindMixing
using BSON
using OceanParameterizations
using OrdinaryDiffEq

# Training data
train_files = ["-1e-3"]

𝒟train = data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=true)
PATH=pwd()
# FILE_PATH = joinpath(PATH, "training_output", "NDE_training_convective_adjustment_1sim_-1e-3_large_long_test.jld2")
FILE_PATH = joinpath("D:\\University Matters\\Massachusetts Institute of Technology\\CLiMA Project\\OceanParameterizations.jl\\training_output", "NDE_training_convective_adjustment_1sim_-1e-3_large_long_test.jld2")

uw_NN = BSON.load(joinpath(PATH, "NDEs", "uw_NN_large.bson"))[:neural_network]
vw_NN = BSON.load(joinpath(PATH, "NDEs", "vw_NN_large.bson"))[:neural_network]
wT_NN = BSON.load(joinpath(PATH, "NDEs", "wT_NN_large.bson"))[:neural_network]

train_epochs = [10, 10, 10, 10, 20]
train_tranges = [1:5:50, 1:5:100, 1:10:200, 1:20:400, 1:20:500]
train_optimizers = [[ADAM(0.01)] for i in 1:length(train_epochs)]
timestepper = ROCK4()

function train(FILE_PATH, train_files, train_epochs, train_tranges, train_optimizers, uw_NN, vw_NN, wT_NN, 𝒟train, timestepper)
    write_metadata_NDE_training(FILE_PATH, train_files, train_epochs, train_tranges, train_optimizers, uw_NN, vw_NN, wT_NN)

    for i in 1:length(train_epochs)
        uw_NN, vw_NN, wT_NN = train_NDE_convective_adjustment(uw_NN, vw_NN, wT_NN, 𝒟train, train_tranges[i], timestepper, train_optimizers[i], train_epochs[i], FILE_PATH, i, 1)
    end

end

train(FILE_PATH, train_files, train_epochs, train_tranges, train_optimizers, uw_NN, vw_NN, wT_NN, 𝒟train, timestepper)


text:
const SIMDFlag = Union{Bool, Symbol, Val{true}, Val{false}, Val{:ivdep}}

"""
    asval(T, x) -> v::Val{x′::T}

`asval(T, x)` validates `x` and returns `Val(x′)` such that `x′ isa T` and
`x == T(x′)`.
"""
function asval(::Type{T}, v::Val{x}) where {T,x}
    x isa T && return v
    error("expected type `$T` in `x` of `Val(x)`; got `x = $x`")
end
asval(::Type{T}, x::T) where {T} = Val(x)
asval(::Type{T}, x) where {T} = Val(T(x))

valueof(::Val{x}) where {x} = x
valueof(x) = x

# @generated foldlval(op, init, ::Val{N}) where {N} =
#     foldl((ex, i) -> :(op($ex, $i)), 1:N; init = :init)

@inline foldlargs(op, x) = x
@inline foldlargs(op, x1, x2, xs...) = foldlargs(op, op(x1, x2), xs...)

function zip_tuples(x::NTuple{N,Any}, xs::NTuple{N,Any}...) where {N}
    tuples = (x, xs...)
    return ntuple(i -> map(t -> t[i], tuples), Val{N}())
end

"""
    static_chunks(xs, ::Val{nchunks}) -> chunks::NTuple{nchunks}

Like `Iterators.partition(xs, cld(length(xs), nchunks))` but eagerly
partition `xs` and return a fixed-size tuple `chunks`. If `xs` is shorter
than `nchunks`, `chunks` contains empty chunks.
"""
function static_chunks end

function static_chunks(xs::AbstractArray, nchunks)
    basesize = cld(length(xs), valueof(nchunks))
    return ntuple(nchunks) do i
        @view xs[min(end+1,begin+(i-1)*basesize):min(end,begin+i*basesize-1)]
    end
end

function static_chunks(iter::Iterators.PartitionIterator{<:AbstractArray}, nchunks)
    xs = iter.c
    n = iter.n
    basesize = cld(length(xs), n * valueof(nchunks))
    return ntuple(nchunks) do i
        ys = @view xs[min(end+1,begin+(i-1)*basesize):min(end,begin+i*basesize-1)]
        Iterators.partition(ys, n)
    end
end

static_chunks(xs::Iterators.Zip, nchunks) =
    map(Base.splat(zip), zip_tuples(map(it -> static_chunks(it, nchunks), xs.is)...))


text:
#=
Copyright 2020 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

function calcNorms(vessels :: Array{Vessel, 1})
    norms = zeros(length(vessels),2)
    @inbounds for (i,v) in enumerate(vessels)
        err = v.P_l[:,4] .- v.P_t[:,4]
        norms[i] = sqrt(sum(err.^2))
    end
    return norms
end

"""
    computeConvError(edge_list, vessels :: Array{Vessel, 1})
"""
function computeConvError(vessels :: Array{Vessel, 1})
    current_norms = calcNorms(vessels)
    maxnorm, ci = findmax(current_norms)
    maxnormloc = vessels[ci[1]].label
    return maxnorm, maxnormloc
end

function printConvError(err::Float64, loc::String)
    err /= 133.332
    if err > 100.0
        @printf(" - Error norm > 100.00 mmHg\n")
    else
        @printf(" - Error norm = %6.2f mmHg @ %s\n", err, loc)
    end

end

checkConvergence(err::Float64, toll::Float64) = err/133.332 <= toll


text:
# test the add_measure_variable extensions
@testset "Expansion Extensions" begin
    # setup the model testing
    m = InfiniteModel()
    @infinite_parameter(m, par in [0, 1], supports = [0, 1])
    @infinite_parameter(m, pars[1:2] in [0, 1], supports = [0, 1])
    @variable(m, x, Infinite(par))
    @variable(m, y, Infinite(par, pars))
    @variable(m, x0, Point(x, 0))
    @variable(m, y0, SemiInfinite(y, 0, pars))
    tm = optimizer_model(m)
    @variable(tm, a)
    @variable(tm, b)
    @variable(tm, c)
    @variable(tm, d)
    data = transcription_data(tm)
    data.infvar_mappings[x] = [a, b, c]
    data.infvar_supports[x] = [(0.,), (0.5,), (1.,)]
    data.infvar_lookup[x] = Dict{Vector{Float64}, Int}([0] => 1, [0.5] => 2, [1] => 3)
    data.infvar_mappings[y] = [a, b, c, d]
    data.infvar_supports[y] = [(0., [0., 0.]), (0., [1., 1.]), (1., [0., 0.]), (1., [1., 1.])]
    data.infvar_lookup[y] = Dict{Vector{Float64}, Int}([0, 0, 0] => 1, [0, 1, 1] => 2,
                                                       [1, 0, 0] => 3, [1, 1, 1] => 4)
    data.infvar_mappings[y0] = [a, b]
    data.infvar_supports[y0] = [(0., [0., 0.]), (0., [1., 1.])]
    data.infvar_lookup[y0] = Dict{Vector{Float64}, Int}([0, 0, 0] => 1, [0, 1, 1] => 2)
    data.finvar_mappings[x0] = a
    key = Val(:TransData)
    IOTO.set_parameter_supports(tm, m)
    # test add_point_variable
    @testset "add_point_variable" begin
        # add one that was already added to the infinite model
        @test isequal(InfiniteOpt.add_point_variable(tm, x, Float64[0], key), x0)
        @test transcription_variable(x0) == a
        # add one that hasn't been added
        vref = GeneralVariableRef(m, -1, PointVariableIndex)
        @test isequal(InfiniteOpt.add_point_variable(tm, x, Float64[1], key), vref)
        @test transcription_variable(vref) == c
        # add one that has been added internally
        @test isequal(InfiniteOpt.add_point_variable(tm, x, Float64[1], key), vref)
        @test transcription_variable(vref) == c
    end
    # test add_semi_infinite_variable
    @testset "add_semi_infinite_variable" begin
        # add one that was already added to the infinite model
        var = SemiInfiniteVariable(y, Dict{Int, Float64}(1 => 0), [2, 3], [2])
        @test isequal(InfiniteOpt.add_semi_infinite_variable(tm, var, key), y0)
        @test transcription_variable(y0) == [a, b]
        # add a new one
        var = SemiInfiniteVariable(y, Dict{Int, Float64}(1 => 1), [2, 3], [2])
        vref = GeneralVariableRef(m, -1, SemiInfiniteVariableIndex)
        @test isequal(InfiniteOpt.add_semi_infinite_variable(tm, var, key), vref)
        @test isequal(data.semi_infinite_vars, [var])
        @test c in transcription_variable(vref)
        @test d in transcription_variable(vref)
        @test sort!(supports(vref)) == [([0., 0.], ), ([1., 1.], )]
        # add one that has already been added internally
        @test isequal(InfiniteOpt.add_semi_infinite_variable(tm, var, key), vref)
        @test isequal(data.semi_infinite_vars, [var])
        @test c in transcription_variable(vref)
        @test d in transcription_variable(vref)
        @test sort!(supports(vref)) == [([0., 0.], ), ([1., 1.], )]
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule pandoc_crossref_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("pandoc_crossref")
JLLWrappers.@generate_main_file("pandoc_crossref", UUID("f96e3c25-ae70-5588-8d8e-4a7ab3ba4c45"))
end  # module pandoc_crossref_jll


text:
"""
    Closed form expression for the wage for labor that produces intermediate goods
"""
function wage_int_goods(Π; par)
    Π_BS = Π[1] # firm revenue
    Π_LS = Π[2] # firm revenue
    Π_LU = Π[3] # firm revenue
    Π_BU = Π[4] # firm revenue

    w_BS = weight_of_type(:BS; par)
    w_LS = weight_of_type(:LS; par)
    w_LU = weight_of_type(:LU; par)
    w_BU = weight_of_type(:BU; par)

    return par.ζ * par.α_L_hat * (w_BS * Π_BS + w_LS * Π_LS + w_LU * Π_LU + w_BU * Π_BU)
end

function wage_gen_data(Π, D, l_G, D_I; par)
    Π_BS = Π[1] # firm revenue
    Π_LS = Π[2] # firm revenue
    Π_LU = Π[3] # firm revenue
    Π_BU = Π[4] # firm revenue

    D_BS = D[1] # data bundle
    D_LS = D[2] # data bundle
    D_LU = D[3] # data bundle
    D_BU = D[4] # data bundle

    D_I_BS = D_I[1] # data bundle
    D_I_LS = D_I[2] # data bundle
    D_I_LU = D_I[3] # data bundle
    D_I_BU = D_I[4] # data bundle

    l_G_BS = l_G[1] # data bundle
    l_G_LS = l_G[2] # data bundle
    l_G_LU = l_G[3] # data bundle
    l_G_BU = l_G[4] # data bundle

    w_BS = weight_of_type(:BS; par) # weights
    w_LS = weight_of_type(:LS; par) # weights
    w_LU = weight_of_type(:LU; par) # weights
    w_BU = weight_of_type(:BU; par) # weights

    D_G_BS = data_gen(l_G_BS; par, type = :BS)
    D_G_LS = data_gen(l_G_LS; par, type = :LS)
    D_G_LU = data_gen(l_G_LU; par, type = :LU)
    D_G_BU = data_gen(l_G_BU; par, type = :BU)

    c_BS = D_G_BS * Π_BS / D_BS * (D_BS / D_I_BS)^(1 / par.ε) # Put some things together
    c_LS = D_G_LS * Π_LS / D_LS * (D_LS / D_I_LS)^(1 / par.ε) # Put some things together
    c_LU = D_G_LU * Π_LU / D_LU * (D_LU / D_I_LU)^(1 / par.ε) # Put some things together
    c_BU = D_G_BU * Π_BU / D_BU * (D_BU / D_I_BU)^(1 / par.ε) # Put some things together

    return par.ζ * par.α_K_hat * (1 - par.ϕ) * ((1 - par.γ_S) * w_BS * c_BS + (1 - par.γ_S) * w_LS * c_LS + (1 - par.γ_U) * w_LU * c_LU + (1 - par.γ_U) * w_BU * c_BU)
end

"""
    Wage for data processing.
"""
function wage_gen_proc(Π, D; par)
    Π_BS = Π[1] # firm revenue
    Π_LS = Π[2] # firm revenue
    Π_LU = Π[3] # firm revenue
    Π_BU = Π[4] # firm revenue

    w_BS = weight_of_type(:BS; par) # weights
    w_LS = weight_of_type(:LS; par) # weights
    w_LU = weight_of_type(:LU; par) # weights
    w_BU = weight_of_type(:BU; par) # weights

    return par.ζ * par.α_K_hat * (par.γ_S * w_BS * Π_BS + par.γ_S * w_LS * Π_LS + par.γ_U * w_LU * Π_LU + par.γ_U * w_BU * Π_BU)
end


text:
using PyCall
using PyCall: hasproperty
using Test, Dates, Serialization

filter(f, itr) = collect(Iterators.filter(f, itr))
filter(f, d::AbstractDict) = Base.filter(f, d)

PYTHONPATH=get(ENV,"PYTHONPATH","")
PYTHONHOME=get(ENV,"PYTHONHOME","")
PYTHONEXECUTABLE=get(ENV,"PYTHONEXECUTABLE","")
@info "Python version $pyversion from $(PyCall.libpython), PYTHONHOME=$(PyCall.PYTHONHOME)\nENV[PYTHONPATH]=$PYTHONPATH\nENV[PYTHONHOME]=$PYTHONHOME\nENV[PYTHONEXECUTABLE]=$PYTHONEXECUTABLE"

roundtrip(T, x) = convert(T, PyObject(x))
roundtrip(x) = roundtrip(PyAny, x)
roundtripeq(T, x) = roundtrip(T, x) == x
roundtripeq(x) = roundtrip(x) == x

const math = pyimport("math")

struct TestConstruct
    x
end

pymodule_exists(s::AbstractString) = !ispynull(pyimport_e(s))

# default integer type for PyAny conversions
const PyInt = pyversion < v"3" ? Int : Clonglong

@testset "PyCall" begin
    # conversion of NumPy scalars before npy_initialized by array conversions (#481)
    np = pyimport_e("numpy")
    if !ispynull(np) # numpy is installed, so test
        let o = get(pycall(np."array", PyObject, 1:3), PyObject, 2)
            @test convert(Int32, o) === Int32(3)
            @test convert(Int64, o) === Int64(3)
            @test convert(Float64, o) === Float64(3)
            @test convert(Complex{Int}, o) === 3+0im
        end
    end

    # test handling of type-tuple changes in Julia 0.4
    import PyCall.pyany_toany
    @test pyany_toany(Int) == Int
    @test pyany_toany(PyAny) == Any
    @test pyany_toany(Tuple{Int,PyAny}) == Tuple{Int,Any}
    @test pyany_toany(Tuple{Int,Tuple{PyAny,Int8}}) == Tuple{Int,Tuple{Any,Int8}}
    @test pyany_toany(Tuple{PyAny,Int,Vararg{PyAny}}) == Tuple{Any,Int,Vararg{Any}}

    @test roundtripeq(17)
    @test roundtripeq(0x39)
    @test roundtripeq(true) && roundtripeq(false)
    @test roundtripeq(3.14159)
    @test roundtripeq(1.3+4.5im)
    @test roundtripeq(nothing)
    @test roundtripeq("Hello world")
    @test roundtripeq("Hëllö")
    @test roundtripeq("Hello \0\0\0world")
    @test roundtripeq("Hël\0\0lö")
    @test roundtripeq(Symbol, :Hello)
    @test roundtripeq(C_NULL) && roundtripeq(convert(Ptr{Cvoid}, 12345))
    @test roundtripeq([1,3,4,5]) && roundtripeq([1,3.2,"hello",true])
    @test roundtripeq([1 2 3;4 5 6]) && roundtripeq([1. 2 3;4 5 6])
    @test roundtripeq((1,(3.2,"hello"),true)) && roundtripeq(())
    @test roundtripeq(Int32)
    @test roundtripeq(Dict(1 => "hello", 2 => "goodbye")) && roundtripeq(Dict())
    @test roundtripeq(UInt8[1,3,4,5])
    @test roundtrip(3 => 4) == (3,4)
    @test roundtrip(Pair{Int,Int}, 3 => 4) == Pair(3,4)
    @test eltype(roundtrip([Ref(1), Ref(2)])) == typeof(Ref(1))

    @test pycall(PyObject(x -> x + 1), PyAny, 314158) == 314159
    @test PyObject(x -> x + 1)(314158) == 314159
    @test PyAny(PyObject(3)) == 3
    @test roundtrip(x -> x + 1)(314158) == 314159

    testkw(x; y=0) = x + 2*y
    @test pycall(PyObject(testkw), PyAny, 314157) == 314157
    @test pycall(PyObject(testkw), PyAny, 314157, y=1) == 314159
    @test roundtrip(testkw)(314157) == 314157
    @test roundtrip(testkw)(314157, y=1) == 314159

    # check type stability of pycall with an explicit return type
    @inferred pycall(PyObject(1).__add__, Int, 2)

    if PyCall.npy_initialized
        @test PyArray(PyObject([1. 2 3;4 5 6])) == [1. 2 3;4 5 6]
        let A = rand(Int, 2,3,4), B = rand(Bool, 2,3,4)
            @test convert(PyAny, PyReverseDims(A)) == permutedims(A, [3,2,1])
            @test convert(PyAny, PyReverseDims(BitArray(B))) == permutedims(B, [3,2,1])
        end
    end
    @test PyVector(PyObject([1,3.2,"hello",true])) == [1,3.2,"hello",true]
    @test PyDict(PyObject(Dict(1 => "hello", 2 => "goodbye"))) == Dict(1 => "hello", 2 => "goodbye")
    @test roundtripeq(BitArray([true, false, true, true]))

    let d = PyDict(Dict(1 => "hello", 34 => "yes" ))
        @test get(d.o, 1) == "hello"
        set!(d.o, 34, "goodbye")
        @test d[34] == "goodbye"
        @test sort!(keys(Int, d)) == sort!(collect(d.o.keys())) == sort!(collect(keys(d))) == [1, 34]
        @test eltype(d) == eltype(typeof(d)) == Pair{Int, String}
    end

    let d = Dict(zip(1:1000, 1:1000)), f
        f(k,v) = iseven(k) # For 0.6
        f(kv) = iseven(kv[1]) # For 0.7
        @test filter(f, d) == filter(f, PyDict(d)) == filter!(f, PyDict(d)) ==
              Dict(zip(2:2:1000, 2:2:1000))
    end

    @test roundtripeq(Any[1 2 3; 4 5 6])
    @test roundtripeq([])
    @test convert(Array{PyAny,1}, PyObject(Any[1 2 3; 4 5 6])) == Any[Any[1,2,3],Any[4,5,6]]
    if PyCall.npy_initialized
        @test roundtripeq(begin A = Array{Int}(undef); A[1] = 3; A; end)
    end
    @test convert(PyAny, PyObject(begin A = Array{Any}(undef); A[1] = 3; A; end)) == 3

    array2py2arrayeq(x) = PyCall.py2array(Float64,PyCall.array2py(x)) == x
    @test array2py2arrayeq(rand(3))
    @test array2py2arrayeq(rand(3,4))
    @test array2py2arrayeq(rand(3,4,5))

    @test roundtripeq(2:10) && roundtripeq(10:-1:2)
    @test roundtrip(2:2.0:10) == convert(Vector{Float64}, 2:2.0:10)

    @test math.sin(3) ≈ sin(3)

    @test collect(PyObject([1,"hello",5])) == [1,"hello",5]

    @test PyObject("hello") == PyObject("hello")
    @test PyObject("hello") != PyObject("hellö")
    @test PyObject(hash) == PyObject(hash)
    @test PyObject(hash) != PyObject(println)
    @test hash(PyObject("hello")) == hash(PyObject("hello"))
    @test hash(PyObject("hello")) != hash(PyObject("hellö"))
    @test hash(PyObject("hello")) != hash("hellö")
    @test hash(PyObject(hash)) == hash(PyObject(hash))
    @test hash(PyObject(hash)) != hash(PyObject(println))
    @test hash(PyObject(hash)) != hash(hash)

    # issue #92:
    let x = PyVector(PyAny[])
        py"lambda x: x.append(\"bar\")"(x)
        @test x == ["bar"]
    end

    @test roundtripeq(Dates.Date(2012,3,4))
    @test roundtripeq(Dates.DateTime(2012,3,4, 7,8,9,11))
    @test roundtripeq(Dates.Millisecond(typemax(Int32)))
    @test roundtripeq(Dates.Millisecond(typemin(Int32)))
    @test roundtripeq(Dates.Second, Dates.Second(typemax(Int32)))
    @test roundtripeq(Dates.Second, Dates.Second(typemin(Int32)))
    @test roundtripeq(Dates.Day, Dates.Day(999999999)) # max allowed day timedelta
    @test roundtripeq(Dates.Day, Dates.Day(-999999999)) # min allowed day timedelta

    # fixme: is there any nontrivial showable test we can do?
    @test !showable("text/html", PyObject(1))

    # in Python 3, we need a specific encoding to write strings or bufferize them
    # (http://stackoverflow.com/questions/5471158/typeerror-str-does-not-support-the-buffer-interface)
    pyutf8(s::PyObject) = pycall(s."encode", PyObject, "utf-8")
    pyutf8(s::String) = pyutf8(PyObject(s))

    # IO (issue #107)
    #@test roundtripeq(stdout) # No longer true since #250
    let buf = IOBuffer(read=false, write=true), obuf = PyObject(buf)
        @test !obuf.isatty()
        @test obuf.writable()
        @test !obuf.readable()
        @test obuf.seekable()
        obuf.write(pyutf8("hello"))
        obuf.flush()  # should be a no-op, since there's no flushing IOBuffer
        @test position(buf) == obuf.tell() == 5
        let p = obuf.seek(-2, 1)
            @test p == position(buf) == 3
        end
        let p = obuf.seek(0, 0)
            @test p == position(buf) == 0
        end
        @test String(take!(buf)) == "hello"
        obuf.writelines(["first\n", "second\n", "third"])
        @test String(take!(buf)) == "first\nsecond\nthird"
        obuf.write(b"möre stuff")
        @test String(take!(buf)) == "möre stuff"
        @test isopen(buf) == !obuf.closed == true
        obuf.close()
        @test isopen(buf) == !obuf.closed == false
    end
    let buf = IOBuffer("hello\nagain"), obuf = PyObject(buf)
        @test !obuf.writable()
        @test obuf.readable()
        @test obuf.readlines() == ["hello\n", "again"]
    end
    let buf = IOBuffer("hello\nagain"), obuf = PyObject(buf)
        @test codeunits(obuf.read(5)) == b"hello"
        @test codeunits(obuf.readall()) == b"\nagain"
    end
    let buf = IOBuffer("hello\nagain"), obuf = PyTextIO(buf)
        @test obuf.encoding == "UTF-8"
        @test obuf.read(3) == "hel"
        @test obuf.readall() == "lo\nagain"
    end
    let nm = tempname()
        open(nm, "w") do f
            # @test roundtripeq(f)  # PR #250
            pf = PyObject(f)
            @test pf.fileno() == fd(f)
            @test pf.writable()
            @test !pf.readable()
            pf.write(pyutf8(nm))
            pf.flush()
        end
        @test read(nm, String) == nm
    end

    # issue #112
    @test roundtripeq(Array, [1,2,3,4])
    @test roundtripeq(Array{Int8}, [1,2,3,4])

    # conversion of numpy scalars
    pyanycheck(x::Any) = pyanycheck(typeof(x), PyObject(x))
    pyanycheck(T, o::PyObject) = isa(convert(PyAny, o), T)
    @test pyanycheck(PyInt, PyVector{PyObject}(PyObject([1]))[1])
    @test pyanycheck(Float64, PyVector{PyObject}(PyObject([1.3]))[1])
    @test pyanycheck(ComplexF64, PyVector{PyObject}(PyObject([1.3+1im]))[1])
    @test pyanycheck(Bool, PyVector{PyObject}(PyObject([true]))[1])

    # conversions of Int128 and BigInt
    let i = 1234567890123456789 # Int64
        @test PyObject(i) - i == 0
    end
    let i = 12345678901234567890 # Int128
        @test PyObject(i) - i == 0
    end
    let i = BigInt(12345678901234567890), o = PyObject(i) # BigInt
        @test o - i == 0
        @test convert(BigInt, o) == i
        if pyversion >= v"3.2"
            @test PyAny(o) == i == convert(Integer, o)
            @test_throws InexactError convert(Int64, o)
        end
    end

    # bigfloat conversion
    if pymodule_exists("mpmath")
        for x in (big(pi), big(pi) + im/big(pi))
            @test pyanycheck(x)
            # conversion may not be exact since it goes through a decimal string
            @test abs(roundtrip(x) - x) < eps(BigFloat) * 1e3 * abs(x)
        end
    end
    @test convert(BigInt, PyObject(1234)) == 1234

    # hasproperty, getproperty, and propertynames
    py"""
    class A:
        class B:
            C = 1
    """
    A = py"A"
    @test hasproperty(A, "B")
    @test getproperty(A, "B") == py"A.B"
    @test :B in propertynames(A)
    @static if VERSION >= v"0.7-"
        @test A.B.C == 1
        @test_throws KeyError A.X
    end
    setproperty!(py"A.B", "C", 2)
    @test py"A.B.C" == 2

    # buffers
    let b = PyCall.PyBuffer(pyutf8("test string"))
        @test ndims(b) == 1
        @test (length(b),) == (length("test string"),) == (size(b, 1),) == size(b)
        @test stride(b, 1) == 1
        @test PyCall.iscontiguous(b) == true
    end

    let o = PyObject(1+2im)
        @test PyCall.hasproperty(o, :real) # replace by Base.hasproperty in the future
        @test :real in keys(o)
        @test o.real == 1
    end

    @testset "get-based sequence access" begin
        a1=[5,8,6]; a2=rand(3,4); a3=rand(3,4,5); o1=PyObject(a1); o2=PyObject(a2); o3=PyObject(a3)
        @test [get(o1,i-1) for i in eachindex(a1)] == a1
        @test [get(o1,-i) for i in eachindex(a1)] == reverse(a1)
        @test all(get(o2,0) == collect(a2[1,:]))
        @test length(o1) == length(o2) == length(o3) == 3
        set!(o1,-1,7)
        @test get(o1,-1) == 7

        # multiple indices are passed as tuples, but this is apparently
        # only supported by numpy arrays.
        if PyCall.npy_initialized
            @test [get(o2,(i-1,j-1)) for i=1:3, j=1:4] == a2
            @test [get(o3,(i-1,j-1,k-1)) for i=1:3, j=1:4, k=1:5] == a3
            @test all(get(o3, (1,2)) == collect(a3[2,3,:]))
            set!(o2, (1,2), 8)
            @test get(o2, (1,2)) == 8
            set!(o3, (1,2,3), 9)
            @test get(o3, (1,2,3)) == 9
        end
    end

    # list operations:
    let o = PyObject(Any[8,3])
        @test collect(push!(o, 5)) == [8,3,5]
        @test pop!(o) == 5 && collect(o) == [8,3]
        @test popfirst!(o) == 8 && collect(o) == [3]
        @test collect(pushfirst!(o, 9)) == [9,3]
        @test collect(prepend!(o, [5,4,2])) == [5,4,2,9,3]
        @test collect(append!(o, [1,6,8])) == [5,4,2,9,3,1,6,8]
        @test isempty(empty!(o))
    end
    let o = PyObject(Any[8,3])
        @test collect(append!(o, o)) == [8,3,8,3]
        push!(o, 1)
        @test collect(prepend!(o, o)) == [8,3,8,3,1,8,3,8,3,1]
    end

    # issue #216:
    @test length(collect(pyimport("itertools").combinations([1,2,3],2))) == 3

    # PyNULL and copy!
    let x = PyNULL(), y = copy!(x, PyObject(314159))
        @test convert(Int, x) == convert(Int, y) == 314159
    end
    @test ispynull(PyNULL())
    @test !ispynull(PyObject(3))
    @test ispynull(pydecref(PyObject(3)))

    ex = try
        pyimport("s p a m")
    catch ex
        ex
    end
    @test ex isa PyCall.PyError
    @test occursin("could not be found by pyimport", ex.msg)
    # Make sure we are testing ModuleNotFoundError here:
    if PyCall.pyversion >= v"3.6"
        @test pyisinstance(ex.val, pybuiltin("ModuleNotFoundError"))
    end

    @test !ispynull(pyimport_conda("inspect", "not a conda package"))
    import Conda
    if PyCall.conda
        # import pyzmq to test PR #294
        let already_installed = "pyzmq" ∈ Conda._installed_packages()
            @test !ispynull(pyimport_conda("zmq", "pyzmq"))
            @test "pyzmq" ∈ Conda._installed_packages()
            if !already_installed
                Conda.rm("pyzmq")
            end
        end
    end

    let x = 7
        py"""
        def myfun(x):
            return x + $x
        """
        @test py"1 + 2" == 3
        @test py"1 + $x" == 8
        @test py"1 + $(x^2)" == 50
        @test py"myfun"(10) == 17
    end

    # issue #352
    let x = "1+1"
        @test py"$x" == "1+1"
        @test py"$$x" == py"$$(x)" == 2
        @test py"7 - $$x - 7" == 0 # evaluates "7 - 1 + 1 - 7"
        @test py"7 - ($$x) - 7" == -2 # evaluates "7 - (1 + 1) - 7"
        @test py"1 + $$(x[1:2]) 3" == 5 # evals 1 + 1+ 3
    end

    # Float16 support:
    if PyCall.npy_initialized
        @test roundtripeq(Float16[17 18 Inf -Inf -0.0 0.0])
        @test isa(roundtrip(Float16[17]), Vector{Float16})
    end

    """
    foobar doc
    """
    foobar(x) = x+1

    # function attributes
    let o = PyObject(foobar)
        @test o.__name__ == o.func_name == string(foobar)
        @test o.__doc__ == o.func_doc == "foobar doc\n"
        @test o.__module__ == o.__defaults__ == o.func_defaults ==
              o.__closure__ == o.func_closure == nothing
    end

    # issue #345
    let weakdict = pyimport("weakref")."WeakValueDictionary"
        # (use weakdict for the value, since Python supports
        #  weak references to type objects)
        @test convert(Dict{Int,PyObject}, weakdict(Dict(3=>weakdict))) == Dict(3=>weakdict)
        @test get(weakdict(Dict(3=>weakdict)),3) == weakdict
    end

    # Weak ref support for pyjlwrap types
    let weakref = pyimport("weakref")
        bar = TestConstruct(1)
        o = PyObject(bar)
        @test PyCall.is_pyjlwrap(o)
        r = weakref.ref(o)
        @test weakref.getweakrefcount(o) == 1
    end

    # Expose python docs to Julia doc system
    py"""
    def foo():
        "foo docstring"
        return 0
    class bar:
        foo = foo
    """
    global foo354 = py"foo"
    global barclass = py"bar"
    # use 'content' since `Text` objects test equality by object identity
    @test @doc(foo354).content == "foo docstring"
    @test @doc(barclass.foo).content == "foo docstring"

    # binary operators
    for b in (4, PyObject(4))
        for op in (+, -, *, /, %, &, |, ^, <<, >>, ⊻)
            let x = op(PyObject(111), b)
                @test isa(x, PyObject)
                @test convert(PyAny, x) == op(111, 4)
            end
            @test convert(PyAny, op(b, PyObject(3))) == op(4, 3)
        end
    end
    @test convert(PyAny, PyObject(3)^4)  == 3^4 # literal integer powers
    @test convert(PyAny, PyObject(3)^0)  == 1   # literal integer powers
    @test convert(PyAny, PyObject(2)^-1) == 0.5 # literal integer powers
    # unary operators
    for op in (+, -, ~, abs)
        let x = op(PyObject(-3))
            @test isa(x, PyObject)
            @test convert(PyAny, x) == op(-3)
        end
    end
    # comparisons
    for x in (3,4,5), y in (3.0,4.0,5.0)
        for op in (<, <=, ==, !=, >, >=, isless, isequal)
            @test op(PyObject(x), PyObject(y)) == op(x, y)
            if op != isequal
                @test op(PyObject(x), y) == op(x, y)
            end
        end
    end

    # updating operators .+= etcetera
    let o = PyObject(Any[1,2]), c = o
        broadcast!(+, o, o, Any[3,4]) # o .+= x doesn't work yet in 0.7
        @test collect(o) == [1,2,3,4]
        GC.@preserve o c @test(PyPtr(o) == PyPtr(c)) # updated in-place
    end

    # more flexible bool conversions, matching Python "truth value testing"
    @test convert(Bool, PyObject(nothing)) === false
    @test convert(Bool, PyObject(0.0)) === false
    @test convert(Bool, PyObject(Any[])) === false
    @test convert(Bool, PyObject(17.3)) === true
    @test convert(Bool, PyObject(Any[0])) === true
    @test convert(Bool, PyVector{PyObject}(PyObject([false]))[1]) === false

    # serialization
    let py_sum_obj = pybuiltin("sum")
        b = IOBuffer()
        serialize(b, py_sum_obj)
        @test py_sum_obj == deserialize(seekstart(b))

        b = IOBuffer()
        serialize(b, PyNULL())
        @test PyNULL() == deserialize(seekstart(b))
    end

    # @pycall macro expands correctly
    _pycall = GlobalRef(PyCall,:pycall)
    @test macroexpand(@__MODULE__, :(@pycall foo(bar)::T)) == :($(_pycall)(foo, T, bar))
    @test macroexpand(@__MODULE__, :(@pycall foo(bar, args...)::T)) == :($(_pycall)(foo, T, bar, args...))
    @test macroexpand(@__MODULE__, :(@pycall foo(bar; kwargs...)::T)) == :($(_pycall)(foo, T, bar; kwargs...))

    # basic @pywith functionality
    fname = tempname()
    try
        @test begin
            @pywith pybuiltin("open")(fname,"w") as f begin
                f.write("test")
            end
            open(io->read(io, String), fname)=="test"
        end
    finally
        rm(fname,force=true)
    end

    @test occursin("integer", Base.Docs.doc(PyObject(1)).content)
    @test occursin("no docstring", Base.Docs.doc(PyObject(py"lambda x: x+1")).content)

    let b = rand(UInt8, 1000)
        @test(convert(Vector{UInt8}, pybytes(b)) == b
              == convert(Vector{UInt8}, pybytes(String(copy(b))))
              == convert(Vector{UInt8}, pybytes(codeunits(String(copy(b))))))
    end

    let t = convert(Tuple, PyObject((3,34)))
        @test isa(t, Tuple{PyObject,PyObject})
        @test t == (PyObject(3), PyObject(34))
    end
    for T in (Tuple{Vararg{PyAny}}, NTuple{2,PyInt}, Tuple{PyInt,PyInt}, Tuple{Vararg{PyInt}}, Tuple{PyInt,Vararg{PyInt}})
        let t = convert(T, PyObject((3,34)))
            @test isa(t, Tuple{PyInt,PyInt})
            @test t == (3,34)
        end
    end
    @test_throws BoundsError convert(NTuple{3,Int}, PyObject((3,34)))

    let p = PyCall.pickle(), buf = IOBuffer()
        p.dump("hello world", buf)
        p.dump(314159, buf)
        p.dump(Any[1,1,2,3,5,8], buf)
        @test p.load(seekstart(buf)) == "hello world"
        @test p.load(buf) == 314159
        @test p.load(buf) == [1,1,2,3,5,8]
    end

    # Test that we can call constructors on the python side
    @test pycall(PyObject(TestConstruct), PyAny, 1).x == 1

    # Test getattr fallback
    @test PyObject(TestConstruct(1)).x == 1
    @test_throws KeyError PyObject(TestConstruct(1)).y

    # iterating over Julia objects in Python:
    @test py"[x**2 for x in $(PyCall.pyjlwrap_new(1:4))]" ==
          py"[x**2 for x in $(x for x in 1:4)]" ==
          py"[x**2 for x in $(PyCall.jlwrap_iterator(1:4))]" ==
          [1,4,9,16]

    let o = PyObject("foo")
        @test pystr(o) == "foo"
        @test pyrepr(o) == "'foo'"
    end

    # pyfunction
    @test pyfunction(factorial, Int)(3) === PyInt(6)
    @test pyfunction(sin, Complex{Int})(3) === sin(3+0im)
    @test pyfunctionret(factorial, Float64, Int)(3) === 6.0
    @test pyfunctionret(factorial, nothing, Int)(3) === nothing
    @test PyCall.is_pyjlwrap(pycall(pyfunctionret(factorial, Any, Int), PyObject, 3))
    @test pyfunctionret(max, Int, Vararg{Int})(3,4,5) === PyInt(5)

    # broadcasting scalars
    let o = PyObject(3) .+ [1,4]
        @test o isa Vector{PyObject}
        @test o == [4,7]
    end

    # issue #533
    @test py"lambda x,y,z: (x,y,z)"(3:6,4:10,5:11) === (PyInt(3):PyInt(6), PyInt(4):PyInt(10), PyInt(5):PyInt(11))
end

######################################################################
#@pydef tests: type declarations need to happen at top level

# issue #389
@pydef mutable struct EmptyClass
end

# @pywith errors correctly handled
@pydef mutable struct IgnoreError
    function __init__(self, ignore)
        self.ignore = ignore
    end
    __enter__(self) = ()
    __exit__(self, typ, value, tb) = self.ignore
end

# @pydef example from README
@pydef mutable struct Doubler <: PyCall.builtin.AssertionError
    __init__(self, x=10) = (self.x = x)
    function my_method(self, arg1::Number)
        return arg1 + 20
    end
    type_str(self, obj::T) where T = string(T)
    x2.get(self) = self.x * 2
    function x2.set!(self, new_val)
        self.x = new_val / 2
    end
end

# @pydef with class variable
@pydef mutable struct ObjectCounter
    obj_count = 1 - 1
    function __init__(::PyObject)
        ObjectCounter.obj_count += 1
    end
end

@testset "pydef" begin
    d = Doubler(5)
    @test d.x == 5
    d.x2 = 30
    @test d.x == 15
    @test d.type_str(10) == string(PyInt)
    @test PyCall.builtin.isinstance(d, PyCall.builtin.AssertionError)

    @test_throws ErrorException @pywith IgnoreError(false) error()
    @test (@pywith IgnoreError(true) error(); true)

    @test ObjectCounter.obj_count == 0
    a = ObjectCounter()
    @test ObjectCounter.obj_count == 1
end

@testset "callback" begin
    # Returning existing PyObject in Julia should not invalidate it.
    # https://github.com/JuliaPy/PyCall.jl/pull/552
    anonymous = Module()
    Base.eval(
        anonymous, quote
            using PyCall
            obj = pyimport("sys")  # get some PyObject
        end)
    py"""
    ns = {}
    def set(name):
        ns[name] = $include_string($anonymous, name)
    """
    py"set"("obj")
    @test anonymous.obj != PyNULL()

    # Test above for pyjlwrap_getattr too:
    anonymous = Module()
    Base.eval(
        anonymous, quote
            using PyCall
            struct S
                x
            end
            obj = S(pyimport("sys"))
        end)
    py"""
    ns = {}
    def set(name):
        ns[name] = $include_string($anonymous, name).x
    """
    py"set"("obj")
    @test anonymous.obj.x != PyNULL()

    # Test above for pyjlwrap_iternext too:
    anonymous = Module()
    Base.eval(
        anonymous, quote
            using PyCall
            sys = pyimport("sys")
            obj = (sys for _ in 1:1)
        end)
    py"""
    ns = {}
    def set(name):
        ns[name] = list(iter($include_string($anonymous, name)))
    """
    py"set"("obj")
    @test anonymous.sys != PyNULL()
end

struct Unprintable end
Base.show(::IO, ::Unprintable) = error("show(::IO, ::Unprintable) called")
Base.show(::IO, ::Type{Unprintable}) = error("show(::IO, ::Type{Unprintable}) called")

py"""
def try_repr(x):
    try:
        return repr(x)
    except Exception as err:
        return err
"""

py"""
def try_call(f):
    try:
        return f()
    except Exception as err:
        return err
"""

@testset "throwing show" begin
    unp = Unprintable()
    @test_throws Exception show(unp)
    @test py"try_repr"("printable") isa String
    @test pyisinstance(py"try_repr"(unp), pybuiltin("Exception"))
    @test pyisinstance(py"try_call"(() -> throw(Unprintable())),
                       pybuiltin("Exception"))
end

@testset "PyIterator" begin
    arr = [1,2]
    o = PyObject(arr)
    c_pyany = collect(PyCall.PyIterator{PyAny}(o))
    @test c_pyany == arr
    @test c_pyany[1] isa Integer
    @test c_pyany[2] isa Integer

    c_f64 = collect(PyCall.PyIterator{Float64}(o))
    @test c_f64 == arr
    @test eltype(c_f64) == Float64

    i1 = PyObject([1])
    i2 = PyObject([2])
    l = PyObject([i1,i2])

    piter = PyCall.PyIterator(l)
    @test length(piter) == 2
    @test length(collect(piter)) == 2
    r1, r2 = collect(piter)
    @test getfield(r1, :o) === getfield(i1, :o)
    @test getfield(r2, :o) === getfield(i2, :o)

    @test Base.IteratorSize(PyCall.PyIterator(PyObject(1))) == Base.SizeUnknown()
    @test Base.IteratorSize(PyCall.PyIterator(PyObject([1]))) == Base.HasLength()

    # 594
    @test collect(zip(py"iter([1, 2, 3])", 1:3)) ==
    [(1, 1), (2, 2), (3, 3)]
    @test collect(zip(PyCall.PyIterator{Int}(py"iter([1, 2, 3])"), 1:3)) ==
    [(1, 1), (2, 2), (3, 3)]
    @test collect(zip(PyCall.PyIterator(py"[1, 2, 3]"o), 1:3)) ==
    [(1, 1), (2, 2), (3, 3)]
end

@test_throws PyCall.PyError py"__testing_pynamespace"

module __isolated_namespace
using PyCall
py"""
__testing_pynamespace = True
"""
get_testing_pynamespace() = py"__testing_pynamespace"
end

@test_throws PyCall.PyError py"__testing_pynamespace"
@test __isolated_namespace.get_testing_pynamespace()

@testset "atexit" begin
    script = """
    Base.load_path_setup_code()

    using PyCall

    pyimport("atexit").register() do
        println("atexit called")
    end
    """
    out = read(`$(Base.julia_cmd()) --startup-file=no -e $script`, String)
    @test occursin("atexit called", out)
end

include("test_pyfncall.jl")
include("testpybuffer.jl")
include("test_venv.jl")


text:
struct WrappedMap{T, A<:Union{AbstractMatrix, LinearMap}} <: LinearMap{T}
    lmap::A
    _issymmetric::Bool
    _ishermitian::Bool
    _isposdef::Bool
end
function WrappedMap(lmap::Union{AbstractMatrix{T}, LinearMap{T}};
    issymmetric::Bool = issymmetric(lmap),
    ishermitian::Bool = ishermitian(lmap),
    isposdef::Bool = isposdef(lmap)) where {T}
    WrappedMap{T, typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef)
end
function WrappedMap{T}(lmap::Union{AbstractMatrix, LinearMap};
    issymmetric::Bool = issymmetric(lmap),
    ishermitian::Bool = ishermitian(lmap),
    isposdef::Bool = isposdef(lmap)) where {T}
    WrappedMap{T, typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef)
end

# properties
Base.size(A::WrappedMap) = size(A.lmap)
LinearAlgebra.issymmetric(A::WrappedMap) = A._issymmetric
LinearAlgebra.ishermitian(A::WrappedMap) = A._ishermitian
LinearAlgebra.isposdef(A::WrappedMap) = A._isposdef

# multiplication with vector
A_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)
Base.:(*)(A::WrappedMap, x::AbstractVector) = *(A.lmap, x)

At_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) =
    (issymmetric(A) || (isreal(A) && ishermitian(A))) ? A_mul_B!(y, A.lmap, x) : At_mul_B!(y, A.lmap, x)

Ac_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) =
    ishermitian(A) ? A_mul_B!(y, A.lmap, x) : Ac_mul_B!(y, A.lmap, x)

# combine LinearMap and Matrix objects: linear combinations and map composition
Base.:(+)(A₁::LinearMap, A₂::AbstractMatrix) = +(A₁, WrappedMap(A₂))
Base.:(+)(A₁::AbstractMatrix, A₂::LinearMap) = +(WrappedMap(A₁), A₂)
Base.:(-)(A₁::LinearMap, A₂::AbstractMatrix) = -(A₁, WrappedMap(A₂))
Base.:(-)(A₁::AbstractMatrix, A₂::LinearMap) = -(WrappedMap(A₁), A₂)

Base.:(*)(A₁::LinearMap, A₂::AbstractMatrix) = *(A₁, WrappedMap(A₂))
Base.:(*)(A₁::AbstractMatrix, A₂::LinearMap) = *(WrappedMap(A₁), A₂)


text:
# pool

export ZeEventPool

mutable struct ZeEventPool
    handle::ze_event_pool_handle_t

    context::ZeContext

    function ZeEventPool(ctx::ZeContext, count::Integer, devs::ZeDevice...;
                         flags=0)
        desc_ref = Ref(ze_event_pool_desc_t(; flags, count))
        handle_ref = Ref{ze_event_pool_handle_t}()
        zeEventPoolCreate(ctx, desc_ref, length(devs), isempty(devs) ? C_NULL : [devs...], handle_ref)
        obj = new(handle_ref[], ctx)
        finalizer(obj) do obj
            zeEventPoolDestroy(obj)
        end
        obj
    end
end

Base.unsafe_convert(::Type{ze_event_pool_handle_t}, pool::ZeEventPool) = pool.handle

Base.getindex(pool::ZeEventPool, i::Integer) = ZeEvent(pool, i)

# event

export ZeEvent, append_wait!, signal, append_signal!, append_reset!, query, kernel_timestamp

mutable struct ZeEvent
    handle::ze_event_handle_t
    pool::ZeEventPool

    function ZeEvent(pool, index::Integer)
        desc_ref = Ref(ze_event_desc_t(; index=index-1))
        handle_ref = Ref{ze_event_handle_t}()
        zeEventCreate(pool, desc_ref, handle_ref)
        obj = new(handle_ref[], pool)
        finalizer(obj) do obj
            zeEventDestroy(obj)
        end
        obj
    end
end

Base.unsafe_convert(::Type{ze_event_handle_t}, event::ZeEvent) = event.handle

signal(event::ZeEvent) = zeEventHostSignal(event)
append_signal!(list::ZeCommandList, event::ZeEvent) = zeCommandListAppendSignalEvent(list, event)

Base.wait(event::ZeEvent, timeout::Number=typemax(UInt32)) =
    zeEventHostSynchronize(event, timeout)
append_wait!(list::ZeCommandList, events::ZeEvent...) =
    zeCommandListAppendWaitOnEvents(list, length(events), [events...])

Base.reset(event::ZeEvent) = zeEventHostReset(event)
append_reset!(list::ZeCommandList, event::ZeEvent) = zeCommandListAppendEventReset(list, event)

function query(event::ZeEvent)
    res = unsafe_zeEventQueryStatus(event)
    if res == RESULT_NOT_READY
        return false
    elseif res == RESULT_SUCCESS
        return true
    else
        throw_api_error(res)
    end
end

function kernel_timestamp(event)
    timestamp_ref = Ref{ze_kernel_timestamp_result_t}()
    zeEventQueryKernelTimestamp(event, timestamp_ref)

    # TODO: convert using ze_device_properties_t.timerResolution
    # TODO: mask by ze_device_properties_t.kernelTimestampValidBits
    # https://spec.oneapi.com/level-zero/latest/core/PROG.html#kernel-timestamp-events
    # but how to get the device?

    timestamp = timestamp_ref[]
    return (;
        :global => (
            start = timestamp._global.kernelStart == -1%UInt32 ? nothing : Int(timestamp._global.kernelStart),
            stop  = timestamp._global.kernelEnd == -1%UInt32 ?   nothing : Int(timestamp._global.kernelEnd)
        ),
        :context => (
            start = timestamp.context.kernelStart == -1%UInt32 ? nothing : Int(timestamp.context.kernelStart),
            stop  = timestamp.context.kernelEnd == -1%UInt32 ?   nothing : Int(timestamp.context.kernelEnd)
        )
    )
end


text:
export multigraph_edge_list

# The following is useful for when we want to plot a path, it takes in a path and outputs
# a Vecot
"""
`multigraph_edge_list(x::Path)`

This is a helper function which can be used in conjunction with `GraphRecipes.graphplot()` to visualise a path.

Output is tuple `(edge_list, label_map)` where
* `edge_list::Vector{Vector{Int}}` describes all edges, whereby edge_list[1] stores all edges starting at node 1, edge_list[2] those starting at 2 and so on...
* `label_map::Dict` this maps the integer labels of nodes back to the of the input data, e.g. label_map[1] = "a" means node 1 represents the vertex "a" in the path `p`.

To plot path with labels call

`graphplot(edge_list, names=label_map)`
"""
function multigraph_edge_list(p::Path{T}) where {T<:Union{Int, String}}
    edge_list = [Vector{Int}() for i in unique(p.data)]
    # We now loop over the vector and store the mapping from integers to strings
    ind = 1
    label_map = Dict{T, Int}(p[1]=>ind)
    ind+=1
    for i in 1:(length(p)-1)
        if ! haskey(label_map, p[i+1])
            label_map[p[i+1]] = ind
            ind += 1
        end
        push!(edge_list[label_map[p[i]]], label_map[p[i+1]])
    end
    # Invert the dictionary
    label_map = Dict{Int, T}(val => key for (key, val) in label_map)
    return edge_list, label_map
end

function multigraph_edge_list(x::Vector{Path{T}}) where {T<:Union{Int, String}}
    V = length(vertices(x))
    # @show V
    edge_lists = Vector{Vector{Vector{Int}}}()
    label_map = Dict{T, Int}()
    ind = 1
    for p in x
        tmp_edge_list = [Vector{Int}() for i in 1:V]
        if ! haskey(label_map, p[1])
            label_map[p[1]] = ind
            ind += 1
        end
        # @show tmp_edge_list, p
        for i in 1:(length(p)-1)
            if ! haskey(label_map, p[i+1])
                label_map[p[i+1]] = ind
                ind += 1
            end
            # @show label_map
            push!(tmp_edge_list[label_map[p[i]]], label_map[p[i+1]])
        end
        # @show tmp_edge_list
        push!(edge_lists, tmp_edge_list)
    end
    # Invert the dictionary
    label_map = Dict{Int, T}(val => key for (key, val) in label_map)
    return edge_lists, label_map
end

max_stars_repo_path: src/Data Processing/Multigraphs.jl

max_issues_repo_path: src/Data Processing/Multigraphs.jl

max_forks_repo_path: src/Data Processing/Multigraphs.jl


text:
using DataFrames
"""
    A_inverse(ped::AbstractString)
---
################################################################################
Given a 2-column pedigree file `ped`, this function returns a SparseArray of `A`
matrix inversed.

# Pedigree file format
- Each line has two `Int` numbers for `sire` and `dam`, respectively and in
  that order.
- Number `0` is for missing parents.
- Line numbers start from `1` are deemed as ID number.

See `dat/henderson1976.ped` for example.
"""
function A_inverse(ped::AbstractString)
    title("Calculate inversed A matrix with pedigree $ped")
    Ai = Any
    open(pipeline(`cat $ped`, `$abgBin/amat D`), "r") do io
        nid = parse(Int, readline(io))
        D = Float64[]
        for _ in 1:nid
            d = parse(Float64, readline(io))
            push!(D, d)
        end
        Ti = Int[]
        Tj = Int[]
        Td = Float64[]
        for line in eachline(io)
            i, j, d = split(line)
            push!(Ti, parse(Int, i))
            push!(Tj, parse(Int, j))
            push!(Td, parse(Float64, d))
        end
        T = sparse(Ti, Tj, Td)
        Ai = T'Diagonal(D)T
    end
    return Ai
end

"""
    inbreeding_coefficient(ped::AbstractString)
---
Given a 2-column pedigree file `ped`, this program calculate inbreeding
coefficient of everybody in the pedigree.
"""
function inbreeding_coefficient(ped::AbstractString)
    title("Calculate inbreeding coefficients with pedigreee $ped")
    ic = Float64[]
    open(pipeline(`cat $ped`, `$abgBin/amat F`), "r") do io
        for line in eachline(io)
            t = parse(Float64, line)
            push!(ic, t)
        end
    end
    return ic
end

"""
    A_matrix(ped::AbstractString, list::AbstractString="")
---
Given a 2-column pedigree file `ped`, this function calculate an **`A`** matrix.
IF a subset of ID `list` of the predigree is given, this gives a partial
**`A`** matrix. Inverse of this partial **`A`** can be done within `Julia`.
"""
function A_matrix(ped::AbstractString, list::AbstractString="")
    title("Calculate A matrix with pedigree $ped")
    ID = Int[]
    A = Float64[]
    println(ped)
    open(pipeline(`cat $ped`, `$abgBin/amat A $list`), "r") do io
        nid = parse(Int, readline(io))
        for _ in 1:nid
            id = parse(Int, readline(io))
            push!(ID, id)
        end
        for line in eachline(io)
            for elm in split(line)
                push!(A, parse(Float64, elm))
            end
        end
        reshape(A, (nid, nid))
    end
    return ID, A
end

"""
    sort_pedigree(pedigree::AbstractString, miss::AbstractString="0")
---
Given a pedigree file, this function will recode the pedigree to consective integers start from `1`.
Input file should have at least 3 columns.  The first 3 columns are `ID`, `sire` and `dam`.
ID in sire and dam columns are coded first.
Then the ID in ID columns are coded if their parents were already coded.
This procedure stops when all ID are coded, or a `loop`, i.e., an offspring is also its ancestor.
In the latter case, an error is reported.

By default, string `0` is for missing.  One can specify other string according to your data.

A 3-column DataFrame is returnd.  Row numbers are the recoded integers.  The 3rd column are their original
names.  The 1st and 2nd columns are their sires and dams.
"""
function sort_pedigree(pedigree::AbstractString, miss::AbstractString="0")
    title("Sort and recode pedigre $pedigree")
    oped = DataFrame(ID = String[], Sire = String[], Dam = String[]) # The pedigree to be sorted
    for line in eachline(pedigree)
        id, pa, ma = split(line)[[1,2,3]] # manual input meaning pedigree can have more columns
        push!(oped, (id, pa, ma))
    end

    coded = Dict()
    push!(coded, miss => 0)     # Name::String and its coded integer

    iid = 1                                                      # The code starts from 1
    sped = DataFrame(Sire = Int[], Dam = Int[], Name = String[]) # The sorted pedigree

    # ID only in the Sire and Dam column
    Parents = Set{String}()
    for id in [oped.Sire; oped.Dam]
        id == miss || push!(Parents, id)
    end
    for id in setdiff(Parents, oped.ID)
        push!(sped, (0, 0, id))
        push!(coded, id => iid) # ID in parent columns only are coded first
        iid += 1
    end

    left = prev = size(oped, 1)
    while left > 1              # got through the pedigree
        for (id, pa, ma) in eachrow(oped)
            if haskey(coded, pa) && haskey(coded, ma)
                push!(sped, (coded[pa], coded[ma], id))
                push!(coded, id => iid)
                iid += 1
                left -= 1
            end
        end
        if left == prev         # left cant't be reduced but >0
            warning("!!! Loop in the pedigree")
            break
        end
        prev = left
    end
    return sped
end


text:
# utilities for power spectrum estimation

function binning_matrix(left_bins, right_bins, weight_function_ℓ; lmax=nothing)
    lmax = isnothing(lmax) ? right_bins[end] : lmax
    bincut = right_bins .≤ lmax
    left_bins = left_bins[bincut]
    right_bings = right_bins[bincut]

    nbins = length(left_bins)
    P = zeros(nbins, lmax+1)

    for b in 1:nbins
        weights = weight_function_ℓ.(left_bins[b]:right_bins[b])
        norm = sum(weights)
        P[b, left_bins[b]+1:right_bins[b]+1] .= weights ./ norm
    end
    return P
end

function read_commented_header(filename; delim=" ", strip_spaces=true)
    header = CSV.read(filename, DataFrame; header=false, delim=delim, ignorerepeated=true,
        limit=1, types=String)
    if strip_spaces
        headers = [String(strip(header[1,"Column$(i)"])) for i in 1:ncol(header)]

    else
        headers = [header[1,"Column$(i)"] for i in 1:ncol(header)]
    end
    if headers[1] == "#"   # skip the #
        headers = headers[2:end]
    elseif headers[1][1] == '#'
        headers[1] = String(strip(headers[1][2:end]))
    end

    table = CSV.read(filename, DataFrame; comment="#", header=headers, delim=delim,
        ignorerepeated=true)
    return table
end

@doc raw"""
    nside2lmax(nside)

Get the Nyquist frequency from nside, ``3n_{\mathrm{side}} - 1``.
"""
nside2lmax(nside) = 3nside - 1

"""A dictionary that always returns one thing, no matter what key."""
struct ConstantDict{K,V} <: AbstractDict{K,V}
    c::V
end
Base.getindex(d::ConstantDict, key) = d.c

@doc raw"""
    function fitdipole(m::HealpixMap{T}, [w::HealpixMap{T}=1]) where T

Fit the monopole and dipole of a map.

# Arguments:
- `m::HealpixMap{T}`: map to fit
- `w::HealpixMap{T}`: weight map. Defaults to a FillArray of ones.

# Returns:
- `Tuple{T, NTuple{3,T}}`: (monopole, (dipole x, dipole y, dipole z))
"""
fitdipole

# basic reference, from  https://healpix.jpl.nasa.gov/html/subroutinesnode86.htm
@refimpl function fitdipole(m::HealpixMap{T}, w::HealpixMap{T}) where T
    upA = zeros(T,4,4)  # upper triangular version of A
    b = zeros(T, 4)
    for p ∈ eachindex(m.pixels)
        x, y, z = pix2vecRing(m.resolution, p)
        s = SA[one(T), x, y, z]
        for i ∈ 1:4
            b[i] += s[i] * w.pixels[p] * m.pixels[p]
            for j ∈ i:4
                upA[i,j] += s[i] * w.pixels[p] * s[j]
            end
        end
    end
    f = Symmetric(upA) \ b
    return f[1], (f[2], f[3], f[4])  # monopole, dipole
end

# more accurate version using carry bits
function fitdipole(m::HealpixMap{T}, w::HealpixMap{T}) where T
    # A and b
    upA = zeros(T,4,4)  # upper triangular version of A
    b = zeros(T, 4)

    # carry bits
    cA = zeros(T,4,4)
    cb = zeros(T,4)

    # using the Kahan-Babuska-Neumaier (KBN) algorithm for additional precision
    for p ∈ eachindex(m.pixels)
        x, y, z = pix2vecRing(m.resolution, p)
        s = SA[one(T), x, y, z]
        for i ∈ 1:4
            # sum the b vector
            inpb = s[i] * w.pixels[p] * m.pixels[p]
            sumb = b[i]
            tb = sumb + inpb
            if abs(sumb) ≥ abs(inpb)
                cb[i] += (sumb - tb) + inpb
            else
                cb[i] += (inpb - tb) + sumb
            end
            b[i] = tb

            for j ∈ i:4
                inpA = s[i] * w.pixels[p] * s[j]
                sumA = upA[i,j]
                tA = sumA + inpA
                if abs(sumA) ≥ abs(inpA)
                    cA[i,j] += (sumA - tA) + inpA
                else
                    cA[i,j] += (inpA - tA) + sumA
                end
                upA[i,j] = tA
            end
        end
    end
    f = Symmetric(upA .+ cA) \ (b .+ cb)
    return f[1], (f[2], f[3], f[4])  # monopole, dipole
end

function fitdipole(m::HealpixMap{T,O}) where {T,O}
    fitdipole(m, HealpixMap{T,O}(Ones(length(m.pixels))))
end

"""
    subtract_monopole_dipole!(map_in, monopole, dipole)

# Arguments:
- `map_in::HealpixMap`: the map to modify
- `monopole::T`: monopole value
- `dipole::NTuple{3,T}`: dipole value
"""
function subtract_monopole_dipole!(map_in::HealpixMap,
        monopole::T, dipole::NTuple{3,T}) where T
    res = map_in.resolution
    for p ∈ eachindex(map_in.pixels)
        x, y, z = pix2vecRing(res, p)
        t = monopole + dipole[1]*x + dipole[2]*y + dipole[3]*z
        map_in.pixels[p] -= t
    end
    map_in
end

"""
    synalm([rng=GLOBAL_RNG], Cl::AbstractArray{T,3}, nside::Int) where T

# Arguments:
- `Cl::AbstractArray{T,3}`: array with dimensions of comp, comp, ℓ
- `nside::Int`: healpix resolution

# Returns:
- `Vector{Alm{T}}`: spherical harmonics realizations for each component

# Examples
```julia
nside = 16
C0 = [3.  2.;  2.  5.]
Cl = repeat(C0, 1, 1, 3nside)  # spectra constant with ℓ
alms = synalm(Cl, nside)
```
"""
function synalm(rng::AbstractRNG, Cl::AbstractArray{T,3}, nside::Int) where T
    ncomp = size(Cl,1)
    @assert ncomp > 0
    alms = [Alm{Complex{T}}(3nside-1, 3nside-1) for i in 1:ncomp]
    synalm!(rng, Cl, alms)
    return alms
end
synalm(Cl::AbstractArray{T,3}, nside::Int) where T = synalm(Random.default_rng(), Cl, nside)

"""
    synalm!([rng=GLOBAL_RNG], Cl::AbstractArray{T,3}, alms) where T

In-place synthesis of spherical harmonic coefficients, given spectra.

# Arguments:
- `Cl::AbstractArray{T,3}`: array with dimensions of comp, comp, ℓ
- `alms::Union{NTuple{N,Alm}, Vector{Alm}}`: array of Alm to fill

# Examples
```julia
nside = 16
C0 = [3.  2.;  2.  5.]
Cl = repeat(C0, 1, 1, 3nside)  # spectra constant with ℓ
alms = [Alm{Complex{Float64}}(3nside-1, 3nside-1) for i in 1:2]
synalm!(Cl, alms)
```
"""
function synalm!(rng::AbstractRNG, Cl::AbstractArray{T,3},
        alms::Union{NTuple{N,Alm}, Vector{A}};
        lmin=0, lmax=-1) where {T, N, A <: Alm}
    # This implementation could be 1.2x faster by storing the cholesky factorization, but
    # typically you also perform two SHTs with each synalm, which dominates the cost.

    ncomp = size(Cl,1)
    @assert ncomp > 0
    @assert size(Cl,1) == size(Cl,2)
    @assert size(alms,1) > 0
    lmax = (lmax < 0) ? first(alms).lmax : lmax

    # first we synthesize just a unit normal for alms. we'll adjust the magnitudes later
    for comp in 1:ncomp
        randn!(rng, alms[comp].alm)
        for ℓ in 0:(lmin-1)
            for m in 0:ℓ
                i_alm = almIndex(alms[comp], ℓ, m)  # compute alm index
                alms[comp].alm[i_alm] = zero(T)
            end
        end
    end
    𝐂 = Array{T,2}(undef, (ncomp, ncomp))  # covariance for this given ℓ
    h𝐂 = Hermitian(Array{T,2}(undef, (ncomp, ncomp)))  # hermitian buffer
    alm_out = zeros(Complex{T}, ncomp)
    alm_in = zeros(Complex{T}, ncomp)

    for ℓ in lmin:lmax
        # build the 𝐂 matrix for ℓ
        for cᵢ in 1:ncomp, cⱼ in 1:ncomp
            𝐂[cᵢ, cⱼ] = Cl[cᵢ, cⱼ, ℓ+1]
        end

        if iszero(𝐂)
            for m in lmin:ℓ
                i_alm = almIndex(alms[1], ℓ, m)  # compute alm index
                for comp in 1:ncomp  # copy buffer back into the alms
                    alms[comp].alm[i_alm] = zero(T)
                end
            end
        else
            h𝐂 .= Hermitian(𝐂)
            cholesky_factorizable = isposdef!(h𝐂)
            if !cholesky_factorizable
                𝐂 .= sqrt(𝐂)
                for m in 0:ℓ
                    i_alm = almIndex(alms[1], ℓ, m)  # compute alm index
                    for comp in 1:ncomp  # copy over the random variates into buffer
                        alm_in[comp] = alms[comp].alm[i_alm]
                    end
                    mul!(alm_out, 𝐂, alm_in)
                    for comp in 1:ncomp  # copy buffer back into the alms
                        alms[comp].alm[i_alm] = alm_out[comp]
                    end
                end
            else
                # cholesky!(h𝐂)  # we already cholesky'd by calling isposdef!
                for m in 0:ℓ
                    i_alm = almIndex(alms[1], ℓ, m)  # compute alm index
                    for comp in 1:ncomp  # copy over the random variates into buffer
                        alm_in[comp] = alms[comp].alm[i_alm]
                    end
                    lmul!(LowerTriangular(h𝐂'), alm_in)  # transform
                    for comp in 1:ncomp  # copy buffer back into the alms
                        alms[comp].alm[i_alm] = alm_in[comp]
                    end
                end
            end
        end
    end
end
function synalm!(Cl::AbstractArray{T,3}, alms::Union{NTuple{N,A}, Vector{A}};
        lmin=0, lmax=-1) where {T, N, A <: Alm}
    synalm!(Random.default_rng(), Cl, alms; lmin=lmin, lmax=lmax)
end

# Healpix parent
# Base.parent(x::HealpixMap) = x.pixels


text:
function InitialGuess(arg0::Vector{jdouble})
    return InitialGuess((Vector{jdouble},), arg0)
end

function get_initial_guess(obj::InitialGuess)
    return jcall(obj, "getInitialGuess", Vector{jdouble}, ())
end


text:
using ScoreDrivenModels, Distributions, BenchmarkTools, Random, Test

include("test/utils.jl")

scaling = 0.0
ω = [0.1, 0.1]
A = [0.2 0; 0 0.2]
B = [0.2 0; 0 0.2]
simulation = simulate_GAS_1_1(LogNormal, scaling, ω, A, B, 1)
gas = Model(1, 1, LogNormal, scaling)
gas.ω = ω
gas.A[1] = A
gas.B[1] = B
@benchmark score_driven_recursion($gas, $simulation)
# BenchmarkTools.Trial:
#   memory estimate:  235.38 KiB
#   allocs estimate:  18
#   --------------
#   minimum time:     457.305 μs (0.00% GC)
#   median time:      467.132 μs (0.00% GC)
#   mean time:        485.643 μs (2.39% GC)
#   maximum time:     52.405 ms (99.08% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1

scaling = 0.5
ω = [0.1, 0.1]
A = [0.2 0; 0 0.2]
B = [0.2 0; 0 0.2]
simulation = simulate_GAS_1_1(LogNormal, scaling, ω, A, B, 1)
gas = Model(1, 1, LogNormal, scaling)
gas.ω = ω
gas.A[1] = A
gas.B[1] = B
@benchmark score_driven_recursion($gas, $simulation)
# BenchmarkTools.Trial:
#   memory estimate:  938.50 KiB
#   allocs estimate:  25018
#   --------------
#   minimum time:     2.018 ms (0.00% GC)
#   median time:      2.045 ms (0.00% GC)
#   mean time:        2.193 ms (3.19% GC)
#   maximum time:     5.779 ms (55.53% GC)
#   --------------
#   samples:          2278
#   evals/sample:     1

scaling = 1.0
ω = [0.1, 0.1]
A = [0.2 0; 0 0.2]
B = [0.2 0; 0 0.2]
simulation = simulate_GAS_1_1(LogNormal, scaling, ω, A, B, 1)
gas = Model(1, 1, LogNormal, scaling)
gas.ω = ω
gas.A[1] = A
gas.B[1] = B
@benchmark score_driven_recursion($gas, $simulation)
# BenchmarkTools.Trial:
#   memory estimate:  391.63 KiB
#   allocs estimate:  5018
#   --------------
#   minimum time:     1.450 ms (0.00% GC)
#   median time:      1.474 ms (0.00% GC)
#   mean time:        1.538 ms (1.69% GC)
#   maximum time:     5.374 ms (58.97% GC)
#   --------------
#   samples:          3247
#   evals/sample:     1


text:
import Distributions.Distribution
import Distributions.rand
import Distributions.support
import StatsBase.quantile
import StatsBase.mean
import StatsBase.median
import StatsBase.middle
import StatsBase.quantile
import StatsBase.std

import Distributions.pdf

import Base: minimum, maximum

abstract type AbstractEmpiricalDistribution end

struct FittedDistribution{D <: Distribution} <: AbstractEmpiricalDistribution
    distribution::D
end

Broadcast.broadcastable(fd::FittedDistribution) = Ref(fd.distribution)

Distributions.rand(fd::FittedDistribution) = rand(fd.distribution)
Distributions.rand(fd::FittedDistribution, n::Int) = rand(fd.distribution, n)
Distributions.support(fd::FittedDistribution) = support(fd.distribution)
Distributions.pdf(fd::FittedDistribution, x) = pdf(fd.distribution, x)
StatsBase.mean(fd::FittedDistribution) = mean(fd.distribution)
StatsBase.median(fd::FittedDistribution) = median(fd.distribution)
StatsBase.middle(fd::FittedDistribution) = middle(fd.distribution)
StatsBase.quantile(fd::FittedDistribution, q) = quantile(fd.distribution, q)
StatsBase.std(fd::FittedDistribution) = std(fd.distribution)

Base.minimum(uv::FittedDistribution) = minimum(uv.distribution)
Base.maximum(uv::FittedDistribution) = maximum(uv.distribution)

export
AbstractEmpiricalDistribution,
FittedDistribution


text:
module LatexTables

using Parameters, Formatting

export Cell, setstyle!, Format, Style, Color, CellColor,
        apply, Table, Tabular, Cells, Rows, Cols, All, table_to_tex, make_tex_file

abstract type LatexStyle end
abstract type Indices end

mutable struct Cell{T}
    val::T
    style::Dict{Symbol, LatexStyle}

    Cell(val::T, style) where T = new{T}(val, convert(Dict{Symbol, LatexStyle}, style))
end

Cell(val) = Cell(val, Dict())
Cell(val::Real) = Cell(val, Dict(:format => Format("3f")))
Cell(val::String) = Cell(val, Dict(:format => Format('s')))

Base.show(io::IO, c::Cell) = show(io, c.val)

const CellMatrix = AbstractMatrix{<:Cell}
const CellArray = AbstractArray{<:Cell}

include("cellstyles.jl")
include("utilities.jl")

function create_table_top!(rows;
               position::String = "!ht",
               caption::String = "",
               label::String = "",
               centering::Bool = true,
               caption_position_top::Bool = true)

    push!(rows, string("\\begin{table}[", position,"] \n"))
    centering && push!(rows, "\\centering \n")
    (caption_position_top && !isempty(caption)) && push!(rows, string("\\caption{", caption, "} \n"))
    (caption_position_top && !isempty(label)) && push!(rows, string("\\label{", label, "} \n"))
end

function find_first_lrc(s::String)

    ii = findfirst.(['l', 'r', 'c'], s)
    minimum(ii[isa.(ii, Number)])
end

function create_table_mid!(rows, body_cells;
                           leading_col = [],
                           header = [],
                           alignment = [],
                           table_type = :booktabs)

    n_row, n_col = size(body_cells)

    if !isempty(leading_col)
        length(leading_col) != n_row && throw(ArgumentError("Wrong leading column size"))
        n_col_mod = n_col + 1
    else
        n_col_mod = n_col
    end

    if isempty(alignment)
        align = "l"^n_col_mod
    elseif length(alignment) == 1 && isempty(leading_col)
        align = alignment^n_col_mod
    elseif length(alignment) == 1 && !isempty(leading_col)
        align = "l" * alignment^(n_col_mod-1)
    elseif length(alignment) == n_col && !isempty(leading_col)
        ii    = find_first_lrc(alignment)
        align = alignment[1:ii-1] * "l" * alignment[ii:end]
    elseif length(alignment) == n_col_mod
        align = alignment
    else
        throw(ArgumentError("Wrong alignment size"))
    end

    if isempty(header)
        header_cells = header
    elseif isa(header, AbstractVector) && length(header) == n_col_mod
        header_cells = Cell.(reshape(string.(header), 1, :))
    elseif isa(header, AbstractVector) && length(header) == n_col && !isempty(leading_col)
        header_cells = Cell.(hcat("", reshape(string.(header), 1, :)))
    elseif isa(header, AbstractMatrix) && size(header,2) == n_col_mod
        header_cells = Cell.(string.(header))
    elseif isa(header, AbstractMatrix) && size(header,2) == n_col && !isempty(leading_col)
        header_cells = Cell.(hcat(repeat([""],size(header,1),1), string.(header)))
    else
        throw(ArgumentError("Wrong header size"))
    end

    if table_type == :booktabs
        push!(rows, string("\\begin{tabular}{@{} ", align...," @{}} \n"))
        push!(rows, "\\toprule \n")
        if !isempty(header_cells)
            push!(rows, apply(header_cells))
            table_type == :booktabs && push!(rows, "\\midrule \n")
        end
        for i in 1:n_row
            add_row = apply(body_cells[i,:])
            isempty(leading_col) ? push!(rows, add_row) : push!(rows, string(leading_col[i]) * " & " * add_row)
        end
        push!(rows, "\\bottomrule \n")
    elseif table_type == :tabular
        push!(rows, string("\\begin{tabular}{@{} ", align...," @{}} \n"))
        push!(rows, "\\hline \n")
        if !isempty(header_cells)
            push!(rows, apply(header_cells))
            table_type == :booktabs && push!(rows, "\\hline \n")
        end
        for i in 1:n_row
            add_row = apply(body_cells[i,:])
            isempty(leading_col) ? push!(rows, add_row) : push!(rows, string(leading_col[i]) * " & " * add_row)
        end
        push!(rows, "\\hline \n")
    else
        throw(ArgumentError("Wrong table type"))
    end
    push!(rows, "\\end{tabular} \n")
end

function create_table_bot!(rows;
                           caption::String = "",
                           label::String = "",
                           caption_position_top::Bool = false)

    (!caption_position_top && !isempty(caption)) && push!(rows, string("\\caption{", caption, "} \n"))
    (!caption_position_top && !isempty(label)) && push!(rows, string("\\label{", label, "} \n"))
    push!(rows, "\\end{table} \n")
end

function table_to_tex(body::AbstractMatrix;
                      col_format = [],
                      row_format = [],
                      table_type = :booktabs,
                      leading_col = [],
                      header = [],
                      position = "!ht",
                      caption = "",
                      label = "",
                      centering = true,
                      alignment = "",
                      floating_table = true,
                      caption_position_top = true,
                      highlight_max_row = false,
                      highlight_min_row = false,
                      highlight_max_col = false,
                      highlight_min_col = false,
                      highlight_max_style = [CellColor(:green!50), Color(:blue), Style(:bold)],
                      highlight_min_style = Color(:orange),
                      other_rules = [],
                  )

    sum(highlight_max_col + highlight_max_row) > 1 && throw(ArgumentError("Both highlight_max_col and highlight_max_row are true"))
    sum(highlight_min_col + highlight_min_row) > 1 && throw(ArgumentError("Both highlight_min_col and highlight_min_row are true"))
    sum(!isempty(col_format) + !isempty(row_format)) > 1 && throw(ArgumentError("Too many formats specified"))

    rules = []

    n_row, n_col = size(body)

    add_rules_col_format!(rules, col_format, n_col)
    add_rules_row_format!(rules, row_format, n_row)

    highlight_max_col && add_rules_max_col!(rules, body, highlight_max_style)
    highlight_min_col && add_rules_min_col!(rules, body, highlight_min_style)
    highlight_max_row && add_rules_max_row!(rules, body, highlight_max_style)
    highlight_min_row && add_rules_min_row!(rules, body, highlight_min_style)

    [push!(rules, rule) for rule in other_rules]

    body_cells = Cell.(body)

    for (inds, style) in rules
        setstyle!(body_cells, inds, style...)
    end

    rows = String[]
    floating_table && create_table_top!(rows; position=position, caption=caption, label=label, centering=centering, caption_position_top=caption_position_top)
    create_table_mid!(rows, body_cells; leading_col=leading_col, header=header, alignment=alignment, table_type=table_type)
    floating_table && create_table_bot!(rows; caption=caption, label=label, caption_position_top=caption_position_top)

    return reduce(string, rows)
end

function add_rules_col_format!(rules, col_format, n_col)

    if !isempty(col_format)
        col_format = repeat_n(col_format, n_col)
        for i in 1:n_col
            if !isempty(col_format[i])
                push!(rules, Cols(i) => (Format(col_format[i]),))
            end
        end
    end
end

function add_rules_row_format!(rules, row_format, n_row)

    if !isempty(row_format)
        row_format = repeat_n(row_format, n_row)
        for i in 1:n_row
            if !isempty(row_format[i])
                push!(rules, Rows(i) => (Format(row_format[i]),))
            end
        end
    end
end

add_rules_max_col!(rules, body, highlight_style::AbstractVector) = add_rules_max_col!.([rules], [body], highlight_style)
add_rules_min_col!(rules, body, highlight_style::AbstractVector) = add_rules_min_col!.([rules], [body], highlight_style)
add_rules_max_row!(rules, body, highlight_style::AbstractVector) = add_rules_max_row!.([rules], [body], highlight_style)
add_rules_min_row!(rules, body, highlight_style::AbstractVector) = add_rules_min_row!.([rules], [body], highlight_style)
add_rules_max_col!(rules, body, highlight_style::LatexStyle) = add_rules_minmax_col!(rules, body, highlight_style, true)
add_rules_min_col!(rules, body, highlight_style::LatexStyle) = add_rules_minmax_col!(rules, body, highlight_style, false)
add_rules_max_row!(rules, body, highlight_style::LatexStyle) = add_rules_minmax_row!(rules, body, highlight_style, true)
add_rules_min_row!(rules, body, highlight_style::LatexStyle) = add_rules_minmax_row!(rules, body, highlight_style, false)

function add_rules_minmax_col!(rules, body, highlight_style::LatexStyle, find_max::Bool)

    for i in 1:size(body,2)
        i_red = findall(isa.(body[:,i], Number))
        if !isempty(i_red)
            x_ext = find_max ? maximum(body[i_red,i]) : minimum(body[i_red,i])
            i_ext = findall(body[i_red,i] .== x_ext)
            for j in 1:length(i_ext)
                push!(rules, Cells((i_red[i_ext[j]],i)) => (highlight_style,))
            end
        end
    end
end

function add_rules_minmax_row!(rules, body, highlight_style::LatexStyle, find_max::Bool)

    for i in 1:size(body,1)
        i_red = findall(isa.(body[i,:], Number))
        if !isempty(i_red)
            x_ext = find_max ? maximum(body[i,i_red]) : minimum(body[i,i_red])
            i_ext = findall(body[i,i_red] .== x_ext)
            for j in 1:length(i_ext)
                push!(rules, Cells((i,i_red[i_ext[j]])) => (highlight_style,))
            end
        end
    end
end

repeat_n(x::String, n) = repeat([x], n)
repeat_n(x::AbstractVector, n) = x

end


text:
using ApproxFun, Base.Test

## Avazzadev et al

# Example 1

x=Fun([0.,1.])
Q=gamma(0.5)*LeftIntegral(0.5)
f=(2/105*sqrt(x)*(105-56x^2+48x^3))
u=Q\f
@test norm(u-(x^3-x^2+1))<100eps()

# Example 2

x=Fun([0.,1.])
Q=gamma(0.5)*LeftIntegral(0.5)
u=Q\(exp(x)-1)
@test norm(u-exp(x)*erf(sqrt(x))/sqrt(π)) < 100eps() # 5.0036177384681187e-14

# Example 3
x=Fun([0.,1.])
Q=gamma(1/5)*LeftIntegral(1/5)
u=Q\(x+1)
@test norm(u-(1+1.25x)*sin(0.8π)/(π*x^(1/5))) < 10eps()

# Example 4
x=Fun([0.,1.])
Q=gamma(1-1/3)*LeftIntegral(1-1/3)
u=Q\x^(7/6)
@test norm(u-7*gamma(1/6)/(18*sqrt(π)*gamma(2/3))*sqrt(x)) < 100eps()

# Example 5

d=Interval([0.,1.])
x=Fun(d)
f=x+4/3*x^(3/2)
S=Legendre(d)⊕JacobiWeight(.5,0.,Jacobi(.5,.5,d))
Q=gamma(.5)*LeftIntegral(S,.5)

@test_approx_eq sum(f/sqrt(1.-x)) last(Q*f)

L=I+Q

@test_approx_eq last(L.ops[2]*f) last(Q*f)

@test_approx_eq last(L*f) last(f)+last(Q*f)

u=L\f
@test norm(u-x)  < 10eps()

# Example 6

d=Interval([0.,1.])
x=Fun(d)
f=x^2+16/15*x^(5/2)
S=Legendre(d)⊕JacobiWeight(.5,0.,Jacobi(.5,.5,d))
Q=gamma(.5)*LeftIntegral(S,.5)
L=I+Q
u=L\f
@test norm(u-x^2) < 10eps()

# Example 7

d=Interval([0.,1.])
x=Fun(d)
f=2sqrt(x)
S=Legendre(d)⊕JacobiWeight(.5,0.,Jacobi(.5,.5,d))
Q=gamma(.5)*LeftIntegral(S,.5)
L=I+Q
u=L\f

@test norm(1-exp(π*x)*erfc(sqrt(π*x))-u) < 100eps()

# Example 8

d=Interval([0.,1.])
x=Fun(d)
f=1/(x+1)+2*Fun(x->asinh(sqrt(x))/sqrt(1+x),JacobiWeight(.5,0.,d))
S=Legendre(d)⊕JacobiWeight(.5,0.,Jacobi(.5,.5,d))
Q=gamma(.5)*LeftIntegral(S,.5)
L=I+Q
u=L\f
norm((u-1/(x+1)).coefficients) < 1000eps()   # 1.2011889731154679e-14

## Test for bug

QL = LeftIntegral(0.5)	: Legendre()					↦	JacobiWeight(0.5,0.,Ultraspherical(1))
QU = LeftIntegral(0.5)	: JacobiWeight(0.5,0.,Ultraspherical(1))	↦	Legendre()

λ=0.25
A=[λ*I QU; QL λ*I]
L=ApproxFun.interlace(A)
@test_approx_eq L[2,5] 0.


text:
# following Joe Fowler's drift_correct from mass

using KernelDensity, Distributions, Optim

abstract PerPulseTransformation

immutable LinearDriftCorrect <: PerPulseTransformation
    indicator_median::Float64
    param::Float64
end
Base.write(h5::HDF5.HDF5File, path::ASCIIString, dc::LinearDriftCorrect) = h5[path] = Float64[dc.indicator_median,dc.param]
uvhist(edges, y) = UnivariateKDE(edges, convert(Array{Float64,1},hist(y,edges)[2]))
density(h::UnivariateKDE) = h.density
edges(h::UnivariateKDE) = h.x
Base.midpoints(h::UnivariateKDE) = midpoints(edges(h))
Base.conv(h::Histogram, y) = conv(UnivariateKDE(edges(h), convert(Vector{Float64},counts(h))), y)
"""drift_correct(indicator, uncorrected, σ_smooth, limit=NaN;sloperange=(-0.01, 0.01))
Compute a drift correction that minimizes the spectral entropy of
`uncorrected` (a filtered pulse height vector) after correction.
The `indicator` vector should be of equal length and should have some
linear correlation with the pulse gain. Generally it will be the
pretrigger mean of the pulses, but you can experiment with other choices.

The entropy will be computed on corrected values only in the range
[0, `limit`], so `limit` should be set to a characteristic large
value of `uncorrected`. If `limit` is NaN (the default), then
in will be compute as 25% larger than the 99%ile point of
`uncorrected`

The model is that the filtered pulse height PH should be scaled by
(1 + param*PTM) where `param` is an arbitrary parameter computed here, and
PTM is the difference between each record's pretrigger mean and the
median value of all pretrigger means. (Or replace pretrigger mean
with whatever quantity you passed in as `indicator`.)

`σ_smooth` is used when calculating the entropy.
`first(sloperange)` and `last(sloperange)` are the extreme low and high values
in which an optimum value for `param` is searched for."""
function drift_correct(indicator, uncorrected, σ_smooth, limit=NaN;sloperange=(-0.01, 0.01))
    median_indicator = median(indicator)
    centered_indicator = indicator-median_indicator
    if isnan(limit)
        limit = 1.25*StatsBase.percentile(uncorrected, 99)
    end
    edges = pow2edges(σ_smooth*0.4, (0, limit))
    function tooptimize_entropy(param)
        corrected = uncorrected.*(1+centered_indicator*param[1])
        uvh = uvhist(edges, corrected)
        entropy(uvh,σ_smooth)
    end
    optimum = optimize(tooptimize_entropy,first(sloperange),last(sloperange)) # optimize(f(x), lower_bound_x, upper_bound_x) minimizes f with respect to x using brent method
    #check for sucess?
    LinearDriftCorrect(median_indicator, optimum.minimum)
end

"""pow2edges(stepsize_guess, limits)
Return a range `first(limits):s:last(limits)`. Stepsize `s` is chosen to such that the range has an even number
of midpoints, and `s` is near `stepsize_guess`."""
function pow2edges(stepsize_guess, limits)
    # calculate a range from first(limits) to last(limits), stepsize is chose to be close to guess but given a power of 2 number of midpoints
    nbins_guess = round(Int, 0.5+(last(limits)-first(limits))/stepsize_guess )
    pow2=1024
    while pow2<nbins_guess
        pow2*=2
    end
    nbins = pow2
    first(limits):(last(limits)-first(limits))/pow2:last(limits)
end
"""entropy(h::UnivariateKDE, σ)
Calculate the entropy in spectrum `h` after first smoothing by convolution with a gaussian with std devaition σ."""
function entropy(h::UnivariateKDE, σ)
    hsmooth=conv(h, Normal(0,σ))
    return -sum(StatsBase.xlogx(density(hsmooth)))
end

"""applycorrection(correction::PerPulseTransformation, indicator, uncorrected)
Aplly the correction stored in `correction`, for example a gain correction based on pretrigger median, to
convert `uncorrected` to `corrected` based on the values in `indicator`."""
function applycorrection(correction::LinearDriftCorrect, indicator, uncorrected)
    centered_indicator = indicator-correction.indicator_median
    corrected = uncorrected.*(1+centered_indicator*correction.param)
end


text:
module Interface

using Compat; import Compat.String

import ..Legacy:

    AbstractEntry,
    Comment,
    Docs,
    Entry,
    Manual,
    Metadata,
    Page

import ..Utilities: Utilities
import ..Cache
import ..Formats

include("legacy.jl")

end