text:
const _lark_defaults = Dict([
        "debug"=> false,
        "keep_all_tokens"=> false,
        "postlex"=> nothing,
        "parser"=> "lalr",
        "lexer"=> "auto",
        "transformer"=> nothing,
        "start"=> "start",
        "priority"=> "auto",
        "ambiguity"=> "auto",
        "propagate_positions"=> false,
        "lexer_callbacks"=> Dict(),
        "maybe_placeholders"=> false,
        "edit_terminals"=> nothing,
        "g_regex_flags"=> 0,
        "import_paths"=> [],
        "source_path"=> nothing,
    ])

"""
    LarkOptions(options_dict)

A LarkOptions instance contains options for running Lark. Options can be accessed both
as a traditional dictionary and as properties using dot notation. The possible options are
described under the Lark constructor.

"""
struct LarkOptions
    options_dict::Dict
    LarkOptions(o::Dict) = begin
        options = Dict()
        for (name,default) in _lark_defaults
            if name in keys(o)
                value = pop!(o,name)
                if default isa Bool
                    value = Bool(value)
                end
            else value = default
            end
            options[name] = value
        end
        if typeof(options["start"]) <: AbstractString
            options["start"] = [options["start"]]
        end
        if !(options["parser"] in ("earley", "lalr", "cyk", nothing))
            error("Option error: Parser must be earley, lalr, cyk or nothing, given $(options["parser"])")
        end
        if options["parser"] == "earley" && !(options["transformer"]==nothing)
            error("Cannot specify an embedded transformer when using the Earley algorithm." *
                  "Please use your transformer on the resulting parse tree, or use a different algorithm (i.e. lalr)")
        end
        new(options)
    end
end

Base.getproperty(lo::LarkOptions,name::Symbol) = begin
    pdict = getfield(lo,:options_dict)
    if String(name) in keys(pdict)
        pdict[String(name)]
    else
        getfield(lo,name)
    end
end

Base.setproperty!(lo::LarkOptions,name::Symbol,value) = begin
    pdict = getfield(lo,:options_dict)
    if String(name) in keys(pdict)
        pdict[String(name)] = value
    else
        throw(KeyError(name))
    end
end

struct Lark
    options::LarkOptions
    source_path
    source_grammar
    terminals
    rules
    ignore_tokens
    lexer_conf
    parser
end

const VALID_PRIORITY_OPTIONS = ("auto", "normal", "invert", nothing)
const VALID_AMBIGUITY_OPTIONS = ("auto", "resolve", "explicit", "forest")

"""
    Lark(grammar::String;options...)

Create a Lark parser based on the Lark grammar expressed in `grammar`. Options are as
follows:

  * `parser` Decides which parser engine to use, "earley" or "lalr". (Default: "lalr")
         Note: "lalr" requires a lexer. "earley" is not currently implemented.

  * `lexer` Decides whether or not to use a lexer stage

     + "standard": Use a standard lexer

     + "contextual": Stronger lexer (only works with parser="lalr")

     + "auto" (default): Choose for me based on grammar and parser

  * `transformer` - Applies the transformer to every parse tree

  * `debug` - Affects verbosity (default: False)

  * `keep_all_tokens` - Don't automagically remove "punctuation" tokens (default: False)

  * `postlex` - Lexer post-processing (Requires standard lexer. Default: None)

  * `start` - The start symbol (Default: start)

  * `propagate_positions` - Propagates [line, column, end line, end column] attributes into all tree branches.

  * `lexer_callbacks` - Dictionary of callbacks for the lexer. May alter tokens during lexing. Use with caution.
"""
Lark(grammar::String;options...) = begin
    source = "<string>"
    cache_file = "larkcache__$(hash(grammar)%(2^32))"
    Lark(grammar,Dict{String,Any}((String(k),v) for (k,v) in options),source,cache_file)
end

"""
    Lark(grammar::IOStream,source;options...)

Create a Lark parser based on the Lark grammar found by reading `grammar`, whose source
is recorded as `source`. See above for options.
"""
Lark(grammar::IOStream,source;options...) = begin
    cache_file = "larkcache_$(basename(source))"
    textgrammar = read(grammar,String)
    Lark(textgrammar,Dict{String,Any}((String(k),v) for (k,v) in options),source,cache_file)
end

"""
       Lark(grammar::String,loptions,source,cache_file)

Create a Lark parser based on the Lark grammar contained in `grammar`, using options
found in dictionary `loptions`, recording the source as `source` and storing the parser information
in `cache_file`. `cache_file` is currently ignored.
"""
Lark(grammar::String,loptions,source,cache_file) = begin
    options = LarkOptions(loptions)
    if options.lexer == "auto"
        if options.parser == "lalr"
            options.lexer = "contextual"
        elseif options.parser == "earley"
            options.lexer = "dynamic"
        elseif options.parser == "cyk"
            options.lexer = "standard"
        else
            @assert !options.parser
        end
    end

    lexer = options.lexer
    @assert lexer in ("standard", "contextual", "dynamic", "dynamic_complete") || (lexer isa Lexer)

    if options.ambiguity == "auto"
        if options.parser == "earley"
            options.ambiguity = "resolve"
        end
    else
        disambig_parsers = ["earley", "cyk"]
        @assert (options.parser in disambig_parsers)
        "Only $(join(", ",disambig_parser)) supports disambiguation right now"
    end

    if !(options.priority in VALID_PRIORITY_OPTIONS)
        throw(ValueError("invalid priority option: $(options.priority). Must be one of $VALID_PRIORITY_OPTIONS"))
    end

    if !(options.ambiguity in VALID_AMBIGUITY_OPTIONS)
        throw(ValueError("invalid ambiguity option: $(options.ambiguity). Must be one of $VALID_AMBIGUITY_OPTIONS"))
    end

    grammar = load_grammar(grammar, source, options.import_paths, options.keep_all_tokens)

    if options.postlex !== nothing
        terminals_to_keep = Set(options.postlex.always_accept)
    else
        terminals_to_keep = Set()
    end

    # Compile the EBNF grammar into BNF
    try
        terminals, rules, ignore_tokens = compile(grammar,options.start,terminals_to_keep)
        if options.edit_terminals !== nothing
            for t in terminals
                options.edit_terminals(t)
            end
        end
        _terminals_dict = Dict([(t.name => t) for t in terminals])

        if options.priority == "invert"
            for rule in rules
                if rule.options.priority !== nothing
                    rule.options.priority = -rule.options.priority
                end
            end
        elseif options.priority ===  nothing
            for rule in rules
                if rule.options.priority !== nothing
                    rule.options.priority = nothing
                end
            end
        end
        lexer_callbacks = Dict() #options.transformer !== nothing ? _get_lexer_callbacks(options.transformer,terminals) : Dict()
        merge!(lexer_callbacks,options.lexer_callbacks)
        lexer_conf = LexerConf(terminals,ignore=ignore_tokens,postlex=options.postlex,
                               callbacks=lexer_callbacks,
                               g_regex_flags=options.g_regex_flags)
        ## Note that lexer is not preserved so will do nothing.
        if options.parser !== nothing
            parser = _build_parser(options,rules,lexer_conf)
        elseif lexer != nothing
            lexer = TraditionalLexer(lexer_conf.tokens,lexer_conf.ignore,lexer_conf.callbacks)
        end
        Lark(options,source,grammar,terminals,rules,ignore_tokens,lexer_conf,parser)
    catch ex
        if ex isa TaskFailedException
            rethrow(ex.task.exception)
        else
            rethrow(ex)
        end
    end
end

#TODO: integrate transformer calls into create_callback
_prepare_callbacks(options,rules) = begin
    parser_class = get_frontend(options.parser,options.lexer)
    _callbacks = nothing
    if options.ambiguity != "forest"
        _parse_tree_builder = ParseTreeBuilder(rules,
                                               propagate_positions = options.propagate_positions,
                                               ambiguous = options.parser != "lalr" && options.ambiguity=="explicit",
                                               maybe_placeholders=options.maybe_placeholders)
        _callbacks = create_callback(_parse_tree_builder) #,transformer=options.transformer)
    end
    return parser_class,_callbacks
end

_build_parser(options,rules,lexer_conf) = begin
    parser_class,callbacks = _prepare_callbacks(options,rules)
    parser_conf = ParserConf(rules,callbacks,options.start)
    return parser_class(lexer_conf,parser_conf,options=options)
end

"""
    open(grammar_filename;rel_to=nothing,options...)

Open `grammar_filename` and read in a Lark grammar, returning a Lark parser. If
`rel_to` is not `nothing`, `grammar_filename` is relative to `rel_to`. Options
are the same as for [`Lark(grammar::String;options...)`](@ref).
"""
open(grammar_filename;rel_to=nothing,options...) = begin
    if rel_to != nothing
        basepath = dirname(rel_to)
        grammar_filename = joinpath(basepath,grammar_filename)
    end
    Lark(Base.open(grammar_filename),grammar_filename,options...)
end

show(io::IOStream,l::Lark) = begin
    print(io,"Lark(open($(l.source)), parser = $(l.options["parser"]), lexer = $(l.options["lexer"]), ...)")
end

lex(l::Lark, text) = begin
    stream = lex(l.lexer,text)
    if l.options.postlex
        return process(l.options.postlex,stream)
    end
    return stream
end

get_terminal(l::Lark,name) = l._terminals_dict[name]

"""
    parse(l::Lark,text;start=nothing,on_error=nothing)

Parse `text` using parser `l`, returning a parse tree. If `start` is
not nothing, it identifies the start symbol. `on_error` is currently
ignored.

If `text` does not conform to the grammar, either
[`UnexpectedToken`](@ref) or [`UnexpectedCharacters`](@ref)
is raised.
"""
parse(l::Lark,text;start=nothing,on_error=nothing) = begin
    try
        parse(l.parser,text,start=start)
    catch ex
        if ex isa TaskFailedException
            rethrow(ex.task.exception)
        else
            rethrow(ex)
        end
    end
    # TODO Add lark 0.11.0 fancy error handling
end


text:
using Revise
using AdFem
using LinearAlgebra
using PyPlot
using DelimitedFiles
using Optim
using MAT

function calc_residual_and_jacobian(θ, u)
    E = θ[1] * constant(ones(get_ngauss(mmesh)))
    nu = 0.3 * constant(ones(get_ngauss(mmesh)))
    μ = E / (2*(1 + nu))
    λ = E*nu/((1 + nu)*(1 - 2*nu))

    ψ, J = neo_hookean(u, μ, λ, mmesh)
    B = compute_fem_source_term(zeros(get_ngauss(mmesh)), -0.5*ones(get_ngauss(mmesh)), mmesh)
    ψ = ψ - B
    J, ψ = impose_Dirichlet_boundary_conditions(J, ψ, bdnode, zeros(length(bdnode)))
    ψ, J
end

mmesh = Mesh(joinpath(PDATA, "twoholes_large.stl"))
mmesh = Mesh(mmesh.nodes * 10, mmesh.elems)

θ = Variable([1.0; 0.3])
bdnode = bcnode((x,y)->(x<1e-5) || (x>0.49-1e-5), mmesh)
bdnode = [bdnode; bdnode .+ mmesh.ndof]
nr = newton_raphson_with_grad(calc_residual_and_jacobian, zeros(2mmesh.ndof), θ)

U = matread("fenics/fwd1.mat")["x"]

loss = sum((nr - U)^2)
g = tf.convert_to_tensor(gradients(loss, θ))

sess = Session(); init(sess)

BFGS!(sess, loss, g, θ)

# close("all")
# gradview(sess, θ, loss, [5.0])
# savefig("gradtest.png")

ext: jl


text:
# Julia wrapper for header: /usr/local/include/libavcodec/avcodec.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_destruct_packet_nofree,
    av_destruct_packet,
    av_init_packet,
    av_new_packet,
    av_shrink_packet,
    av_grow_packet,
    av_dup_packet,
    av_free_packet,
    av_packet_new_side_data,
    av_packet_get_side_data,
    av_audio_resample_init,
    audio_resample,
    audio_resample_close,
    av_resample_init,
    av_resample,
    av_resample_compensate,
    av_resample_close,
    avpicture_alloc,
    avpicture_free,
    avpicture_fill,
    avpicture_layout,
    avpicture_get_size,
    avcodec_get_chroma_sub_sample,
    avcodec_get_pix_fmt_name,
    avcodec_set_dimensions,
    avcodec_pix_fmt_to_codec_tag,
    av_get_codec_tag_string,
    avcodec_get_pix_fmt_loss,
    avcodec_find_best_pix_fmt,
    img_get_alpha_info,
    avpicture_deinterlace,
    av_codec_next,
    #avcodec_version,
    avcodec_configuration,
    avcodec_license,
    avcodec_init,
    avcodec_register,
    avcodec_find_encoder,
    avcodec_find_encoder_by_name,
    avcodec_find_decoder,
    avcodec_find_decoder_by_name,
    avcodec_string,
    av_get_profile_name,
    avcodec_get_context_defaults,
    avcodec_get_context_defaults2,
    avcodec_get_context_defaults3,
    avcodec_alloc_context,
    avcodec_alloc_context2,
    avcodec_alloc_context3,
    avcodec_copy_context,
    avcodec_get_frame_defaults,
    avcodec_alloc_frame,
    avcodec_default_get_buffer,
    avcodec_default_release_buffer,
    avcodec_default_reget_buffer,
    avcodec_get_edge_width,
    avcodec_align_dimensions,
    avcodec_align_dimensions2,
    avcodec_default_get_format,
    avcodec_thread_init,
    avcodec_default_execute,
    avcodec_default_execute2,
    avcodec_open,
    avcodec_open2,
    avcodec_decode_audio3,
    avcodec_decode_audio4,
    avcodec_decode_video2,
    avcodec_decode_subtitle2,
    avsubtitle_free,
    avcodec_encode_audio,
    avcodec_encode_audio2,
    avcodec_fill_audio_frame,
    avcodec_encode_video,
    avcodec_encode_subtitle,
    avcodec_close,
    avcodec_register_all,
    avcodec_flush_buffers,
    avcodec_default_free_buffers,
    av_get_pict_type_char,
    av_get_bits_per_sample,
    av_get_bits_per_sample_format,
    av_parser_next,
    av_register_codec_parser,
    av_parser_init,
    av_parser_parse2,
    av_parser_change,
    av_parser_close,
    av_register_bitstream_filter,
    av_bitstream_filter_init,
    av_bitstream_filter_filter,
    av_bitstream_filter_close,
    av_bitstream_filter_next,
    av_fast_realloc,
    av_fast_malloc,
    av_picture_copy,
    av_picture_crop,
    av_picture_pad,
    av_xiphlacing,
    av_log_missing_feature,
    av_register_hwaccel,
    av_hwaccel_next,
    av_lockmgr_register,
    avcodec_get_type,
    avcodec_get_class,
    avcodec_is_open

