text:
function DiffEqBase.__solve(prob::Union{DiffEqBase.AbstractODEProblem,DiffEqBase.AbstractDAEProblem},
                            alg::Union{OrdinaryDiffEqAlgorithm,DAEAlgorithm}, args...;
                            kwargs...)
  integrator = DiffEqBase.__init(prob, alg, args...; kwargs...)
  solve!(integrator)
  integrator.sol
end

function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem,DiffEqBase.AbstractDAEProblem},
                           alg::Union{OrdinaryDiffEqAlgorithm,DAEAlgorithm},
                           timeseries_init = (),
                           ts_init = (),
                           ks_init = (),
                           recompile::Type{Val{recompile_flag}} = Val{true};
                           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 && !(typeof(alg) <: Union{DAEAlgorithm,FunctionMap}) && isempty(saveat),
                           calck = (callback !== nothing && callback != CallbackSet()) || (dense), # and no dense output
                           dt = alg isa FunctionMap && isempty(tstops) ? eltype(prob.tspan)(1) : eltype(prob.tspan)(0),
                           dtmin = nothing,
                           dtmax = eltype(prob.tspan)((prob.tspan[end]-prob.tspan[1])),
                           force_dtmin = false,
                           adaptive = isadaptive(alg),
                           gamma = gamma_default(alg),
                           abstol = nothing,
                           reltol = nothing,
                           qmin = qmin_default(alg),
                           qmax = qmax_default(alg),
                           qsteady_min = qsteady_min_default(alg),
                           qsteady_max = qsteady_max_default(alg),
                           qoldinit = isadaptive(alg) ? 1//10^4 : 0,
                           fullnormalize = true,
                           failfactor = 2,
                           beta1 = nothing,
                           beta2 = nothing,
                           maxiters = adaptive ? 1000000 : typemax(Int),
                           internalnorm = ODE_DEFAULT_NORM,
                           internalopnorm = LinearAlgebra.opnorm,
                           isoutofdomain = ODE_DEFAULT_ISOUTOFDOMAIN,
                           unstable_check = 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 = "ODE",
                           progress_message = ODE_DEFAULT_PROG_MESSAGE,
                           userdata = nothing,
                           allow_extrapolation = alg_extrapolates(alg),
                           initialize_integrator = true,
                           alias_u0 = false,
                           alias_du0 = false,
                           initializealg = DefaultInit(),
                           kwargs...) where recompile_flag

  if prob isa DiffEqBase.AbstractDAEProblem && alg isa OrdinaryDiffEqAlgorithm
    error("You cannot use an ODE Algorithm with a DAEProblem")
  end

  if prob isa DiffEqBase.AbstractODEProblem && alg isa DAEAlgorithm
    error("You cannot use an DAE Algorithm with a ODEProblem")
  end

  if typeof(prob.f)<:DynamicalODEFunction && typeof(prob.f.mass_matrix)<:Tuple
    if any(mm != I for mm in prob.f.mass_matrix)
      error("This solver is not able to use mass matrices.")
    end
  elseif !(typeof(prob)<:DiscreteProblem) &&
         !(typeof(prob)<:DiffEqBase.AbstractDAEProblem) &&
         !is_mass_matrix_alg(alg) &&
         prob.f.mass_matrix != I
    error("This solver is not able to use mass matrices.")
  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(LogLevel(-1),progress_name,_id=_id = :OrdinaryDiffEq,progress=0)

  tType = eltype(prob.tspan)
  tspan = prob.tspan
  tdir = sign(tspan[end]-tspan[1])

  t = tspan[1]

  if (((!(typeof(alg) <: OrdinaryDiffEqAdaptiveAlgorithm) && !(typeof(alg) <: OrdinaryDiffEqCompositeAlgorithm) && !(typeof(alg) <: DAEAlgorithm)) || !adaptive) && dt == tType(0) && isempty(tstops)) && !(typeof(alg) <: Union{FunctionMap,LinearExponential})
      error("Fixed timestep methods require a choice of dt or choosing the tstops")
  end

  isdae = alg isa DAEAlgorithm || (!(typeof(prob)<:DiscreteProblem) &&
                                     prob.f.mass_matrix != I &&
                                     !(typeof(prob.f.mass_matrix)<:Tuple) &&
                                     ArrayInterface.issingular(prob.f.mass_matrix))
  if alg isa CompositeAlgorithm && alg.choice_function isa AutoSwitch
    auto = alg.choice_function
    _alg = CompositeAlgorithm(alg.algs,
                             AutoSwitchCache(
                                             0,0,
                                             auto.nonstiffalg,
                                             auto.stiffalg,
                                             auto.stiffalgfirst,
                                             auto.maxstiffstep,
                                             auto.maxnonstiffstep,
                                             auto.nonstifftol,
                                             auto.stifftol,
                                             auto.dtfac,
                                             auto.stiffalgfirst,
                                             auto.switch_max
                                            ))
  else
    _alg = alg
  end
  f = prob.f
  p = prob.p

  # Get the control variables

  if alias_u0
    u = prob.u0
  else
    u = recursivecopy(prob.u0)
  end

  if _alg isa DAEAlgorithm
    if alias_du0
      du = prob.du0
    else
      du = recursivecopy(prob.u0)
    end
    duprev = recursivecopy(du)
  else
    du = nothing
    duprev = nothing
  end

  uType = typeof(u)
  uBottomEltype = recursive_bottom_eltype(u)
  uBottomEltypeNoUnits = recursive_unitless_bottom_eltype(u)

  uEltypeNoUnits = recursive_unitless_eltype(u)
  tTypeNoUnits   = typeof(one(tType))

  if typeof(_alg) <: FunctionMap
    abstol_internal = false
  elseif abstol === nothing
    if uBottomEltypeNoUnits == uBottomEltype
      abstol_internal = real(convert(uBottomEltype,oneunit(uBottomEltype)*1//10^6))
    else
      abstol_internal = real.(oneunit.(u).*1//10^6)
    end
  else
    abstol_internal = real.(abstol)
  end

  if typeof(_alg) <: FunctionMap
    reltol_internal = false
  elseif reltol === nothing
    if uBottomEltypeNoUnits == uBottomEltype
      reltol_internal = real(convert(uBottomEltype,oneunit(uBottomEltype)*1//10^3))
    else
      reltol_internal = real.(oneunit.(u).*1//10^3)
    end
  else
    reltol_internal = real.(reltol)
  end

  dtmax > zero(dtmax) && tdir < 0 && (dtmax *= tdir) # Allow positive dtmax, but auto-convert
  # dtmin is all abs => does not care about sign already.

  if !isdae && isinplace(prob) && typeof(u) <: AbstractArray && eltype(u) <: Number && uBottomEltypeNoUnits == uBottomEltype && tType == tTypeNoUnits # Could this be more efficient for other arrays?
    rate_prototype = recursivecopy(u)
  elseif prob isa DAEProblem
    rate_prototype = prob.du0
  else
    if (uBottomEltypeNoUnits == uBottomEltype && tType == tTypeNoUnits) || eltype(u) <: Enum
      rate_prototype = u
    else # has units!
      rate_prototype = u/oneunit(tType)
    end
  end
  rateType = typeof(rate_prototype) ## Can be different if united

  if isdae
    if uBottomEltype == uBottomEltypeNoUnits
      res_prototype = u
    else
      res_prototype = one(u)
    end
    resType = typeof(res_prototype)
  end

  tstops_internal = initialize_tstops(tType, tstops, d_discontinuities, tspan)
  saveat_internal = initialize_saveat(tType, saveat, tspan)
  d_discontinuities_internal = initialize_d_discontinuities(tType, d_discontinuities, tspan)

  callbacks_internal = CallbackSet(callback)

  max_len_cb = DiffEqBase.max_vector_callback_length(callbacks_internal)
  if max_len_cb isa VectorContinuousCallback
    uBottomEltypeReal = real(uBottomEltype)
    if isinplace(prob)
      callback_cache = DiffEqBase.CallbackCache(u,max_len_cb.len,uBottomEltypeReal,uBottomEltypeReal)
    else
      callback_cache = DiffEqBase.CallbackCache(max_len_cb.len,uBottomEltypeReal,uBottomEltypeReal)
    end
  else
    callback_cache = nothing
  end

  ### Algorithm-specific defaults ###
  if save_idxs === nothing
    ksEltype = Vector{rateType}
  else
    ks_prototype = rate_prototype[save_idxs]
    ksEltype = Vector{typeof(ks_prototype)}
  end

  # Have to convert incase passed in wrong.
  if save_idxs === nothing
    timeseries = timeseries_init === () ? uType[] : convert(Vector{uType},timeseries_init)
  else
    u_initial = u[save_idxs]
    timeseries = timeseries_init === () ? typeof(u_initial)[] :
                                          convert(Vector{uType},timeseries_init)
  end

  ts = ts_init === () ? tType[]    : convert(Vector{tType},ts_init)
  ks = ks_init === () ? ksEltype[] : convert(Vector{ksEltype},ks_init)
  alg_choice = typeof(_alg) <: CompositeAlgorithm ? Int[] : ()

  if !adaptive && save_everystep && tspan[2]-tspan[1] != Inf
    if dt == 0
      steps = length(tstops)
    else
      dtmin === nothing && (dtmin = DiffEqBase.prob2dtmin(prob; use_end_time=true))
      abs(dt) < dtmin && throw(ArgumentError("Supplied dt is smaller than dtmin"))
      steps = ceil(Int,internalnorm((tspan[2]-tspan[1])/dt,tspan[1]))
    end
    sizehint!(timeseries,steps+1)
    sizehint!(ts,steps+1)
    sizehint!(ks,steps+1)
  elseif save_everystep
    sizehint!(timeseries,50)
    sizehint!(ts,50)
    sizehint!(ks,50)
  elseif !isempty(saveat_internal)
    sizehint!(timeseries,length(saveat_internal)+1)
    sizehint!(ts,length(saveat_internal)+1)
    sizehint!(ks,length(saveat_internal)+1)
  else
    sizehint!(timeseries,2)
    sizehint!(ts,2)
    sizehint!(ks,2)
  end

  QT = if tTypeNoUnits <: Integer
    typeof(qmin)
  elseif prob isa DiscreteProblem
    # The QT fields are not used for DiscreteProblems
    constvalue(tTypeNoUnits)
  else
    typeof(internalnorm(u, t))
  end

  k = rateType[]

  if uses_uprev(_alg, adaptive) || calck
    uprev = recursivecopy(u)
  else
    # Some algorithms do not use `uprev` explicitly. In that case, we can save
    # some memory by aliasing `uprev = u`, e.g. for "2N" low storage methods.
    uprev = u
  end
  if allow_extrapolation
    uprev2 = recursivecopy(u)
  else
    uprev2 = uprev
  end

  if prob isa DAEProblem
    cache = alg_cache(_alg,du,u,res_prototype,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol_internal,p,calck,Val(isinplace(prob)))
  else
    cache = alg_cache(_alg,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol_internal,p,calck,Val(isinplace(prob)))
  end

  if typeof(_alg) <: OrdinaryDiffEqCompositeAlgorithm
    id = CompositeInterpolationData(f,timeseries,ts,ks,alg_choice,dense,cache)
    beta2 === nothing && ( beta2=_composite_beta2_default(_alg.algs, cache.current, QT) )
    beta1 === nothing && ( beta1=_composite_beta1_default(_alg.algs, cache.current, QT, beta2) )
  else
    id = InterpolationData(f,timeseries,ts,ks,dense,cache)
    beta2 === nothing && ( beta2=beta2_default(_alg) )
    beta1 === nothing && ( beta1=beta1_default(_alg,beta2) )
  end

  dtmin === nothing && (dtmin = DiffEqBase.prob2dtmin(prob; use_end_time=false))

  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 = DEOptions{typeof(abstol_internal),typeof(reltol_internal),QT,tType,
                   typeof(internalnorm),typeof(internalopnorm),typeof(save_end_user),
                   typeof(callbacks_internal),
                   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(failfactor),tType(dtmax),
                       tType(dtmin),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,
                       QT(beta1),QT(beta2),QT(qoldinit),dense,
                       save_on,save_start,save_end,save_end_user,
                       callbacks_internal,isoutofdomain,
                       unstable_check,verbose,
                       calck,force_dtmin,advance_to_tstop,stop_at_next_tstop)

  destats = DiffEqBase.DEStats(0)

  if typeof(_alg) <: OrdinaryDiffEqCompositeAlgorithm
    sol = DiffEqBase.build_solution(prob,_alg,ts,timeseries,
                      dense=dense,k=ks,interp=id,
                      alg_choice=alg_choice,
                      calculate_error = false, destats=destats)
  else
    sol = DiffEqBase.build_solution(prob,_alg,ts,timeseries,
                      dense=dense,k=ks,interp=id,
                      calculate_error = false, destats=destats)
  end

  if recompile_flag == true
    FType = typeof(f)
    SolType = typeof(sol)
    cacheType = typeof(cache)
  else
    FType = Function
    if _alg isa OrdinaryDiffEqAlgorithm
      SolType = DiffEqBase.AbstractODESolution
      cacheType =  OrdinaryDiffEqCache
    else
      SolType = DiffEqBase.AbstractDAESolution
      cacheType =  DAECache
    end
  end

  # rate/state = (state/time)/state = 1/t units, internalnorm drops units
  # we don't want to differentiate through eigenvalue estimation
  eigen_est = inv(one(tType))
  tprev = t
  dtcache = tType(dt)
  dtpropose = tType(dt)
  iter = 0
  kshortsize = 0
  reeval_fsal = false
  u_modified = false
  EEst = QT(1)
  just_hit_tstop = false
  isout = false
  accept_step = false
  force_stepfail = false
  last_stepfail = false
  do_error_check = true
  event_last_time = 0
  vector_event_last_time = 1
  last_event_error = typeof(_alg) <: FunctionMap ? false : zero(uBottomEltypeNoUnits)
  dtchangeable = isdtchangeable(_alg)
  q11 = QT(1)
  success_iter = 0
  erracc = QT(1)
  dtacc = tType(1)
  reinitiailize = true
  saveiter = 0 # Starts at 0 so first save is at 1
  saveiter_dense = 0

  integrator = ODEIntegrator{typeof(_alg),isinplace(prob),uType,typeof(du),
                             tType,typeof(p),
                             typeof(eigen_est),
                             QT,typeof(tdir),typeof(k),SolType,
                             FType,cacheType,
                             typeof(opts),fsal_typeof(_alg,rate_prototype),
                             typeof(last_event_error),typeof(callback_cache),
                             typeof(initializealg)}(
                             sol,u,du,k,t,tType(dt),f,p,uprev,uprev2,duprev,tprev,
                             _alg,dtcache,dtchangeable,
                             dtpropose,tdir,eigen_est,EEst,QT(qoldinit),q11,
                             erracc,dtacc,success_iter,
                             iter,saveiter,saveiter_dense,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,reinitiailize,isdae,
                             opts,destats,initializealg)
  if initialize_integrator
    if isdae
      DiffEqBase.initialize_dae!(integrator)
    end

    if save_start
      integrator.saveiter += 1 # Starts at 1 so first save is at 2
      integrator.saveiter_dense += 1
      copyat_or_push!(ts,1,t)
      if save_idxs === nothing
        copyat_or_push!(timeseries,1,integrator.u)
        copyat_or_push!(ks,1,[rate_prototype])
      else
        copyat_or_push!(timeseries,1,u_initial,Val{false})
        copyat_or_push!(ks,1,[ks_prototype])
      end
    else
      saveiter = 0 # Starts at 0 so first save is at 1
      saveiter_dense = 0
    end

    initialize_callbacks!(integrator, initialize_save)
    initialize!(integrator,integrator.cache)

    if typeof(_alg) <: OrdinaryDiffEqCompositeAlgorithm && save_start
      # Loop to get all of the extra possible saves in callback initialization
      for i in 1:integrator.saveiter
        copyat_or_push!(alg_choice,i,integrator.cache.current)
      end
    end
  end

  handle_dt!(integrator)

  integrator
end

function DiffEqBase.solve!(integrator::ODEIntegrator)
  @inbounds while !isempty(integrator.opts.tstops)
    while integrator.tdir * integrator.t < first(integrator.opts.tstops)
      loopheader!(integrator)
      if integrator.do_error_check && check_error!(integrator) != :Success
        return integrator.sol
      end
      perform_step!(integrator,integrator.cache)
      loopfooter!(integrator)
      if isempty(integrator.opts.tstops)
        break
      end
    end
    handle_tstop!(integrator)
  end
  postamble!(integrator)

  f = integrator.sol.prob.f

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

# Helpers

function handle_dt!(integrator)
  if iszero(integrator.dt) && integrator.opts.adaptive
    auto_dt_reset!(integrator)
    if sign(integrator.dt)!=integrator.tdir && !iszero(integrator.dt) && !isnan(integrator.dt)
      error("Automatic dt setting has the wrong sign. Exiting. Please report this error.")
    end
    if isnan(integrator.dt)
      if integrator.opts.verbose
        @warn("Automatic dt set the starting dt as NaN, causing instability.")
      end
    end
  elseif integrator.opts.adaptive && integrator.dt > zero(integrator.dt) && integrator.tdir < 0
    integrator.dt *= integrator.tdir # Allow positive dt, but auto-convert
  end
end

# time stops
function initialize_tstops(::Type{T}, tstops, d_discontinuities, tspan) where T
  tstops_internal = BinaryMinHeap{T}()

  t0, tf = tspan
  tdir = sign(tf - t0)
  tdir_t0 = tdir * t0
  tdir_tf = tdir * tf

  if isempty(d_discontinuities) && isempty(tstops) # TODO: Specialize more
    push!(tstops_internal, tdir_tf)
  else
    for t in tstops
      tdir_t = tdir * t
      tdir_t0 < tdir_t ≤ tdir_tf && push!(tstops_internal, tdir_t)
    end

    for t in d_discontinuities
      tdir_t = tdir * t
      tdir_t0 < tdir_t ≤ tdir_tf && push!(tstops_internal, tdir_t)
    end

    push!(tstops_internal, tdir_tf)
  end

  return tstops_internal
end

# saving time points
function initialize_saveat(::Type{T}, saveat, tspan) where T
  saveat_internal = BinaryMinHeap{T}()

  t0, tf = tspan
  tdir = sign(tf - t0)
  tdir_t0 = tdir * t0
  tdir_tf = tdir * tf

  if typeof(saveat) <: Number
    directional_saveat = tdir * abs(saveat)
    for t in (t0 + directional_saveat):directional_saveat:tf
      push!(saveat_internal, tdir * t)
    end
  elseif !isempty(saveat)
    for t in saveat
      tdir_t = tdir * t
      tdir_t0 < tdir_t ≤ tdir_tf && push!(saveat_internal, tdir_t)
    end
  end

  return saveat_internal
end

# discontinuities
function initialize_d_discontinuities(::Type{T}, d_discontinuities, tspan) where T
  d_discontinuities_internal = BinaryMinHeap{T}()
  sizehint!(d_discontinuities_internal, length(d_discontinuities))

  t0, tf = tspan
  tdir = sign(tf - t0)

  for t in d_discontinuities
    push!(d_discontinuities_internal, tdir * t)
  end

  return d_discontinuities_internal
end

function initialize_callbacks!(integrator, initialize_save = true)
  t = integrator.t
  u = integrator.u
  callbacks = integrator.opts.callback
  integrator.u_modified = true

  u_modified = initialize!(callbacks,u,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(integrator.sol.prob)
      recursivecopy!(integrator.uprev,integrator.u)
    else
      integrator.uprev = integrator.u
    end

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

    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


text:
###############################################################################
# Callbacks
###############################################################################

function curl_write_cb(buff::Ptr{Cchar}, s::Csize_t, n::Csize_t, p::Ptr{Cvoid})::Csize_t
  ctxt = unsafe_pointer_to_objref(p)
  nbytes = s * n
  write(ctxt.resp.body, unsafe_string(buff, nbytes))
  ctxt.bytes_recd = ctxt.bytes_recd + nbytes

  nbytes
end

function curl_read_cb(out::Ptr{Cchar}, s::Csize_t, n::Csize_t, p::Ptr{Cvoid})::Csize_t
  ctxt = unsafe_pointer_to_objref(p)
  bavail = s * n
  breq = ctxt.rd.sz - ctxt.rd.offset
  b2copy = bavail > breq ? breq : bavail

  if ctxt.rd.typ == :buffer
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
          out, convert(Ptr{UInt8}, ctxt.rd.str) + ctxt.rd.offset, b2copy)
  elseif ctxt.rd.typ == :io
    b_read = read(ctxt.rd.src, b2copy)
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt), out, b_read, b2copy)
  end
  ctxt.rd.offset = ctxt.rd.offset + b2copy

  b2copy
end

function curl_multi_timer_cb(curlm::Ptr{Cvoid}, timeout_ms::Clong, p::Ptr{Cvoid})::Cint
  muctxt = unsafe_pointer_to_objref(p)
  muctxt.timeout = timeout_ms / 1000.0

  @info "Requested timeout value : " * string(muctxt.timeout)

  0
end

null_cb(curl) = nothing

ext: jl


text:
using Dates
using Remarkable.Articles
using Remarkable.Tags

@testset "Live Performance" begin
    tagmatrix = TagMatrix()

    metadata = ArticleMetadata("foo",
                               "The Astronomical Canon",
                               ["Hypatia of Alexandria"],
                               DateTime(350, 1, 1, 12, 0, 0))
    object = LivePerformance(metadata, "Alexandria")

    @test location(object) == "Alexandria"
    @test authors(object) == ["Hypatia of Alexandria"]
    @test title(object) == "The Astronomical Canon"

    tag!(object, tagmatrix, ["math"])

    @test tags(object) == [tagobject(tagmatrix, "math")]
    @test collect(tags(tagmatrix)) == tags(object)
    @test Tags.popularity(tagmatrix, "math") == 1
end


text:
mutable struct Job
    entry_time::Float64
    exit_time::Float64
end

Job(exit_time::Float64) = Job(exit_time, exit_time)

""" Returns the duration of time the job has spent in the system """
duration(job::Job)::Float64 = job.exit_time - job.entry_time

isless(j1::Job, j2::Job) = j1.exit_time < j2.exit_time

function show(io::IO, job::Job)
    print(io, "Job(E: $(job.entry_time), L: $(job.exit_time))")
end


text:
using MomentClosure
using MomentClosure: define_μ, define_M, isvar
using Test
using Catalyst

@parameters t, c₁, c₂, c₃, c₄, Ω
@variables X(t), Y(t)

S_mat = [ 1 -1  1 -1;
         -1  1  0  0]
# NOTE: only holds if combinatoric_ratelaw=false
a = [c₁*X*Y*(X-1)/Ω^2, c₂*X, c₃*Ω, c₄*X]

# checking if ReactionSystem and ReactionSystemMod are consistent

rn1 = @reaction_network begin
    (c₁/Ω^2), 2X + Y → 3X
    (c₂), X → Y
    (Ω*c₃, c₄), 0 ↔ X
end c₁ c₂ c₃ c₄ Ω

rn2 = ReactionSystemMod(t, [X, Y], [c₁, c₂, c₃, c₄, Ω], a, S_mat)

@test isequal(species(rn1), species(rn2))
ps1 = reactionparams(rn1); ps2 = reactionparams(rn2)
@test length(ps1) == length(ps2) == 5
@test all(isvar.(ps1, Ref(ps2)))
@test all(isvar.(ps2, Ref(collect(keys(paramsmap(rn2))))))
@test isequal(speciesmap(rn1), speciesmap(rn2))
@test numspecies(rn1) == numspecies(rn2) == 2
@test numreactions(rn1) == numreactions(rn2) == 4
@test netstoichmat(rn1) == netstoichmat(rn2) == S_mat
@test isequal(propensities(rn1, combinatoric_ratelaw=false), propensities(rn2))


text:
__precompile__(true)
module BernsteinExpansions

using RecursiveArrayTools

export ImplicitForm,
       univariate,
       multivariate

include("implicit_form.jl")
include("tensorial_form.jl")

end # module


text:
using Hydra, DataFrames,Distances
using Test

@testset "Hydra.jl" begin
    # Write your tests here.
    x = rand(1:4,100)
y = randn(100,5) .+ (x .* (x .==3))

preds = [["a","b","c","d"][i] for i in x]
df = DataFrame([preds],[:x])
Hyde = hydra2(df,y,BrayCurtis,@formula(1~x))
Hyde = hydra2(df,y,BrayCurtis,@formula(1~x),pairs =true)
permute(Hyde,df.x,1000)

end

max_stars_repo_name: EvoArt/Hydra

max_issues_repo_name: EvoArt/Hydra

max_forks_repo_name: EvoArt/Hydra


text:
# Autogenerated wrapper script for basiclu_jll for x86_64-w64-mingw32
export libbasiclu

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("basiclu")
JLLWrappers.@declare_library_product(libbasiclu, "libbasiclu.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libbasiclu,
        "bin\\libbasiclu.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#################### Chains ####################

#################### Constructors ####################

function Chains(iters::Integer, params::Integer;
               start::Integer=1, thin::Integer=1, chains::Integer=1,
               names::Vector{T}=AbstractString[]) where {T<:AbstractString}
  value = Array{Float64}(undef, length(start:thin:iters), params, chains)
  fill!(value, NaN)
  Chains(value, start=start, thin=thin, names=names)
end

function Chains(value::Array{T, 3};
               start::Integer=1, thin::Integer=1,
               names::Vector{U}=AbstractString[], chains::Vector{V}=Int[]) where {T<:Real, U<:AbstractString, V<:Integer}
  n, p, m = size(value)

  if isempty(names)
    names = map(i -> "Param$i", 1:p)
  elseif length(names) != p
    throw(DimensionMismatch("size(value, 2) and names length differ"))
  end

  if isempty(chains)
    chains = collect(1:m)
  elseif length(chains) != m
    throw(DimensionMismatch("size(value, 3) and chains length differ"))
  end

  v = convert(Array{Float64, 3}, value)
  Chains(v, range(start, step=thin, length=n), AbstractString[names...], Int[chains...])
end

function Chains(value::Matrix{T};
               start::Integer=1, thin::Integer=1,
               names::Vector{U}=AbstractString[], chains::Integer=1) where {T<:Real, U<:AbstractString}
  Chains(reshape(value, size(value, 1), size(value, 2), 1), start=start,
         thin=thin, names=names, chains=Int[chains])
end

function Chains(value::Vector{T};
               start::Integer=1, thin::Integer=1,
               names::AbstractString="Param1", chains::Integer=1) where {T<:Real}
  Chains(reshape(value, length(value), 1, 1), start=start, thin=thin,
         names=AbstractString[names], chains=Int[chains])
end

#################### Indexing ####################

function Base.getindex(c::Chains, window, names, chains)
  inds1 = window2inds(c, window)
  inds2 = names2inds(c, names)
  Chains(c.value[inds1, inds2, chains],
         start = first(c) + (first(inds1) - 1) * step(c),
         thin = step(inds1) * step(c), names = c.names[inds2],
         chains = c.chains[chains])
end

Base.lastindex(c::AbstractChains, i) = size(c, i)

function Base.setindex!(c::AbstractChains, value, iters, names, chains)
  setindex!(c.value, value, iters2inds(c, iters), names2inds(c, names), chains)
end

macro mapiters(iters, c)
  return esc(quote
    ($iters .- first($c)) ./ step($c) .+ 1.0
  end)
end

window2inds(c::AbstractChains, window) =
  throw(ArgumentError("$(typeof(window)) iteration indexing is unsupported"))
window2inds(c::AbstractChains, ::Colon) = window2inds(c, 1:size(c, 1))
window2inds(c::AbstractChains, window::Int) = window2inds(c, window:window)
window2inds(c::AbstractChains, window::AbstractRange) = begin
  range = @mapiters(window, c)
  a = max(ceil(Int, first(range)), 1)
  b = step(window)
  c = min(floor(Int, last(range)), size(c.value, 1))
  a:b:c
end

iters2inds(c::AbstractChains, iters) = iters
iters2inds(c::AbstractChains, ::Colon) = 1:size(c.value, 1)
iters2inds(c::AbstractChains, iters::AbstractRange) =
  convert(StepRange{Int, Int}, @mapiters(iters, c))
iters2inds(c::AbstractChains, iter::Real) = Int(@mapiters(iter, c))
iters2inds(c::AbstractChains, iters::Vector{T}) where {T<:Real} =
  Int[@mapiters(i, c) for i in iters]

names2inds(c::AbstractChains, names) = names
names2inds(c::AbstractChains, ::Colon) = 1:size(c.value, 2)
names2inds(c::AbstractChains, name::Real) = [name]
names2inds(c::AbstractChains, name::AbstractString) = names2inds(c, [name])
names2inds(c::AbstractChains, names::Vector{T}) where {T<:AbstractString} =
  indexin(names, c.names)

#################### Concatenation ####################

function Base.cat(c1::AbstractChains, args::AbstractChains...; dims::Integer)
  dims == 1 ? cat1(c1, args...) :
  dims == 2 ? cat2(c1, args...) :
  dims == 3 ? cat3(c1, args...) :
    throw(ArgumentError("cannot concatenate along dimension $dim"))
end

function cat1(c1::AbstractChains, args::AbstractChains...)
  range = c1.range
  for c in args
    last(range) + step(range) == first(c) ||
      throw(ArgumentError("noncontiguous chain iterations"))
    step(range) == step(c) ||
      throw(ArgumentError("chain thinning differs"))
    range = first(range):step(range):last(c)
  end

  names = c1.names
  all(c -> c.names == names, args) ||
    throw(ArgumentError("chain names differ"))

  chains = c1.chains
  all(c -> c.chains == chains, args) ||
    throw(ArgumentError("sets of chains differ"))

  value = cat(c1.value, map(c -> c.value, args)..., dims=1)
  Chains(value, start=first(range), thin=step(range), names=names,
         chains=chains)
end

function cat2(c1::AbstractChains, args::AbstractChains...)
  range = c1.range
  all(c -> c.range == range, args) ||
    throw(ArgumentError("chain ranges differ"))

  names = c1.names
  n = length(names)
  for c in args
    names = union(names, c.names)
    n += length(c.names)
    n == length(names) ||
      throw(ArgumentError("non-unique chain names"))
  end

  chains = c1.chains
  all(c -> c.chains == chains, args) ||
    throw(ArgumentError("sets of chains differ"))

  value = cat(c1.value, map(c -> c.value, args)..., dims=2)
  Chains(value, start=first(range), thin=step(range), names=names,
         chains=chains)
end

function cat3(c1::AbstractChains, args::AbstractChains...)
  range = c1.range
  all(c -> c.range == range, args) ||
    throw(ArgumentError("chain ranges differ"))

  names = c1.names
  all(c -> c.names == names, args) ||
    throw(ArgumentError("chain names differ"))

  value = cat(c1.value, map(c -> c.value, args)..., dims=3)
  Chains(value, start=first(range), thin=step(range), names=names)
end

Base.hcat(c1::AbstractChains, args::AbstractChains...) = cat(c1, args..., dims=2)

Base.vcat(c1::AbstractChains, args::AbstractChains...) = cat(c1, args..., dims=1)

#################### Base Methods ####################

function Base.keys(c::AbstractChains)
  c.names
end

function Base.show(io::IO, c::AbstractChains)
  print(io, "Object of type \"$(summary(c))\"\n\n")
  println(io, header(c))
  show(io, c.value)
end

function Base.size(c::AbstractChains)
  dim = size(c.value)
  last(c), dim[2], dim[3]
end

function Base.size(c::AbstractChains, ind)
  size(c)[ind]
end

Base.first(c::AbstractChains) = first(c.range)
Base.step(c::AbstractChains) = step(c.range)
Base.last(c::AbstractChains) = last(c.range)

#################### Auxilliary Functions ####################

function combine(c::AbstractChains)
  n, p, m = size(c.value)
  value = Array{Float64}(undef, n * m, p)
  for j in 1:p
    idx = 1
    for i in 1:n, k in 1:m
      value[idx, j] = c.value[i, j, k]
      idx += 1
    end
  end
  value
end

function header(c::AbstractChains)
  string(
    "Iterations = $(first(c)):$(last(c))\n",
    "Thinning interval = $(step(c))\n",
    "Chains = $(join(map(string, c.chains), ","))\n",
    "Samples per chain = $(length(c.range))\n"
  )
end

function indiscretesupport(c::AbstractChains,
                           bounds::Tuple{Real, Real}=(0, Inf))
  nrows, nvars, nchains = size(c.value)
  result = Array{Bool}(undef, nvars * (nrows > 0))
  for i in 1:nvars
    result[i] = true
    for j in 1:nrows, k in 1:nchains
      x = c.value[j, i, k]
      if !isinteger(x) || x < bounds[1] || x > bounds[2]
        result[i] = false
        break
      end
    end
  end
  result
end

function link(c::AbstractChains)
  cc = copy(c.value)
  for j in 1:length(c.names)
    x = cc[:, j, :]
    if minimum(x) > 0.0
      cc[:, j, :] = maximum(x) < 1.0 ? logit.(x) : log.(x)
    end
  end
  cc
end


text:
function get_n_chains(model::T) where {T <: CmdStanModels}
  model.n_chains[1]
end

function set_n_chains(model::T, n_chains) where {T <: CmdStanModels}
  model.n_chains[1] = n_chains
end


text:
### FlatField types
# these are all just appropriately-parameterized BaseFields
# note: the seemingly-redundant <:AbstractArray{T} in the argument
# (which is enforced in BaseField anyway) is there to help prevent
# method ambiguities

# spin-0
const FlatMap{        M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{Map,        M, T, A}
const FlatFourier{    M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{Fourier,    M, T, A}
# spin-2
const FlatQUMap{      M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{QUMap,      M, T, A}
const FlatQUFourier{  M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{QUFourier,  M, T, A}
const FlatEBMap{      M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{EBMap,      M, T, A}
const FlatEBFourier{  M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{EBFourier,  M, T, A}
# spin-(0,2)
const FlatIQUMap{     M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{IQUMap,     M, T, A}
const FlatIQUFourier{ M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{IQUFourier, M, T, A}
const FlatIEBMap{     M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{IEBMap,     M, T, A}
const FlatIEBFourier{ M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{IEBFourier, M, T, A}

## FlatField unions
# spin-0
const FlatS0{         B<:Union{Map,Fourier},                         M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{B, M, T, A}
# spin-2
const FlatQU{         B<:Union{QUMap,QUFourier},                     M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{B, M, T, A}
const FlatEB{         B<:Union{EBMap,EBFourier},                     M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{B, M, T, A}
const FlatS2Map{      B<:Union{QUMap,EBMap},                         M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{B, M, T, A}
const FlatS2Fourier{  B<:Union{QUFourier,QUFourier},                 M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{B, M, T, A}
const FlatS2{         B<:Union{QUMap,QUFourier,EBMap,EBFourier},     M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{B, M, T, A}
# spin-(0,2)
const FlatIQU{        B<:Union{IQUMap,IQUFourier},                   M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{B, M, T, A}
const FlatIEB{        B<:Union{IEBMap,IEBFourier},                   M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{B, M, T, A}
const FlatS02Map{     B<:Union{IQUMap,IEBMap},                       M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{B, M, T, A}
const FlatS02Fourier{ B<:Union{IQUFourier,IQUFourier},               M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{B, M, T, A}
const FlatS02{        B<:Union{IQUMap,IQUFourier,IEBMap,IEBFourier}, M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{B, M, T, A}
# any flat projection
const FlatField{      B,                                             M<:FlatProj, T, A<:AbstractArray{T} } = BaseField{B, M, T, A}

# A Flat TEB covariance of the form:
#
#    [ΣTT ΣTE  ⋅
#     ΣTE ΣEE  ⋅
#      ⋅   ⋅  ΣBB]
#
# We store the 2x2 block as a 2x2 SMatrix, ΣTE, so that we can easily call sqrt/inv on
# it, and the ΣBB block separately as ΣB.
struct FlatIEBCov{T,F} <: ImplicitOp{T}
    ΣTE :: SMatrix{2,2,Diagonal{T,F},4}
    ΣB :: Diagonal{T,F}
end

### basis-like definitions
LenseBasis(::Type{<:FlatS0})    = Map
LenseBasis(::Type{<:FlatS2})    = QUMap
LenseBasis(::Type{<:FlatS02})   = IQUMap
DerivBasis(::Type{<:FlatS0})    = Fourier
DerivBasis(::Type{<:FlatS2})    = QUFourier
DerivBasis(::Type{<:FlatS02})   = IQUFourier
HarmonicBasis(::Type{<:FlatS0}) = Fourier
HarmonicBasis(::Type{<:FlatQU}) = QUFourier
HarmonicBasis(::Type{<:FlatEB}) = EBFourier

# useful for enumerating some cases below
_sub_components = [
    (:FlatMap,        ("Ix"=>:, "I"=>:)),
    (:FlatFourier,    ("Il"=>:, "I"=>:)),
    (:FlatQUMap,      ("Qx"=>1, "Ux"=>2, "Q" =>1, "U"=>2, "P"=>:)),
    (:FlatQUFourier,  ("Ql"=>1, "Ul"=>2, "Q" =>1, "U"=>2, "P"=>:)),
    (:FlatEBMap,      ("Ex"=>1, "Bx"=>2, "E" =>1, "B"=>2, "P"=>:)),
    (:FlatEBFourier,  ("El"=>1, "Bl"=>2, "E" =>1, "B"=>2, "P"=>:)),
    (:FlatIQUMap,     ("Ix"=>1, "Qx"=>2, "Ux"=>3, "I"=>1, "Q"=>2, "U"=>3, "P"=>2:3, "IP"=>:)),
    (:FlatIQUFourier, ("Il"=>1, "Ql"=>2, "Ul"=>3, "I"=>1, "Q"=>2, "U"=>3, "P"=>2:3, "IP"=>:)),
    (:FlatIEBMap,     ("Ix"=>1, "Ex"=>2, "Bx"=>3, "I"=>1, "E"=>2, "B"=>3, "P"=>2:3, "IP"=>:)),
    (:FlatIEBFourier, ("Il"=>1, "El"=>2, "Bl"=>3, "I"=>1, "E"=>2, "B"=>3, "P"=>2:3, "IP"=>:)),
]

### constructors
# the default constructor uses ProjLambert, but the code in this file
# should be agnostic to the projection type

_reshape_batch(arr::AbstractArray{T,3}) where {T} = reshape(arr, size(arr,1), size(arr,2), 1, size(arr,3))
_reshape_batch(arr) = arr

## constructing from arrays
# spin-0
function FlatMap(Ix::A; kwargs...) where {T, A<:AbstractArray{T}}
    FlatMap(_reshape_batch(Ix), ProjLambert(;Ny=size(Ix,1), Nx=size(Ix,2), T, storage=basetype(A), kwargs...))
end
function FlatFourier(Il::A; Ny, kwargs...) where {T, A<:AbstractArray{T}}
    FlatFourier(_reshape_batch(Il), ProjLambert(;Ny, Nx=size(Il,2), T, storage=basetype(A), kwargs...))
end
# spin-2
function FlatField{B}(X::A, Y::A, metadata)  where {T, A<:AbstractArray{T}, B<:Basis2Prod{<:Union{𝐐𝐔,𝐄𝐁}}}
    FlatField{B}(cat(_reshape_batch(X), _reshape_batch(Y), dims=Val(3)), metadata)
end
function FlatField{B}(X::A, Y::A; kwargs...) where {T, A<:AbstractArray{T}, B<:Basis2Prod{<:Union{𝐐𝐔,𝐄𝐁}}}
    FlatField{B}(X, Y, ProjLambert(;Ny=size(X,1), Nx=size(X,2), T, storage=basetype(A), kwargs...))
end
# spin-(0,2)
function FlatField{B}(X::A, Y::A, Z::A, metadata) where {T, A<:AbstractArray{T}, B<:Basis3Prod{𝐈,<:Union{𝐐𝐔,𝐄𝐁}}}
    FlatField{B}(cat(_reshape_batch(X), _reshape_batch(Y), _reshape_batch(Z), dims=Val(3)), metadata)
end
function FlatField{B}(X::A, Y::A, Z::A; kwargs...) where {T, A<:AbstractArray{T}, B<:Basis3Prod{𝐈,<:Union{𝐐𝐔,𝐄𝐁}}}
    FlatField{B}(X, Y, Z, ProjLambert(;Ny=size(X,1), Nx=size(X,2), T, storage=basetype(A), kwargs...))
end
## constructing from other fields
function FlatField{B}(X::FlatField{B₀}, Y::FlatField{B₀}) where {B₀<:Union{Map,Fourier}, B<:Basis2Prod{<:Union{𝐐𝐔,𝐄𝐁},B₀}}
    FlatField{B}(cat(X.arr, Y.arr, dims=Val(3)), get_metadata_strict(X, Y))
end
function FlatField{B}(X::FlatField{B₀}, Y::FlatField{Basis2Prod{Pol,B₀}}) where {B₀<:Union{Map,Fourier}, Pol<:Union{𝐐𝐔,𝐄𝐁}, B<:Basis3Prod{𝐈,Pol,B₀}}
    FlatField{B}(cat(X.arr, Y.arr, dims=Val(3)), get_metadata_strict(X, Y))
end
function FlatField{B}(X::FlatField{B₀}, Y::FlatField{B₀}, Z::FlatField{B₀}) where {B₀<:Union{Map,Fourier}, B<:Basis3Prod{𝐈,<:Union{𝐐𝐔,𝐄𝐁},B₀}}
    FlatField{B}(cat(X.arr, Y.arr, Z.arr, dims=Val(3)), get_metadata_strict(X, Y, Z))
end

## docstrings
for (F, props) in _sub_components

    B        = basis(@eval($F))
    origin   = B <: BasisProd ? "some `AbstractArrays`" : "an `AbstractArray`"
    dims     = B <: SpatialBasis{Map} ? "(Ny,Nx[,Nbatch])" : "(Ny÷2+1,Nx[,Nbatch])"
    F₀       = B <: SpatialBasis{Map} ? "FlatMap" : "FlatFourier"
    arr_args = join(["$X::AbstractArray" for X in first.(props) if endswith(X,r"[xl]")], ", ")
    f_args   = join(["$(X[1])::$F₀" for X in first.(props) if endswith(X,r"[xl]")], ", ")
    Ny_req   = B <: SpatialBasis{Fourier} ? "`Ny` must be given as keyword argument. " : ""

    doc = """
        $F($arr_args; $(B <: SpatialBasis{Fourier} ? "Ny, " : "")[θpix=1])
        $F($arr_args, proj::FieldMetadata)

    Construct a `$F` object from $origin. The array dimensions should
    be `$dims`, where `Ny` and `Nx` are the number of pixels in the
    map in the y/x direction, and `Nbatch` is an optional batch
    dimension. $Ny_req `θpix` is the angular resolution in arcmin
    (default=1). A second positional argument `proj` can be used
    instead of keyword arguments to specify the projection/metadata
    (currently the only builtin projection is [`ProjLambert`](@ref),
    which is the default if keyword arguments are used).
    """

    if B <: BasisProd
        doc *= """

            $F($f_args)

        Construct a `$F` object from other `$F₀` objects.
        Projection/metadata and batch size is required to be the same
        for the arguments.
        """
    end

    @eval @doc $doc $F
end

### array interface
# most is inherited from BaseField. the main thing we have specify
# here has to do with which dimension is the "batch" dimension
# (dimension 4), since that is not assumed in BaseField
similar(f::FlatField{B}, Nbatch::Int) where {B} = FlatField{B}(similar(f.arr, size(f.arr,1), size(f.arr,2), size(f.arr,3), Nbatch), f.metadata)
batch_axes(f::FlatField{B,M,T,A}) where {B,M,T,A<:AbstractArray{T,4}} = (f.Nbatch,)
batch_axes(f::FlatField{B,M,T,A}) where {B,M,T,A<:AbstractArray{T}} = ()
nonbatch_dims(f::FlatField) = ntuple(identity,min(3,ndims(f.arr)))
require_unbatched(f::FlatField) = (f.Nbatch==1) || error("This function not implemented for batched fields.")

### properties
# generic
getproperty(f::FlatField, ::Val{:Nbatch}) = size(getfield(f,:arr), 4)
getproperty(f::FlatField, ::Val{:Npol})   = size(getfield(f,:arr), 3)
getproperty(f::FlatField, ::Val{:T})      = eltype(f)
getproperty(f::FlatField, ::Val{:proj})   = getfield(f, :metadata)
# sub-components
for (F, props) in _sub_components
    for (k,I) in props
        body = if k[end] in "xl"
            I==(:) ? :(getfield(f,:arr)) : :(view(getfield(f,:arr), :, :, $I, ntuple(_->:,max(0,N-3))...))
        else
            I==(:) ? :f : :($(FlatField{k=="P" ? Basis2Prod{basis(@eval($F)).parameters[end-1:end]...} : basis(@eval($F)).parameters[end]})(_reshape_batch(view(getfield(f,:arr), :, :, $I, ntuple(_->:,max(0,N-3))...)), f.metadata))
        end
        @eval getproperty(f::$F{M,T,A}, ::Val{$(QuoteNode(Symbol(k)))}) where {M<:FlatProj,T,N,A<:AbstractArray{T,N}} = $body
    end
end

### indices
function getindex(f::FlatS0, k::Symbol; full_plane=false)
    maybe_unfold = full_plane ? x->unfold(x,fieldinfo(f).Ny) : identity
    @match k begin
        :I  => f
        :Ix => Map(f).Ix
        :Il => maybe_unfold(Fourier(f).Il)
        _   => throw(ArgumentError("Invalid FlatS0 index: $k"))
    end
end
function getindex(f::FlatS2{Basis2Prod{B₁,B₂}}, k::Symbol; full_plane=false) where {B₁,B₂}
    maybe_unfold = (full_plane && k in [:El,:Bl,:Ql,:Ul]) ? x->unfold(x,fieldinfo(f).Ny) : identity
    B = @match k begin
        (:P)         => identity
        (:E  || :B)  => Basis2Prod{𝐄𝐁,B₂}
        (:Q  || :U)  => Basis2Prod{𝐐𝐔,B₂}
        (:Ex || :Bx) => Basis2Prod{𝐄𝐁,Map}
        (:El || :Bl) => Basis2Prod{𝐄𝐁,Fourier}
        (:Qx || :Ux) => Basis2Prod{𝐐𝐔,Map}
        (:Ql || :Ul) => Basis2Prod{𝐐𝐔,Fourier}
        _ => throw(ArgumentError("Invalid FlatS2 index: $k"))
    end
    maybe_unfold(getproperty(B(f),k))
end
function getindex(f::FlatS02{Basis3Prod{B₁,B₂,B₃}}, k::Symbol; full_plane=false) where {B₁,B₂,B₃}
    maybe_unfold = (full_plane && k in [:Il,:El,:Bl,:Ql,:Ul]) ? x->unfold(x,fieldinfo(f).Ny) : identity
    B = @match k begin
        (:I  || :P)  => identity
        (:E  || :B)  => Basis3Prod{𝐈,𝐄𝐁,B₃}
        (:Q  || :U)  => Basis3Prod{𝐈,𝐐𝐔,B₃}
        (:Ix)        => Basis3Prod{𝐈,B₂,Map}
        (:Il)        => Basis3Prod{𝐈,B₂,Fourier}
        (:Ex || :Bx) => Basis3Prod{𝐈,𝐄𝐁,Map}
        (:El || :Bl) => Basis3Prod{𝐈,𝐄𝐁,Fourier}
        (:Qx || :Ux) => Basis3Prod{𝐈,𝐐𝐔,Map}
        (:Ql || :Ul) => Basis3Prod{𝐈,𝐐𝐔,Fourier}
        _ => throw(ArgumentError("Invalid FlatS02 index: $k"))
    end
    maybe_unfold(getproperty(B(f),k))
end
function getindex(D::DiagOp{<:FlatEBFourier}, k::Symbol)
    @unpack El, Bl, metadata = diag(D)
    @unpack sin2ϕ, cos2ϕ = fieldinfo(diag(D))
    f = @match k begin
        (:QQ)        => FlatFourier((@. Bl*sin2ϕ^2 + El*cos2ϕ^2),   metadata)
        (:QU || :UQ) => FlatFourier((@. (El - Bl) * sin2ϕ * cos2ϕ), metadata)
        (:UU)        => FlatFourier((@. Bl*cos2ϕ^2 + El*sin2ϕ^2),   metadata)
        _            => getproperty(D.diag, k)
    end
    Diagonal(f)
end
function getindex(L::FlatIEBCov, k::Symbol)
    @match k begin
        :IP => L
        :I => L.ΣTE[1,1]
        :E => L.ΣTE[2,2]
        :B => L.ΣB
        :P => Diagonal(FlatEBFourier(L[:E].diag, L[:B].diag))
        (:QQ || :UU || :QU || :UQ) => getindex(L[:P], k)
        _ => throw(ArgumentError("Invalid FlatIEBCov index: $k"))
    end
end

### basis conversion
## spin-0
Fourier(f::FlatMap) = FlatFourier(m_rfft(f.arr, (1,2)), f.metadata)
Fourier(f′::FlatFourier, f::FlatMap) = (m_rfft!(f′.arr, f.arr, (1,2)); f′)
Map(f::FlatFourier) = FlatMap(m_irfft(f.arr, f.Ny, (1,2)), f.metadata)
Map(f′::FlatMap, f::FlatFourier) = (m_irfft!(f′.arr, f.arr, (1,2)); f′)

## spin-2
QUFourier(f::FlatQUMap) = FlatQUFourier(m_rfft(f.arr, (1,2)), f.metadata)
QUFourier(f::FlatEBMap) = f |> EBFourier |> QUFourier
QUFourier(f::FlatEBFourier) = begin
    @unpack El, Bl, sin2ϕ, cos2ϕ = fieldinfo(f)
    Ql = @. - El * cos2ϕ + Bl * sin2ϕ
    Ul = @. - El * sin2ϕ - Bl * cos2ϕ
    FlatQUFourier(Ql, Ul, f.metadata)
end

QUMap(f::FlatQUFourier) = FlatQUMap(m_irfft(f.arr, f.Ny, (1,2)), f.metadata)
QUMap(f::FlatEBMap)      = f |> EBFourier |> QUFourier |> QUMap
QUMap(f::FlatEBFourier)  = f |> QUFourier |> QUMap

EBFourier(f::FlatEBMap) = FlatEBFourier(m_rfft(f.arr, (1,2)), f.metadata)
EBFourier(f::FlatQUMap) = f |> QUFourier |> EBFourier
EBFourier(f::FlatQUFourier) = begin
    @unpack Ql, Ul, sin2ϕ, cos2ϕ = fieldinfo(f)
    El = @. - Ql * cos2ϕ - Ul * sin2ϕ
    Bl = @.   Ql * sin2ϕ - Ul * cos2ϕ
    FlatEBFourier(El, Bl, f.metadata)
end

EBMap(f::FlatEBFourier) = FlatEBMap(m_irfft(f.arr, f.Ny, (1,2)), f.metadata)
EBMap(f::FlatQUMap)     = f |> QUFourier |> EBFourier |> EBMap
EBMap(f::FlatQUFourier) = f |> EBFourier |> EBMap

# in-place
QUMap(f′::FlatQUMap, f::FlatQUFourier) = (m_irfft!(f′.arr, f.arr, (1,2)); f′)
QUFourier(f′::FlatQUFourier, f::FlatQUMap) = (m_rfft!(f′.arr, f.arr, (1,2)); f′)

## spin-(0,2)
IQUFourier(f::FlatIQUMap) = FlatIQUFourier(m_rfft(f.arr, (1,2)), f.metadata)
IQUFourier(f::FlatIEBMap) = f |> IEBFourier |> IQUFourier
IQUFourier(f::FlatIEBFourier) = FlatIQUFourier(f.I, QUFourier(f.P))

IQUMap(f::FlatIQUFourier) = FlatIQUMap(m_irfft(f.arr, f.Ny, (1,2)), f.metadata)
IQUMap(f::FlatIEBMap)      = f |> IEBFourier |> IQUFourier |> IQUMap
IQUMap(f::FlatIEBFourier)  = f |> IQUFourier |> IQUMap

IEBFourier(f::FlatIEBMap) = FlatIEBFourier(m_rfft(f.arr, (1,2)), f.metadata)
IEBFourier(f::FlatIQUMap) = f |> IQUFourier |> IEBFourier
IEBFourier(f::FlatIQUFourier) = FlatIEBFourier(f.I, EBFourier(f.P))

IEBMap(f::FlatIEBFourier) = FlatIEBMap(m_irfft(f.arr, f.Ny, (1,2)), f.metadata)
IEBMap(f::FlatIQUMap)     = f |> IQUFourier |> IEBFourier |> IEBMap
IEBMap(f::FlatIQUFourier) = f |> IEBFourier |> IEBMap

# in-place
IQUMap(f′::FlatIQUMap, f::FlatIQUFourier) = (m_irfft!(f′.arr, f.arr, (1,2)); f′)
IQUFourier(f′::FlatIQUFourier, f::FlatIQUMap) = (m_rfft!(f′.arr, f.arr, (1,2)); f′)

# spin-0 bases applied to spin-2 and spin-(0,2)
Fourier(f::FlatField{B}) where {B<:BasisProd} = Fourier(B)(f)
Map(f::FlatField{B}) where {B<:BasisProd} = Map(B)(f)

### pretty printing
typealias_def(::Type{F}) where {B,M<:FlatProj,T,A,F<:FlatField{B,M,T,A}} = "Flat$(typealias(B)){$(typealias(A)),$(typealias(M))}"
function Base.summary(io::IO, f::FlatField)
    @unpack Nx,Ny,Nbatch,θpix = f
    print(io, "$(length(f))-element $Ny×$Nx$(Nbatch==1 ? "" : "(×$Nbatch)")-pixel $(θpix)′-resolution ")
    Base.showarg(io, f, true)
end

### dot products
# do in Map space (the LenseBasis, Ł) for simplicity
function dot(a::FlatField, b::FlatField)
    z = Ł(a) .* Ł(b)
    batch(sum_dropdims(z.arr, dims=nonbatch_dims(z)))
end

### logdets

function logdet(L::Diagonal{<:Union{Real,Complex},<:FlatField{B}}) where {B<:Union{Fourier,Basis2Prod{<:Any,Fourier},Basis3Prod{<:Any,<:Any,Fourier}}}
    # half the Fourier plane needs to be counted twice since the real
    # FFT only stores half of it
    @unpack Ny, arr = L.diag
    λ = adapt(typeof(arr), rfft_degeneracy_fac(Ny))
    # note: since our maps are required to be real, the logdet of any
    # operator which preserves this property is also guaranteed to be
    # real, hence the `real` and `abs` below are valid
    batch(real.(sum_dropdims(nan2zero.(log.(abs.(arr)) .* λ), dims=nonbatch_dims(L.diag))))
end

function logdet(L::Diagonal{<:Real,<:FlatField{B}}) where {B<:Union{Map,Basis2Prod{<:Any,Map},Basis3Prod{<:Any,<:Any,Map}}}
    batch(
        sum_dropdims(log.(abs.(L.diag.arr)), dims=nonbatch_dims(L.diag))
        .+ dropdims(log.(prod(sign.(L.diag.arr), dims=nonbatch_dims(L.diag))), dims=nonbatch_dims(L.diag))
    )
end

### traces

function tr(L::Diagonal{<:Union{Real,Complex},<:FlatField{B}}) where {B<:Union{Fourier,Basis2Prod{<:Any,Fourier},Basis3Prod{<:Any,<:Any,Fourier}}}
    @unpack Ny, Nx, arr = L.diag
    λ = adapt(typeof(arr), rfft_degeneracy_fac(Ny))
    # the `real` is ok bc the imaginary parts of the half-plane which
    # is stored would cancel with those from the other half-plane
    batch(real.(sum_dropdims(arr .* λ, dims=nonbatch_dims(L.diag))))
end

function tr(L::Diagonal{<:Real,<:FlatField{B}}) where {B<:Union{Map,Basis2Prod{<:Any,Map},Basis3Prod{<:Any,<:Any,Map}}}
    batch(sum_dropdims(L.diag.arr, dims=nonbatch_dims(L.diag)))
end

### simulation
_white_noise(ξ::FlatField, rng::AbstractRNG) =
    (randn!(rng, similar(ξ.arr, real(eltype(ξ)), ξ.Ny, size(ξ.arr)[2:end]...)), ξ.metadata)
white_noise(ξ::FlatS0,  rng::AbstractRNG) = FlatMap(_white_noise(ξ,rng)...)
white_noise(ξ::FlatS2,  rng::AbstractRNG) = FlatEBMap(_white_noise(ξ,rng)...)
white_noise(ξ::FlatS02, rng::AbstractRNG) = FlatIEBMap(_white_noise(ξ,rng)...)

### creating covariance operators
# fixed covariances
Cℓ_to_Cov(pol::Symbol, args...; kwargs...) = Cℓ_to_Cov(Val(pol), args...; kwargs...)
function Cℓ_to_Cov(::Val{:I}, proj::ProjLambert, Cℓ::InterpolatedCℓs; units=proj.Ωpix)
    Diagonal(FlatFourier(Cℓ_to_2D(Cℓ,proj), proj) / units)
end
function Cℓ_to_Cov(::Val{:P}, proj::ProjLambert, CℓEE::InterpolatedCℓs, CℓBB::InterpolatedCℓs; units=proj.Ωpix)
    Diagonal(FlatEBFourier(Cℓ_to_2D(CℓEE,proj), Cℓ_to_2D(CℓBB,proj), proj) / units)
end
function Cℓ_to_Cov(::Val{:IP}, proj::ProjLambert, CℓTT, CℓEE, CℓBB, CℓTE; kwargs...)
    ΣTT, ΣEE, ΣBB, ΣTE = [Cℓ_to_Cov(:I,proj,Cℓ; kwargs...) for Cℓ in (CℓTT,CℓEE,CℓBB,CℓTE)]
    FlatIEBCov(@SMatrix([ΣTT ΣTE; ΣTE ΣEE]), ΣBB)
end
# ParamDependentOp covariances scaled by amplitudes in different ℓ-bins
function Cℓ_to_Cov(::Val{:I}, proj::ProjLambert{T}, (Cℓ, ℓedges, θname)::Tuple; kwargs...) where {T}
    # we need an @eval here since we want to dynamically select a
    # keyword argument name, θname. the @eval happens into Main rather
    # than CMBLensing as a workaround for
    # https://discourse.julialang.org/t/closure-not-shipping-to-remote-workers-except-from-main/38831
    C₀ = diag(Cℓ_to_Cov(:I, proj, Cℓ; kwargs...))
    @eval Main let ℓedges=$((T.(ℓedges))...,), C₀=$C₀
        ParamDependentOp(function (;$θname=ones($T,length(ℓedges)-1),_...)
            _A = $preprocess.(Ref((nothing,C₀.metadata)), $T.($ensure1d($θname)))
            Diagonal(FlatFourier($bandpower_rescale!(ℓedges, C₀.ℓmag, copy(C₀.arr), _A...), C₀.metadata))
        end)
    end
end
function Cℓ_to_Cov(::Val{:P}, proj::ProjLambert{T}, (CℓEE, ℓedges, θname)::Tuple, CℓBB::InterpolatedCℓs; kwargs...) where {T}
    C₀ = diag(Cℓ_to_Cov(:P, proj, CℓEE, CℓBB; kwargs...))
    @eval Main let ℓedges=$((T.(ℓedges))...,), C₀=$C₀
        ParamDependentOp(function (;$θname=ones($T,length(ℓedges)-1),_...)
            _E = $preprocess.(Ref((nothing,C₀.metadata)),      $T.($ensure1d($θname)))
            _B = $preprocess.(Ref((nothing,C₀.metadata)), one.($T.($ensure1d($θname))))
            Diagonal(FlatEBFourier($bandpower_rescale!(ℓedges, C₀.ℓmag, copy(C₀.El), _E...), C₀.Bl .* _B[1], C₀.metadata))
        end)
    end
end
# this is written weird because the stuff inside the broadcast! needs
# to work as a GPU kernel
function bandpower_rescale!(ℓedges, ℓ, Cℓ, A...)
    length(A)==length(ℓedges)-1 || error("Expected $(length(ℓedges)-1) bandpower parameters, got $(length(A)).")
    eltype(A[1]) <: Real || error("Bandpower parameters must be real numbers.")
    if length(A)>30
        # if more than 30 bandpowers, we need to chunk the rescaling
        # because of a maximum argument limit of CUDA kernels
        for p in partition(1:length(A), 30)
            bandpower_rescale!(ℓedges[p.start:(p.stop+1)], ℓ, Cℓ, A[p]...)
        end
    else
        broadcast!(Cℓ, ℓ, Cℓ, A...) do ℓ, Cℓ, A...
            for i=1:length(ℓedges)-1
                (ℓedges[i] < ℓ < ℓedges[i+1]) && return A[i] * Cℓ
            end
            return Cℓ
        end
    end
    Cℓ
end
# cant reliably get Zygote's gradients to work through these
# broadcasts, which on GPU use ForwardDiff, so write the adjoint by
# hand for now. likely more performant, in any case.
@adjoint function bandpower_rescale!(ℓedges, ℓ, Cℓ, A...)
    function back(Δ)
        Ā = map(1:length(A)) do i
            sum(
                real,
                broadcast(Δ, ℓ, Cℓ) do Δ, ℓ, Cℓ
                    (ℓedges[i] < ℓ < ℓedges[i+1]) ? Cℓ*Δ : zero(Cℓ)
                end,
                dims = ndims(Δ)==4 ? (1,2) : (:)
            )
        end
        (nothing, nothing, nothing, Ā...)
    end
    bandpower_rescale!(ℓedges, ℓ, Cℓ, A...), back
end

### spin adjoints
function *(a::SpinAdjoint{F}, b::F) where {B<:Union{Map,Basis2Prod{<:Any,Map},Basis3Prod{<:Any,<:Any,Map}},F<:FlatField{B}}
    FlatMap(dropdims(sum(a.f.arr .* b.arr, dims=3), dims=((a.f.Nbatch>1 || b.Nbatch>1) ? 3 : ())), get_metadata_strict(a, b))
end
function mul!(dst::FlatMap, a::SpinAdjoint{F}, b::F) where {F<:FlatField{<:Union{Map,Basis2Prod{<:Any,Map},Basis3Prod{<:Any,<:Any,Map}}}}
    dst.arr .= reshape(sum(a.f.arr .* b.arr, dims=3), size(dst.arr))
    dst
end

### FlatIEBCov
# applying
*(L::FlatIEBCov, f::FlatS02) =       L * IEBFourier(f)
\(L::FlatIEBCov, f::FlatS02) = pinv(L) * IEBFourier(f)
function *(L::FlatIEBCov, f::FlatIEBFourier)
    (i,e),b = (L.ΣTE * [f.I, f.E]), L.ΣB * f.B
    FlatIEBFourier(i,e,b)
end
# manipulating
size(L::FlatIEBCov) = 3 .* size(L.ΣB)
adjoint(L::FlatIEBCov) = L
sqrt(L::FlatIEBCov) = FlatIEBCov(sqrt(L.ΣTE), sqrt(L.ΣB))
pinv(L::FlatIEBCov) = FlatIEBCov(pinv(L.ΣTE), pinv(L.ΣB))
global_rng_for(::Type{FlatIEBCov{T,F}}) where {T,F} = global_rng_for(F)
diag(L::FlatIEBCov) = FlatIEBFourier(L.ΣTE[1,1].diag, L.ΣTE[2,2].diag, L.ΣB.diag)
similar(L::FlatIEBCov) = FlatIEBCov(similar.(L.ΣTE), similar(L.ΣB))
get_storage(L::FlatIEBCov) = get_storage(L.ΣB)
simulate(rng::AbstractRNG, L::FlatIEBCov; Nbatch=nothing) =
    sqrt(L) * white_noise(similar(diag(L), (isnothing(Nbatch) || Nbatch==1 ? () : (Nbatch,))...), rng)
# arithmetic
*(L::FlatIEBCov, D::DiagOp{<:FlatIEBFourier}) = FlatIEBCov(SMatrix{2,2}(L.ΣTE * [[D[:I]] [0]; [0] [D[:E]]]), L.ΣB * D[:B])
+(L::FlatIEBCov, D::DiagOp{<:FlatIEBFourier}) = FlatIEBCov(@SMatrix[L.ΣTE[1,1]+D[:I] L.ΣTE[1,2]; L.ΣTE[2,1] L.ΣTE[2,2]+D[:E]], L.ΣB + D[:B])
*(La::F, Lb::F) where {F<:FlatIEBCov} = F(La.ΣTE * Lb.ΣTE, La.ΣB * Lb.ΣB)
+(La::F, Lb::F) where {F<:FlatIEBCov} = F(La.ΣTE + Lb.ΣTE, La.ΣB + Lb.ΣB)
+(L::FlatIEBCov, U::UniformScaling{<:Scalar}) = FlatIEBCov(@SMatrix[(L.ΣTE[1,1]+U) L.ΣTE[1,2]; L.ΣTE[2,1] (L.ΣTE[2,2]+U)], L.ΣB+U)
*(L::FlatIEBCov, λ::Scalar) = FlatIEBCov(L.ΣTE * λ, L.ΣB * λ)
*(D::DiagOp{<:FlatIEBFourier}, L::FlatIEBCov) = L * D
+(U::UniformScaling{<:Scalar}, L::FlatIEBCov) = L + U
*(λ::Scalar, L::FlatIEBCov) = L * λ
# copyto!(dst::Σ, src::Σ) where {Σ<:FlatIEBCov} = (copyto!(dst.ΣB, src.ΣB); copyto!.(dst.ΣTE, src.ΣTE); dst)

### batching

batch_length(f::FlatField) = f.Nbatch

"""
    batch(fs::FlatField...)
    batch(fs::Vector{<:FlatField})

Concatenate one of more FlatFields along the "batch" dimension
(dimension 4 of the underlying array). For the inverse operation, see
[`unbatch`](@ref).
"""
batch(fs::FlatField{B}...) where {B} = batch(collect(fs))
batch(fs::AbstractVector{<:FlatField{B}}) where {B} =
    FlatField{B}(cat(getfield.(fs,:arr)..., dims=Val(4)), only(unique(getfield.(fs,:metadata))))

"""
    unbatch(f::FlatField)

Return an Array of FlatFields corresponding to each batch index. For
the inverse operation, see [`batch`](@ref).
"""
unbatch(f::FlatField{B}) where {B} = [f[!,i] for i=1:batch_length(f)]

batch_index(f::FlatField{B}, I) where {B<:Union{Map,Fourier}} = FlatField{B}(f.arr[:,:,1,I], f.metadata)
batch_index(f::FlatField{B}, I) where {B} = FlatField{B}(f.arr[:,:,:,I], f.metadata)

###

make_mask(f::FlatField; kwargs...) = make_mask((f.Ny,f.Nx), f.θpix; kwargs...)

### power spectra

function get_Cℓ(f₁::FlatS0, f₂::FlatS0=f₁; Δℓ=50, ℓedges=0:Δℓ:16000, Cℓfid=ℓ->1, err_estimate=false)
    @unpack Nx, Ny, Δx, ℓmag = fieldinfo(f₁)
    ℓmag = unfold(ℓmag, Ny)
    α = Nx*Ny/Δx^2

    # faster to excise unused parts:
    ℓmask = (ℓmag .> minimum(ℓedges)) .&  (ℓmag .< maximum(ℓedges))
    L = ℓmag[ℓmask]
    CLobs = 1/α .* real.(dot.(
        adapt(Array,f₁)[:Il, full_plane=true][ℓmask],
        adapt(Array,f₂)[:Il, full_plane=true][ℓmask]
    ))
    w = @. nan2zero((2*Cℓfid(L)^2/(2L+1))^-1)

    sum_in_ℓbins(x) = fit(Histogram, L, Weights(x), ℓedges).weights

    local A, Cℓ, ℓ, N, Cℓ²
    Threads.@sync begin
        Threads.@spawn A  = sum_in_ℓbins(w)
        Threads.@spawn Cℓ = sum_in_ℓbins(w .* CLobs)
        Threads.@spawn ℓ  = sum_in_ℓbins(w .* L)
        if err_estimate
            Threads.@spawn N   = sum_in_ℓbins(one.(w)) / 2
            Threads.@spawn Cℓ² = sum_in_ℓbins(w .* CLobs.^2)
        end
    end

    if err_estimate
        σℓ  = sqrt.((Cℓ² ./ A .- Cℓ.^2) ./ N)
        InterpolatedCℓs(ℓ./A,  Cℓ./A .± σℓ)
    else
        InterpolatedCℓs(ℓ./A,  Cℓ./A)
    end
end

function get_Cℓ(f1::FlatS2, f2::FlatS2=f1; which=(:EE,:BB), kwargs...)
    Cℓ = (;[Symbol(x1*x2) => get_Cℓ(getindex(f1,Symbol(x1)),getindex(f2,Symbol(x2)); kwargs...) for (x1,x2) in split.(string.(ensure1d(which)),"")]...)
    which isa Symbol ? Cℓ[1] : Cℓ
end

function get_Cℓ(f1::FlatS02, f2::FlatS02=f1; which=(:II,:EE,:BB,:IE,:IB,:EB), kwargs...)
    Cℓ = (;[Symbol(x1*x2) => get_Cℓ(getindex(f1,Symbol(x1)),getindex(f2,Symbol(x2)); kwargs...) for (x1,x2) in split.(string.(ensure1d(which)),"")]...)
    which isa Symbol ? Cℓ[1] : Cℓ
end

"""
    ud_grade(f::Field, θnew, mode=:map, deconv_pixwin=true, anti_aliasing=true)

Up- or down-grades field `f` to new resolution `θnew` (only in integer steps).
Two modes are available specified by the `mode` argument:

* `:map`     — Up/downgrade by replicating/averaging pixels in map-space
* `:fourier` — Up/downgrade by extending/truncating the Fourier grid

For `:map` mode, two additional options are possible. If `deconv_pixwin` is
true, deconvolves the pixel window function from the downgraded map so the
spectrum of the new and old maps are the same. If `anti_aliasing` is true,
filters out frequencies above Nyquist prior to down-sampling.

"""
function ud_grade(
    f :: FlatField{B},
    θnew;
    mode = :map,
    deconv_pixwin = (mode==:map),
    anti_aliasing = (mode==:map)
) where {B}

    N,θ,T = (f.Ny, f.Nx), f.θpix, f.T
    θnew==θ && return f
    (mode in [:map,:fourier]) || throw(ArgumentError("Available modes: [:map,:fourier]"))

    fac = θnew > θ ? θnew÷θ : θ÷θnew
    (round(Int, fac) ≈ fac) || throw(ArgumentError("Can only ud_grade in integer steps"))
    fac = round(Int, fac)
    Ny_new, Nx_new = @. round(Int, N * θ ÷ θnew)
    proj = ProjLambert(;Ny=Ny_new, Nx=Nx_new, θpix=θnew, T=real(T), f.storage)
    @unpack Δx,ℓx,ℓy,Nx,Ny,nyquist = proj

    PWF = deconv_pixwin ? Diagonal(FlatFourier((@. T((pixwin(θnew,ℓy)*pixwin(θnew,ℓx)')/(pixwin(θ,ℓy)*pixwin(θ,ℓx)'))), proj)) : I

    if θnew > θ
        # downgrade
        if anti_aliasing
            f = Diagonal(FlatFourier(ifelse.((abs.(f.ℓy) .>= nyquist) .| (abs.(f.ℓx') .>= nyquist), 0, 1), f.metadata)) * f
        end
        if mode == :map
            fnew = FlatField{Map(B)}(dropdims(mean(reshape(Map(f).arr, fac, Ny, fac, Nx, size.(Ref(f.arr),nonbatch_dims(f)[3:end])...), dims=(1,3)), dims=(1,3)), proj)
        else
            fnew = FlatField{Fourier(B)}(Fourier(f).arr[1:(Ny_new÷2+1), [1:(isodd(Nx_new) ? Nx_new÷2+1 : Nx_new÷2); (end-Nx_new÷2+1):end], repeated(:, length(nonbatch_dims(f))-2)...], proj)
        end
        if deconv_pixwin
            fnew = Diagonal(FlatFourier((@. T((pixwin(θnew,ℓy)*pixwin(θnew,ℓx)')/(pixwin(θ,ℓy)*pixwin(θ,ℓx)'))), proj)) \ fnew
        end
    else
        error("Not implemented")
        # # upgrade
        # @assert fieldinfo(f).Nside isa Int "Upgrading resolution only implemented for maps where `Nside isa Int`"
        # if mode==:map
        #     fnew = FlatMap{Pnew}(permutedims(hvcat(N,(x->fill(x,(fac,fac))).(f[:Ix])...)))
        #     deconv_pixwin ? FlatFourier{Pnew}(fnew[:Il] .* PWF' .* PWF[1:Nnew÷2+1]) : fnew
        # else
        #     fnew = FlatFourier{Pnew}(zeros(Nnew÷2+1,Nnew))
        #     setindex!.(Ref(fnew.Il), f[:Il], 1:(N÷2+1), [findfirst(fieldinfo(fnew).k .≈ fieldinfo(f).k[i]) for i=1:N]')
        #     deconv_pixwin ? fnew * fac^2 : fnew
        # end

    end
    return fnew
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "PicoSATBuilder"
version = v"9.6.0"

# Collection of sources required to build PicoSATBuilder
sources = [
    "http://fmv.jku.at/picosat/picosat-960.tar.gz" =>
    "edb3184a04766933b092713d0ae5782e4a3da31498629f8bb2b31234a563e817",

]

# Bash recipe for building across all platforms
script = raw"""

cd $WORKSPACE/srcdir
cd picosat-960/
if [[ ${target} == *-apple-* ]]; then
    export CC=/opt/x86_64-linux-gnu/tools/clang
    export CFLAGS="--target=${target} --sysroot=/opt/${target}/${target}/sys-root"
    sed -i 's/-soname/-install_name/' makefile.in
else
    export CC=/opt/${target}/bin/${target}-gcc
fi

if [[ ${nbits} == 32 ]]; then
    ./configure --shared --32
else
    ./configure --shared
fi

make

if [[ ${target} == *-mingw* ]]; then
    mv libpicosat.so libpicosat.dll
elif [[ ${target} == *-apple-* ]]; then
    mv libpicosat.so libpicosat.dylib
fi

mv libpicosat* $prefix
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    MacOS(),
    Windows(:x86_64),
    Windows(:i686),
    Linux(:x86_64, :glibc),
    Linux(:i686, :glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libpicosat", :libpicosat)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using .Makie
using GeometryTypes: Cylinder
import .Makie.AbstractPlotting: plot!, plot, to_value

#######################################################################
# Tools
#######################################################################

function matrixidx(h::AbstractSparseMatrix, row, col)
    for ptr in nzrange(h, col)
        rowvals(h)[ptr] == row && return ptr
    end
    return 0
end

matrixidx(h::DenseMatrix, row, col) = LinearIndices(h)[row, col]

transparent(rgba::RGBAf0, v = 0.5) = RGBAf0(rgba.r, rgba.g, rgba.b, rgba.alpha * v)

function darken(rgba::RGBAf0, v = 0.66)
    r = max(0, min(rgba.r * (1 - v), 1))
    g = max(0, min(rgba.g * (1 - v), 1))
    b = max(0, min(rgba.b * (1 - v), 1))
    RGBAf0(r,g,b,rgba.alpha)
end
function lighten(rgba, v = 0.66)
    darken(rgba, -v)
end

#######################################################################
# plot(::Hamiltonian)
#######################################################################

function plot(h::Hamiltonian{<:Lattice}; resolution = (1000, 1000), kw...)
    scene = hamiltonianplot(h; resolution = resolution, scale_plot = false, kw...)
    plot = scene[end]
    plot[:tooltips][] && addtooltips!(scene, h)
    scale!(scene)
    return scene
end

@recipe(HamiltonianPlot, hamiltonian) do scene
    Theme(
        allintra = false, allcells = true, showsites = true, showlinks = true,
        shadedsites = false, shadedlinks = false, dimming = 0.75,
        siteradius = 0.2, siteborder = 3, siteborderdarken = 1.0, linkdarken = 0.0,
        linkthickness = 6, linkoffset = missing, linkradius = 0.1,
        tooltips = true, digits = 3,
        _tooltips_rowcolhar = Vector{Tuple{Int,Int,Int}}[],
        colors = map(t -> RGBAf0(t...),
            ((0.960,0.600,.327), (0.410,0.067,0.031),(0.940,0.780,0.000),
            (0.640,0.760,0.900),(0.310,0.370,0.650),(0.600,0.550,0.810),
            (0.150,0.051,0.100),(0.870,0.530,0.640),(0.720,0.130,0.250))),
        light = Vec3f0[[0, 0, 10], [0, 10, 0], [10, 0, 0], [10, 10, 10], [-10, -10, -10]]
    )
end

function plot!(plot::HamiltonianPlot)
    h = to_value(plot[1])
    lat = h.lattice
    colors = Iterators.cycle(plot[:colors][])
    sublats = Quantica.sublats(lat)

    mdist = Quantica.meandist(h)
    mdist > 0 || (mdist = 1)
    plot[:siteradius][] *= mdist/2
    plot[:linkradius][] *= mdist/2

    lo = plot[:linkoffset][]
    plot[:linkoffset][] =
        lo === missing ? √max(1-(plot[:linkradius][]/plot[:siteradius][])^2, 0) : lo * mdist

    # plot sites
    plot[:showsites][] &&
        for (n, har) in enumerate(h.harmonics), (ssrc, csrc) in zip(sublats, colors)
            iszero(har.dn) || plot[:allcells][] || break
            csrc´ = iszero(har.dn) ? csrc : transparent(csrc, 1 - plot[:dimming][])
            itr = siterange(lat, ssrc)
            plotsites!(plot, lat, itr, har.dn, n, csrc´)
        end

    # plot links
    plot[:showlinks][] &&
        for (n, har) in enumerate(h.harmonics)
            iszero(har.dn) || plot[:allcells][] || break
            for (ssrc, csrc) in zip(sublats, colors)
                csrc´ = iszero(har.dn) ? csrc : transparent(csrc, 1 - plot[:dimming][])
                csrc´ = darken(csrc´, plot[:linkdarken][])
                for (sdst, cdst) in zip(sublats, colors)
                    itr = Quantica.nonzero_indices(har, siterange(lat, sdst), siterange(lat, ssrc))
                    plotlinks!(plot, lat, itr, har.dn, n, csrc´)
                end
            end
        end

    return plot
end

function plotsites!(plot, lat, srange, dn, n, color)
    allsites = Quantica.sites(lat)
    br = lat.bravais.matrix
    sites = [padright(allsites[i] + br * dn, Val(3)) for i in srange]
    plot[:tooltips][] && (tt = [(site, 0, n) for site in srange])
    if !isempty(sites)
        plot[:shadedsites][] ? plotsites_hi!(plot, sites, color) : plotsites_lo!(plot, sites, color)
        plot[:tooltips][] && push!(plot[:_tooltips_rowcolhar][], tt)
    end
    return plot
end

function plotsites_lo!(plot, sites, color)
    scatter!(plot, sites;
        color = color,
        markersize = 2 * plot[:siteradius][],
        strokewidth = plot[:siteborder][],
        strokecolor = darken(color, plot[:siteborderdarken][]))
    return nothing
end

function plotsites_hi!(plot, sites, color)
    meshscatter!(plot, sites;
        color = color,
        markersize = plot[:siteradius][], light = plot[:light][])
    return nothing
end

function plotlinks!(plot, lat, itr, dn, n, color)
    links = Pair{SVector{3,Float32},SVector{3,Float32}}[]
    plot[:tooltips][] && (tt = Tuple{Int,Int,Int}[])
    sites = Quantica.sites(lat)
    br = lat.bravais.matrix
    for (row, col) in itr
        iszero(dn) && row == col && continue
        rdst = padright(sites[row] + br * dn, Val(3))
        rsrc = padright(sites[col], Val(3))
        rdst = iszero(dn) ? (rdst + rsrc) / 2 : rdst
        rsrc = rsrc + plot[:siteradius][] * plot[:linkoffset][] * normalize(rdst - rsrc)
        push!(links, rsrc => rdst)
        plot[:tooltips][] && push!(tt, (row, col, n))
    end
    if !isempty(links)
        plot[:shadedlinks][] ? plotlinks_hi!(plot, links, color) : plotlinks_lo!(plot, links, color)
        plot[:tooltips][] && push!(plot[:_tooltips_rowcolhar][], tt)
    end
    return plot
end

function plotlinks_lo!(plot, links, color)
    # linesegments!(plot, links; color = darken(color, plot[:siteborderdarken][]), linewidth = 2+plot[:linkthickness][])
    linesegments!(plot, links; color = color, linewidth = plot[:linkthickness][])
    return nothing
end

function plotlinks_hi!(plot, links, color)
    positions = [(r1 + r2) / 2 for (r1, r2) in links]
    rotvectors = [r2 - r1 for (r1, r2) in links]
    radius = plot[:linkradius][]
    scales = [Vec3f0(radius, radius, norm(r2 - r1)/2) for (r1, r2) in links]
    cylinder = Cylinder(Point3f0(0., 0., -1.0), Point3f0(0., 0, 1.0), Float32(1))
    meshscatter!(plot, positions;
        color = color, marker = cylinder, markersize = scales, rotations = rotvectors,
        light = plot[:light][])
    return nothing
end

function addtooltips!(scene, h)
    sceneplot = scene[end]
    visible = Node(false)
    N = Quantica.blockdim(h)
    poprect = lift(scene.events.mouseposition) do mp
        FRect((mp .+ 5), 1,1)
    end
    textpos = lift(scene.events.mouseposition) do mp
        Vec3f0((mp .+ 5 .+ (0, -50))..., 0)
    end
    popup = poly!(campixel(scene), poprect, raw = true, color = RGBAf0(1,1,1,0), visible = visible)
    rect = popup[end]
    translate!(popup, Vec3f0(0, 0, 10000))
    text!(popup, " ", textsize = 30, position = textpos, align = (:center, :center),
        color = :black, strokewidth = 4, strokecolor = :white, raw = true, visible = visible)
    text_field = popup[end]
    on(scene.events.mouseposition) do event
        subplot, idx = mouse_selection(scene)
        layer = findfirst(isequal(subplot), sceneplot.plots)
        if layer !== nothing && idx > 0
            idx´ = fix_linesegments_bug(idx, subplot)
            txt = popuptext(sceneplot, layer, idx´, h)
            text_field[1] = txt
            visible[] = true
        else
            visible[] = false
        end
        return
    end
    return scene
end

# idx returned by mouse_selection seems wrong by a factor 2 in LineSegments subplot
fix_linesegments_bug(idx, subplot::LineSegments) = Int(idx/2) # catches odd idx
fix_linesegments_bug(idx, subplot) = idx

function popuptext(sceneplot, layer, idx, h)
    cache = sceneplot[:_tooltips_rowcolhar][]
    checkbounds(Bool, cache, layer) || return string("Bug in layer: ", layer, " of ", length(cache))
    checkbounds(Bool, cache[layer], idx) || return string("Bug in idx: ", idx, " of ", length(cache[layer]))
    (row, col_or_zero, haridx) = cache[layer][idx]
    if col_or_zero == 0
        col = iszero(col_or_zero) ? row : col_or_zero
        har = h.harmonics[1]
    else
        col = col_or_zero
        har = h.harmonics[haridx]
    end
    element = round.(har.h[row, col], digits = sceneplot[:digits][])
    txt = iszero(col_or_zero) ? matrixstring(col, element) : matrixstring(row, col, element)
    return txt
end

#######################################################################
# plot(::Bandstructure)
#######################################################################

function plot(bs::Bandstructure{1}; kw...)
    scene = bandplot2d(bs; kw...)
    axis = scene[Axis]
    axis[:names, :axisnames] = ("φ", "ε")
    return scene
end

function plot(bs::Bandstructure{2}; kw...)
    scene = bandplot3d(bs; kw...)
    axis = scene[Axis]
    to_value(scene[:show_axis]) && (axis[:names, :axisnames] = ("φ₁", "φ₂", "ε"))
    return scene
end

@recipe(BandPlot2D, bandstructure) do scene
    Theme(
    linewidth = 3,
    wireframe = true,
    colors = map(t -> RGBAf0((0.8 .* t)...),
        ((0.973, 0.565, 0.576), (0.682, 0.838, 0.922), (0.742, 0.91, 0.734),
         (0.879, 0.744, 0.894), (1.0, 0.84, 0.0), (1.0, 1.0, 0.669),
         (0.898, 0.762, 0.629), (0.992, 0.843, 0.93), (0.88, 0.88, 0.88)))
    )
end

function plot!(plot::BandPlot2D)
    bs = to_value(plot[1])
    bands = haskey(plot, :bands) ? to_value(plot[:bands]) : eachindex(bs.bands)
    colors = Iterators.cycle(plot[:colors][])
    for (nb, color) in zip(bands, colors)
        band = bs.bands[nb]
        vertices = band.mesh.vertices
        simplices = band.simplices
        linesegments!(plot, (t -> vertices[first(t)] => vertices[last(t)]).(simplices),
                      linewidth = plot[:linewidth], color = color)
    end
    return plot
 end

@recipe(BandPlot3D, bandstructure) do scene
    Theme(
    linewidth = 1,
    wireframe = false,
    colors = map(t -> RGBAf0(t...),
        ((0.973, 0.565, 0.576), (0.682, 0.838, 0.922), (0.742, 0.91, 0.734),
         (0.879, 0.744, 0.894), (1.0, 0.84, 0.0), (1.0, 1.0, 0.669),
         (0.898, 0.762, 0.629), (0.992, 0.843, 0.93), (0.88, 0.88, 0.88)))
    )
end

function plot!(plot::BandPlot3D)
    bs = to_value(plot[1])
    bands = haskey(plot, :bands) ? to_value(plot[:bands]) : eachindex(bs.bands)
    colors = Iterators.cycle(plot[:colors][])
    for (nb, color) in zip(bands, colors)
        band = bs.bands[nb]
        vertices = band.mesh.vertices
        connectivity = [s[j] for s in band.simplices, j in 1:3]
        if isempty(connectivity)
            scatter!(plot, vertices, color = color)
        else
            mesh!(plot, vertices, connectivity, color = color, transparency = false)
            if plot[:wireframe][]
                edgevertices = collect(Quantica.edgevertices(band.mesh))
                linesegments!(plot, edgevertices, linewidth = plot[:linewidth])
            end
        end
    end
    return plot
 end

function allpairs(simplex::NTuple{D,Int}, vertices) where D
    ntuple(i -> (i))
end


text:
    Pkg.activate(".")
    srcDir  = dirname(@__FILE__)*"\\";
    include(srcDir*"Examples\\Examples.jl")
    include(srcDir*"Examples\\Utils_states.jl")
    include(srcDir*"Moments.jl")
    include(srcDir*"Constraints\\Utils_cons.jl")
    include(srcDir*"Constraints\\R_constraints.jl")
    include(srcDir*"Constraints\\C_constraints.jl")
    include(srcDir*"Model\\Utils_Model.jl")
    include(srcDir*"Model\\R_sep_Model.jl")
    include(srcDir*"Model\\C_sep_Model.jl")
    include(srcDir*"sep_Compute.jl")

    using .Examples      ; const ex = Examples
    using .Utils_states  ; const us = Utils_states
    using .Moments       ; const mom = Moments
    using .Utils_cons    ; const uc = Utils_cons
    using .C_constraints ; const ccon = C_constraints
    using .R_constraints ; const rcon = C_constraints
    using .Utils_Model
    using .R_sep_Model
    using .C_sep_Model
    using .sep_Compute

    using LinearAlgebra ; const la = LinearAlgebra # Pkg.add("LinearAlgebra")
    using JuMP

## Try d = (3,3), t = 3 without block-diag
ρ = ex.get_example("RCai")
d = ex.get_example_meta("RCai").Size[1]
t = 3
ℂ1_sep_mod = sep_Compute.Computeξₜˢᵉᵖ(ρ,d,t;con_list="S1G",noBlock=true)

## Try d = (2,2), t = 4 with block-diag
d = (7,7) ; r = 5 ; t = 2
ρ = Examples.gen_ρ_RAND(d, r, true)
ℂ1_sep_mod = sep_Compute.Computeξₜˢᵉᵖ(ρ,d,t;con_list="S1G",noBlock=false)

## Get variable reduction
d = (3,3) ; t = 3
## Block
MMᴿcoef,MMᴿexᴿ = mom.get_ℂ_block_diag(d,t,noBlock=false)
Unique_moms = unique(vcat([unique(vcat([[b...] for b in [MMᴿexᴿ[B]...]]...)) for B in keys(MMᴿexᴿ)]...))

## No block
binomial(sum(2 .* d) +2*t,2*t)

max_stars_repo_name: JAndriesJ/sep-rank

max_issues_repo_name: JAndriesJ/sep-rank

max_forks_repo_name: JAndriesJ/sep-rank


text:
module Input

export in_grab, in_release
export bind, bindlist, unbind

import GLFW
import Player

in_grabbed = false
function in_grab()
	GLFW.Disable(GLFW.MOUSE_CURSOR)
	GLFW.SetMousePos(0, 0)
	global in_grabbed = true
end
function in_release()
	global in_grabbed = false
	GLFW.Enable(GLFW.MOUSE_CURSOR)
end

const MOUSE_WHEEL_DOWN = (GLFW.KEY_LAST+1)
const MOUSE_WHEEL_UP   = (GLFW.KEY_LAST+2)

bindlist = Dict{Int,Function}()
bind(key::Integer) = get(bindlist, int(key), None)
bind(key::Integer, action::Function) = setindex!(bindlist, action, int(key))
unbind(key::Integer) = delete!(bindlist, int(key))

function event(key::Int, press::Bool)
	action = bind(key)
	if action != None && (in_grabbed || action == in_grab)
		if applicable(action, press)
			action(press)
		elseif press
			action()
		end
	end
end

# GLFW key and mouse button callback
function event(key::Cint, press::Cint)
	event(int(key), press == 1)
	return
end

# GLFW mouse wheel callback
function wheel_event(val::Cint)
	if val < 0
		event(MOUSE_WHEEL_DOWN, true)
	elseif val > 0
		event(MOUSE_WHEEL_UP, true)
	end
	GLFW.SetMouseWheel(0)
	return
end

m_pitch = 0.05
m_yaw = 0.05

function look_event(x::Cint, y::Cint)
	if in_grabbed && (x != 0 || y != 0)
		Player.lookdir!(Player.self, m_pitch * -y, m_yaw * -x)
		GLFW.SetMousePos(0, 0)
	end
	return
end

end