text:
function DiffEqBase.__solve(prob::DiffEqBase.AbstractDDEProblem,
                            alg::AbstractMethodOfStepsAlgorithm, args...;
                            kwargs...)
  integrator = DiffEqBase.__init(prob, alg, args...; kwargs...)
  DiffEqBase.solve!(integrator)
  integrator.sol
end

function DiffEqBase.__init(prob::DiffEqBase.AbstractDDEProblem,
                           alg::AbstractMethodOfStepsAlgorithm,
                           timeseries_init = (),
                           ts_init = (),
                           ks_init = ();
                           saveat = (),
                           tstops = (),
                           d_discontinuities = (),
                           save_idxs = nothing,
                           save_everystep = isempty(saveat),
                           save_on = true,
                           save_start = save_everystep || isempty(saveat) || saveat isa Number || prob.tspan[1] in saveat,
                           save_end = nothing,
                           callback = nothing,
                           dense = save_everystep && isempty(saveat),
                           calck = (callback !== nothing && callback != CallbackSet()) || # Empty callback
                           dense, # and no dense output
                           dt = zero(eltype(prob.tspan)),
                           dtmin = DiffEqBase.prob2dtmin(prob; use_end_time=false),
                           dtmax = eltype(prob.tspan)(prob.tspan[end]-prob.tspan[1]),
                           force_dtmin = false,
                           adaptive = DiffEqBase.isadaptive(alg),
                           gamma = OrdinaryDiffEq.gamma_default(alg.alg),
                           abstol = nothing,
                           reltol = nothing,
                           qmin = OrdinaryDiffEq.qmin_default(alg.alg),
                           qmax = OrdinaryDiffEq.qmax_default(alg.alg),
                           qsteady_min = OrdinaryDiffEq.qsteady_min_default(alg.alg),
                           qsteady_max = OrdinaryDiffEq.qsteady_max_default(alg.alg),
                           qoldinit = DiffEqBase.isadaptive(alg) ? 1//10^4 : 0,
                           controller = nothing,
                           fullnormalize = true,
                           failfactor = 2,
                           beta1 = nothing,
                           beta2 = nothing,
                           maxiters = adaptive ? 1000000 : typemax(Int),
                           internalnorm = DiffEqBase.ODE_DEFAULT_NORM,
                           internalopnorm = LinearAlgebra.opnorm,
                           isoutofdomain = DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN,
                           unstable_check = DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK,
                           verbose = true,
                           timeseries_errors = true,
                           dense_errors = false,
                           advance_to_tstop = false,
                           stop_at_next_tstop = false,
                           initialize_save = true,
                           progress = false,
                           progress_steps = 1000,
                           progress_name = "DDE",
                           progress_message = DiffEqBase.ODE_DEFAULT_PROG_MESSAGE,
                           userdata = nothing,
                           allow_extrapolation = OrdinaryDiffEq.alg_extrapolates(alg),
                           initialize_integrator = true,
                           alias_u0 = false,
                           # keyword arguments for DDEs
                           discontinuity_interp_points::Int = 10,
                           discontinuity_abstol = eltype(prob.tspan)(1//Int64(10)^12),
                           discontinuity_reltol = 0,
                           kwargs... )
  if haskey(kwargs, :initial_order)
    @warn "initial_order has been deprecated. Please specify order_discontinuity_t0 in the DDEProblem instead."
    order_discontinuity_t0::Int = kwargs[:initial_order]
  else
    order_discontinuity_t0 = prob.order_discontinuity_t0
  end

  if alg.alg isa CompositeAlgorithm && alg.alg.choice_function isa AutoSwitch
    auto = alg.alg.choice_function
    alg = MethodOfSteps(CompositeAlgorithm(alg.alg.algs,
                             OrdinaryDiffEq.AutoSwitchCache(
                                             0,0,
                                             auto.nonstiffalg,
                                             auto.stiffalg,
                                             auto.stiffalgfirst,
                                             auto.maxstiffstep,
                                             auto.maxnonstiffstep,
                                             auto.nonstifftol,
                                             auto.stifftol,
                                             auto.dtfac,
                                             auto.stiffalgfirst,
                                             auto.switch_max
                                            )))
  end

  if haskey(kwargs, :minimal_solution)
    @warn "minimal_solution is ignored"
  end

  if !isempty(saveat) && dense
    @warn("Dense output is incompatible with saveat. Please use the SavingCallback from the Callback Library to mix the two behaviors.")
  end

  progress && @logmsg(-1, progress_name, _id = _id = :DelayDiffEq, progress = 0)

  isdae = prob.f.mass_matrix !== I && !(prob.f.mass_matrix isa Tuple) &&
      ArrayInterface.issingular(prob.f.mass_matrix)

  # unpack problem
  @unpack f, u0, h, tspan, p, neutral, constant_lags, dependent_lags = prob

  # determine type and direction of time
  tType = eltype(tspan)
  t0 = first(tspan)
  tdir = sign(last(tspan) - t0)
  tTypeNoUnits = typeof(one(tType))

  # Allow positive dtmax, but auto-convert
  dtmax > zero(dtmax) && tdir < zero(tdir) && (dtmax *= tdir)

  # no fixed-point iterations for constrained algorithms,
  # and thus `dtmax` should match minimal lag
  if isconstrained(alg) && has_constant_lags(prob)
    dtmax = tdir * min(abs(dtmax), minimum(abs, constant_lags))
  end

  # get absolute and relative tolerances
  abstol_internal = get_abstol(u0, tspan, alg.alg; abstol = abstol)
  reltol_internal = get_reltol(u0, tspan, alg.alg; reltol = reltol)

  # get rate prototype
  rate_prototype = rate_prototype_of(u0, tspan)

  # create a history function
  history = build_history_function(prob, alg, rate_prototype, reltol_internal;
                                   dt = dt, dtmin = dtmin, calck = false,
                                   adaptive = adaptive, internalnorm = internalnorm)
  f_with_history = ODEFunctionWrapper(f, history)

  # get states (possibly different from the ODE integrator!)
  u, uprev, uprev2 = u_uprev_uprev2(u0, alg;
                                    alias_u0 = alias_u0,
                                    adaptive = adaptive,
                                    allow_extrapolation = allow_extrapolation,
                                    calck = calck)
  uEltypeNoUnits = recursive_unitless_eltype(u)
  uBottomEltypeNoUnits = recursive_unitless_bottom_eltype(u)

  # initialize output arrays of the solution
  k = typeof(rate_prototype)[]
  ts, timeseries, ks = solution_arrays(u, tspan, rate_prototype;
                                       timeseries_init = timeseries_init,
                                       ts_init = ts_init,
                                       ks_init = ks_init,
                                       save_idxs = save_idxs,
                                       save_start = save_start)

  # build cache
  ode_integrator = history.integrator
  cache = OrdinaryDiffEq.alg_cache(alg.alg, u, rate_prototype, uEltypeNoUnits,
                                   uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2,
                                   f_with_history, t0, zero(tType), reltol_internal, p, calck,
                                   Val(isinplace(prob)))

  # separate statistics of the integrator and the history
  destats = DDEStats(0)

  # create solution
  if iscomposite(alg)
    id = OrdinaryDiffEq.CompositeInterpolationData(f_with_history, timeseries, ts, ks,
                                                   Int[], dense, cache)
    sol = DiffEqBase.build_solution(prob, alg.alg, ts, timeseries;
                                    dense = dense, k = ks, interp = id,
                                    alg_choice = id.alg_choice, calculate_error = false,
                                    destats = destats)
  else
    id = OrdinaryDiffEq.InterpolationData(f_with_history, timeseries, ts, ks, dense, cache)
    sol = DiffEqBase.build_solution(prob, alg.alg, ts, timeseries;
                                    dense = dense, k = ks, interp = id,
                                    calculate_error = false, destats = destats)
  end

  # retrieve time stops, time points at which solutions is saved, and discontinuities
  tstops_internal =
    OrdinaryDiffEq.initialize_tstops(tType, tstops, d_discontinuities, tspan)
  saveat_internal = OrdinaryDiffEq.initialize_saveat(tType, saveat, tspan)
  d_discontinuities_internal = OrdinaryDiffEq.initialize_d_discontinuities(
    Discontinuity{tType,Int}, d_discontinuities, tspan
  )

  maximum_order = OrdinaryDiffEq.alg_maximum_order(alg)
  tstops_propagated, d_discontinuities_propagated =
    initialize_tstops_d_discontinuities_propagated(tType, tstops, d_discontinuities,
                                                   tspan, order_discontinuity_t0,
                                                   maximum_order, constant_lags, neutral)

  # reserve capacity for the solution
  sizehint!(sol, alg, tspan, tstops_internal, saveat_internal;
            save_everystep = save_everystep, adaptive = adaptive, dt = tType(dt),
            dtmin = dtmin, internalnorm = internalnorm)

  # create array of tracked discontinuities
  # used to find propagated discontinuities with callbacks and to keep track of all
  # passed discontinuities
  tracked_discontinuities = Discontinuity{tType,Int}[]
  if order_discontinuity_t0 ≤ maximum_order
    push!(tracked_discontinuities, Discontinuity(tdir * t0, order_discontinuity_t0))
  end

  # Create set of callbacks and its cache
  callback_set, callback_cache = callback_set_and_cache(prob, callback)

  # separate options of integrator and of dummy ODE integrator since ODE integrator always saves
  # every step and every index (necessary for history function)
  QT = tTypeNoUnits <: Integer ? typeof(qmin) : typeof(internalnorm(u, t0))

  # Setting up the step size controller
  if (beta1 !== nothing || beta2 !== nothing) && controller !== nothing
    throw(ArgumentError(
      "Setting both the legacy PID parameters `beta1, beta2 = $((beta1, beta2))` and the `controller = $controller` is not allowed."))
  end

  if (beta1 !== nothing || beta2 !== nothing)
    message = "Providing the legacy PID parameters `beta1, beta2` is deprecated. Use the keyword argument `controller` instead."
    Base.depwarn(message, :init)
    Base.depwarn(message, :solve)
  end

  if controller === nothing
    controller = OrdinaryDiffEq.default_controller(alg.alg, cache, convert(QT, qoldinit)::QT, beta1, beta2)
  end

  save_end_user = save_end
  save_end = save_end === nothing ? save_everystep || isempty(saveat) || saveat isa Number || prob.tspan[2] in saveat : save_end

  opts = OrdinaryDiffEq.DEOptions{typeof(abstol_internal),typeof(reltol_internal),
                                  QT,tType,typeof(controller),
                                  typeof(internalnorm),typeof(internalopnorm),
                                  typeof(save_end_user),
                                  typeof(callback_set),
                                  typeof(isoutofdomain),
                                  typeof(progress_message),typeof(unstable_check),
                                  typeof(tstops_internal),
                                  typeof(d_discontinuities_internal),typeof(userdata),
                                  typeof(save_idxs),
                                  typeof(maxiters),typeof(tstops),
                                  typeof(saveat),typeof(d_discontinuities)}(
                                    maxiters,save_everystep,
                                    adaptive,abstol_internal,reltol_internal,
                                    QT(gamma),QT(qmax),QT(qmin),
                                    QT(qsteady_max),QT(qsteady_min),
                                    QT(qoldinit),QT(failfactor),
                                    tType(dtmax),tType(dtmin),
                                    controller,
                                    internalnorm,internalopnorm,
                                    save_idxs,tstops_internal,saveat_internal,
                                    d_discontinuities_internal,
                                    tstops,saveat,d_discontinuities,
                                    userdata,progress,progress_steps,
                                    progress_name,progress_message,
                                    timeseries_errors,dense_errors,dense,
                                    save_on,save_start,save_end,save_end_user,
                                    callback_set,
                                    isoutofdomain,unstable_check,
                                    verbose,calck,force_dtmin,
                                    advance_to_tstop,stop_at_next_tstop)

  # create fixed point solver
  fpsolver = build_fpsolver(alg, alg.fpsolve, u, uEltypeNoUnits, uBottomEltypeNoUnits,
                            Val(isinplace(prob)))

  # initialize indices of u(t) and u(tprev) in the dense history
  prev_idx = 1
  prev2_idx = 1

  # create integrator combining the new defined problem function with history
  # information, the new solution, the parameters of the ODE integrator, and
  # parameters of fixed-point iteration
  # do not initialize fsalfirst and fsallast
  # rate/state = (state/time)/state = 1/t units, internalnorm drops units
  eigen_est = inv(one(tType))
  tprev = t0
  dtcache = tType(dt)
  dtpropose = tType(dt)
  iter = 0
  kshortsize = 0
  reeval_fsal = false
  u_modified = false
  EEst = QT(1)
  just_hit_tstop = false
  do_error_check = true
  isout = false
  accept_step = false
  force_stepfail = false
  last_stepfail = false
  event_last_time = 0
  vector_event_last_time = 1
  last_event_error = zero(uBottomEltypeNoUnits)
  dtchangeable = OrdinaryDiffEq.isdtchangeable(alg.alg)
  q11 = QT(1)
  success_iter = 0
  erracc = QT(1)
  dtacc = tType(1)

  integrator = DDEIntegrator{typeof(alg.alg),isinplace(prob),typeof(u0),tType,typeof(p),
                             typeof(eigen_est),QT,typeof(tdir),typeof(k),typeof(sol),
                             typeof(f_with_history),typeof(cache),
                             typeof(ode_integrator),typeof(fpsolver),
                             typeof(opts),typeof(discontinuity_abstol),
                             typeof(discontinuity_reltol),typeof(history),
                             typeof(tstops_propagated),typeof(d_discontinuities_propagated),
                             OrdinaryDiffEq.fsal_typeof(alg.alg, rate_prototype),
                             typeof(last_event_error),typeof(callback_cache)}(
                               sol, u, k, t0, tType(dt), f_with_history, p,
                               uprev, uprev2, tprev, prev_idx, prev2_idx,
                               fpsolver,
                               order_discontinuity_t0, tracked_discontinuities,
                               discontinuity_interp_points, discontinuity_abstol,
                               discontinuity_reltol, tstops_propagated,
                               d_discontinuities_propagated, alg.alg,
                               dtcache, dtchangeable, dtpropose, tdir,
                               eigen_est, EEst, QT(qoldinit), q11,
                               erracc, dtacc, success_iter,
                               iter, length(ts), length(ts), cache, callback_cache,
                               kshortsize, force_stepfail, last_stepfail,
                               just_hit_tstop, do_error_check, event_last_time, vector_event_last_time,
                               last_event_error, accept_step,
                               isout, reeval_fsal,
                               u_modified, isdae, opts, destats, history,
                               ode_integrator)

  # initialize DDE integrator
  if initialize_integrator
    initialize_solution!(integrator)
    OrdinaryDiffEq.initialize_callbacks!(integrator, initialize_save)
    OrdinaryDiffEq.initialize!(integrator)
  end

  # take care of time step dt = 0 and dt with incorrect sign
  OrdinaryDiffEq.handle_dt!(integrator)

  integrator
end

function DiffEqBase.solve!(integrator::DDEIntegrator)
  @unpack tdir, opts, sol = integrator
  @unpack tstops = opts

  # step over all stopping time points, similar to solving with ODE integrators
  @inbounds while !isempty(tstops)
    while tdir * integrator.t < first(tstops)
      # apply step or adapt step size
      OrdinaryDiffEq.loopheader!(integrator)

      # abort integration following same criteria as for ODEs:
      # maxiters exceeded, dt <= dtmin, integration unstable
      DiffEqBase.check_error!(integrator) === :Success || return sol

      # calculate next step
      OrdinaryDiffEq.perform_step!(integrator)

      # calculate proposed next step size, handle callbacks, and update solution
      OrdinaryDiffEq.loopfooter!(integrator)

      isempty(tstops) && break
    end

    # remove hit or passed stopping time points
    OrdinaryDiffEq.handle_tstop!(integrator)
  end

  # clean up solution
  DiffEqBase.postamble!(integrator)

  f = sol.prob.f

  if DiffEqBase.has_analytic(f)
    DiffEqBase.calculate_solution_errors!(sol;
                                          timeseries_errors = opts.timeseries_errors,
                                          dense_errors = opts.dense_errors)
  end
  sol.retcode === :Default || return sol

  integrator.sol = DiffEqBase.solution_new_retcode(sol, :Success)
end

function OrdinaryDiffEq.initialize_callbacks!(integrator::DDEIntegrator,
                                              initialize_save = true)
  callbacks = integrator.opts.callback
  prob = integrator.sol.prob

  # set up additional initial values of newly created DDE integrator
  # (such as fsalfirst) and its callbacks

  integrator.u_modified = true

  u_modified = initialize!(callbacks, integrator.u, integrator.t, integrator)

  # if the user modifies u, we need to fix previous values before initializing
  # FSAL in order for the starting derivatives to be correct
  if u_modified

    if isinplace(prob)
      recursivecopy!(integrator.uprev, integrator.u)
    else
      integrator.uprev = integrator.u
    end

    if OrdinaryDiffEq.alg_extrapolates(integrator.alg)
      if isinplace(prob)
        recursivecopy!(integrator.uprev2, integrator.uprev)
      else
        integrator.uprev2 = integrator.uprev
      end
    end

    # update heap of discontinuities
    # discontinuity is assumed to be of order 0, i.e. solution x is discontinuous
    add_next_discontinuities!(integrator, 0, integrator.t)

    # reset this as it is now handled so the integrators should proceed as normal
    reeval_internals_due_to_modification!(integrator, Val{false})

    if initialize_save &&
      (any((c)->c.save_positions[2],callbacks.discrete_callbacks) ||
       any((c)->c.save_positions[2],callbacks.continuous_callbacks))
      savevalues!(integrator, true)
    end
  end

  # reset this as it is now handled so the integrators should proceed as normal
  integrator.u_modified = false
end

function initialize_tstops_d_discontinuities_propagated(
  ::Type{T}, tstops, d_discontinuities, tspan, order_discontinuity_t0, alg_maximum_order,
  constant_lags, neutral,
) where T
  # create heaps for propagated discontinuities and corresponding time stops
  tstops_propagated = BinaryMinHeap{T}()
  d_discontinuities_propagated = BinaryMinHeap{Discontinuity{T,Int}}()

  # add discontinuities and time stops propagated from initial discontinuity
  if constant_lags !== nothing && !isempty(constant_lags) &&
    order_discontinuity_t0 ≤ alg_maximum_order

    sizehint!(tstops_propagated, length(constant_lags))
    sizehint!(d_discontinuities_propagated, length(constant_lags))

    t0, tf = tspan
    tdir = sign(tf - t0)
    maxlag = tdir * (tf - t0)
    next_order = neutral ? order_discontinuity_t0 : order_discontinuity_t0 + 1
    for lag in constant_lags
      if tdir * lag < maxlag
        t = tdir * (t0 + lag)
        push!(tstops_propagated, t)

        d = Discontinuity{T,Int}(t, next_order)
        push!(d_discontinuities_propagated, d)
      end
    end
  end

  return tstops_propagated, d_discontinuities_propagated
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
#Look up Table
include("lut/mc.jl")

function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {T, VertType, FaceType}
    nx, ny, nz = size(sdf)

    # we subtract one from the length along each axis because
    # an NxNxN SDF has N-1 cells on each axis
    s = VertType(widths[1]/(nx-1), widths[2]/(ny-1), widths[3]/(nz-1))

    # arrays for vertices and faces
    vts = VertType[]
    fcs = FaceType[]
    mt = max(nx,ny,nz)
    method.reduceverts && sizehint!(vts, mt*mt*5)
    !method.reduceverts && sizehint!(vts, mt*mt*6)
    sizehint!(fcs, mt*mt*2)
    @inbounds for zi = 1:nz-1, yi = 1:ny-1, xi = 1:nx-1

        iso_vals = (sdf[xi,yi,zi],
                    sdf[xi+1,yi,zi],
                    sdf[xi+1,yi+1,zi],
                    sdf[xi,yi+1,zi],
                    sdf[xi,yi,zi+1],
                    sdf[xi+1,yi,zi+1],
                    sdf[xi+1,yi+1,zi+1],
                    sdf[xi,yi+1,zi+1])

        #Determine the index into the edge table which
        #tells us which vertices are inside of the surface
        cubeindex = method.insidepositive ? _get_cubeindex_pos(iso_vals, method.iso) : _get_cubeindex(iso_vals, method.iso)

        # Cube is entirely in/out of the surface
        _no_triangles(cubeindex) && continue

        points = mc_vert_points(xi,yi,zi,s,origin,VertType)

        # Create the triangle
        method.reduceverts && _mc_unique_triangles!(method, points, iso_vals, vts, fcs, cubeindex, FaceType)
        !method.reduceverts && _mc_create_triangles!(method, points, iso_vals, vts, fcs, cubeindex, FaceType)
    end
    vts,fcs
end

function isosurface(f::Function, method::MarchingCubes, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2), samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T <: Integer, VertType, FaceType}

    nx, ny, nz = samples[1], samples[2], samples[3]

    # we subtract one from the length along each axis because
    # an NxNxN SDF has N-1 cells on each axis
    s = VertType(widths[1]/(nx-1), widths[2]/(ny-1), widths[3]/(nz-1))

    # arrays for vertices and faces
    vts = VertType[]
    fcs = FaceType[]
    mt = max(nx,ny,nz)
    method.reduceverts && sizehint!(vts, mt*mt*5)
    !method.reduceverts && sizehint!(vts, mt*mt*6)
    sizehint!(fcs, mt*mt*2)
    zv = zero(eltype(VertType))
    iso_vals = (zv,zv,zv,zv,zv,zv,zv,zv)
    @inbounds for zi = 1:nz-1, yi = 1:ny-1, xi = 1:nx-1

        points = mc_vert_points(xi,yi,zi,s,origin,VertType)

        if xi == 1
            iso_vals = (f(points[1]),
                        f(points[2]),
                        f(points[3]),
                        f(points[4]),
                        f(points[5]),
                        f(points[6]),
                        f(points[7]),
                        f(points[8]))
        else
            iso_vals = (iso_vals[2],
                        f(points[2]),
                        f(points[3]),
                        iso_vals[3],
                        iso_vals[6],
                        f(points[6]),
                        f(points[7]),
                        iso_vals[7])
        end

        #Determine the index into the edge table which
        #tells us which vertices are inside of the surface
        cubeindex = method.insidepositive ? _get_cubeindex_pos(iso_vals, method.iso) : _get_cubeindex(iso_vals, method.iso)

        # Cube is entirely in/out of the surface
        _no_triangles(cubeindex) && continue

        # Create the triangle
        method.reduceverts && _mc_unique_triangles!(method, points, iso_vals, vts, fcs, cubeindex, FaceType)
        !method.reduceverts && _mc_create_triangles!(method, points, iso_vals, vts, fcs, cubeindex, FaceType)

    end
    vts,fcs