function av_destruct_packet_nofree(pkt)
    ccall((:av_destruct_packet_nofree,libavcodec),Void,(Ptr{AVPacket},),pkt)
end

function av_destruct_packet(pkt)
    ccall((:av_destruct_packet,libavcodec),Void,(Ptr{AVPacket},),pkt)
end

function av_init_packet(pkt)
    ccall((:av_init_packet,libavcodec),Void,(Ptr{AVPacket},),pkt)
end

function av_new_packet(pkt,size::Integer)
    ccall((:av_new_packet,libavcodec),Cint,(Ptr{AVPacket},Cint),pkt,size)
end

function av_shrink_packet(pkt,size::Integer)
    ccall((:av_shrink_packet,libavcodec),Void,(Ptr{AVPacket},Cint),pkt,size)
end

function av_grow_packet(pkt,grow_by::Integer)
    ccall((:av_grow_packet,libavcodec),Cint,(Ptr{AVPacket},Cint),pkt,grow_by)
end

function av_dup_packet(pkt)
    ccall((:av_dup_packet,libavcodec),Cint,(Ptr{AVPacket},),pkt)
end

function av_free_packet(pkt)
    ccall((:av_free_packet,libavcodec),Void,(Ptr{AVPacket},),pkt)
end

function av_packet_new_side_data(pkt,_type::AVPacketSideDataType,size::Integer)
    ccall((:av_packet_new_side_data,libavcodec),Ptr{Uint8},(Ptr{AVPacket},AVPacketSideDataType,Cint),pkt,_type,size)
end

function av_packet_get_side_data(pkt,_type::AVPacketSideDataType,size)
    ccall((:av_packet_get_side_data,libavcodec),Ptr{Uint8},(Ptr{AVPacket},AVPacketSideDataType,Ptr{Cint}),pkt,_type,size)
end

function av_audio_resample_init(output_channels::Integer,input_channels::Integer,output_rate::Integer,input_rate::Integer,sample_fmt_out::AVSampleFormat,sample_fmt_in::AVSampleFormat,filter_length::Integer,log2_phase_count::Integer,linear::Integer,cutoff::Cdouble)
    ccall((:av_audio_resample_init,libavcodec),Ptr{ReSampleContext},(Cint,Cint,Cint,Cint,AVSampleFormat,AVSampleFormat,Cint,Cint,Cint,Cdouble),output_channels,input_channels,output_rate,input_rate,sample_fmt_out,sample_fmt_in,filter_length,log2_phase_count,linear,cutoff)
end

function audio_resample(s,output,input,nb_samples::Integer)
    ccall((:audio_resample,libavcodec),Cint,(Ptr{ReSampleContext},Ptr{Int16},Ptr{Int16},Cint),s,output,input,nb_samples)
end

function audio_resample_close(s)
    ccall((:audio_resample_close,libavcodec),Void,(Ptr{ReSampleContext},),s)
end

function av_resample_init(out_rate::Integer,in_rate::Integer,filter_length::Integer,log2_phase_count::Integer,linear::Integer,cutoff::Cdouble)
    ccall((:av_resample_init,libavcodec),Ptr{AVResampleContext},(Cint,Cint,Cint,Cint,Cint,Cdouble),out_rate,in_rate,filter_length,log2_phase_count,linear,cutoff)
end

function av_resample(c,dst,src,consumed,src_size::Integer,dst_size::Integer,update_ctx::Integer)
    ccall((:av_resample,libavcodec),Cint,(Ptr{AVResampleContext},Ptr{Int16},Ptr{Int16},Ptr{Cint},Cint,Cint,Cint),c,dst,src,consumed,src_size,dst_size,update_ctx)
end

function av_resample_compensate(c,sample_delta::Integer,compensation_distance::Integer)
    ccall((:av_resample_compensate,libavcodec),Void,(Ptr{AVResampleContext},Cint,Cint),c,sample_delta,compensation_distance)
end

function av_resample_close(c)
    ccall((:av_resample_close,libavcodec),Void,(Ptr{AVResampleContext},),c)
end

function avpicture_alloc(picture,pix_fmt::PixelFormat,width::Integer,height::Integer)
    ccall((:avpicture_alloc,libavcodec),Cint,(Ptr{AVPicture},PixelFormat,Cint,Cint),picture,pix_fmt,width,height)
end

function avpicture_free(picture)
    ccall((:avpicture_free,libavcodec),Void,(Ptr{AVPicture},),picture)
end

function avpicture_fill(picture,ptr,pix_fmt::PixelFormat,width::Integer,height::Integer)
    ccall((:avpicture_fill,libavcodec),Cint,(Ptr{AVPicture},Ptr{Uint8},PixelFormat,Cint,Cint),picture,ptr,pix_fmt,width,height)
end

function avpicture_layout(src,pix_fmt::PixelFormat,width::Integer,height::Integer,dest,dest_size::Integer)
    ccall((:avpicture_layout,libavcodec),Cint,(Ptr{AVPicture},PixelFormat,Cint,Cint,Ptr{Cuchar},Cint),src,pix_fmt,width,height,dest,dest_size)
end

function avpicture_get_size(pix_fmt::PixelFormat,width::Integer,height::Integer)
    ccall((:avpicture_get_size,libavcodec),Cint,(PixelFormat,Cint,Cint),pix_fmt,width,height)
end

function avcodec_get_chroma_sub_sample(pix_fmt::PixelFormat,h_shift,v_shift)
    ccall((:avcodec_get_chroma_sub_sample,libavcodec),Void,(PixelFormat,Ptr{Cint},Ptr{Cint}),pix_fmt,h_shift,v_shift)
end

function avcodec_get_pix_fmt_name(pix_fmt::PixelFormat)
    ccall((:avcodec_get_pix_fmt_name,libavcodec),Ptr{Uint8},(PixelFormat,),pix_fmt)
end

function avcodec_set_dimensions(s,width::Integer,height::Integer)
    ccall((:avcodec_set_dimensions,libavcodec),Void,(Ptr{AVCodecContext},Cint,Cint),s,width,height)
end

function avcodec_pix_fmt_to_codec_tag(pix_fmt::PixelFormat)
    ccall((:avcodec_pix_fmt_to_codec_tag,libavcodec),Uint32,(PixelFormat,),pix_fmt)
end

function av_get_codec_tag_string(buf,buf_size::Csize_t,codec_tag::Integer)
    ccall((:av_get_codec_tag_string,libavcodec),Csize_t,(Ptr{Uint8},Csize_t,Uint32),buf,buf_size,codec_tag)
end

function avcodec_get_pix_fmt_loss(dst_pix_fmt::PixelFormat,src_pix_fmt::PixelFormat,has_alpha::Integer)
    ccall((:avcodec_get_pix_fmt_loss,libavcodec),Cint,(PixelFormat,PixelFormat,Cint),dst_pix_fmt,src_pix_fmt,has_alpha)
end

function avcodec_find_best_pix_fmt(pix_fmt_mask::Int64,src_pix_fmt::PixelFormat,has_alpha::Integer,loss_ptr)
    ccall((:avcodec_find_best_pix_fmt,libavcodec),Cint,(Int64,PixelFormat,Cint,Ptr{Cint}),pix_fmt_mask,src_pix_fmt,has_alpha,loss_ptr)
end

function img_get_alpha_info(src,pix_fmt::PixelFormat,width::Integer,height::Integer)
    ccall((:img_get_alpha_info,libavcodec),Cint,(Ptr{AVPicture},PixelFormat,Cint,Cint),src,pix_fmt,width,height)
end

function avpicture_deinterlace(dst,src,pix_fmt::PixelFormat,width::Integer,height::Integer)
    ccall((:avpicture_deinterlace,libavcodec),Cint,(Ptr{AVPicture},Ptr{AVPicture},PixelFormat,Cint,Cint),dst,src,pix_fmt,width,height)
end

function av_codec_next(c)
    ccall((:av_codec_next,libavcodec),Ptr{AVCodec},(Ptr{AVCodec},),c)
end

function avcodec_version()
    ccall((:avcodec_version,libavcodec),Uint32,())
end

function avcodec_configuration()
    ccall((:avcodec_configuration,libavcodec),Ptr{Uint8},())
end

function avcodec_license()
    ccall((:avcodec_license,libavcodec),Ptr{Uint8},())
end

function avcodec_init()
    ccall((:avcodec_init,libavcodec),Void,())
end

function avcodec_register(codec)
    ccall((:avcodec_register,libavcodec),Void,(Ptr{AVCodec},),codec)
end

function avcodec_find_encoder(id::CodecID)
    ccall((:avcodec_find_encoder,libavcodec),Ptr{AVCodec},(CodecID,),id)
end

function avcodec_find_encoder_by_name(name)
    ccall((:avcodec_find_encoder_by_name,libavcodec),Ptr{AVCodec},(Ptr{Uint8},),name)
end

function avcodec_find_decoder(id::CodecID)
    ccall((:avcodec_find_decoder,libavcodec),Ptr{AVCodec},(CodecID,),id)
end

function avcodec_find_decoder_by_name(name)
    ccall((:avcodec_find_decoder_by_name,libavcodec),Ptr{AVCodec},(Ptr{Uint8},),name)
end

function avcodec_string(buf,buf_size::Integer,enc,encode::Integer)
    ccall((:avcodec_string,libavcodec),Void,(Ptr{Uint8},Cint,Ptr{AVCodecContext},Cint),buf,buf_size,enc,encode)
end

function av_get_profile_name(codec,profile::Integer)
    ccall((:av_get_profile_name,libavcodec),Ptr{Uint8},(Ptr{AVCodec},Cint),codec,profile)
end

function avcodec_get_context_defaults(s)
    ccall((:avcodec_get_context_defaults,libavcodec),Void,(Ptr{AVCodecContext},),s)
end

function avcodec_get_context_defaults2(s,arg1::AVMediaType)
    ccall((:avcodec_get_context_defaults2,libavcodec),Void,(Ptr{AVCodecContext},AVMediaType),s,arg1)
end

function avcodec_get_context_defaults3(s,codec)
    ccall((:avcodec_get_context_defaults3,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{AVCodec}),s,codec)
end

function avcodec_alloc_context()
    ccall((:avcodec_alloc_context,libavcodec),Ptr{AVCodecContext},())
end

function avcodec_alloc_context2(arg1::AVMediaType)
    ccall((:avcodec_alloc_context2,libavcodec),Ptr{AVCodecContext},(AVMediaType,),arg1)
end

function avcodec_alloc_context3(codec)
    ccall((:avcodec_alloc_context3,libavcodec),Ptr{AVCodecContext},(Ptr{AVCodec},),codec)
end

function avcodec_copy_context(dest,src)
    ccall((:avcodec_copy_context,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{AVCodecContext}),dest,src)
end

function avcodec_get_frame_defaults(pic)
    ccall((:avcodec_get_frame_defaults,libavcodec),Void,(Ptr{AVFrame},),pic)
end

function avcodec_alloc_frame()
    ccall((:avcodec_alloc_frame,libavcodec),Ptr{AVFrame},())
end

function avcodec_default_get_buffer(s,pic)
    ccall((:avcodec_default_get_buffer,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{AVFrame}),s,pic)
end

function avcodec_default_release_buffer(s,pic)
    ccall((:avcodec_default_release_buffer,libavcodec),Void,(Ptr{AVCodecContext},Ptr{AVFrame}),s,pic)
end

function avcodec_default_reget_buffer(s,pic)
    ccall((:avcodec_default_reget_buffer,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{AVFrame}),s,pic)
end

function avcodec_get_edge_width()
    ccall((:avcodec_get_edge_width,libavcodec),Uint32,())
end

function avcodec_align_dimensions(s,width,height)
    ccall((:avcodec_align_dimensions,libavcodec),Void,(Ptr{AVCodecContext},Ptr{Cint},Ptr{Cint}),s,width,height)
end

function avcodec_align_dimensions2(s,width,height,linesize_align)
    ccall((:avcodec_align_dimensions2,libavcodec),Void,(Ptr{AVCodecContext},Ptr{Cint},Ptr{Cint},Ptr{Cint}),s,width,height,linesize_align)
end

function avcodec_default_get_format(s,fmt)
    ccall((:avcodec_default_get_format,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{PixelFormat}),s,fmt)
end

function avcodec_thread_init(s,thread_count::Integer)
    ccall((:avcodec_thread_init,libavcodec),Cint,(Ptr{AVCodecContext},Cint),s,thread_count)
end

function avcodec_default_execute(c,func,arg,ret,count::Integer,size::Integer)
    ccall((:avcodec_default_execute,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{Void},Ptr{Void},Ptr{Cint},Cint,Cint),c,func,arg,ret,count,size)
end

function avcodec_default_execute2(c,func,arg,ret,count::Integer)
    ccall((:avcodec_default_execute2,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{Void},Ptr{Void},Ptr{Cint},Cint),c,func,arg,ret,count)
end

function avcodec_open(avctx,codec)
    ccall((:avcodec_open,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{AVCodec}),avctx,codec)
end

function avcodec_open2(avctx,codec,options)
    ccall((:avcodec_open2,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{AVCodec},Ptr{Ptr{AVDictionary}}),avctx,codec,options)
