text:
module Gnuplot

using StructC14N, ColorTypes, Printf, StatsBase, ReusePatterns, DataFrames
using ColorSchemes

import Base.reset
import Base.write
import Base.iterate
import Base.convert
import Base.string

export @gp, @gsp, save, linestyles, palette, contourlines, hist

# ╭───────────────────────────────────────────────────────────────────╮
# │                           TYPE DEFINITIONS                        │
# ╰───────────────────────────────────────────────────────────────────╯
# ---------------------------------------------------------------------
mutable struct DataSet
    name::String
    lines::Vector{String}
end

# ---------------------------------------------------------------------
mutable struct SinglePlot
    cmds::Vector{String}
    elems::Vector{String}
    flag3d::Bool
    SinglePlot() = new(Vector{String}(), Vector{String}(), false)
end

# ---------------------------------------------------------------------
@quasiabstract mutable struct DrySession
    sid::Symbol                # session ID
    datas::Vector{DataSet}     # data sets
    plots::Vector{SinglePlot}  # commands and plot commands (one entry for each plot of the multiplot)
    curmid::Int                # current multiplot ID
end

# ---------------------------------------------------------------------
@quasiabstract mutable struct GPSession <: DrySession
    pin::Base.Pipe;
    pout::Base.Pipe;
    perr::Base.Pipe;
    proc::Base.Process;
    channel::Channel{String};
end

# ---------------------------------------------------------------------
Base.@kwdef mutable struct Options
    dry::Bool = false                         # Use "dry" sessions (i.e. without an underlying Gnuplot process)
    cmd::String = "gnuplot"                   # Customizable command to start the Gnuplot process
    default::Symbol = :default                # Default session name
    init::Vector{String} = Vector{String}()   # Commands to initialize the gnuplot session (e.g., to set default terminal)
    verbose::Bool = false                     # verbosity flag (true/false)
    datalines::Int = 4;                       # How many lines of a dataset are printed in log
end
const sessions = Dict{Symbol, DrySession}()
const options = Options()

# ╭───────────────────────────────────────────────────────────────────╮
# │                         LOW LEVEL FUNCTIONS                       │
# ╰───────────────────────────────────────────────────────────────────╯
# ---------------------------------------------------------------------
function string(c::ColorTypes.RGB)
    return string(float(c.r)*255) * " " * string(float(c.g)*255) * " " * string(float(c.b)*255)
end

# ---------------------------------------------------------------------
"""
  # CheckGnuplotVersion

  Check whether gnuplot is runnable with the command given in `cmd`.
  Also check that gnuplot version is >= 4.7 (required to use data
  blocks).
"""
function CheckGnuplotVersion(cmd::AbstractString)
    icmd = `$(cmd) --version`

    proc = open(`$icmd`, read=true)
    s = String(read(proc))
    if !success(proc)
        error("An error occurred while running: " * string(icmd))
    end

    s = split(s, " ")
    ver = ""
    for token in s
        try
            ver = VersionNumber("$token")
            break
        catch
        end
    end

    if ver < v"4.7"
        error("gnuplot ver. >= 4.7 is required, but " * string(ver) * " was found.")
    end
    #@info "  Gnuplot version: " * string(ver)
    return ver
end

# ---------------------------------------------------------------------
function parseKeywords(; kwargs...)
    template = (xrange=NTuple{2, Real},
                yrange=NTuple{2, Real},
                zrange=NTuple{2, Real},
                cbrange=NTuple{2, Real},
                title=AbstractString,
                xlabel=AbstractString,
                ylabel=AbstractString,
                zlabel=AbstractString,
                xlog=Bool,
                ylog=Bool,
                zlog=Bool)

    kw = canonicalize(template; kwargs...)
    out = Vector{String}()
    ismissing(kw.xrange ) || (push!(out, "set xrange  [" * join(kw.xrange , ":") * "]"))
    ismissing(kw.yrange ) || (push!(out, "set yrange  [" * join(kw.yrange , ":") * "]"))
    ismissing(kw.zrange ) || (push!(out, "set zrange  [" * join(kw.zrange , ":") * "]"))
    ismissing(kw.cbrange) || (push!(out, "set cbrange [" * join(kw.cbrange, ":") * "]"))
    ismissing(kw.title  ) || (push!(out, "set title  \"" * kw.title  * "\""))
    ismissing(kw.xlabel ) || (push!(out, "set xlabel \"" * kw.xlabel * "\""))
    ismissing(kw.ylabel ) || (push!(out, "set ylabel \"" * kw.ylabel * "\""))
    ismissing(kw.zlabel ) || (push!(out, "set zlabel \"" * kw.zlabel * "\""))
    ismissing(kw.xlog   ) || (push!(out, (kw.xlog  ?  ""  :  "un") * "set logscale x"))
    ismissing(kw.ylog   ) || (push!(out, (kw.ylog  ?  ""  :  "un") * "set logscale y"))
    ismissing(kw.zlog   ) || (push!(out, (kw.zlog  ?  ""  :  "un") * "set logscale z"))
    return out
end

# ---------------------------------------------------------------------
function data2string(args...)
    @assert length(args) > 0

    # Check types of args
    for iarg in 1:length(args)
        d = args[iarg]

        ok = false
        if typeof(d) <: Number
            ok = true
        elseif typeof(d) <: AbstractArray
            if typeof(d[1]) <: Number
                ok = true
            end
            if typeof(d[1]) <: ColorTypes.RGB
                ok = true
            end
        end
        @assert ok "Invalid argument type at position $iarg"
    end

    # Collect lengths and number of dims
    lengths = Vector{Int}()
    dims = Vector{Int}()
    firstMultiDim = 0
    for i in 1:length(args)
        d = args[i]
        @assert ndims(d) <= 3 "Array dimensions must be <= 3"
        push!(lengths, length(d))
        push!(dims   , ndims(d))
        (firstMultiDim == 0)  &&  (ndims(d) > 1)  &&  (firstMultiDim = i)
    end

    accum = Vector{String}()

    # All scalars
    if minimum(dims) == 0
        #@info "Case 0"
        @assert maximum(dims) == 0 "Input data are ambiguous: either use all scalar or arrays of floats"
        v = ""
        for iarg in 1:length(args)
            d = args[iarg]
            v *= " " * string(d)
        end
        push!(accum, v)
        return accum
    end

    @assert all((dims .== 1)  .|  (dims .== maximum(dims))) "Array size are incompatible"

    # All 1D
    if firstMultiDim == 0
        #@info "Case 1"
        @assert minimum(lengths) == maximum(lengths) "Array size are incompatible"
        for i in 1:lengths[1]
            v = ""
            for iarg in 1:length(args)
                d = args[iarg]
                v *= " " * string(d[i])
            end
            push!(accum, v)
        end
        return accum
    end

    # Multidimensional, no independent indices
    if firstMultiDim == 1
        #@info "Case 2"
        @assert minimum(lengths) == maximum(lengths) "Array size are incompatible"
        i = 1
        for CIndex in CartesianIndices(size(args[1]))
            indices = Tuple(CIndex)
            (i > 1)  &&  (indices[end-1] == 1)  &&  (push!(accum, ""))  # blank line
            v = "" # * join(string.(getindex.(Ref(Tuple(indices)), 1:ndims(args[1]))), " ")
            for iarg in 1:length(args)
                d = args[iarg]
                v *= " " * string(d[i])
            end
            i += 1
            push!(accum, v)
        end
        return accum
    end

    # Multidimensional (independent indices provided in input)
    if firstMultiDim >= 2
        @assert (firstMultiDim-1 == dims[firstMultiDim]) "Not enough independent variables"
        refLength = lengths[firstMultiDim]
        @assert all(lengths[firstMultiDim:end] .== refLength) "Array size are incompatible"

        if lengths[1] < refLength
            #@info "Case 3"
            # Cartesian product of Independent variables
            checkLength = prod(lengths[1:firstMultiDim-1])
            @assert prod(lengths[1:firstMultiDim-1]) == refLength "Array size are incompatible"

            i = 1
            for CIndex in CartesianIndices(size(args[firstMultiDim]))
                indices = Tuple(CIndex)
                (i > 1)  &&  (indices[end-1] == 1)  &&  (push!(accum, ""))  # blank line
                v = ""
                for iarg in 1:firstMultiDim-1
                    d = args[iarg]
                    v *= " " * string(d[indices[iarg]])
                end
                for iarg in firstMultiDim:length(args)
                    d = args[iarg]
                    v *= " " * string(d[i])
                end
                i += 1
                push!(accum, v)
            end
            return accum
        else
            #@info "Case 4"
            # All Independent variables have the same length as the main multidimensional data
            @assert all(lengths[1:firstMultiDim-1] .== refLength) "Array size are incompatible"

            i = 1
            for CIndex in CartesianIndices(size(args[firstMultiDim]))
                indices = Tuple(CIndex)
                (i > 1)  &&  (indices[end-1] == 1)  &&  (push!(accum, ""))  # blank line
                v = ""
                for iarg in 1:length(args)
                    d = args[iarg]
                    v *= " " * string(d[i])
                end
                i += 1
                push!(accum, v)
            end
            return accum
        end
    end

    return nothing
end

# ╭───────────────────────────────────────────────────────────────────╮
# │                SESSION CONSTRUCTORS AND getsession()              │
# ╰───────────────────────────────────────────────────────────────────╯
# ---------------------------------------------------------------------
function DrySession(sid::Symbol)
    global options
    (sid in keys(sessions))  &&  error("Gnuplot session $sid is already active")
    out = DrySession(sid, Vector{DataSet}(), [SinglePlot()], 1)
    sessions[sid] = out
    return out
end

# ---------------------------------------------------------------------
function GPSession(sid::Symbol)
    function readTask(sid, stream, channel)
        global options
        saveOutput = false

        while isopen(stream)
            line = readline(stream)
            if (length(line) >= 1)  &&  (line[1] == Char(0x1b)) # Escape (xterm -ti vt340)
                buf = Vector{UInt8}()
                append!(buf, convert(Vector{UInt8}, [line...]))
                push!(buf, 0x0a)
                c = 0x00
                while c != 0x1b
                    c = read(stream, 1)[1]
                    push!(buf, c)
                end
                c = read(stream, 1)[1]
                push!(buf, c)
                write(stdout, buf)
                continue
            end
            if line == "GNUPLOT_CAPTURE_BEGIN"
                saveOutput = true
            else
                if (line != "")  &&  (line != "GNUPLOT_CAPTURE_END")  &&  (options.verbose)
                    printstyled(color=:cyan, "GNUPLOT ($sid) -> $line\n")
                end
                (saveOutput)  &&  (put!(channel, line))
                (line == "GNUPLOT_CAPTURE_END")  &&  (saveOutput = false)
            end
        end
        delete!(sessions, sid)
        return nothing
    end

    global options

    CheckGnuplotVersion(options.cmd)
    session = DrySession(sid)

    pin  = Base.Pipe()
    pout = Base.Pipe()
    perr = Base.Pipe()
    proc = run(pipeline(`$(options.cmd)`, stdin=pin, stdout=pout, stderr=perr), wait=false)
    chan = Channel{String}(32)

    # Close unused sides of the pipes
    Base.close(pout.in)
    Base.close(perr.in)
    Base.close(pin.out)
    Base.start_reading(pout.out)
    Base.start_reading(perr.out)

    # Start reading tasks
    @async readTask(sid, pout, chan)
    @async readTask(sid, perr, chan)

    out = GPSession(getfield.(Ref(session), fieldnames(concretetype(DrySession)))...,
                    pin, pout, perr, proc, chan)
    sessions[sid] = out

    # Set window title
    term = writeread(out, "print GPVAL_TERM")[1]
    if term in ("aqua", "x11", "qt", "wxt")
        opts = writeread(out, "print GPVAL_TERMOPTIONS")[1]
        if findfirst("title", opts) == nothing
            writeread(out, "set term $term $opts title 'Gnuplot.jl: $(out.sid)'")
        end
    end
    for l in options.init
        writeread(out, l)
    end

    return out
