text:
# Series evaluation of:
#    M_n/(4*b*r)^{n/2} = \int_{-\kappa/2}^{\kappa/2} (k^2-\sin^2{x})^{n/2} dx
# for the top four values of n = n_max-3 to n_max.
function Mn_series!(t::Transit_Struct{T}) where {T <: Real}
# Use series expansion to compute M_n:
k2 = t.k2; n_max = t.n_max
# Computing leading coefficient (n=0):
if k2 < 1
  tol = eps(k2); term = zero(T)
  fac = t.sqonembmr2^(n_max-3)*t.k
# Now, compute higher order terms until desired precision is reached:
  @inbounds for j=1:4
    # Add leading term to M_n:
    t.Mn[n_max-3+j] = t.Mn_coeff[1,j,1]
    k2n = one(T)   # k^{2n}
    @inbounds for n =1:t.jmax-1
      k2n *= k2
      term = k2n*t.Mn_coeff[1,j,n+1]
      t.Mn[n_max-3+j] += term
      if abs(term) < tol
        break
      end
    end
    t.Mn[n_max-3+j] *= fac
    fac *= t.sqonembmr2
  end
  return
else # k^2 >= 1
  t.k2inv = inv(k2)
  tol = eps(t.k2inv); term = zero(T)
  fac = t.sqonembmr2^(n_max-3)
  @inbounds for j=1:4
    t.Mn[n_max-3+j] = t.Mn_coeff[2,j,1]
    k2n = one(T)
    @inbounds for n = 1:t.jmax-1
      k2n *= t.k2inv
      term = k2n*t.Mn_coeff[2,j,n+1]
      t.Mn[n_max-3+j] += term
      if abs(term) < tol
        break
      end
    end
    t.Mn[n_max-3+j] *= fac
    fac *= t.sqonembmr2
  end
  return
end
end

# Recursive computation of M_n starting at m=1 to 4, and raising to n_max:
function Mn_raise!(t::Transit_Struct{T}) where {T <: Real}
r=t.r; b=t.b; n_max=t.n_max; k2 = t.k2; kc = t.kc; t.k2inv = inv(t.k2); sqarea=t.sqarea; onemr2mb2=t.onemr2mb2
# Computes the integrals:
# M_n(r,b) = (4*b*r)^(m/2) \int_{-\kappa/2}^{\kappa/2} (k^2 - \sin^2{x})^{m/2} dx
# where k^2 = (1-(r-b)^2)/(4*b*r), \kappa/2 = \sin^{-1}(k) for k<1; otherwise \kappa=\pi.
# Compute first four integrals:
Mn_four!(t)
@inbounds for m=4:n_max
  t.Mn[m+1]=(2*(m-1)*onemr2mb2*t.Mn[m-1]+(m-2)*sqarea*t.Mn[m-3])*t.ninv[m]
end
return
end

# Compute M_n for m=0 to 3:
function Mn_four!(t::Transit_Struct{T}) where {T <: Real}
r=t.r; b=t.b; n_max=t.n_max; k2 = t.k2; kc = t.kc; t.k2inv = inv(t.k2)
# Computes the integrals:
# M_n(r,b) = (4*b*r)^(m/2) \int_{-\kappa/2}^{\kappa/2} (k^2 - \sin^2{x})^{m/2} dx
# where k^2 = (1-(r-b)^2)/(4*b*r), \kappa/2 = \sin^{-1}(k) for k<1; otherwise \kappa=\pi.
# Compute first four integrals:
if k2 < 1.0
  t.Mn[1] = t.kap0
  t.Mn[2] = 2*t.sqbr*2*k2*t.Em1mKdm
  t.Mn[3] = t.Mn[1]*t.onemr2mb2+t.kite_area2
  t.Mn[4] = (2*t.sqbr)^3*t.twothird*k2*(t.Eofk+(3*k2-2)*t.Em1mKdm)
else
  t.Mn[1] = pi
  t.Mn[2] = 2*t.sqonembmr2*t.Eofk
  t.Mn[3] = pi*t.onemr2mb2
  t.Mn[4] = t.sqonembmr2^3*t.twothird*((3-2*t.k2inv)*t.Eofk+t.k2inv*t.Em1mKdm)
end
return
end

# Recursive computation of M_n starting at m=n_max to n_max-3, and lowering to m=1:
function Mn_lower!(t::Transit_Struct{T}) where {T <: Real}
r=t.r; b=t.b; n_max=t.n_max; k2 = t.k2
# Computes the integrals:
# M_n(r,b) = (4*b*r)^(m/2) \int_{-\kappa/2}^{\kappa/2} (k^2 - \sin^2{x})^{m/2} dx
# where k^2 = (1-(r-b)^2)/(4*b*r), \kappa/2 = \sin^{-1}(k) for k<1; otherwise \kappa=\pi.
# Compute series version:
Mn_series!(t)
# Now iterate downwards:
invsqarea = inv(t.sqarea)
@inbounds for m=n_max-4:-1:4
  t.Mn[m+1]=((m+4)*t.Mn[m+5]-2*(m+3)*t.onemr2mb2*t.Mn[m+3])*invsqarea*t.ninv[m+2]
