text:
@doc raw"""
    CompressibleEulerEquations2D

The compressible Euler equations for an ideal gas in two space dimensions.
"""
struct CompressibleEulerEquations2D{RealT<:Real} <: AbstractCompressibleEulerEquations{2, 4}
  gamma::RealT
end

varnames(::typeof(cons2cons), ::CompressibleEulerEquations2D) = ("rho", "rho_v1", "rho_v2", "rho_e")
varnames(::typeof(cons2prim), ::CompressibleEulerEquations2D) = ("rho", "v1", "v2", "p")

# Set initial conditions at physical location `x` for time `t`
"""
    initial_condition_constant(x, t, equations::CompressibleEulerEquations2D)

A constant initial condition to test free-stream preservation.
"""
function initial_condition_constant(x, t, equations::CompressibleEulerEquations2D)
  rho = 1.0
  rho_v1 = 0.1
  rho_v2 = -0.2
  rho_e = 10.0
  return SVector(rho, rho_v1, rho_v2, rho_e)
end

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

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

  rho = ini
  rho_v1 = ini
  rho_v2 = ini
  rho_e = ini^2

  return SVector(rho, rho_v1, rho_v2, rho_e)
end

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

Source terms used for convergence tests in combination with
[`initial_condition_convergence_test`](@ref)
(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
"""
@inline function source_terms_convergence_test(u, x, t, equations::CompressibleEulerEquations2D)
  # Same settings as in `initial_condition`
  c = 2
  A = 0.1
  L = 2
  f = 1/L
  ω = 2 * pi * f
  γ = equations.gamma

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

  du1 = tmp1
  du2 = tmp5
  du3 = tmp5
  du4 = 2*((tmp6 - 1)*tmp3 + tmp6*γ)*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)
end

"""
    initial_condition_density_pulse(x, t, equations::CompressibleEulerEquations2D)

A Gaussian pulse in the density with constant velocity and pressure; reduces the
compressible Euler equations to the linear advection equations.
"""
function initial_condition_density_pulse(x, t, equations::CompressibleEulerEquations2D)
  rho = 1 + exp(-(x[1]^2 + x[2]^2))/2
  v1 = 1
  v2 = 1
  rho_v1 = rho * v1
  rho_v2 = rho * v2
  p = 1
  rho_e = p/(equations.gamma - 1) + 1/2 * rho * (v1^2 + v2^2)
  return SVector(rho, rho_v1, rho_v2, rho_e)
end

"""
    initial_condition_density_wave(x, t, equations::CompressibleEulerEquations2D)

A sine wave in the density with constant velocity and pressure; reduces the
compressible Euler equations to the linear advection equations.
This setup is the test case for stability of EC fluxes from paper
- Gregor J. Gassner, Magnus Svärd, Florian J. Hindenlang (2020)
  Stability issues of entropy-stable and/or split-form high-order schemes
  [arXiv: 2007.09026](https://arxiv.org/abs/2007.09026)
with the following parameters
- domain [-1, 1]
- mesh = 4x4
- polydeg = 5
"""
function initial_condition_density_wave(x, t, equations::CompressibleEulerEquations2D)
  v1 = 0.1
  v2 = 0.2
  rho = 1 + 0.98 * sinpi(2 * (x[1] + x[2] - t * (v1 + v2)))
  rho_v1 = rho * v1
  rho_v2 = rho * v2
  p = 20
  rho_e = p / (equations.gamma - 1) + 1/2 * rho * (v1^2 + v2^2)
  return SVector(rho, rho_v1, rho_v2, rho_e)
end

"""
    initial_condition_pressure_pulse(x, t, equations::CompressibleEulerEquations2D)

A Gaussian pulse in the pressure with constant velocity and density.
"""
function initial_condition_pressure_pulse(x, t, equations::CompressibleEulerEquations2D)
  rho = 1
  v1 = 1
  v2 = 1
  rho_v1 = rho * v1
  rho_v2 = rho * v2
  p = 1 + exp(-(x[1]^2 + x[2]^2))/2
  rho_e = p/(equations.gamma - 1) + 1/2 * rho * (v1^2 + v2^2)
  return SVector(rho, rho_v1, rho_v2, rho_e)
end

"""
    initial_condition_density_pressure_pulse(x, t, equations::CompressibleEulerEquations2D)

A Gaussian pulse in density and pressure with constant velocity.
"""
function initial_condition_density_pressure_pulse(x, t, equations::CompressibleEulerEquations2D)
  rho = 1 + exp(-(x[1]^2 + x[2]^2))/2
  v1 = 1
  v2 = 1
  rho_v1 = rho * v1
  rho_v2 = rho * v2
  p = 1 + exp(-(x[1]^2 + x[2]^2))/2
  rho_e = p/(equations.gamma - 1) + 1/2 * rho * (v1^2 + v2^2)
  return SVector(rho, rho_v1, rho_v2, rho_e)
end

"""
    initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)

The classical isentropic vortex test case of
- Chi-Wang Shu (1997)
  Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory
  Schemes for Hyperbolic Conservation Laws
  [NASA/CR-97-206253](https://ntrs.nasa.gov/citations/19980007543)
"""
function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)
  # needs appropriate mesh size, e.g. [-10,-10]x[10,10]
  # make sure that the inicenter does not exit the domain, e.g. T=10.0
  # initial center of the vortex
  inicenter = SVector(0.0, 0.0)
  # size and strength of the vortex
  iniamplitude = 0.2
  # base flow
  rho = 1.0
  v1 = 1.0
  v2 = 1.0
  vel = SVector(v1, v2)
  p = 10.0
  vec = SVector(v1, v2)
  rt = p / rho                  # ideal gas equation
  cent = inicenter + vel*t      # advection of center
  cent = x - cent               # distance to centerpoint
  #cent=cross(iniaxis,cent)     # distance to axis, tangent vector, length r
  # cross product with iniaxis = [0,0,1]
  cent = SVector(-cent[2], cent[1])
  r2 = cent[1]^2 + cent[2]^2
  du = iniamplitude/(2*π)*exp(0.5*(1-r2)) # vel. perturbation
  dtemp = -(equations.gamma-1)/(2*equations.gamma*rt)*du^2            # isentrop
  rho = rho * (1+dtemp)^(1\(equations.gamma-1))
  vel = vel + du*cent
  v1, v2 = vel
  p = p * (1+dtemp)^(equations.gamma/(equations.gamma-1))
  prim = SVector(rho, v1, v2, p)
  return prim2cons(prim, equations)
end

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

A 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::CompressibleEulerEquations2D)
  # 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
  rho = r > 0.5 ? 1.0 : 1.1691
  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(rho, v1, v2, p), equations)
end

"""
    initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)

A medium 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_blast_wave(x, t, equations::CompressibleEulerEquations2D)
  # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
  # 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
  rho = r > 0.5 ? 1.0 : 1.1691
  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.0E-3 : 1.245

  return prim2cons(SVector(rho, v1, v2, p), equations)
end

"""
    initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)

The Sedov blast wave setup based on Flash
- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000
"""
function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
  # 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)

  # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000
  r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
  # r0 = 0.5 # = more reasonable setup
  E = 1.0
  p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2)
  p0_outer = 1.0e-5 # = true Sedov setup
  # p0_outer = 1.0e-3 # = more reasonable setup

  # Calculate primitive variables
  rho = 1.0
  v1  = 0.0
  v2  = 0.0
  p   = r > r0 ? p0_outer : p0_inner

  return prim2cons(SVector(rho, v1, v2, p), equations)
end

"""
    initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)

The Sedov blast wave setup based on Flash
- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000
with smaller strength of the initial discontinuity.
"""
function initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
  # 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)

  # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000
  r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
  # r0 = 0.5 # = more reasonable setup
  E = 1.0
  p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2)
  # p0_outer = 1.0e-5 # = true Sedov setup
  p0_outer = 1.0e-3 # = more reasonable setup

  # Calculate primitive variables
  rho = 1.0
  v1  = 0.0
  v2  = 0.0
  p   = r > r0 ? p0_outer : p0_inner

  return prim2cons(SVector(rho, v1, v2, p), equations)
end

"""
    initial_condition_khi(x, t, equations::CompressibleEulerEquations2D)

The classical Kelvin-Helmholtz instability based on
- https://rsaa.anu.edu.au/research/established-projects/fyris/2-d-kelvin-helmholtz-test
"""
function initial_condition_khi(x, t, equations::CompressibleEulerEquations2D)
  # https://rsaa.anu.edu.au/research/established-projects/fyris/2-d-kelvin-helmholtz-test
  # change discontinuity to tanh
  # typical resolution 128^2, 256^2
  # domain size is [-0.5,0.5]^2
  dens0 = 1.0 # outside density
  dens1 = 2.0 # inside density
  velx0 = -0.5 # outside velocity
  velx1 = 0.5 # inside velocity
  slope = 50 # used for tanh instead of discontinuous initial condition
  # pressure equilibrium
  p     = 2.5
  # density
  rho = dens0 + (dens1-dens0) * 0.5*(1+(tanh(slope*(x[2]+0.25)) - (tanh(slope*(x[2]-0.25)) + 1)))
  if iszero(t) # initial condition
    #  y velocity v2 is only white noise
    v2  = 0.01*(rand(Float64,1)[1]-0.5)
    #  x velocity is also augmented with noise
    v1 = velx0 + (velx1-velx0) * 0.5*(1+(tanh(slope*(x[2]+0.25)) - (tanh(slope*(x[2]-0.25)) + 1)))+0.01*(rand(Float64,1)[1]-0.5)
  else # background values to compute reference values for CI
    v2 = 0.0
    v1 = velx0 + (velx1-velx0) * 0.5*(1+(tanh(slope*(x[2]+0.25)) - (tanh(slope*(x[2]-0.25)) + 1)))
  end
  return prim2cons(SVector(rho, v1, v2, p), equations)
end

"""
    initial_condition_blob(x, t, equations::CompressibleEulerEquations2D)

The blob test case taken from
- Agertz et al. (2006)
  Fundamental differences between SPH and grid methods
  [arXiv: astro-ph/0610051](https://arxiv.org/abs/astro-ph/0610051)
"""
function initial_condition_blob(x, t, equations::CompressibleEulerEquations2D)
  # blob test case, see Agertz et al. https://arxiv.org/pdf/astro-ph/0610051.pdf
  # other reference: https://arxiv.org/pdf/astro-ph/0610051.pdf
  # change discontinuity to tanh
  # typical domain is rectangular, we change it to a square, as Trixi can only do squares
  # resolution 128^2, 256^2
  # domain size is [-20.0,20.0]^2
  # gamma = 5/3 for this test case
  R = 1.0 # radius of the blob
  # background density
  dens0 = 1.0
  Chi = 10.0 # density contrast
  # reference time of characteristic growth of KH instability equal to 1.0
  tau_kh = 1.0
  tau_cr = tau_kh/1.6 # crushing time
  # determine background velocity
  velx0 = 2*R*sqrt(Chi)/tau_cr
  vely0 = 0.0
  Ma0 = 2.7 # background flow Mach number Ma=v/c
  c = velx0/Ma0 # sound speed
  # use perfect gas assumption to compute background pressure via the sound speed c^2 = gamma * pressure/density
  p0 = c*c*dens0/equations.gamma
  # initial center of the blob
  inicenter = [-15,0]
  x_rel = x-inicenter
  r = sqrt(x_rel[1]^2 + x_rel[2]^2)
  # steepness of the tanh transition zone
  slope = 2
  # density blob
  dens = dens0 + (Chi-1) * 0.5*(1+(tanh(slope*(r+R)) - (tanh(slope*(r-R)) + 1)))
  # velocity blob is zero
  velx = velx0 - velx0 * 0.5*(1+(tanh(slope*(r+R)) - (tanh(slope*(r-R)) + 1)))
  return prim2cons(SVector(dens, velx, vely0, p0), equations)
end

