text:
###################################################################################
# Antony Schutz 2015, ANR - POLCA - 2016
###################################################################################
## Data Type
###################################################################################
# PAINTER Data Type
#
# Structure containing all user data and oifits data
# OIDATA::PAINTER_Input
#
###################################################################################
# Initialise structure with nothing
###################################################################################
type PAINTER_Input
    PlotFct::Function               # user defined plot function
    Folder::String             # Folder of OIFITS/FITS Files
    FilesName::Array{String}   # name of the files present in Folder
    indfile::Array{Int,1}           # index of file used
    indwvl::Array{Int,1}            # index of wavelength used
    wvl::Array{Float64}             # mean value of anayzed wavelength
    U::Array{Float64}               # U spatial frequencies
    V::Array{Float64}               # V spatial frequencies
    P::Array{Float64}               # Matrix of V2
    W::Array{Float64}               # Matrix of V2err
    Xi::Dict                        # Dictionary of Vector of phases difference
    K::Dict                         # Dictionary of Vector of phases difference error
    Closure_index::Array{Int}       # index of phases closure
    nb::Int                         # number of bases (spatiale frequencies per wavelength)
    nw::Int                         # number of wavelength
    nx::Int                         # side size of image in pixels
    FOV::Real                       # Field Of View in arcsecond
    lambda_spat::Real               # From here to end: admm variables
    lambda_spec::Real
    lambda_L1::Real
    rho_y::Real

    rho_y_gamma::Real
    rho_y_xi::Real

    rho_spat::Real
    rho_spec::Real
    rho_ps::Real
    alpha::Real
    beta::Real
    eps1::Real
    eps2::Real
    epsilon::Real
    mask3D::Array                   # support constraint array
    xinit3D::Array                  # initial estimate array
    Wvlt::Array                     # list of wavelets basis
    T3::Array{Float64}              # Matrix of phases closure
    T3err::Array{Float64}           # Matrix of phases closure error
    DP::Array{Float64}              # Matrix of Differential phases
    DPerr::Array{Float64}           # Matrix of Differential phases error
    isDP::Int                       # To know if there is DP in data
    dptype::String             # way to construct diff phases
    dpprm::Int                      # horizon for sliding dp or index of reference
    baseNb::Dict
    orderedCluster::Dict
    flux::Array
    norm::Real
    # Kbefore::Dict
end
# Structure containing all data which are modified during admm
# PDATA::PAINTER_Data
type PAINTER_Data
    eta::Real
    plan::Array
    F3D::Array
    H::Dict
    M::Array
    x::Array
    z::Array
    tau_s::Array
    w::Array
    v::Array
    r::Array
    tau_w::Array
    tau_v::Array
    tau_r::Array
    Fx::Array
    tau_xc::Array
    tau_pwc::Array
    tau_xic::Array
    yc::Array
    y_v2::Array
    y_phi::Array
    crit1::Array
    crit2::Array
    ind::Int
    CountPlot::Int
    count::Int

end

function painterinputinit()
    return PAINTER_Input( painterplotfct, "", [], [], [], [], [], [], []
                         , [], Dict{}(), Dict{}(), [], 0., 0., 0., 0., 0., 0., 0., 0.
                         , 0., 0.
                         , 0., 0., 0., 0., 0., 0., 0., 0., [], [], []
                         , [], [], [], [], 0, "", 0, Dict{}(), Dict{}(),[],0) # [],Dict{}())
end
function painterdatainit()

  return PAINTER_Data(0., [], [], Dict{}(), [], [], [], [], [], []
                      , [], [], [], [], [], [], [], [], [], []
                      , [], Float64[], Float64[], 0, 0, 1)
end


text:
# ------ Root Finding Methods -----

"""
    brent(f, a, b; args=(), xtol=2e-12, rtol=4*eps(), maxiter=100, full_output=false)

1D root finding using Brent's method.  Based off the brentq implementation in scipy.

**Arguments**
- `f`: scalar function, that optionally takes additional arguments
- `a`::Float, `b::Float`: bracketing interval for a root - sign changes sign between: (f(a) * f(b) < 0)
- `args::Tuple`: tuple of additional arguments to pass to f
- `xtol::Float`: absolute tolerance (positive) for root
- `rtol::Float`: relative tolerance for root
- `maxiter::Int`: maximum number of iterations allowed
- `full_output::Bool`: flag to indicate whether you want just the root, or the root with a
    second argument (tuple) containing the number of iterations, function calls, and a convergence message.
"""
function brent(f, a, b; args=(), xtol=2e-12, rtol=4*eps(), maxiter=100, full_output=false)

    xpre = a; xcur = b
    xblk = 0.0; fblk = 0.0; spre = 0.0; scur = 0.0
    error_num = "INPROGRESS"

    fpre = f(xpre, args...)
    fcur = f(xcur, args...)
    funcalls = 2
    iterations = 0

    if fpre*fcur > 0
        error_num = "SIGNERR"
        return _pack_brent_results(0.0, iterations, funcalls, error_num, full_output)
    end
    if fpre == 0
        error_num = "CONVERGED"
        return _pack_brent_results(xpre, iterations, funcalls, error_num, full_output)
    end
    if fcur == 0
        error_num = "CONVERGED"
        return _pack_brent_results(xcur, iterations, funcalls, error_num, full_output)
    end

    for i = 1:maxiter
        iterations += 1
        if fpre*fcur < 0
            xblk = xpre
            fblk = fpre
            spre = scur = xcur - xpre
        end
        if abs(fblk) < abs(fcur)
            xpre = xcur
            xcur = xblk
            xblk = xpre

            fpre = fcur
            fcur = fblk
            fblk = fpre
        end

        delta = (xtol + rtol*abs(xcur))/2
        sbis = (xblk - xcur)/2
        if fcur == 0 || abs(sbis) < delta
            error_num = "CONVERGED"
            return _pack_brent_results(xcur, iterations, funcalls, error_num, full_output)
        end

        if abs(spre) > delta && abs(fcur) < abs(fpre)
            if xpre == xblk
                # interpolate
                stry = -fcur*(xcur - xpre)/(fcur - fpre)
            else
                # extrapolate
                dpre = (fpre - fcur)/(xpre - xcur)
                dblk = (fblk - fcur)/(xblk - xcur)
                stry = -fcur*(fblk*dblk - fpre*dpre)/(dblk*dpre*(fblk - fpre))
            end
            if 2*abs(stry) < min(abs(spre), 3*abs(sbis) - delta)
                # good short step
                spre = scur
                scur = stry
            else
                # bisect
                spre = sbis
                scur = sbis
            end
        else
            # bisect
            spre = sbis
            scur = sbis
        end

        xpre = xcur; fpre = fcur
        if abs(scur) > delta
            xcur += scur
        else
            xcur += (sbis > 0 ? delta : -delta)
        end

        fcur = f(xcur, args...)
        funcalls += 1
    end
    error_num = "CONVERR"
    return _pack_brent_results(xcur, iterations, funcalls, error_num, full_output)
end

# TODO AN: replace w/ newer Brent method and automatic bracketing?

"""
private method

pack up results from brent method
"""
function _pack_brent_results(x, iter, fcalls, flag, full_output)
    if full_output
        return x, (iter=iter, fcalls=fcalls, flag=flag)
    else
        return x
    end
end

ext: jl


text:
using JAnki
using HTTP

cards = JAnki.anki_find_notes("transfer")["result"]

for card in cards
           card_info = JAnki.anki_get_card(card)["result"][1]
           println(card_info["noteId"])
           for key in keys(card_info["fields"])
               println("$key:")
               println(card_info["fields"][key]["value"])
               println("\n")
           end
           println("\n\n")
       end


text:
ENV["QSG_RENDER_LOOP"] = "basic"
using CxxWrap # for safe_cfunction
using QML
using Observables
using Plots
import GR

ENV["GKSwstype"] = "use_default"
gr(show=true)

const qmlfile = joinpath(dirname(Base.source_path()), "qml_ex.qml")

nbins = Observable(30)
w, h = (600, 450)

# Arguments here need to be the "reference types", hence the "Ref" suffix
function paint(p::CxxPtr{QPainter}, item::CxxPtr{JuliaPaintedItem})
  global w, h

  ENV["GKS_WSTYPE"] = 381
  ENV["GKS_CONID"] = split(repr(p.cpp_object), "@")[2]

  dev = device(p[])[]
  r = effectiveDevicePixelRatio(window(item[])[])
  w, h = width(dev) / r, height(dev) / r

  num_bins = Int64(round(nbins[]))
  histogram(randn(10000), nbins=num_bins, size=(w, h))

  return
end

function mousePosition(eventx, eventy, deltay)
  if w > h
    xn = eventx / w
    yn = (h - eventy) / w
  else
    xn = eventx / h
    yn = (h - eventy) / h
  end
  x, y = GR.ndctowc(xn, yn)
  "($(round(x,digits=4)), $(round(y,digits=4)))"
end

load(qmlfile,
  paint_cfunction = @safe_cfunction(paint, Cvoid, (CxxPtr{QPainter}, CxxPtr{JuliaPaintedItem})),
  nbins = nbins
)

@qmlfunction mousePosition

exec()


text:
using Symbolics, SymbolicUtils
using SymbolicUtils: istree, operation, arguments
using Symbolics: value, get_variables, solve_for, derivative
using SymbolicUtils.Rewriters
using SymbolicUtils.Code

"""
    sym_taylor generates a Taylor series on variable x to order n and at
    point x₀, i.e., Σ aᵢ*(x - x₀)ⁱ
"""
function sym_taylor(x, n, x₀=0)
    x = value(x)
    @variables a[0:n]
    a = value.(a)
    a[1] + sum([a[i+1]*(x-x₀)^i for i = 1:n])
end

"""
    sym_frobenius generates a Frobenius series on variable x to order n and at
    point x₀, i.e., (x - x₀)^α Σ aᵢ*(x - x₀)ⁱ
"""
function sym_frobenius(x, n, x₀=0)
    x = value(x)
    @variables a[0:n]
    a = value.(a)
    @syms α
    (x-x₀)^α * sym_taylor(x, n, x₀)
end

# pox (power-of-x) is a symbolic function to keep track of the powers of x
# pox(k,n) means k*x^n
@syms pox(k, n)

is_pox(x) = istree(x) && operation(x)==pox
is_not_pox(x) = !is_pox(x)

get_coef(p) = is_pox(p) ? arguments(p)[1] : p
get_power(p) = is_pox(p) ? arguments(p)[2] : 0

replace_x(eq, x) = substitute(eq, Dict(x => pox(1,1)))

count_rule1 = @rule ^(pox(~k, ~n1), ~n2) => isequal(~k,1) ? pox(1, ~n1 * ~n2) : pox(^(~k,~n2), ~n1 * ~n2)
count_rule2 = @rule pox(~k1, ~n1) * pox(~k2, ~n2) => pox(~k1 * ~k2, ~n1 + ~n2)
count_rule3 = @acrule pox(~k, ~n) * ~u::is_not_pox => pox(~k * ~u, ~n)

