text:
module MOCNeutronTransport

# Compilation Options
# -------------------------------------------------------------------------------------------------
const path_to_gmsh_api = "/usr/local/lib/gmsh.jl"
const enable_visualization = false
const visualize_ray_tracing = false

# using
# -------------------------------------------------------------------------------------------------
using Logging
using HDF5
using LightXML
using LinearAlgebra
using StaticArrays
using Dates: now, format
using LoggingExtras: TransformerLogger, global_logger

# import
# -------------------------------------------------------------------------------------------------
# import Base: @propagate_inbounds
# import Base: broadcastable, getindex, getproperty, +, -, *, /, in, intersect,
#              isapprox, rand, union

# Optional compilation/local dependencies
# -------------------------------------------------------------------------------------------------
# if enable_visualization
#     using GLMakie: Axis, Figure, LineSegments, Mesh, Scatter, current_axis, record
#     import GLMakie: linesegments!, mesh!, scatter!, convert_arguments
# end
# Gmsh
include(path_to_gmsh_api)

# Setup logger to have time stamps
# -------------------------------------------------------------------------------------------------
const date_format = "HH:MM:SS.sss"
timestamp_logger(logger) = TransformerLogger(logger) do log
  merge(log, (; message = "$(format(now(), date_format)) $(log.message)"))
end

MOCNeutronTransport_timestamps_on = false
function log_timestamps()
    if !MOCNeutronTransport_timestamps_on
        logger = global_logger()
        logger |> timestamp_logger |> global_logger
        global MOCNeutronTransport_timestamps_on = true
    end
end

include("Tree.jl")
# include("constants.jl")
#include("operators.jl")
#include("./gmsh/gmsh_generate_rectangular_grid.jl")
#include("./gmsh/gmsh_group_preserving_fragment.jl")
#include("./gmsh/gmsh_overlay_rectangular_grid.jl")
include("./primitives/Edge.jl")
include("./primitives/Face.jl")
include("./primitives/Point.jl")
include("./primitives/LineSegment.jl")
include("./primitives/AABB.jl")
# include("./primitives/QuadraticSegment_2D.jl")
#include("./primitives/Triangle_2D.jl")
#include("./primitives/Quadrilateral_2D.jl")
#include("./primitives/Triangle6_2D.jl")
#include("./primitives/Quadrilateral8_2D.jl")
#include("L_system.jl")
#include("./mesh/UnstructuredMesh_2D.jl")
#include("./mesh/LinearUnstructuredMesh_2D.jl")
#include("./mesh/QuadraticUnstructuredMesh_2D.jl")
#include("./mesh/HierarchicalRectangularlyPartitionedMesh.jl")
#include("./mesh/IO_abaqus.jl")
##include("./mesh/IO_vtk.jl")
#include("./mesh/IO_xdmf.jl")
# include("gauss_legendre_quadrature.jl")
#include("./ray_tracing/AngularQuadrature.jl")
#include("./ray_tracing/ray_trace.jl")
#include("./ray_tracing/ray_trace_low_level.jl")
#
#
# Convenience operators
#const ∇ = gradient
#const ∇² = laplacian
#const 𝐉 = jacobian
#const × = cross
#const ⋅ = dot

# Structs/Types
export AABB, AABB_2D, AABB_3D, Edge, Face, LineSegment, LineSegment_2D, LineSegment_3D, Point,
       Point_2D, Point_3D, Tree
# export  Edge_2D,
#         Face_2D,
#         LinearUnstructuredMesh_2D,
#         LineSegment_2D,
#         Point,
#         QuadraticSegment_2D,
#         QuadraticUnstructuredMesh_2D,
#         Quadrilateral_2D,
#         Quadrilateral8_2D,
#         QuadrilateralMesh_2D,
#         Rectangle_2D,
#         Tree,
#         Triangle_2D,
#         Triangle6_2D,
#         TriangleMesh_2D,
#         UnstructuredMesh_2D
# Methods
export +, -, ⋅, ×, arclength, area, depth, distance, distance², dot, height, intersect, isapprox,
       isleft, midpoint, norm, norm², rand, sortpoints, sortpoints!, width
# export  +, -, *, /, ×, ⋅, ⪇ , ⪉ , ∇ , ∇²,
#         add_boundary_edges,
#         add_boundary_edges!,
#         add_connectivity,
#         add_connectivity!,
#         add_edges,
#         add_edges!,
#         add_everything,
#         add_everything!,
#         add_edge_face_connectivity,
#         add_face_edge_connectivity,
#         add_materialized_edges,
#         add_materialized_edges!,
#         add_materialized_faces,
#         add_materialized_faces!,
#         arclength,
#         area,
#         boundary_edges,
#         boundingbox,
#         centroid,
#         closest_point,
#         cross,
#         distance,
#         distance²,
#         dot,
#         edges,
#         edge_face_connectivity,
#         edge_points,
#         faces_sharing_vertex,
#         face_edge_connectivity,
#         face_points,
#         find_face,
#         find_segment_faces,
#         generate_angular_quadrature,
#         generate_tracks,
#         height,
#         hilbert_curve,
#         intersect,
#         intersect_edges,
#         intersect_faces,
#         isleft,
#         isstraight,
#         jacobian,
#         J,
#         log_timestamps,
#         materialize_edge,
#         materialize_edges,
#         materialize_face,
#         materialize_faces,
#         midpoint,
#         norm,
#         norm²,
#         partition_rectangularly,
#         plot_track_edge_to_edge,
#         rand,
#         ray_trace,
#         ray_trace_edge_to_edge,
#         ray_trace_angle_edge_to_edge!,
#         ray_trace_track_edge_to_edge,
#         read_abaqus_2d,
#         read_vtk_2d,
#         real_to_parametric,
#         reorder_faces_to_hilbert!,
#         reorder_points_to_hilbert!,
#         reorder_to_hilbert!,
#         segmentize,
#         shared_edge,
#         sortpoints,
#         sortpoints!,
#         sort_intersection_points!,
#         submesh,
#         to_lines,
#         triangulate,
#         union,
#         validate_ray_tracing_data,
#         width,
#         write_vtk_2d,
#         write_xdmf_2d
#
# # Gmsh
# export gmsh,
#        gmsh_generate_rectangular_grid,
#        gmsh_group_preserving_fragment,
#        gmsh_overlay_rectangular_grid
#
# # Plot
# if enable_visualization
#     export Figure, Axis
#     export scatter, linesegments, mesh,
#            scatter!, linesegments!, mesh!
# end

end # module

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
macro eachsubspace(expr)
    args = esc.(expr.args)
    sub_idx = :(eachindex($(args[2]).subspaces))
    return :( [$(args[1])($(args[2:end,]...), k) for k in $sub_idx] )
end

ext: jl


text:
@doc raw"""
    CompressibleEulerEquations3D

The compressible Euler equations for an ideal gas in three space dimensions.
"""
struct CompressibleEulerEquations3D{RealT<:Real} <: AbstractCompressibleEulerEquations{3, 5}
  gamma::RealT
end

get_name(::CompressibleEulerEquations3D) = "CompressibleEulerEquations3D"
varnames(::typeof(cons2cons), ::CompressibleEulerEquations3D) = @SVector ["rho", "rho_v1", "rho_v2", "rho_v3", "rho_e"]
varnames(::typeof(cons2prim), ::CompressibleEulerEquations3D) = @SVector ["rho", "v1", "v2", "v3", "p"]

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

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

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

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

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

  return @SVector [rho, rho_v1, rho_v2, rho_v3, rho_e]
end

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

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

  x1, x2, x3 = x
  si, co = sincos(((x1 + x2 + x3) - t) * ω)
  tmp1 = si * A
  tmp2 = co * A * ω
  tmp3 = ((((((4 * tmp1 * γ - 4 * tmp1) + 4 * c * γ) - 4c) - 3γ) + 7) * tmp2) / 2

  du1 = 2 * tmp2
  du2 = tmp3
  du3 = tmp3
  du4 = tmp3
  du5 = ((((((12 * tmp1 * γ - 4 * tmp1) + 12 * c * γ) - 4c) - 9γ) + 9) * tmp2) / 2

  # Original terms (without performance enhancements)
  # tmp2 = ((((((4 * sin(((x1 + x2 + x3) - t) * ω) * A * γ - 4 * sin(((x1 + x2 + x3) - t) * ω) * A) + 4 * c * γ) - 4c) - 3γ) + 7) * cos(((x1 + x2 + x3) - t) * ω) * A * ω) / 2
  # du1 = 2 * cos(((x1 + x2 + x3) - t) * ω) * A * ω
  # du2 = tmp2
  # du3 = tmp2
  # du4 = tmp2
  # du5 = ((((((12 * sin(((x1 + x2 + x3) - t) * ω) * A * γ - 4 * sin(((x1 + x2 + x3) - t) * ω) * A) + 12 * c * γ) - 4c) - 9γ) + 9) * cos(((x1 + x2 + x3) - t) * ω) * A * ω) / 2

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

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

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::CompressibleEulerEquations3D)
  rho = 1 + exp(-(x[1]^2 + x[2]^2 + x[3]^2))/2
  v1 = 1
  v2 = 1
  v3 = 1
  rho_v1 = rho * v1
  rho_v2 = rho * v2
  rho_v3 = rho * v3
  p = 1
  rho_e = p/(equations.gamma - 1) + 1/2 * rho * (v1^2 + v2^2 + v3^2)
  return @SVector [rho, rho_v1, rho_v2, rho_v3, rho_e]
end

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

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::CompressibleEulerEquations3D)
  # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
  # Set up spherical coordinates
  inicenter = (0, 0, 0)
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  z_norm = x[3] - inicenter[3]
  r = sqrt(x_norm^2 + y_norm^2 + z_norm^2)
  phi   = atan(y_norm, x_norm)
  theta = iszero(r) ? 0.0 : acos(z_norm / r)

  # Calculate primitive variables
  rho = r > 0.5 ? 1.0 : 1.1691
  v1  = r > 0.5 ? 0.0 : 0.1882 * cos(phi) * sin(theta)
  v2  = r > 0.5 ? 0.0 : 0.1882 * sin(phi) * sin(theta)
  v3  = r > 0.5 ? 0.0 : 0.1882 * cos(theta)
  p   = r > 0.5 ? 1.0 : 1.245

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

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

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::CompressibleEulerEquations3D)
  # Calculate radius as distance from origin
  r = sqrt(x[1]^2 + x[2]^2 + x[3]^2)

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

  # Calculate primitive variables
  rho = 1.0
  v1 = 0.0
  v2 = 0.0
  v3 = 0.0

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

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

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

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::CompressibleEulerEquations3D)
  # 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^3, 256^3
  # domain size is [-20.0,20.0]^3
  # gamma = 5/3 for this test case
  R = 1.0 # radius of the blob
  # background density
  rho = 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
  v1 = 2 * R * sqrt(Chi) / tau_cr
  v2 = 0.0
  v3 = 0.0
  Ma0 = 2.7 # background flow Mach number Ma=v/c
  c = v1 / Ma0 # sound speed
  # use perfect gas assumption to compute background pressure via the sound speed c^2 = gamma * pressure/density
  p = c * c * rho / equations.gamma
  # initial center of the blob
  inicenter = [-15, 0, 0]
  x_rel = x - inicenter
  r = sqrt(x_rel[1]^2 + x_rel[2]^2 + x_rel[3]^2)
  # steepness of the tanh transition zone
  slope = 2
  # density blob
  rho = rho + (Chi - 1) * 0.5 * (1 + (tanh(slope * (r + R)) - (tanh(slope *(r - R)) + 1)))
  # velocity blob is zero
  v1 = v1 - v1 * 0.5 * (1 + (tanh(slope *(r + R)) - (tanh(slope *(r - R)) + 1)))
  return prim2cons(SVector(rho, v1, v2, v3, p), equations)
end

"""
    initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D)

The classical inviscid Taylor-Green vortex.
"""
function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D)
  A  = 1.0 # magnitude of speed
  Ms = 0.1 # maximum Mach number

  rho = 1.0
  v1  =  A * sin(x[1]) * cos(x[2]) * cos(x[3])
  v2  = -A * cos(x[1]) * sin(x[2]) * cos(x[3])
  v3  = 0.0
  p   = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms
  p   = p + 1.0/16.0 * A^2 * rho * (cos(2*x[1])*cos(2*x[3]) + 2*cos(2*x[2]) + 2*cos(2*x[1]) + cos(2*x[2])*cos(2*x[3]))

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

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

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::CompressibleEulerEquations3D)
  # 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] + x[3] - t))
  G = 1.0 # gravitational constant

  rho = ini
  v1 = 1.0
  v2 = 1.0
  v3 = 1.0
  p = ini^2 * G * 2 / (3 * pi) # "3" is the number of spatial dimensions

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

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

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::CompressibleEulerEquations3D)
  # 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 = -4 * G / (3 * pi) # "3" is the number of spatial dimensions  # 2D: -2.0*G/pi

  x1, x2, x3 = x
  # TODO: sincospi
  si, co = sincos(pi * (x1 + x2 + x3 - t))
  rhox = A * pi * co
  rho  = c + A * si

  # In "2 * rhox", the "2" is "number of spatial dimensions minus one"
  du1 = 2 * rhox
  du2 = 2 * rhox
  du3 = 2 * rhox
  du4 = 2 * rhox
  du5 = 2 * rhox * (3/2 - C_grav*rho) # "3" in "3/2" is the number of spatial dimensions

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

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

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).

!!! note
    This method is to be used for testing pure Euler simulations with analytic self-gravity.
    If you intend to do coupled Euler-gravity simulations, you need to use
    [`source_terms_eoc_test_coupled_euler_gravity`](@ref) instead.
