text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
using Printf
varname_generic(m::Model, col::Integer) = "VAR$(col)"

function varname_given(m::Model, col::Integer)
    # TODO: deal with non-ascii characters?
    name = getname(m, col)
    for (pat, sub) in [("[", "_"), ("]", ""), (",", "_")]
        name = replace(name, pat => sub)
    end
    name
end

function writeMOP(m, fname::AbstractString, genericnames=false)
    varname = genericnames ? varname_generic : varname_given

    f = open(fname, "w")

    write(f,"NAME   vOptModel\n")

    numRows = length(m.linconstr)

    md = getvOptData(m)
    objlincoef = [prepAffObjective(m, f) for f in md.objs]

    write(f,"OBJSENSE\n")
    if md.objSenses[1] == :Min
        write(f," MIN\n")
    else
        write(f," MAX\n")
    end

    for i = 2:length(md.objs)
        if md.objSenses[i] != md.objSenses[1]
            println("Warning, MOP does not support different optimisation senses. Flipping objective coefficients.")
            objlincoef[i] = -objlincoef[i]
        end
    end

    # Objective and constraint names
    write(f,"ROWS\n")

    for i = 1:length(md.objs)
        @printf(f," N  OBJ%d\n", i)
    end

    hasrange = false
    for c in 1:numRows
        rowsense = JuMP.sense(m.linconstr[c])
        if rowsense == :(<=)
            senseChar = 'L'
        elseif rowsense == :(==)
            senseChar = 'E'
        elseif rowsense == :(>=)
            senseChar = 'G'
        else
            hasrange = true
            senseChar = 'E'
        end
        @printf(f," %c  CON%d\n",senseChar,c)
    end

    A = JuMP.prepConstrMatrix(m)
    colptr = A.colptr
    rowval = A.rowval
    nzval = A.nzval

    # Output each column
    inintegergroup = false
    write(f,"COLUMNS\n")
    for col in 1:m.numCols
        t = m.colCat[col]
        (t == :SemiCont || t == :SemiInt) && error("The MPS file writer does not currently support semicontinuous or semi-integer variables")
        if (t == :Bin || t == :Int) && !inintegergroup
            @printf(f,"    MARKER    'MARKER'                 'INTORG'\n")
            inintegergroup = true
        elseif (t == :Cont || t == :Fixed) && inintegergroup
            @printf(f,"    MARKER    'MARKER'                 'INTEND'\n")
            inintegergroup = false
        end
        for ind in colptr[col]:(colptr[col+1]-1)
            @printf(f,"    %s  CON%d  ",varname(m,col),rowval[ind])
            (Base.Grisu).print_shortest(f,nzval[ind])
            println(f)
        end
        for obj = 1:length(md.objs)
            @printf(f,"    %s  OBJ%d  ",varname(m,col), obj)
            (Base.Grisu).print_shortest(f,objlincoef[obj][col])
            println(f)
        end
    end
    if inintegergroup
        @printf(f,"    MARKER    'MARKER'                 'INTEND'\n")
    end

    # RHSs
    write(f,"RHS\n")
    for c in 1:numRows
        rowsense = JuMP.sense(m.linconstr[c])
        if rowsense != :range
            @printf(f,"    RHS    CON%d    ",c)
            (Base.Grisu).print_shortest(f,JuMP.rhs(m.linconstr[c]))
        else
            @printf(f,"    RHS    CON%d    ",c)
            (Base.Grisu).print_shortest(f,m.linconstr[c].lb)
        end
        println(f)
    end

    # RANGES
    if hasrange
        write(f,"RANGES\n")
        for c in 1:numRows
            rowsense = JuMP.sense(m.linconstr[c])
            if rowsense == :range
                @printf(f,"    rhs    CON%d    ",c)
                (Base.Grisu).print_shortest(f,m.linconstr[c].ub-m.linconstr[c].lb)
                println(f)
            end
        end
    end

    # BOUNDS
    write(f,"BOUNDS\n")
    for col in 1:m.numCols
        if m.colLower[col] == 0
            if m.colUpper[col] != Inf
                # Default lower 0, and an upper
                @printf(f,"  UP BOUND %s ", varname(m,col))
                (Base.Grisu).print_shortest(f, m.colUpper[col])
                println(f)
            else
                # Default bounds. Explicitly state for solvers like Gurobi. See issue #792
                @printf(f,"  PL BOUND %s", varname(m,col))
                println(f)
            end
        elseif m.colLower[col] == -Inf && m.colUpper[col] == +Inf
            # Free
            @printf(f, "  FR BOUND %s\n", varname(m,col))
        elseif m.colLower[col] != -Inf && m.colUpper[col] == +Inf
            # No upper, but a lower
            @printf(f, "  PL BOUND %s\n  LO BOUND %s ",varname(m,col),varname(m,col))
            (Base.Grisu).print_shortest(f,m.colLower[col])
            println(f)
        elseif m.colLower[col] == -Inf && m.colUpper[col] != +Inf
            # No lower, but a upper
            @printf(f,"  MI BOUND %s\n  UP BOUND %s ",varname(m,col),varname(m,col))
            (Base.Grisu).print_shortest(f,m.colUpper[col])
            println(f)
        else
            # Lower and upper
            @printf(f, "  LO BOUND %s ",varname(m,col))
            (Base.Grisu).print_shortest(f,m.colLower[col])
            println(f)
            @printf(f, "  UP BOUND %s ",varname(m,col))
            (Base.Grisu).print_shortest(f,m.colUpper[col])
            println(f)
        end
    end

    write(f,"ENDATA\n")
    close(f)
    nothing
end

nextline(f) = split(chomp(readline(f)), ' ', keepempty=false)
function parseMOP(fname::AbstractString; solver=JuMP.UnsetSolver())

    m = vModel(solver = solver)

    open(fname) do f
        ln = nextline(f)
        while ln[1] != "NAME"
            ln = nextline(f)
        end

        objSense = :Min
        ln = nextline(f)
        if ln[1] == "OBJSENSE"
            ln = nextline(f)
            if ln[1] == "MAXIMIZE" || ln[1] == "MAX"
                objSense = :Max
            end
            ln = nextline(f)
        end

        if ln[1] == "ROWS"
            ln = nextline(f)
            DicoExpr = Dict{AbstractString,Tuple{JuMP.AffExpr, Symbol}}()
            objOrder = Vector{AbstractString}()
            while ln[1] != "COLUMNS"
                push!(DicoExpr, ln[2] => ( zero(JuMP.AffExpr), Symbol(ln[1]) ) )
                ln[1] == "N" && push!(objOrder, ln[2])
                ln = nextline(f)
            end
        end

        if ln[1] == "COLUMNS"
            ln = nextline(f)
            DicoVar = Dict{AbstractString,JuMP.Variable}()
            isInt = false
            while ln[1] != "RHS"

                if ln[3] == "'INTORG'"
                    isInt = true
                elseif ln[3] == "'INTEND'"
                    isInt = false
                else
                    if !(ln[1] in keys(DicoVar))

                        _cat = isInt ? (:Int) : (:Cont)
                        v = @variable(m, category=_cat, basename=ln[1], lowerbound=0.0)

                        # v = @variable(m, basename=ln[1], lowerbound=0.0)
                        # if isInt
                        #     m.colCat[v.col] = :Int
                        # end

                        push!(DicoVar, ln[1] => v)
                    end

                    expr, expr_type = DicoExpr[ln[2]]
                    var = DicoVar[ln[1]]
                    coeff = parse(Float64, ln[3])

                    push!(expr, coeff, var)

                end
                ln = nextline(f)
            end
        end

        DicoCstr = Dict{AbstractString, Tuple{JuMP.LinearConstraint, Symbol}}()
        DicoObj  = Dict{AbstractString, JuMP.QuadExpr}()
        for (k,v) in DicoExpr
            expr, expr_type = v
            if expr_type == :N
                push!(DicoObj, k => convert(QuadExpr, expr))
            else
                push!(DicoCstr, k => (LinearConstraint(expr, -Inf, Inf), expr_type) )
            end
        end

        md = m.ext[:vOpt]
        for key in objOrder
            obj = DicoObj[key]
            push!(md.objs, obj)
            push!(md.objSenses, objSense)
        end

        if ln[1] == "RHS"
            ln = nextline(f)
            while length(ln) > 1

                constraint, constraint_type = DicoCstr[ln[2]]
                if constraint_type == :E
                    constraint.lb = parse(Float64, ln[3])
                    constraint.ub = parse(Float64, ln[3])
                elseif constraint_type == :G
                    constraint.lb = parse(Float64, ln[3])
                elseif constraint_type == :L
                    constraint.ub = parse(Float64, ln[3])
                end
                ln = nextline(f)
            end
        end

        for (k,v) in DicoCstr
            cstr, cstr_type = v
            JuMP.addconstraint(m, cstr)
        end

        if ln[1] == "RANGES"
            ln = nextline(f)
            while length(ln) > 1

                constraint, constraint_type = DicoCstr[ln[2]]
                r = parse(Float64, ln[3])
                if constraint_type == :G
                    constraint.ub = constraint.lb + abs(r)
                elseif constraint_type == :L
                    constraint.lb = constraint.ub - abs(r)
                else
                    if r >= 0
                        constraint.ub = constraint.lb + r
                    else
                        constraint.lb = constraint.ub + r
                    end
                end
                ln = nextline(f)
            end
        end

        fill!(m.colLower, 0.0)

        if ln[1] == "BOUNDS"
            ln = nextline(f)
            while length(ln) > 1

                boundType, varName = ln[1], ln[3]
                var = DicoVar[varName]
                if boundType == "LO"
                    val = parse(Float64, ln[4])
                    JuMP.setlowerbound(var, val)
                elseif boundType == "UP"
                    val = parse(Float64, ln[4])
                    JuMP.setupperbound(var, val)
                elseif boundType == "FX"
                    val = parse(Float64, ln[4])
                    JuMP.setupperbound(var, val)
                    JuMP.setlowerbound(var, val)
                elseif boundType == "FR"
                    JuMP.setupperbound(var, Inf)
                    JuMP.setlowerbound(var, -Inf)
                elseif boundType == "MI"
                    JuMP.setlowerbound(var, -Inf)
                elseif boundType == "PL"
                    JuMP.setupperbound(var, Inf)
                elseif boundType == "BV"
                    setcategory(var, :Bin)
                elseif boundType == "LI"
                    val = parse(Float64, ln[4])
                    setlowerbound(var, val)
                elseif boundType == "UI"
                    val = parse(Float64, ln[4])
                    setupperbound(var, val)
                elseif boundType == "SC"
                    error("Semi-continuous variables are not supported by vOptGeneric")
                end
                ln = nextline(f)
            end
        end
    end

    return m
end


text:
function encode_attachment(filename::String)
    io = IOBuffer()
    iob64_encode = Base64EncodePipe(io)
    open(filename, "r") do f
        write(iob64_encode, f)
    end
    close(iob64_encode)

    filename_ext = split(filename, '.')[end]

    if haskey(mime_types, filename_ext)
        content_type = mime_types[filename_ext]
    else
        content_type = "application/octet-stream"
    end

    if haskey(mime_types, filename_ext) && startswith(mime_types[filename_ext], "image")
        content_disposition = "inline"
    else
        content_disposition = "attachment"
    end

    encoded_str =
        "Content-Disposition: $content_disposition;\r\n" *
        "    filename=\"$(basename(filename))\"\r\n" *
        "Content-Type: $content_type;\r\n" *
        "    name=\"$(basename(filename))\"\r\n" *
        "Content-ID: <$(basename(filename))>\r\n" *
        "Content-Transfer-Encoding: base64\r\n" *
        "\r\n" *
        "$(String(take!(io)))\r\n"
    return encoded_str
end

# See https://www.w3.org/Protocols/rfc1341/7_1_Text.html about charset
function get_mime_msg(message::String, ::Val{:plain}, charset::String = "UTF-8")
    msg =
        "Content-Type: text/plain; charset=\"$charset\"\r\n" *
        "Content-Transfer-Encoding: quoted-printable\r\n\r\n" *
        "$message\r\n"
    return msg
end

get_mime_msg(message::String, ::Val{:utf8}) =
    get_mime_msg(message, Val(:plain), "UTF-8")

get_mime_msg(message::String, ::Val{:usascii}) =
    get_mime_msg(message, Val(:plain), "US-ASCII")

get_mime_msg(message::String) = get_mime_msg(message, Val(:utf8))

function get_mime_msg(message::String, ::Val{:html})
    msg =
        "Content-Type: text/html;\r\n" *
        "Content-Transfer-Encoding: 7bit;\r\n\r\n" *
        "\r\n" *
        "\r\n" *
        message *
        "\r\n"
    return msg
end

get_mime_msg(message::HTML{String}) = get_mime_msg(message.content, Val(:html))

get_mime_msg(message::Markdown.MD) = get_mime_msg(Markdown.html(message), Val(:html))

#Provide the message body as RFC5322 within an IO