"""
    collect_powers separates the powers of x in eq (a polynomial) and returns
    a dictionary of power => term
"""
function collect_powers(eq, x)
    eq = expand(expand_derivatives(eq))
    eq = replace_x(eq, x)
    #eq = Prewalk(PassThrough(count_rule1))(eq)
    eq = Fixpoint(Prewalk(PassThrough(Chain([count_rule1, count_rule2, count_rule3]))))(eq)

    if !istree(eq)
        return Dict{Any, Any}(0 => eq)
    elseif is_pox(eq)
        return Dict{Any, Any}(get_power(eq) => get_coef(eq))
    else
        eqs = Dict{Any, Any}()

        for term in arguments(eq)
            n = get_power(term)
            if haskey(eqs, n)
                eqs[n] = eqs[n] + get_coef(term)
            else
                eqs[n] = get_coef(term)
            end
        end

        return eqs
    end
end

"""
    strip_powers removes the powers from the output of collect_powers
    and returns a sorted list of terms/sub-equations
"""
function strip_powers(eqs)
    @syms α
    d = Dict(α => 0)
    eqs = Dict(substitute(power, d) => eq for (power, eq) in eqs)
    collect(values(sort(eqs)))
end

"""
    solve_rational gets a triangular system of linear equations and solve
    for the variables one by one starting from the first equation

    eqs is a list of terms as prepared by strip_powers

    solve_rational tries to generate rational results but falls back to reals
    if not possible
"""
function solve_rational(eqs)
    eqs = Num.(eqs)
    vals = Dict()

    for eq in eqs
        vars = get_variables(eq)
        vars = filter(x -> !haskey(vals,x), vars)

        if length(vars) > 0
            dep = vars[end]
            eq = substitute(eq, vals)
            # the last variable of each eq is considered the dependent variable
            A = value(substitute(eq, Dict(dep => 0)))
            B = value(substitute(eq, Dict(dep => 1)))
            try
                vals[dep] = A // (B - A)
            catch
                vals[dep] = A / (B - A)
            end
        end
    end

    vals
end

"""
    finalize_taylor substitutes the calculated coefficients back into the
    original Taylor or Frobenius series
"""
function finalize_taylor(x, y, vals)
    u = substitute(y, vals)
    vars = [v for v in get_variables(y) if !isequal(v,x)]
    factor(u, vars)
end

"""
    A main entry point
    solve_taylor solves a Taylor series expansion of a linear differential equation

    x is the independent variable (a Symbolic Sym)
    y is the dependent variable (a Symbolic Sym)
    diffeq is the differential equation, assuming diffeq ~ 0
"""
function solve_taylor(x, y, diffeq)
    eqs = strip_powers(collect_powers(diffeq, x))
    vals = solve_rational(eqs)
    finalize_taylor(x, y, vals)
end

"""
    A main entry point
    solve_frobenius solves a Frobenius series expansion of a linear differential equation

    x is the independent variable (a Symbolic Sym)
    y is the dependent variable (a Symbolic Sym)
    diffeq is the differential equation, assuming diffeq ~ 0
"""
function solve_frobenius(x, y, diffeq; abstol=1e-8, rationalize_coef=false)
    eqs = strip_powers(collect_powers(diffeq, x))
    eq = eqs[1]
    @syms α
    d = Dict(v => 1 for v in get_variables(eq) if !isequal(v, α))

    rs, zs = find_roots(substitute(eq, d), α)
    isempty(rs) && error("The indicial equation has no real root")
    α₁ = rs[end]
    @info "α is $α₁"

    d = Dict(α => α₁)
    eqs = [substitute(eq,d) for eq in eqs]
    vals = solve_rational(eqs[2:end])
    if rationalize_coef
        vals = rationalize(vals)
    end
    vals[α] = α₁
    finalize_taylor(x, y, vals)
end

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

"""
    factor factors a polynomial eq over a list of variables (vars)
    it assumes trivial separatibilty
"""
function factor(eq, vars)
    f = []
    for v in vars
        d = Dict(u => (isequal(v,u) ? 1 : 0) for u in vars)
        push!(f, v * substitute(eq, d))
    end

    +(f...)
end

"""
    solve_newton is a symbolic Newton-Ralphson solver
    f is a symbolic equation to be solved (f ~ 0)
    x is the variable to solve
    x₀ is the initial guess
"""
function solve_newton(f, x, x₀; abstol=1e-8, maxiter=50)
    xₙ = Complex(x₀)
    ∂f = value(derivative(f, x))

    for i = 1:maxiter
        d = Dict(x => xₙ)
        xₙ₊₁ = xₙ - substitute(f, d) / substitute(∂f, d)
        if abs(xₙ₊₁ - xₙ) < abstol
            return xₙ₊₁
        else
            xₙ = xₙ₊₁
        end
    end
    return xₙ
end

"""
    poly_degree returns the degree of a polynomial equation eq
"""
function poly_degree(poly, x)
    eqs = collect_powers(poly, x)
    round(Int, maximum(keys(eqs)))
end

"""
    find_roots returns all the real and complex roots of a polynomial (poly)
    for variable x

    the output is rs, sz, where rs is a list of real roots and zs is a list of
    complex roots

    if poly is not a polynomial, the number of desired roots (n) needs to be
    provided. For example, find_roots(sin(x), x, 10) returns 10 different, but
    not necessarily sequential, multiples of π
"""
function find_roots(poly, x, n=-1; abstol=1e-8)
    n = (n == -1 ? poly_degree(poly, x) : n)

    rs = Float64[]
    zs = Complex[]

    while n > 0
        z = solve_newton(poly, x, Complex(rand(),rand()))
        if abs(imag(z)) < abstol
            r = real(z)
            push!(rs, r)
            poly = poly / (x - r)
            n -= 1
        else
            append!(zs, [z, conj(z)])
            poly = poly / ((x-z)*(x-conj(z)))
            n -= 2
        end
    end
    sort(rs), zs
end

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

function continued_fraction(x, n; abstol=1e-8)
    try
        l = zeros(Int, n)
        for i = 1:n
            l[i] = floor(Int, x)
            Δx = x - l[i]
            Δx < abstol && break
            x = 1 / Δx
        end
        y = 0
        for i = n:-1:1
            if l[i] + y > 0
                y = 1 // (l[i] + y)
            end
        end
        return 1 // y
    catch
        return x
    end
end

function float_to_rational(x; n=10, abstol=1e-8)
    s = x < 0
    y = continued_fraction(abs(x), n) * (s ? -1 : +1)
    if abs(x - y) < abstol
        return y
    else
        return x
    end
end

# function rationalize(vals)
#     r = Dict()
#     for p in vals
#         eq = last(p)
#         vars = get_variables(eq)
#         if length(vars) == 1
#             var = vars[1]
#             eq = eq / var
#             if !(eq isa Int) || !(eq isa Rational)
#                 eq = float_to_rational(eq)
#             end
#             r[first(p)] = eq * var
#         else
#             r[first(p)] = eq
#         end
#     end
#     r
# end

################################### Examples! #################################

function test_harmonic()
    @syms x
    y = sym_taylor(x, 10)
    D = Differential(x)
    solve_taylor(x, y, D(D(y)) + y)
end

function test_airy()
    @syms x
    y = sym_taylor(x, 10)
    D = Differential(x)
    solve_taylor(x, y, D(D(y)) + x*y)
end

function test_bessel(ν = sqrt(2); rationalize_coef=false)
    @syms x
    y = sym_frobenius(x, 10)
    D = Differential(x)
    solve_frobenius(x, y, x^2 * D(D(y)) + x * D(y) - (x^2 + ν^2)*y; rationalize_coef=rationalize_coef)
end


text:
# Automatically generated test file
function white_noise()

g = [
0.297287984535462,
0.382395967790608,
-0.597634476728231,
-0.010445244637376,
-0.839026854388764,
0.311111338498334,
2.295087823837310,
-2.267086348800531,
0.529965576166746,
0.431421526422912,
0.583708287568779,
0.963271605038191,
0.458790955053717,
-0.522336757421508,
0.408395838324752]

ζ_vec = [13.2605, 5.72288, 2.46984, 1.06592, 0.46002, 0.198532, 0.0856812, 0.0369777, 0.0159586, 0.00688728]

I_sols = [
[1, 15],
[1, 3, 15],
[1, 7, 8, 15],
[1, 7, 8, 9, 15],
[1, 5, 7, 8, 9, 15],
[1, 5, 7, 8, 9, 12, 15],
[1, 5, 7, 8, 9, 12, 14, 15],
[1, 5, 6, 7, 8, 9, 12, 14, 15],
[1, 3, 4, 5, 7, 8, 9, 12, 14, 15],
[1, 2, 3, 4, 5, 7, 8, 9, 12, 14, 15],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 12, 14, 15],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 12, 14, 15],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 13, 14, 15],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]]

f_sols = [
13.604680464796299,
13.260506698294682,
11.380933016032856,
6.070779395384360,
1.658799820085756,
1.459921430168823,
0.760578522235333,
0.586871209768404,
0.424250952535800,
0.235164307482329,
0.119283257822485,
0.048351419645437,
0.008609105509041,
-0.000000000000007]

return g, ζ_vec, I_sols, f_sols

end


text:
module DirectMotionPlanning

using LinearAlgebra, ForwardDiff, FiniteDiff, StaticArrays, SparseArrays
using MathOptInterface, Ipopt
using JLD2

include("lqr.jl")
include("unscented.jl")
include("indices.jl")
include("utils.jl")

include("integration.jl")
include("model.jl")

include("problem.jl")

include("objective.jl")
include("objectives/quadratic.jl")
include("objectives/penalty.jl")

include("constraints.jl")
include("constraints/dynamics.jl")

include("moi.jl")
include("solvers/snopt.jl")

end # module

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
"""
    Setup

This module is the first step in getting started with Vulkan. It's purpose is
to create a window using GLFW which supports Vulkan.
"""
module Setup

export main

using GLFW

"""
    main()

The application's entry point.
"""
main() = (
    buildApp()
    |> initWindow
    |> initVulkan
    |> mainLoop
    |> cleanup
)

"""
The application's state. Anything the main app needs to run will be included
here.
"""
Base.@kwdef mutable struct App
    window::Union{Some{GLFW.Window}, Nothing}
end

"""
Construct a new application instance.
"""
function buildApp()
    App(;
        window = nothing
    )
end

"""
    initWindow(app::App) :: App

Create the application's main window. Configure to support Vulkan instead of
OpenGL.
"""
function initWindow(app::App) :: App
    # NO_API disables GLFW's attempt to load an OpenGL context
    GLFW.WindowHint(GLFW.CLIENT_API, GLFW.NO_API)

    # Resizing takes extra management, so don't support it for now
    GLFW.WindowHint(GLFW.RESIZABLE, false)

    app.window = GLFW.CreateWindow(800, 600, "Vulkan") |> Some

    app
end

"""
    initVulkan(app::App) :: App

Initialize Vulkan resources for this application.
"""
function initVulkan(app::App) :: App
    app
end

"""
    mainLoop(app::App) :: App

The application's main loop. This function blocks until the window is closed
by clicking the 'x' or similarly exited.
"""
function mainLoop(app::App) :: App
    window = something(app.window)

    while !GLFW.WindowShouldClose(window)
        GLFW.PollEvents()
    end

    app
end

"""
    cleanup(app::App) :: App

Cleanup any and all graphical resources.
"""
function cleanup(app::App) :: App
    app.window |> something |> GLFW.DestroyWindow
    app.window = nothing
    app
end

end # End Setup Module


text:
using Documenter, diffprog

makedocs(
    modules = [diffprog],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "iamalwaysuncomfortable",
    sitename = "diffprog.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/iamalwaysuncomfortable/diffprog.jl.git",
    push_preview = true
)


text:
module AlgebraOverFields

using Printf: @printf, @sprintf
using ...Interfaces: dimension
using ...Prerequisites: atol, rtol, Float, rawtype
using ...Prerequisites.NamedVectors: NamedVector
using ...Prerequisites.TypeTraits: efficientoperations
using ..Combinatorics: AbstractCombinatorics

import ...Interfaces: rank, add!, sub!, mul!, div!, ⊗, ⋅, sequence, permute

export SimpleID, ID
export idpropertynames, idisless
export Element, ScalarElement, Elements
export scalartype, idtype

"""
    SimpleID <: NamedVector

A simple id is the building block of the id system of an algebra over a field.
"""
abstract type SimpleID <: NamedVector end

"""
    ID{I<:SimpleID, N}

The (composite) id system of an algebra over a field.

Type alias for `NTuple{N, I} where {N, I<:SimpleID}`.
"""
const ID{I<:SimpleID, N} = NTuple{N, I}

"""
    ID(ids::SimpleID...)
    ID(ids::NTuple{N, SimpleID}) where N

Get the composite id from simple ids.
"""
ID(ids::SimpleID...) = ids
ID(ids::NTuple{N, SimpleID}) where {N} = ids

"""
    ID(::Type{SID}, attrs::Vararg{NTuple{N}, M}) where {SID<:SimpleID, N, M}

Get the composite id from the components of simple ids.
"""
@generated function ID(::Type{SID}, attrs::Vararg{NTuple{N, Any}, M}) where {SID<:SimpleID, N, M}
    exprs = []
    for i = 1:N
        args = [:(attrs[$j][$i]) for j = 1:M]
        push!(exprs, :(SID($(args...))))
    end
    return :(ID($(exprs...)))
end

"""
    idpropertynames(::Type{I}) where I<:ID{SimpleID} -> Tuple{Vararg{Symbol}}

Get the property names of a composite id.
"""
@generated function idpropertynames(::Type{I}) where I<:ID{SimpleID}
    exprs = [QuoteNode(Symbol(name, 's')) for name in I|>eltype|>fieldnames]
    return Expr(:tuple, exprs...)
end

"""
    getproperty(cid::ID{SimpleID}, name::Symbol)

Get the property of a composite id.
"""
Base.getproperty(cid::ID{SimpleID}, name::Symbol) = idgetproperty(cid, Val(name), Val(cid|>typeof|>idpropertynames))
@generated function idgetproperty(cid::ID{SimpleID}, ::Val{name}, ::Val{names}) where {name, names}
    index = findfirst(isequal(name), names)::Int
    exprs = [:(getfield(cid[$i], $index)) for i = 1:fieldcount(cid)]
    return Expr(:tuple, exprs...)
end

"""
    show(io::IO, cid::Tuple{SimpleID, Vararg{SimpleID}})

Show a composite id.
"""
Base.show(io::IO, cid::Tuple{SimpleID, Vararg{SimpleID}}) = @printf io "ID(%s)" join(cid, ", ")

"""
    promote_type(::Type{Tuple{}}, ::Type{Tuple{}})
    promote_type(::Type{Tuple{}}, I::Type{<:ID{SimpleID, N}}) where N
    promote_type(I::Type{<:ID{SimpleID, N}}, ::Type{Tuple{}}) where N

Define the promote rule for ID types.
"""
Base.promote_type(::Type{Tuple{}}, ::Type{Tuple{}}) = Tuple{}
Base.promote_type(::Type{Tuple{}}, I::Type{<:ID{SimpleID, N}}) where {N} = ID{I|>eltype}
Base.promote_type(I::Type{<:ID{SimpleID, N}}, ::Type{Tuple{}}) where {N} = ID{I|>eltype}

"""
    idisless(cid1::ID{SimpleID}, cid2::ID{SimpleID}) -> Bool

Compare two ids and judge whether the first is less than the second.

We define that ids with smaller ranks are always less than those with higher ranks. If two ids are of the same rank, the comparison goes just like that between tuples.
"""
function idisless(cid1::ID{SimpleID}, cid2::ID{SimpleID})
    r1, r2 = cid1|>rank, cid2|>rank
    (r1 < r2) ? true : (r1 > r2) ? false : isless(cid1, cid2)
end

"""
    rank(id::ID{SimpleID}) -> Int
    rank(::Type{<:ID{SimpleID}}) -> Any
    rank(::Type{<:ID{SimpleID, N}}) where N -> Int

Get the rank of a composite id.
"""
rank(id::ID{SimpleID}) = id |> typeof |> rank
rank(::Type{<:ID{SimpleID}}) = Any
rank(::Type{<:ID{SimpleID, N}}) where {N} = N

"""
    *(sid1::SimpleID, sid2::SimpleID) -> ID{SimpleID}
    *(sid::SimpleID, cid::ID{SimpleID}) -> ID{SimpleID}
    *(cid::ID{SimpleID}, sid::SimpleID) -> ID{SimpleID}
    *(cid1::ID{SimpleID}, cid2::ID{SimpleID}) -> ID{SimpleID}

Get the product of the id system.
"""
Base.:*(sid1::SimpleID, sid2::SimpleID) = ID(sid1, sid2)
Base.:*(sid::SimpleID, cid::ID{SimpleID}) = ID(sid, cid...)
Base.:*(cid::ID{SimpleID}, sid::SimpleID) = ID(cid..., sid)
Base.:*(cid1::ID{SimpleID}, cid2::ID{SimpleID}) = ID(cid1..., cid2...)

"""
    Element{V, I<:ID{SimpleID}}

An element of an algebra over a field.

The first and second attributes of an element must be
- `value`: the coefficient of the element
- `id::ID{SimpleID}`: the id of the element
"""
abstract type Element{V, I<:ID{SimpleID}} end

"""
    ScalarElement{V}

Identity element.
"""
const ScalarElement{V} = Element{V, Tuple{}}
scalartype(::Type{M}) where M<:Element = rawtype(M){M|>valtype, Tuple{}}

"""
    valtype(m::Element)
    valtype(::Type{<:Element})
    valtype(::Type{<:Element{V}}) where V

Get the type of the value of an element.

The result is also the type of the field over which the algebra is defined.
"""
Base.valtype(m::Element) = m |> typeof |> valtype
Base.valtype(::Type{<:Element}) = Any
Base.valtype(::Type{<:Element{V}}) where {V} = V

"""
    idtype(m::Element)
    idtype(::Type{<:Element{V, I} where V}) where I<:ID{SimpleID}
    idtype(::Type{<:Element})

The type of the id of an element.
"""
idtype(m::Element) = m |> typeof |> idtype
idtype(::Type{<:Element{V, I} where V}) where {I<:ID{SimpleID}} = I
idtype(::Type{<:Element}) = ID{SimpleID}

"""
    rank(m::Element) -> Int
    rank(::Type{M}) where M<:Element -> Int

Get the rank of an element.
"""
rank(m::Element) = m |> typeof |> rank
rank(::Type{M}) where {M<:Element} = M |> idtype |> rank

"""
    ==(m1::Element, m2::Element) -> Bool

Compare two elements and judge whether they are equal to each other.
"""
Base.:(==)(m1::Element, m2::Element) = ==(efficientoperations, m1, m2)

"""
    isequal(m1::Element, m2::Element) -> Bool

Compare two elements and judge whether they are equal to each other.
"""
Base.isequal(m1::Element, m2::Element) = isequal(efficientoperations, m1, m2)

"""
    isapprox(m1::Element, m2::Element; atol::Real=atol, rtol::Real=rtol) -> Bool

Compare two elements and judge whether they are inexactly equivalent to each other.
"""
Base.isapprox(m1::Element, m2::Element; atol::Real=atol, rtol::Real=rtol) = isapprox(efficientoperations, Val((:value,)), m1, m2; atol=atol, rtol=rtol)

"""
    replace(m::Element; kwargs...) -> typeof(m)

Return a copy of a concrete `Element` with some of the field values replaced by the keyword arguments.
"""
Base.replace(m::Element; kwargs...) = replace(efficientoperations, m; kwargs...)

"""
    promote_rule(::Type{M1}, ::Type{M2}) where {M1<:Element, M2<:Element}

Define the promote rule for Element types.
"""
function Base.promote_rule(::Type{M1}, ::Type{M2}) where {M1<:Element, M2<:Element}
    (M1 <: M2) && return M2
    (M2 <: M1) && return M1
    v1, i1, r1 = M1|>valtype, M1|>idtype, M1|>rank
    v2, i2, r2 = M2|>valtype, M2|>idtype, M2|>rank
    V, I = promote_type(v1, v2), promote_type(i1, i2)
    M = (r2 == 0) ? rawtype(M1) : (r1 == 0) ? rawtype(M2) : typejoin(rawtype(M1), rawtype(M2))
    return isconcretetype(I) ? M{V, I} : M{V, <:I}
end

"""
    mul_promote_type(::Type{M}, ::Type{V}) where {M<:Element, V<:Number}
    mul_promote_type(::Type{V}, ::Type{M}) where {M<:Element, V<:Number}

Define the promote rule for the multiplication between an Element and a scalar.
"""
mul_promote_type(::Type{M}, ::Type{V}) where {M<:Element, V<:Number} = mul_promote_type(V, M)
function mul_promote_type(::Type{V}, ::Type{M}) where {M<:Element, V<:Number}
    C, I = promote_type(M|>valtype, V), M|>idtype
    return isconcretetype(I) ? rawtype(M){C, I} : rawtype(M){C, <:I}
end

"""
    one(::Type{M}) where M<:Element

Get the identity operator.
"""
function Base.one(::Type{M}) where M<:Element
    rtype = rawtype(M)
    vtype = isconcretetype(valtype(M)) ? valtype(M) : Int
    @assert fieldnames(rtype) == (:value, :id) "one error: not supproted type($(nameof(rtype)))."
    return rtype(one(vtype), ID())
end

"""
    convert(::Type{M}, m::ScalarElement) where M<:ScalarElement
    convert(::Type{M}, m::Number) where M<:ScalarElement
    convert(::Type{M}, m::Element) where M<:Element

1) Convert a scalar element from one type to another;
2) Convert a scalar to a scalar element;
3) Convert an element from one type to another.
"""
function Base.convert(::Type{M}, m::ScalarElement) where M<:ScalarElement
    (typeof(m) <: M) && return m
    @assert fieldnames(M) == (:value, :id) "convert error: not supported type($(nameof(M)))."
    return rawtype(M)(convert(M|>valtype, m.value), m.id)
end
function Base.convert(::Type{M}, m::Number) where M<:ScalarElement
    @assert fieldnames(M) == (:value, :id) "convert error: not supported type($(nameof(M)))."
    return rawtype(M)(convert(M|>valtype, m), ID())
end
function Base.convert(::Type{M}, m::Element) where M<:Element
    (typeof(m) <: M) && return m
    @assert idtype(m) <: idtype(M) "convert error: dismatched ID type."
    @assert rawtype(typeof(m)) <: rawtype(M) "convert error: dismatched raw Element type."
    return replace(m, value=convert(valtype(M), m.value))
end

"""
    Elements{I<:ID{SimpleID}, M<:Element} <: AbstractDict{I, M}

An set of elements of an algebra over a field.

Type alias for `Dict{I<:ID{SimpleID}, M<:Element}`.
Similar iterms are automatically merged thanks to the id system.
"""
const Elements{I<:ID{SimpleID}, M<:Element} = Dict{I, M}

"""
    Elements(ms)
    Elements(ms::Pair{I, M}...) where {I<:ID{SimpleID}, M<:Element}
    Elements(ms::Element...)

Get the set of elements with similar items merged.
"""
Elements(ms) = Base.dict_with_eltype((K, V)->Dict{K, V}, ms, eltype(ms))
function Elements(ms::Pair{I, M}...) where {I<:ID{SimpleID}, M<:Element}
    result = Elements{I, M}()
    for (id, m) in ms
        result[id] = m
    end
    return result
end
function Elements(ms::Element...)
    result = Elements{ms|>eltype|>idtype, ms|>eltype}()
    for m in ms add!(result, m) end
    return result
end

"""
    show(io::IO, ms::Elements)

Show a set of elements.
"""
function Base.show(io::IO, ms::Elements)
    @printf io "Elements with %s entries:\n" length(ms)
    for m in values(ms)
        @printf io "  %s\n" m
    end
end

"""
    repr(ms::Elements) -> String

Get the repr representation of a set of elements.
"""
function Base.repr(ms::Elements)
    cache = [@sprintf("Elements with %s entries:", length(ms))]
    for m in ms|>values
        push!(cache, @sprintf("  %s", repr(m)))
    end
    return join(cache, "\n")
end

"""
    zero(ms::Elements) -> Nothing
    zero(::Type{<:Elements}) -> Nothing

Get a zero set of elements.

A zero set of elements is defined to be the one with no elements.
"""
Base.zero(ms::Elements) = ms |> typeof |> zero
Base.zero(::Type{<:Elements}) = nothing

"""
    ==(ms::Elements, ::Nothing) -> Bool
    ==(::Nothing, ms::Elements) -> Bool

Judge whether a set of elements is identically empty.
"""
Base.:(==)(ms::Elements, ::Nothing) = length(ms) == 0
Base.:(==)(::Nothing, ms::Elements) = length(ms) == 0

"""
    isequal(ms::Elements, ::Nothing) -> Bool
    isequal(::Nothing, ms::Elements) -> Bool

Judge whether a set of elements is identically empty.
"""
Base.isequal(ms::Elements, ::Nothing) = length(ms) == 0
Base.isequal(::Nothing, ms::Elements) = length(ms) == 0

"""
    add!(ms::Elements) -> typeof(ms)
    add!(ms::Elements, ::Nothing) -> typeof(ms)
    add!(ms::Elements, m::Number) -> typeof(ms)
    add!(ms::Elements, m::Element) -> typeof(ms)
    add!(ms::Elements, mms::Elements) -> typeof(ms)

Get the inplace addition of elements to a set.
"""
add!(ms::Elements) = ms
add!(ms::Elements, ::Nothing) = ms
add!(ms::Elements, m::Number) = add!(ms, convert(scalartype(ms|>valtype), m))
function add!(ms::Elements, m::Element)
    m = convert(ms|>valtype, m)
    old = get(ms, m.id, nothing)
    new = (old === nothing) ? m : replace(old, value=old.value+m.value)
    (abs(new.value) == 0.0) ? delete!(ms, m.id) : (ms[m.id] = new)
    return ms
end
add!(ms::Elements, mms::Elements) = (for m in mms|>values add!(ms, m) end; ms)

"""
    sub!(ms::Elements) -> typeof(ms)
    sub!(ms::Elements, ::Nothing) -> typeof(ms)
    add!(ms::Elements, m::Number) -> typeof(ms)
    sub!(ms::Elements, m::Element) -> typeof(ms)
    sub!(ms::Elements, mms::Elements) -> typeof(ms)

Get the inplace subtraction of elements from a set.
"""
sub!(ms::Elements) = ms
sub!(ms::Elements, ::Nothing) = ms
sub!(ms::Elements, m::Number) = add!(ms, convert(scalartype(ms|>valtype), -m))
function sub!(ms::Elements, m::Element)
    m = convert(ms|>valtype, m)
    old = get(ms, m.id, nothing)
    new = (old === nothing) ? -m : replace(old, value=old.value-m.value)
    (abs(new.value) == 0.0) ? delete!(ms, m.id) : (ms[m.id] = new)
    return ms
end
sub!(ms::Elements, mms::Elements) = (for m in mms|>values sub!(ms, m) end; ms)

"""
    mul!(ms::Elements, factor::ScalarElement) -> Elements
    mul!(ms::Elements, factor::Number) -> Elements

Get the inplace multiplication of elements with a scalar.
"""
mul!(ms::Elements, factor::ScalarElement) = mul!(ms, factor.value)
function mul!(ms::Elements, factor::Number)
    @assert isa(one(ms|>valtype|>valtype)*factor, ms|>valtype|>valtype) "mul! error: dismatched type, $(ms|>valtype) and $(factor|>typeof)."
    for m in values(ms)
        ms[m.id] = replace(m, value=m.value*factor)
    end
    return ms
end

"""
    div!(ms::Elements, factor::ScalarElement) -> Elements
    div!(ms::Elements, factor::Number) -> Elements

Get the inplace division of element with a scalar.
"""
div!(ms::Elements, factor::ScalarElement) = div!(ms, factor.value)
function div!(ms::Elements, factor::Number)
    @assert isa(one(ms|>valtype|>valtype)/factor, ms|>valtype|>valtype) "div! error: dismatched type, $(ms|>valtype) and $(factor|>typeof)."
    for m in values(ms)
        ms[m.id] = replace(m, value=m.value/factor)
    end
    return ms
end

"""
    +(m::Element) -> typeof(m)
    +(m::Element, ::Nothing) -> typeof(m)
    +(::Nothing, m::Element) -> typeof(m)
    +(m::Element, factor::Number) -> Elements
    +(factor::Number, m::Element) -> Elements
    +(m1::Element, m2::Element) -> Elements
    +(ms::Elements) -> typeof(ms)
    +(ms::Elements, ::Nothing) -> typeof(ms)
    +(::Nothing, ms::Elements) -> typeof(ms)
    +(ms::Elements, factor::Number) -> Elements
    +(factor::Number, ms::Elements) -> Elements
    +(ms::Elements, m::Element) -> Elements
    +(m::Element, ms::Elements) -> Elements
    +(ms1::Elements, ms2::Elements) -> Elements

Overloaded `+` operator between elements of an algebra over a field.
"""
Base.:+(m::Element) = m
Base.:+(ms::Elements) = ms
Base.:+(m::Element, ::Nothing) = m
Base.:+(::Nothing, m::Element) = m
Base.:+(ms::Elements, ::Nothing) = ms
Base.:+(::Nothing, ms::Elements) = ms
Base.:+(factor::Number, m::Element) = m + rawtype(m|>typeof)(factor, ID())
Base.:+(m::Element, factor::Number) = m + rawtype(m|>typeof)(factor, ID())
Base.:+(factor::Number, m::ScalarElement) = replace(m, value=m.value+factor)
Base.:+(m::ScalarElement, factor::Number) = replace(m, value=m.value+factor)
Base.:+(m1::ScalarElement, m2::ScalarElement) = replace(m1, value=m1.value+m2.value)
Base.:+(factor::Number, ms::Elements) = ms + rawtype(ms|>valtype)(factor, ID())
Base.:+(ms::Elements, factor::Number) = ms + rawtype(ms|>valtype)(factor, ID())
function Base.:+(m1::Element, m2::Element)
    I = promote_type(m1|>idtype, m2|>idtype)
    M = promote_type(typeof(m1), typeof(m2))
    return add!(Elements{I, M}(m1.id=>m1), m2)
end
Base.:+(m::Element, ms::Elements) = ms + m
function Base.:+(ms::Elements, m::Element)
    I = promote_type(ms|>keytype, m|>idtype)
    M = promote_type(valtype(ms), typeof(m))
    return add!(Elements{I, M}(ms), m)
end
function Base.:+(ms1::Elements, ms2::Elements)
    I = promote_type(ms1|>keytype, ms2|>keytype)
    M = promote_type(valtype(ms1), valtype(ms2))
    return add!(Elements{I, M}(ms1), ms2)
end

"""
    *(m::Element, ::Nothing) -> Nothing
    *(::Nothing, m::Element) -> Nothing
    *(factor::Number, m::Element) -> Element
    *(m::Element, factor::Number) -> Element
    *(m1::Element, m2::Element) -> Element
    *(ms::Elements, ::Nothing) -> Nothing
    *(::Nothing, ms::Elements) -> Nothing
    *(factor::Number, ms::Elements) -> Elements
    *(ms::Elements, factor::Number) -> Elements
    *(m::Element, ms::Elements) -> Elements
    *(ms::Elements, m::Element) -> Elements
    *(ms1::Elements, ms2::Elements) -> Elements

Overloaded `*` operator for element-scalar multiplications and element-element multiplications of an algebra over a field.
"""
Base.:*(m::Element, ::Nothing) = nothing
Base.:*(::Nothing, m::Element) = nothing
Base.:*(ms::Elements, ::Nothing) = nothing
Base.:*(::Nothing, ms::Elements) = nothing
Base.:*(factor::ScalarElement, m::Element) = m * factor.value
Base.:*(m::Element, factor::ScalarElement) = m * factor.value
Base.:*(m1::ScalarElement, m2::ScalarElement) = replace(m1, value=m1.value*m2.value)
Base.:*(factor::Number, m::Element) = m * factor
Base.:*(m::Element, factor::Number) = replace(m, value=factor*m.value)
Base.:*(factor::ScalarElement, ms::Elements) = ms * factor.value
Base.:*(ms::Elements, factor::ScalarElement) = ms * factor.value
Base.:*(factor::Number, ms::Elements) = ms * factor
function Base.:*(ms::Elements, factor::Number)
    (abs(factor) == 0) && return zero(Elements)
    result = Elements{ms|>keytype, mul_promote_type(ms|>valtype, factor|>typeof)}()
    for (id, m) in ms
        result[id] = m * factor
    end
    return result
