text:
#=
    ----------------------------------------------------------------------------
    Implements the main functionality of the package: Markov chain Monte Carlo
    sampler, which allows for exploration of a product space between a path
    space (on which the diffusion of interest is defined) and a parameter space.
    The main function is `mcmc(...)` and it is supported by two main types of
    routines: `impute!(...)` and `update_param(...)`.
    ----------------------------------------------------------------------------
=#

import ForwardDiff: value # currently not needed, but will be

#===============================================================================
                                Main routine
===============================================================================#
"""
    mcmc(setup)

Gibbs sampler alternately imputing unobserved parts of the path and updating
unknown coordinates of the parameter vector. `setup` defines all variables of
the Markov Chain
"""
function mcmc(setup::MCMCSetup)
    adaptive_prop, num_mcmc_steps = setup.adaptive_prop, setup.num_mcmc_steps
    ws, ll, θ = Workspace(setup)
    gibbs = GibbsDefn(setup)
    init_adaptation!(adaptive_prop, ws)

    for i in 1:num_mcmc_steps
        verbose = act(Verbose(), ws, i)
        act(SavePath(), ws, i) && save_path!(ws)
        ws = next_set_of_blocks(ws)
        ll, acc = impute!(ws, ll, verbose, i)
        update!(ws.accpt_tracker, Imputation(), acc)
        act(Readjust(), ws, i) && readjust_pCN!(ws, i)

        if act(ParamUpdate(), ws, i)
            for j in 1:length(gibbs)
                ll, acc, θ = update_param!(gibbs[j], θ, ws, ll, verbose, i)
                update!(ws.accpt_tracker, ParamUpdate(), j, acc)
                update!(ws.θ_chain, θ)
                verbose && print("\n")
            end
            verbose && print("------------------------------------------------",
                             "------\n")
        end
        add_path!(adaptive_prop, ws.XX, i)
        #print_adaptation_info(adaptive_prop, accImpCounter, accUpdtCounter, i)
        adaptive_prop, ll = update!(adaptive_prop, ws, i, ll)
        adaptive_prop = still_adapting(adaptive_prop)
    end
    display_acceptance_rate(ws.blocking)
    ws
end

#===============================================================================
                            Imputation routines       (some also for param updt)
===============================================================================#

"""
    impute!(ws::Workspace{OS,NoBlocking}, ll, verbose=false, it=NaN,
            headstart=false) where OS

Imputation step of the MCMC scheme (without blocking).
...
# Arguments
- `ws`: workspace containing most containers and definitions of the chain
- `ll`: log-likelihood of the old (previously accepted) diffusion path
- `verbose`: whether to print update info while sampling
- `it`: iteration index of the MCMC algorithm
- `headstart`: flag for whether to 'ease into' fpt conditions
...
"""
function impute!(ws::Workspace{OS,NoBlocking}, ll, verbose=false, it=NaN,
                 headstart=false) where OS
    ρ = ws.ρ[1][1]
    # sample proposal starting point
    zᵒ, yᵒ = proposal_start_pt(ws, ws.P[1], ρ)

    # sample proposal path
    m = length(ws.WWᵒ)
    sample_segments!(1:m, ws, yᵒ, ρ, Val{headstart}())

    llᵒ = logpdf(ws.x0_prior, yᵒ)
    llᵒ += path_log_likhd(OS(), ws.XXᵒ, ws.P, 1:m, ws.fpt)
    llᵒ += lobslikelihood(ws.P[1], yᵒ)

    print_info(verbose, it, value(ll), value(llᵒ), "impute")

    if accept_sample(llᵒ-ll, verbose)
        swap!(ws.XX, ws.XXᵒ, ws.WW, ws.WWᵒ, 1:m)
        set!(ws.z, zᵒ)
        return llᵒ, true
    else
        return ll, false
    end
end

"""
    proposal_start_pt(ws::Workspace, P, ρ)

Sample a new proposal starting point using a random walk with pCN step
...
# Arguments
- `ws`: Workspace, contains previously accepted `z` and the definition of a prior
- `P`: diffusion law according to which auxiliary posterior is to be computed
- `ρ`: memory parameter of the preconditioned Crank Nicolson scheme
...
"""
function proposal_start_pt(ws::Workspace, P, ρ)
    x0_prior, z = ws.x0_prior, ws.z.val
    zᵒ = rand(x0_prior, z, ρ)
    yᵒ = start_pt(zᵒ, x0_prior, P)
    zᵒ, yᵒ
end

"""
    sample_segments!(iRange, ws::Workspace{OS}, y, ρ,
                     headstart::Val{false}=Val{false}()
                     ) where OS

Sample paths segments in index range `iRange` using preconditioned
Crank-Nicolson scheme
...
# Arguments
- `iRange`: range of indices of the segments that need to be sampled
- `ws`: Workspace, containing all relevant containers of the Markov chain
- `y`: starting point
- `ρ`: memory parameter of the preconditioned Crank-Nicolson scheme
- `headstart`: whether to ease into first-passage time sampling
...
"""
function sample_segments!(iRange, ws::Workspace{OS}, y, ρ,
                          headstart::Val{false}=Val{false}()) where OS
    for i in iRange
        y = sample_segment!(i, ws, y, ρ)
    end
end

function sample_segments!(iRange, ws::Workspace{OS}, y, ρ,
                          headstart::Val{true}=Val{false}()) where OS
    for i in iRange
        sample_segment!(i, ws, y, ρ)
        while !checkFpt(OS(), ws.XXᵒ[i], ws.fpt[i])
            sample_segment!(i, ws, y, ρ)
        end
        y = ws.XXᵒ[i].yy[end]
    end
end

"""
    sample_segment!(i, ws, y, ρ)

Sample the `i`th path segment using preconditioned Crank-Nicolson scheme
...
# Arguments
- `i`: index of the segment to be sampled
- `ws`: Workspace, containing all relevant containers of the Markov chain
- `y`: starting point
- `ρ`: memory parameter of the preconditioned Crank-Nicolson scheme
...
"""
function sample_segment!(i, ws, y, ρ)
    sample!(ws.WWᵒ[i], ws.Wnr)
    crank_nicolson!(ws.WWᵒ[i].yy, ws.WW[i].yy, ρ)
    solve!(Euler(), ws.XXᵒ[i], y, ws.WWᵒ[i], ws.P[i]) # always according to trgt law
    ws.XXᵒ[i].yy[end]
end

"""
    crank_nicolson!(yᵒ, y, ρ)

Preconditioned Crank-Nicolson update with memory parameter `ρ`, previous vector
`y` and new vector `yᵒ`
"""
crank_nicolson!(yᵒ, y, ρ) = (yᵒ .= √(1-ρ)*yᵒ + √(ρ)*y)

"""
    path_log_likhd(::OS, XX, P, iRange, fpt; skipFPT=false
                   ) where OS <: AbstractObsScheme

Compute likelihood for path `XX` to be observed under `P`. Only segments with
index numbers in `iRange` are considered. `fpt` contains relevant info about
checks regarding adherence to first passage time pattern. `skipFPT` if set to
`true` can skip the step of checking adherence to fpt pattern (used for
conjugate updates, or any updates that keep `XX` unchanged)
"""
function path_log_likhd(::OS, XX, P, iRange, fpt; skipFPT=false
                        ) where OS <: AbstractObsScheme
    ll = 0.0
    for i in iRange
        ll += llikelihood(LeftRule(), XX[i], P[i])
    end
    !skipFPT && (ll = check_full_path_fpt(OS(), XX, iRange, fpt) ? ll : -Inf)
    !skipFPT && (ll += check_domain_adherence(P, XX, iRange) ? 0.0 : -Inf)
    ll
end

"""
    print_info(verbose::Bool, it::Integer, ll, llᵒ, msg="update")

Print information to the console about current likelihood values

...
# Arguments
- `verbose`: flag for whether to print anything at all
- `it`: iteration of the Markov chain
- `ll`: likelihood of the previous, accepted sample
- `llᵒ`: likelihood of the proposal sample
- `msg`: message to start with
...
"""
function print_info(verbose::Bool, it::Integer, ll, llᵒ, msg="update")
    verbose && print(msg, ": ", it, " ll ", round(ll, digits=3), " ",
                     round(llᵒ, digits=3), " diff_ll: ", round(llᵒ-ll,digits=3))
end

"""
    accept_sample(log_threshold, verbose=false)

Make a random MCMC decision for whether to accept a sample or reject it.
"""
function accept_sample(log_threshold, verbose=false)
    if rand(Exponential(1.0)) > -log_threshold # Reject if NaN
        verbose && print("\t ✓\n")
        return true
    else
        verbose && print("\t .\n")
        return false
    end
end

"""
    swap!(A, Aᵒ, iRange)

Swap contents between containers A & Aᵒ in the index range iRange
"""
function swap!(A, Aᵒ, iRange)
    for i in iRange
        A[i], Aᵒ[i] = Aᵒ[i], A[i]
    end
end

"""
    swap!(A, Aᵒ, B, Bᵒ, iRange)

Swap contents between containers A & Aᵒ in the index range iRange, do the same
for containers B & Bᵒ
"""
function swap!(A, Aᵒ, B, Bᵒ, iRange)
    swap!(A, Aᵒ, iRange)
    swap!(B, Bᵒ, iRange)
end

"""
    impute!(ws::Workspace{OS,ChequeredBlocking}, ll, verbose=false, it=NaN,
            headstart=false) where OS

Imputation step of the MCMC scheme (without blocking).
...
# Arguments
- `ws`: workspace containing most containers and definitions of the chain
- `ll`: log-likelihood of the old (previously accepted) diffusion path
- `verbose`: whether to print update info while sampling
- `it`: iteration index of the MCMC algorithm
- `headstart`: flag for whether to 'ease into' fpt conditions
...
"""
function impute!(ws::Workspace{OS,<:ChequeredBlocking}, ll, verbose=false,
                 it=NaN, headstart=false) where OS
    P, XXᵒ, XX = ws.P, ws.XXᵒ, ws.XX

    recompute_accepted_law!(ws)
    # proposal starting point:
    z_prop, y_prop = proposal_start_pt(ws, P[1], ws.ρ[ws.blidx][1])

    ll_total = 0.0
    for (block_idx, block) in enumerate(ws.blocking.blocks[ws.blidx])
        block_flag = Val{block[1]}()
        # previously accepted starting point
        y = XX[block[1]].yy[1]
        # proposal starting point for the block (can be non-y only for the first block)
        yᵒ = choose_start_pt(block_flag, y, y_prop)

        # sample path in block
        sample_segments!(block, ws, yᵒ, ws.ρ[ws.blidx][block_idx])
        set_end_pt_manually!(block_idx, block, ws)

        # starting point, path and observations contribution
        llᵒ = start_pt_log_pdf(block_flag, ws.x0_prior, yᵒ)
        llᵒ += path_log_likhd(OS(), XXᵒ, P, block, ws.fpt)
        llᵒ += lobslikelihood(P[block[1]], yᵒ)

        llPrev = start_pt_log_pdf(block_flag, ws.x0_prior, y)
        llPrev += path_log_likhd(OS(), XX, P, block, ws.fpt; skipFPT=true)
        llPrev += lobslikelihood(P[block[1]], y)

        print_info(verbose, it, value(llPrev), value(llᵒ), "impute")
        if accept_sample(llᵒ-llPrev, verbose)
            swap!(XX, XXᵒ, block)
            register_accpt!(ws, block_idx, true)
            set_z!(block_flag, ws, z_prop)
            ll_total += llᵒ
        else
            register_accpt!(ws, block_idx, false)
            ll_total += llPrev
        end
    end
    # acceptance indicator does not matter for sampling with blocking
    return ll_total, true
end

"""
    recompute_accepted_law!(ws::Workspace)

Recompute the (H, Hν, c) triplet as well as the noise that corresponds to law
`ws.P` that is obtained after switching blocks
"""
function recompute_accepted_law!(ws::Workspace)
    solve_back_rec!(ws, ws.P)         # compute (H, Hν, c) for given blocks
    noise_from_path!(ws, ws.XX, ws.WW, ws.P) # find noise WW that generates XX under 𝔅.P

    # compute white noise generating starting point under 𝔅
    z = inv_start_pt(ws.XX[1].yy[1], ws.x0_prior, ws.P[1])
    set!(ws.z, z)
end

"""
    solve_back_rec!(::NoBlocking, ws::Workspace{OS,B,ST}, P) where {OS,B,ST}

Solve backward recursion to find H, Hν and c, which together define r̃(t,x)
and p̃(x, 𝓓) under the auxiliary law, when no blocking is done
"""
function solve_back_rec!(::NoBlocking, ws::Workspace{OS,B,ST}, P) where {OS,B,ST}
    m = length(P)
    gpupdate!(P[m]; solver=ST())
    for i in (m-1):-1:1
        gpupdate!(P[i], P[i+1].H[1], P[i+1].Hν[1], P[i+1].c[1]; solver=ST())
    end
end

"""
    solve_back_rec!(ws::Workspace{OS,B,ST}, P) where {OS,B,ST}

Solve backward recursion to find H, Hν and c, which together define r̃(t,x)
and p̃(x, 𝓓) under the auxiliary law, when blocking is done
"""
function solve_back_rec!(ws::Workspace{OS,B,ST}, P) where {OS,B,ST}
    for block in reverse(ws.blocking.blocks[ws.blidx])
        gpupdate!(P[block[end]]; solver=ST())
        for i in reverse(block[1:end-1])
            gpupdate!(P[i], P[i+1].H[1], P[i+1].Hν[1], P[i+1].c[1]; solver=ST())
        end
    end
end

"""
    noise_from_path!(ws::Workspace, XX, WW, P)

Compute driving Wiener noise `WW` from path `XX` drawn under law `P`
"""
function noise_from_path!(ws::Workspace, XX, WW, P)
    𝔅 = ws.blocking
    for block in 𝔅.blocks[ws.blidx]
        for i in block
            inv_solve!(Euler(), XX[i], WW[i], P[i])
        end
    end
end

choose_start_pt(::Val{1}, y, yᵒ) = copy(yᵒ)
choose_start_pt(::Any, y, yᵒ) = copy(y)

"""
    set_end_pt_manually!(block_idx, block, ws::Workspace)

Manually set the end-point of the proposal path under blocking so that it agrees
with the end-point of the previously accepted path. If it is the last block,
then do nothing
"""
function set_end_pt_manually!(block_idx, block, ws::Workspace)
    if block_idx < length(ws.blocking.blocks[ws.blidx])
        ws.XXᵒ[block[end]].yy[end] = ws.XX[block[end]].yy[end]
    end
end

"""
    start_pt_log_pdf(::Val{1}, yPr::StartingPtPrior, y)

Compute the log-likelihood contribution of the starting point for a given prior
under a blocking scheme (intended to be used with a first block only)
"""
start_pt_log_pdf(::Val{1}, x0_prior::StartingPtPrior, y) = logpdf(x0_prior, y)

