text:
# Observer
abstract type Observer end

# Lazy Object
abstract type LazyObject <: Observer end

# Process
abstract type StochasticProcess end
abstract type StochasticProcess1D <: StochasticProcess end
abstract type AbstractBlackScholesProcess <: StochasticProcess1D end
abstract type AbstractHestonProcess <: StochasticProcess1D end
abstract type AbstractDiscretization end
abstract type AbstractHestonDiscretization <: AbstractDiscretization end
abstract type BlackScholesType end

# Methods
abstract type Lattice end
abstract type TreeLattice <: Lattice end
abstract type AbstractTree end
abstract type AbstractBinomialTree <: AbstractTree end
abstract type BinomialTreeType end
abstract type EqualProbabilitiesBinomialTreeType <: BinomialTreeType end
abstract type EqualJumpsBinomialTreeType <: BinomialTreeType end
abstract type FdmSchemeDescType end
abstract type FdScheme end
abstract type FdmMesher end
abstract type Fdm1DMesher end
abstract type StepCondition end
abstract type CurveWrapper end
abstract type StepType end
abstract type FdmInnerValueCalculator end
abstract type FdmLinearOpComposite end
abstract type NinePointLinearOp end
abstract type BoundaryConditionType end
abstract type BC_Side end
# Monte Carlo
abstract type AbstractMonteCarloModel end
abstract type MCTrait end
abstract type AbstractPathPricer end
abstract type EarlyExercisePathPricer <: AbstractPathPricer end
abstract type LsmBasisSystemPolynomType end
abstract type LSMBasisSystemFunction <: Function end
# abstract BoundaryCondition

# Exercise
abstract type Exercise end
abstract type EarlyExercise <: Exercise end

# Instruments
abstract type Instrument <: LazyObject end
abstract type PositionType end
abstract type Bond <: Instrument end
abstract type AbstractCallableBond <: Bond end
abstract type AbstractRate <: Instrument end
abstract type Swap <: Instrument end
abstract type AbstractClaim end
abstract type SettlementType end
abstract type AbstractPayoff end
abstract type StrikedTypePayoff <: AbstractPayoff end
abstract type Option{E} <: Instrument end
abstract type OneAssetOption{E} <: Option{E} end
abstract type OptionType end
abstract type SwapType end
abstract type CDSProtectionSide end
abstract type Results end
abstract type CallType end
abstract type AbstractForward <: Instrument end

# Term Structures
abstract type TermStructure <: LazyObject end
abstract type YieldTermStructure <: TermStructure end
abstract type CreditTermStructure <: TermStructure end
# Curves
abstract type Curve <: YieldTermStructure end
abstract type InterpolatedCurve{P} <: Curve end
abstract type AbstractDefaultProbabilityTermStructure <: CreditTermStructure end
abstract type AbstractDefaultProbabilityCurve <: AbstractDefaultProbabilityTermStructure end
abstract type InterpolatedDefaultProbabilityCurve{P} <: AbstractDefaultProbabilityCurve end

abstract type VolatilityTermStructure <: TermStructure end
abstract type OptionletVolatilityStructure <: VolatilityTermStructure end
abstract type SwaptionVolatilityStructure <: VolatilityTermStructure end
abstract type BlackVolTermStructure <: VolatilityTermStructure end
abstract type LocalVolTermStructure <: VolatilityTermStructure end
abstract type CallableBondVolatilityStructure <: TermStructure end
abstract type AbstractSmileSection end
abstract type VolatilityType end
abstract type BootstrapTrait end
abstract type Bootstrap end
abstract type FittingMethod end
abstract type BootstrapHelper <: LazyObject end
abstract type BondHelper <: BootstrapHelper end
abstract type RateHelper <: BootstrapHelper end
abstract type AbstractCDSHelper <: BootstrapHelper end

# Pricing Engines
abstract type PricingEngine end
abstract type DiscretizedAsset end
abstract type DiscretizedOption <: DiscretizedAsset end
abstract type LatticeShortRateModelEngine{S} <: PricingEngine end
abstract type AbstractHestonEngine{HI} <: PricingEngine end
abstract type AbstractFDVanillaEngine <: PricingEngine end
abstract type FDMultiPeriodEngine <: AbstractFDVanillaEngine end
abstract type FDStepConditionEngine <: AbstractFDVanillaEngine end
abstract type AbstractVanillaEngine <: PricingEngine end
abstract type MCEngine{S, RSG} <: PricingEngine end
abstract type MCVanillaEngine{S, RSG} <: MCEngine{S, RSG} end
abstract type MCLongstaffSchwartzEngine{S, P, RSG} <: MCEngine{S, RSG} end
abstract type GaussianProbabilities end

# Cash Flows
abstract type CashFlows end
abstract type Leg <: CashFlows end
abstract type CashFlow end
abstract type Coupon <: CashFlow end
abstract type Duration end
abstract type CouponPricer end
abstract type IborCouponPricer <: CouponPricer end

# Indexes
abstract type InterestRateIndex end

# Models
abstract type Parameter end
abstract type CalibrationErrorType end
abstract type CalibrationHelper <: LazyObject end
abstract type CalibrationBasketType end
abstract type ModelType end
abstract type Model{T <: ModelType} <: LazyObject end
abstract type ShortRateModel{T} <: Model{T} end
abstract type OneFactorModel{T} <: ShortRateModel{T} end
abstract type Gaussian1DModel{T} <: OneFactorModel{T} end
abstract type TwoFactorModel{T} <: ShortRateModel{T} end
abstract type ComplexLogFormula end
abstract type HestonIntegration end
abstract type ShortRateDynamics end
abstract type ShortRateTree <: TreeLattice end
abstract type AbstractMarketModel end
abstract type AbstractMarketModelEvolver end
abstract type ExerciseStrategy end
abstract type MarketModelMultiProduct end
abstract type MarketModelBasisSystem end
abstract type MarketModelPathwiseMultiProduct end
abstract type MultiProductMultiStep <: MarketModelMultiProduct end
abstract type MarketModelExerciseValue end
abstract type PiecewiseConstantCorrelation end
abstract type BrownianGeneratorFactory end
abstract type BrownianGenerator end
abstract type SobolOrdering end
abstract type CurveState end

# Currencies
abstract type AbstractCurrency end

# Interest Rates
abstract type CompoundingType end


text:
using LinearAlgebra: norm
import Base: ==
const BLI = BLInvariants
using JuLIP: decode_dict

BondLengthDesc(transform::String, cutoff) =
         BondLengthDesc(AnalyticTransform(transform), fcut_analyse(cutoff))

# -------------- IO -------------------
Dict(D::BondLengthDesc) = Dict( "__id__"    =>  "BondLengthDesc",
                                "transform" =>  Dict(D.transform),
                                "cutoff"    =>  Dict(D.cutoff) )

BondLengthDesc(D::Dict) = BondLengthDesc( decode_dict(D["transform"]),
                                          decode_dict(D["cutoff"]) )

==(D1::BondLengthDesc, D2::BondLengthDesc) =
      ( (D1.transform == D2.transform) && (D1.cutoff == D2.cutoff) )

Base.convert(::Val{:BondLengthDesc}, D::Dict) = BondLengthDesc(D)

# ------------- Interface Code ---------------

tdegrees(::BondLengthDesc, vN::Val{N}) where {N} = BLI.tdegrees(vN)

@inline ricoords(D::BondLengthDesc, Rs, J) = edge_lengths(Rs, J)

@inline gradri2gradR!(desc::BondLengthDesc, dVsite, dV_dr, Rs, J, r) =
   _grad_len2pos!(dVsite, dV_dr, Rs, J, r)

@inline fcut(D::BondLengthDesc, r) = fcut(D.cutoff, r)
@inline fcut_d(D::BondLengthDesc, r) = fcut_d(D.cutoff, r)

@inline skip_simplex(D::BondLengthDesc, r) = (maximum(r) > cutoff(D.cutoff))

@inline invariants(D::BondLengthDesc, r) = BLI.invariants(transform.(Ref(D), r))

@inline function invariants_ed(D::BondLengthDesc, r)
   x = transform.(Ref(D), r)
   I1, I2, DI1, DI2 = BLI.invariants_ed(x)
   x_d = transform_d.(Ref(D), r)
   return I1, I2, _sdot(x_d, DI1), _sdot(x_d, DI2)
end

# ---------------- inner kernels for the edge length descriptor --------------

"""
`edge_lengths(Rs::AbstractVector, J::SVector{K, Int})`

Take collection `Rs::AbstractVector{JVec}` representing a simplex
and convert into an ordered vector of bondlengths, e.g.,
```
[ r01, r02, r03, r12, r13, r23 ]
[ r01, r02, r03, r04, r12, r13, r14, r23, r24, r34 ]
```

* `Rs` : collection of all R-vectors within a site energy computation
* `J` : indices of the vectors involved in the current nbody computation
"""
@generated function edge_lengths(Rs::AbstractVector, J::SVector{K, Int}) where K
   # note K = N-1 ]
   code = Expr[]
   idx = 0
   for n = 1:K
      idx += 1
      push_str!(code, "s_$idx = norm(Rs[J[$n]])")  # r_0n
   end
   for n = 1:K-1, m = (n+1):K
      idx += 1
      push_str!(code, "s_$idx = norm(Rs[J[$n]] - Rs[J[$m]])")   # r_nm
   end
   s_str = "s = @SVector [s_1"
   if idx > 1; s_str *= prod(", s_$i" for i = 2:idx); end
   push_str!(code, s_str * "]")
   quote
      $(Expr(:meta, :inline))
      @inbounds $(Expr(:block, code...))
      return s
   end