end
Base.:*(m::Element, ms::Elements) = Elements((m * mm for mm in ms|>values)...)
Base.:*(ms::Elements, m::Element) = Elements((mm * m for mm in ms|>values)...)
Base.:*(ms1::Elements, ms2::Elements) = Elements((m1 * m2 for m1 in ms1|>values for m2 in ms2|>values)...)
function Base.:*(m1::Element, m2::Element)
    @assert((m1|>typeof|>nameof == m2|>typeof|>nameof) && (m1|>typeof|>fieldcount == m2|>typeof|>fieldcount == 2),
            "\"*\" error: not implemented between $(m1|>typeof|>nameof) and $(m2|>typeof|>nameof)."
            )
    rawtype(typeof(m1))(m1.value*m2.value, m1.id*m2.id)
end

"""
    -(m::Element) -> typeof(m)
    -(m::Element, ::Nothing) -> typeof(m)
    -(::Nothing, m::Element) -> typeof(m)
    -(m::Element, factor::Number) -> Elements
    -(factor::Number, m::Element) -> Elements
    -(m1::Element, m2::Element) -> Elements
    -(ms::Elements) -> typeof(ms)
    -(ms::Elements, ::Nothing) -> typeof(ms)
    -(::Nothing, ms::Elements) -> typeof(ms)
    -(ms::Elements, factor::Number) -> Elements
    -(factor::Number, ms::Elements) -> Elements
    -(m::Element, ms::Elements) -> Elements
    -(ms::Elements, m::Element) -> Elements
    -(ms1::Elements, ms2::Elements) -> Elements

Overloaded `-` operator between elements of an algebra over a field.
"""
Base.:-(m::Element) = m * (-1)
Base.:-(ms::Elements) = ms * (-1)
Base.:-(m::Element, ::Nothing) = m
Base.:-(::Nothing, m::Element) = -m
Base.:-(ms::Elements, ::Nothing) = ms
Base.:-(::Nothing, ms::Elements) = -ms
Base.:-(factor::Number, m::Element) = rawtype(m|>typeof)(factor, ID()) - m
Base.:-(m::Element, factor::Number) = m + rawtype(m|>typeof)(-factor, ID())
Base.:-(factor::Number, m::ScalarElement) = replace(m, value=factor-m.value)
Base.:-(m::ScalarElement, factor::Number) = replace(m, value=m.value-factor)
Base.:-(m1::ScalarElement, m2::ScalarElement) = replace(m1, value=m1.value-m2.value)
Base.:-(factor::Number, ms::Elements) = rawtype(ms|>valtype)(factor, ID()) - ms
Base.:-(ms::Elements, factor::Number) = ms + rawtype(ms|>valtype)(-factor, ID())
function Base.:-(m1::Element, m2::Element)
    I = promote_type(m1|>idtype, m2|>idtype)
    M = promote_type(typeof(m1), typeof(m2))
    return sub!(Elements{I, M}(m1.id=>m1), m2)
end
function Base.:-(m::Element, ms::Elements)
    I = promote_type(m|>idtype, ms|>keytype)
    M = promote_type(typeof(m), valtype(ms))
    return sub!(Elements{I, M}(m.id=>m), ms)
end
function Base.:-(ms::Elements, m::Element)
    I = promote_type(ms|>keytype, m|>idtype)
    M = promote_type(valtype(ms), typeof(m))
    return sub!(Elements{M|>idtype, M}(ms), m)
end
function Base.:-(ms1::Elements, ms2::Elements)
    I = promote_type(ms1|>keytype, ms2|>keytype)
    M = promote_type(valtype(ms1), valtype(ms2))
    return sub!(Elements{I, M}(ms1), ms2)
end

"""
    /(m::Element, factor::Number) -> Element
    /(m::Element, factor::ScalarElement) -> Element
    /(ms::Elements, factor::Number) -> Elements
    /(ms::Elements, factor::ScalarElement) -> Elements

Overloaded `/` operator for element-scalar division of an algebra over a field.
"""
Base.:/(m::Element, factor::Number) = m * (1/factor)
Base.:/(m::Element, factor::ScalarElement) = m * (1/factor.value)
Base.:/(ms::Elements, factor::Number) = ms * (1/factor)
Base.:/(ms::Elements, factor::ScalarElement) = ms * (1/factor.value)

"""
    //(m::Element, factor::Integer) -> Element
    //(m::Element, factor::ScalarElement) -> Element
    //(ms::Elements, factor::Integer) ->  Elements
    //(ms::Elements, factor::ScalarElement) -> Elements

Overloaded `//` operator for element-scalar division of an algebra over a field.
"""
Base.://(m::Element, factor::Integer) = m * (1//factor)
Base.://(m::Element, factor::ScalarElement) = m * (1//factor.value)
Base.://(ms::Elements, factor::Number) = ms * (1//factor)
Base.://(ms::Elements, factor::ScalarElement) = ms * (1//factor.value)

"""
    ^(m::Element, n::Integer) -> Element
    ^(ms::Elements, n::Integer) -> Elements

Overloaded `^` operator for element-integer power of an algebra over a field.
"""
Base.:^(m::Element, n::Integer) = (@assert n > 0 "^ error: non-positive integers are not allowed."; prod(ntuple(i->m, Val(n))))
Base.:^(ms::Elements, n::Integer) = (@assert n > 0 "^ error: non-positive integers are not allowed."; prod(ntuple(i->ms, Val(n))))

"""
    ⊗(m::Element, ms::Elements) -> Elements
    ⊗(ms::Elements, m::Element) -> Elements
    ⊗(ms1::Elements, ms2::Elements) -> Elements

Overloaded `⊗` operator for element-element multiplications of an algebra over a field.
"""
⊗(m::Element, ms::Elements) = Elements((m ⊗ mm for mm in ms|>values)...)
⊗(ms::Elements, m::Element) = Elements((mm ⊗ m for mm in ms|>values)...)
⊗(ms1::Elements, ms2::Elements) = Elements((m1 ⊗ m2 for m1 in ms1|>values for m2 in ms2|>values)...)

"""
    ⋅(m::Element, ms::Elements) -> Elements
    ⋅(ms::Elements, m::Element) -> Elements
    ⋅(ms1::Elements, ms2::Elements) -> Elements

Overloaded `⋅` operator for element-element multiplications of an algebra over a field.
"""
⋅(m::Element, ms::Elements) = Elements((m ⋅ mm for mm in ms|>values)...)
⋅(ms::Elements, m::Element) = Elements((mm ⋅ m for mm in ms|>values)...)
⋅(ms1::Elements, ms2::Elements) = Elements((m1 ⋅ m2 for m1 in ms1|>values for m2 in ms2|>values)...)

"""
    sequence(m::Element, table=nothing) -> NTuple{rank(m), Int}

Get the sequence of the ids of an element according to a table.
"""
@generated sequence(m::Element, table) = Expr(:tuple, [:(get(table, m.id[$i], nothing)) for i = 1:rank(m)]...)

"""
    split(m::Element) -> Tuple{Any, Vararg{Element}}

Split an element into the coefficient and a sequence of rank-1 elements.
"""
@generated function Base.split(m::Element)
    @assert m|>fieldnames == (:value, :id) "split error: not supported split of $(nameof(typeof(m)))."
    exprs = [:(m.value)]
    for i = 1:rank(m)
        push!(exprs, :(rawtype(typeof(m))(one(m|>valtype), ID(m.id[$i]))))
    end
    return Expr(:tuple, exprs...)
end

"""
    replace(m::Element, pairs::Pair{<:SimpleID, <:Union{Element, Elements}}...) -> Element/Elements
    replace(ms::Elements, pairs::Pair{<:SimpleID, <:Union{Element, Elements}}...) -> Elements

Replace the rank-1 components of an element with new element/elements.
"""
function Base.replace(m::Element, pairs::Pair{<:SimpleID, <:Union{Element, Elements}}...)
    @assert m|>typeof|>fieldnames == (:value, :id) "replace error: not supported replacement of $(nameof(typeof(m)))."
    replacedids = NTuple{length(pairs), idtype(m)|>eltype}(pair.first for pair in pairs)
    ms = split(m)
    result = ms[1]
    for i = 1:rank(m)
        index = findfirst(isequal(m.id[i]), replacedids)
        result = result * (isa(index, Int) ? pairs[index].second : ms[i+1])
    end
    return result
end
function Base.replace(ms::Elements, pairs::Pair{<:SimpleID, <:Union{Element, Elements}}...)
    result = elementstype(pairs...)()
    for m in values(ms)
        add!(result, replace(m, pairs...))
    end
    return result
end
@generated function elementstype(ms::Vararg{Pair{<:SimpleID, <:Union{Element, Elements}}, N}) where N
    mts = ntuple(i->(fieldtype(ms[i], 2) <: Elements) ? fieldtype(ms[i], 2)|>valtype : fieldtype(ms[i], 2), Val(N))
    its = ntuple(i->(fieldtype(ms[i], 2) <: Elements) ? fieldtype(ms[i], 2)|>keytype : fieldtype(ms[i], 1), Val(N))
    @assert mapreduce(mt->(mt|>fieldnames) == (:value, :id), &, mts) "elementstype error: not supported."
    I = mapreduce(eltype, promote_type, its)
    V = mapreduce(valtype, promote_type, mts)
    M = reduce(typejoin, mts)|>rawtype
    @assert isconcretetype(V) @sprintf "elementstype error: not concrete type of coefficients(%s)." V
    return Elements{ID{I}, M{V, <:ID{I}}}
end

"""
    permute!(result::Elements, m::Element, table=nothing) -> Elements
    permute!(result::Elements, ms::Elements, table=nothing) -> Elements

Permute the ids of an-element/a-set-of-elements to the descending order according to a table, and store the permuted elements in result.

!!! note
    To use this function, the user must implement a method of `permute`, which computes the result of the permutation of a rank-2 element and takes the following interface:
    ```julia
    permute(::Type{M}, id1::SimpleID, id2::SimpleID, table=nothing) -> Union{M, Elements}
    ```
    Here, `M` is the type of the input element `m` in `permute!`, `id1` and `id2` are two arbitary simple ids contained in `m.id`.
"""
function Base.permute!(result::Elements, m::Element, table=nothing)
    cache = valtype(result)[m]
    while length(cache) > 0
        current = pop!(cache)
        pos = elementcommuteposition(sequence(current, table))
        if isa(pos, Nothing)
            add!(result, current)
        else
            left, right = elementleft(current, pos), elementright(current, pos)
            for middle in permute(typeof(current), current.id[pos], current.id[pos+1], table)
                temp = left * middle * right
                (temp === nothing) || push!(cache, temp)
            end
        end
    end
    return result