"""
    initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::CompressibleEulerEquations2D)

Setup used for convergence tests of the Euler equations with self-gravity used in
- Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020)
  A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics
  [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593)
in combination with [`source_terms_eoc_test_coupled_euler_gravity`](@ref)
or [`source_terms_eoc_test_euler`](@ref).
"""
function initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::CompressibleEulerEquations2D)
  # OBS! this assumes that γ = 2 other manufactured source terms are incorrect
  if equations.gamma != 2.0
    error("adiabatic constant must be 2 for the coupling convergence test")
  end
  c = 2.0
  A = 0.1
  ini = c + A * sin(pi * (x[1] + x[2] - t))
  G = 1.0 # gravitational constant

  rho = ini
  v1 = 1.0
  v2 = 1.0
  p = ini^2 * G / pi # * 2 / ndims, but ndims==2 here

  return prim2cons(SVector(rho, v1, v2, p), equations)
end

"""
    source_terms_eoc_test_coupled_euler_gravity(u, x, t, equations::CompressibleEulerEquations2D)

Setup used for convergence tests of the Euler equations with self-gravity used in
- Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020)
  A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics
  [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593)
in combination with [`initial_condition_eoc_test_coupled_euler_gravity`](@ref).
"""
@inline function source_terms_eoc_test_coupled_euler_gravity(u, x, t, equations::CompressibleEulerEquations2D)
  # Same settings as in `initial_condition_eoc_test_coupled_euler_gravity`
  c = 2.0
  A = 0.1
  G = 1.0 # gravitational constant, must match coupling solver
  C_grav = -2.0 * G / pi # 2 == 4 / ndims

  x1, x2 = x
  si, co = sincos(pi * (x1 + x2 - t))
  rhox = A * pi * co
  rho  = c + A *  si

  du1 = rhox
  du2 = rhox
  du3 = rhox
  du4 = (1.0 - C_grav*rho)*rhox

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

"""
    source_terms_eoc_test_euler(u, x, t, equations::CompressibleEulerEquations2D)

Setup used for convergence tests of the Euler equations with self-gravity used in
- Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020)
  A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics
  [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593)
in combination with [`initial_condition_eoc_test_coupled_euler_gravity`](@ref).
"""
@inline function source_terms_eoc_test_euler(u, x, t, equations::CompressibleEulerEquations2D)
  # Same settings as in `initial_condition_eoc_test_coupled_euler_gravity`
  c = 2.0
  A = 0.1
  G = 1.0
  C_grav = -2 * G / pi # 2 == 4 / ndims

  x1, x2 = x
  si, co = sincos(pi * (x1 + x2 - t))
  rhox = A * pi * co
  rho  = c + A *  si

  du1 = rhox
  du2 = rhox * (1 -     C_grav * rho)
  du3 = rhox * (1 -     C_grav * rho)
  du4 = rhox * (1 - 3 * C_grav * rho)

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

"""
    initial_condition_sedov_self_gravity(x, t, equations::CompressibleEulerEquations2D)

Adaptation of the Sedov blast wave with self-gravity taken from
- Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020)
  A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics
  [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593)
based on
- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000
Should be used together with [`boundary_condition_sedov_self_gravity`](@ref).
"""
function initial_condition_sedov_self_gravity(x, t, equations::CompressibleEulerEquations2D)
  # Set up polar coordinates
  r = sqrt(x[1]^2 + x[2]^2)

  # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000
  r0 = 0.125 # = 4.0 * smallest dx (for domain length=8 and max-ref=8)
  E = 1.0
  p_inner   = (equations.gamma - 1) * E / (pi * r0^2)
  p_ambient = 1e-5 # = true Sedov setup

  # Calculate primitive variables
  # use a logistic function to tranfer density value smoothly
  L  = 1.0    # maximum of function
  x0 = 1.0    # center point of function
  k  = -150.0 # sharpness of transfer
  logistic_function_rho = L/(1.0 + exp(-k*(r - x0)))
  rho_ambient = 1e-5
  rho = max(logistic_function_rho, rho_ambient) # clip background density to not be so tiny

  # velocities are zero
  v1 = 0.0
  v2 = 0.0

  # use a logistic function to tranfer pressure value smoothly
  logistic_function_p = p_inner/(1.0 + exp(-k*(r - r0)))
  p = max(logistic_function_p, p_ambient)

  return prim2cons(SVector(rho, v1, v2, p), equations)
end

"""
    boundary_condition_sedov_self_gravity(u_inner, orientation, direction, x, t,
                                          surface_flux_function,
                                          equations::CompressibleEulerEquations2D)

Adaptation of the Sedov blast wave with self-gravity taken from
- Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020)
  A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics
  [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593)
based on
- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000
Should be used together with [`initial_condition_sedov_self_gravity`](@ref).
"""
function boundary_condition_sedov_self_gravity(u_inner, orientation, direction, x, t,
                                               surface_flux_function,
                                               equations::CompressibleEulerEquations2D)
  # velocities are zero, density/pressure are ambient values according to
  # initial_condition_sedov_self_gravity
  rho = 1e-5
  v1 = 0.0
  v2 = 0.0
  p = 1e-5

  u_boundary = prim2cons(SVector(rho, v1, v2, p), 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

"""
    boundary_state_slip_wall(u_internal, normal_direction::AbstractVector,
                             equations::CompressibleEulerEquations2D)

Determine the external solution value for a slip wall condition. Sets the normal
velocity of the the exterior fictitious element to the negative of the internal value.

!!! warning "Experimental code"
    This wall function can change any time.
"""
@inline function boundary_state_slip_wall(u_internal, normal_direction::AbstractVector,
                                          equations::CompressibleEulerEquations2D)

  # normalize the outward pointing direction
  normal = normal_direction / norm(normal_direction)

  # compute the normal and tangential components of the velocity
  u_normal  = normal[1] * u_internal[2] + normal[2] * u_internal[3]
  u_tangent = (u_internal[2] - u_normal * normal[1], u_internal[3] - u_normal * normal[2])

  return SVector(u_internal[1],
                 u_tangent[1] - u_normal * normal[1],
                 u_tangent[2] - u_normal * normal[2],
                 u_internal[4])
end

# Calculate 1D flux for a single point
@inline function flux(u, orientation::Integer, equations::CompressibleEulerEquations2D)
  rho, rho_v1, rho_v2, rho_e = u
  v1 = rho_v1/rho
  v2 = rho_v2/rho
  p = (equations.gamma - 1) * (rho_e - 1/2 * rho * (v1^2 + v2^2))
  if orientation == 1
    f1 = rho_v1
    f2 = rho_v1 * v1 + p
    f3 = rho_v1 * v2
    f4 = (rho_e + p) * v1
  else
    f1 = rho_v2
    f2 = rho_v2 * v1
    f3 = rho_v2 * v2 + p
    f4 = (rho_e + p) * v2
  end
  return SVector(f1, f2, f3, f4)
end

# Calculate 1D flux for a single point in the normal direction
# Note, this directional vector is not normalized
@inline function flux(u, normal_direction::AbstractVector, equations::CompressibleEulerEquations2D)
  rho, rho_v1, rho_v2, rho_e = u
  v1 = rho_v1/rho
  v2 = rho_v2/rho
  p = (equations.gamma - 1) * (rho_e - 1/2 * rho * (v1^2 + v2^2))

  v_normal = v1 * normal_direction[1] + v2 * normal_direction[2]
  rho_v_normal = rho * v_normal
  f1 = rho_v_normal
  f2 = rho_v_normal * v1 + p * normal_direction[1]
  f3 = rho_v_normal * v2 + p * normal_direction[2]
  f4 = (rho_e + p) * v_normal
  return SVector(f1, f2, f3, f4)
end

"""
    function flux_shima_etal(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D)

This flux is is a modification of the original kinetic energy preserving two-point flux by
- Kuya, Totani and Kawai (2018)
  Kinetic energy and entropy preserving schemes for compressible flows
  by split convective forms
  [DOI: 10.1016/j.jcp.2018.08.058](https://doi.org/10.1016/j.jcp.2018.08.058)
The modification is in the energy flux to guarantee pressure equilibrium and was developed by
- Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020)
  Preventing spurious pressure oscillations in split convective form discretizations for
  compressible flows
"""
@inline function flux_shima_etal(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations2D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_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
  p_ll  = (equations.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2 + v2_ll^2))
  p_rr  = (equations.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2 + v2_rr^2))

  # Average each factor of products in flux
  rho_avg = 1/2 * (rho_ll + rho_rr)
  v1_avg  = 1/2 * ( v1_ll +  v1_rr)
  v2_avg  = 1/2 * ( v2_ll +  v2_rr)
  p_avg   = 1/2 * (  p_ll +   p_rr)
  kin_avg = 1/2 * (v1_ll*v1_rr + v2_ll*v2_rr)

  # Calculate fluxes depending on orientation
  if orientation == 1
    pv1_avg = 1/2 * (p_ll*v1_rr + p_rr*v1_ll)
    f1 = rho_avg * v1_avg
    f2 = rho_avg * v1_avg * v1_avg + p_avg
    f3 = rho_avg * v1_avg * v2_avg
    f4 = p_avg*v1_avg/(equations.gamma-1) + rho_avg*v1_avg*kin_avg + pv1_avg
  else
    pv2_avg = 1/2 * (p_ll*v2_rr + p_rr*v2_ll)
    f1 = rho_avg * v2_avg
    f2 = rho_avg * v2_avg * v1_avg
    f3 = rho_avg * v2_avg * v2_avg + p_avg
    f4 = p_avg*v2_avg/(equations.gamma-1) + rho_avg*v2_avg*kin_avg + pv2_avg
  end

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

"""
    flux_kennedy_gruber(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D)

Kinetic energy preserving two-point flux by
- Kennedy and Gruber (2008)
  Reduced aliasing formulations of the convective terms within the
  Navier-Stokes equations for a compressible fluid
  [DOI: 10.1016/j.jcp.2007.09.020](https://doi.org/10.1016/j.jcp.2007.09.020)
"""
@inline function flux_kennedy_gruber(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations2D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_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

  # Average each factor of products in flux
  rho_avg = 1/2 * (rho_ll + rho_rr)
  v1_avg = 1/2 * (v1_ll + v1_rr)
  v2_avg = 1/2 * (v2_ll + v2_rr)
  p_avg = 1/2 * ((equations.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2 + v2_ll^2)) +
                 (equations.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2 + v2_rr^2)))
  e_avg = 1/2 * (rho_e_ll/rho_ll + rho_e_rr/rho_rr)

  # Calculate fluxes depending on orientation
  if orientation == 1
    f1 = rho_avg * v1_avg
    f2 = rho_avg * v1_avg * v1_avg + p_avg
    f3 = rho_avg * v1_avg * v2_avg
    f4 = (rho_avg * e_avg + p_avg) * v1_avg
  else
    f1 = rho_avg * v2_avg
    f2 = rho_avg * v2_avg * v1_avg
    f3 = rho_avg * v2_avg * v2_avg + p_avg
    f4 = (rho_avg * e_avg + p_avg) * v2_avg
  end

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

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

Entropy conserving two-point flux by
- Chandrashekar (2013)
  Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes
  for Compressible Euler and Navier-Stokes Equations
  [DOI: 10.4208/cicp.170712.010313a](https://doi.org/10.4208/cicp.170712.010313a)
"""
@inline function flux_chandrashekar(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations2D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_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
  p_ll =  (equations.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2 + v2_ll^2))
  p_rr =  (equations.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2 + v2_rr^2))
  beta_ll = 0.5*rho_ll/p_ll
  beta_rr = 0.5*rho_rr/p_rr
  specific_kin_ll = 0.5*(v1_ll^2 + v2_ll^2)
  specific_kin_rr = 0.5*(v1_rr^2 + v2_rr^2)

  # Compute the necessary mean values
  rho_avg  = 0.5*(rho_ll+rho_rr)
  rho_mean = ln_mean(rho_ll,rho_rr)
  beta_mean = ln_mean(beta_ll,beta_rr)
  beta_avg = 0.5*(beta_ll+beta_rr)
  v1_avg = 0.5*(v1_ll+v1_rr)
  v2_avg = 0.5*(v2_ll+v2_rr)
  p_mean = 0.5*rho_avg/beta_avg
  velocity_square_avg = specific_kin_ll + specific_kin_rr

  # Calculate fluxes depending on orientation
  if orientation == 1
    f1 = rho_mean * v1_avg
    f2 = f1 * v1_avg + p_mean
    f3 = f1 * v2_avg
    f4 = f1 * 0.5*(1/(equations.gamma-1)/beta_mean - velocity_square_avg)+f2*v1_avg + f3*v2_avg
  else
    f1 = rho_mean * v2_avg
    f2 = f1 * v1_avg
    f3 = f1 * v2_avg + p_mean
    f4 = f1 * 0.5*(1/(equations.gamma-1)/beta_mean - velocity_square_avg)+f2*v1_avg + f3*v2_avg
  end

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

"""
    flux_ranocha(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D)

Entropy conserving and kinetic energy preserving two-point flux by
- Ranocha (2018)
  Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods
  for Hyperbolic Balance Laws
  [PhD thesis, TU Braunschweig](https://cuvillier.de/en/shop/publications/7743)
See also
- Ranocha (2020)
  Entropy Conserving and Kinetic Energy Preserving Numerical Methods for
  the Euler Equations Using Summation-by-Parts Operators
  [Proceedings of ICOSAHOM 2018](https://doi.org/10.1007/978-3-030-39647-3_42)
"""
@inline function flux_ranocha(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations2D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_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
  p_ll =  (equations.gamma - 1) * (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2))
  p_rr =  (equations.gamma - 1) * (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2))

  # Compute the necessary mean values
  rho_mean   = ln_mean(rho_ll, rho_rr)
  rho_p_mean = ln_mean(rho_ll / p_ll, rho_rr / p_rr)
  v1_avg = 0.5 * (v1_ll + v1_rr)
  v2_avg = 0.5 * (v2_ll + v2_rr)
  p_avg  = 0.5 * (p_ll + p_rr)
  velocity_square_avg = 0.5 * (v1_ll*v1_rr + v2_ll*v2_rr)

  # Calculate fluxes depending on orientation
  if orientation == 1
    f1 = rho_mean * v1_avg
    f2 = f1 * v1_avg + p_avg
    f3 = f1 * v2_avg
    f4 = f1 * ( velocity_square_avg + 1 / ((equations.gamma-1) * rho_p_mean) ) + 0.5 * (p_ll*v1_rr + p_rr*v1_ll)
  else
    f1 = rho_mean * v2_avg
    f2 = f1 * v1_avg
    f3 = f1 * v2_avg + p_avg
    f4 = f1 * ( velocity_square_avg + 1 / ((equations.gamma-1) * rho_p_mean) ) + 0.5 * (p_ll*v2_rr + p_rr*v2_ll)
  end

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

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the
# maximum velocity magnitude plus the maximum speed of sound
# TODO: This doesn't really use the `orientation` - should it?
@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations2D)
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  # Calculate primitive variables and speed of sound
  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 = (equations.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * v_mag_ll^2)
  c_ll = sqrt(equations.gamma * 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 = (equations.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * v_mag_rr^2)
  c_rr = sqrt(equations.gamma * p_rr / rho_rr)

  λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
end

@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEquations2D)
  return max_abs_speed_naive(u_ll, u_rr, 0, equations) * norm(normal_direction)
end

# Calculate minimum and maximum wave speeds for HLL-type fluxes
@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations2D)
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  # Calculate primitive variables and speed of sound
  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  p_ll = (equations.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2 + v2_ll^2))

  v1_rr = rho_v1_rr / rho_rr
  v2_rr = rho_v2_rr / rho_rr
  p_rr = (equations.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2 + v2_rr^2))

  if orientation == 1 # x-direction
    λ_min = v1_ll - sqrt(equations.gamma * p_ll / rho_ll)
    λ_max = v1_rr + sqrt(equations.gamma * p_rr / rho_rr)
  else # y-direction
    λ_min = v2_ll - sqrt(equations.gamma * p_ll / rho_ll)
    λ_max = v2_rr + sqrt(equations.gamma * p_rr / rho_rr)
  end

  return λ_min, λ_max
end

@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
                                     equations::CompressibleEulerEquations2D)
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  # Calculate primitive variables and speed of sound
  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  p_ll = (equations.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2 + v2_ll^2))

  v1_rr = rho_v1_rr / rho_rr
  v2_rr = rho_v2_rr / rho_rr
  p_rr = (equations.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2 + v2_rr^2))

  v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
  v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]

  λ_min = ( v_normal_ll - sqrt(equations.gamma * p_ll / rho_ll) ) * norm(normal_direction)
  λ_max = ( v_normal_rr + sqrt(equations.gamma * p_rr / rho_rr) ) * norm(normal_direction)

  return λ_min, λ_max
end

# Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction
# has been normalized prior to this rotation of the state vector
@inline function rotate_to_x(u, normal_vector, equations::CompressibleEulerEquations2D)
  # cos and sin of the angle between the x-axis and the normalized normal_vector are
  # the normalized vector's x and y coordinates respectively (see unit circle).
  c = normal_vector[1]
  s = normal_vector[2]

  # Apply the 2D rotation matrix with normal and tangent directions of the form
  # [ 1    0    0   0;
  #   0   n_1  n_2  0;
  #   0   t_1  t_2  0;
  #   0    0    0   1 ]
  # where t_1 = -n_2 and t_2 = n_1

  return SVector(u[1],
                 c * u[2] + s * u[3],
                 -s * u[2] + c * u[3],
                 u[4])
end

# Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction
# has been normalized prior to this back-rotation of the state vector
@inline function rotate_from_x(u, normal_vector, equations::CompressibleEulerEquations2D)
  # cos and sin of the angle between the x-axis and the normalized normal_vector are
  # the normalized vector's x and y coordinates respectively (see unit circle).
  c = normal_vector[1]
  s = normal_vector[2]

  # Apply the 2D back-rotation matrix with normal and tangent directions of the form
  # [ 1    0    0   0;
  #   0   n_1  t_1  0;
  #   0   n_2  t_2  0;
  #   0    0    0   1 ]
  # where t_1 = -n_2 and t_2 = n_1

  return SVector(u[1],
                 c * u[2] - s * u[3],
                 s * u[2] + c * u[3],
                 u[4])
end

 """
    flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D)

Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro
[Lecture slides](http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf)
Signal speeds: [DOI: 10.1137/S1064827593260140](https://doi.org/10.1137/S1064827593260140)
"""
function flux_hllc(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations2D)
  # Calculate primitive variables and speed of sound
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  e_ll  = rho_e_ll / rho_ll
  p_ll = (equations.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2 + v2_ll^2))
  c_ll = sqrt(equations.gamma*p_ll/rho_ll)

  v1_rr = rho_v1_rr / rho_rr
  v2_rr = rho_v2_rr / rho_rr
  e_rr  = rho_e_rr / rho_rr
  p_rr = (equations.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2 + v2_rr^2))
  c_rr = sqrt(equations.gamma*p_rr/rho_rr)

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

  # Compute Roe averages
  sqrt_rho_ll = sqrt(rho_ll)
  sqrt_rho_rr = sqrt(rho_rr)
  sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr
  if orientation == 1 # x-direction
    vel_L = v1_ll
    vel_R = v1_rr
    ekin_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2
  elseif orientation == 2 # y-direction
    vel_L = v2_ll
    vel_R = v2_rr
    ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2
  end
  vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho
  ekin_roe = 0.5 * (vel_roe^2 + ekin_roe / sum_sqrt_rho^2)
  H_ll = (rho_e_ll + p_ll) / rho_ll
  H_rr = (rho_e_rr + p_rr) / rho_rr
  H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho
  c_roe = sqrt((equations.gamma - 1) * (H_roe - ekin_roe))
  Ssl = min(vel_L - c_ll, vel_roe - c_roe)
  Ssr = max(vel_R + c_rr, vel_roe + c_roe)
  sMu_L = Ssl - vel_L
  sMu_R = Ssr - vel_R

  if Ssl >= 0.0
    f1 = f_ll[1]
    f2 = f_ll[2]
    f3 = f_ll[3]
    f4 = f_ll[4]
  elseif Ssr <= 0.0
    f1 = f_rr[1]
    f2 = f_rr[2]
    f3 = f_rr[3]
    f4 = f_rr[4]
  else
    SStar = (p_rr - p_ll + rho_ll*vel_L*sMu_L - rho_rr*vel_R*sMu_R) / (rho_ll*sMu_L - rho_rr*sMu_R)
    if Ssl <= 0.0 <= SStar
      densStar = rho_ll*sMu_L / (Ssl-SStar)
      enerStar = e_ll + (SStar - vel_L) * (SStar + p_ll / (rho_ll * sMu_L))
      UStar1 = densStar
      UStar4 = densStar*enerStar
      if orientation == 1 # x-direction
        UStar2 = densStar*SStar
        UStar3 = densStar*v2_ll
      elseif orientation == 2 # y-direction
        UStar2 = densStar*v1_ll
        UStar3 = densStar*SStar
      end
      f1 = f_ll[1]+Ssl*(UStar1 - rho_ll)
      f2 = f_ll[2]+Ssl*(UStar2 - rho_v1_ll)
      f3 = f_ll[3]+Ssl*(UStar3 - rho_v2_ll)
      f4 = f_ll[4]+Ssl*(UStar4 - rho_e_ll)
    else
      densStar = rho_rr*sMu_R / (Ssr-SStar)
      enerStar = e_rr + (SStar - vel_R) * (SStar + p_rr / (rho_rr * sMu_R))
      UStar1 = densStar
      UStar4 = densStar*enerStar
      if orientation == 1 # x-direction
        UStar2 = densStar*SStar
        UStar3 = densStar*v2_rr
      elseif orientation == 2 # y-direction
        UStar2 = densStar*v1_rr
        UStar3 = densStar*SStar
      end
      f1 = f_rr[1]+Ssr*(UStar1 - rho_rr)
      f2 = f_rr[2]+Ssr*(UStar2 - rho_v1_rr)
      f3 = f_rr[3]+Ssr*(UStar3 - rho_v2_rr)
      f4 = f_rr[4]+Ssr*(UStar4 - rho_e_rr)
    end
  end
  return SVector(f1, f2, f3, f4)
end

@inline function max_abs_speeds(u, equations::CompressibleEulerEquations2D)
  rho, rho_v1, rho_v2, rho_e = u
  v1 = rho_v1 / rho
  v2 = rho_v2 / rho
  p = (equations.gamma - 1) * (rho_e - 1/2 * rho * (v1^2 + v2^2))
  c = sqrt(equations.gamma * p / rho)

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

# Convert conservative variables to primitive
@inline function cons2prim(u, equations::CompressibleEulerEquations2D)
  rho, rho_v1, rho_v2, rho_e = u

  v1 = rho_v1 / rho
  v2 = rho_v2 / rho
  p = (equations.gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2))

  return SVector(rho, v1, v2, p)
end

# Convert conservative variables to entropy
@inline function cons2entropy(u, equations::CompressibleEulerEquations2D)
  rho, rho_v1, rho_v2, rho_e = u

  v1 = rho_v1 / rho
  v2 = rho_v2 / rho
  v_square = v1^2 + v2^2
  p = (equations.gamma - 1) * (rho_e - 0.5 * rho * v_square)
  s = log(p) - equations.gamma*log(rho)
  rho_p = rho / p

  w1 = (equations.gamma - s) / (equations.gamma-1) - 0.5 * rho_p * v_square
  w2 = rho_p * v1
  w3 = rho_p * v2
  w4 = -rho_p

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

@inline function entropy2cons(w, equations::CompressibleEulerEquations2D)
  # See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD
  # [DOI: 10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1)
  @unpack gamma = equations

  # convert to entropy `-rho * s` used by Hughes, France, Mallet (1986)
  # instead of `-rho * s / (gamma - 1)`
  V1, V2, V3, V5 = w .* (gamma-1)

  # s = specific entropy, eq. (53)
  s = gamma - V1 + (V2^2 + V3^2)/(2*V5)

  # eq. (52)
  rho_iota = ((gamma-1) / (-V5)^gamma)^(1/(gamma-1))*exp(-s/(gamma-1))

  # eq. (51)
  rho      = -rho_iota * V5
  rho_v1   =  rho_iota * V2
  rho_v2   =  rho_iota * V3
  rho_e    =  rho_iota * (1-(V2^2 + V3^2)/(2*V5))
  return SVector(rho, rho_v1, rho_v2, rho_e)
end

# Convert primitive to conservative variables
@inline function prim2cons(prim, equations::CompressibleEulerEquations2D)
  rho, v1, v2, p = prim
  rho_v1 = rho * v1
  rho_v2 = rho * v2
  rho_e  = p/(equations.gamma-1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2)
  return SVector(rho, rho_v1, rho_v2, rho_e)
end

@inline function density(u, equations::CompressibleEulerEquations2D)
 rho = u[1]
 return rho
end

@inline function pressure(u, equations::CompressibleEulerEquations2D)
 rho, rho_v1, rho_v2, rho_e = u
 p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho)
 return p
end

@inline function density_pressure(u, equations::CompressibleEulerEquations2D)
 rho, rho_v1, rho_v2, rho_e = u
 rho_times_p = (equations.gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2))
 return rho_times_p
end

# Calculates the entropy flux in direction "orientation" and the entropy variables for a state cons
# NOTE: This method seems to work currently (b82534e) but is never used anywhere. Thus it is
# commented here until someone uses it or writes a test for it.
# @inline function cons2entropyvars_and_flux(gamma::Float64, cons, orientation::Int)
#   entropy = MVector{4, Float64}(undef)
#   v = (cons[2] / cons[1] , cons[3] / cons[1])
#   v_square= v[1]*v[1]+v[2]*v[2]
#   p = (gamma - 1) * (cons[4] - 1/2 * (cons[2] * v[1] + cons[3] * v[2]))
#   rho_p = cons[1] / p
#   # thermodynamic entropy
#   s = log(p) - gamma*log(cons[1])
#   # mathematical entropy
#   S = - s*cons[1]/(gamma-1)
#   # entropy variables
#   entropy[1] = (gamma - s)/(gamma-1) - 0.5*rho_p*v_square
#   entropy[2] = rho_p*v[1]
#   entropy[3] = rho_p*v[2]
#   entropy[4] = -rho_p
#   # entropy flux
#   entropy_flux = S*v[orientation]
#   return entropy, entropy_flux
# end

# Calculate thermodynamic entropy for a conservative state `cons`
@inline function entropy_thermodynamic(cons, equations::CompressibleEulerEquations2D)
  # Pressure
  p = (equations.gamma - 1) * (cons[4] - 1/2 * (cons[2]^2 + cons[3]^2) / cons[1])

  # Thermodynamic entropy
  s = log(p) - equations.gamma*log(cons[1])

  return s
end

# Calculate mathematical entropy for a conservative state `cons`
@inline function entropy_math(cons, equations::CompressibleEulerEquations2D)
  # Mathematical entropy
  S = -entropy_thermodynamic(cons, equations) * cons[1] / (equations.gamma - 1)

  return S
end

# Default entropy is the mathematical entropy
@inline entropy(cons, equations::CompressibleEulerEquations2D) = entropy_math(cons, equations)

# Calculate total energy for a conservative state `cons`
@inline energy_total(cons, ::CompressibleEulerEquations2D) = cons[4]

# Calculate kinetic energy for a conservative state `cons`
@inline function energy_kinetic(u, equations::CompressibleEulerEquations2D)
  rho, rho_v1, rho_v2, rho_e = u
  return (rho_v1^2 + rho_v2^2) / (2 * rho)
end

# Calculate internal energy for a conservative state `cons`
@inline function energy_internal(cons, equations::CompressibleEulerEquations2D)
  return energy_total(cons, equations) - energy_kinetic(cons, equations)
end


text:
export FactoredPolynomial

"""
    FactoredPolynomial{T,X}

A polynomial type that stores its data in a dictionary whose keys are the roots and whose values are the respective multiplicities along with a leading coefficient.

The structure is utilized for scalar multiplication, polynomial multiplication and powers, the finding of roots, and the identification of a greatest common divisor. For other operations, say addition, the operation is done after converting to the `Polynomial` type then converting back. (This requires the identification of the roots, so is subject to numeric issues.)

## Examples

```jldoctest factored_polynomial
julia> using Polynomials

julia> p = FactoredPolynomial(Dict([0=>1, 1=>2, 3=>4]))
FactoredPolynomial(x * (x - 3)⁴ * (x - 1)²)

julia> q = fromroots(FactoredPolynomial, [0,1,2,3])
FactoredPolynomial(x * (x - 2) * (x - 3) * (x - 1))

julia> p*q
FactoredPolynomial(x² * (x - 2) * (x - 3)⁵ * (x - 1)³)

julia> p^1000
FactoredPolynomial(x¹⁰⁰⁰ * (x - 3)⁴⁰⁰⁰ * (x - 1)²⁰⁰⁰)

julia> gcd(p,q)
FactoredPolynomial(x * (x - 3) * (x - 1))

julia> p = Polynomial([24, -50, 35, -10, 1])
Polynomial(24 - 50*x + 35*x^2 - 10*x^3 + x^4)

julia> q = convert(FactoredPolynomial, p) # noisy form of `factor`:
FactoredPolynomial((x - 4.0000000000000036) * (x - 2.9999999999999942) * (x - 1.0000000000000002) * (x - 2.0000000000000018))

julia> map(x->round(x, digits=12), q) # map works over factors and leading coefficient -- not coefficients in the standard basis
FactoredPolynomial((x - 4.0) * (x - 2.0) * (x - 3.0) * (x - 1.0))
```
"""
struct FactoredPolynomial{T <: Number, X} <: StandardBasisPolynomial{T, X}
    coeffs::Dict{T,Int}
    c::T
    function FactoredPolynomial{T, X}(checked::Val{false}, cs::Dict{T,Int}, c::T) where {T, X}
        new{T,X}(cs,T(c))
    end
    function FactoredPolynomial{T, X}(cs::Dict{T,Int}, c=one(T)) where {T, X}
        D = Dict{T,Int}()
        for (k,v) ∈ cs
            v > 0 && (D[k] = v)
        end
        FactoredPolynomial{T,X}(Val(false), D,T(c))
    end
    function FactoredPolynomial(cs::Dict{T,Int}, c::S=1, var::SymbolLike=:x) where {T,S}
        X = Symbol(var)
        R = promote_type(T,S)
        D = convert(Dict{R,Int}, cs)
        FactoredPolynomial{R,X}(D, R(c))
    end
end

# There are idiosyncracies with this type
# * unlike P{T}(...) we allow T to widen here if the roots of the polynomial Polynomial(coeffs) needs
#   a wider type (e.g. Complex{Float64}, not Float64)
# * the handling of Inf and NaN, when a specified coefficient, is to just return a constant poly (Inf or NaN)
function FactoredPolynomial{T,X}(coeffs::AbstractVector{S}) where {T,S,X}
    p = Polynomial{T,X}(T.(coeffs))
    P = FactoredPolynomial{T,X}
    iszero(p)     && return zero(P)
    isconstant(p) && return FactoredPolynomial{T,X}(Val(false), Dict{T,Int}(), T(coeffs[1]))
    any(isnan, p) && return FactoredPolynomial{T,X}(Val(false), Dict{T,Int}(), T(NaN))
    any(isinf, p) && return FactoredPolynomial{T,X}(Val(false), Dict{T,Int}(), T(Inf))

    zs = Multroot.multroot(p)
    c = p[end]
    D = Dict(zip(zs.values, zs.multiplicities))
    FactoredPolynomial(D, c, X)
end

function FactoredPolynomial{T}(coeffs::AbstractVector{S}, var::SymbolLike=:x) where {T,S}
    X = Symbol(var)
    FactoredPolynomial{T,X}(coeffs)
end

function FactoredPolynomial(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T}
    X = Symbol(var)
    FactoredPolynomial{T,X}(coeffs)
end

## ----
# abstract.jl The use of @register FactoredPolynomial didn't quite work, so
# that is replicated here and modified.
Base.convert(::Type{P}, p::P) where {P <: FactoredPolynomial} = p
function Base.convert(P::Type{<:FactoredPolynomial}, p::FactoredPolynomial{T,X}) where {T,X}
    T′ = _eltype(P)
    𝑻 = T′ == nothing ? T : T′
    𝑿 = indeterminate(P, p)
    d = Dict{𝑻,Int}()
    copy!(d, p.coeffs)
    FactoredPolynomial{𝑻,𝑿}(d, p.c)
end
Base.promote(p::P,q::Q) where {X,T,P<:FactoredPolynomial{T,X},Q<:FactoredPolynomial{T,X}} = p,q
Base.promote_rule(::Type{<:FactoredPolynomial{T,X}}, ::Type{<:FactoredPolynomial{S,X}}) where {T,S,X} =
    FactoredPolynomial{promote_type(T,S), X}
Base.promote_rule(::Type{<:FactoredPolynomial{T,X}}, ::Type{S}) where {T,S<:Number,X} =
    FactoredPolynomial{promote_type(T,S), X}
FactoredPolynomial{T,X}(n::S) where {T,X,S<:Number} = T(n) * one(FactoredPolynomial{T,X})
FactoredPolynomial{T}(n::S, var::SymbolLike=:x) where {T,S<:Number} = T(n) * one(FactoredPolynomial{T,X})
FactoredPolynomial(n::S, var::SymbolLike=:x) where {S<:Number} = n * one(FactoredPolynomial{S,Symbol(var)})
FactoredPolynomial(var::SymbolLike=:x) = variable(FactoredPolynomial, Symbol(var))
(p::FactoredPolynomial)(x) = evalpoly(x, p)

function Base.convert(::Type{<:Polynomial}, p::FactoredPolynomial{T,X}) where {T,X}
    x = variable(Polynomial{T,X})
    isconstant(p) && return Polynomial{T,X}(p.c)
    p(x)
end
function Base.convert(P::Type{<:FactoredPolynomial}, p::Polynomial{T,X}) where {T,X}
    isconstant(p) && return ⟒(P)(constantterm(p), X)
    ⟒(P)(coeffs(p), X)
end

## ----
## apply map to factors and the leading coefficient, not the coefficients
function Base.map(fn, p::P, args...)  where {T,X,P<:FactoredPolynomial{T,X}}
    𝒅 = Dict{T, Int}()
    for (k,v) ∈ p.coeffs
        𝒌 = fn(k, args...)
        𝒅[𝒌] = v
    end
    𝒄 = fn(p.c, args...)
    P(𝒅,𝒄)
end

## ----

function printpoly(io::IO, p::FactoredPolynomial{T,X}, mimetype=nothing) where {T,X}
    if iszero(p.c)
        print(io, p.c)
    else
        isone(p.c) && iszero(length(p.coeffs)) && print(io, p.c)
        !isone(p.c) && print(io, p.c)
        x = string(X)
        for (i,(k,v)) ∈ enumerate(p.coeffs)
            (!isone(p.c) || i != 1) && print(io, " * ")
            if iszero(k)
                print(io, x)
            else
                print(io, "(")
                print(io, x)
                if hasneg(T)
                    isneg(k) ?  print(io, " + ", -k) : print(io, " - ", k)
                else
                    print(io, " - ", k)
                end
                print(io, ")")
            end
            v > 1 && unicode_exponent(io, v)
        end
    end
end

## ----
Base.lastindex(p::FactoredPolynomial) = degree(p)
Base.eachindex(p::FactoredPolynomial) = 0:degree(p)

function Base.getindex(p::FactoredPolynomial{T}, idx::Int) where {T <: Number}
    m,M = firstindex(p), lastindex(p)
    idx < m && throw(BoundsError(p, idx))
    idx > M && return zero(T)
    coeffs(p)[idx - m + 1]
end

# pairs,keys, values
Base.keys(p::FactoredPolynomial)   = keys(convert(Polynomial, p))
Base.values(p::FactoredPolynomial) = values(convert(Polynomial, p))
Base.pairs(p::FactoredPolynomial)  = pairs(convert(Polynomial, p))

Base.copy(p::P) where {P<:FactoredPolynomial} = P(copy(p.coeffs), p.c)
Base.:(==)(p1::FactoredPolynomial, p2::FactoredPolynomial) =
    check_same_variable(p1,p2) && p1.c == p2.c && (p1.coeffs == p2.coeffs)

# what does it mean to be approximate here?
# Question pushed off to the Polynomial type
# might be better to compare roots and multiplicities
function Base.isapprox(p1::FactoredPolynomial{T,X},
                       p2::FactoredPolynomial{S,Y};
                       rtol::Real = (Base.rtoldefault(T, S, 0)),
                       atol::Real = 0,) where {T,X,S,Y}
    assert_same_variable(p1, p2)

    # p ≈ q if 𝒑 ≈ 𝒒
    𝑷 = Polynomial
    𝒑,𝒒 = convert(𝑷,p1), convert(𝑷,p2)
    return isapprox(𝒑, 𝒒, atol=atol, rtol=rtol)

end

## ----
Base.iszero(p::FactoredPolynomial) = iszero(p.c)

Base.zero(::Type{FactoredPolynomial{T,X}}) where {T, X} = FactoredPolynomial{T,X}(Dict{T,Int}(),     zero(T))
Base.one(::Type{FactoredPolynomial{T,X}}) where {T, X}  = FactoredPolynomial{T,X}(Dict{T,Int}(),     one(T))
variable(::Type{FactoredPolynomial{T,X}}) where {T, X}  = FactoredPolynomial{T,X}(Dict{T,Int}(0=>1), one(T))

#
function evalpoly(x, p::FactoredPolynomial)
    iszero(length(p.coeffs)) && return p.c * one(x)
    p.c * prod((x - k)^v for (k,v) ∈ p.coeffs)
end

function evalpoly(x::Array, p::FactoredPolynomial)
    iszero(length(p.coeffs)) && return p.c * one(x)
    p.c * prod((x - k*I)^v for (k,v) ∈ p.coeffs)
end

coeffs(p::FactoredPolynomial) = coeffs(convert(Polynomial, p))

function degree(p::P) where {P <: FactoredPolynomial}
    d = length(p.coeffs)
    d > 0 && return sum(values(p.coeffs))
    iszero(p.c)  ? -1 : 0
end

## ----

function fromroots(::Type{P}, r::AbstractVector{T}; var::SymbolLike=:x) where {T <: Number, P<:FactoredPolynomial}
    X = Symbol(var)
    d = Dict{T,Int}()
    for rᵢ ∈ r
        d[rᵢ] = get(d, rᵢ, 0) + 1
    end
    FactoredPolynomial{T, X}(d)
end

roots(p::FactoredPolynomial{T}) where {T} = Base.typed_vcat(T,[repeat([k],v) for (k,v) ∈ p.coeffs]...)

## -----
# unary subtraction
Base.:-(p::P) where {T,X,P<:FactoredPolynomial{T,X}} = (-1)*p

# addition
function Base.:+(p::P, c::S) where {S<:Number, T,X, P<:FactoredPolynomial{T,X}}
    R = promote_type(S,T)
    𝑷 = Polynomial{R,X}
    𝒑,𝒒 = convert(𝑷, p), convert(𝑷, c)
    convert(P, 𝒑 + 𝒒)
end
function Base.:+(p::P, q::P) where {T,X,P<:FactoredPolynomial{T,X}}
    𝑷 = Polynomial{T,X}
    𝒑,𝒒 = convert(𝑷, p), convert(𝑷, q)
    convert(P, 𝒑 + 𝒒 )
end

# multiplication
function Base.:*(p::P, q::P) where {T,X, P<:FactoredPolynomial{T,X}}
    d = copy(p.coeffs)
    for (k,v) ∈ q.coeffs
        d[k] = get(d, k, 0) + v
    end
    P(d, p.c*q.c)
end

# scalar mult
function scalar_mult(p::P, c::S) where {S<:Number, T, X, P <: FactoredPolynomial{T, X}}
    R = promote_type(T,S)
    d = Dict{R, Int}() # wident
    copy!(d, p.coeffs)
    FactoredPolynomial{R,X}(d, c * p.c)
end

# scalar division
function Base.:/(p::P, c::S) where {S, T, X, P <: FactoredPolynomial{T, X}}
    p * (1/c)
end

function Base.:^(p::P, n::Integer) where {T,X, P<:FactoredPolynomial{T,X}}
    n >= 0 || throw(ArgumentError("n must be non-negative"))
    d = Dict{T,Int}()
    for (k,v) ∈ p.coeffs
        d[k] = v*n
    end
    P(d, p.c^n)
end

## ----
## gcd, divrem, uvw
function Base.gcd(p::P, q::P) where {T, X, P<:FactoredPolynomial{T,X}}
    iszero(p) && return q
    iszero(q) && return p
    d = Dict{T,Int}()

    for k ∈ intersect(keys(p.coeffs), keys(q.coeffs))
        d[k] = min(p.coeffs[k], q.coeffs[k])
    end

    P(d)
end

# return u,v,w with p = u*v , q = u*w
function uvw(p::P, q::P; kwargs...) where {T, X, P<:FactoredPolynomial{T,X}}
    du, dv, dw = Dict{T,Int}(), Dict{T,Int}(), Dict{T,Int}()
    dp,dq = p.coeffs, q.coeffs
    kp,kq = keys(dp), keys(dq)

    for k ∈ setdiff(kp, kq)
        dv[k] = dp[k]
    end
    for k ∈ setdiff(kq, kp)
        dw[k] = dq[k]
    end
    for k ∈ intersect(kp, kq)
        pₖ,qₖ = dp[k], dq[k]
        m = min(pₖ, qₖ)
        du[k] = m
        dv[k] = pₖ - m;
        dw[k] = qₖ - m
    end
    P(du), P(dv, p.c), P(dw, q.c)
end

# return a,b with p = a + b*q
function Base.divrem(p::P, q::P) where {T, X, P<:FactoredPolynomial{T,X}}

    u, v, w = uvw(p, q)
    𝑷 = Polynomial
    𝒗,𝒘 = convert(𝑷, v), convert(𝑷, w)
    𝒅,𝒓 = divrem(𝒗, 𝒘)
    d,r = convert(P, 𝒅), convert(P, 𝒓)

    return (d, u*r)

end

## ----

function integrate(p::P) where {P <: FactoredPolynomial}
    𝑷 = Polynomial
    𝒑 = convert(𝑷, p)
    ∫𝒑 = integrate(𝒑)
    convert(P, ∫𝒑)
end

function derivative(p::P,n::Int) where {P <: FactoredPolynomial}
    𝑷 = Polynomial
    𝒑 = convert(𝑷, p)
    𝒑⁽ⁿ⁾ = derivative(𝒑, n)
    convert(P, 𝒑⁽ⁿ⁾)
end

ext: jl


text:
"""
    Result

Stores the results from an experiment.
"""
struct Result
    projectdirname::AbstractString
    experimentdirname::AbstractString
    replicate::Int
    samplename::AbstractString
    cell::Int
    channels::Vector{ChannelData}
    distances::Vector{AbstractFloat}
    mediandistance::AbstractFloat
    percentileranks::Vector{AbstractFloat}
    positivecontrol_distances::Vector{AbstractFloat} # stores the distances obtained when simulating 100% bound from this cell
    negativecontrol_distances::Vector{Vector{AbstractFloat}} # stores the distances obtained when simulating 0% bound from this cell
    positivecontrol_percentileranks::Vector{AbstractFloat} # stores the percentile ranks obtained when simulating 100% bound from this cell
    negativecontrol_percentileranks::Vector{Vector{AbstractFloat}} # stores the percentile ranks obtained when simulating 0% bound from this cell
end

Result(r::Result, newchannels::Vector{ChannelData}) = Result(r.projectdirname, r.experimentdirname, r.replicate, r.samplename, r.cell, newchannels, r.distances, r.mediandistance, r.percentileranks,
                                                             r.positivecontrol_distances, r.negativecontrol_distances, r.positivecontrol_percentileranks, r.negativecontrol_percentileranks)

"""
    ResultOptimizing

Stores the results from an experiment with extra data for optimization purposes.
"""
struct ResultOptimizing
    projectdirname::AbstractString
    experimentdirname::AbstractString
    replicate::Int
    samplename::AbstractString
    cell::Int
    channels::Vector{ChannelDataOptimizing}
    distances::Vector{AbstractFloat}
    mediandistance::AbstractFloat
    percentileranks_by_localdensity::Vector{Vector{Float64}}
end

"""
    ResultSimulate

Stores the results from an experiment, for the simulation method (without positive/negative controls)
"""
struct ResultSimulate
    projectdirname::AbstractString
    experimentdirname::AbstractString
    replicate::Int
    samplename::AbstractString
    cell::Int
    channels::Vector{ChannelData}
    distances::Vector{AbstractFloat}
    mediandistance::AbstractFloat
    percentileranks::Vector{AbstractFloat}
end


text:
function _profile(m::GeneralizedLinearModel, ind::Integer)
    pred = m.pp
    Xorig = pred.X
    remainder = ones(Bool, size(Xorig, 2))
    remainder[ind] = false
    Xreduced = Xorig[:, remainder]
    start = coef(m)[remainder]
    y = m.rr.y
    d = m.rr.d
    l = Link(m)
    β₀ = coef(m)[ind]
    col = view(Xorig, :, ind)
    stderror = StatsBase.stderr(m)[ind]
    dev₀ = deviance(m)
    parvals = [β₀]
    zvals = zeros(parvals)
    newm = fit(GeneralizedLinearModel, Xreduced, y, d, l, offset = β₀ * col, start = start)
    cvals = copy(coef(newm))
    function newval(pv)
        push!(parvals, pv)
        newm.fit = false
        copy!(newm.rr.offset, pv * col)
        push!(zvals, sign(step) * sqrt(deviance(fit!(newm, start = start)) - dev₀))
        append!(cvals, coef(newm))
    end
    fact = step = 0.6
    while abs(zvals[end]) < quantile(Normal(), 0.995)
        newval(β₀ + fact * stderror)
        fact += step
    end
    fact = step = -step
    newval(β₀ + fact * stderror)
    while abs(zvals[end]) < quantile(Normal(), 0.995)
        newval(β₀ + fact * stderror)
        fact += step
    end
    p = sortperm(parvals)
    parvals[p], zvals[p], transpose(reshape(cvals, (length(start), length(p))))[p,:]
end

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
function (rule::UnequalRule)(g::EGraph, match::Match)
    lc = match.id
    rinst = instantiate(g, match.pat_to_inst, match.sub, rule)
    rc, node = addexpr!(g, rinst)

    if find(g, lc) == find(g, rc)
        @log "Contradiction!" rule
        return Contradiction()
    end
    return nothing
end

function (rule::SymbolicRule)(g::EGraph, match::Match)
    rinst = instantiate(g, match.pat_to_inst, match.sub, rule)
    rc, node = addexpr!(g, rinst)
    merge!(g, match.id, rc.id)
    return nothing
end

function (rule::DynamicRule)(g::EGraph, match::Match)
    f = rule.rhs_fun
    actual_params = [instantiate(g, PatVar(v, i), match.sub, rule) for (i, v) in enumerate(rule.patvars)]
    r = f(geteclass(g, match.id), match.sub, g, actual_params...)
    rc, node = addexpr!(g, r)
    merge!(g, match.id, rc.id)
    return nothing
end

using .ReportReasons
function eqsat_apply!(g::EGraph, matches, rep::Report, params::SaturationParams)
    i = 0
    # println.(matches)
    for match ∈ matches
        i += 1

        if params.eclasslimit > 0 && g.numclasses > params.eclasslimit
            @log "E-GRAPH SIZEOUT"
            rep.reason = EClassLimit(params.eclasslimit)
            return
        end

        if reached(g, params.goal)
            @log "Goal reached"
            rep.reason = GoalReached()
            return
        end

        rule=match.rule
        # println("applying $rule")

        halt_reason = rule(g, match)
        if (halt_reason !== nothing)
            rep.reason = halt_reason
            return
        end

        # println(rule)
        # println(sub)
        # println(l); println(r)
        # display(egraph.classes); println()
    end
end


text:
module Tracers

end  # module


text:
@testset "minimal_is_common_interval_test" begin

    @test iscommon(interval(-27.0,-27.0)) === true

    @test iscommon(interval(-27.0, 0.0)) === true

    @test iscommon(interval(0.0,0.0)) === true

    @test iscommon(interval(-0.0,-0.0)) === true

    @test iscommon(interval(-0.0,0.0)) === true

    @test iscommon(interval(0.0,-0.0)) === true

    @test iscommon(interval(5.0, 12.4)) === true

    @test iscommon(interval(-0x1.FFFFFFFFFFFFFp1023, 0x1.FFFFFFFFFFFFFp1023)) === true

    @test iscommon(entireinterval()) === false

    @test iscommon(emptyinterval()) === false

    @test iscommon(interval(-Inf, 0.0)) === false

    @test iscommon(interval(0.0, Inf)) === false

end

@testset "minimal_is_common_interval_dec_test" begin

    @test iscommon(DecoratedInterval(interval(-27.0,-27.0), com)) === true

    @test iscommon(DecoratedInterval(interval(-27.0, 0.0), com)) === true

    @test iscommon(DecoratedInterval(interval(0.0,0.0), com)) === true

    @test iscommon(DecoratedInterval(interval(-0.0,-0.0), com)) === true

    @test iscommon(DecoratedInterval(interval(-0.0,0.0), com)) === true

    @test iscommon(DecoratedInterval(interval(0.0,-0.0), com)) === true

    @test iscommon(DecoratedInterval(interval(5.0, 12.4), com)) === true

    @test iscommon(DecoratedInterval(interval(-0x1.FFFFFFFFFFFFFp1023, 0x1.FFFFFFFFFFFFFp1023), com)) === true

    @test iscommon(DecoratedInterval(interval(-27.0,-27.0), trv)) === true

    @test iscommon(DecoratedInterval(interval(-27.0, 0.0), def)) === true

    @test iscommon(DecoratedInterval(interval(0.0,0.0), dac)) === true

    @test iscommon(DecoratedInterval(interval(-0.0,-0.0), trv)) === true

    @test iscommon(DecoratedInterval(interval(-0.0,0.0), def)) === true

    @test iscommon(DecoratedInterval(interval(0.0,-0.0), dac)) === true

    @test iscommon(DecoratedInterval(interval(5.0, 12.4), def)) === true

    @test iscommon(DecoratedInterval(interval(-0x1.FFFFFFFFFFFFFp1023, 0x1.FFFFFFFFFFFFFp1023), trv)) === true

    @test iscommon(DecoratedInterval(entireinterval(), dac)) === false

    @test iscommon(nai()) === false

    @test iscommon(DecoratedInterval(emptyinterval(), trv)) === false

    @test iscommon(DecoratedInterval(interval(-Inf, 0.0), trv)) === false

    @test iscommon(DecoratedInterval(interval(0.0, Inf), def)) === false

end

@testset "minimal_is_singleton_test" begin

    @test isthin(interval(-27.0,-27.0)) === true

    @test isthin(interval(-2.0, -2.0)) === true

    @test isthin(interval(12.0,12.0)) === true

    @test isthin(interval(17.1, 17.1)) === true

    @test isthin(interval(-0.0,-0.0)) === true

    @test isthin(interval(0.0,0.0)) === true

    @test isthin(interval(-0.0, 0.0)) === true

    @test isthin(interval(0.0, -0.0)) === true

    @test isthin(emptyinterval()) === false

    @test isthin(entireinterval()) === false

    @test isthin(interval(-1.0, 0.0)) === false

    @test isthin(interval(-1.0, -0.5)) === false

    @test isthin(interval(1.0, 2.0)) === false

    @test isthin(interval(-Inf,-0x1.FFFFFFFFFFFFFp1023)) === false

    @test isthin(interval(-1.0,Inf)) === false

end

@testset "minimal_is_singleton_dec_test" begin

    @test isthin(DecoratedInterval(interval(-27.0,-27.0), def)) === true

    @test isthin(DecoratedInterval(interval(-2.0, -2.0), trv)) === true

    @test isthin(DecoratedInterval(interval(12.0,12.0), dac)) === true

    @test isthin(DecoratedInterval(interval(17.1, 17.1), com)) === true

    @test isthin(DecoratedInterval(interval(-0.0,-0.0), def)) === true

    @test isthin(DecoratedInterval(interval(0.0,0.0), com)) === true

    @test isthin(DecoratedInterval(interval(-0.0, 0.0), def)) === true

    @test isthin(DecoratedInterval(interval(0.0, -0.0), dac)) === true

    @test isthin(DecoratedInterval(emptyinterval(), trv)) === false

    @test isthin(nai()) === false

    @test isthin(DecoratedInterval(entireinterval(), def)) === false

    @test isthin(DecoratedInterval(interval(-1.0, 0.0), dac)) === false

    @test isthin(DecoratedInterval(interval(-1.0, -0.5), com)) === false

    @test isthin(DecoratedInterval(interval(1.0, 2.0), def)) === false

    @test isthin(DecoratedInterval(interval(-Inf,-0x1.FFFFFFFFFFFFFp1023), dac)) === false

    @test isthin(DecoratedInterval(interval(-1.0,Inf), trv)) === false

end

@testset "minimal_is_member_test" begin

    @test in(-27.0, interval(-27.0,-27.0)) === true

    @test in(-27.0, interval(-27.0, 0.0)) === true

    @test in(-7.0, interval(-27.0, 0.0)) === true

    @test in(0.0, interval(-27.0, 0.0)) === true

    @test in(-0.0, interval(0.0,0.0)) === true

    @test in(0.0, interval(0.0,0.0)) === true

    @test in(0.0, interval(-0.0,-0.0)) === true

    @test in(0.0, interval(-0.0,0.0)) === true

    @test in(0.0, interval(0.0,-0.0)) === true

    @test in(5.0, interval(5.0, 12.4)) === true

    @test in(6.3, interval(5.0, 12.4)) === true

    @test in(12.4, interval(5.0, 12.4)) === true

    @test in(0.0, entireinterval()) === true

    @test in(5.0, entireinterval()) === true

    @test in(6.3, entireinterval()) === true

    @test in(12.4, entireinterval()) === true

    @test in(0x1.FFFFFFFFFFFFFp1023, entireinterval()) === true

    @test in(-0x1.FFFFFFFFFFFFFp1023, entireinterval()) === true

    @test in(0x1.0p-1022, entireinterval()) === true

    @test in(-0x1.0p-1022, entireinterval()) === true

    @test in(-71.0, interval(-27.0, 0.0)) === false

    @test in(0.1, interval(-27.0, 0.0)) === false

    @test in(-0.01, interval(0.0,0.0)) === false

    @test in(0.000001, interval(0.0,0.0)) === false

    @test in(111110.0, interval(-0.0,-0.0)) === false

    @test in(4.9, interval(5.0, 12.4)) === false

    @test in(-6.3, interval(5.0, 12.4)) === false

    @test in(0.0, emptyinterval()) === false

    @test in(-4535.3, emptyinterval()) === false

    @test in(-Inf, emptyinterval()) === false

    @test in(Inf, emptyinterval()) === false

    @test in(NaN, emptyinterval()) === false

    @test in(-Inf, entireinterval()) === false

    @test in(Inf, entireinterval()) === false

    @test in(NaN, entireinterval()) === false

end

@testset "minimal_is_member_dec_test" begin

    @test in(-27.0, DecoratedInterval(interval(-27.0,-27.0), trv)) === true

    @test in(-27.0, DecoratedInterval(interval(-27.0, 0.0), def)) === true

    @test in(-7.0, DecoratedInterval(interval(-27.0, 0.0), dac)) === true

    @test in(0.0, DecoratedInterval(interval(-27.0, 0.0), com)) === true

    @test in(-0.0, DecoratedInterval(interval(0.0,0.0), trv)) === true

    @test in(0.0, DecoratedInterval(interval(0.0,0.0), def)) === true

    @test in(0.0, DecoratedInterval(interval(-0.0,-0.0), dac)) === true

    @test in(0.0, DecoratedInterval(interval(-0.0,0.0), com)) === true

    @test in(0.0, DecoratedInterval(interval(0.0,-0.0), trv)) === true

    @test in(5.0, DecoratedInterval(interval(5.0, 12.4), def)) === true

    @test in(6.3, DecoratedInterval(interval(5.0, 12.4), dac)) === true

    @test in(12.4, DecoratedInterval(interval(5.0, 12.4), com)) === true

    @test in(0.0, DecoratedInterval(entireinterval(), trv)) === true

    @test in(5.0, DecoratedInterval(entireinterval(), def)) === true

    @test in(6.3, DecoratedInterval(entireinterval(), dac)) === true

    @test in(12.4, DecoratedInterval(entireinterval(), trv)) === true

    @test in(0x1.FFFFFFFFFFFFFp1023, DecoratedInterval(entireinterval(), def)) === true

    @test in(-0x1.FFFFFFFFFFFFFp1023, DecoratedInterval(entireinterval(), dac)) === true

    @test in(0x1.0p-1022, DecoratedInterval(entireinterval(), trv)) === true

    @test in(-0x1.0p-1022, DecoratedInterval(entireinterval(), def)) === true

    @test in(-71.0, DecoratedInterval(interval(-27.0, 0.0), trv)) === false

    @test in(0.1, DecoratedInterval(interval(-27.0, 0.0), def)) === false

    @test in(-0.01, DecoratedInterval(interval(0.0,0.0), dac)) === false

    @test in(0.000001, DecoratedInterval(interval(0.0,0.0), com)) === false

    @test in(111110.0, DecoratedInterval(interval(-0.0,-0.0), trv)) === false

    @test in(4.9, DecoratedInterval(interval(5.0, 12.4), def)) === false

    @test in(-6.3, DecoratedInterval(interval(5.0, 12.4), dac)) === false

    @test in(0.0, DecoratedInterval(emptyinterval(), trv)) === false

    @test in(0.0, DecoratedInterval(emptyinterval(), trv)) === false

    @test in(-4535.3, DecoratedInterval(emptyinterval(), trv)) === false

    @test in(-4535.3, DecoratedInterval(emptyinterval(), trv)) === false

    @test in(-Inf, DecoratedInterval(emptyinterval(), trv)) === false

    @test in(-Inf, nai()) === false

    @test in(Inf, DecoratedInterval(emptyinterval(), trv)) === false

    @test in(Inf, nai()) === false

    @test in(NaN, DecoratedInterval(emptyinterval(), trv)) === false

    @test in(NaN, nai()) === false

    @test in(-Inf, DecoratedInterval(entireinterval(), trv)) === false

    @test in(Inf, DecoratedInterval(entireinterval(), def)) === false

    @test in(NaN, DecoratedInterval(entireinterval(), dac)) === false

end


text:
# Algorithm 2.8
# Approximate the root of a function by constructing a parabola
# between three points and intersect it with y = 0 to find the
# next approximation.

function muller(f, x0, x1, x2, tol, n0)
	h1 = x1 - x0
	h2 = x2 - x1
	δ1 = (f(x1) - f(x0))/h1
	δ2 = (f(x2) - f(x1))/h2
	d = (δ2 - δ1)/(h2 + h1)
	for i ∈ 3:n0
		b = δ2 + h2 * d
		D = (b^2 -4f(x2)d)^(1/2)
		E = 0
		if abs(b - D) < abs(b + D)
			E = b + D
		else
			E = b - D
		end
		h = (-2f(x2))/E
		p = x2 + h
		if abs(h) < tol
			return p
		end
		x0 = x1
		x1 = x2
		x2 = p
		h1 = x1 - x0
		h2 = x2 - x1
		δ1 = (f(x1) - f(x0))/h1
		δ2 = (f(x2) - f(x1))/h2
		d = (δ2 - δ1)/(h2 + h1)
	end
	throw(ArgumentError("Method failed"))
end

println(muller(x -> (16x^4 - 40x^3 + 5x^2 + 20x +6), 0.5, 1.0, 1.5, 10^-5, 50))


text:
module BNBKnapsack

    using LinearAlgebra
    using JuMP, CPLEX
    import Base.print

    mutable struct Problem
        n_objects::Int64        # Also the number of variables.
        values::Array{Int64,1}
        weights::Array{Int64,1}
        max_weight::Int64
        integer_precision::Float64  # Integers are represented as floating-point values: if they are close enough to an integer, they are considered integer.
        gap_precision::Float64
    end

    problem = nothing

    mutable struct BBNode
        # Tree structure.
        parent                              #::Union(Nothing, BBNode)
        children::Array{BBNode, 1}
        # Branching information.
        variable::Int                       # Just an index of the variable for the branching. Indexing starts at 1 (0: no branching; only for root node).
        variable_value::Int
        # Node's LP information.
        lp_value::Float64
        lp_feasibility                      #::Symbol # :Optimal or something else.
        solution::Array{Float64, 1}
        integer_feasible::Bool
    end

    mutable struct BBStatus
        upper_bound::Float64 # Computed from the LPs.
        lower_bound::Float64 # Incumbent.
        incumbent::Array{Int, 1} # Rounded from floating-point solution.
        root::BBNode
        # Optimisations: list of nodes to avoid going through all the tree, sort nodes so the most interesting ones are picked up first, etc.
    end

    function get_lp()
        m = Model(CPLEX.Optimizer)
        @variable(m, 0 <= x[1:problem.n_objects] <= 1)
        @constraint(m, dot(problem.weights, x) <= problem.max_weight)
        @objective(m, Max, dot(problem.values, x))
        return m, x
    end

    function solution_is_boolean(sol::Array{Float64,1})
        for s in sol
            if s > problem.integer_precision && s < 1 - problem.integer_precision
                return false
            end
        end
        return true
    end

    function gap(status::BBStatus)
        # Will return NaN if no integer-feasible solution is found.
        gap = (status.upper_bound - status.lower_bound) / status.upper_bound
        if isnan(gap)
            gap = Inf
        end
        return gap
    end

    function is_pruned(node::BBNode, status::BBStatus; quiet::Bool = true)
        if node.lp_feasibility == :UnfeasibleSolution
            if !quiet
                println("Node pruned by infeasibility.")
            end
            return true, :Feasibility
        end
        if node.integer_feasible
            if !quiet println("Node pruned by optimality.") end
            return true, :Optimality
        end
        if node.lp_value < status.lower_bound
            if !quiet println("Node pruned by bounds.") end
            return true, :Bounds
        end
        return false, :None
    end

    function get_next_node(status)
        # Find a node which has no children and is not pruned.
        # Pruning is then used here, together with tree-exploration strategy!
        function explore_node(node::BBNode)
            pruned, why_pruned = is_pruned(node, status; quiet=false)
            if pruned
                return nothing
            end
            if length(node.children) == 0
                println("self no children: $(node.variable)")
                return node
            end
            # Binary variables: two children or zero, no other case.
            node_left = node.children[1]
            if node_left != nothing
                possible_left = explore_node(node_left)
            else
                possible_left = nothing
            end
            node_right = node.children[2]
            if node_right != nothing
                possible_right = explore_node(node_right)
            else
                possible_right = nothing
            end
            if possible_left != nothing
                return possible_left
            elseif possible_right != nothing
                return possible_right
            end
        end
        return explore_node(status.root)
    end

    function is_variable_branched_on(node::BBNode, var::Int)
        # Depends on the path followed from root to now.
        if node.variable == var
            return true
        end
        if node.variable == 0 # Root node.
            return false
        end
        return is_variable_branched_on(node.parent, var)
    end

    function find_branching_variable(sol::Array{Float64,1}, current_node::BBNode)
        # Only called on non-integer solutions.
        n_vars = length(sol)
        for s in 1:n_vars
            v = sol[s]
            if v > problem.integer_precision && v < 1 - problem.integer_precision && ! is_variable_branched_on(current_node, s)
                return s
            end
        end
        println(" -> Where shall I branch?")
    end

    function branch_lp(current_node::BBNode, new_variable::Int, value::Int)
        lp, vars = get_lp()
        # Retrieve the previous nodes' branching decisions.
        node = current_node
        while true
            if node.variable == 0
                break
            end
            @constraint(lp, vars[node.variable] == node.variable_value)
            node = node.parent
        end
        @constraint(lp, vars[new_variable] == value)
        return lp, vars
    end

    function create_node(parent::BBNode, new_variable::Int, value::Int)
        lp, vars = branch_lp(parent, new_variable, value)
        lp_status = optimize!(lp)           #lp_feasibility = optimize!(lp)

        println(lp)

        lp_feasibility = :UnfeasibleSolution
        if termination_status(lp) == MOI.OPTIMAL
            lp_feasibility = :FeasibleSolution
        end

        lp_value = NaN
        lp_solution = []
        int_feas = false

        if lp_feasibility == :FeasibleSolution
            lp_value = JuMP.objective_value(lp)
            lp_solution = JuMP.value.(vars)
            int_feas = solution_is_boolean(lp_solution)
        end

        return BBNode(parent, [], new_variable, value, lp_value, lp_feasibility, lp_solution, int_feas)
    end

    function round_solution(sol::Array{Float64, 1})
        return [(s < .5) ? 0 : 1 for s in sol]
    end

    function get_upper_bound(status::BBStatus)
        # As the tree is developed such that both children are always created at the same time,
        # retrieve the leaves' LP values, take the minimum of them.
        function get_leaves(node::BBNode)
            if length(node.children) == 0
                pruned, why = is_pruned(node, status)
                if node.integer_feasible || (pruned && why != :Bounds)
                    return Inf
                else
                    return node.lp_value
                end
            else
                left = get_leaves(node.children[1])
                right = get_leaves(node.children[2])
                return min(left, right)
            end
        end
        return get_leaves(status.root)
    end

    function print(tree::BBNode, status::BBStatus)
        function space(level::Int)
            return "  " ^ level
        end
        function sub_print(node::BBNode, level::Int)
            print(" . ")
            print(space(level))

            if node.variable != 0
                print("Variable " * string(node.variable) * "   :   " * string(node.variable_value) * "   ;   " * string(node.lp_value))
            else
                print("Root node:               " * string(node.lp_value))
            end

            if node.integer_feasible
                print("  (I)")
            else
                print("     ")
            end

            pruned, why_pruned = is_pruned(node, status; quiet=true)
            if pruned
                print("  (P: " * string(why_pruned) * ")")
            else
                print("     ")
            end

            print("\n")

            for child in node.children
                sub_print(child, level + 1)
            end
        end
        sub_print(tree, 2)
    end

    function print(status::BBStatus)
        println(" . Complete branch-and-bound status, including tree.")
        println(" . Upper bound: " * string(status.upper_bound))
        println(" . Lower bound: " * string(status.lower_bound))
        println(" . Incumbent: " * string(status.incumbent))
        println(" . Tree: ")
        print(status.root, status)
        print("\n")
    end

    function execute(p::Problem)
        global problem = p

        println("\n=====================\n")

        root_lp, root_vars = get_lp()
        #root_feasibility = optimize!(root_lp)
        root_status = optimize!(root_lp)

        root_feasibility = :UnfeasibleSolution
        if termination_status(root_lp) == MOI.OPTIMAL
            root_feasibility = :FeasibleSolution
        end

        if root_feasibility != :FeasibleSolution
            println("Root node could not be solved to optimality: $root_feasibility")
            exit(0)
        end

        root_solution = JuMP.value.(root_vars)
        root_obj = JuMP.objective_value(root_lp)
        if solution_is_boolean(root_solution)
            println("Solved at root node.")
            println(root_solution)
            println("Objective: $root_obj")
            return
        end

        # parent, children, variable, variable_value, lp_value, lp_feasibility, solution, integer_feasible
        root = BBNode(nothing, [], 0, 0, root_obj, root_feasibility, root_solution, false)
        # upper_bound, lower_bound, incumbent, root
        status = BBStatus(root_obj, -Inf, [], root)

        println("[0] Root node solved. LP value: $root_obj.")

        # Node exploring: do both children at once, while exploring current node. Not what is usually done, but easier
        # (no list of nodes to explore).
        n_nodes_explored = 0
        while gap(status) > problem.gap_precision

            n_nodes_explored += 1

            println("[ $n_nodes_explored ] Exploring one more node.")

            # Find a node to explore. There may be none.
            node_to_explore = 0
            if n_nodes_explored == 1
                node_to_explore = root
            else
                node_to_explore = get_next_node(status)
            end

            if node_to_explore == nothing && n_nodes_explored != 1
                println("No more branches to explore in the tree.")
                break
            end

            # Decide the branching variable.
            var_index = 0
            if n_nodes_explored == 1
                var_index = find_branching_variable(root.solution, root)
            else
                var_index = find_branching_variable(node_to_explore.solution, node_to_explore)
            end

            println("[ $n_nodes_explored ] Branching on variable $var_index")

            # Create new nodes with branching on this *binary* variable (with all constraints from root to this new branch).
            left = create_node(node_to_explore, var_index, 0)
            right = create_node(node_to_explore, var_index, 1)

            # Update the parent with the new children.
            node_to_explore.children = [left, right]

            # Update the bounds and, if necessary, update the incumbent.
            if left.integer_feasible && left.lp_value > status.lower_bound
                println("[ $n_nodes_explored ] Found a new incumbent. Value: $(left.lp_value)")
                status.lower_bound = left.lp_value
                status.incumbent = round_solution(left.solution)
            end

            if right.integer_feasible && right.lp_value > status.lower_bound
                println("[ $n_nodes_explored ] Found a new incumbent. Value: $(right.lp_value)")
                status.lower_bound = right.lp_value
                status.incumbent = round_solution(right.solution)
            end

            status.upper_bound = get_upper_bound(status)

            # End iteration.
            println("[ $n_nodes_explored ] Iteration done. Gap: $(gap(status)). Bound: [ $(status.lower_bound) ; $(status.upper_bound) ].")

        end
        print(status)
        println("\n")
    end

end # module


text:
abstract AbstractInterpolator

export AbstractInterpolator

@lazymod LinTrigInterpMod "lin_trig_interp.jl"
@lazymod LinQuadInterpMod "lin_quad_interp.jl"
@lazymod QuadTrigInterpMod "quad_trig_interp.jl"

dNmatrix() = error("Not implemented")
dNdxmatrix() = error("Not implemented")
Jmatrix() = error("Not implemented")
Nvec() = error("Not implemented")
get_area() = error("Not implemented")
mass_matrix() = error("Not implemented")
mass_matrix_big() = error("Not implemented")
@inline function inv2x2t!(J::Matrix{Float64})
    d = det2x2(J)
    J[1,1], J[2,2] = J[2,2]/d, J[1,1]/d
    J[1,2], J[2,1] = -J[2,1]/d, -J[1,2]/d
    return J
end

@inline function det2x2(J::Matrix{Float64})
    d = J[1,1]*J[2,2] - J[1,2]*J[2,1]
    return d
end

#=
function dNdxmatrix{T <: AbstractInterpolator}(interp::T, local_coords::Point2,
                    vertices::Vector{Int}, nodes::Vector{FENode2})

        dN = dNmatrix(interp, local_coords)
        J = Jmatrix(interp, vertices, nodes, dN)
        dNdx = dN * (inv(J)')

        return dNdx
end

function Jmatrix{T <: AbstractInterpolator}(::T,
                 vertices::Vector{Int}, nodes::Vector{FENode2},
                 dN::Matrix{Float64})

    J = zeros(2, 2)

    for row in 1:size(dN, 1)
        x = nodes[vertices[row]].coords[1]
        y = nodes[vertices[row]].coords[2]

        J[1, 1] += dN[row, 1] * x
        J[1, 2] += dN[row, 1] * y
        J[2, 1] += dN[row, 2] * x
        J[2, 2] += dN[row, 2] * y
    end

    return J
end
=#


text:
function hammingsequence(N::Int)
    if N < 1
        throw("Hamming sequence index must be a positive integer")
    end
    ham = Vector{BigInt}([1])
    base2, base3, base5 = 1, 1, 1
    next2, next3, next5 = BigInt(2), BigInt(3), BigInt(5)
    for _ in 1:N-1
        x = min(next2, next3, next5)
        push!(ham, x)
        next2 <= x && (base2 += 1; next2 = 2ham[base2])
        next3 <= x && (base3 += 1; next3 = 3ham[base3])
        next5 <= x && (base5 += 1; next5 = 5ham[base5])
    end
    ham
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#
# Advent of Code
# Day 1: No Time for a Taxicab
#
# My First Julia Program
#

println("Advent of Code")
println("Day 1: No Time for a Taxicab")
println("")

DIRECTIONS  = ["North", "East", "South", "West"]
DEBUG = false

function travel(inputDataBits, isPart2)
  # lets find where we go!
  location = [0, 0]
  direction = 0
  visitedLocations = Vector{Int}[]
  foundIt = false

  # a strange while loop
  state = start(inputDataBits)
  while !done(inputDataBits, state)
    (item, state) = next(inputDataBits, state)

    if DEBUG println("Instruction: $item") end

    # deal with the item, which should be a string and an integer
    turnDirection = lowercase(item[1])
    distance = parse(Int,item[2:end])

    # work out the new direction
    if turnDirection == 'l'
      if DEBUG println("Turning left") end
      direction -= 1
      if direction < 0
        direction += 4
      end
    else
      if DEBUG println("Turning right") end
      direction += 1
      if direction > 3
        direction -= 4
      end
    end

    if DEBUG @printf("Facing %s, moving %d\n", DIRECTIONS[direction+1], distance) end

    #
    # We want to keep track of *everywhere* we've been,
    # not just destinations
    #

    for counter in range(1, distance)
      if direction == 0
        location[2] += 1
      elseif direction == 1
        location[1] += 1
      elseif direction == 2
        location[2] -= 1
      else
        location[1] -= 1
      end

      # Check our list of places we've been,
      # and if we've already been there, set the found flag
      # otherwise, add the current location to the list of
      # past locations

      if isPart2
        if !in(location, visitedLocations)
          push!(visitedLocations, copy(location))
        else
          foundIt = true
          break
        end
      end
    end

    # we need to break out of the iterating loop too
    if foundIt
      break
    end
  end

  return location
end

# Data file
# Input data has been saved to this file
inputDataFileName = "day01-input.txt"

# Read the input data, parse it down to a list of "commands"
inputDataFile = open(inputDataFileName, "r")
inputData = readstring(inputDataFile)
close(inputDataFile)

# Part 1 test data
#inputData = "R2, L3"
#inputData = "R2, R2, R2"
#inputData=  "R5, L5, R5, R3"

# Part 2 test data
#inputData = "R8, R4, R4, R8"

# like python split, break the string into bits on commas
inputDataBits = split(inputData, ",")

# in-place map function to strip the white spaces off of each inputDataBit
map!(strip,inputDataBits)

part1Location = travel(inputDataBits, false)
println("In part 1, you end up at ", part1Location, " which is at a distance of ", abs(part1Location[1]) + abs(part1Location[2]), " blocks.")

part2Location = travel(inputDataBits, true)
println("In part 2, you end up at ", part2Location, " which is at a distance of ", abs(part2Location[1]) + abs(part2Location[2]), " blocks.")


text:
using PhysicalConstants.CODATA2018: c_0, h, k_B

# Boltzmann constant in cm/K
const k_cm = Float64((k_B * (1 / h) * (1 / (c_0 * 100))).val)

const k_mhz = k_cm * (c_0.val / 1e4)
# c, h, k_B = c_0.val, h.val, k_B.val

export
    c_0,
    k_B,
    h,
    k_cm,
    k_mhz,
    MHz2wavenumber

MHz2wavenumber(ν) = ν / (c_0.val / 1e4)


text:
# See p. 280 of Murphy's Machine Learning
# x_k1 = x_k - alpha * gr + mu * (x - x_previous)

# L should be function or any other callable
immutable MomentumGradientDescent{L} <: Optimizer
    mu::Float64
    linesearch!::L
end

#= uncomment for v0.8.0
MomentumGradientDescent(; mu::Real = 0.01, linesearch = LineSearches.HagerZhang()) =
  MomentumGradientDescent(Float64(mu), linesearch)
=#

Base.summary(::MomentumGradientDescent) = "Momentum Gradient Descent"

function MomentumGradientDescent(; mu::Real = 0.01, linesearch = LineSearches.HagerZhang())
    MomentumGradientDescent(Float64(mu), linesearch)
end

type MomentumGradientDescentState{T,N}
    x::Array{T,N}
    x_previous::Array{T,N}
    f_x_previous::T
    s::Array{T,N}
    @add_linesearch_fields()
end

function initial_state{T}(method::MomentumGradientDescent, options, d, initial_x::Array{T})
    value_gradient!(d, initial_x)

    MomentumGradientDescentState(copy(initial_x), # Maintain current state in state.x
                         copy(initial_x), # Maintain previous state in state.x_previous
                         T(NaN), # Store previous f in state.f_x_previous
                         similar(initial_x), # Maintain current search direction in state.s
                         @initial_linesearch()...) # Maintain a cache for line search results in state.lsr
end

function update_state!{T}(d, state::MomentumGradientDescentState{T}, method::MomentumGradientDescent)
    n = length(state.x)
    # Search direction is always the negative gradient
    @simd for i in 1:n
        @inbounds state.s[i] = -gradient(d, i)
    end

    # Determine the distance of movement along the search line
    lssuccess = perform_linesearch!(state, method, d)

    # Update current position
    @simd for i in 1:n
        # Need to move x into x_previous while using x_previous and creating "x_new"
        @inbounds tmp = state.x_previous[i]
        @inbounds state.x_previous[i] = state.x[i]
        @inbounds state.x[i] = state.x[i] + state.alpha * state.s[i] + method.mu * (state.x[i] - tmp)
    end
    lssuccess == false # break on linesearch error
end


text:
"""
```julia
    ConstantStep()
```

Use a constant step size of ``( √L̅ - √̅μ̅  ) / ( √L̅ + √̅μ̅  )``, where ``L`` is the maximum eigenvalue
of the Hessian and ``μ`` is the minimum eigenvalue of the Hessian.

See also: [`VariableStep`](@ref)
"""
mutable struct ConstantStep <: AbstractStep
    β

    function ConstantStep()
        return new( 1.0 )
    end
end

function init!( step::ConstantStep, L, μ )
    step.β = ( sqrt( L ) - sqrt( μ ) ) / ( sqrt( L ) + sqrt( μ ) )
end

struct ConstantStepIterator{T}
    β::T
end

iterate( cs::ConstantStep ) = (cs.β, ConstantStepIterator{typeof( cs.β )}( cs.β ))
iterate( cs::ConstantStep, step::ConstantStepIterator ) = (cs.β, step)

"""
```julia
    VariableStep()
```

Use a variable step size  computed by finding the roots of a simple polynomial.

See Algorithm II.1 in [1] for more information.

[1] S. Richter, C. N. Jones, and M. Morari, ‘Computational Complexity Certification for Real-Time MPC With Input
    Constraints Based on the Fast Gradient Method’, IEEE Transactions on Automatic Control, vol. 57, no. 6,
    pp. 1391–1403, 2012, doi: 10.1109/TAC.2011.2176389.

See also: [`ConstantStep`](@ref)
"""
mutable struct VariableStep <: AbstractStep
    L
    μ
    func

    function VariableStep()
        return new( 0.0, 0.0, nothing )
    end
end

function init!( step::VariableStep, L, μ )
    step.L = L
    step.μ = μ

    # This let block creates a local copy of the variables for the closure, which
    # speeds up the computations (https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured)
    step.func = let μ = μ, L = L
        (α, αₙ) -> ( 1 - αₙ ) * α^2 + μ/L * αₙ - αₙ^2
    end
end

mutable struct VariableStepIterator{T} <: AbstractStep
    α::T
    αₙ::T
end

function iterate( vs::VariableStep )
    α = sqrt( vs.μ / vs.L )

    step = VariableStepIterator{typeof( α )}( α, α )

    return iterate( vs, step )
end

function iterate( vs::VariableStep, step::VariableStepIterator )
    step.α = step.αₙ

    # Search for the zero crossing inside (0,1)
    step.αₙ = find_zero( αₙ -> vs.func( step.α, αₙ ), (0, 1) )

    β = ( step.α * ( 1 - step.α ) ) / ( step.α^2 + step.αₙ )

    return (β, step)
end


text:
abstract type AbstractPhysicsComponent end
abstract type AbstractTerm end

"""
    Physics

    Essentially a NamedTuple, but one that returns `nothing` by default.
    Example:
        physics = Physics(diffusion = NiceDiffusion())
        `physics.diffusion` returns NiceDiffusion()
        `physics.coriolis` returns nothing

    With great power comes great responsibility:
        Great for flexibility, just like NamedTuple, but potential
        errors pass silently (e.g, typos in keywords)
"""
struct Physics{𝒯} <: AbstractPhysicsComponent
    components::𝒯
end

function Physics(; kwargs...)
    components = (; collect(kwargs)...)
    return Physics(components)
end

"""
    Base module extensions
"""
function Base.getproperty(value::Physics, name::Symbol)
    # Physics default for everything is `nothing`
    components = getfield(value, :components, false)
    return get(components, name, nothing)
end

function Base.propertynames(value::Physics)
    components = getfield(value, :components, false)
    return keys(components)
end

function Base.show(io::IO, physics::Physics)
    color = :white
    printstyled(io, "Physics(\n", color = color)
    for name in propertynames(physics)
        property = getproperty(physics, name)
        print("  ")
        printstyled(property, color = Int(rand(UInt8)))
        print("\n")
    end
    printstyled(io, ")", color = color)
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
"""
    HyperRectangle{N,T}

Hyperrectangle in `N` dimensions given by `low_corner::SVector{N,T}` and
`high_corner::SVector{N,T}`
"""
struct HyperRectangle{N,T}
    low_corner::SVector{N,T}
    high_corner::SVector{N,T}
end

# promote
HyperRectangle(l::Tuple,h::Tuple) = HyperRectangle(SVector(l),SVector(h))
HyperRectangle(l::SVector,h::SVector) = HyperRectangle(promote(l,h)...)
# 1d case
HyperRectangle(a::Number,b::Number) = HyperRectangle(SVector(a),SVector(b))

Base.:(==)(h1::HyperRectangle, h2::HyperRectangle) = (h1.low_corner == h2.low_corner) && (h1.high_corner == h2.high_corner)
Base.isapprox(h1::HyperRectangle,h2::HyperRectangle) = isapprox(h1.low_corner,h2.low_corner) && isapprox(h1.high_corner,h2.high_corner)
Base.in(point,h::HyperRectangle) = all(h.low_corner .<= point .<= h.high_corner)

Base.eltype(h::HyperRectangle{N,T}) where {N,T}     = T
ambient_dimension(h::HyperRectangle{N}) where {N}   = N
geometric_dimension(h::HyperRectangle{N}) where {N} = N

low_corner(r::HyperRectangle) = r.low_corner
high_corner(r::HyperRectangle) = r.high_corner

diameter(cub::HyperRectangle) = norm(cub.high_corner .- cub.low_corner,2)

function bounding_box(data)
    low_corner  = first(data)
    high_corner = first(data)
    for pt in data
        low_corner  = min.(low_corner,pt)
        high_corner = max.(high_corner,pt)
    end
    return HyperRectangle(low_corner,high_corner)
end

center(rec::HyperRectangle) = (rec.low_corner + rec.high_corner) / 2
radius(rec::HyperRectangle) = diameter(rec) / 2

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
module JSON2DataFrame

    # ############################################################################
    #
    # Dependencies
    #
    # ############################################################################

    using JSON, DataFrames

    ##############################################################################
    ##
    ## Exported methods and types
    ##
    ##############################################################################

    export flatten,
        flatten_to_df,
        dicts_to_df

    ##############################################################################
    ##
    ## Load files
    ##
    ##############################################################################

    include("flatten.jl")

end