function get_body(
        to::Vector{String},
        from::String,
        subject::String,
        msg::String;
        cc::Vector{String} = String[],
        replyto::String = "",
        attachments::Vector{String} = String[]
    )

    boundary = "Julia_SMTPClient-" * join(rand(collect(vcat('0':'9','A':'Z','a':'z')), 40))

    tz = mapreduce(
        x -> string(x, pad=2), *,
        divrem( div( ( now() - now(Dates.UTC) ).value, 60000 ), 60 )
    )
    date = join([Dates.format(now(), "e, d u yyyy HH:MM:SS", locale="english"), tz], " ")

    contents =
        "From: $from\r\n" *
        "Date: $date\r\n" *
        "Subject: $subject\r\n" *
        ifelse(length(cc) > 0, "Cc: $(join(cc, ", "))\r\n", "") *
        ifelse(length(replyto) > 0, "Reply-To: $replyto\r\n", "") *
        "To: $(join(to, ", "))\r\n"

    if length(attachments) == 0
        contents *=
            "MIME-Version: 1.0\r\n" *
            "$msg\r\n\r\n"
    else
        contents *=
            "Content-Type: multipart/mixed; boundary=\"$boundary\"\r\n" *
            "MIME-Version: 1.0\r\n" *
            "\r\n" *
            "This is a message with multiple parts in MIME format.\r\n" *
            "--$boundary\r\n" *
            msg *
            "\r\n--$boundary\r\n" *
            join(encode_attachment.(attachments), "\r\n--$boundary\r\n") *
            "\r\n--$boundary--\r\n"
    end
    body = IOBuffer(contents)
    return body
end

ext: jl


text:
export UserId, ChannelId, TeamId, EventTimestamp
export SlackName

@newimmutable UserId <: UTF8String
@stringinterface UserId

@newimmutable ChannelId <: UTF8String
@stringinterface ChannelId

@newimmutable TeamId <: UTF8String
@stringinterface TeamId

@newimmutable GroupId <: UTF8String
@stringinterface GroupId

@newimmutable ImId <: UTF8String
@stringinterface ImId

@newimmutable BotId <: UTF8String
@stringinterface BotId

@newimmutable SlackName <: UTF8String
@stringinterface SlackName

@newimmutable MpimName <: UTF8String
@stringinterface MpimName

@newimmutable BotName <: UTF8String
@stringinterface BotName

@newimmutable Timestamp <: UInt64

@newimmutable Token <: UTF8String
@stringinterface Token

@newimmutable EventTimestamp <: UTF8String
@stringinterface EventTimestamp


text:
COMPANY_OVERVIEW_KEYS =  ["SharesOutstanding",
"ExDividendDate",
"52WeekLow",
"ReturnOnEquityTTM",
"LatestQuarter",
"200DayMovingAverage",
"EVToEBITDA",
"RevenuePerShareTTM",
"Beta",
"Sector",
"ForwardAnnualDividendYield",
"Exchange",
"PercentInsiders",
"QuarterlyEarningsGrowthYOY",
"Currency",
"EBITDA",
"ShortRatio",
"DividendYield",
"AnalystTargetPrice",
"DilutedEPSTTM",
"BookValue",
"LastSplitDate",
"SharesFloat",
"PriceToSalesRatioTTM",
"FullTimeEmployees",
"ShortPercentOutstanding",
"52WeekHigh",
"ReturnOnAssetsTTM",
"PayoutRatio",
"PriceToBookRatio",
"Symbol",
"QuarterlyRevenueGrowthYOY",
"DividendDate",
"ProfitMargin",
"50DayMovingAverage",
"LastSplitFactor",
"RevenueTTM",
"PEGRatio",
"OperatingMarginTTM",
"Industry",
"ForwardPE",
"EVToRevenue",
"ForwardAnnualDividendRate",
"ShortPercentFloat",
"SharesShortPriorMonth",
"TrailingPE",
"SharesShort",
"Country",
"Address",
"EPS",
"GrossProfitTTM",
"Name",
"MarketCapitalization",
"PERatio",
"DividendPerShare",
"Description",
"FiscalYearEnd",
"AssetType",
"PercentInstitutions"]

COMMON_KEYS = ["reportedCurrency", "netIncome"]

COMMON = [[("quarterlyReports", "income_statement", x, "fiscalDateEnding") for x in COMMON_KEYS];
            [("annualReports", "income_statement", x, "fiscalDateEnding") for x in COMMON_KEYS]]

INCOME_STATEMENT_KEYS = [
"incomeTaxExpense",
"otherNonOperatingIncome",
"minorityInterest",
"discontinuedOperations",
"incomeBeforeTax",
"totalOtherIncomeExpense",
"interestIncome",
"researchAndDevelopment",
"grossProfit",
"totalRevenue",
"otherOperatingExpense",
"taxProvision",
"extraordinaryItems",
"ebit",
"otherItems",
"netIncomeApplicableToCommonShares",
"totalOperatingExpense",
"costOfRevenue",
"interestExpense",
"sellingGeneralAdministrative",
"operatingIncome",
"netIncomeFromContinuingOperations",
"netInterestIncome",
"effectOfAccountingCharges",
"nonRecurring",
"preferredStockAndOtherAdjustments"]

INCOME_STATEMENT = [[("quarterlyReports", "income_statement", x, "fiscalDateEnding") for x in INCOME_STATEMENT_KEYS];
                    [("annualReports", "income_statement", x, "fiscalDateEnding") for x in INCOME_STATEMENT_KEYS]]

BALANCE_SHEET_KEYS = ["totalPermanentEquity",
"warrants",
"negativeGoodwill",
"preferredStockTotalEquity",
"accumulatedAmortization",
"inventory",
"additionalPaidInCapital",
"commonStockTotalEquity",
"longTermInvestments",
"netTangibleAssets",
"cashAndShortTermInvestments",
"longTermDebt",
"otherShareholderEquity",
"totalCurrentAssets",
"treasuryStock",
"otherAssets",
"capitalSurplus",
"totalNonCurrentAssets",
"accountsPayable",
"totalNonCurrentLiabilities",
"otherLiabilities",
"totalShareholderEquity",
"liabilitiesAndShareholderEquity",
"otherCurrentAssets",
"totalCurrentLiabilities",
"otherNonCurrrentAssets",
"shortTermDebt",
"commonStockSharesOutstanding",
"capitalLeaseObligations",
"netReceivables",
"retainedEarningsTotalEquity",
"earningAssets",
"totalAssets",
"commonStock",
"cash",
"deferredLongTermLiabilities",
"totalLongTermDebt",
"retainedEarnings",
"shortTermInvestments",
"propertyPlantEquipment",
"goodwill",
"preferredStockRedeemable",
"totalLiabilities",
"otherNonCurrentLiabilities",
"currentLongTermDebt",
"intangibleAssets",
"accumulatedDepreciation",
"otherCurrentLiabilities",
"deferredLongTermAssetCharges"]

BALANCE_SHEET = [[("quarterlyReports", "balance_sheet", x, "fiscalDateEnding") for x in BALANCE_SHEET_KEYS];
                [("annualReports", "balance_sheet", x, "fiscalDateEnding") for x in BALANCE_SHEET_KEYS]]

 CASH_FLOW_KEYS = ["cashflowFromInvestment",
 "changeInInventory",
 "changeInAccountReceivables",
 "changeInCashAndCashEquivalents",
 "otherOperatingCashflow",
 "dividendPayout",
 "changeInReceivables",
 "capitalExpenditures",
 "changeInExchangeRate",
 "operatingCashflow",
 "cashflowFromFinancing",
 "changeInLiabilities",
 "stockSaleAndPurchase",
 "otherCashflowFromFinancing",
 "changeInOperatingActivities",
 "depreciation",
 "changeInCash",
 "netBorrowings",
 "investments",
 "changeInNetIncome",
 "otherCashflowFromInvestment"]

 CASH_FLOW = [[("quarterlyReports", "cash_flow", x, "fiscalDateEnding") for x in CASH_FLOW_KEYS];
              [("annualReports", "cash_flow", x, "fiscalDateEnding") for x in CASH_FLOW_KEYS]]

 EARNINGS_KEYS_Q = ["reportedDate",
 "estimatedEPS",
 "surprise",
 "surprisePercentage",
 "reportedEPS"]

 EARNINGS_KEYS_A = ["reportedEPS"]

 EARNINGS = [[("quarterlyEarnings", "earnings", x, "reportedDate") for x in EARNINGS_KEYS_Q];
             [("annualEarnings", "earnings", x, "fiscalDateEnding") for x in EARNINGS_KEYS_A]]

 FUNDAMENTAL_VALUES = vcat(INCOME_STATEMENT, BALANCE_SHEET, CASH_FLOW, EARNINGS)


text:
let
    verbose = true
    upstream_repo = ""
    local_root = ""

    try

        println("\n", "-"^60)
        println("Testing git utils")
        println("-"^60, "\n")

        # upstream
        branch_name = "test_branch"
        url, upstream_repo = GitLinks.create_local_upstream(tempname(); verbose, branch_name)
        @assert isdir(upstream_repo)

        @test !isempty(GitLinks._remote_HEAD_hash(url))
        @test !isempty(GitLinks._HEAD_hash(upstream_repo))
        @test GitLinks._remote_HEAD_hash(url) == GitLinks._HEAD_hash(upstream_repo)

        @test GitLinks._check_gitdir(upstream_repo)
        @test !GitLinks._check_gitdir(tempname())
        @test GitLinks._curr_branch(upstream_repo) == branch_name

        # local repo
        local_root = tempname()
        GitLinks._rm(local_root)

        gl = GitLinks.GitLink(local_root, url)
        local_repo = GitLinks.repo_dir(gl)
        local_repo_git = joinpath(local_repo, ".git")

        # hard pull
        @test !isdir(local_root)
        pull_ok = GitLinks.hard_pull(gl; verbose, clearwd = true) # clone
        @test pull_ok
        @test isdir(local_root)
        @test isdir(local_repo_git)

        pull_ok = GitLinks.hard_pull(gl; verbose, clearwd = true) # pull
        @test pull_ok
        @test isdir(local_root)
        @test isdir(local_repo_git)

        GitLinks._rm(joinpath(local_repo, ".git")) # break repo
        @test !GitLinks._check_gitdir(local_root)
        pull_ok = GitLinks.hard_pull(gl; verbose, clearwd = true) # must recover
        @test pull_ok
        @test isdir(local_root)
        @test isdir(local_repo_git)

        # soft push
        upsize0 = GitLinks._foldersize(upstream_repo)
        for it in 1:10
            # create 'big' files
            dummy = joinpath(local_repo, "dummy$it")
            write(dummy, GitLinks.rand_str(1000))
            push_ok = GitLinks.soft_push(gl; verbose)
            @test push_ok
        end
        upsize1 = GitLinks._foldersize(upstream_repo)
        @test upsize0 < upsize1

        println("\n", "-"^60)
        println("Before nuking", "\n")
        run(`git -C $(upstream_repo) --no-pager log -l10 --pretty=oneline`)

        # TODO: test that nuking actually reduce size
        # nuke
        # upsize2 = GitLinks._foldersize(upstream_repo)
        # upsize3 = GitLinks._foldersize(local_root)

        nuke_ok = GitLinks.nuke_remote(gl; verbose)
        @test nuke_ok

        GitLinks._rm(local_root)
        pull_ok = GitLinks.hard_pull(gl; verbose, clearwd = true) # clone again
        @test pull_ok
        @test isdir(local_root)
        @test isdir(local_repo_git)

        # upsize4 = GitLinks._foldersize(local_root)
        # upsize5 = GitLinks._foldersize(upstream_repo)
        # GitLinks._run("git -C $(upstream_repo) gc")
        # upsize6 = GitLinks._foldersize(upstream_repo)
        # @test upsize2 < upsize1
        # @show upsize0 upsize1 upsize2 upsize3 upsize4 upsize5 upsize6
        # @show upsize2 upsize5
        # @show upsize3 upsize4
        # @test upsize2 < upsize1

        println("\n", "-"^60)
        println("After nuking", "\n")
        run(`git -C $(upstream_repo) --no-pager log -l10 --pretty=oneline`)
        println()

    finally
        # clear
        GitLinks._rm(local_root)
        GitLinks._rm(upstream_repo)
    end
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
@inline function rot16(x::VECTOR)::VECTOR
    x = reinterpret(Vec{8, UInt16}, x)
    r = shufflevector(x, Val((1, 0, 3, 2, 5, 4, 7, 6)))
    return reinterpret(VECTOR, r)
end

@inline function rot12(x::VECTOR)::VECTOR
    return (x >> 12) ⊻ (x << 20)
end

@inline function rot8(x::VECTOR)::VECTOR
    x = reinterpret(Vec{16, UInt8}, x)
    r = shufflevector(x, Val((1, 2, 3, 0, 5, 6, 7, 4, 9, 10, 11, 8, 13, 14, 15, 12)))
    return reinterpret(VECTOR, r)
end

@inline function rot7(x::VECTOR)::VECTOR
    return (x >> 7) ⊻ (x << 25)
end


text:
using kRPC
using CSV
using LinearAlgebra
using DataFrames
using Plots
pyplot()

conn = kRPCConnect("Aerodata")
using kRPC.Remote.SpaceCenter
totable = 1
outp_a = 1
outp_b = 1
vessel = ActiveVessel()
flight = Flight(vessel)
orbit = Orbit(vessel)
body = Body(orbit)
body_frame = ReferenceFrame(body)

function compute_reference_frame(body, lat, lon, alt::Float64)
    landing_position = SurfacePosition(body, lat, lon, body_frame)
    q_long = (
        0.0,
        sin(-lon * 0.5 * pi / 180),
        0.0,
        cos(-lon * 0.5 * pi / 180))
    q_lat = (
        0.0,
        0.0,
        sin(lat * 0.5 * pi / 180),
        cos(lat * 0.5 * pi / 180))
    cr1 = CreateRelative(body_frame, position=landing_position, rotation=q_long)
    cr2 = CreateRelative(cr1,rotation=q_lat)
    return CreateRelative(cr2,position=(alt,0.0,0.0))