end
# Now, compute lowest four exactly:
Mn_four!(t)
return
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct TrackableObjectGraph_TrackableObject_ObjectReference <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TrackableObjectGraph_TrackableObject_ObjectReference(; kwargs...)
        obj = new(meta(TrackableObjectGraph_TrackableObject_ObjectReference), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TrackableObjectGraph_TrackableObject_ObjectReference
const __meta_TrackableObjectGraph_TrackableObject_ObjectReference = Ref{ProtoMeta}()
function meta(::Type{TrackableObjectGraph_TrackableObject_ObjectReference})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TrackableObjectGraph_TrackableObject_ObjectReference)
            __meta_TrackableObjectGraph_TrackableObject_ObjectReference[] = target = ProtoMeta(TrackableObjectGraph_TrackableObject_ObjectReference)
            allflds = Pair{Symbol,Union{Type,String}}[:node_id => Int32, :local_name => AbstractString]
            meta(target, TrackableObjectGraph_TrackableObject_ObjectReference, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TrackableObjectGraph_TrackableObject_ObjectReference[]
    end
end
function Base.getproperty(obj::TrackableObjectGraph_TrackableObject_ObjectReference, name::Symbol)
    if name === :node_id
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :local_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct TrackableObjectGraph_TrackableObject_SerializedTensor <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TrackableObjectGraph_TrackableObject_SerializedTensor(; kwargs...)
        obj = new(meta(TrackableObjectGraph_TrackableObject_SerializedTensor), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TrackableObjectGraph_TrackableObject_SerializedTensor
const __meta_TrackableObjectGraph_TrackableObject_SerializedTensor = Ref{ProtoMeta}()
function meta(::Type{TrackableObjectGraph_TrackableObject_SerializedTensor})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TrackableObjectGraph_TrackableObject_SerializedTensor)
            __meta_TrackableObjectGraph_TrackableObject_SerializedTensor[] = target = ProtoMeta(TrackableObjectGraph_TrackableObject_SerializedTensor)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :full_name => AbstractString, :checkpoint_key => AbstractString, :optional_restore => Bool]
            meta(target, TrackableObjectGraph_TrackableObject_SerializedTensor, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TrackableObjectGraph_TrackableObject_SerializedTensor[]
    end
end
function Base.getproperty(obj::TrackableObjectGraph_TrackableObject_SerializedTensor, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :full_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :checkpoint_key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :optional_restore
        return (obj.__protobuf_jl_internal_values[name])::Bool
    else
        getfield(obj, name)
    end
end

mutable struct TrackableObjectGraph_TrackableObject_SlotVariableReference <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TrackableObjectGraph_TrackableObject_SlotVariableReference(; kwargs...)
        obj = new(meta(TrackableObjectGraph_TrackableObject_SlotVariableReference), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TrackableObjectGraph_TrackableObject_SlotVariableReference
const __meta_TrackableObjectGraph_TrackableObject_SlotVariableReference = Ref{ProtoMeta}()
function meta(::Type{TrackableObjectGraph_TrackableObject_SlotVariableReference})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TrackableObjectGraph_TrackableObject_SlotVariableReference)
            __meta_TrackableObjectGraph_TrackableObject_SlotVariableReference[] = target = ProtoMeta(TrackableObjectGraph_TrackableObject_SlotVariableReference)
            allflds = Pair{Symbol,Union{Type,String}}[:original_variable_node_id => Int32, :slot_name => AbstractString, :slot_variable_node_id => Int32]
            meta(target, TrackableObjectGraph_TrackableObject_SlotVariableReference, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TrackableObjectGraph_TrackableObject_SlotVariableReference[]
    end
end
function Base.getproperty(obj::TrackableObjectGraph_TrackableObject_SlotVariableReference, name::Symbol)
    if name === :original_variable_node_id
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :slot_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :slot_variable_node_id
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

mutable struct TrackableObjectGraph_TrackableObject <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TrackableObjectGraph_TrackableObject(; kwargs...)
        obj = new(meta(TrackableObjectGraph_TrackableObject), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TrackableObjectGraph_TrackableObject
const __meta_TrackableObjectGraph_TrackableObject = Ref{ProtoMeta}()
function meta(::Type{TrackableObjectGraph_TrackableObject})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TrackableObjectGraph_TrackableObject)
            __meta_TrackableObjectGraph_TrackableObject[] = target = ProtoMeta(TrackableObjectGraph_TrackableObject)
            allflds = Pair{Symbol,Union{Type,String}}[:children => Base.Vector{TrackableObjectGraph_TrackableObject_ObjectReference}, :attributes => Base.Vector{TrackableObjectGraph_TrackableObject_SerializedTensor}, :slot_variables => Base.Vector{TrackableObjectGraph_TrackableObject_SlotVariableReference}]
            meta(target, TrackableObjectGraph_TrackableObject, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TrackableObjectGraph_TrackableObject[]
    end
end
function Base.getproperty(obj::TrackableObjectGraph_TrackableObject, name::Symbol)
    if name === :children
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{TrackableObjectGraph_TrackableObject_ObjectReference}
    elseif name === :attributes
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{TrackableObjectGraph_TrackableObject_SerializedTensor}
    elseif name === :slot_variables
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{TrackableObjectGraph_TrackableObject_SlotVariableReference}
    else
        getfield(obj, name)
    end
end

mutable struct TrackableObjectGraph <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TrackableObjectGraph(; kwargs...)
        obj = new(meta(TrackableObjectGraph), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TrackableObjectGraph
const __meta_TrackableObjectGraph = Ref{ProtoMeta}()
function meta(::Type{TrackableObjectGraph})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TrackableObjectGraph)
            __meta_TrackableObjectGraph[] = target = ProtoMeta(TrackableObjectGraph)
            allflds = Pair{Symbol,Union{Type,String}}[:nodes => Base.Vector{TrackableObjectGraph_TrackableObject}]
            meta(target, TrackableObjectGraph, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TrackableObjectGraph[]
    end
end
function Base.getproperty(obj::TrackableObjectGraph, name::Symbol)
    if name === :nodes
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{TrackableObjectGraph_TrackableObject}
    else
        getfield(obj, name)
    end
end

export TrackableObjectGraph_TrackableObject_ObjectReference, TrackableObjectGraph_TrackableObject_SerializedTensor, TrackableObjectGraph_TrackableObject_SlotVariableReference, TrackableObjectGraph_TrackableObject, TrackableObjectGraph

ext: jl


text:
module HilbertSchmidtIndependenceCriterion

	using StatsFuns, Distances

  # load Base modules
  using Statistics, LinearAlgebra

	# includes
	include("common.jl")
	include("gammaHSIC.jl")

	# function exports
	export gammaHSIC, estimateKernelSize, rbfDotProduct

end # module


text:
# SimpleHttpServer module
#
# Serve HTTP requests and WebSockets in Julia.
#
module SimpleHttpServer

if VERSION < v"0.4.0-dev"
    using Docile
    eval(:(@docstrings(manual = ["../README.md"])))
end

using SimpleHttpCommon
using Compat
using Nettle
using Codecs
using WebSockets
using GnuTLS
using Match
using JSON
import WebSockets: WebSocket, close
import Base: listen, UVError, write
import DataStructures:OrderedDict

const C = SimpleHttpCommon

export HTPHandler,
    Server,
    HTTPServer,
    handle,
    http_handle,
    Request,
    Response,
    headers,
    run,
    listen,
    Iterable,
    readbytes,
    readline,
    readuntil,
    BytesReaderIterator,
    BytesFileIterator,
    LimitedBytesFileIterator,
    FileResponse,
    mh,
    mh_push!

hc_event = SimpleHttpCommon.event

typealias NothingOrUVServer Union(Nothing, Base.UVServer)

type TCPHandler
    handle::Function
    sock::NothingOrUVServer
    events::Dict{String, Function}

    TCPHandler(handle::Function, sock::NothingOrUVServer, events::Dict{String, Function}) = new(handle, sock, events)
    TCPHandler(handle::Function, sock::NothingOrUVServer) = TCPHandler(handle, sock, defaultevents)
    TCPHandler(handle::Function) = TCPHandler(handle, nothing)
end

type Server
    handler::TCPHandler
    error_count::BigInt
    last_client_id::BigInt
    options::Associative{Any, Any}

    Server(handler::TCPHandler, args...) = new(handler, BigInt(0), BigInt(0), Dict{Any,Any}(args))
end

function increase_error(s::Server)
    s.error_count += 1
    return s.error_count
end

function increase_client_id(s::Server)
    s.last_client_id += 1
    return s.last_client_id
end

typealias ClientAddr Union(Nothing, IPv4, IPv6)
typealias ClientHostName Union(Nothing, String)

type Client
    id::BigInt
    server::Server
    sock
    finished::Bool
    addr::ClientAddr
    hostname::ClientHostName

    Client(s::Server, sc, f::Bool) = new(increase_client_id(s), s, sc, f, nothing, nothing)
    Client(s::Server, sc) = Client(s, sc, true)
end

type Error
    id::BigInt
    req::Request
    res
    error
    traceback
end

function print_error(io::IO, e::Error)
    println(io, "[REQUEST_EXCEPTION $(e.id)] BEGIN")
    println(io, "Request error: $(typeof(e.error))")
    println(io, "Request msg: $(e.error)")
    println(io, "Request OPEN: $(isopen(e.req.sock))")
    println(io, "Request: $(e.req)")
    println(io, "")
    showerror(io, e.error, e.traceback)
    println(io, "")
    println(io, "[REQUEST_EXCEPTION $(e.id)] END")
end

function print_respose_iter_error(io::IO, e::Error)
    println(io, "[RESPONSE_ITER_DATA_EXCEPTION $(e.id) BEGIN")
    println(io, "Response: $(e.res)")
    println(io, "")
    showerror(io, e.error, e.traceback)
    println(io, "")
    println(io, "[RESPONSE_ITER_DATA_EXCEPTION $(e.id)] END")
end

function read_ip(sock)
  buffer = Array(Uint8,32)
  bufflen::Int64 = 32

  ccall(:uv_tcp_getpeername,Int64,(Ptr{Void},Ptr{Uint8},Ptr{Int64}), sock.handle,buffer,&bufflen)

  b = buffer[5:8]
  ip = IPv4(b...)
  ip
end

defaultevents = Dict{String, Function}()
defaultevents["error"]  = ( client, err ) -> println( err )
defaultevents["listen"] = ( saddr ) -> println("Listening on $saddr...")
defaultevents["connect"] = ( client ) -> begin
  ip = read_ip(client.sock)
  println("Client connect from: $ip")
  client.addr = ip
end

defaultevents["request_sock_error"] = ( error ) -> begin
    print_error(STDERR, "Sock error: $(error.error)")
end

defaultevents["request_error"] = ( error ) -> begin
    print_error(STDERR, error)
end

defaultevents["post_request_init"] = ( req ) -> begin
    req.events["post_header"] = http_post_header
    req.events["error"] = request_error
    req.events["response_iter_data_error"] = response_iter_data_error
end

defaultevents["response_iter_data_error"] = ( error ) -> begin
    print_respose_iter_error(STDERR, error)
end

function response_iter_data_error(error::Error)
    io = IOBuffer()
    print_respose_iter_error(io, error)
    write_chunk(error.req.sock, takebuf_array(io))
end

function request_error(error::Error, set_res::Function)
    io = IOBuffer()
    print_error(io, error)
    set_res(Response(500, takebuf_array(io)))
end

# copy from: https://github.com/JuliaWeb/WebSockets.jl/blob/master/src/WebSockets.jl
# get key out of request header
get_websocket_key(request::Request) = begin
  return request.headers["Sec-WebSocket-Key"]
end

# copy from: https://github.com/JuliaWeb/WebSockets.jl/blob/master/src/WebSockets.jl
# the protocol requires that a special key
# be processed and sent back with the handshake response
# to prove that received the HTTP request
# and that we *really* know what webSockets means.
function generate_websocket_key(key)
  h = HashState(SHA1)
  update!(h, key*"258EAFA5-E914-47DA-95CA-C5AB0DC85B11")
  bytestring(encode(Base64, digest!(h)))
end

# perform the handshake assuming it's a websocket request
websocket_handshake(req, res) = begin
  key = get_websocket_key(req)
  resp_key = generate_websocket_key(key)

  #TODO: use a proper HTTP response type
  res.status = 101
  res.headers["Upgrade"] = "websocket"
  res.headers["Connection"] = "Upgrade"
  res.headers["Sec-WebSocket-Accept"] = resp_key
  res.data = ""
end

function http_post_header(req::Request, res::Response)
    if !isempty(req.query)
        req.q = parse_qsr(req.query)
    end

    # using web sockets
    if get(req.headers, "Connection", nothing) == "Upgrade"
        hupgrade = get(req.headers, "Upgrade", "")
        if hupgrade == "websocket"
            if get(req.env["server"].options, :ws, true)
                req.websocket = WebSockets.WebSocket(0, req.sock)
                websocket_handshake(req, res)
            else
                throw(HTTPParserHeaderException("Websockts not allowed"))
            end
        else
            throw(HTTPParserHeaderException("Invalid Connection HEADER", "Upgrade: $(hupgrade)"))
        end
    end

    ct = get(req.headers, "Content-Type", "")

    if search(ct, "application/x-www-form-urlencoded").start > 0
        req.p = parse_qsr(C.bodys(req))
    elseif search(ct, "application/json").start > 0
        req.p = JSON.parse(C.bodys(req))
    end
end

handle(h::TCPHandler, c::Client) = h.handle(c)

type HTTPParserHeaderException <: Exception
    message::String
    line::Union(Nothing, String)

    HTTPParserHeaderException(m::String, l::Union(Nothing, String)) = new(m, l)
    HTTPParserHeaderException(m::String) = new(m, nothing)
end

function default_http_processor(client::Client, resource_handler::Function)
    req = Request(client.sock)
    req.env["client"] = client
    req.env["server"] = client.server
    event(client.server, "post_request_init", req)
    res = Response()
    res.req = req
    res_arr = [res]
    write_res = true
    uv_error = false

    try
        try
            parse_http_headers(req)
            hc_event(req, "post_header", req, res)

            # see http://tools.ietf.org/html/rfc7230#section-6.3
            if req.protocol != nothing && get(req.headers, "Connection", nothing) != "close"
                if req.protocol.name == "HTTP" && req.protocol.version >= v"1.1"
                    client.finished = false
                end
            end

            hc_event(req, "pre_handler", req)

            req.finished = client.finished
        catch e
            if isa(e, UVError)
                client.finished = true
                write_res = false
                error = Error(increase_error(client.server), req, res, e, catch_backtrace())
                event(client.server, "request_socket_error", error)
                hc_event(req, "socket_error", error, (new_res) -> begin
                                res_arr[1] = new_res
                            end)
                return
            else
                rethrow()
            end
        end

        if req.websocket != nothing
            write_res = false
            write(req.sock, res)
            resource_handler(req)
            WebSockets.close(req.websocket)
        else
            res2 = resource_handler(req)

            if isa(res2, Response)
                res2.req = req
                res_arr[1] = res2
            end
        end
        hc_event(req, "post_handler", res, (new_res) -> begin
                        res_arr[1] = new_res
                     end)
    catch e
        if isa(e, CloseClientException)
            client.finished = true
            write_res = false
        else
            res = res_arr[1]
            res.status = 500

            error = Error(increase_error(client.server), req, res, e, catch_backtrace())
            event(client.server, "request_error", error)
            hc_event(req, "error", error, (new_res) -> begin
                            res_arr[1] = new_res
                        end)
        end
    finally
        if write_res && isopen(req.sock)
            res = res_arr[1]
            write(req.sock, res)
        end
    end
end

function http_handler(resource_handler::Function)
    callback = client::Client -> default_http_processor(client, resource_handler)
    TCPHandler(callback)
end

Server(handler::Function; args...) = Server(TCPHandler(handler), args...)
HTTPServer(resource_handler::Function; args...) = Server(http_handler(resource_handler), args...)

function event(server::Server, event::String, args...)
    haskey(server.handler.events, event) && server.handler.events[event](args...)
end

@doc "Converts a `Response` to an HTTP response string" ->
function write{T<:IO}(io::T, res::Response)
    write(io, join(["HTTP/1.1", res.status, STATUS_CODES[res.status], "\r\n"], " "))

    it = nothing

    if isa(res.data, Iterable) || isa(res.data, Task)
        it = res.data
    elseif isa(res.data, Function)
        it = Iterable(Task(res.data))
    end

    if it != nothing
        res.headers["Transfer-Encoding"] = "chunked"
    elseif !haskey(res.headers, "Content-Length")
        res.headers["Content-Length"] = string(sizeof(res.data))
    end

    for (header, value) in res.headers
        write(io, @sprintf "%s: %s\r\n" header value)
    end

    for (header, value_list) in res.mheaders
        for value in value_list
            write(io, @sprintf "%s: %s\r\n" header value)
        end
    end

    write(io, "\r\n")

    if it == nothing
        write(io, res.data)
    else
        try
            for data in it.data
                write_chunk(io, data)
            end
        catch e
            bt = catch_backtrace()
            req = res.req
            server = req.env["server"]
            error = Error(increase_error(server), req, res, e, bt)
            event(server, "response_iter_data_error", error)
            hc_event(req, "response_iter_data_error", error)
        end

        write(io, @sprintf "%x\r\n" 0)
        write(io, "\r\n")
        write(io, "\r\n")
    end
end

function listen(server::Server, host::Base.IpAddr, port::Integer)
    server.handler.sock = listen(host, port)
    inet = "$host:$port"
    event(server, "listen", inet)
    return server
end

listen(server::Server, port::Integer) = listen(server, IPv4(0), port)

@unix_only function listen(server::Server, path::String)
    bind(server.handler.sock, path) || throw(ArgumentError("could not listen on path $path"))
    Base.uv_error("listen", Base._listen(server.handler.sock))
    event(server, "listen", path)
    return server
end

function run(server::Server; args...)
    params = Dict(args)

    # parse parameters
    port = get(params, :port, 0)
    host = get(params, :host, IPv4(0))
    use_https = haskey(params, :ssl)
    use_sockets = @unix? (haskey(params, :socket) && !haskey(params, :port)) : false

    server = if use_sockets
        listen(server, params[:socket])  # start server on Unix socket
    else
        listen(server, host, port)       # start server on network socket
    end

    if use_https
        println("HTTPS")
        #handle_https_request(server, params[:ssl])
    else
        handle_http_request(server)
    end
end

@doc """ Handle HTTP request from client """ ->
function handle_http_request(server::Server)
    # handle requests, Base.wait_accept blocks until a connection is made
    while true
        sock = accept(server.handler.sock)
        @async process_sock(server, sock)
    end
end

@doc """ Handle HTTPS request from client """ ->
function handle_https_request(server::Server, cert_store::GnuTLS.CertificateStore)
    while true
        sess = GnuTLS.Session(true)
        set_priority_string!(sess)
        set_credentials!(sess, cert_store)
        client = accept(server.handler.sock)
        try
            associate_stream(sess, client)
            handshake!(sess)
        catch e
            println("Error establishing SSL connection: ", e)
            close(client)
            continue
        end
        sock = accept(server.handler.sock)
        @async process_sock(server, sock)
    end
end

function process_sock(server::Server, sock)
    client = Client(server, sock)
    event(server, "connect", client)

    while isopen(sock)
        try
            handle(server.handler, client)
        catch e
            if isa(e,GnuTLS.GnuTLSException) && e.code == -110
                # Because Chrome is buggy on OS X, ignore E_PREMATURE_TERMINATION
            else
                bt = catch_backtrace()
                err_num = increase_error(server)

                print_error = io -> begin
                    println(io, "[CLIENT_EXCEPTION $err_num] Begin")
                    println(io, "Client: $client")
                    println(io, "Error Type: $(typeof(e))")
                    println(io, "")
                    showerror(io, e, bt)
                    println(io, "")
                    println(io, "[CLIENT_EXCEPTION $err_num] End")
                end

                print_error(STDERR)

                client.finished = true
            end
        finally
            if client.finished
                if isopen(sock)
                    close(sock)
                end

                break
            end
        end
    end

    event(server, "close", client)
end

end # module


text:
using MNPDynamics
using Test


text:
using Genie, Genie.Router
using Genie.Renderer, Genie.Renderer.Html, Genie.Renderer.Json, Genie.Requests

# https://genieframework.github.io/Genie.jl/dev/guides/Interactive_environment.html
function learn_interactive_environment()
    route("/") do
        "Hi there!"
    end

    function hello_world()
        return "Hello World!"
    end

    route("/hello/world", hello_world)

    route("/hello.html") do
        html("Hello World")
    end

    route("/echo/:message") do
        return @params(:message)
    end

    route("/sum/:x::Int/:y::Int") do
        x = @params(:x)
        y = @params(:y)

        return x+y
    end

    route("/sum/:x::Int/:y::Int") do
        x = @params(:x)
        y = @params(:y)
        initial_value = get(@params, :initial_value, "0")
        initial_value = parse(Int, initial_value)

        return x + y + initial_value
    end

    route("/hello.json") do
        json("Hello World")
    end

    route("/hello.txt") do
        respond("Hello World", :text)
    end
end

# https://genieframework.github.io/Genie.jl/dev/guides/Simple_API_backend.html
function learn_simple_api_backend()
    route("/") do
        (:message => "Hi there!") |> json
    end

    route("/echo", method=POST) do
        message = Genie.Requests.jsonpayload()
        Ans = (message["message"] * " ") ^ message["repeat"]
        return (:echo => Ans) |> json
    end
end

#learn_interactive_environment()
learn_simple_api_backend()

up(8001, async=false) # async=false so script does not exit


text:
using LinearAlgebra, Random
using Test
using Quadmath
using IterativeRefinement

const verbose = (get(ENV,"VERBOSITY","0") == "1")

Random.seed!(1101)

include("utils.jl")

function runone(A::Matrix{T},x0::AbstractVector) where {T}
    n = size(A,1)
    DT = widen(T)
    # println("wide type is $DT")
    Ad = DT.(A)
    xd = DT.(x0)
    b = T.(Ad * xd)
    # checkme: Demmel et al. use refined solver here
    xtrue = Ad \ DT.(b)
    xt = T.(xtrue)
    Rv, Cv = equilibrators(A)
    if maximum(abs.(Cv)) > 10
        RA = Diagonal(Rv)*A
    else
        RA = A
    end
    a = opnorm(RA,Inf)
    F = lu(RA,check=false)
    # Some version of OpenBLAS gave exact singularity for one of our "random" cases.
    # Handle gracefully so we can just try again.
    if F.info != 0
        return false
    end
    κnorm = condInfest(RA,F,a)
    RAx = RA*Diagonal(xt)
    a = opnorm(RAx,Inf)
    F = lu(RAx, check=false)
    if F.info != 0
        return false
    end
    κcomp = condInfest(RAx,F,a)
    crit = 1 / (max(sqrt(n),10) * eps(real(T)))
    if verbose
        println("problem difficulty (rel. to convergence criterion):")
        println("normwise: ", κnorm/crit, " componentwise: ", κcomp/crit)
    end
    xhat,Bnorm,Bcomp = @inferred(rfldiv(A,b))
    # xhat,Bnorm,Bcomp = rfldiv(A,b)
    Enorm = norm(xhat-xtrue,Inf)/norm(xtrue,Inf)
    Ecomp = maximum(abs.(xhat-xtrue) ./ abs.(xtrue))
    if verbose
        println("Bounds: $Bnorm $Bcomp")
        println("Errors: $Enorm $Ecomp")
    end
    if Bnorm > 0.1
        @test κcomp > 100 * crit
    else
        γ = max(10,sqrt(n))
        @test Enorm < 1.1*Bnorm
        if κnorm < crit
            @test Bnorm < γ * eps(real(T))
        end
        @test Ecomp < 1.1*Bcomp
        if κcomp < crit
            @test Bcomp < γ * eps(real(T))
        end
    end
    return true
end

# pick log10(condition-number) for various cases
function lkval(class,T)
    if class == :easy
        if real(T) <: Float32
            return 5.0
        elseif real(T) <: Float64
            return 13.0
        elseif real(T) <: Float128
            return 29.0
        end
    elseif class == :moderate
        if real(T) <: Float32
            return 7.5
        elseif real(T) <: Float64
            return 16.0
        end
    elseif class == :painful
        if real(T) <: Float32
            return 9.0
        elseif real(T) <: Float64
            return 18.0
        end
    end
    throw(ArgumentError("undefined lkval"))
end

@testset "matrix rhs $T" for T in (Float32, Float64, ComplexF32, ComplexF64)
    for n in [10]
        A = mkmat(n,lkval(:easy,T),T)
        nrhs = 3
        X = rand(T,n,nrhs)
        B = copy(X)
        X1 = copy(X)
        bn1 = zeros(T,nrhs)
        bc1 = zeros(T,nrhs)
        # check validity w/ view arg (someday maybe more tricky AbstractArrays)
        runone(A,view(X,:,2))
        for j=1:nrhs
            x,bnorm,bcomp = @inferred(rfldiv(A,view(X,:,j)))
            X1[:,j] .= x
            bn1[j] = bnorm
            bc1[j] = bcomp
        end
        X2, bn2, bc2 = @inferred(rfldiv(A,B))
        @test X1 ≈ X2
        @test bn1 ≈ bn2
        @test bc1 ≈ bc2
    end
end

@testset "preprocessed args $T" for T in (Float32, Float128)
    n = 16
    A = mkmat(n,lkval(:easy,T),T)
    # make it badly scaled
    s = 1 / sqrt(floatmax(T))
    A = s * A
    x = rand(T,n)
    b = A * x
    # basic usage for comparison
    x1, bn1, bc1 = rfldiv(A,b)

    # example of use with precomputed factor
    Rv, Cv = equilibrators(A)
    R = Diagonal(Rv)
    As = R * A * Diagonal(Cv)
    bs = R * b
    F = lu(As)
    a = opnorm(As,Inf)
    κnorm = condInfest(As,F,a)
    x2, bn2, bc2 = rfldiv(As,bs; F=F, κ = κnorm, equilibrate = false)
    cx2 = Diagonal(Cv) * x2
    @test cx2 ≈ x1
    @test bn2 ≈ bn1
    @test bc2 ≈ bc1
    # make sure this was not an empty test
    x3, bn3, bc3 = @test_logs (:warn, r"no convergence.*") rfldiv(A,b; F=F, κ = κnorm, equilibrate = false)
    @test ! (x3 ≈ x1)
end

@testset "well-conditioned $T" for T in (Float32, Float64, ComplexF32, ComplexF64)
    for n in [10,30,100]
        A = mkmat(n,lkval(:easy,T),T)
        x = rand(n)
        runone(A,x)
    end
end

@testset "marginally-conditioned $T" for T in (Float32, Float64, ComplexF32, ComplexF64)
    for n in [10,30,100]
        A = mkmat(n,lkval(:moderate,T),T)
        x = rand(n)
        runone(A,x)
    end
end

@info "The next block of tests is expected to produce warnings"

@testset "badly-conditioned $T" for T in (Float32, Float64, ComplexF32, ComplexF64)
    # We don't test for convergence failure here because
    # the method occasionally works in this regime.
    for n in [10,30,100]
        LU_ok = false
        for j in 1:10
            A = mkmat(n,lkval(:painful,T),T)
            x = rand(n)
            LU_ok = runone(A,x)
            if LU_ok
                break
            end
        end
        if !LU_ok
            @info "failed to find nonsingular example for n=$n"
            @test_broken LU_ok
        end
    end
end

include("eigen.jl")


text:
using Porta

TS³(z₀, z₁) = begin
    (X₀, X₁)
end

αⱼ(X₀, X₁) = Xⱼ
α̅ⱼ(X₀, X₁) = X̅ⱼ

A(z₀, z₁) = 1/2 * (z̅₀ * α₀ - z₀ * α̅₀ + z̅₁ * α₁ - z₁ * α̅₁)


text:
module AnsiCodes

const CSI = "\e["

cursor_up!(io::IO, n=1) = print(io, CSI, n, 'A')
cursor_down!(io::IO, n=1) = print(io, CSI, n, 'B')

cursor_horizontal_absolute!(io::IO, n=1) = print(io, CSI, n, 'G')

@enum EraseLineMode begin
    CURSOR_TO_END = 0
    CURSOR_TO_BEGINNING = 1
    ENTIRE_LINE = 2
end
erase_line!(io::IO, mode::EraseLineMode=ENTIRE_LINE) = print(io, CSI, Int(mode), 'K')

end


text:
module MaxHelpingHandGoldenStrawberryCustomConditions

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/GoldenStrawberryCustomConditions" GoldenStrawberryCustomConditions(x::Integer, y::Integer,
    mustNotDieAndVisitFurtherRooms::Bool=true, mustHaveUnlockedCSides::Bool=true, mustHaveCompletedLevel::Bool=true, showGoldenChapterCard::Bool=true)

const placements = Ahorn.PlacementDict(
    "Golden Strawberry (Custom Conditions)\n(max480's Helping Hand)" => Ahorn.EntityPlacement(
        GoldenStrawberryCustomConditions
    )
)

sprite = "collectables/goldberry/idle00"

function Ahorn.selection(entity::GoldenStrawberryCustomConditions)
    x, y = Ahorn.position(entity)
    return Ahorn.getSpriteRectangle(sprite, x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::GoldenStrawberryCustomConditions, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

myFMU = fmi2Load("SpringPendulum1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])
myFMU.executionConfig.assertOnWarning = true

comp = fmi2Instantiate!(myFMU; loggingOn=false)
@test comp != 0

@test fmi2SetupExperiment(comp) == 0
@test fmi2EnterInitializationMode(comp) == 0
@test fmi2ExitInitializationMode(comp) == 0

numStates = length(myFMU.modelDescription.stateValueReferences)
targetValues = [[0.0, -10.0], [1.0, 0.0]]
dir_ders_buffer = zeros(fmi2Real, numStates)
sample_ders_buffer = zeros(fmi2Real, numStates, 1)
for i in 1:numStates

    if fmi2ProvidesDirectionalDerivative(myFMU)
        # multi derivatives calls
        sample_ders = fmi2SampleDirectionalDerivative(comp, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]])
        fmi2SampleDirectionalDerivative!(comp, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]], sample_ders_buffer)

        @test sum(abs.(sample_ders[:,1] - targetValues[i])) < 1e-3
        @test sum(abs.(sample_ders_buffer[:,1] - targetValues[i])) < 1e-3

        dir_ders = fmi2GetDirectionalDerivative(comp, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]])
        fmi2GetDirectionalDerivative!(comp, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]], dir_ders_buffer)

        @test sum(abs.(dir_ders - targetValues[i])) < 1e-3
        @test sum(abs.(dir_ders_buffer - targetValues[i])) < 1e-3

        # single derivative call
        dir_der = fmi2GetDirectionalDerivative(comp, myFMU.modelDescription.derivativeValueReferences[1], myFMU.modelDescription.stateValueReferences[1])
        @test dir_der == targetValues[1][1]
    else
        @warn "Skipping directional derivative testing, FMU from $(ENV["EXPORTINGTOOL"]) doesn't support directional derivatives."
    end