end

"""
    _mc_create_triangles!(vts, fcs, vertlist, cubeindex, FaceType)

Create triangles by adding every point within a triangle to the vertex vector.
"""
@inline function _mc_create_triangles!(method, points, iso_vals, vts, fcs, cubeindex, ::Type{FaceType}) where {FaceType}

    vertlist = find_vertices_interp(points, iso_vals, cubeindex, method.iso, method.eps)

    fct = length(vts) + 3

    push!(vts, vertlist[tri_table[cubeindex][1]],
               vertlist[tri_table[cubeindex][2]],
               vertlist[tri_table[cubeindex][3]])
    push!(fcs, FaceType(fct, fct-1, fct-2))

    iszero(tri_table[cubeindex][4]) && return
    fct += 3
    push!(vts, vertlist[tri_table[cubeindex][4]],
               vertlist[tri_table[cubeindex][5]],
               vertlist[tri_table[cubeindex][6]])
    push!(fcs, FaceType(fct, fct-1, fct-2))

    iszero(tri_table[cubeindex][7]) && return
    fct += 3
    push!(vts, vertlist[tri_table[cubeindex][7]],
               vertlist[tri_table[cubeindex][8]],
               vertlist[tri_table[cubeindex][9]])
    push!(fcs, FaceType(fct, fct-1, fct-2))

    iszero(tri_table[cubeindex][10]) && return
    fct += 3
    push!(vts, vertlist[tri_table[cubeindex][10]],
               vertlist[tri_table[cubeindex][11]],
               vertlist[tri_table[cubeindex][12]])
    push!(fcs, FaceType(fct, fct-1, fct-2))

    iszero(tri_table[cubeindex][13]) && return
    fct += 3
    push!(vts, vertlist[tri_table[cubeindex][13]],
               vertlist[tri_table[cubeindex][14]],
               vertlist[tri_table[cubeindex][15]])
    push!(fcs, FaceType(fct, fct-1, fct-2))
end