"""
    start_pt_log_pdf(::Any, yPr::StartingPtPrior, y)

Default contribution to log-likelihood from the startin point under blocking
"""
start_pt_log_pdf(::Any, x0_prior::StartingPtPrior, y) = 0.0

set_z!(::Val{1}, ws::Workspace, zᵒ) = set!(ws.z, zᵒ)

set_z!(::Any, ::Workspace, ::Any) = nothing

#===============================================================================
                            Parameter update routines
===============================================================================#

"""
    update_param!(pu::ParamUpdtDefn{MetropolisHastingsUpdt,UpdtIdx}, θ,
                  ws::Workspace{OS,NoBlocking}, ll, verbose=false, it=NaN
                  ) where {UpdtIdx,OS}
Update parameters
...
# Arguments
- `pu`: determines transition kernel, priors and which parameters to update
- `θ`: current value of the parameter
- `ws`: workspace with all necessary containers needed by the Markov chain
- `ll`: likelihood of the old (previously accepted) parametrisation
- `verbose`: whether to print updates info while sampling
- `it`: iteration index of the MCMC algorithm
...
"""
function update_param!(pu::ParamUpdtDefn{MetropolisHastingsUpdt,UpdtIdx}, θ,
                       ws::Workspace{OS,NoBlocking}, ll, verbose=false,
                       it=NaN) where {UpdtIdx,OS}
    WW, Pᵒ, P, XXᵒ, XX, fpt = ws.WW, ws.Pᵒ, ws.P, ws.XXᵒ, ws.XX, ws.fpt
    m = length(WW)
    θᵒ = rand(pu.t_kernel, θ, UpdtIdx())               # sample new parameter
    update_laws!(Pᵒ, θᵒ)
    pu.recompute_ODEs && solve_back_rec!(NoBlocking(), ws, Pᵒ) # compute (H, Hν, c)

    # find white noise which for a given θᵒ gives a correct starting point
    y = XX[1].yy[1]
    zᵒ = inv_start_pt(y, ws.x0_prior, Pᵒ[1])

    find_path_from_wiener!(XXᵒ, y, WW, Pᵒ, 1:m)

    llᵒ = logpdf(ws.x0_prior, y)
    llᵒ += path_log_likhd(OS(), XXᵒ, Pᵒ, 1:m, fpt)
    llᵒ += lobslikelihood(Pᵒ[1], y)

    print_info(verbose, it, ll, llᵒ)

    llr = ( llᵒ - ll + prior_kernel_contrib(pu.t_kernel, pu.priors, θ, θᵒ))

    # Accept / reject
    if accept_sample(llr, verbose)
        swap!(XX, XXᵒ, P, Pᵒ, 1:m)
        set!(ws.z, zᵒ)
        return llᵒ, true, θᵒ
    else
        return ll, false, θ
    end
end

"""
    update_laws!(Ps, θᵒ)

Set new parameter `θᵒ` for the laws in vector `Ps`
"""
function update_laws!(Ps, θᵒ)
    m = length(Ps)
    for i in 1:m
        Ps[i] = GuidPropBridge(Ps[i], θᵒ)
    end
end

"""
    find_path_from_wiener!(XX, y, WW, P, iRange)

Find path `XX` (that starts from `y`) that is generated under law `P` from the
Wiener process `WW`. Only segments with indices in range `iRange` are considered
"""
function find_path_from_wiener!(XX, y, WW, P, iRange)
    for i in iRange
        solve!(Euler(), XX[i], y, WW[i], P[i])
        y = XX[i].yy[end]
    end
end

"""
    prior_kernel_contrib(t_kern, priors, θ, θᵒ)

Contribution to the log-likelihood ratio from transition kernel `t_kern` and
`priors`.
"""
function prior_kernel_contrib(t_kern, priors, θ, θᵒ)
    llr = logpdf(t_kern, θᵒ, θ) - logpdf(t_kern, θ, θᵒ)
    for prior in priors
        llr += logpdf(prior, θᵒ) - logpdf(prior, θ)
    end
    llr
end

#NOTE blocking and no-blocking param update should be joined into one function
"""
    update_param!(pu::ParamUpdtDefn{MetropolisHastingsUpdt,UpdtIdx}, θ,
                  ws::Workspace{OS,B,ST}, ll, verbose=false, it=NaN
                  ) where {UpdtIdx,OS,B,ST}
Update parameters
"""
function update_param!(pu::ParamUpdtDefn{MetropolisHastingsUpdt,UpdtIdx},
                       θ, ws::Workspace{OS,B}, ll, verbose=false, it=NaN
                       ) where {UpdtIdx,OS,B}
    WW, Pᵒ, P, XXᵒ, XX, fpt = ws.WW, ws.Pᵒ, ws.P, ws.XXᵒ, ws.XX, ws.fpt
    m = length(P)
    θᵒ = rand(pu.t_kernel, θ, UpdtIdx())               # sample new parameter
    update_laws!(Pᵒ, θᵒ)                   # update law `Pᵒ` accordingly
    solve_back_rec!(ws, Pᵒ)                 # compute (H, Hν, c)

    y = XX[1].yy[1]
    zᵒ = inv_start_pt(y, ws.x0_prior, Pᵒ[1])

    llᵒ = logpdf(ws.x0_prior, y)
    for (block_idx, block) in enumerate(ws.blocking.blocks[ws.blidx])
        y = XX[block[1]].yy[1]
        find_path_from_wiener!(XXᵒ, y, WW, Pᵒ, block)
        set_end_pt_manually!(block_idx, block, ws)

        # Compute log-likelihood ratio
        llᵒ += path_log_likhd(OS(), XXᵒ, Pᵒ, block, ws.fpt)
        llᵒ += lobslikelihood(Pᵒ[block[1]], y)
    end
    print_info(verbose, it, ll, llᵒ)

    llr = ( llᵒ - ll + prior_kernel_contrib(pu.t_kernel, pu.priors, θ, θᵒ))

    # Accept / reject
    if accept_sample(llr, verbose)
        swap!(XX, XXᵒ, P, Pᵒ, 1:m)
        set!(ws.z, zᵒ)
        return llᵒ, true, θᵒ
    else
        return ll, false, θ
    end
end

"""
    update_param!(pu::ParamUpdtDefn{ConjugateUpdt,UpdtIdx}, θ,
                  ws::Workspace{OS,NoBlocking,ST}, ll, verbose=false, it=NaN
                  ) where {UpdtIdx,OS,ST}
Update parameters using conjugate draws
"""
function update_param!(pu::ParamUpdtDefn{ConjugateUpdt,UpdtIdx}, θ,
                       ws::Workspace{OS,NoBlocking}, ll, verbose=false, it=NaN
                       ) where {UpdtIdx,OS}
    WW, Pᵒ, P, XXᵒ, XX, fpt = ws.WW, ws.Pᵒ, ws.P, ws.XXᵒ, ws.XX, ws.fpt
    m = length(WW)
    θᵒ = conjugate_draw(θ, XX, P[1].Target, pu.priors[1], UpdtIdx())   # sample new parameter

    update_laws!(P, θᵒ)
    pu.recompute_ODEs && solve_back_rec!(NoBlocking(), ws, P) # compute (H, Hν, c)

    for i in 1:m    # compute wiener path WW that generates XX
        inv_solve!(Euler(), XX[i], WW[i], P[i])
    end
    # compute white noise that generates starting point
    y = XX[1].yy[1]
    z = inv_start_pt(y, ws.x0_prior, P[1])

    llᵒ = logpdf(ws.x0_prior, y)
    llᵒ += path_log_likhd(OS(), XX, P, 1:m, fpt; skipFPT=true)
    llᵒ += lobslikelihood(P[1], y)
    print_info(verbose, it, value(ll), value(llᵒ))
    set!(ws.z, z)
    return llᵒ, true, θᵒ
end

#NOTE blocking and no-blocking param conjugate update should be joined into one function
"""
    update_param!(pu::ParamUpdtDefn{ConjugateUpdt,UpdtIdx}, θ,
                  ws::Workspace{OS,B}, ll, verbose=false, it=NaN
                  ) where {UpdtIdx,OS,B}
Update parameters
see the definition of  update_param!(…, ::MetropolisHastingsUpdt, …) for the
explanation of the arguments.
"""
function update_param!(pu::ParamUpdtDefn{ConjugateUpdt,UpdtIdx}, θ,
                       ws::Workspace{OS,B}, ll, verbose=false, it=NaN
                       ) where {UpdtIdx,OS,B}
    WW, Pᵒ, P, XXᵒ, XX, fpt = ws.WW, ws.Pᵒ, ws.P, ws.XXᵒ, ws.XX, ws.fpt
    m = length(WW)
    θᵒ = conjugate_draw(θ, XX, P[1].Target, pu.priors[1], UpdtIdx())   # sample new parameter

    update_laws!(P, θᵒ)
    pu.recompute_ODEs && solve_back_rec!(ws, P)
    for i in 1:m    # compute wiener path WW that generates XX
        inv_solve!(Euler(), XX[i], WW[i], P[i])
    end
    # compute white noise that generates starting point
    y = XX[1].yy[1]
    z = inv_start_pt(y, ws.x0_prior, P[1])

    llᵒ = logpdf(ws.x0_prior, y)
    for block in ws.blocking.blocks[ws.blidx]
        llᵒ += path_log_likhd(OS(), XX, P, block, ws.fpt; skipFPT=true)
        llᵒ += lobslikelihood(P[block[1]], XX[block[1]].yy[1])
    end
    print_info(verbose, it, value(ll), value(llᵒ))
    set!(ws.z, z)
    return llᵒ, true, θᵒ
end


text:
"""
    parse_tuple(ps, ret)

`ret` is followed by a comma so tries to parse the rest of the
tuple.
"""
function parse_tuple(ps::ParseState, ret)
    startbyte = ps.nt.startbyte

    # Parsing
    next(ps)
    op = INSTANCE(ps)
    format_comma(ps)

    if isassignment(ps.nt) && ps.nt.kind != Tokens.APPROX
        if ret isa EXPR{TupleH}
            push!(ret.args, op)
            ret.span += op.span
        else
            ret =  EXPR{TupleH}(EXPR[ret, op], ret.span + op.span, Variable[], "")
        end
    elseif closer(ps)
        if ret isa EXPR{TupleH}  #(length(ret.punctuation) == 0 || !(first(ret.args) isa PUNCTUATION{Tokens.LPAREN}))
            push!(ret.args, op)
            ret.span += op.span
        else
            ret = EXPR{TupleH}(EXPR[ret, op], ret.span + op.span, Variable[], "")
        end
    else
        @catcherror ps startbyte nextarg = @closer ps tuple parse_expression(ps)
        if ret isa EXPR{TupleH} && (!(first(ret.args) isa EXPR{PUNCTUATION{Tokens.LPAREN}})) # && length(ret.punctuation) == 0 ||
            push!(ret.args, op)
            push!(ret.args, nextarg)
            ret.span += ps.nt.startbyte - startbyte
        else
            ret = EXPR{TupleH}(EXPR[ret, op, nextarg], ret.span + ps.nt.startbyte - startbyte, Variable[], "")
        end
    end
    return ret
end

ext: jl


text:
#
## Vanilla proximal gradient
#
abstract type ProxGradExtrapolation end
abstract type ProxGradExtrapolationState <: OptimizerState end

struct VanillaProxGrad <: ProxGradExtrapolation end
mutable struct VanillaProxGradState <: ProxGradExtrapolationState end

extrapolation_state(::VanillaProxGrad, x, g) = VanillaProxGradState()
function extrapolation!(::VanillaProxGrad, pgstate, extrapolationstate, pb)
    pgstate.temp .= pgstate.x .- pgstate.γ .* pgstate.∇f_x
    M = prox_αg!(pb, pgstate.x, pgstate.temp, pgstate.γ)
    pgstate.ncalls_proxg += 1
    pgstate.M = M
    return
end

Base.summary(::VanillaProxGrad) = ""

#
## Accelerated proximal gradient
#
## Inertial parameters from FISTA, but with possibilities suggested by Jingwei Liang...
@with_kw struct AcceleratedProxGrad <: ProxGradExtrapolation
    p::Float64 = 1.0
    q::Float64 = 1.0
    r::Float64 = 4.0
end
mutable struct AcceleratedProxGradState{Tx} <: ProxGradExtrapolation
    t::Float64
    y::Tx
    y_old::Tx
end

Base.summary(::AcceleratedProxGrad) = " Accelerated"

function extrapolation_state(::AcceleratedProxGrad, x, g)
    return AcceleratedProxGradState(1.0, zero(x), zero(x))
end
function extrapolation!(apg::AcceleratedProxGrad, pgstate, extrastate, pb)
    @unpack p, q, r = apg

    # y = proxgrad(x)
    pgstate.temp .= pgstate.x .- pgstate.γ .* pgstate.∇f_x
    M = prox_αg!(pb, extrastate.y, pgstate.temp, pgstate.γ)
    pgstate.M = M
    pgstate.ncalls_proxg += 1

    t_next = (p + sqrt(q + r * extrastate.t^2)) / 2

    # x = y + α (y-y_old)
    pgstate.x .=
        extrastate.y .+ (extrastate.t - 1) / t_next * (extrastate.y .- extrastate.y_old)
    extrastate.t = t_next
    extrastate.y_old .= extrastate.y

    return
end

#
## T1
#
@with_kw struct Test1ProxGrad <: ProxGradExtrapolation
    acceleration::AcceleratedProxGrad = AcceleratedProxGrad()
end
mutable struct Test1ProxGradState <: ProxGradExtrapolationState
    M_previousit::Manifold
    acceleration_state::AcceleratedProxGradState
end

function extrapolation_state(o::Test1ProxGrad, x, g::Regularizer)
    return Test1ProxGradState(
        wholespace_manifold(g, x),
        extrapolation_state(o.acceleration, x, g),
    )
end
function extrapolation!(o::Test1ProxGrad, pgstate, extrastate, pb)
    @unpack p, q, r = o.acceleration
    accel_state = extrastate.acceleration_state

    ## Compute proximal gradient iterate
    pgstate.temp .= pgstate.x .- pgstate.γ .* pgstate.∇f_x
    M = prox_αg!(pb, accel_state.y, pgstate.temp, pgstate.γ)
    pgstate.ncalls_proxg += 1

    t_next = (p + sqrt(q + r * accel_state.t^2)) / 2

    ## Compute the possibly accelerated iterate
    if !(M < extrastate.M_previousit)
        pgstate.x .=
            accel_state.y .+
            (accel_state.t - 1) / t_next * (accel_state.y .- accel_state.y_old)
    else
        pgstate.x .= accel_state.y
    end

    pgstate.M = M
    extrastate.M_previousit = M

    accel_state.t = t_next
    accel_state.y_old .= accel_state.y

    return
end

Base.summary(::Test1ProxGrad) = " - Test 1"

