text:
@doc raw"""
    CompressibleEulerMulticomponentEquations2D

!!! warning "Experimental code"
    This system of equations is experimental and can change any time.

Multicomponent version of the compressible Euler equations for an ideal gas in two space dimensions.
"""
struct CompressibleEulerMulticomponentEquations2D{RealT<:Real} <: AbstractCompressibleEulerMulticomponentEquations{2, 5}
  gamma1          ::RealT
  gamma2          ::RealT
  gas_constant1   ::RealT
  gas_constant2   ::RealT
  cv1             ::RealT
  cv2             ::RealT
  cp1             ::RealT
  cp2             ::RealT
end

function CompressibleEulerMulticomponentEquations2D()

  # Set ratio of specific heat of the gas per species
  # (It holds: gamma = cp/cv) (Can be used for consistency check!)
  gamma1  = 1.4
  gamma2  = 1.4

  # Set specific gas constant with respect to the molar mass per species
  # (It holds: gas_constant = R/M, with R = molar gas constant, M = molar mass)
  # (For a calorically and thermally perfect gas Mayer's relation holds: gas_constant = cp - cv) (Can be used for consistency check!)
  gas_constant1     = 0.4
  gas_constant2     = 0.4

  # Set specific heat for a constant volume per species
  cv1    = 1.0
  cv2    = 1.0

  # Set specific heat for a constant pressure per species
  cp1    = 1.4
  cp2    = 1.4

  CompressibleEulerMulticomponentEquations2D(gamma1, gamma2, gas_constant1, gas_constant2, cv1, cv2, cp1, cp2)
end

function CompressibleEulerMulticomponentEquations2D(gamma1, gamma2, gas_constant1, gas_constant2)

  # Set specific heat for a constant volume per species
  cv1    = gas_constant1 / (gamma1 - 1.0)
  cv2    = gas_constant2 / (gamma2 - 1.0)

  # Set specific heat for a constant pressure per species
  cp1    = gas_constant1 + cv1
  cp2    = gas_constant2 + cv2

  CompressibleEulerMulticomponentEquations2D(gamma1, gamma2, gas_constant1, gas_constant2, cv1, cv2, cp1, cp2)
end

get_name(::CompressibleEulerMulticomponentEquations2D) = "CompressibleEulerMulticomponentEquations2D"
varnames(::typeof(cons2cons), ::CompressibleEulerMulticomponentEquations2D) = @SVector ["rho1", "rho2", "rho_v1", "rho_v2", "rho_e"]
varnames(::typeof(cons2prim), ::CompressibleEulerMulticomponentEquations2D) = @SVector ["rho1", "rho2", "v1", "v2", "p"]

# Set initial conditions at physical location `x` for time `t`

"""
    initial_condition_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations2D)

A smooth initial condition used for convergence tests in combination with
[`source_terms_convergence_test`](@ref)
(and [`boundary_condition_convergence_test`](@ref) in non-periodic domains).
"""
function initial_condition_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations2D)
  c       = 2
  A       = 0.1
  L       = 2
  f       = 1/L
  omega   = 2 * pi * f
  ini     = c + A * sin(omega * (x[1] + x[2] - t))

  v1      = 1.0
  v2      = 1.0

  rho     = ini
  rho1    = 0.2 * rho
  rho2    = 0.8 * rho
  rho_v1  = rho * v1
  rho_v2  = rho * v2
  rho_e   = rho^2

  return @SVector [rho1, rho2, rho_v1, rho_v2, rho_e]
end

"""
    source_terms_convergence_test(u, x, t, equations::CompressibleEulerMulticomponentEquations2D)

Source terms used for convergence tests in combination with
[`initial_condition_convergence_test`](@ref)
(and [`boundary_condition_convergence_test`](@ref) in non-periodic domains).
"""
@inline function source_terms_convergence_test(u, x, t, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv1, cv2, gamma1, gamma2 = equations
  # Same settings as in `initial_condition`
  c       = 2
  A       = 0.1
  L       = 2
  f       = 1/L
  omega   = 2 * pi * f

  rho1    = u[1]
  rho2    = u[2]

  gamma   = (rho1*cv1*gamma1 + rho2*cv2*gamma2) / (rho1*cv1 + rho2*cv2)

  x1, x2  = x
  si, co  = sincos((x1 + x2 - t)*omega)
  tmp1    = co * A * omega
  tmp2    = si * A
  tmp3    = gamma - 1
  tmp4    = (2*c - 1)*tmp3
  tmp5    = (2*tmp2*gamma - 2*tmp2 + tmp4 + 1)*tmp1
  tmp6    = tmp2 + c

  du1     = 0.2 * tmp1
  du2     = 0.8 * tmp1
  du3     = tmp5
  du4     = tmp5
  du5     = 2*((tmp6 - 1.0)*tmp3 + tmp6*gamma)*tmp1

  # Original terms (without performanc enhancements)
  # du1 = cos((x1 + x2 - t)*ω)*A*ω
  # du2 = (2*sin((x1 + x2 - t)*ω)*A*γ - 2*sin((x1 + x2 - t)*ω)*A +
  #                             2*c*γ - 2*c - γ + 2)*cos((x1 + x2 - t)*ω)*A*ω
  # du3 = (2*sin((x1 + x2 - t)*ω)*A*γ - 2*sin((x1 + x2 - t)*ω)*A +
  #                             2*c*γ - 2*c - γ + 2)*cos((x1 + x2 - t)*ω)*A*ω
  # du3 = 2*((c - 1 + sin((x1 + x2 - t)*ω)*A)*(γ - 1) +
  #                             (sin((x1 + x2 - t)*ω)*A + c)*γ)*cos((x1 + x2 - t)*ω)*A*ω

  return SVector(du1, du2, du3, du4, du5)
end

"""
    boundary_condition_convergence_test(u_inner, orientation, direction, x, t,
                                        surface_flux_function,
                                        equations::CompressibleEulerMulticomponentEquations2D)

Boundary conditions used for convergence tests in combination with
[`initial_condition_convergence_test`](@ref) and [`source_terms_convergence_test`](@ref).
"""
function boundary_condition_convergence_test(u_inner, orientation, direction, x, t,
                                              surface_flux_function,
                                              equations::CompressibleEulerMulticomponentEquations2D)
  u_boundary = initial_condition_convergence_test(x, t, equations)

  # Calculate boundary flux
  if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
  end

  return flux
end

"""
    initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D)

A shock-bubble testcase for multicomponent Euler equations
- Ayoub Gouasmi, Karthik Duraisamy, Scott Murman
  Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations
  [arXiv: 1904.00972](https://arxiv.org/abs/1904.00972)
"""
function initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D)
  # bubble test case, see Agertz et al. https://arxiv.org/pdf/1904.00972
  # other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf
  # typical domain is rectangular, we change it to a square, as Trixi can only do squares
  @unpack gas_constant1, gas_constant2 = equations

  # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving
  delta   = 0.03

  # Region I
  rho1_1  = delta
  rho2_1  = 1.225 * gas_constant1/gas_constant2 - delta
  v1_1    = 0.0
  v2_1    = 0.0
  p_1     = 101325

  # Region II
  rho1_2  = 1.225-delta
  rho2_2  = delta
  v1_2    = 0.0
  v2_2    = 0.0
  p_2     = 101325

  # Region III
  rho1_3  = 1.6861 - delta
  rho2_3  = delta
  v1_3    = -113.5243
  v2_3    = 0.0
  p_3     = 159060

  # Set up Region I & II:
  inicenter = SVector(0.0, 0.0)
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  r = sqrt(x_norm^2 + y_norm^2)

  if (x[1] > 0.50)
    # Set up Region III
    rho1    = rho1_3
    rho2    = rho2_3
    v1      = v1_3
    v2      = v2_3
    p       = p_3
  elseif (r < 0.25)
    # Set up Region I
    rho1    = rho1_1
    rho2    = rho2_1
    v1      = v1_1
    v2      = v2_1
    p       = p_1
  else
    # Set up Region II
    rho1    = rho1_2
    rho2    = rho2_2
    v1      = v1_2
    v2      = v2_2
    p       = p_2
  end

  return prim2cons(SVector(rho1, rho2, v1, v2, p), equations)
end

"""
    initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerMulticomponentEquations2D)

A for multicomponent adapted weak blast wave taken from
- Sebastian Hennemann, Gregor J. Gassner (2020)
  A provably entropy stable subcell shock capturing approach for high order split form DG
  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
"""
function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerMulticomponentEquations2D)
  # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
  # Set up polar coordinates
  inicenter         = SVector(0.0, 0.0)
  x_norm            = x[1] - inicenter[1]
  y_norm            = x[2] - inicenter[2]
  r                 = sqrt(x_norm^2 + y_norm^2)
  phi               = atan(y_norm, x_norm)
  sin_phi, cos_phi  = sincos(phi)

  # Calculate primitive variables.    !!! NOTE: EC flux does not transfer mass across an interface separating two different species
  rho1              = r > 0.5 ? 0.2*1.0 : 0.2*1.1691
  rho2              = r > 0.5 ? 0.8*1.0 : 0.8*1.1691
  rho               = rho1 + rho2
  v1                = r > 0.5 ? 0.0 : 0.1882 * cos_phi
  v2                = r > 0.5 ? 0.0 : 0.1882 * sin_phi
  p                 = r > 0.5 ? 1.0 : 1.245

  return prim2cons(SVector(rho1, rho2, v1, v2, p), equations)
end

# Calculate 1D flux for a single point
@inline function calcflux(u, orientation, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv1, cv2, gamma1, gamma2 = equations
  rho1, rho2, rho_v1, rho_v2, rho_e = u

  rho   = rho1 + rho2
  v1    = rho_v1/rho
  v2    = rho_v2/rho
  gamma = (rho1*cv1*gamma1 + rho2*cv2*gamma2) / (rho1*cv1 + rho2*cv2)
  p     = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2))
  #p     = pressure(u, equations)

  if orientation == 1
    f1  = rho1 * v1
    f2  = rho2 * v1
    f3  = rho_v1 * v1 + p
    f4  = rho_v2 * v1
    f5  = (rho_e + p) * v1
  else
    f1  = rho1 * v2
    f2  = rho2 * v2
    f3  = rho_v1 * v2
    f4  = rho_v2 * v2 + p
    f5  = (rho_e + p) * v2
  end
  return SVector(f1, f2, f3, f4, f5)
end

"""
    flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D)

Entropy conserving two-point flux by
- Ayoub Gouasmi, Karthik Duraisamy (2020)
  "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations""
  arXiv:1904.00972v3 [math.NA] 4 Feb 2020
"""
@inline function flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D)
  # Unpack left and right state
  @unpack gas_constant1, gas_constant2, cv1, cv2 = equations
  rho1_ll, rho2_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho1_rr, rho2_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  # helpful primitive variables
  rho_ll      = rho1_ll + rho2_ll
  rho_rr      = rho1_rr + rho2_rr
  v1_ll       = rho_v1_ll/rho_ll
  v2_ll       = rho_v2_ll/rho_ll
  v1_rr       = rho_v1_rr/rho_rr
  v2_rr       = rho_v2_rr/rho_rr

  # helpful mean values
  rho1_mean   = ln_mean(rho1_ll,rho1_rr)
  rho2_mean   = ln_mean(rho2_ll,rho2_rr)
  rho1_avg    = 0.5 * (rho1_ll + rho1_rr)
  rho2_avg    = 0.5 * (rho2_ll + rho2_rr)
  v1_avg      = 0.5 * (v1_ll + v1_rr)
  v2_avg      = 0.5 * (v2_ll + v2_rr)
  v1_square   = 0.5 * (v1_ll^2 + v1_rr^2)
  v2_square   = 0.5 * (v2_ll^2 + v2_rr^2)
  v_sum       = v1_avg + v2_avg

  # multicomponent specific values
  enth        = rho1_avg * gas_constant1 + rho2_avg * gas_constant2
  T_ll        = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / (rho1_ll * cv1 + rho2_ll * cv2)
  T_rr        = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / (rho1_rr * cv1 + rho2_rr * cv2)
  T           = 0.5 * (1.0/T_ll + 1.0/T_rr)
  T_log       = ln_mean(1.0/T_ll, 1.0/T_rr)

  # Calculate fluxes depending on orientation
  if orientation == 1
    f1 = rho1_mean * v1_avg
    f2 = rho2_mean * v1_avg
    f3 = (f1 + f2) * v1_avg + enth/T
    f4 = (f1 + f2) * v2_avg
    f5 = (f1 * cv1 + f2 * cv2)/T_log - 0.5 * (v1_square + v2_square) * (f1 + f2) + v1_avg * f3 + v2_avg * f4
  else
    f1 = rho1_mean * v2_avg
    f2 = rho2_mean * v2_avg
    f3 = (f1 + f2) * v1_avg
    f4 = (f1 + f2) * v2_avg + enth/T
    f5 = (f1 * cv1 + f2 * cv2)/T_log - 0.5 * (v1_square + v2_square) * (f1 + f2) + v1_avg * f3 + v2_avg * f4
  end

  return SVector(f1, f2, f3, f4, f5)
end

