text:
# This method is called when a SemidiscretizationHyperbolic is constructed.
# It constructs the basic `cache` used throughout the simulation to compute
# the RHS etc.
function create_cache(mesh::UnstructuredQuadMesh, equations,
                      dg::DG, RealT, uEltype)

  polydeg_ = polydeg(dg.basis)
  nvars = nvariables(equations)

  if polydeg_ > mesh.polydeg && any(mesh.element_is_curved)
    throw(ArgumentError("polynomial degree of DG (= $polydeg_) must be less than or equal to mesh polynomial degree (= $(mesh.polydeg))"))
  end

  elements = init_elements(RealT, uEltype, mesh, dg.basis.nodes, nvars, polydeg_)

  interfaces = init_interfaces(uEltype, mesh, nvars, polydeg_)

  # TODO: if the mesh is periodic we probably can avoid creating this "empty" boundary container
  if isperiodic(mesh)
    boundaries = UnstructuredBoundaryContainer2D{RealT, uEltype, nvars, polydeg_}(0)
  else
    boundaries = init_boundaries(RealT, uEltype, mesh, elements, nvars, polydeg_)
  end

  cache = (; elements, interfaces, boundaries)

  return cache
end

# Note! The mesh also passed to some functions, e.g., calc_volume_integral! to dispatch on the
#       correct version and use the ::UnstructuredQuadMesh variable type below to keep track on it
function rhs!(du, u, t,
              mesh::UnstructuredQuadMesh, equations,
              initial_condition, boundary_conditions, source_terms,
              dg::DG, cache)
  # Reset du
  @timeit_debug timer() "reset ∂u/∂t" du .= zero(eltype(du))

  # Calculate volume integral
  @timeit_debug timer() "volume integral" calc_volume_integral!(
    du, u, mesh,
    have_nonconservative_terms(equations), equations,
    dg.volume_integral, dg, cache)

  # Prolong solution to interfaces
  @timeit_debug timer() "prolong2interfaces" prolong2interfaces!(
    cache, u, mesh, equations, dg)

  # Calculate interface fluxes
  @timeit_debug timer() "interface flux" calc_interface_flux!(
    cache.elements.surface_flux_values, mesh,
    have_nonconservative_terms(equations), equations,
    dg, cache)

  # Prolong solution to boundaries
  @timeit_debug timer() "prolong2boundaries" prolong2boundaries!(
    cache, u, mesh, equations, dg)

  # Calculate boundary fluxes
  # TODO: Meshes, remove initial condition as an input argument here, only needed for hacky BCs
  @timeit_debug timer() "boundary flux" calc_boundary_flux!(
    cache, t, boundary_conditions, mesh, equations, dg, initial_condition)

  # Calculate surface integrals
  @timeit_debug timer() "surface integral" calc_surface_integral!(
    du, mesh, equations, dg, cache)

  # Apply Jacobian from mapping to reference element
  #  Note! this routine is reused from dg_curved/dg_2d.jl
  @timeit_debug timer() "Jacobian" apply_jacobian!(
    du, mesh, equations, dg, cache)

  # Calculate source terms
  @timeit_debug timer() "source terms" calc_sources!(
    du, u, t, source_terms, equations, dg, cache)

  return nothing
end

# compute volume contribution of the DG approximation with the divergence of the contravariant fluxes
function calc_volume_integral!(du, u,
                               mesh::UnstructuredQuadMesh,
                               nonconservative_terms::Val{false}, equations,
                               volume_integral::VolumeIntegralWeakForm,
                               dg::DGSEM, cache)
  @unpack derivative_dhat = dg.basis
  # TODO: Meshes. This is basically the same as the version for the `CurvedMesh`;
  #       the only difference is that `X_xi, X_eta, Y_xi, Y_eta` are stored in
  #       separate arrays here and in one block array `contravariant_vectors`
  #       for the `CurvedMesh`.
  @unpack X_xi, X_eta, Y_xi, Y_eta = cache.elements

  @threaded for element in eachelement(dg, cache)
    for j in eachnode(dg), i in eachnode(dg)
      u_node = get_node_vars(u, equations, dg, i, j, element)

      # compute the physical fluxes in each Cartesian direction
      x_flux = flux(u_node, 1, equations)
      y_flux = flux(u_node, 2, equations)

      # compute the contravariant flux in the x-direction
      flux1  = Y_eta[i, j, element] * x_flux - X_eta[i, j, element] * y_flux
      for ii in eachnode(dg)
        integral_contribution = derivative_dhat[ii, i] * flux1
        add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element)
      end

      # compute the contravariant flux in the y-direction
      flux2  = -Y_xi[i, j, element] * x_flux + X_xi[i, j, element] * y_flux
      for jj in eachnode(dg)
        integral_contribution = derivative_dhat[jj, j] * flux2
        add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element)
      end
    end
  end

  return nothing
end

# prolong the solution into the convenience array in the interior interface container
# Note! this routine is for quadrilateral elements with "right-handed" orientation
function prolong2interfaces!(cache, u,
                             mesh::UnstructuredQuadMesh,
                             equations, dg::DG)
  @unpack interfaces = cache

  @threaded for interface in eachinterface(dg, cache)
    primary_element   = interfaces.element_ids[1, interface]
    secondary_element = interfaces.element_ids[2, interface]

    primary_side   = interfaces.element_side_ids[1, interface]
    secondary_side = interfaces.element_side_ids[2, interface]

    if primary_side == 1
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[1, v, i, interface] = u[v, i, 1, primary_element]
      end
    elseif primary_side == 2
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[1, v, i, interface] = u[v, nnodes(dg), i, primary_element]
      end
    elseif primary_side == 3
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[1, v, i, interface] = u[v, i, nnodes(dg), primary_element]
      end
    else # primary_side == 4
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[1, v, i, interface] = u[v, 1, i, primary_element]
      end
    end

    if secondary_side == 1
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[2, v, i, interface] = u[v, i, 1, secondary_element]
      end
    elseif secondary_side == 2
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[2, v, i, interface] = u[v, nnodes(dg), i, secondary_element]
      end
    elseif secondary_side == 3
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[2, v, i, interface] = u[v, i, nnodes(dg), secondary_element]
      end
    else # secondary_side == 4
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[2, v, i, interface] = u[v, 1, i, secondary_element]
      end
    end
  end

  return nothing
end

# compute the numerical flux interface coupling between two elements on an unstructured quadrilateral mesh
function calc_interface_flux!(surface_flux_values,
                              mesh::UnstructuredQuadMesh,
                              nonconservative_terms::Val{false}, equations, dg::DG, cache)
  @unpack surface_flux = dg
  @unpack u, start_index, index_increment, element_ids, element_side_ids = cache.interfaces
  @unpack normals, scaling = cache.elements

  @threaded for interface in eachinterface(dg, cache)
    # Get neighboring elements
    primary_element   = element_ids[1, interface]
    secondary_element = element_ids[2, interface]

    # Get the local side id on which to compute the flux
    primary_side   = element_side_ids[1, interface]
    secondary_side = element_side_ids[2, interface]

    # initial index for the coordinate system on the secondary element
    secondary_index = start_index[interface]

    # loop through the primary element coordinate system and compute the interface coupling
    for primary_index in eachnode(dg)
      # pull the primary and secondary states from the boundary u values
      u_ll = get_one_sided_surface_node_vars(u, equations, dg, 1, primary_index, interface)
      u_rr = get_one_sided_surface_node_vars(u, equations, dg, 2, secondary_index, interface)

      # TODO: Meshes, performance. The rotation shouldn't be hard-coded here
      #       since it's not generally applicable to all equations. Thus, it
      #       should be moved to the numerical flux, cf. `CurvedMesh`.
      #       This will make this part strictly more general and allow additional
      #       performance optimizations.

      # pull the directional vectors and scaling factors
      #   Note! this assumes a conforming approximation, more must be done in terms of the normals
      #         for hanging nodes and other non-conforming approximation spaces
      normal_vector = get_surface_normal(normals, primary_index, primary_side, primary_element)
      scaling_ll    = scaling[primary_index, primary_side, primary_element]
      scaling_rr    = scaling[secondary_index, secondary_side, secondary_element]

      # rotate states
      u_tilde_ll = rotate_to_x(u_ll, normal_vector, equations)
      u_tilde_rr = rotate_to_x(u_rr, normal_vector, equations)

      # Call pointwise Riemann solver in the rotated direction
      flux_tilde = surface_flux(u_tilde_ll, u_tilde_rr, 1, equations)

      # backrotate the flux into the original direction
      flux = rotate_from_x(flux_tilde, normal_vector, equations)

      # Scale the flux appropriately and copy back to primary/secondary element storage
      # Note the sign change for the normal flux in the secondary element!
      for v in eachvariable(equations)
        surface_flux_values[v, primary_index  , primary_side  , primary_element  ] =  flux[v] * scaling_ll
        surface_flux_values[v, secondary_index, secondary_side, secondary_element] = -flux[v] * scaling_rr
      end

      # increment the index of the coordinate system in the secondary element
      secondary_index += index_increment[interface]
    end
  end

  return nothing
end

# move the approximate solution onto physical boundaries within a "right-handed" element
function prolong2boundaries!(cache, u,
                             mesh::UnstructuredQuadMesh,
                             equations, dg::DG)
  @unpack boundaries = cache

  # TODO: Meshes. Should this become `eachboundary(dg, cache)`?
  @threaded for boundary in eachboundary(boundaries)
    element = boundaries.element_id[boundary]
    side    = boundaries.element_side_id[boundary]

    if side == 1
      for l in eachnode(dg), v in eachvariable(equations)
        boundaries.u[1, v, l, boundary] = u[v, l, 1, element]
      end
    elseif side == 2
      for l in eachnode(dg), v in eachvariable(equations)
        boundaries.u[1, v, l, boundary] = u[v, nnodes(dg), l, element]
      end
    elseif side == 3
      for l in eachnode(dg), v in eachvariable(equations)
        boundaries.u[1, v, l, boundary] = u[v, l, nnodes(dg), element]
      end
    else # side == 4
      for l in eachnode(dg), v in eachvariable(equations)
        boundaries.u[1, v, l, boundary] = u[v, 1, l, element]
      end
    end
  end

  return nothing
end

# TODO: Taal dimension agnostic
function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic,
                             mesh::UnstructuredQuadMesh, equations, dg::DG,
                             initial_condition)
  @assert isempty(eachboundary(cache.boundaries))
end

# TODO: Meshes; `intial_condition` argument for convenience,
#       cleanup later with better boundar condition handling
function calc_boundary_flux!(cache, t, boundary_condition,
                             mesh::UnstructuredQuadMesh, equations, dg::DG,
                             initial_condition)
  @unpack surface_flux = dg
  @unpack normals, scaling, surface_flux_values = cache.elements
  @unpack u, element_id, element_side_id, node_coordinates, name  = cache.boundaries

  @threaded for boundary in eachboundary(cache.boundaries)
    # Get the element and side IDs on the primary element
    primary_element = element_id[boundary]
    primary_side    = element_side_id[boundary]

    for primary_index in eachnode(dg)
      # hacky way to set "exact solution" boundary conditions. Only used to test the orientation
      # for a mesh with flipped elements
      u_external = initial_condition((node_coordinates[1, primary_index, boundary],
                                      node_coordinates[2, primary_index, boundary]),
                                      t, equations)

      # pull the left state from the boundary u values on the primary element as well as the
      # directional vectors and scaling
      #   Note! this assumes a conforming approximation, more must be done in terms of the normals
      #         for hanging nodes and other non-conforming approximation spaces
      u_ll = get_one_sided_surface_node_vars(u, equations, dg, 1, primary_index, boundary)

      normal_vector  = get_surface_normal(normals, primary_index, primary_side, primary_element)
      scaling_ll  = scaling[primary_index, primary_side, primary_element]

      # rotate states
      u_tilde_ll       = rotate_to_x(u_ll, normal_vector, equations)
      u_tilde_external = rotate_to_x(u_external, normal_vector, equations)

      # Call pointwise Riemann solver in the rotated direction
      flux_tilde = surface_flux(u_tilde_ll, u_tilde_external, 1, equations)

      # backrotate the flux into the original direction
      flux = rotate_from_x(flux_tilde, normal_vector, equations)

      # Scale the flux appropriately and copy back to primary element storage
      for v in eachvariable(equations)
        surface_flux_values[v, primary_index, primary_side, primary_element] = flux[v] * scaling_ll
      end
    end
  end

  return nothing
end

# Note! The local side numbering for the unstructured quadrilateral element implementation differs
#       from the structured TreeMesh or CurvedMesh local side numbering:
#
#      TreeMesh/CurvedMesh sides   versus   UnstructuredMesh sides
#                  4                                  3
#          -----------------                  -----------------
#          |               |                  |               |
#          | ^ eta         |                  | ^ eta         |
#        1 | |             | 2              4 | |             | 2
#          | |             |                  | |             |
#          | ---> xi       |                  | ---> xi       |
#          -----------------                  -----------------
#                  3                                  1
# Therefore, we require a different surface integral routine here despite their similar structure.
function calc_surface_integral!(du, mesh::UnstructuredQuadMesh,
                                equations, dg::DGSEM, cache)
  @unpack boundary_interpolation = dg.basis
  @unpack surface_flux_values = cache.elements

  @threaded for element in eachelement(dg, cache)
    for l in eachnode(dg), v in eachvariable(equations)
      # surface contribution along local sides 2 and 4 (fixed x and y varies)
      du[v, 1,          l, element] += ( surface_flux_values[v, l, 4, element]
                                          * boundary_interpolation[1, 1] )
      du[v, nnodes(dg), l, element] += ( surface_flux_values[v, l, 2, element]
                                          * boundary_interpolation[nnodes(dg), 2] )
      # surface contribution along local sides 1 and 3 (fixed y and x varies)
      du[v, l, 1,          element] += ( surface_flux_values[v, l, 1, element]
                                          * boundary_interpolation[1, 1] )
      du[v, l, nnodes(dg), element] += ( surface_flux_values[v, l, 3, element]
                                          * boundary_interpolation[nnodes(dg), 2] )
    end
  end

  return nothing
end


text:
function test_base_function(f)
    F = (f)()

    @testset "Testing properties" begin
        test_properties = get(base_functions_properties, f, (false,false))
        @test MLK.isstationary(F) == test_properties[1]
        @test MLK.isisotropic(F)  == test_properties[2]
    end

    # Test that initial value is set correctly compared to reference
    @testset "Testing $(MLK.base_initiate)" begin
        for T in FloatingPointTypes
            s = MLK.base_initiate(F, T)
            s_tmp = convert(T, get(base_functions_initiate, f, 0))

            @test typeof(s) == T
            @test s == s_tmp
        end
    end

    # Test that aggregation works correctly
    @testset "Testing $(MLK.base_aggregate)" begin
        for T in FloatingPointTypes
            f_tmp = get(base_functions_aggregate, f, (s,x,y) -> convert(T, NaN))
            for x in (T(0), T(1), T(2)), y in (T(0), T(1), T(2))
                s = MLK.base_aggregate(F, T(1), x, y)
                s_tmp = f_tmp(T(1), x, y)

                @test typeof(s) == T
                @test s == s_tmp
            end
        end
    end

    # Test that return works correctly
    @testset "Testing $(MLK.base_return)" begin
        for T in FloatingPointTypes
            f_tmp = get(base_functions_return, f, s -> s)
            for x in (T(0), T(1), T(2))
                s = MLK.base_return(F, x)
                s_tmp = f_tmp(x)

                @test typeof(s) == T
                @test s == s_tmp
            end
        end
    end
end

for f in base_functions
    @testset "Testing $f" begin
        test_base_function(f)
    end
end

ext: jl


text:
# --- Calling Embedded Julia from MATLAB --- #

# entry point for MATLAB calling Julia
jl_mex(plhs::Vector{Ptr{Cvoid}}, prhs::Vector{Ptr{Cvoid}}) = jl_mex_inner(plhs, prhs)

# runs julia function with mex function inputs, catching errors if they occur
function jl_mex_inner(plhs::Vector{Ptr{Cvoid}}, prhs::Vector{Ptr{Cvoid}})

    # get number of outputs
    nlhs = length(plhs)

    for i = 1:nlhs
        # default output is boolean false
        none = MATLAB.mxarray(false);
        # transfer ownership to MATLAB
        none.own = false
        # give pointer to MATLAB
        plhs[i] = none.ptr
    end

    try

        # extract function and arguments (function in first slot, arguments in remaining slots)
        fun = Core.eval(Main, Meta.parse(MATLAB.jvalue(MATLAB.MxArray(prhs[1], false))))
        args = MATLAB.MxArray.(prhs[2:end], false)

        # call Julia function
        vals = fun(args)

        # convert outputs to MxArray type
        for i = 1:length(vals)
            # stop early if max number of outputs is reached
            if i > nlhs-1
                break
            end
            # create MATLAB array for output
            mx = MATLAB.mxarray(vals[i])
            # transfer ownership to MATLAB
            mx.own = false
            # give pointer to MATLAB
            plhs[1+i] = mx.ptr
        end

    catch exn

        # get backtrace
        bt = catch_backtrace()

        # create MATLABException from exception and backtrace
        mexn = MatlabException(exn, bt)

        # return MATLABexception in first output slot
        plhs[1] = mexn.ptr

    end

end

# --- Functions Used by MATLAB when calling embedded Julia --- #

# evaluates a Julia expressions
jl_eval(exprs::Vector{MATLAB.MxArray}) = [Core.eval(Main, Meta.parse(MATLAB.jvalue(e))) for e in exprs]

# Call a julia function, possibly with keyword arguments.
#
# The values in the args array are interpreted as follows:
#      Index = Meaning
# -----------------------------------------------------------
#          1 = the function to call
#          2 = an integer, npos, the number of positional arguments
#   3:2+npos = positional arguments
# 3+npos:end = keyword arguments, in keyword/value pairs
#
# If npos < 0, all arguments are assumed positional.
function jl_call_kw(args::Vector{MATLAB.MxArray})

    # convert arguments to Julia objects
    vals = MATLAB.jvalue.(args)
    nvals = length(args)

    # first argument is the function
    func = Symbol(vals[1])

    # second argument is the number of positional arguments
    npos = vals[2]

    # if npos is negative, all arguments are positional
    if npos < 0
        npos = nvals - 2
    end

    # get number of keyword arguments
    nkw = div(nvals - 2 - npos, 2)

    # initialize arguments to julia expression
    expr_args = Array{Any, 1}(undef, npos+nkw)

    # add positional arguments
    for i = 1:npos
        expr_args[i] = vals[2+i]
    end

    # add keyword arguments
    for i = 1:nkw
        # assemble the key-value pair
        kw = Symbol(vals[2+npos+(2*i-1)])
        val = vals[2+npos+(2*i)]
        expr_args[npos+i] = Expr(:kw, kw, val)
    end

    # construct the expression
    expr = Expr(:call, func, expr_args...)

    # return the evaluated expression
    return [Core.eval(Main, expr)]
end

# used for mimicing a basic Julia repl from the MATLAB console
input(prompt::String="julia> ") = call_matlab(1, "input", prompt, "s")[1]

# -------------- #

# The functions below this point are for interrupt handling and stdout/stderr
# redirection.  These features have various issues that need to be resolved
# before being incorporated with the rest of the code.

# the following function was the original entry point for calling Julia
# from MATLAB which enables interrupt handling and stdout/stderr redirection

# const jl_mex_call_depth = Ref{Int}(0)
# function jl_mex(plhs::Vector{Ptr{Cvoid}}, prhs::Vector{Ptr{Cvoid}})
#
#     # TODO: manually turn stdout and stderr redirection on/off?
#
#     # get stdout and stderr before changing them
#     stdout = Base.stdout
#     stderr = Base.stderr
#
#     # redirect stdout and stderr to matlab terminal
#     mexstdout_rd, mexstdout_wr = redirect_stdout()
#     mexerrout_rd, mexerrout_wr = redirect_stderr()
#
#     # start printing stdout and stderr to matlab console
#     t1 = @async readloop(mexstdout_rd, 1)
#     t2 = @async readloop(mexstderr_rd, 2)
#
#     jl_mex_call_depth[] += 1
#     try
#         if jl_mex_call_depth[] == 1
#             # call with interrupt handling
#             jl_mex_outer(plhs, prhs)
#         else
#             # call without interrupt handling
#             jl_mex_inner(plhs, prhs)
#         end
#     finally
#         jl_mex_call_depth[] -= 1
#     end
#
#     # stop printing stdout and stderr to matlab console
#     @async Base.throwto(t1, InterruptException())
#     @async Base.throwto(t2, InterruptException())
#
#     # restore stdout and stderr
#     redirect_stdout(stdout)
#     redirect_stderr(stderr)
#
#     # TODO: when called from the MATLAB engine, Julia will overwrite the last line
#     # of terminal output. fix this.
#
# end

# this function sets up interrupt handling before calling jl_mex_inner. The
# sync/async loop sometimes causes the output of the MEX function to become
# not be assigned (as least as far as MATLAB can tell).  This is likely due to
# the fact that MATLAB is not thread-safe.  Undocumented matlab references
# claim that the interrupt catching function is thread-safe, so the issue
# therefore likely lies with the use of @async on jl_mex_inner.  I ran into this
# issue while setting up callback functions, but have not taken the time to create
# a MWE which showcases the problem.

# function jl_mex_outer(plhs::Vector{Ptr{Cvoid}}, prhs::Vector{Ptr{Cvoid}})
#     @sync begin
#         # do the desired computation
#         main_task = @async jl_mex_inner(plhs, prhs)
#
#         # check for interrupt
#         @async while(!istaskdone(main_task))
#             if is_interrupt_pending()
#                 Base.throwto(main_task, InterruptException())
#             end
#             yield()
#         end
#     end
# end

# these function is used to poll MATLAB for interrupts, due to the issues
# described above, it is not used in the code

# is_interrupt_pending() = ccall(ut_is_interrupt_pending[], UInt8, ()) != 0
#
# function check_for_interrupt()
#     if is_interrupt_pending()
#         throw(InterruptException())
#     end
# end

# these functions are used for stdout/stderr redirection.  I found that if the
# code tries to print too much to the matlab console in a short period of time
# MATLAB will crash.  Therefore I removed them from the code.  Additionally,
# the redirection is unnecessary (and probably undesirable since it consumes
# some computational resources) if MATLAB is launched from the terminal

# function fwrite(fid, msg)
#     call_matlab(0, "fwrite", convert(Float64, fid), msg, "char")
#     return nothing
# end

# function readloop(stream, fid)
#     try
#         while isopen(stream)
#             fwrite(fid, readavailable(stream))
#         end
#     finally
#     end
# end


text:
import HTTP
import Markdown: htmlesc
import UUIDs: UUID
import Pkg

const found_is_pluto_dev = Ref{Union{Nothing,Bool}}(nothing)
function is_pluto_dev()
    if found_is_pluto_dev[] !== nothing
        return found_is_pluto_dev[]
    end
    found_is_pluto_dev[] = try
        deps = Pkg.dependencies()

        p_index = findfirst(p -> p.name == "Pluto", deps)
        p = deps[p_index]

        p.is_tracking_path
    catch
        false
    end
end

function frontend_directory(; allow_bundled::Bool=true)
    if allow_bundled && isdir(project_relative_path("frontend-dist")) && (get(ENV, "JULIA_PLUTO_FORCE_BUNDLED", "nein") == "ja" || !is_pluto_dev())
        "frontend-dist"
    else
        "frontend"
    end
end

function should_cache(path::String)
    dir, filename = splitdir(path)
    endswith(dir, "frontend-dist") && occursin(r"\.[0-9a-f]{8}\.", filename)
end

# Serve everything from `/frontend`, and create HTTP endpoints to open notebooks.

"Attempts to find the MIME pair corresponding to the extension of a filename. Defaults to `text/plain`."
function mime_fromfilename(filename)
    # This bad boy is from: https://developer.mozilla.org/en-US/docs/Web/HTTP/Basics_of_HTTP/MIME_types/Common_types
    mimepairs = Dict(".aac" => "audio/aac", ".bin" => "application/octet-stream", ".bmp" => "image/bmp", ".css" => "text/css", ".csv" => "text/csv", ".eot" => "application/vnd.ms-fontobject", ".gz" => "application/gzip", ".gif" => "image/gif", ".htm" => "text/html", ".html" => "text/html", ".ico" => "image/vnd.microsoft.icon", ".jpeg" => "image/jpeg", ".jpg" => "image/jpeg", ".js" => "text/javascript", ".json" => "application/json", ".jsonld" => "application/ld+json", ".mjs" => "text/javascript", ".mp3" => "audio/mpeg", ".mp4" => "video/mp4", ".mpeg" => "video/mpeg", ".oga" => "audio/ogg", ".ogv" => "video/ogg", ".ogx" => "application/ogg", ".opus" => "audio/opus", ".otf" => "font/otf", ".png" => "image/png", ".pdf" => "application/pdf", ".rtf" => "application/rtf", ".sh" => "application/x-sh", ".svg" => "image/svg+xml", ".tar" => "application/x-tar", ".tif" => "image/tiff", ".tiff" => "image/tiff", ".ttf" => "font/ttf", ".txt" => "text/plain", ".wav" => "audio/wav", ".weba" => "audio/webm", ".webm" => "video/webm", ".webp" => "image/webp", ".woff" => "font/woff", ".woff2" => "font/woff2", ".xhtml" => "application/xhtml+xml", ".xml" => "application/xml", ".xul" => "application/vnd.mozilla.xul+xml", ".zip" => "application/zip", ".wasm" => "application/wasm")
    file_extension = getkey(mimepairs, '.' * split(filename, '.')[end], ".txt")
    MIME(mimepairs[file_extension])
end

const day = let
    second = 1
    hour = 60second
    day = 24hour
end

function asset_response(path; cacheable::Bool=false)
    if !isfile(path) && !endswith(path, ".html")
        return asset_response(path * ".html"; cacheable)
    end
    if isfile(path)
        data = read(path)
        response = HTTP.Response(200, data)
        m = mime_fromfilename(path)
        push!(response.headers, "Content-Type" => Base.istextmime(m) ? "$(m); charset=UTF-8" : string(m))
        push!(response.headers, "Content-Length" => string(length(data)))
        push!(response.headers, "Access-Control-Allow-Origin" => "*")
        cacheable && push!(response.headers, "Cache-Control" => "public, max-age=$(30day), immutable")
        response
    else
        HTTP.Response(404, "Not found!")
    end
end

function error_response(status_code::Integer, title, advice, body="")
    template = read(project_relative_path(frontend_directory(), "error.jl.html"), String)

    body_title = body == "" ? "" : "Error message:"
    filled_in = replace(replace(replace(replace(replace(template,
        "\$STYLE" => """"""),
        "\$TITLE" => title),
        "\$ADVICE" => advice),
        "\$BODYTITLE" => body_title),
        "\$BODY" => htmlesc(body))

    response = HTTP.Response(status_code, filled_in)
    push!(response.headers, "Content-Type" => string(mime_fromfilename(".html")))
    response
end

function notebook_response(notebook; home_url="./", as_redirect=true)
    if as_redirect
        response = HTTP.Response(302, "")
        push!(response.headers, "Location" => home_url * "edit?id=" * string(notebook.notebook_id))
        return response
    else
        HTTP.Response(200, string(notebook.notebook_id))
    end
end

"""
Return whether the `request` was authenticated in one of two ways:
1. the session's `secret` was included in the URL as a search parameter, or
2. the session's `secret` was included in a cookie.
"""
function is_authenticated(session::ServerSession, request::HTTP.Request)
    (
        secret_in_url = try
            uri = HTTP.URI(request.target)
            query = HTTP.queryparams(uri)
            get(query, "secret", "") == session.secret
        catch e
            @warn "Failed to authenticate request using URL" exception = (e, catch_backtrace())
            false
        end
    ) || (
        secret_in_cookie = try
            cookies = HTTP.cookies(request)
            any(cookies) do cookie
                cookie.name == "secret" && cookie.value == session.secret
            end
        catch e
            @warn "Failed to authenticate request using cookies" exception = (e, catch_backtrace())
            false
        end
    )
    # that ) || ( kind of looks like Krabs from spongebob
end