"""
function source_terms_eoc_test_euler(u, x, t, equations::CompressibleEulerEquations3D)
  # Same settings as in `initial_condition_eoc_test_coupled_euler_gravity`
  c = 2.0
  A = 0.1
  G = 1.0
  C_grav = -4 * G / (3 * pi) # "3" is the number of spatial dimensions

  x1, x2, x3 = x
  # TODO: sincospi
  si, co = sincos(pi * (x1 + x2 + x3 - t))
  rhox = A * pi * co
  rho  = c + A *  si

  du1 = rhox *  2
  du2 = rhox * (2 -     C_grav * rho)
  du3 = rhox * (2 -     C_grav * rho)
  du4 = rhox * (2 -     C_grav * rho)
  du5 = rhox * (3 - 5 * C_grav * rho)

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

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

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::CompressibleEulerEquations3D)
  # Calculate radius as distance from origin
  r = sqrt(x[1]^2 + x[2]^2 + x[3]^2)

  # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000
  r0 = 0.25 # = 4.0 * smallest dx (for domain length=8 and max-ref=7)
  E = 1.0
  p_inner   = (equations.gamma - 1) * E / (4/3 * pi * r0^3)
  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  = -50.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
  v3 = 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, v3, 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::CompressibleEulerEquations3D)
  # velocities are zero, density/pressure are ambient values according to
  # initial_condition_sedov_self_gravity
  rho = 1e-5
  v1 = 0.0
  v2 = 0.0
  v3 = 0.0
  p = 1e-5

  u_boundary = prim2cons(SVector(rho, v1, v2, v3, p), equations)

  # Calculate boundary flux
  if direction in (2, 4, 6) # 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

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

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

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, equations::CompressibleEulerEquations3D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  v3_ll = rho_v3_ll / rho_ll
  v1_rr = rho_v1_rr / rho_rr
  v2_rr = rho_v2_rr / rho_rr
  v3_rr = rho_v3_rr / rho_rr
  p_ll =  (equations.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2 + v2_ll^2 + v3_ll^2))
  p_rr =  (equations.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2 + v2_rr^2 + v3_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)
  v3_avg  = 1/2 * ( v3_ll +  v3_rr)
  p_avg   = 1/2 * (  p_ll +   p_rr)
  kin_avg = 1/2 * (v1_ll*v1_rr + v2_ll*v2_rr + v3_ll*v3_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 = rho_avg * v1_avg * v3_avg
    f5 = p_avg*v1_avg/(equations.gamma-1) + rho_avg*v1_avg*kin_avg + pv1_avg
  elseif orientation == 2
    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 = rho_avg * v2_avg * v3_avg
    f5 = p_avg*v2_avg/(equations.gamma-1) + rho_avg*v2_avg*kin_avg + pv2_avg
  else
    pv3_avg = 1/2 * (p_ll*v3_rr + p_rr*v3_ll)
    f1 = rho_avg * v3_avg
    f2 = rho_avg * v3_avg * v1_avg
    f3 = rho_avg * v3_avg * v2_avg
    f4 = rho_avg * v3_avg * v3_avg + p_avg
    f5 = p_avg*v3_avg/(equations.gamma-1) + rho_avg*v3_avg*kin_avg + pv3_avg
  end

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

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

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, equations::CompressibleEulerEquations3D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll/rho_ll
  v2_ll = rho_v2_ll/rho_ll
  v3_ll = rho_v3_ll/rho_ll
  v1_rr = rho_v1_rr/rho_rr
  v2_rr = rho_v2_rr/rho_rr
  v3_rr = rho_v3_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)
  v3_avg = 1/2 * (v3_ll + v3_rr)
  p_avg = 1/2 * ((equations.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2 + v2_ll^2 + v3_ll^2)) +
                 (equations.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2 + v2_rr^2 + v3_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 * v1_avg * v3_avg
    f5 = (rho_avg * e_avg + p_avg) * v1_avg
  elseif orientation == 2
    f1 = rho_avg * v2_avg
    f2 = rho_avg * v2_avg * v1_avg
    f3 = rho_avg * v2_avg * v2_avg + p_avg
    f4 = rho_avg * v2_avg * v3_avg
    f5 = (rho_avg * e_avg + p_avg) * v2_avg
  else
    f1 = rho_avg * v3_avg
    f2 = rho_avg * v3_avg * v1_avg
    f3 = rho_avg * v3_avg * v2_avg
    f4 = rho_avg * v3_avg * v3_avg + p_avg
    f5 = (rho_avg * e_avg + p_avg) * v3_avg
  end

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

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

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, equations::CompressibleEulerEquations3D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll/rho_ll
  v2_ll = rho_v2_ll/rho_ll
  v3_ll = rho_v3_ll/rho_ll
  v1_rr = rho_v1_rr/rho_rr
  v2_rr = rho_v2_rr/rho_rr
  v3_rr = rho_v3_rr/rho_rr
  p_ll =  (equations.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2 + v2_ll^2 + v3_ll^2))
  p_rr =  (equations.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2 + v2_rr^2 + v3_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 + v3_ll^2)
  specific_kin_rr = 0.5*(v1_rr^2 + v2_rr^2 + v3_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)
  v3_avg = 0.5*(v3_ll+v3_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 * v3_avg
    f5 = f1 * 0.5*(1/(equations.gamma-1)/beta_mean - velocity_square_avg)+ f2*v1_avg + f3*v2_avg + f4*v3_avg
  elseif orientation == 2
    f1 = rho_mean * v2_avg
    f2 = f1 * v1_avg
    f3 = f1 * v2_avg + p_mean
    f4 = f1 * v3_avg
    f5 = f1 * 0.5*(1/(equations.gamma-1)/beta_mean - velocity_square_avg)+ f2*v1_avg + f3*v2_avg + f4*v3_avg
  else
    f1 = rho_mean * v3_avg
    f2 = f1 * v1_avg
    f3 = f1 * v2_avg
    f4 = f1 * v3_avg + p_mean
    f5 = f1 * 0.5*(1/(equations.gamma-1)/beta_mean - velocity_square_avg)+ f2*v1_avg + f3*v2_avg + f4*v3_avg
  end

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

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

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, equations::CompressibleEulerEquations3D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  v3_ll = rho_v3_ll / rho_ll
  v1_rr = rho_v1_rr / rho_rr
  v2_rr = rho_v2_rr / rho_rr
  v3_rr = rho_v3_rr / rho_rr
  p_ll =  (equations.gamma - 1) * (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2 + v3_ll^2))
  p_rr =  (equations.gamma - 1) * (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2 + v3_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)
  v3_avg = 0.5 * (v3_ll + v3_rr)
  p_avg  = 0.5 * (p_ll + p_rr)
  velocity_square_avg = 0.5 * (v1_ll*v1_rr + v2_ll*v2_rr + v3_ll*v3_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 * v3_avg
    f5 = f1 * ( velocity_square_avg + 1 / ((equations.gamma-1) * rho_p_mean) ) + 0.5 * (p_ll*v1_rr + p_rr*v1_ll)
  elseif orientation == 2
    f1 = rho_mean * v2_avg
    f2 = f1 * v1_avg
    f3 = f1 * v2_avg + p_avg
    f4 = f1 * v3_avg
    f5 = f1 * ( velocity_square_avg + 1 / ((equations.gamma-1) * rho_p_mean) ) + 0.5 * (p_ll*v2_rr + p_rr*v2_ll)
  else # orientation == 3
    f1 = rho_mean * v3_avg
    f2 = f1 * v1_avg
    f3 = f1 * v2_avg
    f4 = f1 * v3_avg + p_avg
    f5 = f1 * ( velocity_square_avg + 1 / ((equations.gamma-1) * rho_p_mean) ) + 0.5 * (p_ll*v3_rr + p_rr*v3_ll)
  end

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

function flux_lax_friedrichs(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D)
  # Calculate primitive variables and speed of sound
  rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  v3_ll = rho_v3_ll / rho_ll
  v_mag_ll = sqrt(v1_ll^2 + v2_ll^2 + v3_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
  v3_rr = rho_v3_rr / rho_rr
  v_mag_rr = sqrt(v1_rr^2 + v2_rr^2 + v3_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)

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

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

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

function flux_hll(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D)
  # Calculate primitive variables and speed of sound
  rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = u_rr

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

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

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

  if orientation == 1 # x-direction
    Ssl = v1_ll - sqrt(equations.gamma * p_ll / rho_ll)
    Ssr = v1_rr + sqrt(equations.gamma * p_rr / rho_rr)
  elseif orientation == 2 # y-direction
    Ssl = v2_ll - sqrt(equations.gamma * p_ll / rho_ll)
    Ssr = v2_rr + sqrt(equations.gamma * p_rr / rho_rr)
  else # z-direction
    Ssl = v3_ll - sqrt(equations.gamma * p_ll / rho_ll)
    Ssr = v3_rr + sqrt(equations.gamma * p_rr / rho_rr)
  end

  if Ssl >= 0.0 && Ssr > 0.0
    f1 = f_ll[1]
    f2 = f_ll[2]
    f3 = f_ll[3]
    f4 = f_ll[4]
    f5 = f_ll[5]
  elseif Ssr <= 0.0 && Ssl < 0.0
    f1 = f_rr[1]
    f2 = f_rr[2]
    f3 = f_rr[3]
    f4 = f_rr[4]
    f5 = f_rr[5]
  else
    f1 = (Ssr*f_ll[1] - Ssl*f_rr[1] + Ssl*Ssr*(rho_rr    - rho_ll))    / (Ssr - Ssl)
    f2 = (Ssr*f_ll[2] - Ssl*f_rr[2] + Ssl*Ssr*(rho_v1_rr - rho_v1_ll)) / (Ssr - Ssl)
    f3 = (Ssr*f_ll[3] - Ssl*f_rr[3] + Ssl*Ssr*(rho_v2_rr - rho_v2_ll)) / (Ssr - Ssl)
    f4 = (Ssr*f_ll[4] - Ssl*f_rr[4] + Ssl*Ssr*(rho_v3_rr - rho_v3_ll)) / (Ssr - Ssl)
    f5 = (Ssr*f_ll[5] - Ssl*f_rr[5] + Ssl*Ssr*(rho_e_rr  - rho_e_ll))  / (Ssr - Ssl)
  end

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

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

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, equations::CompressibleEulerEquations3D)
  # Calculate primitive variables and speed of sound
  rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  v3_ll = rho_v3_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 + v3_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
  v3_rr = rho_v3_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 + v3_rr^2))
  c_rr = sqrt(equations.gamma*p_rr/rho_rr)

  # Obtain left and right fluxes
  f_ll = calcflux(u_ll, orientation, equations)
  f_rr = calcflux(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 + (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_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 + (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr)^2
  else # z-direction
    vel_L = v3_ll
    vel_R = v3_rr
    ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2 + (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_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]
    f5 = f_ll[5]
  elseif Ssr <= 0.0
    f1 = f_rr[1]
    f2 = f_rr[2]
    f3 = f_rr[3]
    f4 = f_rr[4]
    f5 = f_rr[5]
  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
      UStar5 = densStar*enerStar
      if orientation == 1 # x-direction
        UStar2 = densStar*SStar
        UStar3 = densStar*v2_ll
        UStar4 = densStar*v3_ll
      elseif orientation == 2 # y-direction
        UStar2 = densStar*v1_ll
        UStar3 = densStar*SStar
        UStar4 = densStar*v3_ll
      else # z-direction
        UStar2 = densStar*v1_ll
        UStar3 = densStar*v2_ll
        UStar4 = 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_v3_ll)
      f5 = f_ll[5]+Ssl*(UStar5 - 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
      UStar5 = densStar*enerStar
      if orientation == 1 # x-direction
        UStar2 = densStar*SStar
        UStar3 = densStar*v2_rr
        UStar4 = densStar*v3_rr
      elseif orientation == 2 # y-direction
        UStar2 = densStar*v1_rr
        UStar3 = densStar*SStar
        UStar4 = densStar*v3_rr
      else # z-direction
        UStar2 = densStar*v1_rr
        UStar3 = densStar*v2_rr
        UStar4 = 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_v3_rr)
      f5 = f_rr[5]+Ssr*(UStar5 - rho_e_rr)
    end
  end
  return SVector(f1, f2, f3, f4, f5)
end

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

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

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

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

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

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

  v1 = rho_v1 / rho
  v2 = rho_v2 / rho
  v3 = rho_v3 / rho
  v_square = v1^2 + v2^2 + v3^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 * v3
  w5 = -rho_p

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

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

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

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

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

# Calculate thermodynamic entropy for a conservative state `cons`
@inline function entropy_thermodynamic(cons, equations::CompressibleEulerEquations3D)
  # Pressure
  p = (equations.gamma - 1) * (cons[5] - 1/2 * (cons[2]^2 + cons[3]^2 + cons[4]^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::CompressibleEulerEquations3D)
  S = -entropy_thermodynamic(cons, equations) * cons[1] / (equations.gamma - 1)
  # Mathematical entropy

  return S
end

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

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

# Calculate kinetic energy for a conservative state `cons`
@inline function energy_kinetic(cons, equations::CompressibleEulerEquations3D)
  return 0.5 * (cons[2]^2 + cons[3]^2 + cons[4]^2)/cons[1]
end

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


text:
using BayesianDataFusion
using Base.Test

using DataFrames

Y  = sprand(15,10, 0.1)
rd = RelationData(Y, class_cut = 0.5, alpha_sample=true)
assignToTest!(rd.relations[1], 2)

result = macau(rd, burnin = 5, psamples = 6, verbose = false)
@test rd.relations[1].model.alpha > 0


text:
if VERSION < v"0.7.0"
    Pkg.clone(pwd())
else
    import Pkg
    Pkg.develop(Pkg.PackageSpec(url=pwd()))
end

Pkg.test("Coverage"; coverage=true)


text:
# Import libraries.
using Soss, Flux, Plots, Random

# Number of points to generate.
N = 80
M = round(Int, N / 4)
Random.seed!(1234)

# Generate artificial data.
x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
xt1s = Array([[x1s[i] + 0.5; x2s[i] + 0.5] for i = 1:M])
x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
append!(xt1s, Array([[x1s[i] - 5; x2s[i] - 5] for i = 1:M]))

x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
xt0s = Array([[x1s[i] + 0.5; x2s[i] - 5] for i = 1:M])
x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
append!(xt0s, Array([[x1s[i] - 5; x2s[i] + 0.5] for i = 1:M]))

# Store all the data for later.
xs = [xt1s; xt0s]
ts = [ones(2*M); zeros(2*M)]

# Plot data points.
function plot_data()
    x1 = map(e -> e[1], xt1s)
    y1 = map(e -> e[2], xt1s)
    x2 = map(e -> e[1], xt0s)
    y2 = map(e -> e[2], xt0s)

    Plots.scatter(x1,y1, color="red", clim = (0,1))
    Plots.scatter!(x2, y2, color="blue", clim = (0,1))
end

plot_data()

# Turn a vector into a set of weights and biases.
function unpack(nn_params::AbstractVector)
    W₁ = reshape(nn_params[1:6], 3, 2);
    b₁ = reshape(nn_params[7:9], 3)

    W₂ = reshape(nn_params[10:15], 2, 3);
    b₂ = reshape(nn_params[16:17], 2)

    Wₒ = reshape(nn_params[18:19], 1, 2);
    bₒ = reshape(nn_params[20:20], 1)
    return W₁, b₁, W₂, b₂, Wₒ, bₒ
end

# Construct a neural network using Flux and return a predicted value.
function nn_forward(xs, nn_params::AbstractVector)
    W₁, b₁, W₂, b₂, Wₒ, bₒ = unpack(nn_params)
    nn = Chain(Dense(W₁, b₁, tanh),
               Dense(W₂, b₂, tanh),
               Dense(Wₒ, bₒ, σ))
    return nn.(xs)
end;

# Create a regularization term and a Gaussain prior variance term.
alpha = 0.09
sig = sqrt(1.0 / alpha)

# Specify the probabalistic model.
bayes_nn = @model xs begin
    # Create the weight and bias vector.
    nn_params ~ Normal(0,3.0) |> iid(20)

    # Calculate predictions for the inputs given the weights
    # and biases in theta.
    preds = nn_forward(xs, nn_params)

    # Observe each prediction.
    ts ~ For(preds) do p
        Bernoulli(p[1])
    end

end;

post = dynamicHMC(bayes_nn(xs=xs), (ts=ts,); ad_backend=Val(:Zygote))

# Perform inference.
N = 5000
ch = simulate(bayes_nn(hcat(xs...), ts), HMC(0.05, 4), N);

# Extract all weight and bias parameters.
theta = ch[:nn_params].value.data;

# Plot the data we have.
plot_data()

# Find the index that provided the highest log posterior in the chain.
_, i = findmax(ch[:lp].value.data)

# Extract the max row value from i.
i = i.I[1]

# Plot the posterior distribution with a contour plot.
x_range = collect(range(-6,stop=6,length=25))
y_range = collect(range(-6,stop=6,length=25))
Z = [nn_forward([x, y], theta[i, :])[1] for x=x_range, y=y_range]
contour!(x_range, y_range, Z)

# Return the average predicted value across
# multiple weights.
function nn_predict(x, theta, num)
    mean([nn_forward(x, theta[i,:])[1] for i in 1:10:num])
end;

# Plot the average prediction.
plot_data()

n_end = 1500
x_range = collect(range(-6,stop=6,length=25))
y_range = collect(range(-6,stop=6,length=25))
Z = [nn_predict([x, y], theta, n_end)[1] for x=x_range, y=y_range]
contour!(x_range, y_range, Z)

# Number of iterations to plot.
n_end = 500

anim = @gif for i=1:n_end
    plot_data()
    Z = [nn_forward([x, y], theta[i,:])[1] for x=x_range, y=y_range]
    contour!(x_range, y_range, Z, title="Iteration $$i", clim = (0,1))
end every 5


text:
module DE

using Random

export DE_evolution

"""
population type defination for Differential Evolutionary Algorithm
"""
mutable struct Population
    x::Array{Float64, 2}
    fit::Array{Float64, 1}
end

"""
    sample(pop::Population)

sample a new population with mutation and recombination in DE
"""
function sample(pop::Population)
    new_pop = Population(zeros(size(pop.x)), zeros(size(pop.fit)))
    cr, f, num, len = rand(), rand(), size(pop.x, 1), size(pop.x, 2)
    for i in 1:num
        ind = shuffle(Vector(1:num))[1:3]
        j = rand(1:len)
        for k in 1:len
            if rand() < cr || k == j
                new_pop.x[i, k] = pop.x[ind[1], k] +
                    f * (pop.x[ind[2], k] - pop.x[ind[3], k])
            else
                new_pop.x[i, k] = pop.x[i, k]
            end
        end
    end
    return new_pop
end

"""
    selection(pop::Population, new_pop::Population)

choose a better individual between 'parent' and 'child'
"""
function selection(pop::Population, new_pop::Population)
    for i in 1:size(pop.x, 1)
        if new_pop.fit[i] <= pop.fit[i]
            pop.x[i, :] = new_pop.x[i, :]
            pop.fit[i] = new_pop.fit[i]
        end
    end
end

"""
    pop_opti(pop::Population)