end

jac = fmi2GetJacobian(comp, myFMU.modelDescription.derivativeValueReferences, myFMU.modelDescription.stateValueReferences)
@test jac ≈ hcat(targetValues...)

jac = fmi2SampleDirectionalDerivative(comp, myFMU.modelDescription.derivativeValueReferences, myFMU.modelDescription.stateValueReferences)
@test jac ≈ hcat(targetValues...)

fmi2Unload(myFMU)


text:
"""
    ControlBenchmarkInfo

A struct to hold various information about a control benchmark problem, such as its name,
id, description and a set of keywords for locating benchmarks problems.
"""
struct ControlBenchmarkInfo
    name::String             # The name of the benchmark
    id::String               # A concise alpha-numeric ID for the benchmark
    desc::String             # A short description of the benchmark
    keywords::Set{Symbol}    # A set of keywords for the benchmark
    configurable::Bool       # True if the problem type contains configurable parameters
end

"""
    benchmarkmap

A dictionary mapping the benchmark ID to the type implementing the benchmark. Benchmarks are added
to this dictionary using the [`@addbenchmark`](@ref) macro.
"""
benchmarkmap = Dict{String,DataType}()

"""
    benchmarkinfo

A dictionary mapping the benchmark type to the struct containing its information. Benchmarks are added
to this dictionary using the [`@addbenchmark`](@ref) macro.
"""
benchmarkinfo = Dict{DataType, ControlBenchmarkInfo}()

