text:
"""
    GridTopology(mesh)

A generic tensor-product topology defined on either an equispaced, regular rectangular
structured grid or irregular one.
"""
struct GridTopology{M <: AbstractMesh, B} <: AbstractTopology
    mesh::M
    boundaries::B
end
function GridTopology(mesh::M) where {M <: AbstractMesh}
    x1boundary = mesh.domain.x1boundary
    x2boundary = mesh.domain.x2boundary
    boundaries = if isnothing(x1boundary)
        if isnothing(x2boundary)
            NamedTuple()
        else
            NamedTuple{x2boundary}((3, 4))
        end
    else
        if isnothing(x2boundary)
            NamedTuple{x1boundary}((1, 2))
        else
            NamedTuple{(x1boundary..., x2boundary...)}((1, 2, 3, 4))
        end
    end
    GridTopology(mesh, boundaries)
end

function Base.show(io::IO, topology::GridTopology{M}) where {M <: AbstractMesh}
    print(io, "GridTopology on ", topology.mesh)
end
domain(topology::GridTopology{M}) where {M <: AbstractMesh} =
    topology.mesh.domain

function nlocalelems(topology::GridTopology{M}) where {M <: AbstractMesh}
    n1 = topology.mesh.n1
    n2 = topology.mesh.n2
    return n1 * n2
end

eachslabindex(topology::GridTopology{M}) where {M <: AbstractMesh} =
    1:nlocalelems(topology)

# InteriorFaceIterator
function Base.length(
    fiter::InteriorFaceIterator{T},
) where {M <: AbstractMesh, T <: GridTopology{M}}
    topology = fiter.topology
    mesh = topology.mesh
    n1 = mesh.n1
    n2 = mesh.n2
    x1periodic = isnothing(mesh.domain.x1boundary)
    x2periodic = isnothing(mesh.domain.x2boundary)
    return (x1periodic ? n1 : n1 - 1) * n2 + n1 * (x2periodic ? n2 : n2 - 1)
end

function Base.iterate(
    fiter::InteriorFaceIterator{T},
    (d, z1, z2) = (1, 0, 0),
) where {M <: AbstractMesh, T <: GridTopology{M}}
    # iteration state (major first)
    #  - d ∈ (1,2): face direction
    #  - z1 ∈ 0:n1-1: 0-based face index in direction 1
    #  - z2 ∈ 0:n2-1: 0-based face index in direction 2

    topology = fiter.topology
    mesh = topology.mesh
    n1 = mesh.n1
    n2 = mesh.n2
    x1periodic = isnothing(mesh.domain.x1boundary)
    x2periodic = isnothing(mesh.domain.x2boundary)

    # skip boundary faces
    if d == 1 && z1 == 0 && !x1periodic
        d = 2
    end
    if d == 2 && z2 == 0 && !x2periodic
        d = 1
        z1 += 1
        if z1 >= n1
            z1 = 0
            z2 += 1
            if !x1periodic
                d = 2
            end
        end
    end

    if z2 >= n2
        return nothing
    end

    if d == 1
        y1 = z1 == 0 ? n1 - 1 : z1 - 1
        y2 = z2
    else
        y1 = z1
        y2 = z2 == 0 ? n2 - 1 : z2 - 1
    end

    elem1 = z2 * n1 + z1 + 1
    elem2 = y2 * n1 + y1 + 1
    if d == 1
        nextstate = (2, z1, z2)
    else
        z1 += 1
        if z1 == n1
            z1 = 0
            z2 += 1
        end
        nextstate = (1, z1, z2)
    end
    if d == 1
        return (elem1, 1, elem2, 2, false), nextstate
    else
        return (elem1, 3, elem2, 4, false), nextstate
    end
end

# BoundaryFaceIterator
function boundary_names(topology::GridTopology{M}) where {M <: AbstractMesh}
    x1boundary = topology.mesh.domain.x1boundary
    x2boundary = topology.mesh.domain.x2boundary
    if isnothing(x1boundary)
        isnothing(x2boundary) ? () : x2boundary
    else
        isnothing(x2boundary) ? x1boundary : (x1boundary..., x2boundary...)
    end
end

function boundary_tag(
    topology::GridTopology{M},
    name::Symbol,
) where {M <: AbstractMesh}
    x1boundary = topology.mesh.domain.x1boundary
    x2boundary = topology.mesh.domain.x2boundary
    if !isnothing(x1boundary)
        x1boundary[1] == name && return 1
        x1boundary[2] == name && return 2
    end
    if !isnothing(x2boundary)
        x2boundary[1] == name && return 3
        x2boundary[2] == name && return 4
    end
    error("Invalid boundary name")
end

function boundaries(topology::GridTopology{M}) where {M <: AbstractMesh}
    return topology.boundaries
end

function Base.length(
    bfiter::BoundaryFaceIterator{T},
) where {M <: AbstractMesh, T <: GridTopology{M}}
    boundary = bfiter.boundary
    topology = bfiter.topology
    if boundary in (1, 2)
        if isnothing(topology.mesh.domain.x1boundary)
            return 0
        else
            return topology.mesh.n2
        end
    end
    if boundary in (3, 4)
        if isnothing(topology.mesh.domain.x2boundary)
            return 0
        else
            return topology.mesh.n1
        end
    end
end

function Base.iterate(
    bfiter::BoundaryFaceIterator{T},
) where {M <: AbstractMesh, T <: GridTopology{M}}
    boundary = bfiter.boundary
    topology = bfiter.topology
    if boundary in (1, 2) && isnothing(topology.mesh.domain.x1boundary)
        return nothing
    end
    if boundary in (3, 4) && isnothing(topology.mesh.domain.x2boundary)
        return nothing
    end
    Base.iterate(bfiter, 0)
end

function Base.iterate(
    bfiter::BoundaryFaceIterator{T},
    z,
) where {M <: AbstractMesh, T <: GridTopology{M}}
    boundary = bfiter.boundary
    topology = bfiter.topology
    mesh = topology.mesh
    n1 = mesh.n1
    n2 = mesh.n2
    if boundary == 1
        z >= n2 && return nothing
        elem = z * n1 + 1
    elseif boundary == 2
        z >= n2 && return nothing
        elem = z * n1 + n1
    elseif boundary == 3
        z >= n1 && return nothing
        elem = z + 1
    elseif boundary == 4
        z >= n1 && return nothing
        elem = (n2 - 1) * n1 + z + 1
    end
    return (elem, boundary), z + 1
end

# VertexIterator
function Base.length(
    viter::VertexIterator{T},
) where {M <: AbstractMesh, T <: GridTopology{M}}
    topology = viter.topology
    mesh = topology.mesh
    n1 = mesh.n1
    n2 = mesh.n2
    x1periodic = isnothing(mesh.domain.x1boundary)
    x2periodic = isnothing(mesh.domain.x2boundary)
    nv1 = x1periodic ? n1 : n1 + 1
    nv2 = x2periodic ? n2 : n2 + 1
    return nv1 * nv2
end

function Base.iterate(
    viter::VertexIterator{T},
    (z1, z2) = (0, 0),
) where {M <: AbstractMesh, T <: GridTopology{M}}
    topology = viter.topology
    mesh = topology.mesh
    n1 = mesh.n1
    n2 = mesh.n2
    x1periodic = isnothing(mesh.domain.x1boundary)
    x2periodic = isnothing(mesh.domain.x2boundary)
    nv1 = x1periodic ? n1 : n1 + 1 # unique vertices in x1 direction
    nv2 = x2periodic ? n2 : n2 + 1 # unique vertices in x2 direction

    if z2 >= nv2
        return nothing
    end
    vertex = Vertex(topology, (z1, z2))
    z1 += 1
    if z1 >= nv1
        nextstate = (0, z2 + 1)
    else
        nextstate = (z1, z2)
    end
    return vertex, nextstate
end

# Vertex
function Base.length(
    vertex::Vertex{T},
) where {M <: AbstractMesh, T <: GridTopology{M}}
    topology = vertex.topology
    mesh = topology.mesh
    n1 = mesh.n1
    n2 = mesh.n2
    x1periodic = isnothing(mesh.domain.x1boundary)
    x2periodic = isnothing(mesh.domain.x2boundary)

    z1, z2 = vertex.num

    k1 = !x1periodic && (z1 == 0 || z1 == n1) ? 1 : 2
    k2 = !x2periodic && (z2 == 0 || z2 == n2) ? 1 : 2
    return k1 * k2
end

function Base.iterate(
    vertex::Vertex{T},
    vert = 0,
) where {M <: AbstractMesh, T <: GridTopology{M}}
    # iterator of (element, vertnum) that share global vertex
    topology = vertex.topology
    mesh = topology.mesh
    n1 = mesh.n1
    n2 = mesh.n2
    x1periodic = isnothing(mesh.domain.x1boundary)
    x2periodic = isnothing(mesh.domain.x2boundary)
    nv1 = x1periodic ? n1 : n1 + 1
    nv2 = x2periodic ? n2 : n2 + 1
    z1, z2 = vertex.num

    vert += 1
    # at the boundary, skip non-existent elements
    if !x1periodic
        if z1 == 0 && (vert == 2 || vert == 4)
            vert += 1
        end
        if z1 == n1 && (vert == 1 || vert == 3)
            vert += 1
        end
    end
    if !x2periodic
        if z2 == 0 && (vert == 3 || vert == 4)
            vert += 2
        end
        if z2 == n2 && (vert == 1 || vert == 2)
            vert += 2
        end
    end

    if vert > 4
        return nothing
    end

    if vert == 2 || vert == 4
        z1 = mod(z1 - 1, nv1)
    end
    if vert == 3 || vert == 4
        z2 = mod(z2 - 1, nv2)
    end
    elem = z2 * n1 + z1 + 1
    return (elem, vert), vert
end

# GridTopology1D
"""
    GridTopology1D(mesh)

A line topology defined on an equispaced linear mesh.
"""

function GridTopology1D(mesh::Meshes.EquispacedLineMesh)
    x1boundary = mesh.domain.x3boundary
    x2boundary = nothing
    boundaries =
        isnothing(x1boundary) ? NamedTuple() : NamedTuple{x1boundary}((1, 2))
    return GridTopology(mesh, boundaries)
end

# Uniform grid implementations, dispatching on EquispacedRectangleMesh
function vertex_coordinates(
    topology::GridTopology{M},
    elem::Integer,
) where {M <: EquispacedRectangleMesh}
    @assert 1 <= elem <= nlocalelems(topology)
    CT = Topologies.coordinate_type(topology)

    # convert to 0-based indices
    mesh = topology.mesh
    n1 = mesh.n1
    n2 = mesh.n2
    range1 = mesh.range1
    range2 = mesh.range2

    z2, z1 = fldmod(elem - 1, n1)

    c1 = CT(range1[z1 + 1], range2[z2 + 1])
    c2 = CT(range1[z1 + 2], range2[z2 + 1])
    c3 = CT(range1[z1 + 1], range2[z2 + 2])
    c4 = CT(range1[z1 + 2], range2[z2 + 2])
    return (c1, c2, c3, c4)
end

function opposing_face(
    topology::GridTopology{M},
    elem::Integer,
    face::Integer,
) where {M <: EquispacedRectangleMesh}
    @assert 1 <= elem <= nlocalelems(topology)
    @assert 1 <= face <= 4

    # convert to 0-based indices
    mesh = topology.mesh
    n1 = mesh.n1
    n2 = mesh.n2
    x1periodic = isnothing(mesh.domain.x1boundary)
    x2periodic = isnothing(mesh.domain.x2boundary)

    z2, z1 = fldmod(elem - 1, n1)
    if face == 1
        z1 -= 1
        if z1 < 0
            if !x1periodic
                return (0, 1, false)
            end
            z1 += n1
        end
        opface = 2
    elseif face == 2
        z1 += 1
        if z1 == n1
            if !x1periodic
                return (0, 2, false)
            end
            z1 -= n1
        end
        opface = 1
    elseif face == 3
        z2 -= 1
        if z2 < 0
            if !x2periodic
                return (0, 3, false)
            end
            z2 += n2
        end
        opface = 4
    elseif face == 4
        z2 += 1
        if z2 == n2
            if !x2periodic
                return (0, 4, false)
            end
            z2 -= n2
        end
        opface = 3
    end
    opelem = z2 * n1 + z1 + 1
    return opelem, opface, false
end

# Non-uniform grid implementations, dispatching on TensorProductMesh

function vertex_coordinates(
    topology::GridTopology{M},
    elem::Integer,
) where {M <: TensorProductMesh}
    @assert 1 <= elem <= nlocalelems(topology)

    # convert to 0-based indices
    mesh = topology.mesh
    n1 = mesh.n1
    n2 = mesh.n2
    coordinates = mesh.coordinates

    z2, z1 = fldmod(elem - 1, n1)

    c1 = coordinates[z1 * (n2 + 1) + (z2 + 1)]
    c2 = coordinates[(z1 + 1) * (n2 + 1) + (z2 + 1)]
    c3 = coordinates[z1 * (n2 + 1) + (z2 + 2)]
    c4 = coordinates[(z1 + 1) * (n2 + 1) + (z2 + 2)]

    return (c1, c2, c3, c4)
end