end
function Base.permute!(result::Elements, ms::Elements, table=nothing)
    for m in values(ms)
        permute!(result, m, table)
    end
    return result
end
function elementcommuteposition(seqs)
    pos = 1
    while pos < length(seqs)
        (seqs[pos] < seqs[pos+1]) && return pos
        pos += 1
    end
    return nothing
end
elementleft(m::Element, i::Int) = rawtype(typeof(m))(m.value, m.id[1:i-1])
elementright(m::Element, i::Int) = rawtype(typeof(m))(one(m.value), m.id[i+2:end])

"""
    permute(m::Element, table=nothing) -> Elements
    permute(ms::Elements, table=nothing) -> Elements

Permute the ids of an-element/a-set-of-elements to the descending order according to a table.
"""
function permute(m::Element, table=nothing)
    M, V, I = m|>typeof|>rawtype, m|>valtype, m|>idtype|>eltype
    permute!(Elements{ID{I}, M{V, <:ID{I}}}(), m, table)
end
function permute(ms::Elements, table=nothing)
    M, V, I = ms|>valtype|>rawtype, ms|>valtype|>valtype, ms|>keytype|>eltype
    permute!(Elements{ID{I}, M{V, <:ID{I}}}(), ms, table)
end

end #module


text:
using DEBmicroTrait
using CSV, DataFrames, Statistics
using DifferentialEquations

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
module GindeleRandomForestBehaviors

using AutomotiveDrivingModels

using RandomForestBehaviors.MvDecisionTrees

import AutomotiveDrivingModels:
    AbstractVehicleBehavior,
    AbstractVehicleBehaviorPreallocatedData,
    AbstractVehicleBehaviorTrainParams,

    preallocate_learning_data,
    select_action,
    calc_action_loglikelihood,
    train

export
    GindeleRandomForestBehavior,
    GRF_TrainParams,

    select_action,
    calc_action_loglikelihood,
    train

type GRF_TrainParams <: AbstractVehicleBehaviorTrainParams

    targets::ModelTargets
    indicators::Vector{AbstractFeature}

    ntrees::Int # number of trees in the forest
    max_tree_depth::Int # max depth for trees in the forest
    min_samples_split::Int # min num samples in a node for split to be considered
    min_samples_leaves::Int # min num samples in resulting leaves for split to be done
    n_split_tries::Int # number of random split attempts
    use_PCA::Bool

    min_split_improvement::Float64 # minimum improvement for split to be accepted
    partial_sampling::Float64 # percentage of dataset to sample

    function GRF_TrainParams(;
        targets::ModelTargets = ModelTargets(Features.FUTUREDESIREDANGLE, Features.FUTUREACCELERATION),
        indicators::Vector{AbstractFeature} = AbstractFeature[],

        ntrees::Integer=10,
        max_tree_depth::Integer=5,
        min_samples_split::Integer=100,
        min_samples_leaves::Integer=20,
        n_split_tries::Integer=10,
        use_PCA::Bool = false,
        min_split_improvement::Float64=0.0,
        partial_sampling::Float64=0.7,
        )

        retval = new()
        retval.targets = targets
        retval.indicators = indicators
        retval.ntrees = ntrees
        retval.max_tree_depth = max_tree_depth
        retval.min_samples_split = min_samples_split
        retval.min_samples_leaves = min_samples_leaves
        retval.n_split_tries = n_split_tries
        retval.use_PCA = use_PCA
        retval.min_split_improvement = min_split_improvement
        retval.partial_sampling = partial_sampling
        retval
    end
end
type GRF_PreallocatedData <: AbstractVehicleBehaviorPreallocatedData

    # the full dataset matrices pre-extracted
    X::Matrix{Float64} # [nfeatures × nframes]
    Y::Matrix{Float64} # [2 × nframes]
    Z::Matrix{Float64} # [n_pca_features × nframes]

    extractor::FeatureSubsetExtractor
    preprocess::DataPreprocessor
    preprocess_pca::DataPreprocessor
    action_clamper::DataClamper

    function GRF_PreallocatedData(dset::ModelTrainingData2, params::GRF_TrainParams)
        retval = new()

        targets = params.targets
        indicators = params.indicators
        trainingframes = dset.dataframe
        nframes = nrow(trainingframes)
        nindicators = length(indicators)

        X = Array(Float64, nindicators, nframes) # NOTE(tim): initially indicators × nframes, we later transpose it
        Y = Array(Float64, 2, nframes)
        pull_design_and_target_matrices!(X, Y, trainingframes, targets, indicators)

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

        extractor = FeatureSubsetExtractor(Array(Float64, nindicators), indicators)
        preprocess = ChainedDataProcessor(extractor)
        push!(preprocess, X, DataNaReplacer, extractor.indicators)
        process!(X, preprocess.processors[end]) # process in place
        push!(preprocess, X, DataScaler) # ensure data is standardized before running PCA
        process!(X, preprocess.processors[end]) # process in place
        retval.X = copy(X)

        # Add a PCA scaler
        preprocess_pca = deepcopy(preprocess)
        n_PCA_features = min(nindicators, 20)
        push!(preprocess_pca, X, DataLinearTransform, n_PCA_features)

        Z = Array(Float64, n_PCA_features, nframes)
        retval.Z = process!(Z, X, preprocess_pca.processors[end])

        retval.Y = Y
        retval.extractor = extractor
        retval.preprocess = preprocess
        retval.preprocess_pca = preprocess_pca
        retval.action_clamper = DataClamper(
                Array(Float64, 2),
                vec(minimum(Y, 2)),
                vec(maximum(Y, 2))
            )

        retval
    end
end
type GindeleRandomForestBehavior <: AbstractVehicleBehavior
    model_μ::Ensemble
    model_Σ::Ensemble

    targets::ModelTargets
    indicators::Vector{AbstractFeature}
    extractor::FeatureSubsetExtractor
    processor::DataPreprocessor
    action_clamper::DataClamper

    # preallocated memory
    μ::Vector{Float64} # [2]
    Σ::Matrix{Float64} # [2×2]

    function GindeleRandomForestBehavior(
        model_μ::Ensemble,
        model_Σ::Ensemble,
        params::GRF_TrainParams,
        pre::GRF_PreallocatedData,
        )

        extractor = pre.extractor
        action_clamper = pre.action_clamper
        processor = params.use_PCA > 0 ? pre.preprocess_pca : pre.preprocess

        retval = new()
        retval.model_μ = model_μ
        retval.model_Σ = model_Σ

        retval.targets = params.targets
        # copy the processor and ensure that the copied io arrays match up
        retval.extractor = FeatureSubsetExtractor(deepcopy(extractor.x), extractor.indicators)
        retval.processor = deepcopy(processor, retval.extractor)
        retval.action_clamper = DataClamper(Array(Float64, 2),
                                   deepcopy(action_clamper.f_lo),
                                   deepcopy(action_clamper.f_hi))

        retval.μ = Array(Float64, 2)
        retval.Σ = Array(Float64, 2, 2)

        retval
    end
end
AutomotiveDrivingModels.get_targets(behavior::GindeleRandomForestBehavior) = behavior.targets
function AutomotiveDrivingModels.get_indicators(behavior::GindeleRandomForestBehavior)

    # run through and pull the indicators that are actually used in the tree

    indicators = Set{AbstractFeature}()

    for index in get_split_index_set(behavior.model_μ)
        push!(indicators, behavior.extractor.indicators[index])
    end
    for index in get_split_index_set(behavior.model_Σ)
        push!(indicators, behavior.extractor.indicators[index])
    end

    return collect(indicators)
end

function preallocate_learning_data(
    dset::ModelTrainingData2,
    params::GRF_TrainParams)

    GRF_PreallocatedData(dset, params)
end

function _calc_mvnormal(behavior::GindeleRandomForestBehavior)

    apply_forest!(behavior.μ, behavior.model_μ, behavior.processor.z)
    apply_forest!(behavior.Σ, behavior.model_Σ, behavior.processor.z)

    MvNormal(behavior.μ, behavior.Σ)
end

function select_action(
    behavior::GindeleRandomForestBehavior,
    runlog::RunLog,
    sn::StreetNetwork,
    colset::UInt,
    frame::Int
    )

    observe!(behavior.extractor, runlog, sn, colset, frame)
    process!(behavior.processor)

    normal = _calc_mvnormal(behavior)

    Distributions._rand!(normal, behavior.action_clamper.x)
    process!(behavior.action_clamper)

    action_lat = behavior.action_clamper.x[1]
    action_lon = behavior.action_clamper.x[2]

    (action_lat, action_lon)
end

function calc_action_loglikelihood(
    behavior::GindeleRandomForestBehavior,
    runlog::RunLog,
    sn::StreetNetwork,
    colset::UInt,
    frame::Int,
    action_lat::Float64,
    action_lon::Float64,
    )

    #=
    Compute the log-likelihood of the action taken during a single frame
    given the VehicleBehaviorGaussian.
    =#

    behavior.action_clamper.x[1] = action_lat
    behavior.action_clamper.x[2] = action_lon

    observe!(behavior.extractor, runlog, sn, colset, frame)
    process!(behavior.processor)

    normal = _calc_mvnormal(behavior)
    logpdf(normal, behavior.action_clamper.x)
end
function calc_action_loglikelihood(
    behavior::GindeleRandomForestBehavior,
    features::DataFrame,
    frameind::Integer,
    )

    behavior.action_clamper.x[1] = features[frameind, symbol(behavior.targets.lat)]::Float64
    behavior.action_clamper.x[2] = features[frameind, symbol(behavior.targets.lon)]::Float64

    observe!(behavior.extractor, features, frameind)
    process!(behavior.processor)

    apply_forest!(behavior.μ, behavior.model_μ, behavior.processor.z)
    apply_forest!(behavior.Σ, behavior.model_Σ, behavior.processor.z)

    mvnorm_model = MvNormal(behavior.μ, behavior.Σ)
    logpdf(mvnorm_model, behavior.action_clamper.x)
end

function train(
    training_data::ModelTrainingData2,
    preallocated_data::GRF_PreallocatedData,
    params::GRF_TrainParams,
    foldset::FoldSet,
    )

    ntrees = params.ntrees
    max_tree_depth = params.max_tree_depth
    min_samples_split = params.min_samples_split
    min_samples_leaves = params.min_samples_leaves
    n_split_tries = params.n_split_tries
    min_split_improvement = params.min_split_improvement
    partial_sampling = params.partial_sampling

    #######

    build_tree_params = BuildTreeParameters(
        ntrees, max_tree_depth, min_samples_split, min_samples_leaves,
        min_split_improvement, n_split_tries, LossFunction_MSE, MeanVecLeaf)

    X_orig = params.use_PCA ? preallocated_data.Z : preallocated_data.X

    Y_orig = preallocated_data.Y

    nframes = length(foldset)
    X = Array(Float64, nframes, size(X_orig, 1)) # NOTE: this is transpose of X_orig
    y = Array(Float64, 2, nframes)

    i = 0
    for frame in foldset
        i += 1

        for j in 1 : size(X_orig, 1)
            X[i,j] = X_orig[j,frame]
        end

        y[1,i] = Y_orig[1,frame]
        y[2,i] = Y_orig[2,frame]
    end
    @assert(i == nframes)

    μ = build_forest(y, X, ntrees, build_tree_params, partial_sampling)

    mean_vector = Array(Float64, 2)
    for i = 1 : size(y, 1)
        y[:,i] -= apply_forest!(mean_vector, μ, vec(X[i,:])) # subtract the predicted mean
    end

    build_tree_params = BuildTreeParameters(
        ntrees, max_tree_depth, min_samples_split, min_samples_leaves,
        min_split_improvement, n_split_tries, LossFunction_MSE, CovLeaf)

    Σ = build_forest(y, X, ntrees, build_tree_params, partial_sampling)

    GindeleRandomForestBehavior(μ, Σ, params, preallocated_data)