end

# ---------------------------------------------------------------------
function getsession(sid::Symbol=options.default)
    global options
    if !(sid in keys(sessions))
        if options.dry
            DrySession(sid)
        else
            GPSession(sid)
        end
    end
    return sessions[sid]
end

# ╭───────────────────────────────────────────────────────────────────╮
# │                       write() and writeread()                     │
# ╰───────────────────────────────────────────────────────────────────╯
# ---------------------------------------------------------------------
"""
  # write

  Send a string to gnuplot's STDIN.

  The commands sent through `write` are not stored in the current
  session (use `newcmd` to save commands in the current session).

  ## Arguments:
  - `gp`: a `DrySession` object;
  - `str::String`: command to be sent;
"""
write(gp::DrySession, str::AbstractString) = nothing
function write(gp::GPSession, str::AbstractString)
    global options
    if options.verbose
        printstyled(color=:light_yellow, "GNUPLOT ($(gp.sid)) $str\n")
    end
    w = write(gp.pin, strip(str) * "\n")
    w <= 0  &&  error("Writing on gnuplot STDIN pipe returned $w")
    flush(gp.pin)
    return w
end

write(gp::DrySession, d::DataSet) = nothing
function write(gp::GPSession, d::DataSet)
    if options.verbose
        v = ""
        printstyled(color=:light_black, "GNUPLOT ($(gp.sid)) $(d.name) << EOD\n")
        n = min(options.datalines, length(d.lines))
        for i in 1:n
            printstyled(color=:light_black, "GNUPLOT ($(gp.sid)) $(d.lines[i])\n")
        end
        if n < length(d.lines)
            printstyled(color=:light_black, "GNUPLOT ($(gp.sid)) ...\n")
        end
        printstyled(color=:light_black, "GNUPLOT ($(gp.sid)) EOD\n")
    end
    write(gp.pin, "$(d.name) << EOD\n")
    write(gp.pin, join(d.lines, "\n") * "\n")
    write(gp.pin, "EOD\n")
    flush(gp.pin)
    return nothing
end

# ---------------------------------------------------------------------
writeread(gp::DrySession, str::AbstractString) = [""]
function writeread(gp::GPSession, str::AbstractString)
    global options
    verbose = options.verbose

    options.verbose = false
    write(gp, "print 'GNUPLOT_CAPTURE_BEGIN'")

    options.verbose = verbose
    write(gp, str)

    options.verbose = false
    write(gp, "print 'GNUPLOT_CAPTURE_END'")
    options.verbose = verbose

    out = Vector{String}()
    while true
        l = take!(gp.channel)
        l == "GNUPLOT_CAPTURE_END"  &&  break
        push!(out, l)
    end
    return out
end

# ╭───────────────────────────────────────────────────────────────────╮
# │              PRIVATE FUNCTIONS TO MANIPULATE SESSIONS             │
# ╰───────────────────────────────────────────────────────────────────╯
# ---------------------------------------------------------------------
function reset(gp::DrySession)
    gp.datas = Vector{DataSet}()
    gp.plots = [SinglePlot()]
    gp.curmid = 1
    exec(gp, "reset session")
    return nothing
end

# ---------------------------------------------------------------------
function setmulti(gp::DrySession, mid::Int)
    @assert mid >= 0 "Multiplot ID must be a >= 0"
    for i in length(gp.plots)+1:mid
        push!(gp.plots, SinglePlot())
    end
    (mid > 0)  &&  (gp.curmid = mid)
end

# ---------------------------------------------------------------------
function newdataset(gp::DrySession, accum::Vector{String}; name="")
    (name == "")  &&  (name = string("data", length(gp.datas)))
    name = "\$$name"
    d = DataSet(name, accum)
    push!(gp.datas, d)
    write(gp, d) # Send now to gnuplot process
    return name
end
newdataset(gp::DrySession, args...; name="") = newdataset(gp, data2string(args...), name=name)

# ---------------------------------------------------------------------
function newcmd(gp::DrySession, v::String; mid::Int=0)
    setmulti(gp, mid)
    (v != "")  &&  (push!(gp.plots[gp.curmid].cmds, v))
    (length(gp.plots) == 1)  &&  (exec(gp, v))  # execute now to check against errors
    return nothing
end

function newcmd(gp::DrySession; mid::Int=0, args...)
    for v in parseKeywords(;args...)
        newcmd(gp, v, mid=mid)
    end
    return nothing
end

# ---------------------------------------------------------------------
function newplot(gp::DrySession, name, opt=""; mid=0)
    setmulti(gp, mid)
    push!(gp.plots[gp.curmid].elems, "$name $opt")
end

# ---------------------------------------------------------------------
function quit(gp::DrySession)
    global options
    delete!(sessions, gp.sid)
    return 0
end

function quit(gp::GPSession)
    close(gp.pin)
    close(gp.pout)
    close(gp.perr)
    wait( gp.proc)
    exitCode = gp.proc.exitcode
    invoke(quit, Tuple{DrySession}, gp)
    return exitCode
end

# ╭───────────────────────────────────────────────────────────────────╮
# │                 execall(), dump() and driver()                    │
# ╰───────────────────────────────────────────────────────────────────╯
# ---------------------------------------------------------------------
execall(gp::DrySession; term::AbstractString="", output::AbstractString="") = nothing
function execall(gp::GPSession; term::AbstractString="", output::AbstractString="")
    if term != ""
        former_term = writeread(gp, "print GPVAL_TERM")[1]
        former_opts = writeread(gp, "print GPVAL_TERMOPTIONS")[1]
        exec(gp, "set term $term")
    end
    (output != "")  &&  exec(gp, "set output '$output'")

    for i in 1:length(gp.plots)
        d = gp.plots[i]
        for j in 1:length(d.cmds)
            exec(gp, d.cmds[j])
        end
        if length(d.elems) > 0
            s = (d.flag3d  ?  "splot "  :  "plot ") * " \\\n  " *
                join(d.elems, ", \\\n  ")
            exec(gp, s)
        end
    end
    (length(gp.plots) > 1)  &&  exec(gp, "unset multiplot")
    (output != "")  &&  exec(gp, "set output")
    if term != ""
        exec(gp, "set term $former_term $former_opts")
    end
    return nothing
end

function savescript(gp::DrySession, filename; term::AbstractString="", output::AbstractString="")
    stream = open(filename, "w")

    println(stream, "reset session")
    if term != ""
        println(stream, "set term $term")
    end
    (output != "")  &&  println(stream, "set output '$output'")

    for i in 1:length(gp.datas)
        d = gp.datas[i]
        println(stream, d.name * " << EOD")
        for j in 1:length(d.lines)
            println(stream, d.lines[j])
        end
        println(stream, "EOD")
    end

    for i in 1:length(gp.plots)
        d = gp.plots[i]
        for j in 1:length(d.cmds)
            println(stream, d.cmds[j])
        end
        if length(d.elems) > 0
            s = (d.flag3d  ?  "splot "  :  "plot ") * " \\\n  " *
                join(d.elems, ", \\\n  ")
            println(stream, s)
        end
    end
    (length(gp.plots) > 1)  &&  println(stream, "unset multiplot")
    println(stream, "set output")
    close(stream)
    return nothing
end

# ---------------------------------------------------------------------
function driver(args...; flag3d=false)
    if length(args) == 0
        gp = getsession()
        execall(gp)
        return nothing
    end

    data = Vector{Any}()
    dataname = ""
    dataplot = nothing

    function dataCompleted()
        if length(data) > 0
            AllArraysAreNotEmpty = true
            for i in 1:length(data)
                if (typeof(data[i]) <: AbstractArray)  &&  (length(data[i]) == 0)
                    @warn "Input array is empty"
                    AllArraysAreNotEmpty = false
                    break
                end
            end
            if AllArraysAreNotEmpty
                last = newdataset(gp, data...; name=dataname)
                (dataplot != nothing)  &&  (newplot(gp, last, dataplot))
            end
        end
        data = Vector{Any}()
        dataname = ""
        dataplot = nothing
    end
    function isPlotCmd(s::String)
        (length(s) >= 2)  &&  (s[1:2] ==  "p "    )  &&  (return (true, false, strip(s[2:end])))
        (length(s) >= 3)  &&  (s[1:3] ==  "pl "   )  &&  (return (true, false, strip(s[3:end])))
        (length(s) >= 4)  &&  (s[1:4] ==  "plo "  )  &&  (return (true, false, strip(s[4:end])))
        (length(s) >= 5)  &&  (s[1:5] ==  "plot " )  &&  (return (true, false, strip(s[5:end])))
        (length(s) >= 2)  &&  (s[1:2] ==  "s "    )  &&  (return (true, true , strip(s[2:end])))
        (length(s) >= 3)  &&  (s[1:3] ==  "sp "   )  &&  (return (true, true , strip(s[3:end])))
        (length(s) >= 4)  &&  (s[1:4] ==  "spl "  )  &&  (return (true, true , strip(s[4:end])))
        (length(s) >= 5)  &&  (s[1:5] ==  "splo " )  &&  (return (true, true , strip(s[5:end])))
        (length(s) >= 6)  &&  (s[1:6] ==  "splot ")  &&  (return (true, true , strip(s[6:end])))
        return (false, false, "")
    end

    gp = nothing
    doDump  = true
    doReset = true

    for loop in 1:2
        if loop == 2
            (gp == nothing)  &&  (gp = getsession())
            doReset  &&  reset(gp)
            gp.plots[gp.curmid].flag3d = flag3d
        end

        for iarg in 1:length(args)
            arg = args[iarg]

            if typeof(arg) == Symbol
                if arg == :-
                    (loop == 1)  &&  (iarg < length(args)) &&  (doReset = false)
                    (loop == 1)  &&  (iarg > 1 )           &&  (doDump  = false)
                else
                    (loop == 1)  &&  (gp = getsession(arg))
                end
            elseif isa(arg, Tuple)  &&  length(arg) == 2  &&  isa(arg[1], Symbol)
                if arg[1] == :term
                    if loop == 1
                        if typeof(arg[2]) == String
                            term = (deepcopy(arg[2]), "")
                        elseif length(arg[2]) == 2
                            term = deepcopy(arg[2])
                        else
                            error("The term tuple must contain at most two strings")
                        end
                    end
                else
                    (loop == 2)  &&  newcmd(gp; [arg]...) # A cmd keyword
                end
            elseif isa(arg, Int)
                (loop == 2)  &&  (@assert arg > 0)
                (loop == 2)  &&  (dataplot = ""; dataCompleted())
                (loop == 2)  &&  setmulti(gp, arg)
                (loop == 2)  &&  (gp.plots[gp.curmid].flag3d = flag3d)
            elseif isa(arg, String)
                # Either a dataname, a plot or a command
                if loop == 2
                    if isa(arg, String)  &&  (length(arg) > 1)  &&  (arg[1] == '$')
                        dataname = arg[2:end]
                        dataCompleted()
                    elseif length(data) > 0
                        dataplot = arg
                        dataCompleted()
                    else
                        (isPlot, is3d, cmd) = isPlotCmd(arg)
                        if isPlot
                            gp.plots[gp.curmid].flag3d = is3d
                            newplot(gp, cmd)
                        else
                            newcmd(gp, arg)
                        end
                    end
                end
            else
                (loop == 2)  &&  push!(data, arg) # a data set
            end
        end
    end

    dataplot = ""
    dataCompleted()
    (doDump)  &&  (execall(gp))

    return nothing
end

#_____________________________________________________________________
#                         EXPORTED FUNCTIONS
#_____________________________________________________________________

