text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#

#-------------------------------------- Default Renderer -------------------------------

initializeVisualization(renderer::Modia3D.AbstractRenderer, allVisuElements::Vector{Object3D}) = error("No renderer defined.")
visualize!(renderer::Modia3D.AbstractRenderer, time::Float64) = error("No renderer defined.")
closeVisualization(renderer::Modia3D.AbstractRenderer)        = error("No renderer defined.")

#-------------------------------------- Default Contact Detection -------------------------------

"""
    ContactPairs(nzmax, collSuperObjs, noCPairs, dummyObject3D)

Generate a new ContactPairs structure used for communication between the Object3D handler and a ContactDetection handler.

- `DummyObject3D::Modia3D.AbstractObject3DFeature`: A dummy Object3D that can be used in the struct as element of a vector of Object3Ds
  to fill the vector with a dummy value of the correct type.
"""
mutable struct ContactPairs
    # Solid shapes used in contact detection (provided by Object3D handler)
    allowedToMove::Vector{Union{Bool,Nothing}}
    dummyObject3D::Modia3D.AbstractObject3DFeature     # Dummy Object3D for non-used elements of z-Vector.

    # Dimensions
    ne::Int                                         # length(collSuperObjs)
    nz::Int                                         # length(z)
    nzContact::Int                                  # length(z | z has contact) length of z where zi has contact

    function ContactPairs(world::Composition.Object3D, AABB::Vector{Vector{Basics.BoundingBox}}, superObjs::Vector{SuperObjsRow},
                        allowedToMove::Vector{Union{Bool,Nothing}}, visualizeBoundingBox::Bool, nVisualContSupPoints::Int,
                        visualizeContactPoints::Bool, visualizeSupportPoints::Bool, defaultContactSphereDiameter::Float64)
        @assert(length(superObjs) > 0)
        @assert(nVisualContSupPoints > 0)
        dummyObject3D = Composition.emptyObject3DFeature

        lengthCollSuperObjs = length(superObjs)
        if lengthCollSuperObjs <= 1
            error("There's nothing to collide. All Object3Ds are rigidly connected.")
        end

        nzContact = 0
        nz = 2

        if lengthCollSuperObjs > 2
            nz = lengthCollSuperObjs
        end

        amountVisuPoints = 1
        amountSuppPoints = 1
        if visualizeContactPoints
            addContactVisuObjToWorld!(world, nVisualContSupPoints, defaultContactSphereDiameter)
        end
        if visualizeSupportPoints
            addSupportVisuObjToWorld!(world, nVisualContSupPoints, defaultContactSphereDiameter)
        end

        if visualizeBoundingBox
            addAABBVisuToWorld!(world, AABB)
        end

        new(allowedToMove, dummyObject3D, lengthCollSuperObjs, nz, nzContact)
    end
end

initializeContactDetection!(ch::Modia3D.AbstractContactDetection, collSuperObjs::Vector{Vector{Modia3D.AbstractObject3DFeature}}, noCPairs::Vector{Vector{Int64}}) = error("No contact detection handler defined.")
selectContactPairsWithEvent!(ch::Modia3D.AbstractContactDetection)   = error("No contact detection handler defined.")
selectContactPairsNoEvent!(ch::Modia3D.AbstractContactDetection)     = error("No contact detection handler defined.")
getDistances!(ch::Modia3D.AbstractContactDetection)                  = error("No contact detection handler defined.")
setComputationFlag(ch::Modia3D.AbstractContactDetection)             = error("No contact detection handler defined.")
closeContactDetection!(ch::Modia3D.AbstractContactDetection)         = error("No contact detection handler defined.")

include(joinpath(Modia3D.path, "src", "contactDetection", "ContactDetectionMPR", "ContactDetectionMPR_handler.jl"))

#-------------------------------------- Gravity field ----------------------------------

"""
    gravityField = NoGravityField()

Generate an instance of type `NoGravityField` that defines no gravity.
"""
struct NoGravityField <: Modia3D.AbstractGravityField
   gvec::SVector{3,Float64} # [m/s^2] Vector of gravity acceleration
   NoGravityField() = new(SVector{3,Float64}(0.0, 0.0, 0.0))
end
gravityAcceleration(grav::NoGravityField, r_abs::SVector{3,Float64})::SVector{3,Float64} = grav.gvec

"""
    UniformGravityField(;g=9.81, n=[0,-1,0])

Generate an instance of type `UniformGravityField` that defines
a uniform gravity field with gravitational acceleration `g` in
direction `n`.

# Arguments
- `g::Float64`: Gravitational acceleration
- `n::AbstractVector`: Direction of gravitational acceleration

# Example
```julia
import Modia3D

grav = Modia3D.UniformGravityField()
   r = Modia3D.EarthRadius
   g = gravityAcceleration(grav,r)  # g is independent of r
```
"""
struct UniformGravityField <: Modia3D.AbstractGravityField
   gvec::SVector{3,Float64} # [m/s^2] Vector of gravity acceleration

   function UniformGravityField(;g=9.81, n=[0,-1,0])
      @assert(g >= 0.0)
      new(g*normalize(n))
   end
end
gravityAcceleration(grav::UniformGravityField, r_abs::SVector{3,Float64})::SVector{3,Float64} = grav.gvec

const G           = 6.67408e-11  # [m3/(kg.s2)]  Newtonian constant of gravitation (https://en.wikipedia.org/wiki/Gravitational_constant)
const EarthMass   = 5.9722e24    # [kg]          Mass of earth (https://en.wikipedia.org/wiki/Earth_mass)
const EarthRadius = 6.3781e6     # [m]           Radius of earth (https://en.wikipedia.org/wiki/Earth_radius)

"""
    PointGravityField([mass|;mue=G*EarthMass])

Generate an instance of type `PointGravityField` that defines
a point gravity field of `mass` or gravity field constant `mue`.
The gravity center is located at the world origin, i.e. the
gravitational acceleration decreases quadratically with distance
from world.

# Example
```julia
import Modia3D

grav = Modia3D.PointGravityField()   # Gravity field of earth
   r = Modia3D.EarthRadius
   g = gravityAcceleration(grav,r)  # g is singular at r=0
```
"""
struct PointGravityField <: Modia3D.AbstractGravityField
   mue::Float64    # [m3/s2] Gravity field constant (default = field constant of earth)

   function PointGravityField(mass)
      @assert(mass >= 0.0)
      mue = G*mass
      new(mue)
   end
   function PointGravityField(;mue = G*EarthMass)
      @assert(mue >= 0.0)
      new(mue)
   end
end
gravityAcceleration(grav::PointGravityField, r_abs::SVector{3,Float64})::SVector{3,Float64} = -(grav.mue/dot(r_abs,r_abs))*normalize(r_abs)

"""
    sceneUpDir = upwardsDirection(gravityField::Modia3D.AbstractGravityField)

Determine upwards direction of scene depending on main gravity field direction.

sceneUpDir: -/+1 = -/+x-direction; -/+2 = -/+y-direction; -/+3 = -/+z-direction
"""
function upwardsDirection(gravityField::Modia3D.AbstractGravityField)
    sceneUpDir = +2
    if typeof(gravityField) == UniformGravityField
        gv = gravityField.gvec
        if norm(gv) > 1e-8
            if (abs(gv[1]) > abs(gv[2])) && (abs(gv[1]) > abs(gv[3]))
                if gv[1] < 0.0
                    sceneUpDir = +1
                else
                    sceneUpDir = -1
                end
            elseif (abs(gv[2]) > abs(gv[3])) && (abs(gv[2]) > abs(gv[1]))
                if gv[2] < 0.0
                    sceneUpDir = +2
                else
                    sceneUpDir = -2
                end
            else
                if gv[3] < 0.0
                    if gv[3] < 0.0
                        sceneUpDir = +3
                    else
                        sceneUpDir = -3
                    end
                end
            end
        end
    end
    return sceneUpDir
end

function camRotation(beta, alpha)
    # beta = longitude angle about y-axis of view ref frame
    # alpha = latitude angle about -x-axis of beta-rotated frame
    # R_cam = transformation matrix v_cam = R_cam * v_view
    cbe = cos(beta)
    sbe = sin(beta)
    cal = cos(alpha)
    sal = sin(alpha)
    R_cam = [cbe  0.0  -sbe; -sal*sbe  cal  -sal*cbe; cal*sbe  sal  cal*cbe]
end

"""
    (r_Camera, R_Camera) = cameraPosition(distance, longitude, latitude, cameraUpDir, sceneUpDir)

Determine camera/light position and orientation.

distance:    Distance between world frame and camera position.
longitude:   Longitude angle of camera position (0 = -y/-z/-x direction).
latitude:    Latitude angle of camera position (0 = horizontal).
cameraUpDir: Upwards direction of camera.
sceneUpDir:  Upwards direction of scene.

r_Camera: Absolute position vector of camera/light.
R_Camera: Absolute rotation matrix of camera/light.
"""
function cameraPosition(distance, longitude, latitude, cameraUpDir, sceneUpDir)
    R_Camera = camRotation(longitude, latitude)
    r_Camera = R_Camera' * [0.0, 0.0, 1.0] * distance
    if sceneUpDir == -1
        R_view = [0 1 0; -1 0 0; 0 0 1]  # scene upwards := -x direction
    elseif sceneUpDir == 1
        R_view = [0 0 1; +1 0 0; 0 1 0]  # scene upwards := +x direction
    elseif sceneUpDir == -2
        R_view = [0 0 1; 0 -1 0; 1 0 0]  # scene upwards := -y direction
    elseif sceneUpDir == -3
        R_view = [1 0 0; 0 0 -1; 0 1 0]  # scene upwards := -z direction
    elseif sceneUpDir == +3
        R_view = [0 1 0; 0 0 +1; 1 0 0]  # scene upwards := +z direction
    else
        R_view = [1 0 0; 0 +1 0; 0 0 1]  # scene upwards := +y direction
    end
    if cameraUpDir == 1
        R_view = [0 1 0; 0 0 1; 1 0 0] * R_view  # camera upwards := +x direction
    elseif cameraUpDir == 3
        R_view = [0 0 1; 1 0 0; 0 1 0] * R_view  # camera upwards := +z direction
    else
        # camera upwards := +y direction
    end
    r_Camera = R_view' * r_Camera
    R_Camera = R_Camera * R_view
    return (r_Camera, R_Camera)
end

#-------------------------------------- Global Scene Options -------------------------------
struct SceneOptions <: Modia3D.AbstractSceneOptions

    # Gravity field
    gravityField::Modia3D.AbstractGravityField

    # Multibody structure
    useOptimizedStructure::Bool    # = true, if the optimized structure (with super objects, and common inertia) is used

    # Contact detection
    contactDetection::Modia3D.AbstractContactDetection
    nVisualContSupPoints::Int               # amount of visual contact or support points
    gap::Float64
    enableContactDetection::Bool            # = true, if contact detection is enabled
    defaultContactSphereDiameter::Float64   # = true, if contact points are visualized
    elasticContactReductionFactor::Float64  # c_res_used = c_res * elasticContactReductionFactor (> 0)

    # Visual defaults
    nominalLength::Float64         # [m]     Nominal length of 3D system
    defaultFrameLength::Float64    # [m]     Default for frame length if visualizeFrames = true (but not world frame)
    defaultJointLength::Float64    # [m]     Default for the fixed length of a shape representing a joint
    defaultJointWidth::Float64     # [m]     Default for the fixed width of a shape representing a joint
    defaultForceLength::Float64    # [m]     Default for the fixed length of a shape representing a force (e.g., damper)
    defaultForceWidth::Float64     # [m]     Default for the fixed width of a shape representing a force (e.g., spring, bushing)
    defaultBodyDiameter::Float64   # [m]     Default for diameter of sphere representing the center of mass of a body
    defaultWidthFraction::Float64  #         Default for shape width as a fraction of shape length
    defaultArrowDiameter::Float64  # [m]     Default for arrow diameter (e.g., of forces, torques, sensors)
    defaultN_to_m::Float64         # [N/m]   Default scaling of force arrows (length = force/defaultN_to_m)
    defaultNm_to_m::Float64        # [N.m/m] Default scaling of torque arrows (length = torque/defaultNm_to_m)

    # Animation
    enableVisualization::Bool             # = true, if online animation is enabled
    animationFile::Union{Nothing,String}  # path&name of animation file
    visualizeGravity::Bool                # = true, if gravity field shall be visualized (acceleration vector or field center)
    visualizeFrames::Bool                 # = true, if all frames shall be visualized
    visualizeConvexHulls::Bool            # = true, if convex hulls (used for contact detection) shall be visualized
    visualizeBoundingBox::Bool            # = true, if AABB's are visualized
    visualizeContactPoints::Bool          # = true, if contact points shall be visualized
    visualizeSupportPoints::Bool          # = true, if support points shall be visualized
    cameraDistance::Float64               # Distance between world frame and camera position
    cameraLongitude::Float64              # Longitude angle of camera position (0 = -y/-z/-x direction)
    cameraLatitude::Float64               # Latitude angle of camera position (0 = horizontal)
    lightDistance::Float64                # Distance between world frame and light position
    lightLongitude::Float64               # Longitude angle of light position (0 = -y/-z/-x direction)
    lightLatitude::Float64                # Latitude angle of light position (0 = horizontal)

    function SceneOptions(;gravityField                  = UniformGravityField(),
                           useOptimizedStructure         = true,
                           contactDetection              = ContactDetectionMPR_handler(),
                           nVisualContSupPoints          = 5,
                           gap                           = 0.001,
                           enableContactDetection        = true,
                           defaultContactSphereDiameter  = 0.1,
                           elasticContactReductionFactor = 1.0,
                           nominalLength                 = 1.0,
                           defaultFrameLength            = 0.2*nominalLength,
                           defaultJointLength            = nominalLength/10,
                           defaultJointWidth             = nominalLength/20,
                           defaultForceLength            = nominalLength/10,
                           defaultForceWidth             = nominalLength/20,
                           defaultBodyDiameter           = nominalLength/9,
                           defaultWidthFraction          = 20,
                           defaultArrowDiameter          = nominalLength/40,
                           defaultN_to_m                 = 1000,
                           defaultNm_to_m                = 1000,
                           enableVisualization           = true,
                           animationFile                 = nothing,
                           visualizeGravity              = true,
                           visualizeFrames               = false,
                           visualizeConvexHulls          = true,
                           visualizeBoundingBox          = false,
                           visualizeContactPoints        = false,
                           visualizeSupportPoints        = false,
                           cameraDistance                = 10.0*nominalLength,
                           cameraLongitude               = 30/180*pi,
                           cameraLatitude                = 15/180*pi,
                           lightDistance                 = 10.0*nominalLength,
                           lightLongitude                = 60/180*pi,
                           lightLatitude                 = 45/180*pi)
        @assert(nVisualContSupPoints > 0)
        @assert(gap > 0.0)
        @assert(defaultContactSphereDiameter > 0.0)
        @assert(nominalLength > 0.0)
        @assert(defaultFrameLength > 0.0)
        @assert(defaultJointLength > 0.0)
        @assert(defaultJointWidth > 0.0)
        @assert(defaultForceLength > 0.0)
        @assert(defaultForceWidth > 0.0)
        @assert(defaultBodyDiameter >= 0.0)
        @assert(defaultArrowDiameter > 0.0)
        @assert(defaultN_to_m > 0.0)
        @assert(defaultNm_to_m  > 0.0)
        @assert(cameraDistance > 0.0)
        @assert(lightDistance > 0.0)

        sceneOptions = new(gravityField,
            useOptimizedStructure,
            contactDetection,
            nVisualContSupPoints,
            gap,
            enableContactDetection,
            defaultContactSphereDiameter,
            elasticContactReductionFactor,
            nominalLength,
            defaultFrameLength,
            defaultJointLength,
            defaultJointWidth,
            defaultForceLength,
            defaultForceWidth,
            defaultBodyDiameter,
            defaultWidthFraction,
            defaultArrowDiameter,
            defaultN_to_m,
            defaultNm_to_m,
            enableVisualization,
            animationFile,
            visualizeGravity,
            visualizeFrames,
            visualizeConvexHulls,
            visualizeBoundingBox,
            visualizeContactPoints,
            visualizeSupportPoints,
            cameraDistance,
            cameraLongitude,
            cameraLatitude,
            lightDistance,
            lightLongitude,
            lightLatitude)

        sceneOptions.contactDetection.visualizeContactPoints       = visualizeContactPoints
        sceneOptions.contactDetection.visualizeSupportPoints       = visualizeSupportPoints
        sceneOptions.contactDetection.defaultContactSphereDiameter = defaultContactSphereDiameter
        return sceneOptions
    end