function opposing_face(
    topology::GridTopology{M},
    elem::Integer,
    face::Integer,
) where {M <: TensorProductMesh}
    @assert 1 <= elem <= nlocalelems(topology)
    @assert 1 <= face <= 4

    return topology.mesh.faces[(elem - 1) * 4 + face][3:5]
end


text:
function _degree_centrality(g::SimpleGraph, gtype::Integer; normalize=true)
   n_v = nv(g)
   c = zeros(n_v)
   for v in 1:n_v
       if gtype == 0    # count both in and out degree if appropriate
           deg = outdegree(g, v) + (typeof(g) == DiGraph? indegree(g, v) : 0.0)
       elseif gtype == 1    # count only in degree
           deg = indegree(g, v)
       else                 # count only out degree
           deg = outdegree(g, v)
       end
       s = normalize? (1.0 / (n_v - 1.0)) : 1.0
       c[v] = deg*s
   end
   return c
end

# TODO avoid repetition of this docstring
"""Calculates the [degree centrality](https://en.wikipedia.org/wiki/Centrality#Degree_centrality)
of the graph `g`, with optional (default) normalization."""
degree_centrality(g::SimpleGraph; all...) = _degree_centrality(g, 0; all...)
"""Calculates the [degree centrality](https://en.wikipedia.org/wiki/Centrality#Degree_centrality)
of the graph `g`, with optional (default) normalization."""
indegree_centrality(g::SimpleGraph; all...) = _degree_centrality(g, 1; all...)
"""Calculates the [degree centrality](https://en.wikipedia.org/wiki/Centrality#Degree_centrality)
of the graph `g`, with optional (default) normalization."""
outdegree_centrality(g::SimpleGraph; all...) = _degree_centrality(g, 2; all...)

ext: jl


text:
module AugmentorAddons

using Augmentor
using Augmentor: ImageOperation
using ImageFiltering
using Random: AbstractRNG, GLOBAL_RNG

import Augmentor.supports_eager
import Augmentor.supports_affineview
import Augmentor.supports_view
import Augmentor.supports_stepview
import Augmentor.applystepview
import Augmentor.applyeager

export PadBorders, RCutOut

include("padborders.jl")
include("rcutout.jl")

end


text:
#! /Users/malcolm/bin/julia -q
#
using ArgParse

const ECHO_PORT = 3000
const ECHO_HOST = "localhost"

function parse_commandline()
  s = ArgParseSettings()
  @add_arg_table s begin
  "--server", "-s"
    help = "hostname of the echo server"
    default = ECHO_HOST
  "--port", "-p"
    help = "port number running the service"
    arg_type = Int
    default = ECHO_PORT
  "--quiet", "-q"
    help = "run quietly, i.e. with no server output"
    action = :store_true
  end
 return parse_args(s)
end

pa = parse_commandline()

ehost = pa["server"]
eport = pa["port"]
vflag = !pa["quiet"]

pp =  (ehost == "localhost" ? "" : "$ehost>")

if vflag println("Listening on port $eport") end
server = listen(eport)
while true
  conn = accept(server)
  @async begin
    try
      while true
        s0 = readline(conn)
        s1 = chomp(s0)
        if length(chomp(s1)) > 0
          s2 = reverse(s1)
          if s2 == "."
            println("Done.")
            close(conn)
            exit(0)
          else
            write(conn,string(pp,s2,"\r\n"))
          end
        end
      end
    catch err
      println("Connection lost:  $err")
      exit(1)
    end
  end
end

max_stars_repo_path: Module 3/Chapter09/echos.jl

max_issues_repo_path: Module 3/Chapter09/echos.jl

max_forks_repo_path: Module 3/Chapter09/echos.jl


text:
@testset "readTopology Hybrid Node isMajor Tests" begin
    global n1, n2, n3
	# n1: see issue #44, both hybrid parent edges were minor
	n1 = readTopology("(A,((B,#H1),(C,(D)#H1:::0.2)));")
	@test writeTopology(n1) == "(A,((B,(D)#H1:::0.8),(C,#H1:::0.2)));"
	hybridParents = [x for x in n1.hybrid[1].edge if x.hybrid]
	# Exclusive or to check exactly one parent edge is major
	@test hybridParents[1].isMajor ⊻ hybridParents[2].isMajor
	# Gammas must sum to 1 and be non-negative for hybrids
	@test hybridParents[1].gamma + hybridParents[2].gamma == 1.0
	@test hybridParents[1].gamma >= 0.0
	@test hybridParents[2].gamma >= 0.0
	# n2 previously failed: both hybrid parent edges were major
	n2 = readTopology("(A,((B,#H1:::0.5),(C,(D)#H1:::0.5)));")
	@test writeTopology(n2) == "(A,((B,#H1:::0.5),(C,(D)#H1:::0.5)));"
	hybridParents = [x for x in n2.hybrid[1].edge if x.hybrid]
	# Exclusive or to check exactly one parent edge is major
	@test hybridParents[1].isMajor ⊻ hybridParents[2].isMajor
	# Gammas must sum to 1 and be non-negative for hybrids
	@test hybridParents[1].gamma + hybridParents[2].gamma == 1.0
	@test hybridParents[1].gamma >= 0.0
	@test hybridParents[2].gamma >= 0.0
	n3 = readTopology("(A,((C,(D)#H1:::0.5),(B,#H1:::0.5)));")
	@test writeTopology(n3) =="(A,((C,(D)#H1:::0.5),(B,#H1:::0.5)));"
end

@testset "parsing extended newick" begin
# second reticulation has minor before major, issue #70:
net1 = (@test_logs readTopology("(((lo,#H3),#H4),((sp)#H3,(mu)#H4));"));
net3 = (@test_logs readTopology("(((sp)#H3,(mu)#H4),((lo,#H3),#H4));"));
@test hardwiredClusterDistance(net1, net3, true) == 0
@test_logs readTopology("((((((((((((((Ae_ca1,Ae_ca2),Ae_ca3))#H1,#H2),(((Ae_um1,Ae_um2),Ae_um3),#H1)),((Ae_co1,Ae_co),(((Ae_un1,Ae_un2),Ae_un3),Ae_un4))),(((Ae_ta1,Ae_ta2),(Ae_ta3,Ae_ta4)),((((((((Ae_lo1,Ae_lo2),Ae_lo3),(Ae_sh1,Ae_sh2)),((Ae_bi1,Ae_bi2),Ae_bi3)),((Ae_se1,Ae_se2),Ae_se3)))#H2,#H3),#H4))),(((T_bo1,(T_bo2,T_bo3)),T_bo4),((T_ur1,T_ur2),(T_ur3,T_ur4)))),(((((Ae_sp1,Ae_sp2),Ae_sp3),Ae_sp4))#H3,((((Ae_mu1,Ae_mu2),Ae_mu3),Ae_mu4))#H4))),Ta_ca),S_va),Er_bo),H_vu);");
@test_throws ErrorException readTopology("(((lo,#H3),#H4);") # "Tree ended while reading in subtree ..."
@test_throws ErrorException readTopology("((lo,#H3),);") # "Expected beginning of subtree but read ..."
end