end

end # end module


text:
## Modification to Mustache and friends to make this easier to work with

## Image substitution: LaTeX, Plot, File
## in Markdown ![alt](url) will display a graphic image by url
## We use mustache and a few files to add other images:
## The pattern: ![alt]({{{:img}}}) can hold
## * a Plots plot, `p`, when `img=Plot(p)`
## * a local file when `img=File(file_name)`.
## * LaTeX includes the ![alt](...) in it
## For each, the image is embedded as a Base64-encode string

# Plots. Workflow
# template include: ![Alt]({{{:figure}}})
# use `figure=Plot(p)` in ther context
# Also File(fname) and LaTeX(snippet) produce figures
function Plot(p)
    io  = IOBuffer()
    show(io, MIME("image/png"), p)
    data = Base64.base64encode(take!(io))
    close(io)

    io = IOBuffer()
    print(io,"data:image/gif;base64,")
    print(io,data)
    String(take!(io))
end

## File ![alt]({{{:img}}})
## where `img = File("imagefile.png"))` is in context
function File(p)
    img = base64encode(read(p, String))
    io = IOBuffer()
    print(io,"data:image/gif;base64,")
    print(io,img)
    String(take!(io))
end

## write p as tikz picture for inclusion into document
## requires pgfplotsx() in Plots
function Tikz(tikzfile)
    tzio = IOBuffer()
    for l in readlines(tikzfile)
        println(tzio, replace(l, ", show background rectangle"=>""))
    end
    String(take!(tzio))
end

abstract type LATEXTemplate end
function question(io::IO, qt::Symbol, q::LATEXTemplate, answers...)
    question(io, qt, LaTeX(q.string[1]), answers...)
end

# indicate latex for conversion to latex
struct LaTeXTemplate <: LATEXTemplate
    tokens
    string::Vector
end

(tpl::LaTeXTemplate)(;kwargs...) = (tpl.string[1] = tpl.tokens(;kwargs...); tpl)

macro ltx_str(s)
    tokens = Mustache.parse(s)
    str = sprint(io -> Mustache.render(io, tokens))
    LaTeXTemplate(tokens, [str])
end

# inicate markdown for conversion to LaTeX
struct MDLaTeXTemplate <: LATEXTemplate
    tokens
    string::Vector
end

(tpl::MDLaTeXTemplate)(;kwargs...) = begin
    str = tpl.tokens(;kwargs...)
    tpl.string[1] = _LaTeX′(str)
    tpl
end

macro mdltx_str(s)
    tokens = Mustache.parse(s)
    str = _LaTeX′(tokens())
    MDLaTeXTemplate(tokens, [str])
end

"""
    MT

Use `<<...>>` or `<<{...}>>` for substitution before randomimization substitution. Useful for plots
"""
macro MT_str(s)
    Mustache.parse(s, ("<<", ">>"))
end

## Templates

# used to include LaTeX into MarkDown
const LaTeX_tpl = mt"""![{{{:alt}}}]({{{:latex}}})"""

# used by tth
const tth_tpl = mt"""
\documentclass{article}
\usepackage{amssymb}
\usepackage{amsmath}
{{#:pkgs}}
\usepackage{ {{.}} }
{{/:pkgs}}

\begin{document}
\thispagestyle{empty}
{{{:txt}}}
\end{document}
"""

# use by LaTeX and preview to create pdf
# uses preview package to frame output
const preview_tpl = mt"""
\documentclass{article}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{float}
\usepackage{ragged2e}
\usepackage{pgfplots}
\pgfplotsset{compat=1.16}
\usepgfplotslibrary{fillbetween}
\usepackage{tikz}

{{#:pkgs}}
\usepackage{ {{.}} }
{{/:pkgs}}

\usepackage[active, tightpage]{preview}
\setlength\PreviewBorder{10pt}%

\begin{document}
\thispagestyle{empty}
\begin{preview}
\{{:fontsize}}
{{{:txt}}}
\end{preview}
\end{document}
"""

## LaTeX to [HTML, png]
# convert latex to png file
latex_to_png(str::Mustache.MustacheTokens;kwargs...) = latex_to_png(str(); kwargs...)
function latex_to_png(str;  fontsize="LARGE", tpl=preview_tpl, pkgs=[])

    fnm = tempname()
    fnmtex = fnm * ".tex"
    fnmpdf = fnm * ".pdf"
    fnmpng = fnm * ".png"

    open(fnmtex, "w") do io
        Mustache.render(io, tpl, (txt=str, fontsize=fontsize, pkgs=collect(pkgs)))
    end
    tectonic() do bin
        run(`$bin -c minimal $fnmtex`)
    end

    run(`cp $fnmtex /tmp/test.tex`)

    if Sys.isapple()
        run(pipeline(`sips -s format png $fnmpdf --out $fnmpng`, stdout=devnull))
    else
        imagemagick_convert() do bin
            run(`$bin -density 150 -depth 8 -quality 100 $fnmpdf $fnmpng`)
        end
    end

    return fnmpng

    # Keep in case useful elsewhere; no longer used as with preview no need to trim
    # hack to trim pdf and write out png
    # function _findfirst(rs)
    #     i = findfirst(iszero, rs)
    #     i == nothing ? 1 : i
    # end
    # function _findlast(rs)
    #     i = findlast(iszero, rs)
    #     i == nothing ? length(rs) : i
    # end
    # WHITE = ImageMagick.RGBA{ImageMagick.N0f16}(1.0,1.0,1.0,0.0)
    # rs = [all(a[i,:] .== WHITE) for i in 1:size(a, 1)]
    # cs = [all(a[:,j] .== WHITE) for j in 1:size(a, 2)]
    # rₘ = _findfirst(rs)
    # rₙ = _findlast(rs)
    # cₘ = _findfirst(cs)
    # cₙ = _findlast(cs)
    # b = a[rₘ:rₙ, cₘ:cₙ]
    # save(fnmpng, b, quality=100)
    # return fnmpng

end

const FONTSIZE = "large"

"""
    LaTeX(str; alt="", tpl=[preview_tpl], fontsize="$FONTSIZE", pkgs=[])

Run LaTeX and produce output for inclusion as a question.

Returns string for a Markdown figure, e.g,
`"![alt](latex_as_encoded_png_file)"` that gets rendered into an image.

* `alt`: alt tag for the image

* `tpl`: what template to use, default uses `preview` package.

* `fontsize`: what fontsize to use, among ( "Huge", "huge", "LARGE", "Large", "large",
"normalsize", "small", "footnotesize", "scriptsize", "tiny)

* `pkgs` a collection of package names for inclusion through `\\usepackage{pkgname}`

To get a png file from a LaTeX snippet, `str`, we have

```
JuliaBlackBoard.latex_to_png(str)
```

The `preview` function can be used to preview the image generated.

"""
function LaTeX(str;
               alt = "",
               tpl=preview_tpl,
               fontsize=FONTSIZE,
               pkgs=[])

    fnm = latex_to_png(str, fontsize=fontsize, tpl=tpl, pkgs=pkgs)
    LaTeX_tpl(latex=File(fnm), alt=create_html(alt, strip=true))

end

# show the png file generated
preview(mt::Mustache.MustacheTokens; kwargs...) = preview(mt(); kwargs...)
preview(ltpl::LATEXTemplate ; kwargs...) = preview(ltpl.tokens; kwargs...)
preview(ltpl::MDLaTeXTemplate ; kwargs...) = preview(_LaTeX′(ltpl.tokens()); kwargs...)
function preview(str; tpl=preview_tpl, fontsize=FONTSIZE, pkgs=[])
    fnm = tempname()
    fnmtex = fnm * ".tex"
    open(fnmtex, "w") do io
        Mustache.render(io, tpl, (txt=str, fontsize=fontsize, pkgs=collect(pkgs)))
    end
    run(`cp $fnmtex /tmp/str.tex`)
    tectonic() do bin
        run(`$bin --print $fnmtex`)
    end
    run(`open $fnm.pdf`)

end

"""
    LaTeX′(str; kwargs...)

Parse `Markdown` string `str` and *then* run `LaTeX` function.

Allows LaTeX questions to use Markdown formatting.
"""
LaTeX′(str::Mustache.MustacheTokens; kwargs...) = LaTeX′(str(); kwargs...)
function LaTeX′(str; kwargs...)
    ast = parser(str)
    out =  sprint(io -> show(io, "text/latex", ast))
    LaTeX(out; kwargs...)
end

function _LaTeX′(str)
    ast = parser(str)
    sprint(io -> show(io, "text/latex", ast))
end

# Take a string of LaTeX code and produce an HTML fragment with tth
function latex_to_html(ltx)

    fnm = tempname() * ".tex"
    open(fnm, "w") do io
        Mustache.render(io, tth_tpl, (txt = ltx, pkgs=[]))
    end

    out = try
        io = IOBuffer()
        tth() do tthbinary
            run(pipeline(`$tthbinary -f5 -i -r  -t -w`, stdin=fnm, stdout=io, stderr=devnull))
        end
        out = String(take!(io))
        out = out[11:end-1]
        out
    catch err
        @warn "Issue running tth on ≪ $ltx ≫"
        ltx
    end

    out

end

# create HTML from string or mustache tokens
create_html(q::Mustache.MustacheTokens; kwargs...) = create_html(q();kwargs...)
function create_html(q; strip=false)

    ast = parser(q)
    qq =  sprint(io -> show(io, "text/html", ast))
    qq = chomp(qq)
    qq = replace(qq, "\n" => " ")
    qq = replace(qq, "\t" => " "^4)
    if strip
        qq = qq[4:end-4]
    end
    qq

end

## IO helpers for managing pools within a single file
## see example-pool.jl for usage
# OPEN = JuliaBlackBoard.OPEN(@__FILE__)
# POOL = JuliaBlackBoard.POOL(@__FILE__)
function OPEN(SCRIPTNAME)
    bnm = replace(SCRIPTNAME, r".jl$" => "")
    fname = bnm * ".txt"
    (f) -> open(f, fname, "w")
end