end

struct animationData
    position::SVector{3,Float32}    # abs. Object3D position
    quaternion::SVector{4,Float32}  # abs. Object3D quaternion
end

struct animationStep
    time::Float32                      # simulation time
    objectData::Vector{animationData}  # animation data of Object3Ds
end

#-------------------------------------- Scene -------------------------------
"""
    Scene(; kwargs...)

Defines global properties of the system, such as the gravity field. Exactly one [Object3D](@ref) must have a `Scene` feature defined. This Object3D is used as inertial system (world, root) and is not allowed to have a parent Object3D.

| Keyword arguments             | defaults                                |
|:------------------------------|:----------------------------------------|
| `gravityField`                  | [`UniformGravityField`](@ref)`()`         |
| `enableVisualization`           | true                                    |
| `visualizeFrames`               | false                                   |
| `nominalLength`                 | 1.0                                     |
| `defaultFrameLength`            | 0.2*`nominalLength`                     |
| `useOptimizedStructure`         | true                                    |
| `enableContactDetection`        | true                                    |
| `contactDetection`              | [`ContactDetectionMPR_handler`](@ref)`()`|
| `elasticContactReductionFactor` | 1.0                                     |
| `visualizeBoundingBox`          | false                                   |
| `visualizeContactPoints`        | false                                   |
| `visualizeSupportPoints`        | false                                   |
| `nVisualContSupPoints`          | 5                                       |
| `defaultContactSphereDiameter`  | 0.1                                     |
| `animationFile`                 | nothing                                 |
| `cameraDistance`                | 10.0*`nominalLength`                      |
| `cameraLongitude`               | 30/180*pi                               |
| `cameraLatitude`                | 15/180*pi                               |
| `lightDistance`                 | 10.0*`nominalLength`                      |
| `lightLongitude`                | 60/180*pi                               |
| `lightLatitude`                 | 45/180*pi                               |

# Arguments
- `gravityField::Modia3D.AbstractGravityField`: Gravity field of the scene. Supported values:
   - [`NoGravityField`](@ref),
   - [`UniformGravityField`](@ref),
   - [`PointGravityField`](@ref)

- `enableVisualization::Bool`: = true, to enable online animation with DLR SimVis. If SimVis is not installed, this flag has no effect.

- `visualizeFrames::Bool`: = true, to visualize the coordinate system of every [Object3D](@ref) that is not explicitly switched off.

- `nominalLength::Float64`: Nominal length in [m].

- `defaultFrameLength::Float64`: Default frame length in [m] for visualizing little [CoordinateSystem](@ref).

- `useOptimizedStructure::Bool`: = true, if pre-processing the whole system. For example, computing the common mass, common center of mass and common inertia tensor of all rigidly connected Object3Ds that have mass properties.

- `enableContactDetection::Bool`: = true, if contact detection is enable, see [Collision Handling](@ref).

- `contactDetection::Modia3D.AbstractContactDetection`: Handler used for contact detection e.g., to determine the smallest distance between two objects.

- `elasticContactReductionFactor::Float64`: (> 0.0)
  - ``usedContactCompliance = contactCompliance * elasticContactReductionFactor``

- `visualizeBoundingBox::Bool`: Flag enabled for visualizing Axis Aligned Bounding Box (AABB) for all solid shapes allowed to collide

- `visualizeContactPoints::Bool`: Flag enabled for visualizing contact points, and `enableVisualization = true`

- `visualizeSupportPoints::Bool`: Flag enabled for visualizing support points, and `enableVisualization = true`

- `nVisualContSupPoints::Int64`: defines how many contact as well as support points are visualized (there is a warning if you'll need to increase it)

- `defaultContactSphereDiameter::Float64`: diameter of [Sphere](@ref) in [m] used for contact and support point visualization

- `animationFile::String`: only if a valid path and name of the animation file is set (it must be a .json file) a json file is exported

- `cameraDistance::Float64`: distance between world frame and camera position

- `cameraLongitude::Float64`: longitude angle of camera position (0 = -y/-z/-x direction)

- `cameraLatitude::Float64`: latitude angle of camera position (0 = horizontal)

- `lightDistance::Float64`: distance between world frame and light position

- `lightLongitude::Float64`: longitude angle of light position (0 = -y/-z/-x direction)

- `lightLatitude::Float64`: latitude angle of light position (0 = horizontal)
"""
mutable struct Scene
    name::String                              # model name
    autoCoordsys::Shapes.CoordinateSystem   # Coordinate system that is automatically included (e.g. due to visualizeFrames=true)
    stack::Vector{Object3D}                   # Stack to traverse objs
    buffer::Vector{Object3D}                  # stores all roots of a super obj

    options::SceneOptions                     # Global options defined for the scene

    visualize::Bool                           # = true, if visualization elements available
    collide::Bool                             # = true, if elements for contact detection available

    # initialization of analysis
    initAnalysis::Bool                        # = true, if analysis is initialized
    initSuperObj::Bool                        # = true, if super objects are initialized
    initMassComp::Bool                        # = true, if mass computation for optimized structure is initialized
    analysis::Modia3D.AnalysisType            # Type of analysis
    superObjs::Vector{SuperObjsRow}          # super objects

    treeAccVelo::Vector{Object3D}
    tree::Vector{Object3D}                    # Spanning tree of the frames in depth-first order (without world)
    treeForComputation::Vector{Object3D}      # Tree that is used for the computation
                                              # treeForComputation = options.useOptimizedStructure ? treeAccVelo : tree

    #cutJoints::Vector{Modia3D.AbstractJoint}  # Vector of all cut-joints
    allVisuElements::Vector{Object3D}         # all Object3Ds that should be visualized (for communiating with SimVis)
    updateVisuElements::Vector{Object3D}      # all Object3Ds that are visibly only e.g. visualization frames (must be updated first)
    allCollisionElements::Vector{Object3D}    # all Object3Ds, which are allowed to collide (no order, no super objects)
    noCPairs::Vector{Vector{Int64}}           # Indices of frames (with respect to collSuperObjs) that can't collide in general (e.g. objects are connected via joints)
    noCPairsHelp::Dict{Modia3D.AbstractJoint,Vector{Int64}}
    allowedToMove::Vector{Union{Bool,Nothing}}
    AABB::Vector{Vector{Basics.BoundingBox}}    # Bounding boxes of elements that can collide
    #forceElements::Vector{Int64}
    exportAnimation::Bool                     # animation file export is enabled
    animation::Vector{animationStep}          # animation data of visible Object3Ds

    # Data specific to a particular joint type
    revolute::Vector{Revolute}
    prismatic::Vector{Prismatic}
    freeMotion::Vector{FreeMotion}

    function Scene(sceneOptions::SceneOptions = SceneOptions())
        exportAnimation = false
        if !isnothing(sceneOptions.animationFile)
            (base, ext) = splitext(sceneOptions.animationFile)
            if ext == ".json"
                exportAnimation = true
            else
                @warn("Extension of animationFile=$(sceneOptions.animationFile) is not 'json'.\n-> Animation export is disabled.")
            end
        end
        new("Scene",
            Shapes.CoordinateSystem(length=sceneOptions.defaultFrameLength),
            Vector{Object3D}[],
            Vector{Object3D}[],
            sceneOptions,
            false,
            false,
            false,
            false,
            false,
            Modia3D.KinematicAnalysis,
            Vector{SuperObjsRow}[],
            Vector{Object3D}[],
            Vector{Object3D}[],
            Vector{Object3D}[],
            Vector{Object3D}[],
            Vector{Object3D}[],
            Vector{Object3D}[],
            Vector{Vector{Int64}}[],
            Dict{Modia3D.AbstractJoint,Vector{Int64}}(),
            Vector{Union{Bool}}[],
            Vector{Vector{Basics.BoundingBox}}[],
            exportAnimation,
            Vector{animationStep}[],

            Revolute[],
            Prismatic[],
            FreeMotion[])
    end
end


text:
"""
    MPIWorkerManager([nprocs])

A [`ClusterManager`](https://docs.julialang.org/en/v1/stdlib/Distributed/#Distributed.ClusterManager)
using the MPI.jl launcher
[`mpiexec`](https://juliaparallel.github.io/MPI.jl/stable/environment/#MPI.mpiexec).

The workers will all belong to an MPI session, and can communicate using MPI
operations. Note that unlike `MPIManager`, the MPI session will not be
initialized, so the workers will need to `MPI.Init()`.

The master process (pid 1) is _not_ part of the session, and will communicate
with the workers via TCP/IP.

# Usage

    using Distributed, MPIClusterManager

    mgr = MPIWorkerManager(4) # launch 4 MPI workers
    mgr = MPIWorkerManager() # launch the default number of MPI workers (determined by `mpiexec`)

    addprocs(mgr; kwoptions...)

The following `kwoptions` are supported:

 - `dir`: working directory on the workers.

 - `mpiexec`: MPI launcher executable (default: use the launcher from MPI.jl)

 - `mpiflags`: additional flags  to pass to `mpiexec`

 - `exename`: Julia executable on the workers.

 - `exeflags`: additional flags to pass to the Julia executable.

 - `threadlevel`: the threading level to initialize MPI. See
 [`MPI.Init()`](https://juliaparallel.github.io/MPI.jl/stable/environment/#MPI.Init)
 for details.

 - `topology`: how the workers connect to each other.

 - `enable_threaded_blas`: Whether the workers should use threaded BLAS.

 - `master_tcp_interface`: Server interface to listen on. This allows direct
   connection from other hosts on same network as specified interface
   (otherwise, only connections from `localhost` are allowed).
"""
mutable struct MPIWorkerManager <: ClusterManager
    "number of MPI processes"
    nprocs::Union{Int, Nothing}
    "map `MPI.COMM_WORLD` rank to Julia pid"
    mpi2j::Dict{Int,Int}
    "map Julia pid to `MPI.COMM_WORLD` rank"
    j2mpi::Dict{Int,Int}
    "are the processes running?"
    launched::Bool
    "have the workers been initialized?"
    initialized::Bool
    "notify this when all workers registered"
    cond_initialized::Condition
    "redirected ios from workers"
    stdout_ios::Vector{IO}

    function MPIWorkerManager(nprocs = nothing)
        mgr = new(nprocs,
                  Dict{Int,Int}(),
                  Dict{Int,Int}(),
                  false,
                  false,
                  Condition(),
                  IO[]
                  )

        return mgr
    end

end