"""
    flux_chandrashekar_stable(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D)

Entropy-Stable two-point flux by
- Ayoub Gouasmi, Karthik Duraisamy (2020)
  "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations""
  arXiv:1904.00972v3 [math.NA] 4 Feb 2020
"""
@inline function flux_chandrashekar_stable(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D)
  # Unpack left and right state
  @unpack gas_constant1, gas_constant2, cv1, cv2, gamma1, gamma2 = equations
  rho1_ll, rho2_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho1_rr, rho2_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  # helpful primitive variables
  rho_ll = rho1_ll + rho2_ll
  rho_rr = rho1_rr + rho2_rr
  v1_ll = rho_v1_ll/rho_ll
  v2_ll = rho_v2_ll/rho_ll
  v1_rr = rho_v1_rr/rho_rr
  v2_rr = rho_v2_rr/rho_rr

  # needed mean values
  rho1_mean = ln_mean(rho1_ll,rho1_rr)
  rho2_mean = ln_mean(rho2_ll,rho2_rr)
  rho1_avg  = 0.5 * (rho1_ll + rho1_rr)
  rho2_avg  = 0.5 * (rho2_ll + rho2_rr)
  rho_avg   = 0.5 * (rho_ll + rho_rr)
  v1_avg    = 0.5 * (v1_ll + v1_rr)
  v2_avg    = 0.5 * (v2_ll + v2_rr)
  v1_square = 0.5 * (v1_ll^2 + v1_rr^2)
  v2_square = 0.5 * (v2_ll^2 + v2_rr^2)
  v_ll_sq   = 0.5 * (v1_ll^2 + v2_ll^2)
  v_rr_sq   = 0.5 * (v1_rr^2 + v2_rr^2)
  v_sum     = v1_avg + v2_avg

  # multicomponent specific values
  enth      = rho1_avg * gas_constant1 + rho2_avg * gas_constant2
  T_ll      = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / (rho1_ll * cv1 + rho2_ll * cv2)
  T_rr      = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / (rho1_rr * cv1 + rho2_rr * cv2)
  T         = 0.5 * (T_ll + T_rr)
  T_log     = ln_mean(1.0/T_ll, 1.0/T_rr)

  gamma     = (rho1_avg*cv1*gamma1 + rho2_avg*cv2*gamma2) / (rho1_avg*cv1 + rho2_avg*cv2)

  # Calculate dissipation operator to be ES
  R_x       = zeros(Float64, (5, 5))
  R_y       = zeros(Float64, (5, 5))
  Lambda_x  = zeros(Float64, (5, 5))
  Lambda_y  = zeros(Float64, (5, 5))
  Tau       = zeros(Float64, (5, 5))
  W         = zeros(Float64, (5))
  Dv_x      = zeros(Float64, (5))
  Dv_y      = zeros(Float64, (5))

  # Help variables
  e1    = equations.cv1 * T
  e2    = equations.cv2 * T
  r     = (rho1_avg / rho_avg) * gas_constant1 + (rho2_avg / rho_avg) * gas_constant2
  a     = sqrt(gamma * T * r)
  k     = 0.5 * (v1_square + v2_square)
  h1    = e1 + gas_constant1 * T
  h2    = e2 + gas_constant2 * T
  h     = (rho1_avg / rho_avg) * h1 + (rho2_avg / rho_avg) * h2
  d1    = h1 - gamma * e1
  d2    = h2 - gamma * e2
  ht    = h + k
  tauh  = sqrt(rho_avg/(gamma * r))

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

  R_x[1,1]  = 1.0
  R_x[1,2]  = 0.0
  R_x[1,3]  = 0.0
  R_x[1,4]  = rho1_avg / rho_avg
  R_x[1,5]  = R_x[1,4]

  R_x[2,1]  = 0.0
  R_x[2,2]  = 1.0
  R_x[2,3]  = 0.0
  R_x[2,4]  = rho2_avg / rho_avg
  R_x[2,5]  = R_x[2,4]

  R_x[3,1]  = v1_avg
  R_x[3,2]  = v1_avg
  R_x[3,3]  = 0.0
  R_x[3,4]  = v1_avg + a
  R_x[3,5]  = v1_avg - a

  R_x[4,1]  = v2_avg
  R_x[4,2]  = v2_avg
  R_x[4,3]  = a
  R_x[4,4]  = v2_avg
  R_x[4,5]  = v2_avg

  R_x[5,1]  = k - (d1/(gamma - 1.0))
  R_x[5,2]  = k - (d2/(gamma - 1.0))
  R_x[5,3]  = a * v2_avg
  R_x[5,4]  = ht + a * v1_avg
  R_x[5,5]  = ht - a * v1_avg

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

  R_y[1,1]  = 1.0
  R_y[1,2]  = 0.0
  R_y[1,3]  = 0.0
  R_y[1,4]  = rho1_avg / rho_avg
  R_y[1,5]  = R_y[1,4]

  R_y[2,1]  = 0.0
  R_y[2,2]  = 1.0
  R_y[2,3]  = 0.0
  R_y[2,4]  = rho2_avg / rho_avg
  R_y[2,5]  = R_y[2,4]

  R_y[3,1]  = v1_avg
  R_y[3,2]  = v1_avg
  R_y[3,3]  = -1.0 * a
  R_y[3,4]  = v1_avg
  R_y[3,5]  = v1_avg

  R_y[4,1]  = v2_avg
  R_y[4,2]  = v2_avg
  R_y[4,3]  = 0.0
  R_y[4,4]  = v2_avg + a
  R_y[4,5]  = v2_avg - a

  R_y[5,1]  = k - (d1/(gamma - 1))
  R_y[5,2]  = k - (d2/(gamma - 1))
  R_y[5,3]  = -1.0 * a * v1_avg
  R_y[5,4]  = ht + a * v2_avg
  R_y[5,5]  = ht - a * v2_avg

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

  Lambda_x[1,1] = abs(v1_avg)
  Lambda_x[2,2] = abs(v1_avg)
  Lambda_x[3,3] = abs(v1_avg)
  Lambda_x[4,4] = abs(v1_avg + a)
  Lambda_x[5,5] = abs(v1_avg - a)

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

  Lambda_y[1,1] = abs(v2_avg)
  Lambda_y[2,2] = abs(v2_avg)
  Lambda_y[3,3] = abs(v2_avg)
  Lambda_y[4,4] = abs(v2_avg + a)
  Lambda_y[5,5] = abs(v2_avg - a)

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

  Tau[1,1]  = tauh * (-1.0 * sqrt((rho1_avg/rho_avg) * (rho2_avg/rho_avg)) * sqrt(gamma * (gas_constant2 / gas_constant1)))
  Tau[1,2]  = tauh * ((rho1_avg / rho_avg) * sqrt(gamma - 1.0))
  Tau[1,3]  = 0.0
  Tau[1,4]  = 0.0
  Tau[1,5]  = 0.0

  Tau[2,1]  = tauh * (sqrt((rho1_avg/rho_avg) * (rho2_avg/rho_avg)) * sqrt(gamma * (gas_constant1 / gas_constant2)))
  Tau[2,2]  = tauh * ((rho2_avg / rho_avg) * sqrt(gamma - 1.0))
  Tau[2,3]  = 0.0
  Tau[2,4]  = 0.0
  Tau[2,5]  = 0.0

  Tau[3,1]  = 0.0
  Tau[3,2]  = 0.0
  Tau[3,3]  = tauh / a
  Tau[3,4]  = 0.0
  Tau[3,5]  = 0.0

  Tau[4,1]  = 0.0
  Tau[4,2]  = 0.0
  Tau[4,3]  = 0.0
  Tau[4,4]  = tauh / sqrt(2)
  Tau[4,5]  = 0.0

  Tau[5,1]  = 0.0
  Tau[5,2]  = 0.0
  Tau[5,3]  = 0.0
  Tau[5,4]  = 0.0
  Tau[5,5]  = tauh / sqrt(2)

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

  w_ll  = cons2entropy(u_ll, equations)
  w_rr  = cons2entropy(u_rr, equations)

  W  = w_rr - w_ll

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

  Dv_x = 0.5 * (R_x * Tau) * Lambda_x * transpose(R_x * Tau) * W
  Dv_y = 0.5 * (R_y * Tau) * Lambda_y * transpose(R_y * Tau) * W

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

  # Calculate fluxes depending on orientation
  if orientation == 1
    f1 = (rho1_mean * v1_avg) - Dv_x[1]
    f2 = (rho2_mean * v1_avg) - Dv_x[2]
    f3 = ((f1 + f2) * v1_avg + enth * T) - Dv_x[3]
    f4 = ((f1 + f2) * v2_avg) - Dv_x[4]
    f5 = ((f1 * cv1 + f2 * cv2)/T_log - 0.5 * (v1_square + v2_square) * (f1 + f2) + v1_avg * f3 + v2_avg * f4) - Dv_x[5]
  else
    f1 = (rho1_mean * v2_avg) - Dv_y[1]
    f2 = (rho2_mean * v2_avg) - Dv_y[2]
    f3 = ((f1 + f2) * v1_avg) - Dv_y[3]
    f4 = ((f1 + f2) * v2_avg + enth * T) - Dv_y[4]
    f5 = ((f1 * cv1 + f2 * cv2)/T_log - 0.5 * (v1_square + v2_square) * (f1 + f2) + v1_avg * f3 + v2_avg * f4) - Dv_y[5]
  end

  return SVector(f1, f2, f3, f4, f5)
end

function flux_lax_friedrichs(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv1, cv2, gamma1, gamma2 = equations
  # Calculate primitive variables and speed of sound
  rho1_ll, rho2_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho1_rr, rho2_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  rho_ll = rho1_ll + rho2_ll
  rho_rr = rho1_rr + rho2_rr

  gamma_ll = (rho1_ll*cv1*gamma1 + rho2_ll*cv2*gamma2)/(rho1_ll*cv1 + rho2_ll*cv2)
  gamma_rr = (rho1_rr*cv1*gamma1 + rho2_rr*cv2*gamma2)/(rho1_rr*cv1 + rho2_rr*cv2)

  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  v_mag_ll = sqrt(v1_ll^2 + v2_ll^2)
  p_ll = (gamma_ll - 1) * (rho_e_ll - 1/2 * rho_ll * v_mag_ll^2)
  c_ll = sqrt(gamma_ll * p_ll / rho_ll)
  v1_rr = rho_v1_rr / rho_rr
  v2_rr = rho_v2_rr / rho_rr
  v_mag_rr = sqrt(v1_rr^2 + v2_rr^2)
  p_rr = (gamma_rr - 1) * (rho_e_rr - 1/2 * rho_rr * v_mag_rr^2)
  c_rr = sqrt(gamma_rr * p_rr / rho_rr)

  # Obtain left and right fluxes
  f_ll = calcflux(u_ll, orientation, equations)
  f_rr = calcflux(u_rr, orientation, equations)

  λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
  f1 = 1/2 * (f_ll[1] + f_rr[1]) - 1/2 * λ_max * (rho1_rr   - rho1_ll)
  f2 = 1/2 * (f_ll[2] + f_rr[2]) - 1/2 * λ_max * (rho2_rr   - rho2_ll)
  f3 = 1/2 * (f_ll[3] + f_rr[3]) - 1/2 * λ_max * (rho_v1_rr - rho_v1_ll)
  f4 = 1/2 * (f_ll[4] + f_rr[4]) - 1/2 * λ_max * (rho_v2_rr - rho_v2_ll)
  f5 = 1/2 * (f_ll[5] + f_rr[5]) - 1/2 * λ_max * (rho_e_rr  - rho_e_ll)

  return SVector(f1, f2, f3, f4, f5)
end

@inline function max_abs_speeds(u, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv1, cv2, gamma1, gamma2 = equations
  rho1, rho2, rho_v1, rho_v2, rho_e = u

  rho   = rho1 + rho2
  v1    = rho_v1 / rho
  v2    = rho_v2 / rho

  gamma = (rho1*cv1*gamma1 + rho2*cv2*gamma2)/(rho1*cv1 + rho2*cv2)
  p     = (gamma - 1) * (rho_e - 1/2 * rho * (v1^2 + v2^2))
  c     = sqrt(gamma * p / rho)

  return abs(v1) + c, abs(v2) + c
end

# Convert conservative variables to primitive
@inline function cons2prim(u, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv1, cv2, gamma1, gamma2 = equations
  rho1, rho2, rho_v1, rho_v2, rho_e = u

  rho   = rho1 + rho2
  v1    = rho_v1 / rho
  v2    = rho_v2 / rho
  gamma = (rho1*cv1*gamma1 + rho2*cv2*gamma2)/(rho1*cv1 + rho2*cv2)
  p     = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2))

  return SVector(rho1, rho2, v1, v2, p)
end

# Convert conservative variables to entropy
@inline function cons2entropy(u, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv1, cv2, gamma1, gamma2, gas_constant1, gas_constant2 = equations
  rho1, rho2, rho_v1, rho_v2, rho_e = u

  rho       = rho1 + rho2
  v1        = rho_v1 / rho
  v2        = rho_v2 / rho
  v_square  = v1^2 + v2^2
  gamma     = (rho1*cv1*gamma1 + rho2*cv2*gamma2)/(rho1*cv1 + rho2*cv2)
  p         = (gamma - 1) * (rho_e - 0.5 * rho * v_square)
  s         = log(p) - gamma*log(rho)
  rho_p     = rho / p

  # Multicomponent stuff
  T         = (rho_e - 0.5 * rho * v_square) / (rho1 * cv1 + rho2 * cv2)
  s1        = cv1 * log(T) - gas_constant1 * log(rho1)
  s2        = cv2 * log(T) - gas_constant2 * log(rho2)

  # Entropy variables
  w1        = -s1 + gas_constant1 + cv1 - (v_square / (2*T)) # + e01/T (bec. e01 = 0 for compressible euler) (DONT confuse it with e1 which is e1 = cv * T)
  w2        = -s2 + gas_constant2 + cv2 - (v_square / (2*T)) # + e02/T (bec. e02 = 0 for compressible euler)
  w3        = (v1)/T
  w4        = (v2)/T
  w5        = (-1.0)/T

  return SVector(w1, w2, w3, w4, w5)
end

# Convert primitive to conservative variables
@inline function prim2cons(prim, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv1, cv2, gamma1, gamma2 = equations
  rho1, rho2, v1, v2, p = prim
  rho     = rho1 + rho2

  gamma  = (rho1*cv1*gamma1 + rho2*cv2*gamma2)/(rho1*cv1 + rho2*cv2)
  rho_v1  = rho * v1
  rho_v2  = rho * v2
  rho_e   = p/(gamma-1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2)

  return SVector(rho1, rho2, rho_v1, rho_v2, rho_e)
end

@inline function density_pressure(u, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv1, cv2, gamma1, gamma2 = equations
 rho1, rho2, rho_v1, rho_v2, rho_e = u

 rho          = rho1 + rho2
 gamma        = (rho1*cv1*gamma1 + rho2*cv2*gamma2)/(rho1*cv1 + rho2*cv2)
 rho_times_p  = (gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2))
 return rho_times_p
end


text:
add_brackets(arr::AbstractArray, vars) = [add_brackets(element, vars) for element in arr]

add_brackets(s::Any, vars) = return s

function add_brackets(ex::Expr, vars)
    ex = postwalk(x -> x in vars ? "\\left[ $(convertSubscript(x)) \\right]" : x, ex)
    return ex
end

"""
    render(::LaTeXString; debug=false, name=tempname(), command="\\Large")

Display a standalone PDF with the given input.
"""
function render(s::LaTeXString; debug=false, name=tempname(), command="\\Large")
    doc = """
    \\documentclass[varwidth]{standalone}
    \\usepackage{amssymb}
    \\usepackage{amsmath}
    \\begin{document}
    {
        $command
        $s
    }
    \\end{document}
    """
    doc = replace(doc, "\\begin{align}"=>"\\[\n\\begin{aligned}")
    doc = replace(doc, "\\end{align}"=>"\\end{aligned}\n\\]")
    open("$(name).tex", "w") do f
        write(f, doc)
    end
    cd(dirname(name)) do
        cmd = `lualatex $(name).tex`
        debug || (cmd = pipeline(cmd, devnull))
        run(cmd)
    end
    if Sys.iswindows()
        run(`cmd /c "start $(name).pdf"`, wait=false)
    elseif Sys.islinux()
        run(`xdg-open $(name).pdf`, wait=false)
    elseif Sys.isapple()
        run(`open $(name).pdf`, wait=false)
    elseif Sys.isbsd()
        run(`xdg-open $(name).pdf`, wait=false)
    end
    return nothing
end

ext: jl


text:
using Knet

"""
Global hyper parameters are defined here for simplicity in
changing their values.
"""
# Image and Data Paths
wf_path         = "/datasets/widerface/WIDER_"
wf_labels_path  = "/datasets/widerface/retinaface/"
wf_eval_path    = "./data/"
fddb_path       = "/datasets/fddb/"
aflw_path       = "/datasets/aflw/"

# Weight Paths
r50_mat_dir     = "./weights/imagenet-resnet-50-dag.mat"

# Image Manipulation Metrics
avg_img         = (123, 117, 104)

# Training Parameters
num_epochs      = 60
weight_decay    = 0.0005
momentum        = 0.9
lrs             = [1e-3, 4e-3, 1e-3, 1e-4]
lr_change_epoch = [10, 20, 30]

# Evaluation Parameters
mAP_ious        = [0.5, 0.05, 0.95] # start, step size, end

# Device Parameters
run_gpu = true
atype = run_gpu ? KnetArray{Float32} : Array{Float32}

# Anchor Metrics
lat_3_anchors   = [
    Dict([("stride", 8), ("anchors", [16, 32])]),
    Dict([("stride", 16), ("anchors", [64, 128])]),
    Dict([("stride", 32), ("anchors", [256, 512])])
]

lat_5_anchors = [
    Dict([("stride", 4), ("anchors", [16, 20.16, 25.40])]),
    Dict([("stride", 8), ("anchors", [32, 40.32, 50.80])]),
    Dict([("stride", 16), ("anchors", [64, 80.63, 101.59])]),
    Dict([("stride", 32), ("anchors", [128, 161.26, 203.19])]),
    Dict([("stride", 64), ("anchors", [256, 322.54, 406.37])])
]

# Loss Calculation Metrics
head1_pos_iou   = 0.7
head1_neg_iou   = 0.3
head2_pos_iou   = 0.5
head2_neg_iou   = 0.4 # 0.3 # --> for no cascaded structure
ohem_ratio      = 3
lambda1         = 2
lambda2         = 1

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


text:
struct ZIBinomialExpert{T<:Real} <: NonZIDiscreteExpert
    p0::T
    n::Int
    p::T
    ZIBinomialExpert{T}(p0, n, p) where {T<:Real} = new{T}(p0, n, p)
end

function ZIBinomialExpert(p0::T, n::Integer, p::T; check_args=true) where {T <: Real}
    check_args && @check_args(ZIBinomialExpert, 0 <= p0 <= 1 && 0 <= p <= 1 && isa(n, Integer))
    return ZIBinomialExpert{T}(p0, n, p)
end

## Outer constructors
# ZIBinomialExpert(p0::Real, n::Integer, p::Real) = ZIBinomialExpert(p0, n, float(p))
ZIBinomialExpert(p0::Integer, n::Integer, p::Integer) = ZIBinomialExpert(float(p0), n, float(p))
ZIBinomialExpert() = ZIBinomialExpert(0.50, 2, 0.50)

## Conversion
function convert(::Type{ZIBinomialExpert{T}}, p0::S, n::Int, p::S) where {T <: Real, S <: Real}
    ZIBinomialExpert(T(p0), n, T(p))
end
function convert(::Type{ZIBinomialExpert{T}}, d::ZIBinomialExpert{S}) where {T <: Real, S <: Real}
    ZIBinomialExpert(T(d.p0), d.n, T(d.p), check_args=false)
end
copy(d::ZIBinomialExpert) = ZIBinomialExpert(d.p0, d.n, d.p, check_args=false)

## Loglikelihood of Expoert
logpdf(d::ZIBinomialExpert, x...) = isinf(x...) ? 0.0 : Distributions.logpdf.(Distributions.Binomial(d.n, d.p), x...)
pdf(d::ZIBinomialExpert, x...) = isinf(x...) ? -Inf : Distributions.pdf.(Distributions.Binomial(d.n, d.p), x...)
logcdf(d::ZIBinomialExpert, x...) = isinf(x...) ? 0.0 : Distributions.logcdf.(Distributions.Binomial(d.n, d.p), x...)
cdf(d::ZIBinomialExpert, x...) = isinf(x...) ? 1.0 : Distributions.cdf.(Distributions.Binomial(d.n, d.p), x...)

## expert_ll, etc
expert_ll_exact(d::ZIBinomialExpert, x::Real) = (x == 0.) ? log.(p_zero(d) + (1-p_zero(d))*exp.(LRMoE.logpdf(d, x))) : log.(1-p_zero(d)) + LRMoE.logpdf(d, x)
function expert_ll(d::ZIBinomialExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_ll_pos = LRMoE.expert_ll(LRMoE.BinomialExpert(d.n, d.p), tl, yl, yu, tu)
    # Deal with zero inflation
    p0 = p_zero(d)
    expert_ll = (yl == 0.) ? log.(p0 + (1-p0)*exp.(expert_ll_pos)) : log.(0.0 + (1-p0)*exp.(expert_ll_pos))
    return expert_ll
end
function expert_tn(d::ZIBinomialExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn_pos = LRMoE.expert_tn(LRMoE.BinomialExpert(d.n, d.p), tl, yl, yu, tu)
    # Deal with zero inflation
    p0 = params(d)[1]
    expert_tn = (tl == 0.) ? log.(p0 + (1-p0)*exp.(expert_tn_pos)) : log.(0.0 + (1-p0)*exp.(expert_tn_pos))
    return expert_tn
end
function expert_tn_bar(d::ZIBinomialExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn_bar_pos = LRMoE.expert_tn_bar(LRMoE.BinomialExpert(d.n, d.p), tl, yl, yu, tu)
    # Deal with zero inflation
    p0 = params(d)[1]
    expert_tn_bar = (tl > 0.) ? log.(p0 + (1-p0)*exp.(expert_tn_bar_pos)) : log.(0.0 + (1-p0)*exp.(expert_tn_bar_pos))
    return expert_tn_bar
end

exposurize_expert(d::ZIBinomialExpert; exposure = 1) = d

## Parameters
params(d::ZIBinomialExpert) = (d.p0, d.n, d.p)
p_zero(d::ZIBinomialExpert) = d.p0
function params_init(y, d::ZIBinomialExpert)
    n_init = Int(maximum(vec(y))) + 2
    μ, σ2 = mean(y), var(y)
    p_init = ( (σ2+μ*μ)/μ - 1 ) / (n_init - 1)
    p0_init = 1 - μ/(n_init*p_init)
    try
        return ZIBinomialExpert(p0_init, n_init, p_init)
    catch;
        return ZIBinomialExpert()
    end
end

## Simululation
sim_expert(d::ZIBinomialExpert) = (1 .- Distributions.rand(Distributions.Bernoulli(d.p0), 1)[1]) .* Distributions.rand(Distributions.Binomial(d.n, d.p), 1)[1]

## penalty
penalty_init(d::ZIBinomialExpert) = []
no_penalty_init(d::ZIBinomialExpert) = []
penalize(d::ZIBinomialExpert, p) = 0.0

## statistics
mean(d::ZIBinomialExpert) = (1-d.p0)*mean(Distributions.Binomial(d.n, d.p))
var(d::ZIBinomialExpert) = (1-d.p0)*var(Distributions.Binomial(d.n, d.p)) + d.p0*(1-d.p0)*(mean(Distributions.Binomial(d.n, d.p)))^2
quantile(d::ZIBinomialExpert, p) = p <= d.p0 ? 0.0 : quantile(Distributions.Binomial(d.n, d.p), p-d.p0)

## EM: M-Step
function EM_M_expert(d::ZIBinomialExpert,
                    tl, yl, yu, tu,
                    exposure,
                    z_e_obs, z_e_lat, k_e;
                    penalty = true, pen_pararms_jk = [2.0 1.0])

    # Old parameters
    p_old = p_zero(d)

    if p_old > 0.999999
        return d
    end

    # Not affected by Exposurize
    tmp_exp = BinomialExpert(d.n, d.p)
    expert_ll_pos = expert_ll.(tmp_exp, tl, yl, yu, tu)
    expert_tn_bar_pos = expert_tn_bar.(tmp_exp, tl, yl, yu, tu)

    # Update zero probability
    z_zero_e_obs = z_e_obs .* EM_E_z_zero_obs(yl, p_old, expert_ll_pos)
    z_pos_e_obs = z_e_obs .- z_zero_e_obs
    z_zero_e_lat = z_e_lat .* EM_E_z_zero_lat(tl, p_old, expert_tn_bar_pos)
    z_pos_e_lat = z_e_lat .- z_zero_e_lat
    p_new = EM_M_zero(z_zero_e_obs, z_pos_e_obs, z_zero_e_lat, z_pos_e_lat, k_e)

    # Update parameters: call its positive part
    tmp_update = EM_M_expert(tmp_exp,
                            tl, yl, yu, tu,
                            exposure,
                            # z_e_obs, z_e_lat, k_e,
                            z_pos_e_obs, z_pos_e_lat, k_e,
                            penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    return ZIBinomialExpert(p_new, tmp_update.n, tmp_update.p)

end

## EM: M-Step, exact observations
function EM_M_expert_exact(d::ZIBinomialExpert,
                    ye, exposure,
                    z_e_obs;
                    penalty = true, pen_pararms_jk = [Inf 1.0 Inf])

    # Old parameters
    p_old = p_zero(d)

    if p_old > 0.999999
        return d
    end

    # Update zero probability
    expert_ll_pos = expert_ll_exact.(LRMoE.BinomialExpert(d.n, d.p), ye)

    z_zero_e_obs = z_e_obs .* EM_E_z_zero_obs(ye, p_old, expert_ll_pos)
    z_pos_e_obs = z_e_obs .- z_zero_e_obs
    p_new = EM_M_zero(z_zero_e_obs, z_pos_e_obs, 0.0, 0.0, 0.0)

    # Update parameters: call its positive part
    tmp_exp = BinomialExpert(d.n, d.p)
    tmp_update = EM_M_expert_exact(tmp_exp,
                            ye, exposure,
                            z_pos_e_obs;
                            penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    return ZIBinomialExpert(p_new, tmp_update.n, tmp_update.p)

end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
"""
    convert(ListRecord, qrec::QueueRecord{E})

Converts a QueueRecord into the corresponding ListRecord.
"""
function Base.convert(::Type{ListRecord{S,D,I}}, qrec::QueueRecord{Entity{S,D,I}}) where {S,D,I}

    frames = Array{RecordFrame}(undef, nframes(qrec))
    states = Array{RecordState{S,I}}(undef, nstates(qrec))
    defs = Dict{I, D}()

    lo = 1
    for (i,pastframe) in enumerate(1-nframes(qrec) : 0)
        frame = qrec[pastframe]

        hi = lo-1
        for entity in frame
            hi += 1
            defs[entity.id] = entity.def
            states[hi] = RecordState{S,I}(entity.state, entity.id)
        end

        frames[i] = RecordFrame(lo, hi)
        lo = hi + 1
    end

    return ListRecord{S,D,I}(get_timestep(qrec), frames, states, defs)
end
Base.convert(::Type{ListRecord}, qrec::QueueRecord{Entity{S,D,I}}) where {S,D,I} = convert(ListRecord{S,D,I}, qrec)

"""
    convert(QueueRecord, lrec::ListRecord)

Converts a ListRecord into the corresponding QueueRecord{Entity{S,D,I}}.
Note that the timesteps for a ListRecord are not necessarily constant timesteps.
"""
function Base.convert(::Type{QueueRecord{Entity{S,D,I}}}, lrec::ListRecord{S,D,I}) where {S,D,I}

    N = nframes(lrec)
    M = maximum(n_objects_in_frame(lrec, i) for i in 1 : N)
    retval = QueueRecord(Entity{S,D,I}, N, get_timestep(lrec), M)

    frame = Frame(Entity{S,D,I}, M)
    for i in 1 : N
        get!(frame, lrec, i)
        update!(retval, frame)
    end

    return retval
end

"""
    get_sparse_lookup(lrec::ListRecord{S,D,I})

Converts a ListRecord into the corresponding sparse matrix containing the states, SparseMatrixCSC{S,Int}.
This requires I <: Integer.

In the sparse array, each column is the state of an entity, and reach row is a frame entry.

A tuple containing the sparse matrix and a dictionary mapping from ids to the column number are returned.
"""
function get_sparse_lookup(rec::ListRecord{S,D,I}) where {S,D,I<:Integer}

    # ids are by time of entry
    # translate id to index on range of 1:n
    id_lookup = Dict{I,Int}(id => index for (index, id) in
        (rec.defs |> keys |> collect |> sort |> enumerate))

    m = nframes(rec) # num rows
    n = nids(rec) # num cols
    n_states = length(rec.states)

    # the row/frame of each state
    Is = Vector{Int}(undef, n_states)
    # the column/car index (not the id, see id_lookup)
    Js = similar(Is)
    # the states themselves
    Vs = similar(Is, S)

    # the index into I, J, and V (on [1:n_states])
    idx = 1
    for (fid, frame) in enumerate(rec.frames)
        for stateid in frame.lo : frame.hi
            recstate = rec.states[stateid]
            Is[idx] = fid
            Js[idx] = id_lookup[recstate.id]
            Vs[idx] = recstate.state

            idx += 1
        end
    end

    sparsemat = sparse(Is, Js, Vs, m, n)
    return (sparsemat, id_lookup)
end

max_stars_repo_name: tawheeler/Records

max_issues_repo_name: tawheeler/Records

max_forks_repo_name: tawheeler/Records


text:
function get_score(stream::IO)
  score = 0
  depth = 0
  is_garbage = false
  garbage_count = 0

  while !eof(stream)
    char = read(stream, Char)

    while char == '!'
      char = read(stream, Char)
      char = read(stream, Char)
    end

    if is_garbage && char != '>'
      garbage_count += 1
    elseif char == '<'
      is_garbage = true
    elseif char == '>'
      is_garbage = false
    elseif char == '{'
      depth += 1
      score += depth
    elseif char == '}'
      depth -= 1
    end
  end

  return score, garbage_count
end

(score, garbage_count) = get_score(STDIN)
println("Score: $score, Garbage count: $garbage_count")


text:
mutable struct ActiveInfo{D} <: InfoDist
    joint::Matrix{Int}
    future::Vector{Int}
    history::Vector{Int}
    B::Int
    bs::NTuple{D,Int}
    k::Int
    N::Int
    function ActiveInfo(bs::NTuple{D, <:Integer}; k::Integer=1) where D
        if any(b -> b < 2, bs)
            throw(ArgumentError("the support of the future must be at least 2"))
        elseif k < 1
            throw(ArgumentError("history length must be at least 1"))
        end
        B = prod(bs)
        new{D}(zeros(Int, B, B^k), zeros(Int, B), zeros(Int, B^k), B, bs, k, 0)
    end
end

ActiveInfo(b::Integer; kwargs...) = ActiveInfo((b,); kwargs...)

function ActiveInfo(xs::AbstractArray{Int,3}; k::Int=1)
    if k < 1
        throw(ArgumentError("history length must be at least 1"))
    elseif length(xs) ≤ k
        throw(ArgumentError("first argument's length must be greater than the history length"))
    end
    xmax = max.(2, maximum(xs; dims=(2,3)))
    observe!(ActiveInfo(tuple(xmax...); k), xs)
end

function ActiveInfo(xs::AbstractArray{Int,2}; k::Int=1)
    if k < 1
        throw(ArgumentError("history length must be at least 1"))
    elseif length(xs) ≤ k
        throw(ArgumentError("first argument's length must be greater than the history length"))
    end
    xmax = max.(2, maximum(xs; dims=2))
    observe!(ActiveInfo(tuple(xmax...); k), xs)
end

function ActiveInfo(xs::AbstractArray{Int,1}; k::Int=1)
    if k < 1
        throw(ArgumentError("history length must be at least 1"))
    elseif length(xs) ≤ k
        throw(ArgumentError("first argument's length must be greater than the history length"))
    end
    xmax = max(2, maximum(xs))
    observe!(ActiveInfo(xmax; k), xs)
end

function estimate(dist::ActiveInfo)
    entropy(dist.future, dist.N) + entropy(dist.history, dist.N) - entropy(dist.joint, dist.N)
end

function observe!(dist::ActiveInfo, xs::AbstractArray{Int,3})
    @views for i in 1:size(xs, 3)
        observe!(dist, xs[:,:,i])
    end
    dist
end

function observe!(dist::ActiveInfo, xs::AbstractArray{Int,2})
    if size(xs, 2) ≤ dist.k
        throw(ArgumentError("data's length must be greater than the history length"))
    end
    dist.N += size(xs, 2) - dist.k
    history, q = 0, 1
    @views for t in 1:dist.k
        q *= dist.B
        history = dist.B * history + index(xs[:,t], dist.bs) - 1
    end
    @views for t in dist.k+1:size(xs, 2)
        x = index(xs[:,t], dist.bs)
        dist.future[x] += 1
        dist.history[history + 1] += 1
        dist.joint[x, history + 1] += 1
        history = dist.B * history - q * (index(xs[:,t-dist.k], dist.bs) - 1) + x - 1
    end
    dist
end

function observe!(dist::ActiveInfo, xs::AbstractArray{Int,1})
    if length(xs) ≤ dist.k
        throw(ArgumentError("data's length must be greater than the history length"))
    end
    dist.N += length(xs) - dist.k
    history, q = 0, 1
    @views for t in 1:dist.k
        q *= dist.bs[1]
        history = dist.bs[1] * history + xs[t] - 1
    end
    @views for t in dist.k+1:length(xs)
        dist.future[xs[t]] += 1
        dist.history[history + 1] += 1
        dist.joint[xs[t], history + 1] += 1
        history = dist.bs[1] * history - q * (xs[t - dist.k] - 1) + xs[t] - 1
    end
    dist
end

@inline function clear!(dist::ActiveInfo)
    dist.joint[:] .= 0
    dist.future[:] .= 0
    dist.history[:] .= 0
    dist.N = 0
    dist
end

function activeinfo!(dist::ActiveInfo, xs::AbstractArray{Int})
    estimate(observe!(dist, xs))
end

activeinfo(xs::AbstractArray{Int}; kwargs...) = estimate(ActiveInfo(xs; kwargs...))

function activeinfo(::Type{Kraskov1}, xs::AbstractMatrix{Float64};
                    k::Int=1, τ::Int=1, nn::Int=1, metric::Metric=Chebyshev())
    hs = history(xs, k, τ, 1)
    fs = @view xs[:, end-size(hs, 2)+1:end]
    mutualinfo(Kraskov1, fs, hs; nn=nn, metric=metric)
end

function activeinfo(::Type{Kraskov1}, xs::AbstractMatrix{Float64},
                    cond::AbstractMatrix{Float64}, conds::AbstractMatrix{Float64}...;
                    k::Int=1, τ::Int=1, nn::Int=1, metric::Metric=Chebyshev())
    hs = history(xs, k, τ, 1)
    N = size(hs, 2)
    fs = @view xs[:, end-N+1:end]
    condview = @view cond[:, end-N:end-1]
    condviews = [@view c[:, end-N:end-1] for c in conds]
    mutualinfo(Kraskov1, fs, hs, condview, condviews...; nn=nn, metric=metric)
end

function activeinfo(::Type{Kraskov}, xs::AbstractMatrix{Float64}, conds::AbstractMatrix{Float64}...; kwargs...)
    activeinfo(Kraskov1, xs, conds...; kwargs...)
end

function activeinfo(xs::AbstractMatrix{Float64}, conds::AbstractMatrix{Float64}...; kwargs...)
    activeinfo(Kraskov1, xs, conds...; kwargs...)
end

function activeinfo(xs::AbstractVector{Float64}, cond::AbstractVector{Float64}...; kwargs...)
    conds = map(c -> reshape(c, 1, length(c)), cond)
    activeinfo(reshape(xs, 1, length(xs)), conds...; kwargs...)
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ b0f0ff46-903b-11eb-3a2f-31e7f1f1ab7a
using ImageDraw, ColorTypes, BenchmarkTools,ColorVectorSpace,ImageMorphology

# ╔═╡ 837217b6-903c-11eb-180f-f14a7e934ff7
begin
	using HistogramThresholding
using TestImages # For the moonsurface image.
using ImageContrastAdjustment # For the build_histogram() function.

img = testimage("moonsurface")
edges, counts = build_histogram(img,256)
t = find_threshold(UnimodalRosin(), counts[1:end], edges)
img_final = zeros(Bool, axes(img))
for i in CartesianIndices(img)
        img_final[i] = img[i] < t ? 0 : 1
end
img_final
end

# ╔═╡ 5bbb0dd4-903e-11eb-0098-9b83c4f34379
using FileIO

# ╔═╡ 5f8fe6da-9040-11eb-086f-53680b92c269
using Noise,ImageShow

# ╔═╡ bb2be408-903b-11eb-3232-99614dda4abc
 begin
vert=CartesianIndex{2}[]
	push!(vert, CartesianIndex(2,2))
	push!(vert, CartesianIndex(3,2))
    push!(vert, CartesianIndex(4,2))
	push!(vert, CartesianIndex(5,2))
    push!(vert, CartesianIndex(6,3))
	push!(vert, CartesianIndex(7,4))
	push!(vert, CartesianIndex(8,4))
	push!(vert, CartesianIndex(9,3))
	push!(vert, CartesianIndex(10,2))
	push!(vert, CartesianIndex(11,3))
	push!(vert, CartesianIndex(12,4))
	push!(vert, CartesianIndex(13,5))
	push!(vert, CartesianIndex(12,6))
	push!(vert, CartesianIndex(11,7))
	push!(vert, CartesianIndex(10,6))#u
	push!(vert, CartesianIndex(9,7))
	push!(vert, CartesianIndex(8,8))
	push!(vert, CartesianIndex(8,9))
	push!(vert, CartesianIndex(8,10))
	push!(vert, CartesianIndex(8,11))
	push!(vert, CartesianIndex(9,11))
	push!(vert, CartesianIndex(10,11))
	push!(vert, CartesianIndex(11,11))
	push!(vert, CartesianIndex(12,11))
	push!(vert, CartesianIndex(13,11))
	push!(vert, CartesianIndex(13,12))
	push!(vert, CartesianIndex(13,13))
	push!(vert, CartesianIndex(12,13))
	push!(vert, CartesianIndex(11,13))
	push!(vert, CartesianIndex(10,13))
	push!(vert, CartesianIndex(9,13))
	push!(vert, CartesianIndex(8,13))
	push!(vert, CartesianIndex(7,13))
	push!(vert, CartesianIndex(6,13))
	push!(vert, CartesianIndex(6,12))
	push!(vert, CartesianIndex(6,11))
	push!(vert, CartesianIndex(5,11))
	push!(vert, CartesianIndex(4,11))
	push!(vert, CartesianIndex(3,11))
	push!(vert, CartesianIndex(2,11))
	push!(vert, CartesianIndex(3,10))
	push!(vert, CartesianIndex(4,9))
	push!(vert, CartesianIndex(5,8))
	push!(vert, CartesianIndex(4,7))
	push!(vert, CartesianIndex(4,6))
	push!(vert, CartesianIndex(5,5))
	push!(vert, CartesianIndex(4,4))
    push!(vert, CartesianIndex(3,3))
	push!(vert, CartesianIndex(2,2))
end

# ╔═╡ a6e92884-903b-11eb-3830-0191de19b831
# begin
# img = zeros(Bool,14,14)
# img[1,1] = 1
# img[13,13] =1
# img[1,13] =1
# img[13,1] =1
# end

# ╔═╡ 8554a5be-903d-11eb-069c-11673c9c9356
img

# ╔═╡ e67457e4-903b-11eb-3b1b-89259950a651
convex_hull_test = convexhull(img_final)

# ╔═╡ 926e9572-903d-11eb-2ec2-e15a01a9524a
# begin
# 	img_final1 = zeros(Gray{Bool}, axes(img))
# for i in CartesianIndices(convex_hull_test)
#         img_final1[i] = convex_hull_test[i][1] < 1 ? 0 : 1
# end
# img_final1
# Gray.(img_final1)
# end

# ╔═╡ 47fe6612-903e-11eb-280d-7f745e6c0612
begin
img_coins = load("../1.jpg")
end

# ╔═╡ a60870b6-903e-11eb-263c-05a406696259
begin
edges1, counts1 = build_histogram(Gray.(img_coins),256)
t1 = find_threshold(UnimodalRosin(), counts1[1:end], edges1)
img_final_coins = zeros(Bool, axes(Gray.(img_coins)))
for i in CartesianIndices(img_coins)
        img_final_coins[i] = Gray.(img_coins[i]) < t1 ? 0 : 1
end
img_final_coins
end

# ╔═╡ 9cd44928-9040-11eb-1f50-8fdfbdb00ab8
Gray.(img_final_coins)

# ╔═╡ 6b951ff6-9040-11eb-287e-6dec9dc1d7fa
img_final_coins1 = salt_pepper(img_final_coins, 0.5)

# ╔═╡ d903f85e-903f-11eb-36ef-d5435d6ed4c6
Gray.(img_final_coins1)

# ╔═╡ 865d4368-903e-11eb-36fc-3745cdadf810
begin
fill_image=imfill(img_final_coins, (1,2))
Gray.(fill_image)
end

# ╔═╡ 24ebf4fa-9041-11eb-245a-31e91bd0f2ee
img_thinning=thinning(img_final_coins,algo=GuoAlgo())

# ╔═╡ 34e54dc0-9041-11eb-325e-43f60e27a486
Gray.(img_thinning)

# ╔═╡ 7236322a-9041-11eb-13f2-1153e9fe7571
Gray.([img_final_coins img_thinning])

# ╔═╡ 007b900c-9042-11eb-27da-611a1b40a7a8
cleared_img = clearborder(img_final_coins, 20)

# ╔═╡ 0b500786-9042-11eb-2b74-79c766338419
Gray.([img_final_coins cleared_img])

# ╔═╡ Cell order:
# ╠═b0f0ff46-903b-11eb-3a2f-31e7f1f1ab7a
# ╠═bb2be408-903b-11eb-3232-99614dda4abc
# ╠═a6e92884-903b-11eb-3830-0191de19b831
# ╠═8554a5be-903d-11eb-069c-11673c9c9356
# ╠═837217b6-903c-11eb-180f-f14a7e934ff7
# ╠═e67457e4-903b-11eb-3b1b-89259950a651
# ╠═926e9572-903d-11eb-2ec2-e15a01a9524a
# ╠═5bbb0dd4-903e-11eb-0098-9b83c4f34379
# ╠═5f8fe6da-9040-11eb-086f-53680b92c269
# ╠═47fe6612-903e-11eb-280d-7f745e6c0612
# ╠═a60870b6-903e-11eb-263c-05a406696259
# ╠═9cd44928-9040-11eb-1f50-8fdfbdb00ab8
# ╠═6b951ff6-9040-11eb-287e-6dec9dc1d7fa
# ╠═d903f85e-903f-11eb-36ef-d5435d6ed4c6
# ╠═865d4368-903e-11eb-36fc-3745cdadf810
# ╠═24ebf4fa-9041-11eb-245a-31e91bd0f2ee
# ╠═34e54dc0-9041-11eb-325e-43f60e27a486
# ╠═7236322a-9041-11eb-13f2-1153e9fe7571
# ╠═007b900c-9042-11eb-27da-611a1b40a7a8
# ╠═0b500786-9042-11eb-2b74-79c766338419

max_stars_repo_name: ashwani-rathee/Demos

max_issues_repo_name: ashwani-rathee/Demos

max_forks_repo_name: ashwani-rathee/Demos


text:
mutable struct Simulation{T <: SimulationModel}
    status::BUILD_STATUS
    problem::Union{Nothing, SciMLBase.DEProblem}
    tspan::NTuple{2, Float64}
    sys::PSY.System
    perturbations::Vector{<:Perturbation}
    x0_init::Vector{Float64}
    initialized::Bool
    tstops::Vector{Float64}
    callbacks::Vector
    simulation_folder::String
    inputs::Union{Nothing, SimulationInputs}
    results::Union{Nothing, SimulationResults}
    console_level::Base.CoreLogging.LogLevel
    file_level::Base.CoreLogging.LogLevel
    multimachine::Bool
end

get_system(sim::Simulation) = sim.sys
get_simulation_inputs(sim::Simulation) = sim.inputs
get_initial_conditions(sim::Simulation) = sim.x0_init
get_tspan(sim::Simulation) = sim.tspan

function Simulation(
    ::Type{T},
    sys::PSY.System;
    tspan,
    initial_conditions,
    initialize_simulation,
    perturbations,
    simulation_folder,
    console_level,
    file_level,
) where {T <: SimulationModel}
    PSY.set_units_base_system!(sys, "DEVICE_BASE")

    return Simulation{T}(
        BUILD_INCOMPLETE,
        nothing,
        tspan,
        sys,
        perturbations,
        initial_conditions,
        !initialize_simulation,
        Vector{Float64}(),
        Vector{SciMLBase.AbstractDiscreteCallback}(),
        simulation_folder,
        nothing,
        nothing,
        console_level,
        file_level,
        false,
    )
end

function Simulation!(
    ::Type{T},
    system::PSY.System,
    simulation_folder::String,
    tspan::NTuple{2, Float64},
    perturbation::Perturbation;
    kwargs...,
) where {T <: SimulationModel}
    return Simulation!(T, system, simulation_folder, tspan, [perturbation]; kwargs...)
end

function Simulation(
    ::Type{T},
    system::PSY.System,
    simulation_folder::String,
    tspan::NTuple{2, Float64},
    perturbation::Perturbation;
    kwargs...,
) where {T <: SimulationModel}
    return Simulation(T, system, simulation_folder, tspan, [perturbation]; kwargs...)
end

"""
    function Simulation!
        ::SimulationModel
        system::PowerSystems.System
        simulation_folder::String
        tspan::NTuple{2, Float64},
        perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
        kwargs...,
    end

Builds the simulation object and conducts the indexing process. The original system is not modified and a copy its created and stored in the Simulation.

# Arguments:
- `::SimulationModel` : Type of Simulation Model. `ResidualModel` or `MassMatrixModel`. See [Models Section](https://nrel-siip.github.io/PowerSimulationsDynamics.jl/stable/models/) for more details
- `system::PowerSystems.System` : System data
- `simulation_folder::String` : Folder directory
- `tspan::NTuple{2, Float64}` : Time span for simulation
- `perturbations::Vector{<:Perturbation}` : Vector of Perturbations for the Simulation. Default: No Perturbations
- `initialize_simulation::Bool` : Runs the initialization routine. If false, simulation runs based on the operation point stored in System
- `initial_conditions::Vector{Float64}` : Allows the user to pass a vector with the initial condition values desired in the simulation. If initialize_simulation = true, these values are used as a first guess and overwritten.
- `frequency_reference` : Default `ReferenceBus`. Determines which frequency model is used for the network. Currently there are two options available:
    - `ConstantFrequency` assumes that the network frequency is 1.0 per unit at all times.
    - `ReferenceBus` will use the frequency state of a Dynamic Generator (rotor speed) or Dynamic Inverter (virtual speed) connected to the Reference Bus (defined in the Power Flow data) as the network frequency. If multiple devices are connected to such bus, the device with larger base power will be used as a reference. If a Voltage Source is connected to the Reference Bus, then a `ConstantFrequency` model will be used.
- `system_to_file::Bool` : Default `false`. Serializes the initialized system
- `console_level::Logging` : Default `Logging.Warn`. Sets the level of logging output to the console. Can be set to `Logging.Error`, `Logging.Warn`, `Logging.Info` or `Logging.Debug`
- `file_level::Logging` : Default `Logging.Info`. Sets the level of logging output to file. Can be set to `Logging.Error`, `Logging.Warn`, `Logging.Info` or `Logging.Debug`
- `disable_timer_output::Bool` : Default `false`. Allows the user to display timer information about the construction and initilization of the Simulation.
"""
function Simulation!(
    ::Type{T},
    system::PSY.System,
    simulation_folder::String,
    tspan::NTuple{2, Float64},
    perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
    kwargs...,
) where {T <: SimulationModel}
    check_folder(simulation_folder)
    # Instantiates the Simulation object
    sim = Simulation(
        T,
        system;
        tspan = tspan,
        initial_conditions = get(kwargs, :initial_conditions, Vector{Float64}()),
        initialize_simulation = get(kwargs, :initialize_simulation, true),
        simulation_folder = simulation_folder,
        perturbations = perturbations,
        console_level = get(kwargs, :console_level, Logging.Warn),
        file_level = get(kwargs, :file_level, Logging.Info),
    )

    build!(sim; kwargs...)
    if get(kwargs, :system_to_file, false)
        PSY.to_json(system, joinpath(simulation_folder, "initialized_system.json"))
    end
    return sim
end

"""
    function Simulation
        ::SimulationModel
        system::PowerSystems.System
        simulation_folder::String
        tspan::NTuple{2, Float64},
        perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
        kwargs...,
    end

Builds the simulation object and conducts the indexing process. The initial conditions are stored in the system.

# Arguments:
- `::SimulationModel` : Type of Simulation Model. `ResidualModel` or `MassMatrixModel`. See [Models Section](https://nrel-siip.github.io/PowerSimulationsDynamics.jl/stable/models/) for more details
- `system::PowerSystems.System` : System data
- `simulation_folder::String` : Folder directory
- `tspan::NTuple{2, Float64}` : Time span for simulation
- `perturbations::Vector{<:Perturbation}` : Vector of Perturbations for the Simulation. Default: No Perturbations
- `initialize_simulation::Bool` : Runs the initialization routine. If false, simulation runs based on the operation point stored in System
- `initial_conditions::Vector{Float64}` : Allows the user to pass a vector with the initial condition values desired in the simulation. If initialize_simulation = true, these values are used as a first guess and overwritten.
- `frequency_reference` : Default `ReferenceBus`. Determines which frequency model is used for the network. Currently there are two options available:
    - `ConstantFrequency` assumes that the network frequency is 1.0 per unit at all times.
    - `ReferenceBus` will use the frequency state of a Dynamic Generator (rotor speed) or Dynamic Inverter (virtual speed) connected to the Reference Bus (defined in the Power Flow data) as the network frequency. If multiple devices are connected to such bus, the device with larger base power will be used as a reference. If a Voltage Source is connected to the Reference Bus, then a `ConstantFrequency` model will be used.
- `system_to_file::Bool` : Default `false`. Serializes the initialized system
- `console_level::Logging` : Default `Logging.Warn`. Sets the level of logging output to the console. Can be set to `Logging.Error`, `Logging.Warn`, `Logging.Info` or `Logging.Debug`
- `file_level::Logging` : Default `Logging.Info`. Sets the level of logging output to file. Can be set to `Logging.Error`, `Logging.Warn`, `Logging.Info` or `Logging.Debug`
- `disable_timer_output::Bool` : Default `false`. Allows the user to display timer information about the construction and initilization of the Simulation.
"""
function Simulation(
    ::Type{T},
    system::PSY.System,
    simulation_folder::String,
    tspan::NTuple{2, Float64},
    perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
    kwargs...,
) where {T <: SimulationModel}
    check_folder(simulation_folder)
    simulation_system = deepcopy(system)
    # Instantiates the Simulation object
    sim = Simulation(
        T,
        simulation_system;
        tspan = tspan,
        initial_conditions = get(kwargs, :initial_conditions, Vector{Float64}()),
        initialize_simulation = get(kwargs, :initialize_simulation, true),
        simulation_folder = simulation_folder,
        perturbations = perturbations,
        console_level = get(kwargs, :console_level, Logging.Warn),
        file_level = get(kwargs, :file_level, Logging.Info),
    )
    build!(sim; kwargs...)
    if get(kwargs, :system_to_file, false)
        PSY.to_json(system, joinpath(simulation_folder, "input_system.json"))
    end
    return sim
end

function reset!(sim::Simulation{T}) where {T <: SimulationModel}
    @info "Rebuilding the simulation after reset"
    sim.inputs = SimulationInputs(T(), get_system(sim), sim.inputs.tspan)
    build!(sim)
    @info "Simulation reset to status $(sim.status)"
    return
end

function configure_logging(sim::Simulation, file_mode; kwargs...)
    return IS.configure_logging(
        console = true,
        console_stream = stderr,
        console_level = get(kwargs, :console_level, sim.console_level),
        file = true,
        filename = joinpath(sim.simulation_folder, SIMULATION_LOG_FILENAME),
        file_level = get(kwargs, :file_level, sim.file_level),
        file_mode = file_mode,
        tracker = nothing,
        set_global = false,
        progress = _prog_meter_enabled(),
    )
end

function _build_inputs!(
    sim::Simulation{T},
    frequency_reference,
) where {T <: SimulationModel}
    simulation_system = get_system(sim)
    sim.inputs = SimulationInputs(T, simulation_system, frequency_reference)
    @debug "Simulation Inputs Created"
    return
end

function _get_flat_start(inputs::SimulationInputs)
    bus_count = get_bus_count(inputs)
    var_count = get_variable_count(inputs)
    initial_conditions = zeros(var_count)
    initial_conditions[1:bus_count] .= 1.0
    return initial_conditions
end

function _initialize_state_space(sim::Simulation{T}) where {T <: SimulationModel}
    simulation_inputs = get_simulation_inputs(sim)
    x0_init = sim.x0_init
    if isempty(x0_init) && sim.initialized
        @warn "Initial Conditions set to flat start"
        sim.x0_init = _get_flat_start(simulation_inputs)
    elseif isempty(x0_init) && !sim.initialized
        sim.x0_init = _get_flat_start(simulation_inputs)
    elseif !isempty(x0_init) && sim.initialized
        if length(sim.x0_init) != get_variable_count(simulation_inputs)
            throw(
                IS.ConflictingInputsError(
                    "The size of the provided initial state space does not match the model's state space.",
                ),
            )
        end
    elseif !isempty(x0_init) && !sim.initialized
        @warn "initial_conditions were provided with initialize_simulation. User's initial_conditions will be overwritten."
        if length(sim.x0_init) != get_variable_count(simulation_inputs)
            sim.x0_init = _get_flat_start(simulation_inputs)
        end
    else
        @assert false
    end
    return
end

function _pre_initialize_simulation!(sim::Simulation)
    _initialize_state_space(sim)
    if sim.initialized != true
        @info("Pre-Initializing Simulation States")
        sim.initialized = precalculate_initial_conditions!(sim.x0_init, sim)
        if !sim.initialized
            error(
                "The simulation failed to find an adequate initial guess for the initialization. Check the intialization routine.",
            )
        end
    else
        @warn(
            "No Pre-initialization conducted. If this is unexpected, check the initialization keywords"
        )
        sim.status = SIMULATION_INITIALIZED
    end
    return
end

function _get_jacobian(sim::Simulation{ResidualModel})
    inputs = get_simulation_inputs(sim)
    x0_init = get_initial_conditions(sim)
    return JacobianFunctionWrapper(
        ResidualModel(inputs, x0_init, JacobianCache),
        x0_init,
        # sparse_retrieve_loop = 0,
    )
end

function _get_jacobian(sim::Simulation{MassMatrixModel})
    inputs = get_simulation_inputs(sim)
    x0_init = get_initial_conditions(sim)
    return JacobianFunctionWrapper(MassMatrixModel(inputs, x0_init, JacobianCache), x0_init)
end

function _build_perturbations!(sim::Simulation)
    @info "Attaching Perturbations"
    if isempty(sim.perturbations)
        @debug "The simulation has no perturbations"
        return SciMLBase.CallbackSet(), [0.0]
    end
    inputs = get_simulation_inputs(sim)
    perturbations = sim.perturbations
    perturbations_count = length(perturbations)
    callback_vector = Vector{SciMLBase.DiscreteCallback}(undef, perturbations_count)
    tstops = Vector{Float64}(undef, perturbations_count)
    for (ix, pert) in enumerate(perturbations)
        @debug pert
        condition = (x, t, integrator) -> t in [pert.time]
        affect = get_affect(inputs, get_system(sim), pert)
        callback_vector[ix] = SciMLBase.DiscreteCallback(condition, affect)
        tstops[ix] = pert.time
    end
    sim.tstops = tstops
    sim.callbacks = callback_vector
    return
end

function _get_diffeq_problem(
    sim::Simulation,
    model::SystemModel{ResidualModel},
    jacobian::JacobianFunctionWrapper,
)
    x0 = get_initial_conditions(sim)
    dx0 = zeros(length(x0))
    simulation_inputs = get_simulation_inputs(sim)
    sim.problem = SciMLBase.DAEProblem(
        SciMLBase.DAEFunction{true}(
            model;
            jac = jacobian,
            tgrad = (dT, u, p, t) -> dT .= false,
            jac_prototype = jacobian.Jv,
        ),
        dx0,
        x0,
        get_tspan(sim),
        simulation_inputs,
        differential_vars = get_DAE_vector(simulation_inputs),
    )
    sim.status = BUILT
    return
end

function _get_diffeq_problem(
    sim::Simulation,
    model::SystemModel{MassMatrixModel},
    jacobian::JacobianFunctionWrapper,
)
    simulation_inputs = get_simulation_inputs(sim)
    sim.problem = SciMLBase.ODEProblem(
        SciMLBase.ODEFunction{true}(
            model,
            mass_matrix = get_mass_matrix(simulation_inputs),
            jac = jacobian,
            jac_prototype = jacobian.Jv,
            # Necessary to avoid unnecessary calculations in Rosenbrock methods
            tgrad = (dT, u, p, t) -> dT .= false,
        ),
        sim.x0_init,
        get_tspan(sim),
        simulation_inputs,
    )
    sim.status = BUILT
    return
end

function _build!(sim::Simulation{T}; kwargs...) where {T <: SimulationModel}
    check_kwargs(kwargs, SIMULATION_ACCEPTED_KWARGS, "Simulation")
    # Branches are a super set of Lines. Passing both kwargs will
    # be redundant.
    TimerOutputs.reset_timer!(BUILD_TIMER)
    if get(kwargs, :disable_timer_outputs, false)
        TimerOutputs.disable_timer!(BUILD_PROBLEMS_TIMER)
    end

    TimerOutputs.@timeit BUILD_TIMER "Build Simulation" begin
        try
            if get(kwargs, :all_branches_dynamic, false)
                TimerOutputs.@timeit BUILD_TIMER "AC Branch Transform to Dynamic" begin
                    sys = get_system(sim)
                    transform_branches_to_dynamic(sys, PSY.ACBranch)
                end
            elseif get(kwargs, :all_lines_dynamic, false)
                TimerOutputs.@timeit BUILD_TIMER "Line Transform to Dynamic" begin
                    sys = get_system(sim)
                    transform_branches_to_dynamic(sys, PSY.Line)
                end
            end
            TimerOutputs.@timeit BUILD_TIMER "Build Simulation Inputs" begin
                f_ref = get(kwargs, :frequency_reference, ReferenceBus)
                _build_inputs!(sim, f_ref)
                # TODO: Update and store f_ref somewhere.
                sim.multimachine =
                    get_global_vars_update_pointers(sim.inputs)[GLOBAL_VAR_SYS_FREQ_INDEX] !=
                    0
            end
            TimerOutputs.@timeit BUILD_TIMER "Pre-initialization" begin
                _pre_initialize_simulation!(sim)
            end
            if sim.status != BUILD_FAILED
                simulation_inputs = get_simulation_inputs(sim)
                TimerOutputs.@timeit BUILD_TIMER "Calculate Jacobian" begin
                    jacobian = _get_jacobian(sim)
                end
                TimerOutputs.@timeit BUILD_TIMER "Make Model Function" begin
                    model = T(simulation_inputs, get_initial_conditions(sim), SimCache)
                end
                TimerOutputs.@timeit BUILD_TIMER "Initial Condition NLsolve refinement" begin
                    refine_initial_condition!(sim, model, jacobian)
                end
                TimerOutputs.@timeit BUILD_TIMER "Build Perturbations" begin
                    _build_perturbations!(sim)
                end
                TimerOutputs.@timeit BUILD_TIMER "Make DiffEq Problem" begin
                    _get_diffeq_problem(sim, model, jacobian)
                end
                @info "Simulations status = $(sim.status)"
            else
                @error "The simulation couldn't be initialized correctly. Simulations status = $(sim.status)"
            end
        catch e
            bt = catch_backtrace()
            @error "$T failed to build" exception = e, bt
            sim.status = BUILD_FAILED
        end
    end
    return
end

function build!(sim; kwargs...)
    logger = configure_logging(sim, "w")
    Logging.with_logger(logger) do
        _build!(sim; kwargs...)
        if sim.status == BUILT
            string_buffer = IOBuffer()
            TimerOutputs.print_timer(
                string_buffer,
                BUILD_TIMER,
                sortby = :firstexec,
                compact = true,
            )
            @info "\n$(String(take!(string_buffer)))\n"
        end
    end
    close(logger)
    return sim.status
end

function simulation_pre_step!(sim::Simulation, reset_sim::Bool)
    if sim.status == BUILD_FAILED
        error(
            "The Simulation status is $(sim.status). Can not continue, correct your inputs and build the simulation again.",
        )
    elseif sim.status != BUILT && !reset_sim
        error(
            "The Simulation status is $(sim.status). Use keyword argument reset_simulation = true",
        )
    end

    reset_sim && reset!(sim)
    return
end

function _prog_meter_enabled()
    return isa(stderr, Base.TTY) &&
           (get(ENV, "CI", nothing) != "true") &&
           (get(ENV, "RUNNING_PSID_TESTS", nothing) != "true")
end

function _execute!(sim::Simulation, solver; kwargs...)
    @debug "status before execute" sim.status
    simulation_pre_step!(sim, get(kwargs, :reset_simulation, false))
    sim.status = SIMULATION_STARTED
    time_log = Dict{Symbol, Any}()
    if get(kwargs, :auto_abstol, false)
        cb = AutoAbstol(true, get(kwargs, :abstol, 1e-9))
        callbacks = SciMLBase.CallbackSet((), tuple(push!(sim.callbacks, cb)...))
    else
        callbacks = SciMLBase.CallbackSet((), tuple(sim.callbacks...))
    end

    solution,
    time_log[:timed_solve_time],
    time_log[:solve_bytes_alloc],
    time_log[:sec_in_gc] = @timed SciMLBase.solve(
        sim.problem,
        solver;
        callback = callbacks,
        tstops = !isempty(sim.tstops) ? [sim.tstops[1] ÷ 2, sim.tstops...] : [],
        progress = get(kwargs, :enable_progress_bar, _prog_meter_enabled()),
        progress_steps = 1,
        advance_to_tstop = !isempty(sim.tstops),
        initializealg = SciMLBase.NoInit(),
        kwargs...,
    )
    if solution.retcode == :Success
        sim.status = SIMULATION_FINALIZED
        sim.results = SimulationResults(
            get_simulation_inputs(sim),
            get_system(sim),
            time_log,
            solution,
        )
    else
        @error("The simulation failed with return code $(solution.retcode)")
        sim.status = SIMULATION_FAILED
    end
end

"""
    execute!(
        sim::Simulation,
        solver;
        kwargs...
    )

Solves the time-domain dynamic simulation model.

# Arguments
- `sim::Simulation` : Initialized simulation object
- `solver` : Solver used for numerical integration. Must be passed correctly depending on the Type of Simulation Model
- `enable_progress_bar::Bool` : Default: `true`. Enables progress bar for the integration routine.
- Additional solver keyword arguments can be included. See [Common Solver Options](https://diffeq.sciml.ai/stable/basics/common_solver_opts/) in the `DifferentialEquations.jl` documentation for more details.
"""
function execute!(sim::Simulation, solver; kwargs...)
    logger = configure_logging(sim, "a"; kwargs...)
    Logging.with_logger(logger) do
        try
            _execute!(sim, solver; kwargs...)
        catch e
            @error "Execution failed" exception = (e, catch_backtrace())
            sim.status = SIMULATION_FAILED
        end
    end
    close(logger)
    return sim.status
end

function read_results(sim::Simulation)
    return sim.results
end

function get_dynamic_wrapper(sim::Simulation, name::AbstractString)
    return get_dynamic_wrapper(get_simulation_inputs(sim), name)
end

"""
    get_setpoints(sim::Simulation)

Function that returns the reference setpoints for all the dynamic devices.

# Arguments

- `sim::Simulation` : Simulation object that contains the initial condition and setpoints.
"""
function get_setpoints(sim::Simulation)
    return get_setpoints(get_simulation_inputs(sim))
end


text:
using Test
using FrankWolfe
using Random
using LinearAlgebra

import FrankWolfe: compute_gradient, compute_value, compute_value_gradient

Random.seed!(123)

@testset "Simple and stochastic function gradient interface" begin
    for n in (1, 5, 20)
        A = rand(Float16, n, n)
        A .+= A'
        b = rand(Float16, n)
        c = rand(Float16)
        @testset "Simple function" begin
            simple_quad(x) = A * x ⋅ x / 2 + b ⋅ x + c
            ∇simple_quad(x) = A * x + b
            f_simple = FrankWolfe.SimpleFunctionObjective(simple_quad, ∇simple_quad)
            for _ in 1:5
                x = randn(n)
                @test FrankWolfe.compute_gradient(f_simple, x) == ∇simple_quad(x)
                @test FrankWolfe.compute_value(f_simple, x) == simple_quad(x)
                @test FrankWolfe.compute_value_gradient(f_simple, x) ==
                      (simple_quad(x), ∇simple_quad(x))
                @test FrankWolfe.compute_gradient(f_simple, convert(Vector{Float32}, x)) isa
                      Vector{Float32}
                @test FrankWolfe.compute_value(f_simple, convert(Vector{Float32}, x)) isa Float32
            end
        end
        @testset "Stochastic function linear regression" begin
            function simple_reg_loss(θ, data_point)
                (xi, yi) = data_point
                (a, b) = (θ[1:end-1], θ[end])
                pred = a ⋅ xi + b
                return (pred - yi)^2 / 2
            end
            function ∇simple_reg_loss(θ, data_point)
                (xi, yi) = data_point
                (a, b) = (θ[1:end-1], θ[end])
                grad_a = xi * (a ⋅ xi + b - yi)
                grad = push!(grad_a, a ⋅ xi + b - yi)
                return grad
            end
            xs = [10 * randn(5) for i in 1:10000]
            params = rand(6) .- 1 # start params in (-1,0)
            bias = 4π
            params_perfect = [1:5; bias]
            @testset "Perfectly representable data" begin
                # Y perfectly representable by X
                data_perfect = [(x, x ⋅ (1:5) + bias) for x in xs]
                f_stoch =
                    FrankWolfe.StochasticObjective(simple_reg_loss, ∇simple_reg_loss, data_perfect)
                @test compute_value(f_stoch, params) > compute_value(f_stoch, params_perfect)
                @test compute_value(f_stoch, params_perfect) ≈ 0 atol = 1e-10
                @test compute_gradient(f_stoch, params_perfect) ≈ zeros(6) atol = 6e-10
                @test !isapprox(compute_gradient(f_stoch, params), zeros(6))
                (f_estimate, g_estimate) = compute_value_gradient(
                    f_stoch,
                    params_perfect,
                    batch_size=length(data_perfect),
                    rng=Random.seed!(33),
                )
                @test f_estimate ≈ 0 atol = 1e-10
                @test g_estimate ≈ zeros(6) atol = 6e-10
                (f_estimate, g_estimate) = compute_value_gradient(
                    f_stoch,
                    params,
                    batch_size=length(data_perfect),
                    rng=Random.seed!(33),
                )
                @test f_estimate ≈ compute_value(
                    f_stoch,
                    params,
                    batch_size=length(data_perfect),
                    rng=Random.seed!(33),
                )
                @test g_estimate ≈ compute_gradient(
                    f_stoch,
                    params,
                    batch_size=length(data_perfect),
                    rng=Random.seed!(33),
                )
            end
            @testset "Noisy data" begin
                data_noisy = [(x, x ⋅ (1:5) + bias + 0.5 * randn()) for x in xs]
                f_stoch_noisy =
                    FrankWolfe.StochasticObjective(simple_reg_loss, ∇simple_reg_loss, data_noisy)
                @test compute_value(f_stoch_noisy, params) >
                      compute_value(f_stoch_noisy, params_perfect)
                # perfect parameters shouldn't have too high of a residual gradient
                @test norm(compute_gradient(f_stoch_noisy, params_perfect)) <=
                      length(data_noisy) * 0.05
                @test norm(compute_gradient(f_stoch_noisy, params_perfect)) <
                      norm(compute_gradient(f_stoch_noisy, params))
                @test !isapprox(compute_gradient(f_stoch_noisy, params), zeros(6))
                (f_estimate, g_estimate) = compute_value_gradient(
                    f_stoch_noisy,
                    params,
                    batch_size=length(data_noisy),
                    rng=Random.seed!(33),
                )
                @test f_estimate ≈ compute_value(
                    f_stoch_noisy,
                    params,
                    batch_size=length(data_noisy),
                    rng=Random.seed!(33),
                )
                @test g_estimate ≈ compute_gradient(
                    f_stoch_noisy,
                    params,
                    batch_size=length(data_noisy),
                    rng=Random.seed!(33),
                )
            end
        end
    end
end


text:
using Dates

effective_date = Date(2020, 01, 15)

days = [369, 462, 553, 735, 1100, 1465]
df = [0.97678403, 0.97153083, 0.96655364, 0.95695995, 0.93829085, 0.91954718]
test_date = effective_date + Day(900)

const dc_act_360 = QuantLib.Time.Actual360()

# test log-linear interpolation
test_curve_loglinear = QuantLib.InterpolatedDiscountCurve(effective_date .+ Day.(days), df, dc_act_360, QuantLib.Math.LogLinear())
disc_factor = QuantLib.discount(test_curve_loglinear, test_date)
x1 = (test_date - test_curve_loglinear.dates[4]).value
dx = (test_curve_loglinear.dates[5]-test_curve_loglinear.dates[4]).value
dy = log(test_curve_loglinear.data[5]) - log(test_curve_loglinear.data[4])
df_recalc = exp(log(test_curve_loglinear.data[4]) + x1*dy/dx)
@test disc_factor == df_recalc

test_curve_linear = QuantLib.InterpolatedDiscountCurve(effective_date .+ Day.(days), df, dc_act_360, QuantLib.Math.LinearInterpolation())
disc_factor = QuantLib.discount(test_curve_linear, test_date)
x1 = (test_date - test_curve_linear.dates[4]).value
dx = (test_curve_linear.dates[5]-test_curve_linear.dates[4]).value
dy = (test_curve_linear.data[5]) - (test_curve_linear.data[4])
df_recalc = ((test_curve_linear.data[4]) + x1*dy/dx)
@test disc_factor == df_recalc

# test joint calendars
const calendar_nyc = QuantLib.Time.USNYSECalendar()
const calendar_target = QuantLib.Time.TargetCalendar()
const calendar_lon = QuantLib.Time.UKLSECalendar()

calendar_nyc_target = QuantLib.Time.JointCalendar(calendar_nyc, calendar_target) # standard constructor
@test calendar_nyc_target isa QuantLib.Time.BusinessCalendar
calendar_nyc_target_lon = QuantLib.Time.JointCalendar(calendar_nyc, calendar_target, calendar_lon) # new constructor
@test calendar_nyc_target_lon isa QuantLib.Time.BusinessCalendar
# at least 2 calendars must be specified for JointCalendar
@test_throws ErrorException QuantLib.Time.JointCalendar(calendar_nyc)


text:
# Autogenerated wrapper script for MibS_jll for i686-linux-musl-cxx03
export libmibs, mibs

using CoinUtils_jll
using Osi_jll
using Clp_jll
using Cgl_jll
using Cbc_jll
using SYMPHONY_jll
using ALPS_jll
using BiCePS_jll
using CHiPPS_BLIS_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("MibS")
JLLWrappers.@declare_library_product(libmibs, "libMibs.so.0")
JLLWrappers.@declare_executable_product(mibs)
function __init__()
    JLLWrappers.@generate_init_header(CoinUtils_jll, Osi_jll, Clp_jll, Cgl_jll, Cbc_jll, SYMPHONY_jll, ALPS_jll, BiCePS_jll, CHiPPS_BLIS_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libmibs,
        "lib/libMibs.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        mibs,
        "bin/mibs",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function [X]=Annulus(DG,P1,P2)
ksi=DG.xw;
n=DG.OrdPoly+1;
X=zeros(n,2);
Rad1=sqrt(P1(1)^2+P1(2)^2);
Rad2=sqrt(P2(1)^2+P2(2)^2);
Rad=0.5*(Rad1*(1-ksi)+...
  Rad2*(1+ksi));
X(:,1)=.5*((1-ksi)*P1(1)+(1+ksi)*P2(1));
X(:,2)=.5*((1-ksi)*P1(2)+(1+ksi)*P2(2));
X=X./sqrt((X(:,1).^2+X(:,2).^2));
X=X.*Rad;
end


text:
using DeferredAcceptance
using Test, Random, LinearAlgebra

@testset "Throws" begin
    @testset "Dim mismatches" begin
        @test_throws AssertionError deferredacceptance([1 2; 2 1], [1 2; 2 1; 3 3], [1, 1, 1])
        @test_throws AssertionError TTC_match([1 2; 2 1; 3 3], [5, 6])
        @test_throws AssertionError choiceaugmented([1 2; 2 1; 3 3], [1, 2])
    end
end

@testset "Discrete input" begin
    @testset "Normal input" begin
        students = [1 2;
                    2 1;
                    3 4;
                    4 3]
        schools = [1 2 2 1;
                   2 1 1 2]
        capacities = [1, 1, 2, 1]

        @test !isempty(deferredacceptance(students, schools, capacities))
        @test !isempty(TTC_match(students, capacities))
    end

    @testset "UInt input" begin
        students = UInt[1 2;
                        2 1;
                        3 4;
                        4 3]
        schools = [1 2 2 1;
                   2 1 1 2]
        capacities = [1, 1, 2, 1]

        @test !isempty(deferredacceptance(students, schools, capacities))
        @test !isempty(deferredacceptance(students, schools, capacities; rev=true))
        @test !isempty(TTC_match(students, capacities))
    end

    @testset "Adjoints" begin
        students = [1 2 3;
                    2 1 3]
        schools = [1 2; 1 2; 2 1]

        capacities = [1, 1, 2]

        @test !isempty(deferredacceptance(students', schools', capacities))
        @test !isempty(TTC_match(students', capacities))
    end

    @testset "1D input" begin
        students = [1, 3, 2, 4]
        schools = [1, 1, 1, 1]
        capacities = [3, 2, 4, 2]

        @test_broken !isempty(deferredacceptance(students, schools, capacities))
        @test_broken !isempty(TTC_match(students, capacities))
    end
end

@testset "Discrete matches" begin
    samp = 5

    @testset "Tiny TTC" begin
        # Basic
        assn = [3, 2, 1]
        students_inv = [1 1 3;
                        2 3 1;
                        3 2 2]
        @test TTC(students_inv, assn) == [1, 2, 3]

        # Lots of cycles
        assn = [1, 3, 2, 4]
        students_inv = [1 2 3 4;
                        2 1 4 3;
                        3 4 1 2;
                        4 3 2 1]
        @test TTC(students_inv, assn) == [1, 2, 3, 4]

        # Underdemanded
        assn = [1, 2]
        students_inv = [4 3;
                        2 1;
                        1 4;
                        3 2]
        @test TTC(students_inv, assn) == [2, 1]

        # Overdemanded
        assn = [1, 2, 2, 3]
        students_inv = [3 1 1 2;
                        1 2 3 1;
                        2 3 2 3]
        @test TTC(students_inv, assn) in ([3, 1, 2, 2], [3, 2, 1, 2])
    end

    @testset "TTC instability" begin
        # About 1% of the time TTC chances upon a stable match. Don't worry.
        res = falses(samp)

        for i in 1:samp
            n = rand(100:200)
            m = rand(20:40)
            students = hcat((randperm(m) for i in 1:n)...)
            schools = hcat((randperm(n) for i in 1:m)...)
            capacities = rand(5:10, m)

            assn, rdist = TTC_match(students, capacities)
            res[i] = isstable(students, schools, capacities, assn) == false
        end

        @test sum(res) / samp ≥ .9
    end

    students = [3 3 4 3 4 3 3 3 3 4;
                4 4 3 4 3 4 4 4 4 3;
                2 1 2 2 2 1 2 2 2 2;
                1 2 1 1 1 2 1 1 1 1]
    schools = [1 5 7 5;
               1 1 1 1;
               5 1 1 1;
               5 1 1 1;
               10 9 7 8;
               1 5 4 5;
               5 1 4 1;
               8 5 7 8;
               8 9 7 8;
               1 5 4 5]

    schools_STB = singletiebreaking(schools)
    schools_MTB = multipletiebreaking(schools)
    schools_WTB = welfaretiebreaking(students, schools, rand(4)')
    capacities = [3, 2, 2, 3]

    @testset "Tiny tiebreaking doesn't fail" begin
        for sch in [schools_STB, schools_MTB, schools_WTB],
            col in eachcol(sch)
            @test all(i in col for i in 1:size(schools)[1])
        end
    end

    @testset "Tiny tiebreaking linear orders agree" begin
        for sch in [schools_STB, schools_MTB, schools_WTB],
            (tb, or) in zip(eachcol(sch), eachcol(schools))
            @test all(diff(or[sortperm(tb)]) .≥ 0)
        end
    end

    @testset "Tiny DA stability" begin
        for sch in [schools_STB, schools_MTB, schools_WTB]
            assn, ranks = deferredacceptance(students, sch, [3, 2, 2, 3])
            @test isstable(students, sch, capacities, assn)
            # Checks linear orders:
            # @test all([all(diff(schools[:, i][sortperm(col)]) .≥ 0)
            #            for (i, col) in enumerate(eachcol(sch))])
        end
    end

    @testset "Big DA stability, fwd" begin
        for _ in 1:samp
            n = rand(100:200)
            m = rand(20:40)
            students = hcat((randperm(m) for i in 1:n)...)
            schools = hcat((randperm(n) for i in 1:m)...)
            capacities = rand(5:10, m)

            assn, ranks = deferredacceptance(students, schools, capacities)
            @test isstable(students, schools, capacities, assn)
        end
    end

    @testset "Big DA stability, rev" begin
        for _ in 1:samp
            n = rand(100:200)
            m = rand(20:40)
            students = hcat((randperm(m) for i in 1:n)...)
            schools = hcat((randperm(n) for i in 1:m)...)
            capacities = rand(5:10, m)

            assn, ranks = deferredacceptance(students, schools, capacities, rev=true)
            @test isstable(students, schools, capacities, assn)
        end
    end

    @testset "CADA stability" begin
        for _ in 1:samp
            m = rand(10:20)
            cap = rand(5:10)
            n = m * cap

            students = repeat(randperm(m), 1, n)
            targets = rand(1:m, n)

            schools = ones(Int, n, m)
            capacities = ones(Int, m) .* cap

            schools_CADA = choiceaugmented(schools, targets)
            assn, rdist = deferredacceptance(students, schools_CADA, capacities)
            @test isstable(students, schools, capacities, assn)
        end
    end

    @testset "DA rank dist" begin
        for _ in 1:samp
            n = rand(100:200)
            m = rand(20:40)
            students = hcat((randperm(m) for i in 1:n)...)
            schools = hcat((randperm(n) for i in 1:m)...)
            capacities = rand(5:10, m)

            out = DA_rank_dist(students, schools, capacities)
            @test out[end] == min(n, sum(capacities))
        end
    end
end

@testset "Nonatomic input" begin
    @testset "Normal input" begin
        students = [1 2;
                    2 1;
                    3 4;
                    4 3]
        students_dist = [0.4, 0.6]
        schools = [1 2 2 1;
                   2 1 1 2]
        capacities = [.2, .2, .2, .3]

        @test !isempty(nonatomicdeferredacceptance(students, students_dist, schools, capacities))
        @test !isempty(nonatomicdeferredacceptance(students, students_dist, nothing, capacities))
    end

    @testset "Adjoints" begin
        students = [1 2 3;
                    2 1 3]
        students_dist = [.1, .9]
        schools = [1 2; 1 2; 2 1]

        capacities = [.1, .6, .1]

        @test !isempty(nonatomicdeferredacceptance(students', students_dist, schools', capacities))
        @test !isempty(nonatomicdeferredacceptance(students', students_dist, nothing, capacities))
    end

    @testset "UInt input" begin
        students = UInt[1 2;
                        2 1;
                        3 4;
                        4 3]
        students_dist = [0.4, 0.6]
        schools = UInt[1 2 2 1;
                       2 1 1 2]
        capacities = [.2, .2, .2, .3]

        @test !isempty(nonatomicdeferredacceptance(students, students_dist, schools, capacities))
        @test !isempty(nonatomicdeferredacceptance(students, students_dist, nothing, capacities))
    end

    @testset "1D input" begin
        #=  Would actually be nice for this to work, as knowing the allocation given
            that everyone has the same preference order is useful.      =#

        students = [1, 3, 2, 4]
        students_dist = [1.]
        schools = [1, 1, 1, 1]
        capacities = [.1, .6, .1, .2]

        @test_broken !isempty(nonatomicdeferredacceptance(students, students_dist, schools, capacities))
        @test_broken !isempty(nonatomicdeferredacceptance(students, students_dist, nothing, capacities))
    end
end

@testset "Nonatomic DA" begin
    samp = 5

    @testset "Azevedo and Leshno (2016)'s example'" begin
        students = [1 2; 2 1]
        students_dist = [0.5, 0.5]
        capacities = [0.25, 0.5]

        assn, rdist, cutoffs = nonatomicdeferredacceptance(students, students_dist, nothing, capacities;
                                            return_cutoffs=true)

        @test cutoffs ≈ [√17 + 1, √17 - 1] ./ 8
        @test ismarketclearing(students, students_dist, capacities, cutoffs)
        @test ismarketclearing(students, students_dist, capacities, [0.1, 0.1])==false
        @test ismarketclearing(students, students_dist, capacities, [0.9, 0.9])==false
    end

    @testset "Assignment and demand operators" begin
        for _ in 1:samp
            n = rand(5:10)    # Number of student profiles in continuum
            m = rand(10:20)    # Number of schools
            α = 0.5 + rand()

            students = hcat((randperm(m) for i = 1:n)...)   # Student profiles
            students_dist = rand(n)                         # Percentage of total student population
            students_dist /= sum(students_dist)             # associated with each profile

            students_inv = mapslices(invperm, students, dims=1)

            capacities = rand(m)
            capacities /= (α * sum(capacities))

            cutoffs = rand(m)       # Should be incorrect

            D_l = demandfrompreflists(students, students_dist, cutoffs)
            assn_d, D_d = assignmentfrompreflists(students_inv, students_dist, cutoffs;
                                              return_demands=true)

            @test D_l ≈ D_d
            @test sum(assn_d, dims=1) ≈ students_dist'
            @test ismarketclearing(students, students_dist, capacities, cutoffs)==false
        end
    end

    @testset "Cutoff algorithm market-clearing" begin
        for _ in 1:samp
            n = rand(5:10)    # Number of student profiles in continuum
            m = rand(5:10)    # Number of schools
            α = 0.5 + rand()       # Proportion by which mkt is overdemanded

            students = hcat((randperm(m) for i = 1:n)...)   # Student profiles
            students_dist = rand(n)                         # Percentage of total student population
            students_dist /= sum(students_dist)             # associated with each profile

            capacities = rand(m)                            # Percentage of total student population
            capacities /= (α * sum(capacities))

            students_inv = mapslices(invperm, students, dims=1)
            cutoffs = nonatomicdeferredacceptance_iid(students, students_dist, capacities)

            D_l = demandfrompreflists(students, students_dist, cutoffs)
            assn_d, D_d = assignmentfrompreflists(students_inv, students_dist, cutoffs;
                                              return_demands=true)

            @test D_l ≈ D_d
            @test sum(assn_d, dims=1) ≈ students_dist'
            @test vec(sum(assn_d[1:end - 1, :], dims=2)) ≤ capacities .+ 1e-8
            @test ismarketclearing(students, students_dist, capacities, cutoffs)
        end
    end

    @testset "No school prefs" begin
        for _ in 1:samp
            n = rand(5:10)    # Number of student profiles in continuum
            m = rand(5:10)    # Number of schools
            α = 0.5 + rand()   # Proportion by which mkt is overdemanded

            students = hcat((randperm(m) for i = 1:n)...)   # Student profiles
            students_dist = rand(n)                         # Percentage of total student population
            students_dist /= sum(students_dist)             # associated with each profile

            capacities = rand(m)                            # Percentage of total student population
            capacities /= (α * sum(capacities))

            assn, rdist, cutoffs = nonatomicdeferredacceptance(students, students_dist, nothing, capacities;
                                                return_cutoffs=true)

            @test sum(assn, dims=1) ≈ students_dist'
            @test sum(assn, dims=2)[1:end - 1] ≤ capacities .+ 1e-4
            @test ismarketclearing(students, students_dist, capacities, cutoffs)
        end
    end

    @testset "With school prefs" begin
        for _ in 1:samp
            n = rand(5:10)    # Number of student profiles in continuum
            m = rand(5:10)    # Number of schools
            α = 0.5 + rand()   # Proportion by which mkt is overdemanded

            students = hcat((randperm(m) for i = 1:n)...)   # Student profiles
            students_dist = rand(n)                         # Percentage of total student population
            students_dist /= sum(students_dist)             # associated with each profile

            capacities = rand(m)                            # Percentage of total student population
            capacities /= (α * sum(capacities))
            schools = hcat((randperm(n) for i = 1:m)...)

            assn = nonatomicdeferredacceptance(students, students_dist, schools, capacities)[1]

            @test sum(assn, dims=1) ≈ students_dist'
            @test sum(assn, dims=2)[1:end - 1] ≤ capacities .+ 1e-8
        end
    end

    @testset "Reverse nonatomic lite" begin
        for _ in 1:samp
            n = rand(5:10)    # Number of student profiles in continuum
            m = rand(5:10)    # Number of schools
            α = 0.5 + rand()   # Proportion by which mkt is overdemanded

            students = hcat((randperm(m) for i = 1:n)...)   # Student profiles
            students_dist = rand(n)                         # Percentage of total student population
            students_dist /= sum(students_dist)             # associated with each profile

            capacities = rand(m)                            # Percentage of total student population
            capacities /= (α * sum(capacities))

            # Could just use DA_nonatomic_lite here but this tests the wrapper too
            _, _, cutoffs_fwd = nonatomicdeferredacceptance(students, students_dist, nothing, capacities;
                                             return_cutoffs=true)
            _, _, cutoffs_rev = nonatomicdeferredacceptance(students, students_dist, nothing, capacities;
                                             return_cutoffs=true, rev=true)

            # Cutoffs should be equal in this idealized case
            @test cutoffs_fwd ≈ cutoffs_rev #atol=1e-6
        end
    end
end

@testset "Funny demand functions" begin
    samp = 5

    @testset "MLN iid" begin

        @testset "WGS" begin
            for _ in 1:samp
                m = rand(5:10)

                qualities = rand(m)

                cutoffs = rand(m)
                delta = copy(cutoffs)

                delta[2:end] .= 1

                demand(cut) = demandfromMNL_iid(qualities, cut)

                # When all schools but 1 increase their cutoffs
                out_orig = demand(cutoffs)
                out_pert = demand(cutoffs + rand() * (delta .- cutoffs))

                # 1's demand should increase
                @test out_orig[1] ≤ out_pert[1]
            end
        end

        @testset "Azevedo and Leshno (2016)'s example" begin
            qualities = [1., 1]         # Schools equally preferable
            capacities = [0.25, 0.5]

            demand(cut) = demandfromMNL_iid(qualities, cut)

            actual = [√17 + 1, √17 - 1] ./ 8

            @test nonatomicdeferredacceptance_iid(demand, capacities) ≈ actual
            @test nonatomicdeferredacceptance_iid(demand, capacities, rev=true) ≈ actual
            @test nonatomictatonnement(demand, capacities) ≈ actual
            @test nonatomicsecant(demand, capacities) ≈ actual
        end

        @testset "DA-lite fwd, rev; secant; tatonnement agree" begin
            for _ in 1:samp
                m = rand(5:10)
                qualities = randexp(m)
                capacities = randexp(m)
                capacities ./= (0.5 + rand()) .* sum(capacities)

                demand(cut) = demandfromMNL_iid(qualities, cut)

                cutoffs = nonatomicdeferredacceptance_iid(demand, capacities)

                @test cutoffs ≈
                      nonatomicdeferredacceptance_iid(demand, capacities; rev=true)
                @test cutoffs ≈
                      nonatomictatonnement(demand, capacities)
                @test cutoffs ≈
                      nonatomicsecant(demand, capacities, verbose=false)

                @test ismarketclearing(qualities, capacities, cutoffs)
                @test ismarketclearing(demand, capacities, cutoffs)
            end
        end

    end

    @testset "MLN one test" begin

        @testset "WGS" begin
            m = rand(5:10)

            qualities = rand(m)

            cutoffs = rand(m)
            delta = copy(cutoffs)

            delta[2:end] .= 1

            demand(cut) = demandfromMNL_singlescore(qualities, cut)

            # When all schools but 1 increase their cutoffs
            out_orig = demand(cutoffs)
            out_pert = demand(cutoffs + rand() * (delta .- cutoffs))

            # 1's demand should increase
            @test out_orig[1] ≤ out_pert[1]
        end

        @testset "Tatonnement clears mkt" begin
            for _ in 1:samp
                m = rand(10:25)
                qualities = rand(m)
                capacities = randexp(m)
                capacities ./= (0.5 + rand()) .* sum(capacities)

                demand(cut) = demandfromMNL_singlescore(qualities, cut)

                cut = nonatomictatonnement(demand, capacities, maxit=800)
                @test ismarketclearing(demand, capacities, cut)
            end
        end
    end

    @testset "p MNL profiles, t tests" begin

        @testset "WGS" begin
            for _ in 1:samp
                m = rand(5:10)
                p = rand(3:5)
                t = rand(3:5)

                qualities = randexp(m, p)
                profile_dist = rand(p)
                profile_dist ./= sum(profile_dist)

                blends = rand(m, t)
                blends ./= sum(blends, dims=2)

                cutoffs = rand(m)

                delta = copy(cutoffs)
                delta[2:end] .= 1

                demand(cut) = demandfrommixedMNL_tscores(qualities, profile_dist, blends, cut)

                # When all schools but 1 increase their cutoffs
                out_orig = demand(cutoffs)
                out_pert = demand(cutoffs + rand() * (delta .- cutoffs))

                # 1's demand should increase
                @test out_orig[1] ≤ out_pert[1] + 1e-3
            end
        end

        @testset "Equivalence with MNL iid" begin
            # We use Monte Carlo integration here, so have to set pretty loose tolerances.

            for _ in 1:samp
                m = rand(5:10)
                qualities = rand(m)
                cutoffs = rand(m)

                out1 = demandfromMNL_iid(qualities, cutoffs)

                # This one hangs when using exact mode, due to too many
                # redundant constraints for Polyhedra.jl to figure out.
                out2 = demandfrommixedMNL_tscores(reshape(qualities, :, 1),
                                           [1.],
                                           Matrix{Float64}(I, m, m)[randperm(m), :],
                                           montecarlo=true,
                                           cutoffs)

                @test isapprox(out1, out2, atol=1e-2)
            end
        end

        @testset "Equivalence with MNL one test" begin
            for _ in 1:samp
                m = rand(5:10)
                qualities = rand(m)
                cutoffs = rand(m)

                out1 = demandfromMNL_singlescore(qualities, cutoffs)

                out2 = demandfrommixedMNL_tscores(reshape(qualities, :, 1),
                                           [1.],
                                           ones(m, 1),
                                           cutoffs)

                @test isapprox(out1, out2, atol=1e-2)
            end
        end

        samp = 3

        @testset "Tatonnement clears mkt" begin
            for _ in 1:samp
                m = rand(5:7)
                p = rand(5:10)
                t = rand(2:4)

                qualities = randexp(m, p)
                profile_dist = rand(p)
                profile_dist ./= sum(profile_dist)

                blends = rand(m, t)
                blends ./= sum(blends, dims=2)

                capacities = randexp(m)
                capacities ./= (0.5 + rand()) .* sum(capacities)

                nit = 1

                # Kinda hacky but works. Use MC for local search then start
                # return precise values later on.
                function demand(cut)
                    if nit < 100
                        nit += 1
                        return demandfrommixedMNL_tscores(qualities, profile_dist, blends, cut,
                                                   montecarlo=true, n_points=1000)
                    else
                        return demandfrommixedMNL_tscores(qualities, profile_dist, blends, cut,
                                                   montecarlo=false)
                    end
                end

                cut = nonatomictatonnement(demand, capacities, maxit=200, β=.0, tol=1e-6)

                @test ismarketclearing(demand, capacities, cut, tol=1e-4)
            end
        end

        @testset "Secant clears mkt" begin
            for _ in 1:samp
                m = rand(5:7)
                p = rand(5:10)
                t = rand(2:4)

                qualities = randexp(m, p)
                profile_dist = rand(p)
                profile_dist ./= sum(profile_dist)

                blends = rand(m, t)
                blends ./= sum(blends, dims=2)

                capacities = randexp(m)
                capacities ./= (0.5 + rand()) .* sum(capacities)

                nit = 1

                function demand(cut)
                    if nit < 50
                        nit += 1
                        return demandfrommixedMNL_tscores(qualities, profile_dist, blends, cut,
                                                   montecarlo=true, n_points=1000)
                    else
                        return demandfrommixedMNL_tscores(qualities, profile_dist, blends, cut,
                                                   montecarlo=false)
                    end
                end

                cut = nonatomicsecant(demand, capacities, maxit=200, verbose=false, tol=1e-6)

                @test ismarketclearing(demand, capacities, cut, tol=1e-2)
            end
        end
    end

end

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
function show_indented(io::IO,x,indent)
    print(io," "^indent,"- ")
    show(io,x)
end

show(io::IO,x::PackageManager) = print(io,"$(typeof(x)) package $(pkg_name(x))")
show(io::IO,x::SimpleBuild) = print(io,"Simple Build Process")
show(io::IO,x::Sources) = print(io,"Sources")
show(io::IO,x::Binaries) = print(io,"Binaries")
show(io::IO,x::Autotools) = print(io,"Autotools Build")