text:
# Routines for extracting segment/face data for tracks (rays) overlaid on a mesh

# Classify a point as on the North, East, South, or West boundary edge of a rectangular mesh
function classify_nesw(p::Point_2D, mesh::UnstructuredMesh_2D)
    ϵ = 1e-4
    # First boundary edge ID -> get the first point in the edge -> Get the desired coordinate
    y_N = mesh.points[mesh.edges[mesh.boundary_edges[1][1]][1]][2]
    x_E = mesh.points[mesh.edges[mesh.boundary_edges[2][1]][1]][1]
    y_S = mesh.points[mesh.edges[mesh.boundary_edges[3][1]][1]][2]
    x_W = mesh.points[mesh.edges[mesh.boundary_edges[4][1]][1]][1]
    if abs(p[2] - y_N) < ϵ
        return 1 # North
    elseif abs(p[1] - x_E) < ϵ
        return 2 # East
    elseif abs(p[2] - y_S) < ϵ
        return 3 # South
    elseif abs(p[1] - x_W) < ϵ
        return 4 # West
    else
        # Used as an index to boundary_edges, so error should be evident
        return 0 # Error
    end
end

# # Get the face indices for all tracks in a single angle
# # Type-stable
# function find_segment_faces_in_angle!(segment_points::Vector{Vector{Point_2D{F}}},
#                                       indices::Vector{Vector{MVector{N, U}}},
#                                       HRPM::HierarchicalRectangularlyPartitionedMesh{F, U}
#                                      ) where {F <: AbstractFloat, U <: Unsigned, N}
#     nt = length(segment_points)
#     nfaces_found = 0
#     # for each track, find the segment indices
#     for it = 1:nt
#         # Returns true if indices were found for all segments in the track
#         nfaces_found += Int64(find_segment_faces_in_track!(segment_points[it], indices[it], HRPM))
#     end
#     return nfaces_found == nt
# end
#
# # Get the face indices for all segments in a single track
# # Type-stable
# function find_segment_faces_in_track!(segment_points::Vector{Point_2D{F}},
#                                       indices::Vector{MVector{N, U}},
#                                       HRPM::HierarchicalRectangularlyPartitionedMesh{F, U}
#                                      ) where {F <: AbstractFloat, U <: Unsigned, N}
#     # Points in the track
#     npoints = length(segment_points)
#     nfaces_found = 0
#     # Test the midpoint of each segment to find the face
#     for iseg = 1:npoints-1
#         p_midpoint = midpoint(segment_points[iseg], segment_points[iseg+1])
#         nfaces_found += Int64(findface!(p_midpoint, indices[iseg], HRPM))
#     end
#     return nfaces_found == npoints - 1
# end
#
# Get the face indices for all tracks in a single angle
function find_segment_faces_in_angle!(segment_points::Vector{Vector{Point_2D}},
                                      segment_faces::Vector{Vector{UInt32}},
                                      mesh::UnstructuredMesh_2D)
    # for each track, find the segment indices
    for it = 1:length(segment_points)
        # Returns true if indices were found for all segments in the track
        find_segment_faces_in_track!(segment_points[it], segment_faces[it], mesh)
    end
    return nothing
end

# Get the face indices for all segments in a single track
function find_segment_faces_in_track!(segment_points::Vector{Point_2D},
                                      segment_faces::Vector{UInt32},
                                      mesh::UnstructuredMesh_2D)
    npoints = length(segment_points)
    # Test the midpoint of each segment to find the face
    for iseg = 1:npoints-1
        p_midpoint = midpoint(segment_points[iseg], segment_points[iseg+1])
        segment_faces[iseg] = findface(p_midpoint, mesh)
    end
    return nothing
end

# # Follows https://mit-crpg.github.io/OpenMOC/methods/track_generation.html
# # Generate tracks with track spacing tₛ for each azimuthal angle in the angular quadrature.
# # These tracks lie within the domain of the mesh.
# # Type-stable
# function generate_tracks(tₛ::F,
#                          ang_quad::ProductAngularQuadrature{nᵧ, nₚ, F},
#                          HRPM::HierarchicalRectangularlyPartitionedMesh{F, U}
#                          ) where {nᵧ, nₚ, F <: AbstractFloat, U <: Unsigned}
#     w = HRPM_width(HRPM)
#     h = HRPM_height(HRPM)
#     # The tracks for each γ
#     tracks = [ generate_tracks(γ, tₛ, w, h) for γ in ang_quad.γ ].data
#     # Shift all tracks if necessary, since the tracks are generated as if the HRPM has a
#     # bottom left corner at (0,0)
#     offset = HRPM.rect.points[1]
#     for iγ = 1:nᵧ
#         for it in 1:length(tracks[iγ])
#             tracks[iγ][it] = LineSegment_2D(tracks[iγ][it].points[1] + offset,
#                                             tracks[iγ][it].points[2] + offset)
#         end
#     end
#     return tracks
# end
#
# Generate tracks with track spacing tₛ for each azimuthal angle in the angular quadrature.
# These tracks lie within the domain of the mesh.
# Follows https://mit-crpg.github.io/OpenMOC/methods/track_generation.html
function generate_tracks(tₛ::Float64,
                         ang_quad::ProductAngularQuadrature{nᵧ, nₚ},
                         mesh::UnstructuredMesh_2D;
                         boundary_shape::String = "Unknown"
                        ) where {nᵧ, nₚ}
    if boundary_shape !== "Rectangle"
        @error "Unsupported boundary shape"
    end
    bb = boundingbox(mesh, boundary_shape = "Rectangle")
    # The tracks for each γ
    tracks = [ generate_tracks(γ, tₛ, width(bb), height(bb)) for γ in ang_quad.γ ].data
    # Shift all tracks if necessary, since the tracks are generated as if the HRPM has a
    # bottom left corner at (0,0)
    offset = bb.bl
    for iγ = 1:nᵧ
        tracks[iγ] .+= offset
    end
    return tracks
end

# Generate tracks for angle γ, with track spacing tₛ for a rectangular domain with width w, height h
# Rectangle has bottom left corner at (0, 0)
function generate_tracks(γ::Float64, tₛ::Float64, w::Float64, h::Float64)
    if γ ≤ π/2
        n_y = ceil(Int64, h*cos(γ)/tₛ)
        # Number of tracks in x direction
        n_x = ceil(Int64, w*sin(γ)/tₛ)
        # Total number of tracks
        nₜ = n_y + n_x
        # Allocate the tracks
        tracks = Vector{LineSegment_2D}(undef, nₜ)
        # Effective angle to ensure cyclic tracks
        γₑ = atan((h*n_x)/(w*n_y))
        t_eff = w*sin(γₑ)/n_x
        Δx = t_eff/sin(γₑ)
        Δy = t_eff/cos(γₑ)
        # Number of tracks in y direction
        # Generate tracks from the bottom edge of the rectangular domain
        for ix = 1:n_x
            x₀ = w - (ix - 0.5)*Δx
            y₀ = 0.0
            # Segment either terminates at the right edge of the rectangle
            # Or on the top edge of the rectangle
            x₁ = min(w, h/tan(γₑ) + x₀)
            y₁ = min((w - x₀)*tan(γₑ), h)
            l = LineSegment_2D(Point_2D(x₀, y₀), Point_2D(x₁, y₁))
            if arclength(l) < minimum_segment_length
                @warn "Small track generated: $l"
            end
            tracks[ix] = l
        end
        # Generate tracks from the left edge of the rectangular domain
        for iy = 1:n_y
            x₀ = 0.0
            y₀ = (iy - 0.5)*Δy
            # Segment either terminates at the right edge of the rectangle
            # Or on the top edge of the rectangle
            x₁ = min(w, (h - y₀)/tan(γₑ))
            y₁ = min(w*tan(γₑ) + y₀, h)
            l = LineSegment_2D(Point_2D(x₀, y₀), Point_2D(x₁, y₁))
            if arclength(l) < minimum_segment_length
                @warn "Small track generated: $l"
            end
            tracks[n_x + iy] = l
        end
    else
        γᶜ = π - γ # Compliment
        # Number of tracks in y direction
        n_y = ceil(Int64, h*cos(γᶜ)/tₛ)
        # Number of tracks in x direction
        n_x = ceil(Int64, w*sin(γᶜ)/tₛ)
        # Total number of tracks
        nₜ = n_y + n_x
        # Allocate the tracks
        tracks = Vector{LineSegment_2D}(undef, nₜ)
        # Effective angle to ensure cyclic tracks
        γₑ = atan((h*n_x)/(w*n_y))
        t_eff = w*sin(γₑ)/n_x
        Δx = t_eff/sin(γₑ)
        Δy = t_eff/cos(γₑ)
        # Generate tracks from the bottom edge of the rectangular domain
        for ix = 1:n_x
            x₀ = (ix - 0.5)*Δx
            y₀ = 0.0
            # Segment either terminates at the left edge of the rectangle
            # Or on the top edge of the rectangle
            x₁ = max(0, -h/tan(γₑ) + x₀)
            y₁ = min(x₀*tan(γₑ), h)
            l = LineSegment_2D(Point_2D(x₀, y₀), Point_2D(x₁, y₁))
            if arclength(l) < minimum_segment_length
                @warn "Small track generated: $l"
            end
            tracks[ix] = l
        end
        # Generate tracks from the right edge of the rectangular domain
        for iy = 1:n_y
            x₀ = w
            y₀ = (iy - 0.5)*Δy
            # Segment either terminates at the left edge of the rectangle
            # Or on the top edge of the rectangle
            x₁ = max(0, w - (h - y₀)/tan(γₑ))
            y₁ = min(w*tan(γₑ) + y₀, h)
            l = LineSegment_2D(Point_2D(x₀, y₀), Point_2D(x₁, y₁))
            if arclength(l) < minimum_segment_length
                @warn "Small track generated: $l"
            end
            tracks[n_x + iy] = l
        end
     end
    return tracks
end

# Get the boundary edge that a point lies on for a rectangular mesh
function get_start_edge_nesw(p::Point_2D, mesh::UnstructuredMesh_2D)
    nesw = classify_nesw(p, mesh)
    if nesw == 1 || nesw == 3
        # On North or South edge. Just check x coordinates
        xₚ = p[1]
        for iedge in mesh.boundary_edges[nesw]
            epoints = edgepoints(mesh.edges[iedge], mesh.points)
            x₁ = epoints[1][1]
            x₂ = epoints[2][1]
            if x₁ ≤ xₚ ≤ x₂ || x₂ ≤ xₚ ≤ x₁
                return iedge
            end
        end
    else # nesw == 2 || nesw == 4
        # On East or West edge. Just check y coordinates
        yₚ = p[2]
        for iedge in mesh.boundary_edges[nesw]
            epoints = edgepoints(mesh.edges[iedge], mesh.points)
            y₁ = epoints[1][2]
            y₂ = epoints[2][2]
            if y₁ ≤ yₚ ≤ y₂ || y₂ ≤ yₚ ≤ y₁
                return iedge
            end
        end
    end
    # Used as index to connectivity, so error should be evident
    return 0x00000000
end

# Get the segment points and the face which the segment lies in for all segments,
# in all tracks in an angle, using the edge-to-edge ray tracing method.
# Assumes a rectangular boundary
function ray_trace_angle_edge_to_edge!(tracks::Vector{LineSegment_2D},
                                       segment_points::Vector{Vector{Point_2D}},
                                       segment_faces::Vector{Vector{UInt32}},
                                       mesh::UnstructuredMesh_2D)
    for it = 1:length(tracks)
        segment_points[it], segment_faces[it] = ray_trace_track_edge_to_edge(tracks[it], mesh)
    end
    return nothing
end

# Get the segment points and the face which the segment lies in for all segments
# in a track, using the edge-to-edge ray tracing method.
# Assumes a rectangular boundary
# Linear edges
function ray_trace_track_edge_to_edge(l::LineSegment_2D, mesh::LinearUnstructuredMesh_2D)
    # Classify line as intersecting north, east, south, or west boundary edge of the mesh
    start_point = l.points[1]; end_point = l.points[2]
    # Find the edges and faces the line starts and ends the mesh on
    start_edge = get_start_edge_nesw(start_point, mesh)
    start_face = mesh.edge_face_connectivity[start_edge][2] # 1st entry should be 0
    end_edge   = get_start_edge_nesw(end_point, mesh)
    end_face   = mesh.edge_face_connectivity[end_edge][2] # 1st entry should be 0
    # Setup for finding the segment points, faces
    segment_points = [start_point]
    segment_faces = UInt32[]
    max_iters = Int64(1E3) # Max iterations of finding intersection points before declaring an error
    current_edge = start_edge
    current_face = start_face
    next_edge = start_edge
    next_face = start_face
    intersection_point = start_point
    last_point = start_point
    end_reached = false
    iters = 0
    # Find the segment points, faces
    while !end_reached && iters < max_iters
        if visualize_ray_tracing
            mesh!(mesh.materialized_faces[current_face], color = (:yellow, 0.2))
        end
        next_edge, next_face, intersection_point = next_edge_and_face(last_point, current_edge,
                                                                        current_face, l, mesh)
        # Could not find next face, or jumping back to a previous face
        if next_face === current_face
            next_edge, next_face = next_edge_and_face_fallback(last_point, current_face,
                                                               segment_faces, l, mesh)
        else
            if visualize_ray_tracing
                mesh!(mesh.materialized_faces[current_face], color = (:green, 0.15))
                scatter!(intersection_point, color = :green)
                linesegments!(mesh.materialized_edges[next_edge], color = :green)
                println("Intersection at point $intersection_point, on face $current_face," *
                        " over edge $next_edge")
            end
            push!(segment_points, intersection_point)
            push!(segment_faces, current_face)
            last_point = intersection_point
        end
        current_edge = next_edge
        current_face = next_face
        # If the furthest intersection is below the minimum segment length to the
        # end point, end here.
        if distance(intersection_point, end_point) < minimum_segment_length ||
                current_face == end_face
            end_reached = true
            if intersection_point ≉  end_point
                push!(segment_points, end_point)
                push!(segment_faces, end_face)
            end
        end
        iters += 1
    end
    if max_iters ≤ iters
        @error "Exceeded max iterations for $l."
    end
    return sort_intersection_points_E2E(l, segment_points, segment_faces)
end