# --------------------------------------------------------------------
"""
`@gp args...`

The `@gp` macro (and its companion `@gsp`, for `splot` operations) allows to exploit all of the **Gnuplot** package functionalities using an extremely efficient and concise syntax.  Both macros accept the same syntax, as described below.

The macros accepts any number of arguments, with the following meaning:
- a symbol: the name of the session to use;
- a string: a command (e.g. "set key left") or plot specification (e.g. "with lines");
- a string starting with a `\$` sign: a data set name;
- an `Int` > 0: the plot destination in a multiplot session;
- a keyword/value pair: a keyword value (see below);
- any other type: a dataset to be passed to Gnuplot.  Each dataset must be terminated by either:
  - a string starting with a `\$` sign (i.e. the data set name);
  - or a string with the plot specifications (e.g. "with lines");
- the `:-` symbol, used as first argument, avoids resetting the Gnuplot session.  Used as last argument avoids immediate execution  of the plot/splot command.  This symbol can be used to split a  single call into multiple ones.

All entries are optional, and there is no mandatory order.  The plot specification can either be:
 - a complete plot/splot command (e.g., "plot sin(x)", both "plot" and "splot" can be abbreviated to "p" and "s" respectively);
 - or a partial specification starting with the "with" clause (if it follows a data set).

The list of accepted keyword is as follows:
- `title::String`: plot title;
- `xlabel::String`: X axis label;
- `ylabel::String`: Y axis label;
- `zlabel::String`: Z axis label;
- `xlog::Bool`: logarithmic scale for X axis;
- `ylog::Bool`: logarithmic scale for Y axis;
- `zlog::Bool`: logarithmic scale for Z axis;
- `xrange::NTuple{2, Number}`: X axis range;
- `yrange::NTuple{2, Number}`: Y axis range;
- `zrange::NTuple{2, Number}`: Z axis range;
- `cbrange::NTuple{2, Number}`: Color box axis range;

The symbol for the above-mentioned keywords may also be used in a shortened form, as long as there is no ambiguity with other keywords.  E.g. you can use: `xr=(1,10)` in place of `xrange=(1,10)`.

# Examples:

## Simple examples with no data:
```
@gp "plot sin(x)"
@gp "plot sin(x)" "pl cos(x)"
@gp "plo sin(x)" "s cos(x)"

# Split a `@gp` call in two
@gp "plot sin(x)" :-
@gp :- "plot cos(x)"

# Insert a 3 second pause between one plot and the next
@gp "plot sin(x)" 2 xr=(-2pi,2pi) "pause 3" "plot cos(4*x)"
```

### Simple examples with data:
```
@gp "set key left" tit="My title" xr=(1,12) 1:10 "with lines tit 'Data'"

x = collect(1.:10)
@gp x
@gp x x
@gp x -x
@gp x x.^2
@gp x x.^2 "w l"

lw = 3
@gp x x.^2 "w l lw \$lw"
```

### A more complex example
```
@gp("set grid", "set key left", xlog=true, ylog=true,
    title="My title", xlab="X label", ylab="Y label",
    x, x.^0.5, "w l tit 'Pow 0.5' dt 2 lw 2 lc rgb 'red'",
    x, x     , "w l tit 'Pow 1'   dt 1 lw 3 lc rgb 'blue'",
    x, x.^2  , "w l tit 'Pow 2'   dt 3 lw 2 lc rgb 'purple'")
```

### Multiplot example:
```
@gp(xr=(-2pi,2pi), "unset key",
    "set multi layout 2,2 title 'Multiplot title'",
    1, "p sin(x)"  ,
    2, "p sin(2*x)",
    3, "p sin(3*x)",
    4, "p sin(4*x)")
```
or equivalently
```
@gp xr=(-2pi,2pi) "unset key" "set multi layout 2,2 title 'Multiplot title'" :-
for i in 1:4
  @gp :- i "p sin(\$i*x)" :-
end
@gp
```

### Multiple gnuplot sessions
```
@gp :GP1 "plot sin(x)"
@gp :GP2 "plot sin(x)"

Gnuplot.quitall()
```

### Further examples
```
x = range(-2pi, stop=2pi, length=100);
y = 1.5 * sin.(0.3 .+ 0.7x) ;
noise = randn(length(x))./2;
e = 0.5 * fill(1, size(x));

name = "\\\$MyDataSet1"
@gp x y name "plot \$name w l" "pl \$name u 1:(2*\\\$2) w l"

@gsp randn(Float64, 30, 50)
@gp randn(Float64, 30, 50) "w image"
@gsp x y y

@gp("set key horizontal", "set grid",
    xrange=(-7,7), ylabel="Y label",
    x, y, "w l t 'Real model' dt 2 lw 2 lc rgb 'red'",
    x, y+noise, e, "w errorbars t 'Data'")

@gp "f(x) = a * sin(b + c*x); a = 1; b = 1; c = 1;"   :-
@gp :- x y+noise e name                               :-
@gp :- "fit f(x) \$name u 1:2:3 via a, b, c;"         :-
@gp :- "set multiplot layout 2,1"                     :-
@gp :- "plot \$name w points" ylab="Data and model"   :-
@gp :- "plot \$name u 1:(f(\\\$1)) w lines"           :-
@gp :- 2 xlab="X label" ylab="Residuals"              :-
@gp :- "plot \$name u 1:((f(\\\$1)-\\\$2) / \\\$3):(1) w errorbars notit"

# Retrieve values for a, b and c
a = Meta.parse(Gnuplot.exec("print a"))
b = Meta.parse(Gnuplot.exec("print b"))
c = Meta.parse(Gnuplot.exec("print c"))

# Save to a PDF file
save(term="pdf", output="gnuplot.pdf")
```

### Display an image
```
using TestImages
img = testimage("lena");
@gp img "w image"
@gp "set size square" img "w rgbimage" # Color image with correct proportions
@gp "set size square" img "u 2:(-\\\$1):3:4:5 with rgbimage" # Correct orientation
```
"""
macro gp(args...)
    out = Expr(:call)
    push!(out.args, :(Gnuplot.driver))
    for iarg in 1:length(args)
        arg = args[iarg]
        if (isa(arg, Expr)  &&  (arg.head == :(=)))
            sym = string(arg.args[1])
            val = arg.args[2]
            push!(out.args, :((Symbol($sym),$val)))
        else
            push!(out.args, arg)
        end
    end
    return esc(out)
end

"""
  # @gsp

  See documentation for `@gp`.
"""
macro gsp(args...)
    out = Expr(:macrocall, Symbol("@gp"), LineNumberNode(1, "Gnuplot.jl"))
    push!(out.args, args...)
    push!(out.args, Expr(:kw, :flag3d, true))
    return esc(out)
end

# ╭───────────────────────────────────────────────────────────────────╮
# │              FUNCTIONS MEANT TO BE INVOKED BY USERS               │
# ╰───────────────────────────────────────────────────────────────────╯
# ---------------------------------------------------------------------
"""
  `quit()`

  Quit the session and the associated gnuplot process (if any).
"""
function quit(sid::Symbol)
    global options
    if !(sid in keys(sessions))
        error("Gnuplot session $sid do not exists")
    end
    return quit(sessions[sid])
end

"""
  `quitall()`

  Quit all the sessions and the associated gnuplot processes.
"""
function quitall()
    global options
    for sid in keys(sessions)
        quit(sid)
    end
    return nothing
end

# --------------------------------------------------------------------
"""
`exec(sid::Symbol, s::Vector{String})`

Directly execute commands on the underlying Gnuplot process, and return the result(s).

## Examples:
```julia
exec("print GPVAL_TERM")
exec("plot sin(x)")
```
"""
exec(gp::DrySession, command::String) = nothing
function exec(gp::GPSession, command::String)
    answer = Vector{String}()
    push!(answer, writeread(gp, command)...)

    verbose = options.verbose
    options.verbose = false
    errno = writeread(gp, "print GPVAL_ERRNO")[1]
    options.verbose = verbose

    if errno != "0"
        printstyled(color=:red, "GNUPLOT ERROR $(gp.sid) -> ERRNO=$errno\n")
        errmsg = writeread(gp, "print GPVAL_ERRMSG")
        write(gp.pin, "reset error\n")
        for line in errmsg
            printstyled(color=:red, "GNUPLOT ERROR $(gp.sid) -> $line\n")
        end
        error("Gnuplot process raised an error")
    end

    return join(answer, "\n")
end
function exec(s::String)
    global options
    exec(getsession(), s)
end
exec(sid::Symbol, s::String) = exec(getsession(sid), s)

# --------------------------------------------------------------------
"""
`setverbose(b::Bool)`

Set verbose flag to `true` or `false` (default: `false`).
"""
function setverbose(b::Bool)
    global options
    options.verbose = b
end

# --------------------------------------------------------------------
"""
`save(...)`

Save the data and commands in the current session to either:
- the gnuplot process (i.e. produce a plot): `save(term="", output="")`;
- an IO stream: `save(stream::IO; term="", output="")`;
- a file: `save(file::AbstractStrings; term="", output="")`.

To save the data and command from a specific session pass the ID as first argument, i.e.:
- `save(sid::Symbol, term="", output="")`;
- `save(sid::Symbol, file::AbstractStrings; term="", output="")`.

In all cases the `term` keyword allows to specify a gnuplot terminal, and the `output` keyword allows to specify an output file.
"""
save(           ; kw...) = execall(getsession()   ; kw...)
save(sid::Symbol; kw...) = execall(getsession(sid); kw...)
save(             file::AbstractString; kw...) = savescript(getsession()   , file, kw...)
save(sid::Symbol, file::AbstractString; kw...) = savescript(getsession(sid), file, kw...)

# ╭───────────────────────────────────────────────────────────────────╮
# │                     HIGH LEVEL FACILITIES                         │
# ╰───────────────────────────────────────────────────────────────────╯
# ---------------------------------------------------------------------
linestyles(s::Symbol) = linestyles(colorschemes[s])
function linestyles(cmap::ColorScheme)
    styles = Vector{String}()
    for i in 1:length(cmap.colors)
        push!(styles, "set style line $i lt 1 lc rgb '#" * Base.hex(cmap.colors[i]))
    end
    return join(styles, "\n")
end

# --------------------------------------------------------------------
palette(s::Symbol) = palette(colorschemes[s])
function palette(cmap::ColorScheme)
    levels = Vector{String}()
    for x in LinRange(0, 1, length(cmap.colors))
        color = get(cmap, x)
        push!(levels, "$x '#" * Base.hex(color) * "'")
    end
    return "set palette defined (" * join(levels, ", ") * ")\nset palette maxcol $(length(cmap.colors))\n"
end

# ╭───────────────────────────────────────────────────────────────────╮
# │                     EXPERIMENTAL FUNCTIONS                        │
# ╰───────────────────────────────────────────────────────────────────╯
# # --------------------------------------------------------------------
# """
#   # repl
#
#   Read/evaluate/print/loop
# """
# function repl(sid::Symbol)
#     verb = options.verbose
#     options.verbose = 0
#     gp = getsession(sid)
#     while true
#         line = readline(stdin)
#         (line == "")  &&  break
#         answer = send(gp, line, true)
#         for line in answer
#             println(line)
#         end
#     end
#     options.verbose = verb
#     return nothing
# end
# function repl()
#     global options
#     return repl(options.default)
# end

# --------------------------------------------------------------------
#=
Example:
v = randn(1000)
h = hist(v, bs=0.2)
@gp h.loc h.counts "w histep" h.loc h.counts "w l"
=#
function hist(v::Vector{T}; range=[NaN,NaN], bs=NaN, nbins=0, pad=true) where T <: Number
    i = findall(isfinite.(v))
    isnan(range[1])  &&  (range[1] = minimum(v[i]))
    isnan(range[2])  &&  (range[2] = maximum(v[i]))
    i = findall(isfinite.(v)  .&  (v.>= range[1])  .&  (v.<= range[2]))
    (nbins > 0)  &&  (bs = (range[2] - range[1]) / nbins)
    if isfinite(bs)
        rr = range[1]:bs:range[2]
        if maximum(rr) < range[2]
            rr = range[1]:bs:(range[2]+bs)
        end
        hh = fit(Histogram, v[i], rr, closed=:left)
        if sum(hh.weights) < length(i)
            j = findall(v[i] .== range[2])
            @assert length(j) == (length(i) - sum(hh.weights))
            hh.weights[end] += length(j)
        end
    else
        hh = fit(Histogram, v[i], closed=:left)
    end
    @assert sum(hh.weights) == length(i)
    x = collect(hh.edges[1])
    x = (x[1:end-1] .+ x[2:end]) ./ 2
    h = hh.weights
    binsize = x[2] - x[1]
    if pad
        x = [x[1]-binsize, x..., x[end]+binsize]
        h = [0, h..., 0]
    end
    return (loc=x, counts=h, binsize=binsize)
end

# --------------------------------------------------------------------
function hist(v1::Vector{T1}, v2::Vector{T2};
              range1=[NaN,NaN], bs1=NaN, nbins1=0,
              range2=[NaN,NaN], bs2=NaN, nbins2=0) where {T1 <: Number, T2 <: Number}
    i = findall(isfinite.(v2))
    isnan(range1[1])  &&  (range1[1] = minimum(v1[i]))
    isnan(range1[2])  &&  (range1[2] = maximum(v1[i]))
    i = findall(isfinite.(v2))
    isnan(range2[1])  &&  (range2[1] = minimum(v2[i]))
    isnan(range2[2])  &&  (range2[2] = maximum(v2[i]))

    i1 = findall(isfinite.(v1)  .&  (v1.>= range1[1])  .&  (v1.<= range1[2]))
    i2 = findall(isfinite.(v2)  .&  (v2.>= range2[1])  .&  (v2.<= range2[2]))
    (nbins1 > 0)  &&  (bs1 = (range1[2] - range1[1]) / nbins1)
    (nbins2 > 0)  &&  (bs2 = (range2[2] - range2[1]) / nbins2)
    if isfinite(bs1) &&  isfinite(bs2)
        hh = fit(Histogram, (v1[i1], v2[i2]), (range1[1]:bs1:range1[2], range2[1]:bs2:range2[2]), closed=:left)
    else
        hh = fit(Histogram, (v1[i1], v2[i2]), closed=:left)
    end
    x1 = collect(hh.edges[1])
    x1 = (x1[1:end-1] .+ x1[2:end]) ./ 2
    x2 = collect(hh.edges[2])
    x2 = (x2[1:end-1] .+ x2[2:end]) ./ 2

    binsize1 = x1[2] - x1[1]
    binsize2 = x2[2] - x2[1]
    return (loc1=x1, loc2=x2, counts=hh.weights, binsize1=binsize1)
end

# --------------------------------------------------------------------
function contourlines(args...; cntrparam="level auto 10")
    tmpfile = Base.Filesystem.tempname()
    sid = Symbol("j", Base.Libc.getpid())
    if !haskey(Gnuplot.sessions, sid)
        gp = getsession(sid)
    end

    Gnuplot.exec(sid, "set term unknown")
    @gsp    sid "set contour base" "unset surface" :-
    @gsp :- sid "set cntrparam $cntrparam" :-
    @gsp :- sid "set table '$tmpfile'" :-
    @gsp :- sid args...
    Gnuplot.exec(sid, "unset table")
    Gnuplot.exec(sid, "reset")

    out = DataFrame()
    curlevel = NaN
    curx = Vector{Float64}()
    cury = Vector{Float64}()
    curid = 1
    elength(x, y) = sqrt.((x[2:end] .- x[1:end-1]).^2 .+
                          (y[2:end] .- y[1:end-1]).^2)
    function dump()
        ((length(curx) < 2)  ||  isnan(curlevel))  &&  return nothing
        tmp = DataFrame([Int, Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}],
                        [:id, :level , :len           , :x             , :y])
        push!(tmp, (curid, curlevel, [0.; elength(curx, cury)], [curx...], [cury...]))
        append!(out, tmp)
        curid += 1
        # d = cumsum(elength(curx, cury))
        # i0 = findall(d .<= width); sort!(i0)
        # i1 = findall(d .>  width)
        # if (length(i0) > 0)  &&  (length(i1) > 0)
        #     rot1 = atan(cury[i0[end]]-cury[i0[1]], curx[i0[end]]-curx[i0[1]]) * 180 / pi
        #     rot = round(mod(rot1, 360))
        #     x = mean(curx[i0])
        #     y = mean(cury[i0])
        #     push!(outl, "set label " * string(length(outl)+1) * " '$curlevel' at $x, $y center front rotate by $rot")
        #     curx = curx[i1]
        #     cury = cury[i1]
        # end
        empty!(curx)
        empty!(cury)
    end

    for l in readlines(tmpfile)
        if length(strip(l)) == 0
            dump()
            continue
        end
        if !isnothing(findfirst("# Contour ", l))
            dump()
            curlevel = Meta.parse(strip(split(l, ':')[2]))
            continue
        end
        (l[1] == '#')  &&  continue

        n = Meta.parse.(split(l))
        @assert length(n) == 3
        push!(curx, n[1])
        push!(cury, n[2])
    end
    rm(tmpfile)

    if nrow(out) > 0
        levels = unique(out.level)
        sort!(levels)
        out[!, :levelcount] .= 0
        for i in 1:length(levels)
            j = findall(out.level .== levels[i])
            out[j, :levelcount] .= i
        end
    end
    return out
end

# --------------------------------------------------------------------
function boxxyerror(x, y; xmin=NaN, ymin=NaN, xmax=NaN, ymax=NaN, cartesian=false)
    @assert length(x) == length(y)
    @assert issorted(x)
    @assert issorted(y)
    xlow  = Vector{Float64}(undef, length(x))
    xhigh = Vector{Float64}(undef, length(x))
    ylow  = Vector{Float64}(undef, length(x))
    yhigh = Vector{Float64}(undef, length(x))
    for i in 2:length(x)-1
        xlow[i]  = (x[i-1] + x[i]) / 2
        ylow[i]  = (y[i-1] + y[i]) / 2
        xhigh[i] = (x[i+1] + x[i]) / 2
        yhigh[i] = (y[i+1] + y[i]) / 2
    end
    xlow[1]    = (isfinite(xmin)  ?  xmin  :  (x[1] - (x[2]-x[1])/2))
    ylow[1]    = (isfinite(ymin)  ?  ymin  :  (y[1] - (y[2]-y[1])/2))
    xlow[end]  = (x[end] - (x[end]-x[end-1])/2)
    ylow[end]  = (y[end] - (y[end]-y[end-1])/2)
    xhigh[1]   = (x[1] + (x[2]-x[1])/2)
    yhigh[1]   = (y[1] + (y[2]-y[1])/2)
    xhigh[end] = (isfinite(xmax)  ?  xmax  :  (x[end] + (x[end]-x[end-1])/2))
    yhigh[end] = (isfinite(ymax)  ?  ymax  :  (y[end] + (y[end]-y[end-1])/2))
    if !cartesian
        return (x, y, xlow, xhigh, ylow, yhigh)
    end
    n = length(x)
    i = repeat(1:n, outer=n)
    j = repeat(1:n, inner=n)
    return (x[i], y[j], xlow[i], xhigh[i], ylow[j], yhigh[j])
end

# --------------------------------------------------------------------
function histo2segments(in_x, counts)
    @assert length(in_x) == length(counts)
    x = Vector{Float64}()
    y = Vector{Float64}()
    push!(x, in_x[1])
    push!(y, counts[1])
    for i in 2:length(in_x)
        xx = (in_x[i-1] + in_x[i]) / 2.
        push!(x, xx)
        push!(y, counts[i-1])
        push!(x, xx)
        push!(y, counts[i])
    end
    push!(x, in_x[end])
    push!(y, counts[end])
    return (x, y)
end

end #module


text:
using TropicalYao.Reversible
using TropicalNumbers
using LinearAlgebra, Test
using Yao

@testset "i-bond tensor" begin
    mat = ones(Tropical{Float64}, 2, 2)
    J = 0.3
    @test i_bondtensor!(mat, J)[1] ≈ i_bondtensor(J)
    mat = Diagonal(ones(Tropical{Float64}, 4))
    @test i_gvb_tensor!(mat, J)[1] ≈ i_gvb_tensor(J)
    mat = ones(Tropical{Float64}, 16, 16)
    Js = randn(16)
    @test i_g16_tensor!(mat, Js)[1] ≈ i_g16_tensor(Js)
end

@testset "tropical block" begin
    reg = ArrayReg(Tropical.(randn(1<<10)))
    S = stack4reg(reg, 10)
    reg1, _, _, S = apply_Ghb!(copy(reg), 2, 0.5, copy(S))
    reg2 = copy(reg) |> put(10, 2=>matblock(i_bondtensor(0.5)))
    @test reg1 ≈ reg2
    @test check_inv(apply_Ghb!, (copy(reg), 2, 0.5, S))
    reg1, _, _, S = apply_Gvb!(copy(reg), (2,4), 0.5, S)
    reg2 = copy(reg) |> put(10, (2,4)=>matblock(i_gvb_tensor(0.5)))
    @test reg1 ≈ reg2
    @test check_inv(apply_Gvb!, (copy(reg), (2,4), 0.5, S))
    Js = randn(16)
    reg1, _, _, S = apply_G16!(copy(reg), (1,2,4,9), Js, S)
    reg2 = copy(reg) |> put(10, (1,2,4,9)=>matblock(i_g16_tensor(Js)))
    @test reg1 ≈ reg2
    @test check_inv(apply_G16!, (copy(reg), (1,2,4,9), Js, S))
    reg1, _, _, S = apply_Gh!(copy(reg), 2, 0.5, S)
    reg2 = copy(reg) |> put(10, 2=>matblock(Diagonal(i_magtensor(0.5))))
    @test reg1 ≈ reg2
    @test check_inv(apply_Gh!, (copy(reg), 2, 0.5, S))
    reg1, _, S = apply_Gcp!(copy(reg), (2, 4), S)
    reg2 = copy(reg) |> put(10, (2, 4)=>matblock(copytensor(Tropical{Float64})))
    @test reg2.state ≈ map(x->x.n<-1000 ? Tropical(-Inf) : x, reg1.state)
    @test check_inv(apply_Gcp!, (copy(reg), (2, 4), S))
    reg1, _, S = apply_Greset!(copy(reg), 2, copy(S))
    reg2 = copy(reg) |> put(10, 2=>matblock(resettensor(Tropical{Float64})))
    @test reg1 ≈ reg2
    @test check_inv(apply_Greset!, (copy(reg), 2, S))
end

ext: jl


text:
# Metida
# Copyright © 2019-2020 Vladimir Arnautov aka PharmCat <mail@pharmcat.net>

__precompile__(true)
module MetidaBase

    using Tables, PrettyTables, TypedTables, StatsBase#, Reexport

    #@reexport using StatsModels

    import StatsModels: StatisticalModel
    import Tables: istable, columnaccess, columns, getcolumn, columnnames, schema, rowaccess, rows
    import CPUSummary: num_cores

    import Base: getindex, length, ht_keyindex, show, pushfirst!, iterate, size

    include("abstracttype.jl")
    include("types.jl")
    include("utils.jl")
    include("iterators.jl")

end


text:
mutable struct ActiveRenderPass
    mCmd::RecordingCommandBuffer
    mBeginInfo::vk.VkRenderPassBeginInfo
    mNextSubpass::UInt32

    function ActiveRenderPass(cmd::RecordingCommandBuffer, info::vk.VkRenderPassBeginInfo)
        this = new()
        this.mCmd = cmd
        this.mBeginInfo = info
        this.mNextSubpass = 0
        return this
    end
end

@class ActiveRenderPass

mutable struct InlineSubpass
    mPass::ActiveRenderPass
    mSubpassIndex::UInt32
    mCommandBuffer::RecordingCommandBuffer

    function InlineSubpass(pass::ActiveRenderPass, subpass::UInt32, commandBuffer::RecordingCommandBuffer)
        this = new()
        this.mPass = pass
        this.mSubpassIndex = subpass
        this.mCommandBuffer = commandBuffer
        return this
    end
end

@class InlineSubpass

function startInlineSubpass(this::ActiveRenderPass)::InlineSubpass
    cmd = handle(getBuffer(this.mCmd))
    if (this.mNextSubpass == 0)
        vk.vkCmdBeginRenderPass(cmd, Ref(this.mBeginInfo), vk.VK_SUBPASS_CONTENTS_INLINE)
    else
        vk.vkCmdNextSubpass(cmd, vk.VK_SUBPASS_CONTENTS_INLINE)
    end
    this.mNextSubpass += 1
    return InlineSubpass(this, this.mNextSubpass, this.mCmd)
end

# When ~ActiveRenderPass is called
function endRenderPass(this::ActiveRenderPass)
    if this.mNextSubpass != 0
        vk.vkCmdEndRenderPass(handle(this.mCmd))
    end
end

function bindPipeline(this::InlineSubpass, pip::GraphicsPipeline)
    setLastLayout(this.mCommandBuffer, getLayout(pip))
    vk.vkCmdBindPipeline(handle(this.mCommandBuffer), vk.VK_PIPELINE_BIND_POINT_GRAPHICS, handleRef(pip)[])
end

function setViewports(this::InlineSubpass, viewports::Vector{vk.VkViewport}, first::UInt32 = UInt32(0))
    vk.vkCmdSetViewport(handle(this.mCommandBuffer), first, length(viewports), pointer(viewports))
end

function bindVertexBuffers(this::InlineSubpass, buffers::Vector{Buffer}, first::UInt32 = UInt32(0))
    vkBuffers = Vector{vk.VkBuffer}()
    for b in buffers
        push!(vkBuffers, handle(b))
    end
    offsets = zeros(vk.VkDeviceSize, length(vkBuffers))

    vk.vkCmdBindVertexBuffers(handle(this.mCommandBuffer), first, length(vkBuffers), pointer(vkBuffers), pointer(offsets))
end

function bindIndexBuffer(this::InlineSubpass, buffer::Buffer, type::vk.VkIndexType = vk.VkIndexType(vk.VK_INDEX_TYPE_UINT32), offset::vk.VkDeviceSize = vk.VkDeviceSize(0))
    vk.vkCmdBindIndexBuffer(handle(this.mCommandBuffer), handle(buffer), offset, type)
end

function pushConstantBlock(this::InlineSubpass, data::T) where T
    ref = Ref(data)
    GC.@preserve ref begin
        pushConstantBlock(this, UInt32(sizeof_obj(data)), ref_to_pointer(Cvoid, ref))
    end
end

function pushConstantBlock(this::InlineSubpass, size::UInt32, data::Ptr{Cvoid})
    pushConstantBlock(this.mCommandBuffer, size, data);
end

function drawIndexed(this::InlineSubpass, indices::UInt32, instances::UInt32 = UInt32(1),
             vertexOffset::UInt32 = UInt32(0), firstIndex::UInt32 = UInt32(0), firstInstance::UInt32 = UInt32(0))
    vk.vkCmdDrawIndexed(handle(this.mCommandBuffer), indices, instances, firstIndex, vertexOffset, firstInstance)
end


text:
using ZipFile

struct ZipPathInfo # Rename to ZipEntry ?
    is_dir::Bool
    path::RelPath
end

struct ZipTreeRoot
    reader::ZipFile.Reader
    file_info::Vector{ZipPathInfo} # same order as reader.files
end

function ZipTreeRoot(path::AbstractString)
    reader = ZipFile.Reader(path)
    ZipTreeRoot(reader, _zip_path_info.(reader.files))
end

Base.show(io::IO, root::ZipTreeRoot) = print(io, "ZipTreeRoot(...)")

function _zip_path_info(f)
    name = f.name
    is_dir = endswith(name, '/')
    name = is_dir ? name[1:end-1] : name
    path = DataSets.RelPath(split(name, '/'))
    ZipPathInfo(is_dir, path)
end

#-------------------------------------------------------------------------------
# FIXME: Factor back together with BlobTree.jl !!

struct ZippedBlobTree <: AbstractBlobTree
    root::ZipTreeRoot
    path::RelPath
end

ZippedBlobTree(root::ZipTreeRoot) = ZippedBlobTree(root, RelPath())

Base.basename(tree::ZippedBlobTree) = basename(tree.path)

function Base.getindex(tree::ZippedBlobTree, path::RelPath)
    newpath = joinpath(tree.path, path)
    i = findfirst(tree.root.file_info) do info
        info.path == newpath
    end
    if i == nothing
        error("Path $newpath doesn't exist in $tree")
    elseif tree.root.file_info[i].is_dir
        ZippedBlobTree(tree.root, newpath)
    else
        Blob(tree.root, newpath)
    end
end

function Base.getindex(tree::ZippedBlobTree, name::AbstractString)
    getindex(tree, joinpath(RelPath(), name))
end

function _tree_children(tree::ZippedBlobTree)
    children = String[]
    for (i,info) in enumerate(tree.root.file_info)
        if dirname(info.path) == tree.path
            push!(children, basename(info.path))
        end
    end
    children
end

Base.IteratorSize(tree::ZippedBlobTree) = Base.SizeUnknown()
function Base.iterate(tree::ZippedBlobTree, state=nothing)
    if state == nothing
        children = _tree_children(tree)
        itr = iterate(children)
    else
        (children, cstate) = state
        itr = iterate(children, cstate)
    end
    if itr == nothing
        return nothing
    else
        (name, cstate) = itr
        (tree[name], (children, cstate))
    end
end

function Base.joinpath(tree::ZippedBlobTree, r::RelPath)
    # Should this AbsPath be rooted at `tree` rather than `tree.root`?
    AbsPath(tree.root, joinpath(tree.path, r))
end

function Base.joinpath(tree::ZippedBlobTree, s::AbstractString)
    AbsPath(tree.root, joinpath(tree.path, s))
end

function Base.open(func::Function, f::Blob{ZipTreeRoot}; write=false, read=!write)
    if write
        error("Error writing file at read-only path $f")
    end
    i = findfirst(info->info.path == f.path, f.root.file_info)
    io = f.root.reader.files[i]
    _seek(io, 0)
    func(io)
end

#=
w = ZipFile.Writer("/tmp/example.zip");
f = ZipFile.addfile(w, "hello.txt");
write(f, "hello world!\n");
f = ZipFile.addfile(w, "julia.txt", method=ZipFile.Deflate);
write(f, "Julia\n"^5);
close(w)
r = ZipFile.Reader("/tmp/example.zip");
for f in r.files
   println("Filename: $(f.name)")
   write(stdout, read(f, String));
end
close(r)
=#

#-------------------------------------------------------------------------------

# Fixes for ZipFile.jl
# TODO: Remove these once https://github.com/fhs/ZipFile.jl/pull/75 is in a release.

# It appears that ZipFile.jl just doesn't have a way to rewind to the start of
# one of the embedded files.
function _seek(io::ZipFile.ReadableFile, n::Integer)
    # Only support seeking to the start
    n == 0 || throw(ArguementError("Cannot efficiently seek zip stream to nonzero offset $n"))
    io._datapos = -1
    io._currentcrc32 = 0
    io._pos = 0
    io._zpos = 0
end

# Needed for use as `src` in `write(dst::IO, src::IO)`.
Base.readavailable(io::ZipFile.ReadableFile) = read(io)


text:
Generate = Command(
   "generate",
   "pangraph generate <options> [input fasta file]",
   "outputs a simulated sequence alignment graph",
   """zero or one fasta file.
      if zero, sequences are assumed to be supplied via stdin.""",
   [
    Arg(
        Float64,
        "mutation rate",
        (short="-m", long="--snp-rate"),
        "rate of mutations per site per genome per generation",
        1e-5,
    ),
    Arg(
        Float64,
        "HGT rate",
        (short="-r", long="--hgt-rate"),
        "rate of horizontal transfer events per genome per generation",
        0,
    ),
    Arg(
        Float64,
        "Length variance prefactor (σ=L/s)",
        (short="-s", long="--sigma-pre"),
        "divisor of mean length that sets the variance",
        10,
    ),
    Arg(
        Float64,
        "Deletion rate",
        (short="-d", long="--delete-rate"),
        "rate of deletion events per genome per generation",
        0,
    ),
    Arg(
        Float64,
        "Inversion rate",
        (short="-i", long="--invert-rate"),
        "rate of inversion events per genome per generation",
        0,
    ),
    Arg(
        String,
        "Graph output path",
        (short="-o", long="--graph-output"),
        "path where to output graph",
        "",
    ),
    Arg(
        Int,
        "Generations to simulate",
        (short="-t", long="--time"),
        "number of generations simulated under WF model",
        35,
    ),
   ],
   function(args)
        if args === nothing || length(args) == 0
            usage(Generate)
            return 2
        end

        input = parse(Generate, args)
        input = if (input === nothing || length(input) == 0)
            "/dev/stdin"
        elseif length(input) == 1
            input[1]
        else
            return 2
        end

        ancestors = open(input) do io
            map(r->r.seq, read_fasta(io))
        end

        N = length(ancestors)
        L = Int(round(sum(length.(ancestors)) / N))
        σ = L ÷ arg(Generate,"-s")

        evolve! = Simulation.model(
            Simulation.Params(;
                N   = N,
                L   = L,
                σ   = σ,
                snp = arg(Generate,"-m"),
                hgt = arg(Generate,"-r"),
                del = arg(Generate,"-d"),
                inv = arg(Generate,"-i"),
            )
        )

        T    = arg(Generate, "-t")
        path = arg(Generate, "-o")

        sequences, graph = Simulation.run(evolve!, T, ancestors; graph=path!="")

        for (i, sequence) in enumerate(sequences)
            write_fasta(stdout, "isolate_$(i)", sequence)
        end

        if length(path) > 0
            open(path,"w") do io
                marshal(io, graph; fmt=:json)
            end
        end
    end
)

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
using Sigma
using Lens
## Non-Negative Matrix Factorization

Sigma.restart_counter!()

dist(A,B; epsilon = 1.0) = sum([abs(A[i] - B[i]) for i = 1:length(A)]) < epsilon
dist3(A,B; epsilon = 1.0) = (&)([abs(A[i] - B[i]) < epsilon for i = 1:length(A)]...)
dist2(A,B; epsilon = 1.0) = sum([abs(A[i] - B[i]) for i = 1:length(A)])

function nmf(V, W, H; epsilon = 1.0)
  dist3(V, W*H; epsilon = epsilon)
end

function benchmark_nmf(m::Integer, k::Integer, n::Integer, nsamples::Integer,
                      epsilon::Float64, precision::Float64)
  Sigma.restart_counter!()
  srand(1234)
  W = Sigma.mvuniform(0.0,1.0,m,k)
  H = Sigma.mvuniform(0.0,1.0,k,n)

  # Generate date
  W_data = rand(m,2)
  H_data = rand(2,n)
  V = W_data*H_data
  # sol = Sigma.point_sample_mc(nmf(V, W, H; epsilon=0.05), nsamples, parallel = true, precision=0.001)
  #
  sol = rand((W,H), nmf(V, W, H; epsilon=epsilon), nsamples; parallel=true, precision=precision, RandVarType=Sigma.SymbolicRandVar)
  for i = 1:length(sol[1])
    W_conc = sol[1][i]
    H_conc = sol[2][i]
    prodWH = W_conc * H_conc
    lens(:distance, dist2(V,prodWH))
  end
  return sol, V
end

# function benchmark_nmf_rj(m::Integer, k::Integer, n::Integer, nsamples::Integer,
#                       epsilon::Float64, precision::Float64)
#   Sigma.restart_counter!()
#   srand(1234)
#   W = Sigma.mvuniform(0.0,1.0,m,k)
#   H = Sigma.mvuniform(0.0,1.0,k,n)
#
#   # Generate date
#   W_data = rand(m,2)
#   H_data = rand(2,n)
#   V = W_data*H_data
#   # sol = Sigma.point_sample_mc(nmf(V, W, H; epsilon=0.05), nsamples, parallel = true, precision=0.001)
#   #
#   nmf(V, W, H; epsilon=epsilon)
# end

function sample_variance{T,N}(samples::Vector{Array{T,N}})
  centerofmass::Array = mean(samples)
  distances = [sqrt(sum((sample - centerofmass).^2)) for sample in samples]
  sum(distances) / (length(samples) - 1)
end

function sample_variance(samples::Vector)
  centerofmass::Array = mean(samples)
  distances = [sqrt(sum((sample - centerofmass).^2)) for sample in samples]
  sum(distances) / (length(samples) - 1)
end

results, stats = capture(()->benchmark_nmf(3,2,3,100,0.01,0.001), [:distance, :sat_check, :post_loop])

timediffs = vcat(get(stats, proc_id=1, lensname=:sat_check), get(stats, proc_id=1, lensname=:sat_check))
timediffs = vcat(get(stats, proc_id=1, lensname=:sat_check), get(stats, proc_id=2, lensname=:sat_check))

timediffs = map(i->float(i)/1e9, timediffs)


text:
function test_index2onehot_layer(backend::Backend, tensor_dim, n_input, T, eps)
  println("-- Testing Index2OnehotLayer on $(typeof(backend)){$T}...")

  expand_dim = max(1, abs(rand(Int)) % tensor_dim)
  println("    > $tensor_dim-dimensional input, expanding along dimension $expand_dim")
  dims = [rand(6:11, tensor_dim) for i in 1:n_input]
  for i = 1:n_input
    dims[i][expand_dim] = 1
  end
  n_class = 6
  input = [convert(Array{T}, abs.(rand(Int, dims[i]...)) .% n_class) for i = 1:n_input]
  input_blob = Blob[make_blob(backend, x) for x in input]
  diff_blob = Blob[NullBlob() for i = 1:n_input]

  layer = Index2OnehotLayer(tops=Array{Symbol}(n_input), bottoms=Array{Symbol}(n_input),
      dim=expand_dim, n_class=n_class)
  state = setup(backend, layer, input_blob, diff_blob)

  forward(backend, state, input_blob)

  for i = 1:n_input
    my_dims  = size(input[i])
    dim_pre  = prod(my_dims[1:expand_dim-1])
    dim_post = prod(my_dims[expand_dim+1:end])
    canonical_input = reshape(input[i], (dim_pre, dim_post))
    canonical_output = zeros(T, dim_pre, n_class, dim_post)
    for x = 1:dim_pre
      for y = 1:dim_post
        canonical_output[x, convert(Int,canonical_input[x,y])+1, y] = 1
      end
    end

    output = reshape(canonical_output, size(state.blobs[i]))
    got_output = to_array(state.blobs[i])

    @test all(-eps .< output - got_output .< eps)
  end

  shutdown(backend, state)
end

function test_index2onehot_layer(backend::Backend, n_input, T, eps)
  for td in [2,4,5]
    test_index2onehot_layer(backend, td, n_input, T, eps)
  end
end
function test_index2onehot_layer(backend::Backend)
  test_index2onehot_layer(backend, 3, Float32, 1e-5)
  test_index2onehot_layer(backend, 3, Float64, 1e-10)
end

if test_cpu
  test_index2onehot_layer(backend_cpu)
end
if test_gpu
  test_index2onehot_layer(backend_gpu)
end


text:
"""
Describes the maximum discharge for a dam output, which can be constant (`Constant*DamOutput`) or depend on the
hydraulic head with a conveyance factor (`Conveyance*DamOutput`); the latter make the assumption that the reservoir level
is the hydraulic head. This discharge can be conditional or not, which is indicated by the `condition` attribute
(whose type must be `DamOutputCondition`).

These releases are considered as degrees of freedom of the operator, including for flood management: the operator
can decide on the exact discharge through each dam output (if the conditions are met). These releases are not
supposed to fulfil any dam purpose (like drinking water), which are either constant or valued (see `Purpose` for these).

Each output must have a unique identifier, whose type is `Symbol`. It is always available with the `getId()` function.
For implementers, a default implementation uses the field `id`. Likewise, the condition is available through
`getCondition()`, and a default implementation retrieves the value from the field `condition`.
"""
abstract DamOutput

getId(dout::DamOutput) = dout.id
getCondition(dout::DamOutput) = dout.condition
hasCondition(dout::DamOutput) = ! isempty(getCondition(dout))

id(dout::DamOutput) = getId(dout)
condition(dout::DamOutput) = getCondition(dout)
isConditional(dout::DamOutput) = hasCondition(dout)

"""
Describes the maximum discharge for a dam output using a conveyance factor, be it an open-channel spillway
(`ConveyanceSpillway`) or a closed conduit (`ConveyancePipe`).

A default implementation of `getConveyanceFactor()` takes its value from the `conveyance` field.
"""
abstract ConveyanceDamOutput <: DamOutput

hasConveyanceFactor(dout::DamOutput) = false
hasConveyanceFactor(dout::ConveyanceDamOutput) = true
getConveyanceFactor(dout::ConveyanceDamOutput) = dout.conveyance
conveyance(dout::DamOutput) = getConveyanceFactor(dout)

"""
Condition at which a dam output is enabled. For example, a sufficient reservoir level is needed for the spillways to be
reachable.
"""
abstract DamOutputCondition

"""
Describes a dam output condition: the dam release is always possible.
"""
immutable NoCondition <: DamOutputCondition; end

isempty(c::DamOutputCondition) = false
isempty(c::NoCondition) = true

"""
Describes a constant dam output, independent of all parameters. Its discharge is given in 10^6 m^3/s.
"""
immutable ConstantDamOutput <: DamOutput
  id::Symbol
  discharge::Float64
  condition::DamOutputCondition

  ConstantDamOutput(discharge::Discharge, condition::DamOutputCondition=NoCondition()) =
    new(:constant, _from_unitful(discharge), condition)
  ConstantDamOutput(id::Symbol, discharge::Discharge, condition::DamOutputCondition=NoCondition()) =
    new(id, _from_unitful(discharge), condition)
end

getDischarge(dout::DamOutput) = error("Discharge has no meaning for output " * string(getId(dout)) * ".")
getDischarge(dout::ConstantDamOutput) = dout.discharge
getDischarge(dout::DamOutput, p::Period) = getDischarge(dout) * Dates.days(p) * 86400
isConstant(dout::DamOutput) = false
isConstant(dout::ConstantDamOutput) = true

discharge(dout::DamOutput) = getDischarge(dout)
discharge(dout::DamOutput, p::Period) = getDischarge(dout, p)

"""
Describes the penstock driving a hydropower unit.

This object is automatically created when creating a `Reservoir` data structure, and thus does not need to be used
outside internal code.
"""
immutable HydropowerDamOutput <: DamOutput
  hp::HydropowerUnit
  condition::DamOutputCondition

  HydropowerDamOutput(hp::HydropowerUnit, cond::DamOutputCondition=NoCondition()) = new(hp, cond)
end

getId(hp::HydropowerDamOutput) = getDamOutputId(hp.hp)
getDischarge(hp::HydropowerDamOutput) = getMaximumDischarge(hp.hp)

isHydropower(dout::DamOutput) = false
isHydropower(hp::HydropowerDamOutput) = true
getHydropowerUnit(hp::HydropowerDamOutput) = hp.hp

hydropower(hp::HydropowerDamOutput) = getHydropowerUnit(hp)

"""
Describes the maximum flow of an open-channel spillway using a conveyance factor. This discharge is limited by the
available hydraulic head.

This kind of dam output obeys the following law:
    Q = conveyance (head)^(3/2)

TODO: How to deal with units? Q is 10^6 m^3/period!
"""
immutable ConveyanceSpillway <: ConveyanceDamOutput
  id::Symbol
  conveyance::Float64
  condition::DamOutputCondition

  ConveyanceSpillway(conveyance::Float64, condition::DamOutputCondition=NoCondition()) =
    new(:spillway, conveyance, condition)
  ConveyanceSpillway(id::Symbol, conveyance::Float64, condition::DamOutputCondition=NoCondition()) =
    new(id, conveyance, condition)
end

"""
Describes the maximum flow of a closed conduit using a conveyance factor. This discharge is limited by the available
hydraulic head.

This kind of dam output obeys the following law:
    Q = conveyance (head)^(1/2)

TODO: Is this law the right one?
TODO: How to deal with units? Q is 10^6 m^3/period, while the coefficient will be given for m^3/s!
"""
immutable ConveyancePipe <: ConveyanceDamOutput
  id::Symbol
  conveyance::Float64
  condition::DamOutputCondition

  ConveyancePipe(conveyance::Float64, condition::DamOutputCondition=NoCondition()) =
    new(:penstock, conveyance, condition)
  ConveyancePipe(id::Symbol, conveyance::Float64, condition::DamOutputCondition=NoCondition()) =
    new(id, conveyance, condition)
end

"""
Describes a dam output condition: the release is possible only if this condition is met, i.e. if the reservoir level is
sufficiently high. This type represents physical constraints on the actual reservoir level (not an approximation for
hydraulic head).
"""
immutable MinimumReservoirLevelCondition <: DamOutputCondition
  min::Float64
  MinimumReservoirLevelCondition(min::Height) = new(_from_unitful(min))
end

getMinimumLevel(c::MinimumReservoirLevelCondition) = c.min

getMinimumLevel(c::DamOutputCondition) = error("Minimum reservoir level has no meaning for the condition.")
minLevel(c::DamOutputCondition) = getMinimumLevel(c)

# Equivalent functions for volumes: see output_fwd.jl (they need the declarations of other stuff).

"""
Builds an array of dam outputs for the most common cases, i.e. one set of bottom outlets and one spillway.
"""
function ConstantDamOutputs(;bottomOutlets::Discharge=0.m^3/s,
                             spillway::Discharge=0.m^3/s, spillway_minLevel::Height=0.m)
  list = DamOutput[] # All data structures expect a generic array of dam outputs, be they constant or not.
  if bottomOutlets != 0.m^3/s
    push!(list, ConstantDamOutput(:bottomOutlets, bottomOutlets))
  end
  if spillway != 0.m^3/s
    if spillway_minLevel != 0.m
      push!(list, ConstantDamOutput(:spillway, spillway, MinimumReservoirLevelCondition(spillway_minLevel)))
    else
      push!(list, ConstantDamOutput(:spillway, spillway))
    end
  end
  return list
end


text:
#day1
input = "818275977931166178424892653779931342156567268946849597948944469863818248114327524824136924486891794739281668741616818614613222585132742386168687517939432911753846817997473555693821316918473474459788714917665794336753628836231159578734813485687247273288926216976992516314415836985611354682821892793983922755395577592859959966574329787693934242233159947846757279523939217844194346599494858459582798326799512571365294673978955928416955127211624234143497546729348687844317864243859238665326784414349618985832259224761857371389133635711819476969854584123589566163491796442167815899539788237118339218699137497532932492226948892362554937381497389469981346971998271644362944839883953967698665427314592438958181697639594631142991156327257413186621923369632466918836951277519421695264986942261781256412377711245825379412978876134267384793694756732246799739464721215446477972737883445615664755923441441781128933369585655925615257548499628878242122434979197969569971961379367756499884537433839217835728263798431874654317137955175565253555735968376115749641527957935691487965161211853476747758982854811367422656321836839326818976668191525884763294465366151349347633968321457954152621175837754723675485348339261288195865348545793575843874731785852718281311481217515834822185477982342271937155479432673815629144664144538221768992733498856934255518875381672342521819499939835919827166318715849161715775427981485233467222586764392783699273452228728667175488552924399518855743923659815483988899924199449721321589476864161778841352853573584489497263216627369841455165476954483715112127465311353411346132671561568444626828453687183385215975319858714144975174516356117245993696521941589168394574287785233685284294357548156487538175462176268162852746996633977948755296869616778577327951858348313582783675149343562362974553976147259225311183729415381527435926224781181987111454447371894645359797229493458443522549386769845742557644349554641538488252581267341635761715674381775778868374988451463624332123361576518411234438681171864923916896987836734129295354684962897616358722633724198278552339794629939574841672355699222747886785616814449297817352118452284785694551841431869545321438468118"
sum=0
for i=1:length(input)
    sum = parse(Int,input[i])==parse(Int,input[Int((i-1+length(input)/2)%length(input)+1)]) ? sum+parse(Int,input[i]) : sum
end
# if input[end] == input[1]
#     sum+=parse(Int,input[end])
# end
dfprintln(sum)

input = readdlm("/Users/ponet/Downloads/input.txt",Int)
i=1
j=0
while 0<i<=length(input)
    t=i
    i+=input[i]
    input[t] = input[t] >=3 ? input[t]-1:input[t]+1
    j+=1
end
dfprintln(j)


text:
module Smooth

#=
"""
! this subroutine returns a zero phase gaussian wi
ndow in time domain
n,& ! n must be either odd or a power of 2
! if n is odd, then the first nlag samples are negative lags
! followed by zero and positive lags
! if n is a power of 2, then lags are in wrap around order
! in this case one less negative lag is present, with gwin(1) as
! zero lag, followed by positive lags and negative lags are towards
! the end of the signal
nwin, & ! this parameter controls the width of the gaussian window in
! signal, the width of an impulse after smoothin will be 2*nwin
gwin, & ! output zerophase gaussain window in time domain
gwin_freq& ! optional output in frequency domain [only when n = npow2]
  ! only real part is returned after cfft
)
"""
function gaussian_filter(nwin::Int64, nlag::Int64=2*nwin)

	window = zeros(2*nlag+1)
	# time domain window
	if(nwin == 0)
		# if nwin = 0, then delta spike in time domain
		window[nlag + 1] =  1.0
	else
		# guassian window
		shrpcf = real(nwin)^(-2.0)
		for iw = 1: 2*nlag+1
			jsq          = ( nlag + 1 - iw )^2
			window[iw] = exp( -1.0 * shrpcf * real(jsq) )
		end
	end

	"""
	normalization of the filter such that the integral equals one
	necessary to preserve amplitudes after smoothing or whatever
	"""
	window ./= (sum(window))

end

"""
! A zerophase gaussian window of length 2*nwin is designed
! and used to smooth datain
! width of an impulsive arrival in datout after smoothin will be 2*nwin
! The window itself can be seen as an the impluse response of the smoothing filter
! higher values of nwin means more smoothing
datin, & ! [n1*n2] input data matrix
datout, &! [n1*n2] output data after smoothing
n1,&     ! first dimension of dat
n2,&     ! second dimension of dat
nwin&    ! half of the width of the gaussian window used for smoothing
"""
function gaussian{N}(x::AbstractArray{Float64,N}, nwin::Vector{Int64})
	xsize=[size(x)...]

	pa=P_conv(dsize=xsize, ssize=ssize, gsize=xsize,)

	# convolution with a Gaussian window
	if(N==1)
		gwin=gaussian_filter(nwin[1])
		xin=deepcopy(x)
		Conv.conv!(x, xin, gwin, :d)
	elseif(N==2)
		gwin1=gaussian_filter(nwin[1])
		gwin2=gaussian_filter(nwin[2])
		xin=deepcopy(x)
		for i in 1:size(x,2)
			xx=view(x,:,i)
			xxin=view(xin,:,i)
			Conv.conv!(xx, xxin, gwin1, :d)
		end
		xin=deepcopy(x)
		for i in 1:size(x,1)
			xx=view(x,i,:)
			xxin=view(xin,i,:)
			Conv.conv!(xx, xxin, gwin2, :d)
		end
	else
		error("only implemented upto 2 dimensions")
	end
	return x
end

function gaussian!(xg::Array{Float64,N},
		      x::Array{Float64,N}, pa) where {N}
	Conv.mod!(pa, :d, d=xg, g=x)
end
=#
end


text:
using Documenter, PowerSystems, DocStringExtensions, PowerFlows

makedocs(
    modules=[PowerFlows],
    format=Documenter.HTML(mathengine=Documenter.MathJax()),
    sitename="PowerFlows.jl",
    pages=Any[
        "Welcome Page" => "index.md",
        "Code Base Developer Guide" => Any["Developer Guide" => "code_base_developer_guide/developer.md",],
        "Public API Reference" => "api/public.md",
        "Internal API Reference" => "api/internal.md",
    ],
)

deploydocs(
    repo="github.com/NREL-SIIP/PowerFlows.jl.git",
    target="build",
    branch="gh-pages",
    devbranch="main",
    devurl="dev",
    push_preview=true,
    versions=["stable" => "v^", "v#.#"],
)


text:
# To create all plots execute this file from within this folder
# "julia create_plots.jl"
# Figures will be created in the 'figures' folder

include("plot_baseline_empirical.jl") # Fig. 1
include("plot_dynamics_no_policy.jl") # Fig. 2
include("plot_variation_policy_mutation.jl") # Fig. 3
include("plot_dynamics_policies.jl") # Fig. 4
include("plot_variation_policy_no_mutation.jl") # Fig. 5a
include("plot_variation_policy_no_mutation_pinf125.jl") # Fig. 5b
include("plot_boxplots_policy_A.jl") # Fig. 6a,b
include("plot_boxplots_policy_D.jl") # Fig. 6c,d
include("plot_vaccination_rollout.jl") # Fig. 7


text:
using Test, QSimulator
using QuadGK: quadgk

EC, EJ1, EJ2 = 0.172, 12.71, 3.69
freq = 0.1
amp = 0.6
harmonics = -50:50
parks = [0.0, 0.25]

@testset "FourierSeries" begin
	times = 0.0:.01:50
	for park in parks
		ϕ(t) = park + amp * sin(2π * freq * t)
		ω(t) = 2π * perturbative_transmon_freq(EC, EJ1, EJ2, ϕ(t))
		fs = FourierSeries(ω, freq, harmonics)
		ans = eval_series(fs, times)
		@test isapprox(ω.(times), ans, rtol=1e-9)
	end
end

@testset "rotating_frame_series" begin
	times = 0.0:.01:1
	for park in parks
		ϕ(t) = park + amp * sin(2π * freq * t)
		ω(t) = 2π * perturbative_transmon_freq(EC, EJ1, EJ2, ϕ(t))
		fs = FourierSeries(ω, freq, harmonics)
		Φ(t) = exp(1im * quadgk(ω, 0.0, t)[1])
		cfs = rotating_frame_series(fs, harmonics)
		# note that we must multiply by the term corresponding to fs.terms[0]
		# as this was ignored in rotating_frame_series
		ans = exp.(1im * fs.terms[0] * times) .* eval_series(cfs, times)
		@test isapprox(Φ.(times), ans, rtol=1e-6)
	end
end


text:
using KernelAbstractions
using KernelGradients
using Enzyme
using ROCKernels
using AMDGPU
using Test

include(joinpath(dirname(pathof(KernelAbstractions)), "..", "test", "testsuite.jl"))
include(joinpath(dirname(pathof(KernelGradients)), "..", "test", "testsuite.jl"))

if parse(Bool, get(ENV, "CI", "false"))
    default = "CPU"
else
    default = "ROCM"
end

backend = get(ENV, "KERNELABSTRACTIONS_TEST_BACKEND", default)

if backend != "ROCM"
    @info "ROCM backend not selected"
    exit()
end

if length(AMDGPU.get_agents(:gpu)) > 0
    AMDGPU.allowscalar(false)
    Testsuite.testsuite(ROCDevice, backend, AMDGPU, ROCArray, AMDGPU.ROCDeviceArray)
    GradientsTestsuite.testsuite(ROCDevice, backend, AMDGPU, ROCArray, AMDGPU.ROCDeviceArray)
else
    error("No AMD GPUs available!")
end


text:
# Randomized PCA

function rpca(X,        # (m,n) array
              k::Int,               # target rank
              p::Int=10,            # oversampling parameter
              q::Int=2)             # power iteration parameter

    m, n = size(X)
    U, Σ, W = rsvd(X, k, p, q)

    return W, (Σ .^ 2) ./ (m - 1), U * Diagonal(Σ)
end


text:
using CellFreeModelGenerationKit

# path to test metabolism file -
path_to_test_metabolism_file = "./test/data/Test.vff"

# need to build a model object -
build_result = build_julia_model_object(path_to_test_metabolism_file,"blank_for_now")
julia_model_object = check(build_result)

# parse -
result = parse_vff_model_document(julia_model_object)
ir_dictionary = check(result)

# look at the reaction table?
rxn_table = ir_dictionary[ir_master_reaction_table_key]
met_symbol_array = ir_dictionary[ir_list_of_molecular_species_key]
rxn_tag_array = ir_dictionary[ir_list_of_reaction_tags_key]

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
module EMLHEX_Engine
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("HEX_Engine/Properties_Average.jl")
	include("HEX_Engine/Properties_In_Out.jl")
	include("HEX_Engine/Properties_Wall.jl")
	include("HEX_Engine/Physical_Properties.jl")
	include("HEX_Engine/Physical_Properties_Heatex.jl")
	include("HEX_Engine/Tube_Pdrop.jl")
	include("HEX_Engine/Tube_Heat_Transfer.jl")
	include("HEX_Engine/Shell_Pdrop.jl")
	include("HEX_Engine/Shell_Heat_Transfer.jl")
	include("HEX_Engine/Baffles_Main.jl")
	include("HEX_Engine/Clearances_Main.jl")
	include("HEX_Engine/NTU_Basic.jl")
	include("HEX_Engine/LMTD_Basic.jl")
	include("HEX_Engine/Details_Main.jl")
	include("HEX_Engine/Tube_Side_Main.jl")
	include("HEX_Engine/Shell_Side_Main.jl")
	include("HEX_Engine/DoublePipe_HeatTransfer.jl")
	include("HEX_Engine/DoublePipe_PressureDrop.jl")
	include("HEX_Engine/Main_DoublePipe.jl")
end


text:
function neighborlist(r, rcut)

dim, N = size(r);

r1 = reshape(r,(dim, N, 1));
r2 = reshape(r,(dim, 1, N));
r12 = r1 .- r2;
d = reshape(sqrt.(sum(r12.*r12,dims=1)), (N,N)); # distances between atoms

A = (d .<= rcut);  # adjacency marix
index = 1:(N+1):N*N; # indices of diagonal entries
A[index] .= 0.0;  # set diagonal of A to zero

# pairs of atoms (i,j) within the cut-off radius rcut
ind = findall(A.==1);
ai = getindex.(ind, 2);
aj = getindex.(ind, 1);

# number of neighbors for each atom i
numneigh = accumarray(ai, ones(length(ai),1));

ai = Int64.(ai);
aj = Int64.(aj);
numneigh = Int64.(numneigh);

return ai, aj, numneigh

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
struct Classic_G2KP{D, S, P} end

is_collection(_::Val{:Classic_G2KP}) = false
is_collection(_::Classic_G2KP{D, S, P}) where {D, S, P} = false

@auto_hash_equals struct G2KP{D, S, P}
	"Length of the knapsack/'original plate'/'large object'."
	L :: S
	"Width of the knapsack/'original plate'/'large object'."
	W :: S
	"Lengths of the pieces."
	l :: Vector{S}
	"Widths of the pieces."
	w :: Vector{S}
	"Demands for the pieces (i.e., maximum number of copies that may be sold)."
	d :: Vector{D}
	"Profits obtained selling the one unit of the corresponding piece."
	p :: Vector{P}
end

