text:
"""
AXLeastSquares(y, X)

An Ordinary Least Squares implementation.
"""

struct AXLeastSquares <: AXEstimator
    β::Array{Float64} # Coefficient
    y::Matrix{Float64} # Response
    X::Matrix{Float64} # Covariates

    # Constructor Function
    function AXLeastSquares(y::Matrix{Float64}, X::Matrix{Float64})
        β = inv(X' * X) * (X' * y)
        new(β, y, X)
    end
end #AXLeastSquares

# Methods ======================================================================
"""
infer(fit::AXLeastSquares; heteroskedastic, print_df)

A method to calculate standard errors of an AXLeastSquares object.
"""
function infer(fit::AXLeastSquares; heteroskedastic::Bool=false, print_df::Bool=true)
    # Retrieve necessary parameters
    N = length(fit.y)
    K = size(fit.X, 2) # number of covariates

    # Calculate the variance-covariance matrix
    resid  = fit.y - predict(fit)
    XX     = fit.X' * fit.X
    residX = fit.X .* resid
    if heteroskedastic
        vcv = inv(XX) * (residX' * residX) * inv(XX) .* (N / (N-K))
    else
        # homoskedasticity
        vcv = inv(XX) .* (resid' * resid) ./ (N-K)
    end

    se = sqrt.(diag(vcv))
    t_stat = fit.β ./ se
    p_val = 2 * cdf.(Normal(), -abs.(t_stat))
    r2 = 1 - sum(resid.^2)/sum((fit.y.-mean(fit.y)).^2)

    # Print estimates
    if print_df
        out_df = DataFrame(hcat(fit.β, se, t_stat, p_val, r2), :auto)
        rename!(out_df, ["coef", "se", "t-stat", "p-val", "R-sqaure"])
        display(out_df)
    end

    # Organize and return output
    output = (β = fit.β, se = se, t = t_stat, p = p_val, r = r2)
    return output
end #infer.AXLeastSquares


text:
let n::Int64=0
	global increment_n() = (n+=1)
	global reset_n() = (n=0)
end

let file::String = ""
	global get_nwk_file() = file
	global set_nwk_file(s) = (file=s)
end

"""
	read_tree(
		nwk_filename::AbstractString;
		node_data_type=DEFAULT_NODE_DATATYPE, force_new_labels=false
	)
	read_tree(
		io::IO;
		node_data_type=DEFAULT_NODE_DATATYPE, force_new_labels=false
	)

Read Newick file and create a `Tree{node_data_type}` object from it. The input file can
contain multiple Newick strings on different lines. The output will then be an array of
`Tree` objects.

`node_data_type` must be a subtype of `TreeNodeData`, and must have a *callable default outer
constructor*: the call `node_data_type()` must exist and return a valid instance of
`node_data_type`. See `?TreeNodeData` for implemented types.

Use `force_new_labels=true` to force the renaming of all internal nodes.

If you have a variable containing a Newick string and want to build a tree from it,
use `parse_newick_string` instead.

## Note on labels
The `Tree` type identifies nodes by their labels. This means that labels have to be unique.
For this reason, the following is done when reading a tree:
- if an internal node does not have a label, a unique one will be created of the form
	`"NODE_i"`
- if a node has a label that was already found before in the tree, a random identifier
	will be appended to it to make it unique. Note that the identifier is created using
	`randstring(8)`, unicity is technically not guaranteed.
- if `force_new_labels` is used, a unique identifier is appended to node labels
- if node labels in the Newick file are identified as confidence/bootstrap values, a random
	identifier is appended to them, even if they're unique in the tree. See
	`?TreeTools.isbootstrap` to see which labels are identified as confidence values.
"""
function read_tree(
	io::IO;
	node_data_type=DEFAULT_NODE_DATATYPE, force_new_labels=false
)
	trees = map(eachline(io)) do line
		t = parse_newick_string(line; node_data_type, force_new_labels)
		check_tree(t)
		t
	end
	return length(trees) == 1 ? trees[1] : trees
end
function read_tree(
	nwk_filename::AbstractString;
	node_data_type=DEFAULT_NODE_DATATYPE, force_new_labels=false
)
	return open(nwk_filename, "r") do io
		read_tree(io; node_data_type, force_new_labels)
	end
end
"""
	parse_newick_string(
		nw::AbstractString;
		node_data_type=DEFAULT_NODE_DATATYPE, force_new_labels=false
	)

Parse newick string into a tree. See `read_tree` for more informations.
"""
function parse_newick_string(
	nw::AbstractString;
	node_data_type=DEFAULT_NODE_DATATYPE, force_new_labels=false,
)
	@assert nw[end] == ';' "Newick string does not end with ';'"

	reset_n()
	root = parse_newick(nw; node_data_type)
	tree = node2tree(root; force_new_labels)
	check_tree(tree)
	return tree
end

"""
	read_newick(nwk_filename::AbstractString)

Read Newick file `nwk_filename` and create a graph of `TreeNode` objects in the process.
  Return the root of said graph.
  `node2tree` or `read_tree` must be called to obtain a `Tree` object.
"""
function read_newick(nwk_filename::AbstractString; node_data_type=DEFAULT_NODE_DATATYPE)
	@assert node_data_type <: TreeNodeData
	set_nwk_file(nwk_filename)

	f = open(nwk_filename)
	nw = readlines(f)
	close(f)
	if length(nw) > 1
		error("File $nwk_filename has more than one line.")
	elseif length(nw) == 0
		error("File $nwk_filename is empty")
	end
	nw = nw[1]
	if nw[end] != ';'
		error("File $nwk_filename does not end with ';'")
	end
	nw = nw[1:end-1]

	reset_n()
	root = parse_newick(nw; node_data_type)
	return root
end

"""
	parse_newick(nw::AbstractString; node_data_type=DEFAULT_NODE_DATATYPE)

Parse newick string into a `TreeNode`.
"""
function parse_newick(nw::AbstractString; node_data_type=DEFAULT_NODE_DATATYPE)
	if isempty(nw)
		error("Cannot parse empty Newick string.")
	end
	reset_n()
	root = TreeNode(node_data_type())
	parse_newick!(nw, root, node_data_type)
	root.isroot = true # Rooting the tree with outer-most node of the newick string
	root.tau = missing
	return root
end

"""
	parse_newick!(nw::AbstractString, root::TreeNode)

Parse the tree contained in Newick string `nw`, rooting it at `root`.
"""
function parse_newick!(nw::AbstractString, root::TreeNode, node_data_type)

	# Setting isroot to false. Special case of the root is handled in main calling function
	root.isroot = false
	# Getting label of the node, after last parenthesis
	parts = map(x->String(x), split(nw, ")"))
	lab, tau = nw_parse_name(String(parts[end]))
	if lab == ""
		lab = "NODE_$(increment_n())"
	elseif length(lab) > 1 && lab[1:2] == "[&" # Deal with extended newick annotations
		lab = "NODE_$(increment_n())" * lab
	end

	root.label, root.tau = (lab,tau)

	if length(parts) == 1 # Is a leaf. Subtree is empty
		root.isleaf = true
	else # Has children
		root.isleaf = false
		if parts[1][1] != '('
			println(parts[1][1])
			@error "Parenthesis mismatch around $(parts[1]). This may be caused by spaces in the newick string."
			error("$(get_nwk_file()): incorrect Newick format.")
		else
			parts[1] = parts[1][2:end] # Removing first bracket
		end
		children = join(parts[1:end-1], ")") # String containing children, now delimited with ','
		l_children = nw_parse_children(children) # List of children (array of strings)

		for sc in l_children
			nc = TreeNode(node_data_type())
			parse_newick!(sc, nc, node_data_type) # Will set everything right for subtree corresponding to nc
			nc.anc = root
			push!(root.child, nc)
		end
	end

end

"""
	nw_parse_children(s::AbstractString)

Idea from http://stackoverflow.com/a/26809037
Split a string of children in newick format to an array of strings.

## Example
`"A,(B,C),D"` --> `["A","(B,C)","D"]`
"""
function nw_parse_children(s::AbstractString)
	parcount = 0
	annotation = false # For reading extended newick grammar, with [&key=val] after label
	l_children = []
	current = ""
	cstart = 1
	cend = 1
	for (i,c) in enumerate("$(s),")
		if c == '['
			annotation = true
		elseif c == ']'
			annotation = false
		elseif c == ',' && parcount == 0 && annotation == false
			cend = i-1
			push!(l_children, s[cstart:cend])
			cstart = i+1
		else
			if c == '('
				parcount +=1
			elseif c == ')'
				parcount -=1
			end
		end
	end
	return l_children
end

"""
	nw_parse_name(s::AbstractString)

Parse Newick string of child into name and time to ancestor. Default value for missing time is `missing`.
"""
function nw_parse_name(s::AbstractString)
	temp = split(s, ":")
	if occursin(':', s) # Node has a time
		if length(temp) == 2 # Node also has a name, return both
			tau = (tau = tryparse(Float64,temp[2]); isnothing(tau) ? missing : tau) # Dealing with unparsable times
			return string(temp[1]), tau
		else # Return empty name
			tau = (tau = tryparse(Float64,temp[1]); isnothing(tau) ? missing : tau) # Dealing with unparsable times
			return "", tau
		end
	else # Node does not have a time, return string as name
		return s, missing
	end
end

ext: jl


text:
using CUDA

using Test

"Dummy kernel doing 100 FMAs."
function kernel_100fma(a, b, c, out)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    @inbounds a_val = a[i]
    @inbounds b_val = b[i]
    @inbounds c_val = c[i]

    for j in 1:33
        a_val = CUDA.fma(a_val, b_val, c_val)
        b_val = CUDA.fma(a_val, b_val, c_val)
        c_val = CUDA.fma(a_val, b_val, c_val)
    end

    @inbounds out[i] = CUDA.fma(a_val, b_val, c_val)

    return
end

function peakflops(n::Integer=10000, dev::CuDevice=CuDevice(0))
    device!(dev) do
        dims = (n, n)
        a = round.(rand(Float32, dims) * 100)
        b = round.(rand(Float32, dims) * 100)
        c = round.(rand(Float32, dims) * 100)
        out = similar(a)

        d_a = CuArray(a)
        d_b = CuArray(b)
        d_c = CuArray(c)
        d_out = CuArray(out)

        len = prod(dims)
        threads = min(len, 1024)
        blocks = len ÷ threads

        # warm-up
        @cuda kernel_100fma(d_a, d_b, d_c, d_out)
        synchronize()

        secs = CUDA.@elapsed begin
            @cuda blocks=blocks threads=threads kernel_100fma(d_a, d_b, d_c, d_out)
        end
        flopcount = 200*len
        flops = flopcount / secs

        return flops
    end
end

f1(x)=CUDA.sqrt(CUDA.sin(x^2))-CUDA.sqrt(CUDA.cos(x^2))

CUDA.@time sum(f1.(cuv1))

myf(x)=sqrt()

v1=fill(1.2345,m*m)

@elapsed @avx sqrt.(v1)

@elapsed @avx for i in 1:m*m
    v1[i]=sqrt(v1[i])
end

v1

v1=rand(cutype,n*n)

CUDA.device!(1)
using CUDA
using LinearAlgebra
BLAS.set_num_threads(1)
const m=5000
const n=5000
const p=5000
tflop=2*m*n*p *10^-12
cutype=Float64
v1=rand(cutype,n,n)
v2=rand(cutype,n,n)
v3=rand(cutype,n,n)

function A_mul_B!(𝐂, 𝐀, 𝐁)
    @simd for m ∈ axes(𝐀,1), n ∈ axes(𝐁,2)
        𝐂ₘₙ = zero(eltype(𝐂))
    @simd for k ∈ axes(𝐀,2)
            𝐂ₘₙ += 𝐀[m,k] * 𝐁[k,n]
        end
        𝐂[m,n] = 𝐂ₘₙ
    end
end

cuv1=CUDA.rand(cutype,m,n);

cuv2=CUDA.rand(cutype,n,p);
cuv3=CUDA.rand(cutype,m,p);

GC.gc()

device_gc()

afv1=rand(AFArray{cutype},m,n);
afv2=rand(AFArray{cutype},n,p);
afv3=rand(AFArray{cutype},m,p);

gput=@elapsed CUDA.@sync afv3=afv1*afv2

finalize(afv1)
finalize(afv2)
finalize(afv3)

1
@elapsed v3.=v1.*v2

gput=CUDA.@elapsed CUDA.@sync mul!(cuv3,cuv1,cuv2)
gput=CUDA.@elapsed CUDA.@sync (cuv3=cuv1*cuv2;)
GC.gc(true)
CUDA.reclaim()

cput=@elapsed A_mul_B!(v3,v1,v2)

finalize(afv1)
finalize(afv2)
finalize(afv3)
CUDA.reclaim()
CUDA.memory_status()

GC.gc(true)
CUDA.reclaim()

cputflops=tflop/cput

gputflops=tflop/gput

max_stars_repo_name: Yangyang-Tan/FRG-Code-Collection

max_issues_repo_name: Yangyang-Tan/FRG-Code-Collection

max_forks_repo_name: Yangyang-Tan/FRG-Code-Collection


text:
using brainflow

# specify logging library to use
brainflow.enable_dev_brainflow_logger(brainflow.BOARD_CONTROLLER)

params = BrainFlowInputParams()
board_shim = brainflow.BoardShim(brainflow.SYNTHETIC_BOARD, params)
sampling_rate = brainflow.get_sampling_rate(brainflow.SYNTHETIC_BOARD)
nfft = brainflow.get_nearest_power_of_two(sampling_rate)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_board_data(board_shim)
brainflow.release_session(board_shim)

eeg_channels = brainflow.get_eeg_channels(brainflow.SYNTHETIC_BOARD)

bands = brainflow.get_avg_band_powers(data, eeg_channels, sampling_rate, true)
println(bands[1])
println(bands[2])

max_stars_repo_name: matthijscox/brainflow

max_issues_repo_name: matthijscox/brainflow

max_forks_repo_name: matthijscox/brainflow


text:
export particlefilter

function particlefilter(hmm::HMM, observations::Matrix{Float}, N::Int,
                        proposal::Proposal;
                        resampling::Function=multinomialresampling,
                        essthresh::Float=0.5
                        )::Tuple{ParticleSet,Vector{Float}}

    K   = size(observations, 2)
    # particle set filter (storage)
    psf = ParticleSet(N, hmm.dimx, K)
    ess = zeros(K)

    (p1,e1) = resample( Particles(
                            [proposal.mu0 + proposal.noise() for i in 1:N],
                            ones(N)/N),
                        essthresh )
    # store
    psf.p[1] = p1
    ess[1]   = e1

    for k=2:K
        pkm1 = psf.p[k-1]
        obsk = observations[:,k]

        logak = zeros(N)
        xk    = similar(pkm1.x)
        # sample (BOOTSTRAP)
        for i in 1:N
            xk[i]    = proposal.mean(k, pkm1.x[i], obsk) + proposal.noise()
            logak[i] = hmm.transloglik(k, pkm1.x[i], xk[i]) +
                        hmm.obsloglik(k, obsk, xk[i]) -
                        proposal.loglik(k, pkm1.x[i], obsk, xk[i])
        end
        Wk  = log.(pkm1.w) + logak
        Wk -= minimum(Wk) # try to avoid underflows
        wk  = exp.(Wk)
        wk /= sum(wk)

        (pk, ek) = resample(Particles(xk,wk), essthresh, resampling)

        psf.p[k] = pk
        ess[k]   = ek
    end
    (psf, ess)
end


text:
struct WebsocketFrame
    config::NamedTuple
    buffers::NamedTuple
    inf::Dict{Symbol, Any}
    function WebsocketFrame(
        config::NamedTuple,
        buffers::NamedTuple,
        binaryPayload::Array{UInt8, 1} = Array{UInt8, 1}();
            ptr::Int = 1,
    )
        inf = Dict{Symbol, Any}(
            :fin => false,
            :mask => false,
            :opcode => 0x00,
            :rsv1 => 0x00,
            :rsv2 => 0x00,
            :rsv3 => 0x00,
            :length => 0x00,
            :parseState => DECODE_HEADER,
            :binaryPayload => binaryPayload,
            :ptr => ptr,
        )
        self = new(
            config,
            buffers,
            inf,
        )
    end
end

function toBuffer(self::WebsocketFrame)

    outBuffer = self.buffers.outBuffer
    !isopen(outBuffer) && return false
    truncate(outBuffer, 0)
    maskBytes = self.buffers.maskBytes
    inf = self.inf
    firstByte = 0x00
    secondByte = 0x00
    inf[:fin] && (firstByte |= 0x80)
    inf[:mask] && (secondByte |= 0x80)
    firstByte |= (inf[:opcode] & 0x0F)

    if inf[:opcode] === CONNECTION_CLOSE_FRAME
        closeStatus = reinterpret(UInt8, [hton(UInt16(inf[:closeStatus]))])
        pushfirst!(inf[:binaryPayload], closeStatus...)
    end

    len = length(inf[:binaryPayload])

    if len <= 125
        secondByte |= (len & 0x7F)
    elseif len > 125 && len <= 0xFFFF
        secondByte |= 126
    elseif len > 0xFFFF
        secondByte |= 127
    end

    write(outBuffer, UInt8(firstByte), UInt8(secondByte))

    if len > 125 && len <= 0xFFFF
        write(outBuffer, hton(UInt16(len)))
    elseif len > 0xFFFF
        write(outBuffer, hton(0x00000000), hton(UInt32(len)))
    end

    if inf[:mask]
        seekstart(maskBytes)
        write(maskBytes, hton(rand(UInt32)))
        seekstart(maskBytes)
        write(outBuffer, maskBytes)
        mask!(maskBytes, inf[:binaryPayload])
    end

    unsafe_write(outBuffer, pointer(inf[:binaryPayload]), len)
    seek(outBuffer, 0)
    return true
end

function addData(self::WebsocketFrame)

    inf = self.inf
    header = self.buffers.frameHeader
    inBuffer = self.buffers.inBuffer
    maskBytes = self.buffers.maskBytes
    seek(inBuffer, inf[:ptr] - 1)

    if inf[:parseState] === DECODE_HEADER && inBuffer.size >= 2
        seekstart(header)
        write(header, read(inBuffer, 2))
        inf[:ptr] = inBuffer.ptr
        firstByte = header.data[1]
        secondByte = header.data[2]
        inf[:fin] = firstByte & WS_FINAL > 0
        inf[:mask] = secondByte & WS_MASK > 0
        inf[:opcode] = firstByte & WS_OPCODE
        inf[:length] = secondByte & WS_LENGTH
        inf[:rsv1] = firstByte & WS_RSV1 > 0
        inf[:rsv2] = firstByte & WS_RSV2 > 0
        inf[:rsv3] = firstByte & WS_RSV3 > 0

        if inf[:rsv1] || inf[:rsv2] || inf[:rsv2]
            throw(error("[$CLOSE_REASON_POLICY_VIOLATION]websocket extensions are not supported"))
        end

        if inf[:opcode] >= 0x08
            inf[:length] > 125 && throw(error("illegal control frame longer than 125 bytes."))
            !inf[:fin] && throw(error("control frames must not be fragmented."))
        end

        if inf[:length] === 0x7e
            inf[:parseState] = WAITING_FOR_16_BIT_LENGTH
        elseif inf[:length] === 0x7f
            inf[:parseState] = WAITING_FOR_64_BIT_LENGTH
        else
            inf[:parseState] = WAITING_FOR_MASK_KEY
        end

    end
    if inf[:parseState] === WAITING_FOR_16_BIT_LENGTH
        if inBuffer.size >= 2
            write(header, read(inBuffer, 2))
            inf[:ptr] = inBuffer.ptr
            ptr = header.ptr
            seek(header, 2)
            inf[:length] = Int(ntoh(read(header, UInt16)))
            seek(header, ptr - 1)
            inf[:parseState] = WAITING_FOR_MASK_KEY
        end
    elseif inf[:parseState] === WAITING_FOR_64_BIT_LENGTH
        if inBuffer.size >= 8
            write(header, read(inBuffer, 8))
            inf[:ptr] = inBuffer.ptr
            ptr = header.ptr
            seek(header, 2)
            lengthPair = [
                Int(ntoh(read(header, UInt32))),
                Int(ntoh(read(header, UInt32)))
            ]
            if lengthPair[1] !== 0
                throw(error("Unsupported 64-bit length frame received"))
            end
            inf[:length] = lengthPair[2]
            inf[:parseState] = WAITING_FOR_MASK_KEY
        end
    end
    if inf[:parseState] === WAITING_FOR_MASK_KEY
        if inf[:mask]
            if inBuffer.size >= 4
                seekstart(maskBytes)
                write(maskBytes, read(inBuffer, 4))
                inf[:ptr] = inBuffer.ptr
                inf[:parseState] = WAITING_FOR_PAYLOAD
            end
        else
            inf[:parseState] = WAITING_FOR_PAYLOAD
        end
    end
    if inf[:parseState] === WAITING_FOR_PAYLOAD

        if inf[:length] > self.config.maxReceivedFrameSize
            throw(error("[$CLOSE_REASON_POLICY_VIOLATION]frame size exceeds maximum of $(self.config.maxReceivedFrameSize) Bytes."))
        end
        if inBuffer.size - (inf[:ptr] - 1) >= inf[:length]
            inf[:binaryPayload] = read(inBuffer, inf[:length])
            inf[:ptr] = inBuffer.ptr
            inf[:mask] && mask!(maskBytes, inf[:binaryPayload])

            if inf[:opcode] === CONNECTION_CLOSE_FRAME
                inf[:length] < 2 && throw(error("close frame too small"))
                status = IOBuffer(splice!(inf[:binaryPayload], 1:2))
                inf[:closeStatus] = Int(ntoh(read(status, UInt16)))
                close(status)
            end

            inf[:parseState] = COMPLETE

            return true
        end
    end
    return false
end

max_stars_repo_name: citkane/SimpleWebsockets

max_issues_repo_name: citkane/SimpleWebsockets

max_forks_repo_name: citkane/SimpleWebsockets


text:
using CSV
using DataFrames

const DIR = @__DIR__
const ROOT_DIR = joinpath(@__DIR__, "..")

res_file_all  = ARGS[1]
res_file_ooqp = ARGS[2]

# Read each results file
df_all = CSV.read(res_file_all)    # Non-OOQP solvers
df_ooqp = CSV.read(res_file_ooqp)  # OOQP results

# Merge results
# This will raise an error if `df_all` already containts `time_OOQP` and/or `success_OOQP`
df_merged = join(df_all, df_ooqp, on = :instance)

# Write results
CSV.write("res_merged.csv", df_merged)


text:
module Dojo

# constants
global REG = 1.0e-10

#TODO: remove
using FiniteDiff

using LinearAlgebra
using Random
using StaticArrays
using SparseArrays
using StaticArrays: SUnitRange
using Quaternions
using Parameters
using Statistics

using Colors
using Colors: RGBA, RGB
using FFMPEG
using LightXML
using MeshCat
import MeshCat: render
using Meshing
using GeometryBasics
using Graphs

using CoordinateTransformations

using JLD2
using DocStringExtensions

# Utilities
include(joinpath("utilities", "methods.jl"))
include(joinpath("utilities", "custom_static.jl"))
include(joinpath("utilities", "normalize.jl"))

# Orientation
include(joinpath("orientation", "quaternion.jl"))
include(joinpath("orientation", "mrp.jl"))
include(joinpath("orientation", "axis_angle.jl"))
include(joinpath("orientation", "mapping.jl"))
include(joinpath("orientation", "rotate.jl"))

# Graph system
include(joinpath("graph", "entry.jl"))
include(joinpath("graph", "system.jl"))
include(joinpath("graph", "linear_system.jl"))
include(joinpath("graph", "adjacency.jl"))
include(joinpath("graph", "depth_first_search.jl"))
include(joinpath("graph", "cycles.jl"))
include(joinpath("graph", "ldu_factorization.jl"))

# Graph objects
include(joinpath("mechanism", "node.jl"))
include(joinpath("mechanism", "edge.jl"))
include(joinpath("mechanism", "id.jl"))

# Bodies
include(joinpath("bodies", "shapes.jl"))
include(joinpath("bodies", "state.jl"))
include(joinpath("bodies", "constructor.jl"))
include(joinpath("bodies", "origin.jl"))
include(joinpath("bodies", "set.jl"))

# Mechanism
include(joinpath("joints", "constraints.jl"))
include(joinpath("contacts", "constructor.jl"))
include(joinpath("contacts", "contact.jl"))

include(joinpath("mechanism", "constructor.jl"))
include(joinpath("mechanism", "gravity.jl"))
include(joinpath("mechanism", "state.jl"))
include(joinpath("mechanism", "system.jl"))
include(joinpath("mechanism", "methods.jl"))
include(joinpath("mechanism", "set.jl"))
include(joinpath("mechanism", "get.jl"))
include(joinpath("mechanism", "urdf.jl"))
include(joinpath("mechanism", "traversal.jl"))

# Simulation
include(joinpath("simulation", "step.jl"))
include(joinpath("simulation", "storage.jl"))
include(joinpath("simulation", "simulate.jl"))

# Mechanics
include(joinpath("mechanics", "momentum.jl"))
include(joinpath("mechanics", "energy.jl"))

# Joints
include(joinpath("joints", "orthogonal.jl"))
include(joinpath("joints", "joint.jl"))
include(joinpath("joints", "translational", "constructor.jl"))
include(joinpath("joints", "translational", "impulses.jl"))
include(joinpath("joints", "translational", "input.jl"))
include(joinpath("joints", "translational", "springs.jl"))
include(joinpath("joints", "translational", "dampers.jl"))
include(joinpath("joints", "translational", "minimal.jl"))
include(joinpath("joints", "rotational", "constructor.jl"))
include(joinpath("joints", "rotational", "impulses.jl"))
include(joinpath("joints", "rotational", "input.jl"))
include(joinpath("joints", "rotational", "springs.jl"))
include(joinpath("joints", "rotational", "dampers.jl"))
include(joinpath("joints", "rotational", "minimal.jl"))
include(joinpath("joints", "limits.jl"))
include(joinpath("joints", "prototypes.jl"))
include(joinpath("joints", "minimal.jl"))
include(joinpath("joints", "impulses.jl"))

# Contacts
include(joinpath("contacts", "constraints.jl"))
include(joinpath("contacts", "cone.jl"))
include(joinpath("contacts", "collisions", "collision.jl"))
include(joinpath("contacts", "collisions", "sphere_halfspace.jl"))
include(joinpath("contacts", "collisions", "sphere_sphere.jl"))
include(joinpath("contacts", "velocity.jl"))
include(joinpath("contacts", "impact.jl"))
include(joinpath("contacts", "linear.jl"))
include(joinpath("contacts", "nonlinear.jl"))
include(joinpath("contacts", "utilities.jl"))

# Solver
include(joinpath("solver", "linear_system.jl"))
include(joinpath("solver", "centering.jl"))
include(joinpath("solver", "complementarity.jl"))
include(joinpath("solver", "violations.jl"))
include(joinpath("solver", "options.jl"))
include(joinpath("solver", "initialization.jl"))
include(joinpath("solver", "correction.jl"))
include(joinpath("solver", "mehrotra.jl"))
include(joinpath("solver", "line_search.jl"))

# Integrator
include(joinpath("integrators", "integrator.jl"))
include(joinpath("integrators", "constraint.jl"))

# Visualizer
include(joinpath("visuals", "visualizer.jl"))
include(joinpath("visuals", "set.jl"))
include(joinpath("visuals", "convert.jl"))
include(joinpath("visuals", "colors.jl"))

# Data
include(joinpath("mechanism", "data.jl"))

# Gradients
include(joinpath("gradients", "finite_difference.jl"))
include(joinpath("gradients", "state.jl"))
include(joinpath("gradients", "data.jl"))
include(joinpath("gradients", "utilities.jl"))

# Environments
include(joinpath("..", "environments", "mechanisms.jl"))
include(joinpath("..", "environments", "environment.jl"))
include(joinpath("..", "environments", "dynamics.jl"))
include(joinpath("..", "environments", "utilities.jl"))
include(joinpath("..", "environments", "include.jl"))

# Bodies
export
    Body,
    Origin,
    Box,
    Cylinder,
    Sphere,
    Pyramid,
    Mesh,
    Shapes,
    get_body,
    get_node

# Joints
export
    Rotational,
    Translational,
    JointConstraint,
    Floating,
    Fixed,
    Prismatic,
    Planar,
    FixedOrientation,
    Revolute,
    Cylindrical,
    PlanarAxis,
    FreeRevolute,
    Orbital,
    PrismaticOrbital,
    PlanarOrbital,
    FreeOrbital,
    Spherical,
    CylindricalFree,
    PlanarFree,
    get_joint,
    nullspace_mask

# Contacts
export
    ContactConstraint,
    ImpactContact,
    LinearContact,
    NonlinearContact,
    get_contact,
    get_sdf,
    contact_location,
    damper_impulses,
    contact_constraint

# Collision
export
    SphereHalfSpaceCollision,
    SphereSphereCollision,
    distance,
    contact_point,
    contact_normal,
    contact_tangent

# Inputs
export
    set_input!,
    add_input!,
    input_dimension

# Mechanism
export
    Mechanism,
    get_mechanism,
    initialize!,
    set_floating_base,
    zero_velocity!

# Maximal
export
    set_maximal_state!,
    set_maximal_configurations!,
    set_maximal_velocities!,
    get_maximal_state,
    maximal_dimension

# Minimal
export
    set_minimal_state!,
    set_minimal_coordinates!,
    set_minimal_velocities!,
    get_minimal_state,
    minimal_coordinates,
    minimal_velocities,
    minimal_dimension

# Maximal <-> Minimal
export
    maximal_to_minimal,
    minimal_to_maximal

# Simulation
export
    simulate!,
    step!,
    generate_storage

# Environment
export
    Environment,
    dynamics,
    dynamics_jacobian_state,
    dynamics_jacobian_input,
    get_environment,
    #step,
    get_observation,
    cost,
    is_done,
    #reset,
    render,
    seed,
    close,
    Space,
    BoxSpace

# Orientation
export
    Quaternion,
    attitude_jacobian

# Data
export
    get_data,
    set_data!,
    get_solution

# Gradients
export
    maximal_to_minimal_jacobian,
    minimal_to_maximal_jacobian

# Mechanics
export
    kinetic_energy,
    potential_energy,
    mechanical_energy,
    momentum

# Solver
export
    mehrotra!,
    SolverOptions

# Linear System "Ax = b"
export
    full_matrix,
    full_data_matrix

# Visuals
export
    Visualizer,
    visualize,
    set_background!,
    set_floor!,
    set_surface!,
    set_light!,
    set_camera!,
    RGBA,
    orange,
    cyan,
    magenta

# Static
export
    szeros,
    sones,
    srand

# Utilities
export
    normalize

end


text:
module TestModifying

using DataTools
using Accessors: @optic
using Test
using Transducers

reduce_bs1(args...; kw...) = reduce(args...; basesize = 1, kw...)

@testset for fold in [foldl, reduce_bs1, reduce]
    @test fold(modifying(a = +), Map(identity), [(a = 1, b = 2), (a = 3, b = 4)]) ==
          (a = 4, b = 2)
end

@testset "map" begin
    @test map(modifying(a = string), [(a = 1, b = 2), (a = 3, b = 4)]) ==
          [(a = "1", b = 2), (a = "3", b = 4)]
    @test map(
        modifying(@optic(_.a[1].b) => x -> 10x),
        [(a = ((b = 1,), 2),), (a = ((b = 3,), 4),)],
    ) == [(a = ((b = 10,), 2),), (a = ((b = 30,), 4),)]
end

end  # module


text:
"""
    ladderize_tree!(root::T, ascending::Bool=true)
        ::Nothing where T<:AbstractNode

This function ladderizes a tree inplace, i.e. sorts the nodes on all levels by the count
of their descendants.

* `root` : root Node of tree.

* `ascending` : Boolean, determines whether to sort in ascending (true) or
                descending (false) order.
"""
function ladderize_tree!(root::T, ascending::Bool=true) where T<:AbstractNode
    root.nchild == 0 && return nothing
    ndescendants = Array{Float64,1}(undef, length(root.children))
    for (index, child) in enumerate(root.children)
        ndescendants[index] = size(child)[1]
    end
    perm = sortperm(ndescendants, rev=!ascending)
    root.children = root.children[perm]

    for child in root.children
        ladderize_tree!(child, ascending)
    end
end

"""
    ladderize_tree(root::T, ascending::Bool=true)::T where T<:AbstractNode

This function returns a ladderized copy of a tree, i.e. a copy with all the
nodes on all levels sorted by the count of their descendants.

* `root` : root Node of tree.

* `ascending` : Boolean, determines whether to sort in ascending (true) or
                descending (false) order.
"""
function ladderize_tree(root::T, ascending::Bool=true)::T where T<:AbstractNode
    copyroot = deepcopy(root)
    ladderize_tree!(copyroot, ascending)
    return copyroot
end


text:
const textDocument_codeAction_request_type = JSONRPC.RequestType("textDocument/codeAction", CodeActionParams, Vector{Command})
const textDocument_completion_request_type = JSONRPC.RequestType("textDocument/completion", CompletionParams, CompletionList)
const textDocument_signatureHelp_request_type = JSONRPC.RequestType("textDocument/signatureHelp", TextDocumentPositionParams, SignatureHelp)
const textDocument_definition_request_type = JSONRPC.RequestType("textDocument/definition", TextDocumentPositionParams, Vector{Location})
const textDocument_formatting_request_type = JSONRPC.RequestType("textDocument/formatting", DocumentFormattingParams, Vector{TextEdit})
const textDocument_references_request_type = JSONRPC.RequestType("textDocument/references", ReferenceParams, Vector{Location})
const textDocument_rename_request_type = JSONRPC.RequestType("textDocument/rename", RenameParams, WorkspaceEdit)
const textDocument_documentSymbol_request_type = JSONRPC.RequestType("textDocument/documentSymbol", DocumentSymbolParams, Vector{SymbolInformation})
const textDocument_documentHighlight_request_type = JSONRPC.RequestType("textDocument/documentHighlight", DocumentHighlightParams, Vector{DocumentHighlight})
const textDocument_hover_request_type = JSONRPC.RequestType("textDocument/hover", TextDocumentPositionParams, Hover)
const textDocument_didOpen_notification_type = JSONRPC.NotificationType("textDocument/didOpen", DidOpenTextDocumentParams)
const textDocument_didClose_notification_type = JSONRPC.NotificationType("textDocument/didClose", DidCloseTextDocumentParams)
const textDocument_didChange_notification_type = JSONRPC.NotificationType("textDocument/didChange", DidChangeTextDocumentParams)
const textDocument_didSave_notification_type = JSONRPC.NotificationType("textDocument/didSave", DidSaveTextDocumentParams)
const textDocument_willSave_notification_type = JSONRPC.NotificationType("textDocument/willSave", WillSaveTextDocumentParams)
const textDocument_willSaveWaitUntil_request_type = JSONRPC.RequestType("textDocument/willSaveWaitUntil", WillSaveTextDocumentParams, Vector{TextEdit})
const textDocument_publishDiagnostics_notification_type = JSONRPC.NotificationType("textDocument/publishDiagnostics", PublishDiagnosticsParams)
const textDocument_selectionRange_request_type = JSONRPC.RequestType("textDocument/selectionRange", SelectionRangeParams, Vector{SelectionRange})

const workspace_executeCommand_request_type = JSONRPC.RequestType("workspace/executeCommand", ExecuteCommandParams, Nothing)
const workspace_symbol_request_type = JSONRPC.RequestType("workspace/symbol", WorkspaceSymbolParams, Vector{SymbolInformation})
const workspace_didChangeWatchedFiles_notification_type = JSONRPC.NotificationType("workspace/didChangeWatchedFiles", DidChangeWatchedFilesParams)
const workspace_didChangeConfiguration_notification_type = JSONRPC.NotificationType("workspace/didChangeConfiguration", DidChangeConfigurationParams)
const workspace_didChangeWorkspaceFolders_notification_type = JSONRPC.NotificationType("workspace/didChangeWorkspaceFolders", DidChangeWorkspaceFoldersParams)
const workspace_applyEdit_request_type = JSONRPC.RequestType("workspace/applyEdit", ApplyWorkspaceEditParams, ApplyWorkspaceEditResponse)
const workspace_configuration_request_type = JSONRPC.RequestType("workspace/configuration", ConfigurationParams, Vector{Any})
const julia_activateenvironment_notification_type = JSONRPC.NotificationType("julia/activateenvironment", NamedTuple{(:envPath,),Tuple{String}})
const julia_refreshLanguageServer_notification_type = JSONRPC.NotificationType("julia/refreshLanguageServer", Nothing)

const initialize_request_type = JSONRPC.RequestType("initialize", InitializeParams, InitializeResult)
const initialized_notification_type = JSONRPC.NotificationType("initialized", InitializedParams)
const shutdown_request_type = JSONRPC.RequestType("shutdown", Nothing, Nothing)
const exit_notification_type = JSONRPC.NotificationType("exit", Nothing)
const client_registerCapability_request_type = JSONRPC.RequestType("client/registerCapability", RegistrationParams, Nothing)

const cancel_notification_type = JSONRPC.NotificationType("\$/cancelRequest", CancelParams)
const setTrace_notification_type = JSONRPC.NotificationType("\$/setTrace", SetTraceParams)
# TODO This seems to not exist in the spec?
const setTraceNotification_notification_type = JSONRPC.NotificationType("\$/setTraceNotification", Nothing)

const window_workDoneProgress_create_request_type = JSONRPC.RequestType("window/workDoneProgress/create", WorkDoneProgressCreateParams, Nothing)
const window_showMessage_notification_type = JSONRPC.NotificationType("window/showMessage", ShowMessageRequestParams)

const progress_notification_type = JSONRPC.NotificationType("\$/progress", ProgressParams)

const telemetry_event_notification_type = JSONRPC.NotificationType("telemetry/event", Any)


text:
using Test, bad; import Distributions

pnull = .25
palt  = .4
p     = Beta(5, 7)
α, β  = .05, .2
problem = Problem(
        minimise( SampleSize(p) ),
        Power(p | pnull) <= α,
        Power(p | palt)  >= 1 - β,
        nmax  = 100,
        n1max = 40
    )
design = optimise(problem; verbosity = 0)

samplesize = SampleSize(p)

@test n1(design) == samplesize(design, 0, 0, .2)
@test n1(design) ≈ samplesize(design, 0, 0)
@test all( n1(design) .≈ samplesize.(design, early_stop_region(design), .2) )
@test all( n1(design) .≈ samplesize.(design,  early_stop_region(design)) )
@test all( n1(design) .< samplesize.(design,  continuation_region(design)) )
@test expectation( p -> sum(
        Distributions.pdf.(Distributions.Binomial(n1(design), p), 0:n1(design)) .*
        n.(design, 0:n1(design))
    ),
    p
) ≈ samplesize(design)

pow = problem.power.score

@test 0 ≈ pow(design, 0, 0, .2)
@test (0 .≈ pow.(design, futility_region(design), pow.bounds[1]) ) |> all
@test (0 .≈ pow.(design, futility_region(design)) ) |> all
@test (1 .≈ pow.(design, efficacy_region(design), pow.bounds[1]) ) |> all
@test (1 .≈ pow.(design, efficacy_region(design)) ) |> all
@test (0 .≈ pow.(design, continuation_region(design), 0.) ) |> all
@test (0.5 .< pow.(design, continuation_region(design)) .< .99) |> all
@test 0 ≈ pow(design, 6, .2)
@test 0 ≈ pow(design, .2)
@test abs(1 - β - pow(design)) < 1e-3

toer = problem.toer.score

@test 0 ≈ toer(design, 0, 0, .2)
@test (0 .≈ toer.(design, futility_region(design), toer.bounds[1]) ) |> all
@test (0 .≈ toer.(design, futility_region(design)) ) |> all
@test (1 .≈ toer.(design, efficacy_region(design), toer.bounds[1]) ) |> all
@test (1 .≈ toer.(design, efficacy_region(design)) ) |> all
@test (0 .≈ toer.(design, continuation_region(design), 0)) |> all
@test (0.0 .< toer.(design, continuation_region(design)) .< .99) |> all
@test 1 ≈ toer(design, 20, pnull)
@test toer(design, pnull) <= α
@test toer(design, pnull/2) < toer(design, pnull)

# check for monotone conditional error/power (for n1partial > 0)
for n1partial in 1:n1(design)
    @test begin
        map(
            x1p -> pow(design; partial_stage_one = (x1p, n1partial) ),
            0:n1partial
        ) |>
        x -> diff(x) |>
        x -> minimum(x) >= -sqrt(eps()) # numerical inaccuracies
    end
    @test begin
        map(
            x1p -> toer(design; partial_stage_one = (x1p, n1partial) ),
            0:n1partial
        ) |>
        x -> diff(x) |>
        x -> minimum(x) >= -sqrt(eps()) # numerical inaccuracies
    end
end

# check for monotone conditional error/power in stage two
for x1 in 0:n1(design)
    for n2partial in 1:n2(design, x1)
        @test begin
            map(
                x2p -> pow(design, x1; partial_stage_two = (x2p, n2partial) ),
                0:n2partial
            ) |>
            x -> diff(x) |>
            x -> minimum(x) >= -sqrt(eps()) # numerical inaccuracies
        end
        @test begin
            map(
                x2p -> toer(design, x1; partial_stage_two = (x2p, n2partial) ),
                0:n2partial
            ) |>
            x -> diff(x) |>
            x -> minimum(x) >= -sqrt(eps()) # numerical inaccuracies
        end
    end
end

@test true


text:
module Scan
include("token_type.jl")

using Printf
using .TokenType: Token, toPreversedType, toPunctuationType

@enum ScannerState begin
    START
    # PREPROCESS # deprecated
    IDENTIFIER
    CONSTANT
    SLASH_MODE
end

mutable struct Scanner
    lc::Int                             # line count
    state::ScannerState                 # state of scanner
    sign::Union{Char,Nothing}           # sign record token
    specifier::Union{Char,Nothing}      # specifier record token
    __stream::String                    # stream from file
    __index::Int                        # index in stream
    Scanner(stream::String) = new(1, START, nothing, nothing, stream, 0)
end

function next(sc::Scanner, ignoreSpace::Bool = false)::Union{Char,Nothing}
    sc.__index += 1
    result = sc.__index > lastindex(sc.__stream) ? nothing : sc.__stream[sc.__index]
    return (result in [' ', '\t'] && ignoreSpace) ? next(sc, true) : result
end

function retract(sc::Scanner, n::Int = 1)::Bool
    isSuccess = ((sc.__index - n) > 0)
    sc.__index = isSuccess ? sc.__index - n : 1
    return isSuccess
end
end

using .Scan: Scanner

function start!(sc::Scanner, c::Char, tokens::Array)
    buf = string(c)
    # new line
    if c == '\r'
        if Scan.next(sc) != '\n'
            Scan.retract(sc)
        end
        sc.lc += 1
        push!(tokens, (TokenType.LN_TOKEN, nothing))
    elseif c == '\n'
        sc.lc += 1
        push!(tokens, (TokenType.LN_TOKEN, nothing))
    elseif c == '#'
        push!(tokens, (TokenType.SHARP, nothing))
        # sc.state = Scan.PREPROCESS # deprecated
    elseif c == ','
        push!(tokens, (TokenType.COMMA, nothing))
    elseif c == ':'
        push!(tokens, (TokenType.COLON, nothing))
    elseif c == ';'
        push!(tokens, (TokenType.SEMICOLON, nothing))
    elseif c == '('
        push!(tokens, (TokenType.LEFT_PARENTHESE, nothing))
    elseif c == ')'
        push!(tokens, (TokenType.RIGHT_PARENTHESE, nothing))
    elseif c == '['
        push!(tokens, (TokenType.LEFT_BRACKET, nothing))
    elseif c == ']'
        push!(tokens, (TokenType.RIGHT_BRACKET, nothing))
    elseif c == '{'
        push!(tokens, (TokenType.LEFT_BRACE, nothing))
    elseif c == '}'
        push!(tokens, (TokenType.RIGHT_BRACE, nothing))
    elseif c == '"'
        if isempty(tokens)
            throw(UnIdentifierError(buf))
        end
        lastTokenType = nothing
        for type in reverse(tokens)
            if type[1] != TokenType.RIGHT_PARENTHESE
                lastTokenType = type[1]
                break
            end
        end
        if lastTokenType == TokenType.INCLUDE
            c = Scan.next(sc)
            while c != '\r' && c != '\n' && c !== nothing
                buf *= c
                if c == '>'
                    push!(tokens, (TokenType.LIBRARY_NAME, buf))
                    buf = ""
                    break
                end
                c = Scan.next(sc)
            end
            if !isempty(buf)
                throw(UnIdentifierError(buf))
            end
        else
            push!(tokens, __string_seq_stream!(sc, c, '"', true))
        end
    elseif c == '\''
        push!(tokens, __string_seq_stream!(sc, c, '"', true))
    elseif c == '+'
        if isempty(tokens)
            throw(UnIdentifierError(buf))
        end
        lastTokenType = nothing
        for type in reverse(tokens)
            if type[1] != TokenType.RIGHT_PARENTHESE
                lastTokenType = type[1]
                break
            end
        end
        # operator id(++|--) +(=|+)
        if lastTokenType in [TokenType.POINTER_IDENTIFIER, TokenType.IDENTIFIER, TokenType.DPLUS, TokenType.DMINUS]
            c = Scan.next(sc)
            # +=
            if c == '='
                push!(tokens, (TokenType.ACC_PLUS, nothing))
            # ++
            elseif c == '+'
                push!(tokens, (TokenType.DPLUS, nothing))
            # only + need retract
            else
                push!(tokens, (TokenType.PLUS, nothing))
                Scan.retract(sc)
            end
        # operator constant +
        elseif lastTokenType in [TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON]
            push!(tokens, (TokenType.PLUS, nothing))
        # positive-sign (+123), += +123, -= +123, + +2, - +2, = +2
        # elseif lastTokenType in [TokenType.LEFT_PARENTHESE, TokenType.ACC_PLUS, TokenType.ACC_MINUS, TokenType.PLUS, TokenType.MINUS, TokenType.ASSIGN]
        #     sc.sign = c
        elseif lastTokenType !== nothing
            sc.sign = c
        else
            throw(UnIdentifierError(buf))
        end
    elseif c == '-'
        if isempty(tokens)
            throw(UnIdentifierError(buf))
        end
        lastTokenType = nothing
        for type in reverse(tokens)
            if type[1] != TokenType.RIGHT_PARENTHESE
                lastTokenType = type[1]
                break
            end
        end
        # operator id(++|--) -(=|-)
        if lastTokenType in [TokenType.POINTER_IDENTIFIER, TokenType.IDENTIFIER, TokenType.DPLUS, TokenType.DMINUS]
            c = Scan.next(sc)
            # -=
            if c == '='
                push!(tokens, (TokenType.ACC_MINUS, nothing))
            # --
            elseif c == '-'
                push!(tokens, (TokenType.DMINUS, nothing))
            # only - need retract
            else
                push!(tokens, (TokenType.MINUS, nothing))
                Scan.retract(sc)
            end
        # operator constant -
        elseif lastTokenType in [TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON]
            push!(tokens, (TokenType.MINUS, nothing))
        # negative-sign (-123), += -123, -= -123, + -2, - -2, = -2
        # elseif lastTokenType in [TokenType.LEFT_PARENTHESE, TokenType.ACC_PLUS, TokenType.ACC_MINUS, TokenType.PLUS, TokenType.MINUS, TokenType.ASSIGN, TokenType.]
        #     sc.sign = c
        elseif lastTokenType !== nothing
            sc.sign = c
        else
            throw(UnIdentifierError(buf))
        end
    elseif c == '*'
        if isempty(tokens)
            throw(UnIdentifierError(buf))
        end
        lastTokenType = nothing
        for type in reverse(tokens)
            if type[1] != TokenType.RIGHT_PARENTHESE
                lastTokenType = type[1]
                break
            end
        end
        # operator id(++)/(--) *(=)
        if lastTokenType in [TokenType.POINTER_IDENTIFIER, TokenType.IDENTIFIER, TokenType.DPLUS, TokenType.DMINUS]
            c = Scan.next(sc)
            # *=
            if c == '='
                push!(tokens, (TokenType.ACC_MULTIPLY, nothing))
            # only * need retract
            else
                push!(tokens, (TokenType.MULTIPLY, nothing))
                Scan.retract(sc)
            end
        # operator constant *
        elseif lastTokenType in [TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON]
            push!(tokens, (TokenType.MULTIPLY, nothing))
        # if target is pointer (allow format *id with space between)
        else
            sc.specifier = c
            c = Scan.next(sc, true)
            if c == '\r' || c == '\n'
                c = nothing
                Scan.retract(sc)
            end
            # if pointer name is valid
            if c !== nothing && (c >= 'a' && c <= 'z' || c >= 'A' && c <= 'Z' || c == '_')
                sc.state = Scan.IDENTIFIER
                Scan.retract(sc)
            else
                throw(UnIdentifierError((sc.sign === nothing ? "" : sc.sign) * buf * string(c === nothing ? "" : c)))
            end
        end
    elseif c == '/'
        sc.state = Scan.SLASH_MODE
    elseif c == '%'
        if isempty(tokens)
            throw(UnIdentifierError(buf))
        end
        lastTokenType = nothing
        for type in reverse(tokens)
            if type[1] != TokenType.RIGHT_PARENTHESE
                lastTokenType = type[1]
                break
            end
        end
        # operator id(++)/(--) %(=)
        if lastTokenType in [TokenType.POINTER_IDENTIFIER, TokenType.IDENTIFIER, TokenType.DPLUS, TokenType.DMINUS]
            c = Scan.next(sc)
            # %=
            if c == '='
                push!(tokens, (TokenType.ACC_MODULO, nothing))
            # only % need retract
            else
                push!(tokens, (TokenType.MODULO, nothing))
                Scan.retract(sc)
            end
        # operator constant %
        elseif lastTokenType in [TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON]
            push!(tokens, (TokenType.MODULO, nothing))
        else
            throw(UnIdentifierError(buf))
        end
    elseif c == '='
        if isempty(tokens)
            throw(UnIdentifierError(buf))
        end
        c = Scan.next(sc)
        lastTokenType = nothing
        for type in reverse(tokens)
            if type[1] != TokenType.RIGHT_PARENTHESE
                lastTokenType = type[1]
                break
            end
        end
        # compare
        if c == '=' && lastTokenType in [TokenType.POINTER_IDENTIFIER, TokenType.IDENTIFIER, TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON, TokenType.DPLUS, TokenType.DMINUS]
            c = Scan.next(sc)
            # === identifier equal
            if c == '='
                push!(tokens, (TokenType.IDENTIFIER_EQUAL, nothing))
            # == equal
            else
                push!(tokens, (TokenType.EQUAL, nothing))
                Scan.retract(sc)
            end
        # = assign
        elseif lastTokenType in [TokenType.POINTER_IDENTIFIER, TokenType.IDENTIFIER, TokenType.DPLUS, TokenType.DMINUS]
            push!(tokens, (TokenType.ASSIGN, nothing))
            Scan.retract(sc)
        else
            throw(UnIdentifierError(buf * c))
        end
    elseif c == '!'
        c = Scan.next(sc)
        # != compare
        if c == '=' && !isempty(tokens)
            if last(tokens)[1] in [TokenType.IDENTIFIER, TokenType.POINTER_IDENTIFIER, TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON, TokenType.DPLUS, TokenType.DMINUS]
                push!(tokens, (TokenType.NOT_EQUAL, nothing))
            else
                throw(UnIdentifierError(buf * c))
            end
        # ! not
        else
            push!(tokens, (TokenType.NOT, nothing))
            Scan.retract(sc)
        end
    elseif c == '&'
        if isempty(tokens)
            throw(UnIdentifierError(buf))
        end
        lastTokenType = nothing
        for type in reverse(tokens)
            if type[1] != TokenType.RIGHT_PARENTHESE
                lastTokenType = type[1]
                break
            end
        end
        # operator id(++)/(--) &(=)
        if lastTokenType in [TokenType.POINTER_IDENTIFIER, TokenType.IDENTIFIER, TokenType.DPLUS, TokenType.DMINUS]
            c = Scan.next(sc)
            # &=
            if c == '='
                push!(tokens, (TokenType.ACC_BIT_AND, nothing))
            # &&
            elseif c == '&'
                push!(tokens, (TokenType.AND, nothing))
            # only & need retract
            else
                push!(tokens, (TokenType.BIT_AND, nothing))
                Scan.retract(sc)
            end
        # operator constant &
        elseif lastTokenType in [TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON]
            push!(tokens, (TokenType.BIT_AND, nothing))
        # if target is address (allow format &id with space between)
        else
            sc.specifier = c
            c = Scan.next(sc, true)
            if c == '\r' || c == '\n'
                c = nothing
                Scan.retract(sc)
            end
            # if address name is valid
            if c !== nothing && (c >= 'a' && c <= 'z' || c >= 'A' && c <= 'Z' || c == '_')
                sc.state = Scan.IDENTIFIER
                Scan.retract(sc)
            else
                throw(UnIdentifierError((sc.sign === nothing ? "" : sc.sign) * buf * (c === nothing ? "" : c)))
            end
        end
    elseif c == '|'
        if isempty(tokens)
            throw(UnIdentifierError(buf))
        end
        lastTokenType = nothing
        for type in reverse(tokens)
            if type[1] != TokenType.RIGHT_PARENTHESE
                lastTokenType = type[1]
                break
            end
        end
        # operator id(++)/(--) |(=)
        if lastTokenType in [TokenType.IDENTIFIER, TokenType.POINTER_IDENTIFIER, TokenType.DPLUS, TokenType.DMINUS]
            c = Scan.next(sc)
            # |=
            if c == '='
                push!(tokens, (TokenType.ACC_BIT_OR, nothing))
            # ||
            elseif c == '|'
                push!(tokens, (TokenType.OR, nothing))
            # only | need retract
            else
                push!(tokens, (TokenType.BIT_OR, nothing))
                Scan.retract(sc)
            end
        # operator constant |
        elseif lastTokenType in [TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON]
            push!(tokens, (TokenType.BIT_OR, nothing))
        else
            throw(UnIdentifierError(buf))
        end
    elseif c == '^'
        if isempty(tokens)
            throw(UnIdentifierError(buf))
        end
        lastTokenType = nothing
        for type in reverse(tokens)
            if type[1] != TokenType.RIGHT_PARENTHESE
                lastTokenType = type[1]
                break
            end
        end
        # operator id(++)/(--) ^(=)
        if lastTokenType in [TokenType.IDENTIFIER, TokenType.POINTER_IDENTIFIER, TokenType.DPLUS, TokenType.DMINUS]
            c = Scan.next(sc)
            # ^=
            if c == '='
                push!(tokens, (TokenType.ACC_BIT_XOR, nothing))
            # only ^ need retract
            else
                push!(tokens, (TokenType.BIT_XOR, nothing))
                Scan.retract(sc)
            end
        # operator constant ^
        elseif lastTokenType in [TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON]
            push!(tokens, (TokenType.BIT_XOR, nothing))
        else
            throw(UnIdentifierError(buf))
        end
    elseif c == '~'
        push!(tokens, (TokenType.BIT_NOT, nothing))
    elseif c == '<'
        if isempty(tokens)
            c = Scan.next(sc)
            while c != '\r' && c != '\n' && c !== nothing
                buf *= c
                if c == '>'
                    push!(tokens, (TokenType.LIBRARY_NAME, buf))
                    buf = ""
                    break
                end
                c = Scan.next(sc)
            end
            if !isempty(buf)
                throw(UnIdentifierError(buf))
            end
        else
            lastTokenType = nothing
            for type in reverse(tokens)
                if type[1] != TokenType.RIGHT_PARENTHESE
                    lastTokenType = type[1]
                    break
                end
            end
            if lastTokenType in [TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON, TokenType.IDENTIFIER, TokenType.POINTER_IDENTIFIER, TokenType.DPLUS, TokenType.DMINUS]
                c = Scan.next(sc)
                buf *= c
                if c == '<'
                    c = Scan.next(sc)
                    # <<=
                    if c == '='
                        if !(lastTokenType in [TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON])
                            push!(tokens, (TokenType.ACC_LEFT_SHIFT, nothing))
                        else
                            throw(UnIdentifierError(buf * c))
                        end
                    # << without = need retract
                    else
                        push!(tokens, (TokenType.LEFT_SHIFT, nothing))
                        Scan.retract(sc)
                    end
                # <=
                elseif c == '='
                    push!(tokens, (TokenType.LESS_EQUAL, nothing))
                # only < need retract
                else
                    push!(tokens, (TokenType.LESS, nothing))
                    Scan.retract(sc)
                end
            elseif lastTokenType == TokenType.INCLUDE
                c = Scan.next(sc)
                while c != '\r' && c != '\n' && c !== nothing
                    buf *= c
                    if c == '>'
                        push!(tokens, (TokenType.LIBRARY_NAME, buf))
                        buf = ""
                        break
                    end
                    c = Scan.next(sc)
                end
                if !isempty(buf)
                    throw(UnIdentifierError(buf))
                end
            else
                throw(UnIdentifierError(buf))
            end
        end
    elseif c == '>'
        if isempty(tokens)
            throw(UnIdentifierError(buf))
        end
        lastTokenType = nothing
        for type in reverse(tokens)
            if type[1] != TokenType.RIGHT_PARENTHESE
                lastTokenType = type[1]
                break
            end
        end
        if lastTokenType in [TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON, TokenType.POINTER_IDENTIFIER, TokenType.IDENTIFIER, TokenType.DPLUS, TokenType.DMINUS]
            c = Scan.next(sc)
            buf *= c
            if c == '>'
                c = Scan.next(sc)
                # >>=
                if c == '='
                    if !(lastTokenType in [TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON])
                        push!(tokens, (TokenType.ACC_RIGHT_SHIFT, nothing))
                    else
                        throw(UnIdentifierError(buf * c))
                    end
                # >> without = need retract
                else
                    push!(tokens, (TokenType.RIGHT_SHIFT, nothing))
                    Scan.retract(sc)
                end
            # >=
            elseif c == '='
                push!(tokens, (TokenType.GREATER_EQUAL, nothing))
            # only > need retract
            else
                push!(tokens, (TokenType.GREATER, nothing))
                Scan.retract(sc)
            end
        else
            throw(UnIdentifierError(buf))
        end
    elseif c >= 'a' && c <= 'z' || c >= 'A' && c <= 'Z' || c == '_'
        sc.state = Scan.IDENTIFIER
        Scan.retract(sc)
    elseif c >= '0' && c <= '9'
        sc.state = Scan.CONSTANT
        Scan.retract(sc)
    elseif c == '.'
        c = Scan.next(sc)
        Scan.retract(sc)
        if c >= '0' && c <= '9'
            sc.state = Scan.CONSTANT
            Scan.retract(sc)
        else
            throw(UnIdentifierError(buf))
        end
    else
        throw(UnIdentifierError(buf))
    end
end

# deprecated
# function preprocess!(sc::Scanner, c::Char, tokens::Array)
#     buffer, keyword, isValid = "", "", true
#     isUserLib = nothing
#     while c !== nothing && c != '\r' && c != '\n'
#         if isempty(keyword)
#             if c >= 'a' && c <= 'z' || c >= 'A' && c <= 'Z'
#                 if !isValid
#                     push!(tokens, (TokenType.UNIDENTIFIER, rstrip(buffer)))
#                     buffer, isValid = "", true
#                 end
#                 buffer *= c
#                 if lowercase(buffer) in ["include", "define"]
#                     keyword, buffer = buffer, ""
#                 end
#             elseif c != ' ' && c != '\t' || !isempty(buffer)
#                 buffer *= c
#                 isValid = false
#             end
#         elseif lowercase(keyword) == "include"
#             if isUserLib === nothing
#                 if c == '"'
#                     isUserLib = true
#                     Scan.retract(sc)
#                 elseif c == '<'
#                     isUserLib = false
#                     Scan.retract(sc)
#                 elseif c != ' ' && c != '\t' || !isempty(buffer)
#                     buffer *= c
#                 end
#             elseif isUserLib
#                 buffer, isValid = __preinclude(buffer, c, tokens, isValid, '"', '"')
#             else
#                 buffer, isValid = __preinclude(buffer, c, tokens, isValid, '<', '>')
#             end
#         elseif lowercase(keyword) == "define"
#         elseif keyword == "finish"
#             if c != ' ' && c != '\t' || !isempty(buffer)
#                 buffer *= c
#             end
#         end
#         c = Scan.next(sc)
#     end
#     buffer = strip(buffer)
#     Scan.retract(sc)
#     if !isempty(buffer)
#         single = findfirst("//", buffer)
#         multiple = findfirst("/*", buffer)
#         if single !== nothing || multiple !== nothing
#             single = (single === nothing ? typemax(Int) : single[1])
#             multiple = (multiple === nothing ? typemax(Int) : multiple[1])
#             minVal = min(single, multiple) - 1
#             push!(tokens, (TokenType.UNIDENTIFIER, rstrip(buffer[1:minVal])))
#             Scan.retract(sc, lastindex(buffer) - minVal)
#         else
#             push!(tokens, (TokenType.UNIDENTIFIER, buffer))
#         end
#     end
#     sc.state = Scan.START
# end

function identifier!(sc::Scanner, c::Char, tokens::Array)
    id, c = string(c), Scan.next(sc)
    while c !== nothing && c != '\r' && c != '\n'
        if c >= 'a' && c <= 'z' || c >= 'A' && c <= 'Z' || c == '_' || c >= '0' && c <= '9'
            id *= c
        # allow multiple pointer
        elseif sc.specifier == '*' && c == '*' && last(id) == '*'
            while c == '*'
                id *= c
                c = Scan.next(sc, true)
            end
            Scan.retract(sc)
        else
            Scan.retract(sc)
            break
        end
        c = Scan.next(sc)
    end
    if c === nothing || c == '\r' || c == '\n'
        Scan.retract(sc)
    end
    pType = TokenType.toPreversedType(id)
    if pType !== nothing
        if sc.specifier !== nothing || sc.sign !== nothing
            throw(UnIdentifierError((sc.sign === nothing ? "" : sc.sign) * (sc.specifier === nothing ? "" : sc.specifier) * id))
        end
        push!(tokens, (pType, nothing))
    else
        tempArr = []
        if sc.specifier == '*'
            if isempty(tokens)
                throw(UnIdentifierError((sc.sign === nothing ? "" : sc.sign) * sc.specifier * id))
            end
            if isempty(findall(x -> x[1] in [TokenType.POINTER_IDENTIFIER] && x[2] == id, tokens))
                for type in reverse(tokens)
                    if type[1] == TokenType.LN_TOKEN
                        throw(UnIdentifierError((sc.sign === nothing ? "" : sc.sign) * sc.specifier * id))
                    elseif type[1] in [TokenType.INT, TokenType.CHAR, TokenType.DOUBLE, TokenType.FLOAT]
                        break
                    end
                end
            end
            push!(tempArr, (TokenType.POINTER_IDENTIFIER, id))
        elseif sc.specifier == '&'
            # if variable exists
            if !isempty(findall(x -> x[1] == TokenType.IDENTIFIER && x[2] == id, tokens))
                push!(tempArr, (TokenType.ADDRESS_IDENTIFIER, id))
            else
                throw(UnIdentifierError((sc.sign === nothing ? "" : sc.sign) * sc.specifier * id))
            end
        # if normal variable
        else
            if isempty(tokens)
                throw(UnIdentifierError((sc.sign === nothing ? "" : sc.sign) * id))
            end
            if isempty(findall(x -> x[1] == TokenType.IDENTIFIER && x[2] == id, tokens))
                for type in reverse(tokens)
                    if type[1] == TokenType.LN_TOKEN
                        throw(UnIdentifierError((sc.sign === nothing ? "" : sc.sign) * id))
                    elseif type[1] in [TokenType.INT, TokenType.CHAR, TokenType.DOUBLE, TokenType.FLOAT]
                        break
                    end
                end
            end
            push!(tempArr, (TokenType.IDENTIFIER, id))
        end
        push!(tokens, tempArr...)
    end
    sc.sign = nothing
    sc.specifier = nothing
    sc.state = Scan.START
end

function constant!(sc::Scanner, c::Char, tokens::Array)
    con, c, conType = string(c), Scan.next(sc), nothing
    # determine whether is decimal or not
    if con == "0"
        if c == 'x' || c == 'X'
            conType = TokenType.HEX_CON
        elseif c >= '0' && c <= '7'
            conType = TokenType.OCT_CON
        elseif c == '.'
            conType = TokenType.FLOAT_CON
        end
        if conType === nothing
            # if constant is zero
            conType = TokenType.DEC_CON
        else
            con *= c
            c = Scan.next(sc)
        end
    # float number
    elseif con == "."
        con *= c
        if c >= '0' && c <= '9'
            conType = TokenType.FLOAT_CON
            c = Scan.next(sc)
        else
            throw(UnIdentifierError(con))
        end
    else
        conType = TokenType.DEC_CON
    end
    while c !== nothing && c != '\r' && c != '\n'
        if conType == TokenType.HEX_CON
            if !occursin(r"^[0-9A-Fa-f]$", string(c))
                break
            end
        elseif conType == TokenType.OCT_CON
            if c < '0' || c > '7'
                break
            end
        elseif conType == TokenType.DEC_CON
            if c == '.'
                conType = TokenType.FLOAT_CON
            elseif c < '0' || c > '9'
                break
            end
        elseif conType == TokenType.FLOAT_CON
            if c < '0' || c > '9'
                break
            end
        end
        con *= c
        c = Scan.next(sc)
    end
    Scan.retract(sc)
    if sc.sign == '+' || sc.sign == '-'
        con = sc.sign * con
    end
    push!(tokens, (conType, con))
    sc.sign = nothing
    sc.specifier = nothing
    sc.state = Scan.START
end

function slash!(sc::Scanner, c::Char, tokens::Array)
    buffer = ""
    # // comment
    if c == '/'
        c = Scan.next(sc)
        while c !== nothing && c != '\r' && c != '\n'
            buffer *= c
            c = Scan.next(sc)
        end
        push!(tokens, (TokenType.SINGLE_COMMENT, buffer))
        Scan.retract(sc)
        buffer = ""
    # /* multiple comment
    elseif c == '*'
        c = Scan.next(sc)
        while c !== nothing && !occursin("*/", buffer)
            buffer *= c
            c = Scan.next(sc)
        end
        push!(tokens, (TokenType.MULTIPLE_COMMENT, buffer[1:(end - 2)]))
        Scan.retract(sc)
        buffer = ""
    elseif !isempty(tokens)
        lastTokenType = last(tokens)[1]
        # operator id(++)/(--) /(=)
        if lastTokenType in [TokenType.IDENTIFIER, TokenType.POINTER_IDENTIFIER, TokenType.DPLUS, TokenType.DMINUS]
            # /=
            if c == '='
                push!(tokens, (TokenType.ACC_DIVIDE, nothing))
            # only / need retract
            else
                push!(tokens, (TokenType.DIVIDE, nothing))
                Scan.retract(sc)
            end
        # operator constant /
        elseif lastTokenType in [TokenType.DEC_CON, TokenType.HEX_CON, TokenType.OCT_CON, TokenType.FLOAT_CON]
            push!(tokens, (TokenType.DIVIDE, nothing))
        else
            throw(UnIdentifierError(string(c)))
        end
    else
        throw(UnIdentifierError(string(c)))
    end
    if !isempty(buffer)
        push!(tokens, (TokenType.UNIDENTIFIER, buffer))
    end
    sc.state = Scan.START
end

# deprecated
# function __preinclude(buffer::String, c::Char, tokens::Array, isValid::Bool, prefix::Char, postfix::Char)::Tuple{String, Bool}
#     if (c >= 'a' && c <= 'z' || c >= 'A' && c <= 'Z' || c == '.' || c == '_') && startswith(buffer, string(prefix))
#         # < (invalid)
#         if !isValid
#             push!(tokens, (TokenType.UNIDENTIFIER, rstrip(buffer[2:end])))
#             buffer, isValid = string(prefix), true
#         end
#         buffer *= c
#     elseif c == postfix && startswith(buffer, string(prefix))
#         # <xx(invalid) >
#         matchString = match(Regex("^" * prefix * "[a-zA-Z._]*"), buffer)
#         matchString = matchString === nothing ? "" : matchString.match
#         remain = strip(buffer[(lastindex(matchString) + 1):end])
#         if !isValid && !isempty(remain)
#             push!(tokens, (TokenType.UNIDENTIFIER, remain))
#             isValid = true
#         end
#         push!(tokens, (TokenType.INCLUDE, nothing))
#         push!(tokens, (TokenType.LIBRARY_NAME, matchString * c))
#         buffer, keyword = "", "finish"
#     elseif c == prefix && !(c in buffer)
#         # (invalid) <
#         if !isValid
#             push!(tokens, (TokenType.UNIDENTIFIER, rstrip(buffer)))
#             buffer, isValid = "", true
#         end
#         buffer *= c
#     elseif c != ' ' && c != '\t' || !isempty(buffer[2:end])
#         buffer *= c
#         isValid = false
#     end
#     return buffer, isValid
# end

function __string_seq_stream!(sc::Scanner, startSymbol::Char, endSymbol::Char, isLimitOneLine::Bool)::Tuple{TokenType.Token, Union{Any,Nothing}}
    startType, endType = TokenType.toPunctuationType(startSymbol), TokenType.toPunctuationType(endSymbol)
    strSeq, charSeq, buf, c = [], [], "", Scan.next(sc)
    while c !== nothing && (c != '\r' && c != '\n' || !isLimitOneLine)
        if c == '\\' || c == '%'
            if !isempty(buf)
                push!(charSeq, split(buf, ' ')...)
                push!(strSeq, (TokenType.CHAR_SEQ, tuple(charSeq...)))
                charSeq = []
                buf = ""
            end
            # \ escape char
            if c == '\\'
                c = Scan.next(sc)
                if c == 'a'
                    push!(strSeq, TokenType.ALERT_ESC)
                elseif c == 'b'
                    push!(strSeq, TokenType.BACKSPACE_ESC)
                elseif c == 'e'
                    push!(strSeq, TokenType.ESCAPE_ESC)
                elseif c == 'f'
                    push!(strSeq, TokenType.FORMFEED_ESC)
                elseif c == 'n'
                    push!(strSeq, TokenType.NEWLINE_ESC)
                elseif c == 'r'
                    push!(strSeq, TokenType.RETURN_ESC)
                elseif c == 't'
                    push!(strSeq, TokenType.HTAB_ESC)
                elseif c == 'v'
                    push!(strSeq, TokenType.VTAB_ESC)
                elseif c == '\\'
                    push!(strSeq, TokenType.BACKSLASH_ESC)
                elseif c == '\''
                    push!(strSeq, TokenType.QUOTE_ESC)
                elseif c == '"'
                    push!(strSeq, TokenType.DQUOTE_ESC)
                elseif c == '?'
                    push!(strSeq, TokenType.QUESTION_ESC)
                elseif c == 'u'
                    for i in 1:4
                        c = Scan.next(sc)
                        if c === nothing || !occursin(r"^[0-9A-Fa-f]$", string(c))
                            Scan.retract(sc)
                            break
                        end
                        buf *= c
                    end
                    if occursin(r"^[0-9A-Fa-f]{4}$", buf)
                        push!(strSeq, (TokenType.HEX_ESC, buf))
                    # hexadecimal format error
                    else
                        push!(strSeq, (TokenType.UNIDENTIFIER, "\\u" * buf))
                    end
                elseif c >= '0' && c <= '7'
                    buf *= c
                    times = (c >= '0' && c <= '3' ? 2 : 1)
                    for i in 1:times
                        c = Scan.next(sc)
                        if c === nothing || c < '0' || c > '9'
                            Scan.retract(sc)
                            break
                        end
                        buf *= c
                    end
                    if occursin(r"^[0-3]?[0-9]{1,2}$", buf)
                        push!(strSeq, (TokenType.OCT_ESC, buf))
                    # octal format error
                    else
                        push!(strSeq, (TokenType.UNIDENTIFIER, "\\" * buf))
                    end
                # unknown escape symbol
                else
                    push!(strSeq, (TokenType.UNIDENTIFIER, "\\" * c))
                end
            # % format char
            elseif c == '%'
                c = Scan.next(sc)
                # fetch the format size and align
                while c !== nothing && c != '\r' && c != '\n'
                    buf *= c
                    if c >= 'a' && c <= 'z' || c >= 'A' && c <= 'Z'
                        if !occursin(r"^[dDfFcCsS]$", string(c))
                            buf *= Scan.next(sc)
                        end
                        break
                    end
                    c = Scan.next(sc)
                end
                if c === nothing || c == '\r' || c == '\n'
                    Scan.retract(sc)
                end
                # %x.x
                format = match(r"^[+-]?(([0-9]+\.?[0-9]*)|([0-9]*\.?[0-9]+))?", buf)
                if format !== nothing
                    format = format.match
                    type = lowercase(buf[(lastindex(format) + 1):end])
                    # format type (ignore case)
                    if type == "d"
                        push!(strSeq, (TokenType.DECIMAL_FORMAT, format))
                    elseif type == "l"
                        push!(strSeq, (TokenType.LONG_FORMAT, format))
                    elseif type == "f"
                        push!(strSeq, (TokenType.FLOAT_FORMAT, format))
                    elseif type == "lf"
                        push!(strSeq, (TokenType.DOUBLE_FORMAT, format))
                    elseif type == "c"
                        push!(strSeq, (TokenType.CHAR_FORMAT, format))
                    elseif type == "s"
                        push!(strSeq, (TokenType.STRING_FORMAT, format))
                    else
                        push!(strSeq, (TokenType.UNIDENTIFIER, '%' * buf))
                    end
                else
                    push!(strSeq, (TokenType.UNIDENTIFIER, '%' * buf))
                end
            end
            buf = ""
        elseif c == endSymbol
            if !isempty(buf)
                push!(charSeq, split(rstrip(buf))...)
                push!(strSeq, (TokenType.CHAR_SEQ, tuple(charSeq...)))
            end
            return TokenType.STRING_SEQ, tuple((startType === nothing ? startSymbol : startType), strSeq..., (endType === nothing ? endSymbol : endType))
        else
            buf *= c
        end
        c = Scan.next(sc)
    end
    Scan.retract(sc)
    if !isempty(buf)
        push!(strSeq, buf)
    end
    return TokenType.UNIDENTIFIER, tuple(startSymbol, strSeq..., endSymbol)
end


text:
struct StateSpace{N}
    basis :: Vector{NTuple{N, Int}}
    energy :: Vector{Float64}
    deg :: Vector{BigInt}
end

function StateSpace(Hm, Nm, Hw, Nw)
    const Lm, Lw = length(Hm), length(Hw)

    if Lm < 2
        throw(ArgumentError("monomer Hamiltonian is too short"))
    elseif Lw < 2
        throw(ArgumentError("water Hamiltonian is too short"))
    elseif Nm < 1
        throw(ArgumentError("space has too few monomers"))
    elseif Nw < Nm
        throw(ArgumentError("space has fewer water molecules than monomers"))
    end

    const Nmf, Nwf = factorial(BigInt(Nm)), factorial(BigInt(Nw))

    basis, energy, degeneracy = NTuple{Nm + Nw, Int}[], Float64[], BigInt[]

    state = Array{Int}(Nm + Nw)
    for partition in Combinatorics.IntegerPartitions(Nm)
        if maximum(partition) <= Lm
            i = 1
            for n in partition
                for _ in 1:n
                    state[i] = n
                    i += 1
                end
            end

            menergy = mapreduce(n -> n*Hm[n], +, partition)
            mdeg = Nmf
            for (n, a) in counter(partition)
                mdeg /= factorial(BigInt(a))
                if n != 1
                    mdeg /= 2^a
                end
            end

            t = length(partition)
            state[i:end-t] = Lw
            state[end-t+1:end] = 1
            push!(basis, tuple(state...))

            wenergy = sum(@view Hw[state[i:end]])
            push!(energy, menergy + wenergy)

            wdeg = Nwf
            for (n, a) in counter(@view state[i:end])
                wdeg ÷= factorial(BigInt(a))
            end
            push!(degeneracy, mdeg * wdeg)

            while state[i] > 2
                for j in (length(state)-t):-1:i
                    if state[j] > 2
                        state[j] -= 1
                        state[j+1:length(state)-t] = state[j]
                        break
                    end
                end
                push!(basis, tuple(state...))

                wenergy = sum(@view Hw[state[i:end]])
                push!(energy, menergy + wenergy)

                wdeg = Nwf
                for (n, a) in counter(@view state[i:end])
                    wdeg ÷= factorial(BigInt(a))
                end
                push!(degeneracy, mdeg * wdeg)
            end
        end
    end

    StateSpace{Nm + Nw}(basis, energy, degeneracy)
end

Base.length(s::StateSpace) = length(s.basis)

function degeneracies(s::StateSpace, t::StateSpace)
    deg_states = NTuple{2, NTuple{2,Int}}[]
    N = BigInt(0)
    A = BigInt(0)
    counts = BigInt[]
    for i in 1:length(s)
        for j in 1:length(t)
            e1 = s.energy[i] + t.energy[j]
            d1 = s.deg[i] * t.deg[j]
            A += 0.5 * d1 * (d1 - 1)
            for l in (j+1):length(t)
                e2 = s.energy[i] + t.energy[l]
                d2 = s.deg[i] * t.deg[l]
                if e1 == e2
                    push!(deg_states, ((i,j), (i,l)))
                    push!(counts, d1*d2)
                    N += counts[end]
                end
            end
            for k in (i+1):length(s)
                for l in 1:length(t)
                    e2 = s.energy[k] + t.energy[l]
                    d2 = s.deg[k] * t.deg[l]
                    if e1 == e2
                        push!(deg_states, ((i,j), (k,l)))
                        push!(counts, d1*d2)
                        N += counts[end]
                    end
                end
            end
        end
    end
    N += A
    deg_states, counts/N, A/N
end


text:
# ---
# title: 403. Frog Jump
# id: problem403
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/frog-jump/description/>
# hidden: true
# ---
#
# A frog is crossing a river. The river is divided into x units and at each unit
# there may or may not exist a stone. The frog can jump on a stone, but it must
# not jump into the water.
#
# Given a list of stones' positions (in units) in sorted ascending order,
# determine if the frog is able to cross the river by landing on the last stone.
# Initially, the frog is on the first stone and assume the first jump must be 1
# unit.
#
# If the frog's last jump was _k_ units, then its next jump must be either _k_
# \- 1, _k_ , or _k_ \+ 1 units. Note that the frog can only jump in the forward
# direction.
#
# **Note:**
#
#   * The number of stones is ≥ 2 and is < 1,100.
#   * Each stone's position will be a non-negative integer < 231.
#   * The first stone's position is always 0.
#
# **Example 1:**
#
#
#
#     **[0,1,3,5,6,8,12,17]**
#
#     There are a total of 8 stones.
#     The first stone at the 0th unit, second stone at the 1st unit,
#     third stone at the 3rd unit, and so on...
#     The last stone at the 17th unit.
#
#     **Return true**. The frog can jump to the last stone by jumping
#     1 unit to the 2nd stone, then 2 units to the 3rd stone, then
#     2 units to the 4th stone, then 3 units to the 6th stone,
#     4 units to the 7th stone, and 5 units to the 8th stone.
#
#
# **Example 2:**
#
#
#
#     **[0,1,2,3,4,8,9,11]**
#
#     **Return false**. There is no way to jump to the last stone as
#     the gap between the 5th and 6th stone is too large.
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function conditional_gaussian_closure(sys::MomentEquations,
                                      binary_vars::Array{Int,1}=Int[])

    if isempty(binary_vars)
      error("conditional closure does not work if there are no binary species")
    end
    N = sys.N

    sys = bernoulli_moment_eqs(sys, binary_vars)
    # define symbolic raw moment expressions
    μ = typeof(sys) == CentralMomentEquations ? define_μ(N, sys.q_order) : copy(sys.μ)
    # express cumulants in terms of raw moments
    K = cumulants_to_raw_moments(N, sys.q_order)

    # closure of higher order raw moments without explicit form of truncated moments
    # e.g. μ₁₄ would still be a function of μ₁₃ even though μ₁₃ is also truncated
    closure_μ = OrderedDict()
    # closure of higher order raw moments explicitly expanding closed lower order moments
    # the form we use in the end when solving the ODEs
    closure_μ_exp = OrderedDict()

    # perform conditional gaussian closure on raw moments μ

    # by nonbernoulli we denote moments which cannot be written in the conditional form
    nonbernoulli_iters = filter(x -> sum(x[binary_vars]) == 0, sys.iter_all)
    for order in sys.m_order+1:sys.q_order

        # building the closed moment expressions order by order (due to such hierarchical functional dependency)
        iter_order = filter(x -> sum(x) == order, sys.iter_q)

        for iter in iter_order

            # Two cases: (i) conditional (⟨gp^j⟩=⟨p^j|g=1⟩⟨g⟩) or (ii) basic non bernoulli (⟨p^j⟩)"
            # Here g is a bernoulli variable (e.g gene state), p is any other stochastic variable (e.g no. of proteins)
            if sum(iter[binary_vars]) > 0
                # Case (i): ⟨gp^j⟩=⟨p^j|g=1⟩⟨g⟩
                # boils down to manipulation on indices:
                # 1. ignore the bernoulli variable state and write ⟨gp^j⟩ as μ₀ⱼ
                # 2. perform gaussian closure on μ₀ⱼ (assuming cumulant κ₀ⱼ = 0) and multiply by ⟨g⟩ (μ₁₀)
                # 3. add back bernoulli variable indices, for example:
                #   moment ⟨gp^3⟩ = ⟨p^3|g=1⟩⟨g⟩
                # i)  ⟨gp^3⟩ = μ₁₃ -> μ₀₃
                # ii) μ₀₃ = 3μ₀₂²*μ₀₁-2μ₀₁^3
                # iii) 3μ₀₂²*μ₀₁-2μ₀₁^3 -> 3μ₁₂²*μ₁₁-2μ₁₁^3

                # step 1
                bernoulli_iter = fill(0, N)
                for ind in binary_vars
                    if iter[ind] > 0
                        bernoulli_iter[ind] = 1
                    end
                end
                bernoulli_iter = Tuple(bernoulli_iter)
                r = iter.-bernoulli_iter

                # step 2
                conditional_μ = -μ[bernoulli_iter]*(K[r] - μ[r])
                conditional_μ = simplify(conditional_μ, expand=true)

                # step 3
                iter_conditional = filter(x -> 0 < sum(x) <= sum(r), nonbernoulli_iters)
                conditional_sub = Dict([Pair(μ[iter], μ[iter.+bernoulli_iter]/μ[bernoulli_iter])
                                       for iter in iter_conditional])
                conditional_μ = substitute(conditional_μ, conditional_sub)
                conditional_μ = simplify(conditional_μ, expand=true)

                closure_μ[μ[iter]] = conditional_μ
                closure_μ_exp[μ[iter]] = substitute(conditional_μ, closure_μ_exp)
                closure_μ_exp[μ[iter]] = simplify(closure_μ_exp[μ[iter]], expand=true)
            else
                # Case (ii)

                # Conditional Gaussian closure has no well-defined approximation for terms
                # such as ⟨p^j⟩ which are independent of the Bernoulli variables
                # we assume that the marginal distribution P(p) also follows a Gaussian
                # so that ⟨p^j⟩ is truncated according to normal closure

                moment = simplify(-(K[iter]-μ[iter]), expand=true)

                closure_μ[μ[iter]] = moment
                closure_μ_exp[μ[iter]] = substitute(moment, closure_μ_exp)
                closure_μ_exp[μ[iter]] = simplify(closure_μ_exp[μ[iter]], expand=true)
            end
        end
    end

    if typeof(sys) == CentralMomentEquations

        central_to_raw = central_to_raw_moments(N, sys.q_order)
        μ_central = Dict()
        for iter in vcat(sys.iter_m, sys.iter_q)
            μ_central[μ[iter]] = central_to_raw[iter]
        end

        μ_M_exp = define_μ(N, 1)
        for i in sys.iter_m
            μ_M_exp[i] = μ_central[μ[i]]
        end
        μ_M = copy(μ_M_exp)
        for i in sys.iter_q
            μ_M[i] = closure_μ[μ[i]]
            μ_M[i] = substitute(μ_M[i], μ_central)
            μ_M[i] = simplify(μ_M[i], expand=true)
            μ_M_exp[i] = closure_μ_exp[μ[i]]
            μ_M_exp[i] = substitute(μ_M_exp[i], μ_central)
            μ_M_exp[i] = simplify(μ_M_exp[i], expand=true)
        end
        closure = OrderedDict()
        closure_exp = OrderedDict()
        # construct the corresponding truncated expressions of higher order
        # central moments from the obtained raw moment expressions
        raw_to_central_exp = raw_to_central_moments(N, sys.q_order, μ_M_exp, bernoulli=true)
        raw_to_central = raw_to_central_moments(N, sys.q_order, μ_M, bernoulli=true)
        for i in sys.iter_q
            closure_exp[sys.M[i]] = simplify(raw_to_central_exp[i], expand=true)
            closure[sys.M[i]] = simplify(raw_to_central[i], expand=true)
        end

    else
        closure_exp = closure_μ_exp
        closure = closure_μ
    end

    close_eqs(sys, closure_exp, closure, false)

end


text:
@testset "Hermite normal form" begin
  for i in 1:10
    r = 20
    c = 10
    A = matrix(FlintZZ, rand([0,0,0,0,0,0,0,0,0,0,1], r, c))
    As = sparse_matrix(A)
    @test hnf(A) == fmpz_mat(hnf(As))

    H, T = Hecke.hnf_kannan_bachem(As, Val{true})
    B = deepcopy(As)

    for t in T
      Hecke.apply_left!(B, t)
    end

    @test fmpz_mat(B) == hnf(A)
  end
end


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

module QCircuits

include("Math.jl")
include("QBase.jl")
include("Registers.jl")
include("Gates.jl")
include("Instructions.jl")
include("OtherGates.jl")
include("ComplexGates.jl")
include("Graph.jl")
include("Qiskit.jl")
include("Circuit.jl")

# using QGen.Quantum.QBase
#
# export QuantumCircuit, QuantumDevice, QuantumGate, add!

end # module


text:
# TODO:
# - test bidirectional rnns
# - test new interface
# - test keepstate

include("header.jl")
using Knet: rnntest

if gpu() >= 0; @testset "rnn" begin

    function rmulti(r,xs,hs...)
        h = Any[hs...]
        y = Any[]
        for x in xs
            push!(y, r(x; hidden=h))
        end
        y = reshape(cat1d(y...), size(y[1],1), size(y[1],2), :)
        return y,h
    end
    rmulti1(x...)=rmulti(x...)[1]
    rnewh(r,x,h...;o...)=r(x; hidden=Any[h...], o...)
    eq(a,b)=all(map((x,y)->(x==y==nothing || isapprox(x,y)),a,b))
    gchk(a...)=gradcheck(a...; rtol=0.2, atol=0.05, args=1)
    gnew(a...; o...)=isdefined(AutoGrad,:gcheck) ? gcheck(a...; kw=o, rtol=0.2, atol=0.05) : true
    rnn1(p,r,b=nothing)=rnnforw(r,p...;batchSizes=b)[1]
    D,X,H,B,T = Float64,32,32,16,8 # Keep X==H to test skipInput

    # rnew=r=w=x1=x2=x3=hx1=cx1=hx2=cx2=hx3=cx3=nothing
    # rcpu=wcpu=x1cpu=x2cpu=x3cpu=hx1cpu=cx1cpu=hx2cpu=cx2cpu=hx3cpu=cx3cpu=nothing

    for M=(:relu,:tanh,:lstm,:gru), L=1:2, I=(:false,:true), BI=(:false,:true)
        # println((:rnninit,X,H,:dataType,D, :rnnType,M, :numLayers,L, :skipInput,I, :bidirectional,BI, :binit,xavier))
        # global rnew,r,w,x1,x2,x3,hx1,cx1,hx2,cx2,hx3,cx3
        # global rcpu,wcpu,x1cpu,x2cpu,x3cpu,hx1cpu,cx1cpu,hx2cpu,cx2cpu,hx3cpu,cx3cpu
        Knet.seed!(2)

        (r,w) = rnninit(X, H; dataType=D, rnnType=M, numLayers=L, skipInput=I, bidirectional=BI, binit=xavier) # binit=zeros does not pass gchk
        (rcpu,wcpu) = rnninit(X, H; dataType=D, rnnType=M, numLayers=L, skipInput=I, bidirectional=BI, binit=xavier, usegpu=false)
        @test eltype(wcpu) == eltype(w)
        @test size(wcpu) == size(w)
        wcpu = Array(w)
        HL = BI ? 2L : L
        BT = B*T

        # rnnforw(r...) calls the cudnn gpu implementation
        # rnntest(r...) calls the manual cpu/gpu implementation (does not support batchSizes)
        # rnnforw(rcpu...) calls rnntest for cpu inputs
        # (::RNN)(x...) is the new interface

        ## Test 1-D x
        x1cpu = randn(D,X); x1 = ka(x1cpu)
        hx1cpu = randn(D,H,1,HL); hx1 = ka(hx1cpu)
        cx1cpu = randn(D,H,1,HL); cx1 = ka(cx1cpu)
        # x
        @test eq(rnnforw(r,w,x1),rnntest(r,w,x1))
        @test eq(rnnforw(r,w,x1),rnnforw(rcpu,wcpu,x1cpu))
        @test gchk(rnn1,[w,x1],r)
        @test gchk(rnn1,[wcpu,x1cpu],rcpu)
        # x,hidden
        @test eq(rnnforw(r,w,x1,hx1,cx1),rnntest(r,w,x1,hx1,cx1))
        @test eq(rnnforw(r,w,x1,hx1,cx1),rnnforw(rcpu,wcpu,x1cpu,hx1cpu,cx1cpu))
        @test gchk(rnn1,[w,x1,hx1,cx1],r)
        @test gchk(rnn1,[wcpu,x1cpu,hx1cpu,cx1cpu],rcpu)
        # x,batchSizes
        @test eq(rnnforw(r,w,x1;batchSizes=[1]),rnntest(r,w,x1))
       #@test eq(rnnforw(r,w,x1;batchSizes=[1]),rnnforw(rcpu,wcpu,x1cpu;batchSizes=[1]))
        @test gchk(rnn1,[w,x1],r,[1])
       #@test gchk(rnn1,[wcpu,x1cpu],rcpu,[1]) # TODO
        # x,hidden,batchSizes
        @test eq(rnnforw(r,w,x1,hx1,cx1;batchSizes=[1]),rnntest(r,w,x1,hx1,cx1))
       #@test eq(rnnforw(r,w,x1,hx1,cx1;batchSizes=[1]),rnnforw(rcpu,wcpu,x1cpu,hx1cpu,cx1cpu;batchSizes=[1]))
        @test gchk(rnn1,[w,x1,hx1,cx1],r,[1])
       #@test gchk(rnn1,[wcpu,x1cpu,hx1cpu,cx1cpu],rcpu,[1]) # TODO

        ## Test 2-D x
        x2cpu =  randn(D,X,B); x2 = ka(x2cpu)
        hx2cpu = randn(D,H,B,HL); hx2 = ka(hx2cpu)
        cx2cpu = randn(D,H,B,HL); cx2 = ka(cx2cpu)
        # x
        @test eq(rnnforw(r,w,x2),rnntest(r,w,x2))
        @test eq(rnnforw(r,w,x2),rnnforw(rcpu,wcpu,x2cpu))
        @test gchk(rnn1,[w,x2],r)
        @test gchk(rnn1,[wcpu,x2cpu],rcpu)
        # x,hidden
        @test eq(rnnforw(r,w,x2,hx2,cx2),rnntest(r,w,x2,hx2,cx2))
        @test eq(rnnforw(r,w,x2,hx2,cx2),rnnforw(rcpu,wcpu,x2cpu,hx2cpu,cx2cpu))
        @test gchk(rnn1,[w,x2,hx2,cx2],r)
        @test gchk(rnn1,[wcpu,x2cpu,hx2cpu,cx2cpu],rcpu)
        # x,hidden,batchSizes
        @test eq(rnnforw(r,w,x2,hx2,cx2;batchSizes=[B]),rnntest(r,w,x2,hx2,cx2))
       #@test eq(rnnforw(r,w,x2,hx2,cx2;batchSizes=[B]),rnnforw(rcpu,wcpu,x2cpu,hx2cpu,cx2cpu;batchSizes=[B]))
        for b in ([B],[B÷2,B÷2],[B÷2,B÷4,B÷4])
            hx2a = ka(randn(D,H,b[1],HL))
            cx2a = ka(randn(D,H,b[1],HL))
            @test gchk(rnn1,[w,x2,hx2a,cx2a],r,b)
            @test gchk(rnn1,[w,x2],r,b)
            #@test gchk(rnn1,[wcpu,x2cpu,hx2cpu,cx2cpu],rcpu,b) # TODO
        end

        ## Test 3-D x
        x3cpu = randn(D,X,B,T); x3 = ka(x3cpu)
        hx3cpu = randn(D,H,B,HL); hx3 = ka(hx3cpu)
        cx3cpu = randn(D,H,B,HL); cx3 = ka(cx3cpu)
        # x
        @test eq(rnnforw(r,w,x3),rnntest(r,w,x3))
        @test eq(rnnforw(r,w,x3),rnnforw(rcpu,wcpu,x3cpu))
        @test gchk(rnn1,[w,x3],r)
        @test gchk(rnn1,[wcpu,x3cpu],rcpu)
        # x,hidden
        @test eq(rnnforw(r,w,x3,hx3,cx3),rnntest(r,w,x3,hx3,cx3))
        @test eq(rnnforw(r,w,x3,hx3,cx3),rnnforw(rcpu,wcpu,x3cpu,hx3cpu,cx3cpu))
        @test gchk(rnn1,[w,x3,hx3,cx3],r)
        @test gchk(rnn1,[wcpu,x3cpu,hx3cpu,cx3cpu],rcpu)
        # x,hidden,batchSizes
        @test eq(rnnforw(r,w,x3,hx3,cx3;batchSizes=[B for t=1:T]),rnntest(r,w,x3,hx3,cx3))
       #@test eq(rnnforw(r,w,x3,hx3,cx3;batchSizes=[B]),rnnforw(rcpu,wcpu,x3cpu,hx3cpu,cx3cpu;batchSizes=[B]))
        for b in ([BT],[BT÷2,BT÷2],[BT÷2,BT÷4,BT÷4])
            hx3a = ka(randn(D,H,b[1],HL))
            cx3a = ka(randn(D,H,b[1],HL))
            @test gchk(rnn1,[w,x3,hx3a,cx3a],r,b)
            @test gchk(rnn1,[w,x3],r,b)
            # @test gchk(rnn1,[wcpu,x3cpu,hx3cpu,cx3cpu],rcpu,b) # TODO
        end

        ## Test new interface in 3-D
        rnew = RNN(X, H; dataType=D, rnnType=M, numLayers=L, skipInput=I, bidirectional=BI, binit=xavier)
        rnew.w = param(w)
        # x
        yold1 = rnnforw(r,w,x3)
        ynew1 = rnew(x3)
        @test isapprox(yold1[1], ynew1)
        @test gnew(rnew, x3)
        # x,hidden
        yold2 = rnnforw(r,w,x3,hx3,cx3)
        hnew2 = Any[hx3,cx3]
        ynew2 = rnew(x3, hidden=hnew2)
        @test isapprox(yold2[1], ynew2)
        @test isapprox(yold2[2], hnew2[1])
        if M == :lstm; @test isapprox(yold2[3], hnew2[2]); end
        @test gnew(rnewh, rnew, Param(x3), Param(hx3), Param(cx3))
        # x,hidden,batchSizes
        bs = [B for t=1:T]
        yold3 = rnnforw(r,w,x3,hx3,cx3; batchSizes=bs)
        hnew3 = Any[hx3,cx3]
        ynew3 = rnew(x3, hidden=hnew3, batchSizes=bs)
        @test isapprox(yold3[1], ynew3)
        @test isapprox(yold3[2], hnew3[1])
        if M == :lstm; @test isapprox(yold3[3], hnew3[2]); end
        @test gnew(rnewh, rnew, Param(x3), Param(hx3), Param(cx3); batchSizes=bs)

        # compare result and diff between multi-step and single-step
        if !BI
            T2 = T÷2
            x3a = x3[:,:,1:T2]
            x3b = x3[:,:,1+T2:end]
            y1,h1 = rmulti(rnew, Any[x3], hx3, cx3)
            y2,h2 = rmulti(rnew, Any[x3a,x3b], hx3, cx3)
            @test isapprox(y1,y2)
            @test isapprox(h1[1],h2[1])
            if M == :lstm; @test isapprox(h1[2],h2[2]); end
            @test gnew(rmulti1, rnew, Any[x3], hx3, cx3)
            @test gnew(rmulti1, rnew, Any[x3a,x3b], hx3, cx3)
        end

        # rnnparam, rnnparams
        for m in (1,2)
            for l in 1:L
                for i in 1:(M==:lstm ? 8 : M==:gru ? 6 : 2)
                    #@show M,L,I,l,i,m
                    @test rnnparam(r,w,l,i,m) == rnnparam(r,wcpu,l,i,m)
                end
            end
        end
        @test all(map(==, rnnparams(r,w), rnnparams(r,wcpu)))
    end # for
end # @testset begin
end # if gpu() >= 0

nothing


text:
# ------------------------------------------------------------------------------------------
# # Character based RNN language model
# (c) Deniz Yuret, 2019. Based on http://karpathy.github.io/2015/05/21/rnn-effectiveness.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# * Objectives: Learn to define and train a character based language model and generate text
# from it. Minibatch blocks of text. Keep a persistent RNN state between updates. Train a
# Shakespeare generator and a Julia programmer using the same type of model.
# * Prerequisites: [RNN basics](60.rnn.ipynb), [Iterators](25.iterators.ipynb)
# * New functions:
# [converge](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.converge)
# ------------------------------------------------------------------------------------------

# Set display width, load packages, import symbols
ENV["COLUMNS"]=72
using Pkg; haskey(Pkg.installed(),"Knet") || Pkg.add("Knet")
using Statistics: mean
using Base.Iterators: cycle
using Knet: Knet, AutoGrad, Data, param, param0, mat, RNN, dropout, value, nll, adam, minibatch, progress!, converge

# ------------------------------------------------------------------------------------------
# ## Define the model
# ------------------------------------------------------------------------------------------

struct Embed; w; end

Embed(vocab::Int,embed::Int)=Embed(param(embed,vocab))

(e::Embed)(x) = e.w[:,x]  # (B,T)->(X,B,T)->rnn->(H,B,T)

struct Linear; w; b; end

Linear(input::Int, output::Int)=Linear(param(output,input), param0(output))

(l::Linear)(x) = l.w * mat(x,dims=1) .+ l.b  # (H,B,T)->(H,B*T)->(V,B*T)

# Let's define a chain of layers
struct Chain
    layers
    Chain(layers...) = new(layers)
end
(c::Chain)(x) = (for l in c.layers; x = l(x); end; x)
(c::Chain)(x,y) = nll(c(x),y)
(c::Chain)(d::Data) = mean(c(x,y) for (x,y) in d)

# The h=0,c=0 options to RNN enable a persistent state between iterations
CharLM(vocab::Int,embed::Int,hidden::Int; o...) =
    Chain(Embed(vocab,embed), RNN(embed,hidden;h=0,c=0,o...), Linear(hidden,vocab))

# ------------------------------------------------------------------------------------------
# ## Train and test utilities
# ------------------------------------------------------------------------------------------

# For running experiments
function trainresults(file,maker,chars)
    if (print("Train from scratch? "); readline()[1]=='y')
        model = maker()
        a = adam(model,cycle(dtrn))
        b = (exp(model(dtst)) for _ in every(100,a))
        c = converge(b, alpha=0.1)
        progress!(c, alpha=1)
        Knet.save(file,"model",model,"chars",chars)
    else
        isfile(file) || download("http://people.csail.mit.edu/deniz/models/tutorial/$file",file)
        model,chars = Knet.load(file,"model","chars")
    end
    Knet.gc() # To save gpu memory
    return model,chars
end

every(n,itr) = (x for (i,x) in enumerate(itr) if i%n == 0);

# To generate text from trained models
function generate(model,chars,n)
    function sample(y)
        p = Array(exp.(y)); r = rand()*sum(p)
        for j=1:length(p); (r -= p[j]) < 0 && return j; end
    end
    x = 1
    reset!(model)
    for i=1:n
        y = model([x])
        x = sample(y)
        print(chars[x])
    end
    println()
end

reset!(m::Chain)=(for r in m.layers; r isa RNN && (r.c=r.h=0); end);

# ------------------------------------------------------------------------------------------
# ## The Complete Works of William Shakespeare
# ------------------------------------------------------------------------------------------

RNNTYPE = :lstm
BATCHSIZE = 256
SEQLENGTH = 100
VOCABSIZE = 84
INPUTSIZE = 168
HIDDENSIZE = 334
NUMLAYERS = 1;

# Load 'The Complete Works of William Shakespeare'
include(Knet.dir("data","gutenberg.jl"))
trn,tst,shakechars = shakespeare()
map(summary,(trn,tst,shakechars))

# Print a sample
println(string(shakechars[trn[1020:1210]]...))

# Minibatch data
function mb(a)
    N = length(a) ÷ BATCHSIZE
    x = reshape(a[1:N*BATCHSIZE],N,BATCHSIZE)' # reshape full data to (B,N) with contiguous rows
    minibatch(x[:,1:N-1], x[:,2:N], SEQLENGTH) # split into (B,T) blocks
end
dtrn,dtst = mb.((trn,tst))
length.((dtrn,dtst))

summary.(first(dtrn))  # each x and y have dimensions (BATCHSIZE,SEQLENGTH)

# 3.30e+00  ┣   /       /       /       /       /    ┫ 122 [04:46, 2.35s/i]
Knet.gc()
shakemaker() = CharLM(VOCABSIZE, INPUTSIZE, HIDDENSIZE; rnnType=RNNTYPE, numLayers=NUMLAYERS)
shakemodel,shakechars = trainresults("shakespeare113.jld2", shakemaker, shakechars);

#exp(shakemodel(dtst))  # Perplexity = 3.30

generate(shakemodel,shakechars,1000)

# ------------------------------------------------------------------------------------------
# ## Julia programmer
# ------------------------------------------------------------------------------------------

RNNTYPE = :lstm
BATCHSIZE = 64
SEQLENGTH = 64
INPUTSIZE = 512
VOCABSIZE = 128
HIDDENSIZE = 512
NUMLAYERS = 2;

# Read julia base library source code
base = joinpath(Sys.BINDIR, Base.DATAROOTDIR, "julia")
text = ""
for (root,dirs,files) in walkdir(base)
    for f in files
        f[end-2:end] == ".jl" || continue
        text *= read(joinpath(root,f), String)
    end
    # println((root,length(files),all(f->contains(f,".jl"),files)))
end
length(text)

# Find unique chars, sort by frequency, assign integer ids.
charcnt = Dict{Char,Int}()
for c in text; charcnt[c]=1+get(charcnt,c,0); end
juliachars = sort(collect(keys(charcnt)), by=(x->charcnt[x]), rev=true)
charid = Dict{Char,Int}()
for i=1:length(juliachars); charid[juliachars[i]]=i; end
hcat(juliachars, map(c->charcnt[c],juliachars))

# Keep only VOCABSIZE most frequent chars, split into train and test
data = map(c->charid[c], collect(text))
data[data .> VOCABSIZE] .= VOCABSIZE
ntst = 1<<19
tst = data[1:ntst]
trn = data[1+ntst:end]
length.((data,trn,tst))

# Print a sample
r = rand(1:(length(trn)-1000))
println(string(juliachars[trn[r:r+1000]]...))

# Minibatch data
function mb(a)
    N = length(a) ÷ BATCHSIZE
    x = reshape(a[1:N*BATCHSIZE],N,BATCHSIZE)' # reshape full data to (B,N) with contiguous rows
    minibatch(x[:,1:N-1], x[:,2:N], SEQLENGTH) # split into (B,T) blocks
end
dtrn,dtst = mb.((trn,tst))
length.((dtrn,dtst))

summary.(first(dtrn))  # each x and y have dimensions (BATCHSIZE,SEQLENGTH)

# 3.25e+00  ┣       /       /       /       /       /┫ 126 [05:43, 2.72s/i]
juliamaker() = CharLM(VOCABSIZE, INPUTSIZE, HIDDENSIZE; rnnType=RNNTYPE, numLayers=NUMLAYERS)
juliamodel,juliachars = trainresults("juliacharlm113.jld2", juliamaker, juliachars);

#exp(juliamodel(dtst))  # Perplexity = 3.27

generate(juliamodel,juliachars,1000)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials