text:
module MoodleQuiz

using LaTeXStrings
using FileIO
using LightXML
using HttpCommon
using Markdown
using OrderedCollections
using Base: @kwdef
using Pipe: @pipe
using Base64

import Base.show

export Question, Category, Quiz, moodle_xml, @cz_str, @en_str

export TrueFalse, True, False, YesNo, Yes, No, MultiChoice, SingleChoice, Choice
export Correct, Wrong, Essay, ShortAnswer, Numerical

## Helper functions
ensurearray(x) = if isa(x, Array); x else [ x ] end

function escape_chars(io::IO, s::AbstractString, chars)
    a = Iterators.Stateful(s)
    for c::AbstractChar in a
        if c in chars
            print(io, '\\', c)
        else
            print(io, c)
        end
    end
end

## Data types & contructors
@kwdef struct Text
    text::Markdown.MD
    lang::String
end
macro cz_str(text) return :( Text(text=Markdown.parse($text), lang="cs") ) end
macro en_str(text) return :( Text(text=Markdown.parse($text), lang="en") ) end

# Answer data types

abstract type Answer end

struct Essay <: Answer end

struct TrueFalse <: Answer
    correct::Bool
end

True() = TrueFalse(true)
False() = TrueFalse(false)

struct YesNo <: Answer
    yes::Bool
end

Yes() = YesNo(true)
No() = YesNo(false)

struct Choice
    correct::Bool
    content::Vector{Any}
end

Choice(content; correct::Bool) = Choice(correct, ensurearray(content))

Correct(content) = Choice(content, correct=true)
Wrong(content) = Choice(content, correct=false)

@enum AnswerNumbering none=0 abc=1 ABCD=2 numeric=3

@kwdef struct MultiChoice <: Answer
    single::Bool = false
    shuffle::Bool = true
    correct_feedback = ""
    partiallycorrectfeedback = ""
    incorrectfeedback = ""
    answernumbering::AnswerNumbering = none

    choices::Vector{Choice}
    # TODO: Enforce single correct choice if single==true
end

MultiChoice(choices...; kwargs...) = MultiChoice(; choices=collect(choices), kwargs...)
SingleChoice(choices...; kwargs...) = MultiChoice(; choices=collect(choices), single=true, kwargs...)

struct Numerical <: Answer
    val::Real
    tolerance::Real
end
Numerical(val) = Numerical(val, 0.0)

struct ShortAnswer <: Answer
    text::String
    fraction::Int               # 0 - 100 (or points for cloze)
end
ShortAnswer(text) = ShortAnswer(text, 100)

# Question, Category, Quiz

@kwdef struct Question
    id::String
    points::Float64
    dates::String
    content::Vector{Any}
    answer::Union{Answer, Nothing} = nothing

    responseformat::String = "editorfilepicker"
    responserequired::Bool = true
    responsefieldlines::Int64 = 5
    attachments::Int64 = 0
    attachmentsrequired::Int64 = 0
end

struct Category
    name::String
    questions::OrderedDict{String,Question}
end
Category(name) = Category(name, OrderedDict())
default_category = Category("")

Question(id, points, dates, content, answer=nothing; kwargs...) =
    Question(default_category, id, points, dates, content, answer, kwargs...)

function Question(category::Category, id, points, dates, content, answer=nothing; kwargs...)
    if isnothing(answer) && !any(typeof.(content) .<: Answer)
        answer = Essay()
    end
    q = Question(; id, points, dates, content, answer, kwargs...)
    category.questions[q.id] = q
    return q
end

@kwdef struct Quiz
    questions::Vector{Question} = [] # questions without category
    categories::Vector{Category} = [] # questions with category
end

Quiz(q::Question) = Quiz(questions=[q])
Quiz(c::Category) = Quiz(categories=[c])

## Preview methods

preview(io::IO, ::Nothing) = nothing
preview(io::IO, ::Essay) = print(io, "Essay")
preview(io::IO, ::ShortAnswer) = print(io, "ShortAnswer")
preview(io::IO, n::Numerical) = print(io, "$(n.val) ± $(n.tolerance)")
preview(io::IO, tf::TrueFalse) = print(io, "($(tf.correct ? "x" : " ")) True  ($(tf.correct ? " " : "x")) False")
preview(io::IO, yn::YesNo) = print(io, "($(yn.yes ? "x" : " ")) Yes  ($(yn.yes ? " " : "x")) No")
function preview(io::IO, choice::Choice)
    for val in choice.content
        print(io, repr("text/html", val))
    end
end