"""
    _mc_unique_triangles!(vts, fcs, vertlist, cubeindex, FaceType)

Create triangles by only adding unique vertices within the voxel.
Each face may share a reference to a vertex with another face.
"""
function _mc_unique_triangles!(method, points, iso_vals, vts, fcs, cubeindex, ::Type{FaceType}) where {FaceType}
    @inbounds begin
        fct = length(vts)

        find_vertices_interp!(vts, points, iso_vals, cubeindex, method.iso, method.eps)

        offsets = _mc_connectivity[_mc_eq_mapping[cubeindex]]

        # There is atleast one face so we can push it immediately
        push!(fcs, FaceType(fct+offsets[3], fct+offsets[2], fct+offsets[1]))

        for i in (4,7,10,13)
            iszero(offsets[i]) && return
            push!(fcs, FaceType(fct+offsets[i+2], fct+offsets[i+1], fct+offsets[i]))
        end
    end
end

"""
    find_vertices_interp(points, iso_vals, cubeindex, iso, eps)

Find the vertices where the surface intersects the cube
"""
@inline function find_vertices_interp(points, iso_vals, cubeindex, iso, eps)
    VT = eltype(points)
    zv = Ref{VT}()[]
    (!iszero(edge_table[cubeindex] & 0x001) ? vertex_interp(iso,points[1],points[2],iso_vals[1],iso_vals[2], eps) : zv,
     !iszero(edge_table[cubeindex] & 0x002) ? vertex_interp(iso,points[2],points[3],iso_vals[2],iso_vals[3], eps) : zv,
     !iszero(edge_table[cubeindex] & 0x004) ? vertex_interp(iso,points[3],points[4],iso_vals[3],iso_vals[4], eps) : zv,
     !iszero(edge_table[cubeindex] & 0x008) ? vertex_interp(iso,points[4],points[1],iso_vals[4],iso_vals[1], eps) : zv,
     !iszero(edge_table[cubeindex] & 0x010) ? vertex_interp(iso,points[5],points[6],iso_vals[5],iso_vals[6], eps) : zv,
     !iszero(edge_table[cubeindex] & 0x020) ? vertex_interp(iso,points[6],points[7],iso_vals[6],iso_vals[7], eps) : zv,
     !iszero(edge_table[cubeindex] & 0x040) ? vertex_interp(iso,points[7],points[8],iso_vals[7],iso_vals[8], eps) : zv,
     !iszero(edge_table[cubeindex] & 0x080) ? vertex_interp(iso,points[8],points[5],iso_vals[8],iso_vals[5], eps) : zv,
     !iszero(edge_table[cubeindex] & 0x100) ? vertex_interp(iso,points[1],points[5],iso_vals[1],iso_vals[5], eps) : zv,
     !iszero(edge_table[cubeindex] & 0x200) ? vertex_interp(iso,points[2],points[6],iso_vals[2],iso_vals[6], eps) : zv,
     !iszero(edge_table[cubeindex] & 0x400) ? vertex_interp(iso,points[3],points[7],iso_vals[3],iso_vals[7], eps) : zv,
     !iszero(edge_table[cubeindex] & 0x800) ? vertex_interp(iso,points[4],points[8],iso_vals[4],iso_vals[8], eps) : zv)
end

"""
    find_vertices_interp(points, iso_vals, cubeindex, iso, eps)

Find the vertices where the surface intersects the cube
"""
function find_vertices_interp!(vts, points, iso_vals, cubeindex, iso, eps)
    @inbounds begin
        vert_to_add = _mc_verts[cubeindex]
        for i = 1:12
            vt = vert_to_add[i]
            iszero(vt) && break
            ed = _mc_edge_list[vt]
            push!(vts, vertex_interp(iso,points[ed[1]],points[ed[2]],iso_vals[ed[1]],iso_vals[ed[2]], eps))
        end
    end
end

"""
    vertex_interp(iso, p1, p2, valp1, valp2, eps = 0.00001)

Linearly interpolate the position where an isosurface cuts
an edge between two vertices, each with their own scalar value
"""
function vertex_interp(iso, p1, p2, valp1, valp2, eps = 0.00001)

    abs(iso - valp1) < eps && return p1
    abs(iso - valp2) < eps && return p2
    abs(valp1-valp2) < eps && return p1
    mu = (iso - valp1) / (valp2 - valp1)
    p = p1 + mu * (p2 - p1)

    return p
end

"""
    mc_vert_points(xi,yi,zi, scale, origin, ::Type{VertType})

Returns a tuple of 8 points corresponding to each corner of a cube
"""
@inline function mc_vert_points(xi,yi,zi, scale, origin, ::Type{VertType}) where VertType
    (VertType(xi-1,yi-1,zi-1) .* scale .+ origin,
     VertType(xi  ,yi-1,zi-1) .* scale .+ origin,
     VertType(xi  ,yi  ,zi-1) .* scale .+ origin,
     VertType(xi-1,yi  ,zi-1) .* scale .+ origin,
     VertType(xi-1,yi-1,zi  ) .* scale .+ origin,
     VertType(xi  ,yi-1,zi  ) .* scale .+ origin,
     VertType(xi  ,yi  ,zi  ) .* scale .+ origin,
     VertType(xi-1,yi  ,zi  ) .* scale .+ origin)
end

ext: jl


text:
module InteractMithril

using WebIO: node, Node, instanceof, props, children, Scope, JSString, @js_str, onimport,
             setobservable!, onjs, WebIO
using Widgets: AbstractWidget, Widget, Widgets
using Observables: on, Observable, AbstractObservable, ObservablePair, Observables
using Dates
using Colors: Colorant, hex

import JSON

export mithril, MithrilComponent

include("component.jl")
include("input.jl")
include("slider.jl")
include("optioninput.jl")

end # module


text:
"""
    LinearModel{M,V}

A continuous-time dynamics model of the form:

```math
\\dot{x} = A x + B u + d
```
"""
struct LinearModel{M,V} <: RD.ContinuousDynamics
    A::M
    B::M
    d::V
    function LinearModel(A::M, B::M, d::V) where {M<:AbstractMatrix, V<:AbstractVector}
        if !(size(A,1) == size(B,1) == size(d,1))
            throw(DimensionMismatch("All inputs must have the same number of rows."))
        elseif size(d,2) != 1
            throw(DimensionMismatch("d must have a single column."))
        end
        new{M,V}(A,B,d)
    end
end

function LinearModel(A::AbstractMatrix, B::AbstractMatrix)
    d = similar(A, size(A,1))
    d .= 0
    LinearModel(A, B, d)
end

function LinearModel(sig::RD.FunctionSignature, diffmethod::RD.DiffMethod,
                     model::RD.ContinuousDynamics, x, u, t=zero(eltype(x)); affine::Bool=true)
    n,m = RD.dims(model)
    J = zeros(n,n+m)
    d = zeros(n)
    z = RD.KnotPoint{n,m}(x, u, t, NaN)
    RD.jacobian!(sig, diffmethod, model, J, d, z)
    if affine
        if sig == RD.InPlace()
            RD.dynamics!(model, d, x, u, t)
        else
            d .= RD.dynamics(model, x, u, t)
        end
    else
        d .= 0
    end
    A = J[:,1:n]
    B = J[:,n+1:n+m]
    LinearModel(A, B, d)
end

# Generates a random controllable linear system
function Base.rand(::Type{LinearModel}, n::Integer, m::Integer)
    A,B = RandomLinearModels.gencontrollable(n, m, :continuous)
    d = randn(n)
    LinearModel(A, B, d)
end

RD.output_dim(model::LinearModel) = length(model.d)
RD.state_dim(model::LinearModel) = size(model.A, 2)
RD.control_dim(model::LinearModel) = size(model.B, 2)
RD.default_signature(model::LinearModel) = RD.InPlace()
RD.default_diffmethod(model::LinearModel) = RD.UserDefined()

getstatecoeffs(model::LinearModel) = model.A
getcontrolcoeffs(model::LinearModel) = model.B
getaffinecoeffs(model::LinearModel) = model.d

function iscontrollable(model::LinearModel)
    A = getstatecoeffs(model)
    B = getcontrolcoeffs(model)
    RandomLinearModels.iscontrollable(A, B)
end

isstable(model::LinearModel) = RandomLinearModels.isstablec(getstatecoeffs(model))

function dynamics(model::LinearModel, x, u)
    model.A * x .+ model.B * u .+ model.d
end

function dynamics!(model::LinearModel, xn, x, u)
    xn .= model.d
    mul!(xn, model.A, x, 1.0, 1.0)
    mul!(xn, model.B, u, 1.0, 1.0)
    nothing
end

function jacobian!(model::LinearModel, J, xn, x, u)
    n,m = RD.dims(model)
    J[:,1:n] .= model.A
    J[:,n+1:n+m] .= model.B
    nothing
end

"""
    DiscreteLinearModel{M,V}

A discrete-time dynamics model of the form:

```math
\\ A x_1 + B u_1 + C x_2 + d = 0
```
"""
struct DiscreteLinearModel{M,V} <: RD.DiscreteDynamics
    A::M
    B::M
    C::M
    d::V
    dt::Float64
    CA::M
    CB::M
    Cd::V
    function DiscreteLinearModel(dt, A::M, B::M, C::M, d::V) where {M<:AbstractMatrix, V<:AbstractVector}
        if !(size(A,1) == size(B,1) == size(C,1) == size(d,1))
            throw(DimensionMismatch("All inputs must have the same number of rows."))
        elseif size(A,2) != size(C,2)
            throw(DimensionMismatch("A and C must have the same number of columns."))
        elseif size(d,2) != 1
            throw(DimensionMismatch("d must have a single column."))
        end
        if C isa SMatrix
            F = lu(C)
        else
            F = factorize(C)
        end
        CA = -(F\A)
        CB = -(F\B)
        Cd = -(F\d)
        new{M,V}(A,B,C,d,dt, CA,CB,Cd)
    end
end

function DiscreteLinearModel(dt, A::AbstractMatrix, B::AbstractMatrix,
                             d::AbstractVector=(similar(A, size(A,1)) .= zero(eltype(A))))
    C = similar(A) * zero(eltype(A))
    for i = 1:size(C,1)
        C[i,i] = -one(eltype(A))
    end
    DiscreteLinearModel(dt, A, B, C, d)
end

function DiscreteLinearModel(sig::RD.FunctionSignature, diffmethod::RD.DiffMethod,
                             model::RD.DiscretizedDynamics,
                             x, u, t, h; affine::Bool=true)
    n,m = RD.dims(model)
    J = zeros(n,n+m)
    d = zeros(n)
    z = RD.KnotPoint{n,m}([x;u], t, h)
    RD.jacobian!(sig, diffmethod, model, J, d, z)
    if affine
        if sig == RD.InPlace()
            RD.discrete_dynamics!(model, d, x, u, t, h)
        else
            d .= RD.discrete_dynamics(model, x, u, t, h)
        end
    else
        d .= 0
    end
    A = J[:,1:n]
    B = J[:,n+1:n+m]
    DiscreteLinearModel(h, A, B, d)
end

function DiscreteLinearModel(sig::RD.FunctionSignature, diffmethod::RD.DiffMethod,
                             model::RD.ImplicitDynamicsModel,
                             x, u, t, h; affine::Bool=true)
    n,m = RD.dims(model)
    J1 = zeros(n,n+m)
    J2 = zeros(n,n+m)
    d = zeros(n)
    y1 = similar(d)
    z1 = RD.StaticKnotPoint{Any,Any}(n,m, [x;u], t, h)
    z2 = RD.StaticKnotPoint{Any,Any}(n,m, [x;u], t, h)
    RD.dynamics_error_jacobian!(sig, diffmethod, model, J2, J1, d, y1, z2, z1)
    if affine
        if sig == RD.InPlace()
            RD.dynamics_error!(model, d, y1, z2, z1)
        else
            d .= RD.dynamics_error(model, d, z2, z1)
        end
    else
        d .= 0
    end
    A = J1[:,1:n]
    B = J1[:,n+1:n+m]
    C = J2[:,1:n]
    if norm(J2[:,n+1:n+m]) > 0
        error("Cannot form a DiscreteLinearModel from an implicit dynamics function dependent on u2.")
    end
    DiscreteLinearModel(h, A, B, C, d)
end