text:
"""
    test_group(
        G,
        g_pts::AbstractVector,
        v_pts::AbstractVector = [],
        ve_pts::AbstractVector = [];
        atol = 1e-10,
        test_mutating = true,
        test_group_exp_log = true,
        test_diff = false,
        test_invariance = false,
        diff_convs = [(), (LeftAction(),), (RightAction(),)],
    )

Tests general properties of the group `G`, given at least three different points
elements of it (contained in `g_pts`).
Optionally, specify `test_diff` to test differentials of translation, using `v_pts`, which
must contain at least one tangent vector at `g_pts[1]`, and the direction conventions
specified in `diff_convs`.
If the group is equipped with an invariant metric, `test_invariance` indicates that the
invariance should be checked for the provided points.
"""
function test_group(
    G,
    g_pts::AbstractVector,
    v_pts::AbstractVector=[],
    ve_pts::AbstractVector=[];
    atol=1e-10,
    test_mutating=true,
    test_group_exp_log=true,
    test_diff=false,
    test_invariance=false,
    diff_convs=[(), (LeftAction(),), (RightAction(),)],
)
    e = make_identity(G, g_pts[1])

    Test.@testset "Basic group properties" begin
        Test.@testset "Closed" begin
            for g1 in g_pts, g2 in g_pts
                g3 = compose(G, g1, g2)
                Test.@test is_point(G, g3, true; atol=atol)
            end
        end

        Test.@testset "Associative" begin
            g12_3 = compose(G, compose(G, g_pts[1], g_pts[2]), g_pts[3])
            g1_23 = compose(G, g_pts[1], compose(G, g_pts[2], g_pts[3]))
            Test.@test isapprox(G, g12_3, g1_23; atol=atol)

            test_mutating && Test.@testset "mutating" begin
                g12, g23, g12_3, g1_23 = allocate.(repeat([g_pts[1]], 4))
                Test.@test compose!(G, g12, g_pts[1], g_pts[2]) === g12
                Test.@test compose!(G, g23, g_pts[2], g_pts[3]) === g23
                Test.@test compose!(G, g12_3, g12, g_pts[3]) === g12_3
                Test.@test compose!(G, g1_23, g_pts[1], g23) === g1_23
                Test.@test isapprox(G, g12_3, g1_23; atol=atol)
            end
        end

        Test.@testset "Identity" begin
            Test.@test isapprox(G, e, e)
            Test.@test identity(G, e) === e
            Test.@test compose(G, e, e) === e
            Test.@test copyto!(e, e) === e

            for g in g_pts
                Test.@test isapprox(G, compose(G, g, e), g)
                Test.@test isapprox(G, compose(G, e, g), g)

                ge = identity(G, g)
                Test.@test isapprox(G, compose(G, g, ge), g)
                Test.@test isapprox(G, compose(G, ge, g), g)
            end

            test_mutating && Test.@testset "mutating" begin
                for g in g_pts
                    h = allocate(g)
                    Test.@test compose!(G, h, g, e) === h
                    Test.@test isapprox(G, h, g)
                    h = allocate(g)
                    Test.@test compose!(G, h, e, g) === h
                    Test.@test isapprox(G, h, g)

                    ge = allocate(g)
                    Test.@test identity!(G, ge, e) === ge
                    Test.@test isapprox(G, compose(G, g, ge), g)
                    Test.@test isapprox(G, compose(G, ge, g), g)

                    ge = allocate(g)
                    Test.@test compose!(G, ge, e, e) === ge
                    Test.@test isapprox(G, ge, e)
                end
            end
        end

        Test.@testset "Inverse" begin
            for g in g_pts
                ginv = inv(G, g)
                Test.@test isapprox(G, compose(G, g, ginv), e; atol=atol)
                Test.@test isapprox(G, compose(G, ginv, g), e; atol=atol)
                Test.@test isapprox(G, e, compose(G, g, ginv); atol=atol)
                Test.@test isapprox(G, e, compose(G, ginv, g); atol=atol)
                Test.@test inv(G, e) === e

                test_mutating && Test.@testset "mutating" begin
                    ginv = allocate(g)
                    Test.@test inv!(G, ginv, g) === ginv
                    Test.@test isapprox(G, compose(G, g, ginv), e; atol=atol)
                    Test.@test isapprox(G, compose(G, ginv, g), e; atol=atol)

                    Test.@test inv(G, e) === e
                    geinv = allocate(g)
                    Test.@test inv!(G, geinv, e) === geinv
                    Test.@test isapprox(G, geinv, e; atol=atol)
                end
            end
        end
    end

    Test.@testset "translation" begin
        convs = ((), (LeftAction(),), (RightAction(),))

        Test.@test isapprox(
            G,
            translate(G, g_pts[1], g_pts[2]),
            compose(G, g_pts[1], g_pts[2]);
            atol=atol,
        )
        Test.@test isapprox(
            G,
            translate(G, g_pts[1], g_pts[2], LeftAction()),
            compose(G, g_pts[1], g_pts[2]);
            atol=atol,
        )
        Test.@test isapprox(
            G,
            translate(G, g_pts[1], g_pts[2], RightAction()),
            compose(G, g_pts[2], g_pts[1]);
            atol=atol,
        )

        for conv in convs
            Test.@test isapprox(
                G,
                inverse_translate(
                    G,
                    g_pts[1],
                    translate(G, g_pts[1], g_pts[2], conv...),
                    conv...,
                ),
                g_pts[2];
                atol=atol,
            )
            Test.@test isapprox(
                G,
                translate(
                    G,
                    g_pts[1],
                    inverse_translate(G, g_pts[1], g_pts[2], conv...),
                    conv...,
                ),
                g_pts[2];
                atol=atol,
            )
        end

        test_mutating && Test.@testset "mutating" begin
            for conv in convs
                g = allocate(g_pts[1])
                Test.@test translate!(G, g, g_pts[1], g_pts[2], conv...) === g
                Test.@test isapprox(
                    G,
                    g,
                    translate(G, g_pts[1], g_pts[2], conv...);
                    atol=atol,
                )

                g = translate(G, g_pts[1], g_pts[2], conv...)
                g2 = allocate(g)
                Test.@test inverse_translate!(G, g2, g_pts[1], g, conv...) === g2
                Test.@test isapprox(G, g2, g_pts[2]; atol=atol)

                g = inverse_translate(G, g_pts[1], g_pts[2], conv...)
                g2 = allocate(g)
                Test.@test translate!(G, g2, g_pts[1], g, conv...) === g2
                Test.@test isapprox(G, g2, g_pts[2]; atol=atol)
            end
        end
    end

    test_diff && Test.@testset "translation differential" begin
        X = v_pts[1]
        g21 = compose(G, g_pts[2], g_pts[1])
        g12 = compose(G, g_pts[1], g_pts[2])
        Test.@test isapprox(
            G,
            g12,
            translate_diff(G, g_pts[2], g_pts[1], X),
            translate_diff(G, g_pts[2], g_pts[1], X, LeftAction());
            atol=atol,
        )
        Test.@test is_vector(
            G,
            g12,
            translate_diff(G, g_pts[2], g_pts[1], X, LeftAction()),
            true;
            atol=atol,
        )
        RightAction() in diff_convs && Test.@test is_vector(
            G,
            g21,
            translate_diff(G, g_pts[2], g_pts[1], X, RightAction()),
            true;
            atol=atol,
        )

        for conv in diff_convs
            g2g1 = translate(G, g_pts[2], g_pts[1], conv...)
            g2invg1 = inverse_translate(G, g_pts[2], g_pts[1], conv...)
            Test.@test isapprox(
                G,
                g_pts[1],
                inverse_translate_diff(
                    G,
                    g_pts[2],
                    g2g1,
                    translate_diff(G, g_pts[2], g_pts[1], X, conv...),
                    conv...,
                ),
                X;
                atol=atol,
            )
            Test.@test isapprox(
                G,
                g_pts[1],
                translate_diff(
                    G,
                    g_pts[2],
                    g2invg1,
                    inverse_translate_diff(G, g_pts[2], g_pts[1], X, conv...),
                    conv...,
                ),
                X;
                atol=atol,
            )
            Xe = inverse_translate_diff(G, g_pts[1], g_pts[1], X, conv...)
            Test.@test isapprox(G, e, Xe, translate_diff(G, e, e, Xe, conv...); atol=atol)
            Test.@test isapprox(
                G,
                e,
                Xe,
                inverse_translate_diff(G, e, e, Xe, conv...);
                atol=atol,
            )
        end

        test_mutating && Test.@testset "mutating" begin
            for conv in diff_convs
                g2g1 = translate(G, g_pts[2], g_pts[1], conv...)
                g2invg1 = inverse_translate(G, g_pts[2], g_pts[1], conv...)
                Y = allocate(X)
                Test.@test translate_diff!(G, Y, g_pts[2], g_pts[1], X, conv...) === Y
                Test.@test isapprox(
                    G,
                    g2g1,
                    Y,
                    translate_diff(G, g_pts[2], g_pts[1], X, conv...);
                    atol=atol,
                )

                Y = translate_diff(G, g_pts[2], g_pts[1], X, conv...)
                Z = allocate(Y)
                Test.@test inverse_translate_diff!(G, Z, g_pts[2], g2g1, Y, conv...) === Z
                Test.@test isapprox(G, g_pts[1], Z, X; atol=atol)

                Y = inverse_translate_diff(G, g_pts[2], g_pts[1], X, conv...)
                Z = allocate(Y)
                Test.@test translate_diff!(G, Z, g_pts[2], g2invg1, Y, conv...) === Z
                Test.@test isapprox(G, g_pts[1], Z, X; atol=atol)
            end
        end
    end

    test_group_exp_log && Test.@testset "group exp/log properties" begin
        Test.@testset "e = exp(0)" begin
            v = group_log(G, identity(G, g_pts[1]))
            g = group_exp(G, v)
            Test.@test isapprox(G, make_identity(G, g_pts[1]), g; atol=atol)

            test_mutating && Test.@testset "mutating" begin
                v = allocate(ve_pts[1])
                Test.@test group_log!(G, v, identity(G, g_pts[1])) === v
                g = allocate(g_pts[1])
                Test.@test group_exp!(G, g, v) === g
                Test.@test isapprox(G, make_identity(G, g_pts[1]), g; atol=atol)
            end
        end

        Test.@testset "v = log(exp(v))" begin
            for v in ve_pts
                g = group_exp(G, v)
                Test.@test is_point(G, g; atol=atol)
                v2 = group_log(G, g)
                Test.@test isapprox(G, make_identity(G, g_pts[1]), v2, v; atol=atol)
            end

            test_mutating && Test.@testset "mutating" begin
                for v in ve_pts
                    g = allocate(g_pts[1])
                    Test.@test group_exp!(G, g, v) === g
                    Test.@test is_point(G, g; atol=atol)
                    Test.@test isapprox(G, g, group_exp(G, v); atol=atol)
                    v2 = allocate(v)
                    Test.@test group_log!(G, v2, g) === v2
                    Test.@test isapprox(G, make_identity(G, g_pts[1]), v2, v; atol=atol)
                end
            end
        end

        Test.@testset "inv(g) = exp(-log(g))" begin
            g = g_pts[1]
            v = group_log(G, g)
            ginv = group_exp(G, -v)
            Test.@test isapprox(G, ginv, inv(G, g); atol=atol)
        end

        Test.@testset "exp(sv)∘exp(tv) = exp((s+t)v)" begin
            g1 = group_exp(G, 0.2 * ve_pts[1])
            g2 = group_exp(G, 0.3 * ve_pts[1])
            g12 = group_exp(G, 0.5 * ve_pts[1])
            g1_g2 = compose(G, g1, g2)
            g2_g1 = compose(G, g2, g1)
            isapprox(G, g1_g2, g12; atol=atol)
            isapprox(G, g2_g1, g12; atol=atol)
        end
    end

    test_group_exp_log &&
        test_diff &&
        Test.@testset "exp/log retract/inverse_retract" begin
            for conv in diff_convs
                y = retract(
                    G,
                    g_pts[1],
                    v_pts[1],
                    Manifolds.GroupExponentialRetraction(conv...),
                )
                Test.@test is_point(G, y; atol=atol)
                v2 = inverse_retract(
                    G,
                    g_pts[1],
                    y,
                    Manifolds.GroupLogarithmicInverseRetraction(conv...),
                )
                Test.@test isapprox(G, g_pts[1], v2, v_pts[1]; atol=atol)
            end

            test_mutating && Test.@testset "mutating" begin
                for conv in diff_convs
                    y = allocate(g_pts[1])
                    Test.@test retract!(
                        G,
                        y,
                        g_pts[1],
                        v_pts[1],
                        Manifolds.GroupExponentialRetraction(conv...),
                    ) === y
                    Test.@test is_point(G, y; atol=atol)
                    v2 = allocate(v_pts[1])
                    Test.@test inverse_retract!(
                        G,
                        v2,
                        g_pts[1],
                        y,
                        Manifolds.GroupLogarithmicInverseRetraction(conv...),
                    ) === v2
                    Test.@test isapprox(G, g_pts[1], v2, v_pts[1]; atol=atol)
                end
            end
        end

    test_invariance && Test.@testset "metric invariance" begin
        if has_invariant_metric(G, LeftAction())
            Test.@testset "left-invariant" begin
                Test.@test has_approx_invariant_metric(
                    G,
                    g_pts[1],
                    v_pts[1],
                    v_pts[end],
                    g_pts,
                    LeftAction(),
                )
            end
        end
        if invariant_metric_dispatch(G, RightAction()) === Val(true)
            Test.@testset "right-invariant" begin
                Test.@test has_approx_invariant_metric(
                    G,
                    g_pts[1],
                    v_pts[1],
                    v_pts[end],
                    g_pts,
                    RightAction(),
                )
            end
        end
    end
    return nothing
end

"""
    test_action(
        A::AbstractGroupAction,
        a_pts::AbstractVector,
        m_pts::AbstractVector,
        v_pts = [];
        atol = 1e-10,
        atol_ident_compose = 0,
        test_optimal_alignment = false,
        test_mutating = true,
        test_diff = false,
        test_switch_direction = true,
    )

Tests general properties of the action `A`, given at least three different points
that lie on it (contained in `a_pts`) and three different point that lie
on the manifold it acts upon (contained in `m_pts`).

# Arguments
- `atol_ident_compose = 0`: absolute tolerance for the test that composition with identity
  doesn't change the group element.
"""
function test_action(
    A::AbstractGroupAction,
    a_pts::AbstractVector,
    m_pts::AbstractVector,
    v_pts=[];
    atol=1e-10,
    atol_ident_compose=0,
    test_optimal_alignment=false,
    test_mutating=true,
    test_diff=false,
    test_switch_direction=true,
)
    G = base_group(A)
    M = g_manifold(A)
    e = make_identity(G, a_pts[1])

    Test.@testset "Basic action properties" begin
        test_switch_direction && Test.@testset "Direction" begin
            Aswitch = switch_direction(A)
            if isa(A, AbstractGroupAction{LeftAction})
                Test.@test direction(A) === LeftAction()
                Test.@test isa(Aswitch, AbstractGroupAction{RightAction})
                Test.@test direction(Aswitch) === RightAction()
            else
                Test.@test direction(A) === RightAction()
                Test.@test isa(Aswitch, AbstractGroupAction{LeftAction})
                Test.@test direction(Aswitch) === LeftAction()
            end
        end

        Test.@testset "Closed" begin
            Test.@testset "over actions" begin
                for a1 in a_pts, a2 in a_pts
                    a3 = compose(A, a1, a2)
                    Test.@test is_point(G, a3, true; atol=atol)
                end
            end
            Test.@testset "over g-manifold" begin
                for a in a_pts, m in m_pts
                    Test.@test is_point(M, apply(A, a, m), true; atol=atol)
                    Test.@test is_point(M, inverse_apply(A, a, m), true; atol=atol)
                end
            end
        end

        Test.@testset "Associative" begin
            a12 = compose(A, a_pts[1], a_pts[2])
            a23 = compose(A, a_pts[2], a_pts[3])

            Test.@testset "over compose" begin
                a12_a3 = compose(A, a12, a_pts[3])
                a1_a23 = compose(A, a_pts[1], a23)
                Test.@test isapprox(G, a12_a3, a1_a23; atol=atol)
            end

            Test.@testset "over apply" begin
                for m in m_pts
                    a12_a3_m = apply(A, a12, apply(A, a_pts[3], m))
                    a1_a23_m = apply(A, a_pts[1], apply(A, a23, m))
                    Test.@test isapprox(M, a12_a3_m, a1_a23_m; atol=atol)
                end
            end

            test_mutating && Test.@testset "mutating" begin
                a12, a23, a12_3, a1_23 = allocate.(repeat([a_pts[1]], 4))
                Test.@test compose!(A, a12, a_pts[1], a_pts[2]) === a12
                Test.@test compose!(A, a23, a_pts[2], a_pts[3]) === a23
                Test.@test compose!(A, a12_3, a12, a_pts[3]) === a12_3
                Test.@test compose!(A, a1_23, a_pts[1], a23) === a1_23
                Test.@test isapprox(G, a12_3, a1_23; atol=atol)

                for m in m_pts
                    a12_a3_m, a1_a23_m = allocate(m), allocate(m)
                    Test.@test apply!(A, a12_a3_m, a12, apply(A, a_pts[3], m)) === a12_a3_m
                    Test.@test apply!(A, a1_a23_m, a_pts[1], apply(A, a23, m)) === a1_a23_m
                    Test.@test isapprox(M, a12_a3_m, a1_a23_m; atol=atol)
                end
            end
        end

        Test.@testset "Identity" begin
            Test.@test compose(A, e, e) === e

            for a in a_pts
                Test.@test isapprox(G, compose(A, a, e), a; atol=atol_ident_compose)
                Test.@test isapprox(G, compose(A, e, a), a; atol=atol_ident_compose)

                ge = identity(G, a)
                Test.@test isapprox(G, compose(A, a, ge), a; atol=atol_ident_compose)
                Test.@test isapprox(G, compose(A, ge, a), a; atol=atol_ident_compose)

                for m in m_pts
                    Test.@test isapprox(M, apply(A, e, m), m)
                    Test.@test isapprox(M, apply(A, ge, m), m)
                    Test.@test isapprox(M, inverse_apply(A, e, m), m)
                    Test.@test isapprox(M, inverse_apply(A, ge, m), m)
                end
            end

            test_mutating && Test.@testset "mutating" begin
                for a in a_pts
                    h = allocate(a)
                    Test.@test compose!(A, h, a, e) === h
                    Test.@test isapprox(G, h, a)
                    h = allocate(a)
                    Test.@test compose!(A, h, e, a) === h
                    Test.@test isapprox(G, h, a)

                    ge = identity(G, a)
                    Test.@test isapprox(G, compose(A, a, ge), a)
                    Test.@test isapprox(G, compose(A, ge, a), a)

                    ge = identity(G, a)
                    Test.@test compose!(A, ge, e, e) === ge
                    Test.@test isapprox(G, ge, e)

                    for m in m_pts
                        em = allocate(m)
                        Test.@test apply!(A, em, e, m) === em
                        Test.@test isapprox(M, em, m)
                        em = allocate(m)
                        Test.@test apply!(A, em, ge, m) === em
                        Test.@test isapprox(M, em, m)
                        em = allocate(m)
                        Test.@test inverse_apply!(A, em, e, m) === em
                        Test.@test isapprox(M, em, m)
                        em = allocate(m)
                        Test.@test inverse_apply!(A, em, ge, m) === em
                        Test.@test isapprox(M, em, m)
                    end
                end
            end
        end

        Test.@testset "Inverse" begin
            for a in a_pts
                ainv = inv(G, a)
                Test.@test isapprox(G, compose(A, a, ainv), e; atol=atol)
                Test.@test isapprox(G, compose(A, ainv, a), e; atol=atol)
                Test.@test isapprox(G, e, compose(A, a, ainv); atol=atol)
                Test.@test isapprox(G, e, compose(A, ainv, a); atol=atol)

                for m in m_pts
                    Test.@test isapprox(
                        M,
                        apply(A, a, m),
                        inverse_apply(A, ainv, m);
                        atol=atol,
                    )
                    Test.@test isapprox(
                        M,
                        apply(A, ainv, m),
                        inverse_apply(A, a, m);
                        atol=atol,
                    )
                end
            end
        end
    end

    test_diff && Test.@testset "apply differential" begin
        for (m, v) in zip(m_pts, v_pts)
            for a in a_pts
                am, av = apply(A, a, m), apply_diff(A, a, m, v)
                ainvm, ainvv = inverse_apply(A, a, m), inverse_apply_diff(A, a, m, v)
                Test.@test is_vector(M, am, av, true; atol=atol)
                Test.@test is_vector(M, ainvm, ainvv, true; atol=atol)
            end

            a12 = compose(A, a_pts[1], a_pts[2])
            a2m = apply(A, a_pts[2], m)
            a12v = apply_diff(A, a12, m, v)
            a2v = apply_diff(A, a_pts[2], m, v)
            Test.@test isapprox(M, a2m, apply_diff(A, a_pts[1], a2m, a2v), a12v; atol=atol)

            Test.@test isapprox(M, m, apply_diff(A, e, m, v), v; atol=atol)
            Test.@test isapprox(M, m, inverse_apply_diff(A, e, m, v), v; atol=atol)
        end

        test_mutating && Test.@testset "mutating" begin
            for (m, v) in zip(m_pts, v_pts)
                for a in a_pts
                    am = apply(A, a, m)
                    av = allocate(v)
                    Test.@test apply_diff!(A, av, a, m, v) === av
                    ainvm = inverse_apply(A, a, m)
                    ainvv = allocate(v)
                    Test.@test inverse_apply_diff!(A, ainvv, a, m, v) === ainvv
                    Test.@test is_vector(M, am, av, true; atol=atol)
                    Test.@test is_vector(M, ainvm, ainvv, true; atol=atol)
                end

                a12 = compose(A, a_pts[1], a_pts[2])
                a2m = apply(A, a_pts[2], m)
                a12m = apply(A, a12, m)
                a12v, a2v, a1_a2v = allocate(v), allocate(v), allocate(v)
                Test.@test apply_diff!(A, a12v, a12, m, v) === a12v
                Test.@test apply_diff!(A, a2v, a_pts[2], m, v) === a2v
                Test.@test apply_diff!(A, a1_a2v, a_pts[1], a2m, a2v) === a1_a2v
                Test.@test isapprox(M, a12m, a1_a2v, a12v; atol=atol)

                ev = allocate(v)
                Test.@test apply_diff!(A, ev, e, m, v) === ev
                Test.@test isapprox(G, m, ev, v; atol=atol)
                ev = allocate(v)
                Test.@test inverse_apply_diff!(A, ev, e, m, v) === ev
                Test.@test isapprox(G, m, ev, v; atol=atol)
            end
        end
    end

    test_optimal_alignment && Test.@testset "Center of orbit" begin
        act = center_of_orbit(A, [m_pts[1]], m_pts[2])
        act2 = center_of_orbit(A, [m_pts[1]], m_pts[2], GradientDescentEstimation())
        act_opt = optimal_alignment(A, m_pts[2], m_pts[1])
        Test.@test isapprox(G, act, act_opt; atol=atol)
        Test.@test isapprox(G, act2, act_opt; atol=atol)

        test_mutating && Test.@testset "mutating" begin
            act_opt2 = allocate(act_opt)
            optimal_alignment!(A, act_opt2, m_pts[2], m_pts[1])
            Test.@test isapprox(G, act_opt, act_opt2; atol=atol)
        end
    end
    return nothing