Distributed.default_addprocs_params(::MPIWorkerManager) =
    merge(Distributed.default_addprocs_params(),
          Dict{Symbol,Any}(
                :mpiexec        => nothing,
                :mpiflags       => ``,
                :master_tcp_interface => nothing,
                :threadlevel    => :serialized,
          ))

# Launch a new worker, called from Base.addprocs
function Distributed.launch(mgr::MPIWorkerManager,
                            params::Dict,
                            instances::Array,
                            cond::Condition)

    mgr.launched && error("MPIWorkerManager already launched. Create a new instance to add more workers")

    master_tcp_interface = params[:master_tcp_interface]

    if mgr.nprocs === nothing
        configs = WorkerConfig[]
    else
        configs = Vector{WorkerConfig}(undef, mgr.nprocs)
    end

    # Set up listener
    port, server = if !isnothing(master_tcp_interface)
        # Listen on specified server interface
        # This allows direct connection from other hosts on same network as
        # specified interface.
        listenany(getaddrinfo(master_tcp_interface), 11000) # port is just a hint
    else
        # Listen on default interface (localhost)
        # This precludes direct connection from other hosts.
        listenany(11000)
    end
    ip = getsockname(server)[1]

    connections = @async begin
        while isnothing(mgr.nprocs) || length(mgr.stdout_ios) < mgr.nprocs
            io = accept(server)
            config = WorkerConfig()
            config.io = io
            config.enable_threaded_blas = params[:enable_threaded_blas]

            # Add config to the correct slot so that MPI ranks and
            # Julia pids are in the same order
            rank = Serialization.deserialize(io)
            config.ident = (rank=rank,)
            nprocs = Serialization.deserialize(io)
            if mgr.nprocs === nothing
                if nprocs === nothing
                    error("Could not determine number of processes")
                end
                mgr.nprocs = nprocs
                resize!(configs, nprocs)
            end
            configs[rank+1] = config
            push!(mgr.stdout_ios, io)
        end
    end

    # Start the workers
    cookie = Distributed.cluster_cookie()
    setup_cmds = "using Distributed; import MPIClusterManagers; MPIClusterManagers.setup_worker($(repr(string(ip))),$(port),$(repr(cookie)); threadlevel=$(repr(params[:threadlevel])))"
    MPI.mpiexec() do mpiexec
        mpiexec = something(params[:mpiexec], mpiexec)
        mpiflags = params[:mpiflags]
        if !isnothing(mgr.nprocs)
            mpiflags = `$mpiflags -n $(mgr.nprocs)`
        end
        exename = params[:exename]
        exeflags = params[:exeflags]
        dir = params[:dir]
        mpi_cmd = Cmd(`$mpiexec $mpiflags $exename $exeflags -e $setup_cmds`, dir=dir)
        open(detach(mpi_cmd))
    end
    mgr.launched = true

    # wait with timeout (https://github.com/JuliaLang/julia/issues/36217)
    launch_timeout = Distributed.worker_timeout()
    timer = Timer(launch_timeout) do t
        schedule(connections, InterruptException(), error=true)
    end
    try
        wait(connections)
    catch e
        error("Could not connect to workers")
    finally
        close(timer)
    end

    # Append our configs and notify the caller
    append!(instances, configs)
    notify(cond)
end

function Distributed.manage(mgr::MPIWorkerManager, id::Integer, config::WorkerConfig, op::Symbol)
    if op == :register
        rank = config.ident.rank
        mgr.j2mpi[id] = rank
        mgr.mpi2j[rank] = id
        if length(mgr.j2mpi) == mgr.nprocs
            # All workers registered
            mgr.initialized = true
            notify(mgr.cond_initialized)
        end
    elseif op == :deregister
        # TODO: Sometimes -- very rarely -- Julia calls this `deregister`
        # function, and then outputs a warning such as """error in running
        # finalizer: ErrorException("no process with id 3 exists")""". These
        # warnings seem harmless; still, we should find out what is going wrong
        # here.
    elseif op == :interrupt
        # TODO: This should never happen if we rmprocs the workers properly
        @assert false
    elseif op == :finalize
        # This is called from within a finalizer after deregistering; do nothing
    else
        @assert false # Unsupported operation
    end
end

ext: jl


text:
using Documenter, JuliaSyntax

makedocs(;
    modules=[JuliaSyntax],
    format=Documenter.HTML(),
    pages=[
        "Overview" => "index.md",
        "API Reference" => "reference.md",
        "Design Discussion" => "design.md",
    ],
    repo="https://github.com/c42f/JuliaSyntax.jl/blob/{commit}{path}#L{line}",
    sitename="JuliaSyntax.jl",
    authors = "Chris Foster and contributors: https://github.com/c42f/JuliaSyntax.jl/graphs/contributors"
)

deploydocs(;
    repo="github.com/c42f/JuliaSyntax.jl",
    push_preview=true
)


text:
# CSCI5573Final/task.jl
# Licensed under the MIT License. See LICENSE.md file in the project root for
# full license information.

mutable struct _Task
    ID::Int64
    Children::Vector{Int64}
    Type::Symbol
    Cost::Int64
    Complexity::Float64

    Dependencies::Vector{Int64}

    _Task(ID::Int64, Children::Vector{Int64}, Type::Symbol, Cost::Int64=0, Complexity::Float64=0.0) =
        new(ID, Children, Type, Cost, Complexity, Vector{Int64}())

end

function working!(task::_Task, n::Int64)
    if task.Cost - n <= 0
        task.Cost = 0
    else
        task.Cost = task.Cost - n
    end
    return task.Cost
end

isDone(task::_Task) = task.Cost === 0

add_dependency!(task::_Task, ID::Int64) = push!(task.Dependencies, ID)
remove_dependency!(task::_Task, ID::Int64) = deleteat!(task.Dependencies, task.Dependencies .== ID)

function ParseTask(cfg::String, final=false)
    config = split(cfg, " ")

    ID = parse(Int, config[2])
    Children = Vector{Int64}()
    Type = Symbol(config[4])

    Cost = Symbol(config[5]) === :nan ? 100 : parse(Float64, config[5])/1_000_000_000 |> round |> Int64
    Cost = Cost > 0 ? Cost : 1

    Complexity = parse(Float64, config[6])

    if !final
        Children = split(config[3], ",") |> s->String.(s) |> i->parse.(Int, i)
    end

    return _Task(ID, Children, Type, Cost, Complexity)
end


text:
using ExtensibleUnions
using Test
using Traceur

macro genstruct!(list)
    x = gensym()
    y = quote
        struct $(x) end
        push!($(list), $(x))
    end
    return y
end

# macro genfunc!(list)
#     x = gensym()
#     y = quote
#         function $(x) end
#         push!($(list), $(x))
#     end
#     return y
# end

@testset "ExtensibleUnions.jl" begin
    @testset "Unit tests" begin
        @testset "test_code_transformation.jl" begin
            include("test_code_transformation.jl")
        end
        @testset "test_extensible_functions.jl" begin
            include("test_extensible_functions.jl")
        end
        @testset "test_extensible_unions.jl" begin
            include("test_extensible_unions.jl")
        end
        @testset "test_inferred.jl" begin
            include("test_inferred.jl")
        end
        @testset "test_traceur.jl" begin
            include("test_traceur.jl")
        end
        @testset "test_update_methods.jl" begin
            include("test_update_methods.jl")
        end
    end
    @testset "Integration tests" begin
        @testset "test_examples.jl" begin
            include("test_examples.jl")
        end
    end
    @testset "reset.jl" begin
        @testset "ExtensibleUnions._unsafe_reset!()" begin
            ExtensibleUnions._unsafe_reset!()
            @test isempty(ExtensibleUnions._registry_extensibleunion_to_genericfunctions)
            @test isempty(ExtensibleUnions._registry_extensibleunion_to_members)
            @test isempty(ExtensibleUnions._registry_genericfunctions_to_extensibleunions)
        end
    end
end


text:
using FinancialMath
using Test

@testset "FinancialMath" begin
  include("core.jl")
end


text:
module BrokenSideEffect

using Matte

const title = "MatteExample :: Broken Side-effects"

function ui()
    sidebar_layout(
        side_panel(
            h1("Push this button"),
            br(),
            button("my_button", "Push me"),
            br(),
            slider("slider", "Choose a number", 1, 10)
        ),
        main_panel(
            p("You have pushed the button: ", text_output("my_count"), " times"),
            br(),
            p("Times the slider is equal to: ", text_output("my_calc"))
        )
    )
end

function register_session_vars()
    Dict(
        :count => 0
    )
end

module Server

function my_calc(slider, session)
    slider * session.count
end

function my_count(my_button, session)
    if my_button
        session.count += 1
    end
    session.count
end

end

end


