text:
module ResultTests

using DiffBase
using Base.Test

using StaticArrays

using DiffBase: DiffResult, GradientResult, JacobianResult, HessianResult,
                value, value!, derivative, derivative!, gradient, gradient!,
                jacobian, jacobian!, hessian, hessian!

k = 4
n0, n1, n2 = rand(), rand(), rand()
x0, x1, x2 = rand(k), rand(k, k), rand(k, k, k)
s0, s1, s2 = SVector{k}(rand(k)), SMatrix{k,k}(rand(k, k)), SArray{Tuple{k,k,k}}(rand(k, k, k))
rn = DiffResult(n0, n1, n2)
rx = DiffResult(x0, x1, x2)
rs = DiffResult(s0, s1, s2)
rsmix = DiffResult(n0, s0, s1)

issimilar(x, y) = typeof(x) == typeof(y) && size(x) == size(y)
issimilar(x::DiffResult, y::DiffResult) = issimilar(value(x), value(y)) && all(issimilar, zip(x.derivs, y.derivs))
issimilar(t::Tuple) = issimilar(t...)

##############
# DiffResult #
##############

@test rn === DiffResult(n0, n1, n2)
@test rx == DiffResult(x0, x1, x2)
@test rs === DiffResult(s0, s1, s2)
@test rsmix === DiffResult(n0, s0, s1)

@test issimilar(GradientResult(x0), DiffResult(first(x0), x0))
@test issimilar(JacobianResult(x0), DiffResult(x0, similar(x0, k, k)))
@test issimilar(JacobianResult(similar(x0, k + 1), x0), DiffResult(similar(x0, k + 1), similar(x0, k + 1, k)))
@test issimilar(HessianResult(x0), DiffResult(first(x0), x0, similar(x0, k, k)))

@test GradientResult(s0) === DiffResult(first(s0), s0)
@test JacobianResult(s0) === DiffResult(s0, zeros(SMatrix{k,k,Float64}))
@test JacobianResult(SVector{k+1}(vcat(s0, 0.0)), s0) === DiffResult(SVector{k+1}(vcat(s0, 0.0)), zeros(SMatrix{k+1,k,Float64}))
@test HessianResult(s0) === DiffResult(first(s0), s0, zeros(SMatrix{k,k,Float64}))

@test eltype(rn) === typeof(n0)
@test eltype(rx) === eltype(x0)
@test eltype(rs) === eltype(s0)

rn_copy = copy(rn)
@test rn == rn_copy
@test rn === rn_copy

rx_copy = copy(rx)
@test rx == rx_copy
@test rx !== rx_copy

rs_copy = copy(rs)
@test rs == rs_copy
@test rs === rs_copy

rsmix_copy = copy(rsmix)
@test rsmix == rsmix_copy
@test rsmix === rsmix_copy

# value/value! #
#--------------#

@test value(rn) === n0
@test value(rx) === x0
@test value(rs) === s0
@test value(rsmix) === n0

rn = value!(rn, n1)
@test value(rn) === n1
rn = value!(rn, n0)

x0_new, x0_copy = rand(k), copy(x0)
rx = value!(rx, x0_new)
@test value(rx) === x0 == x0_new
rx = value!(rx, x0_copy)

s0_new = rand(k)
rs = value!(rs, s0_new)
@test value(rs) == s0_new
@test typeof(value(rs)) === typeof(s0)
rs = value!(rs, s0)

rsmix = value!(rsmix, n1)
@test value(rsmix) === n1
rsmix = value!(rsmix, n0)

rn = value!(exp, rn, n1)
@test value(rn) === exp(n1)
rn = value!(rn, n0)

x0_new, x0_copy = rand(k), copy(x0)
rx = value!(exp, rx, x0_new)
@test value(rx) === x0 == exp.(x0_new)
rx = value!(rx, x0_copy)

s0_new = rand(k)
rs = value!(exp, rs, s0_new)
@test value(rs) == exp.(s0_new)
@test typeof(value(rs)) === typeof(s0)
rs = value!(rs, s0)

rsmix = value!(exp, rsmix, n1)
@test value(rsmix) === exp(n1)
rsmix = value!(rsmix, n0)

ksqrt = Int(sqrt(k))
T = typeof(SMatrix{ksqrt,ksqrt}(rand(ksqrt,ksqrt)))
rs_new = value!(rs, convert(T, value(rs)))
@test rs_new === rs

# derivative/derivative! #
#------------------------#

@test derivative(rn) === n1
@test derivative(rn, Val{2}) === n2

@test derivative(rx) === x1
@test derivative(rx, Val{2}) === x2

@test derivative(rs) === s1
@test derivative(rs, Val{2}) === s2

@test derivative(rsmix) === s0
@test derivative(rsmix, Val{2}) === s1

rn = derivative!(rn, n0)
@test derivative(rn) === n0
rn = derivative!(rn, n1)

x1_new, x1_copy = rand(k, k), copy(x1)
rx = derivative!(rx, x1_new)
@test derivative(rx) === x1 == x1_new
rx = derivative!(rx, x1_copy)

s1_new = rand(k, k)
rs = derivative!(rs, s1_new)
@test derivative(rs) == s1_new
@test typeof(derivative(rs)) === typeof(s1)
rs = derivative!(rs, s1)

s0_new = rand(k)
rsmix = derivative!(rsmix, s0_new)
@test derivative(rsmix) == s0_new
@test typeof(derivative(rsmix)) === typeof(s0)
rsmix = derivative!(rsmix, s0)

rn = derivative!(rn, n1, Val{2})
@test derivative(rn, Val{2}) === n1
rn = derivative!(rn, n2, Val{2})

x2_new, x2_copy = rand(k, k, k), copy(x2)
rx = derivative!(rx, x2_new, Val{2})
@test derivative(rx, Val{2}) === x2 == x2_new
rx = derivative!(rx, x2_copy, Val{2})

s2_new = rand(k, k, k)
rs = derivative!(rs, s2_new, Val{2})
@test derivative(rs, Val{2}) == s2_new
@test typeof(derivative(rs, Val{2})) === typeof(s2)
rs = derivative!(rs, s2, Val{2})

s1_new = rand(k, k)
rsmix = derivative!(rsmix, s1_new, Val{2})
@test derivative(rsmix, Val{2}) == s1_new
@test typeof(derivative(rsmix, Val{2})) === typeof(s1)
rsmix = derivative!(rsmix, s1, Val{2})

rn = derivative!(exp, rn, n0)
@test derivative(rn) === exp(n0)
rn = derivative!(rn, n1)

x1_new, x1_copy = rand(k, k), copy(x1)
rx = derivative!(exp, rx, x1_new)
@test derivative(rx) === x1 == exp.(x1_new)
rx = derivative!(exp, rx, x1_copy)

s1_new = rand(k, k)
rs = derivative!(exp, rs, s1_new)
@test derivative(rs) == exp.(s1_new)
@test typeof(derivative(rs)) === typeof(s1)
rs = derivative!(exp, rs, s1)

s0_new = rand(k)
rsmix = derivative!(exp, rsmix, s0_new)
@test derivative(rsmix) == exp.(s0_new)
@test typeof(derivative(rsmix)) === typeof(s0)
rsmix = derivative!(exp, rsmix, s0)

rn = derivative!(exp, rn, n1, Val{2})
@test derivative(rn, Val{2}) === exp(n1)
rn = derivative!(rn, n2, Val{2})

x2_new, x2_copy = rand(k, k, k), copy(x2)
rx = derivative!(exp, rx, x2_new, Val{2})
@test derivative(rx, Val{2}) === x2 == exp.(x2_new)
rx = derivative!(exp, rx, x2_copy, Val{2})

s2_new = rand(k, k, k)
rs = derivative!(exp, rs, s2_new, Val{2})
@test derivative(rs, Val{2}) == exp.(s2_new)
@test typeof(derivative(rs, Val{2})) === typeof(s2)
rs = derivative!(exp, rs, s2, Val{2})

s1_new = rand(k, k)
rsmix = derivative!(exp, rsmix, s1_new, Val{2})
@test derivative(rsmix, Val{2}) == exp.(s1_new)
@test typeof(derivative(rsmix, Val{2})) === typeof(s1)
rsmix = derivative!(exp, rsmix, s1, Val{2})

# gradient/gradient! #
#--------------------#

x1_new, x1_copy = rand(k, k), copy(x1)
rx = gradient!(rx, x1_new)
@test gradient(rx) === x1 == x1_new
rx = gradient!(rx, x1_copy)

s1_new = rand(k, k)
rs = gradient!(rs, s1_new)
@test gradient(rs) == s1_new
@test typeof(gradient(rs)) === typeof(s1)
rs = gradient!(rs, s1)

x1_new, x1_copy = rand(k, k), copy(x1)
rx = gradient!(exp, rx, x1_new)
@test gradient(rx) === x1 == exp.(x1_new)
rx = gradient!(exp, rx, x1_copy)

s0_new = rand(k)
rsmix = gradient!(exp, rsmix, s0_new)
@test gradient(rsmix) == exp.(s0_new)
@test typeof(gradient(rsmix)) === typeof(s0)
rsmix = gradient!(exp, rsmix, s0)

T = typeof(SVector{k*k}(rand(k*k)))
rs_new = gradient!(rs, convert(T, gradient(rs)))
@test rs_new === rs

# jacobian/jacobian! #
#--------------------#

x1_new, x1_copy = rand(k, k), copy(x1)
rx = jacobian!(rx, x1_new)
@test jacobian(rx) === x1 == x1_new
rx = jacobian!(rx, x1_copy)

s1_new = rand(k, k)
rs = jacobian!(rs, s1_new)
@test jacobian(rs) == s1_new
@test typeof(jacobian(rs)) === typeof(s1)
rs = jacobian!(rs, s1)

x1_new, x1_copy = rand(k, k), copy(x1)
rx = jacobian!(exp, rx, x1_new)
@test jacobian(rx) === x1 == exp.(x1_new)
rx = jacobian!(exp, rx, x1_copy)

s0_new = rand(k)
rsmix = jacobian!(exp, rsmix, s0_new)
@test jacobian(rsmix) == exp.(s0_new)
@test typeof(jacobian(rsmix)) === typeof(s0)
rsmix = jacobian!(exp, rsmix, s0)

T = typeof(SVector{k*k}(rand(k*k)))
rs_new = jacobian!(rs, convert(T, jacobian(rs)))
@test rs_new === rs

# hessian/hessian! #
#------------------#

x2_new, x2_copy = rand(k, k, k), copy(x2)
rx = hessian!(rx, x2_new)
@test hessian(rx) === x2 == x2_new
rx = hessian!(rx, x2_copy)

s2_new = rand(k, k, k)
rs = hessian!(rs, s2_new)
@test hessian(rs) == s2_new
@test typeof(hessian(rs)) === typeof(s2)
rs = hessian!(rs, s2)

x2_new, x2_copy = rand(k, k, k), copy(x2)
rx = hessian!(exp, rx, x2_new)
@test hessian(rx) === x2 == exp.(x2_new)
rx = hessian!(exp, rx, x2_copy)

s1_new = rand(k, k)
rsmix = hessian!(exp, rsmix, s1_new)
@test hessian(rsmix) == exp.(s1_new)
@test typeof(hessian(rsmix)) === typeof(s1)
rsmix = hessian!(exp, rsmix, s1)

T = typeof(SVector{k*k*k}(rand(k*k*k)))
rs_new = hessian!(rs, convert(T, hessian(rs)))
@test rs_new === rs

end # module


text:
"""
# Simple Least Squares fitting

The `CurveFit` module provides functions that
implement a few least squares approximations.

It is simple in an engineering sense. It simply
returns the coefficients and does not do
any error analysis.

It does, however, provide a simple and common interface
to the routines.

The package also includes nonlinear least squares fitting
using a Newton type algorithm

The fitting algorithms include

 * Straight lines
 * Polynomial fitting
 * Power laws
 * log laws
 * exp laws
 * Rational polynopmial fitting
 * A generic non-linear fitting algorithm
 * King's law (used in hotwire anemometry)

"""
module CurveFit
using Polynomials

export linear_fit, log_fit, power_fit, exp_fit, poly_fit
export LinearFit, LogFit, PowerFit, ExpFit
export curve_fit, apply_fit
export nonlinear_fit,  gauss_newton_fit, gauss_newton_generic_fit
export linear_king_fit, king_fit
export LinearKingFit, KingFit
export fit, coef
export linear_rational_fit, RationalPoly, ratval, rational_fit
export Poly

# package code goes here
#"Abstract base class for fitting data"
abstract type AbstractApproxFit end

#"Abstract class for least squares fitting of data"
abstract type AbstractLeastSquares <: AbstractApproxFit end

include("linfit.jl")
include("rationalfit.jl")
include("nonlinfit.jl")
include("king.jl")

"""
#Uses the object created by `curve_fit` to estimate values

The `call` method is overloaded so that the fit object can
be used as a function:

## Example:
x = 1.0:10.0
@. y = 2*x + 1 + randn()

fit = curve_fit(LinearFit, x, y)

y1 = fit(5.1)
y2 = apply_fit(fit, 5.1)
"""
apply_fit(f::T, x) where {T<:AbstractLeastSquares} = f(x)

end # module

ext: jl


text:
"Implicit Euler"
function getTableauImplicitEuler()
    a = ones(Float64, 1, 1)
    b = [1.0]
    c = [1.0]
    o = 1

    TableauFIRK(:implicit_euler, o, a, b, c)
end

"Implicit Midpoint"
function getTableauImplicitMidpoint()
    a = 0.5*ones(Float64, 1, 1)
    b = [1.0]
    c = [0.5]
    o = 2

    TableauFIRK(:implicit_midpoint, o, a, b, c)
end

"Gauss-Legendre Runge-Kutta, s=1"
function getTableauGLRK1()
    a = 0.5*ones(Float64, 1, 1)
    b = [1.0]
    c = [0.5]
    o = 2

    TableauFIRK(:glrk1, o, a, b, c)
end

"Gauss-Legendre Runge-Kutta, s=2"
function getTableauGLRK2()
    fac =
    a = [[0.25       0.25-√3/6]
         [0.25+√3/6  0.25     ]]
    b = [0.5,      0.5     ]
    c = [0.5-√3/6, 0.5+√3/6]
    o = 4

    TableauFIRK(:glrk2, o, a, b, c)