end

rf = compute_reference_frame(body, -0.09720792138114814, -74.55767985150328, 0.0)
gfl = Flight(vessel, referenceFrame=body_frame)

cforce = Vector{Float64}[]
drg = Vector{Float64}[]
lft = Vector{Float64}[]
aoal = Float64[]
machl = Float64[]
p=plot([],[], xlims = (0,180), ylims = (0,1), legend=false, linetype=:scatter)
iter = [0]
laoa = [-10000.0]
lmach = [-10000.0]

rtb = RotationalSpeed(body)

drag_vis = kRPC.Remote.Drawing.AddLine((0.0,0.0,0.0),(0.0,0.0,0.0),body_frame)
lift_vis = kRPC.Remote.Drawing.AddLine((0.0,0.0,0.0),(0.0,0.0,0.0),body_frame)
aerf_vis = kRPC.Remote.Drawing.AddLine((0.0,0.0,0.0),(0.0,0.0,0.0),body_frame)
kRPC.Remote.Drawing.Color!(lift_vis, (1.0,0.0,0.0))
kRPC.Remote.Drawing.Color!(aerf_vis, (0.0,1.0,0.0))
global lvel = Vector{Float64}[[0.0,0.0,0.0]]
global lut = -10000.0
add_multiple_streams([kRPC.Remote.SpaceCenter.Delayed.AerodynamicForce(gfl),
	kRPC.Remote.SpaceCenter.Delayed.Direction(gfl),
	kRPC.Remote.SpaceCenter.Delayed.Velocity(gfl),
	kRPC.Remote.SpaceCenter.Delayed.Position(vessel, body_frame),
	kRPC.Remote.SpaceCenter.Delayed.AtmosphereDensity(gfl),
	kRPC.Remote.SpaceCenter.Delayed.SpeedOfSound(gfl),
	kRPC.Remote.SpaceCenter.Delayed.UT()], function (aerof,direction,velocity,posn,density,spos,iut)
	sos = Float64(spos)
	ifrc,dir,vel,pos = convert(Vector{Float64}, -[aerof...]), convert(Vector{Float64}, -[direction...]), convert(Vector{Float64}, [velocity...]), convert(Vector{Float64}, [posn...])
	aoa = acosd(dot(dir,vel)/norm(vel))
	mach = norm(vel)/sos
	if aoa == laoa[1] && mach == lmach[1]
		return
	end
	accel = [0.0,0.0,0.0]
	if lut > 0
		dt = iut - lut
		accel = (vel-lvel[1])/dt
	end
	laoa[1] = aoa
	lmach[1] = mach
	global lvel[1] = vel
	global lut = iut

	normvel = vel/norm(vel)
	push!(drg, normvel)
	lftv = cross(cross(dir, normvel),normvel)
	lftv = lftv/norm(lftv)
	push!(cforce, ifrc/density)
	push!(lft, lftv)
	push!(aoal, aoa)
	push!(machl, mach)
	iter[1] = iter[1] + 1
	if iter[1] % 10 == 0
		equivalent_vel = norm(pos) * rsb
		base_vector = Float64[0, 0, -equivalent_vel]
		sfrc = Float64[SimulateAerodynamicForceAt(gfl, body, posn, velocity)...]
		lneup = [
		 kRPC.Remote.Drawing.Delayed.Start!(drag_vis, (pos...,)),
		 kRPC.Remote.Drawing.Delayed.End!(drag_vis, ((pos + sfrc/norm(sfrc)*20)...,)),
		 kRPC.Remote.Drawing.Delayed.Start!(lift_vis, (pos...,)),
		 kRPC.Remote.Drawing.Delayed.End!(lift_vis, ((pos + accel/norm(accel)*20)...,)),
		 kRPC.Remote.Drawing.Delayed.Start!(aerf_vis, (pos...,)),
		 kRPC.Remote.Drawing.Delayed.End!(aerf_vis, ((pos + ifrc/norm(ifrc)*20) ...,))]
		kRPC.SendMessage(conn, lneup)
		println(stdout, norm(sfrc-ifrc)/norm(ifrc))
		push!(p, 1, aoa, mach)
		gui()
	end
end)
p
while true
	a = readline(stdin)
	if a == "e"
		drags = [if !isnan(dot(cforce[i], drg[i])) dot(cforce[i], drg[i]) else 0.0 end for i=1:length(cforce)]
		lifts = [if !isnan(dot(cforce[i], lft[i])) dot(cforce[i], lft[i]) else 0.0 end for i=1:length(cforce)]
		other = [if !isnan(dot(cforce[i], cross(lft[i],drg[i]))) dot(cforce[i], cross(lft[i],drg[i])) else 0.0 end for i=1:length(cforce)]
		op = DataFrame(aoa=aoal,mach=machl,drag=drags,lift=lifts,other=other)
		CSV.write("lift_drag_test.csv", op)
		break
	end
end

max_stars_repo_name: BenChung/SuccessiveConvexification

max_issues_repo_name: BenChung/SuccessiveConvexification

max_forks_repo_name: BenChung/SuccessiveConvexification


text:
module BaysorAnalysis

using DataFrames
using DataFramesMeta
using StatsBase

import Baysor as B

include("data_loading.jl")
include("validation.jl")
include("validation_plots.jl")
include("utils.jl")
include("plots.jl")

end

max_stars_repo_name: hiraksarkar/BaysorAnalysis

max_issues_repo_name: hiraksarkar/BaysorAnalysis

max_forks_repo_name: hiraksarkar/BaysorAnalysis


text:
function loadinputs(folder::String)

    p_D_rowdata=CSV.File(folder *"/demand.csv") |> Tables.matrix
    p_D=p_D_rowdata'

    p_D_MW=p_D*max_demand
    D= maximum(sum(p_D_MW, dims=1))

    Τ=CSV.File(folder *"/carbontax_scenarios.csv") |> Tables.matrix

    wind_rowdata= CSV.File(folder *"/wind_existing.csv") |> Tables.matrix
    wind=wind_rowdata'

    wind_opt_rowdata= CSV.File(folder *"/wind_options.csv") |> Tables.matrix

    wind_opt=wind_opt_rowdata'

    Ns_H= CSV.File(folder *"/weights.csv") |> Tables.matrix

return (p_D,D,Τ,wind,wind_opt,Ns_H)
end

#do read headers
#Write a Csv file
#=
https://stackoverflow.com/questions/54410030/read-csv-into-array
df = DataFrame(rand(2,3))
CSV.write("test.csv", df)
CSV.File("test.csv") |> Tables.matrix
=#


text:
module PlotUtils

include(joinpath(@__DIR__, "imports.jl"))
include(joinpath(@__DIR__, "dden_complete.jl"))

# Load python defs
path_to_plot_defs = @__DIR__
pushfirst!(PyCall.PyVector(PyCall.pyimport("sys")."path"), "$path_to_plot_defs")
plot_yz = PyCall.pyimport("plot_yz")
blue2red = PyCall.pyimport("blue2red")
Population = PyCall.pyimport("Population").Population
pyrange(n) = collect(range(0, stop=n-1))

function quantiles(X, q; dims, drop=false)
  Q = mapslices(x -> quantile(x, q), X, dims=dims)
  out = drop ? dropdims(Q, dims=dims) : Q
  return out
end

skipnan(x) = x[.!isnan.(x)]

extract(s, out) = [o[s] for o in out]

function boxplot(x; showmeans=true, whis=[2.5, 97.5], showfliers=false, kw...)
  plt.boxplot(x, showmeans=showmeans, whis=whis, showfliers=showfliers; kw...)
end

function add_gridlines_Z(Z)
  J, K = size(Z)
  for j in pyrange(J)
    plt.axhline(y=j+.5, color="grey", linewidth=.5)
  end

  for k in pyrange(K)
    plt.axvline(x=k+.5, color="grey", linewidth=.5)
  end
end

axhlines(x; kw...) = for xi in x plt.axhline(xi; kw...) end
cm_greys = plt.cm.get_cmap("Greys", 5)

function plot_missmech(beta, i; xlim=[-10, 5],
                      ygrid=range(xlim[1], stop=xlim[2], length=300))
  p = [CytofRepFAM.Model.prob_miss(y, beta[:, i]) for y in ygrid]
  plt.plot(ygrid, p, lw=3)
  plt.title("Missing Mechanism for Sample $(i)")
  plt.xlabel("Probability of Missing")
  plt.ylabel("Expression Level")
end

# TODO
# Steal things from: Cytof5/sims/repfam_fs/test

function plot_Z(Z; colorbar=true)
  J, K = size(Z)
  p = plt.imshow(Z, aspect="auto", vmin=0, vmax=1, cmap=cm_greys)
  add_gridlines_Z(Z)
  plt.yticks(pyrange(J), pyrange(J) .+ 1, fontsize=rcParams["font.size"])
  plt.xticks(pyrange(K), pyrange(K) .+ 1, fontsize=rcParams["font.size"],
             rotation=90)
  if colorbar
    plt.colorbar()
  end
  return p
end

# plot y/z
function make_yz(y, Zs, Ws, lams, imgdir; vlim,
                 w_thresh=.01, lw=3,
                 Z_true=nothing,
                 markernames=[],
                 #
                 feature_names="cell subpopulations",
                 heatmap_ylabel="cells",
                 #
                 fs_y=rcParams["font.size"],
                 fs_z=rcParams["font.size"],
                 fs_ycbar=rcParams["font.size"],
                 fs_zcbar=rcParams["font.size"],
                 rotation=90, ha="center")
  # Make img dir if needed
  mkpath(imgdir)

  # Make img/txt dir
  txtdir = "$(imgdir)/txt"
  mkpath(txtdir)

  I = length(y)
  population = Population()
  for i in 1:I
    idx_best = estimate_ZWi_index(Zs, Ws, i)

    Zi = Int.(Zs[idx_best])
    Wi = Float64.(Ws[idx_best][i, :])
    lami = Int64.(lams[idx_best][i])

    open("$(txtdir)/Z$(i)_best.txt", "w") do io
      writedlm(io, Zi)
    end

    open("$(txtdir)/W$(i)_best.txt", "w") do io
      writedlm(io, Wi)
    end

    open("$(txtdir)/lam$(i)_best.txt", "w") do io
      writedlm(io, lami)
    end

    yi = Float64.(y[i])

    # plot Yi, lami
    plt.figure(figsize=(6, 6))
    plot_yz.plot_y(yi, Wi, lami, vlim=vlim, cm=blue2red.cm(9), lw=lw,
                   fs_xlab=fs_y, fs_ylab=fs_y, fs_lab=fs_y, fs_cbar=fs_ycbar,
                   ylab=heatmap_ylabel,
                   markernames=markernames, rotation=rotation, ha=ha)
    plt.savefig("$(imgdir)/y$(i).pdf", bbox_inches="tight")
    plt.close()

    # plot Zi, Wi
    plt.figure(figsize=(6, 6))
    plot_yz.plot_Z(Zi, Wi, lami, w_thresh=w_thresh, add_colorbar=false,
                   fs_lab=fs_z, fs_celltypes=fs_z, fs_markers=fs_z,
                   fs_cbar=fs_zcbar, markernames=markernames,
                   ylab="$(feature_names) (abundance)",
                   population=population, rotation=rotation, ha=ha)
    plt.savefig("$(imgdir)/Z$(i).pdf", bbox_inches="tight")
    plt.close()

    # plot Zmean
    plot_yz.plot_Z_only(mean(Zs), fs=fs_z,
                        xlab="cell subpopulations", ylab="markers",
                        rotate_xticks=true)
    plt.colorbar()
    plt.savefig("$(imgdir)/Zmean.pdf", bbox_inches="tight")
    plt.close()
  end

  if Z_true != nothing
    # plot Z true
    plt.figure(figsize=(6, 6))
    plot_yz.plot_Z_only(Z_true, fs=fs_z,
                        xlab="cell subpopulations", ylab="markers",
                        rotate_xticks=false)
    plt.savefig("$(imgdir)/Z_true.pdf", bbox_inches="tight")
    plt.close()

    # plot ZT true
    plt.figure(figsize=(6, 6))
    plot_yz.plot_Z_only(Z_true', fs=fs_z,
                        xlab="markers", ylab="cell subpopulations")
    plt.savefig("$(imgdir)/ZT_true.pdf", bbox_inches="tight")
    plt.close()
  end
end

function plot_loglike(loglike, imgdir; fname="loglike.pdf", title=nothing)
  plt.plot(loglike)
  plt.xlabel("iter")
  plt.ylabel("log-likelihood")
  if title != nothing
    plt.title(title)
  end
  plt.savefig("$(imgdir)/$(fname)", bbox_inches="tight")
  plt.close()
end

function plot_alpha(alphas, imgdir; printmean=true)
  plt.hist(alphas, density=true)
  plt.xlabel("alpha")
  plt.ylabel("density")
  plt.savefig("$(imgdir)/alpha.pdf", bbox_inches="tight")
  plt.close()

  # Print mean alpha
  if printmean
    open("$(imgdir)/txt/alpha_mean.txt", "w") do io
      writedlm(io, mean(alphas))
    end
  end
end