text:
module TestData
    using Base.Test
    using DataArrays
    using DataFrames

    #test_group("DataVector creation")
    dvint = @data([1, 2, NA, 4])
    dvint2 = DataArray([5:8])
    dvint3 = DataArray(5:8)
    dvflt = @data([1.0, 2, NA, 4])
    dvstr = @data(["one", "two", NA, "four"])
    dvdict = DataArray(Dict,4)    # for issue #199

    #test_group("constructors")
    df1 = DataFrame({dvint, dvstr}, ["Ints", "Strs"])
    df2 = DataFrame({dvint, dvstr})
    df3 = DataFrame({dvint})
    df4 = DataFrame([1:4 1:4])
    df5 = DataFrame({@data([1,2,3,4]), dvstr})
    df6 = DataFrame({dvint, dvint, dvstr}, ["A", "B", "C"])
    df7 = DataFrame(x = dvint, y = dvstr)
    @assert size(df7) == (4, 2)
    @assert isequal(df7["x"], dvint)

    #test_group("description functions")
    @assert nrow(df6) == 4
    @assert ncol(df6) == 3
    @assert all(names(df6) .== ["A", "B", "C"])
    @assert all(names(df2) .== ["x1", "x2"])
    @assert all(names(df7) .== ["x", "y"])

    #test_group("ref")
    @assert df6[2, 3] == "two"
    @assert isna(df6[3, 3])
    @assert df6[2, "C"] == "two"
    @assert isequal(df6["B"], dvint)
    @assert ncol(df6[[2,3]]) == 2
    @assert nrow(df6[2,:]) == 1
    @assert size(df6[[1, 3], [1, 3]]) == (2, 2)
    @assert size(df6[1:2, 1:2]) == (2, 2)
    @assert size(head(df6,2)) == (2, 3)
    # lots more to do

    #test_group("combining")

    dfc = hcat(df3, df4)
    @assert ncol(dfc) == 3
    @assert all(names(dfc) .== ["x1", "x1_1", "x2"])
    @assert isequal(dfc["x1"], df3["x1"])

    @assert isequal(dfc, [df3 df4])

    dfr = vcat(df4, df4)
    @assert nrow(dfr) == 8
    @assert all(names(df4) .== names(dfr))
    @assert isequal(dfr, [df4, df4])

    dfr = vcat(df2, df3)
    @assert size(dfr) == (8,2)
    @assert all(names(df2) .== names(dfr))
    @assert isna(dfr[8,"x2"])

    #test_group("show")
    # @assert repr(df1) == "4x2 DataFrame:\n        Ints   Strs\n[1,]       1  \"one\"\n[2,]       2  \"two\"\n[3,]      NA     NA\n[4,]       4 \"four\"\n"

    #test_group("assign")
    df6[3] = @data(["un", "deux", "troix", "quatre"])
    @assert df6[1, 3] == "un"
    df6["B"] = [4, 3, 2, 1]
    @assert df6[1,2] == 4
    df6["D"] = [true, false, true, false]
    @assert df6[1,4] == true
    delete!(df6, "D")
    @assert all(names(df6) .== ["A", "B", "C"])
    @assert ncol(df6) == 3

    #test_group("NA handling")
    @assert nrow(df5[complete_cases(df5), :]) == 3

    #test_context("SubDataFrames")

    #test_group("constructors")
    # single index is rows
    sdf6a = sub(df6, 1)
    sdf6b = sub(df6, 2:3)
    sdf6c = sub(df6, [true, false, true, false])
    @assert size(sdf6a) == (1,3)
    sdf6d = sub(df6, [1,3], "B")
    @assert size(sdf6d) == (2,1)

    #test_group("ref")
    @assert sdf6a[1,2] == 4

    #test_context("Within")
    #test_group("Associative")

    srand(1)
    a1 = [:a => [1, 2], :b => [3, 4], :c => [5, 6]]
    a2 = ["a" => [1, 2], "b" => [3, 4], "c" => [5, 6]]
    a3 = {"a" => [1, 2], "b" => [3, 4], :c => [5, 6]}

    @assert isequal(with(a1, :(c + 1)), a1[:c] + 1)
    @assert isequal(with(a2, :(c + 1)), with(a1, :(c + 1)))
    @assert isequal(with(a3, :(c + 1)), with(a1, :(c + 1)))
    @assert isequal(with(a3, :(c + 1 + 0 * b)), with(a1, :(c + 1)))

    a4 = within(a1, :( d = a + b ))
    @assert isequal(a4[:d], a1[:a] + a1[:b])
    @assert isequal(a4[:a], a1[:a])

    a4c = @transform(copy(a1), d => a + b )
    @assert isequal(a4[:d], a4c[:d])

    a4 = within(a2, :( d = a + b ))
    @assert isequal(a4["d"], a2["a"] + a2["b"])
    @assert isequal(a4["a"], a2["a"])

    # a4c = @transform(copy(a2), d => a + b )
    # @assert isequal(a4c["d"], a4["d"])

    a4 = within(a3, :( d = a + b ))
    @assert isequal(a4[:d], a3["a"] + a3["b"])
    @assert isequal(a4["a"], a3["a"])

    # Note: The following won't work.
    #       @transform will only find the keys that are symbols.
    ## a4c = @transform(copy(a3), d => a + b )
    ## @assert isequal(a4c[:d], a4[:d])

    a4 = based_on(a1, :( d = a + b ))
    @assert isequal(a4[:d], a1[:a] + a1[:b])

    a4 = based_on(a2, :( d = a + b ))
    @assert isequal(a4["d"], a2["a"] + a2["b"])

    a4 = based_on(a3, :( d = a + b ))
    @assert isequal(a4[:d], a3["a"] + a3["b"])

    #test_group("DataFrame")

    srand(1)
    N = 20
    d1 = PooledDataArray(rand(1:2, N))
    d2 = (@pdata ["A", "B", NA])[rand(1:3, N)]
    d3 = DataArray(randn(N))
    d4 = DataArray(randn(N))
    df7 = DataFrame({d1, d2, d3}, ["d1", "d2", "d3"])

    @assert isequal(with(df7, :(d3 + d3)), df7["d3"] + df7["d3"])
    @assert isequal(with(df7, :(d3 + $d4)), df7["d3"] + d4)
    x = df7 |> with(:( d3 + d3 ))
    @assert isequal(x, df7["d3"] + df7["d3"])

    df8 = within(df7, :(d4 = d3 + d3 + 1))
    @assert isequal(df7, df8[1:3])
    @assert isequal(df8["d4"], df7["d3"] + df7["d3"] + 1)
    within!(df8, :( d4 = d1 ))
    @assert isequal(df8["d1"], df8["d4"])

    df8 = @transform(copy(df7), d4 => d3 + 1)
    @assert isequal(df7, df8[1:3])

    df8 = based_on(df7, :( d1 = d3 ))
    @assert isequal(df8["d1"], df7["d3"])
    df8 = df7 |> based_on(:( d1 = d3 ))
    @assert isequal(df8["d1"], df7["d3"])
    df8 = based_on(df7, :( sum_d3 = sum(d3) ))
    @assert isequal(df8[1,1], sum(df7["d3"]))

    #@assert all(df7[:( d2 .== "B" )]["d1"] .== PooledDataArray([1,2,1,1]))
    # TODO: Remove all tests that depend upon srand(), which was just changed.
    # TODO: Restore this test
    # @assert all(df7[:( d2 .== "B" ), "d1"] .== PooledDataArray([2,1,1,1,1,2,2,1,2]))

    #test_group("groupby")

    gd = groupby(df7, "d1")
    @assert length(gd) == 2
    # @assert isequal(gd[2]["d2"], PooledDataVector["A", "B", NA, "A", NA, NA, NA, NA])
    @assert sum(gd[2]["d3"]) == sum(df7["d3"][dropna(df7["d1"] .== 2)])

    g1 = groupby(df7, ["d1", "d2"])
    g2 = groupby(df7, ["d2", "d1"])
    @assert sum(g1[1]["d3"]) == sum(g2[1]["d3"])

    res = 0.0
    for x in g1
        res += sum(x["d1"])
    end
    @assert res == sum(df7["d1"])

    df8 = df7 |> groupby(["d2"]) |> :( d3sum = sum(d3); d3mean = mean(dropna(d3)) )
    @assert isequal(df8["d2"], @pdata([NA, "A", "B"]))

    df9 = based_on(groupby(df7, "d2"),
                   :( d3sum = sum(d3); d3mean = mean(dropna(d3)) ))
    @assert isequal(df9, df8)

    df8 = within(groupby(df7, "d2"),
                 :( d4 = d3 + 1; d1sum = sum(d1) ))
    @assert all(df8[:( d2 .== "C" ), "d1sum"] .== 13)

    ## @assert isequal(with(g1, :( sum(d1) )), map(x -> sum(x["d1"]), g1))

    df8 = colwise(df7[[1, 3]], :sum)
    @assert df8[1, "d1_sum"] == sum(df7["d1"])

    df8 = colwise(groupby(df7, "d2"), [:sum, :length])
    @assert nrow(df8) == 3
    @assert ncol(df8) == 5
    #@assert df8[1, "d1_sum"] == 13
    # @assert df8[2, "d1_length"] == 7
    @assert df8[2, "d1_length"] == 8

    df9 = df7 |> groupby(["d2"]) |> [:sum, :length]
    @assert isequal(df9, df8)
    df9 = by(df7, "d2", [:sum, :length])
    @assert isequal(df9, df8)

    #test_group("reshape")

    d1 = DataFrame(quote
        a = [1:3]
        b = [1:4]
        c = randn(12)
        d = randn(12)
    end)

    d1c = @DataFrame(a => [1:3],
                     b => [1:4],
                     c => randn(12),
                     d => randn(12))

    @assert isequal(d1[1:2], d1c[1:2])

    d1s = stack(d1, ["a", "b"])
    d1s2 = stack(d1, ["c", "d"])
    d1s3 = melt(d1, ["c", "d"])
    @assert isequal(d1s[1:12, "c"], d1["c"])
    @assert isequal(d1s[13:24, "c"], d1["c"])
    @assert all(names(d1s) .== ["variable", "value", "c", "d"])
    @assert isequal(d1s, d1s3)
    d1s_df = stack_df(d1, ["a", "b"])
    # TODO: Fix this
    @assert isequal(d1s["variable"], d1s_df["variable"][:])
    @assert isequal(d1s["value"], d1s_df["value"][:])
    @assert isequal(d1s["c"], d1s_df["c"][:])
    @assert isequal(d1s[1,:], d1s_df[1,:])

    d1s["idx"] = [1:12, 1:12]
    d1s2["idx"] = [1:12, 1:12]
    d1us = unstack(d1s, "variable", "idx", "value")
    d1us2 = unstack(d1s2, "variable", "idx", "value")
    @assert isequal(d1us["a"], d1["a"])
    @assert isequal(d1us2["d"], d1["d"])

    d = DataFrame(quote
        a = letters[5:8]
        b = LETTERS[10:11]
        c = LETTERS[13 + [1, 1, 2, 2, 2, 1, 1, 2, 1, 2, 2, 1, 1, 2]]
        d = pi * [1:14]
    end)

    dpv = pivot_table(d, ["a", "b"], "c", "d")
    @assert( dpv[1,"O"] == d[5,"d"])
    @assert( nrow(dpv) == 4 )

    dpv2 = pivot_table(d, ["a"], ["c", "b"], "d")
    @assert( dpv2[1,"O_J"] == d[5,"d"])

    dpv3 = pivot_table(d, ["a"], ["c", "b"], "d", length)
    @assert( dpv3[1,"O_J"] == 1.0)

    #test_group("merge")

    srand(1)
    df1 = DataFrame(quote
        a = shuffle!([1:10])
        b = ["A","B"][rand(1:2, 10)]
        v1 = randn(10)
    end)

    df2 = DataFrame(quote
        a = shuffle!(reverse([1:5]))
        b2 = ["A","B","C"][rand(1:3, 5)]
        v2 = randn(3)    # test unequal lengths in the constructor
    end)

    m1 = join(df1, df2, on = "a")
    @assert isequal(m1["a"], @data([1, 2, 3, 4, 5]))
    # TODO: Re-enable
    # m2 = join(df1, df2, on = "a", kind = :outer)
    # @assert isequal(m2["b2"], DataVector["A", "B", "B", "B", "B", NA, NA, NA, NA, NA])
    # @assert isequal(m2["b2"], DataVector["B", "B", "B", "C", "B", NA, NA, NA, NA, NA])

    df1 = DataFrame({"a" => [1, 2, 3],
                     "b" => ["America", "Europe", "Africa"]})
    df2 = DataFrame({"a" => [1, 2, 4],
                     "c" => ["New World", "Old World", "New World"]})
    m1 = join(df1, df2, on = "a", kind = :inner)
    @assert isequal(m1["a"], @data([1, 2]))
    m2 = join(df1, df2, on = "a", kind = :left)
    @assert isequal(m2["a"], @data([1, 2, 3]))
    m3 = join(df1, df2, on = "a", kind = :right)
    @assert isequal(m3["a"], @data([1, 2, 4]))
    # TODO: Re-enable
    # m4 = join(df1, df2, on = "a", kind = :outer)
    # @assert isequal(m4["a"], DataVector[1, 2, 3, 4])

    # # test with NAs (issue #185)
    # df1 = DataFrame()
    # df1["A"] = DataVector["a", "b", "a", NA]
    # df1["B"] = DataVector[1, 2, 1, 3]

    # df2 = DataFrame()
    # df2["A"] = DataVector["a", NA, "c"]
    # df2["C"] = DataVector[1, 2, 4]

    # m1 = join(df1, df2, on = "A")
    # @assert size(m1) == (3,3)
    # @assert isequal(m1["A"], DataVector[NA,"a","a"])
    # m2 = join(df1, df2, on = "A", kind = :outer)
    # @assert size(m2) == (5,3)
    # @assert isequal(m2["A"], DataVector[NA,"a","a","b","c"])

    srand(1)
    df1 = DataFrame(quote
        a = ["x","y"][rand(1:2, 10)]
        b = ["A","B"][rand(1:2, 10)]
        v1 = randn(10)
    end)

    df2 = DataFrame(quote
        a = ["x","y"][[1,2,1,1,2]]
        b = ["A","B","C"][[1,1,1,2,3]]
        v2 = randn(5)
    end)
    df2[1,"a"] = NA

    # # TODO: Restore this functionality
    # m1 = join(df1, df2, on = ["a","b"])
    # @assert isequal(m1["a"], DataArray(["x", "x", "y", "y", fill("x", 5)]))
    # m2 = join(df1, df2, on = ["a","b"], kind = :outer)
    # @assert isequal(m2[10,"v2"], NA)
    # @assert isequal(m2["a"], DataVector["x", "x", "y", "y", "x", "x", "x", "x", "x", "y", NA, "y"])

    # m1a = join(within(df1, :(key = PooledDataArray(_DF[["a","b"]]))),
    #            based_on(df2, :(key = PooledDataArray(_DF[["a","b"]]); v2 = v2)),
    #            on = "key")
    # m2a = join(within(df1, :(key = PooledDataArray(_DF[["a","b"]]))),
    #            based_on(df2, :(key = PooledDataArray(_DF[["a","b"]]); v2 = v2)),
    #            on = "key",
    #            kind = :outer)
    # @assert isequal(sort(m1["b"]), sort(m1a["b"]))

    srand(1)
    function spltdf(d)
        d["x1"] = map(x -> x[1], d["a"])
        d["x2"] = map(x -> x[2], d["a"])
        d["x3"] = map(x -> x[3], d["a"])
        d
    end
    df1 = DataFrame(quote
        a = ["abc","abx", "axz", "def", "dfr"]
        v1 = randn(5)
    end)
    df1 = spltdf(df1)
    df2 = DataFrame(quote
        a = ["def", "abc","abx", "axz", "xyz"]
        v2 = randn(5)
    end)
    df2 = spltdf(df2)

    # m1 = join(df1, df2, on = "a")
    # m2 = join(df1, df2, on = ["x1", "x2", "x3"])
    # @assert isequal(sort(m1["a"]), sort(m2["a"]))

    #test_group("New DataVector constructors")
    dv = DataArray(Int, 5)
    @assert all(isna(dv))
    dv = DataArray(Float64, 5)
    @assert all(isna(dv))
    dv = @data(zeros(5))
    @assert all(dv .== 0.0)
    dv = @data(ones(5))
    @assert all(dv .== 1.0)

    # No more NA corruption
    dv = @data(ones(10_000))
    @assert !any(isna(dv))

    PooledDataArray(falses(2), falses(2))
    PooledDataArray(falses(2), trues(2))

    # Test vectorized comparisons work for DataVector's and PooledDataVector's
    @data([1, 2, NA]) .== 1
    @pdata([1, 2, NA]) .== 1
    @data(["1", "2", NA]) .== "1"
    @pdata(["1", "2", NA]) .== "1"

    # Test unique()
    #test_group("unique()")
    # TODO: Restore this
    # dv = DataArray(1:4)
    # dv[4] = NA
    # @assert (1 in unique(dv))
    # @assert (2 in unique(dv))
    # @assert (3 in unique(dv))
    # @assert (NA in unique(dv))

    #test_group("find()")
    dv = DataArray([true, false, true])
    @assert isequal(find(dv), [1, 3])

    pdv = PooledDataArray([true, false, true])
    @assert isequal(find(pdv), [1, 3])

    dv[1] = NA
    @assert isequal(find(dv), [3])

    pdv[1] = NA
    @assert isequal(find(pdv), [3])
end


text:
abstract type AbstractNote end

"""
    Note(pitch, velocity, position, duration, channel = 0) <: AbstractNote
Mutable data structure describing a music note. A bundle of many notes results
in the [`Notes`](@ref) struct, which is the output of the [`getnotes`](@ref)
function.

If the `channel` of the note is `0` (default), it is not shown.

You can also create a `Note` with the following keyword constructor:
```julia
Note(pitch, position; velocity = 100, duration = 960, channel = 0)
```

## Fields:
* `pitch::UInt8` : Pitch, starting from C-1 = 0, adding one per semitone.
  Use the functions [`name_to_pitch`](@ref) and
  [`pitch_to_name`](@ref) for integer and string representations.
* `velocity::UInt8` : Dynamic intensity. Cannot be higher than 127 (0x7F).
* `position::UInt` : Position in absolute time (since beginning of track), in ticks.
* `duration::UInt` : Duration in ticks.
* `channel::UInt8 = 0` : Channel of the track that the note is played on.
  Cannot be higher than 127 (0x7F).
"""
mutable struct Note <: AbstractNote
    pitch::UInt8
    velocity::UInt8
    position::UInt
    duration::UInt
    channel::UInt8

    Note(pitch, velocity, position, duration, channel = 0) =
        if channel > 0x7F
            error( "Channel must be less than 128" )
        elseif velocity > 0x7F
            error( "Velocity must be less than 128" )
        else
            new(pitch, velocity, position, duration, channel)
        end