function POOL(SCRIPTNAME)
    bnm = replace(SCRIPTNAME, r".jl$" => "")
    ctr = Ref(0)
    (f, io) -> begin
        num = ctr[] = ctr[] + 1
        question(io, OP, "XXX replace with pool-$num XXX")
        poolnm = bnm * "-pool-" * string(num) * ".txt"
        @show poolnm
        open(f, poolnm, "w")
    end
end

# Use CommonMark -- not Markdown-- for parsing, as we can
# overrule the latex bit easier

# set up common mark parser
parser = Parser()
enable!(parser, MathRule())
enable!(parser, DollarMathRule())

## Math display
## There are two ways to display math marked up in LaTeX within the question's HTML:
## * the default is to use `TtH` to convert LaTeX in to HTML. The drawback is some
##   things convert poorly, such as \bar{x} and \hat{x} type constructs.
## * Setting `ENV["USE_MATHJAX"] = true` will allow BlackBoard to render the formula
##   using MathML. This looks *much* better. **HOWEVER**, to see this one must edit and
##   save the question within BlackBoard.
##
## By overwriting methods from CommonMark  non-fatal warnings are issued when precompiling
function CommonMark.write_html(::CommonMark.Math, rend, node, enter)
    if haskey(ENV, "USE_MATHJAX") && ENV["USE_MATHJAX"] == "true"
        print(rend.buffer, "\$\$" * node.literal * "\$\$")
    else
        tmp = latex_to_html("\$" * node.literal * "\$")
        CommonMark.tag(rend, "span", CommonMark.attributes(rend, node, ["class" => "math"]))
        print(rend.buffer, latex_to_html("\$" * node.literal * "\$"))
        CommonMark.tag(rend, "/span")
    end
end

function CommonMark.write_html(::CommonMark.DisplayMath, rend, node, enter)
    if haskey(ENV, "USE_MATHJAX") && ENV["USE_MATHJAX"] == "true"
        CommonMark.tag(rend, "div")##, CommonMark.attributes(rend, node, ["class" => "display-math"]))
        print(rend.buffer,  "\$\$" * node.literal* "\$\$
")
        CommonMark.tag(rend, "/div")
    else
        CommonMark.tag(rend, "div")##, CommonMark.attributes(rend, node, ["class" => "display-math"]))
        print(rend.buffer,  latex_to_html("\$\$" * node.literal* "\$\$"))
        CommonMark.tag(rend, "/div")
    end
end

## enhance BlackBoards display of code
const CODE_FACE = "Courier New"
function CommonMark.write_html(::CommonMark.Code, r, n, ent)
    CommonMark.tag(r, "code", CommonMark.attributes(r, n))
    CommonMark.tag(r, "font", ["face" => CODE_FACE])
    CommonMark.literal(r, CommonMark.escape_xml(n.literal))
    CommonMark.tag(r, "/font")
    CommonMark.tag(r, "/code")
end

function CommonMark.write_html(::CommonMark.CodeBlock, r, n, ent)

    nliteral = CommonMark.escape_xml(n.literal)
    nliteral = replace(nliteral, "\n" => "
")

    info_words = split(n.t.info === nothing ? "" : n.t.info)
    attrs = CommonMark.attributes(r, n)
    if !isempty(info_words) && !isempty(first(info_words))
        push!(attrs, "class" => "language-$(escape_xml(first(info_words)))")
    end
    CommonMark.cr(r)
    CommonMark.tag(r, "pre")
    CommonMark.tag(r, "code", attrs)
    CommonMark.tag(r, "font", ["face" => CODE_FACE])
    CommonMark.literal(r, nliteral)
    CommonMark.tag(r, "/font")
    CommonMark.tag(r, "/code")
    CommonMark.tag(r, "/pre")
    CommonMark.literal(r, "
")
    CommonMark.cr(r)
end

# write_latex adjustments
function CommonMark.write_latex(image::CommonMark.Image, w, node, ent)
    if ent
        image = CommonMark._smart_link(MIME"text/latex"(), image, node, w.env)
        CommonMark.cr(w)
        CommonMark.literal(w, "\\begin{figure}[H]\n") # add [H] and use `float` pkg
        CommonMark.literal(w, "\\centering\n")
        CommonMark.literal(w, "\\includegraphics{", image.destination, "}\n")
        CommonMark.literal(w, "\\caption{")
    else
        CommonMark.literal(w, "}\n")
        CommonMark.literal(w, "\\end{figure}")
        CommonMark.cr(w)
    end
end

#function CommonMark.write_latex(::CommonMark.Backslash, w, node, ent)
#    CommonMark.literal(w, ent ? raw"" : raw"")
#end

## Adjusments for `latex_escape`ing
let chars = Dict(
    ##        '^'  => "\\^{}",
    ##        '\\' => "{\\textbackslash}", # modify this one
    '~'  => "{\\textasciitilde}",
)
    for c in    "&%#[]" # took out _ \$, {} added []
    chars[c] = "\\$c"
    end
    global function CommonMark.latex_escape(w::CommonMark.Writer, s::AbstractString)
        for ch in s
            CommonMark.literal(w, get(chars, ch, ch))
        end
    end
end


text:
module Media

include("dynamic.jl")
include("system.jl")
include("compat.jl")

function __init__()
  init_compat()
end

end


text:
export Executable

typealias Executable hsa_executable_t

function Executable(;
    profile = HSA_PROFILE_FULL,
    executable_state = HSA_EXECUTABLE_STATE_UNFROZEN,
    options = C_NULL
    )

    res = Ref{hsa_executable_t}(hsa_executable_t(0))

    err = ccall((:hsa_executable_create, libhsa), hsa_status_t, (
    hsa_profile_t,
    hsa_executable_state_t,
    Ptr{Cchar},
    Ptr{hsa_executable_t}),
    profile, executable_state, options, res)

    test_status(err)

    return res.x
end

function executable_load_code_object(
    executable, agent, code_object;
    options = C_NULL)

    err = ccall((:hsa_executable_load_code_object, libhsa), hsa_status_t, (
    hsa_executable_t,
    hsa_agent_t,
    hsa_code_object_t,
    Ptr{Cchar}),
    executable, agent, code_object, options)

    test_status(err)
end

function executable_freeze(
    executable;
    options = C_NULL)

    err = ccall((:hsa_executable_freeze, libhsa), hsa_status_t, (
    hsa_executable_t,
    Ptr{Cchar}),
    executable,
    options)

    test_status(err)
end

function executable_get_symbol(
    executable,
    symbol_name;
    module_name = "",
    agent = hsa_agent_t(0),
    call_convention = 0)

    if module_name == "" || module_name == nothing
        module_name = C_NULL
    end

    res = Ref{hsa_executable_symbol_t}(hsa_executable_symbol_t(0))

    err = ccall((:hsa_executable_get_symbol, libhsa), hsa_status_t, (
    hsa_executable_t,
    Ptr{UInt8},
    Ptr{UInt8},
    hsa_agent_t,
    Int32,
    Ptr{hsa_executable_symbol_t}),
    executable, module_name, symbol_name, agent, call_convention, res)

    test_status(err)

    return res.x
end

const executable_iterate_symbols_cb_ptr = cfunction(iterate_cb, hsa_status_t, (hsa_executable_t, hsa_executable_symbol_t, Ptr{Void}))

function executable_iterate_symbols(executable, cb::Function)
    convert_cb = function(ex, sym)
        name = HSA.executable_symbol_info_name(sym)
        cb(ex, sym, name)
    end
    state = IterState(convert_cb)

    state_ptr = pointer_from_objref(state)

    err = ccall((:hsa_executable_iterate_symbols, libhsa), hsa_status_t, (hsa_executable_t, Ptr{Void}, Ptr{Void}),
                executable, executable_iterate_symbols_cb_ptr, state_ptr)

    if !isnull(state.err)
        rethrow(state.err.value)
    end

    test_status(err)
end

function symbols(executable::Executable)
    res = []

    executable_iterate_symbols(executable, (ex, sym, name) -> push!(res, (sym, name)))

    return res
end


text:
module mpec

    using Ipopt
    using JuMP
    using DataFrames
    using DataFramesMeta

    """
        build_trans(p::Vector{Float64},n::Int)

    Build transition matrix for each of `n` states where each element `p[j]` of `p` is the probability of moving `j-1` states up in the grid defined by `1:n`.
    """
    function build_trans(p::Vector{Float64},n::Int)
        P = spzeros(n,n)
        p = vcat(p,(1-sum(p)))
        N = length(p)
        for i=0:(N-1)
            P =   P + sparse(1:(n-i),(1+i):n,ones(n-i)*p[i+1],n,n)
            P[n-i,n] = 1-sum(p[1:i])
        end
        return P
    end

    type Param
        n :: Int  # num of grid points
        gmax :: Int # upper boudn of mileage grid

        # structural params
        pr :: Vector{Float64}  # transition probs
        RC :: Float64         # replacement cost
        c :: Float64  # cost parameter
        beta :: Float64 # discount factor

        # computation objects
        grid :: LinSpace{Float64}
        P :: SparseMatrixCSC{Float64,Int64}
        # EV0 :: Vector{Float64}   # old value function
        # EV1 :: Vector{Float64}   # new
        pk :: Vector{Float64}  # probability of keeping engine
        cost :: Vector{Float64}   #operating cost

        """
            create default param type
        """
        function Param()
            this = new()
            this.n = 175
            this.gmax = 450
            this.grid = 0.0:(this.n-1)

            this.pr = [0.0937; 0.4475; 0.4459; 0.0127]
            this.RC = 11.7257
            this.c = 2.45569
            this.beta = 0.9999
            this.P = build_trans(this.pr,this.n)
            this.cost = 0.001*this.c*collect(this.grid)
            # this.EV0 = zeros(this.n)
            # this.EV1 = zeros(this.n)
            this.pk = zeros(this.n)
            return this
        end
    end

    function read_busses(p::Param,num_types=4)
        # load data
        d = readtable(joinpath(dirname(@__FILE__),"buses.csv"),header=false)
        # setup data
        rename!(d,[:x1,:x2,:x5,:x7,:x9],[:id,:bus_type,:d1,:x,:dx1])

        # subset to bus type
        d = @where(d,:bus_type .<= num_types)

        # discretize odometer data
        d[:x] = ceil(Int,p.n/(p.gmax*1000) * d[:x])
        d[:dx1] = d[:x] .- vcat(0,d[:x][1:(end-1)])
        d[:dx1] = d[:dx1].*(1-d[:d1]) .+d[:d1].*d[:x]  # replace first diff of x by x if replaced

        # get replacement dummy
        d[:d] = vcat(d[:d1][2:end],0)

        # remove obs with missing lagged mileage
        # this is the first row for each obs
        remove = d[:id] .- vcat(0,d[:id][1:(end-1)])
        d = d[remove.==0,:]

        d = @select(d,:id,:bus_type,:d,:x,:dx1)

        return d
    end

end

max_stars_repo_name: snowdj/CoursePack

max_issues_repo_name: snowdj/CoursePack

max_forks_repo_name: snowdj/CoursePack