function plot_v(vs, imgdir; printmean=true)
  v_mat = hcat(vs...)

  if printmean
    open("$(imgdir)/txt/v_mean.txt", "w") do io
      writedlm(io, mean(vs))
    end
  end

  boxplot(v_mat')
  plt.xlabel("cell phenotypes")
  plt.ylabel("v")
  plt.savefig("$(imgdir)/v.pdf", bbox_inches="tight")
  plt.close()
end

function plot_p(ps, imgdir; printmean=true)
  ps_mat = Matrix(hcat(ps...)')

  # Plot boxplots
  plt.boxplot(ps_mat)
  plt.xlabel("sample")
  plt.ylabel("p")
  plt.savefig("$(imgdir)/p.pdf", bbox_inches="tight")
  plt.close()

  # Traceplot of p
  pcols = size(ps_mat, 2)
  for i in 1:pcols
    plt.plot(ps_mat[:, i])
    plt.savefig("$(imgdir)/p$(i)_trace.pdf", bbox_inches="tight")
    plt.close()
  end

  if printmean
    open("$(imgdir)/txt/p_mean.txt", "w") do io
      writedlm(io, mean(ps))
    end
  end
end

function plot_mus(deltas, imgdir; printmean=true)
  mus0 = Matrix(hcat([-cumsum(d[0]) for d in deltas]...)')
  mus1 = Matrix(hcat([cumsum(d[1]) for d in deltas]...)')

  if printmean
    open("$(imgdir)/txt/mus0_mean.txt", "w") do io
      writedlm(io, mean(mus0, dims=1))
    end

    open("$(imgdir)/txt/mus1_mean.txt", "w") do io
      writedlm(io, mean(mus1, dims=1))
    end
  end

  mus = [mus0 mus1]
  boxplot(mus)
  plt.axhline(0)
  plt.savefig("$(imgdir)/mus.pdf", bbox_inches="tight")
  plt.close()

  # Traceplot of mus
  plt.plot(mus0)
  plt.plot(mus1)
  plt.savefig("$(imgdir)/mus_trace.pdf", bbox_inches="tight")
  plt.close()
end

function plot_sig2(sig2s, imgdir; sig2_true=nothing, printmean=true)
  sig2s = Matrix(hcat(sig2s...)')

  if printmean
    open("$(imgdir)/txt/sig2_mean.txt", "w") do io
      writedlm(io, mean(sig2s, dims=1))
    end
  end

  boxplot(sig2s)
  if sig2_true != nothing
    axhlines(sig2_true)
  end
  plt.savefig("$(imgdir)/sig2.pdf", bbox_inches="tight")
  plt.close()

  # Traceplot of sig2
  plt.plot(sig2s)
  plt.savefig("$(imgdir)/sig2_trace.pdf", bbox_inches="tight")
  plt.close()
end

function plot_W(Ws, imgdir; W_true=nothing, xlabel="cell subpopulations")
  plot_W(Ws, imgdir=imgdir, W_true=W_true, xlabel=xlabel)
end

function plot_W(Ws; imgdir, W_true=nothing, xlabel="cell subpopulations")
  I, K = size(Ws[1])
  # NOTE: Julia bug. Splat (...) throws errors if number of elems is too large.
  # Workaround: use list comprehensions instead.
  # See: https://github.com/JuliaLang/julia/issues/30796
  # Ws_mat = cat(Ws..., dims=3)
  S = length(Ws)

  plt.figure()
  for i in 1:I
    plt.subplot(I, 1, i)
    # boxplot(Ws_mat[i, :, :]')
    Wis = [Ws[s][i, k] for s in 1:S, k in 1:K]
    boxplot(Wis)
    plt.ylabel("W$(i)")
    if W_true != nothing
      K_true = size(W_true, 2)
      for k in 1:K_true
        plt.axhline(W_true[i, k])
      end
    end
  end
  plt.xlabel(xlabel)
  plt.savefig("$(imgdir)/W.pdf", bbox_inches="tight")
  plt.close()
end

grepKmcmc(s) = match(r"(?<=Kmcmc_)\d+", s)

function make_metrics(different_K_runs_dir, outputfname; thresh,
                      grepKmcmcfn=grepKmcmc)
  metrics = Dict{Int, Any}()

  for (root, dirs, files) in walkdir(different_K_runs_dir)
    for file in files
      if file == outputfname
        # Get path to output
        path_to_output = joinpath(root, file)
        println(path_to_output)

        # Parse Kmcmc
        Kmcmc = parse(Int, grepKmcmcfn(splitdir(root)[2]).match)

        # Load output
        output = BSON.load(path_to_output)

        # Append metrics to metrics dictionary
        # metrics[Kmcmc] = output[:metrics]
        metrics[Kmcmc] = Dict{Symbol, Any}()
        metrics[Kmcmc][:LPML] = output[:LPML]
        metrics[Kmcmc][:DIC] = output[:DIC]

        # Posterior samples of W
        Ws = extract(:theta__W, output[:samples][1])

        # Calibration metric
        cmetric = let
          Wmean = mean(Ws)
          metrics[Kmcmc][:I] = size(Wmean, 1)
          num_wik_lt_thresh = sum(Wmean .< thresh)
        end

        # Append calibration metric to metrics dictionary
        metrics[Kmcmc][:cmetric] = cmetric

        # Posterior samples of R (I x K) -- number of active features / sample
        # Rs (I x num_mcmc_samples)
        Rs = hcat([vec(sum(W .> 0, dims=2)) for W in Ws]...)

        # Mean of R, by sample (vector of length I)
        Rmean = vec(mean(Rs, dims=2))

        # Mean of R, by sample (vector of length I)
        R_025 = quantiles(Rs, .025, dims=2, drop=true)

        # Mean of R, by sample (vector of length I)
        R_975 = quantiles(Rs, .975, dims=2, drop=true)

        # Append Rmean to metrics dictionary
        metrics[Kmcmc][:Rmean] = Rmean

        # Append Rmean to metrics dictionary
        metrics[Kmcmc][:R_025] = R_025

        # Append Rmean to metrics dictionary
        metrics[Kmcmc][:R_975] = R_975
      end
    end
  end

  # All Kmcmcs
  Kmcmcs = sort(collect(keys(metrics)))

  # Create dir for metrics output
  metrics_dir = joinpath(different_K_runs_dir, "metrics")
  mkpath(metrics_dir)

  # Plot Calibration metric
  plt.plot([metrics[K][:cmetric] for K in Kmcmcs],
           [metrics[K][:LPML] for K in Kmcmcs], marker="o")
  for K in Kmcmcs
    plt.text(metrics[K][:cmetric], metrics[K][:LPML], K)
  end
  plt.xlabel(L"number of $W_{i,k}$ < 10%")
  plt.ylabel("LPML")
  plt.savefig(joinpath(metrics_dir, "calibration_metric.pdf"),
              bbox_inches="tight")
  plt.close()

  # Plot LPML, DIC
  for metric in [:LPML, :DIC]
    plt.plot(Kmcmcs, map(K -> metrics[K][metric], Kmcmcs), marker="o")
    plt.xlabel("K")
    plt.ylabel(String(metric))
    plt.savefig(joinpath(metrics_dir, "$(metric).pdf"), bbox_inches="tight")
    plt.close()
  end

  # Plot R
  I = metrics[Kmcmcs[1]][:I]
  plt.figure()
  for i in 1:I
    plt.subplot(I, 1, i)
    plt.plot(Kmcmcs, map(K -> metrics[K][:Rmean][i], Kmcmcs), marker="o")
    plt.fill_between(Kmcmcs,
                     map(K -> metrics[K][:R_025][i], Kmcmcs),
                     map(K -> metrics[K][:R_975][i], Kmcmcs), alpha=.5)
    if i == I
      plt.xlabel("K")
    end
    plt.ylabel("R$(i)")
  end
  plt.savefig(joinpath(metrics_dir, "R.pdf"), bbox_inches="tight")
  plt.close()

  return metrics
end

function plot_dden(; ddens, etas, Ws, Zs, sig2s, deltas,
                   ygrid, imgdir, printmean=true,
                   plot_true=true,
                   simdat=nothing, xlabel="expression level", ylabel="density",
                   dfs=nothing,
                   dden_xlim=(-6, 6))
  # Make directories if needed
  mkpath("$(imgdir)/dden")
  mkpath("$(imgdir)/txt")

  # Cache data sizes
  I, J = size(ddens[1])

  # split etas into etas0, etas1
  etas0 = [x[0] for x in etas]
  etas1 = [x[1] for x in etas]

  # Split mus from deltas
  mus0 = Matrix(hcat([-cumsum(d[0]) for d in deltas]...)')
  mus1 = Matrix(hcat([cumsum(d[1]) for d in deltas]...)')

  for i in 1:I
    for j in 1:J
      print("\r i: $i j: $j  ")

      # Get data density for sample i, marker j.
      # NOTE: `dden_ij` will be a vector (of length `length(ddens)`) of vectors
      # (of length `length(ygrid)`), and represents the posterior density at
      # the corresponding ygrid values.
      dden_ij = [ddij[i, j] for ddij in ddens]

      # 95% CI lower bound for dden ij
      dden_ij_lower = [quantile([ddij[g] for ddij in dden_ij], .025)
                       for g in 1:length(ygrid)]

      # 95% CI upper bound for dden ij
      dden_ij_upper = [quantile([ddij[g] for ddij in dden_ij], .975)
                       for g in 1:length(ygrid)]

      # Plot object of the credible interval for dden[i, j]
      p_ci_obs = plt.fill_between(ygrid, dden_ij_lower, dden_ij_upper,
                                  alpha=.5, color="blue")

      plt.xlabel(xlabel)
      plt.ylabel(ylabel)
      if dden_xlim != nothing
        plt.xlim(dden_xlim)
      end

      # Add eta to posterior dden plot
      eta0_ij_mean = mean(eta[0][i, j, :] for eta in etas)
      L0 = length(eta0_ij_mean)
      plt.scatter(mean(mus0, dims=1), zeros(L0),
                  s=eta0_ij_mean * 60 .+ 10, marker="X", color="green")

      eta1_ij_mean = mean(eta[1][i, j, :] for eta in etas)
      L1 = length(eta1_ij_mean)
      p_mu = plt.scatter(mean(mus1, dims=1), zeros(L1),
                         s=eta1_ij_mean * 60 .+ 10, marker="X",
                         color="green")

      if printmean
        mkpath("$(imgdir)/txt/eta")

        open("$(imgdir)/txt/eta/eta0_i$(i)_j$(j)_mean.txt", "w") do io
          writedlm(io, eta0_ij_mean)
        end

        open("$(imgdir)/txt/eta/eta1_i$(i)_j$(j)_mean.txt", "w") do io
          writedlm(io, eta1_ij_mean)
        end
      end

      if simdat != nothing
        # Plot simulated data truth
        p_yobs = sns.kdeplot(skipnan(simdat[:y][i][:, j]),
                             color="red", bw=.1, label="tmp")
        p_yobs = p_yobs.get_legend_handles_labels()[1][1]

        # Histogram of observed data only
        # plt.hist(skipnan(simdat[:y][i][:, j]), color="red",
        #          alpha=.3, label="y (observed)",
        #          density=true, bins=30)

        if :eta in keys(simdat)
          eta_true = simdat[:eta]
        else
          eta_true = Dict(0 => ones(I, J, simdat[:L][0]),
                          1 => ones(I, J, simdat[:L][1]))
        end

        dgrid = dden_complete(ygrid, simdat[:W], eta_true,
                              simdat[:Z], simdat[:mus],
                              simdat[:sig2], i=i, j=j)

        if plot_true
          p_truth_complete, = plt.plot(ygrid, dgrid, color="grey", ls="--")
        end
      else
        # Plot histogram of observed data
        p_yobs = sns.kdeplot(skipnan(dfs.data.y[i][:, j]),
                             color="red", bw=.1, label="tmp")
        p_yobs = p_yobs.get_legend_handles_labels()[1][1]
      end

      # Number of MCMC samples
      num_mcmc_samples = length(Ws)

      # Create mus (vector)
      mus = [Dict(0 => -cumsum(delta[0]),
                  1 => cumsum(delta[1]))
             for delta in deltas]

      # Plot complete posterior density (obs and imputed)
      dd_complete_post = [dden_complete(ygrid,
                                        Ws[b],
                                        etas[b],
                                        Zs[b],
                                        mus[b],
                                        sig2s[b], i=i, j=j)
                          for b in 1:num_mcmc_samples]

      dd_complete_post = hcat(dd_complete_post...)  # legnth(ygird) x num_mcmc_samples
      dcp_lower = quantiles(dd_complete_post, .025, dims=2, drop=true)
      dcp_upper = quantiles(dd_complete_post, .975, dims=2, drop=true)
      p_ci_complete = plt.fill_between(ygrid, dcp_lower, dcp_upper, alpha=.5,
                                       color="orange")

      # TODO: PICK UP HERE
      if plot_true
        plt.legend([p_truth_complete, p_ci_complete,
                    p_yobs, p_ci_obs, p_mu],
                   ["truth (complete)", "95% CI (complete)",
                    "y (obs)", "95% CI (obs)", PyPlot.L"$\mu^\star$"])
      else
        plt.legend([p_ci_complete,
                    p_yobs, p_ci_obs, p_mu],
                   ["95% CI (complete)",
                    "y (obs)", "95% CI (obs)", PyPlot.L"$\mu^\star$"])
      end
      plt.savefig("$(imgdir)/dden/dden_i$(i)_j$(j).pdf",
                  bbox_inches="tight")
      plt.close()
    end
  end
  println()
end

end  # PlotUtils


text:
# This file defines an API that allows clients to access an ORM
# model without scattering XPath strings throughout the code.
#
# Generally speaking:
#  - extract_X operations extract nodes that descend from a model
#      element with no additional context
#  - operations that lookup various kinds of information from
#      some id or model element. Note that these differ from extract
#      operations in that lookups must query the entire model for
#      for context and often do significant calculation.
#  - predicates of various names (all return Booleans) that query
#      whether model elements (or ids of model elements) satisfy
#      some property
#

fact_name(f) = f["_Name"]
id(r) = r["id"]
name(o) = o["Name"]
ref(r) = r["ref"]
refmode(e) = e["_ReferenceMode"]

by_id(i, path, element) = findfirst(path * "[@id='$i']", element)
by_name(n, path, element) = findfirst(path * "[@Name='$n']", element)
by_ref(r, path, element) = findfirst(path * "[@ref='$r']", element)

ref_of(path, element) = element !== nothing ? ref(findfirst(path, element)) : nothing
refs_of(path, element) = findall(path, element) .|> ref

plays_role_path(r_id) = "/orm:PlayedRoles/orm:Role[@ref='$r_id']/../.."

# Functions that extract information from an ORM model element node
#
extract_conceptual_data_type_id(vt) = ref_of("orm:ConceptualDataType", vt)

extract_preferred_id_uc(etype) = ref_of("orm:PreferredIdentifier", etype)

extract_role_by_id(r_id, ft) = by_id(r_id, "orm:FactRoles/orm:Role", ft)
extract_role_ids(ft) = extract_roles(ft) .|> id

"""
Extracts an array of role ids from the lone role sequence of
  any ORM constraint that ranges over a lone role sequence
  (e.g., UC, mandatory)
"""
extract_role_ids_from_simple(c) = refs_of("orm:RoleSequence/orm:Role", c)
extract_roles(ft) = findall("orm:FactRoles/orm:Role", ft)

extract_uc_ids_of(ft) = refs_of("orm:InternalConstraints/orm:UniquenessConstraint", ft)
extract_mc_ids_of(ft) = refs_of("orm:InternalConstraints/orm:MandatoryConstraint", ft)

# Functions that consult the ORM model to lookup information of various
# kinds about ORM model elements or ids
#

"""
Given a diagram D that should depict a central concept
   and its constellation and model M, verify correct
   applicaiton of the pattern and return the entity type
   that defines the concept

   NOTE: We designate a schema subgraph centered by some concept
   with a constellation of relations that are keyed by that
   concept using ORM Diagrams that are named according to the
   convention "X:concept" where X is the name of the entity
   type that represents the central concept (e.g.,
   "Location:concept" or "Subscription:concept").
"""
function central_concept(D, M)

   concept_name_pattern = Regex("(?<concept>[^(:| )]+):concept")

   concept_match = match(concept_name_pattern, name(D))
   if concept_match !== nothing
      concept_name = concept_match[:concept]
      etype = entity_by_name(concept_name, M)
      if is_drawn_on(etype, D) return etype end
   end
   return nothing
end

constraint_by_id(c_id, M) = by_id(c_id, "orm:Constraints/*", M)

data_type_of_player(r_id, M) =
   extract_conceptual_data_type_id(
      findfirst("orm:Objects/orm:ValueType" * plays_role_path(r_id), M))

entity_by_name(n, M) = by_name(n, "orm:Objects/orm:EntityType", M)

"""
Find the entity for which uc_id is the preferred identifier, or nothing
"""
function entity_identified_by(uc_id, M)
   uc = uc_by_id(uc_id, M)
   uc === nothing && return nothing
   reft = findfirst("orm:PreferredIdentifierFor", uc)
   reft === nothing ? nothing : by_id(ref(reft), "orm:Objects/orm:EntityType", M)
end

entity_that_plays(r_id, M) =
   findfirst("orm:Objects/orm:EntityType" * plays_role_path(r_id), M)

fact_type_with_role(role_id, M) =
   findfirst("orm:Facts/orm:Fact[orm:FactRoles/orm:Role[@id='$role_id']]", M)

"""
Create sequence of ORM fact types that are drawn on diagram D
"""
function fact_types_in_diagram(D, M)
   drawn_fact_type_path = "ormDiagram:Shapes/ormDiagram:FactTypeShape/ormDiagram:Subject"
   global_relations = findfirst("orm:Facts", M)
   # outputs = EzXML.Node[]
   outputs = []
   for shape in findall(drawn_fact_type_path, D)
      ft = by_id(ref(shape), "orm:Fact", global_relations)
      if ft !== nothing
         push!(outputs, ft)
      end
   end
   outputs
end

"""
Given an entity type, looks up the fact type that identifies
  instances of that entity type, if it exists, or nothing
"""
function identifying_fact_type_for(etype, M)
   uc_id = extract_preferred_id_uc(etype)
   uc = uc_by_id(uc_id, M)
   role_ids = extract_role_ids_from_simple(uc)
   if size(role_ids)[1] == 1
      return fact_type_with_role(role_ids[1], M)
   else
      return nothing
   end
end

"""
Retrieves the lone role spanned by the constraint identified by
   c_id, or nothing if the c_id is invalid or if it spans more
   than one role
"""
function lone_role_of(c_id, M)
   c = constraint_by_id(c_id, M)
   c !== nothing || return nothing

   rid_seq = extract_role_ids_from_simple(c)
   size(rid_seq)[1] == 1 ? rid_seq[1] : nothing
end

mc_by_id(mc_id, M) = by_id(mc_id, "orm:Constraints/orm:MandatoryConstraint", M)

role_player_by_id(r_id, M) = findfirst("orm:Objects/*" * plays_role_path(r_id), M)

"""
This function picks out that role id of a fact type that is not covered
   by the given uc. We return an array that will be empty if the uc
   is spanning and will contain exactly one element otherwise
"""
function roles_excluded_by_uc_in_fact_type(uc, ft, M)
   role_seq = String[]
   key_role_id_seq = extract_role_ids_from_simple(uc)
   for role_id in extract_role_ids(ft)
      if findfirst((r_id -> r_id === role_id), key_role_id_seq) !== nothing
         continue
      else
         push!(role_seq, role_id)
      end
   end
   return role_seq
end

"""
Returns the sequence of entity supertypes of etype in M
"""
function super_types_of(etype, M)
   path_prefix = "orm:Facts/orm:SubtypeFact/orm:FactRoles"
   path_suffix = "/orm:SupertypeMetaRole/orm:RolePlayer"
   refs_of("orm:PlayedRoles/orm:SubtypeMetaRole", etype) .|>
      (oid -> findfirst(path_prefix *
                          "[orm:SubtypeMetaRole[@id='$oid']]" *
                          path_suffix,
                        M)) .|>
      ref .|>
      (reft -> by_id(reft, "orm:Objects/orm:EntityType", M))
end

"""
Returns the sequence of entity subtypes of etype in M
"""
function sub_types_of(etype, M)
   path_prefix = "orm:Facts/orm:SubtypeFact/orm:FactRoles"
   path_suffix = "/orm:SubtypeMetaRole/orm:RolePlayer"
   refs_of("orm:PlayedRoles/orm:SupertypeMetaRole", etype) .|>
      (oid -> findfirst(path_prefix *
                          "[orm:SupertypeMetaRole[@id='$oid']]" *
                          path_suffix,
                        M)) .|>
      ref .|>
      (reft -> by_id(reft, "orm:Objects/orm:EntityType", M))
end

uc_by_id(uc_id, M) = by_id(uc_id, "orm:Constraints/orm:UniquenessConstraint", M)

#
# An ORM fact type often gives rise to multiple Rel relations because
#   roles are unordered in a fact type and may be ordered differently
#   to support different readings, while in Rel the ordering of
#   roles in a relation is fixed. We use the internal uniqueness
#   constraints declared over the roles of a fact type to generate
#   relations with definite column orderings. When a fact type contains
#   non-spanning UCs, we often need to associate a UC with the role
#   that it does not span and vice versa.
#

"""
This function extracts that UC of a fact type that does not cover
   the given role.
"""
function uc_that_excludes_role_in_fact_type(role, ft, M)
   role_id = id(role)
   for uc_id in extract_uc_ids_of(ft)
      uc = by_id(uc_id, "orm:Constraints/orm:UniquenessConstraint", M)
      occurrence = by_ref(role_id, "orm:RoleSequence/orm:Role", uc)
      if occurrence === nothing return uc end
   end
   return nothing
end

#
# ORM model element or id predicates
#

function has_named_roles(ft)::Bool
   for role in extract_roles(ft)
      if name(role) != ""
         return true
      end
   end
   return false
end

"""
Check whether the model element is drawn on diagram D
"""
function is_drawn_on(elem, D)::Bool
   elem === nothing && return false

   shape_subject = "ormDiagram:Shapes/*/ormDiagram:Subject"
   by_ref(id(elem), "$shape_subject", D) !== nothing
end

"""
Check whether a fact type is used to provide the preferred identifier
  for some entity type
"""
function is_refmode_fact_type(ft, M)::Bool
   role_ids = extract_role_ids(ft)
   if size(role_ids)[1] == 2
      uc_ids = extract_uc_ids_of(ft)
      if size(uc_ids)[1] == 2
         for uc_id in uc_ids
            # Choose the UC that is the preferred identifier for some
            # entity type, and note that in a well-formed refmode
            # pattern, that UC spans the role that is *not* played
            # by that entity type
            etype = entity_identified_by(uc_id, M)
            if etype !== nothing
               # Find the role spanned by that UC
               identifying_role_id = lone_role_of(uc_id, M)
               identifying_role_id !== nothing || continue

               # Verify that the other role is played by etype and is mandatory
               for rid in role_ids
                  if rid != identifying_role_id
                     rid_player = entity_that_plays(rid, M)
                     if id(rid_player) == id(etype)
                        for mc_id in extract_mc_ids_of(ft)
                           if rid == lone_role_of(mc_id, M)
                              return true
                           end
                        end
                     end
                  end
               end
            end
         end
      end
   end
   return false
end

function is_top_level_concept(etype, M)::Bool
   size(super_types_of(etype, M), 1) == 0
end

identifies_cdt(dt_id, cdt, M)::Bool = by_id(dt_id, "orm:DataTypes/$cdt", M)  !== nothing

function maps_to_rel_string(dt_id, M)::Bool
   identifies_cdt(dt_id, "orm:VariableLengthTextDataType", M) ||
   identifies_cdt(dt_id, "orm:FixedLengthTextDataType", M) ||
   identifies_cdt(dt_id, "orm:AutoCounterNumericDataType", M)
end

function maps_to_rel_int(dt_id, M)::Bool
   identifies_cdt(dt_id, "orm:SignedIntegerNumericDataType", M) ||
   identifies_cdt(dt_id, "orm:TrueOrFalseLogicalDataType", M) ||
   identifies_cdt(dt_id, "orm:UnsignedIntegerNumericDataType", M) ||
   identifies_cdt(dt_id, "orm:UnsignedSmallIntegerNumericDataType", M) ||
   identifies_cdt(dt_id, "orm:TimeTemporalDataType", M)
end

function maps_to_rel_date(dt_id, M)::Bool
   identifies_cdt(dt_id, "orm:DateTemporalDataType", M)
end

function maps_to_rel_datetime(dt_id, M)::Bool
   identifies_cdt(dt_id, "orm:DateAndTimeTemporalDataType", M)
end

max_stars_repo_name: RelationalAI-oss/kse-tools

max_issues_repo_name: RelationalAI-oss/kse-tools

max_forks_repo_name: RelationalAI-oss/kse-tools


text:
module ClipData

using CSV, Tables

using InteractiveUtils: clipboard

export cliptable, cliparray, mwetable, mwearray, @mwetable, @mwearray

"""
    cliptable(; kwargs...)

Make a table from the clipboard. Returns
a `CSV.File`, which can then be transformed
into a `DataFrame` or other Tables.jl-compatible
object.

`cliptable` auto-detects delimiters, and all keyword
arguments are passed directly to the `CSV.File`
constructor.

## Examples

```julia-repl
julia> # Send string to the clipboard
       \"\"\"
       a, b
       1, 2
       100, 200
       \"\"\" |> clipboard

julia> cliptable()
2-element CSV.File{false}:
 CSV.Row: (a = 1,  b = 2)
 CSV.Row: (a = 100,  b = 200)
```

"""
function cliptable(; kwargs...)
    CSV.File(IOBuffer(clipboard()); kwargs...)
end

"""
    cliparray(; kwargs...)

Make a `Vector` or `Matrix` from the clipboard.
Auto-detects delimiters, and all keyword arguments are passed
directly to a `CSV.File` constructor. If the returned `CSV.File`
has one column, `cliparray` returns a `Vector`. Otherwise, it returns
a `Matrix`.

# Examples

```julia-repl
julia> # Send string to clipboard
       \"\"\"
       1 2
       3 4
       \"\"\" |> clipboard

julia> cliparray()
2×2 Matrix{Int64}:
 1  2
 3  4

julia> \"\"\"
       1
       2
       3
       4
       \"\"\" |> clipboard

julia> cliparray()
4-element Vector{Int64}:
 1
 2
 3
 4
```
"""
function cliparray(; kwargs...)
    t = CSV.File(IOBuffer(clipboard()); header=false, kwargs...)
    mat = Tables.matrix(t)
    if size(mat, 2) == 1 || size(mat, 1) == 1
        return vec(mat)
    else
        return mat
    end
end

"""
    cliptable(t; delim = '\t', kwargs...)

Send a Tables.jl-compatible object to the clipboard.
Default delimiter is tab. Accepts all keyword arguments
that can be pased to `CSV.write`.

# Example

```julia-repl
julia> t = (a = [1, 2, 3], b = [100, 200, 300])
(a = [1, 2, 3], b = [100, 200, 300])

julia> cliptable(t)
a   b
1   100
2   200
3   300
```
"""
function cliptable(t; returnstring = false, delim = '\t', kwargs...)
    io = IOBuffer()
    CSV.write(io, t, delim = delim, kwargs...)
    s = chop(String(take!(io)), tail = 1)
    clipboard(s)

    if returnstring == true
      return s
    else
      return nothing
    end
end

"""
    cliparray(t::AbstractVecOrMat; kwargs...)

Send a `Vector` or `Matrix` to the clipboard.
Default delimiter is tab and with no header.
Accepts all keyword arguments that can be passed
to `CSV.write`.

# Examples

```julia-repl
julia> \"\"\"
       1 2
       3 4
       \"\"\" |> clipboard

julia> cliparray()
2×2 Matrix{Int64}:
 1  2
 3  4
```
"""
function cliparray(t::AbstractVecOrMat; returnstring = false, delim='\t',
                   header=false, kwargs...)
    if t isa AbstractVector
        t = reshape(t, :, 1)
    end
    io = IOBuffer()
    CSV.write(io, Tables.table(t); delim=delim, header=header, kwargs...)
    s = chop(String(take!(io)), tail = 1)
    clipboard(s)

    if returnstring == true
      return s
    else
      return nothing
    end
end

"""
    mwetable(; name="df")

Create a Minimum Working Example (MWE) using
the clipboard. `tablmwe` prints out a multi-line
comma-separated string and provides the necessary
code to read that string using `CSV.File`.
The object is assigned the name given by
`name` (default `"df"`).

# Examples

```julia-repl
julia> \"\"\"
       a b
       1 2
       100 200
       \"\"\" |> clipboard

julia> mwetable()
df = \"\"\"
a,b
1,2
100,200
\"\"\" |> IOBuffer |> CSV.File
```

"""
function mwetable(; returnstring=false, name="df")
    t = cliptable()
    mwetable(t, returnstring=returnstring, name=name)
end

"""
    mwetable(t; name="df")

Create a Minimum Working Example (MWE) from
an existing Tables.jl-compatible object.
`tablmwe` prints out a multi-line
comma-separated string and provides the necessary
code to read that string using `CSV.File`.
The object is assigned the name given by
`name` (default `:df`).

# Examples

```julia-repl
julia> t = (a = [1, 2, 3], b = [100, 200, 300])
(a = [1, 2, 3], b = [100, 200, 300])

julia> mwetable(t)
df = \"\"\"
a,b
1,100
2,200
3,300
\"\"\" |> IOBuffer |> CSV.File
```
"""
function mwetable(t; returnstring=false, name="df")
    main_io = IOBuffer()
    table_io = IOBuffer()

    start_str = """
$name = \"\"\"
"""
    print(main_io, start_str)

    CSV.write(table_io, t;)
    print(main_io, String(take!(table_io)))

    end_str = """
\"\"\" |> IOBuffer |> CSV.File"""
    print(main_io, end_str)
    s = String(take!(main_io))

    if returnstring == true
      return s
    else
      return nothing
    end
end

function mwetable_helper(t::Symbol)
    t_name = QuoteNode(t)
    :(mwetable($t, name = $t_name))
end

"""
    @mwetable(t)

Create a Minimum Working Example (MWE) from
an existing Tables.jl-compatible object.
`tablmwe` prints out a multi-line
comma-separated string and provides the necessary
code to read that string using `CSV.File`. The name
assigned to the object in the MWE is the
same as the name of the input object.

# Examples

```julia-repl
julia> my_special_table = (a = [1, 2, 3], b = [100, 200, 300])
(a = [1, 2, 3], b = [100, 200, 300])

julia> @mwetable my_special_table
my_special_table = \"\"\"
a,b
1,100
2,200
3,300
\"\"\" |> IOBuffer |> CSV.File
```
"""
macro mwetable(t)
    esc(mwetable_helper(t))
end

"""
    mwearray(; name=:X)

Create a Minimum Working Example (MWE) from
the clipboard to create an array. `mwearray`
returns the a multi-line string with the
code necessary to read the string stored in
clipboard as a `Vector` or `Matrix`.

# Examples

```julia-repl
julia> \"\"\"
       1 2
       3 4
       \"\"\" |> clipboard

julia> mwearray()
X = \"\"\"
1,2
3,4
\"\"\" |> IOBuffer |> CSV.File |> Tables.matrix
```
"""
function mwearray(; returnstring=false, name=nothing)
    t = cliparray()
    name= t isa AbstractVector ? :x : :X
    mwearray(t, returnstring=returnstring, name=name)
end

"""
    mwearray(t::AbstractMatrix; name=:X)

Create a Minimum Working Example (MWE) from
a `Matrix`. `mwearray`
returns the a multi-line string with the
code necessary to recreate `t`.

# Examples

```julia-repl
julia> X = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> mwearray(X)
X = \"\"\"
1,2
3,4
\"\"\" |> IOBuffer |> CSV.File |> Tables.matrix
```
"""
function mwearray(t::AbstractMatrix; returnstring=false, name=:X)
    main_io = IOBuffer()
    array_io = IOBuffer()

    start_str = """
$name = \"\"\"
"""
    print(main_io, start_str)

    CSV.write(array_io, Tables.table(t); header=false)
    print(main_io, String(take!(array_io)))

    end_str = """
\"\"\" |> IOBuffer |> CSV.File |> Tables.matrix"""
    print(main_io, end_str)
    s = String(take!(main_io))

    if returnstring == true
      return s
    else
      return nothing
    end
end

"""
    mwearray(t::AbstractVector; name=:x)

Create a Minimum Working Example (MWE) from
a `Vector`. `mwearray`
returns the a multi-line string with the
code necessary to recreate `t`.

# Example

```julia-repl
julia> mwearray(x; name=:x)
x = \"\"\"
1
2
3
4
\"\"\" |> IOBuffer |> CSV.File |> Tables.matrix |> vec
```
"""
function mwearray(t::AbstractVector; returnstring=false, name=:x)
    main_io = IOBuffer()
    array_io = IOBuffer()

    t = reshape(t, :, 1)

    start_str = """
$name = \"\"\"
"""
    print(main_io, start_str)

    CSV.write(array_io, Tables.table(t); header=false)
    print(main_io, String(take!(array_io)))

    end_str = """
\"\"\" |> IOBuffer |> CSV.File |> Tables.matrix |> vec"""
    print(main_io, end_str)
    s = String(take!(main_io))

    if returnstring == true
      return s
    else
      return nothing
    end
end

function mwearray_helper(t::Symbol)
    t_name = QuoteNode(t)
    :(mwearray($t, name=$t_name))
end

"""
    @mwearray(t)

Create a Minimum Working Example (MWE)
from a `Vector` or `Matrix` with the same
name as the object in the Julia session.

# Examples

```julia-repl
julia> my_special_matrix = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> @mwearray my_special_matrix
my_special_matrix = \"\"\"
1,2
3,4
\"\"\" |> IOBuffer |> CSV.File |> Tables.matrix
```
"""
macro mwearray(t)
    esc(mwearray_helper(t))
end

end

max_stars_repo_name: pdeffebach/DataClip

max_issues_repo_name: pdeffebach/DataClip

max_forks_repo_name: pdeffebach/DataClip


text:
"""
**Representation of a GBIF taxon**

All taxonomic level fields can either be `missing`, or a pair linking the name
of the taxon/level to its unique key in the GBIF database.

`name` - the vernacular name of the taxon

`scientific` - the accepted scientific name of the species

`status` - the status of the taxon

`match` - the type of match

`kingdom` - a `Pair` linking the name of the kingdom to its unique ID

`phylum` - a `Pair` linking the name of the phylum to its unique ID

`class` - a `Pair` linking the name of the class to its unique ID

`order` - a `Pair` linking the name of the order to its unique ID

`family` - a `Pair` linking the name of the family to its unique ID

`genus` - a `Pair` linking the name of the genus to its unique ID

`species` - a `Pair` linking the name of the species to its unique ID

`confidence` - an `Int64` to note the confidence in the match

`synonym` - a `Boolean` indicating whether the taxon is a synonym
"""
struct GBIFTaxon
   name::AbstractString
   scientific::AbstractString
   status::Union{Missing, Symbol}
   match::Union{Missing, Symbol}
   kingdom::Union{Missing, Pair{String, Int64}}
   phylum::Union{Missing, Pair{String, Int64}}
   class::Union{Missing, Pair{String, Int64}}
   order::Union{Missing, Pair{String, Int64}}
   family::Union{Missing, Pair{String, Int64}}
   genus::Union{Missing, Pair{String, Int64}}
   species::Union{Missing, Pair{String, Int64}}
   confidence::Union{Missing, Int64}
   synonym::Bool
end

function GBIFTaxon(o::Dict{String, Any})
   levels = ["kingdom", "phylum", "class", "order", "family", "genus", "species"]
   r = Dict{Any,Any}()
   for l in levels
      if haskey(o, l)
         r[l] = Pair(o[l], o[l*"Key"])
      else
         r[l] = missing
      end
   end
   return GBIFTaxon(
      o["canonicalName"],
      o["scientificName"],
      haskey(o, "status") ? Symbol(o["status"]) : missing,
      haskey(o, "matchType") ? Symbol(o["matchType"]) : missing,
      r["kingdom"],
      r["phylum"],
      r["class"],
      r["order"],
      r["family"],
      r["genus"],
      r["species"],
      get(o, "confidence", missing),
      o["synonym"]
   )
end


text:
"""
    plot_surfaces(femesh, compartments = 1:ncompartment)

Plot surfaces.
"""
function plot_surfaces(surfaces, boundaries = 1:maximum(surfaces.facetmarkers))
    (; points, facets, facetmarkers) = surfaces
    scene = nothing
    first = true
    for b ∈ boundaries
        f = facets[:, facetmarkers .== b]
        color = points[3, :]
        if first
            scene = poly(points', f'; color, strokewidth = 1, shading = false)
            first = false
        else
            poly!(points', f'; color, strokewidth = 1, shading = false)
        end
    end
    scene
end


text:
const SCALING_FACTOR = 1_000_000

# use specialty type for "read counts"
# with basic function: 'sum' that can be
# overridden for more complex count bias models
#const ReadCount = JointBiasCounter #DefaultCounter
const DEFAULT_DEF   = (v = zero(DefaultCounter); v.isadjusted = true; v)
const DEFAULT_JOINT = (v = zero(JointBiasCounter); v.isadjusted = true; v)

const DEF_READCOUNT = 0.0
const DEF_READVALUE = 1.0

default(::Type{DefaultCounter}) = DEFAULT_DEF
default(::Type{JointBiasCounter}) = DEFAULT_JOINT

# Use this abstraction to store SGAlignPaths before
# quantification or assignment
abstract type SGAlignContainer end

struct SGAlignSingle <: SGAlignContainer
   fwd::SGAlignPath
end

struct SGAlignPaired <: SGAlignContainer
   fwd::SGAlignPath
   rev::SGAlignPath
end

# Store unique path's for gene/node but ignore
# SGAlignScore which is irrelevant at this stage
Base.isequal( a::SGAlignSingle, b::SGAlignSingle ) = a.fwd == b.fwd
Base.isequal( a::SGAlignPaired, b::SGAlignPaired ) = a.fwd == b.fwd && a.rev == b.rev

Base.:(==)( a::SGAlignSingle, b::SGAlignSingle ) = a.fwd == b.fwd
Base.:(==)( a::SGAlignPaired, b::SGAlignPaired ) = a.fwd == b.fwd && a.rev == b.rev

Base.hash( v::SGAlignSingle ) = hash( v.fwd )
Base.hash( v::SGAlignSingle, h::UInt64 ) = hash( v.fwd, h )
Base.hash( v::SGAlignPaired ) = hash( v.fwd, hash( v.rev ) )
Base.hash( v::SGAlignPaired, h::UInt64 ) = hash( v.fwd, hash( v.rev, h ) )

ispaired( cont::SGAlignSingle ) = false
ispaired( cont::SGAlignPaired ) = true

# check if one alignment path is a subset of another
function Base.in( first::SGAlignPath, last::SGAlignPath )
   for i in 1:(length(last)-length(first))+1
      offset = i - 1
      match = true
      for j in 1:length(first)
         if !(first[j] == last[offset+j])
            match = false
            break
         end
      end
      if match
         return true
      end
   end
   return false
end

# reimplementation of to edge code in events.jl
# we require not just the nodes to exist in the intset
# but also that they are adjacent in the intset
function Base.in( first, last, is::IntSet )
   if first in is &&
      last in is
      for i in first+1:last-1
         if i in is
            return false
         end
      end
      return true
   end
   return false
end

Base.in( edge::I, is::IntSet ) where I <: IntervalValue = in( edge.first, edge.last, is )

# This looks for an edge in a Vector of IntSets using ^
function Base.in( edge::I, viset::Vector{IntSet} ) where I <: IntervalValue
   for iset in viset
      if edge in iset
         return true
      end
   end
   false
end

# For events.jl:
function inall( edge::I, viset::Vector{IntSet} ) where I <: IntervalValue
   inone = false
   for iset in viset
      if edge.first >= first(iset) && edge.last <= last(iset)
         if !(edge in iset)
            return false
         end
         inone = true
      end
   end
   inone ? true : false
end

# each edge in the path must exist in the int set to be true
function Base.in( path::SGAlignPath, is::IntSet )
   if length(path) == 1
      if !(path[1].node in is)
         return false
      end
   else
      for i in 1:length(path)-1
         if !in( path[i].node, path[i+1].node, is )
            return false
         end
      end
   end
   return true
end

Base.in( aln::SGAlignSingle, is::IntSet ) = in( aln.fwd, is )
Base.in( aln::SGAlignPaired, is::IntSet ) = in( aln.fwd, is ) && in( aln.rev, is )

# This is where we count reads for nodes/edges/circular-edges/effective_lengths
# bias is an adjusting multiplier for nodecounts to bring them to the same level
# as junction counts, which should always be at a lower level
mutable struct SpliceGraphQuant{C <: SGAlignContainer, R <: ReadCounter}
   node::Vector{R}
   edge::IntervalMap{ExonInt,R}
   long::Dict{C,R}
   circ::Dict{Tuple{ExonInt,ExonInt},R}
   leng::Vector{Float64}
   bias::Float64

   function SpliceGraphQuant{C,R}( sg::SpliceGraph ) where {C <: SGAlignContainer, R <: ReadCounter}
      return new( R[zero(R) for x in 1:length(sg.nodelen)],
             IntervalMap{ExonInt,R}(),
             Dict{C,R}(),
             Dict{Tuple{ExonInt,ExonInt},R}(),
             zeros( length(sg.nodelen) ), 1.0 )
   end
end

abstract type EquivalenceClass end

const ClassInt  = Int32
const ClassType = Vector{ClassInt}

# type-stable isoform compatibility class
mutable struct IsoCompat <: EquivalenceClass
   class::ClassType
   prop::Vector{Float64}
   prop_sum::Float64
   count::Float64
   isdone::Bool

   function IsoCompat( arr::Vector{Int}, count::Float64 )
      return new( arr, ones( length(arr) ) / length(arr), 1.0, count, false )
   end
end

# Here we store whole graphome quantification
struct GraphLibQuant{C <: SGAlignContainer, R <: ReadCounter}
   tpm::Vector{Float64}                 # isoform tpm
   count::Vector{Float64}               # isoform counts
   length::Vector{Int64}                # isoform lengths
   geneidx::Vector{Int64}               # 0-based offset of isoforms for gene
   genetpm::Vector{Float64}             # sum of isoform tpms for each gene
   quant::Vector{SpliceGraphQuant{C,R}} # splice graph quant structs for gene
   classes::Vector{IsoCompat}           # isoform compatibility classes

   function GraphLibQuant{C,R}( lib::GraphLib ) where {C <: SGAlignContainer, R <: ReadCounter}
      isonum = zeros(Int, length(lib.graphs))
      for i in 1:length(lib.graphs)
         isonum[i] = length(lib.graphs[i].annoname)
      end
      isolen  = sum( isonum )
      tpm     = zeros( isolen )
      count   = zeros( isolen )
      leng    = zeros( isolen )
      geneoff = zeros(Int, length(lib.graphs))
      genetpm = zeros( length(lib.graphs) )
      quant   = Vector{SpliceGraphQuant}( length(lib.graphs) )
      classes = Vector{IsoCompat}()
      cumul_i = 0
      for i in 1:length(lib.graphs)
         geneoff[i] = cumul_i
         for j in 1:length(lib.graphs[i].annopath)
            leng[cumul_i+j] = sum(map( y->lib.graphs[i].nodelen[y], collect(lib.graphs[i].annopath[j]) ))
         end
         cumul_i += isonum[i]
         quant[i] = SpliceGraphQuant{C,R}( lib.graphs[i] )
      end
      new( tpm, count, leng, geneoff, genetpm, quant, classes )
   end
end

# Adjust each MultiCompat raw value and assign its adjusted to .count
function adjust!( sgquant::SpliceGraphQuant{C,R}, mod::B, func=adjust! ) where {C <: SGAlignContainer, R <: ReadCounter, B <: BiasModel}
   for node in sgquant.node
      func( node, mod )
   end
   for edg in sgquant.edge
      func( edg.value, mod )
   end
   for lng in values(sgquant.long)
      func( lng, mod )
   end
   for cir in values(sgquant.circ)
      func( cir, mod )
   end
   sgquant
end

primer_adjust!( sgquant::SpliceGraphQuant{C,R}, mod::B ) where {C <: SGAlignContainer, R <: ReadCounter, B <: BiasModel} = adjust!( sgquant, mod, primer_adjust! )
gc_adjust!( sgquant::SpliceGraphQuant{C,R}, mod::B ) where {C <: SGAlignContainer, R <: ReadCounter, B <: BiasModel} = adjust!( sgquant, mod, gc_adjust! )

# Adjust each sgquant storage value
function adjust!( gquant::GraphLibQuant{C,R}, mod::B, func=adjust! ) where {C <: SGAlignContainer, R <: ReadCounter, B <: BiasModel}
   for sgq in gquant.quant
      func( sgq, mod )
   end
   gquant
end

primer_adjust!( gquant::GraphLibQuant{C,R}, mod::B ) where {C <: SGAlignContainer, R <: ReadCounter, B <: BiasModel} = adjust!( gquant, mod, primer_adjust! )
gc_adjust!( gquant::GraphLibQuant{C,R}, mod::B ) where {C <: SGAlignContainer, R <: ReadCounter, B <: BiasModel} = adjust!( gquant, mod, gc_adjust! )

function gc_normalize!( mod::JointBiasMod, lib::GraphLib, quant::GraphLibQuant )
   for i in 1:length(lib.graphs)
      sg  = lib.graphs[i]
      idx = quant.geneidx[i]
      for j in 1:length(sg.annobias)
         add_to_back!( mod.gcback, sg.annobias[j], quant.tpm[j+idx] )
      end
   end
   foresum = sum(mod.gcfore)
   backsum = sum(mod.gcback)
   @inbounds for i in 1:length(mod.gcfore)
      @fastmath mod.gcfore[i] /= foresum
      @fastmath mod.gcback[i] /= backsum
   end
   mod
end

@inbounds function set_gene_tpm!( quant::GraphLibQuant, lib::GraphLib )
   for i in 1:length(lib.graphs)
      tpm = 0.0
      curidx = quant.geneidx[i]
      for j in 1:length( lib.graphs[i].annoname )
         tpm += quant.tpm[ curidx + j ]
      end
      quant.genetpm[i] = tpm
   end
end

assign_path!( graphq::GraphLibQuant{SGAlignSingle}, path::SGAlignSingle,
              value, temp_iset::IntSet ) = assign_path!( graphq, path, value )

# This function re-count!'s SGAlignPaths which were not count!'ed originally
@inbounds function assign_path!( graphq::GraphLibQuant{SGAlignSingle}, path::SGAlignSingle, value )
   const init_gene = path.fwd[1].gene
   const sgquant = graphq.quant[ init_gene ]

   if length(path.fwd) == 1
      push!( sgquant.node[ path.fwd[1].node ], value )
   else
      for i in 1:length(path.fwd)-1
         const lnode = path.fwd[i].node
         const rnode = path.fwd[i+1].node
         if lnode < rnode
            interv = Interval{ExonInt}( lnode, rnode )
            pushzero!( sgquant.edge, interv, value )
         elseif lnode >= rnode
            pushzero!( sgquant.circ, (lnode,rnode), value )
         end
      end
   end
end

# This function re-counts paired SGAlignPaths that weren't counted because they were too long, disjointed, or both
@inbounds function assign_path!( graphq::GraphLibQuant{SGAlignPaired}, path::SGAlignPaired, value, temp_iset::IntSet=IntSet() )
   const init_gene = path.fwd[1].gene
   const sgquant = graphq.quant[ init_gene ]

   if length(path.fwd) == 1
      push!( sgquant.node[ path.fwd[1].node ], value )
      push!( temp_iset, path.fwd[1].node )
   else
      for i in 1:length(path.fwd)-1
         const lnode = path.fwd[i].node
         const rnode = path.fwd[i+1].node
         push!(temp_iset, lnode)
         push!(temp_iset, rnode)
         if lnode < rnode
            interv = Interval{ExonInt}( lnode, rnode )
            pushzero!( sgquant.edge, interv, value )
         elseif lnode >= rnode
            pushzero!( sgquant.circ, (lnode,rnode), value )
         end
      end
   end

   if length(path.rev) == 1
      if !(path.rev[1].node in temp_iset)
         push!( sgquant.node[ path.rev[1].node ], value )
      end
   else
      for i in 1:length(path.rev)-1
         const lnode = path.rev[i].node
         const rnode = path.rev[i+1].node
         (lnode in temp_iset && rnode in temp_iset) && continue
         if lnode < rnode
            interv = Interval{ExonInt}( lnode, rnode )
            pushzero!( sgquant.edge, interv, value )
         elseif lnode >= rnode
            pushzero!( sgquant.circ, (lnode,rnode), value )
         end
      end
   end
   empty!(temp_iset)
end

# build compatibility classes from SpliceGraphQuant node and edge "counts"
@inbounds function build_equivalence_classes!( graphq::GraphLibQuant, lib::GraphLib; assign_long=true )
   # initialize re-used data structures
   iset = IntSet()
   resize!(iset.bits, 64)
   temp = Dict{Vector{Int},Float64}()

   @inline function handle_iset_value!( idx, value )
      if length(iset) == 1
         graphq.count[first(iset) + idx] += value
      elseif length(iset) > 1
         arr = collect(iset)
         if haskey(temp, arr)
            temp[arr] += value
         else
            temp[arr] = value
         end
      end
   end

   # go through splice graphs (to re-use temp data structures)
   for i in 1:length(lib.graphs)
      # go through single node mapping reads
      for n in 1:length(graphq.quant[i].node)
         for p in 1:length(lib.graphs[i].annopath)
            if n in lib.graphs[i].annopath[p]
               push!(iset, p)
            end
         end
         handle_iset_value!( graphq.geneidx[i], get(graphq.quant[i].node[n]) )
         empty!(iset)
      end
      # go through edges
      for edg in graphq.quant[i].edge
         for p in 1:length(lib.graphs[i].annopath)
            if edg in lib.graphs[i].annopath[p]
               push!(iset, p)
            end
         end
         handle_iset_value!( graphq.geneidx[i], get(edg.value) )
         empty!(iset)
      end
      # go through multi-edge crossing paths
      for align in keys(graphq.quant[i].long)
         for p in 1:length(lib.graphs[i].annopath)
            if align in lib.graphs[i].annopath[p]
               push!(iset, p)
            end
         end
         handle_iset_value!( graphq.geneidx[i], get(graphq.quant[i].long[align]) )
         empty!(iset)
         if assign_long
            assign_path!( graphq, align, get(graphq.quant[i].long[align]), iset )
         end
      end
      # create compatibility classes for isoforms
      for arr in keys(temp)
         if temp[arr] > 0.0
            push!( graphq.classes, IsoCompat(arr + graphq.geneidx[i], temp[arr]) )
         end
      end
      empty!(temp)
   end
end

# This function takes a vector of alignments and returns the compatibility class
# and a boolean describing whether or not we can assign the read through EM
@inbounds function equivalence_class( graphq::GraphLibQuant,
                                      lib::GraphLib,
                                      cont::Vector{C},
                                      temp_iset::IntSet=IntSet() ) where C <: SGAlignContainer
   matches_all = true
   for path in cont
      const g = path.fwd[1].gene
      has_match = false
      for p in 1:length(lib.graphs[g].annopath)
         if path in lib.graphs[g].annopath[p]
            push!(temp_iset, p + graphq.geneidx[g])
            has_match = true
         end
      end
      if !has_match
         matches_all = false
         empty!(temp_iset)
         return ClassType(),false
      end
   end
   retval = collect(ClassInt, temp_iset)
   empty!(temp_iset)
   retval,matches_all
end

## MultiAlignment containers & Equivalence classes

mutable struct MultiCompat{R} <: EquivalenceClass
   class::ClassType
   prop::Vector{Float64}
   prop_sum::Float64
   count::Float64
   isdone::Bool
   postassign::Bool
   raw::R

   function MultiCompat{R}( graphq::GraphLibQuant{C,R}, lib::GraphLib, cont::Vector{C}, temp_iset::IntSet=IntSet() ) where {C <: SGAlignContainer, R <: ReadCounter}
      class,matches = equivalence_class( graphq, lib, cont, temp_iset )
      new( class, ones(length(class)) / length(class), 1.0, 1.0, !matches, !matches, zero(R) )
   end
end

# This hash structure stores multi-mapping
# equivalence classes
struct MultiMapping{C <: SGAlignContainer, R <: ReadCounter}
   map::Dict{Vector{C},MultiCompat{R}}
   iset::IntSet

   function MultiMapping{C,R}() where {C <: SGAlignContainer, R <: ReadCounter}
      new( Dict{Vector{C},MultiCompat{R}}(), IntSet() )
   end
end

function total_multi( mm::MultiMapping )
   total = 0.0
   for v in values(mm.map)
      total += v.count
   end
   total
end

function Base.push!( ambig::MultiMapping{SGAlignSingle,R}, alns::Vector{SGAlignment}, value,
                     graphq::GraphLibQuant, lib::GraphLib ) where R <: ReadCounter
   cont = Vector{SGAlignSingle}( length(alns) )
   for i in 1:length(alns)
      cont[i] = SGAlignSingle( alns[i].path )
   end
   if haskey( ambig.map, cont )
      push!( ambig.map[cont].raw, value )
   else
      mc = MultiCompat{R}( graphq, lib, cont, ambig.iset )
      push!( mc.raw, value )
      ambig.map[cont] = mc
   end
end

function Base.push!( ambig::MultiMapping{SGAlignPaired,R}, fwd::Vector{SGAlignment}, rev::Vector{SGAlignment}, value,
                     graphq::GraphLibQuant, lib::GraphLib ) where R <: ReadCounter
   cont = Vector{SGAlignPaired}( length(fwd) )
   for i in 1:length(fwd)
      cont[i] = SGAlignPaired( fwd[i].path, rev[i].path )
   end
   if haskey( ambig.map, cont )
      push!( ambig.map[cont].raw, value )
   else
      mc = MultiCompat{R}( graphq, lib, cont, ambig.iset )
      push!( mc.raw, value )
      ambig.map[cont] = mc
   end
end

# Adjust each MultiCompat raw value and assign its adjusted to .count
function adjust!( multi::MultiMapping{C,R}, mod::B, func=adjust! ) where {C <: SGAlignContainer, R <: ReadCounter, B <: BiasModel}
   for mc in values(multi.map)
      func( mc.raw, mod )
      mc.count = get(mc.raw)
   end
   multi
end

primer_adjust!( multi::MultiMapping{C,R}, mod::B ) where {C <: SGAlignContainer, R <: ReadCounter, B <: BiasModel} = adjust!( multi, mod, primer_adjust! )
gc_adjust!( multi::MultiMapping{C,R}, mod::B ) where {C <: SGAlignContainer, R <: ReadCounter, B <: BiasModel} = adjust!( multi, mod, gc_adjust! )

# for a given alignment container, set the temp_iset in multi to the compatibility class
@inbounds function set_equivalence_class!( multi::MultiMapping{C,R}, graphq::GraphLibQuant,
                                           lib::GraphLib, aln::C ) where {C <: SGAlignContainer, R <: ReadCounter}
   empty!(multi.iset)
   const g = aln.fwd[1].gene
   for p in 1:length(lib.graphs[g].annopath)
      if aln in lib.graphs[g].annopath[p]
         push!(multi.iset, p + graphq.geneidx[g])
      end
   end
end

function get_class_value( iset::IntSet, compat::MultiCompat )
   value = 0.0
   @inbounds for i in 1:length(compat.class)
      if compat.class[i] in iset
         value += compat.prop[i]
      end
   end
   value
end

function assign_ambig!( graphq::GraphLibQuant{C,R}, lib::GraphLib, ambig::MultiMapping{C,R} ) where {C <: SGAlignContainer, R <: ReadCounter}
   @inbounds for (vp,mc) in ambig.map
      # store relative weight of each alignment path
      vals = zeros(length(vp))

      if mc.postassign  # not all paths match annotated transcripts, use gene level values to assign
         total_tpm = 0.0
         for i in 1:length(vp)
            vals[i] = graphq.genetpm[ vp[i].fwd[1].gene ]
            total_tpm += vals[i]
         end
         @fastmath vals /= sum(vals)
      else # assign based on transcripts
         # since two alignments for a read can match the same transcript
         # lets sum the prop values of the MultiCompat for each alignment
         # then we normalize by the total sum of those, and assign reads accordingly.
         # Example: imagine a single read matches the same transcript in two positions
         # that creates a single MultiCompat object that has two alignments hitting the same
         # transcript class.  This method would assign that prop to both, and then re-normalize
         # so that the total count of the MultiCompat is assigned for the sum of the alignments
         vals = zeros(length(vp))
         for i in 1:length(vp)
            set_equivalence_class!( ambig, graphq, lib, vp[i] )
            vals[i] = get_class_value( ambig.iset, mc )
         end
         @fastmath vals /= sum(vals)
      end

      # now assign alignment paths proportional to their weight
      for i in 1:length(vp)
         assign_path!( graphq, vp[i], vals[i] * mc.count, ambig.iset )
      end
   end
end

# Count a single alignment in quant
function count!( graphq::GraphLibQuant{C,R}, align::SGAlignment, value ) where {C <: SGAlignContainer, R <: ReadCounter}
   align.isvalid == true || return
   init_gene = align.path[1].gene
   sgquant   = graphq.quant[ init_gene ]

   # single node support
   if length(align.path) == 1
      # access node -> [ SGNode( gene, *node* ) ]
      push!( sgquant.node[ align.path[1].node ], value )
   else # multi-node mapping read
      # if any node in path maps to other gene, ignore read
      for n in 1:length(align.path)
         if align.path[n].gene != init_gene
            return
         end
      end

      # add a single edge support
      if length(align.path) == 2
         lnode = align.path[1].node
         rnode = align.path[2].node
         if lnode < rnode
            interv = Interval{ExonInt}( lnode, rnode )
            pushzero!( sgquant.edge, interv, value )
         elseif lnode >= rnode
            pushzero!( sgquant.circ, (lnode,rnode), value )
         end
      else # or we have long path, add long count
         curkey = SGAlignSingle(align.path)
         if haskey(sgquant.long, curkey)
            push!( sgquant.long[curkey], value )
         else
            entry = zero(R)
            push!( entry, value )
            sgquant.long[curkey] = entry
         end
      end
   end
end

## Extension of count! for paired end counting
function count!( graphq::GraphLibQuant{C,R}, fwd::SGAlignment, rev::SGAlignment, value ) where {C <: SGAlignContainer, R <: ReadCounter}
   (fwd.isvalid == true && rev.isvalid == true) || return
   init_gene = fwd.path[1].gene
   rev_gene  = rev.path[1].gene
   (init_gene == rev_gene) || return
   sgquant = graphq.quant[ init_gene ]
   singlebool = false

   if length(fwd.path) <= length(rev.path)
      if fwd.path in rev.path
         singlebool = true
         single = rev.path
      end
   else
      if rev.path in fwd.path
         singlebool = true
         single = fwd.path
      end
   end

   # single node support
   if singlebool
      if length(single) == 1
         # access node -> [ SGNode( gene, *node* ) ]
         push!( sgquant.node[ single[1].node ], value )
      else # multi-node mapping read
         # if any node in path maps to other gene, ignore read
         for n in 1:length(single)
            if single[n].gene != init_gene
               return
            end
         end

         # add a single edge support
         if length(single) == 2
            lnode = single[1].node
            rnode = single[2].node
            if lnode < rnode
               interv = Interval{ExonInt}( lnode, rnode )
               pushzero!( sgquant.edge, interv, value )
            elseif lnode >= rnode
               pushzero!( sgquant.circ, (lnode,rnode), value )
            end
         end
      end
   else
       # or we have two paths, add long count
       curkey = SGAlignPaired(fwd.path, rev.path)
       if haskey(sgquant.long, curkey)
          push!( sgquant.long[curkey], value )
       else
          entry = zero(R)
          push!( entry, value )
          sgquant.long[curkey] = entry
       end
   end
end

@inline function calculate_tpm!( quant::GraphLibQuant, counts::Vector{Float64}=quant.count; readlen::Int64=50, sig::Int64=1 )
   for i in 1:length(counts)
      @fastmath quant.tpm[ i ] = counts[i] / max( (quant.length[i] - readlen), 1.0 )
   end
   rpk_sum = max( sum( quant.tpm ), 1.0 )
   for i in 1:length(quant.tpm)
      if sig > 0
         @fastmath quant.tpm[i] = round( quant.tpm[i] * SCALING_FACTOR / rpk_sum, sig )
      else
         @fastmath quant.tpm[i] = ( quant.tpm[i] * SCALING_FACTOR / rpk_sum )
      end
   end
end

calculate_bias!( sgquant::SpliceGraphQuant ) = 1.0
global_bias( graphq::GraphLibQuant ) = 1.0,0.0

# Every exon-exon junction has an effective mappable space of readlength - K*2.
# Since we only count nodes that don't map over edges, we can give
# give each junction an effective_length of one, and now the space
# inside of a node for which a read could map without overlapping a junction
# is put into units of junction derived effective_length
# kadj is minimum of (readlength - minalignlen) and k-1
@inline function eff_length( node, sg::SpliceGraph, eff_len::Int, kadj::Int )
   len = sg.nodelen[node] + (istxstart( sg.edgetype[node] ) ? 0 : kadj) +
                            (istxstop( sg.edgetype[node+1] ) ? 0 : kadj)
   @fastmath len / eff_len
end

function eff_lengths!( sg::SpliceGraph, sgquant::SpliceGraphQuant, eff_len::Int, kadj::Int )
   for i in 1:length( sg.nodelen )
      sgquant.leng[i] = sg.nodelen[i] #eff_length( i, sg, eff_len, kadj )
   end
end

function effective_lengths!( lib::GraphLib, graphq::GraphLibQuant, eff_len::Int, kadj::Int )
   for i in 1:length( lib.graphs )
      eff_lengths!( lib.graphs[i], graphq.quant[i], eff_len, kadj )
   end
end

function Base.unsafe_copy!( dest::Vector{T}, src::Vector{T}; indx_shift=0 ) where T <: Number
   @inbounds for i in 1:length(src)
      dest[i+indx_shift] = src[i]
   end
end

function copy_isdone!( dest::Vector{T}, src::Vector{T}, denominator::Float64=1.0, sig::Int=4 ) where T <: AbstractFloat
   isdifferent = false
   for i in 1:length(dest)
      val = round( src[i] / denominator, sig )
      if val != dest[i]
         isdifferent = true
      end
      dest[i] = isnan(val) ? 0.0 : val
   end
   !isdifferent || sum(denominator) == 0
end

# This function performs expectation maximization
# and sets the quant.tpm array as the proposed expression set
# at each iteration.   It also requires that calculate_tpm! be
# run initially once prior to gene_em!() call.

function gene_em!( quant::GraphLibQuant, ambig::MultiMapping{C,R};
                   it::Int64=1, maxit::Int64=1000,
                   sig::Int64=4, readlen::Int64=50 ) where {C <: SGAlignContainer, R <: ReadCounter}

   count_temp = ones(length(quant.count))
   tpm_temp   = ones(length(quant.count))
   prop_temp  = zeros( 500 )
   uniqsum    = sum(quant.count)
   ambigsum   = length(ambig.map)

   @inline function maximize_and_assign!( eq::E, maximize::Bool=true ) where E <: EquivalenceClass
      eq.isdone && return
      if maximize
         eq.prop_sum = 0.0
         c = 1
         for i in eq.class
            cur_tpm = quant.tpm[i]
            prop_temp[c] = cur_tpm
            eq.prop_sum += cur_tpm
            c += 1
         end
         eq.isdone = copy_isdone!( eq.prop, prop_temp, eq.prop_sum )
      end
      c = 1
      for i in eq.class
         if eq.isdone
            quant.count[i] += eq.prop[c] * eq.count
         end
         count_temp[i] += eq.prop[c] * eq.count
         c += 1
      end
   end

   while tpm_temp != quant.tpm && it < maxit

      unsafe_copy!( count_temp, quant.count )
      unsafe_copy!( tpm_temp, quant.tpm )

      for eq in values(ambig.map)
         maximize_and_assign!( eq, (it > 1) )
      end
      for eq in quant.classes
         maximize_and_assign!( eq, (it > 1) )
      end

      calculate_tpm!( quant, count_temp, sig=sig, readlen=readlen )

      it += 1
   end
   it
end

function output_tpm( file, lib::GraphLib, gquant::GraphLibQuant )
   geneio = open( file * ".gene.tpm.gz", "w" )
   genestream = ZlibDeflateOutputStream( geneio )
   isoio = open( file * ".isoform.tpm.gz", "w" )
   isostream = ZlibDeflateOutputStream( isoio )
   output_tpm_header( genestream, "Gene" )
   output_tpm_header( isostream, "Isoform" )
   for i in 1:length(lib.names)
      idx = gquant.geneidx[i]
      genetpm = 0.0
      genecnt = 0.0
      for j in 1:length(lib.graphs[i].annoname)
         genetpm += gquant.tpm[j+idx]
         genecnt += gquant.count[j+idx]
         tab_write( isostream, lib.graphs[i].annoname[j] )
         tab_write( isostream, string(gquant.tpm[j+idx]) )
         tab_write( isostream, string(gquant.count[j+idx]) )
         write( isostream, '\n' )
      end
      tab_write( genestream, lib.names[i] )
      tab_write( genestream, string(genetpm) )
      tab_write( genestream, string(genecnt) )
      write( genestream, '\n' )
   end
   close( isostream )
   close( isoio )
   close( genestream )
   close( geneio )
end


text:
@testset "GaussianBlur" begin
    @testset "constructor" begin
        @test_throws ArgumentError GaussianBlur(0)
        @test_throws ArgumentError GaussianBlur(3, 0)
        @test_throws ArgumentError GaussianBlur([1, 3], 0)
        @test_throws ArgumentError GaussianBlur([1, 0], 1)
        @test_throws ArgumentError GaussianBlur(3, -1:1)
    end
    @testset "randparam" begin
        img = testpattern()
        ks = [3, 1:2:5, [1, 3, 5]]
        σs = [1, 1:0.1:2, [1, 2]]
        for k in ks, σ in σs
            op = GaussianBlur(k, σ)
            p = @inferred Augmentor.randparam(op, img)
            @test p[1] in k
            @test p[2] in σ
        end
    end
    @testset "eager" begin
        @test Augmentor.supports_eager(GaussianBlur)

        k = 3
        σ = 3.2

        imgs = [testpattern(), camera]

        for img in imgs
            img = convert.(RGB, img)
            ref = imfilter(img, KernelFactors.gaussian((σ, σ), (k, k)))
            res = Augmentor.applyeager(GaussianBlur(k, σ), img)
            @test assess_psnr(res, ref) > 25
        end
    end
end