#
## T2
#
@with_kw struct Test2ProxGrad <: ProxGradExtrapolation
    acceleration::AcceleratedProxGrad = AcceleratedProxGrad()
end
mutable struct Test2ProxGradState{Tx} <: ProxGradExtrapolationState
    M_previousit::Manifold
    acceleration_state::AcceleratedProxGradState{Tx}
    candidate_pg::Tx
    candidate_apg::Tx
end

function extrapolation_state(o::Test2ProxGrad, x, g::Regularizer)
    return Test2ProxGradState(
        wholespace_manifold(g, x),
        extrapolation_state(o.acceleration, x, g),
        zero(x),
        zero(x),
    )
end
function extrapolation!(o::Test2ProxGrad, pgstate, extrastate, pb)
    @unpack p, q, r = o.acceleration
    accel_state = extrastate.acceleration_state

    ## Compute proximal gradient iterate
    pgstate.temp .= pgstate.x .- pgstate.γ .* pgstate.∇f_x
    M_accelerated = prox_αg!(pb, accel_state.y, pgstate.temp, pgstate.γ)
    pgstate.ncalls_proxg += 1

    t_next = (p + sqrt(q + r * accel_state.t^2)) / 2

    extrastate.candidate_pg .= accel_state.y
    extrastate.candidate_apg .=
        accel_state.y .+ (accel_state.t - 1) / t_next * (accel_state.y .- accel_state.y_old)

    # Compute image of candidates by PG, only manifold is relevant.
    ## 1. proxgrad
    pgstate.temp .= extrastate.candidate_pg .- pgstate.γ .* ∇f(pb, extrastate.candidate_pg)
    M_PG_pg = prox_αg!(pb, pgstate.temp, pgstate.temp, pgstate.γ)
    pgstate.ncalls_∇f += 1; pgstate.ncalls_proxg += 1

    pgstate.temp .=
        extrastate.candidate_apg .- pgstate.γ .* ∇f(pb, extrastate.candidate_apg)
    M_PG_apg = prox_αg!(pb, pgstate.temp, pgstate.temp, pgstate.γ)
    pgstate.ncalls_∇f += 1; pgstate.ncalls_proxg += 1

    ## Here, if apg has less structure than pg (in image by proxgrad), no acceleration
    if !(M_PG_pg < M_PG_apg)
        pgstate.x .= extrastate.candidate_apg
        pgstate.M = M_PG_pg
    else
        pgstate.x .= extrastate.candidate_pg
        pgstate.M = M_PG_apg
    end

    extrastate.M_previousit = pgstate.M

    accel_state.t = t_next
    accel_state.y_old .= accel_state.y

    return
end

Base.summary(::Test2ProxGrad) = " - Test 2"

#
## MFISTA
#
@with_kw struct MFISTA <: ProxGradExtrapolation
    acceleration::AcceleratedProxGrad = AcceleratedProxGrad()
end
mutable struct MFISTAState{Tx} <: ProxGradExtrapolationState
    acceleration_state::AcceleratedProxGradState{Tx}
    z::Tx
    z_old::Tx
end

function extrapolation_state(o::MFISTA, x, g::Regularizer)
    return MFISTAState(extrapolation_state(o.acceleration, x, g), copy(x), copy(x))
end

"""
    extrapolation!(o::MFISTA, pgstate, extrastate, pb)

MFISTA algorithm.

NOTE: implementaiton might be improved.
"""
function extrapolation!(o::MFISTA, pgstate, extrastate, pb)
    @unpack p, q, r = o.acceleration
    accel_state = extrastate.acceleration_state

    ## Compute proximal gradient iterate
    pgstate.temp .= pgstate.x .- pgstate.γ .* pgstate.∇f_x
    M = prox_αg!(pb, accel_state.y, pgstate.temp, pgstate.γ)
    pgstate.ncalls_proxg += 1

    t_next = (p + sqrt(q + r * accel_state.t^2)) / 2

    if F(pb, accel_state.y) <= F(pb, extrastate.z_old)
        extrastate.z .= accel_state.y
    else
        extrastate.z .= extrastate.z_old
    end

    pgstate.x .=
        extrastate.z .+ accel_state.t / t_next .* (accel_state.y .- extrastate.z) .+
        (accel_state.t - 1) / t_next .* (extrastate.z .- extrastate.z_old)
    pgstate.M = M

    extrastate.z_old .= extrastate.z
    accel_state.t = t_next

    return
end

Base.summary(::MFISTA) = "MFISTA"

#
## Restarted Accelerated proximal gradient
#
## Inertial parameters from FISTA, but with possibilities suggested by Jingwei Liang...
@with_kw struct RestartedAPG <: ProxGradExtrapolation
    acceleration::AcceleratedProxGrad = AcceleratedProxGrad()
    period::Int64 = 15
end
mutable struct RestartedAPGState{Tx} <: ProxGradExtrapolation
    acceleration_state::AcceleratedProxGradState{Tx}
    nbAPGiterations::Int64
end

Base.summary(::RestartedAPG) = " Restarted A."

function extrapolation_state(o::RestartedAPG, x, g::Regularizer)
    return RestartedAPGState(extrapolation_state(o.acceleration, x, g), 0)
end

function extrapolation!(rapg::RestartedAPG, pgstate, extrastate, pb)
    @unpack p, q, r = rapg.acceleration
    accel_state = extrastate.acceleration_state

    ## Compute proximal gradient iterate
    pgstate.temp .= pgstate.x .- pgstate.γ .* pgstate.∇f_x
    M = prox_αg!(pb, accel_state.y, pgstate.temp, pgstate.γ)
    pgstate.ncalls_proxg += 1

    t_next = (p + sqrt(q + r * accel_state.t^2)) / 2

    pgstate.x .=
        accel_state.y .+
        (accel_state.t - 1) / t_next * (accel_state.y .- accel_state.y_old)

    pgstate.M = M

    accel_state.t = t_next
    accel_state.y_old .= accel_state.y

    extrastate.nbAPGiterations += 1
    if extrastate.nbAPGiterations == rapg.period
        extrastate.nbAPGiterations = 0
        extrastate.acceleration_state.t = 1
        extrastate.acceleration_state.y .= pgstate.x
        extrastate.acceleration_state.y_old .= pgstate.x
    end
    return
end


text:
struct WeberSoundHooks <: TimedSound.Hooks
end

function TimedSound.show_latency_warnings(::WeberSoundHooks)
  show_sound_latency_warnings()
end

function show_sound_latency_warnings()
  if in_experiment()
    !info(get_experiment()).warn_on_trials_only || Weber.trial() > 0
  else
    true
  end
end

function TimedSound.on_play(::WeberSoundHooks)
  if in_experiment() && !experiment_running()
    error("You cannot call `play` during experiment `setup`. During `setup`",
          " you should add play to a trial (e.g. ",
          "`addtrial(moment(play,my_sound))`).")
  end
end

TimedSound.tick(::WeberSoundHooks) = Weber.tick()
function TimedSound.on_high_latency(::WeberSoundHooks,latency)
  record("high_latency",value=latency)
end

function TimedSound.on_no_timing(::WeberSoundHooks)
  warn("Cannot guarantee the timing of a sound. Add a delay before playing the",
       " sound if precise timing is required.",moment_trace_string())
end

function TimedSound.get_streamers(::WeberSoundHooks)
  if in_experiment()
    data(get_experiment()).streamers
  else
    error("Cannot stream sounds during setup. Stream a sound within a moment.",
          "(e.g. moment(play,tone(1kHz))")
  end
end


text:
using RecipesBase

@recipe function f(ga::GeoArray; band=1)
    xflip --> false
    yflip --> false
    aspect_ratio --> 1
    seriestype := :heatmap
    color := :viridis

    is_rotated(ga) && (ga = straighten(ga))

    c = GeoArrays.coords(ga, Vertex())
    x = map(x -> x[1], c[:, 1])
    y = map(x -> x[2], c[end, :])
    z = ga.A[:, :, band]'

    # Can't use x/yflip as x/y coords
    # have to be sorted for Plots
    if ga.f.linear[1] < 0
        z = reverse(z, dims=2)
        reverse!(x)
    end
    if ga.f.linear[4] < 0
        z = reverse(z, dims=1)
        reverse!(y)
    end

    xlims --> (extrema(x))
    ylims --> (extrema(y))

    x, y, z
end


text:
module TestSelect

using Julio
using Test

include("../../../examples/select.jl")

function test_select_take()
    output = nothing
    selected = nothing
    ie1, _ = Julio.channel(Int)
    ie2, oe2 = Julio.stack(Int)
    begin
        begin
            put!(ie2, 222)
            Julio.select(
                (put!, ie1, 111) => _ -> begin
                    selected = :put!
                end,
                (take!, oe2) => x -> begin
                    output = x
                    selected = :take!
                end,
            )
        end
    end
    @test selected === :take!
    @test output == 222
end

function test_select_cmd()
    Julio.open(`echo "hello"`) do io
        _, o1 = Julio.channel()
        begin
            selected = Julio.select(
                Events.readline(io; keep = true),  # should win
                (take!, o1),
            )
            @test selected == "hello\n"
        end
    end
end

end  # module


text:
include("../Preprocessing/Encoder.jl")

"""
    SVC(; alpha=0.01, ni=1000)
Support Vector Machine Classifier.

This struct learns classifiers using One-Vs-Rest.
One-Vs-Rest generates two-class classifiers divided into one class and the other classes using Logistic Regression, adopting the most likely one among all classifiers.

Parameter `α` indicates the learning rate, and `ni` indicates the number of learnings.

# Example
```jldoctest classification
julia> model = SVC()
SVC(0.01, 1000, Logistic[])

julia> fit!(model, ct)
3-element Vector{Logistic}:
 Logistic(0.01, 1000, [0.8116709490679518 1.188329050932049; 1.7228257190036231 0.2771742809963788; -0.1519960725403138 2.1519960725403116])
 Logistic(0.01, 1000, [0.9863693439936144 1.0136306560063886; 0.8838433946106077 1.11615660538939; 1.4431044559203794 0.5568955440796174])
 Logistic(0.01, 1000, [1.262510641510418 0.7374893584895849; 0.5242383002319192 1.4757616997680822; 1.864635796779504 0.135364203220495])

julia> println(predict(model, x))
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3]
```
"""
mutable struct SVC
    α::Float64
    n_iter::Int64
    classifiers::Array{Logistic, 1}
    SVC(; alpha=0.01, ni=1000) = new(alpha, ni, Array{Logistic}(undef, 0))
end

function fit!(model::SVC, x, t)
    alpha, n_iter = model.α, model.n_iter
    check_size(x, t)
    c = size(t, 2)
    w = Array{Float64}(undef, 0, 0)
    classifiers = Array{Logistic}(undef, 0)
    for i in 1 : c
        classifier = Logistic(alpha = alpha, ni = n_iter)
        OHE = OneHotEncoder()
        fit!(classifier, x, OHE(t[:, i]))
        push!(classifiers, classifier)
    end
    model.classifiers = classifiers
end

function predict(model::SVC, x)
    p = Array{Float64}(undef, 0, 0)
    for i in 1 : length(model.classifiers)
        if i == 1
            p = forecast(model.classifiers[i], x)[:, 2]
        else
            p = hcat(p, forecast(model.classifiers[i], x)[:, 2])
        end
    end
    return [findfirst(p[i, :] .== maximum(p[i, :])) for i in 1:size(p, 1)]
end


text:
"""
    Base.copy(m::StandardModel)

Shallow copy of a [`StandardModel`](@ref)
"""
Base.copy(m::StandardModel) = StandardModel(m.id, m.reactions, m.metabolites, m.genes)

"""
    Base.copy(r::Reaction)

Shallow copy of a [`Reaction`](@ref)
"""
Base.copy(r::Reaction) = Reaction(
    r.id;
    name = r.name,
    metabolites = r.metabolites,
    lb = r.lb,
    ub = r.ub,
    grr = r.grr,
    subsystem = r.subsystem,
    notes = r.notes,
    annotations = r.annotations,
    objective_coefficient = r.objective_coefficient,
)

"""
    Base.copy(m::Metabolite)

Shallow copy of a [`Metabolite`](@ref)
"""
Base.copy(m::Metabolite) = Metabolite(
    m.id;
    name = m.name,
    formula = m.formula,
    charge = m.charge,
    compartment = m.compartment,
    notes = m.notes,
    annotations = m.annotations,
)

"""
    Base.copy(g::Gene)

Shallow copy of a [`Gene`](@ref)
"""
Base.copy(g::Gene) = Gene(g.id; name = g.name, notes = g.notes, annotations = g.annotations)

"""
    atom_exchange(flux_dict::Dict{String, Float64}, model::StandardModel)

Return a dictionary mapping the flux of atoms across the boundary of the model
given `flux_dict` (the solution of a constraint based analysis) of reactions in `model`.
"""
function atom_exchange(flux_dict::Dict{String,Float64}, model::StandardModel)
    atom_flux = Dict{String,Float64}()
    for (rxn_id, flux) in flux_dict
        if is_boundary(model.reactions[rxn_id])
            for (met, stoich_rxn) in model.reactions[rxn_id].metabolites
                adict = get_atoms(model.metabolites[met])
                for (atom, stoich_molecule) in adict
                    atom_flux[atom] =
                        get(atom_flux, atom, 0.0) + flux * stoich_rxn * stoich_molecule
                end
            end
        end
    end
    return atom_flux
end

"""
    atom_exchange(rxn_id::String, model::StandardModel)

Return a dictionary mapping the flux of atoms through a reaction in `model`.
"""
function atom_exchange(rxn_id::String, model::StandardModel)
    atom_flux = Dict{String,Float64}()
    for (met, stoich_rxn) in model.reactions[rxn_id].metabolites
        adict = get_atoms(model.metabolites[met])
        for (atom, stoich_molecule) in adict
            atom_flux[atom] = get(atom_flux, atom, 0.0) + stoich_rxn * stoich_molecule
        end
    end
    return atom_flux
end

"""
    metabolite_fluxes(flux_dict::Dict{String, Float64}, model::StandardModel)

Return two dictionaries of metabolite `id`s mapped to reactions that consume or
produce them given the flux distribution supplied in `fluxdict`.
"""
function metabolite_fluxes(flux_dict::Dict{String,Float64}, model::StandardModel)
    S = stoichiometry(model)
    rxnids = reactions(model)
    metids = metabolites(model)

    producing = Dict{String,Dict{String,Float64}}()
    consuming = Dict{String,Dict{String,Float64}}()
    for (row, metid) in enumerate(metids)
        for (col, rxnid) in enumerate(rxnids)
            mf = flux_dict[rxnid] * S[row, col]
            # ignore zero flux
            if mf < -_constants.tolerance # consuming rxn
                if haskey(consuming, metid)
                    consuming[metid][rxnid] = mf
                else
                    consuming[metid] = Dict(rxnid => mf)
                end
            elseif mf > _constants.tolerance
                if haskey(producing, metid)
                    producing[metid][rxnid] = mf
                else
                    producing[metid] = Dict(rxnid => mf)
                end
            end
        end
    end
    return consuming, producing