end


text:
import CuArrays
import Zygote

include("./DataTools.jl")
include("./Transformer.jl")

const RootDir = "./models/AttentionEmbed/"

sentences = vcat(TextToArray(RootDir * "data/plots.txt"), TextToArray(RootDir * "data/quotes.txt"))
cvData, trainData = DatasetSplit(sentences, cvPercent = 0.2)

const NumEpochs = 5
const BatchSize = 32
const NumBatches = 20

const Vocab = Int32(1e3)
const ModelVocab = Vocab + 3
const EmbDims = 64
const TimeSteps = GetMaxSequenceLength(trainData)
const Layers = 8
const Stacks = 1

const Input = SequenceFeed(EmbDims, ModelVocab, TimeSteps)
const Encode = Encoding(EmbDims, TimeSteps, Layers, Stacks)
const Decode = Decoding(EmbDims, TimeSteps, Layers, Stacks)
const Output = DenseProjection(ModelVocab)

const θ = Flux.params(Input, Encode, Decode, Output)
const opt = Flux.ADAM(0.01)

const Lookup = TopNLookup(TokenFrequencies(sentences), Vocab)
const cvX = TokenizeLookup(cvData, Lookup)

function MakeDecodeTarget(x::AbstractArray{String, 1})
    x = Tokenize(x)
    target = TagEnding(x)
    decode = TagBeginning(x)

    target = TokenizeLookup(map(s -> join(s, " "), target), Lookup)
    decode = TokenizeLookup(map(s -> join(s, " "), decode), Lookup)
    return decode, target
end

function Model(x::AbstractArray{T, 1}, y::AbstractArray{T, 1}) where T
    EncSeqLens = SequenceLengths(x)
    x = Input(x)
    x, EncMask, context = Encode(x, EncSeqLens)

    DecSeqLens = SequenceLengths(y)
    y = Input(y)
    y = Decode(x, y, DecSeqLens, context, EncMask, Encode.Attn)
    y = Output(y, Input.Emb.E)
    return y
end

function Loss(x::AbstractArray{T, 1}, xdecode::AbstractArray{T, 1}, y::AbstractArray{T, 1}) where T
    ϵ = Float32(1e-8)
    ŷ = Model(x, xdecode)
    (Steps, Labels, Batches) = size(ŷ)

    y = map(seq -> Float32.(permutedims(Flux.onehotbatch(seq, 1:Labels), [2, 1])) |> Flux.param, y)
    y = ArrayPad(y, Steps)

    s = TensorSoftmax(ŷ, dims=2)
    cost = y .* log.(s .+ ϵ) + (1 .- y) .* log.(1 .- s .+ ϵ)
    return - sum(cost) / prod(size(cost))
end

function callback(x, xdecode, y)
    GC.gc()
    println("\t mini-batch loss: $(Loss(x, xdecode, y))")
end

for epoch in 1:NumEpochs
    println("[INFO] Training on epoch $epoch")
    _sent = Random.shuffle(trainData)[1:(BatchSize * NumBatches)]

    for minibatch in 1:NumBatches
        start = minibatch * BatchSize
        stop = (minibatch + 1) * BatchSize

        if start > length(_sent) || stop > length(_sent)
            break
        end

        sents = _sent[start:stop]
        X = TokenizeLookup(sents, Lookup)
        XDecode, Y = MakeDecodeTarget(sents)
        # Flux.Optimise.train!(Loss, θ, [(X, XDecode, Y)], opt; cb = () -> callback(X, XDecode, Y))

        loss = Loss(X, XDecode, Y)
        # grads = Flux.Tracker.gradient(() -> loss, θ)
        grads = Zygote.gradient(() -> loss, θ)
        Flux.Tracker.update!(opt, θ, grads)
        println("\t mini-batch loss: $loss")
        GC.gc()

        if isnan(loss)
            Ŷ = Model(X, XDecode)
            s = TensorSoftmax(Ŷ, dims=2)
            println("model output check: $(size(Ŷ[isnan.(Ŷ)]))")
            println("nan check: $(size(Ŷ[isnan.(s)]))")
        end
    end
end


text:
# Autogenerated wrapper script for libpng_jll for aarch64-apple-darwin
export libpng

using Zlib_jll
JLLWrappers.@generate_wrapper_header("libpng")
JLLWrappers.@declare_library_product(libpng, "@rpath/libpng16.16.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_library_product(
        libpng,
        "lib/libpng16.16.37.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module ConjugateGradient

function conjgrad(A::Array, b::Array, N::Integer, x1::Array)

    issym(A) || error("Must provide symmetric matrix")

    ∘ = dot                                 # convenience notation

    x = x1
    p = r = b - A*x
    rrold = r∘r

    for i in 1:N
        Ap = A*p
        α = rrold / (p∘Ap)                  # scale correction
        x += α*p                            # add correction
        r -= α*Ap
        rr = r∘r
        √rr < 1e-16 && break                # stop if answer is close
        p *= (rr / rrold)
        p += r
        rrold = rr
    end

    return x

end

# convenience methods
conjgrad(A, b, N::Integer) = conjgrad(A, b, N, ones(size(b)))
conjgrad(A, b, x1::Array) =  conjgrad(A, b, length(b), x1)
conjgrad(A, b) = conjgrad(A, b, length(b), ones(size(b)))

end


text:
#=
Converts requests for files in certain folders and subfolders into responses.
May pick the right folder based on request extension.
Generates index.html on request.
=#
# defines MIMETYPES
include("mimetypes.jl")
"Root of LOCFLDRS"
const FOLDERHOME = replace(joinpath(@__DIR__, "..", "assets") |> realpath, "\\" => "/")
const RGX_FOLDERHOME = Regex("^"* FOLDERHOME)

"Contains the local path to folders to expose. Subfolders will also be exposed."
const LOCFLDRS = [ "/js/", "/html/", "/svg/", "/img/", "/mjs/", "/css/", "/fonts/"]
const IMGEXT = r"jpg|jpeg|png|gif|tiff|tif"
const TINDEX ="""<!DOCTYPE html>
<meta charset="UTF-8"><meta name="viewport" content="width=device-width, initial-scale=1.0">
<meta http-equiv="x-ua-compatible" content="ie=edge">
<title>Index.html</title>"""

"""
Index in String / html format to 1st level in the given folder.
Note, 'Mustache.jl' does these things more elegantly.
"""
function _htmlindex(di::String)
    pubdi = _pubfilename(di)
    if !isdir(di)
        @info "Invalid path:"* di
        error(Base.ArgumentError("Invalid path"))
    end
    shnams =  readdir(di)
    if pubdi *"/" ∈ LOCFLDRS
        # Include links to the other LOCFLDRS."
        foreach(LOCFLDRS) do fol
                              if fol != pubdi *"/"
                                  push!(shnams, ".." * fol[1:end-1])
                              end
                        end
    end
    longnams = di * "/" .* shnams
    timstamps = map(mtime, longnams)
    dtgs = map(timstamps) do ts
        Dates.format(Dates.unix2datetime(ts), Dates.RFC1123Format)
    end
    pubnams = pubdi * "/" .* shnams
    pubnams = map(shnams, longnams) do sn, lna
        isdir(lna) ? pubdi * "/" * sn * "/" : pubdi * "/" * sn
    end
    shownams = map(shnams, longnams) do sn, lna
                            isdir(lna) ? "/" * sn * "/" : sn
                        end
    matr = Array{Any,2}(undef, length(pubnams), 4)
    matr[:, 1] = timstamps
    matr[:, 2] = shownams
    matr[:, 3] = pubnams
    matr[:, 4] = dtgs
    smatr = sortslices(matr, dims =1, rev = true)
    # HTML Table rows with folders /../ in di
    tablcont = join(mapslices(smatr, dims = 2) do rw
        shownam, pubnam, dtg = rw[2], rw[3], rw[4]
        if shownam[1:4] == "/../"
            """
$shownam	$dtg
"""
        else
            ""
        end
    end)
    # HTML Table rows with subordinate folders in di
    tablcont *= join(mapslices(smatr, dims = 2) do rw
                        shownam, pubnam, dtg = rw[2], rw[3], rw[4]
                        if shownam[1] == '/' && shownam[1:4] != "/../"
                            """
$shownam	$dtg
"""
                        else
                            ""
                        end
                    end
                )
    # HTML table rows with files in di
    tablcont *= join(mapslices(smatr, dims = 2) do rw
                        shownam, pubnam, dtg = rw[2], rw[3], rw[4]
                        if shownam[1] != '/'
                            """
$shownam	$dtg
"""
                        else
                            ""
                        end
                    end
                )
    bodyc = "
Index $pubdi

" * tablcont * "

"
    replace(TINDEX, "" => bodyc)
end

_bck2fwdslash(pth) = replace(pth, "\\" => "/")
_doubleslash2slash(pth) = replace(pth, "//" => "/")
_nobadchars(pth) = replace(pth, r"^\s*(#|$<>~!)" => "")
_noschemechars(pth) = replace(pth, r";|:|@|=|&" => "")
_noquery(pth) = replace(pth, "?" => "")
_nojumpup(pth) = replace(pth,r"\.{2,}" => "")
_nokeywords(pth) = replace(pth, r"prn"i=>"")
_maybeindexhtml(pth) = pth =="" || pth[end]=='/' ? pth * "index.html" : pth
_startslash(pth) = startswith(pth, "/") ? pth : "/" * pth
"""
Converts e.g. %20 to space
From URIParser.jl, added: misformed -> ""
"""
function _unescapeuri(str)
    r = UInt8[]
    l = length(str)
    i = 1
    while i <= l
        c = str[i]
        i += 1
        if c == '%'
            try
                c = parse(UInt8, str[i:i+1], base=16)
                i += 2
            catch
                return ""
            end
        end
        push!(r, c)
    end
   return String(r)
end

"
Makes an attempt at safeguarding against malicious requests or mishandled requests before file access.
Does not failsafe against all strings, e.g. passing 'ø' may crash the handler and send no response.
"
_censoredrequest(pubreq::String) = pubreq |>  _bck2fwdslash |> _unescapeuri |> _nobadchars |>
                                   _noschemechars |> _noquery |> lowercase |> _nojumpup |> _nokeywords |>
                                    _maybeindexhtml |> _doubleslash2slash |> _startslash

"""
Returns _locfolder ∈ LOCFLDRS from a requested resource.
If /folder/ is missing, but file extension corresponds to one
of LOCFLDRS,
Defaults to /html/
# Examples
```julia-repl
_locfolder("/cat.svg") -> "/svg/cat.svg"
_locfolder("/cat.unknown") -> "/html/cat.unknown"
_locfolder("/secret/file.txt") -> "/html/"
```
Assumes input is already cleaned of maliciuos requests like "/svg/../../secret/big.txt"
and is lowercase
"""
function _locfolder(creq::String)
    fofirst ="/" * get(split(creq, "/"), 2, "") *"/"
    flast = split(creq, "/")[end]
    fileextension =  splitext(flast)[end][2:end]
    fofilextension = "/" * replace(fileextension, IMGEXT => "img") * "/"
    if flast != "index.html" && fofilextension ∈ LOCFLDRS
        return fofilextension
    end
    if fofirst ∈ LOCFLDRS
        return fofirst
    end
    return "/html/"
end

"
public request - >  (top folder, modified request)
Picks a path ∈ LOCFLDRS, based on
    1) file extension (/img/fig.svg -> /svg/fig.svg)
    2) default  (/X.abc -> /html/X.abc)
"
function _topfolder_modreq(pubreq::String)
    censoredreq = pubreq |> _censoredrequest
    topfolder = _locfolder(censoredreq)
    if startswith(censoredreq, topfolder)
        # Make checks against /html/html/doc.html?
        return topfolder, censoredreq
    end
    # top folder was not taken directly from the request.
    if topfolder != "/html/" && occursin("/html/", censoredreq)
        #= Document link:           "./pic.jpg"
        -> request:                 "/html/pic.jpg"
        -> topfolder, censoredreq:  "/img/", "/html/pic.jpg"
        -> topfolder, modreq:       "/img/", "/img/pic.jpg"
        =#
        return topfolder, replace(censoredreq, "/html/" => topfolder)
    end
    #= Document link:           "./holiday/pic.jpg"
    -> request:                 "/html/holidaypic.jpg"
    -> topfolder, censoredreq:  "/img/", "/html/holidaypic.jpg"
    -> topfolder, modreq:       "/img/", "/img/pic.jpg"
    =#
    return topfolder, topfolder * censoredreq |> _doubleslash2slash
end

"Full path to requested file"
_longfilename(modreq::String) = occursin(RGX_FOLDERHOME, modreq) ? modreq : FOLDERHOME * modreq

"Visible name, for index.html"
_pubfilename(longfnam::String) = replace(longfnam, RGX_FOLDERHOME => "")

"Returns (status{Int}, body{Vector{UInt8} , mime{String}, lastmodified{String}:<day-name>, <day> <month> <year> <hour>:<minute>:<second> GMT
Used for building HTTP.Response"
function _rawresponse(pubreq::String)
    topfolder, modreq = pubreq |>  _topfolder_modreq
    longreq = modreq |> _longfilename
    #@debug("Public request:\t" * pubreq *
    #    "\n\tModified request:\t" * modreq *
    #    "\n\tTop folder:\t" * topfolder *
    #    "\n\tLong request:\t" * longreq)
    if occursin(r"index.htm"i, splitdir(modreq)[2]) && isdir(splitdir(longreq)[1])
        status = 200
        body = _htmlindex(splitdir(longreq)[1])
        tstamp =  time()
        (_, ext) = splitext(modreq)
        mime = get(MIMETYPES, ext[2:end], "application/octet-stream")
    elseif isfile(longreq)
        status = 200
        body = read(longreq)
        tstamp = mtime(longreq)
        (_, ext) = splitext(modreq)
        mime = get(MIMETYPES, ext[2:end], "application/octet-stream")
    else
        status = 404
        body = Vector{UInt8}("Not Found (404):  $modreq could not be found.")
        mime = "text/html"
        tstamp = time()
        @info(modreq * " could not be found. Looked in " * topfolder)
    end
    return status, body, mime, Dates.format(Dates.unix2datetime(tstamp), Dates.RFC1123Format) * " GMT"
end

"To be called by handler. -> Response is the corresponding file with defined mime type info."
function _handle_file(pubreq::WebSockets.Request, resp::WebSockets.Response)
    status, body, mime, tstamp = _rawresponse(pubreq.target)
    if typeof(body) == String
        resp.body = codeunits(body)
    else
        resp.body = body
    end
    WebSockets.setheader(resp, WebSockets.Header("Content-Type", mime))
    WebSockets.setheader(resp, WebSockets.Header("Content-Length",  string(length(body))))
    WebSockets.setheader(resp, WebSockets.Header("Last-Modified", tstamp))
    resp.status = status
    if status != 200
        @info "Failure to handle $(pubreq.target)"
    else
        @debug "Served $(pubreq.target)"
    end
    resp
end
handle_file(pubreq::WebSockets.Request) = _handle_file(pubreq, WebSockets.Response())
nothing

max_stars_repo_name: hustf/EquationEditor

max_issues_repo_name: hustf/EquationEditor

max_forks_repo_name: hustf/EquationEditor


text:
export Closed, Partial, Open, closure

"""
Trait to indicate that a binary operation • is closed over set S. Only methods
of • with the signature •(x::S, y::S) are to be considered.

The definition of closed is that •(x::S, y::S) shall not throw an error, and
•(x::S, y::S) shall return a result of type S.
"""
abstract Closed

"""
Trait to indicate that a binary operation • is a partial function over set S.
Only methods of • with the signature •(x::S, y::S) are to be considered.

The definition of a partial function in this context is that •(x::S, y::S) shall
either throw an error or return a result of type S.
"""
abstract Partial

"""
Trait to indicate that a binary operation • is not known to be closed or to be a
partial function over set S. Note that this trait does not imply the existence
of a counterexample to closure.
"""
abstract Open

"""
Return the closure trait of the set S under the binary operation •.
"""
closure(::Type, ::Function) = Open

# numbers should be closed under + and *
closure{N<:Number}(::Type{N}, ::typeof(+)) = Closed
closure{N<:Number}(::Type{N}, ::typeof(*)) = Closed

# Strings are closed under *
closure(::Type{String}, ::typeof(*)) = Closed


text:
function calc_clamped_piecewise(x::T, x_1::Float64, x_2::Float64, y_1::Float64, y_2::Float64) where {T <: Real}
	# The function varies linearly between y_1 at x_1 to y_2 at x_2.
	# It is assumed that y_2 < y_1.
	# The output is bounded between y_1 and y_2.

	k = (y_2 - y_1) / (x_2 - x_1)
	y = y_1 + (x - x_1) * k
	return clamp(y, y_2, y_1)
end

### Traction
function traction(ins::Regularized, vel_t::SVector{3,T}, p_dA::T) where {T}
	# TODO write unit test to ensure that this function cannot produce a NaN in partials

	mag_vel_t² = dot(vel_t, vel_t)
	μs = ins.μs
	μd = ins.μd
	v_c = ins.v_c

	if mag_vel_t² < (ins.v_c^2)
		T̄c = -μs * vel_t / v_c
	else
		mag_vel_t = sqrt(mag_vel_t²)
		μ = calc_clamped_piecewise(mag_vel_t, ins.v_μs, ins.v_μd, μs, μd)
		T̄c = -μ * vel_t / mag_vel_t
	end

	return T̄c * p_dA
end

function traction(ins::Bristle, T̄s::SVector{3,T}, p_dA::T) where {T}
	# TODO write unit test to ensure that this function cannot produce a NaN in partials

	mag_T̄s² = dot(T̄s, T̄s)
	μs = ins.μs
	μd = ins.μd

	if mag_T̄s² < (μs^2)
		T̄c = T̄s
	else
		mag_T̄s = sqrt(mag_T̄s²)
		μ = calc_clamped_piecewise(mag_T̄s, ins.T̄s_μs, ins.T̄s_μd, μs, μd)
		T̄c = μ * T̄s / mag_T̄s
	end

	return T̄c * p_dA
end

function yes_contact!(fric_type::Regularized, tm::TypedMechanismScenario{T}, c_ins::ContactInstructions) where {T}
    b = tm.bodyBodyCache
    frame = b.mesh_2.FrameID
    twist = b.twist_r²_r¹_r²
    vʳᵉˡ = as_static_vector(twist)
    wrench_lin = zeros(SVector{3,T})
    wrench_ang = zeros(SVector{3,T})
    for k_trac = 1:length(b.TractionCache)
        trac = b.TractionCache[k_trac]
        r = trac.r_cart
        n̂ = trac.n̂
        cart_vel = spatial_vel_formula(vʳᵉˡ, r)
        cart_vel_t = vec_sub_vec_proj(cart_vel, n̂)
        p_dA = calc_p_dA(trac)

		T_c = traction(fric_type, cart_vel_t, p_dA)
        traction_k = p_dA * n̂ + T_c

        wrench_lin += traction_k
        wrench_ang += cross(r, traction_k)
    end
    return Wrench(frame, wrench_ang, wrench_lin)
end

##########################

no_contact!(fric_type::Regularized, tm::TypedMechanismScenario{T}, c_ins::ContactInstructions) where {T} = nothing
function no_contact!(fric_type::Bristle, tm::TypedMechanismScenario{T}, c_ins::ContactInstructions) where {T}
    BF = c_ins.FrictionModel
    bristle_id = BF.BristleID
    get_bristle_d1(tm, bristle_id) .= -(1 / BF.τ) * get_bristle_d0(tm, bristle_id)
end

#########################################################

function calc_K̄_sqrt_inv(s::spatialStiffness{T}) where {T}
    # TODO: make this allocate less

    EF = eigen!(s.K̄)
    s.σ_sqrt.diag .= EF.values
    max_σ = maximum(s.σ_sqrt.diag)
    for k = 1:6
        s.σ_sqrt.diag[k] = 1.0 / sqrt(max(s.σ_sqrt.diag[k], max_σ * 1.0e-16))
    end
    mul!(s.mul_pre, EF.vectors, s.σ_sqrt)
    mul!(s.K̄⁻¹_sqrt, s.mul_pre, EF.vectors')
end

@inline function calc_trace_12(K)
    trace_1 = K[1]  + K[8]  + K[15]
    trace_2 = K[22] + K[29] + K[36]
    return trace_1, trace_2
end

function decompose_K!(s::spatialStiffness{T}, magic::Float64) where {T}
	t_1, t_2 = calc_trace_12(s.K)
	s_1 = 1 / sqrt(t_1)
	for k = 1:3
        s.S⁻¹.diag[k] = s_1 * magic
	end
	s_2 = 1 / sqrt(t_2)
    for k = 4:6
        s.S⁻¹.diag[k] = s_2
	end

    s.K̄.data .= s.S⁻¹ * s.K * s.S⁻¹
    calc_K̄_sqrt_inv(s)
end

function bristle_wrench_in_world(tm::TypedMechanismScenario{T}, c_ins::ContactInstructions) where {T}
    # cop: patch center of pressure

    BF = c_ins.FrictionModel
    bristle_id = BF.BristleID
    b = tm.bodyBodyCache
    spatial_stiffness = b.spatialStiffness
    τ⁻¹ = 1 / BF.τ
    cop, wrench²_normal = normal_wrench_cop(b)
    calc_patch_spatial_stiffness!(tm, BF, cop.v)
	s = get_bristle_d0(tm, bristle_id)
    decompose_K!(spatial_stiffness, BF.magic)
	Δ² = spatial_stiffness.S⁻¹ * (spatial_stiffness.K̄⁻¹_sqrt * s)
	Δ² = SVector{6,T}(Δ²[1], Δ²[2], Δ²[3], Δ²[4], Δ²[5], Δ²[6])
    wrenchᶜᵒᵖ_fric, wrench²_fric = calc_spatial_bristle_force(tm, c_ins, Δ², b.twist_r²_r¹_r², cop.v)
	get_bristle_d1(tm, bristle_id) .= -τ⁻¹ * ( spatial_stiffness.K̄⁻¹_sqrt * (spatial_stiffness.S⁻¹ * wrenchᶜᵒᵖ_fric) + s)
    return wrench²_normal, wrench²_fric
end

#########################################################

function yes_contact!(fric_type::Bristle, tm::TypedMechanismScenario{T}, c_ins::ContactInstructions) where {T}
    wrench²_normal, wrench²_fric = bristle_wrench_in_world(tm, c_ins)
    return wrench²_normal + wrench²_fric
end

spatial_vel_formula(v::SVector{6,T}, b::SVector{3,T}) where {T} = last_3_of_6(v) + cross(first_3_of_6(v), b)

function calc_patch_spatial_stiffness!(tm::TypedMechanismScenario{T}, BF, cop::SVector{3,T}) where {T}
    b = tm.bodyBodyCache
    tc = b.TractionCache
    K_11_sum = zeros(SMatrix{3,3,T,9})
    K_12_sum = zeros(SMatrix{3,3,T,9})
    K_22_sum = zeros(SMatrix{3,3,T,9})
    for k = 1:length(tc)
        trac = tc.vec[k]
        n̂ = trac.n̂
        p_dA = calc_p_dA(trac)
        r = trac.r_cart - cop
        K_22_sum += p_dA * (I - n̂ * n̂')
        r_x_n̂     = cross(r, n̂)
        K_12_sum += p_dA * (vector_to_skew_symmetric(r) - r_x_n̂ * n̂')
        r_skew²   = vector_to_skew_symmetric_squared(r)
        K_11_sum -= p_dA * (r_skew² + r_x_n̂ * r_x_n̂')
    end
    b.spatialStiffness.K.data[1:3, 1:3] .= K_11_sum
    b.spatialStiffness.K.data[4:6, 1:3] .= K_12_sum'
    b.spatialStiffness.K.data[1:3, 4:6] .= K_12_sum
    b.spatialStiffness.K.data[4:6, 4:6] .= K_22_sum
    b.spatialStiffness.K.data .*= BF.k̄
end

function calc_spatial_bristle_force(tm::TypedMechanismScenario{T}, c_ins::ContactInstructions, Δ²::SVector{6,T},
    twist, cop::SVector{3,T}) where {T}

    b = tm.bodyBodyCache
    tc = b.TractionCache
    BF = c_ins.FrictionModel
    τ = BF.τ
    k̄ = BF.k̄
    vʳᵉˡ = as_static_vector(twist)
    wrench_lin = zeros(SVector{3,T})
    wrench_ang = zeros(SVector{3,T})
    for k = 1:length(tc)
        trac = tc.vec[k]
        n̂ = trac.n̂
        r = trac.r_cart
        x² = r - cop
        δ² = spatial_vel_formula(Δ², x²)  # TODO: the stated frame is not correct
        ṙ²_perp = spatial_vel_formula(vʳᵉˡ, r)
        p_dA = calc_p_dA(trac)

		T̄s = -k̄ * (δ² + τ * ṙ²_perp)
        T̄s = vec_sub_vec_proj(T̄s, n̂)
		T_c = traction(BF, T̄s, p_dA)

        wrench_lin += T_c
        wrench_ang += cross(x², T_c)
    end
    wrenchᶜᵒᵖ_fric = vcat(wrench_ang, wrench_lin)
    wrench²_fric = vcat(wrench_ang + cross(cop, wrench_lin), wrench_lin)
    return wrenchᶜᵒᵖ_fric, Wrench(b.mesh_2.FrameID, wrench_ang + cross(cop, wrench_lin), wrench_lin)
end

# function transform_stiffness!(s::spatialStiffness{T}, xform) where {T}
#     R = rotation(xform)
#     t = translation(xform)
#     rx = vector_to_skew_symmetric(t)
#     X = vcat(hcat(R, (rx * R)), hcat(zeros(SMatrix{3,3,T,9}), R))
#     mul!(s.mul_pre, X, s.K)
#     mul!(s.K.data, s.mul_pre, X')
# end

# function calc_s(s::spatialStiffness, Δ)
#     # this function should not be called in tight loops
#     return inv(s.K⁻¹_sqrt) * Δ
# end

# function transform_δ(v::SVector{6,T}, x) where {T}
#     ang = SVector{3,T}(v[1], v[2], v[3])
#     lin = SVector{3,T}(v[4], v[5], v[6])
#     ang, lin = RigidBodyDynamics.Spatial.transform_spatial_motion(ang, lin, rotation(x), translation(x))
#     return vcat(ang, lin)
# end

# function set_s_from_Δʷ(mech_scen, c_ins::ContactInstructions, δʷ)
#     force_single_elastic_intersection!(mech_scen, mech_scen.float, c_ins)
#     δ² = transform_δ(δʷ, mech_scen.float.bodyBodyCache.x_rʷ_r²)
#     Δ = calc_s(mech_scen.float.bodyBodyCache.spatialStiffness, δ²)
#     b_id = c_ins.FrictionModel.BristleID
#     mech_scen.float.s.segments[b_id] .= Δ
# end


text:
using MeshArrays, MITgcmTools, OceanStateEstimation
using CSV, DataFrames, Statistics, Plots
#using FortranFiles,

import Plots: heatmap

"""
    heatmap(x::MeshArray; args...)

Apply heatmap to each subdomain in a MeshArray
"""
function heatmap(x::MeshArray; args...)
    n=x.grid.nFaces
    p=()
    for i=1:n; p=(p...,heatmap(x[i]; args...)); end
    plot(p...)
end

#Convert Velocity (m/s) to transport (m^3/s)
function convert_velocities(U::MeshArray,V::MeshArray,G::NamedTuple)
    for i in eachindex(U)
        tmp1=U[i]; tmp1[(!isfinite).(tmp1)] .= 0.0
        tmp1=V[i]; tmp1[(!isfinite).(tmp1)] .= 0.0
        U[i]=G.DRF[i[2]]*U[i].*G.DYG[i[1]]
        V[i]=G.DRF[i[2]]*V[i].*G.DXG[i[1]]
    end
    return U,V
end

##

"""
    trsp_read(myspec::String,mypath::String)

Function that reads files that were generated by `trsp_prep`
"""
function trsp_read(myspec::String,mypath::String)
    γ=GridSpec(myspec,mypath)
    TrspX=γ.read(mypath*"TrspX.bin",MeshArray(γ,Float32))
    TrspY=γ.read(mypath*"TrspY.bin",MeshArray(γ,Float32))
    TauX=γ.read(mypath*"TauX.bin",MeshArray(γ,Float32))
    TauY=γ.read(mypath*"TauY.bin",MeshArray(γ,Float32))
    SSH=γ.read(mypath*"SSH.bin",MeshArray(γ,Float32))
    return TrspX, TrspY, TauX, TauY, SSH
end

"""
    trsp_prep(γ,Γ,dirOut)

Function that generates small binary files (2D) from large netcdf ones (4D).

```
using FortranFiles, MeshArrays
!isdir("nctiles_climatology") ? error("missing files") : nothing
include(joinpath(dirname(pathof(MeshArrays)),"gcmfaces_nctiles.jl"))
(TrspX, TrspY, TauX, TauY, SSH)=trsp_prep(γ,Γ,MeshArrays.GRID_LLC90);
```
"""
function trsp_prep(γ::gcmgrid,Γ::NamedTuple,dirOut::String="")

    #wind stress
    fileName="nctiles_climatology/oceTAUX/oceTAUX"
    oceTAUX=read_nctiles(fileName,"oceTAUX",γ)
    fileName="nctiles_climatology/oceTAUY/oceTAUY"
    oceTAUY=read_nctiles(fileName,"oceTAUY",γ)
    oceTAUX=mask(oceTAUX,0.0)
    oceTAUY=mask(oceTAUY,0.0)

    #sea surface height anomaly
    fileName="nctiles_climatology/ETAN/ETAN"
    ETAN=read_nctiles(fileName,"ETAN",γ)
    fileName="nctiles_climatology/sIceLoad/sIceLoad"
    sIceLoad=read_nctiles(fileName,"sIceLoad",γ)
    rhoconst=1029.0
    myssh=(ETAN+sIceLoad./rhoconst)
    myssh=mask(myssh,0.0)

    #seawater transports
    fileName="nctiles_climatology/UVELMASS/UVELMASS"
    U=read_nctiles(fileName,"UVELMASS",γ)
    fileName="nctiles_climatology/VVELMASS/VVELMASS"
    V=read_nctiles(fileName,"VVELMASS",γ)
    U=mask(U,0.0)
    V=mask(V,0.0)

    #time averaging and vertical integration
    TrspX=similar(Γ.DXC)
    TrspY=similar(Γ.DYC)
    TauX=similar(Γ.DXC)
    TauY=similar(Γ.DYC)
    SSH=similar(Γ.XC)

    for i=1:γ.nFaces
        tmpX=mean(U.f[i],dims=4)
        tmpY=mean(V.f[i],dims=4)
        for k=1:length(Γ.RC)
            tmpX[:,:,k]=tmpX[:,:,k].*Γ.DYG.f[i]
            tmpX[:,:,k]=tmpX[:,:,k].*Γ.DRF[k]
            tmpY[:,:,k]=tmpY[:,:,k].*Γ.DXG.f[i]
            tmpY[:,:,k]=tmpY[:,:,k].*Γ.DRF[k]
        end
        TrspX.f[i]=dropdims(sum(tmpX,dims=3),dims=(3,4))
        TrspY.f[i]=dropdims(sum(tmpY,dims=3),dims=(3,4))
        TauX.f[i]=dropdims(mean(oceTAUX.f[i],dims=3),dims=3)
        TauY.f[i]=dropdims(mean(oceTAUY.f[i],dims=3),dims=3)
        SSH.f[i]=dropdims(mean(myssh.f[i],dims=3),dims=3)
    end

    if !isempty(dirOut)
        write_bin(TrspX,dirOut*"TrspX.bin")
        write_bin(TrspY,dirOut*"TrspY.bin")
        write_bin(TauX,dirOut*"TauX.bin")
        write_bin(TauY,dirOut*"TauY.bin")
        write_bin(SSH,dirOut*"SSH.bin")
    end

    return TrspX, TrspY, TauX, TauY, SSH
end

"""
    trsp_prep(γ,Γ,dirOut)

Function that writes a `MeshArray` to a binary file using `FortranFiles`.
"""
function write_bin(inFLD::MeshArray,filOut::String)
    recl=prod(inFLD.grid.ioSize)*4
    tmp=Float32.(convert2gcmfaces(inFLD))
    println("saving to file: "*filOut)
    f =  FortranFile(filOut,"w",access="direct",recl=recl,convert="big-endian")
    write(f,rec=1,tmp)
    close(f)
end

##

"""
    rotate_uv(uv,γ)

    1. Convert to `Sv` units and mask out land
    2. Interpolate `x/y` transport to grid cell center
    3. Convert to `Eastward/Northward` transport
    4. Display Subdomain Arrays (optional)
"""
function rotate_uv(uv::Dict,G::NamedTuple)
    u=1e-6 .*uv["U"]; v=1e-6 .*uv["V"];
    u[findall(G.hFacW[:,1].==0)].=NaN
    v[findall(G.hFacS[:,1].==0)].=NaN;

    nanmean(x) = mean(filter(!isnan,x))
    nanmean(x,y) = mapslices(nanmean,x,dims=y)
    (u,v)=exch_UV(u,v); uC=similar(u); vC=similar(v)
    for iF=1:u.grid.nFaces
        tmp1=u[iF][1:end-1,:]; tmp2=u[iF][2:end,:]
        uC[iF]=reshape(nanmean([tmp1[:] tmp2[:]],2),size(tmp1))
        tmp1=v[iF][:,1:end-1]; tmp2=v[iF][:,2:end]
        vC[iF]=reshape(nanmean([tmp1[:] tmp2[:]],2),size(tmp1))
    end

    cs=G.AngleCS
    sn=G.AngleSN
    u=uC.*cs-vC.*sn
    v=uC.*sn+vC.*cs;

    return u,v,uC,vC
end

"""
    interp_uv(u,v)
"""
function interp_uv(u,v)
    mypath=MeshArrays.GRID_LLC90
    SPM,lon,lat=read_SPM(mypath) #interpolation matrix (sparse)
    uI=MatrixInterp(write(u),SPM,size(lon)) #interpolation itself
    vI=MatrixInterp(write(v),SPM,size(lon)); #interpolation itself
    return transpose(uI),transpose(vI),vec(lon[:,1]),vec(lat[1,:])
end

max_stars_repo_name: gaelforget/Notebooks

max_issues_repo_name: gaelforget/Notebooks

max_forks_repo_name: gaelforget/Notebooks


text:
type ClosestDispatch <: DispatchModel
    drivetime::DataFrame
    candidates::Vector{Vector{Int}}
end

function ClosestDispatch(p::DeploymentProblem, drivetime::DataFrame)
    candidates = Array(Vector{Int}, p.nregions)
    I = 1:p.nlocations
    for region in 1:p.nregions
        candidates[region] = I[vec(p.coverage[region,:])]
    end
    ClosestDispatch(drivetime, candidates)
end

function ClosestDispatch(p::DeploymentProblem, problem::DispatchProblem)
    stn_names = [Symbol("stn$(i)_min") for i in 1:size(p.coverage,2)]
    ClosestDispatch(p, problem.emergency_calls[:, stn_names])
end

update_ambulances!(model::ClosestDispatch, i::Int, delta::Int) = nothing

function available_for(dispatch::ClosestDispatch, id::Int, problem::DispatchProblem)
    location = 0
    min_time = Inf
    for i in dispatch.candidates[problem.emergency_calls[id, :neighborhood]]
        if problem.available[i] > 0 && dispatch.drivetime[id, i] < min_time
            location = i
            min_time = dispatch.drivetime[id, i]
        end
    end
    location
end


text:
"""
    BasisFun{T,F} <: ParametricFun

A type representing a linear combination of basis functions.
"""
struct BasisFun{T<:AbstractVector{<:Real},F} <: ParametricFun
    a::T
    basis::F
end

function Base.show(io::IO, ::MIME"text/plain", x::BasisFun{T,F}) where {T,F}
    println(io, "$(length(x.a))-element BasisFun with basis $F and coefficients")
    println(io, x.a)
end

"""
    eval_fun(x, a, basis)

Evaluate a linear combination of basis functions at a point `x`.

## Arguments
* `x`: the point at which to evaluate the function
* `a`: a vector of coefficients representing the function
* `basis`: a function of the form `b(k,x)` that evaluates the `k`-th basis element at `x`.
"""
function eval_fun(x, a, basis)
    y = zero(x * a[1])
    for k in 1:length(a)
        y += a[k] * basis(k, x)
    end
    return y
end

function ChainRulesCore.rrule(::typeof(eval_fun), x, a, basis)
    y = eval_fun(x, a, basis)
    function eval_fun_pullback(δy)
        δf = NoTangent()
        δx = @thunk ForwardDiff.derivative(t -> eval_fun(t, a, basis), x) * δy
        δa = [basis(k,x) for k in 1:length(a)] * δy
        return δf, δx, δa, NoTangent()
    end
    return y, eval_fun_pullback
end

function (f::BasisFun)(x)
    eval_fun(x, f.a, f.basis)
end

function to_grad(f::BasisFun, x)
    [f.basis(k,x) for k in 1:length(f.a)]
end

function from_grad(f::BasisFun, ∇f)
    return BasisFun(∇f, f.basis)
end

params(f::BasisFun) = f.a

Base.:+(f::BasisFun{T,F}, g::BasisFun{T,F}) where {T,F} = BasisFun{T,F}(f.a+g.a, f.basis)
Base.:-(f::BasisFun{T,F}, g::BasisFun{T,F}) where {T,F} = BasisFun{T,F}(f.a-g.a, f.basis)
Base.zero(f::BasisFun) = BasisFun(zero(f.a), f.basis)

Base.length(f::BasisFun) = length(f.a)

Base.:*(α::Number, f::BasisFun) = BasisFun(α * f.a, f.basis)


text:
module DataAnalyses

function start()
    root =  normpath(joinpath(Base.find_package("DataAnalyses"), "..", ".."))
    indexpath = joinpath(root, "index.jl")
    exe = joinpath(Sys.BINDIR, "julia")
    script = :(using Pkg;
                Pkg.activate($root);
                using Pluto;
                Pluto.run(notebook = $indexpath)
    )
    run(`$exe -e $script`)
end

end


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

name = "mimalloc"
version = v"2.0.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/microsoft/mimalloc.git",
              "8e35ccc43be293a9bfd6e63da310a79c235d25d9"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/mimalloc/
if [[ "${target}" == i686-*-mingw* ]]; then
    atomic_patch -p1 ../patches/mimalloc-redirect32.patch
fi
mkdir -p build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DMI_BUILD_STATIC=OFF \
    -DMI_BUILD_OBJECT=OFF \
    -DMI_INSTALL_TOPLEVEL=ON \
    -DMI_BUILD_TESTS=OFF \
    ..
make -j ${nproc}
make -j ${nproc} install
# Manually install the mimalloc-redirect lib for Windows
if [[ "${target}" == i686-*-mingw* ]]; then
    cp "mimalloc-redirect32.${dlext}" "${libdir}/."
elif [[ "${target}" == x86_64-*-mingw* ]]; then
    cp "mimalloc-redirect.${dlext}" "${libdir}/."
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"6.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@testset "define_NumParam" begin

    @test_throws AssertionError HA_solver.define_NumParam(0, 0.015, 0.01,0.0, 0.001)
    @test_throws AssertionError HA_solver.define_NumParam(100, -0.015, 0.01,0.0, 0.001)
    @test_throws AssertionError HA_solver.define_NumParam(100, 0.015, -0.01,0.0, 0.001)
    a_n, a_grid, a_min = HA_solver.define_NumParam(100, 0.015, 0.01,0.0, 0.001)
    @test a_grid == sort(a_grid)
end

@testset "DiscretizeAR" begin

    @test_throws AssertionError HA_solver.DiscretizeAR(-0.1, 0.2, 10, "T")
    @test_throws AssertionError HA_solver.DiscretizeAR(0.9, -0.2, 10, "T")
    @test_throws AssertionError HA_solver.DiscretizeAR(0.9, 0.2, 0, "T")
    P,ϵ =  HA_solver.DiscretizeAR(0.9, 0.2, 10, "T")
    check_P =1

    for i in 1:size(P,1)
        if abs(sum(P[i,:])- 1.0)>1e-9
            check_P = 0
        end
    end
    @test check_P>0
end


text:
@testset "wrinting_padding" begin
    for i ∈ 1:N_PAD_TESTS
        base = rand(0:256)*Arrow.ALIGNMENT
        r = rand(1:(Arrow.ALIGNMENT-1))
        n = base + r
        @test padding(base) == base
        @test padding(n) == base + Arrow.ALIGNMENT
    end
end

@testset "writing_Primitive" begin
    for i ∈ 1:N_WRITE_TESTS
        dtype = rand(PRIMITIVE_ELTYPES)
        v = reinterpret(UInt8, rand(dtype, rand(1:MAX_VECTOR_LENGTH)))
        p = Primitive(v)
        pad = zeros(UInt8, padding(length(v)) - length(v))
        buff = vcat(convert(Vector{UInt8}, v), pad)
        io = IOBuffer()
        @test writepadded(io, p, values) % 8 == 0
        @test rawpadded(p, values) == buff
        @test take!(io) == buff
    end
end

@testset "writing_NullablePrimitive" begin
    for i ∈ 1:N_WRITE_TESTS
        v, buff = rand_nullableprimitive_buffer2()
        p = NullablePrimitive(v)
        io = IOBuffer()
        @test writepadded(io, p, bitmask, values) % Arrow.ALIGNMENT == 0
        @test rawpadded(p, bitmask, values) == buff
        @test take!(io) == buff
    end
end

@testset "writing_List" begin
    for i ∈ 1:N_WRITE_TESTS
        len = rand(1:MAX_VECTOR_LENGTH)
        offstype = rand([Int32,Int64])
        strs = String[randstring(rand(0:MAX_STRING_LENGTH)) for i ∈ 1:(len-1)]
        # won't work if all strings are empty
        push!(strs, "a")
        offs = Arrow.offsets(offstype, UInt8, strs)
        offs_reint = reinterpret(UInt8, offs)
        vals = convert(Vector{UInt8}, codeunits(reduce(*, strs)))
        valpad = makepad(length(vals))
        buff = vcat(offs_reint, makepad(length(offs_reint)), vals, valpad)
        l = List{String,offstype}(UInt8, strs)
        io = IOBuffer()
        @test writepadded(io, l, offsets, values) % Arrow.ALIGNMENT == 0
        @test rawpadded(l, offsets, values) == buff
        @test take!(io) == buff
    end
end


text:
links = map(x -> split(x, '-'), eachline("input.txt"))
caverns = links |> Base.Iterators.flatten |> collect |> sort! |> unique!
const graph = Dict(c => Set{typeof(c)}() for c in caverns)
for c in caverns
    for link in links
        c in link && push!(graph[c], (only∘setdiff)(link, (c, )))
    end
    graph[c] = setdiff(graph[c], ("start",))
end

const smalls = setdiff(filter(islowercase ∘ first, caverns), ("start", "end"))
function count_small(path) # return true if no small cave is visited twice
    for s in smalls
        count(==(s), path) > 1 && return false
    end
    return true
end

function find_paths(path)
    paths = Vector{typeof(path)}()
    for adj ∈ graph[last(path)]
        if adj ∉ path || isuppercase(first(adj)) || count_small(path)
            newpath = copy(path)
            push!(newpath, adj)
            if adj == "end"
                push!(paths, newpath)
            else
                newpaths = find_paths(newpath)
                for np in newpaths; push!(paths, np); end
            end
        end
    end
    return paths
end

paths = find_paths(["start"])
@info "Part 1: $(length(filter(count_small, paths)))"
@info "Part 2: $(length(paths))"


text:
# # Figures neither in the paper nor in the supplement
#

using Rhone2020  #hide
const R = Rhone2020  #hide

using PyPlot #hide

include("../read_measurements.jl"); #hide
include("../derived_quantities.jl"); #hide
use_progress_meter = false #hide
include("../size_evolution_models.jl"); #hide
include("../heat_transfer.jl"); #hide

# ## Moody plot
#
# The Moody plot shows the relation between the Reynolds number and the friction factor.
# Between AM15 on 9 August and AM13 on 21 August the friction factor does not change even though the Reynolds numbers are considerably different.
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams") #hide
rcParams["font.size"] = 9 # general #hide
figure(figsize=(9,5)) #hide
for date in keys(pick) #hide
    errorbar(mean.(mid_309_265[date][:Re][pick[date]]),mean.(mid_309_265[date][:f][pick[date]]),xerr=std.(mid_309_265[date][:Re][pick[date]]),yerr=std.(mid_309_265[date][:f][pick[date]]),fmt="_",label=join(["Aug-",date[1:2]])) #hide
    xlabel("Reynolds number") #hide
    ylabel("Friction coefficient") #hide
end #hide
legend() #hide
xscale("log") #hide
yscale("log") #hide
gcf() #hide

# ### Water offset-temperature as a histogram
#
# This plot shows the same data as Figure 5d in the paper but plotted as a histogram instead of an errorbar plot.
R.plot_tauw_hist(tau_w, tauw_measured) #hide


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))
hist(res, density=true)

savefig(joinpath(@OUTPUT, "ISL-lab-3-res2.svg")) # hide


text:
type Clock
    ptr::Ptr{Void}

    function Clock(ptr::Ptr{Void})
        c = new(ptr)
        finalizer(c, destroy)
        c
    end
end

function Clock()
    Clock(ccall((:sfClock_create, libcsfml_system), Ptr{Void}, ()))
end

function copy(clock::Clock)
    return Clock(ccall((:sfClock_copy, libcsfml_system), Ptr{Void}, (Ptr{Void},), clock.ptr))
end

function destroy(clock::Clock)
    ccall((:sfClock_destroy, libcsfml_system), Void, (Ptr{Void},), clock.ptr)
end

function get_elapsed_time(clock::Clock)
    return ccall((:sfClock_getElapsedTime, libcsfml_system), Time, (Ptr{Void},), clock.ptr)
end

function restart(clock::Clock)
    return ccall((:sfClock_restart, libcsfml_system), Time, (Ptr{Void},), clock.ptr)
end


text:
using Test, TestSetExtensions
using Logging
using Suppressor

output = @capture_out begin
    @testset ExtendedTestSet "top-level tests" begin
        @testset "2nd-level tests 1" begin
            @test true
            @test 1 == 1
        end
        @testset "2nd-level tests 2" begin
            @test true
            @test 1 == 1
        end
    end
end

try
    @info "You should see 3 failing tests with pretty diffs..."
    include(joinpath("..", "diffdemo.jl"))
catch
end
try
    @info "These 4 failing tests don't have pretty diffs to display"
    @testset ExtendedTestSet "not-pretty" begin
        @testset "No pretty diff for matrices" begin
            @test [1 2; 3 4] == [1 4; 3 4]
        end
        @testset "don't diff non-equality" begin
            @test 1 > 2
        end
        @testset "don't diff non-comparisons" begin
            @test iseven(7)
        end
        @testset "errors don't have diffs either" begin
            throw(ErrorException("This test is supposed to throw an error"))
        end
    end
catch
end

@testset ExtendedTestSet "TestSetExtensions Tests" begin
    @testset "check output dots" begin
        @test split(output, '\n')[1] == "...."
    end

    @testset "Auto-run test files" begin
        global file1_run = false
        global file2_run = false
        global file3_run = false

        @includetests

        @test file1_run
        @test file2_run
        @test file3_run
    end

    @testset "run selected test files" begin
        global file1_run = false
        global file2_run = false
        global file3_run = false

        @includetests ["file1", "file3"]

        @test file1_run
        @test !file2_run
        @test file3_run
    end

    @testset "more than one arg to @includetests is an error" begin
        @test_throws LoadError macroexpand(@__MODULE__, :(@includetests one two))
    end
end

@info "ExtendedTestSet{FallbackTestSet} test sets should exit when the first test fails"
@testset "ExtendedTestSet{FallbackTestSet} Tests" begin
    ets_fallback = ExtendedTestSet{Test.FallbackTestSet}

    # Single-level test set
    err = nothing
    try
        @testset ets_fallback "top-level tests" begin
            @test 1 == 2
            @test 1 == 1
        end
    catch err
    end

    @test err isa TestSetExtensions.ExtendedTestSetException
    @test err.msg == "FallbackTestSetException occurred"

    # Nested test sets
    err = nothing
    try
        @testset ets_fallback "top-level tests" begin
            @testset "Test set with failing test" begin
                @test 1 == 2
                @test 1 == 1
            end
            @testset "Test set with no failing tests" begin
                @test 2 == 2
                @test 3 == 3
            end
        end
    catch err
    end

    @test err isa TestSetExtensions.ExtendedTestSetException
    @test err.msg == "FallbackTestSetException occurred"

    # --- Nested DefaultTestSet tests
    #
    # * Tests Test.record(ExtendedTestSet{FallbackTestSet}, DefaultTestSet) needed for
    #   backward compatibility with Julia<=1.3.

    default_test_set = Test.DefaultTestSet

    # ------ DefaultTest nested under single ExtendedTestSet{FallbackTestSet} test set

    # With failing tests
    err = nothing
    try
        @testset ets_fallback "top-level tests" begin
            @testset default_test_set "Failing test" begin
                @test 1 == 2
                @test 1 == 1
            end
            @testset default_test_set "No failing tests" begin
                @test 2 == 2
                @test 3 == 3
            end
        end
    catch err
    end

    @test err isa Test.FallbackTestSetException

    # With no failing tests
    err = nothing
    try
        @testset ets_fallback "top-level tests" begin
            @testset default_test_set "No failing tests" begin
                @test 2 == 2
                @test 3 == 3
            end
        end
    catch err
    end

    @test err === nothing  # Note: isnothing() is not used for backward compatibility
                           # with Julia 1.0

    # ------ DefaultTest nested under multiple ExtendedTestSet{FallbackTestSet} test sets

    # With failing tests
    try
        @testset ets_fallback "top-level tests" begin
            @testset ets_fallback "2nd-level tests" begin
                @testset default_test_set "Failing test" begin
                    @test 1 == 2
                    @test 1 == 1
                end
                @testset default_test_set "No failing tests" begin
                    @test 2 == 2
                    @test 3 == 3
                end
            end
        end
    catch err
    end

    @test err isa TestSetExtensions.ExtendedTestSetException
    @test err.msg == "FallbackTestSetException occurred"

    # With no failing tests
    err = nothing
    try
        @testset ets_fallback "top-level tests" begin
            @testset ets_fallback "2nd-level tests" begin
                @testset default_test_set "No failing tests" begin
                    @test 2 == 2
                    @test 3 == 3
                end
            end
        end
    catch err
    end

    @test err === nothing  # Note: isnothing() is not used for backward compatibility
                           # with Julia 1.0

end


text:
abstract type AbstractRunningVariable{T,C,VT} <: AbstractVector{T} end

"""
    RunningVariable(Zs; cutoff = 0.0, treated = :≥)

Represents the running variable values for data in a regression
discontinuity setting. The discontinuity is at `cutoff`, and `treated` is one of
`[:>; :>=; :≥; :≧; :<; :<=; :≤; :≦ ]' and determines the treatment based on the
running variable value compared to the cutoff.
"""
struct RunningVariable{T,C,VT} <: AbstractRunningVariable{T,C,VT}
    Zs::VT
    cutoff::C
    treated::Symbol
    Ws::BitArray{1}
    ZsC::VT
    function RunningVariable{T,C,VT}(Zs::VT, cutoff::C, treated) where {T,C,VT}
        treated = Symbol(treated)
        if treated ∉ [:>; :>=; :≥; :≧; :<; :<=; :≤; :≦]
            error("treated should be one of [:>; :>=; :≥; :≧; :<; :<=; :≤; :≦ ]")
        elseif treated ∈ [:>=; :≥; :≧]
            treated = :≥
        elseif treated ∈ [:<=; :≤; :≦]
            treated = :≤
        end
        Ws = broadcast(getfield(Base, treated), Zs, cutoff)
        ZsC = Zs .- cutoff
        new(Zs, cutoff, treated, Ws, ZsC)
    end
end

function RunningVariable(Zs::VT, cutoff::C, treated) where {C,VT}
    RunningVariable{eltype(VT),C,VT}(Zs, cutoff, treated)
end

RunningVariable(Zs; cutoff = 0.0, treated = :≥) = RunningVariable(Zs, cutoff, treated)

Base.size(ZsR::AbstractRunningVariable) = Base.size(ZsR.Zs)
Base.maximum(ZsR::AbstractRunningVariable) = Base.maximum(ZsR.Zs)
Base.minimum(ZsR::AbstractRunningVariable) = Base.minimum(ZsR.Zs)
StatsBase.nobs(ZsR::AbstractRunningVariable) = length(ZsR)

Base.@propagate_inbounds function Base.getindex(ZsR::AbstractRunningVariable, x::Int)
    @boundscheck checkbounds(ZsR.Zs, x)
    @inbounds ret = getindex(ZsR.Zs, x)
    return ret
end

Base.@propagate_inbounds function Base.getindex(
    ZsR::AbstractRunningVariable,
    i::AbstractArray,
)
    @boundscheck checkbounds(ZsR, i)
    @inbounds Zs = ZsR.Zs[i]
    RunningVariable(Zs, ZsR.cutoff, ZsR.treated)
end

struct DiscretizedRunningVariable{T,C,VT} <: AbstractRunningVariable{T,C,VT}
    Zs::VT
    cutoff::C
    treated::Symbol
    Ws::BitArray{1}
    ZsC::VT
    weights::Array{Int,1}
    h::Array{Float64,1}
    binmap::Array{Int,1}

    function DiscretizedRunningVariable{T,C,VT}(
        ZsR::RunningVariable{T,C,VT},
        nbins::Int,
    ) where {T,C,VT}
        hist = fit(Histogram, ZsR; nbins = nbins)
        Zs = midpoints(hist.edges...)
        Ws = broadcast(getfield(Base, ZsR.treated), Zs, ZsR.cutoff)
        weights = hist.weights
        h = Zs[2:length(Zs)] .- Zs[1:(length(Zs)-1)]
        binmap = StatsBase.binindex.(Ref(hist), ZsR.Zs)
        new(Zs, ZsR.cutoff, ZsR.treated, Ws, Zs .- ZsR.cutoff, weights, h, binmap)
    end
end

function DiscretizedRunningVariable(ZsR::RunningVariable{T,C,VT}, nbins::Int) where {T,C,VT}
    DiscretizedRunningVariable{T,C,VT}(ZsR, nbins)
end

# Tables interface

Tables.istable(ZsR::AbstractRunningVariable) = true
Tables.columnaccess(ZsR::AbstractRunningVariable) = true
Tables.columns(ZsR::AbstractRunningVariable) = (Ws = ZsR.Ws, Zs = ZsR.Zs, ZsC = ZsR.ZsC)
function Tables.schema(ZsR::AbstractRunningVariable)
    Tables.Schema((:Ws, :Zs, :ZsC), (eltype(ZsR.Ws), eltype(ZsR.Zs), eltype(ZsR.ZsC)))
end

function fit(
    ::Type{Histogram{T}},
    ZsR::AbstractRunningVariable;
    nbins = StatsBase.sturges(length(ZsR)),
) where {T}
    @unpack cutoff, Zs, treated, ZsC = ZsR
    if treated in [:<; :≥]
        closed = :left
    else
        closed = :right
    end
    nbins = iseven(nbins) ? nbins : nbins + 1

    min_Z, max_Z = extrema(Zs)

    bin_width = (max_Z - min_Z) * 1.01 / nbins

    prop_right = (max_Z - cutoff) * 1.01 / (max_Z - min_Z)
    prop_left = (cutoff - min_Z) * 1.01 / (max_Z - min_Z)
    nbins_right = ceil(Int, nbins * prop_right)
    nbins_left = ceil(Int, nbins * prop_left)

    breaks_left = reverse(range(cutoff; step = -bin_width, length = nbins_left))
    breaks_right = range(cutoff; step = bin_width, length = nbins_right)

    breaks = collect(sort([breaks_left; breaks_right]))

    fit(Histogram{T}, ZsR, breaks; closed = closed)
end

@recipe function f(ZsR::AbstractRunningVariable)

    nbins = get(plotattributes, :bins, StatsBase.sturges(length(ZsR)))

    fitted_hist = fit(Histogram, ZsR; nbins = nbins)

    yguide --> "Frequency"
    xguide --> "Running variable"
    grid --> false
    label --> nothing
    fillcolor --> :lightgray
    thickness_scaling --> 1.7
    linewidth --> 0.3
    ylims --> (0, 1.5 * maximum(fitted_hist.weights))

    @series begin
        fitted_hist
    end

    @series begin
        seriestype := :vline
        linecolor := :purple
        linestyle := :dash
        linewidth := 1.7
        [ZsR.cutoff]
    end
end

"""
    RDData(Ys, ZsR::RunningVariable)

A dataset in the regression discontinuity setting. `Ys` is a vector of outcomes.
"""
struct RDData{V,R<:AbstractRunningVariable}
    Ys::V
    ZsR::R
end

StatsBase.nobs(rdd_data::RDData) = nobs(rdd_data.ZsR)

Base.@propagate_inbounds function Base.getindex(rdd_data::RDData, i::AbstractArray)
    @boundscheck checkbounds(rdd_data.Ys, i)
    @boundscheck checkbounds(rdd_data.ZsR, i)

    @inbounds Ys = rdd_data.Ys[i]
    @inbounds ZsR = rdd_data.ZsR[i]
    RDData(Ys, ZsR)
end

# Tables interface

Tables.istable(::RDData) = true
Tables.columnaccess(::RDData) = true
Tables.columns(rdd_data::RDData) = merge((Ys = rdd_data.Ys,), Tables.columns(rdd_data.ZsR))
function Tables.schema(rdd_data::RDData)
    Tables.Schema(
        (:Ys, :Ws, :Zs, :cutoff, :ZsC),
        (
            eltype(rdd_data.Ys),
            eltype(rdd_data.ZsR.Ws),
            eltype(rdd_data.ZsR.Zs),
            typeof(cutoff),
            eltype(rdd_data.ZsR.ZsC),
        ),
    )
end

function Base.getproperty(obj::RDData, sym::Symbol)
    if sym in (:Ys, :ZsR)
        return getfield(obj, sym)
    else
        return getfield(obj.ZsR, sym)
    end
end

function Base.propertynames(obj::RDData)
    (Base.fieldnames(typeof(obj))..., Base.fieldnames(typeof(obj.ZsR))...)
end

abstract type RDDIndexing end

struct Treated <: RDDIndexing end
struct Untreated <: RDDIndexing end

function Base.getindex(ZsR::R, i::Interval) where {R<:Union{RunningVariable,RDData}}
    idx = in.(ZsR.Zs, i)
    Base.getindex(ZsR, idx)
end

@recipe function f(ZsR::AbstractRunningVariable, Ys::AbstractVector)
    RDData(Ys, ZsR)
end

#-------------------------------------------------------
# include this temporarily, from Plots.jl codebase
# until issue 2360 is resolved
error_tuple(x) = x, x
error_tuple(x::Tuple) = x
_cycle(v::AbstractVector, idx::Int) = v[mod(idx, axes(v, 1))]
_cycle(v, idx::Int) = v

nanappend!(a::AbstractVector, b) = (push!(a, NaN); append!(a, b))

function error_coords(errorbar, errordata, otherdata...)
    ed = Vector{Float64}(undef, 0)
    od = [Vector{Float64}(undef, 0) for odi in otherdata]
    for (i, edi) in enumerate(errordata)
        for (j, odj) in enumerate(otherdata)
            odi = _cycle(odj, i)
            nanappend!(od[j], [odi, odi])
        end
        e1, e2 = error_tuple(_cycle(errorbar, i))
        nanappend!(ed, [edi - e1, edi + e2])
    end
    return (ed, od...)
end
#----------------------------------------------------------
@recipe function f(rdd_data::RDData)

    ZsR = rdd_data.ZsR
    nbins = get(plotattributes, :bins, StatsBase.sturges(length(ZsR)))

    fitted_hist = fit(Histogram, ZsR; nbins = nbins)
    zs = StatsBase.midpoints(fitted_hist.edges[1])
    err_length = (zs[2] - zs[1]) / 2
    binidx = StatsBase.binindex.(Ref(fitted_hist), ZsR)

    tmp_df = rdd_data |> DataFrame
    tmp_df.binidx = binidx
    tmp_df = combine(groupby(tmp_df, :binidx), [:Ys => mean, :Ys => length])

    zs = zs[tmp_df.binidx]
    ys = tmp_df.Ys_mean

    yguide --> "Response"
    xguide --> "Running variable"
    grid --> false
    linecolor --> :grey
    linewidth --> 1.2
    thickness_scaling --> 1.7
    legend --> :outertop
    background_color_legend --> :transparent
    foreground_color_legend --> :transparent

    zs_scatter, ys_scatter = error_coords(err_length, zs, ys)

    #ylims --> extrema(ys)

    @series begin
        label --> "Regressogram"
        seriestype := :path
        zs_scatter, ys_scatter
    end

    @series begin
        label := nothing
        seriestype := :vline
        linecolor := :purple
        linestyle := :dash
        linewidth := 1.7
        [ZsR.cutoff]
    end
end