end

"Gauss-Legendre Runge-Kutta, s=3"
function getTableauGLRK3()
    a = [[5/36         2/9-√15/15  5/36-√15/30]
         [5/36+√15/24  2/9         5/36-√15/24]
         [5/36+√15/30  2/9+√15/15  5/36       ]]
    b = [5/18, 4/9, 5/18]
    c = [0.5-√15/10, 0.5, 0.5+√15/10]
    o = 6

    TableauFIRK(:glrk3, o, a, b, c)
end


text:
using DataFrames, MySQL, DBInterface, ProgressMeter, DataStructures

export DBconnect, DBclose, DBtransaction, DBprepare, DB, @DB_str, DBsource, upload_table!

const DBConn = Vector{DBInterface.Connection}()
DBconnect() = DBConn[end]

function DBclose()
    DBInterface.close!.(DBConn)
    empty!(DBConn)
end

function DBconnect(host; user=nothing, pass=nothing, dbname=nothing)
    DBclose()

    if !isnothing(user)  &&  isnothing(pass)
        pass = askpass("Enter password for DB user $user")
    end
    conn = DBInterface.connect(MySQL.Connection, host, user, pass)

    # "Driver={MariaDB};SERVER=127.0.0.1"
    #conn = ODBC.Connection(host *
    #                       (isnothing(user)  ?  ""  :  ";USER=" * user * (
    #                           isnothing(pass)  ?  ""  :  ";PWD=" * pass)))
    push!(DBConn, conn)

    if !isnothing(dbname)
        DB("USE $dbname")
    end
    nothing
end

DBtransaction(f) = MySQL.transaction(f, DBconnect())
DBprepare(sql::AbstractString) = DBInterface.prepare(DBconnect(), string(sql))

DB(stmt, params...) = DBInterface.execute(stmt, params)
function DB(stmt, df::DataFrame)
    desc = split(stmt.sql)[1] * join(fill(" ", 9))
    desc = desc[1:9] * " "
    barlen = ProgressMeter.tty_width(desc, stderr, false)
    (barlen > 50)  &&  (barlen = 50)
    prog = Progress(nrow(df), desc=desc, dt=0.5, color=:light_black, barlen=barlen,
                    barglyphs=BarGlyphs('|','█', ['▏','▎','▍','▌','▋','▊','▉'],' ','|',))
    DBtransaction() do
        for (i, row) in enumerate(Tables.rows(df))
            ProgressMeter.update!(prog, i)
            DBInterface.execute(stmt, Tables.Row(row))
        end
    end
    nothing
end

function DB(sql::AbstractString)
    out = DataFrame(DBInterface.execute(DBconnect(), string(sql)))
    if  (ncol(out) == 0)  &&
        (nrow(out) == 0)
        return nothing
    end
    if  (ncol(out) == 1)  &&
        (nrow(out) == 1)
        return out[1,1]
    end
    return out
end

macro DB_str(sql)
    return :(DB($sql))
end

function DBsource(file::AbstractString, subst::Vararg{Pair{String,String}, N}) where N
    delim = ";"
    sql = ""
    for line in readlines(file)
        if !isnothing(match(r"delimiter"i, line))
            if line == "DELIMITER ;"
                delim = ";"
            elseif line == "DELIMITER //"
                delim = "//"
            else
                error("Unexpected delimiter statement: $line")
            end
        else
            if !isnothing(match(Regex(delim), line))
                s = string.(strip.(split(line, delim)))
                @assert length(s) == 2 "Too many delimiters in one line: $line"
                sql *= s[1]
                for r in subst
                    sql = replace(sql, r)
                end
                println(sql)
                DB(sql)
                sql = ""
                (s[2] != "")  &&  (sql = s[2] * "\n")
            else
                sql *= line * "\n"
            end
        end
    end
    if strip(sql) != ""
        for r in subst
            sql = replace(sql, r)
        end
        println(sql)
        DB(sql)
    end
end

mutable struct DBColumn{T}
    sql::String
    null::Bool
    length::Int

    function DBColumn(data::AbstractVector{T}, sql) where T
        new{T}(sql, false, 0)
    end

    function DBColumn(data::AbstractVector{Union{Missing, T}}, sql) where T
        @assert !isa(T, Union)
        new{T}(sql, true, 0)
    end

    function DBColumn(data::AbstractVector{Vector{T}}, sql) where T
        @assert !isa(T, Union)
        new{T}(sql, false, length(data[1]))
    end
end

DBColumn(data::Vector{Union{Missing, Float32}}) = DBColumn(data, "FLOAT")
function DBColumn(data::Vector{Float32})
    if count(isnan.(data)) > 0
        data = convert(Vector{Union{Missing, Float32}}, data)
    end
    return DBColumn(data, "FLOAT")
end

DBColumn(data::Vector{Union{Missing, Float64}}) = DBColumn(data, "DOUBLE")
function DBColumn(data::Vector{Float64})
    if count(isnan.(data)) > 0
        data = convert(Vector{Union{Missing, Float64}}, data)
    end
    return DBColumn(data, "DOUBLE")
end

DBColumn(data::Vector{Union{Missing, UInt8 }}) = DBColumn(data, "TINYINT  UNSIGNED")
DBColumn(data::Vector{               UInt8  }) = DBColumn(data, "TINYINT  UNSIGNED")
DBColumn(data::Vector{Union{Missing,  Int8 }}) = DBColumn(data, "TINYINT    SIGNED")
DBColumn(data::Vector{                Int8  }) = DBColumn(data, "TINYINT    SIGNED")
DBColumn(data::Vector{Union{Missing, UInt16}}) = DBColumn(data, "SMALLINT UNSIGNED")
DBColumn(data::Vector{               UInt16 }) = DBColumn(data, "SMALLINT UNSIGNED")
DBColumn(data::Vector{Union{Missing,  Int16}}) = DBColumn(data, "SMALLINT   SIGNED")
DBColumn(data::Vector{                Int16 }) = DBColumn(data, "SMALLINT   SIGNED")
DBColumn(data::Vector{Union{Missing, UInt32}}) = DBColumn(data, "     INT UNSIGNED")
DBColumn(data::Vector{               UInt32 }) = DBColumn(data, "     INT UNSIGNED")
DBColumn(data::Vector{Union{Missing,  Int32}}) = DBColumn(data, "     INT   SIGNED")
DBColumn(data::Vector{                Int32 }) = DBColumn(data, "     INT   SIGNED")
DBColumn(data::Vector{Union{Missing, UInt64}}) = DBColumn(data, "  BIGINT UNSIGNED")
DBColumn(data::Vector{               UInt64 }) = DBColumn(data, "  BIGINT UNSIGNED")
DBColumn(data::Vector{Union{Missing,  Int64}}) = DBColumn(data, "  BIGINT   SIGNED")
DBColumn(data::Vector{                Int64 }) = DBColumn(data, "  BIGINT   SIGNED")
DBColumn(data::Vector{Union{Missing,   Bool}}) = DBColumn(data, "TINYINT    SIGNED")
DBColumn(data::Vector{                 Bool }) = DBColumn(data, "TINYINT    SIGNED")

function DBColumn(data::Vector{Union{Missing, Symbol}})
    i = findall(.!ismissing.(data))
    @assert length(i) > 0 "All values are missing"
    return DBColumn(data, "ENUM(" * join("'" .* sort(unique(string.(data[i]))) .* "'", ", ") * ")")
end

DBColumn(data::Vector{Symbol}) =
    return DBColumn(data, "ENUM(" * join("'" .* sort(unique(string.(data))) .* "'", ", ") * ")")

function DBColumn(data::Vector{Union{Missing, String}})
    i = findall(.!ismissing.(data))
    if length(i) == 0
        @warn "All values are missing, assuming a length of 20"
        maxlen = 20
    else
        maxlen = maximum(length.(data[i]))
    end
    return DBColumn(data, "VARCHAR($(maxlen))")
end

function DBColumn(data::Vector{String})
    maxlen = maximum(length.(data))
    return DBColumn(data, "VARCHAR($(maxlen))")
end

prepare_column!(data::DataFrame, col::DBColumn, name::Symbol) =
    col.null  &&  allowmissing!(data, name)

function prepare_column!(data::DataFrame, col::Union{DBColumn{Float32}, DBColumn{Float64}}, name::Symbol)
    i = findall(.!ismissing.(data[:, name])  .&  isnan.(data[:, name]))
    @assert (length(i) == 0)  ||  col.null
    if col.null
        allowmissing!(data, name)
        data[i, name] .= missing
    end
end

function prepare_column!(data::DataFrame, col::DBColumn{Bool}, name::Symbol)
    data[!, name] = fill(Int8(0), nrow(data))
    col.null  &&  allowmissing!(data, name)
    i = findall(  ismissing.(data[:, name]));  data[i, name] .= missing
    i = findall(.!ismissing.(data[:, name]));  data[i, name] .= Int8(data[:, name])
end

function DBColumns(data::DataFrame)
    out = OrderedDict{Symbol, DBColumn}()
    for name in Symbol.(names(data))
        out[name] = DBColumn(data[:, name])
    end
    return out
end

upload_table!(data::DataFrame, tbl_name::String; kw...) =
    upload_table!(data, DBColumns(data), tbl_name; kw...)

function upload_table!(data::DataFrame, meta::OrderedDict{Symbol, DBColumn}, tbl_name::String; drop=true, temp=false, memory=false)
    coldefs = Vector{String}()
    for name in Symbol.(names(data))
        @info "Preparing column $name"
        col = meta[name]
        if col.length == 0
            prepare_column!(data, col, name)
            push!(coldefs, "`$(name)` $(col.sql) " * (col.null  ?  ""  :  " NOT NULL"))
        else
            for i in 1:col.length
                data[!, Symbol(name, i)] .= getindex.(data[:, name], i)
                prepare_column!(data, col, Symbol(name, i))
                push!(coldefs, "`$(name)$(i)` $(col.sql) " * (col.null  ?  ""  :  " NOT NULL"))
            end
            select!(data, Not(name))
        end
    end

    if drop
        DB("DROP TABLE IF EXISTS $tbl_name")
        sql = "CREATE " * (temp ? "TEMPORARY" : "") * " TABLE $tbl_name"
        sql *= " ( " * join(coldefs, ", ") * ")"
        memory  &&  (sql *= " ENGINE=MEMORY")
        println(sql)
        DB(sql)
    end
    params = join(repeat("?", ncol(data)), ",")
    stmt = DBprepare("INSERT INTO $tbl_name VALUES ($params)")
    DB(stmt, data)
    nothing
end


text:
module PGFPlots

export LaTeXString, @L_str, @L_mstr
export plot, ErrorBars, Axis, Axes, PolarAxis, SmithAxis, TernaryAxis, GroupPlot, Plots, ColorMaps, save, define_color
export pushPGFPlotsOptions, popPGFPlotsOptions, resetPGFPlotsOptions, pgfplotsoptions
export pushPGFPlotsPreamble, popPGFPlotsPreamble, resetPGFPlotsPreamble, pgfplotspreamble
export pushPGFPlots, popPGFPlots
import Colors: RGB
import Contour: contours, levels

using Discretizers
using DataFrames
using Printf

mutable struct ErrorBars
    data::AbstractMatrix{Real}
    style
    mark
    ErrorBars(data::AbstractMatrix{T}; style=nothing, mark=nothing) where {T <: Real} = new(data, style, mark)
end

mylength(x) = (x == nothing) ? 0 : length(x)

function ErrorBars(; x=nothing, y=nothing, xplus=nothing, xminus=nothing, yplus=nothing, yminus=nothing, style=nothing, mark=nothing)
    if x != nothing
        xplus = x
        xminus = x
    end
    if y != nothing
        yplus = y
        yminus = y
    end
    n = maximum(map(mylength, Any[xplus, xminus, yplus, yminus]))
    xplus = mylength(xplus) == n ? xplus : zeros(n)
    xminus = mylength(xminus) == n ? xminus : zeros(n)
    yplus = mylength(yplus) == n ? yplus : zeros(n)
    yminus = mylength(yminus) == n ? yminus : zeros(n)
    data = hcat(xplus, xminus, yplus, yminus)'
    ErrorBars(data, style=style, mark=mark)
end

include("colormaps.jl")
include("plots.jl")

Plots.Linear(df::DataFrame; x::Symbol=:x, y::Symbol=:y, kwargs...) = Plots.Linear(df[!,x], df[!,y]; kwargs...)
Plots.Linear3(df::DataFrame; x::Symbol=:x, y::Symbol=:y, z::Symbol=:z, kwargs...) = Plots.Linear3(df[!,x], df[!,y], df[!,z]; kwargs...)
Plots.Histogram(df::DataFrame; x::Symbol=:x, kwargs...) = Plots.Histogram(df[!,x]; kwargs...)
Plots.Scatter(df::DataFrame; x::Symbol=:x, y::Symbol=:y, kwargs...) = Plots.Scatter(df[!,x], df[!,y]; kwargs...)
Plots.Quiver(df::DataFrame; x::Symbol=:x, y::Symbol=:y, u::Symbol=:u, v::Symbol=:v, kwargs...) = Plots.Quiver(convert(Vector, df[!,x]), convert(Vector, df[!,y]), convert(Vector, df[!,u]), convert(Vector, df[!,v]); kwargs...)
Plots.Histogram2(df::DataFrame; x::Symbol=:x, y::Symbol=:y, kwargs...) = Plots.Histogram2(convert(Vector, df[!,x]), convert(Vector, df[!,y]); kwargs...)
Plots.Histogram2(df::DataFrame, edges_x::AbstractVector{C}, edges_y::AbstractVector{C}; x::Symbol=:x, y::Symbol=:y, kwargs...) where {C<:Real} = Plots.Histogram2(convert(Vector, df[!,x]), convert(Vector, df[!,y]), edges_x, edges_y; kwargs...)

import TikzPictures: TikzPicture, PDF, TEX, TIKZ, SVG, save, LaTeXString, @L_str, @L_mstr

_pgfplotsoptions = Any[""]

pushPGFPlotsOptions(options::AbstractString) = push!(_pgfplotsoptions, options)
function popPGFPlotsOptions()
    if length(_pgfplotsoptions) > 1
        pop!(_pgfplotsoptions)
    end