end

"""
convert ∇V (where ∇ is the gradient w.r.t. bond-lengths) into forces, i.e., into
∇Vsite (where ∇ is the gradient w.r.t. positions)

J : neighbour (sub-) indices
"""
@generated function _grad_len2pos!(dVsite, dV, Rs, J::SVector{K, Int}, r) where {K}
   # K is the number of neighbours, i.e. N = K+1 counting also the center atom
   # length(dV) == length(s) == length(S) == K * (K+1)/2
   # ------
   code = Expr[]
   idx = 0
   # the first K entries of dV, s, S are the |Ri - 0|
   for k = 1:K
      idx += 1   # idx == k of course
      push!(code, :( dVsite[J[$k]] += dV[$idx] * (Rs[J[$k]]/r[$idx]) ))
   end
   # the remaining ones are |R_i - R_j|
   for n = 1:K-1, m = (n+1):K
      idx += 1
      push!(code, :( S = (Rs[J[$n]]-Rs[J[$m]]) / r[$idx] ))
      push!(code, :( dVsite[J[$n]] += dV[$idx] * S ))
      push!(code, :( dVsite[J[$m]] -= dV[$idx] * S ))
   end
   quote
      $(Expr(:meta, :inline))
      @inbounds $(Expr(:block, code...))
      return dVsite
   end
end

# =============== Experimental Code for Cluster-based summation!!!!!!

ClusterBLDesc(transform::String, cutoff) =
         ClusterBLDesc(AnalyticTransform(transform), fcut_analyse(cutoff))

# -------------- IO -------------------
Dict(D::ClusterBLDesc) = Dict( "__id__"    =>  "ClusterBLDesc",
                                "transform" =>  Dict(D.transform),
                                "cutoff"    =>  Dict(D.cutoff) )

ClusterBLDesc(D::Dict) = ClusterBLDesc( decode_dict(D["transform"]),
                                        decode_dict(D["cutoff"]) )

==(D1::ClusterBLDesc, D2::ClusterBLDesc) =
      ( (D1.transform == D2.transform) && (D1.cutoff == D2.cutoff) )

Base.convert(::Val{:ClusterBLDesc}, D::Dict) = ClusterBLDesc(D)

ClusterBLDesc(D::BondLengthDesc) = ClusterBLDesc(D.transform, D.cutoff)

# ------------- Interface Code ---------------

tdegrees(::ClusterBLDesc, vN::Val{N}) where {N} = BLI.tdegrees(vN)

@inline ricoords(D::ClusterBLDesc, Rs, J) = edge_lengths(Rs, J)

@inline gradri2gradR!(desc::ClusterBLDesc, dVsite, dV_dr, Rs, J, r) =
   _grad_len2pos!(dVsite, dV_dr, Rs, J, r)

@inline fcut(D::ClusterBLDesc, r) = fcut(D.cutoff, r)
@inline fcut_d(D::ClusterBLDesc, r) = fcut_d(D.cutoff, r)

@inline skip_simplex(D::ClusterBLDesc, r) = (maximum(r) > cutoff(D.cutoff))

@inline invariants(D::ClusterBLDesc, r) = BLI.invariants(transform.(Ref(D), r))

@inline function invariants_ed(D::ClusterBLDesc, r)
   x = transform.(Ref(D), r)
   I1, I2, DI1, DI2 = BLI.invariants_ed(x)
   x_d = transform_d.(Ref(D), r)
   return I1, I2, _sdot(x_d, DI1), _sdot(x_d, DI2)
end

ext: jl

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
using LaTeXStrings
using Plots
using Plots.PlotMeasures
pyplot()

wp = JvWorker(grid_size=25)
v_init = wp.x_grid .* 0.5

f(x) = bellman_operator(wp, x)
V = compute_fixed_point(f, v_init, max_iter=300)

s_policy, ϕ_policy = bellman_operator(wp, V, ret_policies=true)

# === plot solution === #
p = plot(wp.x_grid, [ϕ_policy s_policy V],
         title=["ϕ policy" "s policy" "value function"],
         color=[:orange :blue :green],
         xaxis=("x", (0.0, maximum(wp.x_grid))),
         yaxis=((-0.1, 1.1)), size=(800, 800),
         legend=false, layout=(3, 1),
         bottom_margin=Measures.Length(:mm, 20))


text:
# ## Objectives
#
# Here we will show you how to perform an aerodynamic analysis of an airfoil.
# > **Recipe**
# > 1. Compute the coordinates of a NACA 4-digit series airfoil.
# > 2. Plot its camber, thickness, upper and lower surface representations.
# > 3. Perform an aerodynamic analysis at a given angle of attack.
# > 4. Plot its aerodynamic properties.
# For this, we will need to import some packages which will be convenient for plotting.
using AeroMDAO      # Main package
using Plots         # Plotting library
gr(dpi = 300)       # Plotting backend
using LaTeXStrings  # For LaTeX printing in plots

# ## Your First Airfoil
#
# You can define a NACA 4-digit series airfoil using the following function.
airfoil = naca4((2,4,1,2), 60)

# You can access the $x$- and $y$-coordinates as the following fields.
airfoil.x, airfoil.y

# You can obtain the coordinates as an array by calling the following function.
coordinates(airfoil)

# ### Geometric Representations
# You can convert these coordinates into the camber-thickness representation.
xcamthick = camber_thickness(airfoil, 60)

# You can split the coordinates into their upper and lower surfaces.
upper, lower = split_surface(airfoil);

x_upper, y_upper = @views upper[:,1], upper[:,2]
x_lower, y_lower = @views lower[:,1], lower[:,2];

# ### Plotting
# Let's plot everything!

## Plot object
af_plot = plot(aspect_ratio = 1, xlabel=L"(x/c)", ylabel = L"y")

## Upper surface
plot!(upper[:,1], upper[:,2], label = "$(airfoil.name) Upper",
      ls = :solid, lw = 2, c = :cornflowerblue)

## Lower surface
plot!(lower[:,1], lower[:,2], label = "$(airfoil.name) Lower",
      ls = :solid, lw = 2, c = :orange)

## Camber
plot!(xcamthick[:,1], xcamthick[:,2], label = "$(airfoil.name) Camber",
      ls = :dash, lw = 2, c = :burlywood3)

## Thickness
plot!(xcamthick[:,1], xcamthick[:,3], label = "$(airfoil.name) Thickness",
      ls = :dash, lw = 2, c = :grey)

# ## Your First Doublet-Source Analysis
#
# Now we have an airfoil, and we would like to analyze its aerodynamic characteristics. The potential flow panel method for inviscid analyses of airfoils, which you may have studied in your course in aerodynamics, provides decent estimations of the lift generated by the airfoil.

# Our analysis also requires boundary conditions, which is the freestream flow defined by a magnitude ``V_\infty`` and angle of attack ``\alpha``. We provide these to the analysis by defining variables and feeding them to a `Uniform2D` type, corresponding to uniform flow in 2 dimensions.
V       = 1.0
alpha   = 4.0 # degrees
uniform = Uniform2D(V, alpha)

# You can analyze this airfoil for the given freestream flow as shown.
system  = solve_case(
                     airfoil, uniform;
                     num_panels = 80
                    );

# This will run the analysis and return a system which can be used to obtain the aerodynamic quantities of interest and post-processing. The panels used for the analysis can be accessed as follows.
panels = system.surface_panels

# You can compute the lift coefficient from the system.
cl = lift_coefficient(system)

# You can also compute the sectional lift, moment and pressure coefficients.
cls, cms, cps = surface_coefficients(system)

# Note the difference between the lift coefficient computed and the sum of the sectional lift coefficients; this is due to numerical errors in the solution procedure and modeling.
cl, sum(cls)

#md # !!! info
#md #     Support for drag prediction with boundary layer calculations will be added soon. For now, try out the amazing [Webfoil](http://webfoil.engin.umich.edu/) developed by the MDOLab at University of Michigan -- Ann Arbor!

# ### Visualization
#
# Let's see what the pressure and lift distribution curves look like over the airfoil. AeroMDAO provides more helper functions for post-processing data. For example, you can make your fancy plots by segregating the values depending on the locations of the panels by defining the following function.
cp_upper, cp_lower = get_surface_values(panels, cps)

# Now let's plot the results!

## Pressure coefficients
plot(yflip = true, xlabel = L"(x/c)", ylabel = L"C_p", lw = 2)
plot!(cp_upper, label = L"$C_p$ Upper",
      ls = :dash, lw = 2, c = :cornflowerblue)
plot!(cp_lower, label = L"$C_p$ Lower",
      ls = :dash, lw = 2, c = :orange)
plot!(x_upper, -y_upper, label = "$(airfoil.name) Upper",
      ls = :solid, lw = 2, c = :cornflowerblue)
plot!(x_lower, -y_lower, label = "$(airfoil.name) Lower",
      ls = :solid, lw = 2, c = :orange)

#
## Lift coefficients
cl_upper, cl_lower = get_surface_values(panels, cls)

cl_plot = plot(xlabel = L"(x/c)", ylabel = L"C_l", lw = 2)
plot!(cl_upper, label = L"$C_l$ Upper",
      ls = :dash, lw = 2, c = :cornflowerblue)
plot!(cl_lower, label = L"$C_l$ Lower",
      ls = :dash, lw = 2, c = :orange)
plot!(x_lower, y_lower, label = "$(airfoil.name) Lower",
      ls = :solid, lw = 2, c = :orange)
plot!(x_upper, y_upper, label = "$(airfoil.name) Upper",
      ls = :solid, lw = 2, c = :cornflowerblue)

#
## Moment coefficients
cm_upper, cm_lower = get_surface_values(panels, cms)

cm_plot = plot(xlabel = L"(x/c)", ylabel = L"C_m")
plot!(cm_upper, label = L"$C_m$ Upper",
      ls = :dash, lw = 2, c = :cornflowerblue)
plot!(cm_lower, label = L"$C_m$ Lower",
      ls = :dash, lw = 2, c = :orange)
plot!(x_upper, y_upper, label = L"$C_m$ Upper",
      ls = :solid, lw = 2, c = :cornflowerblue)
plot!(x_lower, y_lower, label = L"$C_m$ Lower",
      ls = :solid, lw = 2, c = :orange)

# Great! We've created our first airfoil and run an aerodynamic analysis in 2 dimensions. Now we can move on to analyzing an aircraft configuration in the [Aircraft Aerodynamic Analysis](tutorials-aircraft.md) tutorial.


text:
module hello

using Libdl

lib = dlopen("../libs/libeosio_native.dylib") # Open the library explicitly.
sym_is_account = dlsym(lib, :is_account)
sym_action_data_size = dlsym(lib, :action_data_size)
sym_read_action_data = dlsym(lib, :read_action_data)
sym_string_to_uint64 = dlsym(lib, :string_to_uint64)
sym_uint64_to_string = dlsym(lib, :uint64_to_string)

function is_account(account::UInt64)
  ret = ccall(sym_is_account, Cuchar, (Culonglong,), account)
  return ret != 0
end

function action_data_size()
  ret = ccall(sym_action_data_size, Cuint, ())
  return ret
end

function read_action_data()
  size = action_data_size()
  buf = Array{UInt8, 1}(undef, size)
  ret = ccall(sym_read_action_data, Cuint, (Ptr{UInt8}, Cuint), buf, size)
  return String(buf)
end

function N(account::String)
  ret = ccall(sym_string_to_uint64, Culonglong, (Cstring,), account)
  return ret
end

function n2s(n::UInt64)
  buf = Array{UInt8, 1}(undef, 13)
  ret = ccall(sym_uint64_to_string, Cuint, (Culonglong, Ptr{UInt8}, Cuint), n, buf, sizeof(buf))
  return String(buf[1:ret])
end

function apply(receiver::UInt64, code::UInt64, action::UInt64)
  if action == N("sayhello")
    println(n2s(receiver))
    println(n2s(code))
    println(n2s(action))

    println(is_account(receiver))
    println(is_account(action))
    println(N("hello"))

    r = read_action_data()
    println(r)
  end
  return 1
end

end

max_stars_repo_name: learnforpractice/pyeos

max_issues_repo_name: openchatproject/safeos

max_forks_repo_name: learnforpractice/pyeos


text:
function like(
    state::Tuple{AbstractGraph, AbstractArray}, agent_idx::Integer,
    config::Config
)
    agent_list = state[2]
    this_agent = agent_list[agent_idx]
    inclin_interact = deepcopy(this_agent.inclin_interact)
    i = 1
    while inclin_interact > rand()
        if i < length(this_agent.feed)
            if ((abs(this_agent.feed[i].opinion - this_agent.opinion) < config.opinion_treshs.like)
                && !(this_agent.feed[i] in this_agent.liked_Tweets))
                this_agent.feed[i].like_count += 1
                this_agent.feed[i].weight *= 1.01
                push!(this_agent.liked_Tweets, this_agent.feed[i])
            end
        else
            break
        end
        i += 1
        inclin_interact -= 1
    end
    return state
end

max_stars_repo_name: digitalemuendigkeit/twitter-abm

max_issues_repo_name: digitalemuendigkeit/twitter-abm

max_forks_repo_name: digitalemuendigkeit/twitter-abm


text:
using Statistics
using ..Transformers.Basic: gather

using ZygoteRules: @adjoint, pullback

using AbstractTrees
import AbstractTrees: children

const ACT2FN = (gelu = gelu, relu = relu, swish = swish, gelu_new = gelu, mish = mish)

abstract type THModule end
abstract type HGFPreTrainedModel <: THModule end

# torch.nn.LayerNorm

struct FakeTHLayerNorm{T<:AbstractArray} <: THModule
  eps::Float32
  weight::T
  bias::T
end

Functors.functor(::Type{<:FakeTHLayerNorm}, layernorm) = (weight = layernorm.weight, bias = layernorm.bias), y -> FakeTHLayerNorm(layernorm.eps, y...)

(ln::FakeTHLayerNorm)(x) = ln.weight .* Flux.normalise(x, dims=1, ϵ=ln.eps) .+ ln.bias

function load_state!(layer::FakeTHLayerNorm, state)
  for k in keys(state)
    if k == :gamma
      nk = :weight
    elseif k == :beta
      nk = :bias
    else
      nk = k
    end
    load_state!(getfield(layer, nk), getfield(state, k))
  end
end

# torch.nn.Linear

struct FakeTHLinear{W<:AbstractArray, B<:Union{Nothing, AbstractArray}} <: THModule
  weight::W
  bias::B
end

FakeTHLinear(w) = FakeTHLinear(w, nothing)

_has_bias(::FakeTHLinear{W, Nothing}) where {W<:AbstractArray}= false
_has_bias(::FakeTHLinear) = true

function Functors.functor(::Type{<:FakeTHLinear}, linear)
  (_has_bias(linear) ? (weight = linear.weight, bias = linear.bias) : (weight = linear.weight,)),
  y -> FakeTHLinear(y...)
end

(l::FakeTHLinear)(x::AbstractMatrix) = _has_bias(l) ? l.weight * x .+ l.bias : l.weight * x
function (l::FakeTHLinear)(x::AbstractArray)
  old_size = size(x)
  new_size = Base.setindex(old_size, size(l.weight, 1), 1)

  new_x = reshape(x, old_size[1], :)
  y = l(new_x)
  return reshape(y, new_size)
end

# torch.nn.Embedding

struct FakeTHEmbedding{T<:AbstractArray} <: THModule
  pad_idx::Int
  weight::T
end

function FakeTHEmbedding(pad_idx, weight)
  @assert pad_idx <= size(weight, 2)
  if !iszero(pad_idx)
    @view(weight[:, pad_idx]) .= 0
  end

  FakeTHEmbedding(pad_idx, weight)
end

Functors.functor(::Type{<:FakeTHEmbedding}, embedding) = (weight = embedding.weight,), y -> FakeTHEmbedding(embedding.pad_idx, y...)

#TODO: refactor basic/embeds/gather.jl
_padded_gather(w, x, padded_idx) = gather(w, x)
@adjoint function _padded_gather(w, x, padded_idx)
  y, back = pullback(gather, w, x)
  return y, Δ -> begin
    Δ′, _ = back(Δ)
    if !iszero(padded_idx)
      @view(Δ′[:, padded_idx]) .= 0
    end
    return (Δ′, nothing, nothing)
  end
end

(e::FakeTHEmbedding)(x) = _padded_gather(e.weight, x, e.pad_idx)

function load_state!(layer::FakeTHEmbedding, state)
  load_state!(layer.weight, state.weight')
end

function get_state_dict(state, prefix, embedding::FakeTHEmbedding)
  param = Functors.functor(embedding)[1]
  k = :weight
  cprefix = isnothing(prefix) ? String(k) : join((prefix, k), '.')
  state[cprefix] = param.weight'
end

# torch.nn.ModuleList

struct FakeTHModuleList{N, T<:Tuple} <: THModule
  _modules::T
  FakeTHModuleList(ms) = new{length(ms), typeof(ms)}(ms)
end

FakeTHModuleList(ms...) = FakeTHModuleList(ms)
FakeTHModuleList(ms::Vector) = FakeTHModuleList(Tuple(ms))

Functors.functor(::Type{<:FakeTHModuleList}, modulelist) = modulelist._modules, y -> FakeTHModuleList(y)

Base.iterate(modulelist::FakeTHModuleList) = iterate(modulelist._modules)
Base.iterate(modulelist::FakeTHModuleList, i...) = iterate(modulelist._modules, i...)
Base.length(::FakeTHModuleList{N}) where N = N
Base.getindex(modulelist::FakeTHModuleList, i) = modulelist._modules[i]

function get_state_dict(state, prefix, modulelist::FakeTHModuleList)
  param = Functors.functor(modulelist)[1]
  for (i, layer) in enumerate(modulelist)
    cprefix = join((prefix, i-1), '.')
    get_state_dict(state, cprefix, layer)
  end
end

function load_state!(layer::FakeTHModuleList, state)
  for (layerᵢ, stateᵢ) in zip(layer, state)
    load_state!(layerᵢ, stateᵢ)
  end
end

# displaying layer in tree structure

function Base.show(io::IO, x::M; depth=10) where {M<:THModule}
  print_tree(_printnode, io, x, depth)
  io
end

children(m::THModule) = Tuple(pairs(first(Functors.functor(m))))
function children(::THModule, x::Pair{Symbol})
  v = x[2]
  v isa AbstractArray ? () : children(v)
end
children(::THModule, x::Pair{Int}) = children(x[2])

_printnode(io, x) = print(io, Base.typename(typeof(x)))
function _printnode(io::IO, p::Pair)
  print(io, p[1])
  print(io, " => ")
  v = p[2]
  if v isa AbstractArray
    summary(io, p[2])
  else
    print(io, Base.typename(typeof(v)))
  end
end

# transformers.modeling_utils.Conv1D

struct FakeHGFConv1D{W<:AbstractArray, B<:AbstractArray} <: THModule
  weight::W
  bias::B
end

@functor FakeHGFConv1D

(c::FakeHGFConv1D)(x::AbstractMatrix) = c.weight' * x .+ c.bias

function (c::FakeHGFConv1D)(x::AbstractArray)
  old_size = size(x)
  new_size = Base.setindex(old_size, size(c.weight, 2), 1)

  new_x = reshape(x, old_size[1], :)
  y = c(new_x)
  return reshape(y, new_size)
end

# transformers.modeling_utils.SequenceSummary

abstract type AbstractSummaryType end
struct LastSummary <: AbstractSummaryType end
struct FirstSummary <: AbstractSummaryType end
struct MeanSummary <: AbstractSummaryType end
struct CLSindexSummary <: AbstractSummaryType end

struct FakeHGFSequenceSummary{T<:AbstractSummaryType, S, F} <: THModule
  summary::S
  activation::F

  FakeHGFSequenceSummary{T}(s, a) where {T} = new{T, typeof(s), typeof(a)}(s, a)
end

SummaryType(s::Symbol) = (LastSummary, FirstSummary, MeanSummary, CLSindexSummary)[
  findfirst(isequal(s), (:last, :first, :mean, :cls_index))
]

summarytype(s::FakeHGFSequenceSummary{T}) where {T<:AbstractSummaryType} = T

Functors.functor(::Type{<:FakeHGFSequenceSummary}, ss) = (summary = ss.summary, activation = ss.activation), y -> FakeHGFSequenceSummary{summarytype(ss)}(y...)

summary_input(s::FakeHGFSequenceSummary{LastSummary}, x) = x[:, end, :, :]
summary_input(s::FakeHGFSequenceSummary{FirstSummary}, x) = x[:, 1, :, :]
summary_input(s::FakeHGFSequenceSummary{MeanSummary}, x) = mean(x, dims=2)

summary_input(s::FakeHGFSequenceSummary, x, _) = summary_input(s, x)
summary_input(s::FakeHGFSequenceSummary{CLSindexSummary}, x, ::Nothing) = x[:, end, :, :]
summary_input(s::FakeHGFSequenceSummary{CLSindexSummary}, x, cls_index) = gather(x, cls_index)

function (s::FakeHGFSequenceSummary)(hidden_states, cls_index)
  output = summary_input(s, hidden_states, cls_index)
  output = s.summary(output)
  output = s.activation(output)
  return output
end


text:
export ec2download, ec2stream, getXML

"""
  ec2download(bucket,filelist,OUTDIR)

Download files using pmap from S3 to EC2.
# Arguments
- `bucket::String`: S3 bucket to download from.
- `filelist::Array{String}`: Filepaths to download in `bucket`.
- `OUTDIR::String`: The output directory on EC2 instance.

# Keywords
- `v::Int=0`: Verbosity level. Set v = 1 for download progress.
- `XML::Bool=false`: Download StationXML files for request. Downloads StationXML to
    `joinpath(OUTDIR,"XML")`.

"""
function ec2download(
	aws::AWSConfig,
	bucket::String,
	filelist::Array{String},
	OUTDIR::String;
	v::Int=0,
	XML::Bool=false,
)

	# check being run on AWS
	tstart = now()
	!AWS.localhost_is_ec2() && @warn("Running locally. Run on EC2 for maximum performance.")

	println("Starting Download...      $(now())")
	println("Using $(nworkers()) cores...")

	# send files everywhere
	@eval @everywhere filelist=$filelist
	@eval @everywhere aws=$aws

	# create output files
	OUTDIR = expanduser(OUTDIR)
	outfiles = [joinpath(OUTDIR,f) for f in filelist]
	filedir = unique([dirname(f) for f in outfiles])
	for ii = 1:length(filedir)
		if !isdir(filedir[ii])
			mkpath(filedir[ii])
		end
	end

	# get XML
	if XML
		XMLDIR = joinpath(OUTDIR,"XML")
		getXML(aws,bucket,filelist,XMLDIR,v=v)
	end

    # send outfiles everywhere
	@eval @everywhere outfiles=$outfiles
	# do transfer to ec2
    startsize = diskusage(OUTDIR)
    if v > 0
	    pmap(
			s3_file_map,
			fill(aws,length(outfiles)),
            fill(bucket,length(outfiles)),
            filelist,
            outfiles
        )
    else
        pmap(
			s3_get_file,
			fill(aws,length(outfiles)),
            fill(bucket,length(outfiles)),
            filelist,
            outfiles,
        )
    end

	println("Download Complete!        $(now())          ")
	tend = now()
	# check data in directory
	endsize = diskusage(OUTDIR)
	downloadsize = endsize - startsize
	downloadseconds = (tend - tstart).value / 1000
	println("Download took $(Dates.canonicalize(Dates.CompoundPeriod(tend - tstart)))")
	println("Download size $(formatbytes(downloadsize))")
	println("Download rate $(formatbytes(Int(round(downloadsize / downloadseconds))))/s")
	return nothing
end
ec2download(a...;b...) = ec2download(global_aws_config(region="us-west-2"), a...; b...)
ec2download(d::Dict,a...;b...) = ec2download(global_aws_config(region=d[:region]), a...; b...)

"""
  ec2stream(bucket,filelist)

Stream files using pmap from S3 to EC2.
# Arguments
- `bucket::String`: S3 bucket to download from.
- `filelist::Array{String}`: Filepaths to stream from `bucket`.

# Keywords
- `demean::Bool=false`: Demean data after streaming.
- `detrend::Bool=false`: Detrend data after streaming.
- `msr::Bool=false`: Get multi-stage response after streaming.
- `prune::Bool=false`: Prune empty channels after streaming.
- `rr::Bool=false`: Remove instrument response after streaming.
- `taper::Bool=false`: Taper data after streaming.
- `ungap::Bool=false`: Ungap data after streaming.
- `resample::Bool=false`: Resample data after streaming.
- `fs::Float64=0`: New sampling rate.
- `rtype`: Return requested data as `SeisData` or `Array` of `SeisData`. Defaults
    to `SeisData`. Use `Array` to return an `Array` of `SeisData`.
"""
function ec2stream(
	aws::AWSConfig,
	bucket::String,
	filelist::Array{String};
    demean::Bool = false,
    detrend::Bool = false,
    msr::Bool = false,
    prune::Bool = false,
    rr::Bool = false,
    taper::Bool = false,
    ungap::Bool = false,
    unscale::Bool = false,
    resample::Bool = false,
	fs::Real = Float64(0),
    rtype = SeisData,
)

	# check being run on AWS
	!AWS.localhost_is_ec2() && @warn("Running locally. Run on EC2 for maximum performance.")

	# send files everywhere
	@eval @everywhere aws=$aws
    @eval @everywhere bucket=$bucket
	@eval @everywhere filelist=$filelist

	# do transfer to ec2
	Sarray =  pmap(
		s3_get_seed,
		fill(aws,length(filelist)),
        fill(bucket,length(filelist)),
        filelist,
        fill(demean,length(filelist)),
        fill(detrend,length(filelist)),
        fill(msr,length(filelist)),
        fill(prune,length(filelist)),
        fill(rr,length(filelist)),
        fill(taper,length(filelist)),
        fill(ungap,length(filelist)),
        fill(unscale,length(filelist)),
        fill(resample,length(filelist)),
        fill(fs,length(filelist)),
    )
    if rtype == SeisData
	    return SeisData(Sarray...)
    end
    return Sarray
end
ec2stream(a...;b...) = ec2stream(global_aws_config(region="us-west-2"), a...; b...)
ec2stream(d::Dict,a...;b...) = ec2stream(global_aws_config(region=d[:region]), a...; b...)

"""
  getXML(bucket,filelist,XMLDIR)

Download XML files using pmap from S3 to EC2.
# Arguments
- `bucket::String`: S3 bucket to download from.
- `filelist::Array{String}`: Filepaths to download in `bucket`.
- `XMLDIR::String`: The output directory for stationXML files on EC2 instance.

# Keywords
- `v::Int=0`: Verbosity level. Set v = 1 for download progress.
- `getall::Bool=false`: Download all available StationXML files from `scedc-pds`.
"""
function getXML(
	aws::AWSConfig,
	bucket::String,
	filelist::AbstractArray,
	XMLDIR::String;
	v::Int=0,
	getall::Bool=false,
)

	if !isdir(XMLDIR)
		mkpath(XMLDIR)
	end

	# get all files if needed
	prefix = "FDSNstationXML/CI/" # change this if more nets added
	if getall
		filelist = collect(s3_list_keys(bucket,prefix))
		filelist = [replace(f,"_"=>"") for f in filelist]
		filelist = [replace(f,".xml"=>"__") for f in filelist]
	end

	basenames = basename.(filelist)
	nets = [b[1:2] for b in basenames]
	stas = [replace(b[3:7],"_"=>"") for b in basenames]

	# only get response for CI station
	ind = findall(nets .== "CI")
	if length(ind) > 0
		nets = nets[ind]
		stas = stas[ind]
	else
		@warn("No XML for request.")
	end

	# get input/output files names
	infiles = [joinpath(prefix,"$(nets[ii])_$(stas[ii]).xml") for ii = 1:length(nets)]
	infiles = unique(infiles)
	outfiles = [joinpath(XMLDIR,basename(f)) for f in infiles]

	# download files
	if v > 0
	    pmap(
			s3_file_map,
			fill(aws,length(infiles)),
			fill(bucket,length(infiles)),
			infiles,
			outfiles,
        )
    else
		pmap(
			s3_get_file,
			fill(aws,length(infiles)),
			fill(bucket,length(infiles)),
			infiles,
			outfiles,
		)
    end
	return nothing
end
getXML(a...;b...) = getXML(global_aws_config(region="us-west-2"),a...;b...)
getXML(d::Dict,a...;b...) = getXML(global_aws_config(region=d[:region]),a...;b...)

function s3_file_map(aws::AWSConfig,bucket::String,filein::String,fileout::String)
    s3_get_file(aws,bucket, filein, fileout)
    println("Downloading file: $filein       \r")
	return nothing
end
s3_file_map(a...) = s3_file_map(global_aws_config(region="us-west-2"),a...)
s3_file_map(d::Dict,a...) = s3_file_map(global_aws_config(region=d[:region]),a...)

function s3_get_seed(
	aws::AWSConfig,
	bucket::String,
	filein::String,
    demean::Bool,
    detrend::Bool,
    msr::Bool,
    prune::Bool,
    rr::Bool,
    taper::Bool,
    ungap::Bool,
    unscale::Bool,
    resample::Bool,
    fs::Real,
)
    stream = S3.get_object(bucket, filein,aws_config=aws)
	S = parseseed(stream)

	# remove empty channels
	if prune == true
		prune!(S)
	end

	# Get list of channels with sane instrument codes
	CC = get_seis_channels(S)

	if msr == true
		@warn("Getting response not implemented yet.")
	end

	# unscale
	if unscale == true
  		unscale!(S, chans=CC)
	end

	# Demean
	if demean == true
	  demean!(S, chans=CC)
	end

	# Taper
	if taper == true
	  taper!(S, chans=CC)
	end

	# Ungap
	if ungap == true
	  ungap!(S, chans=CC)
	end

	# resample data
	if resample == true && fs != 0
		resample!(S, chans=CC, fs=fs)
	end

	# Remove response
	# need to implement attaching response
	if rr == true
	  @warn("Removing response not implemented yet.")
	end

	return S
end
s3_get_seed(a...;b...) = s3_get_seed(global_aws_config(region="us-west-2"),a...;b...)
s3_get_seed(d::Dict,a...;b...) = s3_get_seed(global_aws_config(region=d[:region]),a...;b...)

function formatbytes(bytes::Real, digits::Int=1)
	units = ["B", "KB", "MB", "GB", "TB","PB"]
	bytes = max(bytes,0)
	pow = Int(floor((bytes > 0 ? log(bytes) : 0) / log(1024)))
	pow = min(pow,length(units))
	powind = pow < length(units) ? pow + 1 : pow
	return string(round(bytes / 1024 ^ pow,digits=digits)) * units[powind]
end

function diskusage(dir)
	s = read(`du -k $dir`, String)
	kb = parse(Int, split(s)[1])
	return 1024 * kb
end

"""
  parseseed(f)

Convert uint8 data to SeisData.
"""
function parseseed(f::AbstractArray)
    S = SeisData()
    SeisIO.SEED.parsemseed!(S,IOBuffer(f),SeisIO.KW.nx_new,SeisIO.KW.nx_add,false,0)
    return S
end


text:
#######################################################################
# Utilities for working with a file system
#######################################################################

# Returns a list of (dPath, dName) for every subdirectory of `dirPath`
# ASSUMPTION: `dirPath` is a path to a directory
subdirPathsWithNames(dirPath :: String) :: Vector{Tuple{String, String}} = begin
    pathsWithNames = map(
        name -> (joinpath(dirPath, name), name),
        readdir(dirPath)
    )
    filter(pn -> isdir(pn[1]), pathsWithNames)
end

max_stars_repo_name: prl-julia/julia-sub

max_issues_repo_name: prl-julia/julia-sub

max_forks_repo_name: prl-julia/julia-sub


text:
using DataFrames

simulationid = :simulation01

DATADIR = "../../data"
OUTDIR = joinpath("../../output", string(simulationid))

### Load data

DATAFILE = joinpath(DATADIR, "apobec1.txt")

df = readtable(DATAFILE, separator = '\t')

sites = DataFrame(chr = ["chr10", "chr11"], coord = [57515985, 100476998])
# sites = DataFrame(chr = ["chr2", "chr16"], coord = [122152902, 84954513])
# sites = DataFrame(chr = ["chrX", "chr4", "chr5"], coord = [109162174, 6395365, 104440840])
nsites = nrow(sites)

df = df[intersect(findin(df[:chr], sites[:chr]), findin(df[:coord], sites[:coord])), :]

data = Dict{Symbol, Any}(:m=>df[:, :Bulk_er],
                         :coverage=>float(array(df[:, 5:2:51])),
                         :rate=>float(array(df[:, 6:2:52])))
data[:edited] = float(data[:rate].*data[:coverage])

cells = Dict{Int, Vector{Int}}()
for i in 1:nsites
  cells[i] = find(data[:coverage][i, :] .!= 0)
end

ncells = Dict{Int, Int}()
for i in 1:nsites
  ncells[i] = length(cells[i])
end


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2021 Varnerlab
# Robert Frederick School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #

"""
    update_is_minimum_flag(data_dictionary, vargs...)

    Updates where the problem is setup as either a min (is_minimum_flag = true) or max (is_minimum_flag = false)
"""
function update_is_minimum_flag(data_dictionary::Dict{String,Any}, vargs...)::Bool
    @warn("Ooops! You've called an autoegenerated dummy implementation that uses the default problem parameters and setup. Replace this implementation with your code to customize")

    # return the default -
    return true
end

"""
    update_species_bounds_array(data_dictionary, vargs...)

    Update the species bounds array with current state or measurements
"""
function update_species_bounds_array(data_dictionary::Dict{String,Any}, vargs...)::Array{Float64,2}
    @warn("Ooops! You've called an autoegenerated dummy implementation that uses the default problem parameters and setup. Replace this implementation with your code to customize")

    # return the default -
    return data_dictionary["species_bounds_array"]
end

"""
    update_flux_bounds_array(data_dictionary, vargs...)
"""
function update_flux_bounds_array(data_dictionary::Dict{String,Any}, vargs...)::Array{Float64,2}
    @warn("Ooops! You've called an autoegenerated dummy implementation that uses the default problem parameters and setup. Replace this implementation with your code to customize")

    # return the default -
    return data_dictionary["flux_bounds_array"]
end

"""
    update_objective_coefficient_array(data_dictionary, vargs...)
"""
function update_objective_coefficient_array(data_dictionary::Dict{String,Any}, vargs...)::Array{Float64,1}
    @warn("Ooops! You've called an autoegenerated dummy implementation that uses the default problem parameters and setup. Replace this implementation with your code to customize")

    # return the default -
    return data_dictionary["objective_coefficient_array"]
end

"""
    update_species_concentration_array(data_dictionary, vargs...)
"""
function update_species_concentration_array(data_dictionary::Dict{String,Any}, vargs...)::(Union{T, Nothing} where T<:Any)
    @warn("Ooops! You've called an autoegenerated dummy implementation that uses the default problem parameters and setup. Replace this implementation with your code to customize")

    # return the default -
    return data_dictionary["species_concentration_array"]
end

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
function isisogram(s)
    replace(s, r"[\s-]" => "") |> lowercase |> allunique
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
using TensorNetworkCodes
using Test

@testset "pauli_are_commuting" begin
    # 5-qubit code stabilizers: XZZXI, IXZZX, XIXZZ, ZXIXZ
    stabilizers = [[1, 3, 3, 1, 0], [0, 1, 3, 3, 1], [1, 0, 1, 3, 3], [3, 1, 0, 1, 3]]
    logical_x = [1, 1, 1, 1, 1]  # XXXXX
    logical_z = [3, 3, 3, 3, 3]  # ZZZZZ
    @test pauli_are_commuting(stabilizers)
    @test pauli_are_commuting(vcat(stabilizers, [logical_x]))
    @test pauli_are_commuting(vcat(stabilizers, [logical_z]))
    @test !pauli_are_commuting([logical_x, logical_z])
    @test pauli_are_commuting([Int[]])  # empty operator
    @test pauli_are_commuting([])  # empty operator list
end

@testset "pauli_are_independent" begin
    # 5-qubit code stabilizers: XZZXI, IXZZX, XIXZZ, ZXIXZ
    stabilizers = [[1, 3, 3, 1, 0], [0, 1, 3, 3, 1], [1, 0, 1, 3, 3], [3, 1, 0, 1, 3]]
    logical_x = [1, 1, 1, 1, 1]  # XXXXX
    logical_z = [3, 3, 3, 3, 3]  # ZZZZZ
    @test pauli_are_independent(stabilizers)
    @test pauli_are_independent(vcat(stabilizers, [logical_x, logical_z]))
    extra_stabilizer = pauli_product.(stabilizers[1], stabilizers[2])
    @test !pauli_are_independent(vcat(stabilizers, [extra_stabilizer]))
    @test pauli_are_independent([Int[]])  # empty operator
    @test pauli_are_independent([])  # empty operator list
end

@testset "pauli_commutation" begin
    # single-qubit
    @test pauli_commutation.(0, [0, 1, 2, 3]) == [0, 0, 0, 0]  # I with IXYZ
    @test pauli_commutation.(1, [0, 1, 2, 3]) == [0, 0, 1, 1]  # X with IXYZ
    @test pauli_commutation.(2, [0, 1, 2, 3]) == [0, 1, 0, 1]  # Y with IXYZ
    @test pauli_commutation.(3, [0, 1, 2, 3]) == [0, 1, 1, 0]  # Z with IXYZ
    # multi-qubit
    # 5-qubit code stabilizers: XZZXI, IXZZX, XIXZZ, ZXIXZ
    stabilizers = [[1, 3, 3, 1, 0], [0, 1, 3, 3, 1], [1, 0, 1, 3, 3], [3, 1, 0, 1, 3]]
    @test pauli_commutation.(Ref([1, 1, 1, 1, 1]), stabilizers) == zeros(4) # XXXXX
    @test pauli_commutation.(Ref([3, 3, 3, 3, 3]), stabilizers) == zeros(4) # ZZZZZ
    @test pauli_commutation.(Ref([1, 1, 0, 0, 0]), stabilizers) == [1, 0, 0, 1] # XXIII
    @test pauli_commutation.(Ref([0, 0, 0, 3, 3]), stabilizers) == [1, 1, 0, 1] # IIIZZ
    @test pauli_commutation.(Ref([0, 0, 2, 0, 0]), stabilizers) == [1, 1, 1, 0] # IIYII
    @test pauli_commutation(Int[], Int[]) == 0  # empty operator
end

@testset "pauli_pow" begin
    @test pauli_pow.(0, [-1, 0, 1, 2, 3]) == [0, 0, 0, 0, 0]  # I to powers
    @test pauli_pow.(1, [-1, 0, 1, 2, 3]) == [1, 0, 1, 0, 1]  # X to powers
    @test pauli_pow.(2, [-1, 0, 1, 2, 3]) == [2, 0, 2, 0, 2]  # Y to powers
    @test pauli_pow.(3, [-1, 0, 1, 2, 3]) == [3, 0, 3, 0, 3]  # Z to powers
end

@testset "pauli_product" begin
    @test pauli_product.(0, [0, 1, 2, 3]) == [0, 1, 2, 3]  # I with IXYZ -> IXYZ
    @test pauli_product.(1, [0, 1, 2, 3]) == [1, 0, 3, 2]  # X with IXYZ -> XIZY
    @test pauli_product.(2, [0, 1, 2, 3]) == [2, 3, 0, 1]  # Y with IXYZ -> YZIX
    @test pauli_product.(3, [0, 1, 2, 3]) == [3, 2, 1, 0]  # Z with IXYZ -> ZYXI
end

@testset "pauli_product" begin
    operators = [[0, 1, 2, 3], [1, 0, 3, 2], [2, 3, 0, 1]]
    operators_copy = deepcopy(operators)
    @test pauli_product(operators) == [3, 2, 1, 0]
    @test operators == operators_copy  # operators not modified
    @test pauli_product([Int[]]) == []  # empty operator
end

@testset "pauli_product_pow" begin
    operators = [[0, 1, 2, 3], [1, 0, 3, 2], [2, 3, 0, 1]]
    operators_copy = deepcopy(operators)
    @test pauli_product_pow(operators, [0, 0, 0]) == [0, 0, 0, 0]
    @test pauli_product_pow(operators, [0, 0, 1]) == [2, 3, 0, 1]
    @test pauli_product_pow(operators, [0, 1, 0]) == [1, 0, 3, 2]
    @test pauli_product_pow(operators, [0, 1, 1]) == [3, 3, 3, 3]
    @test pauli_product_pow(operators, [1, 0, 0]) == [0, 1, 2, 3]
    @test pauli_product_pow(operators, [1, 0, 1]) == [2, 2, 2, 2]
    @test pauli_product_pow(operators, [1, 1, 0]) == [1, 1, 1, 1]
    @test pauli_product_pow(operators, [1, 1, 1]) == [3, 2, 1, 0]
    @test pauli_product_pow(operators, [-1, 3, 1]) == [3, 2, 1, 0]  # non-binary powers
    @test operators == operators_copy  # operators not modified
    @test pauli_product_pow([Int[]], [1]) == []  # empty operator
end

@testset "pauli_rep_change" begin
    @test pauli_rep_change.([0, 1, 2, 3]) == ['I', 'X', 'Y', 'Z']
    @test pauli_rep_change.(['I', 'X', 'Y', 'Z']) == [0, 1, 2, 3]
end

@testset "pauli_weight" begin
    @test pauli_weight([0, 0, 0, 0, 0]) == 0
    @test pauli_weight([0, 0, 1, 2, 0]) == 2
    @test pauli_weight([3, 0, 2, 0, 1]) == 3
    @test pauli_weight([3, 2, 1, 2, 2]) == 5
    @test pauli_weight(Int[]) == 0  # empty operator
end


text:
export Manager, Node
export CUDD_UNIQUE_SLOTS, CUDD_CACHE_SLOTS
export initialize_cudd
export add_var, add_ith_var, add_level_var, add_const, add_constraint
export ref, deref, recursive_deref
export add_constraint, add_restrict, evaluate
export get_value, is_const, is_nonconst, read_index, node_count
export output_dot, output_stats

mutable struct Manager
end

mutable struct Node
end

mutable struct FILE
end

const CUDD_UNIQUE_SLOTS = 256
const CUDD_CACHE_SLOTS  = 262144

function initialize_cudd()
    dd_manager = ccall((:Cudd_Init, libcudd),
        Ptr{Manager}, (Cuint, Cuint, Cuint, Cuint, Csize_t), 0,0,CUDD_UNIQUE_SLOTS,CUDD_CACHE_SLOTS,0) # default in CUDD
    if dd_manager == C_NULL # Could not allocate memory
        throw(OutOfMemoryError())
    end
    return dd_manager
end

function add_var(mgr::Ptr{Manager})
    new_var = ccall((:Cudd_addNewVar, libcudd),
        Ptr{Node}, (Ptr{Manager},), mgr)
    if new_var == C_NULL # Could not allocate memory
        throw(OutOfMemoryError())
    end
    return new_var
end

function add_ith_var(mgr::Ptr{Manager}, i::Int)
    new_var = ccall((:Cudd_addIthVar, libcudd),
        Ptr{Node}, (Ptr{Manager}, Cuint), mgr, i)
    if new_var == C_NULL # Could not allocate memory
        throw(OutOfMemoryError())
    end
    return new_var
end

function add_level_var(mgr::Ptr{Manager}, level::Int)
    new_var = ccall((:Cudd_addNewVarAtLevel, libcudd),
        Ptr{Node}, (Ptr{Manager}, Cuint), mgr, level)
    if new_var == C_NULL
        throw(OutOfMemoryError())
    end
    return new_var
end

function add_const(mgr::Ptr{Manager}, c::Real)
    const_node = ccall((:Cudd_addConst, libcudd),
        Ptr{Node}, (Ptr{Manager}, Cdouble), mgr, c)
    if const_node == C_NULL
        throw(OutOfMemoryError())
    end
    return const_node
end

function add_constraint(mgr::Ptr{Manager}, f::Ptr{Node}, c::Ptr{Node})
    res = ccall((:Cudd_addConstrain, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Node}, Ptr{Node}), mgr, f, c)
    if res == C_NULL
        throw(OutOfMemoryError())
    end
    return res
end

function add_restrict(mgr::Ptr{Manager}, f::Ptr{Node}, c::Ptr{Node})
    res = ccall((:Cudd_addRestrict, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Node}, Ptr{Node}), mgr, f, c)
    if res == C_NULL
        throw(OutOfMemoryError())
    end
    return res
end

# WARNING: the assignment array starts with the variable indexed 0;
# hence, be sure to create the ith variable from index 0.
function evaluate(mgr::Ptr{Manager}, f::Ptr{Node}, assignment::Array{Cint})
    res = ccall((:Cudd_Eval, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Node}, Ptr{Cint}), mgr, f, assignment)
    if res == C_NULL
        throw(OutOfMemoryError())
    end
    return res
end

function ref(n::Ptr{Node})
    ccall((:Cudd_Ref, libcudd), Cvoid, (Ptr{Node},), n)
end

function deref(n::Ptr{Node})
    ccall((:Cudd_Deref, libcudd), Cvoid, (Ptr{Node},), n)
end

function recursive_deref(table::Ptr{Manager}, n::Ptr{Node})
    ccall((:Cudd_RecursiveDeref, libcudd), Cvoid, (Ptr{Manager}, Ptr{Node}), table, n)
end

function get_value(c::Ptr{Node})
    res = ccall((:Cudd_V, libcudd), Cdouble, (Ptr{Node},), c)
    if res == C_NULL
        throw(OutOfMemoryError())
    end
    return res
end

function is_const(n::Ptr{Node})
    return ccall((:Cudd_IsConstant, libcudd), Cint, (Ptr{Node},), n)
end

function is_nonconst(n::Ptr{Node})
    return ccall((:Cudd_IsNonConstant, libcudd), Cint, (Ptr{Node},), n)
end

function read_index(n::Ptr{Node})
    return ccall((:Cudd_NodeReadIndex, libcudd), Cuint, (Ptr{Node},), n)
end

function node_count(mgr::Ptr{Manager})
    return ccall((:Cudd_ReadNodeCount, libcudd), Clong, (Ptr{Manager},), mgr)
end

function output_dot(mgr::Ptr{Manager}, f::Ptr{Node}, filename::String)
    outfile = ccall(:fopen, Ptr{FILE}, (Cstring, Cstring), filename, "w")
    res = ccall((:Cudd_DumpDot, libcudd), Cint, (Ptr{Manager}, Cint, Ref{Ptr{Node}},
        Ptr{Ptr{UInt8}}, Ptr{Ptr{UInt8}}, Ptr{FILE}), mgr, 1, f, C_NULL, C_NULL, outfile)
    ccall(:fclose, Cint, (Ptr{FILE},), outfile)  # here returns 0 suggesting the file is successfully closed
    return res
end

function output_stats(mgr::Ptr{Manager}, filename::String)
    outfile = ccall(:fopen, Ptr{FILE}, (Cstring, Cstring), filename, "w")
    res = ccall((:Cudd_PrintInfo, libcudd), Cint, (Ptr{Manager}, Ptr{FILE}), mgr, outfile)
    ccall(:fclose, Cint, (Ptr{FILE},), outfile)
    return res
end


text:
"""
Revise.jl tracks source code changes and incorporates the changes to a running Julia session.

Revise.jl works behind-the-scenes. To track a package, e.g. `Example`:
```julia
(@v1.6) pkg> dev Example        # make a development copy of the package
[...pkg output omitted...]

julia> using Revise             # this must come before the package under development

julia> using Example

[...develop the package...]     # Revise.jl will automatically update package functionality to match code changes

```

Functions in Revise.jl that may come handy in special circumstances:
- `Revise.track`: track updates to `Base` Julia itself or `Core.Compiler`
- `includet`: load a file and track future changes. Intended for small, quick works
- `entr`: call an additional function whenever code updates
- `revise`: evaluate any changes in `Revise.revision_queue` or every definition in a module
- `Revise.retry`: perform previously-failed revisions. Useful in cases of order-dependent errors
- `Revise.errors`: report the errors represented in `Revise.queue_errors`
"""
module Revise

# We use a code structure where all `using` and `import`
# statements in the package that load anything other than
# a Julia base or stdlib package are located in this file here.
# Nothing else should appear in this file here, apart from
# the `include("packagedef.jl")` statement, which loads what
# we would normally consider the bulk of the package code.
# This somewhat unusual structure is in place to support
# the VS Code extension integration.

using OrderedCollections, CodeTracking, JuliaInterpreter, LoweredCodeUtils

using CodeTracking: PkgFiles, basedir, srcfiles, line_is_decl, basepath
using JuliaInterpreter: whichtt, is_doc_expr, step_expr!, finish_and_return!, get_return,
                        @lookup, moduleof, scopeof, pc_expr, is_quotenode_egal,
                        linetable, codelocs, LineTypes, is_GotoIfNot, isassign, isidentical
using LoweredCodeUtils: next_or_nothing!, trackedheads, structheads, callee_matches

include("packagedef.jl")

end # module


text:
module CayleyMengerDeterminant
using Base: OneTo
using LinearAlgebra
import ArrayInterface
import Static: StaticInt, SOneTo
using StaticArrays
import StaticArrays: Dynamic
import InverseFunctions

### utility functions

"""
    binomial2(x::Union{Int,StaticInt,Dynamic}) -> Union{Int,StaticInt,Dynamic}

Compute the second binomial coefficient of `x`.

This function is a convenience shorthand for `binomial(x, 2)`.

```julia-repl
julia> binomial2(4)
6
```
"""
function binomial2 end

@inline binomial2(x::Int)::Int = binomial(x, 2)

@inline binomial2(::StaticInt{X}) where {X} = StaticInt(binomial2(X))

@inline binomial2(::Dynamic) = Dynamic()

"""
    inverse_binomial2(x::Union{Int,StaticInt,Dynamic) -> Union{Int,StaticInt,Dynamic}

Compute the integer value that `x` is the second binomial coefficient of.

This function is the inverse of `binomial2`.

```julia-repl
julia> inverse_binomial2(6)
4
```
"""
function inverse_binomial2 end

@inline inverse_binomial2(x::Int)::Int = (convert(Int, sqrt(8x + 1)) + 1) ÷ 2

@inline inverse_binomial2(::StaticInt{X}) where {X} = StaticInt(inverse_binomial2(X))

@inline inverse_binomial2(::Dynamic) = Dynamic()

InverseFunctions.inverse(::typeof(binomial2)) = inverse_binomial2

InverseFunctions.inverse(::typeof(inverse_binomial2)) = binomial2

"""
    CayleyMengerDeterminant.index_triangular_nodiag(ixA::Int, ixB::Int) -> Int

For row `ixA` and column `ixB`, compute the linear index into the flat vector storage for a symmetric or
lower-triangular matrix with a zero diagonal.

`ixA` and `ixB` are both 1-based and contiguous, and `ixA` is assumed to be greater than or equal to `ixB`.

```julia-repl
julia> CayleyMengerDeterminant.index_triangular_nodiag(4,2)
5
```
"""
@inline index_triangular_nodiag(ixA::Int, ixB::Int)::Int = binomial2(ixA - 1) + ixB

export binomial2, inverse_binomial2

### main matrix type

"""
    CayleyMengerDistanceMatrix{T,Sz}(simplex_dimensions::Int, square_distances::Vector{T}) where {T<:Real,Sz::Union{Int,Dynamic}}
    CayleyMengerDistanceMatrix(points::Vararg{NTuple{N,T},N+1})::CayleyMengerDistanceMatrix{T,N} where {T<:Real,N::Int}
    CayleyMengerDistanceMatrix(P::AbstractMatrix{T})::CayleyMengerDistanceMatrix{T,Dynamic} where {T<:Real}

The zero-diagonal symmetric matrix of square distances among the points of an `N`-simplex, backed by efficient linear storage.

When providing `points`, it must be N+1 tuples of N values each, which are the coordinates of the points of the `N`-simplex.
When providing `P`, it must be a tall near-square matrix with N+1 rows and N columns, where the rows are the coordinates of the
points of the `N`-simplex.
When providing `simplex_dimensions` and `square_distances` directly, `simplex_dimensions` must be the integer number of dimension
`N` of the `N-simplex, and `square_distances` must be the precomputed backing linear storage of the square distances among the
points of the `N`-simplex and so must have length `binomial2(simplex_dimensions + 1)`.
"""
struct CayleyMengerDistanceMatrix{T<:Real,Sz} <: ArrayInterface.AbstractArray2{T,2}
    "The number of dimensions `N` of the `N`-simplex whose points are being calculated with."
    simplex_dimensions::Int

    "The squared distances of points in the simplex, in natural iteration order as initially given, stored flat and triangular for efficiency."
    square_distances::Vector{T}

    function CayleyMengerDistanceMatrix{T,Sz}(
        simplex_dimensions::Int,
        square_distances::Vector{T},
    ) where {T<:Real,Sz}
        @assert Sz == Dynamic || Sz == simplex_dimensions
        @assert length(square_distances) == binomial2(simplex_dimensions + 1)
        new{T,Sz}(simplex_dimensions, square_distances)
    end
end

CayleyMengerDistanceMatrix(onlyPoint::Tuple{}) =
    CayleyMengerDistanceMatrix{Union{},0}(0, (Union{})[])

function CayleyMengerDistanceMatrix(
    first_point::NTuple{N,T},
    other_points::Vararg{NTuple{N,T},N},
) where {N,T<:Real}
    points = (first_point, other_points...)

    num_distances = binomial2(N + 1)
    sqr_dists = Vector{T}(undef, num_distances)
    for pointIxA = 1:(N+1)
        for pointIxB = 1:(pointIxA-1)
            pos = index_triangular_nodiag(pointIxA, pointIxB)
            sqr_dists[pos] = sum((points[pointIxA] .- points[pointIxB]) .^ 2)
        end
    end

    CayleyMengerDistanceMatrix{T,N}(N, sqr_dists)
end