end
Note(pitch, position; velocity = 100, duration = 960, channel = 0) =
Note(pitch, velocity, position, duration, channel)

@inline Note(n::Note) = n

"""
    DrumNote(pitch, position, duration = 960; velocity = 100)
Shorthand constructor for a [`Note`](@ref) that is always on channel `9`.
It is possible to specify pitch as a `String` (e.g. "Acoustic Snare"), which
will be converted to actual pitch using [`DRUMKEY`](@ref).
"""
DrumNote(pitch, position, duration = 960; velocity = 100, gmmap = DRUMKEY) =
Note(_drumpitch(pitch, gmmap), velocity, position, duration, 9)

_drumpitch(pitch::Real, gmmap) = pitch
_drumpitch(pitch::String, gmmap) = name_to_pitch(gmmap[pitch])

import Base.==
==(n1::AbstractNote, n2::AbstractNote) =
    n1.pitch == n2.pitch &&
    n1.duration == n2.duration &&
    n1.position == n2.position &&
    n1.channel == n2.channel &&
    n1.velocity == n2.velocity

Base.copy(n::N) where {N<:AbstractNote} =
N(n.pitch, n.velocity, n.position, n.duration, n.channel)

"""
    Notes(note_vector, tpq = 960) -> Notes
A data structure describing a collection of music notes, bundled with the ticks
per quarter note (so that the notes can be attributed rhythmic value).

`Notes` can be iterated and accessed as the given `note_vector`.
"""
struct Notes{N <: AbstractNote}
    notes::Vector{N}
    tpq::Int16
end

# Constructors for Notes:
function Notes(notes::Vector{N}, tpq::Int = 960) where {N <: AbstractNote}
    if tpq < 1
        throw(ArgumentError("Ticks per quarter note (tpq) must be >= 1"))
    end
    Notes{N}(notes, tpq)
end

Notes(; tpq = 960) = Notes{Note}(Vector{Note}[], tpq)

# Iterator Interface for notes:
Base.iterate(n::Notes, i = 1) = iterate(n.notes, i)
# Base.iterate(n::Notes, i) = iterate(n.notes, i)
Base.eltype(::Type{Notes{N}}) where {N} = N

# Indexing
Base.length(n::Notes) = length(n.notes)
Base.lastindex(n::Notes) = lastindex(n.notes)
Base.firstindex(n::Notes) = firstindex(n.notes)
Base.getindex(n::Notes, i::Int) = n.notes[i]
Base.getindex(n::Notes, r) = Notes(n.notes[r], n.tpq)
Base.view(n::Notes, r) = view(n.notes, r)

# Pushing
Base.push!(no::Notes{N}, n::N) where {N <: AbstractNote} = push!(no.notes, n)
function Base.append!(n1::Notes{N}, n2::Notes{N}) where {N}
    n1.tpq == n2.tpq || throw(ArgumentError("The Notes do not have same tpq."))
    append!(n1.notes, n2.notes)
    return n1
end

Base.copy(notes::Notes) = Notes([copy(n) for n in notes], notes.tpq)

#######################################################
# string name <-> midi pitch
#######################################################
const PITCH_TO_NAME = Dict(
0=>"C", 1=>"C♯", 2=>"D", 3=>"D♯", 4=>"E", 5=>"F", 6=>"F♯", 7=>"G", 8=>"G♯", 9=>"A",
10 =>"A♯", 11=>"B")
const NAME_TO_PITCH = Dict(
(v => k for (v, k) in zip(values(PITCH_TO_NAME), keys(PITCH_TO_NAME)))...,
"D♭"=>1, "E♭"=>3, "G♭"=>6, "A♭"=>8, "B♭"=>10)
const SHARP_TO_FLAT = Dict(
"C♯"=>"D♭", "D♯"=>"E♭", "F♯"=>"G♭", "G♯"=>"A♭", "A♯"=>"B♭"
)
"""
    pitch_to_name(pitch; flat=false) -> String
Return the name of the pitch, e.g. `F5`, `A♯3` etc. in modern notation given the
pitch value in integer. When `flat=true`, accidentals are printed as flats,
e.g. `A♯3` is printed as `B♭3`.

Reminder: middle C has pitch `60` and is displayed as `C4`.
"""
function pitch_to_name(pitch; flat::Bool=false)
    i = Int(pitch)
    notename = PITCH_TO_NAME[mod(i, 12)]
    if flat
        notename = get(SHARP_TO_FLAT, notename, notename)
    end
    octave = (i÷12)-1
    return notename*string(octave)
end

"""
    name_to_pitch(name::String) -> Int
Return the pitch value of the given note name, which can be of the form
`capital_letter*sharp*octave` where:

* `capital_letter` : from `"A"` to `"G"`.
* `sharp` : one of `"#"` `"♯"` `"b"` `"♭"` or `""`.
* `octave` : any integer (as a string), the octave number (an octave is 12 pitches).
  If not given it is assumed `"5"`.

We define E.g. `name_to_pitch("C4") === 60` (i.e. string
`"C4"`, representing the middle-C, corresponds to pitch `60`).

See http://newt.phys.unsw.edu.au/jw/notes.html
and https://en.wikipedia.org/wiki/C_(musical_note) .
"""
function name_to_pitch(name)
    pe = collect(Unicode.graphemes(name))
    pitch = NAME_TO_PITCH[pe[1]]
    x = 0
    if length(pe) >= 2
        if pe[2] == "#" || pe[2] == "♯"
            x = 1
            deleteat!(pe,2)
        elseif pe[2] == "b" || pe[2] == "♭"
            x = -1
            deleteat!(pe,2)
        end
    end
    if length(pe) == 1
        push!(pe, string(4))
    end

    octave = parse(Int, join(pe[2:end]))

    return pitch + x + 12(octave+1) # lowest possible octave is -1 but pitch starts from 0
end

"""
    pitch_to_hz(pitch::Integer, A4::Real = 440) -> hz::Real
Return the frequency value of the given midi note, optionally given the reference for middle A.
See https://en.wikipedia.org/wiki/Piano_key_frequencies
and https://librosa.org/doc/main/_modules/librosa/core/convert.html#midi_to_hz.
"""
function pitch_to_hz(pitch, A4 = 440)
    return A4 * (2^ ((pitch-69) / 12))
end

"""
    hz_to_pitch(hz::Real, A4::Real = 440) -> pitch::Int
Inverse of [`pitch_to_hz`](@ref).
"""
function hz_to_pitch(freq, A4 = 440)
    return 12 * (log2(freq) - log2(A4)) + 69
end

#######################################################
# pretty printing
#######################################################
function Base.show(io::IO, note::N) where {N<:AbstractNote}
    mprint = nameof(N)
    nn = rpad(pitch_to_name(note.pitch), 3)
    chpr = note.channel == 0 ? "" : " | channel $(note.channel)"
    velprint = rpad("vel = $(Int(note.velocity))", 9)
    print(io, "$(mprint) $nn | $velprint | "*
    "pos = $(Int(note.position)), "*
    "dur = $(Int(note.duration))"*chpr)
end

function Base.show(io::IO, ::MIME"text/plain", notes::Notes{N}) where {N}
    mprint = nameof(N)
    print(io, "$(length(notes)) $(mprint)s with tpq=$(notes.tpq)")
    _notevectorprint(io, notes)
end
function Base.show(io::IO, notes::Vector{N}) where {N<:AbstractNote}
    mprint = nameof(N)
    print(io, "$(length(notes))-element Vector{$mprint}")
    _notevectorprint(io, notes)
end

function _notevectorprint(io, notes)
    s = 7
    if length(notes) > 2s
        for note in (@view notes[1:s])
            print(io, "\n", " ", note)
        end
        print(io, "\n", "  ⋮")
        for note in (@view notes[end-s+1:end])
            print(io, "\n", " ", note)
        end
    else
        for note in notes
            print(io, "\n", " ", note)
        end
    end
end

function Base.show(io::IO, notes::Notes{N}) where {N}
    mprint = nameof(N)
    print(io, "Notes{$(mprint)} with $(length(notes)) notes")
end


text:
"""
Truncate a cutoff vector to the probability simplex. Useful because some methods
produce negative cutoffs while iterating.
"""
function cutbox(cut)
    return max.(0, min.(1, cut))
end

"""
    assignmentfrompreflists(students_inv, students_dist, cutoffs;
                        return_demands=false)

Return assignment associated with given cutoffs and, if `return_demands=true`,
the demands. For demands only, `demandfrompreflists()` is faster. Ignores
capacity constraints.
"""
function assignmentfrompreflists(students_inv     ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                             students_dist    ::AbstractArray{<:AbstractFloat, 1},
                             cutoffs          ::AbstractArray{<:AbstractFloat, 1};
                             return_demands   ::Bool=false)
    (m, n) = size(students_inv)
    @assert size(cutoffs) == (m, )        "Dim mismatch between students_inv and cutoffs"
    @assert size(students_dist) == (n, )  "Dim mismatch between students_inv and students_dist"

    cutoffs = cutbox(cutoffs)

    assn = zeros(m + 1, n)
    unassigned = copy(students_dist)

    for s in 1:n, c in 1:m
        got_in = unassigned[s] * (1 - cutoffs[students_inv[c, s]])
        assn[students_inv[c, s], s] += got_in
        unassigned[s] -= got_in
    end

    assn[end, :] = unassigned

    if return_demands    # For testing purposes
        demands = sum(assn[1:end - 1, :], dims=2)
        return assn, demands
    else
        return assn
    end
end

"""
    demandfrompreflists(students, students_dist, cutoffs)

Return demand for each school given a set of cutoffs and ignoring capacity, using
given preference lists and assuming iid scores at each school.

Satifies WGS and score independence, so equilibrium can be computed using
`DA_nonatomic_lite()`.
"""
function demandfrompreflists(students      ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                           students_dist ::AbstractArray{<:AbstractFloat, 1},
                           cutoffs       ::AbstractArray{<:AbstractFloat, 1},
                           )::AbstractArray{<:AbstractFloat, 1}
    (m, n) = size(students)
    @assert size(cutoffs) == (m, )        "Dim mismatch between students and cutoffs"
    @assert size(students_dist) == (n, )  "Dim mismatch between students and students_dist"

    cutoffs = cutbox(cutoffs)

    demands = [(1 - cutoffs[c]) * sum(students_dist[s] *
               prod(cutoffs[students[:, s] .< students[c, s]]) for s in 1:n)
               for c in 1:m]

    return demands
end

"""
    demandfromMNL_iid(qualities, cutoffs)

Return demand for each school given a set of cutoffs and ignoring capacity, using
multinomial logit choice model and assuming scores are iid uniform (equivalent
to MTB when schools have no preferences).

Satifies WGS and score independence, so equilibrium can be computed using
`DA_nonatomic_lite()`.
"""
function demandfromMNL_iid(qualities   ::AbstractArray{<:AbstractFloat, 1},
                         cutoffs     ::AbstractArray{<:AbstractFloat, 1},
                         )::AbstractArray{<:AbstractFloat, 1}
    (m, ) = size(qualities)
    @assert (m, ) == size(cutoffs) "Dim mismatch"
    cutoffs = cutbox(cutoffs)

    demands = zeros(m)

    γ = exp.(qualities)

    for c in 1:m
        C_minus_c = setdiff(1:m, c)
        for C♯ ∈ powerset(C_minus_c)
            demands[c] += prod(e in C♯ ? 1 - cutoffs[e] : cutoffs[e] for e in C_minus_c) /
                          (γ[c] + sum(AbstractFloat[γ[d] for d in C♯]))
        end
        demands[c] *= (1 - cutoffs[c]) * γ[c]
    end

    return demands
end

"""
    demandfromMNL_singlescore(qualities, cutoffs)

Return demand for each school given a set of cutoffs and ignoring capacity, using
multinomial logit choice model and assuming all schools use a single score (equivalent
to STB when schools have no preferences).

Satisfies WGS but not score independence, so equilibrium must be computed using
`nonatomic_tatonnement()`.
"""
function demandfromMNL_singlescore(qualities   ::AbstractArray{<:AbstractFloat, 1},
                             cutoffs   ::AbstractArray{<:AbstractFloat, 1},
                             )::AbstractArray{<:AbstractFloat, 1}
    (m, ) = size(qualities)
    @assert (m, ) == size(cutoffs) "Dim mismatch"

    cutoffs = cutbox(cutoffs)

    demands = zeros(m)

    sort_order = sortperm(cutoffs)
    cutoffs[sort_order]

    γ = exp.(qualities)
    demands = zeros(m)

    consideration_set_probabilities = diff([cutoffs[sort_order]; 1])

    for c in 1:m, d in c:m     # For each score threshold
        demands[sort_order[c]] += consideration_set_probabilities[d] *
                                  γ[sort_order[c]] / sum(γ[sort_order[1:d]])
    end

    return demands
end