# Get the segment points and the face which the segment lies in for all segments
# in a track, using the edge-to-edge ray tracing method.
# Assumes a rectangular boundary
# Quadratic edges
function ray_trace_track_edge_to_edge(l::LineSegment_2D, mesh::QuadraticUnstructuredMesh_2D)
    # Classify line as intersecting north, east, south, or west boundary edge of the mesh
    start_point = l.points[1]; end_point = l.points[2]
    # Find the edges and faces the line starts and ends the mesh on
    start_edge = get_start_edge_nesw(start_point, mesh)
    start_face = mesh.edge_face_connectivity[start_edge][2] # 1st entry should be 0
    end_edge   = get_start_edge_nesw(end_point, mesh)
    end_face   = mesh.edge_face_connectivity[end_edge][2] # 1st entry should be 0
    # Setup for finding the segment points, faces
    segment_points = [start_point]
    segment_faces = UInt32[]
    max_iters = Int64(1E3) # Max iterations of finding the next point before declaring an error
    current_face = start_face
    next_face = start_face
    intersection_point = start_point
    last_point = start_point
    end_reached = false
    iters = 0
    # Find the segment points, faces
    while !end_reached && iters < max_iters
        if visualize_ray_tracing
            mesh!(mesh.materialized_faces[current_face], color = (:yellow, 0.2))
        end
        next_face, intersection_point = next_face_E2E(last_point, current_face, l, mesh)
        # Could not find next face, or jumping back to a previous face
        if next_face === current_face
            next_face = next_face_fallback(last_point, current_face, segment_faces, l, mesh)
        else
            if visualize_ray_tracing
                mesh!(mesh.materialized_faces[current_face], color = (:green, 0.15))
                scatter!(intersection_point, color = :green)
                println("Intersection at point $intersection_point, on face $current_face")
            end
            push!(segment_points, intersection_point)
            push!(segment_faces, current_face)
            last_point = intersection_point
        end
        current_face = next_face
        # If the furthest intersection is below the minimum segment length to the
        # end point, end here.
        if distance(intersection_point, end_point) < minimum_segment_length
            end_reached = true
            if intersection_point ≉  end_point
                push!(segment_points, end_point)
                push!(segment_faces, end_face)
            end
        end
        iters += 1
    end
    if max_iters ≤ iters
        @error "Exceeded max iterations for $l."
    end
    return sort_intersection_points_E2E(l, segment_points, segment_faces)
end

# Return the next edge, next face, and intersection point on the next edge
# This is for linear, materialized edges
function next_edge_and_face(last_point::Point_2D, current_edge::UInt32, current_face::UInt32,
                            l::LineSegment_2D, mesh::LinearUnstructuredMesh_2D)
    if visualize_ray_tracing # Compile time constant. Compiler will prune branch if not visualizing
        ax = current_axis()
    end
    next_edge = current_edge
    next_face = current_face
    start_point = l.points[1]
    min_distance = distance(start_point, last_point) + minimum_segment_length
    intersection_point = Point_2D(1e10)
    # For each edge in this face, intersect the track with the edge
    for edge in mesh.face_edge_connectivity[current_face]
        # If we used this edge to get to this face, skip it.
        if edge === current_edge
            continue
        end
        if visualize_ray_tracing
            lplot = linesegments!(mesh.materialized_edges[edge], color = :orange)
        end
        # Edges are linear, so 1 intersection point max
        npoints, point = l ∩ mesh.materialized_edges[edge]
        # If there's an intersection
        if npoints === 0x00000001
            # If the intersection point on this edge is further along the ray than
            # the last_point, then we want to leave the face from this edge
            if min_distance < distance(start_point, point)
                intersection_point = point
                # Make sure not to pick the current face for the next face
                if mesh.edge_face_connectivity[edge][1] == current_face
                    next_face = mesh.edge_face_connectivity[edge][2]
                else
                    next_face = mesh.edge_face_connectivity[edge][1]
                end
                next_edge = edge
                break
            end
        end
        if visualize_ray_tracing
            readline()
            delete!(ax.scene, lplot)
        end
    end
    if visualize_ray_tracing
        readline()
    end
    return next_edge, next_face, intersection_point
end

# Return the next face and intersection point on the next edge
# This is for quadratic, materialized edges
function next_face_E2E(last_point::Point_2D, current_face::UInt32, l::LineSegment_2D,
                       mesh::QuadraticUnstructuredMesh_2D)
    if visualize_ray_tracing # Compile time constant. Prune branch if not visualizing
        ax = current_axis()
    end
    next_face = current_face
    start_point = l[1]
    min_distance = distance(start_point, last_point) + minimum_segment_length
    intersection_point = Point_2D(1e10)
    # For each edge in this face, intersect the track with the edge
    for edge in mesh.face_edge_connectivity[current_face]
        if visualize_ray_tracing
            lplot = linesegments!(mesh.materialized_edges[edge], color = :orange)
        end
        npoints, ipoints = l ∩ mesh.materialized_edges[edge]
        # If there's an intersection
        for i = 1:npoints
            point = ipoints[i]
            # If the intersection point on this edge is further along the ray than
            # the last_point, then we want to leave the face from this edge
            if min_distance < distance(start_point, point) < distance(start_point, intersection_point)
                intersection_point = point
                # Make sure not to pick the current face for the next face
                if mesh.edge_face_connectivity[edge][1] == current_face
                    next_face = mesh.edge_face_connectivity[edge][2]
                else
                    next_face = mesh.edge_face_connectivity[edge][1]
                end
            end
        end
        if visualize_ray_tracing
            readline()
            delete!(ax.scene, lplot)
        end
    end
    return next_face, intersection_point
end

# Return the next face, and edge to skip, for the edge-to-edge algorithm to check in the event
# that they could not be determined simply by checking edges of the current face. This typically
# means there was an intersection around a vertex, and floating point error is
# causing an intersection to erroneously return 0 points intersected.
# Requires materialized faces
function next_edge_and_face_fallback(last_point::Point_2D, current_face::UInt32,
                                     segment_faces::Vector{UInt32}, l::LineSegment_2D,
                                     mesh::LinearUnstructuredMesh_2D)
    # If the next face could not be determined, or the ray is jumping back to the
    # previous face, this means either:
    # (1) The ray is entering or exiting through a vertex, and floating point error
    #       means the exiting edge did not register an intersection.
    # (2) You're supremely unlucky and a fallback method kicked you to another face
    #       where the next face couldn't be determined
    next_face = current_face
    start_point = l[1]
    # The closest point along l intersected in this iteration
    closest_point = Point_2D(1e10)

    # Check adjacent faces first to see if that is sufficient to solve the problem
    next_face, closest_point = adjacent_faces_fallback(last_point, current_face, l, mesh)
    # If adjacent faces were not sufficient, try all faces sharing the vertices of this face
    if next_face === current_face || next_face ∈  segment_faces
        next_face, closest_point = shared_vertex_fallback(last_point, current_face, l, mesh)
        if visualize_ray_tracing
            readline()
        end
    end

    # If the next face STILL couldn't be determined, you're screwed
    if next_face == current_face
        @error "Could not find next face, even using fallback methods, for segment $l."
    end
    # Determine the edge that should be skipped by choosing the edge with intersection point
    # closest to the start of the line.
    next_edge = skipped_edge_fallback(next_face, l, mesh)
    return next_edge, next_face
end

# Return the next face for the edge-to-edge algorithm to check in the event
# that it could not be determined simply by checking edges of the current face. This typically
# means there was an intersection around a vertex, and floating point error is
# causing an intersection to erroneously return 0 points intersected.
# Requires materialized faces
# Quadratic edges
function next_face_fallback(last_point::Point_2D, current_face::UInt32, segment_faces::Vector{UInt32},
            l::LineSegment_2D, mesh::QuadraticUnstructuredMesh_2D)
    # If the next face could not be determined, this means either:
    # (1) The ray is entering or exiting through a vertex, and floating point error
    #       means the exiting edge did not register an intersection.
    # (2) You're supremely unlucky and a fallback method kicked you to another face
    #       where the next face couldn't be determined
    next_face = current_face
    start_point = l[1]
    # The closest point along l intersected in this iteration
    closest_point = Point_2D(1e10)
    next_face, closest_point = quadratic_shared_vertex_fallback(last_point, current_face, l, mesh)
    if visualize_ray_tracing
        readline()
    end
    # If the next face STILL couldn't be determined, you're screwed
    if next_face === current_face
        @error "Could not find next face, even using fallback methods, for segment $l."
    end
    return next_face
end

# Check to see if one of the adjacent faces should be the next face in edge-to-edge ray tracing
function adjacent_faces_fallback(last_point::Point_2D, current_face::UInt32,
                                 l::LineSegment_2D, mesh::UnstructuredMesh_2D)
    next_face = current_face
    start_point = l.points[1]
    min_distance = distance(start_point, last_point) + minimum_segment_length
    intersection_point = Point_2D(1e10)
    for face in adjacent_faces(current_face, mesh)
        npoints, ipoints = l ∩ mesh.materialized_faces[face]
        if 0 < npoints
            for point in ipoints[1:npoints]
                if min_distance < distance(start_point, point) < distance(start_point, intersection_point)
                    intersection_point = point
                    next_face = face
                end
            end
        end
    end
    return next_face, intersection_point
end

# Check to see if one of the faces sharing a vertex with the current face
# should be the next face in edge-to-edge ray tracing
function shared_vertex_fallback(last_point::Point_2D, current_face::UInt32,
                                l::LineSegment_2D, mesh::UnstructuredMesh_2D)
    next_face = current_face
    start_point = l.points[1]
    min_distance = distance(start_point, last_point) + minimum_segment_length
    intersection_point = Point_2D(1e10)
    # Get the vertex ids for each vertex in the face
    already_tested = adjacent_faces(current_face, mesh)
    push!(already_tested, current_face)
    vertex_ids = mesh.faces[current_face]
    faces_OI = UInt32[] # faces of interest
    for vertex in vertex_ids
        candidate_faces = faces_sharing_vertex(vertex, mesh)
        for face in candidate_faces
            if face ∉ already_tested && face ∉ faces_OI
                push!(faces_OI, face)
            end
        end
    end
    if visualize_ray_tracing
        mesh_vec = []
    end
    for face in faces_OI
        npoints, ipoints = l ∩ mesh.materialized_faces[face]
        if visualize_ray_tracing
            println("Face: $face")
            push!(mesh_vec, mesh!(mesh.materialized_faces[face], color = (:black, 0.2)))
            readline()
        end
        if 1 < npoints
            contains_last_point = false
            for point in ipoints[1:npoints]
                if distance(point, last_point) < minimum_segment_length
                    contains_last_point = true
                end
            end
            for point in ipoints[1:npoints]
                distance_to_point = distance(start_point, point)
                # If this is a valid intersection point
                if min_distance < distance_to_point
                    # If this point is the closest point, use this
                    if distance_to_point ⪉  distance(start_point, intersection_point)
                        intersection_point = point
                        next_face = face
                    # If this face contains the last point, we want to prioritize this face
                    elseif contains_last_point
                        next_face = face
                    end
                end
                if visualize_ray_tracing
                    scatter!(point, color = :yellow)
                    readline()
                end
            end
        end
    end

    if visualize_ray_tracing
        ax = current_axis()
        for m in mesh_vec
            delete!(ax.scene, m)
        end
    end
    return next_face, intersection_point
end

# Check to see if one of the faces sharing a vertex with the current face
# should be the next face in edge-to-edge ray tracing
function quadratic_shared_vertex_fallback(last_point::Point_2D, current_face::UInt32,
                                          l::LineSegment_2D, mesh::QuadraticUnstructuredMesh_2D)
    next_face = current_face
    start_point = l[1]
    min_distance = distance(start_point, last_point) + minimum_segment_length
    intersection_point = Point_2D(1e10)
    # Get the vertex ids for each vertex in the face
    vertex_ids = mesh.faces[current_face]
    faces_OI = UInt32[] # faces of interest
    for vertex in vertex_ids
        candidate_faces = faces_sharing_vertex(vertex, mesh)
        for face in candidate_faces
            if face !== current_face && face ∉ faces_OI
                push!(faces_OI, face)
            end
        end
    end
    if visualize_ray_tracing
        mesh_vec = []
    end
    next_face_contains_last_point = false
    for face in faces_OI
        npoints, ipoints = l ∩ mesh.materialized_faces[face]
        if visualize_ray_tracing
            println("Face: $face")
            push!(mesh_vec, mesh!(mesh.materialized_faces[face], color = (:black, 0.2)))
            readline()
        end
        if 1 < npoints
            contains_last_point = false
            for point in ipoints[1:npoints]
                if distance(point, last_point) < minimum_segment_length
                    contains_last_point = true
                end
            end
            if visualize_ray_tracing
                println(contains_last_point ? "Contains last point" : "Does not contain last point")
            end
            for point in ipoints[1:npoints]
                distance_to_point = distance(start_point, point)
                # If this is a valid intersection point
                if min_distance < distance_to_point
                    # If this point is the closest point, use this
                    # Could be a problem for Float32.
                    if distance_to_point ⪉ distance(start_point, intersection_point)
                        if visualize_ray_tracing
                            println("New intersection point")
                            println("New face: $face")
                        end
                        intersection_point = point
                        next_face = face
                        if contains_last_point
                            next_face_contains_last_point = true
                        end
                    # If this face contains the last point, we want to prioritize this face,
                    # but only if the face with the closest point doesn't also contain
                    # the last point
                    elseif contains_last_point && !next_face_contains_last_point
                        next_face = face
                        if visualize_ray_tracing
                            println("New face: $face")
                        end
                    end
                end
                if visualize_ray_tracing
                    println("Point: $point")
                    scatter!(point, color = :yellow)
                    readline()
                end
            end
        end
    end

    if visualize_ray_tracing
        ax = current_axis()
        for m in mesh_vec
            delete!(ax.scene, m)
        end
        println("Next face: $next_face, Intersection: $intersection_point")
    end
    return next_face, intersection_point
end

# Determine the edge that should be skipped by choosing the edge with intersection point
# closest to the start of the line.
# Linear edges
function skipped_edge_fallback(next_face::UInt32, l::LineSegment_2D, mesh::LinearUnstructuredMesh_2D)
    start_point = l.points[1]
    closest_point = Point_2D(1.0e10)
    edges_OI = mesh.face_edge_connectivity[next_face]
    next_edge = 0x00000000
    for iedge in edges_OI
        npoints, point = l ∩ mesh.materialized_edges[iedge]
        if 0 < npoints
            if  distance(start_point, point) < distance(start_point, closest_point)
                closest_point = point
                next_edge = iedge
            end
        end
    end
    return next_edge
end

function sort_intersection_points_E2E(l::LineSegment_2D, segment_points::Vector{Point_2D},
                                      segment_faces::Vector{UInt32})
    if 2 <= length(segment_points)
        npoints = length(segment_points)
        perm = sortperm(distance².(l[1], view(segment_points, 2:npoints)))
        permute!(view(segment_points, 2:npoints), perm)
        permute!(segment_faces, perm)
        # Eliminate any points and faces for which the distance between consecutive points
        # is less than the minimum segment length
        delete_ids = Int64[]
        id_start = 1
        for id_stop ∈ 2:npoints
            if distance²(segment_points[id_start], segment_points[id_stop]) < minimum_segment_length^2
                push!(delete_ids, id_stop)
            else
                id_start = id_stop
            end
        end
        deleteat!(segment_points, delete_ids)
        deleteat!(segment_faces, delete_ids .- 1)
        return segment_points, segment_faces
    else
        return segment_points, segment_faces
    end
end

function use_E2E_raytracing(mesh::UnstructuredMesh_2D)
    return 0 < length(mesh.boundary_edges)
end


text:
struct SMCState{ioT<:SMCIO, modelT<:SMCModel,
    lMT<:Union{Function,Nothing}, setT<:Function, refT}
    model::modelT
    io::ioT
    lM::lMT
    set_param!::setT
    ref::refT
end

"""
   state = SMCState(T, N, ParticleType, ScratchType, set_param!, lG, M!,
                    [lM=nothing]; kwargs...)

Data structures for Sequential Monte Carlo.

# Arguments:
- `T`: Length of the data record (or the hidden Markov model).
- `N`: Number of particles.
- `ParticleType`: Data type of particles.
- `ScratchType`: Data type of scratch (auxiliary data) object.
- `set_param!`: Function which sets the parameters of the model:
  `set_param!(scratch, par)` sets `scratch` to match parameters `par`.
- `lG`: Log-potential function: `lG(k, x, scratch)` returns value of potential
  at time `k` for particle `x` given `scratch`.
- `M!`: Function which simulates from the transition density:
  `M!(x, rng, k, x_, scratch)` evolves `x_` from time `k-1` to `x` at time `k`
  given random number generator `rng` and `scratch`.
- `lM`: Function which returns transition density values:
  `lM(k, x_, x, scratch)` where arguments are as above. (Only required for PG)

# Keyword arguments:
- `nthreads`: Number of threads used within the particle filter; default `1`.
- `essThreshold`: Effective sample size threshold for resampling; default `2.0`.

For details, see documentation of `SequentialMonteCarlo.jl`
"""
function SMCState(T::Int, N::Int, ParticleType, ScratchType,
    set_param!::Function, lG::Function, M!::Function,
    lM::Union{Function,Nothing}=nothing; nthreads=1, essThreshold=2.0)
    model = SMCModel(M!, lG, T, ParticleType, ScratchType)
    io = SMCIO{ParticleType,ScratchType}(N, T, nthreads, true, essThreshold)
    ref = [model.particle() for i=1:T]
    SMCState(model, io, lM, set_param!, ref)
end

# Helper which calculates model log-likelihood for given reference
@inline function _reference_log_likelihood(state::SMCState)
    ref = state.ref; lG = state.model.lG; lM = state.lM;
    pScratch = state.io.internal.particleScratch
    L = 0.0
    @inbounds for k=1:length(ref)
        k_ = max(1, k-1); x = ref[k]; x_ = ref[k_]
        L += lG(k, x, pScratch)
        L += lM(k, x_, x, pScratch)
    end
    L
end

# Wrappers to SequentialMonteCarlo interface
@inline function _set_model_param!(state::SMCState, theta)
    state.set_param!(state.io.internal.particleScratch, theta)
    nothing
end
@inline function _run_smc!(state::SMCState)
    smc!(state.model, state.io)
end
@inline function _log_likelihood(state::SMCState)
    state.io.logZhats[end]
end
@inline function _pick_particle!(state::SMCState)
    SequentialMonteCarlo.pickParticle!(state.ref, state.io)
end
@inline _save_reference!(state::SMCState) = nothing
@inline function _init_path_storage(state::SMCState, nsim)
    [[state.model.particle() for i=1:state.io.n] for j=1:nsim]
end
@inline function _copy_reference!(out, state::SMCState)
    SequentialMonteCarlo._copyParticles!(out, state.ref)
end
@inline function _run_csmc!(state::SMCState, backward_sampling::Bool)
    csmc!(state.model, state.io, state.ref, state.ref)
    if backward_sampling
        SequentialMonteCarlo.pickParticleBS!(state.ref, state.io, state.lM)
    end
    nothing
end

ext: jl


text:
using MNIST

type DataLoader
    cur_id::Int
    order::Vector{Int}
end

DataLoader() = DataLoader(1, shuffle(1:60000))

function next_batch(loader::DataLoader, batch_size)
    x = zeros(Float32, batch_size, 784)
    y = zeros(Float32, batch_size, 10)
    for i in 1:batch_size
        x[i, :] = trainfeatures(loader.order[loader.cur_id])
        label = trainlabel(loader.order[loader.cur_id])
        y[i, Int(label)+1] = 1.0
        loader.cur_id += 1
        if loader.cur_id > 60000
            loader.cur_id = 1
        end
    end
    x, y
end

function load_test_set(N=10000)
    x = zeros(Float32, N, 784)
    y = zeros(Float32, N, 10)
    for i in 1:N
        x[i, :] = testfeatures(i)
        label = testlabel(i)
        y[i, Int(label)+1] = 1.0
    end
    x,y
end


text:
using Test
using Yota
using Yota: Tape, Input, Call, Constant, trace, ctrace, irtrace, play!, transform, binarize_ops
using Yota: ∇mean, setfield_nested!, copy_with, simplegrad, remove_unused
using Yota: find_field_source_var, unwind_iterate, eliminate_common
using Yota: unvcat, unhcat, uncat

include("test_trace_cassette.jl")
include("test_trace_irtools.jl")
# include("test_trace_interp.jl")  # doesn't work for Julia 1.4, maybe will drop support at all
include("gradcheck.jl")
include("test_helpers.jl")
include("test_simple.jl")
include("test_grad.jl")
include("test_distributions.jl")
include("test_dynamic.jl")
include("test_simplegrad.jl")
include("test_update.jl")
include("test_transform.jl")
include("test_examples.jl")


text:
# Playing Cards: Classic Game of War!
using PlayingCards

function game()
    rounds = 0
    deck = shuffle(Deck());
    player1 = Hand(cards = deal(deck, div(length(deck), 2)))
    player2 = Hand(cards = deal(deck, length(deck)))
    while true
        if isempty(player1)
            println("Player 2 wins in $rounds rounds!")
            break;
        end
        if isempty(player2)
            println("Player 1 wins in $rounds rounds!")
            break;
        end
        rounds += 1
        play!(player1, player2)
        if rounds == 1500
            println("We have played $rounds rounds, let's stop!")
            break;
        end
        if rounds % 100 == 0
            println("We have played $rounds rounds, let's shuffle our hands!")
            shuffle(player1); shuffle(player2);
        end
    end
    return nothing;
end

function play!(player1::Hand, player2::Hand; cards::Array{Card, 1} = Card[])
    card1 = player1[1]
    card2 = player2[1]
    if card1 > card2
        append!(cards, [card1, card2])
        player1 - card1; player2 - card2;
        player1 + cards
    elseif card2 > card1
        append!(cards, [card1, card2])
        player1 - card1; player2 - card2;
        player2 + cards
    else # card1 == card2
        println("WAR!")
        if length(player1) > 1
            player1 - card1
            push!(cards, card1)
        end
        if length(player2) > 1
            player2 - card2
            push!(cards, card2)
        end
        war!(player1, player2, cards)
    end
    return nothing;
end

function war!(player1::Hand, player2::Hand, cards::Array{Card, 1})
    a = (length(player1) > 3) ? 3 : length(player1) - 1
    b = (length(player2) > 3) ? 3 : length(player2) - 1
    append!(cards, player1[1:a])
    append!(cards, player2[1:b])
    player1 - player1[1:a]
    player2 - player2[1:b]
    play!(player1, player2, cards = cards)
    return nothing;
end

max_stars_repo_name: Joshua-Robison/PlayingCards

max_issues_repo_name: Joshua-Robison/PlayingCards

max_forks_repo_name: Joshua-Robison/PlayingCards


text:
@testset "Testing method dispatch..." begin

# LinearGaussianSSM
#    - LinearKalmanFilter
@test hasmethod(predict, (LinearGaussianSSM, AbstractMvNormal))
@test hasmethod(update, (LinearGaussianSSM, AbstractMvNormal, Vector))
@test hasmethod(update,
    (LinearGaussianSSM, AbstractMvNormal, Vector, Int))
@test hasmethod(update,
    (LinearGaussianSSM, AbstractMvNormal, Vector, KalmanFilter, Int))

# NonlinearGaussianSSM
#    - NonlinearKalmanFilter
@test hasmethod(update,
    (NonlinearGaussianSSM, AbstractMvNormal, Vector))

for T in [NonlinearKalmanFilter, ExtendedKalmanFilter, UnscentedKalmanFilter]
    println(T)
    @test hasmethod(update,
    (NonlinearGaussianSSM, AbstractMvNormal, Vector, T))
    @test hasmethod(update,
        (NonlinearGaussianSSM, AbstractMvNormal, Vector, T, Int))
end

# Special predict methods that don't return just a Distribution
@test hasmethod(predict, (NonlinearGaussianSSM, AbstractMvNormal, UKF))
@test hasmethod(predict, (NonlinearGaussianSSM, Matrix, EnKF))

@test hasmethod(update,
    (NonlinearGaussianSSM, Matrix, Vector, EnKF))

# NonlinearSSM             NonlinearFilter

end


text:
include("plotting_utils.jl")

files = ["local", "nonlocal", "nonlocal_symmetric", "nonlocal_less_diffusivity", "nonlocal_more_velocity", "nonlocal_less_transition", "semilocal", "nonlocal_more_diffusivity_less_transition", "nonlocal_reverse"]

total_time = 0
for file in files
    println("------")
    println("currently on file ", file)
    println("-----")
    tic = Base.time()
    fig = file_kernel_plot(file * ".jld2")
    save("data/" * file * "_kernel.png", fig)
    fig = file_local_diffusivity_plot(file * ".jld2")
    save("data/" * file * "_local_diffusivity.png", fig)
    toc = Base.time()

    global total_time += toc - tic
end
println("total time to generate plots ", total_time, " seconds ")

max_stars_repo_name: sandreza/StatisticalNonlocality

max_issues_repo_name: sandreza/StatisticalNonlocality

max_forks_repo_name: sandreza/StatisticalNonlocality


text:
import HTTP
#import JSON3

HTTP.request("POST", "http://localhost:8000/jsonpayload",
    [("Content-Type", "application/json")], """{"name":"Adrian"}""")

#
# import HTTP
# HTTP.request("POST", "http://rwresearch.ca/jsonpayload",
#     [("Content-Type", "application/json")], """{"name":"Adrian"}""")

max_stars_repo_name: RoyCCWang/WebAppRKHS

max_issues_repo_name: RoyCCWang/WebAppRKHS

max_forks_repo_name: RoyCCWang/WebAppRKHS


text:
export makeSampler, Sampler, sample, sample!, call

type Sampler
    image::Array{Float32,2}
    ind::Array{Int,1}
    mi::Array{Int,2}
    ma::Array{Int,2}
    clampedpos::Array{Int,2}
    buf::Array{Float32,2}
    normmeanstd::Bool
end

function makeSampler(image, bsize::Int, scale = 1.; centered::Bool = true, col::Bool = false, normmeanstd::Bool = false)
    patchsize = tuple(repeat(bsize, ndims(image))...)
    if centered
        grid = @p centeredmeshgrid patchsize | times scale | plus 1 |  asint
    else
        grid = @p meshgrid patchsize | times scale | asint
    end
    ind = @p map grid subtoind image | minus 1

    buf = zeros(eltype(image), patchsize)
    if col buf = FunctionalData.col(buf) end

    mi = ones(Int,ndims(image),1) - minimum(grid,2) + 1
    ma = siz(image) - maximum(grid,2) + 1
    clampedpos = zeros(Int, size(mi))

    Sampler(image, ind, mi, ma, clampedpos, buf, normmeanstd)
end

@static if VERSION < v"0.5-"
    call(a::Sampler, pos) = sample(pos, a)
else
    (a::Sampler)(pos) = sample(pos, a)
end

sample(a::Sampler,b::Sampler) = error("Only one of the params can be a sampler.")
sample(pos, a::Sampler) = sample!(a.buf, pos, a::Sampler)
sample(a::Sampler, pos) = sample!(a.buf, pos, a::Sampler)

function sample!(buf, pos::Array{Int,2}, a::Sampler)
    c = @p clamp! a.clampedpos pos a.mi a.ma | subtoind a.image
    for i = 1:length(a.buf)
        buf[i] = a.image[c+a.ind[i]]
    end
    a.normmeanstd && normmeanstd!(buf)
    buf
end


text:
using StochasticDiffEq, DiffEqProblemLibrary, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
prob = DiffEqProblemLibrary.SDEProblemLibrary.oval2ModelExample(largeFluctuations=true,useBigs=false)

Random.seed!(250)
prob = remake(prob,tspan = (0.0,500.0))
sol = solve(prob,SRIW1(),dt=(1/2)^(18),progress=true,qmax=1.125,
          saveat=0.1,abstol=1e-5,reltol=1e-3,maxiters=1e7);
Random.seed!(250)
prob = remake(prob,tspan = (0.0,500.0))
@time sol = solve(prob,SRIW1(),dt=(1/2)^(18),progress=true,qmax=1.125,
    saveat=0.1,abstol=1e-5,reltol=1e-3,maxiters=1e7);

println(maximum(sol[:,2]))
using Plots; gr()
lw = 2
lw2 = 3
p1 = plot(sol,vars=(0,16),
          title="(A) Timeseries of Ecad Concentration",xguide="Time (s)",
          yguide="Concentration",guidefont=font(16),tickfont=font(16),
          linewidth=lw,leg=false)

p2 = plot(sol,vars=(0,17),
          title="(B) Timeseries of Vim Concentration",xguide="Time (s)",
          yguide="Concentration",guidefont=font(16),
          tickfont=font(16),linewidth=lw,leg=false)

prob = remake(prob,tspan = (0.0,1.0))
## Little Run
sol = solve(prob,EM(),dt=(1/2)^(20),
          progressbar=true,saveat=0.1)
println("EM")
@time sol = solve(prob,EM(),dt=(1/2)^(20),
          progressbar=true,saveat=0.1)

sol = solve(prob,SRI(),dt=(1/2)^(18),adaptive=false,
            progressbar=true,save_everystep=false)
println("SRI")
@time sol = solve(prob,SRI(),dt=(1/2)^(18),adaptive=false,
          progressbar=true,save_everystep=false)

sol = solve(prob,SRIW1(),dt=(1/2)^(18),adaptive=false,
          adaptivealg=:RSwM3,progressbar=false,qmax=4,saveat=0.1)
println("SRIW1")
@time sol = solve(prob,SRIW1(),dt=(1/2)^(18),adaptive=false,
          adaptivealg=:RSwM3,progressbar=false,qmax=4,saveat=0.1)

sol = solve(prob,SRI(),dt=(1/2)^(18),
          adaptivealg=:RSwM3,progressbar=false,qmax=1.125,
          saveat=0.1,abstol=1e-6,reltol=1e-4)
println("SRI Adaptive")
@time sol = solve(prob,SRI(),dt=(1/2)^(18),
          adaptivealg=:RSwM3,progressbar=false,qmax=1.125,
          saveat=0.1,abstol=1e-6,reltol=1e-4)
@show length(sol.t)

sol = solve(prob,SRIW1(),dt=(1/2)^(18),
          adaptivealg=:RSwM3,progressbar=false,qmax=1.125,
          saveat=0.1,abstol=1e-6,reltol=1e-4)
println("SRIW1 Adaptive")
@time sol = solve(prob,SRIW1(),dt=(1/2)^(18),
          adaptivealg=:RSwM3,progressbar=false,qmax=1.125,
          saveat=0.1,abstol=1e-6,reltol=1e-4)
@show length(sol.t)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
module Vars

jets = mapreduce(vcat, 0:5) do i
    ["Jet_Px[$i]", "Jet_Py[$i]", "Jet_Pz[$i]", "Jet_E[$i]"]
end
tight_lepton = "TightLepton_" .* ["Px", "Py", "Pz", "E"] .* "[0]"
jet_csv = string.("Jet_CSV[", 0:5, "]")