# Exponential integration
function DiscreteLinearModel(model::LinearModel, h)
    if norm(getaffinecoeffs(model)) > 0
        throw(ArgumentError("Cannot use exponential integrator for an integrator with an affine term."))
    end
    n,m = RD.dims(model)
    A = getstatecoeffs(model)
    B = getcontrolcoeffs(model)
    expAB = exp([A*h B*h; zeros(m,n+m)])
    A = expAB[1:n,1:n]
    B = expAB[1:n,n+1:n+m]
    DiscreteLinearModel(h, A, B)
end

ouput_dim(model::DiscreteLinearModel) = size(model.A, 1)
state_dim(model::DiscreteLinearModel) = size(model.A, 2)
control_dim(model::DiscreteLinearModel) = size(model.B, 2)
RD.default_signature(model::DiscreteLinearModel) = RD.InPlace()
RD.default_diffmethod(model::DiscreteLinearModel) = RD.UserDefined()
getstatecoeffs(model::DiscreteLinearModel) = model.A
getcontrolcoeffs(model::DiscreteLinearModel) = model.B
getnextstatecoeffs(model::DiscreteLinearModel) = model.C
getaffinecoeffs(model::DiscreteLinearModel) = model.d

isstable(model::DiscreteLinearModel) = RandomLinearModels.isstabled(getstatecoeffs(model))

iscontrollable(model::DiscreteLinearModel) = RandomLinearModels.iscontrollable(model.CA, model.CB)

function checktimestep(model::DiscreteLinearModel, dt)
    if dt ≉ model.dt
        throw(ArgumentError("Invalid time step. Passed in $dt, but the model is defined using a time step of $(model.dt)."))
    end
end

function RD.discrete_dynamics(model::DiscreteLinearModel, x, u, t, h)
    checktimestep(model, h)
    model.CA * x + model.CB * u + model.Cd
end

function RD.discrete_dynamics!(model::DiscreteLinearModel, xn, x, u, t, h)
    checktimestep(model, h)
    xn .= model.Cd
    mul!(xn, model.CA, x, true, true)
    mul!(xn, model.CB, u, true, true)
    nothing
end

function RD.jacobian!(model::DiscreteLinearModel, J, xn, x, u, t, h)
    checktimestep(model, h)
    n,m = RD.dims(model)
    J[:,1:n] .= model.CA
    J[:,n+1:n+m] .= model.CB
    nothing
end

function RD.dynamics_error(model::DiscreteLinearModel,
                           z2::RD.AbstractKnotPoint, z1::RD.AbstractKnotPoint)
    checktimestep(model, RD.timestep(z1))
    x1,u1 = RD.state(z1), RD.control(z1)
    x2 = RD.state(z1)
    model.A * x1 + model.B * u1 + model.C * x2 + model.d
end

function RD.dynamics_error!(model::DiscreteLinearModel, y2, y1,
                            z2::RD.AbstractKnotPoint, z1::RD.AbstractKnotPoint)
    checktimestep(model, RD.timestep(z1))
    x1,u1 = RD.state(z1), RD.control(z1)
    x2 = RD.state(z1)
    y2 .= model.d
    mul!(y2, model.A, x1, true, true)
    mul!(y2, model.B, u1, true, true)
    mul!(y2, model.C, x2, true, true)
    nothing
end

function RD.dynamics_error_jacobian!(model::DiscreteLinearModel, J2, J1, y2, y1,
                                     z2::RD.AbstractKnotPoint, z1::RD.AbstractKnotPoint)
    checktimestep(model, RD.timestep(z1))
    n,m = RD.dims(model)
    J2[:,1:n] .= model.C
    J1[:,1:n] .= model.A
    J1[:,n+1:n+m] .= model.B
    nothing
end


text:
@testset "execution" begin

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

dummy() = return

@testset "cufunction" begin
    @test_throws UndefVarError cufunction(exec_undefined_kernel, ())

    cufunction(dev, dummy, Tuple{})

    # NOTE: other cases are going to be covered by tests below,
    #       as @cuda internally uses cufunction
end

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

@testset "@cuda" begin

@test_throws UndefVarError @cuda undefined()
@test_throws MethodError @cuda dummy(1)

@testset "compilation params" begin
    @cuda dummy()

    @test_throws CuError @cuda threads=2 maxthreads=1 dummy()
    @cuda threads=2 dummy()
end

@testset "reflection" begin
    CUDAnative.code_lowered(dummy, Tuple{})
    CUDAnative.code_typed(dummy, Tuple{})
    CUDAnative.code_warntype(devnull, dummy, Tuple{})
    CUDAnative.code_llvm(devnull, dummy, Tuple{})
    CUDAnative.code_ptx(devnull, dummy, Tuple{})
    CUDAnative.code_sass(devnull, dummy, Tuple{})

    @device_code_lowered @cuda dummy()
    @device_code_typed @cuda dummy()
    @device_code_warntype io=devnull @cuda dummy()
    @device_code_llvm io=devnull @cuda dummy()
    @device_code_ptx io=devnull @cuda dummy()
    @device_code_sass io=devnull @cuda dummy()

    mktempdir() do dir
        @device_code dir=dir @cuda dummy()
    end

    @test_throws ErrorException @device_code_lowered nothing

    # make sure kernel name aliases are preserved in the generated code
    @test occursin("ptxcall_dummy", sprint(io->(@device_code_llvm io=io @cuda dummy())))
    @test occursin("ptxcall_dummy", sprint(io->(@device_code_ptx io=io @cuda dummy())))
    @test occursin("ptxcall_dummy", sprint(io->(@device_code_sass io=io @cuda dummy())))
end

@testset "shared memory" begin
    @cuda shmem=1 dummy()
end

@testset "streams" begin
    s = CuStream()
    @cuda stream=s dummy()
end

@testset "external kernels" begin
    @eval module KernelModule
        export external_dummy
        external_dummy() = return
    end
    import ...KernelModule
    @cuda KernelModule.external_dummy()
    @eval begin
        using ...KernelModule
        @cuda external_dummy()
    end

    @eval module WrapperModule
        using CUDAnative
        @eval dummy() = return
        wrapper() = @cuda dummy()
    end
    WrapperModule.wrapper()
end

@testset "calling device function" begin
    @noinline child(i) = sink(i)
    function parent()
        child(1)
        return
    end

    @cuda parent()
end

end

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

@testset "argument passing" begin

dims = (16, 16)
len = prod(dims)

@testset "manually allocated" begin
    function kernel(input, output)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x

        val = unsafe_load(input, i)
        unsafe_store!(output, val, i)

        return
    end

    input = round.(rand(Float32, dims) * 100)
    output = similar(input)

    input_dev = Mem.upload(input)
    output_dev = Mem.alloc(input)

    @cuda threads=len kernel(Base.unsafe_convert(Ptr{Float32}, input_dev),
                             Base.unsafe_convert(Ptr{Float32}, output_dev))
    Mem.download!(output, output_dev)
    @test input ≈ output
end

@testset "scalar through single-value array" begin
    function kernel(a, x)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        max = gridDim().x * blockDim().x
        if i == max
            _val = unsafe_load(a, i)
            unsafe_store!(x, _val)
        end
        return
    end

    arr = round.(rand(Float32, dims) * 100)
    val = [0f0]

    arr_dev = Mem.upload(arr)
    val_dev = Mem.upload(val)

    @cuda threads=len kernel(Base.unsafe_convert(Ptr{Float32}, arr_dev),
                             Base.unsafe_convert(Ptr{Float32}, val_dev))
    @test arr[dims...] ≈ Mem.download(eltype(val), val_dev)[1]
end

@testset "scalar through single-value array, using device function" begin
    function child(a, i)
        return unsafe_load(a, i)
    end
    @noinline function parent(a, x)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        max = gridDim().x * blockDim().x
        if i == max
            _val = child(a, i)
            unsafe_store!(x, _val)
        end
        return
    end

    arr = round.(rand(Float32, dims) * 100)
    val = [0f0]

    arr_dev = Mem.upload(arr)
    val_dev = Mem.upload(val)

    @cuda threads=len parent(Base.unsafe_convert(Ptr{Float32}, arr_dev),
                             Base.unsafe_convert(Ptr{Float32}, val_dev))
    @test arr[dims...] ≈ Mem.download(eltype(val), val_dev)[1]
end

@testset "tuples" begin
    # issue #7: tuples not passed by pointer

    function kernel(keeps, out)
        if keeps[1]
            unsafe_store!(out, 1)
        else
            unsafe_store!(out, 2)
        end
        return
    end

    keeps = (true,)
    d_out = Mem.alloc(Int)

    @cuda kernel(keeps, Base.unsafe_convert(Ptr{Int}, d_out))
    @test Mem.download(Int, d_out) == [1]
end

@testset "ghost function parameters" begin
    # bug: ghost type function parameters are elided by the compiler

    len = 60
    a = rand(Float32, len)
    b = rand(Float32, len)
    c = similar(a)

    d_a = Mem.upload(a)
    d_b = Mem.upload(b)
    d_c = Mem.alloc(c)

    @eval struct ExecGhost end

    function kernel(ghost, a, b, c)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        unsafe_store!(c, unsafe_load(a,i)+unsafe_load(b,i), i)
        return
    end
    @cuda threads=len kernel(ExecGhost(),
                             Base.unsafe_convert(Ptr{Float32}, d_a),
                             Base.unsafe_convert(Ptr{Float32}, d_b),
                             Base.unsafe_convert(Ptr{Float32}, d_c))
    Mem.download!(c, d_c)
    @test a+b == c

    # bug: ghost type function parameters confused aggregate type rewriting

    function kernel(ghost, out, aggregate)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        unsafe_store!(out, aggregate[1], i)
        return
    end
    @cuda threads=len kernel(ExecGhost(), Base.unsafe_convert(Ptr{Float32}, d_c), (42,))

    Mem.download!(c, d_c)
    @test all(val->val==42, c)
end

@testset "immutables" begin
    # issue #15: immutables not passed by pointer

    function kernel(ptr, b)
        unsafe_store!(ptr, imag(b))
        return
    end

    buf = Mem.upload([0f0])
    x = ComplexF32(2,2)

    @cuda kernel(Base.unsafe_convert(Ptr{Float32}, buf), x)
    @test Mem.download(Float32, buf) == [imag(x)]
end

@testset "automatic recompilation" begin
    buf = Mem.alloc(Int)

    function kernel(ptr)
        unsafe_store!(ptr, 1)
        return
    end

    @cuda kernel(Base.unsafe_convert(Ptr{Int}, buf))
    @test Mem.download(Int, buf) == [1]

    function kernel(ptr)
        unsafe_store!(ptr, 2)
        return
    end

    @cuda kernel(Base.unsafe_convert(Ptr{Int}, buf))
    @test Mem.download(Int, buf) == [2]
end

@testset "non-isbits arguments" begin
    function kernel1(T, i)
        sink(i)
        return
    end
    @cuda kernel1(Int, 1)

    function kernel2(T, i)
        sink(unsafe_trunc(T,i))
        return
    end
    @cuda kernel2(Int, 1.)
end

@testset "splatting" begin
    function kernel(out, a, b)
        unsafe_store!(out, a+b)
        return
    end

    out = [0]
    out_dev = Mem.upload(out)
    out_ptr = Base.unsafe_convert(Ptr{eltype(out)}, out_dev)

    @cuda kernel(out_ptr, 1, 2)
    @test Mem.download(eltype(out), out_dev)[1] == 3

    all_splat = (out_ptr, 3, 4)
    @cuda kernel(all_splat...)
    @test Mem.download(eltype(out), out_dev)[1] == 7

    partial_splat = (5, 6)
    @cuda kernel(out_ptr, partial_splat...)
    @test Mem.download(eltype(out), out_dev)[1] == 11
end

@testset "object invoke" begin
    # this mimics what is generated by closure conversion

    @eval struct KernelObject{T} <: Function
        val::T
    end

    function (self::KernelObject)(a)
        unsafe_store!(a, self.val)
        return
    end

    function outer(a, val)
       inner = KernelObject(val)
       @cuda inner(a)
    end

    a = [1.]
    a_dev = Mem.upload(a)

    outer(Base.unsafe_convert(Ptr{Float64}, a_dev), 2.)

    @test Mem.download(eltype(a), a_dev) ≈ [2.]
end

@testset "closures" begin
    function outer(a_dev, val)
       function inner(a)
            # captures `val`
            unsafe_store!(a, val)
            return
       end
       @cuda inner(Base.unsafe_convert(Ptr{Float64}, a_dev))
    end

    a = [1.]
    a_dev = Mem.upload(a)

    outer(a_dev, 2.)

    @test Mem.download(eltype(a), a_dev) ≈ [2.]
end