# Currently unused.
"""
    numericaladmitrate(B, cutoffs; n_points=1000)

Numerically compute the admissions rate for each school when school's
blends are the rows of B and scores are iid uniform. That is, compute
the probability that `B[c, :] * x >= cutoffs[c]`, where x is iid uniform.
Uses a meshgrid iterator with (about) `n_points` of evaluation.
"""
function numericaladmitrate(B         ::AbstractArray{<:AbstractFloat, 2},
                            cutoffs   ::AbstractArray{<:AbstractFloat, 1};
                            n_points  ::Int=1000,
                            )::AbstractArray{<:AbstractFloat, 1}

    (m, n_tests) = size(B)
    density = Int(ceil(n_points ^ (1 / n_tests)))

    ran = range(0, stop=1, length=density)
    admits = zeros(Int, m)

    for x in Iterators.product(repeat([ran], n_tests)...), j in 1:m
        admits[j] += sum(B[j, k] * x[k] for k in 1:n_tests) >= cutoffs[j]
    end

    return admits ./ (density ^ n_tests)
end

"""
    demandfrommixedMNL_tscores(qualities, profile_dist, blends, cutoffs;
                        n_points=10000, verbose=false, montecarlo=false)

Return demand for each school given a set of cutoffs and ignoring capacity, using
mixed multinomial logit choice model with `p` student profiles and `t` test scores
shared among schools.

`qualities` is a `C` by `p` matrix, where `C` is the number of schools. Each column
of `qualities` gives the quality vector with respect to a student profile. This
generalizes `demandfromMNL_iid()`.

If an applicant's test score vector is `x`, then her composite score at school
`c` is `blend[c, :] * x`, where `blend[c, :]` is a school-specific convex
combination of the `t` test scores, which we assume are iid uniform. This
generalizes `demandfromMNL_singlescore()`.

`n_points` is the approximate number of test points to use if `montecarlo` evaluation
is selected. Which method is more effecient depends on the problem dimensions, but
in general, for problems with many schools or tests, consider MC.

Satisfies WGS but not score independence, so equilibrium must be computed using
`nonatomic_tatonnement()`.
"""
function demandfrommixedMNL_tscores(qualities      ::AbstractArray{<:AbstractFloat, 2},
                             profile_dist   ::AbstractArray{<:AbstractFloat, 1},
                             blends         ::AbstractArray{<:AbstractFloat, 2},
                             cutoffs        ::AbstractArray{<:AbstractFloat, 1};
                             montecarlo     ::Bool=false,
                             n_points       ::Int=10000,
                             verbose        ::Bool=false,
                             )::AbstractArray{<:AbstractFloat, 1}
    (m, p) = size(qualities)

    @assert (p, ) == size(profile_dist)     "Dim mismatch between qualities and profile_dist"
    @assert m == size(blends)[1]            "Dim mismatch between qualities and blends"
    @assert (m, ) == size(cutoffs)          "Dim mismatch between qualities and cutoffs"
    @assert sum(profile_dist) ≈ 1           "profile_dist doesn't sum to 1"
    @assert all(sum(blends, dims=2) .≈ 1)   "rows of blends don't sum to 1"

    cutoffs = cutbox(cutoffs)

    (m, t) = size(blends)

    demands = zeros(m)

    γ = exp.(qualities)

    if montecarlo
        density = Int(ceil(n_points ^ (1 / t)))
        sample_size = density ^ t   # approx. n_points

        # To use a grid instead of Monte Carlo. Tends to be less accurate because
        # of bias toward the edges of the cube.
        # ran = range(0, stop=1, length=density)

        # for x in Iterators.product(repeat([ran], t)...)
        for _ in 1:sample_size
            x = rand(t)
            verbose ? println("Student with scores:   ", x) : nothing
            # Consideration set
            C♯ = blends * [x...] .>= cutoffs
            verbose ? println("Has consideration set: ", C♯) : nothing

            if true in C♯
                for s in 1:p
                    coef = profile_dist[s] / sum(γ[d, s] * C♯[d] for d in 1:m)
                    verbose ? println("  For profile $s, coef is $coef") : nothing

                    for c in 1:m
                        if C♯[c]
                            demands[c] += γ[c, s] * coef
                            verbose ? println("    Go to school $c wp $(γ[c, s]) * coef = ",
                                              γ[c, s] * coef) : nothing
                        end
                    end
                end
            end
        end

        return demands ./ sample_size

    else        # Exact method
        bounds = vcat([HalfSpace(-col, 0.) for col in eachcol(I(t))],
                      [HalfSpace(1col, 1.) for col in eachcol(I(t))]) # Need 1 so I allocates

        for C♯ in powerset(1:m)
            if !isempty(C♯)
                hspaces = HalfSpace{Float64,Array{Float64,1}}[]

                # Probability of having this choice set is the volume of
                # this m-dimensional polyhedron.
                for c in C♯
                    if cutoffs[c] ≥ 1
                        @goto volume_is_zero
                    elseif cutoffs[c] > 0
                        push!(hspaces, HalfSpace(-blends[c, :], -cutoffs[c]))
                    end
                end

                for c in setdiff(1:m, C♯)
                    if cutoffs[c] ≤ 0
                        @goto volume_is_zero
                    elseif cutoffs[c] < 1
                        push!(hspaces, HalfSpace( blends[c, :],  cutoffs[c]))
                    end
                end

                poly = polyhedron(hrep(vcat(bounds, hspaces)))
                vol = volume(poly)

                for s in 1:p
                    mult = vol * profile_dist[s] / sum(γ[C♯, s])

                    for c in C♯
                        demands[c] += mult * γ[c, s]
                    end
                end

                @label volume_is_zero
            end
        end

        return demands
    end
end

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
# https://adventofcode.com/2020/day/10

# Each of your joltage adapters is rated for a specific output joltage (your puzzle input).
# Any given adapter can take an input 1, 2, or 3 jolts lower than its rating and still produce
# its rated output joltage.

# In addition, your device has a built-in joltage adapter rated for 3 jolts higher than the
# highest-rated adapter in your bag. (If your adapter list were 3, 9, and 6, your device's
# built-in adapter would be rated for 12 jolts.)

# Treat the charging outlet near your seat as having an effective joltage rating of 0.

function read_array(string::AbstractString)::Vector{Int}
    return parse.(Int, split(string, '\n', keepempty=false))
end