scalar_features = [
    "Evt_CSV_avg_tagged",
    "Evt_CSV_avg",
    "Evt_M_minDrLepTag",
    "Evt_Pt_JetsAverage",
    "Evt_Pt_minDrTaggedJets",
    "Evt_Deta_TaggedJetsAverage",
    "Evt_Deta_JetsAverage",
    #"N_Jets",
    "Reco_JABDT_ttbar_Jet_CSV_whaddau2",
    "Reco_ttbar_toplep_m",
    "Reco_tHq_bestJABDToutput",
    "Reco_JABDT_tHq_abs_ljet_eta",
    "Reco_JABDT_tHq_Jet_CSV_hdau1",
    "Reco_JABDT_tHW_Jet_CSV_btop",
    #"memDBp",
]

all_features = [jets; tight_lepton; jet_csv; scalar_features]

weights = [
    "Weight_XS",
    #"Weight_btagSF",
    "Weight_CSV",
    "Weight_GEN_nom",
    "lumiWeight",
]

classes = [:ttH, :ttbb, :ttcc, :ttlf]

names_lola_out = [
    mapreduce(vcat, 0:5) do i
        ["Jet_m²[$i]"; "Jet_p_T[$i]"; "Jet_weighted_E[$i]"; string.("Jet_weighted_d", 1:4, "[$i]")]
    end
    ["Lepton_m²"; "Lepton_p_T"; "Lepton_weighted_E"; string.("Lepton_weighted_d", 1:4)]
]

tier1 = [
    "Evt_CSV_avg_tagged",
    "Evt_CSV_avg",
    "Evt_M_minDrLepTag",
    "Evt_Pt_JetsAverage",
    "Evt_Pt_minDrTaggedJets",
    "Evt_Deta_TaggedJetsAverage",
]
tier2 = [
    "Evt_Deta_JetsAverage",
]
tier3 = [
    "Reco_JABDT_ttbar_Jet_CSV_whaddau2",
    "Reco_ttbar_toplep_m",
    "Reco_tHq_bestJABDToutput",
    "Reco_JABDT_tHq_abs_ljet_eta",
    "Reco_JABDT_tHq_Jet_CSV_hdau1",
    "Reco_JABDT_tHW_Jet_CSV_btop",
]

end


text:
using Documenter
using RenoiseOSC

makedocs(sitename="RenoiseOSC.jl", modules=[RenoiseOSC], pages=["index.md", "api.md", "config.md"])
deploydocs(repo="github.com/stellartux/RenoiseOSC.jl.git")

max_stars_repo_name: stellartux/RenoiseOSC

max_issues_repo_name: stellartux/RenoiseOSC

max_forks_repo_name: stellartux/RenoiseOSC


text:
using NestedTuples: lazymerge
import MeasureTheory

import MeasureBase: insupport
export insupport

# function MeasureBase.insupport(c::AbstractConditionalModel{A,B,M}, x=NamedTuple()) where {A,B,M}
#     _insupport(M, model(c), argvals(c), observations(c), x)
# end

# export sourceInsupport

# sourceInsupport(m::AbstractModel) = sourceInsupport()(model(m))

# function sourceInsupport()
#     function(_m::AbstractModel)
#         proc(_m, st :: Assign)     = :($(st.x) = $(st.rhs))
#         proc(_m, st :: Return)     = nothing
#         proc(_m, st :: LineNumber) = nothing
#         function proc(_m, st :: Sample)
#             x = st.x
#             rhs = st.rhs
#             @q begin
#                 Tilde.dynamic(Tilde.insupport($rhs, $x)) || return false
#                 $x = Tilde.predict($rhs, $x)
#             end
#         end

#         wrap(kernel) = @q begin
#             $kernel
#             return true
#         end

#         buildSource(_m, proc, wrap) |> MacroTools.flatten
#     end
# end

# # MeasureTheory.insupport(d::Distribution, val, tr) = logdensityof(d, val)


text:
@doc raw"""
    Segment3

An object `s` of the data type [`Segment3`](@ref) is a directed straight line
segment in the three-dimensional Euclidean space ``Ε^3``, i.e. a straight line
segment ``[p, q]`` connecting two points ``p, q ∈ ℝ^3``.

The segment is topologically closed, i.e. the end points belong to it. Point `p`
is called the *source* and `q` is called the *target* of `s`. The length of `s`
is the Euclidean distance between `p` and `q`. Note that there is only a
function to compute the square of the length, because otherwise we had to
perform a square root operation which is expensive, and may not be exact.
"""
Segment3

"""
    Segment3(p::Point3, q::Point3)

Introduces a segment `s` with source `p` and target `q`.

The segment is directed from the source towards the target.
"""
Segment3(p::Point3, q::Point3)

"""
    ==(s::Segment3, q::Segment3)

Test for equality: Two segments are equal, iff their sources and targets are
equal.
"""
==(s::Segment3, q::Segment3)

"""
    source(s::Segment3)

The source of `s`.
"""
source(s::Segment3)

"""
    target(s::Segment3)

The target of `s`.
"""
target(s::Segment3)

"""
    min(s::Segment3)

Returns the point of `s` with lexicographically smallest coordinate.
"""
min(s::Segment3)

"""
    max(s::Segment3)

Returns the point of `s` with lexicographically largest coordinate.
"""
max(s::Segment3)

@doc raw"""
    vertex(s::Segment3, i::Integer)

Returns source of target or `s`: `vertex(s, 0)` returns the source of `s`,
`vertex(s, 1)` returns the target of `s`

The parameter `i` is taken module 2, which gives easy access to the other
vertex.
"""
vertex(s::Segment3, i::Integer)

"""
    point(s::Segment3, i::Integer)

Returns `vertex(s, i)`.
"""
point(s::Segment3, i::Integer)

"""
    squared_length(s::Segment3)

Returns the squared length of `s`.
"""
squared_length(s::Segment3)

"""
    to_vector(s::Segment3)

Returns the vector `target(s) - source(s)`.
"""
to_vector(s::Segment3)

"""
    direction(s::Segment3)

Returns the direction from source to target.
"""
direction(s::Segment3)

"""
    opposite(s::Segment3)

Returns a segment with source and target point interchanged.
"""
opposite(s::Segment3)

"""
    supporting_line(s::Segment3)

Returns the line `l` passing through `s`.

Line `l` has the same orientation as segment `s`, that is from the source to the
target of `s`.
"""
supporting_line(s::Segment3)

"""
    is_degenerate(s::Segment3)

Segment `s` is degenerate, if source and target fall together.
"""
is_degenerate(s::Segment3)

"""
    has_on(s::Segment3, p::Point3)

A point is on `s`, iff it is equal to the source or target of `s` or if it is
in the interior of `s`.
"""
has_on(s::Segment3, p::Point3)

"""
    bbox(s::Segment3)

Returns a bounding box containing `s`.
"""
bbox(s::Segment3)

"""
    transform(s::Segment3, t::AffTransformation3)

Returns the segment obtained by applying `t` on the source and the target of
`s`.
"""
transform(s::Segment3, t::AffTransformation3)


text:
# To view in browser start a server in the build dir:
# python -m http.server --bind localhost

using Documenter
using QuasinormalModes

makedocs(sitename = "QuasinormalModes.jl",
    modules = [QuasinormalModes],
    pages = [
        "index.md",
        "intro.md",
        "org.md",
        "schw.md",
        "sho.md",
        "api_ref.md"
    ]
)

deploydocs(
    repo = "github.com/lucass-carneiro/QuasinormalModes.jl.git",
)