function CayleyMengerDistanceMatrix(P::AbstractMatrix{T}) where {T<:Real}
    ps, n = size(P)
    @assert (ps == n + 1) "P should be a tall near-square matrix (N+1 by N)"

    sqr_dists = Vector{T}(undef, binomial2(n + 1))
    for pointIxA = 1:(n+1)
        for pointIxB = 1:(pointIxA-1)
            pos = index_triangular_nodiag(pointIxA, pointIxB)
            sqr_dists[pos] = sum((P[pointIxA, :] .- P[pointIxB, :]) .^ 2)
        end
    end

    CayleyMengerDistanceMatrix{T,Dynamic}(n, sqr_dists)
end

export CayleyMengerDistanceMatrix

### equality and hashing

Base.:(==)(A::CayleyMengerDistanceMatrix, B::CayleyMengerDistanceMatrix) =
    (A.simplex_dimensions == B.simplex_dimensions) &&
    all(A.square_distances .== B.square_distances)

Base.isequal(A::CayleyMengerDistanceMatrix, B::CayleyMengerDistanceMatrix) =
    (A.simplex_dimensions == B.simplex_dimensions) &&
    all(isequal.(A.square_distances, B.square_distances))

Base.hash(A::CayleyMengerDistanceMatrix, h::UInt) =
    hash(A.square_distances, hash(A.simplex_dimensions, h))

# conversion and promotion

Base.widen(A::CayleyMengerDistanceMatrix{T,Sz}) where {T<:Real,Sz} =
    CayleyMengerDistanceMatrix{T,Sz}(A.simplex_dimensions, widen.(A.square_distances))

Base.convert(::Type{U}, A::CayleyMengerDistanceMatrix{T,Sz}) where {U<:Real,T<:Real,Sz} =
    CayleyMengerDistanceMatrix{U,Sz}(A.simplex_dimensions, convert.(U, A.square_distances))

"Computes `only` on all non-`Dynamic()` arguments unless all arguments are `Dynamic()`, in which case it returns `Dynamic()`."
only_or_dynamic(vals...) =
    all(vals .== Dynamic()) ? Dynamic() : only([v for v in vals if v != Dynamic()])

Base.promote_rule(
    ::Type{<:CayleyMengerDistanceMatrix{T,SzT}},
    ::Type{<:CayleyMengerDistanceMatrix{U,SzU}},
) where {T<:Real,U<:Real,SzT,SzU} =
    CayleyMengerDistanceMatrix{promote_type(T, U),only_or_dynamic(SzT, SzU)}

### AbstractArray interface

Base.IndexStyle(::Type{<:CayleyMengerDistanceMatrix}) = IndexCartesian()

Base.size(A::CayleyMengerDistanceMatrix) =
    (A.simplex_dimensions + 2, A.simplex_dimensions + 2)

Base.length(A::CayleyMengerDistanceMatrix) = (A.simplex_dimensions + 2)^2

function Base.getindex(A::CayleyMengerDistanceMatrix{T}, i::Int, j::Int) where {T<:Real}
    n = A.simplex_dimensions + 2
    i == j && return zero(T)
    (i == n || j == n) && return oneunit(T)
    return A.square_distances[index_triangular_nodiag(max(i, j), min(i, j))]
end

Base.axes(A::CayleyMengerDistanceMatrix{T}) where {T<:Real} =
    (OneTo(A.simplex_dimensions + 2), OneTo(A.simplex_dimensions + 2))

Base.axes(A::CayleyMengerDistanceMatrix{T}, i::Int) where {T<:Real} =
    OneTo(A.simplex_dimensions + 2)

### LinearAlgebra interface

LinearAlgebra.issymmetric(::CayleyMengerDistanceMatrix) = true

LinearAlgebra.ishermitian(::CayleyMengerDistanceMatrix) = true

### StaticArrays integration

@inline function StaticArrays.Size(A::CayleyMengerDistanceMatrix{T,Sz}) where {T<:Real,Sz}
    N = Sz isa Dynamic ? Dynamic() : Sz + 2
    N, N
end

### ArrayInterface integration

ArrayInterface.can_change_size(::Type{<:CayleyMengerDistanceMatrix}) = false

ArrayInterface.can_setindex(::Type{<:CayleyMengerDistanceMatrix}) = false

ArrayInterface.has_sparsestruct(::CayleyMengerDistanceMatrix) = true

ArrayInterface.ismutable(::CayleyMengerDistanceMatrix) = false

ArrayInterface.isstructured(::CayleyMengerDistanceMatrix) = true

ArrayInterface.known_first(::CayleyMengerDistanceMatrix{T}) where {T<:Real} = zero(T)

ArrayInterface.known_last(::CayleyMengerDistanceMatrix{T}) where {T<:Real} = zero(T)

ArrayInterface.parent(A::CayleyMengerDistanceMatrix{T}) where {T<:Real} =
    convert.(T, A.square_distances)

function ArrayInterface.findstructralnz(A::CayleyMengerDistanceMatrix{T}) where {T<:Real}
    n = A.simplex_dimensions + 1
    v = n * n - n
    I = 1:n
    J = Vector{Int}(undef, n)
    k = 0
    for i = 1:n
        for j = 1:(i-1)
            k += 1
            J[k] = j
        end
    end
    I, J
end

ArrayInterface.getindex(A::CayleyMengerDistanceMatrix, i::Int, j::Int) = A[i, j]

ArrayInterface.unsafe_getindex(A::CayleyMengerDistanceMatrix, i::Int, j::Int) = A[i, j]

ArrayInterface.size(::CayleyMengerDistanceMatrix{T,N}) where {T<:Real,N} =
    (StaticInt(N + 2), StaticInt(N + 2))

ArrayInterface.size(A::CayleyMengerDistanceMatrix{T,Dynamic}) where {T<:Real} =
    (A.simplex_dimensions + 2, A.simplex_dimensions + 2)

ArrayInterface.axes_types(::Type{CayleyMengerDistanceMatrix{T,N}}) where {T<:Real,N} =
    NTuple{2,SOneTo{N + 2}}

ArrayInterface.axes_types(::Type{CayleyMengerDistanceMatrix{T,Dynamic}}) where {T<:Real} =
    NTuple{2,OneTo{Int}}

ArrayInterface.axes(::CayleyMengerDistanceMatrix{T,N}) where {T<:Real,N} =
    (SOneTo{N + 2}(), SOneTo{N + 2}())

ArrayInterface.axes(A::CayleyMengerDistanceMatrix{T,Dynamic}) where {T<:Real} =
    (OneTo(A.simplex_dimensions + 2), OneTo(A.simplex_dimensions + 2))

ArrayInterface.axes(::CayleyMengerDistanceMatrix{T,N}, i) where {T<:Real,N} =
    (SOneTo{N + 2}(), SOneTo{N + 2}())

ArrayInterface.axes(
    A::CayleyMengerDistanceMatrix{T,Dynamic},
    i::Union{Int,StaticInt},
) where {T<:Real} = OneTo(A.simplex_dimensions + 2)

# the operation

"""
    simplex_volume(distances::CayleyMengerDistanceMatrix{T,Sz}; distance_type::Union{Nothing,Type{<:Real}} = nothing) where {N,T<:Real}
    simplex_volume(points::Vararg{NTuple{N,T},N+1}; distance_type::Union{Nothing,Type{<:Real}} = nothing) where {N,T<:Real}
    simplex_volume(P::AbstractMatrix{T}; distance_type::Union{Nothing,Type{<:Real}} = nothing) where {T<:Real}

Calculates the (interior) measure of an `N`-simplex given the coordinates of its points, using the Cayley-Menger
determinant involving the squared distances between the points.

For example, if `N == 2`, then this function calculates the area of a triangle given the 2-dimensional coordinates
of its 3 points; likewise, if `N == 3`, then then this function calculates the volume of a tetrahedron given the
3-dimensional coordinates of its 4 points.

When providing `points`, it must be N+1 tuples of N values each, which are the coordinates of the points of the `N`-simplex.
When providing `P`, it must be a tall near-square matrix with N+1 rows and N columns, where the rows are the coordinates of the
points of the `N`-simplex.

If `distance_type` is provided and is not `nothing`, then the internal calculations on the squared distances will be done using
`distance_type` as the type; otherwise, the type used for the internal calculations will be automatically derived from `T`.
"""
function simplex_volume end

function simplex_volume(
    distances::CayleyMengerDistanceMatrix{T};
    distance_type::Union{Nothing,Type{<:Real}} = nothing,
) where {T<:Real}
    Ty = isnothing(distance_type) ? T : distance_type
    n::Int = distances.simplex_dimensions
    n == 0 && return isnothing(distance_type) ? false : zero(distance_type)
    distances_converted = convert(Ty, distances)
    sqrt(abs(det(distances_converted)) / (2^n)) / factorial(n)
end

function simplex_volume(
    points::NTuple{N,T}...;
    distance_type::Union{Nothing,Type{<:Real}} = nothing,
) where {N,T<:Real}
    distances = CayleyMengerDistanceMatrix(points...)
    simplex_volume(distances; distance_type)
end

simplex_volume(::Tuple{}; distance_type::Union{Nothing,Type{<:Real}} = nothing) =
    isnothing(distance_type) ? false : zero(distance_type)

function simplex_volume(
    P::AbstractMatrix{T};
    distance_type::Union{Nothing,Type{<:Real}} = nothing,
) where {T<:Real}
    distances = CayleyMengerDistanceMatrix(P)
    simplex_volume(distances; distance_type)
end

export simplex_volume

end