end
resetPGFPlotsOptions() = deleteat!(_pgfplotsoptions, 2:length(_pgfplotsoptions))

function pushPGFPlots(options::AbstractString)
    warn("Use pushPGFPlotsOptions")
    pushPGFPlotsOptions(options)
end

function popPGFPlots()
    warn("Use popPGFPlotsOptions")
    popPGFPlotsOptions()
end

pgfplotsoptions() = join(_pgfplotsoptions, ",\n")

_pgfplotspreamble = Any[read(joinpath(dirname(@__FILE__), "preamble.tex"),String)]

pushPGFPlotsPreamble(preamble::AbstractString) = push!(_pgfplotspreamble, preamble)
function popPGFPlotsPreamble()
    if length(_pgfplotspreamble) > 1
        pop!(_pgfplotspreamble)
    end
end
resetPGFPlotsPreamble() = deleteat!(_pgfplotspreamble,2:length(_pgfplotspreamble))

pgfplotspreamble() = join(_pgfplotspreamble, "\n")

function define_color(name::AbstractString, color::RGB)
    _pgfplotspreamble[end] = _pgfplotspreamble[end] * "\n\\definecolor{$name}{rgb}{$(color.r), $(color.g), $(color.b)}"
end

define_color(name::AbstractString, color) = define_color(name, convert(RGB, color))

function define_color(name::AbstractString, color::Vector{T}) where {T<:AbstractFloat}
    if length(color) == 3
        r, g, b = color[1], color[2], color[3]
        _pgfplotspreamble[end] = _pgfplotspreamble[end] * "\n\\definecolor{$name}{rgb}{$r, $g, $b}"
    elseif length(color) == 4
        c, m, y, k = color[1], color[2], color[3], color[4]
        _pgfplotspreamble[end] = _pgfplotspreamble[end] * "\n\\definecolor{$name}{cmyk}{$c, $m, $y, $k}"
    else
        error("Color must have three or four components")
    end
end

function define_color(name::AbstractString, color::Vector{T}) where {T<:Integer}
    @assert(length(color) == 3, "Color must have three components")
    r, g, b = color[1], color[2], color[3]
    _pgfplotspreamble[end] = _pgfplotspreamble[end] * "\n\\definecolor{$name}{RGB}{$r, $g, $b}"
end

function define_color(name::AbstractString, color::UInt32)
    _pgfplotspreamble[end] = _pgfplotspreamble[end] * "\n\\definecolor{$name}{HTML}{$(uppercase(string(color, base=16)))}"
end

function define_color(name::AbstractString, color::Real)
    _pgfplotspreamble[end] = _pgfplotspreamble[end] * "\n\\definecolor{$name}{gray}{$color}"
end

histogramMap = Dict(
    :bins => "bins",
    :density => "density",
    :cumulative => "cumulative"
    )

linearMap = Dict(
    :mark => "mark",
    :markSize => "mark size",
    :style => "",
    :onlyMarks => "only marks"
    )

barMap = Dict(
    :style => ""
    )

scatterMap = Dict(
    :mark => "mark",
    :markSize => "mark size",
    :style => "",
    :onlyMarks => "only marks",
    :scatterClasses => "scatter/classes"
    )

errorbarsMap = Dict(
    :mark => "error mark",
    :style => "error bar style"
    )

quiverMap = Dict(
    :style => ""
    )

contourMap = Dict(
    :number => "number",
    :levels => "levels",
    :style => "",
    :labels => "labels"
    )

patch2DMap = Dict(
    :style => "",
    :patch_type => "patch type",
    :shader => "shader",
    :legendentry => "legendentry"
    )

using .Plots
using .ColorMaps

const IntegerRange = Tuple{Integer,Integer}
const RealRange = Tuple{Real,Real}

mutable struct Axis
    plots::Vector{Plot}
    title
    xlabel
    xlabelStyle
    ylabel
    ylabelStyle
    zlabel
    zlabelStyle
    xmin
    xmax
    ymin
    ymax
    axisEqual
    axisEqualImage
    enlargelimits
    axisOnTop
    view
    width
    height
    style
    legendPos
    legendStyle
    xmode
    ymode
    colorbar
    colorbarStyle
    hideAxis
    axisLines
    axisKeyword

    Axis(plots::Vector{P};title=nothing, xlabel=nothing, xlabelStyle=nothing, ylabel=nothing, ylabelStyle=nothing, zlabel=nothing, zlabelStyle=nothing, xmin=nothing, xmax=nothing,
                    ymin=nothing, ymax=nothing, axisEqual=nothing, axisEqualImage=nothing, enlargelimits=nothing, axisOnTop=nothing, view=nothing, width=nothing,
                    height=nothing, style=nothing, legendPos=nothing, legendStyle=nothing, xmode=nothing, ymode=nothing, colorbar=nothing, colorbarStyle=nothing, hideAxis=nothing, axisLines=nothing, axisKeyword="axis") where {P <: Plot} =
        new(plots, title, xlabel, xlabelStyle, ylabel, ylabelStyle, zlabel, zlabelStyle, xmin, xmax, ymin, ymax, axisEqual, axisEqualImage, enlargelimits, axisOnTop, view, width, height, style, legendPos, legendStyle, xmode, ymode, colorbar, colorbarStyle, hideAxis, axisLines, axisKeyword
            )

    Axis(;kwargs...) = Axis(Plot[]; kwargs...)
    Axis(plot::Plot; kwargs...) = Axis(Plot[plot]; kwargs...)

    Axis(plot::Contour; kwargs...) = Axis(Plot[plot]; kwargs..., view="{0}{90}")
    Axis(plots::Vector{Contour}; kwargs...) = Axis(Plot[plots...]; kwargs..., view="{0}{90}")
end

function PolarAxis(args...; kwargs...)
    ternary_axis_string = "\\usepgfplotslibrary{polar}"
    if ternary_axis_string ∉ _pgfplotspreamble
        pushPGFPlotsPreamble(ternary_axis_string)
    end
    return Axis(args...; kwargs..., axisKeyword = "polaraxis")
end

function TernaryAxis(args...; kwargs...)
    ternary_axis_string = "\\usepgfplotslibrary{ternary}"
    if ternary_axis_string ∉ _pgfplotspreamble
        pushPGFPlotsPreamble(ternary_axis_string)
    end
    return Axis(args...; kwargs..., axisKeyword = "ternaryaxis")
end

function SmithAxis(args...; kwargs...)
    ternary_axis_string = "\\usepgfplotslibrary{smithchart}"
    if ternary_axis_string ∉ _pgfplotspreamble
        pushPGFPlotsPreamble(ternary_axis_string)
    end
    return Axis(args...; kwargs..., axisKeyword = "smithchart")
end

const Axes = Vector{Axis}

function Base.push!(g::Axis, p::Plot)
    push!(g.plots, p)
end

axisMap = Dict(
    :title => "title",
    :xlabel => "xlabel",
    :xlabelStyle => "xlabel style",
    :ylabel => "ylabel",
    :ylabelStyle => "ylabel style",
    :zlabel => "zlabel",
    :zlabelStyle => "zlabel style",
    :xmin => "xmin",
    :xmax => "xmax",
    :ymin => "ymin",
    :ymax => "ymax",
    :axisEqual => "axis equal",
    :axisEqualImage => "axis equal image",
    :enlargelimits => "enlargelimits",
    :axisOnTop => "axis on top",
    :view => "view",
    :width => "width",
    :height => "height",
    :style => "",
    :legendPos => "legend pos",
    :legendStyle => "legend style",
    :xmode => "xmode",
    :ymode => "ymode",
    :colorbar => "colorbar",
    :colorbarStyle => "colorbar style",
    :hideAxis => "hide axis",
    :axisLines => "axis lines"
    )

mutable struct GroupPlot
    axes::AbstractArray{Axis,1}
    dimensions::IntegerRange
    style
    groupStyle
    GroupPlot(;style=nothing,groupStyle=nothing) = new(AbstractArray{Axis,1}[], (0,0), style, groupStyle)
    GroupPlot(axes::AbstractArray{Axis,1};style=nothing,groupStyle=nothing) = new(axes, (length(axes), 1), style, groupStyle)
    GroupPlot(plots;style=nothing,groupStyle=nothing) = new([Axis(p) for p in plots], (length(plots), 1), style, groupStyle)
    GroupPlot(dimensions::IntegerRange;style=nothing,groupStyle=nothing) = new(AbstractArray{Axis,1}[], dimensions, style, groupStyle)
    GroupPlot(rows::Integer, columns::Integer;style=nothing,groupStyle=nothing) = new(AbstractArray{Axis,1}[], (rows, columns), style, groupStyle)
end

function Base.push!(g::GroupPlot, a::Axis)
    push!(g.axes, a)
    if length(g.axes) > prod(g.dimensions)
        g.dimensions = (length(g.axes), 1)
    end
    g
end
Base.push!(g::GroupPlot, p::Plot) = push!(g, Axis(p))
function Base.append!(g::GroupPlot, arr::AbstractVector{T}) where {T<:Union{Plot,Axis}}
    for a in arr
        push!(g, a)
    end
    g
end

function printList(o::IO, a::AbstractVector; brackets=false)
    first = true
    for elem in a
        if first
            first = false
            if brackets
                print(o, "{")
            end
            print(o, "$elem")
        else
            print(o, ", $elem")
        end
    end
    if !first && brackets
        print(o, "}")
    end
end

printObject(o::IO, object) = print(o, "$(object)")
printObject(o::IO, object::AbstractVector) = printList(o, object, brackets=true)

function optionHelper(o::IO, m, object; brackets=false, otherOptions=Dict{AbstractString,AbstractString}[], otherText=nothing)
    first = true
    for (sym, str) in m
        if getfield(object,sym) != nothing
            if first
                first = false
                if brackets
                    print(o, "[")
                end
            else
                print(o, ", ")
            end
            if length(str) > 0
                print(o, "$str = {")
            end
            printObject(o, getfield(object,sym))
            if length(str) > 0
                print(o, "}")
            end
        end
    end
    for (k, v) in otherOptions
        if first
            first = false
            if brackets
                print(o, "[")
            end
        else
            print(o, ", ")
        end
        print(o, "$k = $v")
    end
    if otherText != nothing
        for t in otherText
            if t != nothing
                if first
                    first = false
                    if brackets
                        print(o, "[")
                    end
                else
                    print(o, ", ")
                end
                print(o, "$t")
            end
        end
    end
    if !first && brackets
        print(o, "]")
    end
end

function plotHelper(o::IO, p::Plots.Histogram)
    if (p.discretization == :pgfplots && p.bins > 0) ||
       (p.discretization == :default && length(p.data) ≤ Plots.THRESHOLD_NSAMPLES_DISC_OURSELVES)

        print(o, "\\addplot+ [mark=none, $(p.style), hist={")
        optionHelper(o, histogramMap, p)
        print(o, "}] table [row sep=\\\\, y index = 0] {")
        println(o, "data\\\\")
        for d in p.data
            println(o, "$d \\\\ ")
        end
        println(o, "};")

    else
        # discretize using Discretizers.jl

        if p.discretization == :specified && p.bins > 0
            edges = binedges(DiscretizeUniformWidth(p.bins), convert(Vector{Float64}, p.data))
        else
            discretization = p.discretization == :default ? :auto : p.discretization
            edges = binedges(DiscretizeUniformWidth(discretization), convert(Vector{Float64}, p.data))
        end

        linear = Plots._construct_histogram_linear_data(p.data, edges, p.density, p.cumulative)
        if p.style != "fill=blue!10"
            linear.style = p.style
            if !occursin("ybar interval", linear.style)
                linear.style = "ybar interval," * linear.style
            end
        end
        plotHelper(o, linear)
    end
end

function plotHelper(o::IO, p::BarChart)
    print(o, "\\addplot+ ")
    if p.errorBars == nothing
        optionHelper(o, barMap, p, brackets=true)
        println(o, "coordinates {")
        for (k,v) in zip(p.keys, p.values)
            println(o, "($k, $v)")
        end
        println(o, "};")
    else
        plotHelperErrorBars(o, p)
    end
    plotLegend(o, p.legendentry)
end
function PGFPlots.Axis(p::BarChart; kwargs...)
    # TODO : What's a better way to do this?
    style = "ybar=0pt, bar width=18pt, xtick=data, symbolic x coords={$(Plots.symbolic_x_coords(p))},"
    kwargs_unsplat = isempty(kwargs) ? Array{Pair{Symbol,Any}}(undef,0) : convert(Array{Pair{Symbol,Any}},collect(kwargs))
    i = findfirst(tup->tup[1] == :style, kwargs_unsplat)
    if i != nothing
        kwargs_unsplat[i] = :style=>(style * kwargs_unsplat[i][2])
    else
        push!(kwargs_unsplat,:style=>style)
    end
    i = findfirst(tup->tup[1] == :ymin, kwargs_unsplat)
    if i == nothing
        push!(kwargs_unsplat, :ymin=>0)
    end

    return Axis(Plots.Plot[p]; kwargs_unsplat...)
end

plotLegend(o::IO, entry) = nothing
plotLegend(o::IO, entry::AbstractString) = println(o, "\\addlegendentry{$entry}")
function plotLegend(o::IO, entries::Vector{T}) where {T <: AbstractString}
    for entry in entries
        plotLegend(o, entry)
    end
end

# todo: add error bars style
function plotHelperErrorBars(o::IO, p::Linear)
    println(o, "[")
    optionHelper(o, linearMap, p, brackets=false)
    println(o, ", error bars/.cd, ")
    optionHelper(o, errorbarsMap, p.errorBars, otherText=["x dir=both", "x explicit", "y dir=both", "y explicit"])
    println(o, "]")
    x = vcat(p.data, p.errorBars.data)
    println(o, "table [")
    println(o, "x error plus=ex+, x error minus=ex-, y error plus=ey+, y error minus=ey-")
    println(o, "] {")
    println(o, "x y ex+ ex- ey+ ey-")
    for i = 1:size(p.data, 2)
        println(o, "$(x[1,i]) $(x[2,i]) $(x[3,i]) $(x[4,i]) $(x[5,i]) $(x[6,i])")
    end
    if p.closedCycle
        println(o, "} \\closedcycle;")
    else
        println(o, "};")
    end