end

function avcodec_decode_audio3(avctx,samples,frame_size_ptr,avpkt)
    ccall((:avcodec_decode_audio3,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{Int16},Ptr{Cint},Ptr{AVPacket}),avctx,samples,frame_size_ptr,avpkt)
end

function avcodec_decode_audio4(avctx,frame,got_frame_ptr,avpkt)
    ccall((:avcodec_decode_audio4,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{AVFrame},Ptr{Cint},Ptr{AVPacket}),avctx,frame,got_frame_ptr,avpkt)
end

function avcodec_decode_video2(avctx,picture,got_picture_ptr,avpkt)
    ccall((:avcodec_decode_video2,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{AVFrame},Ptr{Cint},Ptr{AVPacket}),avctx,picture,got_picture_ptr,avpkt)
end

function avcodec_decode_subtitle2(avctx,sub,got_sub_ptr,avpkt)
    ccall((:avcodec_decode_subtitle2,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{AVSubtitle},Ptr{Cint},Ptr{AVPacket}),avctx,sub,got_sub_ptr,avpkt)
end

function avsubtitle_free(sub)
    ccall((:avsubtitle_free,libavcodec),Void,(Ptr{AVSubtitle},),sub)
end

function avcodec_encode_audio(avctx,buf,buf_size::Integer,samples)
    ccall((:avcodec_encode_audio,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{Uint8},Cint,Ptr{Int16}),avctx,buf,buf_size,samples)
end

function avcodec_encode_audio2(avctx,avpkt,frame,got_packet_ptr)
    ccall((:avcodec_encode_audio2,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{AVPacket},Ptr{AVFrame},Ptr{Cint}),avctx,avpkt,frame,got_packet_ptr)
end

function avcodec_fill_audio_frame(frame,nb_channels::Integer,sample_fmt::AVSampleFormat,buf,buf_size::Integer,align::Integer)
    ccall((:avcodec_fill_audio_frame,libavcodec),Cint,(Ptr{AVFrame},Cint,AVSampleFormat,Ptr{Uint8},Cint,Cint),frame,nb_channels,sample_fmt,buf,buf_size,align)
end

function avcodec_encode_video(avctx,buf,buf_size::Integer,pict)
    ccall((:avcodec_encode_video,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{Uint8},Cint,Ptr{AVFrame}),avctx,buf,buf_size,pict)
end

function avcodec_encode_subtitle(avctx,buf,buf_size::Integer,sub)
    ccall((:avcodec_encode_subtitle,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{Uint8},Cint,Ptr{AVSubtitle}),avctx,buf,buf_size,sub)
end

function avcodec_close(avctx)
    ccall((:avcodec_close,libavcodec),Cint,(Ptr{AVCodecContext},),avctx)
end

function avcodec_register_all()
    ccall((:avcodec_register_all,libavcodec),Void,())
end

function avcodec_flush_buffers(avctx)
    ccall((:avcodec_flush_buffers,libavcodec),Void,(Ptr{AVCodecContext},),avctx)
end

function avcodec_default_free_buffers(s)
    ccall((:avcodec_default_free_buffers,libavcodec),Void,(Ptr{AVCodecContext},),s)
end

function av_get_pict_type_char(pict_type::Integer)
    ccall((:av_get_pict_type_char,libavcodec),Uint8,(Cint,),pict_type)
end

function av_get_bits_per_sample(codec_id::CodecID)
    ccall((:av_get_bits_per_sample,libavcodec),Cint,(CodecID,),codec_id)
end

function av_get_bits_per_sample_format(sample_fmt::AVSampleFormat)
    ccall((:av_get_bits_per_sample_format,libavcodec),Cint,(AVSampleFormat,),sample_fmt)
end

function av_parser_next(c)
    ccall((:av_parser_next,libavcodec),Ptr{AVCodecParser},(Ptr{AVCodecParser},),c)
end

function av_register_codec_parser(parser)
    ccall((:av_register_codec_parser,libavcodec),Void,(Ptr{AVCodecParser},),parser)
end

function av_parser_init(codec_id::Integer)
    ccall((:av_parser_init,libavcodec),Ptr{AVCodecParserContext},(Cint,),codec_id)
end

function av_parser_parse2(s,avctx,poutbuf,poutbuf_size,buf,buf_size::Integer,pts::Int64,dts::Int64,pos::Int64)
    ccall((:av_parser_parse2,libavcodec),Cint,(Ptr{AVCodecParserContext},Ptr{AVCodecContext},Ptr{Ptr{Uint8}},Ptr{Cint},Ptr{Uint8},Cint,Int64,Int64,Int64),s,avctx,poutbuf,poutbuf_size,buf,buf_size,pts,dts,pos)
end

function av_parser_change(s,avctx,poutbuf,poutbuf_size,buf,buf_size::Integer,keyframe::Integer)
    ccall((:av_parser_change,libavcodec),Cint,(Ptr{AVCodecParserContext},Ptr{AVCodecContext},Ptr{Ptr{Uint8}},Ptr{Cint},Ptr{Uint8},Cint,Cint),s,avctx,poutbuf,poutbuf_size,buf,buf_size,keyframe)
end

function av_parser_close(s)
    ccall((:av_parser_close,libavcodec),Void,(Ptr{AVCodecParserContext},),s)
end

function av_register_bitstream_filter(bsf)
    ccall((:av_register_bitstream_filter,libavcodec),Void,(Ptr{AVBitStreamFilter},),bsf)
end

function av_bitstream_filter_init(name)
    ccall((:av_bitstream_filter_init,libavcodec),Ptr{AVBitStreamFilterContext},(Ptr{Uint8},),name)
end

function av_bitstream_filter_filter(bsfc,avctx,args,poutbuf,poutbuf_size,buf,buf_size::Integer,keyframe::Integer)
    ccall((:av_bitstream_filter_filter,libavcodec),Cint,(Ptr{AVBitStreamFilterContext},Ptr{AVCodecContext},Ptr{Uint8},Ptr{Ptr{Uint8}},Ptr{Cint},Ptr{Uint8},Cint,Cint),bsfc,avctx,args,poutbuf,poutbuf_size,buf,buf_size,keyframe)
end

function av_bitstream_filter_close(bsf)
    ccall((:av_bitstream_filter_close,libavcodec),Void,(Ptr{AVBitStreamFilterContext},),bsf)
end

function av_bitstream_filter_next(f)
    ccall((:av_bitstream_filter_next,libavcodec),Ptr{AVBitStreamFilter},(Ptr{AVBitStreamFilter},),f)
end

function av_fast_realloc(ptr,size,min_size::Csize_t)
    ccall((:av_fast_realloc,libavcodec),Ptr{Void},(Ptr{Void},Ptr{Uint32},Csize_t),ptr,size,min_size)
end

function av_fast_malloc(ptr,size,min_size::Csize_t)
    ccall((:av_fast_malloc,libavcodec),Void,(Ptr{Void},Ptr{Uint32},Csize_t),ptr,size,min_size)
end

function av_picture_copy(dst,src,pix_fmt::PixelFormat,width::Integer,height::Integer)
    ccall((:av_picture_copy,libavcodec),Void,(Ptr{AVPicture},Ptr{AVPicture},PixelFormat,Cint,Cint),dst,src,pix_fmt,width,height)
end

function av_picture_crop(dst,src,pix_fmt::PixelFormat,top_band::Integer,left_band::Integer)
    ccall((:av_picture_crop,libavcodec),Cint,(Ptr{AVPicture},Ptr{AVPicture},PixelFormat,Cint,Cint),dst,src,pix_fmt,top_band,left_band)
end

function av_picture_pad(dst,src,height::Integer,width::Integer,pix_fmt::PixelFormat,padtop::Integer,padbottom::Integer,padleft::Integer,padright::Integer,color)
    ccall((:av_picture_pad,libavcodec),Cint,(Ptr{AVPicture},Ptr{AVPicture},Cint,Cint,PixelFormat,Cint,Cint,Cint,Cint,Ptr{Cint}),dst,src,height,width,pix_fmt,padtop,padbottom,padleft,padright,color)
end

function av_xiphlacing(s,v::Integer)
    ccall((:av_xiphlacing,libavcodec),Uint32,(Ptr{Cuchar},Uint32),s,v)
end

function av_log_missing_feature(avc,feature,want_sample::Integer)
    ccall((:av_log_missing_feature,libavcodec),Void,(Ptr{Void},Ptr{Uint8},Cint),avc,feature,want_sample)
end

function av_register_hwaccel(hwaccel)
    ccall((:av_register_hwaccel,libavcodec),Void,(Ptr{AVHWAccel},),hwaccel)
end

function av_hwaccel_next(hwaccel)
    ccall((:av_hwaccel_next,libavcodec),Ptr{AVHWAccel},(Ptr{AVHWAccel},),hwaccel)
end

function av_lockmgr_register(cb)
    ccall((:av_lockmgr_register,libavcodec),Cint,(Ptr{Void},),cb)
end

function avcodec_get_type(codec_id::CodecID)
    ccall((:avcodec_get_type,libavcodec),Cint,(CodecID,),codec_id)
end

function avcodec_get_class()
    ccall((:avcodec_get_class,libavcodec),Ptr{AVClass},())
end

function avcodec_is_open(s)
    ccall((:avcodec_is_open,libavcodec),Cint,(Ptr{AVCodecContext},),s)
end


text:
let
  # testing constructor and pdf
  d = BoolDistribution(0.3)
  @test pdf(d, true) == 0.3
  @test pdf(d, false) == 0.7

  # testing support
  @test support(d) == [true, false]

  # testing ==
  d2 = BoolDistribution(0.3)
  @test d == d2
  @test BoolDistribution(0.4) != BoolDistribution(0.1)

  # testing hash
  @test hash(d) == hash(d2)

  @test sprint((io,d)->show(io,MIME("text/plain"),d), d) == sprint((io,d)->showdistribution(io,d,title="BoolDistribution"), d)
end


text:
writing_output = false
@everywhere Random.seed!(42)

if VERSION < v"1.5"
    ver = "111"
else
    ver = "150"
end

weights = rand(400)
weights = weights ./ sum(weights)

test_sys_resample    = SMC.resample(weights, method = :systematic)
test_multi_resample  = SMC.resample(weights, method = :multinomial)
test_poly_resample   = SMC.resample(weights, method = :polyalgo)

saved_filename = string("reference/resample_version=", ver, ".jld2")
if writing_output
    jldopen(saved_filename, true, true, true, IOStream) do file
        write(file, "sys", test_sys_resample)
        write(file, "multi", test_multi_resample)
        write(file, "poly", test_poly_resample)
    end
end

saved_sys_resample   = load(saved_filename, "sys")
saved_multi_resample = load(saved_filename, "multi")
saved_poly_resample  = load(saved_filename, "poly")

####################################################################

@testset "Resampling methods" begin
    @test test_sys_resample   == saved_sys_resample
    @test test_multi_resample == saved_multi_resample
    @test test_poly_resample  == saved_poly_resample
end


text:
"""
```
transform_transition_matrices(m::AbstractDSGEModel, T::Matrix{T},
                              R::Matrix{T}, C::Matrix{T}) where {T<:AbstractFloat}
```

### Arguments

-`m`: the m object
-`T`, `R`, and `C`: matrices of the state transition equation

### Return values

- `TTT`, `RRR`, and `CCC`: transform transition equation matrices to accomodate states simulated
    between observations

### Description
"""
function transform_transition_matrices(m::AbstractDSGEModel, TT::Matrix{T},
                                       R::Matrix{T}, C::Vector{T};
                                       track_lag::Bool = true) where {T<:AbstractFloat}

    freq = get_setting(m, :state_simulation_freq)
    if track_lag
        TTT = zeros(size(TT, 1) * (freq+1), size(TT, 2) * (freq+1)) # pre-allocate
        RRR = zeros(size(R, 1) * (freq+1), size(R, 2) * freq) # pre-allocate
    else
        TTT = zeros(size(TT, 1) * freq, size(TT, 2) * freq) # pre-allocate
        RRR = zeros(size(R, 1) * freq, size(R, 2) * freq) # pre-allocate
    end

    # Construct powers of TT iteratively
    TT_powers = Dict{Int64, Matrix{T}}()
    for i = 1:freq
        if i == 1
            TT_powers[i] = TT
        else
            TT_powers[i] = TT_powers[i - 1] * TT
        end
    end

    # Insert directly into TTT
    if track_lag
        for i = 1:freq+1
            if i == 1
                TTT[1:size(TT,1), 1 + size(TT,2) * freq:end] = eye(TT) # track last lag
            else
                TTT[1 + size(TT,1)*(i-1):size(TT,1)*i, 1 + size(TT,2)*freq:end] = TT_powers[i-1]
            end
        end
    else
        for i = 1:freq
            TTT[1 + size(TT,1)*(i-1):size(TT,1)*i, 1 + size(TT,2)*(freq-1):end] = TT_powers[i]
        end
    end

    # Create RRR
    TR_powers = Dict{Int64, Matrix{T}}()
    TR_powers[0] = R

    # Construct T^i * R iteratively
    for i = 1:(freq - 1)
        TR_powers[i] = TT * TR_powers[i - 1]
    end

    # Fill in RRR and CCC
    if track_lag
        RRR[1+size(R,1):size(R, 1)*2, 1:size(R,2)] = TR_powers[0]
        for i in 3:freq+1
            for j = (i-2):-1:0
                RRR[1+size(R,1)*(i-1):size(R,1)*i, 1+size(R,2)*j:size(R,2)*(j+1)] = TR_powers[i-j-2]
            end
        end
        CCC = vcat(zeros(C), repeat(C, freq))
    else
        RRR[1:size(R, 1), 1:size(R,2)] = TR_powers[0]
        for i in 2:freq
            for j = (i-1):-1:0
                # inserts values from right nonzero blocks to left
                RRR[1+size(R,1)*(i-1):size(R,1)*i, 1+size(R,2)*j:size(R,2)*(j+1)] = TR_powers[i-j-1]
            end
        end
        CCC = repeat(C, freq) #vec(repeat(C, 1, freq))
    end

    return TTT, RRR, CCC
end


text:
#!/usr/bin/env julia
using Test

try
  using BenchmarkTools
catch
  import Pkg
  Pkg.add("BenchmarkTools")
  Pkg.update("BenchmarkTools")
  using BenchmarkTools
end

if length(ARGS) < 1
  N = 1000000
else
  N = parse(Int, ARGS[1]);
end

println("--> Julia $VERSION    N: $N")

function g(N::Int)
  s = 0.
  for k = 1:N+1
    s = s + (-1.)^(k+1) / (2*k-1)
  end

  x::Real=4.0*s
return x
end

@test isapprox(g(N), pi, atol=1e-4)
o = @benchmark g(N)
println(minimum(o).time/1e9, " seconds")

max_stars_repo_name: scienceopen/numba-examples

max_issues_repo_name: scienceopen/numba-examples

max_forks_repo_name: scienceopen/numba-examples


text:
# technically, this is a product kernel of a continuous and discrete input kernel ...
struct SeparableKernel{T, K, M<:AbstractMatOrFac{T}} <: MultiKernel{T}
    k::K # kernel for input covariances
    B::M # matrix for output covariances
    # potentially initialize with different type
end
function SeparableKernel(B::AbstractMatrix, k)
    T, K, M = promote_type(eltype(B), eltype(k)), typeof(k), typeof(B)
    SeparableKernel{T, K, M}(k, B)
end
const Separable = SeparableKernel
# outputs matrix-valued covariance
(k::Separable)(x, y) = k.B * k.k(x, y)

Base.size(K::Separable) = size(K.B)
Base.getindex(K::Separable, i::Integer, j::Integer) = K.B[i,j] * K.k

# NOTE: maybe this is not worth it for small d!
function evaluate_block!(Gij::AbstractMatrix, k::SeparableKernel, x, y, T = input_trait(k))
    kxy = k.k(x, y)
    @. Gij = B * kxy
end

# this is isotopic marginalization, where we observe all outputs for all inputs
# function gramian(k::Separable, x::AbstractVector, y::AbstractVector = x)
#     k.B ⊗ gramian(k.k, x, y)
# end

function factorize(G::Gramian{<:AbstractMatrix, <:Separable})
    factorize(kronecker(G))
end

function KroneckerProducts.kronecker(G::Gramian{<:Any, <:Separable})
    gramian(G.k.k, G.x, G.y) ⊗ G.k.B
end

# also, need special multiplication
function LinearAlgebra.mul!(y::AbstractVecOfVec, G::Gramian{<:AbstractMatrix, <:Separable},
                            x::AbstractVecOfVec)
    K = gramian(G.k.k, G.x, G.y) ⊗ G.k.B # convert to kronecker # order?
    return mul!(y, kronecker(G), x)
end

# simply applies kernel in d dimensions, useful to implement T-SNE acceleration via matrix-valued Barnes-Hut
# struct DiagonalKernel{T, KT} <: MultiKernel{T}
#     k::KT
#     d::Int
# end
#
# (D::DiagonalKernel)(x, y) = (D.k(x, y)*I)(d)
#
# function evaluate_block!(Gij::Diagonal, k::DiagonalKernel, x, y, T = input_trait(k))
#     Gij.diag .= k.k(x, y)
# end


text:
# test_generic_domain.jl

# We test the generic functionality of a domain.
# These tests check whether the given domain correctly implements the
# interface of a domain.
function test_generic_domain(d::Domain)
    @test eltype(eltype(d)) == subeltype(d)
    @test isreal(d) == isreal(subeltype(d))

    if !isempty(d)
        x = point_in_domain(d)
        @test x ∈ d
    else
        try
            x = point_in_domain(d)
            @test false
        catch
        end
    end
end

@testset "generic domains" begin
    domains = [
        0..1,
        UnitInterval()^3
    ]

    @testset "$(rpad("Generic domains",80))" begin
        for domain in domains
            test_generic_domain(domain)
        end
    end
end


text:
# BBI R-Tree
# ==========

const RTREE_MAGIC = 0x2468ACE0

# Supplemental Table 14.
struct RTreeHeader
    magic::UInt32
    block_size::UInt32
    item_count::UInt64
    lo::NTuple{2,UInt32}
    up::NTuple{2,UInt32}
    end_file_offset::UInt64
    items_per_slot::UInt32
    reserved::UInt32
end

const RTREE_HEADER_SIZE = 48

function Base.read(stream::IO, ::Type{RTreeHeader})
    return RTreeHeader(
        read(stream, UInt32),
        read(stream, UInt32),
        read(stream, UInt64),
        readbound(stream),
        readbound(stream),
        read(stream, UInt64),
        read(stream, UInt32),
        read(stream, UInt32))
end

function Base.write(stream::IO, header::RTreeHeader)
    return write(
        stream,
        header.magic,
        header.block_size,
        header.item_count,
        header.lo[1], header.lo[2],
        header.up[1], header.up[2],
        header.end_file_offset,
        header.items_per_slot,
        header.reserved)
end

# Supplemental Table 15.
struct RTreeNodeFormat
    isleaf::UInt8
    reserved::UInt8
    count::UInt16
end

const RTREE_NODE_FORMAT_SIZE = 4

function Base.write(io::IO, node::RTreeNodeFormat)
    return write(io, node.isleaf, node.reserved, node.count)
end

function isleaf(node::RTreeNodeFormat)
    return node.isleaf != 0x00
end

function Base.read(io::IO, ::Type{RTreeNodeFormat})
    return RTreeNodeFormat(read(io, UInt8), read(io, UInt8), read(io, UInt16))
end

# Supplemental Table 16.
struct RTreeLeafNode
    lo::NTuple{2,UInt32}
    up::NTuple{2,UInt32}
    offset::UInt64
    size::UInt64
end

const RTREE_LEAF_NODE_SIZE = 32

function Base.write(stream::IO, node::RTreeLeafNode)
    return write(
        stream,
        node.lo[1], node.lo[2],
        node.up[1], node.up[2],
        node.offset, node.size)
end

# Supplemental Table 17
struct RTreeInternalNode
    lo::NTuple{2,UInt32}
    up::NTuple{2,UInt32}
    offset::UInt64
end

const RTREE_INTERNAL_NODE_SIZE = 32

function Base.write(stream::IO, node::RTreeInternalNode)
    return write(
        stream,
        node.lo[1], node.lo[2],
        node.up[1], node.up[2],
        node.offset)
end

# TODO: Merge this with RTreeLeafNode and SectionSummary
# Data block indexed by R-tree (0-origin, left-closed and right-open).
struct Block
    # lower bound of genomic interval (chromid, chromstart)
    lo::Tuple{UInt32,UInt32}

    # upper bound of genomic interval (chromid, chromend)
    up::Tuple{UInt32,UInt32}

    # offset to data block
    offset::UInt64

    # size of data block (likely compressed, disk dize)
    size::UInt64
end

# disk-serialized R tree
struct RTree{T<:IO}
    stream::T
    offset::UInt64
    header::RTreeHeader
end

function RTree(stream::IO, offset::Integer)
    seek(stream, offset)
    header = read(stream, RTreeHeader)
    if header.magic != RTREE_MAGIC
        error("invalid R tree magic bytes")
    end
    return RTree(stream, offset, header)
end

# Find overlapping blocks.
function find_overlapping_blocks(tree::RTree, chromid::UInt32, chromstart::UInt32, chromend::UInt32)
    root_offset = tree.offset + (RTREE_HEADER_SIZE % UInt64)
    stack = [root_offset]
    ret = Block[]
    while !isempty(stack)
        offset = pop!(stack)
        seek(tree.stream, offset)
        node = read(tree.stream, RTreeNodeFormat)
        for i in 1:node.count
            lo = readbound(tree.stream)
            up = readbound(tree.stream)
            if overlaps(chromid, chromstart, chromend, lo, up)
                offset = read(tree.stream, UInt64)
                if isleaf(node)
                    datasize = read(tree.stream, UInt64)
                    push!(ret, Block(lo, up, offset, datasize))
                else
                    push!(stack, offset)
                end
            else
                skip(tree.stream, isleaf(node) ? 16 : 8)
            end
        end
    end
    # Leaves may be already sorted?
    sort!(ret, by=n->n.offset)
    return ret
end

function overlaps(chromid, chromstart, chromend, lo, up)
    if chromid < lo[1] || (chromid == lo[1] && chromend ≤ lo[2])
        # query is strictly left
        return false
    end

    if chromid > up[1] || (chromid == up[1] && chromstart ≥ up[2])
        # query is strictly right
        return false
    end

    return true
end

struct InMemoryRTree
    lo::NTuple{2,UInt32}
    up::NTuple{2,UInt32}
    offset::UInt64
    size::UInt64
    children::Vector{InMemoryRTree}

    function InMemoryRTree(lo, up, offset, size, children=nothing)
        if children == nothing
            return new(lo, up, offset, size)
        end

        return new(lo, up, offset, size, children)
    end
end

function isleaf(rtree::InMemoryRTree)
    return rtree.offset == 1
end

function InMemoryRTree(children::Vector{InMemoryRTree})
    @assert !isempty(children)
    # TODO: assert this
    #@assert issorted(n.lo for n in children)
    return InMemoryRTree(
        #children[1].lo,
        Base.minimum(n.lo for n in children),
        Base.maximum(n.up for n in children),
        0, 0,
        children)
end

function write_rtree(stream::IO, blocks::Vector{Block})
    blocksize = 64
    items_per_slot = 1
    root = build_inmemory_rtree(blocks, blocksize)
    header = RTreeHeader(
        RTREE_MAGIC,
        blocksize,
        root.up[1] - root.lo[1] + 1,
        root.lo,
        root.up,
        position(stream),
        items_per_slot,
        # reserved
        0)
    n = write(stream, header)

    # write R-tree recursively
    function rec(node, offset)
        # compute the size of node
        nodesize = RTREE_NODE_FORMAT_SIZE + (isleaf(node) ? RTREE_LEAF_NODE_SIZE : RTREE_INTERNAL_NODE_SIZE) * length(node.children)

        # write children
        offsets = UInt64[]
        datasizes = UInt64[]
        child_offset = offset + nodesize
        for child in node.children
            if isleaf(node)
                push!(offsets, child.offset)
                push!(datasizes, child.size)
            else
                push!(offsets, child_offset)
                child_offset = rec(child, child_offset)
            end
        end

        # write node itself
        truncate(stream, child_offset)  # This is necessary to assert `position(stream) == offset` after seeking.
        seek(stream, offset)
        @assert position(stream) == offset
        write(stream, RTreeNodeFormat(isleaf(node) ? 0x01 : 0x00, 0x00, length(node.children)))
        for i in 1:lastindex(node.children)
            child = node.children[i]
            if isleaf(node)
                write(stream, RTreeLeafNode(child.lo, child.up, offsets[i], datasizes[i]))
            else
                write(stream, RTreeInternalNode(child.lo, child.up, offsets[i]))
            end
        end

        return child_offset
    end

    root_offset = position(stream)
    return n + (rec(root, root_offset) - root_offset)
end

# Build an in-memory B-Tree from leaves to root (bottom-up).
function build_inmemory_rtree(blocks::Vector{Block}, blocksize::Int)
    if !issorted(blocks, by=b->b.lo)
        blocks = sort(blocks, by=b->b.lo)
    end
    function rec(blocks)
        if length(blocks) ≤ blocksize
            # store indexes in leaves
            children = InMemoryRTree[]
            if isempty(blocks)
                lo = up = (UInt32(0), UInt32(0))
            else
                lo = blocks[1].lo
                up = blocks[1].up
                for block in blocks
                    up = max(up, block.up)
                    push!(children, InMemoryRTree(block.lo, block.up, block.offset, block.size))
                end
            end
            return InMemoryRTree(lo, up, 1, 0, children)
        else
            d = cld(length(blocks), blocksize)
            children = InMemoryRTree[]
            for i in 1:blocksize
                idx = (i-1)*d+1:min(i*d,lastindex(blocks))
                push!(children, rec(view(blocks, idx)))
            end
            return InMemoryRTree(children)
        end
    end
    return rec(view(blocks, 1:lastindex(blocks)))
end

function readbound(stream::IO)
    return read(stream, UInt32), read(stream, UInt32)
end


text:
@testset "Bayes" begin
    m=NaiveBayesClassifier([:spam, :ham])
    TextAnalysis.fit!(m, "this is ham", :ham);
    TextAnalysis.fit!(m, "this is spam", :spam);
    r = TextAnalysis.predict(m, "is this spam?")
    @test r[:spam] > r[:ham]

    n = NaiveBayesClassifier([:spam, :ham])
    TextAnalysis.fit!(n, StringDocument("this is ham"), :ham)
    TextAnalysis.fit!(n, StringDocument("this is spam"), :spam)
    p = TextAnalysis.predict(n, TokenDocument("is this spam?"))
    @test p[:spam] > p[:ham]
    @test p[:spam] == r[:spam]
end


text:
@testset "Generating ZhuXian data file..." begin
        home = "D:\\multi-criteria-cws\\data\\other\\zx"
        zx = ZhuXian(home)
        @show zx

        file = KongYiji.dir("zhuxian.jld2")
        @time Knet.save(file, "zx", zx)
        @time zfile = KongYiji.compress(file)
        rm(file)
        @time file = KongYiji.decompress(zfile)
        @time zx2 = Knet.load(file, "zx")
        @test zx == zx2
end


text:
using Images
using Downloads
import BM3DDenoise

σ = [5, 10, 15, 20, 25, 30, 35, 40, 50, 60, 70, 75, 80, 90, 100]
PSNR_data = [38, 35, 33, 32, 31, 30, 29, 28, 27, 26, 25, 25, 25, 24, 23]

# download file
if !isfile(joinpath(tempdir(), "Lena512.png"))
	println("Downloading Lena512.png")
	Downloads.download(
		"http://www.cs.tut.fi/~foi/GCF-BM3D/images/Lena512.png",
		joinpath(tempdir(), "Lena512.png"),
	)