function preview(io::IO, mc::MultiChoice)
    print(io, "
")
    for choice in mc.choices
        x = choice.correct ? "x" : " "
        print(io, "
" * (mc.single ? "($x)" : "[$x]") * " ")
        preview(io, choice)
        print(io, "
")
    end
    println(io, "
")
end

function preview(io::IO, q::Question)
    println(io, "
Question $(q.id)
")
    println(io, "
Points: $(q.points)
")
    println(io, to_html(q))
    println(io, "
Answer: ")
    preview(io, q.answer)
    println(io, "
")
end

function preview(io::IO, cat::Category)
    println(io, "
Category $(cat.name)
")
    for q in values(cat.questions)
        preview(io, q)
    end
end

function preview(io::IO, quiz::Quiz)
    println(io, "
Quiz
")
    for q in quiz.questions
        preview(io, q)
    end
    for cat in quiz.categories
        preview(io, cat)
    end
end

function preview(object)
    open("questions.html", "w") do io
        # TODO: Remove next line - temporary addition
        println(io, """""")
        preview(io, object)
    end
    run(`firefox-reload`)
end

function Base.show(io::IO, ::MIME"text/html", t::Text)
    print(io, """""" * repr("text/html", t.text) * "")
end

png_cache = Dict{String, Any}()

function render_png(s::LaTeXString; dpi=150, debug=false, name=tempname(cleanup=false))
    doc = """
    \\documentclass[varwidth=100cm]{standalone}
    \\usepackage{tikz}
    \\usetikzlibrary{arrows.meta}
    \\usepackage{amssymb}
    \\usepackage{amsmath}
    \\begin{document}
    {
        $s
    }
    \\end{document}
    """
    doc = replace(doc, "\\begin{align}"=>"\\[\n\\begin{aligned}")
    doc = replace(doc, "\\end{align}"=>"\\end{aligned}\n\\]")

    global png_cache
    if doc in keys(png_cache)
        return png_cache[doc]
    end
    try
        open("$(name).tex", "w") do f
            write(f, doc)
        end
        cd(dirname(name)) do
            cmd = `lualatex --interaction=$(debug ? "nonstopmode" : "batchmode") $(name).tex`
            debug || (cmd = pipeline(cmd, devnull))
            run(cmd)
            pdftoppm = run(`pdftoppm -r $(dpi) -png $(name).pdf $(name)`)
        end
        png = load("$(name)-1.png")
        png_cache[doc] = png
        return png
    finally
        Base.Filesystem.rm("$(name).tex")
        Base.Filesystem.rm("$(name).pdf")
        Base.Filesystem.rm("$(name)-1.png")
    end
end

to_html(x) = repr("text/html", x)

# Cloze sub-questions
function cloze(type::String,
               answers::String;
               weight::Int = 1) # Why UInt doesn't work?
    " {$(weight):$(type):$(answers)} "
end
function cloze(type::String, answers::Vector{Choice}; kwargs...)
    ansstr = join(map(answers) do a
                  (a.correct ? "=" : "") *
                  # https://docs.moodle.org/310/en/Embedded_Answers_(Cloze)_question_type#Detailed_syntax_explanations
                  sprint(escape_chars, join(map(to_html, a.content)), "#~/\"\\")
                  end, "~")
    cloze(type, ansstr; kwargs...)
end
cloze_answer(answer, feedback) = "$answer#$feedback"
to_html(yn::YesNo) = cloze("MULTICHOICE", [Choice(yn.yes, ["Yes"]), Choice(!yn.yes, ["No"])])
to_html(sa::ShortAnswer) = cloze("SA", [Correct(sa.text)])
to_html(num::Numerical) = cloze("NUMERICAL", "=$(num.val):$(num.tolerance)")
to_html(mc::MultiChoice) = cloze("MULTICHOICE", "=TODO")

function to_html(x::LaTeXString)
    png = render_png(x)
    io = IOBuffer()
    iob64_encode = Base64EncodePipe(io);
    save(Stream(format"PNG", iob64_encode), png)
    """"""
end

to_html(a::Vector{Any}) = join(to_html.(a))

function Base.show(io::IO, ::MIME"text/html", q::Question)
    for el in q.content
        print(io, to_html(el))
    end
end

function Base.show(io::IO, ::MIME"text/plain", q::Question)
    println(io, "Question($(q.id))")

    # FIXME: Not very generic way of quickly previewing the question
    preview(q)
end

moodle_xml_type(::TrueFalse) = "truefalse"
moodle_xml_type(::ShortAnswer) = "shortanswer"
moodle_xml_type(::Numerical) = "numerical"
moodle_xml_type(::MultiChoice) = "multichoice"
moodle_xml_type(::YesNo) = "multichoice"
moodle_xml_type(::Essay) = "essay"
moodle_xml_type(x) = nothing

function moodle_xml_type(q::Question)
    type = moodle_xml_type(q.answer)
    if isnothing(type) && any(typeof.(q.content) .<: Answer)
        type = "cloze"
    end
   return type
end

function moodle_xml_answer(xquestion, text::String, fraction::Int; kwargs...)
    xans = new_child(xquestion, "answer")
    set_attribute(xans, "fraction", string(fraction))
    xtext = new_child(xans, "text")
    add_text(xtext, text)
    for (tag, val) in kwargs
        xtag = new_child(xans, string(tag))
        add_text(xtag, string(val))
    end
end

moodle_xml(tf::TrueFalse, xquestion) = begin
    moodle_xml_answer(xquestion, "true", tf.correct ? 100 : 0)
    moodle_xml_answer(xquestion, "false", !tf.correct ? 100 : 0)
end
moodle_xml(yn::YesNo, xquestion) = begin
    add(elem::String, val) = add_text(new_child(xquestion, elem), string(val))
    add("single", true)
    moodle_xml_answer(xquestion, "true", yn.yes ? 100 : 0)
    moodle_xml_answer(xquestion, "false", !yn.yes ? 100 : 0)
end
moodle_xml(na::Numerical, xquestion) = moodle_xml_answer(xquestion, string(na.val), 100; na.tolerance)
moodle_xml(sa::ShortAnswer, xquestion) = moodle_xml_answer(xquestion, sa.text, sa.fraction)
# TODO: Multiple short and numerical answers
moodle_xml(mc::MultiChoice, xquestion) = begin
    for c in mc.choices
        moodle_xml_answer(xquestion, to_html(c.content), c.correct ? 100 : 0)
    end

    add(elem::String, val) = add_text(new_child(xquestion, elem), string(val))

    add("single", mc.single)
    add("shuffleanswers", mc.shuffle)
    add("answernumbering", "none")
end
moodle_xml(a::Answer, xquestion) = nothing

function moodle_xml(q::Question, xroot)
    xquestion = new_child(xroot, "question")
    # Set the type
    set_attribute(xquestion, "type", moodle_xml_type(q))

    name = new_child(new_child(xquestion, "name"), "text")
    add_text(name, q.id)

    # Put the question text
    questiontext = new_child(xquestion, "questiontext")
    set_attribute(questiontext, "format", "html")
    text = new_child(questiontext, "text")
    add_text(text, repr("text/html", q))

    if q.answer != nothing
        moodle_xml(q.answer, xquestion)
    end

    add(elem::String, val) = add_text(new_child(xquestion, elem), string(val))
    add(elem::String, val::Bool) = add(elem, Int(val))

    add("defaultgrade", q.points)
    add("responseformat", q.responseformat)
    add("responserequired", q.responserequired)
    add("responsefieldlines", q.responsefieldlines)
    add("attachments", q.attachments)
    add("attachmentsrequired", q.attachmentsrequired)
end

function moodle_xml(cat::Category, xroot)
    xquestion = new_child(xroot, "question")
    # Set the type
    set_attribute(xquestion, "type", "category")
    @pipe xquestion |> new_child(_, "category") |> new_child(_, "text") |> add_text(_, "\$course\$/top/MoodleQuiz.jl/$(cat.name)")
    for q in values(cat.questions)
        moodle_xml(q, xroot)
    end
end

moodle_xml(q::Question) = moodle_xml(Quiz(q))
moodle_xml(c::Category) = moodle_xml(Quiz(c))

function moodle_xml(quiz::Quiz)::XMLDocument
    xdoc = XMLDocument()
    # Create test
    xroot = create_root(xdoc, "quiz")
    for q in quiz.questions
        moodle_xml(q, xroot)
    end
    if default_category in quiz.categories
        for q in values(default_category.questions)
            moodle_xml(q, xroot)
        end
    end
    for cat in quiz.categories
        cat === default_category && continue
        moodle_xml(cat, xroot)
    end
    return xdoc
end

end


text:
using Dash, DashHtmlComponents, DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_input(id = "input-1", type = "text", value = "Montreal"),
    dcc_input(id = "input-2", type = "text", value = "Canada"),
    html_div(id = "output-keywords")
end

callback!(
    app,
    Output("output-keywords", "children"),
    Input("input-1", "value"),
    Input("input-2", "value"),
) do input_1, input_2
    return "Input 1 is \"$input_1\" and Input 2 is \"$input_2\""
end

run_server(app, "0.0.0.0", debug=true)

ext: jl

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -3.428579430434188 <= q <= 8.991868736935128)

                     add_NL_constraint(m, :(gelu(0.758903986598098 + 0.8106379052679875*gelu(0.47487042435595583 + 0.4488775132618512*gelu(0.1716391672055586 + 0.7301849089692007*gelu(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + 0.6276077641275091*gelu(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))) + 0.015398554931854491*gelu(0.8366427516313624 + 0.8654420505988725*gelu(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + -0.07448940361379819*gelu(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3])))) + -0.6668215058198093*gelu(0.8216722410540171 + 0.5667551555087149*gelu(0.1716391672055586 + 0.7301849089692007*gelu(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + 0.6276077641275091*gelu(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))) + 0.5912184082121308*gelu(0.8366427516313624 + 0.8654420505988725*gelu(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + -0.07448940361379819*gelu(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))))) + gelu(-0.09143860592972564 + 0.3407854070498746*gelu(0.47487042435595583 + 0.4488775132618512*gelu(0.1716391672055586 + 0.7301849089692007*gelu(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + 0.6276077641275091*gelu(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))) + 0.015398554931854491*gelu(0.8366427516313624 + 0.8654420505988725*gelu(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + -0.07448940361379819*gelu(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3])))) + 0.24612879529032217*gelu(0.8216722410540171 + 0.5667551555087149*gelu(0.1716391672055586 + 0.7301849089692007*gelu(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + 0.6276077641275091*gelu(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))) + 0.5912184082121308*gelu(0.8366427516313624 + 0.8654420505988725*gelu(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + -0.07448940361379819*gelu(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Base.LinAlg.CHOLMOD

for Ti in (:Int32,:Int64)
    @eval begin
        function chm_analyze_ordering{Tv<:CHOLMOD.CHMVTypes}(A::CHOLMOD.CholmodSparse{Tv,$Ti},
                                                             ordering::Cint,
                                                             Perm::Vector{$Ti})
            Parent   = Array($Ti,A.c.m)
            Post     = Array($Ti,A.c.m)
            ColCount = Array($Ti,A.c.m)
            First    = Array($Ti,A.c.m)
            Level    = Array($Ti,A.c.m)

            ccall((@CHOLMOD.chm_nm("analyze_ordering",$Ti),:libcholmod), Cint,
                  (Ptr{CHOLMOD.c_CholmodSparse{Tv,$Ti}}, Cint, Ptr{$Ti}, Ptr{$Ti}, Csize_t,
                   Ptr{$Ti}, Ptr{$Ti}, Ptr{$Ti},
                   Ptr{$Ti}, Ptr{$Ti}, Ptr{Uint8}),
                   &A.c, ordering, Perm, C_NULL, zero(Csize_t), Parent, Post, ColCount, First, Level, CHOLMOD.cmn($Ti))
            Parent, Post, ColCount
        end

        function chm_analyze_p{Tv<:CHOLMOD.CHMVTypes}(A::CholmodSparse{Tv,$Ti}, Perm::Vector{$Ti})
            F = ccall((@CHOLMOD.chm_nm("analyze_p",$Ti),:libcholmod), Ptr{CHOLMOD.c_CholmodFactor{Tv,$Ti}},
                      (Ptr{CHOLMOD.c_CholmodSparse{Tv,$Ti}}, Ptr{$Ti}, Ptr{$Ti}, Csize_t, Ptr{Uint8}),
                      &A.c, Perm, C_NULL, zero(Csize_t), CHOLMOD.cmn($Ti))
            CholmodFactor(F)
        end

        function chm_factorize_p!{Tv<:CHOLMOD.CHMVTypes}(A::CholmodSparse{Tv,$Ti}, beta::Tv, L::CholmodFactor{Tv,$Ti})
             ccall((@CHOLMOD.chm_nm("factorize_p",$Ti),:libcholmod), Cint,
                  (Ptr{CHOLMOD.c_CholmodSparse{Tv,$Ti}}, Ptr{Tv}, Ptr{$Ti}, Csize_t,
                   Ptr{CHOLMOD.c_CholmodFactor{Tv,$Ti}}, Ptr{Uint8}),
                   &A.c, &beta, C_NULL, zero(Csize_t), &L.c, CHOLMOD.cmn($Ti))
        end

    end
end

function pothen_sun{T<:Integer}(par::Vector{T}, post::Vector{T}, colcount::Vector{T})
    n = length(par)
    flag = -ones(T,n)
    ns = n
    for j=post
        mdeg = colcount[j] - 1
        if par[j] != 0 && mdeg == colcount[par[j]] && flag[par[j]] == -1
            # par[j] not assigned to supernode
            ns -= 1
            if flag[j] < 0    # j is a repr. vertex
                flag[par[j]] = j
                flag[j] -= 1
            else              # j is not a repr. vertex
                flag[par[j]] = flag[j]
                flag[flag[j]] -= 1
            end
        end
    end
    ns, flag
end

function chordal_embedding{Tv<:Number,Ti<:Int}(A::SparseMatrixCSC{Tv,Ti}, Perm::Vector{Ti})
    m, n = size(A)
    S = CholmodSparse(SparseMatrixCSC(m,n,A.colptr,A.rowval,zeros(length(A.nzval))))
    F = chm_analyze_p(S, length(Perm)>0 ? Perm-1 : Perm)
    chm_factorize_p!(S, 1.0, F)

    par, post, colcount = chm_analyze_ordering(S, F.c.ordering, F.Perm)
    ns, flag = pothen_sun(par+1, post+1, colcount)

    # extract cliques
    L = sparse(F)
    [ sort(1+F.Perm[L.rowval[L.colptr[i]:L.colptr[i+1]-1]]) for i = find(flag .< 0) ]
end

# If no permutation is specified, use the CHOLMOD ordering
chordal_embedding{Tv<:Number,Ti<:Int}(A::SparseMatrixCSC{Tv,Ti}) = chordal_embedding(A, Array(Int,0))

function clique_index(I::Array{Array{Int,1},1}, S::Array{Int,1})
    s = IntSet(S...)
    for k=1:length(I)
        if s <= IntSet(I[k]...) return k end
    end
    throw(ArgumentError)
end


text:
export set_options!, set_parameters!, get_option, get_parameter,
       set_limits!, set_setpoint!, set_user_parameters!

"""
    set_options!( model::GRAMPC.Model; options... )

Set solver options in the specified GRAMPC model. The options are specified as keyword arguments,
where the keyword is the option name and the value is the desired value.
"""
set_options!( model::GRAMPC.Model; options... ) = _set_settings!( model; options... )

"""
    set_parameters!( model::GRAMPC.Model; params... )

Set solver parameters in the specified GRAMPC model. The parameters are specified as keyword arguments,
where the keyword is the parameter name and the value is the desired value.
"""
set_parameters!( model::GRAMPC.Model; params... ) = _set_settings!( model; params... )

"""
    get_option( model::GRAMPC.Model, option::Symbol )

Get the current value of the solver option `option` in the GRAMPC model `model`.
"""
get_option( model::GRAMPC.Model, option::Symbol ) = _get_setting( model, option )

"""
    get_parameter( model::GRAMPC.Model, param::Symbol )

Get the current value of the parameter `param` in GRAMPC model `model`.
"""
get_parameter( model::GRAMPC.Model, param::Symbol ) = _get_setting( model, param )

"""
    set_user_parameters!( model::GRAMPC.Model{T}, userparam::T ) where T

Set the user-defined parameters that are passed into every function callback.
"""
set_user_parameters!( model::GRAMPC.Model{T}, userparam::T ) where T = ( model.param= userparam )

"""
    set_limits!( model::GRAMPC.Model, type::Symbol; lb=nothing, ub=nothing )

Set the upper and lower bounds for the inputs and parameters. The valid values for `type` are
`:input` or `:u` for input bounds and `:parameters` or `:p` for parameter bounds. If a bound is
set to `nothing`, it will default to ``∞/-∞`` for the upper/lower bounds respectively.
"""
function set_limits!( model::GRAMPC.Model, type::Symbol; lb=nothing, ub=nothing )
    settingsKey = ""
    len = 0

    if type == :input || type == :u
        settingsKey = "u"
        len = model.Nu
    elseif type == :parameters || type == :p
        settingsKey = "p"
        len = model.Np
    else
        error( "Unknown type $(type)")
    end

    # If a bound isn't specified, it is infinite
    isa( lb, Nothing ) || ( lb = fill( -Inf, len ) )
    isa( ub, Nothing ) || ( ub = fill(  Inf, len ) )

    d = Dict( Symbol( "$(settingsKey)min" ) => lb,
              Symbol( "$(settingsKey)max" ) => ub )

    set_parameters!( model; d... )
end

"""
    set_setpoint!( model::GRAMPC.Model; input=nothing, state=nothing )

Set the desired setpoint for the states and inputs. If it is set to `nothing`, then
no change is made to the setpoint in the problem.
"""
function set_setpoint!( model::GRAMPC.Model; input=nothing, state=nothing )
    d = Dict()

    # If a setpoint isn't specified, it isn't set
    isa( state, Nothing ) || ( d[:xdes] = state )
    isa( input, Nothing ) || ( d[:udes] = input )

    set_parameters!( model; d... )
end


text:
using AlgebraResultTypes: result_field, result_ring, number_type
using Test, Random
import ForwardDiff

####
#### add to this list to test for more types
####

TEST_TYPES = (# Base
              Float64, Int, Rational{Int}, Complex{Float64}, Rational{Int16}, Int8,
              Float32,
              # packages -- ADD TESTS HERE
              ForwardDiff.Dual{:foo,Float64,3})

####
#### deterministic tests
####

function test_field(a)
    abs2(a + a + one(a))/(-a)
end

function test_field(a, b)
    abs2(a + a + one(a))/(-a) + (a/b - b/a + a/a + b/b)
end

function test_ring(a)
    abs2(a + a + one(a))*(-a)
end

function test_ring(a, b)
    abs2(a + a + one(a))*(-a) + (a*b - b*a + a - b)
end

@testset "result_field deterministic tests" begin
    @test result_field(Real) ≡ Number
    for T in TEST_TYPES
        @test typeof(test_field(one(T))) ≡ result_field(T)
        @test result_field(T, Real) ≡ Number
        for S in TEST_TYPES
            @test typeof(test_field(one(T), one(S))) ≡ result_field(T, S)
            @test result_field(T, S, Real) ≡ Number
            for Z in TEST_TYPES
                @test typeof(test_field(one(T), test_field(one(S), one(Z)))) ≡ result_field(T, S, Z)
            end
        end
    end
end

@testset "result_ring deterministic tests" begin
    @test result_ring(Real) ≡ Number
    for T in TEST_TYPES
        @test typeof(test_ring(one(T))) ≡ result_ring(T)
        @test result_ring(T, Real) ≡ Number
        for S in TEST_TYPES
            @test typeof(test_ring(one(T), one(S))) ≡ result_ring(T, S)
            @test result_ring(T, S, Real) ≡ Number
            for Z in TEST_TYPES
                @test typeof(test_ring(one(T), test_ring(one(S), one(Z)))) ≡ result_ring(T, S, Z)
            end
        end
    end
end

####
#### random tests
####

function all_combinations(ops, types)
    A = [op(one(T), one(S)) for T in types, S in types, op in ops]
    reduce(-, A)                # NOTE: avoid `sum` as it widens
end

RING_OPS = (+, -, *)
FIELD_OPS = (RING_OPS..., /)

randsub(x) = shuffle([x...])[1:rand(1:length(x))]

@testset "ring random tests" begin
    for _ in 1:10000
        Ts = randsub(TEST_TYPES)
        y = all_combinations(RING_OPS, Ts)
        T = result_ring(Ts...)
        @test y isa T
        y isa T || @info "test failure ring" Ts y T
    end
end

@testset "field random tests" begin
    for _ in 1:10000
        Ts = randsub(TEST_TYPES)
        y = all_combinations(FIELD_OPS, Ts)
        T = result_field(Ts...)
        @test y isa T
        y isa T || @info "test failure field" Ts y T
    end
end

####
#### utilities
####

@testset "number type" begin
    @test number_type(1) ≡ Int
    @test number_type([2.0, 4.0]) ≡ Float64
    @test number_type(Real[3.0]) ≡ Real
    @test_throws MethodError number_type(Any[1.0])
end


text:
# Deprecated functionality

@deprecate CuDevice(ctx::CuContext) device(ctx)
@deprecate CuCurrentDevice() current_device()
@deprecate CuCurrentContext() current_context()
@deprecate CuContext(ptr::Union{Ptr,CuPtr}) context(ptr)
@deprecate CuDevice(ptr::Union{Ptr,CuPtr}) device(ptr)


text:
using DifferentialEquations
using Plots
using Markdown

description = md"""

Solve the differential equation
``du/dt = f(u,p,t) = α u``
where t∈[0,1] and u₀= 1/2

"""
display(description)

α = 1.01
f(u, p, t) = α*u

# define IC's and time span
u₀ = 1/2
tspan = (0.0, 1.0)

prob = ODEProblem(f, u₀, tspan)
sol = solve(prob)  # type ?solve at repl to see list of solver options (methods, tolerance setting, etc...)

# visualize the solution
plot(sol, linewidth=5, title="Solution to the linear ODE", xaxis="Time (t)", yaxis="u(t) (μm)", label="Numerical Solution")
plot!(sol.t, t->0.5*exp(1.01t), lw=3, ls=:dash, label="True Solution")

max_stars_repo_name: john-waczak/SciML-examples

max_issues_repo_name: john-waczak/SciML-examples

max_forks_repo_name: john-waczak/SciML-examples


text:
include("white_list.jl")

const THRESHOLD = 30
const MIN_READ_SUPPORT = 2
const MAX_CLUSTER_DISTANCE = 100
const MATCH_EDIT_DISTANCE_LIMIT = 5

const FUSION_ON_MERGED_READ = 0
const FUSION_ON_READ1 = 1
const FUSION_ON_READ2 = 2
const FUSION_ON_CROSS_READS = 3

# ref_path is a folder contains fasta files by chromosomes
# like chr1.fa, chr2.fa ...
function detect(ref_path::AbstractString, bed_file::AbstractString, r1fq::AbstractString, r2fq::AbstractString="")
    index = index_bed(ref_path, bed_file)
    panel = index["panel"]
    panel_seq = index["seq"]
    panel_kmer_coord = index["kmer_coord"]
    ref_kmer_coords = index["ref_kmer_coords"]
    # pair end sequencing
    if r2fq != ""
        io = fastq_open_pair(r1fq, r2fq)
        fusion_pairs = Dict()
        while (pair = fastq_read_pair(io)) != false
            matches = detect_pair(panel_kmer_coord, pair)
            if length(matches)>1
                fusion_left, fusion_right, fusion_site, seg_result = verify_fusion_pair(ref_kmer_coords, panel_kmer_coord, panel, panel_seq, pair)
                if fusion_left!=false && fusion_right!=false && distance(fusion_left, fusion_right)>1000
                    # drop the low quality reads
                    if pass_quality_filter(pair.read1.quality) && pass_quality_filter(pair.read2.quality)
                        add_to_fusion_pair(fusion_pairs, fusion_left, fusion_right, fusion_site, seg_result, pair)
                    end
                end
            end
        end
        printed = print_fusion_pair(fusion_pairs, panel_seq, panel)
        if printed == 0
            println("# no fusion detected")
        end
    else
        error("FusionDirect only supports pair-end sequencing data")
    end
end

export pass_quality_filter
function pass_quality_filter(quality)
    # '0' = q15 as threshold
    const QUAL_T = '0'
    lowqual = 0
    for q in quality.qual
        if q<QUAL_T
            lowqual += 1
        end
    end
    # allow up to 8% low quality bases
    if lowqual/length(quality) > 0.08
        return false
    else
        return true
    end
end

function print_fusion_pair(fusion_pairs, panel_seq, panel)
    printed = 0
    gencode = gencode_load("GRCh37")
    for (fusion_key, all_fusion_reads) in fusion_pairs
        contig1, contig2 = fusion_key
        name1 = panel[contig1]["name"]
        name2 = panel[contig2]["name"]

        # clustering of fusion pairs
        clustered_fusion_reads = fusion_clustering(all_fusion_reads, panel)

        for fusion_reads in clustered_fusion_reads
            unique_fusion_reads, read_support = get_unique_fusion_pairs(fusion_reads)
            total_num = length(fusion_reads)
            unique_num = length(unique_fusion_reads)

            # filter out those fusion with very few reads support (like only 1 unique reads) or have no consistent fusion readsd
            if unique_num< MIN_READ_SUPPORT || !have_consistency(unique_fusion_reads)
                # check if this fusion is in the white list
                if !is_important_fusion(unique_fusion_reads, panel, gencode)
                    continue
                end
            end

            # give the fusion as a fasta comment line
            println("#Fusion:$name1-$name2 (total: $total_num, unique: $unique_num)")
            # display all reads support this fusion
            for read in unique_fusion_reads
                fusion_left, fusion_right, fusion_site, conjunct, pair = read
                support = read_support[read]
                name = ">$support\_"
                name = name * get_fusion_site_string(conjunct, fusion_site, pair) * "_"
                name = name * get_gene_coord_string(fusion_left, panel, gencode) * "_"
                name = name * get_gene_coord_string(fusion_right, panel, gencode) * "/"
                print(name, "1\n",pair.read1.sequence.seq,"\n")
                # println("# ", pair.read1.quality.qual)
                print(name, "2\n",pair.read2.sequence.seq,"\n")
                # println("# ", pair.read2.quality.qual)

                if fusion_site == FUSION_ON_MERGED_READ
                    offset, overlap_len, distance = overlap(pair)
                    merged_seq = simple_merge(pair.read1.sequence, pair.read2.sequence, overlap_len, offset)
                    print(name, "merged\n",merged_seq.seq,"\n")
                end

                printed += 1
            end
        end
    end
    return printed
end

function is_important_fusion(fusion_reads, panel, gencode)
    if length(IMPORTANT_FUSIONS) == 0
        return false
    end
    for read in fusion_reads
        fusion_left, fusion_right, fusion_site, conjunct, pair = read

        find1 = get_coord_gene(panel, fusion_left, gencode)
        find2 = get_coord_gene(panel, fusion_right, gencode)

        if find1[1] == false || find2[1] == false
            continue
        end

        for (fusion1, fusion2) in IMPORTANT_FUSIONS
            if is_fusion_identical(find1, fusion1) && is_fusion_identical(find2, fusion2)
                return true
            elseif is_fusion_identical(find1, fusion2) && is_fusion_identical(find2, fusion1)
                return true
            end
        end
    end
    return false
end

function is_fusion_identical(find, fusion)
    if find[1] != fusion[1]
        return false
    elseif fusion[2] != "*" && find[2] != fusion[2]
        return false
    elseif fusion[3] != "*" && find[3] != parse(Int, fusion[3])
        return false
    else
        return true
    end
end

function get_coord_gene(panel, coord, gencode)
    genename = panel[coord.contig]["name"]
    from = panel[coord.contig]["from"]
    chr = panel[coord.contig]["chr"]
    pos = from + abs(coord.pos)

    str = genename
    genelocs = gencode_locate(gencode, chr, pos)
    if length(genelocs) > 0
        for loc in genelocs
            if genename == loc["gene"]
                kind = loc["type"]
                number = loc["number"]
                return genename, kind, number
            end
        end
    end
    return false, false, false
end

function have_consistency(unique_fusion_reads)
    for i in 1:length(unique_fusion_reads) - 1
        for j in i+1:length(unique_fusion_reads)
            fusion_left1, fusion_right1, fusion_site1, conjunct1, pair1 = unique_fusion_reads[i]
            fusion_left2, fusion_right2, fusion_site2, conjunct2, pair2 = unique_fusion_reads[j]
            if fusion_site1 == FUSION_ON_CROSS_READS || fusion_site2 == FUSION_ON_CROSS_READS
                if is_bad_pair(pair1, pair2)
                    continue
                else
                    return true
                end
            else
                r11, r12 = get_fusion_seqs(unique_fusion_reads[i])
                r21, r22 = get_fusion_seqs(unique_fusion_reads[j])
                if is_consistent(r11, r21) && is_consistent(r12, r22)
                    return true
                elseif is_consistent(r11, r22) && is_consistent(r12, r21)
                    return true
                end
            end
        end
    end
    return false
end

function is_bad_pair(pair1, pair2)
    # filter the case that one read is found as dup, but the other read has lots of wrong bases
    dup, ed = is_dup(pair1.read1.sequence.seq, pair2.read1.sequence.seq)
    if dup
        if hamming(pair1.read2.sequence.seq[1:20], pair2.read2.sequence.seq[1:20]) <= 5
            return true
        end
    end
    dup, ed = is_dup(pair1.read1.sequence.seq, pair2.read2.sequence.seq)
    if dup
        if hamming(pair1.read2.sequence.seq[1:20], pair2.read1.sequence.seq[1:20]) <= 5
            return true
        end
    end
    return false
end

function is_consistent(r1, r2)
    const CONSISTENT_T = 5
    overlap_len = min(length(r1), length(r2))
    overlap_r1 = r1[1:overlap_len]
    overlap_r2 = r2[1:overlap_len]
    ed = edit_distance(overlap_r1.seq, overlap_r2.seq)
    return ed <= CONSISTENT_T
end

function get_fusion_seqs(fusion_read)
    fusion_left, fusion_right, fusion_site, conjunct, pair = fusion_read
    seq = nothing
    if fusion_site == FUSION_ON_READ1
        seq = pair.read1.sequence
    elseif fusion_site == FUSION_ON_READ2
        seq = pair.read2.sequence
    else
        offset, overlap_len, distance = overlap(pair)
        seq = simple_merge(pair.read1.sequence, pair.read2.sequence, overlap_len, offset)
    end
    len = length(seq)
    return ~(seq[1:conjunct]), seq[conjunct:len]
end

function get_gene_coord_string(coord, panel, gencode)
    genename = panel[coord.contig]["name"]
    from = panel[coord.contig]["from"]
    chr = panel[coord.contig]["chr"]
    pos = from + abs(coord.pos)

    str = genename
    genelocs = gencode_locate(gencode, chr, pos)
    if length(genelocs) > 0
        for loc in genelocs
            if genename == loc["gene"]
                kind = loc["type"]
                number = loc["number"]
                str *= ":$kind"
                str *= ":$number"
                break
            end
        end
    end
    chrloc = chr * ":" * string(pos)
    str *= "|" * strand_name(coord) * chrloc
end

function fusion_clustering(all_fusion_reads, panel)
    const FUSION_THRESHOLD = 100
    fusion_clusters = []
    for r in all_fusion_reads
        found = false
        for cluster in fusion_clusters
            for read in cluster
                left = read[1]
                right = read[2]
                left1 = r[1]
                right1 = r[2]
                if (abs_distance(left,left1)< MAX_CLUSTER_DISTANCE && abs_distance(right, right1)<MAX_CLUSTER_DISTANCE) || (abs_distance(left, right1)< MAX_CLUSTER_DISTANCE && abs_distance(left1, right)<MAX_CLUSTER_DISTANCE)
                    found = true
                    push!(cluster, r)
                    break
                end
            end
            if found == true
                break
            end
        end
        if found == false
            new_cluster = [r]
            push!(fusion_clusters, new_cluster)
        end
    end
    return fusion_clusters
end

function get_fusion_site_string(conjunct, fusion_site, pair)
    if fusion_site == FUSION_ON_MERGED_READ
        return "merged_$conjunct"
    elseif fusion_site == FUSION_ON_READ1
        return "read1_$conjunct"
    elseif fusion_site == FUSION_ON_READ2
        return "read2_$conjunct"
    else
        return "crosspair_0"
    end
end

function get_unique_fusion_pairs(fusion_reads)
    unique_fusion_reads=[]
    read_support = Dict()
    # display all reads support this fusion
    for read in fusion_reads
        fusion_left, fusion_right, fusion_site, conjunct, pair = read
        unique = true
        for uread in unique_fusion_reads
            uleft, uright, usite, uconjunct, upair = uread
            if is_dup_pair(upair, pair)
                unique = false
                read_support[uread] += 1
                break
            end
        end
        if unique
            push!(unique_fusion_reads, read)
            read_support[read] = 1
        end
    end
    return unique_fusion_reads, read_support
end

function is_dup_pair(pair1::FastqPair, pair2::FastqPair)
    # edit distance threshold
    const ED_T = length(pair1.read1) * 0.1

    seq11 = pair1.read1.sequence.seq
    seq12 = pair1.read2.sequence.seq
    seq21 = pair2.read1.sequence.seq
    seq22 = pair2.read2.sequence.seq

    dup1, ed1 = is_dup(seq11, seq21)
    dup2, ed2 = is_dup(seq12, seq22)
    if dup1 == true && (dup2 == true || ed2 <= ED_T || share_start_bases(seq12, seq22, 5))
        return true
    elseif dup2 == true && (dup1 == true || ed1 <= ED_T || share_start_bases(seq11, seq21, 5))
        return true
    end
    dup1, ed1 = is_dup(seq11, seq22)
    dup2, ed2 = is_dup(seq12, seq21)
    if dup1 == true && (dup2 == true || ed2 <= ED_T || share_start_bases(seq12, seq21, 5))
        return true
    elseif dup2 == true && (dup1 == true || ed1 <= ED_T || share_start_bases(seq11, seq22, 5))
        return true
    end

    # a work around for bad sequence quality
    if share_start_bases(seq11, seq21, 10) && share_start_bases(seq11, seq21, 10)
        return true
    elseif share_start_bases(seq12, seq21, 10) && share_start_bases(seq11, seq22, 10)
        return true
    end

    return false
end

function share_start_bases(s1::ASCIIString, s2::ASCIIString, len = 5)
    return hamming(s1[1:len], s2[1:len]) <= 1
end

function is_dup(s1::ASCIIString, s2::ASCIIString)
    ed = edit_distance(s1, s2)
    hm = hamming(s1, s2)
    if ed <= 3
        return true, ed
    elseif hm <= length(s1) * 0.1
        return true, ed
    elseif edit_distance(s1, s2) <=5 && ( hamming(s1[1:5], s2[1:5])<=1 || hamming(s1[length(s1)-4:length(s1)], s2[length(s2)-4:length(s2)])<=1 )
        return true, ed
    else
        return false, ed
    end
end

# hamming distance
function hamming(s1::ASCIIString, s2::ASCIIString)
    d = 0
    len = min(length(s1), length(s2))
    for i in 1:len
        if s1[i] != s2[i]
            d += 1
        end
    end
    return d
end

function strand_name(coord)
    return coord.pos>0?"+":"-"
end

function add_to_fusion_pair(fusion_pairs, fusion_left, fusion_right, fusion_site, seg_result, pair)
    key = (min(fusion_left.contig, fusion_right.contig), max(fusion_left.contig, fusion_right.contig))
    if (key in keys(fusion_pairs)) == false
        fusion_pairs[key]=[]
    end
    push!(fusion_pairs[key], (fusion_left, fusion_right, fusion_site, seg_result, pair))
end

function verify_fusion_se(ref_kmer_coords, panel_kmer_coord, panel, panel_seq, sequence)
    seg1, coords1 = segment(ref_kmer_coords, panel_kmer_coord, sequence, panel_seq)
    if length(seg1) > 1
        return make_connected_fusion(seg1, coords1)
    end
end

function found_seq_near_coord(panel_seq, seq, coord, pos_in_seq, threshold = 10)
    search_pos = abs(coord.pos -  pos_in_seq)
    seqlen = length(seq)
    # in reverse strand
    if coord.pos < 0
        seq = ~seq
        search_pos -= seqlen
    end
    const SEARCH_WINDOW = 10

    geneseq = panel_seq[coord.contig]
    genelen = length(geneseq)

    for start = max(search_pos-SEARCH_WINDOW, 1) : min(search_pos+SEARCH_WINDOW, genelen-seqlen-SEARCH_WINDOW)
        genepart = uppercase(geneseq.seq[start : start+seqlen-1])
        ed = edit_distance(genepart, seq.seq)
        if ed < threshold
            return true
        end
    end
    return false
end

# check if the splitted two sequences are consistent with ref
function splitted_consistent_with_ref(panel_seq, seq, fusion_left, fusion_right, left, right, conjunct)
    const CONJUNCT_WINDOW = 5
    seqlen = length(seq)

    #println("splitted_consistent_with_ref", " ", conjunct)
    #println(seq)
    for conj = -CONJUNCT_WINDOW:CONJUNCT_WINDOW
        break_point = conjunct + conj
        left_seq = seq[left:break_point]
        right_seq = seq[break_point+1:right]
        # coord of the starts of these two seqs
        left_coord = Coord(fusion_left.contig, fusion_left.pos - (conjunct-left))
        right_coord = Coord(fusion_right.contig, fusion_right.pos + conj)
        if found_seq_near_coord(panel_seq, left_seq, left_coord, 1, MATCH_EDIT_DISTANCE_LIMIT) && found_seq_near_coord(panel_seq, right_seq, right_coord, 1, MATCH_EDIT_DISTANCE_LIMIT)
            return true
        end
    end
    return false

end

function verify_fusion_pair(ref_kmer_coords, panel_kmer_coord, panel, panel_seq, pair)
    offset, overlap_len, distance = overlap(pair)
    # this pair is overlapped, so merged it and segment the merged sequence
    if overlap_len>0 && distance<5
        merged_seq = simple_merge(pair.read1.sequence, pair.read2.sequence, overlap_len, offset)
        seg_result, coords = segment(ref_kmer_coords, panel_kmer_coord, merged_seq, panel_seq)
        if length(seg_result) > 1
            if !is_seq_connected_on_ref(seg_result, ref_kmer_coords, merged_seq)
                fusion_left, fusion_right, conjunct = make_connected_fusion(seg_result, coords)
                l = seg_result[1][1]
                r = seg_result[2][2]
                if !found_seq_near_coord(panel_seq, merged_seq[l:r], fusion_left, conjunct-l+1) && !found_seq_near_coord(panel_seq, merged_seq[l:r], fusion_right, conjunct-l+1)
                    if splitted_consistent_with_ref(panel_seq, merged_seq, fusion_left, fusion_right, l, r, conjunct)
                        return (fusion_left, fusion_right, FUSION_ON_MERGED_READ, conjunct)
                    end
                end
            end
        end
    else
        seg_result1, coords1 = segment(ref_kmer_coords, panel_kmer_coord, pair.read1.sequence, panel_seq)
        if length(seg_result1) > 1
            if !is_seq_connected_on_ref(seg_result1, ref_kmer_coords, pair.read1.sequence)
                fusion_left, fusion_right, conjunct = make_connected_fusion(seg_result1, coords1)
                l = seg_result1[1][1]
                r = seg_result1[2][2]
                if !found_seq_near_coord(panel_seq, pair.read1.sequence[l:r], fusion_left, conjunct-l+1) && !found_seq_near_coord(panel_seq, pair.read1.sequence, fusion_right, conjunct-l+1)
                    if splitted_consistent_with_ref(panel_seq, pair.read1.sequence, fusion_left, fusion_right, l, r, conjunct)
                        return (fusion_left, fusion_right, FUSION_ON_READ1, conjunct)
                    end
                end
            end
        end
        seg_result2, coords2 = segment(ref_kmer_coords, panel_kmer_coord, pair.read2.sequence, panel_seq)
        if length(seg_result2) > 1
            if !is_seq_connected_on_ref(seg_result2, ref_kmer_coords, pair.read2.sequence)
                fusion_left, fusion_right, conjunct = make_connected_fusion(seg_result2, coords2)
                l = seg_result2[1][1]
                r = seg_result2[2][2]
                if !found_seq_near_coord(panel_seq, pair.read2.sequence[l:r], fusion_left, conjunct-l+1) && !found_seq_near_coord(panel_seq, pair.read2.sequence, fusion_right, conjunct-l+1)
                    if splitted_consistent_with_ref(panel_seq, pair.read2.sequence, fusion_left, fusion_right, l, r, conjunct)
                        return (fusion_left, fusion_right, FUSION_ON_READ2, conjunct)
                    end
                end
            end
        end
        if length(seg_result1)==0 || length(seg_result2) == 0 || seg_result1[1][1] == 0 || seg_result2[1][1] == 0
            return false, false, false, 0
        end
        # fusion of different contigs on a pair
        if coords1[seg_result1[1][1]].contig != coords2[seg_result2[1][1]].contig
            if is_pair_connected_on_ref(pair, ref_kmer_coords)
                return false, false, false, 0
            end

            l1 = seg_result1[1][1]
            l2 = seg_result2[1][1]
            len1 = length(coords1)
            len2= length(coords2)
            fusion_left = Coord(coords1[l1].contig, coords1[l1].pos + (len1 - l1))
            # reverse the direction of fusion right to align the pair in same direction
            fusion_right = Coord(coords2[l2].contig, coords2[l2].pos + (len2 - l2), -1)

            # check if consistent with ref
            left_coord = Coord(coords1[l1].contig, coords1[l1].pos - l1)
            right_coord = Coord(coords2[l2].contig, coords2[l2].pos - l2)
            if found_seq_near_coord(panel_seq, pair.read1.sequence, left_coord, 1, MATCH_EDIT_DISTANCE_LIMIT) && found_seq_near_coord(panel_seq, pair.read2.sequence, right_coord, 1, MATCH_EDIT_DISTANCE_LIMIT)
                return (fusion_left, fusion_right, FUSION_ON_CROSS_READS, 0)
            end
        end
    end
    return false, false, false, 0
end

function is_seq_connected_on_ref(seg_result, ref_kmer_coords, original_seq)
    l = seg_result[1][1]
    r = seg_result[2][2]
    seq = original_seq[l:r]
    coords_list = stat_ref(ref_kmer_coords, seq)
    clusters_on_ref = clustering_ref(coords_list)
    if length(clusters_on_ref) == 0
        # the cluster set is too big
        # so it is considered in a heavy repeating region
        # treat it as not fusion directly
        # println("seq not clustered")
        return true
    end
    for cluster in clusters_on_ref
        (left, right) = span_ref(cluster, length(seq))
        # nearly cover whole sequence
        if (right - left) > length(seq) - 20
            # println("whole seq on ref")
            return true
        end
        # nearly covers the two regions
        if left <  15 && right > length(seq) - 15
            # println("two regions on ref")
            return true
        end
    end
    # println("seq not on ref")
    return false
end

function is_pair_connected_on_ref(pair, ref_kmer_coords)
    coords_list1 = stat_ref(ref_kmer_coords, pair.read1.sequence)
    clusters1 = clustering_ref(coords_list1)
    coords_list2 = stat_ref(ref_kmer_coords, pair.read2.sequence)
    clusters2 = clustering_ref(coords_list2)
    # TODO
    # one or two reads are not clustered because the cluster set is too big
    # so they are considered in a heavy repeating region
    # treat it as not fusion directly
    if length(clusters1) == 0 || length(clusters2) == 0
        # println("pair not clustered")
        return true
    end

    # try to find if there is any intersection between clusters_on_ref1 and clusters_on_ref2
    for c1 in clusters1
        (left, right) = span_ref(c1, length(pair.read1.sequence))
        # for any cluster>30 in the read1
        if (right - left) > 30
            for c2 in clusters2
                # if we find a cluster from read2, which is >30 and near c1
                # it means they are connected
                (left2, right2) = span_ref(c2, length(pair.read2.sequence))
                if (right2 - left2) > 30 && min_distance(c1, c2) < 500
                    # two reads are close
                    # println("\npair too close, distance:", min_distance(c1, c2))
                    return true
                end
            end
        end
    end

    return false
end

function min_distance(c1, c2)
    dis = Inf
    for (pos1, coord1) in c1
        for (pos2, coord2) in c2
            # the pair should be in different strand
            d = abs(distance(coord1, Coord(coord2.contig, -coord2.pos)))
            if d < dis
                dis = d
            end
        end
    end
    return dis
end

function make_connected_fusion(seg_result, coords)
    l1 = seg_result[1][1]
    r1 = seg_result[1][2]
    l2 = seg_result[2][1]
    r2 = seg_result[2][2]
    conjunct = div(r1+l2, 2)
    #display_coords(coords)
    fusion_left = Coord(coords[l1].contig, coords[l1].pos + (conjunct - l1))
    fusion_right = Coord(coords[l2].contig, coords[l2].pos + (conjunct - l2))
    return fusion_left, fusion_right, conjunct
end

# detect fusion and find the segmentations
function segment(ref_kmer_coords, panel_kmer_coord::KmerCoord, seq::Sequence, panel_seq)
    counts, coords = stat(panel_kmer_coord, seq)
    clusters = clustering(coords)
    regions = Dict()
    for cluster in clusters
        # skip small clusters
        if length(cluster) < THRESHOLD
            continue
        end
        left, right = span(cluster, length(coords))
        regions[left]=right
    end
    # filter out those regions which are mostly inside other regions
    seg_result = filter_region(regions, coords)
    return seg_result, coords
end

# remove some inside regions
function filter_region(regions, coords)
    len = length(coords)
    #find the biggest region
    left = 0
    right = 0
    for (l, r) in regions
        if r-l>right-left
            left = l
            right = r
        end
    end
    clean_region = Dict(left=>right)
    # if some region is more or less a sub-region of the biggest region
    # then skip it
    for (l, r) in regions
        if r-right>10 || left-l>10
            clean_region[l]=r
        end
    end
    # sort the regions by the left
    sorted = sort(collect(clean_region))
    return sorted
end

# span the cluster on the panel_seq to find the region
function span(cluster, seqlen)
    left = seqlen
    right = 1
    for point in cluster
        left = min(left, point)
        right = max(right, point)
    end
    right = min(right+KMER-1, seqlen)
    return left, right
end

# span the cluster on the reference to find the region
function span_ref(cluster, seqlen)
    left = seqlen
    right = 1
    for point in cluster
        left = min(left, point.first)
        right = max(right, point.first)
    end
    right = min(right+KMER-1, seqlen)
    return left, right
end

# clustering by distance
function clustering(coords::Array{Coord, 1})
    points = Set(1:length(coords))
    clusters = []
    while !isempty(points)
        # create a new cluster
        # get a valid seed
        seed = pop!(points)
        while(!valid(coords[seed]) && !isempty(points))
            seed = pop!(points)
        end
        # we don't find a valid seed, so break
        if !valid(coords[seed])
            break
        end
        # create a new cluster
        cluster = Set(seed)
        # create a new working set
        working = Set(seed)
        while !isempty(working)
            current = pop!(working)
            for p in points
                # if it is not valid, just pop it out
                if !valid(coords[p])
                    pop!(points, p)
                elseif consistent(coords, current, p)
                    pop!(points, p)
                    push!(working, p)
                end
            end
            # current is done, push it to current cluster
            push!(cluster, current)
        end
        push!(clusters, cluster)
    end
    return clusters
end

function consistent(coords::Array{Coord, 1}, p1, p2)
    dis = coords[p2] - coords[p1]
    if dis > 1000
        return false
    end

    # the line of these two points be near 45 degree
    if abs(p2-p1)<40 && abs((p2-p1) - dis * sign(coords[p1].pos)) < 4
        return true
    end

    return false
end

# detect fusion in a pair of reads
function detect_pair(panel_kmer_coord::KmerCoord, pair::FastqPair)
    match1 = detect_read(panel_kmer_coord, pair.read1.sequence)
    match2 = detect_read(panel_kmer_coord, pair.read2.sequence)
    # merge them
    for m in match2
        if (m in match1) == false
            push!(match1, m)
        end
    end
    return match1
end

# simple fusion detection in a sequence
# this function is very fast and is not accurate
function detect_read(panel_kmer_coord::KmerCoord, seq::Sequence)
    counts, coords = stat(panel_kmer_coord, seq)
    matches = Array{Int16, 1}()
    for (k,v) in counts
        if v > THRESHOLD
            push!(matches, k)
        end
    end
    return matches
end

# get the coordinations on the panel for given sequence
function stat(panel_kmer_coord::KmerCoord, seq::Sequence)
    len = length(seq)
    coords = [unknown_coord() for i in 1:len]
    for i in 1:len-KMER+1
        kmer = seq[i:i+KMER-1]
        key = kmer2key(kmer)
        if haskey(panel_kmer_coord, key)
            coords[i] = panel_kmer_coord[key]
        end
    end

    counts = Dict{Int16, Int}()
    for coord in coords
        if valid(coord)
            contig = coord.contig
            if haskey(counts, contig)
                counts[contig] += 1
            else
                counts[contig] = 1
            end
        end
    end
    return counts, coords
end

# get the coordinations on the whole reference for given sequence
function stat_ref(ref_kmer_coords::KmerCoordList, seq::Sequence)
    len = length(seq)
    coord_lists = Dict()
    for i in 1:len-KMER+1
        kmer = seq[i:i+KMER-1]
        key = kmer2key(kmer)
        if haskey(ref_kmer_coords, key)
            coord_lists[i] = ref_kmer_coords[key]
        end
    end
    return coord_lists
end

# clustering of the kmer_coords
function clustering_ref(coord_lists)
    total = Set()
    # covert the coord list to a set of pos=>coord
    for (pos, arr) in coord_lists
        # if two many entries on this kmer
        # we consider this kmer_coords as less informative
        # for performance considering, we skip it
        if length(arr) < 300
            for coord in arr
                push!(total, pos=>coord)
            end
        end
    end

    clusters = []
    # println(length(total))
    # TODO: handle big cluster set
    # the cluster set is too big
    # so it is considered in a heavy repeating region
    # treat it as not fusion directly for performance considering
    if(length(total) > 10000)
        return clusters
    end
    while !isempty(total)
        # create a new cluster
        seed = pop!(total)
        # create a new cluster
        cluster = Set()
        push!(cluster, seed)
        # create a new working set
        working = Set()
        push!(working, seed)
        while !isempty(working)
            current = pop!(working)
            for p in total
                # if it is not valid, just pop it out
                if consistent_on_ref(current, p)
                    pop!(total, p)
                    push!(working, p)
                end
            end
            # current is done, push the pos to current cluster
            push!(cluster, current)
        end
        push!(clusters, cluster)
    end
    return clusters
end

function consistent_on_ref(p1, p2)
    dis = p2.second - p1.second
    if dis > 1000
        return false
    end

    # the line of these two points be near 45 degree
    if abs(p2.first-p1.first)<40 && abs((p2.first-p1.first) - dis * sign(p1.second.pos)) < 4
        return true
    end

    return false
end


text:
"""
```
brisk_params = BRISK([pattern_scale = 1.0])
```

| Argument | Type | Description |
| :--------| :--- | :---------- |
| `pattern_scale` | `Float64` | Scaling factor for the sampling window |
"""
mutable struct BRISK <: Params
    threshold::Float64
    octaves::Int
    pattern_scale::Float64
    pattern_table::Vector{Vector{SamplePair}}
    smoothing_table::Vector{Vector{Float16}}
    orientation_weights::Vector{OrientationWeights}
    short_pairs::Vector{OrientationPair}
    long_pairs::Vector{OrientationPair}
end

function BRISK(; threshold::Float64 = 0.25, octaves::Int = 4, pattern_scale = 1.0)
    pattern_table, smoothing_table = _brisk_tables(pattern_scale)
    orientation_weights = OrientationWeights[]
    short_pairs = OrientationPair[]
    long_pairs = OrientationPair[]
    dminsq = (brisk_dmin * pattern_scale) ^ 2
    dmaxsq = (brisk_dmax * pattern_scale) ^ 2
    for i in 2:brisk_points
        for j in 1:i - 1
            dy = pattern_table[1][j][1] - pattern_table[1][i][1]
            dx = pattern_table[1][j][2] - pattern_table[1][i][2]
            norm = dy ^ 2 + dx ^ 2
            if norm > dminsq
                push!(long_pairs, OrientationPair((j, i)))
                push!(orientation_weights, OrientationWeights((dy / norm, dx / norm)))
            elseif norm < dmaxsq
                push!(short_pairs, OrientationPair((j, i)))
            end
        end
    end
    BRISK(threshold, octaves, pattern_scale, pattern_table, smoothing_table, orientation_weights, short_pairs, long_pairs)
end

function _brisk_orientation(int_img::AbstractArray{T, 2}, keypoint::Keypoint, pattern::Array{SamplePair},
                               orientation_weights::Array{OrientationWeights}, sigmas::Array{Float16}, long_pairs::Array{OrientationPair}) where T<:Gray
    direction_sum_y = 0.0
    direction_sum_x = 0.0
    for (i, o) in enumerate(long_pairs)
        offset_1 = pattern[o[1]]
        offset_2 = pattern[o[2]]
        intensity_diff = ImageFeatures._freak_mean_intensity(int_img, keypoint, offset_1, sigmas[o[1]]) - ImageFeatures._freak_mean_intensity(int_img, keypoint, offset_2, sigmas[o[2]])
        direction_sum_y += orientation_weights[i][1] * intensity_diff
        direction_sum_x += orientation_weights[i][2] * intensity_diff
    end
    angle = atan(direction_sum_y, direction_sum_x)
    scaled_angle = ceil(Int, (angle + pi) * brisk_orientation_steps / (2 * pi))
    scaled_angle
end

function _brisk_tables(pattern_scale::Float64)
    pattern_table = Vector{SamplePair}[]
    smoothing_table = Vector{Float16}[]
    for ori in 0:brisk_orientation_steps - 1
        theta = ori * 2 * pi / brisk_orientation_steps
        pattern = SamplePair[]
        sigmas = Float16[]
        for (i, n) in enumerate(brisk_num_circular_pattern)
            for circle_number in 0:n - 1
                angle = (circle_number * 2 * pi / n) + theta

                push!(pattern, SamplePair((brisk_radii[i] * sin(angle) * pattern_scale * 0.85,
                                            brisk_radii[i] * cos(angle) * pattern_scale * 0.85)))
                push!(sigmas, brisk_sigma[i] * pattern_scale * 0.85)
            end
        end
        push!(pattern_table, pattern)
        push!(smoothing_table, sigmas)
    end
    pattern_table, smoothing_table
end

function create_descriptor(img::AbstractArray{T, 2}, features::Features, params::BRISK) where T<:Gray
    int_img = integral_image(img)
    descriptors = BitArray{1}[]
    ret_features = Feature[]
    window_size = ceil(Int, (brisk_radii[end] + brisk_sigma[end]) * params.pattern_scale * 0.85) + 1
    lim = CartesianIndex(window_size, window_size)
    for feature in features
        keypoint = Keypoint(feature)
        checkbounds(Bool, img, keypoint - lim) && checkbounds(Bool, img, keypoint + lim) || continue
        orientation = _brisk_orientation(int_img, keypoint, params.pattern_table[1], params.orientation_weights, params.smoothing_table[1], params.long_pairs)
        sampled_intensities = T[]
        for (i, p) in enumerate(params.pattern_table[orientation])
            push!(sampled_intensities, ImageFeatures._freak_mean_intensity(int_img, keypoint, p, params.smoothing_table[orientation][i]))
        end
        descriptor = falses(size(params.short_pairs, 1))
        for (i, f) in enumerate(params.short_pairs)
            point_1 = sampled_intensities[f[1]]
            point_2 = sampled_intensities[f[2]]
            descriptor[i] = point_1 < point_2
        end
        push!(descriptors, descriptor)
        push!(ret_features, Feature(keypoint, orientation))
    end
    descriptors, ret_features
end


text:
using LinearAlgebra
import SpecialFunctions: gamma

# Fractional Brownian motion
struct FBM
  t::Vector{Float64}
  n::Int64
  h::Float64 # Hurst index

  function FBM(t::Vector{Float64}, n::Int64, h::Float64)
    t[1] == 0. || error("Starting time point must be equal to 0.0.")
    issorted(t, lt=<=) || error("The time points must be strictly sorted.")
    length(t) == n || error("Number of time points must be equal to the vector holding the time points.")
    0 < h < 1 || error("Hurst index must be between 0 and 1.")
    new(t, n, h)
  end
end

FBM(t::Vector{Float64}, h::Float64) = FBM(t, Int64(length(t)), h)
FBM(t::AbstractRange, h::Float64) = FBM(collect(t), Int64(length(t)), h)
FBM(t::Float64, n::Int64, h::Float64) = FBM(0.0:t/n:t-t/n, h)
FBM(t::Float64, h::Float64) = FBM([t], 1, h)

FBM(t::Matrix{Float64}, h::Float64) = FBM[FBM(t[:, i], h) for i = 1:size(t, 2)]
FBM(t::AbstractRange, np::Int, h::Float64) = FBM[FBM(t, h) for i = 1:np]
FBM(t::Float64, n::Int64, np::Int, h::Float64) = FBM[FBM(t, n, h) for i = 1:np]

# Fractional Gaussian noise
struct FGN
  σ::Float64
  h::Float64 # Hurst index

  function FGN(σ::Float64, h::Float64)
    σ > 0. || error("Standard deviation must be positive.")
    0. < h < 1. || error("Hurst index must be between 0 and 1.")
    new(σ, h)
  end
end

FGN(h::Float64) = FGN(1., h)

convert(::Type{FBM}, p::BrownianMotion) = FBM(p.t, p.n, 0.5)
convert(::Type{FGN}, p::BrownianMotion) = FGN((p.t[end]/(p.n-1))^0.5, 0.5)
convert(::Type{FGN}, p::FBM) = FGN((p.t[end]/(p.n-1))^p.h, p.h)

function autocov(p::FGN, i::Int64, j::Int64)
  twoh::Float64 = 2*p.h
  0.5*abs2(p.σ)*(abs(j-i+1)^twoh+abs(j-i-1)^twoh-2*abs(j-i)^twoh)
end

function autocov!(c::Matrix{Float64}, p::FGN)
  n::Int64 = size(c, 1)

  for i = 1:n
    for j = 1:i
      c[i, j] = autocov(p, i, j)
    end
  end

  for i = 1:n
    for j = (i+1):n
      c[i, j] = c[j, i]
    end
  end

  c
end

autocov(p::FGN, maxlag::Int64) = autocov!(zeros(Float64, maxlag, maxlag), p)

function autocov!(y::Vector{Float64}, p::FGN, lags::IntegerVector)
  nlags = length(lags)
  twoh::Float64 = 2*p.h

  for i = 1:nlags
    y[i] = 0.5*abs2(p.σ)*(abs(lags[i]+1)^twoh+abs(lags[i]-1)^twoh-2*abs(lags[i])^twoh)
  end

  y
end

autocov(p::FGN, lags::IntegerVector) = autocov!(zeros(Float64, length(lags)), p, lags)

function autocov(p::FBM, i::Int64, j::Int64)
  twoh::Float64 = 2*p.h
  0.5*((p.t[i])^twoh+(p.t[j])^twoh-abs(p.t[i]-p.t[j])^twoh)
end

function autocov!(c::Matrix{Float64}, p::FBM)
  n::Int64 = p.n-1

  for i = 1:n
    for j = 1:i
      c[i, j] = autocov(p, i+1, j+1)
    end
  end

  for i = 1:n
    for j = (i+1):n
      c[i, j] = c[j, i]
    end
  end

  c
end

function autocov(p::FBM)
  n::Int64 = p.n-1
  autocov!(zeros(Float64, n, n), p)
end

"""
rand_chol

generates FBM using the method based on Cholesky decomposition.
T. Dieker, Simulation of Fractional Brownian Motion, master thesis, 2004.
The complexity of the algorithm is O(n^3), where n is the number of FBM samples.
"""
function rand_chol(p::FBM, fbm::Bool=true)

  w = (cholesky(autocov(p)).U') * randn(p.n-1)

  # If fbm is true return FBM, otherwise return FGN
  insert!(w, 1, 0.0)
  if !fbm
    w = diff(w)
  end

  w
end

"""
rand_fft

generates FBM using fast Fourier transform (FFT).
The time interval of FBM is [0, 1] with a stepsize of 2^p, where p is a natural number.
The algorithm is known as the Davies-Harte method or the method of circular embedding.
R.B. Davies and D.S. Harte, Tests for Hurst Effect, Biometrika, 74 (1), 1987, pp. 95-102.
For a more recent publication, see P.F. Craigmile, Simulating a Class of Stationary Gaussian Processes Using the
Davies–Harte Algorithm, With Application to Long Memory Processes, Journal of Time Series Analysis, 24 (5), 2003,
pp. 505-511.
The complexity of the algorithm is O(n*log(n)), where n=2^p is the number of FBM samples.
"""
function rand_fft(p::FBM, fbm::Bool=true)
  # Determine number of points of simulated FBM
  pnmone = p.n-1
  n = 1 << Int(ceil(log2(pnmone)))

  # Compute autocovariance sequence of underlying FGN
  c = zeros(Float64, n+1)
  autocov!(c, FGN((1/n)^p.h, p.h), 0:n)

  # Compute square root of eigenvalues of circular autocovariance sequence
  l = real(fft(vcat(1, c, c[end-1:-1:2])))
  all(i->(i>0), l) || error("Non-positive eigenvalues encountered.")
  lsqrt = sqrt.(l)

  # Simulate standard random normal variables
  twon::Int64 = 2n
  z = randn(twon)

  # Generate fractional Gaussian noise (retain only the first p.n-1 values)
  x = sqrt(0.5) .* lsqrt[2:n] .* complex.(z[2*(2:n) .- 2], z[2*(2:n) .- 1])
  w = real(bfft(vcat(1, lsqrt[1]*z[1], x, lsqrt[n+1]*z[twon], conj(reverse(x)))))[1:pnmone]/sqrt(twon)

  # If fbm is true return FBM, otherwise return FGN
  if fbm
    w = insert!(cumsum(w), 1, 0.0)
  end

  w
end

"""
    rand_rl

Implementation of "exact discrete" method for simulating Riemann-Liouville fBm
Based on Muniandy, S. & Lim, S. Modeling of locally self-similar processes using multifractional
Brownian motion of Riemann-Liouville type. Physical Review E 63, 046104. ISSN: 1063-651X (2001).
Specifically, Eqn. 17, 18 and 19.
From Muniandy (2001), Reimann-Liouville fBm is defined as,

```math
B_H(t)=\\frac{1}{\\Gamma(H+0.5)}\\int^{t}_{0}(t-s)^(H-0.5)dB(s), t\\geq0.
```

Here we focus on the discrete time ``t_j=j\\Delta t`` approximation,

```math
B_H(t_j)=\\frac{1}{\\Gamma(H+0.5)}\\sum^{j}_{i=1}\\int^{i\\Delta t}_{(i-1)\\Delta t}(t_j-\\tau)^(H-0.5)dB(\\tau)
```

The exact solution to the interior integral results in a weighting function, implemented here with the :Exact key.
An "improved" weighting function was proposed by
Rambaldi, S. & Pinazza, O. An accurate fractional Brownian motion generator. Physica A 208, 21–30 (1994).
We have implemented this here with the :Improved key.
"""
function rand_rl(p::FBM, fbm::Bool=true; wtype::Symbol=:exact)
  # We are going to require that the time points are evenly spaced to keep things simple
  dt = unique(diff(p.t))
  if length(dt) != 1
	  error("For this simulation technique, the fBm process must be sampled across a uniform temporal grid")
  end

  dt = dt[1]

  w_mat = zeros(p.n-1,p.n-1) # Prealocate a weight matrix. Note that this will be upper triangular

  for j=1:p.n-1
    for i=1:j
      w_mat[i,j] = weight(p.h, (j-i+1)*dt, dt, wtype)
    end
  end

  # Multiply a vector of white noise by the weight matrix and scale appropriately.
	X = dropdims(sqrt(2)*(randn(p.n-1)'*w_mat),dims=1).*sqrt(dt)

	insert!(X, 1, 0.0)  # Snap to zero. Not clear if this causes a discontinuity in the correlation structure.
  #It is worth considering alternative for the above line: generate path X of length p.n, then let X = X-X[1].

  return X
end

weight(H::Float64, t::Float64, dt::Float64, method::Symbol=:exact) = weight(H, t, dt, Val{method})

function weight(H::Float64, t::Float64, dt::Float64, ::Type{Val{:exact}})
	nom = t^(H+.5)-(t-dt)^(H+.5)
	denom = gamma(H+.5)*(H+.5)*dt
	return nom/denom
end

function weight(H::Float64, t::Float64, dt::Float64, ::Type{Val{:improved}})
	nom = t^(2*H)-(t-dt)^(2*H)
	denom = 2*H*dt
	return sqrt(nom/denom)/gamma(H+.5)
end

### API functions

function rand(p::FBM; fbm::Bool=true, method::Symbol=:fft, args...)
  if method == :fft
    rand_fft(p, fbm)
  elseif method == :rl
    rand_rl(p, fbm; args...)
  elseif method == :chol
    rand_chol(p, fbm)
  else
    error("Non-existing method $method")
  end
end

function rand(p::Vector{FBM}; fbm::Bool=true, method::Symbol=:fft, args...)
  n::Int64 = fbm ? p[1].n : p[1].n-1
  np = length(p)
  x = zeros(Float64, n, np)

  for i = 1:np
    x[:, i] = rand(p[i]; fbm=fbm, method=method, args...)
  end

  x
end

function chol_update(r0::Array{Float64,2},A::Array{Float64,2})
# r0 is a lower triangular matrix found via cholesky decomposition of
# A(1:end-1,1:end-1), therefore r0 is a matrix of size (N-1)-by-(N-1) when A
# is of size N-by-N. A is a positive symetric semidefinite matrix by construction.

    m = size(A,1); # Assumes A is square
    r1 = zeros(m,m)
    r1[1:end-1,1:end-1] = r0

  @inbounds  for i=1:m-1
        r1[m,i] =  (1/r1[i,i])*(A[m,i]-sum(r1[m,1:i].*r1[i,1:i]))[1]
    end

    r1[m,m] = sqrt(A[m,m]-sum(r1[m,1:m-1] .* r1[m,1:m-1]))[1]

    return r1
end


text:
"""
    const INTERFACE

Methods comprising the package's interface.
"""
const INTERFACE = Vector{Symbol}()

# Maybe just write the text into a file instead of using a macro here?
"""
    macro import_interface()

Import all symbols inside [`INTERFACE`](@ref).
"""
macro import_interface()
    ex = Expr(:block)
    ex.args = [:(import WavePropBase: $f) for f in INTERFACE]
    return ex
end

"""
    macro export_interface()

Export all symbols inside [`INTERFACE`](@ref).
"""
macro export_interface()
    ex = Expr(:block)
    ex.args = [:(export $f) for f in INTERFACE]
    return ex
end

"""
    ambient_dimension(x)

Dimension of the ambient space where `x` lives. For geometrical objects this can
differ from its [`geometric_dimension`](@ref); for example a triangle in `ℝ³` has
ambient dimension `3` but geometric dimension `2`, while a curve in `ℝ³` has
ambient dimension 3 but geometric dimension 1.
"""
function ambient_dimension end
push!(INTERFACE,:ambient_dimension)

"""
    geometric_dimension(x)
    geometric_dimension(Ω::Domain)

Number of degrees of freedom necessary to locally represent the geometrical
object. For example, lines have geometric dimension of 1 (whether in `ℝ²` or in
`ℝ³`), while surfaces have geometric dimension of 2.

When the argument is a `Domain`, return the largest geometric dimension
encoutered.
"""
function geometric_dimension end
push!(INTERFACE,:geometric_dimension)

"""
    dimension(space)

The length of a basis for `space`; i.e. the number of linearly independent elements
required to span `space`.
"""
function dimension end
push!(INTERFACE,:dimension)

"""
    boundary(ω)

Return the boundary of `ω`. For a mesh element gives the `d-1` dimensional
elements composing its boundary, while for an entity gives the corresponding
`d-1` dimensional entities.
"""
function boundary end
push!(INTERFACE,:boundary)

"""
    diameter(Ω)

Largest distance between `x` and `y` for `x,y ∈ Ω`.
"""
function diameter end
push!(INTERFACE,:diameter)

"""
    distance(Ω₁,Ω₂)

Minimal Euclidean distance between `Ω₁` and `Ω₂`.
"""
function distance end
push!(INTERFACE,:distance)

"""
    radius(Ω)

Half the [`diameter`](@ref).
"""
function radius end
push!(INTERFACE,:radius)

"""
    center(Ω)

Center of the smallest possible ball containing `Ω`.
"""
function center end
push!(INTERFACE,:center)

"""
    return_type(f)

The type returned by function-like objects.
"""
function return_type end
push!(INTERFACE,:return_type)

"""
    jacobian(f,x)

Given a (possibly vector-valued) function `f : 𝐑ᵐ → 𝐅ᵐ`, return the `m × n`
matrix `Aᵢⱼ = ∂fᵢ/∂xⱼ`.
"""
function jacobian end
push!(INTERFACE,:jacobian)

"""
    normal(el,x̂)
    normal(jac::SMatrix)

The unit normal vector at coordinate `x̂`, guaranteed to be orthogonal to all
columns of `jacobian(el,x)`.
"""
function normal end
push!(INTERFACE,:normal)

"""
    domain(f)

The domain of the function `f`; i.e. if `f: Ω → R`, return `Ω`.
"""
function domain end
push!(INTERFACE,:domain)

"""
    image(f)

The image of the function `f`; i.e. if `f: Ω → R`, return `R`.
"""
function image end
push!(INTERFACE,:image)

"""
    parametrization(el)

Return the underlying parametrization of `el`.
"""
function parametrization end
push!(INTERFACE,:parametrization)

function key end
push!(INTERFACE,:key)

"""
    coords(x)

Return an `SVector` giving a cartesian coordinate used to sort `x`. For points,
this is simply the coordinates of the point, while for elements this may be the
center of the element. You must overload this function for your own type if you
want e.g. the clustering algorithms to work for `x`.
"""
function coords end
push!(INTERFACE,:coords)

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
for i in 1:100
    msg = "Fizz" ^ (i % 3 == 0) * "Buzz" ^ (i % 5 == 0)
    println(isempty(msg) ? i : msg)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using BrainSim
using Documenter

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

makedocs(;
    modules=[BrainSim],
    authors="terrypang <terrypang@aliyun.com> and contributors",
    repo="https://github.com/THBI-BrainSim/BrainSim.jl/blob/{commit}{path}#{line}",
    sitename="BrainSim.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://thbi-brainsim.github.io/BrainSim.jl/dev/",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/THBI-BrainSim/BrainSim.jl",
)


text:
using Pkg
Pkg.API.develop(Pkg.PackageSpec(name="MrPhelps", path="/home/caseykneale/Desktop/Playground/MrPhelps/"))
using MrPhelps, Distributed, Dates

localonly = true
#                        Connect Machines!
#First a user needs to define machines it has access too...
MrPhelps.greet()
!localonly && run(`scp -r /home/caseykneale/Desktop/Playground optics@192.168.0.14:/home/optics/`)
#Define NodeManager - to keep track of Distributed Hooks to all machines
nm = NodeManager()
Local  = "127.0.0.1"
Remote = "optics@192.168.0.14"
#Connect to a remote box and use a core/thread on the local machine!
LocalNode   = @async addprocs( 2; restrict = true )
if !localonly
    RemoteNode  = @async addprocs( [ ( Remote, 2 ) ], tunnel = true, max_parallel = 4,
                                    exename = "/home/optics/julia-1.3.0/bin/julia") #, sshflags = "-vvv" )
end

@everywhere begin
    using Pkg
    pth = "/home/caseykneale/Desktop/Playground/MrPhelps/"
    pth = isdir( "/home/caseykneale/" ) ? pth : "/home/optics/Playground/MrPhelps/"
    Pkg.API.develop(Pkg.PackageSpec(name="MrPhelps", path=pth))
end
@everywhere using MrPhelps
#Update our node manager, because we've added connections
#Or we could make the node manager after adding the connections
update!(nm)
#Grok the machines we have available to our node manager...
println( nm.machinenodemap )

#                        Define Some Tasks
mission = MissionGraph()
#Add an unconnected node to the graph
data_file_fomat = "data_{year}-{month}-{day}.csv"
expansion = Dict( "year"  => [2020], "month" => [7,11], "day" => collect(1:30) )

add_node!(mission, StashIterator( Expand(data_file_fomat, expansion),
                                 @thunk( string ), [ Local ], 1 ) )
# add_node!(mission, Stash(   "/home/caseykneale/Desktop/megacsv.csv",
#                             @thunk( string ), [ Local ], 1 ) )
#Add a new node to the graph but connect it to the last node laid down
attach_node!(mission, Agent( @thunk( uppercase ), [ Local ] ) )
#Add another new node, but give it a bookmark so we can find it later!
attach_node!(mission, :prod => Agent( @thunk( lowercase ), [ Local ] ) )
#Look we can add another new node to the graph unattached to anything
add_node!(mission, :final => Agent( @thunk( println ), [Local] ) )
#And now we can connect it to something else we bookmarked!
connect!(mission, :prod, :final)
#We made a very simple linear chain. Yay!
add_node!(mission, :references => Stash("/home/caseykneale/Desktop/refcsv.csv",
                                    @thunk( string ), [ Local ], 1 ) )
connect!( mission, :references, :prod )

#begin tasking...
sc = Scheduler( nm, mission )
execute_mission( sc )

@async spawn_listeners( sc )

mission.meta[1]
next!(mission.meta[1] )

#for debugging
#sc.worker_communications
isready(sc.worker_communications[2])
isready(sc.worker_communications[3])
isready(sc.worker_channels[2])
isready(sc.worker_channels[3])

# Testing StashIterators
# mission = MissionGraph()
# data_file_fomat = "data_{year}-{month}-{day}.csv"
# expansion = Dict(   "year"  => [2020],
#                     "month" => [7,11],
#                     "day"   => collect(1:30))
#
# add_node!(mission, StashIterator( Expand(data_file_fomat, expansion),
#                                  @thunk( string ), [ Local ], 1 ) )
# attach_node!(mission, Agent( @thunk( uppercase ), [ Local ] ) )
# attach_node!(mission, :prod => Agent( @thunk( lowercase ), [ Local ] ) )
# add_node!(mission, :final => Agent( @thunk( println ), [Local] ) )
# connect!(mission, :prod, :final)
# add_node!(mission, :references => Stash("/home/caseykneale/Desktop/refcsv.csv",
#                                     @thunk( string ), [ Local ], 1 ) )
# connect!( mission, :references, :prod )
##########################################################################
# Most naive scheme
#   Worker starts process at a source
#   Finishes process, alerts, scheduler, and gets next task (stay or move)
#       After each alert, the scheduler decides what to do next.
#       Are there available nodes for open tasks?
#       Do we need to wait for a worker to open to complete a task?
##########################################################################
# ===============================================
#               Below is all WIP
# ===============================================
# Graph     <->     Tasks       <->     Machine
# Node      <->     Machine
# Assign:
#   - Find parent task
#   - Find nodes to meet requirements
# Execute:
#   - (async) Dispatch task thunks to nodes and run
#   - (async) Result is stored to a channel & statistics sent back to local
#   - Local decides the next task for the jobs available
# More fleshed out...
#   - Every worker gets a communication remotechannel!
#   - Dispatch thunks
#   - Thunks wrapped in a function which
#       :wraps thunks and deploys args
#       :computes runtime stats
#       :ferries results to remote channel

macro ift(ex...)
    @assert( length(ex) == 3, "malformed if/then statement.")
    @assert( ex[ 2 ] == :then, "malformed if/then statement.")
    quote
        local conditional   = $(esc( ex[1] ))
        local trueresponse  = $(esc( ex[3] ))
        return conditional  ? trueresponse : nothing
    end
end

#abc() = "dfg"
#innerfun1() = @ift 1==1 then abc()
#innerfun3(a, b) = @ift a ==1 then b

max_stars_repo_name: caseykneale/MrPhelps

max_issues_repo_name: caseykneale/MrPhelps

max_forks_repo_name: caseykneale/MrPhelps


text:
using Mocha

backend = CPUBackend()
init(backend)

include("mnist-layers.jl")

train_net = Net("train-net", backend, [train_dl, non_data_layers...])
test_net = Net("test-net", backend, [test_dl, non_data_layers...])

using JLD
base_dir = "snapshots_sun"
jdl_file = jldopen(joinpath(base_dir, "snapshot-150000.jld"))
load_network(jdl_file, train_net)
close(jdl_file)

forward(train_net)
backward(train_net)

forward(test_net)
backward(test_net)

data = Array(Float32, 28, 28, 100)
copy!(data, get_layer_state(train_net, "train-data").blobs[1])
recon = Array(Float32, 28, 28, 100)
copy!(recon, get_layer_state(train_net, "x_mean").blobs[1])
open("train_datarecon", "w") do f
    serialize(f, (data, recon))
end

data2 = Array(Float32, 28, 28, 100)
copy!(data2, get_layer_state(test_net, "test-data").blobs[1])
recon2 = Array(Float32, 28, 28, 100)
copy!(recon2, get_layer_state(test_net, "x_mean").blobs[1])
open("test_datarecon", "w") do f
    serialize(f, (data2, recon2))
end


text:
using UnregisteredWithProject
using Test
using Example

@test f(2) == 2


text:
#---------------------------------------------------------------------
#	Biconnected components
#---------------------------------------------------------------------

"""
	verts2verts(EV::Lar.Cells)::Lar.Cells

Adjacency lists of vertices of a cellular 1-complex.

# Example

```julia
julia> V,(VV,EV,FV) = Lar.cuboidGrid([3,3],true);

julia> verts2verts(EV::Lar.Cells)::Lar.Cells
16-element Array{Array{Int64,1},1}:
 [2, 5]
 [1, 3, 6]
 [2, 4, 7]
 [3, 8]
 [1, 6, 9]
 [2, 5, 7, 10]
 [3, 6, 8, 11]
 [4, 7, 12]
 [5, 10, 13]
 [6, 9, 11, 14]
 [7, 10, 12, 15]
 [8, 11, 16]
 [9, 14]
 [10, 13, 15]
 [11, 14, 16]
 [12, 15]
```
"""
function verts2verts(EV::Lar.Cells)::Lar.Cells
    cscEV = Lar.characteristicMatrix(EV)
    cscVE = convert(Lar.ChainOp, cscEV')
    cscVV = cscVE * cscEV
    rows,cols,data = findnz(cscVV)
    VV = [ findnz(cscVV[k,:])[1] for k=1:size(cscVV,2) ]
    return [setdiff(vs,[k]) for (k,vs) in enumerate(VV)] # remove reflexive pairs
end

"""
	DFV_visit( VV::cells, out::Array, count::Int, visited::Array, parent::Array, d::Array, low::Array, stack::Array, u::Int )::Array

Hopcroft-Tarjan algorithm. Depth-First-Visit recursive algorithm.
"""
function DFV_visit( VV::Lar.Cells, out::Array, count::Int, visited::Array, parent::Array, d::Array, low::Array, stack::Array, u::Int )::Array

    visited[u] = true
    count += 1
    d[u] = count
    low[u] = d[u]
    for v in VV[u]
        if ! visited[v]
            push!(stack, [(u,v)])
            parent[v] = u
            DFV_visit( VV,out,count,visited,parent,d,low,stack, v )
            if low[v] >= d[u]
                push!(out, [outputComp(stack,u,v)])
            end
            low[u] = min( low[u], low[v] )
        else
            if ! (parent[u]==v) && (d[v] < d[u])
                push!(stack, [(u,v)])
            end
            low[u] = min( low[u], d[v] )
        end
    end
    out
end

visited=[]; count=0; out=[]; parent=[]; d=[];low=[]; stack=[]; u=1

"""
	outputComp(stack::Array, u::Int, v::Int)::Array

Internal output of biconnected components. Utility function for DFV_visit graph algorithm (Hopcrof, Tarjan (1973)).
"""
function outputComp(stack::Array, u::Int, v::Int)::Array
    out = []
    while true
        e = pop!(stack)[1]
        push!(out,e)
        if e == (u,v)
        	break
        end
    end
    return [out]
end

"""
	biconnectedComponent(model)

Main procedure for computation of biconnected components of a graph
represented by a LAR unsigned pair (Points, Cells).

# Example

```julia
julia> V =  [3.0  6.0  6.0  6.0  6.0  0.0  3.0  10.0  10.0  10.0  0.0  3.0  0.0;
 			 2.0  2.0  0.0  4.0  8.0  4.0  4.0   4.0   0.0   8.0  8.0  0.0  0.0]

julia> EV = [[1, 2], [2, 3], [2, 4], [5, 10], [3, 12], [5, 11], [6, 7], [8, 9],
			 [6, 13], [12, 13], [1, 7], [1, 12], [4, 7], [4, 8], [4, 5], [8, 10]]

julia> VV = [[k] for k=1:size(V,2)]

julia> Plasm.view( Plasm.numbering(3)((V,[VV, EV])) )
```

```julia
model = (V,EV)
V,EVs = biconnectedComponent(model)
EW = convert(Lar.Cells, cat(EVs))
Plasm.view( Plasm.numbering(3)((V,[VV,EW])) )
```

Another example:

```julia
V = [0. 1  0  0 -1;
	 0  0  1 -1  0]
EV = [[1,2],[2,3],[1,4],[4,5],[5,1],[1,3]  ,[2,4]]
model = (V,EV)
Plasm.view( Plasm.numbering(1)((V,[[[k] for k=1:size(V,2)],EV])) )
V,EVs = Lar.biconnectedComponent(model)

cscEV = Lar.characteristicMatrix(EV)
biconcomp = Lar.Arrangement.biconnected_components(cscEV)

Matrix(Lar.characteristicMatrix(EV))

V,EVs = biconnectedComponent(model)
EW = convert(Lar.Cells, cat(EVs))
Plasm.view( Plasm.numbering(.3)((V,[VV,EW])) )
```

"""
function biconnectedComponent(model)
    W,EV = model
    V = collect(1:size(W,2))
    count = 0
    stack,out = [],[]
    visited = [false for v in V]
    parent = Union{Int, Array{Any,1}}[[] for v in V]
    d = Any[0 for v in V]
    low = Any[0 for v in V]
    VV = Lar.verts2verts(EV)
    out = Any[]
    for u in V
        if ! visited[u]
            out = DFV_visit( VV,out,count,visited,parent,d,low,stack, u )
        end
    end
    out = [component for component in out if length(component) >= 1]
    EVs = [[map(sort∘collect,edges) for edges in cat(comp) if length(edges)>1]
    		for comp in out]
    EVs = cat(filter(x->!isempty(x), EVs))
    EVs = sort(EVs, lt=(x,y)->length(x)>length(y))
    return W, EVs
end

"""
	depth_first_search(EV::Lar.Cells [, start::Int=1])::Tuple{Lar.Cells, Lar.Cells}

*Depth-First-Search algorithm* (Tarjan, 1972) to compute a *spanning tree* of an undirected graph. The graph is entered as an array of edges, given as an array `EV` of 1-cells. Complexity linear in the number of nodes and edges: O(V, E).

Return a partition of `EV` into two arrays `spanningtree` and `fronds` which define a ``palm tree``.  The algorithm may `start` from any vertex index.

# Example

```
julia> EV = [[1,2],[2,3],[3,4],[4,5],[5,6],[1,6],
       [6,7],[4,7],[2,7],[5,8],[1,8],[3,8]];

julia> depth_first_search(EV,8)
(Array{Int64,1}[[8, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]], Array{Int64,1}[[6, 1], [7, 2], [7, 4], [5, 8], [3, 8]])

julia> depth_first_search(EV)
(Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [5, 8]], Array{Int64,1}[[6, 1], [7, 2], [7, 4], [8, 1], [8, 3]])
```
"""
function depth_first_search(EV::Lar.Cells, start::Int=1)::
		Tuple{Lar.Cells, Lar.Cells}
	# initialize the data structures and start DFS
    VV = Lar.verts2verts(EV)
    spanningtree = Array{Int,1}[];
    fronds = Array{Int,1}[]
	number = zeros(Int, length(VV))

	function DFS(v::Int, u::Int)
	# vertex u is the father of
	# vertex v in the spanning tree being constructed
		number[v] = i; i += 1
		for w in VV[v]
			if number[w] == 0
				# w is not yet numbered
				push!(spanningtree, [v,w])
				DFS(w, v)
			elseif (number[w] < number[v]) && (w != u)
				push!(fronds, [v,w])
			end
		end
	end

	# initialization of DFS algorithm
	i = 1
	DFS(start, 1)
	return spanningtree::Lar.Cells, fronds::Lar.Cells
end

"""
	edge_biconnect(EV::Lar.Cells)

Compute maximal 2-edge-connected components of an undirected graph.

# Examples

```julia
V = hcat([[1.,2],[1,3],[0,3],[1,4],[2,3],[2,2],[1,1],[1,0],[0,1]]...)
EV = [[3,4],[2,4],[2,3],[2,5],[1,2],[5,6],[1,6],[1,7],[7,9],[8,9],[7,8]]
using Plasm
Plasm.view( Plasm.numbering(1.3)((V,[[[k] for k=1:size(V,2)], EV])) )

julia> edge_biconnect(EV::Lar.Cells)
4-element Array{Array{Array{Int64,1},1},1}:
 [[4, 2], [3, 4], [2, 3]]
 [[6, 1], [5, 6], [2, 5], [1, 2]]
 [[9, 7], [8, 9], [7, 8]]
 [[1, 7]]
```

```julia
julia> V = hcat([[0,4],[4,4],[4,0],[2,0],[0,0],[2,2],[1,2],[3,2],[1,3],[3,3]]...)
julia> EV = [[1,2],[9,10],[1,5],[7,9],[8,10],[2,3],[6,7],[6,8],[4,6],[4,5],[3,4]]
julia> Plasm.view( Plasm.numbering(1.3)((V,[[[k] for k=1:size(V,2)], EV])) )

julia> EVs = edge_biconnect(EV::Lar.Cells)
3-element Array{Array{Array{Int64,1},1},1}:
 [[8, 6], [10, 8], [9, 10], [7, 9], [6, 7]]
 [[4, 6]]
 [[5, 1], [4, 5], [3, 4], [2, 3], [1, 2]]

julia> map(Lar.depth_first_search, EVs)
```
"""
function edge_biconnect(EV::Lar.Cells)
    VV = Lar.verts2verts(EV)
	number = zeros(Int, length(VV))
	lowpt = zeros(Int, length(VV))
	components = Array{Array{Int,1},1}[]

	function biconnect(v, u)
		number[v] = i; i +=  1
		lowpt[v] = number[v]
		for w in VV[v]
			if number[w] == 0  # w is not yet numbered
				push!(stack, [v, w]) # add [v, w] to edges stack

				biconnect(w, v)

				lowpt[v] = min(lowpt[v], lowpt[w])
				if lowpt[w] ≥ number[v]
					# start new biconnected component
					component = Array{Int,1}[]
					while number[stack[end][1]] ≥ number[w]
						u1,u2 = pop!(stack)
						push!(component, [u1,u2])
					end
					v,w = pop!(stack)
					push!(component, [v,w])
					push!(components, component)
				end
			else
				if number[w] < number[v] && w ≠ u
					push!(stack, [v,w]) # add (v, w) to edge stack;
					lowpt[v] = min(lowpt[v], number[w]);
				end
			end
		end
	end

	i = 1
	stack = Array{Int64,1}[]
	for w=1:length(VV)
		if number[w] == 0 # w is not yet numbered
			biconnect(w, 0)
		end
	end
	return components
end


text:
# Digital signal plots

import .LibCImPlot: PlotDigital

function PlotDigital(x::AbstractArray{T}, y::AbstractArray{T}; label_id::String = "",
                     count::Integer = min(length(x),length(y)), offset::Integer = 0,
                     stride::Integer = 1) where {T<:ImPlotData}

    LibCImPlot.PlotDigital(label_id, x, y, count, offset, stride * sizeof(T))
end

function PlotDigital(x::AbstractArray{T1}, y::AbstractArray{T2}; kwargs...) where {T1<:Real,T2<:Real}
    PlotDigital(promote(x,y)...; kwargs...)
end


text:
export  formats,
        dataless_support,
        mseed_support,
        parsemseed!,
        read_dataless,
        read_mseed_file!,
        read_seed_resp!,
        resp_wont_read,
        scan_seed,
        seed_support


text:
isdefined(Main, :__precompile__) && __precompile__(true)

module Debug

include("AST.jl")
include("Meta.jl")
include("Analysis.jl")
include("Runtime.jl")
include("Graft.jl")
include("Eval.jl")
include("Flow.jl")
include("UI.jl")
using Debug.AST, Debug.Meta, Debug.Analysis, Debug.Graft, Debug.Eval
using Debug.Flow, Debug.UI

# It seems that @instrument has to be exported in order not to be deleted
export @debug, @bp, @instrument, @localscope, debug_eval, @debug_analyze, @notrap

is_trap(::Event)    = false
is_trap(::LocNode)  = false
is_trap(node::Node) = isblocknode(parentof(node))

macro debug(ex)
    code_debug(UI.instrument(ex))
end
never_trap(ex) = false
macro debug_analyze(ex)
    code_debug(instrument(never_trap, nothing, ex))
end
macro instrument(trap_ex, ex)
    @gensym trap_var
    code_debug(quote
        const $trap_var = $trap_ex
        $(instrument(is_trap, trap_var, ex))
    end)
end

function code_debug(ex)
    globalvar = gensym("globalvar")
    esc(quote
        $globalvar = nothing
        if !isdefined($(quot(globalvar)))
            error("@debug: must be applied in global (i.e. module) scope!")
        end
        $ex
    end)
end

end # module


text:
"
    _ftIdx(ft, grp)

Finds the index of grp (string) in ft.labels which corresponds to ft.mat.
For Index i for the grp  returned by this function ft[i, :, :] returns the
2D array [[TP, FP], [FN, TN]] for that group.
"
function _ftIdx(ft::FairTensor, grp)
    idx = findfirst(x->x==string.(grp), ft.labels)
    if idx==nothing throw(ArgumentError("$grp not found in the fairness tensor")) end
    return idx
end

# Helper function for calculating TruePositive, FalseNegative, etc.
# If grp is :, it calculates combined value for all groups, else the specified group
_calcmetric(ft::FairTensor, grp, inds...) = typeof(grp)==Colon ? sum(ft[:, inds...]) : sum(ft[_ftIdx(ft, grp), inds...])

(::TruePositive)(ft::FairTensor; grp=:) = _calcmetric(ft, grp, 1, 1)
(::FalsePositive)(ft::FairTensor; grp=:) = _calcmetric(ft, grp, 1, 2)
(::FalseNegative)(ft::FairTensor; grp=:) = _calcmetric(ft, grp, 2, 1)
(::TrueNegative)(ft::FairTensor; grp=:) = _calcmetric(ft, grp, 2, 2)

# The functions true_positive, false_negative, etc are instances of TruePositive, FalseNegative, etc.
# So on using true_positive(fair_tensor) will use the above defined functions
(::TPR)(ft::FairTensor; grp=:) = 1/(1+false_negative(ft; grp=grp)/(1e-15 + true_positive(ft; grp=grp)))
(::TNR)(ft::FairTensor; grp=:) = 1/(1+false_positive(ft; grp=grp)/(1e-15 + true_negative(ft; grp=grp)))
(::FPR)(ft::FairTensor; grp=:) = 1-true_negative_rate(ft; grp=grp)
(::FNR)(ft::FairTensor; grp=:) = 1-true_positive_rate(ft; grp=grp)

(::FDR)(ft::FairTensor; grp=:) = 1/(1+false_positive(ft; grp=grp)/(1e-15 + true_positive(ft; grp=grp)))
(::Precision)(ft::FairTensor; grp=:) = 1 - false_discovery_rate(ft; grp=grp)
(::NPV)(ft::FairTensor; grp=:) =  1/(1+false_negative(ft; grp=grp)/(1e-15 + true_negative(ft; grp=grp)))

(::MLJBase.Accuracy)(ft::FairTensor; grp=:) = (true_positive(ft; grp = grp) + true_negative(ft; grp = grp))/(
        true_positive(ft; grp = grp) + true_negative(ft; grp = grp) +
        false_positive(ft; grp = grp) + false_negative(ft; grp = grp))