example1 = read_array("
16
10
15
5
1
11
7
19
6
12
4
")

example2 = read_array("
28
33
18
42
31
14
46
20
48
47
24
23
49
45
19
38
39
11
1
32
25
35
8
17
7
9
4
2
34
10
3
")

# What is the number of 1-jolt differences multiplied by the number of 3-jolt differences?

function part1(numbers)
    sorted = sort(numbers)
    sorted = [0; sorted; sorted[end] + 3]
    differences = diff(sorted)
    @assert minimum(differences) == 1
    @assert maximum(differences) == 3

    return sum(differences .== 1) * sum(differences .== 3)
end

@assert part1(example2) == 220

# What is the total number of distinct ways you can arrange the adapters to connect
# the charging outlet to your device?

"""
Cache for holding the partial solution, it maps: node => solutions_count.
"""
Cache = Dict{Int,Int}

"""
Use memoized shortcuts when traversing the tree. On the way, update the cache
in place. Returns the number of paths from the beggining of the sorted list
of nodes.
"""
function count_solutions!(adapters, cache::Cache=Cache())
    socket = adapters[1]
    if socket in keys(cache)
        return cache[socket]
    end

    solutions_found = 0
    for i = 2:4
        if i > length(adapters)
            break
        end
        if adapters[i] <= (socket + 3)
            if i == length(adapters)
                solutions_found += 1
            else
                # updates cache in place, will be available for next iteration
                solutions_found += count_solutions!(adapters[i:end], cache)
            end
        end
    end
    cache[socket] = solutions_found
    return solutions_found
end

@assert count_solutions!([0; sort(example1); 22], Cache()) == 8
@assert count_solutions!([0; sort(example2); 52], Cache()) == 19208

"""
Traverse the tree starting from the back, to build the memoization
cache.
"""
function init_cache(adapters, step_size=0.05)
    n = length(adapters)
    step = Int(max(1, round(n * step_size)))
    cache = Cache()

    pos = n
    while pos >= n / 3
        pos -= step
        count_solutions!(adapters[pos:end], cache)
    end
    return cache
end

function part2(numbers)
    sorted = sort(numbers)
    full_sequence = [0; sorted; sorted[end] + 3]
    cache = init_cache(full_sequence)
    return count_solutions!(full_sequence, cache)
end

@assert part2(example1) == 8
@assert part2(example2) == 19208

test = read_array(read("data/day-10.txt", String))
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

part1(test) == 2738
part2(test) == 74049191673856

println()
print("Example 1: ")
@time @assert part2(example1) == 8
print("Example 2: ")
@time @assert part2(example2) == 19208
print("Part 2:    ")
@time @assert result2 == 74049191673856

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
# These tests are closely based on the server tests in HTTP.jl:
# https://github.com/JuliaWeb/HTTP.jl/blob/master/test/server.jl

using HTTP
import HTTP.IOExtras.bytes

# Function to make a simple request that can be used in tests
function sample_request()
    body = """{"x": 5, "y": 3}"""
    encoded_body = bytes(body)
    return HTTP.Request("GET", "", [], encoded_body)
end

include("utilities-tests.jl")
include("middleware-tests.jl")
include("responders-tests.jl")
include("integration-tests.jl")
include("dockerfile-tests.jl")


text:
"""
Is there zero intersection between two simplices?
"""
function zerointersection(S1::AbstractArray{Float64, 2}, S2::AbstractArray{Float64, 2})
    if toofaraparattointersect(S1, S2) == true
      return(true)
    end

    # If some simplex is singular, they cannot intersect
    if somesimplexissingular(S1, S2) == true;
      return(true)
    end
end

"""
Checks whether two simplices are too far apart to intersect.
"""
function toofaraparattointersect(S1::AbstractArray{Float64, 2}, S2::AbstractArray{Float64, 2})

      r1 = Geometry.squaredradius(S1)
      r2 = Geometry.squaredradius(S2)
      c1 = Geometry.centroid(S1)
      c2 = Geometry.centroid(S2)

      dist = norm(c1 - c2) - (sqrt(r1) + sqrt(r2))
      if dist >= 0
        return(true)
      else
        return(false)
      end
end

"""
Checks whether at least one of the two simplices defined by the vertices
in  'X' and 'Y' are singular. X and Y are matrices of dimension nx(n+1).
"""
function somesimplexissingular(X::Array{Float64, 2}, Y::Array{Float64, 2})
    n = size(X)[1]
    if (vcat(ones(1, n+1), X) == vcat(ones(1, n + 1), Y))
      print("At least one of the simplices is singular!")
      return true
    else
      return false
    end
end

"""
Make sure simplices are of the same dimension
"""
function validatesimplices(S1::Array{Float64, 2}, S2::Array{Float64, 2})
    assert(size(S1)[1] == size(S1)[1])
end


text:
using NumCME
using Test
using BenchmarkTools
using SparseArrays

𝕊 = [[1] [-1]]
λ = 10.0
γ = 5.0
a1 = propensity() do x, p
    p[1]
end
a2 = propensity() do x, p
    p[2]
end
sensmodel = CmeModelWithSensitivity(CmeModel(𝕊, [a1, a2], [λ, γ]))
test_space = StateSpaceSparse(𝕊, [[i] for i in 1:100])
sensmat = ForwardSensFspMatrixSparse{Float64}(sensmodel, test_space)
v = ones(3*get_rowcount(sensmat.fspmatrix))
out = similar(v)
matvec!(out, 0.0, sensmat, v)
@test sum(out) ≈ 0.0


text:
function solveKnapHeur(data::KnapData)
    perm = sortperm(data.items, by = x -> x.profit / x.weight)

    solution = Int64[]
    best_val = 0
    capacity = data.capacity
    for index in perm
        if data.items[index].weight <= capacity
            push!(solution, index)
            best_val += data.items[index].profit
            capacity -= data.items[index].weight
        end
    end

    return best_val, solution
end


text:
"Normalize( data::SingleCellExpr, reads::Int64 )

Normalize Single Cell Data the hard way - resampling

This will add/update the samples column ResampleN where it will store the
final re-scale value used to obtain the total amount of requested reads in the sample.
Lets see - possibly this value will be of use later on.

This function will add -1 values for all sample/gene touples that had a value in the raw data,
but were set to 0 due to the resampling.
"
function Normalize( data::SingleCellExpr, reads::Int64 )

  if (typeof(data.raw) != Int64)
    println( "The data has been normalized before - reverting to raw data")
    data.data = data.raw
    data.raw = 0
    data.zscored=0
  end

  println(string("the data will be re-sampled to ", reads, " reads"))

  f = function(x)
    sum(nonzeros(x))
  end
  csum = mapslices( f, data.data, 1)

  nm = length( nonzeros(data.data) )
  pos = 0
  #define max result values
  v = zeros(nm)
  r = zeros(nm)
  c = zeros(nm)
  ns = zeros(size(data.data)[2])
  for i = 1:size(data.data)[2]
    n = reads
    p = data.data[1:size(data.data)[1],i] / csum[i]
    if  sum(map(round, p*reads)) > reads
       while  sum(map(round, p*n)) > reads
         n -= 1
       end
    elseif sum(map(round, p*reads)) < reads
       while  sum(map(round, p*n)) < reads
         n += 1
       end
    end
    val = map(round, p*n)
    ns[i] = n
    ## now I likely have lost values
    ## which I want to set to -1
    t = countmap(vcat(p.nzind, val.nzind))
    bad = collect(keys(t))[find(collect(values(t)) .== 1)]
    ## collect the Int64 > 0 values
    ok = (pos+1):(pos+length(val.nzind)+length(bad))
    v[ok] = vcat(collect(val.nzval), fill( -1, length(bad) ) )
    r[ok] = vcat(val.nzind, bad)
    c[ok] = fill(i, length(val.nzind) + length(bad))
    pos = pos+length(val.nzind)+ length(bad)
  end
  data.raw = data.data
  data.data = sparse( r[1:pos], c[1:pos] , v[1:pos] )
  data.samples[:ResampleN] = ns
  #= build the new sparse matrix =#
return (data)
end


text:
using FooBar
using Documenter

DocMeta.setdocmeta!(FooBar, :DocTestSetup, :(using FooBar); recursive=true)

makedocs(;
    modules=[FooBar],
    authors="hkauhanen <nenahuak@gmail.com> and contributors",
    repo="https://github.com/hkauhanen/FooBar.jl/blob/{commit}{path}#{line}",
    sitename="FooBar.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://hkauhanen.github.io/FooBar.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/hkauhanen/FooBar.jl",
    devbranch="main",
)


text:
function get_parameter(obj::ParametersController, arg0::JString)
    return jcall(obj, "getParameter", jdouble, (JString,), arg0)
end

function set_parameter(obj::ParametersController, arg0::JString, arg1::jdouble)
    return jcall(obj, "setParameter", void, (JString, jdouble), arg0, arg1)
end


text:
using ConstructiveSolidGeometry
using Base.Test

# Unit tests for ray <-> plane intersection

# Parallel
test_ray = Ray(Coord(0.0, 0.0, 0.0), Coord(0.0, 1.0, 0.0) )
test_plane = Plane(Coord(1.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0))
@test raytrace(test_ray, test_plane)[1] == false

# Hit
test_ray = Ray(Coord(0.0, 0.0, 0.0), unitize(Coord(1.0, 1.0, 0.0)) )
test_plane = Plane(Coord(1.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0))
@test raytrace(test_ray, test_plane) == (true, sqrt(2))

# Miss
test_ray = Ray(Coord(0.0, 0.0, 0.0), unitize(Coord(-1.0, -1.0, 0.0)) )
test_plane = Plane(Coord(1.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0))
@test raytrace(test_ray, test_plane)[1] == false

# Inside Plane
test_ray = Ray(Coord(0.0, 0.0, 0.0), Coord(0.0, 1.0, 0.0) )
test_plane = Plane(Coord(0.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0))
@test raytrace(test_ray, test_plane)[1] == true

# Unit tests for ray <-> sphere intersection

# Hit from outside
test_ray = Ray(Coord(-2.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0) )
test_sphere = Sphere(Coord(0.0, 0.0, 0.0), 1.0)
@test raytrace(test_ray, test_sphere) == (true, 1.0)

# Hit from inside (backwards is closer)
test_ray = Ray(Coord(-0.9, 0.0, 0.0), Coord(1.0, 0.0, 0.0) )
test_sphere = Sphere(Coord(0.0, 0.0, 0.0), 1.0)
@test raytrace(test_ray, test_sphere) == (true, 1.9)

# Hit from inside (backwards is farther)
test_ray = Ray(Coord(-0.5, 0.0, 0.0), Coord(-1.0, 0.0, 0.0) )
test_sphere = Sphere(Coord(0.0, 0.0, 0.0), 1.0)
@test raytrace(test_ray, test_sphere) == (true, 0.5)

# Miss from outside
test_ray = Ray(Coord(-2.0, 0.0, 0.0), Coord(-1.0, 0.0, 0.0) )
test_sphere = Sphere(Coord(0.0, 0.0, 0.0), 1.0)
@test raytrace(test_ray, test_sphere)[1] == false

# Ray origin is on sphere
test_ray = Ray(Coord(-1.0, 0.0, 0.0), Coord(0.0, 1.0, 0.0) )
test_sphere = Sphere(Coord(0.0, 0.0, 0.0), 1.0)
@test raytrace(test_ray, test_sphere) == (true, 0.0)

# Unit tests for ray <-> Cone intersection

# Hit from inside
test_ray = Ray(Coord(0.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0) )
test_cone = Cone(Coord(0.0, 0.0, 1.0), Coord(0.0, 0.0, -1.0), pi/4.0 )
@test raytrace(test_ray, test_cone)[1] == true
@test raytrace(test_ray, test_cone)[2] ≈ 1.0

# Hit from inside
test_ray = Ray(Coord(-0.5, 0.0, 0.0), Coord(1.0, 0.0, 0.0) )
test_cone = Cone(Coord(0.0, 0.0, 1.0), Coord(0.0, 0.0, -1.0), pi/4.0 )
@test raytrace(test_ray, test_cone)[1] == true
@test raytrace(test_ray, test_cone)[2] ≈ 1.5

# Hit from outside
test_ray = Ray(Coord(-5.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0) )
test_cone = Cone(Coord(0.0, 0.0, 1.0), Coord(0.0, 0.0, -1.0), pi/4.0 )
@test raytrace(test_ray, test_cone)[1] == true
@test raytrace(test_ray, test_cone)[2] ≈ 4.0

# Hit at angle from outside
test_ray = Ray(Coord(-5.0, 1.0, 1.0), unitize(Coord(1.0, 0.0, -1.0)) )
test_cone = Cone(Coord(0.0, 0.0, 1.0), unitize(Coord(0.0, 0.1, -1.0)), pi/4.0 )
@test raytrace(test_ray, test_cone)[1] == true
@test raytrace(test_ray, test_cone)[2] ≈ 3.5511721769144815

# Hit at 45 at angle from outside, for rotated cylinder
test_ray = Ray(Coord(1.0, 1.0, 2.0), unitize(Coord(1.0, 1.0, 0.0)) )
test_cone = Cone(Coord(1.0, 1.0, 1.0), unitize(Coord(1.0, 1.0, 0.0)), pi/4.0 )
@test raytrace(test_ray, test_cone)[1] == true
@test raytrace(test_ray, test_cone)[2] ≈ 1.0

# Hit of shadow of cone
test_ray = Ray(Coord(0.0, 0.0, 2.0), Coord(1.0, 0.0, 0.0) )
test_cone = Cone(Coord(0.0, 0.0, 1.0), Coord(0.0, 0.0, -1.0), pi/4.0 )
@test raytrace(test_ray, test_cone)[1] == false

# Valid hit after shadow of cone
test_ray = Ray(Coord(0.5, 0.0, 2.0), Coord(0.0, 0.0, -1.0) )
test_cone = Cone(Coord(0.0, 0.0, 1.0), Coord(0.0, 0.0, -1.0), pi/4.0 )
@test raytrace(test_ray, test_cone)[1] == true
@test raytrace(test_ray, test_cone)[2] ≈ 1.5

# Valid hit before shadow of cone
test_ray = Ray(Coord(0.5, 0.0, 0.0), Coord(0.0, 0.0, 1.0) )
test_cone = Cone(Coord(0.0, 0.0, 1.0), Coord(0.0, 0.0, -1.0), pi/4.0 )
@test raytrace(test_ray, test_cone)[1] == true
@test raytrace(test_ray, test_cone)[2] ≈ 0.5

# Unit tests for ray <-> Infinite Cylinder intersection
# Z Cylinder Tests

test_cylinder = InfCylinder(Coord(-2.0, 0.0, 0.0), Coord(0.0, 0.0, 1.0), 1.0)

# Hit from outside
test_ray = Ray(Coord(0, 0.0, 0.0), Coord(-1.0, 0.0, 0.0) )
@test raytrace(test_ray, test_cylinder) == (true, 1.0)

# Hit from inside (backwards is closer)
test_ray = Ray(Coord(-1.5, 0.0, 0.0), Coord(-1.0, 0.0, 0.0) )
@test raytrace(test_ray, test_cylinder) == (true, 1.5)

# Hit from inside (forwards is closer)
test_ray = Ray(Coord(-2.5, 0.0, 0.0), Coord(-1.0, 0.0, 0.0) )
@test raytrace(test_ray, test_cylinder) == (true, 0.5)

# Miss from outside
test_ray = Ray(Coord(0.0, 2.5, 0.0), Coord(-1.0, 0.0, 0.0) )
@test raytrace(test_ray, test_cylinder)[1] == false

# Tilted Cylinder tests
test_cylinder = InfCylinder(Coord(0.0, 0.0, 0.0), unitize(Coord(1.0, 0.0, 1.0)), 1.0)

# Hit from outside
test_ray = Ray(Coord(-2.0, 0.0, 2.0), unitize(Coord(1.0, 0.0, -1.0)) )
@test raytrace(test_ray, test_cylinder)[1] == true
@test raytrace(test_ray, test_cylinder)[2] ≈ sqrt(8)-1.0

# Hit from inside (backwards is closer)
test_ray = Ray(Coord(-.5, 0.0, 0.5), unitize(Coord(1.0, 0.0, -1.0)) )
@test raytrace(test_ray, test_cylinder)[1] == true
@test raytrace(test_ray, test_cylinder)[2] ≈ sqrt(0.5^2+0.5^2) + 1

# Hit from inside (forwards is closer)
test_ray = Ray(Coord(0.5, 0.0, -.5), unitize(Coord(1.0, 0.0, -1.0)) )
@test raytrace(test_ray, test_cylinder)[2] ≈ 1 - sqrt(0.5^2+0.5^2)

# Miss from outside
test_ray = Ray(Coord(0.0, 2.5, 0.0), Coord(-1.0, 0.0, 0.0) )
@test raytrace(test_ray, test_cylinder)[1] == false

# Reflection testing
test_plane = Plane(Coord(1, 0, 0), Coord(-1, 0, 0), "reflective")

# Straight on
test_ray = Ray(Coord(0, 0, 0), Coord(1, 0, 0))
@test reflect(test_ray, test_plane).direction.x ≈ -1.0
@test reflect(test_ray, test_plane).direction.y ≈ 0.0
@test reflect(test_ray, test_plane).direction.z ≈ 0.0

# 45 degree angle
test_ray = Ray(Coord(0, 0, 0), unitize(Coord(1, 1, 0)))
@test reflect(test_ray, test_plane).direction.x ≈ -1.0/sqrt(2.0)
@test reflect(test_ray, test_plane).direction.y ≈ 1.0/sqrt(2.0)
@test reflect(test_ray, test_plane).direction.z ≈ 0.0

# 45 degree plane
test_plane = Plane(Coord(1, 1, 0), unitize(Coord(-1, -1, 0)), "reflective")

# Straight on
test_ray = Ray(Coord(0, 0, 0), unitize(Coord(1, 1, 0)))
@test reflect(test_ray, test_plane).direction.x ≈ -1.0/sqrt(2.0)
@test reflect(test_ray, test_plane).direction.y ≈ -1.0/sqrt(2.0)
@test reflect(test_ray, test_plane).direction.z ≈ 0.0

# At an angle
test_ray = Ray(Coord(0, 0, 0), Coord(1, 0, 0))
@test reflect(test_ray, test_plane).direction.x <= 1e-10
@test reflect(test_ray, test_plane).direction.y ≈ -1.0
@test reflect(test_ray, test_plane).direction.z ≈ 0.0

# Test plane halfspace

# Point is below a z-plane with normal pointing up
test_coord = Coord(0.0, 0.0, 0.0)
test_plane = Plane(Coord(0.0, 0.0, 1.0), unitize(Coord(0.0, 0.0, 1.0)))
@test halfspace(test_coord, test_plane) == -1

# Point is below a z-plane with normal pointing down
test_coord = Coord(0.0, 0.0, 0.0)
test_plane = Plane(Coord(0.0, 0.0, 1.0), unitize(Coord(0.0, 0.0, -1.0)))
@test halfspace(test_coord, test_plane) == 1

# Point is above a z-plane with normal pointing up
test_coord = Coord(0.0, 0.0, 2.0)
test_plane = Plane(Coord(0.0, 0.0, 1.0), unitize(Coord(0.0, 0.0, 1.0)))
@test halfspace(test_coord, test_plane) == 1

# Point is above a z-plane with normal pointing down
test_coord = Coord(0.0, 0.0, 2.0)
test_plane = Plane(Coord(0.0, 0.0, 1.0), unitize(Coord(0.0, 0.0, -1.0)))
@test halfspace(test_coord, test_plane) == -1

# Point is on a z-plane
test_coord = Coord(0.0, 0.0, 1.0)
test_plane = Plane(Coord(0.0, 0.0, 1.0), unitize(Coord(0.0, 0.0, -1.0)))
@test halfspace(test_coord, test_plane) == -1

# Test sphere halfspace

# Point is inside a sphere
test_coord = Coord(-1.9, -1.9, -2.5)
test_sphere = Sphere(Coord(-2.0, -2.0, -2.0), 1.0)
@test halfspace(test_coord, test_sphere) == -1

# Point is outside a sphere
test_coord = Coord(-1.9, -1.9, -2.5)
test_sphere = Sphere(Coord(-2.0, -2.0, -2.0), 0.1)
@test halfspace(test_coord, test_sphere) == 1

# Point on a sphere
test_coord = Coord(1.0, 0.0, 0.0)
test_sphere = Sphere(Coord(0.0, 0.0, 0.0), 1.0)
@test halfspace(test_coord, test_sphere) == -1

# Test Cone halfspace

# Inside cone
test_coord = Coord(0.0, 0.0, 0.0)
test_cone = Cone(Coord(0.0, 0.0, 1.0), Coord(0.0, 0.0, -1.0),  pi/8.0)
@test halfspace(test_coord, test_cone) == 1

# outside cone
test_coord = Coord(5.0, 0.0, 0.0)
@test halfspace(test_coord, test_cone) == -1

# outside cone (but inside shadow)
#test_coord = Coord(0.0, 0.0, 2.0)
@test halfspace(test_coord, test_cone) == -1

# Inside rotated cylinder
test_coord = Coord(1.5, 1.5, 1.1)
test_cone = Cone(Coord(1.0, 1.0, 1.0), unitize(Coord(1.0, 1.0, 0.0)), pi/4.0 )
@test halfspace(test_coord, test_cone) == 1

# Outside rotated cylinder
test_coord = Coord(1.5, 1.5, 2.1)
test_cone = Cone(Coord(1.0, 1.0, 1.0), unitize(Coord(1.0, 1.0, 0.0)), pi/4.0 )
@test halfspace(test_coord, test_cone) == -1

# Test InfCylinder Halfspace

# Inside cylinder
test_coord = Coord(2.0, 2.0, 2.0)
test_cylinder = InfCylinder(Coord(1.9, 1.9, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 1.0)
@test halfspace(test_coord, test_cylinder) == -1

# Outside cylinder
test_coord = Coord(-2.0, 2.0, 2.0)
test_cylinder = InfCylinder(Coord(1.9, 1.9, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 1.0)
@test halfspace(test_coord, test_cylinder) == 1

# On cylinder
test_coord = Coord(2.9, 1.9, 10.0)
test_cylinder = InfCylinder(Coord(1.9, 1.9, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 1.0)
@test halfspace(test_coord, test_cylinder) == -1

# Bounding Box for the Problem
bounding_box = Box(Coord(-.63, -.63, -150), Coord(.63, .63, 150))

# Containers to store cells as we make them
cells = Array{Cell}(0)

# Define all the surfaces we will use up front
# Box
top = Plane(Coord(0.0, 0.0, 150.0), unitize(Coord(0.0, 0.0, 1.0)),   "reflective")
bot = Plane(Coord(0.0, 0.0, -150.0), unitize(Coord(0.0, 0.0, -1.0)), "reflective")
left = Plane(Coord(-.63, 0.0, 0.0), unitize(Coord(-1.0, 0.0, 0.0)),  "reflective")
right = Plane(Coord(0.63, 0.0, 0.0), unitize(Coord(1.0, 0.0, 0.0)),  "reflective")
up = Plane(Coord(0.0, 0.63, 0.0), unitize(Coord(0.0, 1.0, 0.0)),     "reflective")
down = Plane(Coord(0.0, -0.63, 0.0), unitize(Coord(0.0, -1.0, 0.0)), "reflective")
clad_outer = InfCylinder(Coord(0.0, 0.0, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 0.4750)
clad_inner = InfCylinder(Coord(0.0, 0.0, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 0.4180)
fuel = InfCylinder(Coord(0.0, 0.0, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 0.4096)

# Make the Water Region
regions = Array{Region}(0)

push!(regions, Region(top, -1))
push!(regions, Region(bot, -1))
push!(regions, Region(left, -1))
push!(regions, Region(right, -1))
push!(regions, Region(up, -1))
push!(regions, Region(down, -1))
push!(regions, Region(clad_outer, 1))

ex = :(1 ^ 2 ^ 3 ^ 4 ^ 5 ^ 6 ^ 7)

push!(cells, Cell(regions, ex))

# Make the Cladding
regions = Array{Region}(0)

push!(regions, Region(top, -1))
push!(regions, Region(bot, -1))
push!(regions, Region(clad_outer, -1))
push!(regions, Region(clad_inner, 1))

ex = :(1 ^ 2 ^ 3 ^ 4)

push!(cells, Cell(regions, ex))

# Make the Gap
regions = Array{Region}(0)

push!(regions, Region(top, -1))
push!(regions, Region(bot, -1))
push!(regions, Region(clad_inner, -1))
push!(regions, Region(fuel, 1))

ex = :(1 ^ 2 ^ 3 ^ 4)

push!(cells, Cell(regions, ex))

# Make the Fuel
regions = Array{Region}(0)

push!(regions, Region(top, -1))
push!(regions, Region(bot, -1))
push!(regions, Region(fuel, -1))

ex = :(1 ^ 2 ^ 3)

push!(cells, Cell(regions, ex))

geometry = Geometry(cells, bounding_box)

# Some tests for our test geometry
@test find_cell_id(Coord(0.0,0.1,0), geometry) == 4
@test find_cell_id(Coord(0.0,0.41,0), geometry) == 3
@test find_cell_id(Coord(0.0,0.42,0), geometry) == 2
@test find_cell_id(Coord(0.0,0.55,0), geometry) == 1

# Some tests for weirder definitions
# Bounding Box for the Problem
bounding_box = Box(Coord(-.63, -.63, -150), Coord(.63, .63, 150))

# Containers to store cells as we make them
cells = Array{Cell}(0)

# Define all the surfaces we will use up front
# Box
top = Plane(Coord(0.0, 0.0, 150.0), unitize(Coord(0.0, 0.0, 1.0)), "reflective")
bot = Plane(Coord(0.0, 0.0, -150.0), unitize(Coord(0.0, 0.0, -1.0)), "reflective")
left = Plane(Coord(-.63, 0.0, 0.0), unitize(Coord(-1.0, 0.0, 0.0)), "reflective")
right = Plane(Coord(0.63, 0.0, 0.0), unitize(Coord(1.0, 0.0, 0.0)), "reflective")
up = Plane(Coord(0.0, 0.63, 0.0), unitize(Coord(0.0, 1.0, 0.0)), "reflective")
down = Plane(Coord(0.0, -0.63, 0.0), unitize(Coord(0.0, -1.0, 0.0)), "reflective")
left_cyl = InfCylinder(Coord(-0.175, 0.0, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 0.3 )
right_cyl = InfCylinder(Coord(0.1, 0.0, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 0.2)

# Make the outer Region
regions = Array{Region}(0)

push!(regions, Region(top, -1))
push!(regions, Region(bot, -1))
push!(regions, Region(left, -1))
push!(regions, Region(right, -1))
push!(regions, Region(up, -1))
push!(regions, Region(down, -1))
push!(regions, Region(left_cyl, -1))
push!(regions, Region(right_cyl, -1))

ex = :(1 ^ 2 ^ 3 ^ 4 ^ 5 ^ (~(7 | 8)) ^ 6 )

push!(cells, Cell(regions, ex))

# Make the inner Region
regions = Array{Region}(0)

push!(regions, Region(top, -1))
push!(regions, Region(bot, -1))
push!(regions, Region(left_cyl, -1))
push!(regions, Region(right_cyl, -1))

ex = :(1 ^ 2 ^ (3 | 4))

push!(cells, Cell(regions, ex))

geometry = Geometry(cells, bounding_box)

@test find_cell_id(Coord(0.41,0.0,0), geometry) == 1
@test find_cell_id(Coord(0.05,0.0,0), geometry) == 2

# Test random ray generation
ray = generate_random_ray(geometry.bounding_box)
@test ray.origin.x >= geometry.bounding_box.lower_left.x
@test ray.origin.y >= geometry.bounding_box.lower_left.y
@test ray.origin.z >= geometry.bounding_box.lower_left.z
@test ray.origin.x <= geometry.bounding_box.upper_right.x
@test ray.origin.y <= geometry.bounding_box.upper_right.y
@test ray.origin.z <= geometry.bounding_box.upper_right.z
ray2 = generate_random_ray(geometry.bounding_box)
@test ray.origin != ray2.origin
@test ray.direction != ray2.direction

# Test intersection routines
ray = Ray(Coord(0.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0) )
new_ray, id, bc_type = find_intersection(ray, geometry)
@test new_ray.origin.x ≈ 0.125
@test new_ray.origin.y ≈ 0.0
@test new_ray.origin.z ≈ 0.0
@test new_ray.direction.x ≈ 1.0
@test new_ray.direction.y ≈ 0.0
@test new_ray.direction.z ≈ 0.0
@test id == 3
@test bc_type == "transmission"

println("Tests complete")


text:
@testset "MutableBinaryHeapWithReduction" begin
    @testset for reduce = (+,*)
        @testset for order = (Forward,Reverse)
            @testset for n = 0
                h = MutableBinaryHeapWithReduction([]; reduce, order)
                @test h.v == []
                @test h.r == []
                @test h.k2i == []
                @test h.i2k == []
            end
            @testset for n = 1
                h = MutableBinaryHeapWithReduction([0]; reduce, order)
                @test h.v == [0]
                @test h.r == [0]
                @test h.k2i == [1]
                @test h.i2k == [1]
            end
            @testset for n = 2:10
                Random.seed!(42)
                for _ = 1:10
                    v̂ = rand(0:99,n)
                    h = MutableBinaryHeapWithReduction(
                        copy(v̂);
                        reduce, order
                    )
                    @test isheap(h.v,order)
                    @test isreduction(h.v,h.r,reduce)
                    @test h.v[h.k2i] == v̂
                    @test h.k2i[h.i2k] == 1:n
                end
            end
        end
    end
end

@testset "update!(::MutableBinaryHeapWithReduction)" begin
    @testset for reduce = (+,#=*=#)
        @testset for order = (Forward,Reverse)
            @testset for n = 0
                h = MutableBinaryHeapWithReduction([],[]; reduce, order)
                update!(h, 1, 0)
                @test h.v == []
                @test h.r == []
                @test h.k2i == []
                @test h.i2k == []
            end
            @testset for n = 1
                h = MutableBinaryHeapWithReduction([0],[-1]; reduce, order)
                update!(h, 1, 1)
                @test h.v == [1]
                @test h.r == [1]
                @test h.k2i == [1]
                @test h.i2k == [1]
            end
            @testset for n = 2:10
                Random.seed!(42)
                v̂ = rand(0:99,n)
                h = MutableBinaryHeapWithReduction(
                    copy(v̂),
                    Vector{Int}(undef, n);
                    reduce, order
                )
                for _ = 1:10
                    i = rand(1:n)
                    vi = rand(0:99)
                    update!(h, i, vi)
                    v̂[i] = vi
                    @test isheap(h.v,order)
                    @test isreduction(h.v,h.r,reduce)
                    @test h.v[h.k2i] == v̂
                    @test h.k2i[h.i2k] == 1:n
                end
            end
        end
    end
end


text:
promote_rule(::Type{DoubleFloat{Float64}}, ::Type{DoubleFloat{Float32}}) = DoubleFloat{Float64}
promote_rule(::Type{DoubleFloat{Float64}}, ::Type{DoubleFloat{Float16}}) = DoubleFloat{Float64}
promote_rule(::Type{DoubleFloat{Float32}}, ::Type{DoubleFloat{Float16}}) = DoubleFloat{Float32}

promote_rule(::Type{DoubleFloat{Float64}}, ::Type{T}) where {T<:IEEEFloat} = DoubleFloat{Float64}
promote_rule(::Type{DoubleFloat{Float32}}, ::Type{T}) where {T<:IEEEFloat} = DoubleFloat{Float32}
promote_rule(::Type{DoubleFloat{Float16}}, ::Type{T}) where {T<:IEEEFloat} = DoubleFloat{Float16}

promote_rule(::Type{DoubleFloat{Float64}}, ::Type{T}) where {T<:Integer} = DoubleFloat{Float64}
promote_rule(::Type{DoubleFloat{Float32}}, ::Type{T}) where {T<:Integer} = DoubleFloat{Float32}
promote_rule(::Type{DoubleFloat{Float16}}, ::Type{T}) where {T<:Integer} = DoubleFloat{Float16}

promote_rule(::Type{DoubleFloat{Float64}}, ::Type{BigInt}) = BigInt
promote_rule(::Type{DoubleFloat{Float32}}, ::Type{BigInt}) = BigInt
promote_rule(::Type{DoubleFloat{Float16}}, ::Type{BigInt}) = BigInt
promote_rule(::Type{BigInt}, ::Type{DoubleFloat{Float64}}) = BigInt
promote_rule(::Type{BigInt}, ::Type{DoubleFloat{Float32}}) = BigInt
promote_rule(::Type{BigInt}, ::Type{DoubleFloat{Float16}}) = BigInt

promote_rule(::Type{DoubleFloat{Float64}}, ::Type{BigFloat}) = BigFloat
promote_rule(::Type{DoubleFloat{Float32}}, ::Type{BigFloat}) = BigFloat
promote_rule(::Type{DoubleFloat{Float16}}, ::Type{BigFloat}) = BigFloat
promote_rule(::Type{BigFloat}, ::Type{DoubleFloat{Float64}}) = BigFloat
promote_rule(::Type{BigFloat}, ::Type{DoubleFloat{Float32}}) = BigFloat
promote_rule(::Type{BigFloat}, ::Type{DoubleFloat{Float16}}) = BigFloat

promote_rule(::Type{AbstractIrrational}, ::Type{DoubleFloat{Float16}}) = DoubleFloat{Float16}
promote_rule(::Type{AbstractIrrational}, ::Type{DoubleFloat{Float32}}) = DoubleFloat{Float32}
promote_rule(::Type{AbstractIrrational}, ::Type{DoubleFloat{Float64}}) = DoubleFloat{Float64}

widen(::Type{DoubleFloat{Float16}}) = DoubleFloat{Float32}
widen(::Type{DoubleFloat{Float32}}) = DoubleFloat{Float64}
widen(::Type{DoubleFloat{Float64}}) = BigFloat