end
if !isfile(joinpath(tempdir(), "image_Lena512rgb_noi_s100.png"))
	println("Downloading image_Lena512rgb_noi_s100.png")
	Downloads.download(
		"http://www.cs.tut.fi/~foi/GCF-BM3D/images/image_Lena512rgb_noi_s100.png",
		joinpath(tempdir(), "image_Lena512rgb_noi_s100.png"),
	)
end
if !isfile(joinpath(tempdir(), "image_Lena512rgb.png"))
	println("Downloading image_Lena512rgb.png")
	Downloads.download(
		"http://www.cs.tut.fi/~foi/GCF-BM3D/images/image_Lena512rgb.png",
		joinpath(tempdir(), "image_Lena512rgb.png"),
	)
end
Threads.@threads for i in σ
	if !isfile(joinpath(tempdir(), "Lena512_noi_s" * string(i) * ".png"))
		println("Downloading Lena512_noi_s" * string(i) * ".png")
		Downloads.download(
			"http://www.cs.tut.fi/~foi/GCF-BM3D/images/Lena512_noi_s" *
			string(i) *
			".png",
			joinpath(tempdir(), "Lena512_noi_s" * string(i) * ".png"),
		)
	end
end

Lena = load(joinpath(tempdir(), "Lena512.png"))
config = BM3DDenoise.bm3d_config()
config.show_info = true
@testset "BM3DDenoise_gray" begin
	for i in 1:length(σ)
		@testset "σ=$(σ[i])" begin
			img = load(
				joinpath(
					tempdir(),
					"Lena512_noi_s" * string(σ[i]) * ".png",
				),
			)
			bm3d_psnr = assess_psnr(
				Lena,
				BM3DDenoise.bm3d(img, σ[i] / 255, config),
			)
			println("σ=$(σ[i]), $(bm3d_psnr) > $(PSNR_data[i])")
			@test bm3d_psnr >= PSNR_data[i]
		end
	end
end

Lena = load(joinpath(tempdir(), "image_Lena512rgb.png"))
@testset "BM3DDenoise_color" begin
	img = load(joinpath(tempdir(), "image_Lena512rgb_noi_s100.png"))
	bm3d_psnr = assess_psnr(Lena, BM3DDenoise.bm3d(img, 100 / 255, config))
	println("σ=100(rgb), $(bm3d_psnr)")
	@test bm3d_psnr >= 22
end


text:
using Pkg
cd( "/home/graham_s/julia/vw/Visualisations/")
Pkg.activate(".")
Pkg.develop( path="/home/graham_s/julia/DashBootstrapComponents/" )
# Pkg.update()
include( "/home/graham_s/julia/vw/Visualisations/src/scotbudg.jl")


text:
"""
温度流的Bubble积分
"""

"""
获取所有的Bubble(pp, fs, nfs, ex, nex)
qpp = k1+k2; qfs = k3-k2; qex = k1-k3
"""
function all_bubble_ec_mt(Γ4::Gamma4, lval)
    #
    lamb = Γ4.λ_0 * exp(-lval)
    brlu_area = area(Γ4.model.brillouin)
    #首先得到需要被积分的等能线
    posimat = Matrix{Vector{Segment}}(undef, Γ4.model.bandnum, Γ4.patchnum)
    negamat = Matrix{Vector{Segment}}(undef, Γ4.model.bandnum, Γ4.patchnum)
    #bidx-能带的指标，nidx-patch的指标
    for idxs in CartesianIndices(posimat)
        (bidx, nidx) = Tuple(idxs)
        posimat[bidx, nidx] = []
        negamat[bidx, nidx] = []
    end
    #
    for bidx in 1:1:Γ4.model.bandnum
        posi, peidx = const_energy_line_in_patches(
            Γ4.ltris, Γ4.ladjs, Γ4.lpats, lamb, Γ4.model.dispersion[bidx]
        )
        #将正能量的线加到对应的位置上
        for (edg, pidx) in zip(posi, peidx)
            push!(posimat[bidx, pidx], edg)
        end
        #
        nega, neidx = const_energy_line_in_patches(
            Γ4.ltris, Γ4.ladjs, Γ4.lpats, -lamb, Γ4.model.dispersion[bidx]
        )
        #将负能量的线加到对应的位置上
        for (edg, pidx) in zip(nega, neidx)
            push!(negamat[bidx, pidx], edg)
        end
    end
    #计算pp散射的bubble
    bubbval = Array{Float64, 7}(
        undef,
        #alpha, beta, b1, b2
        Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum,
        #nidx(alpha), k1(b1), k2(b2)
        Γ4.patchnum, Γ4.patchnum, Γ4.patchnum
    )
    Threads.@threads for idxs in CartesianIndices(bubbval)
        alpha, beta, b1, b2, i_n, i1, i2 = Tuple(idxs)
        k1, k2 = Γ4.patches[b1][i1], Γ4.patches[b2][i2]
        qpp = Γ4.model.kadd(k1, k2)
        bubbres = pi_αβ_minus_ec(
            posimat[alpha, i_n], negamat[alpha, i_n],
            brlu_area, lamb,
            qpp, Γ4.model.dispersion[beta]
        )
        bubbval[alpha, beta, b1, b2, i_n, i1, i2] = bubbres
    end
    bubb_qpp = Bubble{:ec, :minus}(lval, bubbval)
    #计算fs和-fs
    bubbval_fs = Array{Float64, 7}(
        undef,
        #alpha, beta, b2, b3
        Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum,
        #nidx(alpha), k2(b2), k3(b3)
        Γ4.patchnum, Γ4.patchnum, Γ4.patchnum
    )
    bubbval_nfs = Array{Float64, 7}(
        undef,
        #alpha, beta, b2, b3
        Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum,
        #nidx(alpha), k2(b2), k3(b3)
        Γ4.patchnum, Γ4.patchnum, Γ4.patchnum
    )
    Threads.@threads for idxs in CartesianIndices(bubbval_fs)
        alpha, beta, b2, b3, i_n, i2, i3 = Tuple(idxs)
        k2, k3 = Γ4.patches[b2][i2], Γ4.patches[b3][i3]
        qfs = Γ4.model.kadd(k3, -k2)
        nqfs = -qfs
        bubbres_fs = pi_αβ_plus_ec(
            posimat[alpha, i_n], negamat[alpha, i_n],
            brlu_area, lamb,
            qfs, Γ4.model.dispersion[beta]
        )
        bubbres_nfs = pi_αβ_plus_ec(
            posimat[alpha, i_n], negamat[alpha, i_n],
            brlu_area, lamb,
            nqfs, Γ4.model.dispersion[beta]
        )
        bubbval_fs[alpha, beta, b2, b3, i_n, i2, i3] = bubbres_fs
        bubbval_nfs[alpha, beta, b2, b3, i_n, i2, i3] = bubbres_nfs
    end
    bubb_qfs = Bubble{:ec, :plus}(lval, bubbval_fs)
    bubb_nqfs = Bubble{:ec, :plus}(lval, bubbval_nfs)
    #计算ex和-ex
    bubbval_ex = Array{Float64, 7}(
        undef,
        #alpha, beta, b1, b3
        Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum,
        #nidx(alpha), k1(b1), k3(b3)
        Γ4.patchnum, Γ4.patchnum, Γ4.patchnum
    )
    bubbval_nex = Array{Float64, 7}(
        undef,
        #alpha, beta, b1, b3
        Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum,
        #nidx(alpha), k1(b1), k3(b3)
        Γ4.patchnum, Γ4.patchnum, Γ4.patchnum
    )
    Threads.@threads for idxs in CartesianIndices(bubbval_ex)
        alpha, beta, b1, b3, i_n, i1, i3 = Tuple(idxs)
        k1, k3 = Γ4.patches[b1][i1], Γ4.patches[b3][i3]
        qex = Γ4.model.kadd(k1, -k3)
        nqex = -qex
        bubbres_ex = pi_αβ_plus_ec(
            posimat[alpha, i_n], negamat[alpha, i_n],
            brlu_area, lamb,
            qex, Γ4.model.dispersion[beta]
        )
        bubbres_nex = pi_αβ_plus_ec(
            posimat[alpha, i_n], negamat[alpha, i_n],
            brlu_area, lamb,
            nqex, Γ4.model.dispersion[beta]
        )
        bubbval_ex[alpha, beta, b1, b3, i_n, i1, i3] = bubbres_ex
        bubbval_nex[alpha, beta, b1, b3, i_n, i1, i3] = bubbres_nex
    end
    bubb_qex = Bubble{:ec, :plus}(lval, bubbval_ex)
    bubb_nqex = Bubble{:ec, :plus}(lval, bubbval_nex)
    return bubb_qpp, bubb_qfs, bubb_nqfs, bubb_qex, bubb_nqex
end

"""
使用能量cutoff作为flow parameter的bubble\n
posi是能量为+LAMBDA的边，nega是能量为-LAMBDA的边, lamb是LAMBDA\n
disp是色散关系，注意这个需要是beta带的，qval是需要平移的大小，应该用一个Point来包装，\n
kshf是动量相加的函数, 这个函数应该能处理好到第一布里渊区的映射\n
area是第一布里渊区的面积\n
```(10.112)本身已经处理好了动量守恒，k, k-q是需要满足动量守恒的关系的，而处理好```
```k-q到第一布里渊区的映射就处理好了Umklapp```
"""
function pi_αβ_plus_ec(
    posi::Vector{Basics.Segment}, nega::Vector{Basics.Segment},
    area::Float64, lamb::Float64,
    qval::Point2D,
    disp::Function
)
    """
    10.112中的 PI^+(n, q) = +LAMBDA (2pi)^-2 beta^-1 Int_{k in k_n} G'(k)G(k - Q)
    其中有一个beta是频率积分带来的，2pi^2是动量积分带来的
    G(k)=CITA(LAMBDA < abs(disp(k))) / i*omega - disp(k)
    G'(k)=-DELTA(abs(disp(k))-LAMBDA) / i*omege - disp(k)
    在零温的情况下10.112中的频率部分可以积分出来，此后的k都是不包含频率的
    = +LAMBDA (2pi)^-2 Int_{k in k_n} CITA() -DELTA()
    { beta^-1 sum_{omega} [(i*omega-disp(k))(i*omega-disp(k - q))]^-1 }
    花括号中的内容求和完之后等于 - CITA(-disp(k)disp(k-q)) / (abs(disp(k)) + abs(disp(k-p)))
    积分会变成
    = +LAMBDA (2pi)^-2 Int_{k in k_n} DELTA(abs(disp(k))-LAMBDA) CITA(LAMBDA<abs(disp(k-q)))
    CITA(-disp(k)disp(k-q)) / (abs(disp(k)) + abs(disp(k-p)))
    因为采用的能量cutoff中有一个 DELTA(abs(disp(k))-LAMBDA)，disp(k)等于正的或者负的LAMBDA
    而CITA(-disp(k)disp(k-q))限制了disp(k)和disp(k-q)符号相反
    所以上式变成
    (第一项disp(k)=LAMBDA>0，于是disp(k-q)<0，而且abs(disp(k))=-disp(k)>LAMBDA)
    (第二项类似，分子中的abs(disp(k))都可以直接换成LAMBDA，abs(disp(k-q))也都知道符号)
    = +LAMBDA (2pi)^-2 Int_{k in kn} {
        DELTA(disp(k)-LAMBDA)CITA(-disp(k-q)-LAMBDA) / (LAMBDA - disp(k - q))
        DELTA(disp(k)+LAMBDA)CITA(disp(k-q)-LAMBDA) / (LAMBDA + disp(k - q)) }
    还可以从积分里面把DELTA给积分掉，这样对于二维平面的积分也会变成对
    disp(k) = LAMBDA 或者 -LAMBDA的线的积分
    = +LAMBDA (2pi)^-2 *
    [Int_{disp(k) = +LAMBDA} CITA(-disp(k-q)-LAMBDA) / (LAMBDA - disp(k - q))]
    +[Int_{disp(k) = -LAMBDA} CITA(disp(k-q)-LAMBDA) / (LAMBDA + disp(k - q))  ]
    """
    nega_q = -qval
    #积分正LAMBDA的线
    intposi = 0.
    for edg in posi
        kval = middle_point(edg.pt1, edg.pt2)
        kprim = kval + nega_q
        #kprim = kadd(kval, nega_q)
        #CITA
        disp_kprim = disp(kprim.x, kprim.y)
        if -disp_kprim < lamb
            continue
        end
        #线积分，计算线元的长度
        intposi += edg.length / (lamb - disp_kprim)
    end
    #积分负LAMBDA的线
    intnega = 0.
    for edg in nega
        kval = middle_point(edg.pt1, edg.pt2)
        kprim = kval + nega_q
        #kprim = kadd(kval, nega_q)
        #CITA
        disp_kprim = disp(kprim.x, kprim.y)
        if disp_kprim < lamb
            continue
        end
        intnega += edg.length / (lamb + disp_kprim)
    end
    #乘上系数
    result = lamb * (intposi + intnega) / area
    return result
end