@testset "conversions" begin
    @eval struct Host   end
    @eval struct Device end

    CUDAnative.cudaconvert(a::Host) = Device()

    Base.convert(::Type{Int}, ::Host)   = 1
    Base.convert(::Type{Int}, ::Device) = 2

    out = [0]

    # convert arguments
    out_dev = Mem.upload(out)
    let arg = Host()
        @test Mem.download(eltype(out), out_dev) ≈ [0]
        function kernel(arg, out)
            unsafe_store!(out, convert(Int, arg))
            return
        end
        @cuda kernel(arg, Base.unsafe_convert(Ptr{Int}, out_dev))
        @test Mem.download(eltype(out), out_dev) ≈ [2]
    end

    # convert tuples
    out_dev = Mem.upload(out)
    let arg = (Host(),)
        @test Mem.download(eltype(out), out_dev) ≈ [0]
        function kernel(arg, out)
            unsafe_store!(out, convert(Int, arg[1]))
            return
        end
        @cuda kernel(arg, Base.unsafe_convert(Ptr{Int}, out_dev))
        @test Mem.download(eltype(out), out_dev) ≈ [2]
    end

    # convert named tuples
    out_dev = Mem.upload(out)
    let arg = (a=Host(),)
        @test Mem.download(eltype(out), out_dev) ≈ [0]
        function kernel(arg, out)
            unsafe_store!(out, convert(Int, arg.a))
            return
        end
        @cuda kernel(arg, Base.unsafe_convert(Ptr{Int}, out_dev))
        @test Mem.download(eltype(out), out_dev) ≈ [2]
    end

    # don't convert structs
    out_dev = Mem.upload(out)
    @eval struct Nested
        a::Host
    end
    let arg = Nested(Host())
        @test Mem.download(eltype(out), out_dev) ≈ [0]
        function kernel(arg, out)
            unsafe_store!(out, convert(Int, arg.a))
            return
        end
        @cuda kernel(arg, Base.unsafe_convert(Ptr{Int}, out_dev))
        @test Mem.download(eltype(out), out_dev) ≈ [1]
    end
end

@testset "argument count" begin
    val = [0]
    val_dev = Mem.upload(val)
    ptr = Base.unsafe_convert(Ptr{Int}, val_dev)
    for i in (1, 10, 20, 35)
        variables = ('a':'z'..., 'A':'Z'...)
        params = [Symbol(variables[j]) for j in 1:i]
        # generate a kernel
        body = quote
            function kernel($(params...))
                unsafe_store!($ptr, $(Expr(:call, :+, params...)))
                return
            end
        end
        eval(body)
        args = [j for j in 1:i]
        call = Expr(:call, :kernel, args...)
        cudacall = :(@cuda $call)
        eval(cudacall)
        @test Mem.download(eltype(val), val_dev)[1] == sum(args)
    end
end

end

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

end


text:
module CUFFT
# FFT interface for CuArrays

using ..CuArrays: libcufft, configured, CuArray
import AbstractFFTs: plan_fft, plan_fft!, plan_bfft, plan_bfft!,
    plan_rfft, plan_brfft, plan_inv, normalization, fft, bfft, ifft, rfft,
    Plan, ScaledPlan
import Base: show, *, convert, unsafe_convert, size, strides, ndims
import LinearAlgebra: mul!
import Base.Sys: WORD_SIZE

using LinearAlgebra

include("libcufft_types.jl")
include("error.jl")

const cufftNumber = Union{cufftDoubleReal,cufftReal,cufftDoubleComplex,cufftComplex}
# note trailing s to deconflict w/ header file
const cufftReals = Union{cufftDoubleReal,cufftReal}
const cufftComplexes = Union{cufftDoubleComplex,cufftComplex}
const cufftDouble = Union{cufftDoubleReal,cufftDoubleComplex}
const cufftSingle = Union{cufftReal,cufftComplex}
const cufftTypeDouble = Union{Type{cufftDoubleReal},Type{cufftDoubleComplex}}
const cufftTypeSingle = Union{Type{cufftReal},Type{cufftComplex}}

include("genericfft.jl")

# K is a flag for forward/backward
# also used as an alias for r2c/c2r

abstract type CuFFTPlan{T<:cufftNumber, K, inplace} <: Plan{T} end

mutable struct cCuFFTPlan{T<:cufftNumber,K,inplace,N} <: CuFFTPlan{T,K,inplace}
    plan::cufftHandle_t
    sz::NTuple{N,Int} # Julia size of input array
    osz::NTuple{N,Int} # Julia size of output array
    xtype::Int
    region::Any
    pinv::ScaledPlan # required by AbstractFFT API

    function cCuFFTPlan{T,K,inplace,N}(plan::cufftHandle_t, X::CuArray{T,N},
                                       sizey::Tuple, region, xtype::Integer
                                       ) where {T<:cufftNumber,K,inplace,N}
        # maybe enforce consistency of sizey
        p = new(plan, size(X), sizey, xtype, region)
        finalizer(destroy_plan, p)
        p
    end
end

cCuFFTPlan(plan,X,region,xtype::Integer) = cCuFFTPlan(plan,X,size(X),region,xtype)

mutable struct rCuFFTPlan{T<:cufftNumber,K,inplace,N} <: CuFFTPlan{T,K,inplace}
    plan::cufftHandle_t
    sz::NTuple{N,Int} # Julia size of input array
    osz::NTuple{N,Int} # Julia size of output array
    xtype::Int
    region::Any
    pinv::ScaledPlan # required by AbstractFFT API

    function rCuFFTPlan{T,K,inplace,N}(plan::cufftHandle_t, X::CuArray{T,N},
                                       sizey::Tuple, region, xtype::Integer
                                       ) where {T<:cufftNumber,K,inplace,N}
        # maybe enforce consistency of sizey
        p = new(plan, size(X), sizey, xtype, region)
        finalizer(destroy_plan, p)
        p
    end
end

rCuFFTPlan(plan,X,region,xtype::Integer) = rCuFFTPlan(plan,X,size(X),region,xtype)

const xtypenames = Dict{cufftType,String}(CUFFT_R2C => "real-to-complex",
                                          CUFFT_C2R => "complex-to-real",
                                          CUFFT_C2C => "complex",
                                          CUFFT_D2Z => "d.p. real-to-complex",
                                          CUFFT_Z2D => "d.p. complex-to-real",
                                          CUFFT_Z2Z => "d.p. complex")

function showfftdims(io, sz, T)
    if isempty(sz)
        print(io,"0-dimensional")
    elseif length(sz) == 1
        print(io, sz[1], "-element")
    else
        print(io, join(sz, "×"))
    end
    print(io, " CuArray of ", T)
end

function show(io::IO, p::CuFFTPlan{T,K,inplace}) where {T,K,inplace}
    print(io, inplace ? "CUFFT in-place " : "CUFFT ",
          xtypenames[p.xtype],
          K == CUFFT_FORWARD ? " forward" : " backward",
          " plan for ")
    showfftdims(io, p.sz, T)
end

# Note: we don't implement padded storage dimensions
function _mkplan(xtype, xdims, region)
    nrank = length(region)
    sz = [xdims[i] for i in region]
    csz = copy(sz)
    csz[1] = div(sz[1],2) + 1
    batch = prod(xdims) ÷ prod(sz)

    pp = Ref{cufftHandle_t}()
    if (nrank == 1) && (batch == 1)
        @check ccall((:cufftPlan1d,libcufft),cufftStatus_t,
                     (Ptr{cufftHandle_t}, Cint, cufftType, Cint),
                     pp, sz[1], xtype, 1)
    elseif (nrank == 2) && (batch == 1)
        @check ccall((:cufftPlan2d,libcufft),cufftStatus_t,
                     (Ptr{cufftHandle_t}, Cint, Cint, cufftType),
                     pp, sz[2], sz[1], xtype)
    elseif (nrank == 3) && (batch == 1)
        @check ccall((:cufftPlan3d,libcufft),cufftStatus_t,
                     (Ptr{cufftHandle_t}, Cint, Cint, Cint, cufftType),
                     pp, sz[3], sz[2], sz[1], xtype)

    else
        rsz = (length(sz) > 1) ? rsz = reverse(sz) : sz
        if ((region...,) == ((1:nrank)...,))
            # handle simple case ... simply! (for robustness)
            @check ccall((:cufftPlanMany,libcufft),cufftStatus_t,
                         (Ptr{cufftHandle_t}, Cint, Ptr{Cint}, # rank, dims
                          Ptr{Cint}, Cint, Cint, # nembed,stride,dist (input)
                          Ptr{Cint}, Cint, Cint, # nembed,stride,dist (output)
                          cufftType, Cint),
                         pp, nrank, Cint[rsz...], C_NULL, 1, 1, C_NULL, 1, 1,
                         xtype, batch)
        else
            if nrank==1 || all(diff(collect(region)) .== 1)
                # _stride: successive elements in innermost dimension
                # _dist: distance between first elements of batches
                if region[1] == 1
                    istride = 1
                    idist = prod(sz)
                    cdist = prod(csz)
                else
                    if region[end] != length(xdims)
                        throw(ArgumentError("batching dims must be sequential"))
                    end
                    istride = prod(xdims[1:region[1]-1])
                    idist = 1
                    cdist = 1
                end
                inembed = Cint[rsz...]
                cnembed = (length(csz) > 1) ? Cint[reverse(csz)...] : Cint[csz[1]]
                ostride = istride
                if xtype == CUFFT_R2C || xtype == CUFFT_D2Z
                    odist = cdist
                    onembed = cnembed
                else
                    odist = idist
                    onembed = inembed
                end
                if xtype == CUFFT_C2R || xtype == CUFFT_Z2D
                    idist = cdist
                    inembed = cnembed
                end
            else
                if any(diff(collect(region)) .< 1)
                    throw(ArgumentError("region must be an increasing sequence"))
                end
                cdims = collect(xdims)
                cdims[region[1]] = div(cdims[region[1]],2)+1

                if region[1] == 1
                    istride = 1
                    ii=1
                    while (ii < nrank) && (region[ii] == region[ii+1]-1)
                        ii += 1
                    end
                    idist = prod(xdims[1:ii])
                    cdist = prod(cdims[1:ii])
                    ngaps = 0
                else
                    istride = prod(xdims[1:region[1]-1])
                    idist = 1
                    cdist = 1
                    ngaps = 1
                end
                nem = ones(Int,nrank)
                cem = ones(Int,nrank)
                id = 1
                for ii=1:nrank-1
                    if region[ii+1] > region[ii]+1
                        ngaps += 1
                    end
                    while id < region[ii+1]
                        nem[ii] *= xdims[id]
                        cem[ii] *= cdims[id]
                        id += 1
                    end
                    @assert nem[ii] >= sz[ii]
                end
                if region[end] < length(xdims)
                    ngaps += 1
                end
                # CUFFT represents batches by a single stride (_dist)
                # so we must verify that region is consistent with this:
                if ngaps > 1
                    throw(ArgumentError("batch regions must be sequential"))
                end

                inembed = Cint[reverse(nem)...]
                cnembed = Cint[reverse(cem)...]
                ostride = istride
                if xtype == CUFFT_R2C || xtype == CUFFT_D2Z
                    odist = cdist
                    onembed = cnembed
                else
                    odist = idist
                    onembed = inembed
                end
                if xtype == CUFFT_C2R || xtype == CUFFT_Z2D
                    idist = cdist
                    inembed = cnembed
                end
            end
            @check ccall((:cufftPlanMany,libcufft),cufftStatus_t,
                         (Ptr{cufftHandle_t}, Cint, Ptr{Cint}, # rank, dims
                          Ptr{Cint}, Cint, Cint, # nembed,stride,dist (input)
                          Ptr{Cint}, Cint, Cint, # nembed,stride,dist (output)
                          cufftType, Cint),
                         pp, nrank, Cint[rsz...],
                         inembed, istride, idist, onembed, ostride, odist,
                         xtype, batch)
        end
    end
    pp[]
end

# this is used implicitly in the unsafe_execute methods below:
unsafe_convert(::Type{cufftHandle_t}, p::CuFFTPlan) = p.plan

convert(::Type{cufftHandle_t}, p::CuFFTPlan) = p.plan

destroy_plan(plan::CuFFTPlan) =
    ccall((:cufftDestroy,libcufft), Nothing, (cufftHandle_t,), plan.plan)

function assert_applicable(p::CuFFTPlan{T,K}, X::CuArray{T}) where {T,K}
    (size(X) == p.sz) ||
        throw(ArgumentError("CuFFT plan applied to wrong-size input"))
end

function assert_applicable(p::CuFFTPlan{T,K}, X::CuArray{T}, Y::CuArray{Ty}) where {T,K,Ty}
    assert_applicable(p, X)
    (size(Y) == p.osz) ||
        throw(ArgumentError("CuFFT plan applied to wrong-size output"))
    # type errors should be impossible by dispatch, but just in case:
    if p.xtype ∈ [CUFFT_C2R, CUFFT_Z2D]
        (Ty == real(T)) ||
            throw(ArgumentError("Type mismatch for argument Y"))
    elseif p.xtype ∈ [CUFFT_R2C, CUFFT_D2Z]
        (Ty == complex(T)) ||
            throw(ArgumentError("Type mismatch for argument Y"))
    else
        (Ty == T) ||
            throw(ArgumentError("Type mismatch for argument Y"))
    end
end

function unsafe_execute!(plan::cCuFFTPlan{cufftComplex,K,true,N},
                         x::CuArray{cufftComplex,N}) where {K,N}
    @assert plan.xtype == CUFFT_C2C
    @check ccall((:cufftExecC2C,libcufft), cufftStatus_t,
                 (cufftHandle_t, Ptr{cufftComplex}, Ptr{cufftComplex},
                  Cint),
                 plan, x, x, K)
end
function unsafe_execute!(plan::rCuFFTPlan{cufftComplex,K,true,N},
                         x::CuArray{cufftComplex,N}) where {K,N}
    @assert plan.xtype == CUFFT_C2R
    @check ccall((:cufftExecC2R,libcufft), cufftStatus_t,
                 (cufftHandle_t, Ptr{cufftComplex}, Ptr{cufftComplex}),
                 plan, x, x)
end

function unsafe_execute!(plan::cCuFFTPlan{cufftComplex,K,false,N},
                         x::CuArray{cufftComplex,N}, y::CuArray{cufftComplex}
                         ) where {K,N}
    @assert plan.xtype == CUFFT_C2C
    @check ccall((:cufftExecC2C,libcufft), cufftStatus_t,
                 (cufftHandle_t, Ptr{cufftComplex}, Ptr{cufftComplex}, Cint),
                 plan, x, y, K)
end
function unsafe_execute!(plan::rCuFFTPlan{cufftComplex,K,false,N},
                         x::CuArray{cufftComplex,N}, y::CuArray{cufftReal}
                         ) where {K,N}
    @assert plan.xtype == CUFFT_C2R
    @check ccall((:cufftExecC2R,libcufft), cufftStatus_t,
                 (cufftHandle_t, Ptr{cufftComplex}, Ptr{cufftReal}),
                 plan, x, y)
end

function unsafe_execute!(plan::rCuFFTPlan{cufftReal,K,false,N},
                         x::CuArray{cufftReal,N}, y::CuArray{cufftComplex,N}
                         ) where {K,N}
    @assert plan.xtype == CUFFT_R2C
    @check ccall((:cufftExecR2C,libcufft), cufftStatus_t,
                 (cufftHandle_t, Ptr{cufftReal}, Ptr{cufftComplex}),
                 plan, x, y)
end

# double prec.
function unsafe_execute!(plan::cCuFFTPlan{cufftDoubleComplex,K,true,N},
                         x::CuArray{cufftDoubleComplex,N}) where {K,N}
    @assert plan.xtype == CUFFT_Z2Z
    @check ccall((:cufftExecZ2Z,libcufft), cufftStatus_t,
                 (cufftHandle_t, Ptr{cufftDoubleComplex}, Ptr{cufftDoubleComplex},
                  Cint),
                 plan, x, x, K)
end
function unsafe_execute!(plan::rCuFFTPlan{cufftDoubleComplex,K,true,N},
                         x::CuArray{cufftDoubleComplex,N}) where {K,N}
    @assert plan.xtype == CUFFT_Z2D
    @check ccall((:cufftExecZ2D,libcufft), cufftStatus_t,
                 (cufftHandle_t, Ptr{cufftDoubleComplex}, Ptr{cufftDoubleComplex}),
                 plan, x, x)
end

function unsafe_execute!(plan::cCuFFTPlan{cufftDoubleComplex,K,false,N},
                         x::CuArray{cufftDoubleComplex,N}, y::CuArray{cufftDoubleComplex}
                         ) where {K,N}
    @assert plan.xtype == CUFFT_Z2Z
    @check ccall((:cufftExecZ2Z,libcufft), cufftStatus_t,
                 (cufftHandle_t, Ptr{cufftDoubleComplex}, Ptr{cufftDoubleComplex}, Cint),
                 plan, x, y, K)
end
function unsafe_execute!(plan::rCuFFTPlan{cufftDoubleComplex,K,false,N},
                         x::CuArray{cufftDoubleComplex,N}, y::CuArray{cufftDoubleReal}
                         ) where {K,N}
    @assert plan.xtype == CUFFT_Z2D
    @check ccall((:cufftExecZ2D,libcufft), cufftStatus_t,
                 (cufftHandle_t, Ptr{cufftDoubleComplex}, Ptr{cufftDoubleReal}),
                 plan, x, y)
end

function unsafe_execute!(plan::rCuFFTPlan{cufftDoubleReal,K,false,N},
                         x::CuArray{cufftDoubleReal,N}, y::CuArray{cufftDoubleComplex,N}
                         ) where {K,N}
    @assert plan.xtype == CUFFT_D2Z
    @check ccall((:cufftExecD2Z,libcufft), cufftStatus_t,
                 (cufftHandle_t, Ptr{cufftDoubleReal}, Ptr{cufftDoubleComplex}),
                 plan, x, y)
end

###############
# Public API  #
###############

# region is an iterable subset of dimensions
# spec. an integer, range, tuple, or array

# inplace complex
function plan_fft!(X::CuArray{T,N}, region) where {T<:cufftComplexes,N}
    K = CUFFT_FORWARD
    inplace = true
    xtype = (T == cufftComplex) ? CUFFT_C2C : CUFFT_Z2Z

    pp = _mkplan(xtype, size(X), region)

    cCuFFTPlan{T,K,inplace,N}(pp, X, size(X), region, xtype)
end

function plan_bfft!(X::CuArray{T,N}, region) where {T<:cufftComplexes,N}
    K = CUFFT_INVERSE
    inplace = true
    xtype =  (T == cufftComplex) ? CUFFT_C2C : CUFFT_Z2Z

    pp = _mkplan(xtype, size(X), region)

    cCuFFTPlan{T,K,inplace,N}(pp, X, size(X), region, xtype)
end

# out-of-place complex
function plan_fft(X::CuArray{T,N}, region) where {T<:cufftComplexes,N}
    K = CUFFT_FORWARD
    xtype =  (T == cufftComplex) ? CUFFT_C2C : CUFFT_Z2Z
    inplace = false

    pp = _mkplan(xtype, size(X), region)

    cCuFFTPlan{T,K,inplace,N}(pp, X, size(X), region, xtype)
end

function plan_bfft(X::CuArray{T,N}, region) where {T<:cufftComplexes,N}
    K = CUFFT_INVERSE
    inplace = false
    xtype =  (T == cufftComplex) ? CUFFT_C2C : CUFFT_Z2Z

    pp = _mkplan(xtype, size(X), region)

    cCuFFTPlan{T,K,inplace,N}(pp, X, size(X), region, xtype)
end

# out-of-place real-to-complex
function plan_rfft(X::CuArray{T,N}, region) where {T<:cufftReals,N}
    K = CUFFT_FORWARD
    inplace = false
    xtype =  (T == cufftReal) ? CUFFT_R2C : CUFFT_D2Z

    pp = _mkplan(xtype, size(X), region)

    ydims = collect(size(X))
    ydims[region[1]] = div(ydims[region[1]],2)+1

    rCuFFTPlan{T,K,inplace,N}(pp, X, (ydims...,), region, xtype)
end

function plan_brfft(X::CuArray{T,N}, d::Integer, region::Any) where {T<:cufftComplexes,N}
    K = CUFFT_INVERSE
    inplace = false
    xtype =  (T == cufftComplex) ? CUFFT_C2R : CUFFT_Z2D
    ydims = collect(size(X))
    ydims[region[1]] = d

    pp = _mkplan(xtype, (ydims...,), region)

    rCuFFTPlan{T,K,inplace,N}(pp, X, (ydims...,), region, xtype)
end

# FIXME: plan_inv methods allocate needlessly (to provide type parameters)
# Perhaps use FakeArray types to avoid this.

function plan_inv(p::cCuFFTPlan{T,CUFFT_FORWARD,inplace,N}) where {T,N,inplace}
    X = CuArray{T}(undef, p.sz)
    pp = _mkplan(p.xtype, p.sz, p.region)
    ScaledPlan(cCuFFTPlan{T,CUFFT_INVERSE,inplace,N}(pp, X, p.sz, p.region,
                                                     p.xtype),
               normalization(X, p.region))
end

function plan_inv(p::cCuFFTPlan{T,CUFFT_INVERSE,inplace,N}) where {T,N,inplace}
    X = CuArray{T}(undef, p.sz)
    pp = _mkplan(p.xtype, p.sz, p.region)
    ScaledPlan(cCuFFTPlan{T,CUFFT_FORWARD,inplace,N}(pp, X, p.sz, p.region,
                                                     p.xtype),
               normalization(X, p.region))
end

function plan_inv(p::rCuFFTPlan{T,CUFFT_INVERSE,inplace,N}
                  ) where {T<:cufftComplexes,N,inplace}
    X = CuArray{real(T)}(undef, p.osz)
    Y = CuArray{T}(undef, p.sz)
    xtype = p.xtype == CUFFT_C2R ? CUFFT_R2C : CUFFT_D2Z
    pp = _mkplan(xtype, p.osz, p.region)
    ScaledPlan(rCuFFTPlan{real(T),CUFFT_FORWARD,inplace,N}(pp, X, p.sz, p.region,
                                                     xtype),
               normalization(X, p.region))
end

function plan_inv(p::rCuFFTPlan{T,CUFFT_FORWARD,inplace,N}
                  ) where {T<:cufftReals,N,inplace}
    X = CuArray{complex(T)}(undef, p.osz)
    Y = CuArray{T}(undef, p.sz)
    xtype = p.xtype == CUFFT_R2C ? CUFFT_C2R : CUFFT_Z2D
    pp = _mkplan(xtype, p.sz, p.region)
    ScaledPlan(rCuFFTPlan{complex(T),CUFFT_INVERSE,inplace,N}(pp, X, p.sz,
                                                              p.region, xtype),
               normalization(Y, p.region))
end

# The rest of the standard API

size(p::CuFFTPlan) = p.sz

function mul!(y::CuArray{Ty}, p::CuFFTPlan{T,K,false}, x::CuArray{T}
                  ) where {Ty,T,K}
    assert_applicable(p,x,y)
    unsafe_execute!(p,x,y)
    return y
end

function *(p::cCuFFTPlan{T,K,true,N}, x::CuArray{T,N}) where {T,K,N}
    assert_applicable(p,x)
    unsafe_execute!(p,x)
    x
end

function *(p::rCuFFTPlan{T,CUFFT_FORWARD,false,N}, x::CuArray{T,N}
           ) where {T<:cufftReals,N}
    @assert p.xtype ∈ [CUFFT_R2C,CUFFT_D2Z]
    y = CuArray{complex(T),N}(undef, p.osz)
    mul!(y,p,x)
    y
end

function *(p::rCuFFTPlan{T,CUFFT_INVERSE,false,N}, x::CuArray{T,N}
           ) where {T<:cufftComplexes,N}
    @assert p.xtype ∈ [CUFFT_C2R,CUFFT_Z2D]
    y = CuArray{real(T),N}(undef, p.osz)
    mul!(y,p,x)
    y
end

function *(p::cCuFFTPlan{T,K,false,N}, x::CuArray{T,N}) where {T,K,N}
    y = CuArray{T,N}(undef, p.osz)
    mul!(y,p,x)
    y
end

end # module


text:
function diagoSlide(square::UInt64)
    visited = EMPTY

    #NW
    for shift = 9:9:54
        square & (MASK_FILES[1] | MASK_RANKS[8]) != EMPTY ? break : visit = square << shift
        visit & (MASK_FILES[1] | MASK_RANKS[8]) != EMPTY ? break : visited |= visit
    end
    #SE
    for shift = 9:9:54
        square & (MASK_FILES[8] | MASK_RANKS[1]) != EMPTY ? break : visit = square >> shift
        visit & (MASK_FILES[8] | MASK_RANKS[1]) != EMPTY ? break : visited |= visit
    end
    #NE
    for shift = 7:7:42
        square & (MASK_FILES[8] | MASK_RANKS[8]) != EMPTY ? break : visit = square << shift
        visit & (MASK_FILES[8] | MASK_RANKS[8]) != EMPTY ? break : visited |= visit
    end
    #SW
    for shift = 7:7:42
        square & (MASK_FILES[1] | MASK_RANKS[1]) != EMPTY ? break : visit = square >> shift
        visit & (MASK_FILES[1] | MASK_RANKS[1]) != EMPTY ? break : visited |= visit
    end

    return visited
end

function bishopMasksGen()
    masks = Dict{UInt64,UInt64}()
    [push!(masks, s=>diagoSlide(s)) for s in values(PGN2UINT)]
    return masks
end
#diagonal masks
const BISHOP_MASKS = bishopMasksGen()

function diagoAttack(square::UInt64, occupancy::UInt64)
    visited = EMPTY

    #NW
    for shift = 9:9:63
        square & (MASK_FILES[1] | MASK_RANKS[8]) != EMPTY ? break : visit = square << shift
        visited |= visit
        if visit & (MASK_FILES[1] | MASK_RANKS[8]) != EMPTY || visit & occupancy != EMPTY; break end
    end
    #SE
    for shift = 9:9:63
        square & (MASK_FILES[8] | MASK_RANKS[1]) != EMPTY ? break : visit = square >> shift
        visited |= visit
        if visit & (MASK_FILES[8] | MASK_RANKS[1]) != EMPTY || visit & occupancy != EMPTY; break end
    end
    #NE
    for shift = 7:7:63
        square & (MASK_FILES[8] | MASK_RANKS[8]) != EMPTY ? break : visit = square << shift
        visited |= visit
        if visit & (MASK_FILES[8] | MASK_RANKS[8]) != EMPTY || visit & occupancy != EMPTY; break end
    end
    #SW
    for shift = 7:7:63
        square & (MASK_FILES[1] | MASK_RANKS[1]) != EMPTY ? break : visit = square >> shift
        visited |= visit
        if visit & (MASK_FILES[1] | MASK_RANKS[1]) != EMPTY || visit & occupancy != EMPTY; break end
    end

    return visited
end


text:
#include("../dataset_reader.jl")
const datasetDemo_labels = Float64[-1, 1]
const datasetDemo_nfeatures = 3
const datasetDemo_filename = "dados1.txt"

function get_numerical_datasetDemo(path::String="")
    data = get_data(joinpath(path, datasetDemo_filename), dlm=" ", delete_first_line=false)

    conversion_pattern = [Float64, Float64, Float64, Float64]
    convert_values(data, conversion_pattern)
end

function get_datasetDemo(path::String=""; shuffle::Bool=false, test_perc::Float64=0.2,
                         normalise::Bool=false, separate_test_data::Bool=true, balance_classes::Bool=false)
    data = get_numerical_datasetDemo(path)
    separate_dataset(data, test_percentage=test_perc, shuffle=shuffle, normalise=normalise,
                     separate_test_data=separate_test_data, balance_classes=balance_classes)
end

function get_dataset0_setup(path::String=""; normalise::Bool=false, separate::Bool=false,
                            balance::Bool=false)
    dataset = get_datasetDemo(path, normalise=normalise, separate_test_data=separate, balance_classes=balance)
    if separate
        SeparatedDataset(dataset..., datasetDemo_labels, datasetDemo_nfeatures, 0)
    else
        Dataset(dataset..., datasetDemo_labels, datasetDemo_nfeatures, 0)
    end
end

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
# Factory method - create objects that conforms to certain interface

# Sample interface

"""
 format(::Formatter, x::T) where {T <: Number}

Format a number `x` using the specified formatter.
Returns a string.
"""
function format end

# Example

module FactoryExample

abstract type Formatter end
struct IntegerFormatter <: Formatter end
struct FloatFormatter <: Formatter end

# constructor (the "factory")
# choose implementation using regular dispatch mechanism
formatter(::Type{T}) where {T <: Integer} = IntegerFormatter()
formatter(::Type{T}) where {T <: AbstractFloat} = FloatFormatter()
formatter(::Type{T}) where T = error("No formatter defined for type $T")

# sample formatters
# (concrete implementation for the interface with "format" function)
using Printf
format(nf::IntegerFormatter, x) = @sprintf("%d", x)
format(nf::FloatFormatter, x) = @sprintf("%.2f", x)

# client usage
# (no direct coupling with the concrete types)
function test()
    nf = formatter(Int)
    println(format(nf, 1234))
    nf = formatter(Float64)
    println(format(nf, 1234))
end

end # module

using .FactoryExample
FactoryExample.test()


text:
export fit, predict, fitted_params, JLBoostMLJModel, JLBoostClassifier, JLBoostRegressor, JLBoostCount

import MLJBase
import MLJBase: Probabilistic, Deterministic, clean!, fit, predict, fitted_params, load_path, Table
import MLJBase: package_name, package_uuid, package_url, is_pure_julia, package_license
import MLJBase: input_scitype, target_scitype, docstring, UnivariateFinite

using ScientificTypes: Continuous, OrderedFactor, Count, Multiclass, Finite

using LossFunctions: PoissonLoss, L2DistLoss
using JLBoost: LogitLogLoss, jlboost, AUC, gini, feature_importance

using DataFrames: DataFrame, nrow, levels, categorical

# supervised determinstinistic model
#abstract type JLBoostMLJModel <: Supervised end

mutable struct JLBoostClassifier <: Probabilistic
    loss
    nrounds
    subsample
    eta
    max_depth
    min_child_weight
    lambda
    gamma
    colsample_bytree
end

mutable struct JLBoostRegressor <: Deterministic
    loss
    nrounds
    subsample
    eta
    max_depth
    min_child_weight
    lambda
    gamma
    colsample_bytree
end
mutable struct JLBoostCount <: Deterministic
    loss
    nrounds
    subsample
    eta
    max_depth
    min_child_weight
    lambda
    gamma
    colsample_bytree
end

"""
    JLBoostClassifier(;
        loss = LogitLogLoss(),
        nrounds = 1,
        subsample = 1,
        eta = 1,
        max_depth = 6,
        min_child_weight = 1,
        lambda = 0,
        gamma = 0,
        colsample_bytree = 1)

Return an MLJ.jl compatible Model. The parameters are the same as `jlboost`. See `?jlboost`
"""
JLBoostClassifier(;
    loss = LogitLogLoss(),
    nrounds = 1,
    subsample = 1.0,
    eta = 1.0,
    max_depth = 6,
    min_child_weight = 1.0,
    lambda = 0.0,
    gamma = 0.0,
    colsample_bytree = 1) = JLBoostClassifier(loss, nrounds, subsample, eta, max_depth, min_child_weight, lambda, gamma, colsample_bytree)

"""
    JLBoostRegressor(;
        loss = L2DistLoss(),
        nrounds = 1,
        subsample = 1,
        eta = 1,
        max_depth = 6,
        min_child_weight = 1,
        lambda = 0,
        gamma = 0,
        colsample_bytree = 1)

Return an MLJ.jl compatible Model. The parameters are the same as `jlboost`. See `?jlboost`
"""
JLBoostRegressor(;
    loss = L2DistLoss(),
    nrounds = 1,
    subsample = 1.0,
    eta = 1.0,
    max_depth = 6,
    min_child_weight = 1.0,
    lambda = 0.0,
    gamma = 0.0,
    colsample_bytree = 1) = JLBoostRegressor(loss, nrounds, subsample, eta, max_depth, min_child_weight, lambda, gamma, colsample_bytree)

"""
    JLBoostCount(;
        loss = L2DistLoss(),
        nrounds = 1,
        subsample = 1,
        eta = 1,
        max_depth = 6,
        min_child_weight = 1,
        lambda = 0,
        gamma = 0,
        colsample_bytree = 1)

Return an MLJ.jl compatible Model. The parameters are the same as `jlboost`. See `?jlboost`
"""
JLBoostCount(;
    loss = PoissonLoss(),
    nrounds = 1,
    subsample = 1.0,
    eta = 1.0,
    max_depth = 6,
    min_child_weight = 1.0,
    lambda = 0.0,
    gamma = 0.0,
    colsample_bytree = 1) = JLBoostCount(loss, nrounds, subsample, eta, max_depth, min_child_weight, lambda, gamma, colsample_bytree)

const JLBoostMLJModel = Union{JLBoostClassifier, JLBoostRegressor, JLBoostCount}

# see https://alan-turing-institute.github.io/MLJ.jl/stable/adding_models_for_general_use/#The-fit-method-1
fit(model::Union{JLBoostRegressor, JLBoostCount}, verbosity::Int, X, y::AbstractVector) = begin
    y = DataFrame(__y__ = y)
    df = hcat(X, y)

    target = names(y)[1]
    features = setdiff(names(X), names(y))
    warm_start = fill(0.0, nrow(X))
    fitresult = jlboost(df, target, features, warm_start , model.loss;
        nrounds = model.nrounds, subsample = model.subsample, eta = model.eta,
        colsample_bytree = model.colsample_bytree, max_depth = model.max_depth,
        min_child_weight = model.min_child_weight, lambda = model.lambda,
        gamma = model.gamma, verbose = verbosity >= 1
     )

     return (
            fitresult = fitresult,
            cache = nothing,
            report = (
                placeholder = true
                # feature_importance3 = feature_importance(fitresult, df)
            )
        )
end

fit(model::JLBoostClassifier, verbosity::Int, X, y::AbstractVector) = begin
    ydf = DataFrame(__y__ = y)
    df = hcat(X, ydf)

    target = names(ydf)[1]
    features = setdiff(names(X), names(ydf))
    warm_start = fill(0.0, nrow(X))
    treemodel = jlboost(df, target, features, warm_start , model.loss;
        nrounds = model.nrounds, subsample = model.subsample, eta = model.eta,
        colsample_bytree = model.colsample_bytree, max_depth = model.max_depth,
        min_child_weight = model.min_child_weight, lambda = model.lambda,
        gamma = model.gamma, verbose = verbosity >= 1
     )

     fitresult = (treemodel = treemodel, target_levels = levels(y))

     if length(levels(y)) == 2
         res = (
            fitresult = fitresult,
            cache = nothing,
            report = (
                AUC = abs(AUC(JLBoost.predict(fitresult.treemodel, X), y)),
                feature_importance = feature_importance(fitresult.treemodel, df)
            )
        )
        return res
    else
         return (
            fitresult = fitresult,
            cache = nothing,
            report = (
                feature_importance = feature_importance(fitresult.treemodel, df),
            )
        )
    end
end

# see https://alan-turing-institute.github.io/MLJ.jl/stable/adding_models_for_general_use/#The-fitted_params-method-1
fitted_params(model::JLBoostMLJModel, fitresult) = (fitresult = fitresult.treemodel, trees = trees(fitresult.treemodel))

#  seehttps://alan-turing-institute.github.io/MLJ.jl/stable/adding_models_for_general_use/#The-predict-method-1
function MLJBase.predict(model::JLBoostClassifier, fitresult, Xnew)
    res = JLBoost.predict(fitresult.treemodel, Xnew)
    p = 1 ./ (1 .+ exp.(-res))
    levels_cate = categorical(fitresult.target_levels)
    [UnivariateFinite(levels_cate, [p, 1-p]) for p in p]
end

function MLJBase.predict(model::JLBoostMLJModel, fitresult, Xnew)
    JLBoost.predict(fitresult, Xnew)
end

# see https://alan-turing-institute.github.io/MLJ.jl/stable/adding_models_for_general_use/#Trait-declarations-1
input_scitype(::Type{<:JLBoostMLJModel}) = Table(Union{Continuous, OrderedFactor, Count})

target_scitype(::Type{<:JLBoostClassifier}) = AbstractVector{<:Union{Finite{2}, Continuous, Count}} #AbstractVector{<:Multiclass{2}}
target_scitype(::Type{<:JLBoostRegressor}) = AbstractVector{<:Continuous}
target_scitype(::Type{<:JLBoostCount}) = AbstractVector{<:Count}

# Misc see https://alan-turing-institute.github.io/MLJ.jl/stable/adding_models_for_general_use/
load_path(::Type{JLBoostClassifier}) = "JLBoostMLJ.JLBoostClassifier"
load_path(::Type{<:JLBoostRegressor}) = "JLBoostMLJ.JLBoostRegressor"
load_path(::Type{<:JLBoostCount}) = "JLBoostMLJ.JLBoostCount"