return the best individual from a population
"""
function pop_opti(pop::Population)
    index = argmin(pop.fit)
    return pop.x[index, :], pop.fit[index]
end

"""
    evolution(eval, pops::Int64, dim::Int64, gen::Int64)

Optimization for function 'eval' uses algorithm 'ea'.

# Argument

- 'pops': population size
- 'dim': dimension of the problem
- 'gen': max generation number
"""
function DE_evolution(f, pops::Int64, dim::Int64, gen::Int64)
    pop = Population(rand(pops, dim), zeros(pops))
    f(pop)
    println("start optimization:")
    println("gen: ", 0, " fitness: ", pop_opti(pop)[2])
    for i in 1:gen
        new_pop = sample(pop)
        f(new_pop)
        selection(pop, new_pop)
        println("gen: ", i, " fitness: ", pop_opti(pop)[2])
    end
    return pop_opti(pop)[1], pop_opti(pop)[2]
end

end  # module DE


text:
using Test: @test, @testset
using LogicCircuits
using ProbabilisticCircuits

@testset "Probabilistic circuits learner tests" begin
    data = dataset(twenty_datasets("nltcs"); do_shuffle=false, batch_size=-1)
    train_x = train(data)
    pc = learn_probabilistic_circuit(train_x; pseudocount = 1.0, algo = "chow-liu", algo_kwargs=(α=1.0, clt_root="graph_center"))

    # simple test
    @test pc isa ProbΔ
    @test check_parameter_integrity(pc)
    @test num_parameters(pc) == 62
    @test pc[26].log_thetas[1] ≈ -0.023528423773273476 atol=1.0e-7

    # all evidence sums to 1
    N = num_features(train_x);
    data_all = XData(generate_data_all(N))
    fc = FlowΔ(pc, max_batch_size(train_x), Bool, opts_accumulate_flows)
    calc_prob_all = log_likelihood_per_instance(fc, data_all)
    calc_prob_all = exp.(calc_prob_all)
    sum_prob_all = sum(calc_prob_all)
    @test sum_prob_all ≈ 1 atol = 1.0e-7;
end


text:
"""
Based on an opaque `RoutingData` object `rd` and a solver (for instance, `CPLEX.Optimizer`), computes a routing.
It returns one object, a `RoutingSolution`, containing all intermediary results.
"""
function compute_routing(rd::RoutingData)
    # Only perform dispatch on the type of model.
    # TODO: everywhere if iterative, add parameter to remember all generated routings or not (store it in RoutingData?).
    mt = rd.model_type
    return compute_routing(rd, mt.edge_obj, mt.agg_obj, mt.type, Val(mt.cg), mt.algo, mt.unc, mt.uncparams)
end

function compute_routing(rd::RoutingData, edge_obj::EdgeWiseObjectiveFunction, agg_obj::AggregationObjectiveFunction,
                         type::FormulationType, cg::Val, algo::AlgorithmChoice, unc::UncertaintyHandling, uncparams::UncertainParameters)
    # When dispatch fails, use this fallback.
    msg = "Model type not yet implemented: ModelType($edge_obj, $agg_obj, $type, $cg, $algo, $unc, $uncparams).\n"
    msg *= "These models are available:\n"
    msg *= string(methods(compute_routing))
    error(msg)
end


text:
using HistogramThresholding, Images, TestImages
using Test

function gaussian_mixture(z;P₁,μ₁,σ₁,P₂,μ₂,σ₂)
    (P₁ / (σ₁*√(2*π))) * exp( -(z-μ₁)^2 / (2*σ₁^2) )  +  (P₂ / (σ₂*√(2*π))) * exp( -(z-μ₂)^2 / (2*σ₂^2) )
end

function surge(t;A,k)
    A*t*exp(-k*t)
end

@testset "HistogramThresholding.jl" begin
    include("otsu.jl")
    include("minimum.jl")
    include("intermodes.jl")
    include("unimodal.jl")
    include("moments.jl")
    include("minimum_error.jl")
end


text:
using Test
using ParserCombinator
using GraphIO.GML

@testset "GML" begin
    for g in values(allgraphs)
        readback_test(GMLFormat(), g, testfail=true)
    end
    write_test(GMLFormat(), allgraphs)
end


text:
###########################################################
# LOAD MODULES
###########################################################
using Combinatorics
using Random
using CSV
using DataFrames
###########################################################
# DEFINE CONSTANTS
###########################################################
const combns2 = combinations(1:4,2)|>collect
const combns3 = combinations(1:4,3)|>collect
const full_hand = combinations(1:4,4)|>collect
const suit_names = ["♣", "♦", "♥", "♠"]
const rank_names = ["A", "2", "3", "4", "5", "6", "7", "8", "9", "10", "J", "Q", "K"]
###########################################################
# SET UP STRUCTS
###########################################################
# define card struct
struct card
    suit :: Int64
    rank :: Int64
    value :: Int64
    function card(suit::Int64, rank::Int64, value::Int64)
        @assert(1 ≤ suit ≤ 4, "suit is not between 1 and 4")
        @assert(1 ≤ rank ≤ 13, "rank is not between 1 and 13")
        @assert(1 ≤ value ≤ 10, "value is not between 1 and 10")
        new(suit, rank, value)
    end
end
# define print method
function Base.show(io::IO, card::card)
    print(io, rank_names[card.rank], suit_names[card.suit])
end
#define deck struct
abstract type card_set end

struct deck <: card_set
    cards :: Array{card, 1}
end
function deck()
        Deck = deck(card[])
        for suit in 1:4
            for rank in 1:13
                if rank ≤ 10
                push!(Deck.cards, card(suit, rank, rank))
            else
                push!(Deck.cards, card(suit, rank, 10))
            end
            end
        end
        Deck
    end
end
function Base.show(io::IO, deck::deck)
    for card in deck.cards
        print(io, card, " ")
    end
    println()
end
struct hand <: card_set
    cards :: Array{card, 1}
    label :: String
end
function hand(label::String="")
    hand(card[], label)
end
function Base.show(io::IO, cs::card_set)
    for card in cs.cards
        print(io, card, " ")
    end
end
function Base.pop!(cs::card_set)
    pop!(cs.cards)
end
function Base.deleteat!(cs::card_set, i)
    deleteat!(cs.cards, i::Int64)
end
function Base.push!(cs::card_set, card::card)
    push!(cs.cards, card)
    nothing
end
function move!(cs1::card_set, cs2::card_set, n::Int)
    @assert 1 ≤ n ≤ length(cs1.cards)
    for i in 1:n
        card = pop!(cs1)
        push!(cs2, card)
    end
    nothing
end
function Random.shuffle!(deck::deck)
    shuffle!(deck.cards)
    deck
end
function Base.isless(c1::card, c2::card)
    c1.rank < c2.rank
end
###########################################################
# DEFINE FUNCTIONS
###########################################################
# score hand needs some work--a lot of the garbage collection overhead comes from this function
# what are the alternatives to creating storage vectors in every iteration?
# using temporary vectors is admittedly inelegant and the major speed bottleneck in the simulations (~50% of the random pegging time is GC overhead)
# they may be fast enough to answer the question but 30-50% gc time suggests I need an alternative...
function score_hand(hand_combn)
    hand_scores = Int[]
    for sel_hand in hand_combn
        fifteens = Int[]
        pairs = Int[]
        runs = Int[]
        flush = Int[]
        for combn in combns2
            # check 2 card combinations for 15s
            if sel_hand[combn[1]].value + sel_hand[combn[2]].value == 15
                push!(fifteens, 2::Int)
            end
            # check 2 card combinations for pairs
            if sel_hand[combn[1]].rank == sel_hand[combn[2]].rank
                push!(pairs, 2::Int)
            end
        end
        for combn in combns3
            # put combination in order
            sorted_combn = sort([sel_hand[combn[1]].rank, sel_hand[combn[2]].rank, sel_hand[combn[3]].rank])
            # run of 3: if the absolute difference between cards 1 and 2 is one and the abs diff of cards 2 and 3 is one
            if abs(diff([sorted_combn[1], sorted_combn[2]])[1]) == 1 && abs(diff([sorted_combn[2], sorted_combn[3]])[1]) == 1
                push!(runs, 3::Int)
            end
            # check for 3 card 15s
            if sel_hand[combn[1]].value + sel_hand[combn[2]].value + sel_hand[combn[3]].value == 15
                push!(fifteens, 2::Int)
            end
        end
        for combn in full_hand
            # put combination in order
            sorted_combn = sort([sel_hand[combn[1]].rank, sel_hand[combn[2]].rank, sel_hand[combn[3]].rank, sel_hand[combn[4]].rank])
            # run of 4: same logic as run of 3
            if abs(diff([sorted_combn[1], sorted_combn[2]])[1]) == 1 && abs(diff([sorted_combn[2],sorted_combn[3]])[1]) == 1 && abs(diff([sorted_combn[3],sorted_combn[4]])[1]) == 1
                push!(runs, 4::Int)
            end
            # check for flush
            if sel_hand[combn[1]].suit == sel_hand[combn[2]].suit && sel_hand[combn[2]].suit == sel_hand[combn[3]].suit && sel_hand[combn[3]].suit == sel_hand[combn[4]].suit
                push!(flush, 4::Int)
            end
            # check for 4 card 15s
            if sel_hand[combn[1]].value + sel_hand[combn[2]].value + sel_hand[combn[3]].value + sel_hand[combn[4]].value == 15
                push!(fifteens, 2::Int)
            end
        end
        push!(hand_scores, sum(vcat(fifteens, pairs, runs, flush)))
    end
    # grab the first hand in the list with the maximum number of points
    # future versions could include a "smarter" selection process, e.g.
    # given the hand Q 5 8 3 9 2, the algorithm currently picks Q 5 8 2 (4 pts) vs Q 5 2 3 (4 pts, ~10% potential for run)
    # Desirable to implement a modifier that gives additional weight to logical combinations of cards
    # Not worthwhile to implement probability maximization--not realistic to most cribbage games
    scored_hand = hand(hand_combn[findfirst(x->x == maximum(hand_scores), hand_scores)], "Scored Hand")
    return scored_hand
end

function peg(scored_hand::hand, peg_pile::hand, pile_score::Vector{Int64}, s_card, s_reason, h_idx, passed_p1::Vector{Bool}, passed_p2::Vector{Bool}, turn::Vector{Int}, idx::Int)

    player = sum(turn) % 2 == 0 ? 2::Int : 1::Int
    # get all legal plays
    pos_plays = Bool[]
    for c in 1:length(scored_hand.cards)
        push!(pos_plays, scored_hand.cards[c].value + sum(pile_score) <= 31)
    end
    # pick pegging choice that yields max points
    if any(pos_plays)
        card_scores = Int[]
        score_reasons = String[]
        # check--are top two cards of pegging stack a pair?
        if length(peg_pile.cards) >= 2
            top_pair = peg_pile.cards[length(peg_pile.cards)].rank == peg_pile.cards[length(peg_pile.cards)-1].rank
        else
            top_pair = false
        end
        # check--if top cards are a pair, is it a triple?
        if top_pair && length(peg_pile.cards) >= 3
            top_trip = peg_pile.cards[length(peg_pile.cards)-1].rank == peg_pile.cards[length(peg_pile.cards)-2].rank
        else
            top_trip = false
        end
        # iterate over possible plays
        for play in 1:length(pos_plays)
            if pos_plays[play] == true
                fifteens = Int[]
                pairs = Int[]
                runs = Int[]
                thirty_one = Int[]
                reason = String[]

                if scored_hand.cards[play].value + sum(pile_score) == 15
                    push!(fifteens, 2::Int)
                    push!(reason, "f_")
                end

                if scored_hand.cards[play].value + sum(pile_score) == 31
                    push!(thirty_one, 2::Int)
                    push!(reason, "t_")
                end

                if length(peg_pile.cards) >= 1 && !top_pair && !top_trip && scored_hand.cards[play].rank == peg_pile.cards[length(peg_pile.cards)].rank
                    push!(pairs, 2::Int)
                    push!(reason, "p_")
                end

                if length(peg_pile.cards) >= 2 && top_pair && !top_trip && scored_hand.cards[play].rank == peg_pile.cards[length(peg_pile.cards)].rank
                    push!(pairs, 6::Int)
                    push!(reason, "p6_")
                end

                if top_trip && scored_hand.cards[play].rank == peg_pile.cards[length(peg_pile.cards)].rank
                    push!(pairs, 12::Int)
                    push!(reason, "p12_")
                end
                # checks for runs up to length 4--it would be possible to go deeper, but for the time being I arrived at 4 as a realistic "max" run depth
                # after validation, doable to add checks for greater depth
                if !top_pair && length(peg_pile.cards) >= 3
                    sorted_stack = sort([scored_hand.cards[play].rank, peg_pile.cards[length(peg_pile.cards)].rank, peg_pile.cards[length(peg_pile.cards)-1].rank, peg_pile.cards[length(peg_pile.cards)-2].rank])
                    if abs(diff([sorted_stack[1], sorted_stack[2]])[1]) == 1 && abs(diff([sorted_stack[2],sorted_stack[3]])[1]) == 1 && abs(diff([sorted_stack[3],sorted_stack[4]])[1]) == 1
                        push!(runs, 4::Int)
                        push!(reason, "r4_")
                        is_r4 = true
                    else
                        is_r4 = false
                    end
                else
                    is_r4 = false
                end

                if !top_pair && !is_r4 && length(peg_pile.cards) >= 2
                    sorted_stack = sort([scored_hand.cards[play].rank, peg_pile.cards[length(peg_pile.cards)].rank, peg_pile.cards[length(peg_pile.cards)-1].rank])
                    if abs(diff([sorted_stack[1], sorted_stack[2]])[1]) == 1 && abs(diff([sorted_stack[2],sorted_stack[3]])[1]) == 1
                        push!(runs, 3::Int)
                        push!(reason, "r3_")
                    end
                end

                push!(card_scores, sum(vcat(fifteens, pairs, runs, thirty_one)))

                full_reason = ""
                # concatenate reasons into single string
                for r in reason
                    full_reason = r*full_reason
                end
                push!(score_reasons, full_reason)
            else
                push!(card_scores, -1)
                push!(score_reasons, "")
            end
        end
        push!(peg_pile, scored_hand.cards[findfirst(x->x == maximum(card_scores), card_scores)])
        push!(pile_score, scored_hand.cards[findfirst(x->x == maximum(card_scores), card_scores)].value)
        if maximum(card_scores) > 0
            push!(s_card, scored_hand.cards[findfirst(x->x == maximum(card_scores), card_scores)].rank)
            push!(s_reason, score_reasons[findfirst(x->x == maximum(card_scores), card_scores)])
            push!(h_idx, idx)
        end
        deleteat!(scored_hand, findfirst(x->x == maximum(card_scores), card_scores))
        if sum(pile_score) == 31
            for i in 1:length(pile_score)
                deleteat!(pile_score, 1)
            end
            for i in 1:length(peg_pile.cards)
                deleteat!(peg_pile.cards, 1)
            end

            deleteat!(passed_p1, 1)
            push!(passed_p1, false)
            deleteat!(passed_p2, 1)
            push!(passed_p2, false)
        end
        push!(turn, 1)
    elseif !any(pos_plays)
        if player == 1 && passed_p2[1] == true
            # reset pile score, pegging pile, passed tracker, advance turn
            for i in 1:length(pile_score)
                deleteat!(pile_score, 1)
            end
            for i in 1:length(peg_pile.cards)
                deleteat!(peg_pile.cards, 1)
            end
            deleteat!(passed_p2, 1)
            push!(passed_p2, false)
            push!(turn, 1)

        elseif player == 2 && passed_p1[1] == true
            # reset pile score, pegging pile, passed tracker, advance turn
            for i in 1:length(pile_score)
                deleteat!(pile_score, 1)
            end
            for i in 1:length(peg_pile.cards)
                deleteat!(peg_pile.cards, 1)
            end
            deleteat!(passed_p1, 1)

            push!(passed_p1, false)
            push!(turn, 1)

        elseif player == 1 && passed_p2[1] == false
            deleteat!(passed_p1, 1)
            push!(passed_p1, true)
            push!(turn, 1)
        else player == 2 && passed_p1[1] == false
            deleteat!(passed_p2, 1)
            push!(passed_p2, true)
            push!(turn, 1)
        end
    end
end

function peg_random(scored_hand::hand, peg_pile::hand, pile_score::Vector{Int64}, s_card, s_reason, h_idx, passed_p1::Vector{Bool}, passed_p2::Vector{Bool}, turn::Vector{Int}, idx::Int)
    player = sum(turn) % 2 == 0 ? 2::Int : 1::Int
    # get all legal plays
    pos_plays = Bool[]
    for c in 1:length(scored_hand.cards)
        push!(pos_plays, scored_hand.cards[c].value + sum(pile_score) <= 31)
    end
    if any(pos_plays)
        play = findfirst(x->x == true, pos_plays)

        reason = String[]
            # check--are top two cards of pegging stack a pair?
        if length(peg_pile.cards) >= 2
            top_pair = peg_pile.cards[length(peg_pile.cards)].rank == peg_pile.cards[length(peg_pile.cards)-1].rank
        else
            top_pair = false
        end
            # check--if top cards are a pair, is it a triple?
        if top_pair && length(peg_pile.cards) >= 3
            top_trip = peg_pile.cards[length(peg_pile.cards)-1].rank == peg_pile.cards[length(peg_pile.cards)-2].rank
        else
            top_trip = false
        end

        reason = String[]

        if scored_hand.cards[play].value + sum(pile_score) == 15
            push!(reason, "f_")
        end

        if scored_hand.cards[play].value + sum(pile_score) == 31
            push!(reason, "t_")
        end

        if length(peg_pile.cards) >= 1 && !top_pair && !top_trip && scored_hand.cards[play].rank == peg_pile.cards[length(peg_pile.cards)].rank
            push!(reason, "p_")
        end

        if length(peg_pile.cards) >= 2 && top_pair && !top_trip && scored_hand.cards[play].rank == peg_pile.cards[length(peg_pile.cards)].rank
            push!(reason, "p6_")
        end

        if top_trip && scored_hand.cards[play].rank == peg_pile.cards[length(peg_pile.cards)].rank
            push!(reason, "p12_")
        end
        # checks for runs up to length 4--it would be possible to go deeper, but for the time being I arrived at 4 as a realistic "max" run depth
        # after validation, doable to add checks for greater depth
        if !top_pair && length(peg_pile.cards) >= 3
            sorted_stack = sort([scored_hand.cards[play].rank, peg_pile.cards[length(peg_pile.cards)].rank, peg_pile.cards[length(peg_pile.cards)-1].rank, peg_pile.cards[length(peg_pile.cards)-2].rank])
            if abs(diff([sorted_stack[1], sorted_stack[2]])[1]) == 1 && abs(diff([sorted_stack[2],sorted_stack[3]])[1]) == 1 && abs(diff([sorted_stack[3],sorted_stack[4]])[1]) == 1
                push!(reason, "r4_")
                is_r4 = true
            else
                is_r4 = false
            end
        else
            is_r4 = false
        end

        if !top_pair && !is_r4 && length(peg_pile.cards) >= 2
            sorted_stack = sort([scored_hand.cards[play].rank, peg_pile.cards[length(peg_pile.cards)].rank, peg_pile.cards[length(peg_pile.cards)-1].rank])
            if abs(diff([sorted_stack[1], sorted_stack[2]])[1]) == 1 && abs(diff([sorted_stack[2],sorted_stack[3]])[1]) == 1
                push!(reason, "r3_")
            end
        end

        if length(reason) > 0
            full_reason = ""
            # concatenate reasons into single string
            for r in reason
                full_reason = r*full_reason
            end
            push!(s_reason, full_reason)
            push!(s_card, scored_hand.cards[play].rank)
            push!(h_idx, idx)
        end

        push!(peg_pile, scored_hand.cards[play])
        push!(pile_score, scored_hand.cards[play].value)
        deleteat!(scored_hand, play)

        if sum(pile_score) == 31
            for i in 1:length(pile_score)
                deleteat!(pile_score, 1)
            end
            for i in 1:length(peg_pile.cards)
                deleteat!(peg_pile.cards, 1)
            end
            deleteat!(passed_p1, 1)
            push!(passed_p1, false)
            deleteat!(passed_p2, 1)
            push!(passed_p2, false)
        end
        push!(turn, 1)
    elseif !any(pos_plays)
        if player == 1 && passed_p2[1] == true
            # reset pile score, pegging pile, passed tracker, advance turn
            for i in 1:length(pile_score)
                deleteat!(pile_score, 1)
            end
            for i in 1:length(peg_pile.cards)
                deleteat!(peg_pile.cards, 1)
            end
            deleteat!(passed_p2, 1)
            push!(passed_p2, false)
            push!(turn, 1)

        elseif player == 2 && passed_p1[1] == true
            # reset pile score, pegging pile, passed tracker, advance turn
            for i in 1:length(pile_score)
                deleteat!(pile_score, 1)
            end
            for i in 1:length(peg_pile.cards)
                deleteat!(peg_pile.cards, 1)
            end
            deleteat!(passed_p1, 1)
            push!(passed_p1, false)
            push!(turn, 1)

        elseif player == 1 && passed_p2[1] == false
            deleteat!(passed_p1, 1)
            push!(passed_p1, true)
            push!(turn, 1)
        else player == 2 && passed_p1[1] == false
            deleteat!(passed_p2, 1)
            push!(passed_p2, true)
            push!(turn, 1)
        end
    end
end
###########################################################
# RANDOM HAND, RANDOM PEGGING
###########################################################
full_rand_s_card = Int[]
full_rand_s_reason = String[]
full_rand_h_idx = Int[]
fr_p1_hands = String[]
fr_p2_hands = String[]
all_fr_h_idx = Int[]
for i in 1:10000000
    cribbage_deck = deck()
    shuffle!(cribbage_deck)

    hand1 = hand()
    hand2 = hand()
    # random 4 card hands
    move!(cribbage_deck, hand1, 4)
    move!(cribbage_deck, hand2, 4)

    push!(fr_p1_hands, join(hand1.cards))
    push!(fr_p2_hands, join(hand2.cards))
    push!(all_fr_h_idx, i)

    peg_pile = hand()
    pile_score = Int[]

    passed_p1 = Bool[false]
    passed_p2 = Bool[false]
    # init pegging
    push!(peg_pile, hand1.cards[1])
    push!(pile_score, hand1.cards[1].value)
    deleteat!(hand1, 1)

    turn = Int[2]

    while length(hand1.cards) + length(hand2.cards) > 0
        if sum(turn) % 2 == 0
            peg_random(hand2, peg_pile, pile_score, full_rand_s_card, full_rand_s_reason, full_rand_h_idx, passed_p1, passed_p2, turn, i)
        else
            peg_random(hand1, peg_pile, pile_score, full_rand_s_card, full_rand_s_reason, full_rand_h_idx, passed_p1, passed_p2, turn, i)
        end
    end
end

###########################################################
# RANDOM HAND, PEGGING STRATEGY
###########################################################
rh_s_card = Int[]
rh_s_reason = String[]
for i in 1:1000000
    cribbage_deck = deck()
    shuffle!(cribbage_deck)

    hand1 = hand()
    hand2 = hand()

    move!(cribbage_deck, hand1, 4)
    move!(cribbage_deck, hand2, 4)

    peg_pile = hand()
    pile_score = Int[]

    passed_p1 = Bool[false]
    passed_p2 = Bool[false]

    # init pegging
    lead_plays = Bool[]
    for c in 1:length(hand1.cards)
        push!(lead_plays, hand1.cards[c].value != 5)
    end
    if any(lead_plays)
        push!(peg_pile, hand1.cards[findfirst(x->x == true, lead_plays)])
        push!(pile_score, hand1.cards[findfirst(x->x == true, lead_plays)].value)
        deleteat!(hand1, findfirst(x->x == true, lead_plays))
    else
        push!(peg_pile, hand1.cards[1])
        push!(pile_score, hand1.cards[1].value)
        deleteat!(hand1, 1)
    end
    # sort hands high to low
    sort!(hand1.cards, rev = true)
    sort!(hand2.cards, rev = true)

    turn = Int[2]

    while length(hand1.cards) + length(hand2.cards) > 0
        if sum(turn) % 2 == 0
            peg(hand2, peg_pile, pile_score, rh_s_card, rh_s_reason, passed_p1, passed_p2, turn)
        else
            peg(hand1, peg_pile, pile_score, rh_s_card, rh_s_reason, passed_p1, passed_p2, turn)
        end
    end
end
###########################################################
# STRATEGY: PICK HIGHEST SCORING HAND & PEG FOR MOST POINTS
###########################################################
s_card = Int[]
s_reason = String[]
h_idx = Int[]
p1_hands = String[]
p2_hands = String[]
all_h_idx = Int[]
@time for i in 1:10000000
    cribbage_deck = deck()
    shuffle!(cribbage_deck)

    hand1 = hand()
    hand2 = hand()

    move!(cribbage_deck, hand1, 6)
    move!(cribbage_deck, hand2, 6)

    h1_combns = combinations(hand1.cards, 4) |> collect
    h2_combns = combinations(hand2.cards, 4) |> collect
    hand1 = score_hand(h1_combns)
    hand2 = score_hand(h2_combns)

    push!(p1_hands, join(hand1.cards))
    push!(p2_hands, join(hand2.cards))
    push!(all_h_idx, i)

    peg_pile = hand()
    pile_score = Int[]

    passed_p1 = Bool[false]
    passed_p2 = Bool[false]

    # init pegging
    lead_plays = Bool[]
    for c in 1:length(hand1.cards)
        push!(lead_plays, hand1.cards[c].value != 5)
    end
    if any(lead_plays)
        push!(peg_pile, hand1.cards[findfirst(x->x == true, lead_plays)])
        push!(pile_score, hand1.cards[findfirst(x->x == true, lead_plays)].value)
        deleteat!(hand1, findfirst(x->x == true, lead_plays))
    else
        push!(peg_pile, hand1.cards[1])
        push!(pile_score, hand1.cards[1].value)
        deleteat!(hand1, 1)
    end
    # sort hands high to low
    sort!(hand1.cards, rev = true)
    sort!(hand2.cards, rev = true)

    turn = Int[2]

    while length(hand1.cards) + length(hand2.cards) > 0
        if sum(turn) % 2 == 0
            peg(hand2, peg_pile, pile_score, s_card, s_reason, h_idx, passed_p1, passed_p2, turn, i)
        else
            peg(hand1, peg_pile, pile_score, s_card, s_reason, h_idx, passed_p1, passed_p2, turn, i)
        end
    end
end

###########################################################
# WRITE OUT RESULTS
###########################################################
CSV.write("D:/Projects/julia_plays_cribbage/full_random_card_scores.csv", DataFrame(score_card = full_rand_s_card, score_reason = full_rand_s_reason, h_idx = full_rand_h_idx))
CSV.write("D:/Projects/julia_plays_cribbage/full_random_hands.csv", DataFrame(p1_hand = fr_p1_hands, p2_hand = fr_p2_hands, h_idx = all_fr_h_idx))
#CSV.write("random_hand_control_data_prelim.csv", DataFrame(score_card = rh_s_card, score_reason = rh_s_reason))
CSV.write("D:/Projects/julia_plays_cribbage/strategy_score_data.csv", DataFrame(score_card = s_card, score_reason = s_reason, h_idx = h_idx))
CSV.write("D:/Projects/julia_plays_cribbage/strategy_hands.csv", DataFrame(p1_hand = p1_hands, p2_hand = p2_hands, h_idx = all_h_idx))


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

function configure(version)
    name = "GMP"

    hash = Dict(
        v"6.1.2" => "5275bb04f4863a13516b2f39392ac5e272f5e1bb8057b18aec1c9b79d73d8fb2",
        v"6.2.0" => "f51c99cb114deb21a60075ffb494c1a210eb9d7cb729ed042ddb7de9534451ea",
        v"6.2.1" => "eae9326beb4158c386e39a356818031bd28f3124cf915f8c5b1dc4c7a36b4d7c",
    )

    # Collection of sources required to complete build
    sources = [
        ArchiveSource("https://gmplib.org/download/gmp/gmp-$(version).tar.bz2", hash[version]),
        DirectorySource("./bundled"; follow_symlinks=true),
    ]

    # Bash recipe for building across all platforms
    script = raw"""