"""
使用能量cutoff作为flow parameter的bubble\n
posi是能量为+LAMBDA的边，nega是能量为-LAMBDA的边, lamb是LAMBDA\n
disp是色散关系，注意这个需要是beta带的，qval是需要平移的大小，应该用一个Point来包装，\n
kshf是动量相加的函数, 这个函数应该能处理好到第一布里渊区的映射\n
area是第一布里渊区的面积\n
```(10.112)本身已经处理好了动量守恒，k, k-q是需要满足动量守恒的关系的，而处理好```
```k-q到第一布里渊区的映射就处理好了Umklapp```
因为算能量的时候，超出第一布里渊区也没有关系，所以不再需要kadd
"""
function pi_αβ_minus_ec(
    posi::Vector{Basics.Segment}, nega::Vector{Basics.Segment},
    area::Float64, lamb::Float64,
    qval::Point2D,
    disp::Function
)
    """
    10.112中的 PI^-(n, q) = -LAMBDA (2pi)^-2 beta^-1 Int_{k in k_n} G'(k)G(- k + Q)
    = -LAMBDA (2pi)^-2 Int_{k in k_n} CITA() -DELTA()
    { beta^-1 sum_{omega} [(i*omega-disp(k))(-i*omega-disp(-k + q))]^-1 }
    在零温下这个频率积分等于，注意-k那里把频率也给反过来了
    +CITA(+disp(k)disp(-k+q)) / (abs(disp(k)) + abs(disp(-k+q)))
    原式就等于
    = LAMBDA (2pi)^-2 Int_{k in k_n} {
        DELTA(abs(disp(k))-LAMBDA) CITA(abs(disp(-k+q)-LAMBDA))
        CITA(disp(k)disp(-k+q)) / (abs(disp(k)) + abs(disp(-k+q))) }
    第二个CITA限制了disp(k)和disp(-k+q)同号，积分积掉DELTA，分类讨论正负
    = LAMBDA (2pi)^-2 {
        Int_{disp(k) = +LAMBDA} CITA(disp(-k+q) - LAMBDA) / (LAMBDA + disp(-k+q)) +
        Int_{disp(k) = -LAMBDA} CITA(-disp(-k+q) -LAMBDA) / (LAMBDA - disp(-k+q))
    }
    """
    #积分正LAMBDA的线
    intposi = 0.
    for edg in posi
        kval = middle_point(edg.pt1, edg.pt2)
        kprim = qval - kval
        #kprim = kadd(-kval, qval)
        #CITA
        disp_kprim = disp(kprim.x, kprim.y)
        if disp_kprim < lamb
            continue
        end
        #线积分，计算线元的长度
        intposi += edg.length / (lamb + disp_kprim)
    end
    #积分负LAMBDA的线
    intnega = 0.
    for edg in nega
        kval = middle_point(edg.pt1, edg.pt2)
        kprim = qval - kval
        #kprim = kadd(-kval, qval)
        #CITA
        disp_kprim = disp(kprim.x, kprim.y)
        if -disp_kprim < lamb
            continue
        end
        intnega += edg.length / (lamb - disp_kprim)
    end
    #乘上系数
    result = lamb * (intposi + intnega) / area
    return result
end


text:
function create(::Type{GaussianCurveFitter})
    return jcall(GaussianCurveFitter, "create", GaussianCurveFitter, ())
end

function with_max_iterations(obj::GaussianCurveFitter, arg0::jint)
    return jcall(obj, "withMaxIterations", GaussianCurveFitter, (jint,), arg0)
end

function with_start_point(obj::GaussianCurveFitter, arg0::Vector{jdouble})
    return jcall(obj, "withStartPoint", GaussianCurveFitter, (Vector{jdouble},), arg0)
end


text:
# Extract a column as a Dagger array
export getindexcol, getdatacol, dindex, ddata,
       DColumns, column, columns, rows, pairs, as

import Base: keys, values
import IndexedTables: DimName, Columns, column, columns,
       rows, pairs, as, As, Tup, namedtuple, itable

import Dagger: DomainBlocks, ArrayDomain, DArray,
                ArrayOp, domainchunks, chunks, Distribute

const DColumns = DArray{<:Tup}

function DColumns(arrays::Tup)
    if length(arrays) == 0
        error("""DColumns must be constructed with at least
                 one column.""")
    end

    i = findfirst(x->isa(x, ArrayOp), arrays)
    wrap = isa(arrays, Tuple) ? tuple :
                                namedtuple(fieldnames(arrays)...)
    if i == 0
        error("""At least 1 array passed to
                 DColumns must be a DArray""")
    end

    darrays = asyncmap(arrays) do x
        isa(x, ArrayOp) ? compute(get_context(), x) : x
    end

    dist = domainchunks(darrays[i])
    darrays = map(darrays) do x
        if isa(x, DArray)
            domainchunks(x) == dist ?
                x : error("Distribution incompatible")
        else
            Distribute(dist, x)
        end
    end

    darrays = asyncmap(darrays) do x
        compute(get_context(), x)
    end

    if length(darrays) == 1
        cs = chunks(darrays[1])
        chunkmatrix = reshape(cs, length(cs), 1)
    else
        chunkmatrix = reduce(hcat, map(chunks, darrays))
    end
    cs = mapslices(x -> delayed((c...) -> Columns(wrap(c...)))(x...), chunkmatrix, 2)[:]
    T = isa(arrays, Tuple) ? Tuple{map(eltype, arrays)...} :
        wrap{map(eltype, arrays)...}

    DArray(T, domain(darrays[1]), domainchunks(darrays[1]), cs, (i, x...)->vcat(x...))
end

function itable(keycols::DColumns, valuecols::DColumns)
    cs = map(delayed(itable), chunks(keycols), chunks(valuecols))
    cs1 = compute(get_context(),
                  delayed((xs...) -> [xs...]; meta=true)(cs...))
    fromchunks(cs1)
end

function extractarray(t::Union{DNDSparse,DColumns}, accessor)
    arraymaker = function (cs_tup...)
        cs = [cs_tup...]
        lengths = length.(domain.(cs))
        dmnchunks = DomainBlocks((1,), (cumsum(lengths),))
        T = promote_eltype_chunktypes(cs)
        DArray(T, ArrayDomain(1:sum(lengths)), dmnchunks, [cs...], (i, x...)->vcat(x...))
    end

    cs = map(delayed(accessor), t.chunks)
    compute(delayed(arraymaker; meta=true)(cs...))
end

isas(d) = isa(d, As) && d.f !== identity

function columns(t::Union{DNDSparse, DColumns}, which::Tuple...)
    if !isempty(which) && any(isas, which[1])
        return _columns_as(t, which...)
    end

    cs = map(delayed(x->columns(x, which...)), t.chunks)
    f = delayed() do c
        map(tochunk, c)
    end

    tuples = collect(get_context(), treereduce(delayed(vcat), map(f, cs)))

    if isa(tuples, Tuple)
        tuples = [tuples]
    end

    # tuples is a vector of tuples
    map(tuples...) do cstup...
        cs = [cstup...]
        T = chunktype(cs[1])
        ls = length.(domain.(cs))
        d = ArrayDomain((1:sum(ls),))
        dchunks = DomainBlocks((1,), (cumsum(ls),))
        DArray(eltype(T), d, dchunks, cs, (i, x...) -> vcat(x...))
    end
end

function _columns_as(t, which)
    stripas(w) = isa(w, As) ? w.src : w
    which_ = ntuple(i->as(stripas(which[i]), i), length(which))
    cs = columns(t, which_)
    asvecs = find(isas, which)
    outvecs = Any[cs...]
    outvecs[asvecs] = map((w,x) -> w.f(x), [which[asvecs]...], outvecs[asvecs])
    tup = IndexedTables._output_tuple(which)
    tup(outvecs...)
end

for f in [:rows, :keys, :values]
    @eval function $f(t::Union{DNDSparse, ArrayOp}, which::Tuple)
        if !any(isas, which)
            # easy
            extractarray(t, x -> $f(x, which))
        else
            DColumns(columns($f(t), which))
        end
    end
end

for f in [:rows, :keys, :values]
    @eval function $f(t::DNDSparse)
        extractarray(t, x -> $f(x))
    end

    @eval function $f(t::DNDSparse, which::Union{Int, Symbol})
        extractarray(t, x -> $f(x, which))
    end

    @eval function $f(t::DNDSparse, which::As)
        which.f($f(t, which.src))
    end
end

function column(t::DNDSparse, name)
    extractarray(t, x -> column(x, name))
end

columns(t::DNDSparse, which::Union{Int,Symbol,As}) = column(t, which)

function pairs(t::DNDSparse)
    extractarray(t, x -> map(Pair, x.index, x.data))
end

Base.@deprecate getindexcol(t::DNDSparse, dim) keys(t, dim)
Base.@deprecate getdatacol(t::DNDSparse, dim)  values(t, dim)
Base.@deprecate dindex(t::DNDSparse) keys(t)
Base.@deprecate ddata(t::DNDSparse)  values(t)


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ 5ccca702-c59f-11eb-34b6-011ca82d5945
md"""
# Neural Network
Only support classification tasks
"""

# ╔═╡ 79e17e30-a87c-450c-ac6f-bd7849c5c9c3
md"""
### Linear Layers (The Basics)
"""

# ╔═╡ e1d2cdc9-2628-40e3-9ca8-e63341dc1f2e
module tools include("../tools.jl") end

# ╔═╡ a4ab2ea6-8467-4995-b093-890f288404b2
JuTools = tools.JuTools

# ╔═╡ 14ff27dd-2926-4219-acdf-39fe6f1168e3
import Random

# ╔═╡ b0005ec6-5119-4786-b3d0-042677c61ee7
import Plots

# ╔═╡ bb1b6558-737a-4121-92a7-a78d1f1798a4
begin
X_train, Y_train = JuTools.data_generate_linear_2d(random_scale=10.0)
I = X_train
Y = Y_train
size(I), size(Y)
end

# ╔═╡ f19893ad-dc1d-40ee-8022-ee9535cc5a27
let
# plot the data
plot1 = X_train[Y_train .== 0, :]
plot2 = X_train[Y_train .== 1, :]
Plots.gr()
Plots.scatter(plot1[:, 1], plot1[:, 2], leg=false, c="red", background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot2[:, 1], plot2[:, 2], leg=false, c="blue")
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ b3aebb54-d447-400d-8291-473c5abb606a
begin
W₁ = Random.randn((size(X_train)[2], 10)) ./ 1000 # input layer (weights)
b₁ = zeros((1, 10))                               # input layer (biases)
W₂ = Random.randn((10, 5)) ./ 1000                # hidden layer (weights)
b₂ = zeros((1, 5))                                # hidden layer (biases)
W₃ = Random.randn((5, 1)) ./ 1000                 # output layer (weights)
b₃ = zeros((1, 1))                                # output layer (biases)
size(W₁), size(W₂), size(W₃)
end

# ╔═╡ ea42ea92-5794-462d-ba7e-0adad9848697
md"""
$$\text{MSE}=\frac{1}{2n}\sum^n_{i=1}(Y_i-\hat{Y}_i)^2$$
"""

# ╔═╡ da41e33b-ad12-4c23-b992-68bdbead07e4
# define a loss function (MSE)
function loss(Y::AbstractArray, Ŷ::AbstractArray)
    @assert size(Y) == size(Ŷ)
    N = size(Y)[1]
    L = (Y .- Ŷ) .^ 2
    return sum(L) / (N * 2)
end

# ╔═╡ 0cde5d94-41e6-4133-adb2-597b632ca863
md"""
$$C=\frac{1}{2n}\sum^n_{i=1}(\hat{Y} - Y)^2$$\

$$\hat{Y} = I_3\cdot W_3+b_3$$\

$$I_3 = I_2\cdot W_2+b_2$$\

$$I_2 = I\cdot W_1+b_1$$\

where $C$ is cost, $Y$ is true target, $W_i$ is the weights of layer $i$, $b_i$ is the biases of layer $i$, and $I$ is input data
"""

# ╔═╡ 50cc214c-3399-4cc0-a039-d7354f490bd6
begin
# forward pass
I₂ = I  * W₁ .+ b₁
I₃ = I₂ * W₂ .+ b₂
Ŷ  = I₃ * W₃ .+ b₃
Ŷ  = reshape(Ŷ, :)
loss(Y, Ŷ)
end

# ╔═╡ 0a10a7ff-22c2-4a70-a9ac-8b85026dbeda
md"""
Then we can compute gradient for each variable:\

$$\begin{aligned}
\frac{\partial C}{\partial W_3} &= \frac{\partial C}{\partial \hat{Y}}\cdot\frac{\partial \hat{Y}}{\partial W_3} \\
&= \frac{1}{n}\sum^n_{i=1} (\hat{Y}-Y) \cdot (I_3) \\
\end{aligned}$$

$$\begin{aligned}
\frac{\partial C}{\partial b_3} &= \frac{\partial C}{\partial \hat{Y}}\cdot\frac{\partial \hat{Y}}{\partial b_3} \\
&= \frac{1}{n}\sum^n_{i=1}(\hat{Y}-Y) \cdot (1) \\
\end{aligned}$$

$$\begin{aligned}
\frac{\partial C}{\partial W_2} &= \frac{\partial C}{\partial \hat{Y}}\cdot\frac{\partial \hat{Y}}{\partial I_3}\cdot\frac{\partial I_3}{\partial W_2} \\
&= \frac{1}{n}\sum^n_{i=1}(\hat{Y}-Y) \cdot (W_3) \cdot (I_2) \\
\end{aligned}$$

$$\begin{aligned}
\frac{\partial C}{\partial b_2} &= \frac{\partial C}{\partial \hat{Y}}\cdot\frac{\partial \hat{Y}}{\partial I_3}\cdot\frac{\partial I_3}{\partial b_2} \\
&= \frac{1}{n}\sum^n_{i=1}(\hat{Y}-Y) \cdot (W_3) \cdot (1) \\
\end{aligned}$$

$$\begin{aligned}
\frac{\partial C}{\partial W_1} &= \frac{\partial C}{\partial \hat{Y}}\cdot\frac{\partial \hat{Y}}{\partial I_3}\cdot\frac{\partial I_3}{\partial I_2}\cdot\frac{\partial I_2}{\partial W_1} \\
&= \frac{1}{n}\sum^n_{i=1}(\hat{Y}-Y) \cdot (W_3) \cdot (W_2) \cdot (I) \\
\end{aligned}$$

$$\begin{aligned}
\frac{\partial C}{\partial b_1} &= \frac{\partial C}{\partial \hat{Y}}\cdot\frac{\partial \hat{Y}}{\partial I_3}\cdot\frac{\partial I_3}{\partial I_2}\cdot\frac{\partial I_2}{\partial b_1} \\
&= \frac{1}{n}\sum^n_{i=1}(\hat{Y}-Y) \cdot (W_3) \cdot (W_2) \cdot (1) \\
\end{aligned}$$
"""