end


text:
include("../coin-or-common.jl")

sources = [
    GitSource("https://github.com/coin-or/Bonmin.git", Bonmin_gitsha),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/Bonmin

if [[ ${target} == *mingw* ]]; then
    sed -i s/dllimport/dllexport/ "${includedir}/coin-or/IpoptConfig.h"
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/fix_dllexport.patch
fi

# Remove misleading libtool files
rm -f ${libdir}/*.la
update_configure_scripts

# old and custom autoconf
sed -i s/elf64ppc/elf64lppc/ configure

export CPPFLAGS="${CPPFLAGS} -I${includedir} -I${includedir}/coin -I${includedir}/coin-or"
export CXXFLAGS="${CXXFLAGS} -std=c++11"

if [[ ${target} == *mingw* ]]; then
    export LDFLAGS="-L${libdir}"
    export LT_LDFLAGS="-no-undefined"
fi

./configure \
    --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --disable-pkg-config \
    --disable-debug \
    --enable-shared \
    lt_cv_deplibs_check_method=pass_all \
    --with-asl-lib="-lasl" \
    --with-coindepend-lib="-lCbc -lCgl -lOsiClp -lClp -lOsi -lCoinUtils -lipoptamplinterface -lipopt"

make -j${nproc}
make install

# Prevent IpoptConfig.h from being included in the tarball after editing it above.
rm "${includedir}/coin-or/IpoptConfig.h"
"""

filter!(x -> cxxstring_abi(x) != "cxx03", platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libbonmin", :libbonmin),
    LibraryProduct("libbonminampl", :libbonminampl),
    ExecutableProduct("bonmin", :amplexe),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("ASL_jll", ASL_version),
    Dependency("Cbc_jll", compat="$(Cbc_version)"),
    Dependency("Cgl_jll", compat="$(Cgl_version)"),
    Dependency("Clp_jll", compat="$(Clp_version)"),
    Dependency("Osi_jll", compat="$(Osi_version)"),
    Dependency("CoinUtils_jll", compat="$(CoinUtils_version)"),
    Dependency("Ipopt_jll", compat="$(Ipopt_version)"),
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(
    ARGS,
    "Bonmin",
    Bonmin_version,
    sources,
    script,
    platforms,
    products,
    dependencies,
    julia_compat = "1.6",
)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 01-2018
#--------------------------------------------------------------------

using ScalarWave
using Base.Test

libraries = ["SpecCalculus",
             "Utilities",
             "Physics",
             "Patch",
             "Grid",
             "Projection",
             "Dispatch",
             "Futures",
             "Visualization"]

for file in libraries
    info("Testing $file")
    include("$(file)Test.jl")
end


text:
#=
Test Resynthesizer by comparing out images to reference good result image.

Using testing framework.
=#

# See "Code Loading"
# Either include, or import package
include("../../src/resynthesizer.jl")
# using Resynthesizer

# the Julia testing framework
using ReferenceTests
using Test

# profiling framework
using Profile

# for "testimage" of packaged, remote  test files ?
using TestImages

# For "load" local  files
using FileIO

using ColorTypes

#using Pkg
#Pkg.instantiate()

# This is executable, not a definition
@testset "Resynthesizer" begin

  # donkey image

  # Don't understand why I can't use relative path here
  inPath  = "/home/bootch/git/resynthesizerJulia/test/data/in/lighthouse.png"
  # Reference image in this directory
  outPath = "data/out/lighthouse.png"

  image = load(inPath)

  mask = falses(size(image))
  # upright rectangular mask of the lighthouse itself
  mask[50:350,300:400] .= true

  out = resynthesize(image, mask)

  # Julia macro that compares file and expression
  @test_reference outPath out
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
# OPERATORS defines single characters that should be converted to an MML
# operator element: <mo></mo>.  There are other MML operators that can not be
# represented by a single ASCII character in TeX, e.g, \{, \}, \subset, \vee,
# \cup, \in, \forall, \approx, etc. They are processed in the command reader.
# chars in BRACKETS are converted to TK_PAREN type which will also be
# converted to <mo> but could carry size attributes.
const OPERATOR_IDS = "+-*/!:=,.\'"
const BRACKETS = "|<>()[]"

# types
# NOTE: Most are not types in TeX, e.g., objects in TK_OPID just represent
# regular symbols in TeX, but they will be converted to operator elements,
# i.e., <mo> elements in MathML
const TK_MATH = 0               # to the top elelment <math>
const TK_ID = 1                 # to <mi>, single-char identifiers
const TK_NUM = 2                # to <mn>
#const TK_SYM = 3               # to <mi> with UCS codepoint literal
const TK_FUN = 4                # to <mi>, multi-char identifiers
#const TK_CMD = 5
const TK_ENV = 6
const TK_OPID = 7               # convert to <mo>
# parenthesis' are also OPIDs. They are a little more special than single-char
# OPIDs in that they are paired, and may be rendered in a specified size when
# a \big TeX command is present (in this case, they are converted to <mo> with
# `maxsize` or `minsize` attributes)
const TK_PAREN = 8
const TK_UNKNOWN = 9
# types that need special layout treatment when sub-sup operators are present.
# TK_LIMs are parsed to TK_FUN, which are ultimately just special identifiers.
# TK_INTs and TK_BIGOPs are parsed to TK_OPID (they are math operator identifiers,
# not true TeX language operators)
const TK_LIM = 33               # \lim, \liminf, ... only subscripts
const TK_INT = 34               # \int, ..., sub-super-scripts
const TK_BIGOP = 35             # \sum, \prod, \bigcup, ..., sub-super-scripts

# TeX language commands, each with its own type
const TK_LBRACE = 11            # '{': initiate an expression group
const TK_RBRACE = 12            # '}': end an expression group
const TK_SUB = 13               # '-': subscript
const TK_SUP = 14               # '^': superscript
const TK_AMP = 15               # '&': ampersand
const TK_NEWLINE = 16           # '\\': new line
const TK_OPNAME = 70            # \operatorname
const TK_TEXT = 41
const TK_SQRT = 42
const TK_FRAC = 43
const TK_LEFT = 44
const TK_RIGHT = 45
const TK_MIDDLE = 46
const TK_BEGIN = 47
const TK_END = 48
const TK_BINOM = 49
const TK_OVERSET = 50
const TK_UNDERSET = 51
const TK_FRACOVER = 52
# TK_OEVRACCENT(TK_UNDERACCENT) are tokens for TeX command: \dot,
# \tilde, \hat, etc. The <mover> (<munder>) element carries an attrubute
# `accent` which is set to be "true", e.g.,
# "\hat{A}" is converted to
#   <mover accent="true">
#     <mi>A</mi>
#     <mo> ^<!--CIRCUMFLEX ACCENT--></mo>
#   </mover>
const TK_OVERACCENT = 54
const TK_UNDERACCENT = 55
const TK_SPACE = 60
const TK_BIG = 61
# \overbrace (\underbrace) are converted to <mover> (<munder>). They could
# carry a super(sub)-script.
const TK_OVERBRACE = 52
const TK_UNDERBRACE = 53

# other MML operators

# an expression group with specified style
# It seems MathJax does not support mathvariant, we may not implement
# style for now... need to do more experiments.
const TK_STYLE = 100            # \mathrm, \mathcal, ... convert to <mstyle>
const TK_IGSTYLE = 101          # in-group style: \bf, \it, \red, ...
const TK_ENV_NAME = 110         # env-names that follow \begin(\end) commands
const TK_ENV_ATTR = 111         # env-attributes that follow \begin(\end) commands

struct Token
    type::Int
    val::String
end
Token(t::Int, v::Char) = Token(t, string(v))

# symbols are copied from `latex_symbols.jl` in the `REPL` module.
const command_to_token = Dict(
    # styles
    "mathbf"     => Token(TK_STYLE, "mathbf"),
    "mathcal"    => Token(TK_STYLE, "mathcal"),
    "mathbb"     => Token(TK_STYLE, "mathbb"),
    "mathfrak"   => Token(TK_STYLE, "mathfrak"),
    "bigl"       => Token(TK_BIG, "1.2em"),
    "bigr"       => Token(TK_BIG, "1.2em"),
    "Bigl"       => Token(TK_BIG, "1.623em"),
    "Bigr"       => Token(TK_BIG, "1.623em"),
    "biggl"      => Token(TK_BIG, "2.047em"),
    "biggr"      => Token(TK_BIG, "2.047em"),
    "Biggl"      => Token(TK_BIG, "2.470em"),
    "Biggr"      => Token(TK_BIG, "2.470em"),
    "it"         => Token(TK_IGSTYLE, "it"),
    "bf"         => Token(TK_IGSTYLE, "bf"),
    "rm"         => Token(TK_IGSTYLE, "rf"),
    # commands
    "text"          => Token(TK_TEXT, ""),
    "sqrt"          => Token(TK_SQRT, ""),
    "frac"          => Token(TK_FRAC, ""),
    "over"          => Token(TK_FRACOVER, ""),
    "binom"         => Token(TK_BINOM, ""),
    "left"          => Token(TK_LEFT, ""),
    "right"         => Token(TK_RIGHT, ""),
#    "middle"        => Token(TK_MIDDLE, ""),
    "begin"         => Token(TK_BEGIN, ""),
    "end"           => Token(TK_END, ""),
    "overset"       => Token(TK_OVERSET, ""),
    "underset"      => Token(TK_UNDERSET, ""),
    # if \overbrace (\underspace) carries a super-script (sub-script)
    # operator, it has to be dealt with specially as the result is NOT
    # a usual super-scripted (sub-scripted) glyph.
    "overbrace"     => Token(TK_OVERBRACE, '\u23de'),
    "underbrace"    => Token(TK_UNDERBRACE, '\u23df'),
    # operator name
    "operatorname"  => Token(TK_OPNAME, ""),
    # over/under-script operator as accent
    "dot"               => Token(TK_OVERACCENT, '˙'),
    "ddot"              => Token(TK_OVERACCENT, '¨'),
    "bar"               => Token(TK_OVERACCENT, '¯'),
    "hat"               => Token(TK_OVERACCENT, '^'),
    "check"             => Token(TK_OVERACCENT, 'ˇ'),
    "breve"             => Token(TK_OVERACCENT, '˘'),
    "acute"             => Token(TK_OVERACCENT, '´'),
    "grave"             => Token(TK_OVERACCENT, '`'),
    "tilde"             => Token(TK_OVERACCENT, '~'),
    "vec"               => Token(TK_OVERACCENT, '→'),
    "overline"          => Token(TK_OVERACCENT, '_'),
    "underline"         => Token(TK_OVERACCENT, '_'),
    "widehat"           => Token(TK_OVERACCENT, '^'),
    "widetilde"         => Token(TK_OVERACCENT, '~'),
    "overrightarrow"    => Token(TK_OVERACCENT, '→'),
    "overleftarrow"     => Token(TK_OVERACCENT, '←'),
    # spaces and new-line
    # Space tokens are conveted to <mspace> elements
    # new-line is also converted to <mspace> with `linebreak="newline"`
    "\\"         => Token(TK_NEWLINE, ""),
    "!"          => Token(TK_SPACE, "-0.167em"),    # -3./18.
    ","          => Token(TK_SPACE, "0.167em"),     #  3./18.
    ":"          => Token(TK_SPACE, "0.222em"),     #  4./18.
    ";"          => Token(TK_SPACE, "0.278em"),     #  5./18"
    " "          => Token(TK_SPACE, "1.0em"),       #  1.
    "quad"       => Token(TK_SPACE, "1.0em"),       #  1.
    "qquad"      => Token(TK_SPACE, "2.0em"),       #  2.
    # parenthesis
    # parenthesis token are converted to either <mo> or <mo> with sizes if
    # \left (\right) command is present
    "langle"     => Token(TK_PAREN, "⟨"),
    "rangle"     => Token(TK_PAREN, "⟩"),
    "{"          => Token(TK_PAREN, "{"),
    "}"          => Token(TK_PAREN, "}"),
    "lceil"      => Token(TK_PAREN, "⌈"),
    "rceil"      => Token(TK_PAREN, "⌉"),
    "lfloor"     => Token(TK_PAREN, "⌊"),
    "rfloor"     => Token(TK_PAREN, "⌋"),
    "|"          => Token(TK_PAREN, '∥'),
    # \lim, \int, and \bigop, we might need some special treatment?
    "lim"        => Token(TK_LIM, "lim"),
    "liminf"     => Token(TK_LIM, "lim inf"),
    "limsup"     => Token(TK_LIM, "lim sup"),
    "min"        => Token(TK_LIM, "min"),
    "max"        => Token(TK_LIM, "max"),
    "inf"        => Token(TK_LIM, "inf"),
    "sup"        => Token(TK_LIM, "sup"),
    "int"        => Token(TK_INT, '∫'),
    "iint"       => Token(TK_INT, '∬'),
    "iiint"      => Token(TK_INT, '∭'),
    "oint"       => Token(TK_INT, '∮'),
    "sum"        => Token(TK_BIGOP, '∑'),
    "prod"       => Token(TK_BIGOP, '∏'),
    "coprod"     => Token(TK_BIGOP, '∐'),
    "bigcap"     => Token(TK_BIGOP, '⋂'),
    "bigcup"     => Token(TK_BIGOP, '⋃'),
    "bigsqcup"   => Token(TK_BIGOP, '⨆'),
    "bigvee"     => Token(TK_BIGOP, '⋁'),
    "bigwedge"   => Token(TK_BIGOP, '⋀'),
    "bigodot"    => Token(TK_BIGOP, '⨀'),
    "bitotimes"  => Token(TK_BIGOP, '⨂'),
    "bigoplus"   => Token(TK_BIGOP, '⨁'),
    "biguplus"   => Token(TK_BIGOP, '⨄'),

    # environments
    "matrix"     => Token(TK_ENV, "matrix"),
    "pmatrix"    => Token(TK_ENV, "pmatrix"),
    "bmatrix"    => Token(TK_ENV, "bmatrix"),
    "vmatrix"    => Token(TK_ENV, "vmatrix"),

    # symbols
    "Alpha"      => Token(TK_ID, 'Α'),
    "alpha"      => Token(TK_ID, 'α'),
    "Beta"       => Token(TK_ID, 'Β'),
    "beta"       => Token(TK_ID, 'β'),
    "Gamma"      => Token(TK_ID, 'Γ'),
    "gamma"      => Token(TK_ID, 'γ'),
    "digamma"    => Token(TK_ID, 'ϝ'),
    "Delta"      => Token(TK_ID, 'Δ'),
    "delta"      => Token(TK_ID, 'δ'),
    "Epsilon"    => Token(TK_ID, 'Ε'),
    "epsilon"    => Token(TK_ID, 'ϵ'),
    "varepsilon" => Token(TK_ID, 'ε'),
    "Zeta"       => Token(TK_ID, 'Ζ'),
    "zeta"       => Token(TK_ID, 'ζ'),
    "Eta"        => Token(TK_ID, 'Η'),
    "eta"        => Token(TK_ID, 'η'),
    "Theta"      => Token(TK_ID, 'Θ'),
    "theta"      => Token(TK_ID, 'θ'),
    "vartheta"   => Token(TK_ID, 'ϑ'),
    "Iota"       => Token(TK_ID, 'Ι'),
    "iota"       => Token(TK_ID, 'ι'),
    "Kappa"      => Token(TK_ID, 'Κ'),
    "kappa"      => Token(TK_ID, 'κ'),
    "Lambda"     => Token(TK_ID, 'Λ'),
    "lambda"     => Token(TK_ID, 'λ'),
    "Mu"         => Token(TK_ID, 'Μ'),
    "mu"         => Token(TK_ID, 'μ'),
    "Nu"         => Token(TK_ID, 'Ν'),
    "nu"         => Token(TK_ID, 'ν'),
    "Xi"         => Token(TK_ID, 'Ξ'),
    "xi"         => Token(TK_ID, 'ξ'),
    "Omicron"    => Token(TK_ID, 'Ο'),
    "omicron"    => Token(TK_ID, 'ο'),
    "Pi"         => Token(TK_ID, 'Π'),
    "pi"         => Token(TK_ID, 'π'),
    "varpi"      => Token(TK_ID, 'ϖ'),
    "Rho"        => Token(TK_ID, 'Ρ'),
    "rho"        => Token(TK_ID, 'ρ'),
    "varrho"     => Token(TK_ID, 'ϱ'),
    "Sigma"      => Token(TK_ID, 'Σ'),
    "sigma"      => Token(TK_ID, 'σ'),
    "varsigma"   => Token(TK_ID, 'ς'),
    "Tau"        => Token(TK_ID, 'Τ'),
    "tau"        => Token(TK_ID, 'τ'),
    "Upsilon"    => Token(TK_ID, 'Υ'),
    "upsilon"    => Token(TK_ID, 'υ'),
    "Phi"        => Token(TK_ID, 'Φ'),
    "phi"        => Token(TK_ID, 'ϕ'),
    "varphi"     => Token(TK_ID, 'φ'),
    "Chi"        => Token(TK_ID, 'Χ'),
    "chi"        => Token(TK_ID, 'χ'),
    "Psi"        => Token(TK_ID, 'Ψ'),
    "psi"        => Token(TK_ID, 'ψ'),
    "Omega"      => Token(TK_ID, 'Ω'),
    "omega"      => Token(TK_ID, 'ω'),
    "aleph"      => Token(TK_ID, 'ℵ'),
    "O"          => Token(TK_ID, 'Ø'),
    "o"          => Token(TK_ID, 'ø'),
    "ss"         => Token(TK_ID, 'ß'),
    "imath"      => Token(TK_ID, 'ı'),
    "jmath"      => Token(TK_ID, 'ȷ'),
    "ell"        => Token(TK_ID, 'ℓ'),
    "hbar"       => Token(TK_ID, 'ℏ'),
    "hslash"     => Token(TK_ID, 'ℏ'),
    "infty"      => Token(TK_ID, '∞'),
    "Re"         => Token(TK_ID, 'ℜ'),
    "Im"         => Token(TK_ID, 'ℑ'),
    "complement" => Token(TK_ID, '∁'),
    "emptyset"   => Token(TK_ID, '∅'),
    "therefore"  => Token(TK_ID, '∴'),
    "because"    => Token(TK_ID, '∵'),
    "Diamond"    => Token(TK_ID, '◊'),
    "Box"        => Token(TK_ID, '◻'),
    "triangle"   => Token(TK_ID, '△'),
    "angle"      => Token(TK_ID, '∠'),
    # function names
    "sin"        => Token(TK_FUN, "sin"),
    "cos"        => Token(TK_FUN, "cos"),
    "tan"        => Token(TK_FUN, "tan"),
    "csc"        => Token(TK_FUN, "csc"),
    "sec"        => Token(TK_FUN, "sec"),
    "cot"        => Token(TK_FUN, "cot"),
    "arcsin"     => Token(TK_FUN, "arcsin"),
    "arccos"     => Token(TK_FUN, "arccos"),
    "arctan"     => Token(TK_FUN, "arctan"),
    "sinh"       => Token(TK_FUN, "sinh"),
    "cosh"       => Token(TK_FUN, "cosh"),
    "tanh"       => Token(TK_FUN, "tanh"),
    "coth"       => Token(TK_FUN, "coth"),
    "exp"        => Token(TK_FUN, "exp"),
    "ln"         => Token(TK_FUN, "ln"),
    "log"        => Token(TK_FUN, "log"),
    "erf"        => Token(TK_FUN, "erf"),
    "erfc"       => Token(TK_FUN, "erfc"),
    "arg"        => Token(TK_FUN, "arg"),
    "ker"        => Token(TK_FUN, "ker"),
    "dim"        => Token(TK_FUN, "dim"),
    "det"        => Token(TK_FUN, "det"),
    # operators, mostly infix, some are prefix
    "times"              => Token(TK_OPID, '×'),
    "oplus"              => Token(TK_OPID, '⊕'),
    "ominus"             => Token(TK_OPID, '⊖'),
    "otimes"             => Token(TK_OPID, '⊗'),
    "oslash"             => Token(TK_OPID, '⊘'),
    "odot"               => Token(TK_OPID, '⊙'),
    "bigcirc"            => Token(TK_OPID, '◯'),
    "amalg"              => Token(TK_OPID, '⨿'),
    "pm"                 => Token(TK_OPID, '±'),
    "mp"                 => Token(TK_OPID, '∓'),
    "cdot"               => Token(TK_OPID, '·'),
    "cdots"              => Token(TK_OPID, '⋯'),
    "vdots"              => Token(TK_OPID, '⋮'),
    "ldots"              => Token(TK_OPID, '…'),
    "ddots"              => Token(TK_OPID, '⋱'),
    "circ"               => Token(TK_OPID, '∘'),
    "bullet"             => Token(TK_OPID, '∙'),
    "star"               => Token(TK_OPID, '⋆'),
    "div"                => Token(TK_OPID, '÷'),
    "lnot"               => Token(TK_OPID, '¬'),
    "land"               => Token(TK_OPID, '∧'),
    "lor"                => Token(TK_OPID, '∨'),
    "sim"                => Token(TK_OPID, '∼'),
    "simeq"              => Token(TK_OPID, '≃'),
    "nsim"               => Token(TK_OPID, '≁'),
    "cong"               => Token(TK_OPID, '≅'),
    "approx"             => Token(TK_OPID, '≈'),
    "ne"                 => Token(TK_OPID, '≠'),
    "neq"                => Token(TK_OPID, '≠'),
    "equiv"              => Token(TK_OPID, '≡'),
    "nequiv"             => Token(TK_OPID, '≢'),
    "prec"               => Token(TK_OPID, '≺'),
    "succ"               => Token(TK_OPID, '≻'),
    "preceq"             => Token(TK_OPID, '⪯'),
    "succeq"             => Token(TK_OPID, '⪰'),
    "dashv"              => Token(TK_OPID, '⊣'),
    "asymp"              => Token(TK_OPID, '≍'),
    "doteq"              => Token(TK_OPID, '≐'),
    "propto"             => Token(TK_OPID, '∝'),
    "barwedge"           => Token(TK_OPID, '⊼'),
    "ltimes"             => Token(TK_OPID, '⋉'),
    "rtimes"             => Token(TK_OPID, '⋊'),
    "Join"               => Token(TK_OPID, '⋈'),
    "lhd"                => Token(TK_OPID, '⊲'),
    "rhd"                => Token(TK_OPID, '⊳'),
    "unlhd"              => Token(TK_OPID, '⊴'),
    "unrhd"              => Token(TK_OPID, '⊵'),
    "vee"                => Token(TK_OPID, '∨'),
    "uplus"              => Token(TK_OPID, '⊎'),
    "wedge"              => Token(TK_OPID, '∧'),
    "veebar"             => Token(TK_OPID, '⊻'),
    "cap"                => Token(TK_OPID, '∩'),
    "cup"                => Token(TK_OPID, '∪'),
    "mid"                => Token(TK_OPID, '∣'),
    "nmid"               => Token(TK_OPID, '∤'),
    "parallel"           => Token(TK_OPID, '∥'),
    "perp"               => Token(TK_OPID, '⊥'),
    "forall"             => Token(TK_OPID, '∀'),
    "exists"             => Token(TK_OPID, '∃'),
    "nexists"            => Token(TK_OPID, '∄'),
    "lt"                 => Token(TK_OPID, '<'),
    "gt"                 => Token(TK_OPID, '>'),
    "leq"                => Token(TK_OPID, '≤'),
    "geq"                => Token(TK_OPID, '≥'),
    "ll"                 => Token(TK_OPID, '≪'),
    "gg"                 => Token(TK_OPID, '≫'),
    "lessapprox"         => Token(TK_OPID, '⪅'),
    "lesssim"            => Token(TK_OPID, '≲'),
    "eqslantless"        => Token(TK_OPID, '⪕'),
    "leqslant"           => Token(TK_OPID, '⩽'),
    "leqq"               => Token(TK_OPID, '≦'),
    "geqq"               => Token(TK_OPID, '≧'),
    "geqslant"           => Token(TK_OPID, '⩾'),
    "eqslantgtr"         => Token(TK_OPID, '⪖'),
    "gtrsim"             => Token(TK_OPID, '≳'),
    "gtrapprox"          => Token(TK_OPID, '⪆'),
    "approxeq"           => Token(TK_OPID, '≊'),
    "lessdot"            => Token(TK_OPID, '⋖'),
    "lll"                => Token(TK_OPID, '⋘'),
    "lessgtr"            => Token(TK_OPID, '≶'),
    "lesseqgtr"          => Token(TK_OPID, '⋚'),
    "lesseqqgtr"         => Token(TK_OPID, '⪋'),
    "doteqdot"           => Token(TK_OPID, '≑'),
    "leftarrow"          => Token(TK_OPID, '←'),
    "gets"               => Token(TK_OPID, '←'),
    "rightarrow"         => Token(TK_OPID, '→'),
    "to"                 => Token(TK_OPID, '→'),
    "nleftarrow"         => Token(TK_OPID, '↚'),
    "nrightarrow"        => Token(TK_OPID, '↛'),
    "leftrightarrow"     => Token(TK_OPID, '↔'),
    "nleftrightarrow"    => Token(TK_OPID, '↮'),
    "longleftarrow"      => Token(TK_OPID, '⟵'),
    "longrightarrow"     => Token(TK_OPID, '⟶'),
    "longleftrightarrow" => Token(TK_OPID, '⟷'),
    "Leftarrow"          => Token(TK_OPID, '⇐'),
    "Rightarrow"         => Token(TK_OPID, '⇒'),
    "nLeftarrow"         => Token(TK_OPID, '⇍'),
    "nRightarrow"        => Token(TK_OPID, '⇏'),
    "uparrow"            => Token(TK_OPID, '↑'),
    "downarrow"          => Token(TK_OPID, '↓'),
    "updownarrow"        => Token(TK_OPID, '↕'),
    "in"                 => Token(TK_OPID, '∈'),
    "ni"                 => Token(TK_OPID, '∉'),
    "notin"              => Token(TK_OPID, '∉'),
    "subset"             => Token(TK_OPID, '⊂'),
    "supset"             => Token(TK_OPID, '⊃'),
    "subseteq"           => Token(TK_OPID, '⊆'),
    "supseteq"           => Token(TK_OPID, '⊇'),
    "nsubseteq"          => Token(TK_OPID, '⊈'),
    "nsupseteq"          => Token(TK_OPID, '⊉'),
    "subsetneq"          => Token(TK_OPID, '⊊'),
    "supsetneq"          => Token(TK_OPID, '⊋'),
    "sqsubset"           => Token(TK_OPID, '⊏'),
    "sqsubseteq"         => Token(TK_OPID, '⊑'),
    "sqsupset"           => Token(TK_OPID, '⊐'),
    "sqsupseteq"         => Token(TK_OPID, '⊒'),
    "sqcap"              => Token(TK_OPID, '⊓'),
    "sqcup"              => Token(TK_OPID, '⊔'),
    "setminus"           => Token(TK_OPID, '∖'),
    "partial"            => Token(TK_OPID, '∂'),
    "nabla"              => Token(TK_OPID, '∇'),
    "__final__"  => Token(TK_ID, "__final__")
   )

const unicode_to_optoken = Dict(
    # integrals and big operators
    '∫' => Token(TK_INT, '∫'),
    '∬' => Token(TK_INT, '∬'),
    '∭' => Token(TK_INT, '∭'),
    '∮' => Token(TK_INT, '∮'),
    '∑' => Token(TK_BIGOP, '∑'),
    '∏' => Token(TK_BIGOP, '∏'),
    '∐' => Token(TK_BIGOP, '∐'),
    '⋂' => Token(TK_BIGOP, '⋂'),
    '⋃' => Token(TK_BIGOP, '⋃'),
    '⨆' => Token(TK_BIGOP, '⨆'),
    '⋁' => Token(TK_BIGOP, '⋁'),
    '⋀' => Token(TK_BIGOP, '⋀'),
    '⨀' => Token(TK_BIGOP, '⨀'),
    '⨂' => Token(TK_BIGOP, '⨂'),
    '⨁' => Token(TK_BIGOP, '⨁'),
    '⨄' => Token(TK_BIGOP, '⨄'),
    # operators
    '×' => Token(TK_OPID, '×'),
    '⊕' => Token(TK_OPID, '⊕'),
    '⊖' => Token(TK_OPID, '⊖'),
    '⊗' => Token(TK_OPID, '⊗'),
    '⊘' => Token(TK_OPID, '⊘'),
    '⊙' => Token(TK_OPID, '⊙'),
    '◯' => Token(TK_OPID, '◯'),
    '⨿' => Token(TK_OPID, '⨿'),
    '±' => Token(TK_OPID, '±'),
    '∓' => Token(TK_OPID, '∓'),
    '·' => Token(TK_OPID, '·'),
    '⋯' => Token(TK_OPID, '⋯'),
    '⋮' => Token(TK_OPID, '⋮'),
    '…' => Token(TK_OPID, '…'),
    '⋱' => Token(TK_OPID, '⋱'),
    '∘' => Token(TK_OPID, '∘'),
    '∙' => Token(TK_OPID, '∙'),
    '⋆' => Token(TK_OPID, '⋆'),
    '÷' => Token(TK_OPID, '÷'),
    '¬' => Token(TK_OPID, '¬'),
    '∧' => Token(TK_OPID, '∧'),
    '∨' => Token(TK_OPID, '∨'),
    '∼' => Token(TK_OPID, '∼'),
    '≃' => Token(TK_OPID, '≃'),
    '≁' => Token(TK_OPID, '≁'),
    '≅' => Token(TK_OPID, '≅'),
    '≈' => Token(TK_OPID, '≈'),
    '≠' => Token(TK_OPID, '≠'),
    '≠' => Token(TK_OPID, '≠'),
    '≡' => Token(TK_OPID, '≡'),
    '≢' => Token(TK_OPID, '≢'),
    '≺' => Token(TK_OPID, '≺'),
    '≻' => Token(TK_OPID, '≻'),
    '⪯' => Token(TK_OPID, '⪯'),
    '⪰' => Token(TK_OPID, '⪰'),
    '⊣' => Token(TK_OPID, '⊣'),
    '≍' => Token(TK_OPID, '≍'),
    '≐' => Token(TK_OPID, '≐'),
    '∝' => Token(TK_OPID, '∝'),
    '⊼' => Token(TK_OPID, '⊼'),
    '⋉' => Token(TK_OPID, '⋉'),
    '⋊' => Token(TK_OPID, '⋊'),
    '⋈' => Token(TK_OPID, '⋈'),
    '⊲' => Token(TK_OPID, '⊲'),
    '⊳' => Token(TK_OPID, '⊳'),
    '⊴' => Token(TK_OPID, '⊴'),
    '⊵' => Token(TK_OPID, '⊵'),
    '∨' => Token(TK_OPID, '∨'),
    '⊎' => Token(TK_OPID, '⊎'),
    '∧' => Token(TK_OPID, '∧'),
    '⊻' => Token(TK_OPID, '⊻'),
    '∩' => Token(TK_OPID, '∩'),
    '∪' => Token(TK_OPID, '∪'),
    '∣' => Token(TK_OPID, '∣'),
    '∤' => Token(TK_OPID, '∤'),
    '∥' => Token(TK_OPID, '∥'),
    '⊥' => Token(TK_OPID, '⊥'),
    '∀' => Token(TK_OPID, '∀'),
    '∃' => Token(TK_OPID, '∃'),
    '∄' => Token(TK_OPID, '∄'),
    '<' => Token(TK_OPID, '<'),
    '>' => Token(TK_OPID, '>'),
    '≤' => Token(TK_OPID, '≤'),
    '≥' => Token(TK_OPID, '≥'),
    '≪' => Token(TK_OPID, '≪'),
    '≫' => Token(TK_OPID, '≫'),
    '⪅' => Token(TK_OPID, '⪅'),
    '≲' => Token(TK_OPID, '≲'),
    '⪕' => Token(TK_OPID, '⪕'),
    '⩽' => Token(TK_OPID, '⩽'),
    '≦' => Token(TK_OPID, '≦'),
    '≧' => Token(TK_OPID, '≧'),
    '⩾' => Token(TK_OPID, '⩾'),
    '⪖' => Token(TK_OPID, '⪖'),
    '≳' => Token(TK_OPID, '≳'),
    '⪆' => Token(TK_OPID, '⪆'),
    '≊' => Token(TK_OPID, '≊'),
    '⋖' => Token(TK_OPID, '⋖'),
    '⋘' => Token(TK_OPID, '⋘'),
    '≶' => Token(TK_OPID, '≶'),
    '⋚' => Token(TK_OPID, '⋚'),
    '⪋' => Token(TK_OPID, '⪋'),
    '≑' => Token(TK_OPID, '≑'),
    '←' => Token(TK_OPID, '←'),
    '←' => Token(TK_OPID, '←'),
    '→' => Token(TK_OPID, '→'),
    '→' => Token(TK_OPID, '→'),
    '↚' => Token(TK_OPID, '↚'),
    '↛' => Token(TK_OPID, '↛'),
    '↔' => Token(TK_OPID, '↔'),
    '↮' => Token(TK_OPID, '↮'),
    '⟵' => Token(TK_OPID, '⟵'),
    '⟶' => Token(TK_OPID, '⟶'),
    '⟷' => Token(TK_OPID, '⟷'),
    '⇐' => Token(TK_OPID, '⇐'),
    '⇒' => Token(TK_OPID, '⇒'),
    '⇍' => Token(TK_OPID, '⇍'),
    '⇏' => Token(TK_OPID, '⇏'),
    '↑' => Token(TK_OPID, '↑'),
    '↓' => Token(TK_OPID, '↓'),
    '↕' => Token(TK_OPID, '↕'),
    '∈' => Token(TK_OPID, '∈'),
    '∉' => Token(TK_OPID, '∉'),
    '∉' => Token(TK_OPID, '∉'),
    '⊂' => Token(TK_OPID, '⊂'),
    '⊃' => Token(TK_OPID, '⊃'),
    '⊆' => Token(TK_OPID, '⊆'),
    '⊇' => Token(TK_OPID, '⊇'),
    '⊈' => Token(TK_OPID, '⊈'),
    '⊉' => Token(TK_OPID, '⊉'),
    '⊊' => Token(TK_OPID, '⊊'),
    '⊋' => Token(TK_OPID, '⊋'),
    '⊏' => Token(TK_OPID, '⊏'),
    '⊑' => Token(TK_OPID, '⊑'),
    '⊐' => Token(TK_OPID, '⊐'),
    '⊒' => Token(TK_OPID, '⊒'),
    '⊓' => Token(TK_OPID, '⊓'),
    '⊔' => Token(TK_OPID, '⊔'),
    '∖' => Token(TK_OPID, '∖'),
    '∂' => Token(TK_OPID, '∂'),
    '∇' => Token(TK_OPID, '∇')
)

const _to_mathbf = Dict(
    "A" => "𝐀",  # mathematical bold capital a
    "B" => "𝐁",  # mathematical bold capital b
    "C" => "𝐂",  # mathematical bold capital c
    "D" => "𝐃",  # mathematical bold capital d
    "E" => "𝐄",  # mathematical bold capital e
    "F" => "𝐅",  # mathematical bold capital f
    "G" => "𝐆",  # mathematical bold capital g
    "H" => "𝐇",  # mathematical bold capital h
    "I" => "𝐈",  # mathematical bold capital i
    "J" => "𝐉",  # mathematical bold capital j
    "K" => "𝐊",  # mathematical bold capital k
    "L" => "𝐋",  # mathematical bold capital l
    "M" => "𝐌",  # mathematical bold capital m
    "N" => "𝐍",  # mathematical bold capital n
    "O" => "𝐎",  # mathematical bold capital o
    "P" => "𝐏",  # mathematical bold capital p
    "Q" => "𝐐",  # mathematical bold capital q
    "R" => "𝐑",  # mathematical bold capital r
    "S" => "𝐒",  # mathematical bold capital s
    "T" => "𝐓",  # mathematical bold capital t
    "U" => "𝐔",  # mathematical bold capital u
    "V" => "𝐕",  # mathematical bold capital v
    "W" => "𝐖",  # mathematical bold capital w
    "X" => "𝐗",  # mathematical bold capital x
    "Y" => "𝐘",  # mathematical bold capital y
    "Z" => "𝐙",  # mathematical bold capital z
    "a" => "𝐚",  # mathematical bold small a
    "b" => "𝐛",  # mathematical bold small b
    "c" => "𝐜",  # mathematical bold small c
    "d" => "𝐝",  # mathematical bold small d
    "e" => "𝐞",  # mathematical bold small e
    "f" => "𝐟",  # mathematical bold small f
    "g" => "𝐠",  # mathematical bold small g
    "h" => "𝐡",  # mathematical bold small h
    "i" => "𝐢",  # mathematical bold small i
    "j" => "𝐣",  # mathematical bold small j
    "k" => "𝐤",  # mathematical bold small k
    "l" => "𝐥",  # mathematical bold small l
    "m" => "𝐦",  # mathematical bold small m
    "n" => "𝐧",  # mathematical bold small n
    "o" => "𝐨",  # mathematical bold small o
    "p" => "𝐩",  # mathematical bold small p
    "q" => "𝐪",  # mathematical bold small q
    "r" => "𝐫",  # mathematical bold small r
    "s" => "𝐬",  # mathematical bold small s
    "t" => "𝐭",  # mathematical bold small t
    "u" => "𝐮",  # mathematical bold small u
    "v" => "𝐯",  # mathematical bold small v
    "w" => "𝐰",  # mathematical bold small w
    "x" => "𝐱",  # mathematical bold small x
    "y" => "𝐲",  # mathematical bold small y
    "z" => "𝐳"   # mathematical bold small z
    )
const _to_mathcal = Dict(
    "A" => "𝒜",  # mathematical script capital a
    "B" => "ℬ",  # mathematical script capital b
    "C" => "𝒞",  # mathematical script capital c
    "D" => "𝒟",  # mathematical script capital d
    "E" => "ℰ",  # mathematical script capital e
    "F" => "ℱ",  # mathematical script capital f
    "G" => "𝒢",  # mathematical script capital g
    "H" => "ℋ",  # mathematical script capital h
    "I" => "ℐ",  # mathematical script capital i
    "J" => "𝒥",  # mathematical script capital j
    "K" => "𝒦",  # mathematical script capital k
    "L" => "ℒ",  # mathematical script capital l
    "M" => "ℳ",  # mathematical script capital m
    "N" => "𝒩",  # mathematical script capital n
    "O" => "𝒪",  # mathematical script capital o
    "P" => "𝒫",  # mathematical script capital p
    "Q" => "𝒬",  # mathematical script capital q
    "R" => "ℛ",  # mathematical script capital r
    "S" => "𝒮",  # mathematical script capital s
    "T" => "𝒯",  # mathematical script capital t
    "U" => "𝒰",  # mathematical script capital u
    "V" => "𝒱",  # mathematical script capital v
    "W" => "𝒲",  # mathematical script capital w
    "X" => "𝒳",  # mathematical script capital x
    "Y" => "𝒴",  # mathematical script capital y
    "Z" => "𝒵",  # mathematical script capital z
    "a" => "𝒶",  # mathematical script small a
    "b" => "𝒷",  # mathematical script small b
    "c" => "𝒸",  # mathematical script small c
    "d" => "𝒹",  # mathematical script small d
    "e" => "ℯ",  # mathematical script small e
    "f" => "𝒻",  # mathematical script small f
    "g" => "ℊ",  # mathematical script small g
    "h" => "𝒽",  # mathematical script small h
    "i" => "𝒾",  # mathematical script small i
    "j" => "𝒿",  # mathematical script small j
    "k" => "𝓀",  # mathematical script small k
    "l" => "𝓁",  # mathematical script small l
    "m" => "𝓂",  # mathematical script small m
    "n" => "𝓃",  # mathematical script small n
    "o" => "ℴ",  # mathematical script small o
    "p" => "𝓅",  # mathematical script small p
    "q" => "𝓆",  # mathematical script small q
    "r" => "𝓇",  # mathematical script small r
    "s" => "𝓈",  # mathematical script small s
    "t" => "𝓉",  # mathematical script small t
    "u" => "𝓊",  # mathematical script small u
    "v" => "𝓋",  # mathematical script small v
    "w" => "𝓌",  # mathematical script small w
    "x" => "𝓍",  # mathematical script small x
    "y" => "𝓎",  # mathematical script small y
    "z" => "𝓏"   # mathematical script small z
    )

const _to_mathfrak = Dict(
    "A" => "𝔄",  # mathematical fraktur capital a
    "B" => "𝔅",  # mathematical fraktur capital b
    "C" => "ℭ",  # mathematical fraktur capital c
    "D" => "𝔇",  # mathematical fraktur capital d
    "E" => "𝔈",  # mathematical fraktur capital e
    "F" => "𝔉",  # mathematical fraktur capital f
    "G" => "𝔊",  # mathematical fraktur capital g
    "H" => "ℌ",  # mathematical fraktur capital h
    "I" => "𝕴",  # mathematical fraktur capital i (bold)
    "J" => "𝔍",  # mathematical fraktur capital j
    "K" => "𝔎",  # mathematical fraktur capital k
    "L" => "𝔏",  # mathematical fraktur capital l
    "M" => "𝔐",  # mathematical fraktur capital m
    "N" => "𝔑",  # mathematical fraktur capital n
    "O" => "𝔒",  # mathematical fraktur capital o
    "P" => "𝔓",  # mathematical fraktur capital p
    "Q" => "𝔔",  # mathematical fraktur capital q
    "R" => "𝕽",  # mathematical fraktur capital r (bold)
    "S" => "𝔖",  # mathematical fraktur capital s
    "T" => "𝔗",  # mathematical fraktur capital t
    "U" => "𝔘",  # mathematical fraktur capital u
    "V" => "𝔙",  # mathematical fraktur capital v
    "W" => "𝔚",  # mathematical fraktur capital w
    "X" => "𝔛",  # mathematical fraktur capital x
    "Y" => "𝔜",  # mathematical fraktur capital y
    "Z" => "ℨ",  # mathematical fraktur capital z
    "a" => "𝔞",  # mathematical fraktur small a
    "b" => "𝔟",  # mathematical fraktur small b
    "c" => "𝔠",  # mathematical fraktur small c
    "d" => "𝔡",  # mathematical fraktur small d
    "e" => "𝔢",  # mathematical fraktur small e
    "f" => "𝔣",  # mathematical fraktur small f
    "g" => "𝔤",  # mathematical fraktur small g
    "h" => "𝔥",  # mathematical fraktur small h
    "i" => "𝔦",  # mathematical fraktur small i
    "j" => "𝔧",  # mathematical fraktur small j
    "k" => "𝔨",  # mathematical fraktur small k
    "l" => "𝔩",  # mathematical fraktur small l
    "m" => "𝔪",  # mathematical fraktur small m
    "n" => "𝔫",  # mathematical fraktur small n
    "o" => "𝔬",  # mathematical fraktur small o
    "p" => "𝔭",  # mathematical fraktur small p
    "q" => "𝔮",  # mathematical fraktur small q
    "r" => "𝔯",  # mathematical fraktur small r
    "s" => "𝔰",  # mathematical fraktur small s
    "t" => "𝔱",  # mathematical fraktur small t
    "u" => "𝔲",  # mathematical fraktur small u
    "v" => "𝔳",  # mathematical fraktur small v
    "w" => "𝔴",  # mathematical fraktur small w
    "x" => "𝔵",  # mathematical fraktur small x
    "y" => "𝔶",  # mathematical fraktur small y
    "z" => "𝔷"   # mathematical fraktur small z
    )
const _to_mathbb = Dict(
    "A" => "𝔸",  # mathematical double-struck capital a
    "B" => "𝔹",  # mathematical double-struck capital b
    "C" => "ℂ",  # mathematical double-struck capital c
    "D" => "𝔻",  # mathematical double-struck capital d
    "E" => "𝔼",  # mathematical double-struck capital e
    "F" => "𝔽",  # mathematical double-struck capital f
    "G" => "𝔾",  # mathematical double-struck capital g
    "H" => "ℍ",  # mathematical double-struck capital h
    "I" => "𝕀",  # mathematical double-struck capital i
    "J" => "𝕁",  # mathematical double-struck capital j
    "K" => "𝕂",  # mathematical double-struck capital k
    "L" => "𝕃",  # mathematical double-struck capital l
    "M" => "𝕄",  # mathematical double-struck capital m
    "N" => "ℕ",  # mathematical double-struck capital o
    "O" => "𝕆",  # mathematical double-struck capital o
    "P" => "ℙ",  # mathematical double-struck capital p
    "Q" => "ℚ",  # mathematical double-struck capital q
    "R" => "ℝ",  # mathematical double-struck capital r
    "S" => "𝕊",  # mathematical double-struck capital s
    "T" => "𝕋",  # mathematical double-struck capital t
    "U" => "𝕌",  # mathematical double-struck capital u
    "V" => "𝕍",  # mathematical double-struck capital v
    "W" => "𝕎",  # mathematical double-struck capital w
    "X" => "𝕏",  # mathematical double-struck capital x
    "Y" => "𝕐",  # mathematical double-struck capital y
    "Z" => "ℤ",  # mathematical double-struck capital z
    "a" => "𝕒",  # mathematical double-struck small a
    "b" => "𝕓",  # mathematical double-struck small b
    "c" => "𝕔",  # mathematical double-struck small c
    "d" => "𝕕",  # mathematical double-struck small d
    "e" => "𝕖",  # mathematical double-struck small e
    "f" => "𝕗",  # mathematical double-struck small f
    "g" => "𝕘",  # mathematical double-struck small g
    "h" => "𝕙",  # mathematical double-struck small h
    "i" => "𝕚",  # mathematical double-struck small i
    "j" => "𝕛",  # mathematical double-struck small j
    "k" => "𝕜",  # mathematical double-struck small k
    "l" => "𝕝",  # mathematical double-struck small l
    "m" => "𝕞",  # mathematical double-struck small m
    "n" => "𝕟",  # mathematical double-struck small n
    "o" => "𝕠",  # mathematical double-struck small o
    "p" => "𝕡",  # mathematical double-struck small p
    "q" => "𝕢",  # mathematical double-struck small q
    "r" => "𝕣",  # mathematical double-struck small r
    "s" => "𝕤",  # mathematical double-struck small s
    "t" => "𝕥",  # mathematical double-struck small t
    "u" => "𝕦",  # mathematical double-struck small u
    "v" => "𝕧",  # mathematical double-struck small v
    "w" => "𝕨",  # mathematical double-struck small w
    "x" => "𝕩",  # mathematical double-struck small x
    "y" => "𝕪",  # mathematical double-struck small y
    "z" => "𝕫"   # mathematical double-struck small z
    )

const LETTER_TRANSLATION = Dict(
    "mathbf"    => _to_mathbf,
    "mathcal"   => _to_mathcal,
    "mathbb"    => _to_mathbb,
    "mathfrak"  => _to_mathfrak
    )

# this dictionary is for printing token names in debug
const token_names = Dict(
    TK_ID           =>   "TK_ID         ",
    TK_NUM          =>   "TK_NUM        ",
    #TK_SYM         =>   "TK_SYM        ",
    TK_FUN          =>   "TK_FUN        ",
    #TK_CMD         =>   "TK_CMD        ",
    TK_ENV          =>   "TK_ENV        ",
    TK_OPID         =>   "TK_OPID       ",
    #TK_OPID2       =>   "TK_OPID2      ",
    TK_LIM          =>   "TK_LIM        ",
    TK_INT          =>   "TK_INT        ",
    TK_BIGOP        =>   "TK_BIGOP      ",
    TK_LBRACE       =>   "TK_LBRACE     ",
    TK_RBRACE       =>   "TK_RBRACE     ",
    TK_SUB          =>   "TK_SUB        ",
    TK_SUP          =>   "TK_SUP        ",
    TK_AMP          =>   "TK_AMP        ",
    TK_NEWLINE      =>   "TK_NEWLINE    ",
    TK_TEXT         =>   "TK_TEXT       ",
    TK_SQRT         =>   "TK_SQRT       ",
    TK_FRAC         =>   "TK_FRAC       ",
    TK_LEFT         =>   "TK_LEFT       ",
    TK_RIGHT        =>   "TK_RIGHT      ",
    TK_MIDDLE       =>   "TK_MIDDLE     ",
    TK_BEGIN        =>   "TK_BEGIN      ",
    TK_END          =>   "TK_END        ",
    TK_BINOM        =>   "TK_BINOM      ",
    TK_OVERSET      =>   "TK_OVERSET    ",
    TK_UNDERSET     =>   "TK_UNDERSET   ",
    TK_OVERBRACE    =>   "TK_OVERBRACE  ",
    TK_UNDERBRACE   =>   "TK_UNDERBRACE ",
    TK_OVERACCENT   =>   "TK_OVERACCENT ",
    TK_UNDERACCENT  =>   "TK_UNDERACCENT",
    TK_SPACE        =>   "TK_SPACE      ",
    TK_BIG          =>   "TK_BIG        ",
    TK_PAREN        =>   "TK_PAREN      ",
    TK_ENV_NAME     =>   "TK_ENV_NAME   ",
    TK_ENV_ATTR     =>   "TK_ENV_ATTR   ",
    TK_OPNAME       =>   "TK_OPNAME     ",
    TK_STYLE        =>   "TK_STYLE      ",
    TK_IGSTYLE      =>   "TK_IGSTYLE    ",
    TK_UNKNOWN      =>   "TK_UNKNOWN    "
   )


text:
immutable Multinomial <: DiscreteMultivariateDistribution
    n::Int
    p::Vector{Float64}

    function Multinomial(n::Integer, p::Vector{Float64})
        if n <= 0
            throw(ArgumentError("n must be a positive integer."))
        end
        if !isprobvec(p)
            throw(ArgumentError("p = $p is not a probability vector."))
        end
        @compat new(round(Int, n), p)
    end

    @compat Multinomial(n::Integer, p::Vector{Float64}, ::NoArgCheck) = new(round(Int, n), p)
    @compat Multinomial(n::Integer, k::Integer) = new(round(Int, n), fill(1.0 / k, k))
end

# Parameters

ncategories(d::Multinomial) = length(d.p)
length(d::Multinomial) = ncategories(d)
probs(d::Multinomial) = d.p
ntrials(d::Multinomial) = d.n

params(d::Multinomial) = (d.n, d.p)

# Statistics

mean(d::Multinomial) = d.n .* d.p

function var(d::Multinomial)
    p = probs(d)
    k = length(p)
    n = ntrials(d)

    v = Array(Float64, k)
    for i = 1:k
        @inbounds pi = p[i]
        v[i] = n * pi * (1.0 - pi)
    end
    v
end

function cov(d::Multinomial)
    p = probs(d)
    k = length(p)
    n = ntrials(d)

    C = Array(Float64, k, k)
    for j = 1:k
        pj = p[j]
        for i = 1:j-1
            @inbounds C[i,j] = - n * p[i] * pj
        end

        @inbounds C[j,j] = n * pj * (1.0-pj)
    end

    for j = 1:k-1
        for i = j+1:k
            @inbounds C[i,j] = C[j,i]
        end
    end
    C
end

function mgf(d::Multinomial, t::AbstractVector)
    p = probs(d)
    n = ntrials(p)
    s = 0.0
    for i in 1:length(p)
        s += p[i] * exp(t[i])
    end
    return s^n
end

function cf(d::Multinomial, t::AbstractVector)
    p = probs(d)
    n = ntrials(d)
    s = 0.0 + 0.0im
    for i in 1:length(p)
        s += p[i] * exp(im * t[i])
    end
    return s^n
end

# Evaluation

function insupport{T<:Real}(d::Multinomial, x::AbstractVector{T})
    k = length(d)
    length(x) == k || return false
    s = 0.0
    for i = 1:k
        @inbounds xi = x[i]
        if !(isinteger(xi) && xi >= 0)
            return false
        end
        s += xi
    end
    return s == ntrials(d)  # integer computation would not yield truncation errors
end

function _logpdf{T<:Real}(d::Multinomial, x::AbstractVector{T})
    p = probs(d)
    n = ntrials(d)
    s = lgamma(n + 1.0)
    t = 0
    for i = 1:length(p)
        @inbounds xi = x[i]
        @inbounds pi = p[i]
        t += xi
        s -= lgamma(xi + 1.0)
        @inbounds s += xi * log(pi)
    end
    return ifelse(t == n, s, -Inf)::Float64
end

# Sampling

_rand!{T<:Real}(d::Multinomial, x::AbstractVector{T}) = multinom_rand!(ntrials(d), probs(d), x)

sampler(d::Multinomial) = MultinomialSampler(ntrials(d), probs(d))

## Fit model

immutable MultinomialStats <: SufficientStats
    n::Int  # number of trials in each experiment
    scnts::Vector{Float64}  # sum of counts
    tw::Float64  # total sample weight

    @compat MultinomialStats(n::Int, scnts::Vector{Float64}, tw::Real) = new(n, scnts, Float64(tw))
end

function suffstats{T<:Real}(::Type{Multinomial}, x::Matrix{T})
    K = size(x, 1)
    n::T = zero(T)
    scnts = zeros(K)

    for j = 1:size(x,2)
        nj = zero(T)
        for i = 1:K
            @inbounds xi = x[i,j]
            @inbounds scnts[i] += xi
            nj += xi
        end

        if j == 1
            n = nj
        elseif nj != n
            error("Each sample in X should sum to the same value.")
        end
    end
    MultinomialStats(n, scnts, size(x,2))
end

function suffstats{T<:Real}(::Type{Multinomial}, x::Matrix{T}, w::Array{Float64})
    length(w) == size(x, 2) || throw(ArgumentError("Inconsistent argument dimensions."))

    K = size(x, 1)
    n::T = zero(T)
    scnts = zeros(K)
    tw = 0.

    for j = 1:size(x,2)
        nj = zero(T)
        @inbounds wj = w[j]
        tw += wj
        for i = 1:K
            @inbounds xi = x[i,j]
            @inbounds scnts[i] += xi * wj
            nj += xi
        end

        if j == 1
            n = nj
        elseif nj != n
            error("Each sample in X should sum to the same value.")
        end
    end
    MultinomialStats(n, scnts, tw)
end

fit_mle(::Type{Multinomial}, ss::MultinomialStats) = Multinomial(ss.n, ss.scnts * inv(ss.tw * ss.n))

function fit_mle{T<:Real}(::Type{Multinomial}, x::Matrix{T})
    ss = suffstats(Multinomial, x)
    Multinomial(ss.n, multiply!(ss.scnts, inv(ss.tw * ss.n)))
end

function fit_mle{T<:Real}(::Type{Multinomial}, x::Matrix{T}, w::Array{Float64})
    ss = suffstats(Multinomial, x, w)
    Multinomial(ss.n, multiply!(ss.scnts, inv(ss.tw * ss.n)))
end


text:
# Autogenerated wrapper script for libsodium_jll for x86_64-apple-darwin
export libsodium

JLLWrappers.@generate_wrapper_header("libsodium")
JLLWrappers.@declare_library_product(libsodium, "@rpath/libsodium.23.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libsodium,
        "lib/libsodium.23.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# State: 4 continuous variables
#   cart position:      ∈ [-4.8, 4.8]
#   cart speed:         ∈ [-Inf, Inf]
#   pole angle:         ∈ [-24°, 24°]
#   pole rotation rate: ∈ [-Inf, Inf]

# Action: Push cart right (1) or left (2)

@with_kw struct CartPole
    γ::Float64 = 1.0
    force_magnitude::Float64 = 10.0
    pole_length::Float64 = 1.0  # Distance from pole to mass at end of pole
    mass_cart::Float64 = 1.0
    mass_pole::Float64 = 0.1
    gravity::Float64 = 9.8
    Δt::Float64 = 0.02  # Time to advance by with each step
    Δx_max::Float64 = 4.8  # Max position deviation
    Δθ_max::Float64 = deg2rad(24)  # Max angular deviation
end

# const CART_POLE = CartPole(1.0, 10.0, 1.0, 1.0, 0.1, 9.8, 0.02, 4.8, deg2rad(24))

struct CartPoleState
    x::Float64 # cart pos [m]
    v::Float64 # cart speed [m/s]
    θ::Float64 # pole angle [rad]
    ω::Float64 # pole rotation range [rad/s]
end

n_actions(mdp::CartPole) = 2
discount(mdp::CartPole) = mdp.γ
ordered_actions(mdp::CartPole) = [1,2]

Base.vec(s::CartPoleState) = [s.x, s.v, s.θ, s.ω]

function generate_start_state(mdp::CartPole)
    U = Uniform(-0.05,0.05)
    return CartPoleState(rand(U), rand(U), rand(U), rand(U))
end

function is_terminal(mdp::CartPole, s::CartPoleState)
    return abs(s.x) > mdp.Δx_max || abs(s.θ) > mdp.Δθ_max
end

function cart_pole_transition(mdp::CartPole, s::CartPoleState, a::Int)
    @assert(1 <= a <= 2)

    if is_terminal(mdp, s)
        return s
    end

    x, v, θ, ω = s.x, s.v, s.θ, s.ω
    force = a == 1 ? mdp.force_magnitude : -mdp.force_magnitude

    cθ, sθ = cos(θ), sin(θ)
    half_length = mdp.pole_length / 2
    mass_tot = mdp.mass_cart + mdp.mass_pole
    temp = (force + half_length*sθ*ω^2) / mass_tot

    # Pole acceleration
    pole_α = (mdp.gravity*sθ - temp*cθ) / (half_length * (4.0/3.0 - mdp.mass_pole * cθ^2 / mass_tot))

    # Cart acceleration
    cart_a = temp - half_length * mdp.mass_pole * pole_α * cθ / mass_tot

    # Euler integration
    Δt = mdp.Δt
    x += Δt*v
    v += Δt*cart_a
    θ += Δt*ω
    ω += Δt*pole_α

    return CartPoleState(x,v,θ,ω)
end

function reward(mdp::CartPole, s::CartPoleState, a::Int)
    @assert(1 <= a <= 2)
    if is_terminal(mdp, s)
        return 0.0
    else
        return 1.0
    end
end

function MDP(mdp::CartPole; γ::Float64=mdp.γ)
    return MDP(
            γ,
            nothing, # no ordered states
            ordered_actions(mdp),
            nothing, # no probabilistic transition function
            (s,a) -> reward(mdp, s, a),
            (s, a)->begin
                s′ = cart_pole_transition(mdp, s, a)
                r = reward(mdp, s, a)
                return (s′, r)
            end
        )
end


text:
abstract Struct

# Base.keys(t::Struct)       = fieldnames(t)
# Base.values(t::Struct)     = [getfield(t,i) for i in fieldnames(t)]
# Base.length(t::Struct)     = length(fieldnames(t))
# Base.start(t::Struct)      = 1
# Base.done(t::Struct, iter) = iter == length(fieldnames(t))+1
# Base.next(t::Struct, iter) = ((fieldnames(t)[iter], getfield(t,iter)), iter += 1)
# Base.endof(t::Struct)      = length(t)
# Base.last(t::Struct)       = (fieldnames(t)[end], t[end])

# Base.getindex(t::Struct, i::Int)              = getfield(t,i)
# Base.getindex(t::Struct, i::UnitRange{Int64}) = slice(t,i)
# Base.getindex(t::Struct, i::Symbol)           = getfield(t,i)
# Base.getindex(t::Struct, i::Symbol, default)  = get(t,i,default)
# Base.get(t::Struct, i::Symbol, default)       = i in keys(t) ? t[i] : default

Base.writemime(io::IO, mime::MIME"text/plain", t::Struct) = show(io,t)

function Base.show(io::IO, t::Struct; cols=100)
    shown_set = get(task_local_storage(), :SHOWNSET, nothing)
    if shown_set === nothing
        shown_set = ObjectIdDict()
        task_local_storage(:SHOWNSET, shown_set)
    end
    t in keys(shown_set) && (print(io, "#= circular reference =#"); return)

    try
        shown_set[t] = true
        fields = fieldnames(t)
        ks = map(string,fields)
        keylen = maximum(map(length,ks))

        for (i,k) in enumerate(fields)
            print(io, Base.rpad(ks[i], keylen))
            print(io, " = ")
            val = Base.with_output_limit(()->sprint(show, getfield(t,k)))
            val = Base._truncate_at_width_or_chars(val, cols - keylen, "\r\n")
            print(io, val)
            print(io, '\n')
        end
    finally
        delete!(shown_set, t)
    end
end

# function show(io::IO, t::Struct)
#     shown_set = get(task_local_storage(), :SHOWNSET, nothing)
#     if shown_set === nothing
#         shown_set = ObjectIdDict()
#         task_local_storage(:SHOWNSET, shown_set)
#     end
#     t in keys(shown_set) && (print(io, "#= circular reference =#"); return)

#     try
#         shown_set[t] = true
#         # if compact
#         #     # show in a Julia-syntax-like form: Dict(k=>v, ...)
#         #     if isempty(t)
#         #         print(io, typeof(t), "()")
#         #     else
#         #         if isleaftype(K) && isleaftype(V)
#         #             print(io, typeof(t).name)
#         #         else
#         #             print(io, typeof(t))
#         #         end
#         #         print(io, '(')
#         #         first = true
#         #         n = 0
#         #         for (k, v) in t
#         #             first || print(io, ',')
#         #             first = false
#         #             show(io, k)
#         #             print(io, "=>")
#         #             show(io, v)
#         #             n+=1
#         #             limit && n >= 10 && (print(io, "…"); break)
#         #         end
#         #         print(io, ')')
#         #     end
#         #     return
#         # end

#         # Otherwise show more descriptively, with one line per key/value pair
#         # rows, cols = sz
#         print(io, summary(t))
#         # isempty(t) && return
#         print(io, ":")
#         # if limit
#         #     rows < 2   && (print(io, " …"); return)
#         #     cols < 12  && (cols = 12) # Minimum widths of 2 for key, 4 for value
#         #     cols -= 6 # Subtract the widths of prefix "  " separator " => "
#         #     rows -= 2 # Subtract the summary and final ⋮ continuation lines

#         #     # determine max key width to align the output, caching the strings
#         #     ks = Array(AbstractString, min(rows, length(t)))
#         #     keylen = 0
#         #     for (i, k) in enumerate(keys(t))
#         #         i > rows && break
#         #         ks[i] = sprint(show, k)
#         #         keylen = clamp(length(ks[i]), keylen, div(cols, 3))
#         #     end
#         # end

#         for (i, (k, v)) in enumerate(t)
#             print(io, "\n  ")
#             # limit && i > rows && (print(io, rpad("⋮", keylen), " => ⋮"); break)

#             # if limit
#                 # key = rpad(_truncate_at_width_or_chars(ks[i], keylen, "\r\n"), keylen)
#             # else
#             key = sprint(show, k)
#             # end
#             print(io, key)
#             print(io, " = ")

#             # if limit
#             #     val = with_output_limit(()->sprint(show, v))
#             #     val = _truncate_at_width_or_chars(val, cols - keylen, "\r\n")
#             #     print(io, val)
#             # else
#             show(io, v)
#             # end
#         end
#     finally
#         delete!(shown_set, t)
#     end
# end

function Base.(:(==))(a::Struct, b::Struct)
    a === b                && return true
    typeof(a) != typeof(b) && return false
    for k in fieldnames(a)
        if getfield(a,k) != getfield(b,k)
            return false
        end
    end
    true
end

function Base.hash(t::Struct, h::UInt)
    h_ = 17
    for k in fieldnames(t)
        h_ = 23*h_ + hash(getfield(t,k), h)
    end
    h_
end

trans(expr::Expr)                     = trans(Val{expr.head}, expr)                 #
trans(sym::Symbol)                    = (sym, :Any, nothing)                        # x
trans(::Type{Val{:(::)}}, expr::Expr) = (expr.args[1], expr.args[2], nothing)       # x::T

function trans(::Type{Val{:kw}}, expr::Expr)                                        # x::T = v
    sym, typ = trans(expr.args[1])
    (sym, typ, expr.args[2])
end

function make_type(base::Symbol, fields::Vector{Symbol})
    name = symbol(base,'_',join(fields,'_'))
    if !isdefined(Element,name)
        types = [symbol('T',i) for i in 1:length(fields)]
        tfields = [Expr(:(::), fields[i], types[i]) for i in 1:length(fields)]
        def = Expr(:type, true, Expr(:(<:), Expr(:curly, name, types...), base), Expr(:block, tfields...))
        eval(def)
    end
    return name
end

function make_tuple(base::Symbol, exprs::Union{Symbol,Expr}...)
    local flag
    n = length(exprs)

    fields = Array{Symbol}(n)
    values = Array{Any}(n)
    types  = Array{Any}(n)

    for i in 1:n
        fields[i],T,v = trans(exprs[i])
        if i == 1; flag = is(v,nothing) end
        is(v,nothing) != flag && error("all values must be specified during construction")
        types[i] = esc(T)
        values[i] = is(T,:Any) ? esc(v) : esc(Expr(:call, :convert, T, v))
    end

    # sort fields alphabetically
    ind = sortperm(fields)
    fields = fields[ind]
    types = types[ind]
    values = values[ind]

    name = make_type(base, fields)  # create type if nonexistent

    if flag
        return Expr(:curly, name, types...)
    else
        return Expr(:call, name, values...)
    end
end

abstract Atom <: Struct

macro atom(exprs::Expr...)
    make_tuple(:Atom, exprs...)
end

macro struct(exprs::Expr...)
    make_tuple(:Struct, exprs...)
end

@generated function merge(a::Struct, b::Struct)
    merge_impl(a,b)
end

function merge_impl{T1<:Struct,T2<:Struct}(a::Type{T1}, b::Type{T2})
    fields = sort(unique([fieldnames(a);fieldnames(b)]))
    T = symbol(last(split(string(super(a)),'.')))
    name = make_type(T, fields)
    vals = [k in fieldnames(b) ? Expr(:call,:getfield,:b,QuoteNode(k)) : Expr(:call,:getfield,:a,QuoteNode(k)) for k in fields]
    Expr(:call, name, vals...)
end

# @doc doc"""
# Create a new Struct with the new value set on it, either overwriting
# the old value or appending a new value.
# This copies the underlying data.
# """ ->
# function setindex{V}( t::Struct, key::Symbol, val::V)
#     nt = eval( make_tuple( super(typeof(t)), [ Expr( :(=>), key, val  )]))
#     return merge( t, nt )
# end

# @doc doc"""
# Create a new Struct with the specifed element removed.
# This copies the underlying data.
# """ ->
# function delete( t::Struct, key::Symbol )
#     nms = filter( x->x!=key, fieldnames( t ) )
#     name = create_tuple( super(typeof(t)), nms )
#     vals = [ getindex( t, nm ) for nm in nms ]
#     return eval( Expr( :call, name, vals... ) )
# end

max_stars_repo_name: johncsnyder/Element


text:
######################################################################
# fov.jl
#
# field of view sensor
######################################################################

struct FOV <: Sensor
    cone_width::Float64
    alpha::Float64          # mistake rate

    blind_distance::Float64     # distance in which any obs is equal likely
end
function FOV(cone_width::Real, alpha::Float64)
    return FOV(float(cone_width), alpha, 0.0)
end
FOV() = FOV(120.0, 0.1, 0.0)

# 1 means it is in the field of view of front antenna
# 0 means it is not
function observe(theta::TargetTuple, s::FOV, p::Pose)

    # ensure bearing is reflected across
    rel_bearing = fit_180(p[3] - true_bearing(p, theta))
    if rel_bearing < 0.0
        rel_bearing = -1.0 * rel_bearing
    end

    a1 = s.cone_width / 2.0
    a2 = 180.0 - a1

    prob_in_view = 0.0
    if rel_bearing < a1
        prob_in_view = 1.0 - s.alpha
    elseif rel_bearing < a2
        prob_in_view = 0.5
    else
        prob_in_view = s.alpha
    end

    # if we are too close, then prob_in_view = 0.5
    dx = p[1] - theta[1]
    dy = p[2] - theta[2]
    if (dx*dx + dy*dy) < s.blind_distance*s.blind_distance
        prob_in_view = 0.5
    end

    o = (rand() < prob_in_view) ? 1 : 0
    return o
end

# o technically has to be a float, but we only have two obs:
#  1 = in field of view
#  0 = not in field of view
function O(s::FOV, theta::LocTuple, p::Pose, o)

    # determine relative bearing and fix it in 0 to 180
    rel_bearing = fit_180(p[3] - true_bearing(p, theta))
    if rel_bearing < 0.0
        rel_bearing = -1.0 * rel_bearing
    end

    a1 = s.cone_width / 2.0
    a2 = 180.0 - a1

    prob_in_view = 0.0  # probability in view of front antenna

    if rel_bearing < a1
        prob_in_view = 1.0 - s.alpha
    elseif rel_bearing < a2
        prob_in_view = 0.5
    else
        prob_in_view = s.alpha
    end

    # if we are too close, then prob_in_view = 0.5
    dx = p[1] - theta[1]
    dy = p[2] - theta[2]
    if (dx*dx + dy*dy) < s.blind_distance*s.blind_distance
        prob_in_view = 0.5
    end

    # prob that it's one (in view of front antenna)
    ret_val = (o == 1) ? prob_in_view : (1.0 - prob_in_view)

    return ret_val
end

get_prob(s::FOV, theta::LocTuple, xp::Pose, o) = O(s, theta, xp, o)
export get_prob


text:
function rtuples(T::Type{<:Number}, D::Int, k::Int)::Vector{NTuple{D,T}}
    return [tuple(rand(T,D)...) for i in 1:k];
end
function rarrays(T::Type{<:Number}, D::Int, k::Int)::Vector{Vector{T}}
    return [rand(T,D) for i in 1:k];
end
function totuples(data)
    return [tuple(e...) for e in data];
end
function get_perfect_data(; n = 100)
    cs = [(-5,-5),(0,5),(5,-5)]
    dt = Vector{Tuple{Float64,Float64}}();
    xr = randn(n);
    yr = randn(n);
    for rd in 1:n
        for c in cs
            push!(dt, (xr[rd] + c[1], yr[rd] + c[2]))
        end
    end
    return dt;
end
function center_coords(ps)
    s = ps[1];
    for i in 2:length(ps)
         s = s .+ ps[i]
    end
    return s ./ length(ps);
end
function assign(p::T, cs::Vector{T}) where T
    findmin(distance(p,cs))[2];
end
function assign(ps::Vector{T}, cs::Vector{T}) where T
    return [findmin(distance(p,cs))[2] for p in ps];
end
function distance(p1, p2)
    return sqrt(sum((p1 .- p2).^2))
end
function distance(p1::T, ps::Vector{T}) where T
    return [distance(p1,p) for p in ps];
end
function ave_distance(p, ps)
    length(ps) == 0 && return zero(eltype(p))
    ds = distance(p,ps)
    return sum(ds)/length(ps);
end
function generate_test_2d_data(k; r = 35, n = 1000)
    cs = [(rand()*r,rand()*r) for i in 1:k];
    dt = Vector{Tuple{Float64,Float64}}();
    xr = randn(n);
    yr = randn(n);
    for rd in 1:n
        for c in cs
            push!(dt, (xr[rd] + c[1], yr[rd] + c[2]))
        end
    end
    return dt;
end
export fake_data
function fake_data()
    return get_perfect_data();
end
function remove(col::Vector, td)
    col[[1:(td-1);(td+1):end]];
end
function remove!(col::Vector, td)
    temp = remove(col,td);
    pop!(col);
    col .= temp;
    return col;
end


text:
using jlmie
using Plots

# settings
# structure
nmat = 4  # refractive index of a sphere
radius = 100*1e-9

# environment
nenv = 1.00

# wavelenth
lbd0 = (400:1000)*1e-9; lbdp = lbd0.*1e9;

# other settings
nmax = -1  # -1: namx large enough (determined from x)

# main
m,x = jlmie_mx(nmat,radius,lbd0,nenv)
Qsca = jlmie_Qsca(m,x,nmax)
println("calculation completed")

# show results
# validation can be done by comparing to ITMO Mie Calculator
# https://physics.itmo.ru/en/mie
plt = plot(lbdp,Qsca,
    xlabel       = "Wavelength (nm)",
    ylabel       = "Scattering efficiency",
    legend       = false,
    # label        ="LABEL",
    # xlims        =(-3,3),
    # ylims        =(-3,3),
    # aspect_ratio =0.5,
    # title        ="TITLE",
    # linecolor    =:blue,
    # linewidth    =5,
    # linestyle    =:dot,
    # size         =(400,300),
) # https://qiita.com/I_ppp/items/dca3552affa6a672e4bd
display(plt)

# save
if true
    runningfilename = splitext(splitpath(@__FILE__)[end])[1]
    outfilename = ".\\result\\" * runningfilename
    savefig(plt,outfilename)
    println("Graph has been saved in result directory")
end

max_stars_repo_name: Hinamoooon/jlmie

max_issues_repo_name: Hinamoooon/jlmie

max_forks_repo_name: Hinamoooon/jlmie


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test

@bench @testset "NLEIGS: Naive Gun" begin
    nep = nep_gallery("nlevp_native_gun")
    square = [-1-1im; -1+1im; 1+1im; 1-1im]
    λ, X, _ = nleigs(nep, 150^2 .+ 200.0*square, logger=displaylevel)
    verify_lambdas(1, nep, λ, X)
end