"""
    @addbenchmark( T, id::String, name::String, desc::String, keywords )

Add a benchmark to the main benchmarks list.

`T` is the type for the benchmark. `id` is an alpha-numeric string used for identifying the benchmark
and easily locating it using either [`controlbenchmark`](@ref) or [`getbenchmark`](@ref).
`name` and `description` are short strings naming and describing the problem.
`keywords` is a set containing symbols that describe the system and are used by the
[`listbenchmarks`](@ref) function.

# Example
```
struct TestSystem <: ControlBenchmarkProblem end

@addbenchmark TestSystem "test01" "Sample system" "A sample system" Set([:sample, :linear])
```
"""
macro addbenchmark( T, id::String, name::String, desc::String, keywords )
    quote
        benchmarkmap[$(esc(id))] = $(esc(T))
        obj = $(T)()
        benchmarkinfo[$(esc(T))] = ControlBenchmarkInfo( $(esc(name)), $(esc(id)), $(esc(desc)), $(esc(keywords)), !isimmutable(obj) )
    end
end

"""
    @checkbounds( T, params )

Performs a bounds check on all the elements of the struct `params`, which is of type `T`. The bounds
are given inside a dictionary returned by the function `parameterbounds(::T)`, where the dictionary keys
must be the fieldnames. Every field must be included in the dictionary for this macro to be used. The comparison
is done using the `∈` operator, so the dictionary must return a type which implements that operator (such as [`Bound`](@ref))
"""
macro checkbounds( T, params )
    block = Expr(:block)
    push!(block.args, :(b = parameterbounds( $(esc(params)) ) ) )

    for f in eval( :(fieldnames( $(T) ) ) )
        str = String(f)
        push!(block.args, :($(esc(params)).$(f) ∈ b[$(str)]  || throw( DomainError( $(esc(params)).$(f), "$($(str)) is outside the allowed domain" ) ) ) )
    end

    return block
end


text:
using OrdinaryDiffEq, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear
@testset "Composite Algorithm Tests" begin
  prob = prob_ode_2Dlinear
  choice_function(integrator) = (Int(integrator.t<0.5) + 1)
  alg_double = CompositeAlgorithm((Tsit5(),Tsit5()),choice_function)
  alg_double2 = CompositeAlgorithm((Vern6(),Vern6()),choice_function)
  alg_switch = CompositeAlgorithm((Tsit5(),Vern7()),choice_function)

  @time sol1 = solve(prob_ode_linear,alg_double)
  @time sol2 = solve(prob_ode_linear,Tsit5())
  @test sol1.t == sol2.t
  @test sol1(0.8) == sol2(0.8)

  integrator1 = init(prob,alg_double2)
  integrator2 = init(prob,Vern6())
  solve!(integrator1)
  solve!(integrator2)

  @test integrator1.sol.t == integrator2.sol.t

  sol = solve(prob,alg_switch)
end


text:
input_1 = """
          30 1 0 120
          50 1 150 90
          60 1 60 90
          70 1 180 180
          90 1 180 60
          """
input_2 = """
        30 1 350 350
        180 1 10 100
        35 1 67 23
        31 1 30 4
        32 1 50 7
        """
input = input_1
type Disk
  speed::Int
  openings::Set{Int}
  function Disk(s::Int, o::Vector{(Int,Int)})
    opens = Set{Int}()
    for (start,length) in o
      for i in start:(start + length)
        push!(opens,i % 360)
      end
    end
    new(s,opens)
  end
end

function first_range!(numbers)
  start = numbers[1]
  stop = shift!(numbers)
  while !isempty(numbers) && stop + 1 == numbers[1]
    stop = shift!(numbers)
  end
  return start:stop
end

function Base.show(io::IO, disk::Disk)
  print(io,"Disk: $(disk.speed)\n\t")
  numbers = sort!([o for o in disk.openings])
  ranges = Any[]
  while !isempty(numbers)
   push!(ranges, first_range!(numbers))
  end
  println(io,ranges)
end

lines = split(input,"\n")
disks = Disk[]
for line in lines
  if line == ""
    continue
  end
  nums = [int(n) for n in split(line)]
  speed = shift!(nums)
  num_openings = shift!(nums)
  openings = (Int,Int)[]
  for i in 1:num_openings
    push!(openings, (shift!(nums),shift!(nums)))
  end
  push!(disks, Disk(speed, openings))
end

# Each argument is vector of Set{Int}, one Set per disk.
function allows_light(ds...)
  !isempty(intersect(ds...))
end

# Advance the disk by one second.
function advance!(disk::Disk)
  next_set = Set{Int}()
  for n in disk.openings
    m = n + disk.speed # advance by one second
    push!(next_set, m % 360)
  end
  disk.openings = next_set
end

found_light = -1
for i = 1:(360*5)
  if allows_light([disk.openings for disk in disks]...)
    found_light = i - 1
    break
  end
  for disk in disks
    advance!(disk)
  end
end

if found_light == -1
  println("none")
else
  println(found_light)
end

max_stars_repo_name: astrieanna/usaco

max_issues_repo_name: astrieanna/usaco

max_forks_repo_name: astrieanna/usaco


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(ROT),GVar{GVar{Float64, Float64}, GVar{Float64, Float64}},GVar{GVar{Float64, Float64}, GVar{Float64, Float64}},GVar{GVar{Float64, Float64}, GVar{Float64, Float64}}})   # time: 0.08222882
    Base.precompile(Tuple{MinusEq{typeof(complex)},GVar{Float64, Float64},Union{GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}},Union{GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}}})   # time: 0.06912181
    Base.precompile(Tuple{typeof(wrap_bcastgrad),Type{AutoBcast{Int64, 4}},Vector{Float64}})   # time: 0.05026503
    Base.precompile(Tuple{Type{GVar},Matrix{Float64}})   # time: 0.044000912
    Base.precompile(Tuple{MinusEq{typeof(sincos)},Tuple{GVar{Float64, Float64}, GVar{Float64, Float64}},GVar{Float64, Float64}})   # time: 0.042384867
    Base.precompile(Tuple{MinusEq{typeof(identity)},GVar{Float64, Float64},Any})   # time: 0.03920624
    Base.precompile(Tuple{typeof(NiLangCore.loaddata),Vector{Any},Vector{Any}})   # time: 0.03768901
    Base.precompile(Tuple{typeof(grad),Tuple{PlusEq{typeof(-)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.036892496
    Base.precompile(Tuple{typeof(NiLangCore.loaddata),Type{Vector{Float64}},Vector{Float64}})   # time: 0.035732098
    Base.precompile(Tuple{MinusEq{typeof(log)},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.03535623
    Base.precompile(Tuple{PlusEq{typeof(complex)},GVar{Float64, Float64},Union{GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}},Union{GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}}})   # time: 0.03526533
    Base.precompile(Tuple{typeof(grad),NTuple{4, GVar{Float64, Float64}}})   # time: 0.032369606
    Base.precompile(Tuple{typeof(wrap_jacobian),Type{AutoBcast{Float64, 1}},Vector{Float64}})   # time: 0.032170866
    Base.precompile(Tuple{Type{GVar},Vector{Float64}})   # time: 0.029704165
    Base.precompile(Tuple{typeof(IROT),GVar{Float64, Float64},GVar{Float64, Float64},GVar{Float64, Float64}})   # time: 0.025739137
    Base.precompile(Tuple{PlusEq{typeof(^)},GVar{GVar{Float64, Float64}, GVar{Float64, Float64}},GVar{GVar{Float64, Float64}, GVar{Float64, Float64}},GVar{GVar{Float64, Float64}, GVar{Float64, Float64}}})   # time: 0.024799693
    Base.precompile(Tuple{MinusEq{typeof(identity)},GVar{Float64, Float64},GVar{Float64, Float64}})   # time: 0.023590498
    Base.precompile(Tuple{typeof(grad),Vector{GVar{Float64, Float64}}})   # time: 0.02317364
    Base.precompile(Tuple{typeof(grad),Tuple{Vector{GVar{Float64, Float64}}, Vector{GVar{Float64, Float64}}, Vector{GVar{Float64, Float64}}, GVar{Float64, Float64}}})   # time: 0.022934826
    Base.precompile(Tuple{PlusEq{typeof(^)},GVar{GVar{Float64, Float64}, GVar{Float64, Float64}},GVar{GVar{Float64, Float64}, GVar{Float64, Float64}},Int64})   # time: 0.022084313
    Base.precompile(Tuple{typeof(NiLangCore.loaddata),Vector{GVar{Float64, Float64}},Vector{Float64}})   # time: 0.02116515
    Base.precompile(Tuple{MinusEq{typeof(/)},GVar{Float64, AutoBcast{Float64, 3}},GVar{Float64, AutoBcast{Float64, 3}},GVar{Float64, AutoBcast{Float64, 3}}})   # time: 0.020155573
    Base.precompile(Tuple{PlusEq{typeof(*)},GVar{Float64, Float64},Any,Float64})   # time: 0.016616046
    Base.precompile(Tuple{typeof(IROT),GVar{Float64, AutoBcast{Float64, 3}},GVar{Float64, AutoBcast{Float64, 3}},GVar{Float64, AutoBcast{Float64, 3}}})   # time: 0.016194412
    Base.precompile(Tuple{Type{GVar},Tuple{Float64, ComplexF64},Tuple{ComplexF64, ComplexF64}})   # time: 0.01619109
    Base.precompile(Tuple{PlusEq{typeof(^)},GVar{Float64, Float64},Any,Int64})   # time: 0.015741551
    Base.precompile(Tuple{Type{GVar},AbstractArray,AbstractArray})   # time: 0.013830282
    Base.precompile(Tuple{typeof(grad),Tuple{MinusEq{typeof(inv)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.011192932
    Base.precompile(Tuple{typeof(wrap_jacobian),Type{AutoBcast{Float64, 3}},Vector{Float64}})   # time: 0.011075067
    Base.precompile(Tuple{typeof(wrap_jacobian),Type{AutoBcast{Float64, 2}},Vector{Float64}})   # time: 0.011053615
    Base.precompile(Tuple{typeof(grad),Tuple{MinusEq{typeof(^)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}, GVar{Float64, Float64}}})   # time: 0.010716249
    Base.precompile(Tuple{typeof(grad),Tuple{MinusEq{typeof(complex)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, GVar{Float64, Float64}, GVar{Float64, Float64}}})   # time: 0.010689858
    Base.precompile(Tuple{typeof(grad),Tuple{MinusEq{typeof(/)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.010517918
    Base.precompile(Tuple{typeof(grad),Tuple{MinusEq{typeof(*)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.01049392
    Base.precompile(Tuple{typeof(grad),Tuple{MinusEq{typeof(-)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.010452602
    Base.precompile(Tuple{typeof(grad),Tuple{PlusEq{typeof(*)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.01041573
    Base.precompile(Tuple{typeof(grad),Tuple{PlusEq{typeof(^)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}, GVar{Float64, Float64}}})   # time: 0.010397522
    Base.precompile(Tuple{typeof(grad),Tuple{MinusEq{typeof(+)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.010353689
    Base.precompile(Tuple{typeof(grad),Tuple{PlusEq{typeof(complex)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, GVar{Float64, Float64}, GVar{Float64, Float64}}})   # time: 0.010343891
    Base.precompile(Tuple{typeof(grad),Tuple{PlusEq{typeof(/)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.010335974
    Base.precompile(Tuple{typeof(grad),Tuple{PlusEq{typeof(exp)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.010250594
    Base.precompile(Tuple{typeof(grad),Tuple{PlusEq{typeof(+)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.010101405
    Base.precompile(Tuple{typeof(grad),Tuple{MinusEq{typeof(log)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.009795744
    Base.precompile(Tuple{typeof(grad),Tuple{MinusEq{typeof(identity)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.009747697
    Base.precompile(Tuple{typeof(grad),Tuple{MinusEq{typeof(exp)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.009696897
    Base.precompile(Tuple{typeof(grad),Tuple{MinusEq{typeof(^)}, Union{GVar{Float64, Float64}, Float64}, Vararg{Any, N} where N}})   # time: 0.009630851
    Base.precompile(Tuple{typeof(grad),Tuple{PlusEq{typeof(inv)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.009559213
    Base.precompile(Tuple{typeof(grad),Tuple{PlusEq{typeof(log)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.009503618
    Base.precompile(Tuple{MinusEq{typeof(identity)},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.009323083
    Base.precompile(Tuple{typeof(grad),Tuple{PlusEq{typeof(identity)}, GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}, Complex{GVar{Float64, Float64}}}})   # time: 0.0093146
    Base.precompile(Tuple{MinusEq{typeof(cos)},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.009097533
    Base.precompile(Tuple{typeof(grad),Tuple{PlusEq{typeof(^)}, Union{GVar{Float64, Float64}, Float64}, Vararg{Any, N} where N}})   # time: 0.009068273
    Base.precompile(Tuple{typeof(grad),Tuple{MinusEq{typeof(complex)}, Union{GVar{Float64, Float64}, Float64}, Complex{GVar{Float64, Float64}}, GVar{Float64, Float64}, GVar{Float64, Float64}}})   # time: 0.008931454
    Base.precompile(Tuple{PlusEq{typeof(^)},GVar{Float64, Float64},GVar,Int64})   # time: 0.008520359
    Base.precompile(Tuple{typeof(grad),Tuple{PlusEq{typeof(complex)}, Union{GVar{Float64, Float64}, Float64}, Complex{GVar{Float64, Float64}}, GVar{Float64, Float64}, GVar{Float64, Float64}}})   # time: 0.00843281
    Base.precompile(Tuple{typeof(wrap_bcastgrad),Type{AutoBcast{Int64, 4}},Tuple{Float64}})   # time: 0.007901799
    Base.precompile(Tuple{typeof(grad),Tuple{typeof(NEG), GVar{Float64, Float64}, Complex{GVar{Float64, Float64}}}})   # time: 0.007657704
    Base.precompile(Tuple{PlusEq{typeof(/)},GVar{Float64, Float64},GVar{Float64, Float64},GVar{Float64, Float64}})   # time: 0.007580364
    Base.precompile(Tuple{typeof(grad),Tuple{GVar{Float64, Float64}, Matrix{GVar{Float64, Float64}}}})   # time: 0.007370161
    Base.precompile(Tuple{typeof(grad),Tuple{GVar{Float64, Float64}, GVar{Float64, Float64}, Vector{GVar{Float64, Float64}}}})   # time: 0.007349056
    Base.precompile(Tuple{MinusEq{typeof(*)},GVar{Float64, Float64},GVar{Float64, Float64},GVar{Float64, Float64}})   # time: 0.007118214
    Base.precompile(Tuple{PlusEq{typeof(identity)},Vector{GVar{Float64, Float64}},Vector{GVar{Float64, Float64}}})   # time: 0.007062924
    Base.precompile(Tuple{MinusEq{typeof(^)},GVar{Float64, AutoBcast{Float64, 3}},GVar{Float64, AutoBcast{Float64, 3}},GVar{Float64, AutoBcast{Float64, 3}}})   # time: 0.007059243
    Base.precompile(Tuple{typeof(grad),Tuple{Float64, GVar{Float64, Float64}}})   # time: 0.006490685
    Base.precompile(Tuple{MinusEq{typeof(identity)},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.005542749
    Base.precompile(Tuple{MinusEq{typeof(exp)},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.004707843
    Base.precompile(Tuple{typeof(zero),GVar{Int64, AutoBcast{Int64, 3}}})   # time: 0.004680167
    Base.precompile(Tuple{typeof(grad),Tuple{GVar{Float64, Float64}, GVar{Float64, Float64}, GVar{Float64, Float64}}})   # time: 0.00466041
    Base.precompile(Tuple{MinusEq{typeof(abs)},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.004496061
    Base.precompile(Tuple{MinusEq{typeof(log)},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.003723256
    Base.precompile(Tuple{PlusEq{typeof(/)},GVar{Float64, Float64},GVar,GVar{Float64, Float64}})   # time: 0.003611694
    Base.precompile(Tuple{typeof(grad),Tuple{GVar{Float64, Float64}, GVar{Float64, Float64}}})   # time: 0.003529446
    Base.precompile(Tuple{MinusEq{typeof(/)},GVar{Float64, Float64},GVar{Float64, Float64},GVar{Float64, Float64}})   # time: 0.003528577
    Base.precompile(Tuple{MinusEq{typeof(tan)},GVar{Float64, Float64},GVar{Float64, Float64}})   # time: 0.003438413
    Base.precompile(Tuple{MinusEq{typeof(identity)},GVar{Float64, Float64},GVar})   # time: 0.003361321
    Base.precompile(Tuple{typeof(grad),Tuple{Float64, Float64}})   # time: 0.003360894
    Base.precompile(Tuple{Type{Inv{GVar}},AbstractArray})   # time: 0.003287998
    Base.precompile(Tuple{typeof(SWAP),GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.003202025
    Base.precompile(Tuple{MinusEq{typeof(*)},GVar{Float64, AutoBcast{Float64, 3}},GVar{Float64, AutoBcast{Float64, 3}},GVar{Float64, AutoBcast{Float64, 3}}})   # time: 0.003030533
    Base.precompile(Tuple{MinusEq{typeof(*)},GVar{Float64, AutoBcast{Float64, 3}},GVar{Float64, AutoBcast{Float64, 3}}})   # time: 0.002981777
    Base.precompile(Tuple{MinusEq{typeof(/)},GVar{Float64, Float64},GVar{Float64, Float64},Int64})   # time: 0.002951113
    Base.precompile(Tuple{MinusEq{typeof(log)},GVar{Float64, Float64},GVar{Float64, Float64},GVar{Float64, Float64}})   # time: 0.002914175
    Base.precompile(Tuple{Type{GVar},Float64,AutoBcast{Float64, _A} where _A})   # time: 0.002864071
    Base.precompile(Tuple{PlusEq{typeof(*)},GVar{Float64, Float64},GVar,Float64})   # time: 0.00271298
    Base.precompile(Tuple{PlusEq{typeof(/)},GVar{Float64, Float64},Real,GVar{Float64, Float64}})   # time: 0.002609891
    Base.precompile(Tuple{typeof(NiLangCore.loaddata),Type{GVar{Float64, AutoBcast{Float64, 3}}},Float64})   # time: 0.002460321
    Base.precompile(Tuple{typeof(-),GVar{Float64, AutoBcast{Float64, 1}}})   # time: 0.002408763
    Base.precompile(Tuple{typeof(NEG),GVar{Float64, AutoBcast{Float64, 1}}})   # time: 0.00236964
    Base.precompile(Tuple{PlusEq{typeof(^)},GVar{Float64, Float64},Float64,Int64})   # time: 0.002257135
    Base.precompile(Tuple{PlusEq{typeof(identity)},GVar{Float64, Float64},GVar{Float64, Float64}})   # time: 0.002224634
    Base.precompile(Tuple{PlusEq{typeof(*)},GVar{Float64, Float64},Real,Float64})   # time: 0.002007978
    Base.precompile(Tuple{SubConst{Float64},GVar{Float64, Float64}})   # time: 0.001776433
    Base.precompile(Tuple{PlusEq{typeof(^)},GVar{Float64, Float64},Real,Int64})   # time: 0.00169386
    Base.precompile(Tuple{MinusEq{typeof(*)},GVar{Float64, Float64},Float64,GVar{Float64, Float64}})   # time: 0.001672073
    Base.precompile(Tuple{MinusEq{typeof(*)},GVar{Float64, Float64},Int64,GVar{Float64, Float64}})   # time: 0.001642821
    Base.precompile(Tuple{PlusEq{typeof(*)},GVar{Float64, Float64},Float64,GVar{Float64, Float64}})   # time: 0.001626397
    Base.precompile(Tuple{AddConst{Float64},GVar{Float64, Float64}})   # time: 0.001618372
    Base.precompile(Tuple{PlusEq{typeof(*)},GVar{Float64, Float64},Int64,GVar{Float64, Float64}})   # time: 0.001491792
    Base.precompile(Tuple{PlusEq{typeof(*)},GVar{Float64, Float64},GVar{Float64, Float64},Float64})   # time: 0.001471461
    Base.precompile(Tuple{MinusEq{typeof(*)},GVar{Float64, Float64},GVar{Float64, Float64},Float64})   # time: 0.001428104
    Base.precompile(Tuple{Type{GVar},ComplexF64})   # time: 0.001348701
    Base.precompile(Tuple{MinusEq{typeof(cos)},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.001196095
    Base.precompile(Tuple{MinusEq{typeof(sin)},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.001177611
    Base.precompile(Tuple{MinusEq{typeof(identity)},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.001173775
    Base.precompile(Tuple{MinusEq{typeof(abs)},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.001170435
    Base.precompile(Tuple{typeof(chfield),Complex{GVar{Float64, Float64}},typeof(grad),ComplexF64})   # time: 0.00116481
    Base.precompile(Tuple{MinusEq{typeof(abs)},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.001158102
    Base.precompile(Tuple{MinusEq{typeof(exp)},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.001156032
    Base.precompile(Tuple{MinusEq{typeof(*)},GVar{Float64, Float64},GVar{Float64, Float64}})   # time: 0.001134847
    Base.precompile(Tuple{MinusEq{typeof(sin)},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.001125385
    Base.precompile(Tuple{MinusEq{typeof(cos)},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.001105154
    Base.precompile(Tuple{MinusEq{typeof(exp)},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.001097152
    Base.precompile(Tuple{MinusEq{typeof(log)},GVar{Float64, AutoBcast{Float64, 2}}})   # time: 0.001096638
    Base.precompile(Tuple{MinusEq{typeof(/)},GVar{Float64, AutoBcast{Float64, 3}}})   # time: 0.001095843
    Base.precompile(Tuple{MinusEq{typeof(^)},GVar{Float64, AutoBcast{Float64, 3}}})   # time: 0.001095656
    Base.precompile(Tuple{typeof(chfield),GVar{Float64, Float64},typeof(grad),Float64})   # time: 0.001003963
    Base.precompile(Tuple{MinusEq{typeof(identity)},GVar{Float64, Float64},Real})   # time: 0.00100021
end


text:
maxwell_lep(sim::Simulation{1};ky::Number=0,kz::Number=0) = Maxwell_LEP(sim;ky=ky,kz=kz)

function Maxwell_LEP(
            sim::Simulation{1};
            ky::Number=0,
            kz::Number=0)

    αε = _compute_αε(sim)
    return Maxwell_LEP(sim,
            _compute_Σs(sim,ky,kz),
            sim.curlcurl(ky,kz),
            αε,
            spdiagm(0=>zeros(ComplexF64,size(αε,1))),
            0,ky,kz,0,0,0)
end
function (ml::Maxwell_LEP{1})(ω::Number)
    sim = ml.sim
    n::Int = length(sim.x)
    ka, ky, kz = ml.ka, ml.ky, ml.kz
    f = ml.sim.self_energy.f
    Σ01,Σ02,Σ11,Σ12,Σ2 = ml.Σs

    Σ0 = Σ01*f[1](ω,ka,ky,kz) + Σ02*f[2](ω,ka,ky,kz)
    Σ1 = Σ11*f[1](ω,ka,ky,kz) + Σ12*f[2](ω,ka,ky,kz)

    χs = map(d->d.χ,sim.domains)
    Fs = map(d->d.pump,sim.domains)
    χχ = map(x->susceptability(x,ω),χs)
    foreach((χ,d)->begin
                    for i ∈ eachindex(ml.αχ.nzval)
                        j = mod1(i,n)
                        if d==ml.sim.domain_indices[j]
                            ml.αχ.nzval[i]=0
                            for x ∈ χ
                                ml.αχ.nzval[i] += Fs[d](ml.sim.x[j])*ml.sim.α[1][j]*x
                            end
                        end
                    end
                end, χχ, eachindex(χχ))
    return Σ0 + Σ1 + Σ2 + ml.curlcurl, ml.αε + ml.αχ
end
function (ml::Maxwell_LEP{1})(ω::Number,ωs::Array,ψs::Array)
    sim = ml.sim
    n::Int = length(sim.x)
    ka, ky, kz = ml.ka, ml.ky, ml.kz
    f = ml.sim.self_energy.f
    Σ01,Σ02,Σ11,Σ12,Σ2 = ml.Σs

    Σ0 = Σ01*f[1](ω,ka,ky,kz) + Σ02*f[2](ω,ka,ky,kz)
    Σ1 = Σ11*f[1](ω,ka,ky,kz) + Σ12*f[2](ω,ka,ky,kz)

    χs = map(d->d.χ,sim.domains)
    Fs = map(d->d.pump,sim.domains)
    foreach(x->susceptability(x,ω,ωs,ψs),χs)
    foreach((χ,d)->begin
                    for i ∈ eachindex(ml.αχ.nzval)
                        j = mod1(i,n)
                        if d==ml.sim.domain_indices[j]
                            ml.αχ.nzval[i]=0
                            for x ∈ χ
                                ml.αχ.nzval[i] += Fs[d](ml.sim.x[j])*ml.sim.α[1][j]*x.chi[i]
                            end
                        end
                    end
                end, χs, eachindex(χs))
    return Σ0 + Σ1 + Σ2 + ml.curlcurl, ml.αε + ml.αχ
end

maxwell(sim::Simulation{1};ky::Number=0,kz::Number=0) = Maxwell(maxwell_lep(sim))
function (m::Maxwell{1})(ω::Number,args...)
    A,B =m.lep(ω,args...)
    return A-B*ω^2
end

function maxwell_nep(
            sim::Simulation{1,C,T};
            ky::Number = 0,
            kz::Number = 0,
            FType::DataType = T,
            kwargs...
            ) where {C,T}

    ka = 0
    f = map(f->(ω->f(ω,ka,ky,kz)),sim.self_energy.f)
    Σ01,Σ02,Σ11,Σ12,Σ2 = _compute_Σs(sim,ky,kz)
    αε = _compute_αε(sim)
    αF, fχ = _compute_αχ(sim)
    As = vcat([Σ01,Σ02,Σ11,Σ12,Σ2,sim.curlcurl(ky,kz),αε],αF)
    fs = [f[1],f[2],f[1],f[2],one,one,ω->-ω^2,map(ϝ->(ω->-ω^2*ϝ(ω)),fχ)...]
    return SPMF_NEP(As,fs;kwargs...)
end

function _compute_Σs(sim::Simulation{1},ky::Number,kz::Number)

    ky², kz², kykz = ky^2, kz^2, ky*kz

    Σ0 = sim.self_energy.Σ0
    Σ1 = sim.self_energy.Σ1
    Σ2 = sim.self_energy.Σ2

    I1 = sparse([2,3],[2,3],[1,1],3,3)
    I2 = sparse([2,3,1,1],[1,1,2,3],[ky,kz,ky,kz],3,3)

    Σ01 = kron(I1,Σ0[1])
    Σ02 = kron(I1,Σ0[2])
    Σ11 = kron(I2,Σ1[1])
    Σ12 = kron(I2,Σ1[2])
    Σ2 = kron(sparse([1,2,3,2,3],[1,2,2,3,3],[ky²+kz²,kz²,kykz,kykz,ky²],3,3),Σ2)
    return Σ01,Σ02,Σ11,Σ12,Σ2
end

_compute_αε(sim::Simulation{1})::SparseMatrixCSC{ComplexF64,Int} =
    kron(I3,spdiagm(0=>sim.α[1].*sim.ε))

function _compute_αχ(sim::Simulation{1})
    χs = flatten(map(d->d.χ,sim.domains),AbstractDispersion)
    domain_inds = flatten(map(d->ntuple(identity,length(d.χ)),sim.domains))
    fχs = map(x->(ω->susceptability(x,ω)),χs)

    n = length(sim.x)
    Fs = map(d->sim.domains[d].pump,domain_inds)
    αFs = Array{SparseMatrixCSC{ComplexF64,Int},1}(undef,length(domain_inds))
    for d ∈ eachindex(αFs)
        αFs[d] = spdiagm(0=>zeros(ComplexF64,3n))
        for i ∈ 1:3n
            j = mod1(i,n)
            if domain_inds[d]==sim.domain_indices[j]
                αFs[d].nzval[i] = Fs[d](sim.x[j])*sim.α[1][j]
            end
        end
    end
    return αFs,fχs
end


text:
# Create a default kernel from two parameters k[1] and k[2]
make_kernel(k) = softplus(k[1]) * (SqExponentialKernel() ∘ ScaleTransform(softplus(k[2])))

# Computes the optimal closed form solution for the variational posterior
# q(u) (e.g. # https://krasserm.github.io/2020/12/12/gaussian-processes-sparse/
# equations (11) & (12)). Assumes a ZeroMean function.
function optimal_variational_posterior(fu, fx, y)
    fu.f.mean isa AbstractGPs.ZeroMean ||
        error("The exact posterior requires a GP with ZeroMean.")
    σ² = fx.Σy[1]
    Kuf = cov(fu, fx)
    Kuu = Symmetric(cov(fu))
    Σ = (Symmetric(cov(fu) + (1 / σ²) * Kuf * Kuf'))
    m = ((1 / σ²) * Kuu * (Σ \ Kuf)) * y
    S = Symmetric(Kuu * (Σ \ Kuu))
    return MvNormal(m, S)
end


text:
using Statistics
using LinearAlgebra
using Flux
using LaTeXStrings
using KrylovMethods
using Plots
pyplot()

include("../src/multigrid/helmholtz_methods.jl");

# Grid

n = 256
n_cells = [n,n];
h = 1.0./n;
src = [128,128]

b = zeros(ComplexF64,n-1,n-1)
b[src[1],src[2]] = 1.0 ./mean(h.^2);

# Parameters

f = 20.0
kappa = ones(Float64,tuple(n_cells .-1...))
omega = 2*pi*f;

# Gamma

gamma_val = 0.00001
gamma = gamma_val*2*pi * ones(Float64,size(kappa));
pad_cells = [20;20]
gamma = absorbing_layer!(gamma, pad_cells, omega);

# V cycle

x = x0 = zeros(ComplexF64,n-1,n-1)
x = fgmres_v_cycle_helmholtz!(n, h, b, kappa, omega, gamma; restrt=30, maxIter=1)

v1_iter = 1
v2_iter = 20
iterations = 1
use_gmres = 1
level = 3

if use_gmres == 1
    # GMRES V cycle

    iterations = 1
    restrt = 30
    max_iter = 5
    shifted_laplacian_matrix, helmholtz_matrix = get_helmholtz_matrices!(kappa, omega, gamma; alpha=0.5)
    A(v) = vec(helmholtz_chain!(reshape(v, n-1, n-1, 1, 1), helmholtz_matrix; h=h))
    function M(v)
        v = reshape(v, n-1, n-1)
        x = zeros(ComplexF64,n-1,n-1)
        for i = 1:iterations
            x, = v_cycle_helmholtz!(n, h, x, v, kappa, omega, gamma; u=1,
                    v1_iter = v1_iter, v2_iter = v2_iter, alpha=0.5, log = 0, level = level)
        end
        return vec(x)
    end
    x,flag,err,iter,resvec = KrylovMethods.fgmres(A, vec(b), restrt, tol=1e-10, maxIter=max_iter,
                                                    M=M, x=vec(x), out=2, flexible=true)

    residual = resvec
    x = reshape(x,n-1,n-1)
else
    # V cycle
    v_cycle_iter = 100
    residual = zeros(v_cycle_iter)

    res = zeros(n-1,n-1)
    bnorm = norm(b)
    for i = 1:v_cycle_iter
        global x, helmholtz_matrix = v_cycle_helmholtz!(n, h, x, b, kappa, omega, gamma; u = 1,
                                        v1_iter = v1_iter, v2_iter = v2_iter, alpha=0.5, log = 0, level = level)
        global res = helmholtz_chain!(reshape(x, n-1, n-1, 1, 1), helmholtz_matrix; h=h)
        residual[i] = norm(b .- res[:,:,1,1]) / bnorm
        println(residual[i])
    end

end

q=M(vec(b))

heatmap(reshape(real(q),n-1,n-1), color=:grays)
print(residual)

iter = range(1, length=length(residual), iterations*length(residual))
p = plot(iter,residual,label="V cycle")
yaxis!(L"\Vert b - Hx \Vert_2", :log10)
xlabel!("Iterations")
savefig(replace("test/multigrid_helmholtz/results/f=$(f) gamma_val=$(gamma_val) residual graph","."=>"_"))

heatmap(reshape(real(x), n-1, n-1), color=:jet,size=(240,160))
savefig(replace("test/multigrid_helmholtz/results/without gamma n=$(n) f=$(f) gamma_val=$(gamma_val) result","."=>"_"))

show_last_residual = 0
if show_last_residual == 1
    last_res = helmholtz_chain!(reshape(x, n-1, n-1, 1, 1), helmholtz_matrix; h=h)
    heatmap(real(last_res[:,:,1,1] - b), color=:grays)
end


text:
using ODLCommonTools
using SummationByParts
using PumiInterface
using PdePumiInterface

dmg_name =  ".null"
smb_name = "cube8f.smb"
order = 4
shape_type = 1

sbp = TetSBP{Float64}(degree=order)

mesh = PumiMesh3{Float64}(dmg_name, smb_name, order, sbp, dofpernode=4, shape_type=shape_type)

#=
for i=1:mesh.numBoundaryFaces
  println("Element ", mesh.boundary_nums[i,1] -1, " face ", mesh.boundary_nums[i,2])
end

el = mesh.elements[1]
downwards, numdown = apf.getDownward(mesh.m_ptr, el, 2)
face = downwards[4]

#=
for i=1:numdown
  numup = apf.countAdjacent(mesh.m_ptr, downwards[i], 3)
  println("i = ", i, " numup = ", numup)
end
=#

which, flip, rotate  = apf.getAlignment(mesh.m_ptr, el, face)
=#


text:
function fpi(f, x0, k)
    x = zeros(k+1)
    x[1] = x0
    for i = 1:k
        x[i+1] = f(x[i])
    end

    return x[k+1]
end


text:
import MittagLefflerFunctions.MLF.mlf3
import MittagLefflerFunctions.MLF.radius
import MittagLefflerFunctions.Γ

α = 0.5
β = 2.0
z = 8.7 * exp(Complex(0.0, 3π/4))
p = 10
@assert (   (  abs(angle(z)) > π*α )
         && ( abs(abs(angle(z)) - (π*α)) > 10.0^(-p) ) )

r0 = radius(α, β, z, p)
val = mlf3(α, β, z, p, r0)
ref = ( erfcx(-z) - 1.0 - z / Γ(3/2) ) / z^2

@test abs(val-ref) ≤ abs(ref)


text:
struct MarsBarycenter <: abstractBarycenter end

naif_id(::Type{MarsBarycenter}) = 4
parent(::Type{MarsBarycenter}) = SSB


text:
@testcase :DICT_MERGE begin
	mydict1 = ["asdfasdf"=>1]
	mydict2 = ["asdfasdf"=>:asdfasdf]
	merge(mydict1, mydict2)
end

@testcase :STR_CONCAT begin
	asdfasdf = "asdfasdf"
	qwerty = "qwerty"
	asdfasdf + qwerty
end

@testcase :REGEX_ITER begin
	dataset = {1, 2, 3, 4, 5}
	query = "12345"
	for m in eachmatch(r"\d", query)
		dataset[m]
	end
end

@testcase :REGEX_ITER begin
	dataset = {"a"=>1, "b"=>2}
	query = "ab"
	for m in eachmatch(r"[a-z]", query)
		dataset[m]
	end
end

@testcase :INVARIANTS begin
	function alpha(asdfwef::Dict{ASCIIString, Number})
		asdfwef
	end

	alpha(["test"=>255])
end

@testcase :NOT_PYTHON begin
	eval(parse("for i in 1:12"))
end

@testcase :NOT_PYTHON begin
	eval(parse("for i in 1:12; identity(i); end; end"))
end

@testcase :BASE_ANDOR begin
	eval(parse("function Base.&&(a,b); a&&b; end"))
end

@testcase :BASE_ANDOR begin
	eval(parse("function Base.||(a,b); a||b; end"))
end

@testcase :ARRAY_SIZE begin
	function jwngwoijenfoer(asdfwef::Array{Float32,8})
		asdfwef
	end

	jwngwoijenfoer(Float32[1:8])
end

@testcase :FUNC_COLON begin
	eval(parse("function ooewjdbwoedjfvobwe11(): nkjsndf=[1:8];end; ooewjdbwoedjfvobwe11()"))
end

@testcase :FUNC_COLON begin
	eval(parse("function ooewjdbwoedjfvobwe12(): identity(9);end; ooewjdbwoedjfvobwe12()"))
end

@testcase :FUNC_COLON begin
	eval(parse("function ooewjdbwoedjfvobwe13(): Base.identity(9);end; ooewjdbwoedjfvobwe13()"))
end

@testcase :DICTS_ITER begin
	oqwjej = {"sdfre"=>3, "jkndivnf"=>"grewerg", "sdfwerf"=>{}}
	for i in oqwjej
		"Hello, " * oqwjej
	end
end

@testcase :IMPORT_LVL begin
	function nuinkjnbk()
		import Base
	end

	nuinkjnbk()
end

@testcase :IS_DEFINED begin
	function nuinkjnbk()
		info("Hello World!")
	end

	isdefined(nuinkjnbk)
end

@testcase :IS_DEFINED begin
	function nuinkjnbk()
		info("Hello World!")
	end

	isdefined("nuinkjnbk")
end

@testcase :IS_DEFINED begin
	isdefined(nuinkjnbk)
end

@testcase :IF_BOOLEAN begin
	myarray = {}
	if length(myarray)
		print(myarray)
	end
end

@testcase :IF_BOOLEAN begin
	mystring = "Hello World!"
	if mystring[1:0]
		print(mystring)
	end
end

@testcase :IF_BOOLEAN begin
	myarray = Number[]
	if myarray
		print(myarray)
	end
end

@example :IF_BOOLEAN begin
	argaegfawerg = 4313
	nvijenvierng = 12342
	if argaegfawerg & nvijenvierng
	elseif argaegfawerg | nvijenvierng
	end
end

@example :IF_BOOLEAN begin
	myarrray = rand(-10:10, (4, 5, 6))
	if length(myarrray)
	end
end

@testcase :PUSHAPPEND begin
	asdfgsdg = Float32[]
	push!(asdfgsdg, [3.14])
end

@testcase :PUSHAPPEND begin
	asdfgsdg = {}
	append!(asdfgsdg, 1)
end

@testcase :INDEX_BASE begin
	myarray = Float32[1:1000]
	for i in 0:1000
		println(myarray[i])
	end
end

@testcase :STRCONVERT begin
	mystring = "Hello, I'm test #"
	mynumber = 25
	println(mystring * mynumber)
end

@testcase :RANGEARRAY begin
	function jansdifnqwiljebf(asdfqwef::Array{Int})
		return asdfqwef
	end

	print(jansdifnqwiljebf(1:1000))
end

@testcase :KEYWORDARG begin
	function asdgergawerg(dfqefqwfqwr, qfregqwregr=Dict())
		print(dfqefqwfqwr)
		dump(qfregqwregr)
	end

	asdgergawerg("Hello!", qfregqwregr=Dict())
end

@testcase :TYPEASSERT begin
	aaefwregqre = 6.0
	egqrevqwefq = aaefwregqre::Int16
end

@testcase :ISOPERATOR begin
	if string(:Hello) <: Symbol
		print("World!")
	end
end

@testcase :SYMBOLEXPR begin
	ijansdvjiowneoijfnew(asdfwqe::Expr) = eval(asdfwqe)
	ijansdvjiowneoijfnew(:x)
end

@testcase :SYMBOLEXPR begin
	ijansdvjiowneoijfnew(asdfwqe::Symbol) = Expr(asdfwqe)
	ijansdvjiowneoijfnew(:(-x))
end

@testcase :DOT_EQUALS begin
	eval(parse("myarray=Int[];mymatches.===234"))
end

@testcase :TYPE_PARAM begin
	eval(parse("function asdasdfqwef{T::Any}(x::T, y::T)x+y;end"))
end

@testcase :INDIRECTFN begin
	function sdfgerfqefwerg(i::Int, f::Number)
		print() = print(i, f)
		print()
	end
	sdfgerfqefwerg(5, 12.3)
end

@testcase :COLLECTSUM begin
	mydict = Dict{Int32,Int32}()
	v = values(mydict)
	sum(v)
end

@testcase :DIMS_MATCH begin
	myarrray = rand(-10:10, (4, 5, 6))
	myarrray[:,:,1] * (myarrray[:,:,1] * 0.5)
end

@example :BOOLASSIGN begin
	eval(parse("while myexamplefoobarvariable='this is a string';end"))
end

@example :MAKESTRING begin
	mystrrrrring = String("Hello, World!")
end

@example :MAKESYMBOL begin
	mysyyyyyymbol = Symbol(:FOOBAR)
end


text:
# NOTE:
# in principal we only allow one to input a tuple of address to instruct
# but most of the instruct function in this file defines a forward method
# for single qubit version with an int as input.
#
#
# In order to make multi threading work, the state vector MUST be named as state

using YaoBase, BitBasis, LuxurySparse, StaticArrays
export instruct!

function YaoBase.instruct!(
    r::ArrayReg,
    op,
    locs::Tuple,
    control_locs::Tuple,
    control_bits::Tuple,
)
    instruct!(r.state, op, locs, control_locs, control_bits)
end

function YaoBase.instruct!(
    r::ArrayReg{1},
    op,
    locs::Tuple,
    control_locs::Tuple,
    control_bits::Tuple,
)
    instruct!(vec(r.state), op, locs, control_locs, control_bits)
end

function YaoBase.instruct!(r::ArrayReg, op, locs::Tuple)
    instruct!(r.state, op, locs)
end

function YaoBase.instruct!(r::ArrayReg{1}, op, locs::Tuple)
    instruct!(vec(r.state), op, locs)
end

function YaoBase.instruct!(
    r::ArrayReg,
    op,
    locs::Tuple,
    control_locs::Tuple,
    control_bits::Tuple,
    theta::Number,
)
    instruct!(r.state, op, locs, control_locs, control_bits, theta)
end

function YaoBase.instruct!(r::ArrayReg, op, locs::Tuple, theta::Number)
    instruct!(r.state, op, locs, theta)
end

"""
    SPECIALIZATION_LIST::Vector{Symbol}

A list of symbol for specialized gates/operators.
"""
const SPECIALIZATION_LIST =
    Symbol[:X, :Y, :Z, :S, :T, :Sdag, :Tdag, :H, :SWAP, :PSWAP, :CPHASE]

const THREAD_THRESHOLD = 10

# generates the threading expression
macro threads(ex)
    esc(quote
        if (Threads.nthreads() == 1) || log2dim1(state) < $THREAD_THRESHOLD
            @inbounds $ex
        else
            @inbounds Threads.@threads $ex
        end
    end)
end

# to avoid potential ambiguity, we limit them to tuple for now
# but they only has to be an iterator over integers
const Locations{T} = NTuple{N,T} where {N}
const BitConfigs{T} = NTuple{N,T} where {N}

function YaoBase.instruct!(
    state::AbstractVecOrMat{T1},
    operator::AbstractMatrix{T2},
    locs::Tuple{},
    control_locs::NTuple{C,Int} = (),
    control_bits::NTuple{C,Int} = (),
) where {T1,T2,M,C}
    return state
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T1},
    operator::AbstractMatrix{T2},
    locs::NTuple{M,Int},
    control_locs::NTuple{C,Int} = (),
    control_bits::NTuple{C,Int} = (),
) where {T1,T2,M,C}

    @warn "Element Type Mismatch: register $(T1), operator $(T2). Converting operator to match, this may cause performance issue"
    return instruct!(
        state,
        copyto!(similar(operator, T1), operator),
        locs,
        control_locs,
        control_bits,
    )
end

function _prepare_instruct(
    state,
    U,
    locs::NTuple{M},
    control_locs,
    control_bits::NTuple{C},
) where {M,C}
    N, MM = log2dim1(state), size(U, 1)

    locked_bits = MVector(control_locs..., locs...)
    locked_vals = MVector(control_bits..., (0 for k = 1:M)...)
    locs_raw_it = (b + 1 for b in itercontrol(N, setdiff(1:N, locs), zeros(Int, N - M)))
    locs_raw = SVector(locs_raw_it...)
    ic = itercontrol(N, locked_bits, locked_vals)
    return locs_raw, ic
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    operator::AbstractMatrix{T},
    locs::Tuple{},
    control_locs::NTuple{C,Int} = (),
    control_bits::NTuple{C,Int} = (),
) where {T,M,C}
    return state
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    operator::AbstractMatrix{T},
    locs::NTuple{M,Int},
    control_locs::NTuple{C,Int} = (),
    control_bits::NTuple{C,Int} = (),
) where {T,M,C}
    operator = sort_unitary(operator, locs)
    locs_raw, ic = _prepare_instruct(state, operator, locs, control_locs, control_bits)
    return _instruct!(state, autostatic(operator), locs_raw, ic)
end

function _instruct!(
    state::AbstractVecOrMat{T},
    U::AbstractMatrix{T},
    locs_raw::SVector,
    ic::IterControl,
) where {T}
    @threads for j = 1:length(ic)
        @inbounds i = ic[j]
        @inbounds unrows!(state, locs_raw .+ i, U)
    end
    return state
end

function _instruct!(
    state::AbstractVecOrMat{T},
    U::SDSparseMatrixCSC{T},
    locs_raw::SVector,
    ic::IterControl,
) where {T}
    work =
        ndims(state) == 1 ? similar(state, length(locs_raw)) :
        similar(state, length(locs_raw), size(state, 2))
    for j = 1:length(ic)
        @inbounds i = ic[j]
        @inbounds unrows!(state, locs_raw .+ i, U, work)
    end
    return state
end

YaoBase.instruct!(state::AbstractVecOrMat, U::IMatrix, locs::NTuple{N,Int}) where {N} =
    state
YaoBase.instruct!(state::AbstractVecOrMat, U::IMatrix, locs::Tuple{Int}) = state

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    U1::AbstractMatrix{T},
    (loc,)::Tuple{Int},
) where {T}
    a, c, b, d = U1
    instruct_kernel(state, loc, 1 << (loc - 1), 1 << loc, a, b, c, d)
    return state
end

@inline function instruct_kernel(state::AbstractVecOrMat, loc, step1, step2, a, b, c, d)
    @threads for j = 0:step2:size(state, 1)-step1
        @inbounds for i = j+1:j+step1
            u1rows!(state, i, i + step1, a, b, c, d)
        end
    end
    return state
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    U1::SDPermMatrix{T},
    (loc,)::Tuple{Int},
) where {T}
    U1.perm[1] == 1 && return instruct!(state, Diagonal(U1), (loc,))
    mask = bmask(loc)
    b, c = U1.vals
    step = 1 << (loc - 1)
    step_2 = 1 << loc

    @threads for j = 0:step_2:size(state, 1)-step
        @inbounds for i = j+1:j+step
            swaprows!(state, i, i + step, c, b)
        end
    end
    return state
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    U1::SDDiagonal{T},
    (loc,)::Tuple{Int},
) where {T}
    mask = bmask(loc)
    a, d = U1.diag
    step = 1 << (loc - 1)
    step_2 = 1 << loc
    @threads for j = 0:step_2:size(state, 1)-step
        @inbounds for i = j+1:j+step
            mulrow!(state, i, a)
            mulrow!(state, i + step, d)
        end
    end
    return state
end

# specialization
# paulis
function YaoBase.instruct!(
    state::AbstractVecOrMat,
    ::Val{:X},
    locs::NTuple{N,Int},
) where {N}
    mask = bmask(locs)
    do_mask = bmask(first(locs))
    @threads for b in basis(state)
        @inbounds if anyone(b, do_mask)
            i = b + 1
            i_ = flip(b, mask) + 1
            swaprows!(state, i, i_)
        end
    end
    return state
end

function YaoBase.instruct!(
    state::AbstractVecOrMat,
    ::Val{:H},
    locs::NTuple{N,Int},
) where {N}
    instruct!(state, YaoBase.Const.H, locs)
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    ::Val{:Y},
    locs::NTuple{N,Int},
) where {T,N}
    mask = bmask(Int, locs)
    do_mask = bmask(Int, first(locs))
    bit_parity = iseven(length(locs)) ? 1 : -1
    factor = T(-im)^length(locs)

    @threads for b in basis(Int, state)
        if anyone(b, do_mask)
            i = b + 1
            i_ = flip(b, mask) + 1
            factor1 = isodd(count_ones(b & mask)) ? -factor : factor
            factor2 = factor1 * bit_parity
            swaprows!(state, i, i_, factor2, factor1)
        end
    end
    return state
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    ::Val{:Z},
    locs::NTuple{N,Int},
) where {T,N}
    mask = bmask(Int, locs)
    @threads for b in basis(Int, state)
        if isodd(count_ones(b & mask))
            mulrow!(state, b + 1, -1)
        end
    end
    return state
end

for (G, FACTOR) in zip(
    [:S, :T, :Sdag, :Tdag],
    [:(im), :($(exp(im * π / 4))), :(-im), :($(exp(-im * π / 4)))],
)
    @eval function YaoBase.instruct!(
        state::AbstractVecOrMat{T},
        ::Val{$(QuoteNode(G))},
        locs::NTuple{N,Int},
    ) where {T,N}
        mask = bmask(Int, locs)
        @threads for b in basis(Int, state)
            mulrow!(state, b + 1, $FACTOR^count_ones(b & mask))
        end
        return state
    end
end

for (G, FACTOR) in zip(
    [:Z, :S, :T, :Sdag, :Tdag],
    [:(-1), :(im), :($(exp(im * π / 4))), :(-im), :($(exp(-im * π / 4)))],
)
    # no effect (to fix ambiguity)
    @eval YaoBase.instruct!(st::AbstractVecOrMat, ::Val{$(QuoteNode(G))}, ::Tuple{}) = st

    @eval function YaoBase.instruct!(
        state::AbstractVecOrMat{T},
        ::Val{$(QuoteNode(G))},
        (loc,)::Tuple{Int},
    ) where {T}
        mask = bmask(loc)
        step = 1 << (loc - 1)
        step_2 = 1 << loc
        @threads for j = 0:step_2:size(state, 1)-step
            for i = j+step+1:j+step_2
                mulrow!(state, i, $FACTOR)
            end
        end
        state
    end
end

# Specialized
import YaoBase: rot_mat

rot_mat(::Type{T}, ::Val{:Rx}, theta::Number) where {T} =
    T[cos(theta / 2) -im*sin(theta / 2); -im*sin(theta / 2) cos(theta / 2)]
rot_mat(::Type{T}, ::Val{:Ry}, theta::Number) where {T} =
    T[cos(theta / 2) -sin(theta / 2); sin(theta / 2) cos(theta / 2)]
rot_mat(::Type{T}, ::Val{:Rz}, theta::Number) where {T} =
    Diagonal(T[exp(-im * theta / 2), exp(im * theta / 2)])
rot_mat(::Type{T}, ::Val{:CPHASE}, theta::Number) where {T} =
    Diagonal(T[1, 1, 1, exp(im * theta)])
rot_mat(::Type{T}, ::Val{:PSWAP}, theta::Number) where {T} = rot_mat(T, Const.SWAP, theta)

for G in [:Rx, :Ry, :Rz, :CPHASE]
    # forward single gates
    @eval function YaoBase.instruct!(
        state::AbstractVecOrMat{T},
        g::Val{$(QuoteNode(G))},
        locs::NTuple{N3,Int},
        control_locs::NTuple{N1,Int},
        control_bits::NTuple{N2,Int},
        theta::Number,
    ) where {T,N1,N2,N3}
        m = rot_mat(T, g, theta)
        instruct!(state, m, locs, control_locs, control_bits)
        return state
    end
end # for

@inline function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    ::Val{:Rx},
    (loc,)::Tuple{Int},
    theta::Number,
) where {T,N}
    b, a = sincos(theta / 2)
    instruct_kernel(state, loc, 1 << (loc - 1), 1 << loc, a, -im * b, -im * b, a)
    return state
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    ::Val{:Ry},
    (loc,)::Tuple{Int},
    theta::Number,
) where {T,N}
    b, a = sincos(theta / 2)
    instruct_kernel(state, loc, 1 << (loc - 1), 1 << loc, a, -b, b, a)
    return state
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    ::Val{:Rz},
    (loc,)::Tuple{Int},
    theta::Number,
) where {T,N}
    a = exp(-im * theta / 2)
    instruct_kernel(state, loc, 1 << (loc - 1), 1 << loc, a, zero(T), zero(T), a')
    return state
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    ::Val{:CPHASE},
    locs::NTuple{N,Int},
    theta::Number,
) where {T,N}
    m = rot_mat(T, Val(:CPHASE), theta)
    instruct!(state, m, locs)
    return state
end

# forward single gates
function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    g::Val,
    locs::Union{Int,NTuple{N1,Int}},
    theta::Number,
) where {T,N1}
    instruct!(state, g, locs, (), (), theta)
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    ::Val{:X},
    locs::NTuple{N1,Int},
    control_locs::NTuple{N2,Int},
    control_bits::NTuple{N3,Int},
) where {T,N1,N2,N3}

    ctrl = controller((control_locs..., locs[1]), (control_bits..., 0))
    mask2 = bmask(locs)
    @threads for b in basis(state)
        if ctrl(b)
            i = b + 1
            i_ = flip(b, mask2) + 1
            swaprows!(state, i, i_)
        end
    end
    return state
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    ::Val{:Y},
    locs::NTuple{N1,Int},
    control_locs::NTuple{N2,Int},
    control_bits::NTuple{N3,Int},
) where {T,N1,N2,N3}

    ctrl = controller((control_locs..., locs[1]), (control_bits..., 0))
    mask2 = bmask(locs)
    @threads for b in basis(state)
        if ctrl(b)
            i = b + 1
            i_ = flip(b, mask2) + 1
            swaprows!(state, i, i_, im, -im)
        end
    end
    return state
end

for (G, FACTOR) in zip(
    [:Z, :S, :T, :Sdag, :Tdag],
    [:(-1), :(im), :($(exp(im * π / 4))), :(-im), :($(exp(-im * π / 4)))],
)
    @eval function YaoBase.instruct!(
        state::AbstractVecOrMat{T},
        ::Val{$(QuoteNode(G))},
        locs::NTuple{N1,Int},
        control_locs::NTuple{N2,Int},
        control_bits::NTuple{N3,Int},
    ) where {T,N1,N2,N3}

        ctrl = controller([control_locs..., locs[1]], [control_bits..., 1])
        @threads for b in basis(state)
            if ctrl(b)
                mulrow!(state, b + 1, $FACTOR)
            end
        end
        return state
    end
end

## single controlled paulis
for G in [:X, :Y, :Z, :S, :T, :Sdag, :Tdag]
    # forward single controlled
    @eval YaoBase.instruct!(
        state::AbstractVecOrMat,
        g::Val{$(QuoteNode(G))},
        locs::NTuple{N,Int},
        control_locs::Tuple{Int},
        control_bits::Tuple{Int},
    ) where {N} = instruct!(state, g, locs, control_locs..., control_bits...)

    # forward single location
    @eval YaoBase.instruct!(
        state::AbstractVecOrMat,
        g::Val{$(QuoteNode(G))},
        locs::Tuple{Int},
        control_locs::Int,
        control_bits::Int,
    ) where {N} = instruct!(state, g, locs..., control_locs, control_bits)

end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    ::Val{:X},
    loc::Int,
    control_locs::Int,
    control_bits::Int,
) where {T}

    mask2 = bmask(loc)
    mask = bmask(control_locs, loc)
    step = 1 << (control_locs - 1)
    step_2 = 1 << control_locs
    start = control_bits == 1 ? step : 0
    @threads for j = start:step_2:size(state, 1)-step+start
        for b = j:j+step-1
            @inbounds if allone(b, mask2)
                i = b + 1
                i_ = flip(b, mask2) + 1
                swaprows!(state, i, i_)
            end
        end
    end
    return state
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    ::Val{:Y},
    loc::Int,
    control_locs::Int,
    control_bits::Int,
) where {T}

    mask2 = bmask(loc)
    mask = bmask(control_locs, loc)
    step = 1 << (control_locs - 1)
    step_2 = 1 << control_locs
    start = control_bits == 1 ? step : 0
    @threads for j = start:step_2:size(state, 1)-step+start
        for b = j:j+step-1
            @inbounds if allone(b, mask2)
                i = b + 1
                i_ = flip(b, mask2) + 1
                if allone(b, mask2)
                    factor = T(im)
                else
                    factor = T(-im)
                end
                swaprows!(state, i, i_, -factor, factor)
            end
        end
    end
    return state
end

for (G, FACTOR) in zip(
    [:Z, :S, :T, :Sdag, :Tdag],
    [:(-1), :(im), :($(exp(im * π / 4))), :(-im), :($(exp(-im * π / 4)))],
)
    @eval function YaoBase.instruct!(
        state::AbstractVecOrMat{T},
        ::Val{$(QuoteNode(G))},
        loc::Int,
        control_locs::Int,
        control_bits::Int,
    ) where {T}

        mask2 = bmask(loc)
        step = 1 << (control_locs - 1)
        step_2 = 1 << control_locs
        start = control_bits == 1 ? step : 0
        @threads for j = start:step_2:size(state, 1)-step+start
            for i = j+1:j+step
                if allone(i - 1, mask2)
                    mulrow!(state, i, $FACTOR)
                end
            end
        end
        return state
    end
end

# 2-qubit gates
function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    ::Val{:SWAP},
    locs::Tuple{Int,Int},
) where {T}

    mask1 = bmask(locs[1])
    mask2 = bmask(locs[2])
    mask12 = mask1 | mask2
    @threads for b in basis(state)
        if b & mask1 == 0 && b & mask2 == mask2
            i = b + 1
            i_ = b ⊻ mask12 + 1
            swaprows!(state, i, i_)
        end
    end
    return state
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    ::Val{:PSWAP},
    locs::Tuple{Int,Int},
    theta::Number,
) where {T}
    mask1 = bmask(locs[1])
    mask2 = bmask(locs[2])
    mask12 = mask1 | mask2
    a = T(cos(theta / 2))
    c = T(-im * sin(theta / 2))
    e = T(exp(-im / 2 * theta))
    @threads for b in basis(state)
        if b & mask1 == 0
            i = b + 1
            i_ = b ⊻ mask12 + 1
            if b & mask2 == mask2
                u1rows!(state, i, i_, a, c, c, a)
            else
                mulrow!(state, i, e)
                mulrow!(state, i_, e)
            end
        end
    end
    return state
end

function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    ::Val{:PSWAP},
    locs::Tuple{Int,Int},
    control_locs::NTuple{C,Int},
    control_bits::NTuple{C,Int},
    theta::Number,
) where {T,C}
    mask1 = bmask(locs[1])
    mask2 = bmask(locs[2])
    mask12 = mask1 | mask2
    a = T(cos(theta / 2))
    c = T(-im * sin(theta / 2))
    e = T(exp(-im / 2 * theta))
    @threads for b in itercontrol(
        log2i(size(state, 1)),
        collect(Int, control_locs),
        collect(Int, control_bits),
    )
        if b & mask1 == 0
            i = b + 1
            i_ = b ⊻ mask12 + 1
            if b & mask2 == mask2
                u1rows!(state, i, i_, a, c, c, a)
            else
                mulrow!(state, i, e)
                mulrow!(state, i_, e)
            end
        end
    end
    return state
end