# ╔═╡ a0173ce1-8a16-4493-85f6-aaa138ef83c3
# define learning rate r
r = 0.1

# ╔═╡ 083f9702-bbfa-405e-bc18-f01f9c23d195
n = size(I)[1]

# ╔═╡ 4cb51df1-73d8-44b8-beba-9cbae56c3db2
gW₃ = I₃' * ((Ŷ .- Y) ./ n)

# ╔═╡ 91de2648-115e-4128-a98c-c69712ada92b
gb₃ = sum(Ŷ .- Y, dims=1) ./ n

# ╔═╡ e2b9b837-7b27-4543-be24-b1f261f2d824
gW₂ = I₂' * (((Ŷ .- Y) ./ n) * W₃')

# ╔═╡ 8fa62798-fc61-4b11-82f1-f5302bf574fb
gb₂ = (sum(Ŷ .- Y, dims=1) ./ n) * W₃'

# ╔═╡ fb07afc6-abfb-4434-a502-f3b67926ff3b
gW₁ = I' * ((((Ŷ .- Y) ./ n) * W₃') * W₂')

# ╔═╡ 0deb5d6b-991b-40dc-9cef-be3797d41ff5
gb₁ = ((sum(Ŷ .- Y, dims=1) ./ n) * W₃') * W₂'

# ╔═╡ bd0756cb-00ee-48c3-b663-d122e9e03553
begin
# backpropagate
W₁ .= W₁ .- r .* gW₁
b₁ .= b₁ .- r .* gb₁
W₂ .= W₂ .- r .* gW₂
b₂ .= b₂ .- r .* gb₂
W₃ .= W₃ .- r .* gW₃
b₃ .= b₃ .- r .* gb₃
end

# ╔═╡ 2af547e8-960e-4542-9cea-fcb7e97e06c2
let
# forward pass again
I₂ = I  * W₁ .+ b₁
I₃ = I₂ * W₂ .+ b₂
Ŷ  = I₃ * W₃ .+ b₃
Ŷ  = reshape(Ŷ, :)
loss(Y, Ŷ)
end

# ╔═╡ a8adef33-608f-4fdb-b482-0cac3a084220
# now ensemble everything into structs and functions
mutable struct ComputeGraph
    W_array::AbstractArray{AbstractArray} # weights
    b_array::AbstractArray{AbstractArray} # biases
    gW_array::AbstractArray{AbstractArray} # gradients of weights
    gb_array::AbstractArray{AbstractArray} # gradients of biases
end

# ╔═╡ 2b3e8e94-580a-4796-9218-2d16dff110a7
function createGraph(weightShapes::AbstractArray; randomize=false)::ComputeGraph
    # check validity of weight shapes
    for (shape1, shape2) in zip(weightShapes, weightShapes[2:end])
        @assert length(shape1) == length(shape2)
        @assert shape1[2] == shape2[1]
    end
    # initialize data
    W_array = []
    b_array = []
    gW_array = []
    gb_array = []
    for shape in weightShapes
        if randomize
            push!(W_array, Random.randn(shape) .* ((2.0 / sum(shape)) ^ 0.5))
        else
            push!(W_array, zeros(shape))
        end
        push!(b_array, zeros((1, shape[2])))
        push!(gW_array, zeros(shape))
        push!(gb_array, zeros((1, shape[2])))
    end
    # create graph
    graph = ComputeGraph(W_array, b_array, gW_array, gb_array)
    return graph
end

# ╔═╡ 75f3884d-540d-469d-a3b8-15245503224c
# this function compute output at layer defined by index recursively
function computeI(index::Integer, X::AbstractArray, graph::ComputeGraph)
    @assert index >= 1
    @assert index <= length(graph.W_array) + 1
    if index == 1
        return X
    end
    I = computeI(index-1, X, graph) * graph.W_array[index-1] .+ graph.b_array[index-1]
    return I
end

# ╔═╡ 63a5d479-ca9a-4f70-8c1b-5ce0d840fef1
# this function compute forward pass
# and return intermediate outputs at each layer
# this function will create huge arrays if input data is large
# so I will not use it in practice
function forward(X::AbstractArray, graph::ComputeGraph)::AbstractArray
    I_array = []
    for (W,b) in zip(graph.W_array, graph.b_array)
        if length(I_array <= 0)
            push!(I_array, (W * X .+ b))
        else
            push!(I_array, (W * I_array[end] .+ b))
        end
    end
    return I_array
end

# ╔═╡ a4e21ead-df90-4e58-9659-db9074199ff9
# this function compute gradients for each W,b and set values inplace
# it includes forward pass computation as well
function backpropagate!(X::AbstractArray, Y::AbstractArray, graph::ComputeGraph)
    tmp = (reshape(computeI(length(graph.W_array)+1, X, graph), :) .- Y) ./ length(Y)
    for index in length(graph.W_array):-1:1
        graph.gW_array[index] .= (computeI(index, X, graph)' * tmp)
        graph.gb_array[index] .= sum(tmp, dims=1)
        tmp = tmp * graph.W_array[index]'
    end
end

# ╔═╡ 312200d0-87f5-4591-91ca-035d6a1dcf04
# this function updates W,b based on their gradients
function step!(lr::AbstractFloat, graph::ComputeGraph)
    for index in 1:length(graph.W_array)
        graph.W_array[index] .= graph.W_array[index] .- lr .* graph.gW_array[index]
        graph.b_array[index] .= graph.b_array[index] .- lr .* graph.gb_array[index]
    end
end

# ╔═╡ 4cf70783-3486-4ab4-bf78-f862bc996d1e
# this function maps continuous prediction to categorical values (0, 1)
function predict(X::AbstractArray, graph::ComputeGraph)
    out = reshape(computeI(length(graph.W_array)+1, X, graph), :)
    out[out .>= 0.5] .= 1.0
    out[out .< 0.5]  .= 0.0
    return out
end

# ╔═╡ 8bfbf0e3-ffff-4e89-a5bd-aaca59ba7555
function predict_proba(X::AbstractArray, graph::ComputeGraph)
    out = reshape(computeI(length(graph.W_array)+1, X, graph), :)
    return out
end

# ╔═╡ a15c1a7d-15d8-44b9-8cff-b2c65afc9371
function train!(X::AbstractArray, Y::AbstractArray, graph::ComputeGraph, lr::AbstractFloat; max_iter::Integer=20, batch_size=50, verbose=false)
    @assert max_iter >= 0
    @assert 0 <= batch_size <= length(Y)
    for epoch in 1:max_iter
        X, Y = JuTools.shuffle_data(X, Y)
        for i in 1:batch_size:length(Y)
            if i + batch_size - 1 > length(Y)
                backpropagate!(X[i:end, :], Y[i:end, :], graph)
            else
                backpropagate!(X[i:(i+batch_size-1), :], Y[i:(i+batch_size-1), :], graph)
            end
            step!(lr, graph)
        end
        if verbose
            Ŷ = predict_proba(X, graph)
            L = loss(Y, Ŷ)
            println(string("Epoch: ", epoch, "\tLoss: ", L))
        end
    end
end

# ╔═╡ 046c8937-7acf-4e87-a8dd-421dfedc7ea2
begin
shapes = [(size(X_train)[2], 5), (5, 1)]
graph = createGraph(shapes, randomize=true)
train!(X_train, Y_train, graph, 0.0001, max_iter=10000, batch_size=50, verbose=false)
end

# ╔═╡ a5a626c0-ae4b-47fd-820a-aeef016add10
let
# plot decision probs
plot1 = X_train[Y_train .== 0, :]
plot2 = X_train[Y_train .== 1, :]
plot3 = collect(0:100)
plot_contour = hcat(vec(repeat(reshape(plot3, :, 1), 1, 101)'), vec(repeat(plot3, 101, 1)))
plot_contour = reshape(predict_proba(plot_contour, graph), 101, :)
Plots.gr()
Plots.contour(plot3, plot3, plot_contour, leg=true, fill=true, c=Plots.cgrad(:matter, rev=false, categorical=false),
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot1[:, 1], plot1[:, 2], leg=false, c="red")
Plots.scatter!(plot2[:, 1], plot2[:, 2], leg=false, c="blue")
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ eed53a70-abbe-4fbc-b105-da35cc237c23
let
# plot decision boundary
plot1 = X_train[Y_train .== 0, :]
plot2 = X_train[Y_train .== 1, :]
plot3 = collect(0:100)
plot_contour = hcat(vec(repeat(reshape(plot3, :, 1), 1, 101)'), vec(repeat(plot3, 101, 1)))
plot_contour = reshape(predict(plot_contour, graph), 101, :)
Plots.gr()
Plots.contour(plot3, plot3, plot_contour, leg=true, fill=true, c=Plots.cgrad(:matter, rev=false, categorical=false),
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot1[:, 1], plot1[:, 2], leg=false, c="red")
Plots.scatter!(plot2[:, 1], plot2[:, 2], leg=false, c="blue")
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ 35b5af9f-a410-47e5-bf77-09d9ab60b3be
JuTools.compute_accuracy(predict(X_train, graph), Y_train)

# ╔═╡ 49d0407a-224f-4200-83a5-c7564bcf1a8c
md"""
### About weights initialization

In current Neural Network setup, weights should __never__ be initialized to all zeros

For example:\
$$\begin{aligned}
\frac{\partial C}{\partial W_2} &= \frac{\partial C}{\partial \hat{Y}}\cdot\frac{\partial \hat{Y}}{\partial I_3}\cdot\frac{\partial I_3}{\partial W_2} \\
&= \frac{1}{n}\sum^n_{i=1}(\hat{Y}-Y) \cdot (W_3) \cdot (I_2)
\end{aligned}$$\

where gradient of $W_2$ is dependent on $W_3$\

However, since:\
$$\begin{aligned}
\frac{\partial C}{\partial W_3} &= \frac{\partial C}{\partial \hat{Y}}\cdot\frac{\partial \hat{Y}}{\partial W_3} \\
&= \frac{1}{n}\sum^n_{i=1} (\hat{Y}-Y) \cdot (I_3)
\end{aligned}$$\

where gradient of $W_3$ is dependent on input of layer 3, which is computed from $W_2$\

Both $W_2$ and $W_3$ will remain zeros during training\

Therefore, we need to initialize weights randomly\
But there are 2 problems:
* If weights are too large, gradients will go infinity, which is called __exploding gradients__
* If weights are too small, gradients will be near zero, which is called __vanishing gradients__

In both cases, the network fails to learn\

Assume we want a weight $w$ of shape $(m,n)$\
To solve this issue, we have 2 methods to initialize a weight:
* He Initialization
$$w=\text{randn}\cdot\sqrt{\frac{2}{n}}$$
* Xavier Initialization
$$w=\text{randn}\cdot\sqrt{\frac{1}{n}}$$

And I use the following method:

$$w=\text{randn}\cdot\sqrt{\frac{2}{m+n}}$$
"""

# ╔═╡ e3afc7de-f6ea-44e4-9be6-78ee952f62ab
md"""
### About learning rate

A large learning rate makes weights go infinity, while a small learning rate makes learning very slow

A better solution than fixed learning rate is to adjust learning rate based on gradient updates (Adaptive Learning Rate)\
Popular algorithms are:
* RMSProp
* Adagrad
* Adam

Other gradient descent optimization algorithms can be found [here](https://ruder.io/optimizing-gradient-descent/index.html#gradientdescentoptimizationalgorithms)
"""

# ╔═╡ Cell order:
# ╟─5ccca702-c59f-11eb-34b6-011ca82d5945
# ╟─79e17e30-a87c-450c-ac6f-bd7849c5c9c3
# ╠═e1d2cdc9-2628-40e3-9ca8-e63341dc1f2e
# ╠═a4ab2ea6-8467-4995-b093-890f288404b2
# ╠═14ff27dd-2926-4219-acdf-39fe6f1168e3
# ╠═b0005ec6-5119-4786-b3d0-042677c61ee7
# ╠═bb1b6558-737a-4121-92a7-a78d1f1798a4
# ╠═f19893ad-dc1d-40ee-8022-ee9535cc5a27
# ╠═b3aebb54-d447-400d-8291-473c5abb606a
# ╟─ea42ea92-5794-462d-ba7e-0adad9848697
# ╠═da41e33b-ad12-4c23-b992-68bdbead07e4
# ╟─0cde5d94-41e6-4133-adb2-597b632ca863
# ╠═50cc214c-3399-4cc0-a039-d7354f490bd6
# ╟─0a10a7ff-22c2-4a70-a9ac-8b85026dbeda
# ╠═a0173ce1-8a16-4493-85f6-aaa138ef83c3
# ╠═083f9702-bbfa-405e-bc18-f01f9c23d195
# ╠═4cb51df1-73d8-44b8-beba-9cbae56c3db2
# ╠═91de2648-115e-4128-a98c-c69712ada92b
# ╠═e2b9b837-7b27-4543-be24-b1f261f2d824
# ╠═8fa62798-fc61-4b11-82f1-f5302bf574fb
# ╠═fb07afc6-abfb-4434-a502-f3b67926ff3b
# ╠═0deb5d6b-991b-40dc-9cef-be3797d41ff5
# ╠═bd0756cb-00ee-48c3-b663-d122e9e03553
# ╠═2af547e8-960e-4542-9cea-fcb7e97e06c2
# ╠═a8adef33-608f-4fdb-b482-0cac3a084220
# ╠═2b3e8e94-580a-4796-9218-2d16dff110a7
# ╠═75f3884d-540d-469d-a3b8-15245503224c
# ╠═63a5d479-ca9a-4f70-8c1b-5ce0d840fef1
# ╠═a4e21ead-df90-4e58-9659-db9074199ff9
# ╠═312200d0-87f5-4591-91ca-035d6a1dcf04
# ╠═4cf70783-3486-4ab4-bf78-f862bc996d1e
# ╠═8bfbf0e3-ffff-4e89-a5bd-aaca59ba7555
# ╠═a15c1a7d-15d8-44b9-8cff-b2c65afc9371
# ╠═046c8937-7acf-4e87-a8dd-421dfedc7ea2
# ╠═a5a626c0-ae4b-47fd-820a-aeef016add10
# ╠═eed53a70-abbe-4fbc-b105-da35cc237c23
# ╠═35b5af9f-a410-47e5-bf77-09d9ab60b3be
# ╟─49d0407a-224f-4200-83a5-c7564bcf1a8c
# ╟─e3afc7de-f6ea-44e4-9be6-78ee952f62ab

max_stars_repo_name: teamclouday/JuliaLearn

max_issues_repo_name: teamclouday/JuliaLearn

max_forks_repo_name: teamclouday/JuliaLearn


text:
#=
An IP implementation of the Router for comparison against the Pathfinder
algorithm implemented in Mapper2.

I have no plans of making a full fledged router for full models like the one in
Mapper2. Instead, this model is going to be specialized for the network style
of KiloCore/KiloCore2 in that it will be an equivalent abstraction of the
circuit switched network, but not contain all of the details for the routing
muxes.

This is because a lot of these details are unnecessary to get a valid route, and
the fewer unnecessary variables we introduce to the IP forumlation, the better
it will perform (given the usual caveats that implies to IPs)

The basic idea is that the layered network of KC will be represented as two
connected mesh grids like the diagram below

                Circuit Layer 0

              o ----- o ----- o
             /       /       /|
            /       /       / |
           o ----- o ----- o  |
          /       /       /   |
         /       /       /    |
        o ----- o ----- o     o Sink Node
        |                     |
        | Source nodes will   |
        | connect both layers |
 Source |                     |
 Node   |                     |
        o     o ----- o ----- o
        |    /       /       /
        |   /       /       /
        |  o ----- o ----- o
        | /       /       /     Circuit Layer 1
        |/       /       /
        o ----- o ----- o

Both layers will not be connected to eachother except at:

1. Source nodes, to model processor's ability to write to any circuit net.
2. Destination nodes if the "preserve_dest" metadata is set to "false".
   Otherwise, only the destination node on the appropriate layer will be marked
   as a suitable destination.

-------------------
Problem Formulation
-------------------

The problem formulation is that of a constrained multi-commodity maximum network
flow. Communication channels from one processor to another are equivalend to
flows through a network, with the start of the flow coming from the vertex
representing the communicating processor and the end of the flow stopping at the
destination processor.

The primary objective is to minimize the total sum of links used to satisfy all
communication channels.

Constraints are:

1. Flow conservation constraints. The flow into a vertex must be equal to the
   flow out of a vertex. This ensures a continuous, non-branching path from
   source to sink and ensures we can't have spontaneous generation of flows.

   The equivalent to this constraint for terminal nodes is that flow out of a
   source must be 1, and flow into a destination must be 1.

2. Capacity constraints. We must ensure that a flow going through any node does
   not exceed 1 to model the physical constraints of the Kilocore communication
   architecture.

------------------------
Objective Interpretation
------------------------
The final objective is the number of global links used in the final circuit
switched network. We will just use this number and not bother trying to route
the actual Map data structure.

-------------
Optimizations
-------------
To reduce model size, mesh generation for source-sink pairs will be limited to
3 or 4 addresses around the bounding box generated by the source-sink pair.
=#

using JuMP
using Gurobi
using DataStructures
using LightGraphs
using ProgressMeter

import Mapper2.Routing
using Mapper2.Helper

# This is not necessarily scalable to an arbitrary number of layers, or
# interconnect topology.
@enum Network layer_0 layer_1 memory_request memory_return

struct AddressCircuit
    address :: Address{2}
    network :: Network
end

address(a::AddressCircuit) = a.address
getnetwork(a::AddressCircuit) = a.network

function iproute(m::Map{A,2}) where {A}
    # Constuct model
    ip_model, struct_time, struct_bytes, _, _ = @timed build_ip_model(m)
    # Solve the model
    status, solve_time, solve_bytes, _, _ = @timed solve(ip_model)

    m.metadata["ip_route_struct_time"] = struct_time
    m.metadata["ip_route_struct_bytes"] = struct_bytes
    m.metadata["ip_route_total_time"] = solve_time
    m.metadata["ip_route_solve_bytes"] = solve_bytes
    m.metadata["ip_route_objective"] = getobjectivevalue(ip_model)
    m.metadata["ip_route_solve_time"] = getsolvetime(ip_model)

    if status == :Optimal
        m.metadata["ip_optimal"] = true
    else
        m.metadata["ip_opcimal"] = false
    end

    m.metadata["ip_status"] = status

    return m
end

function build_ip_model(m::Map{2})
    # Build the network data type.
    network = build_network(m)

    # Instantiate an empty IP model
    ip_model = Model(solver = GurobiSolver(Threads=1))

    # To create variables for JuMP, we need to do it all in one shot. So, we
    # must iterate over the network to get all the variables
    # (nodes in the network) that we need.
    vertices_for_taskgraph_edge = Vector{Int}[]

    # Amount to grow the bounding box around the extrema of the source-sink
    # pair. This helps keep the number of generated variables down without
    # hurting routing. Keep somewhere between 2 and 4.
    box_growth = 3

    # Collection of vertices used for each edge in the taskgraph.
    vertex_collection = Vector{Int}[]
    # Indices of start and stop nodes for each edge of the graph.
    source_collection = Vector{Int}[]
    sink_collection   = Vector{Int}[]

    num_edges = Mapper2.num_edges(m.taskgraph)

    for (edge_index, taskgraph_edge) in enumerate(getedges(m.taskgraph))
        # Get the source and destination tasks, as well as the addresses
        # for these tasks.
        #
        # Since there are no taskgraphs with fanout, just get the first
        # of the sources/sinks as there should only be 1.
        source_name = first(getsources(taskgraph_edge))
        sink_name = first(getsinks(taskgraph_edge))

        source_address = getaddress(m, source_name)
        sink_address = getaddress(m, sink_name)

        source_task = getnode(m.taskgraph, source_name)
        sink_task = getnode(m.taskgraph, sink_name)

        # Create a bounding box for this address
        box = BoundingBox(source_address, sink_address, box_growth)

        # Get the network for this edge.

        # If the sink is a memory, this uses the "memory_request" link
        if ismemory(sink_task)
            links = [memory_request]

        # Otherwise if the source is a memory, it is a "memory_return" link.
        # Note that there should be no links where both the sources and the
        # sinks are memories.
        elseif ismemory(source_task)
            links = [memory_return]

        # Default - use the circuit links.
        else
            links = [layer_0, layer_1]
        end

        # Get the source and sink vertices for this layer.
        local_source = getvertices(network, source_address, links)
        local_sink = getvertices(network, sink_address, links)
        push!(source_collection, local_source)
        push!(sink_collection, local_sink)

        # Now, collect all vertices that are within the bounding box
        local_vertices = Int[]
        for addr in addresses(box)
            append!(local_vertices, getvertices(network, addr, links))
        end
        push!(vertex_collection, local_vertices)
    end

    # We are now ready to generate the IP model.
    @variable(ip_model,
        x[i = 1:num_edges,
          j = vertex_collection[i],
          k = outneighbors(network, j);
          in(k, vertex_collection[i])
        ], Bin)

    @objective(ip_model, Min, sum(x))

    # Add flow constraints. For start and stop edges, make sure the
    # outgoing/incoming sum is 1.
    #
    # Must also ensure that flow out of a sink node and into a source node are
    # zero to avoid out->in loops.
    #
    # For all other edges, incoming flow must equal outgoing flow.

    for i in 1:num_edges
        valid_vertices = vertex_collection[i]
        @constraint(ip_model, sum(x[i, j, k] for j in source_collection[i], k in outneighbors(network, j) if k in valid_vertices) == 1)
        @constraint(ip_model, sum(x[i, k, j] for j in source_collection[i], k in inneighbors(network, j) if k in valid_vertices) == 0)

        @constraint(ip_model, sum(x[i, j, k] for j in sink_collection[i], k in outneighbors(network, j) if k in valid_vertices) == 0)
        @constraint(ip_model, sum(x[i, k, j] for j in sink_collection[i], k in inneighbors(network, j) if k in valid_vertices) == 1)

        for vertex in valid_vertices
            # Skip source and sink vertices
            (vertex in source_collection[i] || vertex in sink_collection[i]) && continue

            # Flow in and flow out must be equal.
            @constraint(ip_model,
                sum(x[i,vertex,k] for k in outneighbors(network, vertex) if k in valid_vertices) -
                sum(x[i,k,vertex] for k in  inneighbors(network, vertex) if k in valid_vertices) == 0)
        end
    end

    # Add capacity constraints.
    for edge in edges(network)
        source = src(edge)
        dest   = dst(edge)

        # Collect all channels that have this edge.
        indices = Int[]
        for i in 1:num_edges
            if source in vertex_collection[i] && dest in vertex_collection[i]
                push!(indices, i)
            end
        end

        # Add capacity constraint.
        @constraint(ip_model, sum(x[i, source, dest] for i in indices) <= 1)
    end

    return ip_model
end

function build_network(m::Map{2})
    # First, build a light graph for the network given above.
    network = SimpleNetwork()

    # Iterate through all addresses. Add circuit and memory nodes as necessary.
    for addr in addresses(m.architecture)
        # Get the component at this address, and then get the mappable component
        # child of this component.
        component = m.architecture[addr]
        mappables = collect(Iterators.filter(ismappable, component[p] for p in walk_children(component)))
        @assert length(mappables) == 1

        mappable_component = first(mappables)

        # Create vertices for this component
        if isinput(mappable_component) || isoutput(mappable_component)
            # Add just one vertex because the input/output handlers only connect
            # to circuit net 0
            #
            # Since we just added a vertex, the index of that vertex is simply
            # the number of vertices.
            add_vertex!(network, addr, layer_0)
        elseif ismemoryproc(mappable_component)
            # Add both layers for memory processor as well as memory request and
            # and memory return nodes.
            add_vertex!(network, addr, layer_0)
            add_vertex!(network, addr, layer_1)
            add_vertex!(network, addr, memory_return)
            add_vertex!(network, addr, memory_request)
        elseif isproc(mappable_component)
            # Just add the two circuit layers
            add_vertex!(network, addr, layer_0)
            add_vertex!(network, addr, layer_1)
        elseif ismemory(mappable_component)
            # Add the two memory networks
            add_vertex!(network, addr, memory_return)
            add_vertex!(network, addr, memory_request)
        else
            error("Unrecognized component type")
        end
    end

    # Now, extract adjacency information to start connecting nodes together
    connected_components = MapperCore.connected_components(m.architecture)
    for (source,destinations) in connected_components
        for dest in destinations
            add_edge!(network, source, dest)
        end
    end

    return network
end

################################################################################
# Abstraction of the network for convenience purposes.

struct SimpleNetwork
    graph               :: LightGraphs.SimpleGraphs.SimpleDiGraph{Int64}
    vertex_data         :: Dict{Int, AddressCircuit}
    address_to_vertex   :: Dict{Address{2}, Vector{Int}}
end

SimpleNetwork() = SimpleNetwork(
    DiGraph(),
    Dict{Int, AddressCircuit}(),
    Dict{Address{2},Vector{Int}}()
)

getdata(n::SimpleNetwork, i::Int) = n.vertex_data[i]
getvertices(n::SimpleNetwork, a::Address{2}) = get(n.address_to_vertex, a, Int[])

function getvertices(n::SimpleNetwork, a::Address{2}, networks)
    return collect(Iterators.filter(
        x -> getdata(n, x).network in networks,
        i for i in getvertices(n, a)
    ))
end

function LightGraphs.add_vertex!(n::SimpleNetwork, address::Address{2}, network)
    add_vertex!(n.graph)
    n.vertex_data[nv(n.graph)] = AddressCircuit(address, network)

    # Register a vertex at this address.
    if haskey(n.address_to_vertex, address)
        push!(n.address_to_vertex[address], nv(n.graph))
    else
        n.address_to_vertex[address] = [nv(n.graph)]
    end
    return nothing
end

function LightGraphs.add_edge!(n::SimpleNetwork, source::Address, dest::Address)
    # Get the lists of nodes at the source and destination vertices.
    source_vertices = n.address_to_vertex[source]
    dest_vertices   = n.address_to_vertex[dest]

    # Add edges from source to destination by matching circuit types.
    # Note that for things like memory links, this will generate a few extra
    # variables that aren't strictly needed, but this is minimal and won't
    # affect the final routing.
    for sv in source_vertices
        source_data = n.vertex_data[sv]
        for dv in dest_vertices
            dest_data = n.vertex_data[dv]
            if getnetwork(source_data) == getnetwork(dest_data)
                add_edge!(n.graph, sv, dv)
            end
        end
    end
end

# Simple forwards.
LightGraphs.nv(n::SimpleNetwork) = nv(n.graph)
LightGraphs.ne(n::SimpleNetwork) = ne(n.graph)
LightGraphs.outneighbors(n::SimpleNetwork, i) = outneighbors(n.graph, i)
LightGraphs.inneighbors(n::SimpleNetwork, i) = inneighbors(n.graph, i)
LightGraphs.edges(n::SimpleNetwork) = edges(n.graph)

################################################################################
struct BoundingBox{D}
    bounds::NTuple{D,UnitRange{Int}}
end

function BoundingBox(a::Address{D}, b::Address{D}, growth = 0) where D
    return BoundingBox(Tuple(map(1:D) do i
        minval, maxval = minmax(a.I[i], b.I[i])
        return (minval - growth):(maxval + growth)
    end))
end

Mapper2.addresses(b::BoundingBox) = (CartesianIndex(i) for i in Iterators.product(b.bounds...))

function BoundingBox(min::NTuple{D}, max::NTuple{D}) where D
    return BoundingBox(Tuple(min[i]:max[i] for i in 1:D))
end

function Base.in(addr::Address{D}, box::BoundingBox{D}) where D
    for i in 1:D
        in(addr.I[i], box.bounds[i]) || return false
    end
    return true
end