cd $WORKSPACE/srcdir/gmp-*

# Include Julia-carried patches
for f in $WORKSPACE/srcdir/patches/*.patch; do
    echo "Applying patch ${f}"
    atomic_patch -p1 ${f}
done

flags=(--enable-cxx --enable-shared --disable-static)

# On x86_64 architectures, build fat binary
if [[ ${proc_family} == intel ]]; then
    flags+=(--enable-fat)
fi

autoreconf
./configure --prefix=$prefix --build=${MACHTYPE} --host=${target} ${flags[@]}

make -j${nproc}
make install

# On Windows, we need to make sure that the non-versioned dll names exist too
if [[ ${target} == *mingw* ]]; then
    cp -v ${libdir}/libgmp-*.dll "${libdir}/libgmp.dll"
    cp -v ${libdir}/libgmpxx-*.dll "${libdir}/libgmpxx.dll"
fi

# GMP is dual-licensed, install all license files
install_license COPYING*
"""

    # We enable experimental platforms as this is a core Julia dependency
    platforms = expand_cxxstring_abis(supported_platforms(;experimental=true))

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

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

    return name, version, sources, script, platforms, products, dependencies
end

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
import Base: exp, log, sqrt, ^
import Base: cos, sin, tan, cosh, sinh, tanh, acos, asin, atan, acosh, asinh, atanh

function exp(q::Quaternion)
    z = sqrt(q.q1*q.q1 + q.q2*q.q2 + q.q3*q.q3)
    return exp(q.q0)*(cos(z) + imag(q)*(z != zero(q.q0) ? sin(z)/z : one(q.q0)))
end

function log(q::Quaternion)
    z = abs(imag(q))
    return log(abs(q)) + ((imag(q)!=zero(q)) ? imag(q)/z*atan(z,q.q0) :
    q.q0 < zero(q.q0) ? quaternion(pi*im) : zero(q))
end

function cos(q::Quaternion)
    z = abs(imag(q))
    return cos(q.q0)*cosh(z)-imag(q)*sin(q.q0)*(z != zero(q.q0) ? sinh(z)/z : one(q.q0))
end

function sin(q::Quaternion)
    z = abs(imag(q))
    return sin(q.q0)*cosh(z)+imag(q)*cos(q.q0)*(z != zero(q.q0) ? sinh(z)/z : one(q.q0))
end

function tan(q::Quaternion)
    return sin(q)/cos(q)
end

function cosh(q::Quaternion)
    return (exp(q)+exp(-q))/2
end

function sinh(q::Quaternion)
    return (exp(q)-exp(-q))/2
end

function tanh(q::Quaternion)
    return sinh(q)/cosh(q)
end

function sqrt(q::Quaternion)
    return exp(log(q)/2)
end

function (^)(q::Quaternion, p::AbstractFloat)
    return exp(p*log(q))
end

function (^)(q::Quaternion, p::Quaternion)
    return exp(p*log(q))
end

function acos(q::Quaternion)
    x = imag(q)
    y = abs(x)
    z = y!=zero(q.q0) ? x/y : quaternion(one(q.q0)*im)

    return -z*log(q + z*sqrt(one(q.q0) - q*q))
    # return -z*log(q - z*sqrt(one(q.q0) - q*q))
end

function asin(q::Quaternion)
    x = imag(q)
    y = abs(x)
    z = y!=zero(q.q0) ? x/y : quaternion(one(q.q0)*im)

    return -z*log(z*q + sqrt(one(q.q0) - q*q))
    #= return -z*log(z*q - sqrt(one(q.q0) - q*q)) =#
end

function atan(q::Quaternion)
    x = imag(q)
    y = abs(x)
    z = y!=zero(q.q0) ? x/y : quaternion(one(q.q0)*im)

    return -z*log((one(q.q0)+z*q)/(one(q.q0)-z*q))/2
end


text:
module ArbiterBase

import ..ArbiterGraphs: ArbiterGraph, NodeSet, ImmutableNodeSet
import ..ArbiterSchedulers: Scheduler, start_task!, end_task!, isfinished

export task_loop

immutable Results
    completed::ImmutableNodeSet
    failed::ImmutableNodeSet
end

Results(completed, failed) = Results(ImmutableNodeSet(completed), ImmutableNodeSet(failed))

immutable TaskResult
    name::Symbol
    successful::Bool
end

"""
The inner task loop for a task runner.

tasks: An iterable of tasks to run
execute: A function that runs a task. It should take a task as its
    sole argument, and returns a Nullable{TaskResult} which may be null.
wait: (optional, None) A function to run whenever there aren't any
    runnable tasks (but there are still tasks listed as running).
    If given, this function should take no arguments, and should
    return an iterable of TaskResults.
"""
function task_loop(tasks, execute::Function, wait::Nullable{Function})
    completed_tasks = NodeSet()
    failed_tasks = NodeSet()

    Scheduler(tasks, completed_tasks, failed_tasks) do scheduler
        while !isfinished(scheduler)
            task = start_task!(scheduler)

            while !isnull(task)
                let result = execute(get(task))

                    # result exists iff execute is synchronous
                    if !isnull(result)
                        let task_result = get(result)
                            end_task!(scheduler, task_result.name, task_result.successful)
                        end
                    end
                end

                task = start_task!(scheduler)
            end

            if !isnull(wait)
                for task_result in wait()
                    end_task!(scheduler, task_result.name, task_result.successful)
                end
            end
        end
    end

    return Results(completed_tasks, failed_tasks)
end

task_loop(tasks, execute::Function, wait::Function) = task_loop(tasks, execute, Nullable(wait))

task_loop(tasks, execute::Function) = task_loop(tasks, execute, Nullable{Function}())

end


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function syminit(app)

nd = app.dim;
ncq = app.ncq;
# ncmu1a = app.ncmu1a;
# ncmu1b = app.ncmu1b;
# ncmu2a = app.ncmu2a;
# ncmu2b = app.ncmu2b;
# ncmu2c = app.ncmu2c;
# ncmu3a = app.ncmu3a;
# ncmu3b = app.ncmu3b;
# ncmu3c = app.ncmu3c;
# ncmu4a = app.ncmu4a;
# ncmu4b = app.ncmu4b;
# ncmu4c = app.ncmu4c;

xij = [SymPy.symbols("xij$i") for i=1:nd];
xik = [SymPy.symbols("xik$i") for i=1:nd];
xil = [SymPy.symbols("xil$i") for i=1:nd];

xi = [SymPy.symbols("xi$i") for i=1:nd];
xj = [SymPy.symbols("xj$i") for i=1:nd];
xk = [SymPy.symbols("xk$i") for i=1:nd];
xl = [SymPy.symbols("xl$i") for i=1:nd];

qi = [SymPy.symbols("qi$i") for i=1:ncq];
qj = [SymPy.symbols("qj$i") for i=1:ncq];
qk = [SymPy.symbols("qk$i") for i=1:ncq];
ql = [SymPy.symbols("ql$i") for i=1:ncq];

ti = [SymPy.symbols("ti$i") for i=1:1];
tj = [SymPy.symbols("tj$i") for i=1:1];
tk = [SymPy.symbols("tk$i") for i=1:1];
tl = [SymPy.symbols("tl$i") for i=1:1];

ai = [SymPy.symbols("ai$i") for i=1:1];
aj = [SymPy.symbols("aj$i") for i=1:1];
ak = [SymPy.symbols("ak$i") for i=1:1];
al = [SymPy.symbols("al$i") for i=1:1];

# ti = SymPy.symbols("ti");
# tj = SymPy.symbols("tj");
# tk = SymPy.symbols("tk");
# tl = SymPy.symbols("tl");

# ai = SymPy.symbols("ai");
# aj = SymPy.symbols("aj");
# ak = SymPy.symbols("ak");
# al = SymPy.symbols("al");

# mu1a = [SymPy.symbols("mu1a$i") for i=1:ncmu1a];
# mu1b = [SymPy.symbols("mu1b$i") for i=1:ncmu1a];

# mu2a = [SymPy.symbols("mu2a$i") for i=1:ncmu2a];
# mu2b = [SymPy.symbols("mu2b$i") for i=1:ncmu2b];
# mu2c = [SymPy.symbols("mu2c$i") for i=1:ncmu2c];

# mu3a = [SymPy.symbols("mu3a$i") for i=1:ncmu3a];
# mu3b = [SymPy.symbols("mu3b$i") for i=1:ncmu3b];
# mu3c = [SymPy.symbols("mu3c$i") for i=1:ncmu3c];

# mu4a = [SymPy.symbols("mu4a$i") for i=1:ncmu4a];
# mu4b = [SymPy.symbols("mu4b$i") for i=1:ncmu4b];
# mu4c = [SymPy.symbols("mu4c$i") for i=1:ncmu4c];

return xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al
#, mu1a, mu1b, mu2a, mu2b, mu2c, mu3a, mu3b, mu3c, mu4a, mu4b, mu4c

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
module BabyNeuralODE

export read_halfmoons, train

using Flux
using Random

function ode_solve(fun, u, t1, N)
    δt = t1/N
    for _ in 1:N
        u += fun(u)*δt
    end
    u
end

function read_halfmoons(fname, n)
    to_float(x) = parse(Float64, x)
    open(fname) do f
        map(readlines(f)) do l
            (x1, x2, y) = split(l) .|> to_float
            ([x1, x2], y)
        end
    end
end

function test(m, data)
    acc = 0
    for (x, y) in data
        (p0, p1) = m(x)
        yy = p0>p1 ? 0 : 1
        if yy == y
            acc += 1
        end
    end
    acc / length(data)
end

function train(data, typ)
    Random.shuffle!(data)
    ndata = length(data)
    ntrain = Int(round(0.8*ndata))
    ntest = ndata - ntrain

    nhidden = 16

    w2 = param(zeros(nhidden, nhidden))
    b2 = param(zeros(nhidden))
    odelayer(x) = ode_solve(x->(w2*x .+ b2), x, 1, 10)

    if typ == 1
        layer2 = Dense(nhidden, nhidden, σ)
    elseif typ == 2
        layer2 = odelayer
    else
        layer2 = x->x
    end

    m = Chain(
        Dense(2, nhidden, σ),
        layer2,
        Dense(nhidden, 2, σ),
        softmax
    )

    function loss(x, y)
        y = Flux.onehot(y, 0:1)
        Flux.crossentropy(m(x), y)
    end

    pm = Flux.params(m, w2, b2)
    opt = Flux.Optimise.ADAM()

    println("epoch: ",  0,
            "\tloss: ", sum(loss(d...) for d in data),
            "\ttest: ", test(m, data))

    batchsize = 400

    acc = Float64[]
    for epoch in 1:100
        Random.shuffle!(@view data[1:ntrain])

        i = 1
        while i<ntrain
            Flux.train!(loss, pm, data[i:(i+batchsize-1)], opt)
            i += batchsize
        end

        t = test(m, data[ntrain+1:end])
        println("epoch: ",  epoch,
                "\tloss: ", sum(loss(data[i]...) for i in 1:ntrain),
                "\ttest: ", t)
        push!(acc, t)
        t > 0.999 && break
    end

    acc
end

end # module


text:
#real number/dedekind cut formulas.
@universal function lub(x::Ubound)
  is_exact(x.upper) && return x.upper
  is_positive(x.upper) && return outer_exact(x.upper)
  return inner_exact(x.upper)
end

@universal function glb(x::Ubound)
  is_exact(x.lower) && return x.lower
  is_positive(x.lower) && return inner_exact(x.lower)
  return outer_exact(x.lower)
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
# This is for testing using Symata code from Julia

## These just call apprules, so they are in general no more efficient than calling from Symata directly.
## In contrast, the trig functions below filter numeric arguments and are exactly as efficient as pure Julia
## when called from within Julia functions for which the compiler can infer types.
## Similar filters could be written for the functions wrapped here.
##
for f in (:Expand, :Factor, :Head, :Take, :Simplify, :Integrate, :DirichletEta, :Times, :List, :Plus, :Exp, :Power, :Length)
    @eval ($f)(mx...) = apprules(mxpr($(QuoteNode(f)),mx...))
    @eval export $f
end

## Use Julia inference and dispatch to call numeric functions directly with no overhead when possible.
## TODO: Handle two arg functions, etc.
for f in (:cos, :sin, :abs, :tan, :exp, :log, (:acos, :ArcCos) , (:asin, :ArcSin), (:atan, :ArcTan ),
          :cot, :cosh, :sinh, :tanh, :sqrt, :erf, :erfc, :gamma, :zeta)
    local uf
    if isa(f, Tuple)
        uf = f[2]
        f = f[1]
    else
        uf = Symbol(uppercasefirst(string(f)))
    end
    @eval ($uf)(x::AbstractFloat) = ($f)(x)
    @eval ($uf)(x::AbstractArray{T}) where {T<:AbstractFloat} = ($f)(x)
    @eval ($uf)(x::Complex{T}) where {T<:AbstractFloat} = ($f)(x)
    @eval ($uf)(x) = apprules(mxpr($(QuoteNode(uf)),x))
    @eval export $uf
end

# Integrate(mx::Mxpr,symorlist) = apprules(mxpr(:Integrate,mx,symorlist))
# export Integrate

const Pi = :Pi
export Pi

const E = Base.MathConstants.e     # could make this the symbol ?
export E

#### Arithmetic operators

## methods for Julia math functions that operate
## on symbols and Mxpr. Some of these are used
## in Symata code. They can also be used at the Julia
## repl.

# SJSym is an alias of Symbol

"""
    _symatamath()

define math methods in the Symata module that operate on symbols. This allows Julia expressions such as
```
:a + :b
:a - 3
:c^2
```

Note: we now call syamtamath() after its definition. It defines these operators in the Symata module scope.

obsolete: These methods are disabled by default because they extend `Base` methods for `Base` types. These methods are
in general reserved for definition in future versions of Julia. There is some chance that these will be given conflicting
definitions in the Julia base language in the future. However, we know of no such plans at present.

These are all binary methods between Julia `Number`s and `Symbol`s extending `*`, `+`, `-`, and `^`.

There are similar methods for arithmetic operators between numbers or symbols and Symata expressions, but
these are defined by default.
"""
function _symatamath()
    @eval begin
        Base.:*(a::Number,b::SJSym) = mxpr(:Times,a,b)  # why not mmul ? TODO: try using mmul, etc. here
        Base.:*(a::SJSym,b::SJSym) = mxpr(:Times,a,b)
        Base.:*(a::SJSym,b::Number) = mxpr(:Times,b,a)
        Base.:+(a::SJSym,b::Number) = mxpr(:Plus,b,a)
        Base.:+(a::SJSym,b::SJSym) = mxpr(:Plus,a,b)
        Base.:+(a::Number,b::SJSym) = mxpr(:Plus,a,b)
        Base.:-(a::Number,b::SJSym) = mplus(a, mxpr(:Times,-1,b))
        Base.:-(a::SJSym,b::Number) = mplus(a, -b)
        Base.:-(a::SJSym,b::SJSym) = mplus(a, mxpr(:Times,-1,b))
        Base.:/(a::SJSym,b::SJSym) = mmul(a, mmul(-1,b))
        Base.:/(a::Number,b::SJSym) = mmul(a, mmul(-1,b))
        Base.:/(a::SJSym,b::Number) = mmul(a, mmul(-1,b))
        Base.:^(base::SJSym,expt::Integer) = mxpr(:Power,base,expt)
        Base.:^(base::SJSym,expt) = mxpr(:Power,base,expt)
    end
    nothing
end

_symatamath()

"""
    symatamath()

defines Base methods for arithmetic operators `*`, `+`, `-`, `^` between `Symbols` and between `Symbols` and `Numbers`. `symatamath()` is *not*
necessary when running Symata and using J(expr) to evaluate Julia code.  `J()` evaluates in the Symata
module where these methods are already defined.  `symatamath()` is useful in Julia mode where expressions
are evaluated in  `Main`.

These methods are not defined in `Main` by default because defining `Base` methods between core Julia objects could
conflict with future versions of Julia.
"""
function symatamath()
    @eval begin
        Base.:*(a::Number,b::SJSym) = mxpr(:Times,a,b)  # why not mmul ?
        Base.:*(a::SJSym,b::SJSym) = mxpr(:Times,a,b)
        Base.:*(a::SJSym,b::Number) = mxpr(:Times,b,a)
        Base.:+(a::SJSym,b::Number) = mxpr(:Plus,b,a)
        Base.:+(a::SJSym,b::SJSym) = mxpr(:Plus,a,b)
        Base.:+(a::Number,b::SJSym) = mxpr(:Plus,a,b)
        Base.:-(a::Number,b::SJSym) = mplus(a, mxpr(:Times,-1,b))
        Base.:-(a::SJSym,b::Number) = mplus(a, -b)
        Base.:-(a::SJSym,b::SJSym) = mplus(a, mxpr(:Times,-1,b))
        Base.:^(base::SJSym,expt::Integer) = mxpr(:Power,base,expt)
        Base.:^(base::SJSym,expt) = mxpr(:Power,base,expt)
    end
    nothing
end

## Arithmetic methods involving annotated Symata types. These will never conflict with Base Julia,
## so they are safe to define.

Base.:*(a::Mxpr,b::Mxpr) = mxpr(:Times,a,b)
Base.:*(a::Mxpr,b) = mxpr(:Times,a,b)
Base.:*(a,b::Mxpr) = mxpr(:Times,a,b)

# Base.:*(a::Mxpr,b::Mxpr) = mxpr(:Times,a,b)
# Base.:*(a::Mxpr,b) = mxpr(:Times,a,b)
# Base.:*(a,b::Mxpr) = mxpr(:Times,a,b)

Base.:+(a::Mxpr,b::Mxpr) = mxpr(:Plus,a,b)
Base.:+(a::Mxpr,b) = mxpr(:Plus,a,b)
Base.:+(a,b::Mxpr) = mxpr(:Plus,a,b)

# Base.:+(a::Mxpr,b::Mxpr) = mxpr(:Plus,a,b)
# Base.:+(a::Mxpr,b) = mxpr(:Plus,a,b)
# Base.:+(a,b::Mxpr) = mxpr(:Plus,a,b)

Base.:-(a,b::Mxpr) = mxpr(:Plus,a,mxpr(:Times,-1,b))
Base.:-(a::Mxpr) = mxpr(:Times,-1,a)

Base.:^(base::Mxpr,expt::Integer) = mxpr(:Power,base,expt)
Base.:^(base::Mxpr,expt) = mxpr(:Power,base,expt)

Base.:/(a::Mxpr,b) = mxpr(:Times,a,mxpr(:Power,b,-1))

## (Aug 16, 2018) Upgrading for Julia v1.0 from v0.7 Comment these
## out. I think this code is very old
# Base.:-(a,b::Mxpr) = mxpr(:Plus,a,mxpr(:Times,-1,b))
# Base.:-(a::Mxpr) = mxpr(:Times,-1,a)
# Base.:^(base::Mxpr,expt::Integer) = mxpr(:Power,base,expt)
# Base.:^(base::Mxpr,expt) = mxpr(:Power,base,expt)
# Base.:/(a::Mxpr,b) = mxpr(:Times,a,mxpr(:Power,b,-1))

## Symata uses module-local functions * + - ^ /
## For anything not defined in Symata, the Base methods are called.
##
## (Aug 16, 2018) Upgrading for Julia v1.0 from v0.7
## Julia v1.0 says I am extending Base on the LHS
## and that I must import the Symbols. That is:
## Base.:*(args...) = Base.:*(args...) Makes no sense.
## What's going on here ?
## I will comment thest out for now,.

# *(args...) = Base.:*(args...)
# +(args...) = Base.:+(args...)
# -(args...) = Base.:-(args...)
# ^(args...) = Base.:^(args...)
# /(args...) = Base.:/(args...)

## This was commented out long before (Aug 16, 2018)
# Already defined elsewhere (... where ?)
# I = im
# export I

# TODO: make an infix assignment operator... hm or  macro


text:
#-----------------------------------------------------------------------------------------------------
# Tensor Basis
#-----------------------------------------------------------------------------------------------------
"""
    TensorBasis

Basis without any symmetries.

Properties:
-----------
- dgt : Vector{Int}, Digits.
- B   : Int, levels on each site.
"""
struct TensorBasis <: AbstractBasis
    dgt::Vector{Int}
    B::Int
end
eltype(::TensorBasis) = Int
#-----------------------------------------------------------------------------------------------------
export tensorbasis
tensorbasis(L::Integer; base::Integer=2) = TensorBasis(zeros(Int, L), base)
#-----------------------------------------------------------------------------------------------------
function change!(b::TensorBasis, i::Integer)::Int
    change!(b.dgt, i, base=b.B)
    1
end
#-----------------------------------------------------------------------------------------------------
function index(b::TensorBasis)::Tuple{Int, Int}
    1, index(b.dgt, base=b.B)
end
#-----------------------------------------------------------------------------------------------------
size(b::TensorBasis, i::Integer) = (i == 2 || i == 1) ? b.B ^ length(b.dgt) : 1
size(b::TensorBasis) = (l=b.B^length(b.dgt); (l, l))
#-----------------------------------------------------------------------------------------------------
# Helper Functions
#-----------------------------------------------------------------------------------------------------
"""
    index(dgt::AbstractVector{<:Integer})

Digit ⟹ index.

The method compute the polynomial

number = bits[i] * base^(L-i) + 1

in the most efficient way.
"""
function index(dgt::AbstractVector{Int}; base::Int=2)::Int
    N = 0
    for i = 1:length(dgt)
        N = muladd(base, N, dgt[i])
    end
    N + 1
end

function index(dgt::AbstractVector{Int}, sites::Vector{Int}; base::Int=2)::Int
    N = 0
    for i in sites
        N = muladd(base, N, dgt[i])
    end
    N + 1
end
#-----------------------------------------------------------------------------------------------------
"""
    change!(dgt::Digit, index::Integer)

Index ⟹ Digit.

The method compute the bits vector and write to bits.
"""
function change!(dgt::Vector{Int}, ind::Int; base::Int=2)
    N = ind - 1
    for i = length(dgt):-1:1
        N, dgt[i] = divrem(N, base)
    end
end

function change!(dgt::Vector{Int}, sites::Vector{Int}, ind::Int; base::Int=2)
    N = ind - 1
    for i = length(sites):-1:1
        N, dgt[sites[i]] = divrem(N, base)
    end
end


text:
"""
This function returns a vector of the products provided by the project
based on the data provided in the input files.
"""
function create_products(simulation_data::AgentSimulationData,
                         projectdata::DataFrames.DataFrameRow)
    products = Product[]

    markets = get_markets(simulation_data)

    reserve_up_mkt_df = read_data(joinpath(get_data_dir(get_case(simulation_data)), "markets_data", "reserve_up_mkt_param.csv"))
    reserve_down_mkt_df = read_data(joinpath(get_data_dir(get_case(simulation_data)), "markets_data", "reserve_down_mkt_param.csv"))

    reserve_up_eligible_categories = reserve_up_mkt_df[findfirst(x-> x == projectdata["Zone"], reserve_up_mkt_df[:, "zones"]), "eligible categories"]
    reserve_down_eligible_categories = reserve_down_mkt_df[findfirst(x-> x == projectdata["Zone"], reserve_down_mkt_df[:, "zones"]), "eligible categories"]

    # Energy market product
    variable_cost = projectdata["Fuel Price \$/MMBTU"] * projectdata["HR_avg_0"] / 1000
        push!(products, Energy(:Energy,
                               Dict{String, Array{Float64, 2}}(),
                               variable_cost))

    if markets[:Reserves] && occursin(projectdata["Category"], reserve_up_eligible_categories) > 0
        push!(products, OperatingReserve{ReserveUpEMIS}(:ReserveUp,
                                        projectdata["Ramp Rate pu/Hr"],
                                        0.0))
    end

    if markets[:Reserves] && occursin(projectdata["Category"], reserve_down_eligible_categories) > 0
        push!(products, OperatingReserve{ReserveDownEMIS}(:ReserveDown,
                                        projectdata["Ramp Rate pu/Hr"],
                                        0.0))
    end

    if markets[:Capacity] && projectdata["Capacity Eligible"]
        push!(products, Capacity(:Capacity,
                                 0.0,
                                 Dict{String, Array{Float64, 1}}(),
                                 0.0))
    end

    if markets[:REC] && projectdata["REC Eligible"]
        push!(products, REC(:REC,
                            0.0,
                            0.0))
    end

    return products
end

"""
This function returns a vector of the products provided by the project
based on the data provided in the input files and PSY Device data.
"""
function create_products(simulation_data::AgentSimulationData,
                         device::T,
                         projectdata::DataFrames.DataFrameRow) where T <: Union{PSY.Generator, PSY.Storage}

    products = Product[]

    markets = get_markets(simulation_data)

    service_types = typeof.(PSY.get_services(device))

    # Energy market product
    variable_cost = projectdata["Fuel Price \$/MMBTU"] * projectdata["HR_avg_0"] / 1000.0

    push!(products, Energy(:Energy,
                               Dict{String, Array{Float64, 2}}(),
                               variable_cost))

    fields = collect(fieldnames(typeof(device)))

    reserve_up_lim = 1.0
    reserve_down_lim = 1.0

    if in(:ramp_limits, fields)
        ramp_limits = PSY.get_ramp_limits(device)
        rating = get_device_size(device)
        if !isnothing(ramp_limits)
            reserve_up_lim = ramp_limits[:up] * 60 / rating
            reserve_down_lim = ramp_limits[:down] * 60 / rating
        end
    end

    if in(PSY.VariableReserve{PSY.ReserveUp}, service_types)
        push!(products, OperatingReserve{ReserveUpEMIS}(:ReserveUp,
                                        reserve_up_lim,
                                        0.0))
    end

    if in(PSY.VariableReserve{PSY.ReserveDown}, service_types)
        push!(products, OperatingReserve{ReserveDownEMIS}(:ReserveDown,
                                        reserve_down_lim,
                                        0.0))
    end

    if markets[:Capacity] && projectdata["Capacity Eligible"]
        push!(products, Capacity(:Capacity,
                                 0.0,
                                 Dict{String, Array{Float64, 1}}(),
                                 0.0))
    end

    if markets[:REC] && projectdata["REC Eligible"]
        push!(products, REC(:REC,
                            0.0,
                            0.0))
    end

    return products
end