end
# todo there is a lot of code redundancy
function plotHelperErrorBars(o::IO, p::BarChart)
    println(o, "[")
    optionHelper(o, barMap, p, brackets=false)
    println(o, ", error bars/.cd, ")
    optionHelper(o, errorbarsMap, p.errorBars, otherText=["x dir=both", "x explicit", "y dir=both", "y explicit"])
    println(o, "]")
    println(o, "table [")
    println(o, "x error plus=ex+, x error minus=ex-, y error plus=ey+, y error minus=ey-")
    println(o, "] {")
    println(o, "x y ex+ ex- ey+ ey-")
    x = p.errorBars.data
    for i = 1:length(p.values)
        println(o, "$(p.keys[i]) $(p.values[i]) $(x[1,i]) $(x[2,i]) $(x[3,i]) $(x[4,i])")
    end
    println(o, "};")
end

function plotHelper(o::IO, p::Linear)
    print(o, "\\addplot+ ")
    if p.errorBars == nothing
        optionHelper(o, linearMap, p, brackets=true)
        println(o, "coordinates {")
        for i in 1:size(p.data,2)
            println(o, "($(p.data[1,i]), $(p.data[2,i]))")
        end
        if p.closedCycle
            println(o, "} \\closedcycle;")
        else
            println(o, "};")
        end
    else
        plotHelperErrorBars(o, p)
    end
    plotLegend(o, p.legendentry)
end

function plotHelper(o::IO, p::SmithData)
    print(o, "\\addplot+ ")
    optionHelper(o, linearMap, p, brackets=true)
    println(o, "coordinates {")
    for item in p.data
        print(o, "($(real(item)), $(imag(item))) ")
    end
    println(o, "};")
    plotLegend(o, p.legendentry)
end

function plotHelper(o::IO, p::SmithCircle)
    if p.style != nothing
        println(o, "\\draw[$(p.style)] ($(p.xc), $(p.yc)) circle ($(p.radius)cm);")
    else
        println(o, "\\draw ($(p.xc), $(p.yc)) circle ($(p.radius)cm);")
    end
end

function plotHelper(o::IO, p::Scatter)
    if size(p.data,1) == 2
        print(o, "\\addplot+[draw=none, ")
        p.onlyMarks = nothing
    elseif p.scatterClasses == nothing
        print(o, "\\addplot+[scatter, scatter src=explicit, ")
    else
        print(o, "\\addplot+[scatter, scatter src=explicit symbolic, ")
    end
    optionHelper(o, scatterMap, p)
    println(o, "] coordinates {")
    if size(p.data,1) == 2
        for i in 1:size(p.data,2)
            println(o, "($(p.data[1,i]), $(p.data[2,i]))")
        end
    else
        for i in 1:size(p.data,2)
            println(o, "($(p.data[1,i]), $(p.data[2,i])) [$(p.data[3,i])]")
        end
    end
    println(o, "};")
    plotLegend(o, p.legendentry)
end

# Specific version for Linear3 mutable struct
# Changes are addplot3 (vs addplot) and iterate over all 3 columns
function plotHelper(o::IO, p::Linear3)
    print(o, "\\addplot3+ ")
    optionHelper(o, linearMap, p, brackets=true)
    println(o, "coordinates {")
    for i = 1:size(p.data,2)
        println(o, "($(p.data[1,i]), $(p.data[2,i]), $(p.data[3,i]))")
    end
    println(o, "};")
    plotLegend(o, p.legendentry)
end

function plotHelper(o::IO, p::Node)

    axis = p.axis != nothing ? p.axis : "axis cs"

    if p.style != nothing
        println(o, "\\node at ($(axis):$(p.x), $(p.y)) [$(p.style)] {$(p.data)};")
    else
        println(o, "\\node at ($(axis):$(p.x), $(p.y)) {$(p.data)};")
    end
end

function plotHelper(o::IO, p::ErrorBars)
    print(o, "\\addplot+ [")
    optionHelper(o, errorbarsMap, p)
    print(o, ",error bars/.cd, x dir=both, x explicit, y dir=both, y explicit")
    println(o,"]")
    println(o, "coordinates {")
    for i = 1:size(p.data,2)
        println(o, "($(p.data[1,i]), $(p.data[2,i])) +=($(p.data[3,i]),$(p.data[4,i])) -=($(p.data[5,i]),$(p.data[6,i]))")
    end
    println(o, "};")
    plotLegend(o, p.legendentry)
end

function plotHelper(o::IO, p::Quiver)
    print(o, "\\addplot+ ")
    optionHelper(o, quiverMap, p, brackets=true, otherOptions=Dict("quiver"=>"{u=\\thisrow{u},v=\\thisrow{v}}"))
    println(o, "table {")
    println(o, "x y u v")
    for i = 1:size(p.data,2)
        println(o, "$(p.data[1,i]) $(p.data[2,i]) $(p.data[3,i]) $(p.data[4,i])")
    end
    println(o, "};")
    plotLegend(o, p.legendentry)
end

function plotHelper(o::IO, p::Contour)
    arg = 5
    if p.number != nothing
        arg = p.number
    elseif p.levels != nothing
        arg = p.levels
    end
    C = contours(convert(Vector{Float64}, p.xbins), convert(Vector{Float64}, p.ybins), convert(Matrix{Float64}, p.data), arg)

    print(o, "\\addplot3[contour prepared")
    if p.contour_style != nothing
        print(o, "={$(p.contour_style)}")
    elseif p.labels == false
        print(o, "={labels=false}")
    end
    if p.style != nothing
        print(o, ", $(p.style)")
    end
    print(o, "] table {")

    for c in levels(C)
        level = c.level
        for l in c.lines
            for v in l.vertices
                println(o, "$(v[1]) $(v[2]) $level")
            end
            println(o)
        end
    end
    println(o, "};")
end

function plotHelper(o::IO, p::Circle)
    if p.style != nothing
        println(o, "\\draw[$(p.style)] (axis cs:$(p.xc), $(p.yc)) circle[radius=$(p.radius)];")
    else
        println(o, "\\draw (axis cs:$(p.xc), $(p.yc)) circle[radius=$(p.radius)];")
    end
end

function plotHelper(o::IO, p::Ellipse)
    if p.style != nothing
        println(o, "\\draw[$(p.style)] (axis cs:$(p.xc), $(p.yc)) ellipse[x radius=$(p.xradius), y radius=$(p.yradius)];")
    else
        println(o, "\\draw (axis cs:$(p.xc), $(p.yc)) ellipse[x radius=$(p.xradius), y radius=$(p.yradius)];")
    end
end

function plotHelper(o::IO, p::Command)
    println(o, p.cmd*";") #PGFPlots expects commands to be terminated with a ;
end

function plotHelper(o::IO, p::Image)
    if p.zmin == p.zmax
        error("Your colorbar range limits must not be equal to each other.")
    end
    if p.style != nothing
        println(o, "\\addplot [$(p.style), point meta min=$(p.zmin), point meta max=$(p.zmax)] graphics [xmin=$(p.xmin), xmax=$(p.xmax), ymin=$(p.ymin), ymax=$(p.ymax)] {$(p.filename)};")
    else
        println(o, "\\addplot [point meta min=$(p.zmin), point meta max=$(p.zmax)] graphics [xmin=$(p.xmin), xmax=$(p.xmax), ymin=$(p.ymin), ymax=$(p.ymax)] {$(p.filename)};")
    end
end

function plotHelper(o::IO, p::Patch2D)
    print(o, "\\addplot")

    optionHelper(o, patch2DMap, p, brackets=true)
    println(o)

    m = p.patch_type == "rectangle" ? 4 : 3
    if size(p.data, 1) == 3 # include color
        println(o, "table[point meta=\\thisrow{c}] {")
        println(o, "\tx y c")
    elseif size(p.data, 1) == 2
        println(o, "table {")
        println(o, "\tx y")
    else
        error("Incorrect dimensions of matrix passed to Patch2D")
    end
    for i in 1 : size(p.data, 2)
        @printf(o, "\t %g %g", p.data[1,i], # x
                                   p.data[2,i]) # y
        if size(p.data, 1) > 2
            @printf(o, " %g\n", p.data[3,i]) # color
        else
            println(o, "")
        end

        if mod(i, m) == 0
            println(o) # extra space between patches (trainges or rectangles)
        end
    end
    println(o, "};")
    plotLegend(o, p.legendentry)
end

function plotHelper(o::IO, p::MatrixPlot)
    if p.zmin >= p.zmax
        error("Your colorbar range limits must not be equal to each other.")
    end
    if p.raster
        if p.style != nothing
            println(o, "\\addplot [$(p.style), point meta min=$(p.zmin), point meta max=$(p.zmax)] graphics[xmin=$(p.xmin), xmax=$(p.xmax), ymin=$(p.ymin), ymax=$(p.ymax)] {$(p.filename)};")
        else
            println(o, "\\addplot [point meta min=$(p.zmin), point meta max=$(p.zmax)] graphics[xmin=$(p.xmin), xmax=$(p.xmax), ymin=$(p.ymin), ymax=$(p.ymax)]{$(p.filename)};")
        end
    else
        if p.style != nothing
            println(o, "\\addplot [matrix plot*, $(p.style), point meta=explicit, point meta min=$(p.zmin), point meta max=$(p.zmax), mesh/cols=$(p.cols), mesh/rows=$(p.rows)] table[meta=data] {$(p.filename)};")
        else
            println(o, "\\addplot [matrix plot*, point meta=explicit, point meta min=$(p.zmin), point meta max=$(p.zmax), mesh/cols=$(p.cols), mesh/rows=$(p.rows)] table[meta=data] {$(p.filename)};")
        end
    end
end

# plot option string and contents; no \begin{axis} or \nextgroupplot
function plotHelper(o::IO, axis::Axis)
    optionHelper(o, axisMap, axis, brackets=true, otherText=[axisOptions(p) for p in axis.plots])
    for p in axis.plots
        plotHelper(o, p)
        try
            if p.texlabel !== nothing
                println(o, "\\label{$(p.texlabel)}")
            end
        catch
        end
    end
end

function plot(axis::Axis)
    o = IOBuffer()
    print(o, "\\begin{$(axis.axisKeyword)}")
    plotHelper(o, axis)
    println(o, "\\end{$(axis.axisKeyword)}")
    TikzPicture(String(take!(o)), options=pgfplotsoptions(), preamble=pgfplotspreamble())
end

function plot(axes::Axes)
    o = IOBuffer()

    for axis in axes
        print(o, "\\begin{$(axis.axisKeyword)}")
        plotHelper(o, axis)
        println(o, "\\end{$(axis.axisKeyword)}")
    end
    TikzPicture(String(take!(o)), options=pgfplotsoptions(), preamble=pgfplotspreamble())
end

function plot(p::GroupPlot)
    o = IOBuffer()
    style = ""
    if p.style != nothing
        style = p.style * ", "
    end
    groupStyle = ""
    if p.groupStyle != nothing
        groupStyle = p.groupStyle * ", "
    end
    println(o, "\\begin{groupplot}[$(style)group style={$(groupStyle)group size=$(p.dimensions[1]) by $(p.dimensions[2])}]")
    for a in p.axes
        print(o, "\\nextgroupplot ")
        plotHelper(o, a)
    end
    println(o, "\\end{groupplot}")
    mypreamble = pgfplotspreamble() * "\\usepgfplotslibrary{groupplots}"
    TikzPicture(String(take!(o)), options=pgfplotsoptions(), preamble=mypreamble)
end

const Plottable = Union{Plot,GroupPlot,Axis,Axes,TikzPicture}

plot(p::Plot) = plot(Axis(p))

function plot(p::Contour)
    plot(Axis(p, xmin=p.xbins[1], xmax=p.xbins[end], ymin=p.ybins[1], ymax=p.ybins[end]))
end

function plot(
    x::AbstractArray{A,1},
    y::AbstractArray{B,1};
    kwargs...,
    ) where {A<:Real,B<:Real}

    return plot(Linear(x, y; kwargs...))
end

function plot(
    x::AbstractVector{A},
    y::AbstractVector{B},
    z::AbstractVector{C};
    kwargs...,
    ) where {A<:Real,B<:Real,C<:Real}

    return plot(Linear3(x, y, z; kwargs...))
end

function Plots.Linear(f::Function, arg_range::RealRange; xbins=100, mark="none", kwargs...)
    x = range(arg_range[1], stop=arg_range[2], length=xbins)
    Linear(x, map(f, x); mark=mark, kwargs...)
end

plot(f::Function, arg_range::RealRange; kwargs...) = plot(Linear(f, arg_range; kwargs...))

plot(tkz::TikzPicture) = tkz # tikz pic doesn't need plot, here for convenience

Base.showable(::MIME"image/svg+xml", p::Plottable) = true

cleanup(p::Axis) = map(cleanup, p.plots)
cleanup(axes::Axes) = map(cleanup, axes)
cleanup(p::GroupPlot) = map(cleanup, p.axes)
cleanup(p::Plot) = nothing
cleanup(p::Circle) = nothing
cleanup(p::Ellipse) = nothing
cleanup(p::Command) = nothing
cleanup(p::Image) = rm(p.filename)
cleanup(p::MatrixPlot) = rm(p.filename)
cleanup(p::Contour) = nothing
cleanup(p::TikzPicture) = nothing

axisOptions(p::Plot) = nothing

function colormapOptions(cm::ColorMaps.GrayMap)
    if cm.invert
        return "colormap={wb}{gray(0cm)=(1); gray(1cm)=(0)}"
    else
        return "colormap={wb}{gray(0cm)=(0); gray(1cm)=(1)}"
    end
end

function colormapOptions(cm::ColorMaps.RGBArrayMap)
    o = IOBuffer()
    print(o, "colormap={mycolormap}{ ")
    n = length(cm.colors)
    for i = 1:n
        c = cm.colors[i]
        print(o, "rgb($(i-1)cm)=($(c.r),$(c.g),$(c.b)) ")
    end
    print(o, "}")
    String(take!(o))
end

function axisOptions(p::Image)
    if p.colorbar
        cmOpt = colormapOptions(p.colormap)
        if p.colorbarStyle == nothing
            return "enlargelimits = false, axis on top, $cmOpt, colorbar"
        else
            return "enlargelimits = false, axis on top, $cmOpt, colorbar, colorbar style = {$(p.colorbarStyle)}"
        end
    else
        return "enlargelimits = false, axis on top"
    end
end

function axisOptions(p::MatrixPlot)
    if p.colorbar
        cmOpt = colormapOptions(p.colormap)
        option_string = "enlargelimits = false, axis on top, $cmOpt, colorbar, xmin=$(p.xmin), xmax=$(p.xmax), ymin=$(p.ymin), ymax=$(p.ymax)"
        if p.colorbarStyle != nothing
            option_string *= ", style = {$(p.colorbarStyle)}"
        end
        return option_string
    else
        return "enlargelimits = false, axis on top"
    end
end

canPlot(p::Axis) = all(map(canPlot, p.plots))
canPlot(axes::Axes) = all(map(canPlot, axes))
canPlot(p::GroupPlot) = all(map(canPlot, p.axes))
canPlot(p::Plot) = true
canPlot(p::Circle) = true
canPlot(p::Ellipse) = true
canPlot(p::Command) = true
canPlot(p::Image) = isfile(p.filename)
canPlot(p::MatrixPlot) = isfile(p.filename)
canPlot(p::Contour) = true
canPlot(p::TikzPicture) = true

function Base.show(f::IO, a::MIME"image/svg+xml", p::Plottable)
    if canPlot(p)
        r = Base.show(f, a, plot(p))
        cleanup(p)
        r
    end
end

function save(filename::AbstractString, o::Plottable; include_preamble::Bool=true)
    _, ext = splitext(filename)
    ext = lowercase(ext)
    if ext == ".pdf"
        save(PDF(filename), plot(o))
        cleanup(o)
    elseif ext == ".svg"
        save(SVG(filename), plot(o))
        cleanup(o)
    elseif ext == ".tex"
        save(TEX(filename, include_preamble=include_preamble), plot(o))
    elseif ext == ".tikz"
        save(TIKZ(filename), plot(o))
    elseif ext == "." || ext == ""
        error("You must specify a file extension.")
    else
        error("Unsupported file extensions: $ext")
    end
end

end # module


text:
module ChemistryFeaturization

export AtomFeat, atom_data_df, build_atom_feats, make_feature_vectors, decode_feature_vector, default_nbins
include("atomfeat.jl")

export AtomGraph, normalized_laplacian, add_features!, visualize_graph
include("atomgraph.jl")

export inverse_square, exp_decay, build_graph, build_graphs_batch, read_graphs_batch
include("pmg_graphs.jl")

# Sean: add which functions/variables to be exported from your files; I've just added the ones to make the tests pass
export weave_featurize, smiles_atom_features, smiles_bond_features
include("smiles_rdkit_features.jl")

end


text:
"""
Functions to return MaterialPropertyOfX for X appropriate to type
MaterialPropertyOfX
"""

"""
interp1_bounded:
bounded linear interpolation of values xarr, yarr for value v
"""
function interp1_bounded(xarr::Vector{<:Real},yarr::Vector{<:Real} , v::Real)
    #xarr musst be sorted low to hight
    tol = 1E-8
    res = 0.0

    if v <= xarr[1]
        res =  yarr[1]
    elseif v >= xarr[end]
        res =  yarr[end]
    else
        xlow = findfirst(f -> f .> (v+tol), xarr)

        res = ((yarr[xlow] - yarr[xlow-1])/(xarr[xlow] - xarr[xlow-1])
            * (v-xarr[xlow-1]) + yarr[xlow-1])
    end

    return res
end

"""
interp1_extra:
linear extrapolated linear interpolation of values xarr, yarr for value v
"""
function interp1_extra(xarr::Vector{<:Real},yarr::Vector{<:Real} , v::Real)
    #xarr musst be sorted low to hight
    tol = 1E-8

    if v <= xarr[1]
        xlow = 2
        res = ((yarr[xlow] - yarr[xlow-1])/(xarr[xlow] - xarr[xlow-1])
        * (v-xarr[xlow-1]) + yarr[xlow-1])

    elseif v >= xarr[end]
        xlow = length(xarr)
        res = ((yarr[xlow] - yarr[xlow-1])/(xarr[xlow] - xarr[xlow-1])
        * (v-xarr[xlow-1]) + yarr[xlow-1])
    else
        xlow = findfirst(f -> f .> (v+tol), xarr)

        res = ((yarr[xlow] - yarr[xlow-1])/(xarr[xlow] - xarr[xlow-1])
            * (v-xarr[xlow-1]) + yarr[xlow-1])
    end

    return res
end

"""
In the follwowing the implementations for the different MaterialPropertyOfX are given:
"""
# import Base

Base.iterate(x::MaterialPropertyOfX) = (x, nothing)
Base.iterate(x::MaterialPropertyOfX, ::Any) = nothing
Base.isempty(x::MaterialPropertyOfX) = false
Base.length(x::MaterialPropertyOfX) = 1

"""
Sepcific vector versions of functions allow for little optimization for some types
TODO: Add Variantes
"""
function pX(p::T1, X::T2)::T2 where T2 <: Vector{<:Real} where T1 <: MaterialPropertyOfX
    return pX.(p, X)
end

function pX_bounded(p::T1, X::T2)::T2 where T2 <: Vector{<:Real} where T1 <: MaterialPropertyOfX
    return pX.(p, X)
end

function pX_bounded(p::T1, X::T2)::T2 where T2 <: Vector{<:Real} where T1 <: BoundedMaterialPropertyOfX
    return  pX_bounded.(p, X)
end

function pX(p::T1, X::T2)::T2 where T2 <: Vector{<:Real} where T1 <: ConstantMaterialPropertyOfX
    return fill( pX(p, X[1]), length(X))
end

function pX_bounded(p::T1, X::T2)::T2 where T2 <: Vector{<:Real} where T1 <: ConstantMaterialPropertyOfX
    return pX(p, X)
end

"""
MaterialPropertyOfX types which have no bounds property should be able
to bound anyhow
"""
function pX_bounded(p::T1, X::T2)::T2 where T2 <: Real where T1 <:MaterialPropertyOfX
    p_val = pX(p, X)

    return p_val
end

function pX_bounded(p::T1, X::T2)::T2 where T2 <: Real where T1 <:BoundedMaterialPropertyOfX
    p_val = pX(p,X)

    if !isnothing(p.bounds)
        if p_val < p.bounds[1]
            p_val = p.bounds[1]
        end

        if p_val > p.bounds[2]
            p_val = p.bounds[2]
        end
    end

    return p_val
end

function pX(p::T1, X::T2)::T2 where T2 <: Real where T1 <:MaterialPropertyOfX
    error("This function should not be called ...")
    return X
end

function pX(p::ConstantMaterialPropertyOfX{T}, X::T)::T where T <: Real
    return p.val
end

function pX(p::TabulatedMaterialPropertyOfX, X::T)::T where T <: Real
    index = floor(Int64, ((X - p.range.start) + p.rangestep + 1E-16) / p.rangestep)

    return @boundscheck p.val_arr[index]
end

function pX(p::LinearInterpMaterialPropertyOfX{T}, X::T)::T where T <: Real
    if p.bounded
        return interp1_bounded.(p.x,p.val, X)
    else
        return interp1_extra.(p.x,p.val, X)
    end
end

function pX(p::LinearMaterialPropertyOfX{T}, X::T)::T where T <: Real
    return p.m * X + p.y1
end

function pX(p::PolynomialMaterialPropertyOfX{T1, T2}, X::T1)::T1 where {T1 <: Real, T2 <: Integer}
    p_val = 0.0
    for i = 1:1:p.n
        p_val += p.a[i] * X^(i-1)
    end
    return p_val
end

function pX(p::ReciprocalMaterialPropertyOfX{T}, X::T)::T where T <: Real
    return p.m * 1.0 /X + p.y1
end

function pX(p::ExponentialMaterialPropertyOfX{T}, X::T)::T where T <: Real
    return p.a * exp(p.b * x) + p.y1
end

function pX(p::LogarithmicMaterialPropertyOfX{T}, X::T)::T where T <: Real
    return p.a * ln.(X) + p.y1
end

"""
num_integrate_pX:
in general use numerical integration
(uses QuadGK.jl)
"""
function num_integrate_pX(p::MaterialPropertyOfX,
        X_from::T, X_to::T)::T where T <: Real

    H, err = quadgk.( x -> pX(p, x), X_from, X_to)

    if err > 1e-3
        error("Error integration of x is high")
    end
    return H
end

"""
integrate_pX(p::PolynomialMaterialPropertyOfX, X_from::Real,
        X_to::Real):

Analytical integration of polynomial material properties
"""
function integrate_pX(p::PolynomialMaterialPropertyOfX{T, <:Integer}, X_from::T,
        X_to::T)::T where T <: Real

    local  p_val_from::T
    local p_val_to::T
    p_val_from = 0.0
    p_val_to = 0.0

    for i = 1:1:p.n
        p_val_from += 1.0/(i) * p.a[i] * X_from^(i)
        p_val_to += 1.0/(i) * p.a[i] * X_to^(i)
    end

    return (p_val_to - p_val_from)
end

"""
integrate_pX(p::ConstantMaterialPropertyOfX, X_from::Real,
        X_to::Real):

Analytical integration of constant material properties
"""
function integrate_pX(p::ConstantMaterialPropertyOfX{T}, X_from::T,
        X_to::T)::T where T <: Real

    return (p.val * (X_to - X_from))
end

"""
integrate_pX(p::TabulatedMaterialPropertyOfX, X_from::Real,
        X_to::Real):

Analytical integration of Tabulated material properties
"""
function integrate_pX(p::TabulatedMaterialPropertyOfX, X_from::T,
        X_to::T)::T where T <: Real

    ifrom = floor(Int64, (X_from - (p.range.start - p.rangestep)) / p.rangestep)
    ito = floor(Int64, (X_to - (p.range.start - p.rangestep)) / p.rangestep)

    int = 0.0
    if ito > ifrom
        @boundscheck for i = ifrom:1:ito
            int += p.val_arr[i] * p.rangestep
        end
    else
        @boundscheck for i = ito:-1:ifrom
            int -= p.val_arr[i] * p.rangestep
        end
    end

    return int
end

"""
integrate_pX(p::MaterialPropertyOfX, X_from::Real,
        X_to::Real):

Fall back to numerical integration until specific version is available!

TODO: Impelement integration functions for other MaterialPropertyOfX
"""
function integrate_pX(p::MaterialPropertyOfX, X_from::T,
    X_to::T)::T where T <: Real

    return num_integrate_pX(p, X_from, X_to)
end


text:
### QuantumBit.jl
# Convenience functionality for working with qubit-based systems
###

"""
    qubit([name="Qubit"::QName])

Create qubit in computational basis.

**Note:** Rows are permuted from usual matrix representation.
This is so the 0 state can represent the ground state of energy
when constructing a qubit Hamiltonian.

**Example:**
```julia
q = qubit()
q('z') * groundvec(q) == sparse([QComp(-1); 0])  # True
H = (ħ*ωq/2)*q('z') # Correct energy structure
```

### Default ops for QFactor:
  - 'i' : identity operator
  - 'd' : lowering operator (``σ_-``)
  - 'u' : raising operator  (``σ_+``)
  - 'x' : Pauli x operator  (``σ_x``)
  - 'y' : Pauli y operator  (``σ_y``)
  - 'z' : Pauli z operator  (``σ_z``)

"""
function qubit(name="Qubit")
    q = osc(2, name)
    # Replace the number operator with a rescaling to get Pauli z
    merge!(q.ops, Dict('z'=> 2 .* q('n') .- q('i')))
    delete!(q.ops,'n')
    q
end


text:
# This file should contain site-specific commands to be executed on Julia startup;
# Users may store their own personal commands in `~/.julia/config/startup.jl`.

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
# Events

function parseevent(window::Window, ::Val{SDL.KEYDOWN}, data::Vector{UInt8})
    return geteventdata(:key_down, data, :type => UInt32, :timestamp => UInt32, :window_id => UInt32,
                        :state => UInt8, :repeat => UInt8, :empty => UInt16, :scancode => UInt32,
                        :keycode => UInt32, :mod => UInt16)
end

function parseevent(window::Window, ::Val{SDL.KEYUP}, data::Vector{UInt8})
    return geteventdata(:key_up, data, :type => UInt32, :timestamp => UInt32, :window_id => UInt32,
                        :state => UInt8, :repeat => UInt8, :empty => UInt16, :scancode => UInt32,
                        :keycode => UInt32, :mod => UInt16)
end

"""
    iskey(e::Event, key::String; repeat=false)

Check if the event `e` is a key event with scancode corresponding to
`key`.  If the keyword argument `repeat` is true, this function also
returns true if the event was triggered by a key repeat (i. e., a held
key).

"""
iskey(e::Event, key::String; repeat=false) = e.scancode == getscancode(key) && (repeat || convert(Bool, e.repeat) == repeat)

"""
    iskeycode(e::Event, key::String; repeat=false)

Check if the event `e` is a key event with keycode corresponding to
`key`.  If the keyword argument `repeat` is true, this function also
returns true if the event was triggered by a key repeat (i. e., a held
key).

"""
iskeycode(e::Event, key::String; repeat=false) = e.keycode == getkeycode(key) && (repeat || convert(Bool, e.repeat) == repeat)

getkeyname(e::Event) = unsafe_string(SDL.GetScancodeName(convert(Int32, e.scancode)))
getkeycodename(e::Event) = unsafe_string(SDL.GetKeycodeName(convert(Int32, e.keycode)))
getscancode(key::String) = SDL.GetScancodeFromName(key)
getkeycode(key::String) = SDL.GetKeycodeFromName(key)

# Keyboard state

struct KeyboardState
    state::Base.ReinterpretArray{Bool, 1, UInt8, Array{UInt8, 1}}
    KeyboardState() = new(reinterpret(Bool, unsafe_wrap(Array{UInt8,1}, SDL.GetKeyboardState(C_NULL), 285)))
end

ispressed(keyboard::KeyboardState, key::String) = keyboard.state[getscancode(key) + 1]


text:
function ShowFGENWindow(fgen, fgen_conf, rev_state_dict, refresh_cnt, base)
	CImGui.Begin("GFG3015 Function Generator")
	# DRAW DISPLAY INFO

	CImGui.Text("Func:"), CImGui.SameLine(), CImGui.Text(fgen_conf.crt_func)
	CImGui.Text("Freq:")
	CImGui.SameLine(), CImGui.Text(fgen_conf.freq)
	CImGui.SameLine(), CImGui.Text("Hz")
	CImGui.Text("Ampl:")
	CImGui.SameLine(), CImGui.Text(fgen_conf.amplit)
	CImGui.SameLine(), CImGui.Text(fgen_conf.amplit_unit)
	CImGui.Text("Offs:")
	CImGui.SameLine(), CImGui.Text(fgen_conf.offs)
	CImGui.SameLine(), CImGui.Text(fgen_conf.amplit_unit)
	CImGui.Text("Duty:")
	CImGui.SameLine(), CImGui.Text(fgen_conf.duty)
	CImGui.SameLine(), CImGui.Text("%")

	CImGui.Text("EDIT ")
	@c CImGui.InputDouble("", &fgen_conf.crt_value, 0.01, 1.0, "%.4f")

	CImGui.Button(" FUNC ") && begin
		# toggle between the three functions
		toggle_func!(fgen, fgen_conf)
	end
	CImGui.SameLine(),CImGui.Button(" AMPL ") && begin
		set_amplit(fgen, fgen_conf.crt_chan, fgen_conf.crt_value)
	end
	CImGui.SameLine(),CImGui.Button(" DUTY ") && begin
		set_duty(fgen, fgen_conf.crt_chan, fgen_conf.crt_value)
	end

	CImGui.Button(" FREQ ") && begin
		# toggle between the three functions
		set_freq(fgen, fgen_conf.crt_chan, fgen_conf.crt_value)
	end
	CImGui.SameLine(),CImGui.Button("OFFSET") && begin
		set_offs(fgen, fgen_conf.crt_chan, fgen_conf.crt_value)
	end
	CImGui.SameLine(),CImGui.Button("MOD ON") && begin
		# don't need
	end

	CImGui.SameLine()

	# Draw REFRESH
	CImGui.SameLine()
	fgen_conf.refresh_btn, pressed = draw_toggle_button(fgen_conf.refresh_btn)
	fgen_conf.refresh_btn.state && @async update_fgen_conf!(fgen_conf, fgen, refresh_cnt, base)

	CImGui.End()
	return fgen_conf
end

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
const TermTuple = Tuple{Vararg{<:Terms.Term,N}} where N

"""
    nc(t::Type{T}) where T <: Terms.Term

Return the number of columns generated by a term of type `t`.  For `Continuous`,
`Intercept`, and `FunctionTerm`, this is 1.  For `Categorical`, it's the size of
the contrasts matrix.  For an `Interaction` term, it's the product of the number
of columns generated by the children terms.  For an `Eval` term, an error is
thrown (since it can't be determined).

For a custom term type, you need to provide an `nc` method as well as a
`rowval` method.

"""
nc(t::Type{Terms.Eval}) =
    throw(ArgumentError("can't compute number of columns generated by un-typed " *
                        "Eval term $t. Did you forget to `set_schema`?"))
nc(::Type{<:Terms.Term}) = 0
nc(::Type{Terms.Intercept}) = 1
nc(::Type{Terms.Continuous}) = 1
nc(::Type{Terms.Categorical{N}}) where N = N
nc(::Type{Terms.Interaction{T}}) where T = mapreduce(nc, *, T.parameters)
# TODO: is this right? should we handle n-ary functions?
nc(::Type{Terms.FunctionTerm{F1,F2}}) where {F1,F2} = 1
nc(::Type{T}) where T<:TermTuple where N = mapreduce(nc, +, T.parameters)

"""
    (term::<:Terms.Term)(data::NamedTuple)

Calling a term with a named tuple should realize the data corresponding to that
named tuple as a numeric vector

"""
(t::Terms.Intercept)(data::NamedTuple) = 1
(t::Terms.Continuous)(data::NamedTuple) = data[t.name]
(t::Terms.Categorical)(data::NamedTuple) =
    t.contrasts.matrix[t.invindex[data[t.name]], :]
(t::Terms.FunctionTerm)(data::NamedTuple) = t.f(data)
# TODO: does this really need to be @generated??
@generated function (t::Terms.Interaction)(data::NamedTuple)
    nterms = length(t.parameters[1].parameters)
    Expr(:call, :kron, [:(t.terms[$ti](data)) for ti in 1:nterms]...)
end

# TODO: add methods for named tuple of vectors (e.g. Data.Table): I think it's
# just the interaction and function terms that need special cases for this.

################################################################################
# Now we have a Formula with a vector of Terms that are typed.  We just need to
# have a generated function that generates the code for filling in a single
# model matrix row (similar to what AnonFactory does).  This requires knowing
# the number of columns generated by each term (`nc`).

function modelmatrix(source, f::Formula)
    f.schema_set || set_schema!(f, summarize(source, f))
    modelmatrix(source, f.term)
end

function modelmatrix(source, term::Terms.FormulaTerm)
    # get a tuple of response, matrix, and other terms
    terms = Terms.flatten_formula(term)
    modelmatrix(source, terms...)
end

nrows(t::Data.Table) = length(t[1])
nrows(r::Data.Rows) = length(r)

init_matrix(source, term::T) where T<:Union{Terms.Term, TermTuple} =
    Array{Float64}(nrows(source), nc(T))

function modelmatrix(source, terms::T) where T<:Tuple{Vararg{<:Union{Terms.Term, TermTuple},N}} where N
    mats = tuple(init_matrix(source, t) for t in terms...)
    mat_terms = zip(mats, terms)

    iter = Data.rows(source)

    for (i, data) in enumerate(iter)
        for (mat, term) in mat_terms
            modelmatrow!(view(mat, i, :), data, term)
        end
    end

    mats
end

"""
    modelmatrow!(row::AbstractVector, data::NamedTuple, terms::Tuple)
    modelmatrow(data::NamedTuple, terms::Tuple)

Fill one model matrix row based on a data named tuple and terms.

"""
@generated function modelmatrow!(row::AbstractVector, data::NamedTuple, terms::Tuple)
    Ts = collect(terms.parameters)
    ncols = mapreduce(nc, +, Ts)
    @debug "$ncols columns"
    @debug "Terms:"
    for T in Ts
        @debug T
    end

    func_body = Expr(:block)
    ci = 0
    for ti in eachindex(Ts)
        starti, ci = ci+1, ci+nc(Ts[ti])
        push!(func_body.args, :(@inbounds row[$starti:$ci] = terms[$ti](data)))
    end
    push!(func_body.args, :(row))
    @debug func_body
    return func_body
end

@generated function modelmatrow(data::NamedTuple, terms::Tuple)
    Ts = collect(terms.parameters)
    cols = mapreduce(nc, +, Ts)
    :(modelmatrow!(Array{Float64}($cols), data, terms))
end

modelmatrow!(row::AbstractVector, data::NamedTuple, term::T) where T<:Terms.Term =
    modelmatrow!(row, data, (term,))


text:
function triple03()
#
#  test polynomial suggested by Goedecker
#
    p = triple(18, 10, 16)
    z = [[0.9;1;1.1]  [18; 10; 16]]

    p, PolyZeros(z)
end


text:
module ComplexVisualGtk

using Cairo
import Gtk
import ComplexVisual
@ComplexVisual.import_huge

include("monkeypatch.jl")

"""macro for importing the *huge* set of symbols."""
macro import_huge()
    :(
        using ComplexVisualGtk:
            cvg_wait_for_destroy, cvg_create_win_for_canvas,
            cvg_handler_for_canvas, cvg_visualize, cvg_close
    )
end

"""
wait for destroy signal on widget.
"""
function cvg_wait_for_destroy(widget)
    cond = Gtk.Condition()
    Gtk.signal_connect(widget, :destroy) do widget
        Gtk.notify(cond)
    end
    Gtk.wait(cond)
end

"""
create GtkWindow with GtkDrawArea as single child.
"""
function cvg_create_win_with_draw_canvas(
        title; width=-1, height=-1, resizeable=false, toplevel=true)

    win = Gtk.GtkWindow(title, width, height, resizeable, toplevel)
    canvas = Gtk.GtkCanvas(width, height)
    push!(win, canvas)
    return (win, canvas)
end

"""
Callable handler to copy a 2DCanvas to a Gtk-widget (typically a GtkDrawArea)
"""
struct CVG_CopyCanvas2Widget{canvasT<:CV_2DCanvas}
    canvas :: canvasT
end
function (cc2widget::CVG_CopyCanvas2Widget)(widget)
    try
        ctx = Gtk.getgc(widget)
        set_operator(ctx, Cairo.OPERATOR_SOURCE)
        set_source(ctx, cc2widget.canvas.surface)
        paint(ctx)
    catch err
        @warn("Error in (Gtk) draw CVG_CopyCanvas2Widget",
              exception=(err, catch_backtrace()))
    end
    return nothing
end

"""
create GtkWindows vor 2DCanvas and attach as draw-method CVG_CopyCanvas2Widget
"""
function cvg_create_win_for_canvas(canvas::T,
        title; resizeable=false, toplevel=true) where {T<:CV_2DCanvas}

    win, gtk_canvas = cvg_create_win_with_draw_canvas(
        title; width=canvas.pixel_width, height=canvas.pixel_height,
        resizeable=resizeable, toplevel=toplevel)
    cc2widget = CVG_CopyCanvas2Widget{T}(canvas)
    gtk_canvas.draw = widget -> cc2widget(widget)
    Gtk.showall(win)
    return win, gtk_canvas
end

"""
Stores a Point.
"""
mutable struct CVG_Point{T<:Number}
    x   :: T
    y   :: T
end

"""
A callable struct which can be used as a gtk-eventhandler.
If called (with `event` data) `update_action_cb` with the mouse
coordinates of the event. The seen coordinate is saved in
`last_action_point`. This is used as a cache to only propagate
events with a different position.
After the `update_action_cb` call the draw-method of the `gtk_canvas`
is called.
"""
struct CVG_CanvasActionHandler{winT, gtkcanT, actionCB}
    window              :: winT
    gtk_canvas          :: gtkcanT
    last_action_point   :: CVG_Point{Int32}
    update_action_cb    :: actionCB
end

function (cah::CVG_CanvasActionHandler)(widget, event)
    gx, gy = round(Int32, event.x), round(Int32, event.y)  # Pos in gtk_canvas
    last = cah.last_action_point
    if gx != last.x  ||  gy != last.y || !(event isa Gtk.GdkEventMotion)
        last.x, last.y = gx, gy
        cah.update_action_cb(gx, gy)

        canvas = cah.gtk_canvas
        canvas.draw(canvas)
        Gtk.reveal(canvas, false)
    end
    return true
end

function cvg_close(cah::CVG_CanvasActionHandler)
    close_window(cah.window)
    return nothing
end

"""
constructs window for a `CV_2DCanvas`.
"""
function cvg_handler_for_canvas(canvas::canT, title;
        update_action_cb=(pixel_x, pixel_y)->nothing,
        update_state_cb=(pixel_x, pixel_y)->nothing) where {canT<:CV_2DCanvas}
    win, gtk_canvas = cvg_create_win_for_canvas(canvas, title)
    handler = CVG_CanvasActionHandler(
        win, gtk_canvas,
        CVG_Point(Int32(-1), Int32(-1)),
        update_action_cb)
    handler_state = (widget, event) -> begin
        gx, gy = round(Int32, event.x), round(Int32, event.y)
        update_state_cb(gx, gy)
        handler.last_action_point.x = -1
        handler(widget, event)
        return nothing
    end
    gtk_canvas.mouse.button1motion = (widget,event) -> handler(widget,event)
    gtk_canvas.mouse.button1press = (widget,event) -> handler(widget,event)
    gtk_canvas.mouse.button3press = handler_state
    return handler
end

function cvg_visualize(scene::CV_2DScene;
        title::AbstractString="Complex Visual")
    handler = cvg_handler_for_canvas(
        cv_get_can_layout(scene), title,
        update_action_cb = cv_get_actionpixel_update(scene),
        update_state_cb = cv_get_statepixel_update(scene))
    return handler
end

end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
function percolate_vertices(graph::AbstractGraph, p)::AbstractGraph
    return percolate_vertices(graph, p, [0])
end

"""
    percolate_vertices(graph, p, exclude)

Perform vertex percolation on graph with probability p, excluding specified vertices.
"""
function percolate_vertices(graph::AbstractGraph, p, exclude::Array{Int64})::AbstractGraph
    perc_graph = deepcopy(graph)

    for vertex in vertices(graph)
        # if !(vertex in exclude)
            if rand(Float64) <= p
                rem_vertex!(perc_graph, vertex)
            end
        # end
    end

    return perc_graph
end

# TODO: Modify this so it removes edges in both directions
function percolate_edges(graph::AbstractGraph, p)::AbstractGraph
    perc_graph = deepcopy(graph)

    for edge in edges(graph)
        if rand(Float64) <= p
            rem_edge!(perc_graph, edge)
        end
    end

    return perc_graph
end

function percolate_edges(network::QNetwork, p)::QNetwork
    perc_network = deepcopy(network)

    for channel in perc_network.channels
        if rand(Float64) <= p
            channel.active = false
        end
    end
    refresh_graph!(perc_network)
    return perc_network
end

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
"""
This module defines extensions to the `Base.IO` interface to support:
 - an `unread!` function for pushing excess bytes back into a stream,
 - `startwrite`, `closewrite`, `startread` and `closeread` for streams
    with transactional semantics.
"""

module IOExtras

export IOError, isioerror,
       unread!,
       startwrite, closewrite, startread, closeread,
       tcpsocket, localport, peerport

"""
    isioerror(exception)

Is `exception` caused by a possibly recoverable IO error.
"""

isioerror(e) = false
isioerror(::Base.EOFError) = true
isioerror(::Base.UVError) = true
isioerror(e::ArgumentError) = e.msg == "stream is closed or unusable"

"""
The request terminated with due to an IO-related error.

Fields:
 - `e`, the error.
"""

struct IOError <: Exception
    e
    message
end

Base.show(io::IO, e::IOError) = print(io, "IOError(", e.e, " ", e.message, ")\n")

"""
    unread!(::IO, bytes)

Push bytes back into a connection (to be returned by the next read).
"""

function unread!(io::IOBuffer, bytes)
    l = length(bytes)
    if l == 0
        return
    end

    @assert bytes == io.data[io.ptr - l:io.ptr-1]

    if io.seekable
        seek(io, io.ptr - (l + 1))
        return
    end

    println("WARNING: Can't unread! non-seekable IOBuffer")
    println("         Discarding $(length(bytes)) bytes!")
    @assert false
    return
end

function unread!(io, bytes)
    if length(bytes) == 0
        return
    end
    println("WARNING: No unread! method for $(typeof(io))!")
    println("         Discarding $(length(bytes)) bytes!")
    return
end

start_close_read_write_doc = """
    startwrite(::IO)
    closewrite(::IO)
    startread(::IO)
    closeread(::IO)

Signal start/end of write or read operations.
"""

@doc start_close_read_write_doc -> startwrite(io) = nothing
@doc start_close_read_write_doc -> closewrite(io) = nothing
@doc start_close_read_write_doc -> startread(io) = nothing
@doc start_close_read_write_doc -> closeread(io) = nothing

using MbedTLS.SSLContext
tcpsocket(io::SSLContext)::TCPSocket = io.bio
tcpsocket(io::TCPSocket)::TCPSocket = io

localport(io) = try !isopen(tcpsocket(io)) ? 0 :
                    VERSION > v"0.7.0-DEV" ?
                    getsockname(tcpsocket(io))[2] :
                    Base._sockname(tcpsocket(io), true)[2]
                catch
                    0
                end

peerport(io) = try !isopen(tcpsocket(io)) ? 0 :
                  VERSION > v"0.7.0-DEV" ?
                  getpeername(tcpsocket(io))[2] :
                  Base._sockname(tcpsocket(io), false)[2]
               catch
                   0
               end

end


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the `EnhancedTestSet` type.

This set of unit tests checks the behavior of `EnhancedTestSet` for passing tests.
"""

# --- Imports

# Standard library
using Test

# External packages
using Suppressor

# Local modules
using TestTools.jltest

# --- Tests

@testset EnhancedTestSet "EnhancedTestSet: Check output dots" begin
    output = @capture_out begin
        @testset EnhancedTestSet "top-level tests" begin
            @testset "2nd-level tests 1" begin
                @test true
                @test 1 == 1
            end
            @testset "2nd-level tests 2" begin
                @test true
                @test 1 == 1
            end
        end
    end

    @test output == "...."
end


text:
using ChebyshevPicardIteration
using DifferentialEquations
using Plots

# First Test
ϵ = 0.001
f(t,y) = cos(t .+ ϵ*y)

t0 = 0
tf = 256*π
y0 = 1

# Integrate with Chebyshev-Picard Iteration
(ts,ys) = Integrate(f, y0, (t0, tf), 200, 200, 1e-20, 1000)

# Integrate with DifferentialEquations.jl
prob = ODEProblem((y,p,t) -> f(t,y), y0, (t0, tf))
sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8)

# Analytic Solution
γ = ϵ/(1 + sqrt(1 - ϵ^2))
α = 2*ϵ/(1 - ϵ + sqrt(1 - ϵ^2))
β = tan(ϵ*y0/2)/(1 + α)
ϕ(t) = 0.5*(1 - γ*ϵ)*t
σ(t) = α*(sin(ϕ(t)) + β*cos(ϕ(t)))
y(t) = -γ*t + (2/ϵ)*atan((β + σ(t)*cos(ϕ(t)))/(1 + σ(t)*sin(ϕ(t))))

plot(sol.t,sol.u - map(y,sol.t))
plot!(ts,ys .- map(y,ts))

# Second Test
#g = 9.81
#k = 0.1
#m = 0.2
#f2(t,y) = [y[2], -k*y[1]/m + g]

#y0 = [5, 0]

# Integrate with Chebyshev-Picard Iteration
#(ts,ys) = Integrate(f2, y0, (t0, tf), 80, 80, 1e-14, 1000)

# Integrate with DifferentialEquations.jl
#prob = ODEProblem((y,p,t) -> f2(t,y), y0, (t0, tf))
#sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8)

#plot(sol)
#plot!(ts,ys)


text:
"""
    InfExtendedTime{T} <: TimeType

The type `T` extended with positive and negative infinity.
"""
struct InfExtendedTime{T<:TimeType} <: TimeType
    flag :: InfFlag
    finitevalue :: T

    InfExtendedTime{T}(x::T) where {T<:TimeType} = new{T}(FINITE, x)
    InfExtendedTime{T}(x::Infinite) where {T<:TimeType} = new{T}(x==PosInf ? POSINF : NEGINF)
end

function Base.getproperty(x::InfExtendedTime, s::Symbol)
    if s === :val
        return x.flag != FINITE ? (x.flag == NEGINF ? -∞ : ∞) : x.finitevalue
    elseif s === :instant
        return x.flag != FINITE ? (x.flag == NEGINF ? -∞ : ∞) : x.finitevalue.instant
    else
        return getfield(x, s)
    end
end

InfExtendedTime{T}(x::TimeType) where {T<:TimeType} = InfExtendedTime{T}(convert(T, x))
InfExtendedTime{T}(x::InfExtendedTime) where {T<:TimeType} = InfExtendedTime{T}(convert(T, x.finitevalue))

"""
    InfExtendedTime(T)

Return `InfExtendedTime{T}` for any TimeType `T`.
"""
InfExtendedTime(::Type{T}) where {T<:TimeType} = InfExtendedTime{T}

"""
    InfExtendedTime(x)

Converts `x` to `InfExtendedTime(typeof(x))`.
"""
InfExtendedTime(x::T) where {T<:TimeType} = InfExtendedTime{T}(x)

Utils.isposinf(x::InfExtendedTime) = x.flag == POSINF
Utils.isneginf(x::InfExtendedTime) = x.flag == NEGINF


text:
module Physics

using Colors
using DataStructures: CircularBuffer, OrderedSet, SortedSet
using LinearAlgebra
using StaticArrays

using Gloria: AbstractObject, AbstractLayer, AbstractScene, RenderTask,
    add!, kill!, populate!
using ..Shapes: Shapes, AbstractShape, BoundingBox, Vertex,
    edges, rotate, trace, translate, vertices

export Physical
    # collide!, timeintersects

import Gloria: after_update!, before_update!, render!, update!, visible
import ..Shapes: inside, intersects

const flatten = Iterators.flatten

const COLLISION_TIME_OFFSET = 1e-5
const MAX_COLLISION_DEPTH = 500

abstract type AbstractPhysicalState{D,R,T} end

mutable struct PhysicalState{D,T,R} <: AbstractPhysicalState{D,T,R}
    t::T
    static::Bool
    resting::Bool
    mass::SMatrix{D,D,T}
    angularmass::SMatrix{R,R,T}
    position::SVector{D,T}
    angle::SVector{R,T}
    velocity::SVector{D,T}
    angularvelocity::SVector{R,T}
    impulse::SVector{D,T}
    angularimpulse::SVector{R,T}

    function PhysicalState{D,T}(
        t::T,
        static::Bool,
        resting::Bool,
        mass::SMatrix{D,D,T},
        angularmass::SMatrix{R,R,T},
        position::SVector{D,T},
        angle::SVector{R,T},
        velocity::SVector{D,T},
        angularvelocity::SVector{R,T},
        impulse::SVector{D,T},
        angularimpulse::SVector{R,T}
    ) where {D,T,R}
        D*(D-1)/2 != R && error("Angular and spatial dimensions do not match: $D*($D-1)/2 != $R")
        return new{D,T,R}(t, static, resting, mass, angularmass, position, angle, velocity, angularvelocity, impulse, angularimpulse)
    end
end
function PhysicalState(t::T, static::Bool, resting::Bool, mass::SMatrix{D,D,T}, angularmass::SMatrix{R,R,T}, position::SVector{D,T}, angle::SVector{R,T}, velocity::SVector{D,T}, angularvelocity::SVector{R,T}) where {D,T,R}
    return PhysicalState{D,T}(t, static, resting, mass, angularmass, position, angle, velocity, angularvelocity, SVector{D}(fill(zero(T), D)), SVector{R}(fill(zero(T), R)))
end

function PhysicalState(t::T, static::Bool, resting::Bool, mass::T, angularmass::T, position::SVector{D,T}, angle::SVector{R,T}, velocity::SVector{D,T}, angularvelocity::SVector{R,T}) where {D,T,R}
    return PhysicalState(t, static, resting, SMatrix{D,D}(Matrix(mass*I, D, D)), SMatrix{R,R}(Matrix(angularmass*I, R, R)), position, angle, velocity, angularvelocity)
end

function PhysicalState(t::T, static::Bool, mass::T, angularmass::T, position::SVector{D,T}, angle::SVector{R,T}, velocity::SVector{D,T}, angularvelocity::SVector{R,T}) where {D,T,R}
    return PhysicalState(t, static, false, SMatrix{D,D}(Matrix(mass*I, D, D)), SMatrix{R,R}(Matrix(angularmass*I, R, R)), position, angle, velocity, angularvelocity)
end

Base.copy(state::PhysicalState{D,T}) where {D,T} = PhysicalState{D,T}(
    state.t,
    state.static,
    state.resting,
    state.mass,
    state.angularmass,
    state.position,
    state.angle,
    state.velocity,
    state.angularvelocity,
    state.impulse,
    state.angularimpulse)

abstract type AbstractPhysicalObject <: AbstractObject end

"""
    Physical{WRAPPED_TYPE, SHAPE_TYPE, PHYSICAL_STATE_TYPE}
"""
mutable struct Physical{T1<:AbstractObject, T2<:AbstractShape, T3<:AbstractPhysicalState} <: AbstractPhysicalObject
    wrapped::T1
    shape::T2
    state_history::CircularBuffer{T3}
end
function Physical(wrapped, shape, state::PhysicalState)
    self = Physical(wrapped, shape, CircularBuffer{PhysicalState}(3))
    fill!(self.state_history, state)
    return self
end
Physical(wrapped, shape, static, mass, angularmass, position, angle, velocity, angularvelocity) =
    Physical(wrapped, shape, PhysicalState(0., static, mass, angularmass, position, angle, velocity, angularvelocity))

function Physical(
    wrapped,
    shape;
    D=2,
    R=floor(Int, D*(D-1)/2),
    static=false,
    mass=1.,
    angularmass=1.,
    position=SVector{D}(fill(0., D)),
    angle=SVector{R}(fill(0., R)),
    velocity=SVector{D}(fill(0., D)),
    angularvelocity=SVector{R}(fill(0., R))
)
    return Physical(wrapped, shape, static, mass, angularmass, position, angle, velocity, angularvelocity)
end

@inline currentstate(obj::Physical) = Core.getfield(obj, :state_history)[end]
@inline previousstate(obj::Physical) = Core.getfield(obj, :state_history)[end-1]
@inline previousstate(obj::Physical, n::Integer) = Core.getfield(obj, :state_history)[end-n]

@inline Base.getproperty(obj::Physical, v::Symbol) = _getproperty(obj, Val(v))
@inline _getproperty(obj::Physical, ::Val{V}) where V = Core.getfield(obj.wrapped, V)
for v in [:(:wrapped), :(:shape), :(:state_history)]
    @eval @inline _getproperty(obj::Physical, ::Val{$v}) = Core.getfield(obj, $v)
end
# for v in [:(:static), :(:m), :(:I), :(:x), :(:y), :(:θ), :(:vx), :(:vy), :(:ω)]
#     @eval @inline _getproperty(obj::Physical, ::Val{$v}) = Core.getfield(currentstate(obj), $v)
# end

@inline Base.setproperty!(obj::Physical, v::Symbol, x) = _setproperty!(obj, Val(v), x)
@inline _setproperty!(obj::Physical, ::Val{V}, x) where V = Core.setfield!(obj.wrapped, V, x)
for v in [:(:wrapped), :(:shape), :(:state_history)]
    @eval @inline _setproperty!(obj::Physical, ::Val{$v}, x) = Core.setfield!(obj, $v, x)
end
# for v in [:(:static)]
#     @eval @inline _setproperty!(obj::Physical, ::Val{$v}, x) = Core.setfield!(currentstate(obj), $v, x)
# end
# for v in [:(:m), :(:I), :(:x), :(:y), :(:θ), :(:vx), :(:vy), :(:ω)]
#     @eval @inline _setproperty!(obj::Physical, ::Val{$v}, x) = Core.setfield!(currentstate(obj), $v, convert(Float64, x))
# end

tocoordinates(obj::Physical, p::SVector) = p |> translate(-position(obj)...) |> rotate(-angle(obj)...)

shape(obj::Physical) = obj.shape |> rotate(angle(obj)...) |> translate(position(obj)...)
isstatic(obj::Physical) = currentstate(obj).static

for name in [:mass, :angularmass, :position, :angle, :velocity, :angularvelocity, :impulse, :angularimpulse]
    @eval function $name(obj::Physical)
        return currentstate(obj).$name
    end

    @eval function $(Symbol("set$(name)!"))(obj::Physical, val)
        any(isnan, val) && error("value contains NaN: $val")
        currentstate(obj).$name = val
        return obj
    end
end

function translate!(obj::Physical, val::SVector)
    any(isnan, val) && error("value contains NaN: $val")
    setposition!(obj, position(obj) + val)
end

function rotate!(obj::Physical, val::SVector)
    any(isnan, val) && error("value contains NaN: $val")
    setangle!(obj, angle(obj) + val)
end

function boost!(obj::Physical, val::SVector)
    any(isnan, val) && error("value contains NaN: $val")
    setvelocity!(obj, velocity(obj) + val)
end

function angularboost!(obj::Physical, val::SVector)
    any(isnan, val) && error("value contains NaN: $val")
    setangularvelocity!(obj, angularvelocity(obj) + val)
end

function resetmotion!(obj::Physical)
    setvelocity!(obj, zero(velocity(obj)))
    setangularvelocity!(obj, zero(angularvelocity(obj)))
    return obj
end

function applyimpulse!(obj::Physical, j::SVector, p::SVector)
    r = p - position(obj)
    setimpulse!(obj, impulse(obj) + j)
    setangularimpulse!(obj, angularimpulse(obj) + SVector((r × j)*180/pi))
    return obj
end

function resetimpulse!(obj::Physical)
    setimpulse!(obj, zero(impulse(obj)))
    setangularimpulse!(obj, zero(angularimpulse(obj)))
    return obj
end

momentum(obj::Physical) = mass(obj)*velocity(obj)
angularmomentum(obj::Physical) = angularmass(obj)*angularvelocity(obj)
kineticenergy(obj::Physical) = 0.5*velocity(obj)'mass(obj)*velocity(obj)
angularkineticenergy(obj::Physical) = 0.5*angularvelocity(obj)'angularmass(obj)*angularvelocity(obj)
totalkineticenergy(obj::Physical) = kineticenergy(obj) + angularkineticenergy(obj)

##################################################
# Update
##################################################

function before_update!(obj::Physical, ::AbstractLayer, t::Float64, dt::Float64)
    savestate!(obj, t)
    if !isstatic(obj)
        boost!(obj, impulse(obj)/mass(obj)[1])
        angularboost!(obj, angularimpulse(obj)/angularmass(obj)[1])
        timeevolve!(obj, dt)
    end
    resetimpulse!(obj)
    return obj
end

##################################################
# State functions
##################################################

"""
    savestate!(obj, t)