"""
    read_from_string(_ :: Classic_G2KP{D, S, P}, s :: AbstractString)

Expect a string in the following format:
```
<L> <W>
<N>
<l_1> <w_1> <p_1> <d_1>
<l_2> <w_2> <p_2> <d_2>
...
```
`L` and `W` are the original plate length and width. `N` is the number of
pieces in the instance. The `l`, `w`, `p`, and `d`, are the length, width,
profit, and demand, of a piece.

Returns a `G2KP{D, S, P}` object.

# Notes

- If `N` is smaller than the number of lines after it, just the first `N`
  lines will be parsed and passed to the returned object.
- Empty lines and multiples spaces/tabs between columns may happen, but each
  non-empty line should contain the expected number of columns (i.e., the
  first non-empty line should have one number, the second non-empty line two
  numbers, all remaining non-empty lines four numbers).

See also: [`read_from_file`](@ref)
"""
@timeit TIMER function read_from_string(
	_ :: Classic_G2KP{D, S, P}, s :: AbstractString
) :: G2KP{D, S, P} where {D, S, P}
	N = zero(D)
	L = W = zero(S)
	(l, w, p, d) = (
		Vector{S}(), Vector{S}(), Vector{P}(), Vector{D}()
	)
	lnc = one(D)
	for ln in split(s, isequal('\n'); keepempty = false)
		if lnc - 2 == N + 1
			@warn "There are more piece lines than the specified number of pieces." *
				" The extra lines will not be used."
			break
		end
		if      lnc == (1::D)
			(L, W) = map(x->parse(S, x), split(ln))
		elseif  lnc == (2::D)
			N = parse(D, ln)
		else    # all other lines
			# cp_ == current plate
			piece_cols = split(ln)
			cp_l = parse(S, piece_cols[1])
			cp_w = parse(S, piece_cols[2])
			cp_p = parse(P, piece_cols[3])
			cp_d = parse(D, piece_cols[4])
			push!.((l, w, p, d), (cp_l, cp_w, cp_p, cp_d))
		end
		# The instances of paper "Improved state space relaxation for
		# constrained two-dimensional guillotine cutting problems" have a N
		# different of the real number of lines because the instances with 50
		# and 25 items are the same, only changing the value (they could
		# have removed the unused lines...).
		lnc = lnc + (1::D)
	end
	return G2KP{D, S, P}(L, W, l, w, d, p)
end

function read_from_string(
	_ :: Val{:Classic_G2KP}, s :: AbstractString
)
	return read_from_string(Classic_G2KP{Int, Int, Int}(), s)
end

# === `GuillotineModels.solution_value` implementations

import ..CutPattern
function solution_value(
	:: Val{:G2KP}, instance :: G2KP{D, S, P}, cp :: CutPattern{D, S}
) where {D, S, P}
	p = instance.p
	if iszero(cp.piece_idx)
		z = zero(eltype(p))
		for sp in cp.subpatterns
			z += solution_value(Val(:G2KP), instance, sp)
		end
		return z
	else
		return p[cp.piece_idx]
	end
end


text:
using CurveProximityQueries
using ConvexBodyProximityQueries
using IntervalArithmetic
using Plots
using StaticArrays

pgfplots()

include("EllipseBounds.jl")
using .EllipseBounds

curve = Bernstein([[-0.6, 0.0], [-0.6, 0.0], [0.1, -0.4], [0.4, 1.6], [-0.2, -1.6], [-0.1, 0.4], [0.8, 0.0], [0.8, 0.0]])
ell = Ellipse(curve, 0.2..0.8)

curve2 = Bernstein([[0.6, 0.6], [0.4, 1.0], [0.2, 0.2], [0.0, 0.6]])
ell2 = Ellipse(curve2, 0.2..0.9)

pts = closest_points(ell, ell2, @SVector(ones(2)))

plot(ell,
    linecolor=RGBA(102/255,130/225,223/255,0.5),
    fillcolor=RGBA(102/255,130/225,223/255,0.5),
    )
plot!(ell2,
     linecolor=RGBA(102/255,130/225,223/255,0.5),
     fillcolor=RGBA(102/255,130/225,223/255,0.5),
     )
plot!(curve, 0.0..0.2,
      linewidth=2.0,
      linestyle=:dash,
      linecolor=Gray(0.5),
     )
plot!(curve, 0.2..0.8,
      linewidth=2.0,
      linecolor=:black,
     )
plot!(curve, 0.8..1.0,
      linewidth=2.0,
      linestyle=:dash,
      linecolor=Gray(0.5),
     )
scatter!(curve(0.2), markersize=4, markercolor=:black)
scatter!(curve(0.8), markersize=4, markercolor=:black)

plot!(curve2, 0.0..0.2,
      linewidth=2.0,
      linestyle=:dash,
      linecolor=Gray(0.5),
     )
plot!(curve2, 0.2..0.9,
      linewidth=2.0,
      linecolor=:black,
     )
plot!(curve2, 0.9..1.0,
      linewidth=2.0,
      linestyle=:dash,
      linecolor=Gray(0.5),
     )
scatter!(curve2(0.2), markersize=4, markercolor=:black)
scatter!(curve2(0.9), markersize=4, markercolor=:black)

plot!(pts,
      linestyle=:dash,
      markersize=8.0,
      linewidth=2.0,
     )
savefig("out/fig4b.pdf")


text:
export Ensemble, CEnsemble, GCEnsemble

abstract type Ensemble end

struct CEnsemble <: Ensemble
  "coupling parameter"
  λ :: Float64
  "inverse temperature"
  β :: Float64
  "particle number"
  N :: Int
end

struct GCEnsemble <: Ensemble
  "coupling parameter"
  λ :: Float64
  "inverse temperature"
  β :: Float64
  "chemical potential"
  µ :: Float64
end


text:
module PyPlotly

using REPL
using Reexport: @reexport
using JSON
using PyCall: PyNULL, pyimport_conda

export plotly
const plotly = PyNULL()

include("graph_objects.jl")
@reexport using .GraphObjects

include("express.jl")
@reexport using .Express

Base.show(io, mime::MIME"text/html", fig::GraphObjects.Figure) = show(io, mime, fig.pyobj)

# Call fig.show() when one is trying to display `fig` object from Julia REPL
Base.display(d::REPL.REPLDisplay, fig::GraphObjects.Figure) = fig.show()

function __init__()
    copy!(plotly, pyimport_conda("plotly", "plotly", "plotly"))
end

end


text:
println("=======================================================================\n")
println("===============  TestInversion using Julia solver   ===================\n");
println("=======================================================================\n")

#############################################################################################################
modelDir ="../examples/";

dataDir = pwd();
include("../../examples/FWI/drivers/prepareFWIDataFiles.jl");
include("../../examples/FWI/drivers/setupFWI.jl");

plotting = false;

if plotting
	using  jInvVis
	using  PyPlot
	close("all");
end
#############################################################################################################
dataDir = pwd();
resultsDir = pwd();
########################################################################################################
dim     = 2;
pad     = 10;
jumpSrc = 15;
newSize = [128,64];
offset  = ceil(Int64,(newSize[1]*(13.5/13.5)));
println("Offset is: ",offset)
(m,Minv,mref,boundsHigh,boundsLow) = readModelAndGenerateMeshMref(modelDir,"SEGmodel2Dsalt.dat",dim,pad,[0.0,13.5,0.0,4.2],newSize,1.752,2.7);
maxBatchSize = 5;
omega = [0.5,1.0]*2*pi;

# ###################################################################################################################
dataFilenamePrefix = string(dataDir,"/DATA_SEG",tuple((Minv.n.+1)...));
resultsFilenamePrefix = string(resultsDir,"/SEG");
#######################################################################################################################

if plotting
	limits = [1.5,4.8];
	figure()
	plotModel(m,true,true,Minv,pad*0,limits);
	figure()
	plotModel(mref,true,true,Minv,pad*0,limits);
end

######################## DIRECT SOLVER #################################################

# Ainv = getMUMPSsolver([],0,0,2);
Ainv = getJuliaSolver();
# Ainv = getPARsolver([],0,0,6);

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

println("omega*maximum(h): ",omega*maximum(Minv.h)*sqrt(maximum(1.0./(boundsLow.^2))));

# This is a list of workers for FWI. Ideally they should be on different machines.
workersFWI = [workers()[1]];
println("The workers that we allocate for FWI are:");
println(workersFWI)

prepareFWIDataFiles(m,Minv,mref,boundsHigh,boundsLow,dataFilenamePrefix,omega,one(ComplexF64)*ones(size(omega)), pad,pad + 10,jumpSrc,
					offset,workersFWI,maxBatchSize,Ainv);

(Q,P,pMis,SourcesSubInd,contDiv,Iact,sback,mref,boundsHigh,boundsLow,resultsFilename) =
   setupFWI(m,dataFilenamePrefix,resultsFilenamePrefix,plotting,workersFWI,maxBatchSize,Ainv,SSDFun);

########################################################################################################
# Setting up the inversion for slowness
########################################################################################################
function dump(mc,Dc,iter,pInv,PMis,resultsFilename)
	fullMc = slowSquaredToVelocity(reshape(Iact*pInv.modelfun(mc)[1] .+ sback,tuple((pInv.MInv.n)...)))[1];
	# Temp = splitext(resultsFilename);
	# if iter>0
		# Temp = string(Temp[1],iter,Temp[2]);
	# else
		# Temp = resultsFilename;
	# end
	# if resultsFilename!=""
		# writedlm(Temp,convert(Array{Float16},fullMc));
	# end
	# if plotting
		# close(888);
		# figure(888);
		# plotModel(fullMc,true,false,[],0,[1.5,4.8],splitdir(Temp)[2]);
	# end
end

mref 		= velocityToSlow(mref)[1];
t    		= copy(boundsLow);
boundsLow 	= velocityToSlow(boundsHigh)[1];
boundsHigh 	= velocityToSlow(t)[1]; t = 0;
modfun1		= slowToSlowSquared;

########################################################################################################
# Set up Inversion #################################################################################
########################################################################################################

maxStep=0.05*maximum(boundsHigh);

regparams = [1.0,1.0,1.0,1e-4];
regfun(m,mref,M) = wdiffusionReg(m,mref,M,Iact=Iact,C=[]);
cgit = 2;
alpha = 1e-10;
pcgTol = 1e-1;
maxit = 2;

HesPrec=getExactSolveRegularizationPreconditioner();

pInv = getInverseParam(Minv,modfun1,regfun,alpha,mref[:],boundsLow,boundsHigh,
                         maxStep=maxStep,pcgMaxIter=cgit,pcgTol=pcgTol,
						 minUpdate=1e-3, maxIter = maxit,HesPrec=HesPrec);
mc = copy(mref[:]);
mc,Dc = freqCont(mc, pInv, pMis,contDiv, 3, resultsFilename,dump,"Joint",1,1,"projGN");

println("======================================================================================\n")
println("===============  TestInversion using shifted Laplacian multigrid   ===================\n")
println("======================================================================================\n")

######################## ITERATIVE SOLVER #############################################
levels      = 2;
numCores 	= 4;
BLAS.set_num_threads(numCores);
maxIter     = 50;
relativeTol = 1e-5;
relaxType   = "SPAI";
relaxParam  = 1.0;
relaxPre 	= 2;
relaxPost   = 2;
cycleType   ='W';
coarseSolveType = "NoMUMPS";
MG = getMGparam(ComplexF64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType,0.0,0.0);
shift = 0.2;
Hparam = HelmholtzParam(Minv,zeros(0),zeros(0),0.0,true,true);
AinvMG = getShiftedLaplacianMultigridSolver(Hparam, MG,shift);

(Q,P,pMis,SourcesSubInd,contDiv,Iact,sback,mref,boundsHigh,boundsLow,resultsFilename) =
   setupFWI(m,dataFilenamePrefix,resultsFilenamePrefix,plotting,workersFWI,maxBatchSize,AinvMG,SSDFun);

mref 		= velocityToSlow(mref)[1];
pInv.mref = mref[:];
mc = copy(mref[:]);

mc,Dc = freqCont(mc, pInv, pMis,contDiv, 3, resultsFilename,dump,"Joint",1,0,"projGN");

clear!(pMis);
##############################################################################################
rm("DATA_SEG(149, 75)_freq0.5.dat");
rm("DATA_SEG(149, 75)_freq1.0.dat");
rm("DATA_SEG(149, 75)_rcvMap.dat");
rm("DATA_SEG(149, 75)_srcMap.dat");
rm("DATA_SEG(149, 75)_PARAM.mat");
rm("SEG(149, 75)_FC2_HisGN.mat");
rm("SEG(149, 75)_FC1_HisGN.mat");
rm("SEG(149, 75)_Cyc1_FC2_HisGN.mat");
rm("SEG(149, 75)_Cyc1_FC1_HisGN.mat");
rm("jInv.out");