text:
function _gen_temp_env(envdir)

    mkpath(envdir)

    # proj
    gwtproj = joinpath(envdir, "Project.toml")
    write(gwtproj,
        string(
            "[deps]", "\n",
            "GitWorkers = \"9e5d5e8c-ba81-11ea-39e8-15323c797f61\"", "\n"
        )
    )

    # manif (The same as GitWorkers)
    gw_pkgdir = pkgdir(GitWorkers)
    gwmanf = joinpath(gw_pkgdir, "Manifest.toml")
    gwtmanf = joinpath(envdir, "Manifest.toml")
    _cp(gwmanf, gwtmanf)

    return envdir
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
module superoperators

import Base: ==, *, /, +, -
import ..operators.check_samebases

using ..bases, ..operators, ..operators_sparse

export DenseSuperOperator, SparseSuperOperator, spre, spost, liouvillian

"""
Base class for all super operator classes.

Super operators are bijective mappings from operators given in one specific
basis to operators, possibly given in respect to another, different basis.
To embed super operators in an algebraic framework they are defined with a
left hand basis ``basis_l`` and a right hand basis ``basis_r`` where each of
them again consists of a left and right hand basis.

.. math::

    A_{bl_1,bl_2} &= S_{(bl_1,bl_2)<->(br_1,br_2)} B_{br_1,br_2}
    \\
    A_{br_1,br_2} &= B_{bl_1,bl_2} S_{(bl_1,bl_2)<->(br_1,br_2)}
"""
abstract SuperOperator

type DenseSuperOperator <: SuperOperator
    basis_l::Tuple{Basis, Basis}
    basis_r::Tuple{Basis, Basis}
    data::Matrix{Complex128}
    function DenseSuperOperator(basis_l::Tuple{Basis, Basis}, basis_r::Tuple{Basis, Basis}, data::Matrix{Complex128})
        if length(basis_l[1])*length(basis_l[2]) != size(data, 1) || length(basis_r[1])*length(basis_r[2]) != size(data, 2)
            throw(DimensionMismatch())
        end
        new(basis_l, basis_r, data)
    end
end

type SparseSuperOperator <: SuperOperator
    basis_l::Tuple{Basis, Basis}
    basis_r::Tuple{Basis, Basis}
    data::SparseMatrixCSC{Complex128}
    function SparseSuperOperator(basis_l::Tuple{Basis, Basis}, basis_r::Tuple{Basis, Basis}, data::SparseMatrixCSC{Complex128})
        if length(basis_l[1])*length(basis_l[2]) != size(data, 1) || length(basis_r[1])*length(basis_r[2]) != size(data, 2)
            throw(DimensionMismatch())
        end
        new(basis_l, basis_r, data)
    end
end

Base.full(a::SparseSuperOperator) = DenseSuperOperator(a.basis_l, a.basis_r, full(a.data))
Base.full(a::DenseSuperOperator) = deepcopy(a)

=={T<:SuperOperator}(a::T, b::T) = (a.basis_l == b.basis_l) && (a.basis_r == b.basis_r) && (a.data == b.data)

operators.check_samebases(a::SuperOperator, b::SuperOperator) = ((a.basis_l!=b.basis_l) || (a.basis_r!=b.basis_r) ? throw(IncompatibleBases()) : nothing)

function *{T<:SuperOperator}(a::T, b::DenseOperator)
    if a.basis_r[1] != b.basis_l || a.basis_r[2] != b.basis_r
        throw(DimensionMismatch())
    end
    data = a.data*reshape(b.data, length(b.data))
    return DenseOperator(a.basis_l[1], a.basis_l[2], reshape(data, length(a.basis_l[1]), length(a.basis_l[2])))
end

function *{T<:SuperOperator}(a::T, b::T)
    if a.basis_r != b.basis_l
        throw(DimensionMismatch())
    end
    return T(a.basis_l, b.basis_r, a.data*b.data)
end

/{T<:SuperOperator}(a::T, b::Number) = T(a.basis_l, a.basis_r, a.data/complex(b))

+{T<:SuperOperator}(a::T, b::T) = (operators.check_samebases(a, b); T(a.basis_l, a.basis_r, a.data+b.data))

-{T<:SuperOperator}(a::T, b::T) = (operators.check_samebases(a, b); T(a.basis_l, a.basis_r, a.data-b.data))
-{T<:SuperOperator}(a::T) = T(a.basis_l, a.basis_r, -a.data)

"""
Create a super-operator equivalent for right side operator multiplication.

For operators :math:`A`, :math:`B` the relation

.. math::

    \\mathrm{spre}(A) B = A B

holds.
"""
spre(op::DenseOperator) = DenseSuperOperator((op.basis_l, op.basis_r), (op.basis_l, op.basis_r), tensor(identity(op), op).data)
spre(op::SparseOperator) = SparseSuperOperator((op.basis_l, op.basis_r), (op.basis_l, op.basis_r), tensor(identity(op), op).data)

"""
Create a super-operator equivalent for left side operator multiplication.

For operators :math:`A`, :math:`B` the relation

.. math::

    \\mathrm{spost}(A) B = B A

holds.
"""
spost(op::DenseOperator) = DenseSuperOperator((op.basis_l, op.basis_r), (op.basis_l, op.basis_r), kron(transpose(op.data), identity(op).data))
spost(op::SparseOperator) = SparseSuperOperator((op.basis_l, op.basis_r), (op.basis_l, op.basis_r),  kron(transpose(op.data), identity(op).data))

function _check_input(H::Operator, J::Vector, Jdagger::Vector, Gamma::Union{Vector{Float64}, Matrix{Float64}})
    for j=J
        @assert typeof(j) <: Operator
        operators.check_samebases(H, j)
    end
    for j=Jdagger
        @assert typeof(j) <: Operator
        operators.check_samebases(H, j)
    end
    @assert length(J)==length(Jdagger)
    if typeof(Gamma) == Matrix{Float64}
        @assert size(Gamma, 1) == size(Gamma, 2) == length(J)
    elseif typeof(Gamma) == Vector{Float64}
        @assert length(Gamma) == length(J)
    else
        error()
    end
end

"""
Create a super-operator equivalent to the master equation.

The  so that :math:`\\dot \\rho = S \\rho`

The super-operator :math:`S` is defined by

.. math::

    \\dot \\rho = S \\rho = -\\frac{i}{\\hbar} [H,\\rho]
            + 2 J \\rho J^\\dagger - J^\\dagger J \\rho - \\rho J^\\dagger J

Arguments
---------
H
    Hamiltonian
J
    Vector of jump operators

Keyword Arguments
-----------------
Gamma (optional)
    Vector or matrix specifying the coefficients for the jump operators.
Jdagger (optional)
    Vector containing the hermitian conjugates of the jump operators. If they
    are not given they are calculated automatically.
"""
function liouvillian{T<:Operator}(H::T, J::Vector{T};
            Gamma::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)),
            Jdagger::Vector{T}=map(dagger, J))
    _check_input(H, J, Jdagger, Gamma)
    L = spre(-1im*H) + spost(1im*H)
    if typeof(Gamma) == Matrix{Float64}
        for i=1:length(J), j=1:length(J)
            jdagger_j = Gamma[i,j]/2*Jdagger[j]*J[i]
            L -= spre(jdagger_j) + spost(jdagger_j)
            L += spre(Gamma[i,j]*J[i]) * spost(Jdagger[j])
        end
    elseif typeof(Gamma) == Vector{Float64}
        for i=1:length(J)
            jdagger_j = Gamma[i]/2*Jdagger[i]*J[i]
            L -= spre(jdagger_j) + spost(jdagger_j)
            L += spre(Gamma[i]*J[i]) * spost(Jdagger[i])
        end
    else
        error()
    end
    return L
end

Base.expm(op::DenseSuperOperator) = DenseSuperOperator(op.basis_l, op.basis_r, expm(op.data))

end # module


text:
using Documenter
using SnakeIterator
using Test

@testset "SnakeIterator.jl" begin
    doctest(SnakeIterator; manual = false)
end