Save the current state of `obj` and push a copy of it with a new
timestamp, `t`, onto the state history buffer.

"""
function savestate!(obj::Physical, t::Float64)
    push!(obj.state_history, copy(currentstate(obj)))
    currentstate(obj).t = t
    return obj
end

##################################################
# Time translation
##################################################

function timeevolve!(obj::Physical, dt::Float64)
    translate!(obj, dt*velocity(obj))
    rotate!(obj, dt*angularvelocity(obj))
end

function timeinterpolate!(obj::Physical, dt::Float64)
    state = currentstate(obj)
end

##################################################
# Time projection
##################################################

function timecapture(shape1::AbstractShape, shape2::AbstractShape, (state1, prevstate1)::Tuple{PhysicalState, PhysicalState}, (state2, prevstate2)::Tuple{PhysicalState, PhysicalState})
    shape = shape1 |> rotate(prevstate1.θ) |> translate(prevstate1.x, prevstate1.y)
    lines = unique!(collect(vertices(shape2))) |> rotate(prevstate2.θ) |>
        #!extrude(rotate(state2.θ - prevstate2.θ) |> translate(state2.x - prevstate2.x, state2.y - prevstate2.y)) |>
        translate(prevstate2.x, prevstate2.y)
    return shape, lines
end

function timetrace(shape1::AbstractShape, shape2::AbstractShape, (state1, prevstate1)::Tuple{PhysicalState, PhysicalState}, (state2, prevstate2)::Tuple{PhysicalState, PhysicalState})
    dt = state1.t - prevstate1.t
    ts = trace(timecapture(shape1, shape2, (state1, prevstate1), (state2, prevstate2))...)
    return (dt * (t - 1.) for t in ts)
end

function timecapture(obj1::Physical, obj2::Physical)
    state1, state2 = currentstate(obj1), currentstate(obj2)
    prevstate1, prevstate2 = obj1.state_history[end-1], obj2.state_history[end-1]
    return timecapture(obj1.shape, obj2.shape, (state1, prevstate1), (state2, prevstate2))
end

function timetrace(obj1::Physical, obj2::Physical)
    state1, state2 = currentstate(obj1), currentstate(obj2)
    prevstate1, prevstate2 = obj1.state_history[end-1], obj2.state_history[end-1]
    return timetrace(obj1.shape, obj2.shape, (state1, prevstate1), (state2, prevstate2))
end

timeintersects(obj1::Physical, obj2::Physical) =
    length(timetrace(obj1, obj2)) > 0 || length(timetrace(obj2, obj1)) > 0

##################################################
# Springs
##################################################

# abstract type AbstractHinge{T} <: AbstractPhysicalObject end
# abstract type AbstractSpring{T} <: AbstractPhysicalObject end

# struct Hinge{T} <: AbstractHinge{T}
#     obj::Physical
#     p::Vertex{T}
# end

# struct Spring{T} <: AbstractSpring{T}
#     hinge1::Hinge{T}
#     hinge2::Hinge{T}
#     k::T
#     l0::T
# end

# position(h::AbstractHinge) = h.p |> rotate(h.obj.θ) |> translate(h.obj.x, h.obj.y)
# previousposition(h::AbstractHinge) = h.p |> rotate(previousstate(h.obj).θ) |> translate(previousstate(h.obj).x, previousstate(h.obj).y)
# angle(v::Vertex) = 180*atan(v.y, v.x)/π

intersects(obj1::Physical, obj2::Physical) = intersects(shape(obj1), shape(obj2))
inside(obj::Physical, p::Vertex, θ) = inside(shape(obj), p, θ)

##################################################
# Collisions
##################################################

struct CollisionSolver{D} <: AbstractPhysicalObject
    contact_cache::Dict{Physical,Int}
    displacement_cache::Dict{Physical,SVector{D}}
end
CollisionSolver{D}() where {D} = CollisionSolver{D}(Dict{Physical,Int}(), Dict{Physical,SVector{D}}())

struct Collision{D, T1<:Physical, T2<:Physical, T} <: AbstractPhysicalObject
    t::Float64
    obj1::T1
    obj2::T2
    p::SVector{D,T}
    n::SVector{D,T}
    d::T
end
Base.isless(col1::Collision, col2::Collision) = col1.d < col2.d

function render!(layer::AbstractLayer, self::Collision, frame::Int, fps::Float64)
    render!(layer, Shapes.Polyline(Shapes.Vertex(self.p...), Shapes.Vertex((self.p |> translate((2*self.d)*self.n...))...)), 0., 0., 0., color=colorant"#F00")
    render!(layer, Shapes.Polyline(Shapes.Vertex(self.p...), Shapes.Vertex((self.p |> translate((self.d)*self.n...))...)), 0., 0., 0., color=colorant"#00F")
end

"""
    collisions(layer1, layer2, t, dt)

Return a sorted array of collisions that happened during the last
update step between objects in `layer1` and `layer2`

"""
function collisions(layer::AbstractLayer, t::Float64, dt::Float64)
    # TODO: broad phase / narrow phase separation
    cs = Collision[]
    for (i, obj1) in enumerate(layer.objects), (j, obj2) in enumerate(layer.objects)
        (i <= j || !(obj1 isa Physical) || !(obj2 isa Physical) || (isstatic(obj1) && isstatic(obj2))) && continue
        shape1 = shape(obj1)
        shape2 = shape(obj2)
        !intersects(BoundingBox(shape1), BoundingBox(shape2)) && continue
        d1, n1 = Shapes.separatingaxis(shape1, shape2)
        d2, n2 = Shapes.separatingaxis(shape2, shape1)
        # Shapes are not intersecting
        (d1 < 0 || d2 < 0) && continue
        # Find the shape
        if !(shape2 isa Shapes.HalfSpace) && (d1 <= d2 || shape1 isa Shapes.HalfSpace)
            d, n, p = d1, n1, Shapes.vector(Shapes.supportvertex(shape2, n1))
            push!(cs, Collision(t, obj1, obj2, p, n, d))
        else
            d, n, p = d2, n2, Shapes.vector(Shapes.supportvertex(shape1, n2))
            push!(cs, Collision(t, obj2, obj1, p, n, d))
        end
    end
    return cs
end

function before_update!(self::CollisionSolver, layer::AbstractLayer, t::Float64, dt::Float64)
    for obj in layer
        if obj isa Physical
            self.contact_cache[obj] = get(self.contact_cache, obj, 0)
            self.displacement_cache[obj] = get(self.displacement_cache, obj, zero(position(obj)))
        end
    end
end

function after_update!(self::CollisionSolver, layer::AbstractLayer, t::Float64, dt::Float64)
    cs = collisions(layer, t, dt)

    # Solve positions
    for c in cs
        self.contact_cache[c.obj1] += 1
        self.contact_cache[c.obj2] += 1
    end

    for _ in 1:10
        for c in cs
            solveposition!(c, self)
        end
    end

    for obj in keys(self.displacement_cache)
        self.contact_cache[obj] = 0
        translate!(obj, self.displacement_cache[obj])
        self.displacement_cache[obj] = zero(self.displacement_cache[obj])
    end

    # Solve collisions iteratively
    for _ in 1:1
        for c in cs
            solvevelocity!(c, self)
        end
    end
end

function solveposition!(self::Collision, solver::CollisionSolver)
    Δ = self.d*self.n + solver.displacement_cache[self.obj1] - solver.displacement_cache[self.obj2]
    j = (self.n ⋅ Δ)
    f = isstatic(self.obj1) || isstatic(self.obj2) ? 1. : 0.5
    if !isstatic(self.obj1)
        solver.displacement_cache[self.obj1] -= self.n*j*f/solver.contact_cache[self.obj1]
    end
    if !isstatic(self.obj2)
        solver.displacement_cache[self.obj2] += self.n*j*f/solver.contact_cache[self.obj2]
    end
end

function solvevelocity!(self::Collision, solver::CollisionSolver)
    vr, j, Δ = collide(self.obj1, self.obj2, self.n, self.p, self.d)
    vn = vr ⋅ self.n
    if vn <= 0
        if !isstatic(self.obj1)
            # x = Δ/mass(self.obj1)[1]
            # translate!(self.obj1, -x)
            #x = solver.displacement_cache[self.obj1]
            applyimpulse!(self.obj1, -j, self.p)
            # if norm(impulse(self.obj1)) < 30
            #     setimpulse!(self.obj1, zero(impulse(self.obj1)))
            #     setvelocity!(self.obj1, zero(velocity(self.obj1)))
            # end
        end
        if !isstatic(self.obj2)
            # x = Δ/mass(self.obj2)[1]
            # translate!(self.obj2, x)
            #x = solver.displacement_cache[self.obj2]
            applyimpulse!(self.obj2, j, self.p)
            # if norm(impulse(self.obj2)) < 30
            #     setimpulse!(self.obj2, zero(impulse(self.obj2)))
            #     setvelocity!(self.obj2, zero(velocity(self.obj2)))
            # end
        end
        # setvelocity!(self.obj1, SVector(0., 0.))
        # setvelocity!(self.obj2, SVector(0., 0.))
        # setangularvelocity!(self.obj1, SVector(0.))
        # setangularvelocity!(self.obj2, SVector(0.))
    end
end

"""
    collide(obj1, obj2, n, p, d, CR)

Return a tuple, `vr, j, Δ`, where `vr` is the relative velocity at the
contact point, `p`, `j` is the resulting impulse to be applied at `p`,
and `Δ` is the mass–weighted displacement for resolving the collision.

"""
function collide(obj1::Physical, obj2::Physical, normal::SVector, p::SVector, d::Real)
    m1, m2 = isstatic(obj1) ? Inf : mass(obj1)[1], isstatic(obj2) ? Inf : mass(obj2)[1]
    I1, I2 = isstatic(obj1) ? Inf : angularmass(obj1)[1], isstatic(obj2) ? Inf : angularmass(obj2)[1]
    p1, p2 = position(obj1), position(obj2)
    θ1, θ2 = angle(obj1), angle(obj2)
    v1 = isinf(m1) ? velocity(obj1) : velocity(obj1) + impulse(obj1)/m1
    v2 = isinf(m2) ? velocity(obj2) : velocity(obj2) + impulse(obj2)/m2
    ω1 = isinf(I1) ? angularvelocity(obj1)[1] : angularvelocity(obj1)[1] + angularimpulse(obj1)[1]/I1
    ω2 = isinf(I2) ? angularvelocity(obj2)[1] : angularvelocity(obj2)[1] + angularimpulse(obj2)[1]/I2

    # Setup
    r1 = p - p1
    r2 = p - p2
    vr = (v2 + ω2*(r2 |> rotate(90))*pi/180) - (v1 + ω1*(r1 |> rotate(90))*pi/180)
    vn = vr ⋅ normal

    vperp = vr - vn*normal
    N = norm(vperp)
    tangent = iszero(N) ? vperp : vperp / N
    vt = vr ⋅ tangent

    A1 = r1 × normal
    A2 = r2 × normal
    m = inv(inv(m1) + inv(m2) + inv(I1)*A1^2 + inv(I2)*A2^2)

    # Collision impulse
    R = restitution(obj1, obj2)
    jc = -(1 + R)*vn*m
    if vn^2 < 1000
        jc = max(0., -vn*m)
    end

    # Friction impulse
    F = norm(vr) > 0 ? dynamicfriction(obj1, obj2) : staticfriction(obj1, obj2)
    jf = min(F*jc, vt*m)

    # Total impulse
    j = jc*normal - jf*tangent

    # Displacement share
    Δ = inv(inv(m1) + inv(m2))*normal*d

    return vr, j, Δ
end

restitution(::Physical, ::Physical) = 0.
dynamicfriction(::Physical, ::Physical) = 0.
staticfriction(::Physical, ::Physical) = 0.

end #module


text:
"""
Polynomial optimization example from [Bertsimas, Nohadani, 2010]
"""
function poly1()
    m = JuMP.Model()
    @variable(m, -1 <= x <= 4)
    @variable(m, -1 <= y <= 4)
    @variable(m, obj)
    @objective(m, Min, obj)
    gm = GlobalModel(model = m)
    set_param(gm, :ignore_accuracy, true)
    set_param(gm, :ignore_feasibility, true)
    set_param(gm, :abstol, 1e-3)

    add_nonlinear_constraint(gm, :((x,y) -> 2*x^6 - 12.2*x^5 + 21.2*x^4 + 6.2*x - 6.4*x^3 - 4.7*x^2 +
    y^6 - 11*y^5 + 43.3*y^4 - 10*y - 74.8*y^3 + 56.9*y^2 - 4.1*x*y - 0.1*y^2*x^2 + 0.4*y^2*x + 0.4*x^2*y),
    name = "objective", dependent_var = obj)
    add_nonlinear_constraint(gm, :((x,y) -> 10.125 - (x-1.5)^4 - (y-1.5)^4), name = "h1")
    add_nonlinear_constraint(gm, :((x,y) -> (2.5 - x)^3 + (y+1.5)^3 - 15.75), name = "h2")
    return gm
end