function http_router_for(session::ServerSession)
    router = HTTP.Router()
    security = session.options.security

    function add_set_secret_cookie!(response::HTTP.Response)
        push!(response.headers, "Set-Cookie" => "secret=$(session.secret); SameSite=Strict; HttpOnly")
        response
    end

    """
        with_authentication(f::Function)

    Returns a function `HTTP.Request → HTTP.Response` which does three things:
    1. Check whether the request is authenticated (by calling `is_authenticated`), if not, return a 403 error.
    2. Call your `f(request)` to create the response message.
    3. Add a `Set-Cookie` header to the response with the session's `secret`.
    """
    function with_authentication(f::Function; required::Bool)::Function
        return function (request::HTTP.Request)
            if !required || is_authenticated(session, request)
                response = f(request)
                add_set_secret_cookie!(response)
                if !required
                    filter!(p -> p[1] != "Access-Control-Allow-Origin", response.headers)
                    push!(response.headers, "Access-Control-Allow-Origin" => "*")
                end
                response
            else
                error_response(403, "Not yet authenticated", "Open the link that was printed in the terminal where you launched Pluto. It includes a secret, which is needed to access this server.

If you are running the server yourself and want to change this configuration, have a look at the keyword arguments to Pluto.run.

Please report this error if you did not expect it!")
            end
        end
    end

    function create_serve_onefile(path)
        return request::HTTP.Request -> asset_response(normpath(path))
    end

    # / does not need security.require_secret_for_open_links
    # because this is how we handle the case:
    #    require_secret_for_open_links == true
    #    require_secret_for_access == false
    # Access to all 'risky' endpoints is still restricted to requests that have the secret cookie, but visiting `/` is allowed, and it will set the cookie. From then on the security situation is identical to
    #    secret_for_access == true
    HTTP.@register(router, "GET", "/", with_authentication(
        create_serve_onefile(project_relative_path(frontend_directory(), "index.html"));
        required=security.require_secret_for_access
        ))
    HTTP.@register(router, "GET", "/edit", with_authentication(
        create_serve_onefile(project_relative_path(frontend_directory(), "editor.html"));
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links,
    ))
    # the /edit page also uses with_authentication, but this is not how access to notebooks is secured: this is done by requiring the WS connection to be authenticated.
    # we still use it for /edit to do the cookie stuff, and show a more helpful error, instead of the WS never connecting.

    HTTP.@register(router, "GET", "/ping", r -> HTTP.Response(200, "OK!"))
    HTTP.@register(router, "GET", "/possible_binder_token_please", r -> session.binder_token === nothing ? HTTP.Response(200,"") : HTTP.Response(200, session.binder_token))

    function try_launch_notebook_response(action::Function, path_or_url::AbstractString; title="", advice="", home_url="./", as_redirect=true, action_kwargs...)
        try
            nb = action(session, path_or_url; action_kwargs...)
            notebook_response(nb; home_url, as_redirect)
        catch e
            if e isa SessionActions.NotebookIsRunningException
                notebook_response(e.notebook; home_url, as_redirect)
            else
                error_response(500, title, advice, sprint(showerror, e, stacktrace(catch_backtrace())))
            end
        end
    end

    serve_newfile = with_authentication(;
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links
    ) do request::HTTP.Request
        notebook_response(SessionActions.new(session); as_redirect=(request.method == "GET"))
    end
    HTTP.@register(router, "GET", "/new", serve_newfile)
    HTTP.@register(router, "POST", "/new", serve_newfile)

    # This is not in Dynamic.jl because of bookmarks, how HTML works,
    # real loading bars and the rest; Same for CustomLaunchEvent
    serve_openfile = with_authentication(;
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links
    ) do request::HTTP.Request
        try
            uri = HTTP.URI(request.target)
            query = HTTP.queryparams(uri)
            as_sample = haskey(query, "as_sample")
            if haskey(query, "path")
                path = tamepath(query["path"])
                if isfile(path)
                    return try_launch_notebook_response(SessionActions.open, path, as_redirect=(request.method == "GET"), as_sample=as_sample, title="Failed to load notebook", advice="The file $(htmlesc(path)) could not be loaded. Please report this error!")
                else
                    return error_response(404, "Can't find a file here", "Please check whether $(htmlesc(path)) exists.")
                end
            elseif haskey(query, "url")
                url = query["url"]
                return try_launch_notebook_response(SessionActions.open_url, url, as_redirect=(request.method == "GET"), as_sample=as_sample, title="Failed to load notebook", advice="The notebook from $(htmlesc(url)) could not be loaded. Please report this error!")
            else
                # You can ask Pluto to handle CustomLaunch events
                # and do some magic with how you open files.
                # You are responsible to keep this up to date.
                # See Events.jl for types and explanation
                #
                maybe_notebook_response = try_event_call(session, CustomLaunchEvent(query, request, try_launch_notebook_response))
                isnothing(maybe_notebook_response) && return error("Empty request")
                return maybe_notebook_response
            end
        catch e
            return error_response(400, "Bad query", "Please report this error!", sprint(showerror, e, stacktrace(catch_backtrace())))
        end
    end

    HTTP.@register(router, "GET", "/open", serve_openfile)
    HTTP.@register(router, "POST", "/open", serve_openfile)

    serve_sample = with_authentication(;
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links
    ) do request::HTTP.Request
        uri = HTTP.URI(request.target)
        sample_filename = split(HTTP.unescapeuri(uri.path), "sample/")[2]
        sample_path = project_relative_path("sample", sample_filename)

        try_launch_notebook_response(SessionActions.open, sample_path; as_redirect=(request.method == "GET"), home_url="../", as_sample=true, title="Failed to load sample", advice="Please report this error!")
    end
    HTTP.@register(router, "GET", "/sample/*", serve_sample)
    HTTP.@register(router, "POST", "/sample/*", serve_sample)

    notebook_from_uri(request) = let
        uri = HTTP.URI(request.target)
        query = HTTP.queryparams(uri)
        id = UUID(query["id"])
        session.notebooks[id]
    end
    serve_notebookfile = with_authentication(;
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links
    ) do request::HTTP.Request
        try
            notebook = notebook_from_uri(request)
            response = HTTP.Response(200, sprint(save_notebook, notebook))
            push!(response.headers, "Content-Type" => "text/julia; charset=utf-8")
            push!(response.headers, "Content-Disposition" => "inline; filename=\"$(basename(notebook.path))\"")
            response
        catch e
            return error_response(400, "Bad query", "Please report this error!", sprint(showerror, e, stacktrace(catch_backtrace())))
        end
    end
    HTTP.@register(router, "GET", "/notebookfile", serve_notebookfile)

    serve_statefile = with_authentication(;
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links
    ) do request::HTTP.Request
        try
            notebook = notebook_from_uri(request)
            response = HTTP.Response(200, Pluto.pack(Pluto.notebook_to_js(notebook)))
            push!(response.headers, "Content-Type" => "application/octet-stream")
            push!(response.headers, "Content-Disposition" => "inline; filename=\"$(without_pluto_file_extension(basename(notebook.path))).plutostate\"")
            response
        catch e
            return error_response(400, "Bad query", "Please report this error!", sprint(showerror, e, stacktrace(catch_backtrace())))
        end
    end
    HTTP.@register(router, "GET", "/statefile", serve_statefile)

    serve_notebookexport = with_authentication(;
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links
    ) do request::HTTP.Request
        try
            notebook = notebook_from_uri(request)
            response = HTTP.Response(200, generate_html(notebook))
            push!(response.headers, "Content-Type" => "text/html; charset=utf-8")
            push!(response.headers, "Content-Disposition" => "inline; filename=\"$(basename(notebook.path)).html\"")
            response
        catch e
            return error_response(400, "Bad query", "Please report this error!", sprint(showerror, e, stacktrace(catch_backtrace())))
        end
    end
    HTTP.@register(router, "GET", "/notebookexport", serve_notebookexport)

    serve_notebookupload = with_authentication(;
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links
    ) do request::HTTP.Request
        save_path = SessionActions.save_upload(request.body)
        try_launch_notebook_response(
            SessionActions.open,
            save_path,
            as_redirect=false,
            as_sample=false,
            title="Failed to load notebook",
            advice="Make sure that you copy the entire notebook file. Please report this error!"
        )
    end
    HTTP.@register(router, "POST", "/notebookupload", serve_notebookupload)

    function serve_asset(request::HTTP.Request)
        uri = HTTP.URI(request.target)
        filepath = project_relative_path(frontend_directory(), relpath(HTTP.unescapeuri(uri.path), "/"))
        asset_response(filepath; cacheable=should_cache(filepath))
    end
    HTTP.@register(router, "GET", "/*", serve_asset)
    HTTP.@register(router, "GET", "/favicon.ico", create_serve_onefile(project_relative_path(frontend_directory(allow_bundled=false), "img", "favicon.ico")))

    return router
end


text:
##############################################################################
#
# TF
#
##############################################################################

tf(dtm::Matrix{T}) where {T <: Real} = tf!(dtm, Array{Float64}(undef, size(dtm)...))

tf(dtm::SparseMatrixCSC{T}) where {T <: Real} =  tf!(dtm, similar(dtm, Float64))

tf!(dtm::AbstractMatrix{T}) where {T <: Real} = tf!(dtm, dtm)

tf!(dtm::SparseMatrixCSC{T}) where {T <: Real} = tf!(dtm, dtm)

tf(dtm::DocumentTermMatrix) = tf(dtm.dtm)

# The second Matrix will be overwritten with the result
# Will work correctly if dtm and tfidf are the same matrix
function tf!(dtm::AbstractMatrix{T1}, tf::AbstractMatrix{T2}) where {T1 <: Real, T2 <: AbstractFloat}
    n, p = size(dtm)

    # TF tells us what proportion of a document is defined by a term
    for i in 1:n
        words_in_document = 0
        for j in 1:p
            words_in_document += dtm[i, j]
        end
        tf[i, :] = dtm[i, :] ./ max(words_in_document, one(T1))
    end

    return tf
end

# assumes second matrix has same nonzeros as first one
function tf!(dtm::SparseMatrixCSC{T}, tf::SparseMatrixCSC{F}) where {T <: Real, F <: AbstractFloat}
    rows = rowvals(dtm)
    dtmvals = nonzeros(dtm)
    tfvals = nonzeros(tf)
    @assert size(dtmvals) == size(tfvals)

    # TF tells us what proportion of a document is defined by a term
    words_in_documents = sum(dtm,dims=2)

    n, p = size(dtm)
    for i = 1:p
       for j in nzrange(dtm, i)
          row = rows[j]
          tfvals[j] = dtmvals[j] / max(words_in_documents[row], one(T))
       end
    end
    tf
end

##############################################################################
#
# TF-IDF
#
##############################################################################

tf_idf(dtm::Matrix{T}) where {T <: Real} = tf_idf!(dtm, Array{Float64}(undef, size(dtm)...))

tf_idf(dtm::SparseMatrixCSC{T}) where {T <: Real} =  tf_idf!(dtm, similar(dtm, Float64))

tf_idf!(dtm::AbstractMatrix{T}) where {T <: Real} = tf_idf!(dtm, dtm)

tf_idf!(dtm::SparseMatrixCSC{T}) where {T <: Real} = tf_idf!(dtm, dtm)

tf_idf(dtm::DocumentTermMatrix) = tf_idf(dtm.dtm)

# This does not make sense, since DocumentTermMatrix is based on an array of integers
#tf_idf!(dtm::DocumentTermMatrix) = tf_idf!(dtm.dtm)

# The second Matrix will be overwritten with the result
# Will work correctly if dtm and tfidf are the same matrix
function tf_idf!(dtm::AbstractMatrix{T1}, tfidf::AbstractMatrix{T2}) where {T1 <: Real, T2 <: AbstractFloat}
    n, p = size(dtm)

    # TF tells us what proportion of a document is defined by a term
    tf!(dtm, tfidf)

    # IDF tells us how rare a term is in the corpus
    documents_containing_term = vec(sum(dtm .> 0, dims=1))
    idf = log.(n ./ documents_containing_term)

    # TF-IDF is the product of TF and IDF
    for i in 1:n
        for j in 1:p
           tfidf[i, j] = tfidf[i, j] * idf[j]
        end
    end

    return tfidf
end

# sparse version
function tf_idf!(dtm::SparseMatrixCSC{T}, tfidf::SparseMatrixCSC{F}) where {T <: Real, F <: AbstractFloat}
    rows = rowvals(dtm)
    dtmvals = nonzeros(dtm)
    tfidfvals = nonzeros(tfidf)
    @assert size(dtmvals) == size(tfidfvals)

    n, p = size(dtm)

    # TF tells us what proportion of a document is defined by a term
    words_in_documents = F.(sum(dtm, dims=2))
    oneval = one(F)

    # IDF tells us how rare a term is in the corpus
    documents_containing_term = vec(sum(dtm .> 0, dims=1))
    idf = log.(n ./ documents_containing_term)

    for i = 1:p
       for j in nzrange(dtm, i)
          row = rows[j]
          tfidfvals[j] = dtmvals[j] / max(words_in_documents[row], oneval) * idf[i]
       end
    end

    tfidf
end


text:
using CairoMakie, Random
Random.seed!(145)
n = 10
x = LinRange(1,20,n)
y = 5*rand(n)
yerr = 0.4*abs.(randn(n))
set_theme!(theme_dark())
fig, ax,  = barplot(x, y; color = 1:n, width = 1,
    colormap = Reverse(:inferno), strokewidth = 1, strokecolor = :white,
    figure = (;resolution = (600,400)))
errorbars!(x, y, yerr, color = 1:n, whiskerwidth = 10, linewidth = 2,
    colormap = Reverse(:linear_wcmr_100_45_c42_n256))
ax.xticks = (x, string.('A':'J'))
hideydecorations!(ax; grid = false)
save("barplot_and_errorbars.png", fig, px_per_unit = 2)
display(fig)


text:
# All the types supported by mshadow.
typealias DType Union{Float32, Float64, Float16, UInt8, Int32}
@enum TypeFlag kFloat32 kFloat64 kFloat16 kUint8 kInt32
typealias DEFAULT_DTYPE Float32

function toTypeFlag{T <: DType}(:: Type{T})
  if T == Float32
    return kFloat32
  elseif T == Float64
    return kFloat64
  elseif T == Float16
    return kFloat16
  elseif T == UInt8
    return kUint8
  elseif T == Int32
    return kInt32
  else
    throw(ArgumentError("Can't convert $T to DType."))
  end
end

function fromTypeFlag(T :: TypeFlag)
  if T == kFloat32
    return Float32
  elseif T == kFloat64
    return Float64
  elseif T == kFloat16
    return Float16
  elseif T == kUint8
    return UInt8
  elseif T == kInt32
    return Int32
  else
    throw(ArgumentError("Can't convert DType $T."))
  end
end

# create a NDArray handle of specific shape
function _ndarray_alloc{N}(shape :: NTuple{N, Int}, ctx :: Context, delay_alloc :: Bool)
  h_ref  = Ref{MX_handle}(0)
  shape  = flipdim(MX_uint[shape...],1)
  @mxcall(:MXNDArrayCreate, (Ptr{MX_uint}, MX_uint, Cint, Cint, Cint, Ref{MX_handle}),
      shape, length(shape), ctx.device_type, ctx.device_id, delay_alloc, h_ref)
  handle = MX_NDArrayHandle(h_ref[])
  return handle
end

# create a NDArray handle of specific shape type
function _ndarray_alloc{T <: DType,N}(:: Type{T}, shape :: NTuple{N, Int}, ctx :: Context, delay_alloc :: Bool)
  h_ref  = Ref{MX_handle}(0)
  shape  = flipdim(MX_uint[shape...],1)
  dtype  = toTypeFlag(T)
  @mxcall(:MXNDArrayCreateEx, (Ptr{MX_uint}, MX_uint, Cint, Cint, Cint, Cint, Ref{MX_handle}),
      shape, length(shape), ctx.device_type, ctx.device_id, delay_alloc, dtype, h_ref)
  handle = MX_NDArrayHandle(h_ref[])
  return handle
end

# create a handle to an empty NDArray, this handle can be used to hold
# results returned by libmx API calls
function _ndarray_alloc()
  h_ref = Ref{MX_handle}(0)
  @mxcall(:MXNDArrayCreateNone, (Ref{MX_handle},), h_ref)
  return MX_NDArrayHandle(h_ref[])
end

################################################################################
# NDArray Type
################################################################################
"""
    NDArray

Wrapper of the `NDArray` type in `libmxnet`. This is the basic building block
of tensor-based computation.

!!! note
      since C/C++ use row-major ordering for arrays while Julia follows a
      column-major ordering. To keep things consistent, we keep the underlying data
      in their original layout, but use *language-native* convention when we talk
      about shapes. For example, a mini-batch of 100 MNIST images is a tensor of
      C/C++/Python shape (100,1,28,28), while in Julia, the same piece of memory
      have shape (28,28,1,100).
"""
type NDArray
  handle   :: MX_NDArrayHandle
  writable :: Bool

  function NDArray(handle, writable=true)
    new(handle, writable)
  end
end

function Base.show(io :: IO, arr :: NDArray)
  print(io, "mx.NDArray{$(eltype(arr))}$(size(arr))")
end

function NDArray{T<:Real}(data :: Array{T})
  copy(data, cpu())
end

function Base.unsafe_convert(::Type{MX_handle}, obj::NDArray)
  Base.unsafe_convert(MX_handle, obj.handle)
end
Base.convert(t::Type{MX_handle}, obj::NDArray) = Base.unsafe_convert(t, obj)
Base.cconvert(t::Type{MX_handle}, obj::NDArray) = Base.unsafe_convert(t, obj)

################################################################################
# NDArray functions exported to the users
################################################################################
"""
    context(arr :: NDArray)

Get the context that this `NDArray` lives on.
"""
function context(arr :: NDArray)
  ref_typeid = Ref{Cint}(0)
  ref_devid  = Ref{Cint}(0)
  @mxcall(:MXNDArrayGetContext, (MX_handle, Ref{Cint}, Ref{Cint}),
          arr, ref_typeid, ref_devid)
  return Context(ref_typeid[], ref_devid[])
end

"""
    empty(DType, shape :: Tuple, ctx :: Context)
    empty(DType, shape :: Tuple)
    empty(DType, dim1, dim2, ...)

Allocate memory for an uninitialized `NDArray` with a specified type.
"""
function empty{N,T<:DType}(::Type{T}, shape :: NTuple{N, Int})
  empty(T, shape, cpu())
end
function empty{N,T<:DType}(:: Type{T}, shape :: NTuple{N, Int}, ctx :: Context)
  NDArray(_ndarray_alloc(T, shape, ctx, false))
end
function empty{T<:DType}(:: Type{T}, shape :: Int...)
  empty(T, shape)
end

"""
    empty(shape :: Tuple, ctx :: Context)
    empty(shape :: Tuple)
    empty(dim1, dim2, ...)

Allocate memory for an uninitialized `NDArray` with specific shape of type Float32.
"""
function empty{N}(shape :: NTuple{N, Int})
  empty(shape, cpu())
end
function empty{N}(shape :: NTuple{N, Int}, ctx :: Context)
  NDArray(_ndarray_alloc(shape, ctx, false))
end
function empty(shape :: Int...)
  empty(shape)
end

"""
    zeros(DType, shape :: Tuple, ctx :: Context)
    zeros(DType, shape :: Tuple)
    zeros(DType, dim1, dim2, ...)

Create zero-ed `NDArray` with specific shape and type
"""
function zeros{N,T<:DType}(:: Type{T}, shape :: NTuple{N, Int})
  zeros(T, shape, cpu())
end
function zeros{N,T<:DType}(:: Type{T}, shape :: NTuple{N, Int}, ctx :: Context)
  arr = empty(T, shape, ctx)
  arr[:] = zero(T)
  return arr
end
function zeros{T<:DType}(:: Type{T}, shape :: Int...)
  zeros(T, shape)
end

"""
    zeros(shape :: Tuple, ctx :: Context)
    zeros(shape :: Tuple)
    zeros(dim1, dim2, ...)

Create zero-ed `NDArray` with specific shape.
"""
function zeros{N}(shape :: NTuple{N, Int})
  zeros(shape, cpu())
end
function zeros{N}(shape :: NTuple{N, Int}, ctx :: Context)
  arr = empty(shape, ctx)
  arr[:] = 0
  return arr
end
function zeros(shape :: Int...)
  zeros(shape)
end

"""
    ones(DType, shape :: Tuple, ctx :: Context)
    ones(DType, shape :: Tuple)
    ones(DType, dim1, dim2, ...)

Create an `NDArray` with specific shape & type, and initialize with 1.
"""
function ones{N,T<:DType}(:: Type{T}, shape :: NTuple{N, Int})
  ones(T, shape, cpu())
end
function ones{N,T<:DType}(:: Type{T}, shape :: NTuple{N, Int}, ctx :: Context)
  arr = empty(T, shape, ctx)
  arr[:] = one(T)
  return arr
end
function ones{T<:DType}(:: Type{T}, shape :: Int...)
  ones(T, shape)
end

"""
    ones(shape :: Tuple, ctx :: Context)
    ones(shape :: Tuple)
    ones(dim1, dim2, ...)

Create an `NDArray` with specific shape and initialize with 1.
"""
function ones{N}(shape :: NTuple{N, Int})
  ones(shape, cpu())
end
function ones{N}(shape :: NTuple{N, Int}, ctx :: Context)
  arr = empty(shape, ctx)
  arr[:] = 1
  return arr
end
function ones(shape :: Int...)
  ones(shape)
end

import Base: size, length, ndims, eltype

"""
    size(arr :: NDArray)
    size(arr :: NDArray, dim :: Int)

Get the shape of an `NDArray`. The shape is in Julia's column-major convention. See
also the notes on NDArray shapes [`NDArray`](@ref).
"""
function size(arr :: NDArray)
  ref_ndim  = Ref{MX_uint}(0)
  ref_shape = Ref{Ptr{MX_uint}}(0)
  @mxcall(:MXNDArrayGetShape, (MX_handle, Ref{MX_uint}, Ref{Ptr{MX_uint}}),
          arr, ref_ndim, ref_shape)
  tuple(map(Int, flipdim(unsafe_wrap(Array, ref_shape[], ref_ndim[]),1))...)
end
function size(arr :: NDArray, dim :: Int)
  size(arr)[dim]
end

"""
    length(arr :: NDArray)

Get the number of elements in an `NDArray`.
"""
function length(arr :: NDArray)
  prod(size(arr))
end

"""
    ndims(arr :: NDArray)

Get the number of dimensions of an `NDArray`. Is equivalent to `length(size(arr))`.
"""
function ndims(arr :: NDArray)
  length(size(arr))
end

"""
    eltype(arr :: NDArray)

Get the element type of an `NDArray`.
"""
function eltype{T <: Union{NDArray, MX_NDArrayHandle}}(arr :: T)
  dtype_ref = Ref{Cint}(0)
  @mxcall(:MXNDArrayGetDType, (MX_handle, Ptr{Cint}), arr, dtype_ref)

  if dtype_ref[] == -1 # arr->is_none()
    warn("Eltype of $arr is not defined")
    Base.show_backtrace(STDOUT,backtrace())
    println()
    return Float32
  else
    return fromTypeFlag(TypeFlag(dtype_ref[]))
  end
end

import Base: slice
"""
    slice(arr :: NDArray, start:stop)

Create a view into a sub-slice of an `NDArray`. Note only slicing at the slowest
changing dimension is supported. In Julia's column-major perspective, this is the last
dimension. For example, given an `NDArray` of shape (2,3,4), `slice(array, 2:3)` will create
a `NDArray` of shape (2,3,2), sharing the data with the original array. This operation is
used in data parallelization to split mini-batch into sub-batches for different devices.
"""
function slice(arr :: NDArray, ::Colon)
  arr
end
function slice(arr :: NDArray, slice::UnitRange{Int})
  dim1 = size(arr)[end]
  @assert(1 <= slice.start <= slice.stop <= dim1)
  if slice.start == 1 && slice.stop == dim1
    return arr
  end

  hdr_ref = Ref{MX_handle}(0)
  # note Julia is 1-based, inclusive-inclusive indexing, while C++ is
  # 0-based, inclusive-exclusive indexing. So 1:3 in Julia should
  # translates into 0:3 in C++.
  @mxcall(:MXNDArraySlice, (MX_handle, MX_uint, MX_uint, Ref{MX_handle}),
          arr, slice.start-1, slice.stop, hdr_ref)
  return NDArray(MX_NDArrayHandle(hdr_ref[]), arr.writable)
end

import Base: setindex!

"""
    setindex!(arr :: NDArray, val, idx)

Assign values to an `NDArray`. Elementwise assignment is not implemented, only the following
scenarios are supported

* `arr[:] = val`: whole array assignment, `val` could be a scalar or an array (Julia `Array`
  or `NDArray`) of the same shape.
* `arr[start:stop] = val`: assignment to a *slice*, `val` could be a scalar or an array of
  the same shape to the slice. See also [`slice`](@ref).
"""
function setindex!(arr :: NDArray, val :: Real, ::Colon)
  @assert(arr.writable)
  _set_value(convert(eltype(arr), val), arr)
  return arr
end
function setindex!{T<:Real}(arr :: NDArray, val :: Array{T}, ::Colon)
  copy!(arr, val)
end
function setindex!(arr :: NDArray, val :: NDArray, ::Colon)
  copy!(arr, val)
end
function setindex!{T<:Real}(arr :: NDArray, val :: Union{T,Array{T},NDArray}, idx::UnitRange{Int})
  setindex!(slice(arr, idx), val, Colon())
end

import Base: getindex
"""
    getindex(arr :: NDArray, idx)

Shortcut for [`slice`](@ref). A typical use is to write

```julia
  arr[:] += 5
```

which translates into

```julia
  arr[:] = arr[:] + 5
```

which furthur translates into

```julia
  setindex!(getindex(arr, Colon()), 5, Colon())
```

!!! note
    The behavior is quite different from indexing into Julia's `Array`. For example, `arr[2:5]`
    create a **copy** of the sub-array for Julia `Array`, while for `NDArray`, this is
    a *slice* that shares the memory.
"""
function getindex(arr :: NDArray, ::Colon)
  return arr
end

"""
Shortcut for [`slice`](@ref). **NOTE** the behavior for Julia's built-in index slicing is to create a
copy of the sub-array, while here we simply call `slice`, which shares the underlying memory.
"""
function getindex(arr :: NDArray, idx::UnitRange{Int})
  slice(arr, idx)
end

import Base: copy!, copy, convert
"""
    copy!(dst :: Union{NDArray, Array}, src :: Union{NDArray, Array})

Copy contents of `src` into `dst`.
"""
function copy!(dst :: NDArray, src :: NDArray)
  @assert(dst.writable)
  if dst.handle == src.handle
    warn("Copying an NDArray to itself")
    return
  end

  _copyto(src, dst)
  return dst
end

function copy!{T<:DType}(dst :: Array{T}, src :: NDArray)
  @assert T == eltype(src)
  @assert size(dst) == size(src)
  @mxcall(:MXNDArraySyncCopyToCPU, (MX_handle, Ptr{Void}, Csize_t),
          src, pointer(dst), length(dst))
  return dst
end
function copy!{T<:Real}(dst :: Array{T}, src :: NDArray)
  copy!(dst, copy(src))
end

function copy!{T<:Real}(dst :: NDArray, src :: Array{T})
  @assert dst.writable
  @assert size(dst) == size(src)
  src = convert(Array{eltype(dst)}, src) # this might involve copying
  @mxcall(:MXNDArraySyncCopyFromCPU, (MX_handle, Ptr{Void}, Csize_t),
          dst.handle, pointer(src), length(src))
  return dst
end

function copy_ignore_shape!{T<:Real}(dst :: NDArray, src :: Array{T})
  @assert dst.writable
  @assert length(dst) == length(src)
  src = convert(Array{eltype(dst)}, src) # this might involve copying
  @mxcall(:MXNDArraySyncCopyFromCPU, (MX_handle, Ptr{Void}, Csize_t),
          dst.handle, pointer(src), length(src))
  return dst
end

"""
    copy(arr :: NDArray)
    copy(arr :: NDArray, ctx :: Context)
    copy(arr :: Array, ctx :: Context)

Create a copy of an array. When no `Context` is given, create a Julia `Array`.
Otherwise, create an `NDArray` on the specified context.
"""
# Create copy: NDArray -> Julia Array
function copy(arr :: NDArray)
  j_arr = Array{eltype(arr)}(size(arr))
  copy!(j_arr, arr)
end

# Create copy: NDArray -> NDArray in a given context
function copy(arr :: NDArray, ctx :: Context)
  dst = NDArray(_ndarray_alloc(eltype(arr), size(arr), ctx, true))
  copy!(dst, arr)
end

# Create copy: Julia Array -> NDArray in a given context
function copy{T<:DType}(arr :: Array{T}, ctx :: Context)
  dst = empty(T, size(arr), ctx)
  copy!(dst, arr)
end

"""
    convert(::Type{Array{T}}, arr :: NDArray)

Convert an `NDArray` into a Julia `Array` of specific type. Data will be copied.
"""
function convert{T<:Real}(t::Type{Array{T}}, arr :: NDArray)
  convert(t, copy(arr))
end

"""
    @inplace

Julia does not support re-definiton of `+=` operator (like `__iadd__` in python),
When one write `a += b`, it gets translated to `a = a+b`. `a+b` will allocate new
memory for the results, and the newly allocated `NDArray` object is then assigned
back to a, while the original contents in a is discarded. This is very inefficient
when we want to do inplace update.

This macro is a simple utility to implement this behavior. Write

```julia
  @mx.inplace a += b
```

will translate into

```julia
  mx.add_to!(a, b)
```

which will do inplace adding of the contents of `b` into `a`.
"""
macro inplace(stmt)
  if stmt.head == :+= || stmt.head == :.+=
    Expr(:call, :add_to!, esc(stmt.args[1]), esc(stmt.args[2]))
  elseif stmt.head == :-= || stmt.head == :.-=
    Expr(:call, :sub_from!, esc(stmt.args[1]), esc(stmt.args[2]))
  elseif stmt.head == :.*=
    Expr(:call, :mul_to!, esc(stmt.args[1]), esc(stmt.args[2]))
  elseif stmt.head == :./=
    Expr(:call, :div_from!, esc(stmt.args[1]), esc(stmt.args[2]))
  else
    error("unsupported inplace translation for $stmt")
  end
end

"""
    add_to!(dst :: NDArray, args :: Union{Real, NDArray}...)

Add a bunch of arguments into `dst`. Inplace updating.
"""
function add_to!(dst :: NDArray, args :: Union{Real, NDArray}...)
  @assert dst.writable
  for arg in args
    if isa(arg, Real)
      _plus_scalar(dst, convert(eltype(dst), arg), dst)
    else
      _plus(dst, arg, dst)
    end
  end
  return dst
end

import Base: +, .+

"""
    +(args...)
    .+(args...)

Summation. Multiple arguments of either scalar or `NDArray` could be
added together. Note at least the first or second argument needs to be an `NDArray` to
avoid ambiguity of built-in summation.
"""
function +(arg0 :: NDArray, args :: Union{Real, NDArray}...)
  ret = copy(arg0, context(arg0))
  add_to!(ret, args...)
end
function .+(arg0 :: NDArray, args :: Union{Real, NDArray}...)
  +(arg0, args...)
end
function +(arg0 :: Real, arg1 :: NDArray, args :: Union{Real, NDArray}...)
  +(arg1, arg0, args...)
end
function .+(arg0 :: Real, arg1 :: NDArray, args :: Union{Real, NDArray}...)
  .+(arg1, arg0, args...)
end

"""
    sub_from!(dst :: NDArray, args :: Union{Real, NDArray}...)

Subtract a bunch of arguments from `dst`. Inplace updating.
"""
function sub_from!(dst :: NDArray, arg :: Union{Real, NDArray})
  @assert dst.writable
  if isa(arg, Real)
    _minus_scalar(dst, convert(eltype(dst), arg), dst)
  else
    _minus(dst, arg, dst)
  end
end

import Base: -, .-

"""
    -(arg0, arg1)
    -(arg0)
    .-(arg0, arg1)

Subtraction `arg0 - arg1`, of scalar types or `NDArray`. Or create
the negative of `arg0`.
"""
function -(arg0 :: NDArray, arg1 :: Union{Real, NDArray})
  ret = copy(arg0, context(arg0))
  sub_from!(ret, arg1)
end
function .-(arg0 :: NDArray, arg1 :: Union{Real, NDArray})
  -(arg0, arg1)
end
function -(arg0 :: Real, arg1 :: NDArray)
  ret = -arg1
  add_to!(ret, arg0)
  return ret
end
function .-(arg0 :: Real, arg1 :: NDArray)
  -(arg0, arg1)
end

function -(arg0 :: NDArray)
  _mul_scalar(arg0, -one(eltype(arg0)))
end

"""
    mul_to!(dst :: NDArray, arg :: Union{Real, NDArray})

Elementwise multiplication into `dst` of either a scalar or an `NDArray` of the same shape.
Inplace updating.
"""
function mul_to!(dst :: NDArray, arg :: Union{Real, NDArray})
  @assert dst.writable
  if isa(arg, Real)
    _mul_scalar(dst, convert(eltype(dst), arg), dst)
  else
    _mul(dst, arg, dst)
  end
  return dst
end

import Base: .*, *

"""
    .*(arg0, arg1)

Elementwise multiplication of `arg0` and `arg`, could be either scalar or `NDArray`.
"""
function .*(arg0 :: NDArray, arg :: Union{Real, NDArray})
  ret = copy(arg0, context(arg0))
  mul_to!(ret, arg)
end
function .*(arg0 :: Real, arg :: NDArray)
  .*(arg, arg0)
end

"""
    *(arg0, arg1)

Currently only multiplication a scalar with an `NDArray` is implemented. Matrix multiplication
is to be added soon.
"""
function *(arg0 :: NDArray, arg :: Real)
  ret = copy(arg0, context(arg0))
  mul_to!(ret, arg)
end
function *(arg0 :: Real, arg :: NDArray)
  *(arg, arg0)
end

"""
    div_from!(dst :: NDArray, arg :: Union{Real, NDArray})

Elementwise divide a scalar or an `NDArray` of the same shape from `dst`. Inplace updating.
"""
function div_from!(dst :: NDArray, arg :: Union{Real, NDArray})
  @assert dst.writable
  if isa(arg, Real)
    _div_scalar(dst, convert(eltype(dst), arg), dst)
  else
    _div(dst, arg, dst)
  end
end

import Base: ./, /
"""
    ./(arg0 :: NDArray, arg :: Union{Real, NDArray})

Elementwise dividing an `NDArray` by a scalar or another `NDArray` of the same shape.
"""
function ./(arg0 :: NDArray, arg :: Union{Real, NDArray})
  ret = copy(arg0, context(arg0))
  div_from!(ret, arg)
end

"""
    /(arg0 :: NDArray, arg :: Real)

Divide an `NDArray` by a scalar. Matrix division (solving linear systems) is not implemented yet.
"""
function /(arg0 :: NDArray, arg :: Real)
  ./(arg0, arg)
end

"""
Manipulating as Julia Arrays
----------------------------

    @nd_as_jl(captures..., statement)

A convenient macro that allows to operate `NDArray` as Julia Arrays. For example,

```julia
  x = mx.zeros(3,4)
  y = mx.ones(3,4)
  z = mx.zeros((3,4), mx.gpu())

  @mx.nd_as_jl ro=(x,y) rw=z begin
    # now x, y, z are just ordinary Julia Arrays
    z[:,1] = y[:,2]
    z[:,2] = 5
  end
```

Under the hood, the macro convert all the declared captures from `NDArray` into Julia
Arrays, by using `try_get_shared`. And automatically commit the modifications back into
the `NDArray` that is declared as `rw`. This is useful for fast prototyping and when
implement non-critical computations, such as `AbstractEvalMetric`.

!!! note
* Multiple `rw` and / or `ro` capture declaration could be made.
* The macro does **not** check to make sure that `ro` captures are not modified. If the
  original `NDArray` lives in CPU memory, then it is very likely the corresponding
  Julia Array shares data with the `NDArray`, so modifying the Julia Array will also
  modify the underlying `NDArray`.
* More importantly, since the `NDArray` is
  asynchronized, we will wait for *writing* for `rw` variables but wait only for *reading*
  in `ro` variables. If we write into those `ro` variables, **and** if the memory is
  shared, racing condition might happen, and the behavior is undefined.
* When an `NDArray` is declared to be captured as `rw`, its contents is always sync
  back in the end.
* The execution results of the expanded macro is always `nothing`.
* The statements are wrapped in a `let`, thus locally introduced new variables will not be
  available after the statements. So you will need to declare the variables before calling the
  macro if needed.
"""
macro nd_as_jl(m_args...)
  @assert(length(m_args) > 0)
  stmts = m_args[end]
  @assert(isa(stmts, Expr) && stmts.head == :block,
          "The last argument should be a statement block (begin-end); but get $stmts")
  stmts = esc(stmts)

  dclrs  = m_args[1:end-1]
  nd_ro  = []
  nd_rw  = []
  nd_all = []
  for declr in dclrs
    @assert(isa(declr, Expr) && declr.head == :(=) && length(declr.args)==2 && declr.args[1] ∈ (:ro,:rw),
            "Invalid declaration, should be rw=(x,y) or ro=z; but get $declr")

    declr_vars = declr.args[2]
    if isa(declr_vars, Symbol)
      declr_vars = (declr_vars,)
    elseif isa(declr_vars, Expr)
      @assert(declr_vars.head ∈ (:tuple, :vect),
              "Capture declaration should be a variable or a tuple of variables; but got $declr_vars")
      declr_vars = declr_vars.args
    else
      @assert(false, "Capture declaration should be a variable or a tuple of variables; but got $declr_vars")
    end
    for declr_var in declr_vars
      @assert(isa(declr_var, Symbol),
              "Captured ndarrays in ro/rw declaration should be variables, but get $(declr_var)")
    end
    append!(nd_all, [declr_vars...])
    if declr.args[1] == :ro
      append!(nd_ro, [declr_vars...])
    else
      append!(nd_rw, [declr_vars...])
    end
  end

  nd_ro    = map(esc, nd_ro)
  nd_rw    = map(esc, nd_rw)
  nd_all   = map(esc, nd_all)
  rw_origs = [gensym() for _ in nd_rw]

  save_statements  = Expr(:block, [:($v_orig = $v) for (v_orig, v) in zip(rw_origs, nd_rw)]...)
  wait_statements  = Expr(:block, [:(_wait_to_read($v)) for v in nd_ro]...,
                                  [:(_wait_to_write($v)) for v in nd_rw]...)
  clear_statements = Expr(:block, [:($v_orig = nothing) for v_orig in rw_origs]...)
  let_assignments  = [:($v = try_get_shared($v)) for v in nd_all]
  sync_statements  = map(rw_origs, nd_rw) do v_orig, v
    quote
      if !is_shared($v, $v_orig)
        # copy data back if not or no longer sharing data
        copy!($v_orig, $v)
      end
    end
  end
  sync_statements  = Expr(:block, sync_statements...)

  let_statement = Expr(:let, quote
    $stmts
    $sync_statements
  end, let_assignments...)
  m_body = quote
    $wait_statements
    $save_statements
    $let_statement
    $clear_statements
    nothing # the final results is always nothing
  end

  m_body
end

# NOTE: internal use only. Accessing pointers on a different device (e.g. accessing GPU
# pointers from CPU) leads to undefined behavior.
import Base.pointer
function pointer(arr :: NDArray)
  pdata = Ref{Ptr{MX_float}}(0)
  @mxcall(:MXNDArrayGetData, (MX_handle, Ref{Ptr{MX_float}}), arr, pdata)
  return convert(Ptr{eltype(arr)}, pdata[])
end
function _wait_to_read(arr :: NDArray)
  @mxcall(:MXNDArrayWaitToRead, (MX_handle,), arr)
end
function _wait_to_write(arr :: NDArray)
  @mxcall(:MXNDArrayWaitToWrite, (MX_handle,), arr)
end

"""
    try_get_shared(arr)

Try to create a Julia array by sharing the data with the underlying `NDArray`.

# Arguments:
* `arr::NDArray`: the array to be shared.

!!! note
    The returned array does not guarantee to share data with the underlying `NDArray`.
    In particular, data sharing is possible only when the `NDArray` lives on CPU.
"""
function try_get_shared(arr :: NDArray)
  if context(arr).device_type == CPU
    # try to do data sharing
    return unsafe_wrap(Array, pointer(arr), size(arr))
  else
    # impossible to share, just copying
    return copy(arr)
  end
end

"""
    is_shared(j_arr, arr)

Test whether `j_arr` is sharing data with `arr`.

# Arguments:
* Array j_arr: the Julia Array.
* NDArray arr: the `NDArray`.
"""
function is_shared(j_arr :: Array, arr :: NDArray)
  false
end
function is_shared{T<:DType}(j_arr :: Array{T}, arr :: NDArray)
  if length(j_arr) != length(arr)
    return false
  end
  if context(arr).device_type != CPU
    return false
  end
  return pointer(j_arr) == pointer(arr)
end

"""
    load(filename, ::Type{NDArray})

Load NDArrays from binary file.

# Arguments:
* `filename::String`: the path of the file to load. It could be S3 or HDFS address.

Returns either `Dict{Symbol, NDArray}` or `Vector{NDArray}`.

`filename` can point to `s3` or `hdfs` resources if the `libmxnet` is built with the
corresponding components enabled. Examples:
* `s3://my-bucket/path/my-s3-ndarray`
* `hdfs://my-bucket/path/my-hdfs-ndarray`
* `/path-to/my-local-ndarray`
"""
function load(filename::AbstractString, ::Type{NDArray})
  out_size      = Ref{MX_uint}(0)
  out_hdrs      = Ref{Ptr{MX_handle}}(0)
  out_name_size = Ref{MX_uint}(0)
  out_names     = Ref{char_pp}(0)
  @mxcall(:MXNDArrayLoad, (char_p, Ref{MX_uint}, Ref{Ptr{MX_handle}}, Ref{MX_uint}, Ref{char_pp}),
          filename, out_size, out_hdrs, out_name_size, out_names)
  out_name_size = out_name_size[]
  out_size      = out_size[]
  if out_name_size == 0
    return [NDArray(MX_NDArrayHandle(hdr)) for hdr in unsafe_wrap(Array, out_hdrs[], out_size)]
  else
    @assert out_size == out_name_size
    return Dict([(Symbol(unsafe_wrap(String, k)), NDArray(MX_NDArrayHandle(hdr))) for (k,hdr) in
                 zip(unsafe_wrap(Array, out_names[], out_size), unsafe_wrap(Array, out_hdrs[], out_size))])
  end
end

"""
    save(filename :: AbstractString, data)

Save NDarrays to binary file. Filename could be S3 or HDFS address, if `libmxnet` is built
with corresponding support (see `load`).

* `filename::String`: path to the binary file to write to.
* `data`: data to save to file. Data can be a`NDArray`, a `Vector{NDArray}`, or a `Dict{Base.Symbol, NDArray}`.
"""
function save(filename::String, data::NDArray)
  save(filename, [data])
end
function save(filename::String, data::Vector{NDArray})
  @mxcall(:MXNDArraySave, (char_p, MX_uint, Ptr{MX_handle}, char_pp),
          filename, length(data), MX_handle[data...], char_pp(0))
end
function save(filename::String, data::Dict{Base.Symbol,NDArray})
  names  = [k for k in keys(data)]
  arrays = MX_handle[data[k] for k in names]
  names  = String[string(k) for k in names]

  @mxcall(:MXNDArraySave, (char_p, MX_uint, Ptr{MX_handle}, char_pp),
          filename, length(names), arrays, names)
end

################################################################################
# NDArray functions dynamically imported from libmxnet
################################################################################
function _invoke_mxfunction(func_handle::MX_handle, use_vars, scalars, mut_vars; kwargs...)
  names = String[string(entry[1]) for entry in kwargs]
  args = String[string(entry[2]) for entry in kwargs]
  @mxcall(:MXFuncInvokeEx,
          (MX_handle, Ptr{MX_handle}, Ptr{MX_float}, Ptr{MX_handle}, Cint, char_pp, char_pp),
          func_handle, use_vars, scalars, mut_vars, length(names), names, args)
end

@enum(LIBMX_FUNC_TYPE_MASK,
  NDARRAY_ARG_BEFORE_SCALAR = 1,
  ACCEPT_EMPTY_MUTATE_TARGET = (1 << 2)
)

# Import corresponding math functions from base so the automatically defined libmxnet
# functions can overload them
import Base: sqrt

"""
The libxmnet APIs are automatically imported from `libmxnet.so`. The functions listed
here operate on `NDArray` objects. The arguments to the functions are typically ordered
as

```julia
  func_name(arg_in1, arg_in2, ..., scalar1, scalar2, ..., arg_out1, arg_out2, ...)
```

unless `NDARRAY_ARG_BEFORE_SCALAR` is not set. In this case, the scalars are put before the input arguments:

```julia
  func_name(scalar1, scalar2, ..., arg_in1, arg_in2, ..., arg_out1, arg_out2, ...)
```

If `ACCEPT_EMPTY_MUTATE_TARGET` is set. An overloaded function without the output arguments will also be defined:

```julia
  func_name(arg_in1, arg_in2, ..., scalar1, scalar2, ...)
```

Upon calling, the output arguments will be automatically initialized with empty NDArrays.

Those functions always return the output arguments. If there is only one output (the typical situation), that
object (`NDArray`) is returned. Otherwise, a tuple containing all the outputs will be returned.
"""

function _get_ndarray_functions()
  n = Ref{MX_uint}(0)
  handles = Ref{Ptr{MX_handle}}(0)

  @mxcall(:MXListFunctions, (Ref{MX_uint}, Ref{Ptr{MX_handle}}), n, handles)

  funcs = unsafe_wrap(Array, handles[], n[])
  return funcs
end

const _function_cache = Dict{Symbol, MX_handle}()
function _get_function(name :: Symbol)
  if !haskey(_function_cache, name)
    handle = Ref{MX_handle}(0)

    @mxcall(:MXGetFunction, (Cstring, Ref{MX_handle}), name, handle)
    _function_cache[name] = handle[]
    return handle[]
  else
    return _function_cache[name]
  end
end

function _get_function_description(handle :: MX_handle)
  # get function information (human readable)
  ref_name = Ref{char_p}(0)
  ref_desc = Ref{char_p}(0)
  ref_narg = Ref{MX_uint}(0)

  ref_arg_names = Ref{char_pp}(0)
  ref_arg_types = Ref{char_pp}(0)
  ref_arg_descs = Ref{char_pp}(0)

  ref_ret_type  = Ref{char_p}(0)

  @mxcall(:MXFuncGetInfo,
         (MX_handle, Ref{char_p}, Ref{char_p}, Ref{MX_uint}, Ref{char_pp},
          Ref{char_pp}, Ref{char_pp}, Ref{char_p}),
          handle, ref_name, ref_desc, ref_narg, ref_arg_names,
          ref_arg_types, ref_arg_descs, ref_ret_type)

  name = Symbol(unsafe_wrap(String, ref_name[]))
  signature = _format_signature(Int(ref_narg[]), ref_arg_names)
  desc = "    " * string(name) * "(" * signature * ")\n\n"
  desc *= unsafe_wrap(String, ref_desc[]) * "\n\n"
  desc *= "# Arguments\n"
  desc *= _format_docstring(Int(ref_narg[]), ref_arg_names, ref_arg_types, ref_arg_descs)
  return name, desc
end

function _get_function_expressions(handle :: MX_handle, name)
  # get function specification
  ref_n_use_vars = Ref{MX_uint}(0)
  ref_n_scalars  = Ref{MX_uint}(0)
  ref_n_mut_vars = Ref{MX_uint}(0)
  ref_type_mask  = Ref{Cint}(0)
  @mxcall(:MXFuncDescribe,
         (MX_handle, Ref{MX_uint}, Ref{MX_uint}, Ref{MX_uint}, Ref{Cint}),
          handle, ref_n_use_vars, ref_n_scalars, ref_n_mut_vars, ref_type_mask)

  n_used_vars   = ref_n_use_vars[]
  n_scalars     = ref_n_scalars[]
  n_mutate_vars = ref_n_mut_vars[]
  type_mask     = ref_type_mask[]
  accept_empty_mutate = (type_mask & convert(Cint,ACCEPT_EMPTY_MUTATE_TARGET)) != 0
  arg_before_scalar   = (type_mask & convert(Cint,NDARRAY_ARG_BEFORE_SCALAR)) != 0

  # general ndarray function
  if arg_before_scalar
    args = vcat([Expr(:(::), Symbol("in$i"), NDArray) for i=1:n_used_vars],
                [Expr(:(::), Symbol("sca$i"), Real) for i=1:n_scalars],
                [Expr(:(::), Symbol("out$i"), NDArray) for i=1:n_mutate_vars])
  else
    args = vcat([Expr(:(::), Symbol("sca$i"), Real) for i=1:n_scalars],
                [Expr(:(::), Symbol("in$i"), NDArray) for i=1:n_used_vars],
                [Expr(:(::), Symbol("out$i"), NDArray) for i=1:n_mutate_vars])
  end

  _use_vars = Expr(:ref, :MX_handle, [Symbol("in$i") for i=1:n_used_vars]...)
  _scalars  = Expr(:ref, :MX_float, [Symbol("sca$i") for i=1:n_scalars]...)
  _mut_vars = Expr(:ref, :MX_handle, [Symbol("out$i") for i=1:n_mutate_vars]...)

  # XXX: hacky way of solving the problem that the arguments of `dot` should be swapped
  # See https://github.com/dmlc/MXNet.jl/issues/55
  if name == :dot
    _use_vars.args[2:end] = flipdim(_use_vars.args[2:end], 1)
  end

  # XXX: hacky way of solving the semantic difference of the axes parameter in Julia
  # and in libmxnet.
  # See https://github.com/dmlc/MXNet.jl/pull/123
  if name == :transpose
    transform = quote
      kwargs = Any[key != :axes ? (key, arg) : (key, reverse(map(i->length(arg)-i, arg))) for (key, arg) in kwargs]
    end
  else
    transform = :()
  end

  stmt_call = quote
    local handle = _get_function($(QuoteNode(name)))
    _invoke_mxfunction(handle, $_use_vars, $_scalars, $_mut_vars; kwargs...)
  end
  if n_mutate_vars == 1
    stmt_ret = :(return out1)
  else
    stmt_ret = Expr(:return, Expr(:tuple, [Symbol("out$i") for i=1:n_mutate_vars]...))
  end

  func_def = quote
    function $name($(args...); kwargs...)
      $transform
      $stmt_call
      $stmt_ret
    end
  end

  if accept_empty_mutate
    args0      = args[1:n_used_vars+n_scalars]
    _mut_vars0 = [:(NDArray(_ndarray_alloc())) for i=1:n_mutate_vars]

    func_def0 = quote
      function $name($(args0...); kwargs...)
        $name($(args0...), $(_mut_vars0...); kwargs...)
      end
    end
    return func_def, func_def0
  else
    return func_def, :()
  end
end

macro _import_ndarray_functions()
  funcs = _get_ndarray_functions()
  func_exprs = Expr[]

  for i = 1:length(funcs)
    handle = funcs[i]

    name, desc = _get_function_description(handle)
    func_def, func_def0 = _get_function_expressions(handle, name)

    expr = quote
      $(isdefined(Base, name) ? :(import Base.$name) : :())
      @doc $desc ->
      $func_def
      $func_def0
    end

    push!(func_exprs, expr)
  end

  esc(quote
    $(func_exprs...)
  end)
end

@_import_ndarray_functions()


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    int{L} x = 7;
    final int z = m();
    static final Condition c = new Condition();

    int m() {
	int y = ({L /c T}) 42;
	return y;
    }

    public static void main(String[] args) {
	C obj = new C();
	obj.c.set();
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
using Abaco
using Test

interval=10
ts = nowts()

sn1 = "tn01"
sn2 = "tn02"
sn3 = "tn03"

count = 0
function onresult(ts, ne, name, value, inputs)
    global count
    count += 1
    @debug "timestamp [$ts]: ne: [$ne] $name = $value"
    if count == 1
        @test value.contribs == 1
        @test value.value == 250
    elseif count == 2
        @test value.contribs == 2
        @test value.value == 950
        @test value.expected == 3
    elseif count == 3
        @test value.contribs == 3
        @test value.value == 1950
    end
end

abaco = abaco_init(onresult, interval=interval)

@debug "starting demo with time span: $interval, values timestamp: $ts"

#setup_settings(abaco, "sensor", oncomplete=onresult)
setup_settings(abaco, "hub", oncomplete=onresult)

## If you have a lot of formulas use a csv file
#using DelimitedFiles
#datadir=joinpath(@__DIR__,"data")
#df = readdlm(joinpath(datadir,"formulas_co2.csv"), ';', header=true)
#add_formulas(abaco, df)

formula(abaco, "total_footprint", "sum(sensor.footprint)", "hub")

city = node(abaco, "trento", "hub")

node(abaco, city, sn1, "sensor")
node(abaco, city, sn2, "sensor")
node(abaco, city, sn3, "sensor")

current_footprint = get_collected(abaco, city.ne, "sensor.footprint")
@debug "initial carbon footprints: $current_footprint"

ingest(abaco, ts, sn1, Dict("footprint" => 250.0))
ingest(abaco, ts, sn2, Dict("footprint" => 700.0))

current_footprint = get_collected(abaco, city.ne, "sensor.footprint", ts)
@debug "current carbon footprint: $current_footprint"

footprint_sum = sum_collected(abaco, city.ne, "sensor.footprint", ts)
@debug "footprint sum: $footprint_sum"
@test footprint_sum.contribs == 2
@test footprint_sum.expected == 3
@test footprint_sum.value == 950

ingest(abaco, ts, sn3, Dict("footprint" => 1000.0))

footprint_sum = sum_collected(abaco, city.ne, "sensor.footprint", ts)
@debug "footprint sum: $footprint_sum"
@test footprint_sum.contribs == 3
@test footprint_sum.value == 1950

@test count == 3


text:
using Test
using TestSetExtensions
using Qaintessent

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "registers" begin
    @testset "quantum register construction" begin
        N = 3
        q1 = qreg(3)

        cgc = Circuit(q1)
        add!(X, q1(3))
        add!(Y, q1(1))
        add!(Z, q1(2))

        cgc_ref = Circuit{N}([
            circuit_gate(3, X),
            circuit_gate(1, Y),
            circuit_gate(2, Z)
        ])
        @test sparse_matrix(cgc) ≈ sparse_matrix(cgc_ref)
    end

    @testset "quantum register array operations " begin
        N = 3
        q1 = qreg(3)

        cgc = Circuit(q1)

        @test all(q1[:] .== [1,2,3])
        @test all(q1[1:2] .== [1,2])
    end

    @testset "quantum register wire numbers with multiple registers" begin
        N = 12
        q1 = qreg(3)
        q2 = qreg(2)
        q3 = qreg(5)
        q4 = qreg(2)
        cgc = Circuit(q1,q2,q3,q4)

        add!(X, q1(3))
        add!(Y, q1(1))
        add!(Z, q1(2))
        add_control!(Y, q2(1), q3(3))
        add_control!(Z, q1(2), q2(2))
        add_control!(TGate(), q3(2), q3(4))
        add_control!(SGate(), q1(2), q2(1))
        add!(X, q1)
        add!(SwapGate(), q2)
        add!(SwapGate(), q2(1), q4(2))
        add_control!(X, q2, q4)

        cgc_ref = Circuit{N}([
            circuit_gate(3, X),
            circuit_gate(1, Y),
            circuit_gate(2, Z),
            circuit_gate(4, Y, 8),
            circuit_gate(2, Z, 5),
            circuit_gate(7, TGate(), 9),
            circuit_gate(2, SGate(), 4),
            circuit_gate(1, X),
            circuit_gate(2, X),
            circuit_gate(3, X),
            circuit_gate(4, 5, SwapGate()),
            circuit_gate(4, 12, SwapGate()),
            circuit_gate(4, X, 11),
            circuit_gate(5, X, 12)
            ])

        @test sparse_matrix(cgc) ≈ sparse_matrix(cgc_ref)
    end

    @testset "quantum register exceptions" begin
        N = 3
        q1 = qreg(3)
        q2 = qreg(3)
        q3 = qreg(4)
        @test_throws ErrorException("Register object has yet to be used in a Circuit object") q1[1]
        @test_throws ErrorException("Register object has yet to be used in a Circuit object") q1(1)
        @test_throws ErrorException("Register object has yet to be used in a Circuit object") add!(X, q1(1))
        @test_throws ErrorException("Register object has yet to be used in a Circuit object") add!(X, q1)

        @test_throws ErrorException("Target quantum register is not included in any circuit") add_control!(X, q1, q2)
        cgc1 = Circuit(q1)
        @test_throws ErrorException("Control quantum register is not included in any circuit") add_control!(X, q1, q2)
        cgc2 = Circuit(q2)

        @test_throws ErrorException("AbstractGate SwapGate() is applied to 2 wires, quantum register of size 3 provided. Unable to determine gates to be applied.") add!(SwapGate(), q1)
        @test_throws ErrorException("Targetted Registers belong to different Circuit objects") add!(SwapGate(), q1(1), q2(1))
        @test_throws ErrorException("Target quantum register and control quantum register are used in different circuits") add_control!(X, q1, q2)
    end
end


text:
"""
    MaxPool(k::NTuple; stride = k, padding = 0)
This is a layer for max pooling with kernel size `k`.

Expects as input an array with `ndims(x) == N+2`, i.e. channel and batch dimensions, after the `N` feature dimensions, where `N = length(out)`.

The default stride is the same as kernel size `k`.

# Example
```
julia> N = NetWork(Conv((2, 2), 5=>2, relu), MaxPool(2, 2))
Layer1 : Convolution(k:(2, 2), IO:5=>2, σ:relu)
Layer2 : MaxPool(k:(2, 2), stride:(2, 2) padding:(0, 0, 0, 0))

julia> x = rand(Float64, 10, 10, 5, 5) |> size
(10, 10, 5, 5)

julia> N(x) |> size
(4, 4, 5, 2)
```
"""
mutable struct MaxPool{N, M} <: NParam
    k::NTuple{N, Int}
    stride::NTuple{N, Int}
    padding::NTuple{M, Int}
    record::Array
    shape::IOShape
end

function MaxPool(k::NTuple{N, Int}; stride = k, padding = 0) where N
    if N == 1
        k = (k[1], 1)
    end
    stride = convert(2, stride)
    padding = convert(4, padding)
    return MaxPool(k, stride, padding, [], IOShape())
end

function Base.show(io::IO, MP::MaxPool)
    print("MaxPool(k:"*string(MP.k)*", stride:"*string(MP.stride)*", padding:"*string(MP.padding)*")")
end

function (mp::MaxPool)(x::Array{T, N}) where {T, N}
    if N != 4
        x = x[:, :, :, :]
    end
    IC = Im2Col(x, mp.k, mp.stride, mp.padding, trans = "Pool")
    return reshape(maximum(IC.x, dims = 2), IC.Oh, IC.Ow, IC.b, IC.c)
end

function (mp::MaxPool)(x::Array{T, N}, record::Dict, i) where {T, N}
    if N != 4
        x = x[:, :, :, :]
    end
    Ih, Iw = size(x)[1:2]
    IC = Im2Col(x, mp.k, mp.stride, mp.padding, trans = "Pool")
    mp.shape(Ih, Iw, IC.Oh, IC.Ow)
    m = argmax(IC.x, dims = 2)
    r = zero(IC.x)
    for i in m r[i] = 1 end
    mp.record = r
    m = reshape(IC.x[m], IC.Oh, IC.Ow, IC.b, IC.c)
    record[i+1] = (x, m)
    return m
end

function (mp::MaxPool)(Δ, z, back::Bool)
    k = mp.k
    x = reshape(permutedims(Δ, [1, 3, 4, 2]), length(Δ), 1)
    C, B = size(z)[3:4]
    Ih, Iw = mp.shape.IShape
    Oh, Ow = mp.shape.OShape
    Δ = reshape(x .* mp.record, B*Oh*Ow, C*k[1]*k[2])'
    return Col2Im(Δ, k, mp.shape, mp.stride, mp.padding).x
end


text:
"""
    ReductionOperation

A reduction operation on a set of values (e.g. maximum).
"""
abstract type ReductionOperation end

struct MaxReduction  <: ReductionOperation end
struct MeanReduction <: ReductionOperation end

_reduce(op::MaxReduction, x)  = maximum(x)
_reduce(op::MeanReduction, x) = sum(x) / length(x)

"""
    GenericHausdorff(inner_op, outer_op)

The generalized Hausdorff distance with inner reduction `inner_op`
and outer reduction `outer_op`.

## References

Dubuisson, M-P; Jain, A. K., 1994. *A Modified Hausdorff Distance for Object-Matching*.

See also: [`Hausdorff`](@ref), [`ModifiedHausdorff`](@ref)
"""
struct GenericHausdorff{I<:ReductionOperation, O<:ReductionOperation} <: Metric
    inner_op::I
    outer_op::O
end

@doc raw"""
    Hausdorff <: GenericHausdorff
    hausdorff(x::BoolImage, y::BoolImage)

Hausdorff distance between two point sets ``\mathcal{A}`` and ``\mathcal{B}``, which
consist of point positions of 1s of `x` and `y`, respectively. The distance is computed
using the following formula:

```math
    \text{hausdorff}(x,y) = max(d(\mathcal{A}, \mathcal{B}), d(\mathcal{B}, \mathcal{A}))
```
where
```math
    d(\mathcal{A}, \mathcal{B}) = max_{a\in\mathcal{A}}d(a, \mathcal{B}) \\
    d(a, \mathcal{B}) = max_{b\in\mathcal{B}}d(a, b)
```
and the point distance is calcuated using [`Euclidean`](@ref) distance.

## References

Dubuisson, M-P; Jain, A. K., 1994. *A Modified Hausdorff Distance for Object-Matching*.

See also: [`modified_hausdorff`](@ref)
"""
const Hausdorff = GenericHausdorff{MaxReduction, MaxReduction}
Hausdorff() = Hausdorff(MaxReduction(), MaxReduction())

@doc raw"""
    ModifiedHausdorff <: GenericHausdorff
    modified_hausdorff(x::BoolImage, y::BoolImage)

Modified Hausdorff distance between two point sets ``\mathcal{A}`` and ``\mathcal{B}``,
which consist of point positions of 1s of `x` and `y`, respectively. The distance is
computed using the following formula:

```math
    \text{modified_hausdorff}(x,y) = max(d(\mathcal{A}, \mathcal{B}), d(\mathcal{B}, \mathcal{A}))
```
where
```math
    d(\mathcal{A}, \mathcal{B}) = \frac{1}{N_a}\sum_{a\in\mathcal{A}}d(a, \mathcal{B}) \\
    d(a, \mathcal{B}) = max_{b\in\mathcal{B}}d(a, b)
```
and the point distance is calcuated using [`Euclidean`](@ref) distance.

## References

Dubuisson, M-P; Jain, A. K., 1994. *A Modified Hausdorff Distance for Object-Matching*.

See also: [`hausdorff`](@ref)
"""
const ModifiedHausdorff = GenericHausdorff{MeanReduction, MaxReduction}
ModifiedHausdorff() = ModifiedHausdorff(MeanReduction(), MaxReduction())

# convert binary image to a point set format
function img2pset(img::AbstractArray{T}) where T<:Union{Gray{Bool}, Bool}
    inds = findall(x->x==true, img)
    [inds[j][i] for i=1:ndims(img), j=1:length(inds)]
end
function img2pset(img::GenericGrayImage)
    try
        return img2pset(of_eltype(Bool, img))
    catch e
        e isa InexactError && throw(ArgumentError("Binary image is needed."))
        rethrow(e)
    end
end

function evaluate_pset(d::GenericHausdorff, psetA, psetB)
    # trivial cases
    psetA == psetB && return 0.
    (isempty(psetA) || isempty(psetA)) && return Inf

    D = pairwise(Euclidean(), psetA, psetB, dims=2)

    dAB = _reduce(d.inner_op, minimum(D, dims=2))
    dBA = _reduce(d.inner_op, minimum(D, dims=1))

    _reduce(d.outer_op, (dAB, dBA))
end

(d::GenericHausdorff)(imgA::GenericGrayImage, imgB::GenericGrayImage) =
    evaluate_pset(d, img2pset(imgA), img2pset(imgB))

# helper functions
@doc (@doc Hausdorff)
hausdorff(imgA::GenericGrayImage, imgB::GenericGrayImage) =
    Hausdorff()(imgA, imgB)

@doc (@doc ModifiedHausdorff)
modified_hausdorff(imgA::GenericGrayImage, imgB::GenericGrayImage)  =
    ModifiedHausdorff()(imgA, imgB)

# precalculate psets to accelerate computing
function pairwise(d::GenericHausdorff,
                  imgsA::AbstractVector{<:GenericGrayImage},
                  imgsB::AbstractVector{<:GenericGrayImage})
    psetsA = [img2pset(imgA) for imgA in imgsA]
    psetsB = [img2pset(imgB) for imgB in imgsB]

    m, n = length(imgsA), length(imgsB)
    D = zeros(m, n)

    for j=1:n
      psetB = psetsB[j]
      for i=min(m, j+1):m
        psetA = psetsA[i]
        D[i,j] = evaluate_pset(d, psetA, psetB)
      end
      for i=1:min(m, j+1)
        psetA = psetsA[i]
        D[i,j] = evaluate_pset(d, psetA, psetB)
      end
    end

    D
end

function pairwise(d::GenericHausdorff, imgs::AbstractVector{<:GenericGrayImage})
    psets = [img2pset(img) for img in imgs]

    n = length(imgs)
    D = zeros(n, n)
    for j=1:n
      psetB = psets[j]
      for i=j+1:n
        psetA = psets[i]
        D[i,j] = evaluate_pset(d, psetA, psetB)
      end
      # nothing to be done to the diagonal (always zero)
      for i=1:j-1
        D[i,j] = D[j,i] # leverage the symmetry
      end
    end

    D
end


text:
module WindowAbstractions

using DocStringExtensions

@template (FUNCTIONS, METHODS, MACROS) =
    """
    $(DOCSTRING)
    $(TYPEDSIGNATURES)
    """

@template TYPES =
    """
    $(DOCSTRING)
    $(TYPEDEF)
    $(TYPEDSIGNATURES)
    $(TYPEDFIELDS)
    """

include("utils.jl")
include("types.jl")
include("keys.jl")
include("exceptions.jl")
include("inputs.jl")
include("window.jl")
include("window_manager.jl")
include("events.jl")

export
        # Windows
        AbstractWindow,
        AbstractWindowManager,
        WindowException,
        InvalidWindow,
        CloseWindow,

        # Event types
        Event,
        EventData,
        KeyEvent,
        KeyAction,
        KeyPressed,
        KeyReleased,
        MouseEvent,
        MouseAction,
        ButtonPressed,
        ButtonReleased,
        PointerEvent,
        PointerMoves,
        PointerEntersWindow,
        PointerLeavesWindow,
        Resize,
        ResizeEvent,
        Expose,
        ExposeEvent,
        PointerEntersWindowEvent,
        PointerLeavesWindowEvent,
        PointerMovesEvent,

        # Keyboard
        KeyModifierState,
        KeyContext,
        KeyCombination,
        no_input,
        KeySymbol,
        key_symbol,
        is_key_defined,

        # Mouse
        MouseButton,
        ButtonLeft,
        ButtonMiddle,
        ButtonRight,
        ButtonScrollDown,
        ButtonScrollUp,
        MouseState,

        # Windowing API functions
        extent,
        dpi,
        set_extent,
        terminate_window!,
        map_window,
        unmap_window,
        set_title,
        set_icon,
        set_icon_title,
        attach_graphics_context!,
        wait_for_event,
        poll_for_event,
        get_window,
        get_window_symbol,
        pressed_buttons,

        # Events and callbacks
        EventDetails,
        action,
        Callbacks,
        WindowCallbacks,
        callback_symbol,
        execute_callback,

        # Execution modes
        ExecutionMode,
        Synchronous,
        Asynchronous,

        # Macros
        @key_str

end


text:
# `NT` for number of tensors
abstract type EinRule{NT} end

struct Sum <: EinRule{1} end
struct Tr <: EinRule{1} end
struct PairWise <: EinRule{Any} end
struct Permutedims <: EinRule{1} end
struct Hadamard <: EinRule{Any} end
struct PTrace <: EinRule{1} end
struct MatMul <: EinRule{2} end
struct Identity <: EinRule{Any} end
struct BatchedContract <: EinRule{2} end
struct DefaultRule <: EinRule{Any} end

"""
a einsum code is trace
"""
function match_rule(::Type{Tr}, ixs, iy)
    iy == () &&
    length.(ixs) === (2,) &&
    ixs[1][1] == ixs[1][2]
end

"""
a einsum code is a pairwise graph.
"""
function match_rule(::Type{PairWise}, ixs::NTuple{N, NTuple{X,T} where X}, iy::NTuple{M, T}) where {N, M, T}
    allinds = TupleTools.vcat(ixs..., iy)
    all(i -> count(==(i), allinds) == 2, allinds) && allunique(iy)
end

"""
a einsum code is sum.
"""
function match_rule(::Type{Sum}, ixs, iy)
    length(ixs) != 1 && return false
    (ix,) = ixs
    allunique(ix) && allunique(iy) && all(i -> i in ix, iy)
end

"""
permutation rule
"""
function match_rule(::Type{Permutedims}, ixs, iy)
    length(ixs) == 1 || return false
    (ix,) = ixs
    length(ix) == length(iy) && allunique(ix) && allunique(iy) &&
     all(i -> i in iy, ix)
end

"""
Hadamard
"""
function match_rule(::Type{Hadamard}, ixs, iy)
    allunique(iy) && all(allunique, ixs) &&  all(==(iy), ixs)
end

"""
Ptrace rule if all indices of one ix in ixs all appear in iy or
appear twice and don't appear in iy
"""
function match_rule(::Type{PTrace}, ixs, iy)
    length(ixs) == 1 || return false
    (ix,) = ixs
    all(iy) do i
        count(==(i), ix) == 1 && count(==(i), iy) == 1
    end || return false
    all(ix) do i
        ciy, cix = count.(==(i), (iy,ix))
        (ciy == 0 && cix == 2) ||
        (ciy == 1 && cix == 1) ||
        (ciy > 1)
    end || return false
end

function match_rule(::Type{MatMul}, ixs, iy)
    length.(ixs) == (2,2) && length(iy) == 2 || return false
    ix1, ix2 = ixs
    length(intersect(ix1,ix2)) == 1 && Set(symdiff(ix1,ix2)) == Set(iy)
end

function match_rule(::Type{Identity}, ixs, iy)
    ixs === (iy,) && allunique(iy)
end

function match_rule(::Type{BatchedContract}, ixs::NTuple{X, NTuple{N,T} where N}, iy::NTuple{M,T}) where {T,X,M}
    X == 2 || return false
    iA, iB = ixs
    all(allunique, (iA, iB, iy)) || return false
    allinds = (iA..., iB..., iy...)
    all(x -> count(==(x), allinds) > 1, allinds)
end

const einsum_rules = [
    Identity,
    MatMul,
    Permutedims,
    Hadamard,
    Tr,
    PTrace,
    Sum,
    BatchedContract,
    PairWise,
    ]

@doc raw"
    match_rule(ixs, iy)
    match_rule(code::EinCode{ixs, iy})
    match_rule(code::NestedEinCode)

go through all operations specified in the `einsum_rules`-vector and return
the first `T` for which `match_rule(T, ixs, iy)` returns true.
"
function match_rule(ixs, iy)
    # the first rule with the higher the priority
    for T in einsum_rules
        match_rule(T, ixs, iy) && return T()
    end
    return DefaultRule()
end

match_rule(code::EinCode{ixs, iy}) where {ixs, iy} = match_rule(ixs, iy)