package_name(::Type{<:JLBoostMLJModel}) = "JLBoostMLJ"
package_uuid(::Type{<:JLBoostMLJModel}) = "cb937e20-20f2-4cea-8a28-54eef8bab285"
package_url(::Type{<:JLBoostMLJModel}) = "https://github.com/xiaodaigh/JLBoostMLJ.jl"
is_pure_julia(::Type{<:JLBoostMLJModel}) = true
package_license(::Type{<:JLBoostMLJModel}) = "MIT"

docstring(::Type{JLBoostClassifier}) =
    "The JLBoost gradient boosting method, for use with "*
    "`Binary=Finite{2}` univariate targets"

docstring(::Type{JLBoostRegressor}) =
    "The JLBoost gradient boosting method, for use with "*
    "`Continuous` univariate targets"

docstring(::Type{JLBoostCount}) =
    "The JLBoost gradient boosting method, for use with "*
    "`Count` univariate targets, using a Poisson loss function"


text:
f(x, y) = cos.(x.^2 .+ transpose(y.^2))

xmin, xmax, nx = -1, 1, 16
ymin, ymax, ny = -1, 1, 16

# Evaluation points
nxe, nye = 64, 128
xe = collect(range(xmin, stop=xmax, length=nxe))
ye = collect(range(ymin, stop=ymax, length=nye))
fe = f(xe, ye)
dx = abs(xe[2] - xe[1])
dy = abs(ye[2] - ye[1])

ϵ = 0.1
γ = 3.0

hermite_x = Hermite(Chebyshev(xmax, xmin, nx), ϵ, γ)
hermite_y = Hermite(Chebyshev(ymax, ymin, ny), ϵ, γ)

xk = hermite_x.nodes
yk = hermite_y.nodes
fk = f(xk, yk)

@testset "Hermite 2D" begin

    @time s = interpolate( hermite_x, hermite_y, fk, xe, ye )

    max_error = maximum(abs.(s .- f(xe, ye)))
    l2_error = sqrt(trapz((s .- f(xe, ye)).^2, dx, dy))

    println(" Max error $max_error ")
    println(" L2  error $l2_error  ")

    @test l2_error < 1.0e-7

end

radial_x  = Radial(Chebyshev(xmax, xmin, nx), ϵ)
radial_y  = Radial(Chebyshev(ymax, ymin, ny), ϵ)

xk = radial_x.nodes
yk = radial_y.nodes
fk = f(xk, yk)

@testset "Radial 2D" begin

    s = interpolate( radial_x, radial_y, fk, xe, ye )

    max_error = maximum(abs.(s .- f(xe, ye)))
    l2_error = sqrt(trapz((s .- f(xe, ye)).^2, dx, dy))

    println(" Max error $max_error ")
    println(" L2  error $l2_error  ")

    @test true

end


text:
export dilation

"""
function dilation(mc,M::AbstractMesh,ssize=2)

dilated the binary-valued image mc by ssize cells.
The code works in 2D and 3D

"""
function dilation(mc::BitArray,M::AbstractMesh,ssize=2,boundLow=0.1,boundHigh=0.9)

    n   = M.n
    mc  = reshape(copy(mc),tuple(n...))

    A   = ones(1+2*ssize)
    A /= sum(vec(A))
    rId = Array{Float64}(mc)

    perm = zeros(Int64,M.dim)
    perm[1:M.dim-1] = collect(2:M.dim)
    perm[end]=1
    perm
    for k=1:M.dim
        szrId = size(rId)
        rId = reshape(rId,size(rId,1),:)
        for j=1:size(rId,2)
            rId[:,j] = conv(A,rId[:,j])[1+ssize:end-ssize]
        end
        rId = reshape(rId,szrId)
        rId = permutedims(rId,perm)
    end
    rId = (vec(rId).>=boundLow) .& (vec(rId) .<= boundHigh)
    return rId
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
using Base.Sort, Base.Order

refs(v::PooledArray) = v.refs
refs(v::AbstractArray) = v
refs(v::StructArray) = StructArray(map(refs, fieldarrays(v)))

function Base.permute!!(c::StructArray, p::AbstractVector{<:Integer})
    Base.invoke(Base.permute!!, Tuple{AbstractArray, AbstractVector{<:Integer}}, refs(c), p)
    return c
end

struct GroupPerm{V<:AbstractVector, P<:AbstractVector{<:Integer}, U<:AbstractUnitRange}
    vec::V
    perm::P
    within::U
end

GroupPerm(vec, perm=sortperm(vec)) = GroupPerm(vec, perm, axes(vec, 1))

Base.sortperm(g::GroupPerm) = g.perm
Base.parent(g::GroupPerm) = g.vec

function Base.iterate(g::GroupPerm, i = first(g.within))
    vec, perm = g.vec, g.perm
    l = last(g.within)
    i > l && return nothing
    @inbounds pi = perm[i]
    i1 = i+1
    @inbounds while i1 <= l && roweq(vec, pi, perm[i1])
        i1 += 1
    end
    return (i:(i1-1), i1)
end

Base.IteratorSize(::Type{<:GroupPerm}) = Base.SizeUnknown()

Base.eltype(::Type{<:GroupPerm}) = UnitRange{Int}

@inline function roweq(x::AbstractVector, i, j)
    r = refs(x)
    @inbounds eq = isequal(r[i], r[j])
    return eq
end

@generated function roweq(c::StructVector{D,C}, i, j) where {D,C}
    N = fieldcount(C)
    ex = :(roweq(getfield(fieldarrays(c),1), i, j))
    for n in 2:N
        ex = :(($ex) && (roweq(getfield(fieldarrays(c),$n), i, j)))
    end
    ex
end

function uniquesorted(keys, perm=sortperm(keys))
    (keys[perm[idxs[1]]] for idxs in GroupPerm(keys, perm))
end

function finduniquesorted(keys, perm=sortperm(keys))
    func = function (idxs)
        p_idxs = perm[idxs]
        return keys[p_idxs[1]] => p_idxs
    end
    (func(idxs) for idxs in GroupPerm(keys, perm))
end

function Base.sortperm(c::StructVector{T}) where {T<:Union{Tuple, NamedTuple}}
    cols = fieldarrays(c)
    x = cols[1]
    p = sortperm(x)
    if length(cols) > 1
        y = cols[2]
        refine_perm!(p, cols, 1, x, y, 1, length(x))
    end
    return p
end

Base.sort!(c::StructArray{<:Union{Tuple, NamedTuple}}) = (Base.permute!!(c, sortperm(c)); c)
Base.sort(c::StructArray{<:Union{Tuple, NamedTuple}}) = c[sortperm(c)]

# Given an ordering `p`, return a vector `v` such that `Perm(Forward, v)` is
# equivalent to `p`. Return `nothing` if such vector is not found.
forward_vec(p::Perm{ForwardOrdering}) = p.data
forward_vec(::Ordering) = nothing

# Methods from IndexedTables to refine sorting:
# # assuming x[p] is sorted, sort by remaining columns where x[p] is constant
function refine_perm!(p, cols, c, x, y′, lo, hi)
    temp = similar(p, 0)
    order = Perm(Forward, y′)
    y = something(forward_vec(order), y′)
    nc = length(cols)
    for idxs in GroupPerm(x, p, lo:hi)
        i, i1 = extrema(idxs)
        if i1 > i
            sort_sub_by!(p, i, i1, y, order, temp)
            if c < nc-1
                z = cols[c+2]
                refine_perm!(p, cols, c+1, y, z, i, i1)
            end
        end
    end
end

# sort the values in v[i0:i1] in place, by array `by`
Base.@noinline function sort_sub_by!(v, i0, i1, by, order, temp)
    empty!(temp)
    sort!(v, i0, i1, MergeSort, order, temp)
end

Base.@noinline function sort_sub_by!(v, i0, i1, by::AbstractVector{T}, order, temp) where T<:Integer
    min = max = by[v[i0]]
    @inbounds for i = i0+1:i1
        val = by[v[i]]
        if val < min
            min = val
        elseif val > max
            max = val
        end
    end
    rangelen = max-min+1
    n = i1-i0+1
    if rangelen <= n
        sort_int_range_sub_by!(v, i0-1, n, by, rangelen, min, temp)
    else
        empty!(temp)
        sort!(v, i0, i1, MergeSort, order, temp)
    end
    v
end

# in-place counting sort of x[ioffs+1:ioffs+n] by values in `by`
function sort_int_range_sub_by!(x, ioffs, n, by, rangelen, minval, temp)
    offs = 1 - minval

    where = fill(0, rangelen+1)
    where[1] = 1
    @inbounds for i = 1:n
        where[by[x[i+ioffs]] + offs + 1] += 1
    end
    cumsum!(where, where)

    length(temp) < n && resize!(temp, n)
    @inbounds for i = 1:n
        xi = x[i+ioffs]
        label = by[xi] + offs
        wl = where[label]
        temp[wl] = xi
        where[label] = wl+1
    end

    @inbounds for i = 1:n
        x[i+ioffs] = temp[i]
    end
    x
end


text:
isnanmissing(x) =  ismissing(x) || isnan(x)

function nantomissing!(d::AbstractDataFrame,x::Union{Integer,Symbol,String})
    if any(isnan,d[!,x])
        allowmissing!(d,x)
	    replace!(d[!,x], NaN => missing)
    end
end

function nantomissing!(d::AbstractDataFrame)
	for nm in names(d)
		if eltype(d[!,nm]) <: Float64
			nantomissing!(d,nm)
		end
	end
end

showallrows(x) = show(x; allrows=true)
showallrows(io,x) = show(io,x; allrows=true)


text:
# tests metadata

@generator MDGeneratesGen() begin # prefix MD (for Meta Data) to avoid type name clashes
    generates: ["a roll of a die", "a roll of a dice"]
    start() = choose(Int, 1, 6)
end

@generator MDBlockGen begin
begin
    generates: ["a month number"]
end
begin
    start() = choose(Int, 1, 12)
end
end

@testset "metadata" begin

	@testset "generator with 'generates' metadata" begin

		gn = MDGeneratesGen()

	    md = meta(gn)
	    @test md[:generates] == ["a roll of a die", "a roll of a dice"]

	    td = choose(gn)
	    @test 1 <= td <= 6

	end

	@testset "generator with 'generates' metadata inside a block" begin

		gn = MDBlockGen()

	    md = meta(gn)
	    @test md[:generates] == ["a month number"]

	end

end


text:
using Revise
using GF256
using Test

function check_inv()::Bool
    for i = 1:255
        if GF(1) != (GF(i) * inv(GF(i)))
            return false
        end
    end
    return true
end

a = GF(2)
b = GF(2)
c = GF(3)

@test a == GF(2)
@test (a == b) == true
@test (a == c) == false
@test bvec(c) == [1,1,0,0,0,0,0,0]
@test (a + c) == GF(1)
@test (c + c) == GF(0)
@test log(c) == 25
@test (a * a) == GF(4)
@test (GF(7) * GF(22)) == GF(98)
@test check_inv()
@test inv(GF(0)) == GF(0)
@test (c ^ 5) == GF(51)
@test (a ^ 300) == GF(193)
@test (a ^ -257) == GF(71)


text:
import CPLEX

function run_rmp_cplex(frmp)

    CPXENV = CPLEX.Env()

    # Set parameters
    CPLEX.set_param!(CPXENV, "CPXPARAM_ScreenOutput", 1)
    CPLEX.set_param!(CPXENV, "CPXPARAM_TimeLimit", 3600)
    CPLEX.set_param!(CPXENV, "CPXPARAM_Threads", 1)  # Single threads
    CPLEX.set_param!(CPXENV, "CPXPARAM_SolutionType", 2)  # No crossover
    CPLEX.set_param!(CPXENV, "CPXPARAM_LPMethod", 4)  # barrier
    CPLEX.set_param!(CPXENV, "CPXPARAM_Preprocessing_Presolve", 0)  # No presolve
    cplex = CPLEX.Model(CPXENV)

    # Read file
    CPLEX.read_model(cplex, frmp)

    # Solve instance
    t = @elapsed CPLEX.optimize!(cplex)

    println("Solve time: $t")
    println("Barrier iterations: ", CPLEX.c_api_getbaritcnt(cplex))
    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    run_rmp_cplex(ARGS[1])
end

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks