text:
using GmshTools

function meshgeo(geofile; physical=nothing, dim=2, kwargs...)

    io = open(geofile)
    str = read(io, String)
    close(io)

    for (key,val) in kwargs
        key_str = String(key)
        pat = Regex("$key_str\\s*=\\s*[0-9]*.?[0-9]*;")
        sub = SubstitutionString("$key_str = $val;")
        str = replace(str, pat => sub)
    end

    println(str)
    # return

    temp_geo = tempname()
    open(temp_geo, "w") do io
        print(io, str)
    end

    temp_msh = tempname() * ".msh"
    gmsh.initialize()
    gmsh.option.setNumber("Mesh.MshFileVersion",2)
    gmsh.open(temp_geo)
    gmsh.model.mesh.generate(dim)
    gmsh.write(temp_msh)
    gmsh.finalize()

    # @show temp_msh

    # run(`gmsh $temp_geo -2 -format msh2 -o $temp_msh`)
    if dim == 2
        m = read_gmsh_mesh(temp_msh, physical=physical)
    elseif dim == 3
        m = read_gmsh3d_mesh(temp_msh, physical=physical)
    else
        error("gmsh files of dimension $(dim) cannot be read.")
    end

    rm(temp_msh)
    rm(temp_geo)

    return m
end

function meshsegment(L::T, delta::T, udim=2) where T<:Real

  PT = SVector{udim, T}
  CT = SVector{2,Int}

  num_segments = ceil(Int, L/delta)
  actual_delta = L/num_segments
  x = collect(0:num_segments) * actual_delta

  vertices = zeros(PT, num_segments+1)
  for i in 1 : length(vertices)
    a = zeros(T, udim)
    a[1] = x[i]
    vertices[i] = PT(a)
  end

  faces = Array{CT}(undef,num_segments)
  for i in 1 : length(faces)
    faces[i] = SVector(i, i+1)
  end

  Mesh(vertices, faces)
end

function meshcircle(radius::T, delta::T, udim=2) where T<:Real

  PT = SVector{udim,T}
  CT = SVector{2,Int}

    circumf = 2 * pi *radius
    num_segments = ceil(Int, circumf / delta)
    delta = circumf / num_segments
    dalpha = delta / radius
    alpha = collect(0 : num_segments-1) * dalpha

    vertices = Array{PT}(undef,num_segments)
    for i in 1 : num_segments
        a = zeros(T, udim)
        a[1] = radius * cos(alpha[i])
        a[2] = radius * sin(alpha[i])
        vertices[i] = PT(a)
    end

    faces = Array{CT}(undef,num_segments)
    for i in 1 : length(faces)-1
        faces[i] = SVector{2,Int}(i, i+1)
    end
    faces[end] = SVector{2,Int}(num_segments, 1)

    return Mesh(vertices, faces)
end

"""
    meshrectangle(width, height, delta, udim)

Create a mesh for a rectangle of width (along the x-axis) `width` and height (along
    the y-axis) `height`.

The target edge size is `delta` and the dimension of the
    embedding universe is `udim` (>= 2).

The mesh is oriented such that the normal is pointing down. This is subject to change.
"""
function meshrectangle(width::T, height::T, delta::T, udim=3; structured=true) where T
  if !structured
	  @assert udim==3 "Only 3D Unstructured mesh currently supported"
	  return meshrectangle_unstructured(width, height, delta)
  end

  PT = SVector{udim,T}
  CT = SVector{3,Int}

    @assert 2 <= udim

    nx = round(Int, ceil(width/delta));  nx = max(nx,1); dx = width/nx
    ny = round(Int, ceil(height/delta)); ny = max(ny,1); dy = height/ny

    xs = (0:nx) * dx
    ys = (0:ny) * dy

    vertices = zeros(PT, (nx+1)*(ny+1))
    k = 1
    for x in xs
        for y in ys
            p = zeros(T, udim)
            p[1] = x
            p[2] = y
            vertices[k] = PT(p)
            k += 1
        end
    end

    faces = zeros(CT, 2*nx*ny)
    k = 1
    for i in 1 : nx
        for j in 1 : ny
            v11 = (i-1)*(ny+1) + j
            v12 = i*(ny+1) + j
            v21 = (i-1)*(ny+1) + (j+1)
            v22 = i*(ny+1) + (j+1)
            faces[k]   = CT(v11,v21,v12)
            faces[k+1] = CT(v21,v22,v12)
            k += 2
        end
    end

    Mesh(vertices, faces)
end

"""
    meshsphere(radius, delta)
    meshsphere(;radius, h)

Create a mesh of a sphere of radius `radius` by parsing a .geo script
    incorporating these parameters into the GMSH mesher.

The target edge size is `delta`.
"""
function meshsphere(radius, delta)
    s = """
lc = $delta;

Point(1)={0,0,0,lc};
Point(2)={$radius,0,0,lc};
Point(3)={0,$radius,0,lc};
Point(4)={-$radius,0,0,lc};
Point(5)={0,-$radius,0,lc};
Point(6)={0,0,$radius,lc};
Point(7)={0,0,-$radius,lc};

Circle(1)={2,1,3};
Circle(2)={3,1,4};
Circle(3)={4,1,5};
Circle(4)={5,1,2};
Circle(5)={6,1,3};
Circle(6)={3,1,7};
Circle(7)={7,1,5};
Circle(8)={5,1,6};

Line Loop(1)={-1,-6,-7,-4};
Line Loop(2)={1,-5,-8,4};
Line Loop(3)={-2,-3,7,6};
Line Loop(4)={2,3,8,5};

Ruled Surface(1)={1} In Sphere{1};
Ruled Surface(2)={2} In Sphere{1};
Ruled Surface(3)={3} In Sphere{1};
Ruled Surface(4)={4} In Sphere{1};
"""

    fn = tempname()
    io = open(fn, "w")
    try
        print(io, s)
    finally
        close(io)
    end
    fno = tempname() * ".msh"

    gmsh.initialize()
    gmsh.option.setNumber("Mesh.MshFileVersion",2)
    gmsh.open(fn)
    gmsh.model.mesh.generate(2)
    gmsh.write(fno)
    gmsh.finalize()

    m = read_gmsh_mesh(fno)

    rm(fno)
    rm(fn)

    return m

end

meshsphere(;radius, h) = meshsphere(radius, h)

"""
not working yet
"""
function tetmeshsphere(radius,delta)
    s = """
    lc = $delta;

    Point(1)={0,0,0,lc};
    Point(2)={$radius,0,0,lc};
    Point(3)={-$radius,0,0,lc};
    Point(4)={0,$radius,0,lc};
    Point(5)={0,-$radius,0,lc};
    Point(6)={0,0,$radius,lc};
    Point(7)={0,0,-$radius,lc};

    Circle(1)={7,1,3};
    Circle(2)={3,1,6};
    Circle(3)={6,1,2};
    Circle(4)={2,1,7};
    Circle(5)={7,1,4};
    Circle(6)={4,1,6};
    Circle(7)={6,1,5};
    Circle(8)={5,1,7};
    Circle(9)={3,1,5};
    Circle(10)={5,1,2};
    Circle(11)={2,1,4};
    Circle(12)={4,1,3};

    Curve Loop(1)={8,-4,-10};
    Curve Loop(2)={11,-5,-4};
    Curve Loop(3)={5,12,-1};
    Curve Loop(4)={1,9,8};
    Curve Loop(5)={6,-2,-12};
    Curve Loop(6)={9,-7,-2};
    Curve Loop(7)={7,10,-3};
    Curve Loop(8)={6,3,11};

    Surface(1)={1};
    Surface(2)={2};
    Surface(3)={3};
    Surface(4)={4};
    Surface(5)={5};
    Surface(6)={6};
    Surface(7)={7};
    Surface(8)={8};

    Surface Loop(1) = {4, 3, 2, 8, 5, 6, 7, 1};
    Volume(1) = {1};
    Physical Volume(1) = {1};
    """

    fn = tempname()
    io = open(fn, "w")
    try
        print(io, s)
    finally
        close(io)
    end

    fno = tempname() * ".msh"
    gmsh.initialize()
    gmsh.option.setNumber("Mesh.MshFileVersion",2)
    gmsh.open(fn)
    gmsh.model.mesh.generate(3)
    gmsh.write(fno)
    gmsh.finalize()

    m = read_gmsh3d_mesh(fno)

    rm(fno)
    rm(fn)

    return m
end

function meshball(;radius, h)

    fn = joinpath(@__DIR__,"geos/structured_sphere.geo")
    io = open(fn)
    str = read(io, String)
    close(io)

    str = replace(str, "lc = 1.0;" => "lc = $h;")
    str = replace(str, "r = 1.0;" => "r = $radius;")

    temp_geo = tempname()
    open(temp_geo, "w") do io
        print(io, str)
    end

    temp_msh = tempname() * ".msh"
    gmsh.initialize()
    gmsh.option.setNumber("Mesh.MshFileVersion",2)
    gmsh.open(temp_geo)
    gmsh.model.mesh.generate(3)
    gmsh.write(temp_msh)
    gmsh.finalize()

    m = read_gmsh3d_mesh(temp_msh)

    rm(temp_msh)
    rm(temp_geo)

    return m
end

function meshcylinder(;radius, height, h)

    fn = joinpath(@__DIR__,"geos/cylinder3.geo")
    io = open(fn)
    str = read(io, String)
    close(io)

    str = replace(str, "r = 1.0;" => "r = $radius;")
    str = replace(str, "z = 1.0;" => "z = $height;")
    str = replace(str, "h = 1.0;" => "h = $h;")

    temp_geo = tempname()
    open(temp_geo, "w") do io
        print(io, str)
    end

    temp_msh = tempname() * ".msh"
    gmsh.initialize()
    gmsh.option.setNumber("Mesh.MshFileVersion",2)
    gmsh.open(temp_geo)
    gmsh.model.mesh.generate(3)
    gmsh.write(temp_msh)
    gmsh.finalize()
    m = read_gmsh3d_mesh(temp_msh)

    rm(temp_msh)
    rm(temp_geo)

    return m
end

"""
    meshcuboid(width, height, length, delta)

Creates a mesh for a cuboid of width (along the x-axis) `width`, height (along
    the y-axis) `height` and length (along the z-axis) `length` by parsing a .geo script
    incorporating these parameters into the GMSH mesher.

The target edge size is `delta`.
physical => in {"TopPlate", "BottomPlate", "SidePlates", "OpenBox"} extracts and
returns only the specified part of the cuboid

"""
function meshcuboid(width, height, length, delta;physical="ClosedBox")
    s =
"""
lc = $delta;

Point(1)={0,0,0,lc};
Point(2)={0,0,$length,lc};
Point(3)={$width,0,$length,lc};
Point(4)={$width,0,0,lc};
Point(5)={0,$height,0,lc};
Point(6)={0,$height,$length,lc};
Point(7)={$width,$height,$length,lc};
Point(8)={$width,$height,0,lc};

Line(1)={1,2};
Line(2)={2,3};
Line(3)={3,4};
Line(4)={4,1};
Line(5)={5,6};
Line(6)={6,7};
Line(7)={7,8};
Line(8)={8,5};
Line(9)={1,5};
Line(10)={2,6};
Line(11)={3,7};
Line(12)={4,8};

Line Loop(1)={-1,-2,-3,-4};
Line Loop(2)={1,-9,-5,10};
Line Loop(3)={2,-10,-6,11};
Line Loop(4)={3,-11,-7,12};
Line Loop(5)={4,-12,-8,9};
Line Loop(6)={5,6,7,8};

Plane Surface(1)={1};
Plane Surface(2)={2};
Plane Surface(3)={3};
Plane Surface(4)={4};
Plane Surface(5)={5};
Plane Surface(6)={6};

//classify parts of geometry
Physical Surface("TopPlate") = {3};
Physical Surface("BottomPlate") = {5};
Physical Surface("SidePlates") = {1,2,4,6};
Physical Surface("OpenBox") = {1,2,4,5,6};
Physical Surface("ClosedBox") = {1,2,3,4,5,6};

Surface Loop(1)={1,2,3,4,5,6};
Volume(1)={1};

"""

    fn = tempname()
    io = open(fn, "w")
    try
        print(io, s)
    finally
        close(io)
    end

    # feed the file to gmsh
    fno = tempname() * ".msh"

    gmsh.initialize()
    gmsh.option.setNumber("Mesh.MshFileVersion",2)
    gmsh.open(fn)
    gmsh.model.mesh.generate(2)
    gmsh.write(fno)
    gmsh.finalize()

    # m = read_gmsh_mesh(fno)
    m = read_gmsh_mesh(fno,physical=physical)

    rm(fno)
    rm(fn)

    return m

end

"""
    meshrectangle_unstructured(width, height, delta)
	Meshes unstructured rectangle (Delaunay Triangulation)
"""
function meshrectangle_unstructured(width, height, delta)
    s =
		"""
		lc = $delta;

		Point(1)={0,0,0,lc};
		Point(4)={$width,0,0,lc};
		Point(5)={0,$height,0,lc};
		Point(8)={$width,$height,0,lc};

		Line(4)={4,1};
		Line(8)={8,5};
		Line(9)={1,5};
		Line(12)={4,8};

		Line Loop(5)={4,-12,-8,9};

		Plane Surface(5)={5};
		"""

    fn = tempname()
    io = open(fn, "w")
    try
        print(io, s)
    finally
        close(io)
    end

    # feed the file to gmsh
    fno = tempname() * ".msh"

    gmsh.initialize()
    gmsh.option.setNumber("Mesh.MshFileVersion",2)
    gmsh.open(fn)
    gmsh.model.mesh.generate(2)
    gmsh.write(fno)
    gmsh.finalize()

    m = read_gmsh_mesh(fno)

    rm(fno)
    rm(fn)
    return m

end


text:
"""
    DynamicsExpansion

Stores the expansion of the dynamics about a point. Stores the jacobians for
both the full state and error state. The full state Jacobians are accessed via
`D.A` and `D.B`. The error state Jacobians are accessed via `D.fx` and `D.fu`.

# Constructor

    DynamicsExpansion{T}(n, e, m)

where `n` is the state dimension, `e` is the error state dimension, and `m` is
the control dimension.
"""
struct DynamicsExpansion{T}
    f::Vector{T}   # (Nx,)
    ∇f::Matrix{T}  # (Nx, Nx+Nu)
    A::SubArray{T,2,Matrix{T},Tuple{ColonSlice, UnitRange{Int}}, true}
    B::SubArray{T,2,Matrix{T},Tuple{ColonSlice, UnitRange{Int}}, true}
    ∇e::Matrix{T}
    fx::SubArray{T,2,Matrix{T},Tuple{ColonSlice, UnitRange{Int}}, true}
    fu::SubArray{T,2,Matrix{T},Tuple{ColonSlice, UnitRange{Int}}, true}
    tmp::Matrix{T}  # (Nx,Ne)
    function DynamicsExpansion{T}(n, e, m; errstate::Bool=n!=e) where T
        f = zeros(T,n)
        ∇f = zeros(T,n,n+m)
        A = view(∇f, :, 1:n)
        B = view(∇f, :, n+1:n+m)
        if errstate
            ∇e = zeros(T,e,e+m)
        else  # alias the error expansion to avoid extra work
            ∇e = ∇f
        end
        fx = view(∇e, :, 1:e)
        fu = view(∇e, :, e+1:e+m)
        tmp = zeros(T,n,e)
        new{T}(f,∇f, A,B, ∇e,fx,fu, tmp)
    end
end

@inline function RD.jacobian!(sig::InPlace, diff::DiffMethod, model::DiscreteDynamics,
    D::DynamicsExpansion, z::KnotPoint{n,m}
) where {n,m}
    RD.jacobian!(sig, diff, model, D.∇f, D.f, z)
end

@inline function RD.jacobian!(sig::StaticReturn, diff::DiffMethod, model::DiscreteDynamics,
    D::DynamicsExpansion, z::KnotPoint{n,m}
) where {n,m}
    _z = SVector{n+m}(RD.getdata(z))
    RD.jacobian!(sig, diff, model, D.∇f, D.f, RD.StaticKnotPoint(z, _z))
end

"""
    errstate_jacobian!(model, G, Z)

Calculate the error state Jacobians for the trajectory `Z`, storing the
    result in `G`.
"""
function errstate_jacobians!(model::Vector{<:DiscreteDynamics}, G, Z)
    # NOTE: assumes all models have the same statevectortype
    errstate_jacobians!(RD.statevectortype(model[1]), model, G, Z)
end

function errstate_jacobians!(::RD.EuclideanState, model::Vector{<:DiscreteDynamics}, G, Z)
    return nothing
end

function errstate_jacobians!(::RD.StateVectorType, models::Vector{<:DiscreteDynamics}, G, Z)
    N = length(Z)
	for k in eachindex(Z)
        model = models[min(k, N-1)]
		G[k] .= 0
		RD.errstate_jacobian!(RD.statevectortype(model), model, G[k], Z[k])
	end
end

# """
#     dynamics_expansion!(sig, diff, model, D, Z)

# Evaluate the expansion of the dynamics about the trajectory `Z`, storing the
# result in `D`.
# """
# function dynamics_expansion!(sig::FunctionSignature, diff::DiffMethod,
#                              model::DiscreteDynamics, D::Vector{<:DynamicsExpansion}, Z::Traj)
#     for k in eachindex(D)
#         RobotDynamics.jacobian!(sig, diff, model, D[k].∇f, D[k].f, Z[k])
#     end
# end

"""
    error_expansion!(model, D, G)

Evaluate the dynamics error state expansion. Assumes the error state Jacobians `G`
have already been calculated using [`errstate_jacobians!`](@ref) and that the
full state Jacobians have been calculated and stored in `D` using
[`dynamics_expansion`](@ref).
"""
error_expansion!(model::Vector{<:DiscreteDynamics}, D::Vector{<:DynamicsExpansion}, G) =
    _error_expansion!(RD.statevectortype(model[1]), D, model, G)

_error_expansion!(::RD.EuclideanState, D::Vector{<:DynamicsExpansion},
                 model::Vector{<:DiscreteDynamics}, G) = nothing

function _error_expansion!(::RD.RotationState, D::Vector{<:DynamicsExpansion},
    model::Vector{<:DiscreteDynamics}, G
)
    for k in eachindex(D)
        _error_expansion!(D[k], G[k], G[k+1])
    end
end

function _error_expansion!(D, G1, G2)
    matmul!(D.tmp, D.A, G1)  # D.tmp = A * G1
    matmul!(D.fx, G2', D.tmp)   # D.fx = G2'A*G1
    matmul!(D.fu, G2', D.B)
end

ext: jl


text:
# Tested in Julia 1.5.1
using DelimitedFiles, Combinatorics
data = readdlm("./data/adventInput01.txt", Int)
function ReportRepairPerm(data; nums = 2, sumto = 2020)
    for p in permutations(data, nums)
        if sum(p) == sumto
            print(p, " Sum: $(sum(p)),  Product : $(prod(p))")
            break
        end
    end
end
ReportRepairPerm(data)
ReportRepairPerm(data, nums = 3)


text:
#=
We want to simulate a number of warehouses that store the same product. At a regular
interval, a product is required. But the origin of the product does not matter.

We create our own type `Warehouse` with two field that allow to identify the
warehouse and that allow to track its stock by by means of a `Store`.

The production process works adds a random quantity to a random warehouse
(the first available) and works as follows:
* generate a product
* generate the requests for all resources
* yield the requests the will occur first. *note*: if two events occur at the same time,
they will both happen. We deal with this later.
* cancel all the other requests that have not occured yet. For those that have occured,
we decrement the store with the value it was increased by.

The simulation stops when either all warehouses are full (or cannot handle the
produced quantity).

A similar approach can be followed when dealing with `get` requests instead of
`put` requests.

This application illustrates how you can deal with resource concurrency,
i.e. taking whatever resource(s) come(s) available first without blocking the
other ones or introducing unwanted artifacts.
=#

using Logging
using SimJulia
using ResumableFunctions
import Base.show

struct Warehouse
    name::String
    stock::Container
    function Warehouse(env::Environment, name::String; capacity::Int=10)
        return new(name, Container(env,capacity))
    end
end

show(io::Core.IO, w::Warehouse) = print(io::Core.IO, "$(w.name) - capacity $(w.stock.level)/$(w.stock.capacity)")

@resumable function production(env::Environment, warehouses::Array{Warehouse,1})
    while true
        @info "$(env.time) - warehouse status: $(status(warehouses))"
        # Start making the product (variable production time)
        @debug "$(env.time) - Production starts"
        @yield timeout(env, rand(6:10))
        quant = rand(1:5)
        @info "$(env.time) - Production finished ($(quant) units produced)"
        # get an overview of available resources -> Array{Store, 1}
        resources = map(w::Warehouse -> w.stock, warehouses)
        @debug "$(env.time) - resources: $(resources)"
        # Signal a product needs storage -> Array{SimJulia.Put,1}
        signal = map(r -> put(r, quant), resources)
        @debug "$(env.time) - signal: $(signal)"
        # Yield the different events -> Dict{AbstractEvent,SimJulia.StateValue}
        delivery_requests = @yield Operator(SimJulia.eval_or, signal...) # only one required
        @debug "$(env.time) - delivery requests: $(delivery_requests)"
        # Get the delivery status
        delivery_status = [delivery_requests[s].state for s in signal]
        @debug "$(env.time) - delivery status: $(delivery_status)"
        # Select "winning" warehouse from the status
        wh = findfirst(x -> x == SimJulia.processed, delivery_status)
        @debug "$(env.time) - winning warehouse : $(warehouses[wh])"
        @debug "$(env.time) - BAD warehouse status: $(status(warehouses))"
        # Cancel other requests or reduce value if they happened at the same time
        for i in 1:length(warehouses)
            if i ≠ wh
                @debug "$(env.time) - cancelling warehouse request for $(warehouses[i])"
                if haskey(resources[i].put_queue, signal[i])
                    cancel(resources[i], signal[i])
                else
                    warehouses[i].stock.level -= quant
                end
            end
        end
    end
end

function mysim()
    # setup simulation
    @info "\n$("-"^70)\nwarehouse application\n$("-"^70)\n"
    sim = Simulation()
    # warehouses
    names = ["Brussels", "Amsterdam", "Paris"]
    warehouses = [Warehouse(sim, name, capacity=10) for name in names]
    # start production
    @process production(sim, warehouses)
    # start client
    run(sim)
    @info "$(sim.time) - FINAL warehouse status: $(status(warehouses))"
end

# helper function for printing
function status(v::Array{Warehouse,1})
    string("\n",["\t- $(w)\n" for w in v]..., "\n")
end

# logging settings
Logging.disable_logging(LogLevel(-1001)) # (de-)activate debug messages
logger = Logging.ConsoleLogger(stdout, Logging.Debug)

# run simulation
with_logger(logger) do
    mysim()
end


text:
module GradientTest

import Calculus

using Base.Test
using ForwardDiff

include(joinpath(dirname(@__FILE__), "utils.jl"))

#############################
# rosenbrock hardcoded test #
#############################

x = [0.1, 0.2, 0.3]
v = rosenbrock(x)
g = [-9.4, 15.6, 52.0]

for c in (Chunk{1}(), Chunk{2}(), Chunk{3}())

    # single-threaded #
    #-----------------#
    @test_approx_eq g ForwardDiff.gradient(rosenbrock, x, c)

    out = similar(x)
    ForwardDiff.gradient!(out, rosenbrock, x, c)
    @test_approx_eq out g

    out = GradientResult(x)
    ForwardDiff.gradient!(out, rosenbrock, x, c)
    @test_approx_eq ForwardDiff.value(out) v
    @test_approx_eq ForwardDiff.gradient(out) g

    # multithreaded #
    #---------------#
    if ForwardDiff.IS_MULTITHREADED_JULIA
        @test_approx_eq g ForwardDiff.gradient(rosenbrock, x, c; multithread = true)

        out = similar(x)
        ForwardDiff.gradient!(out, rosenbrock, x, c; multithread = true)
        @test_approx_eq out g

        out = GradientResult(x)
        ForwardDiff.gradient!(out, rosenbrock, x, c; multithread = true)
        @test_approx_eq ForwardDiff.value(out) v
        @test_approx_eq ForwardDiff.gradient(out) g
    end
end

########################
# test vs. Calculus.jl #
########################

for f in VECTOR_TO_NUMBER_FUNCS
    v = f(X)
    g = ForwardDiff.gradient(f, X)
    @test_approx_eq_eps g Calculus.gradient(f, X) FINITEDIFF_ERROR
    for c in CHUNK_SIZES
        for usecache in (true, false)
            println("  ...testing $f with (chunk size = $c) and (usecache = $usecache)")
            chunk = Chunk{c}()

            # single-threaded #
            #-----------------#
            out = ForwardDiff.gradient(f, X, chunk; usecache = usecache)
            @test_approx_eq out g

            out = similar(X)
            ForwardDiff.gradient!(out, f, X, chunk; usecache = usecache)
            @test_approx_eq out g

            out = GradientResult(X)
            ForwardDiff.gradient!(out, f, X, chunk; usecache = usecache)
            @test_approx_eq ForwardDiff.value(out) v
            @test_approx_eq ForwardDiff.gradient(out) g

            # multithreaded #
            #---------------#
            if ForwardDiff.IS_MULTITHREADED_JULIA
                out = ForwardDiff.gradient(f, X, chunk; multithread = true, usecache = usecache)
                @test_approx_eq out g

                out = similar(X)
                ForwardDiff.gradient!(out, f, X, chunk; multithread = true, usecache = usecache)
                @test_approx_eq out g

                out = GradientResult(X)
                ForwardDiff.gradient!(out, f, X, chunk; multithread = true, usecache = usecache)
                @test_approx_eq ForwardDiff.value(out) v
                @test_approx_eq ForwardDiff.gradient(out) g
            end
        end
    end
end

end # module


text:
mutable struct SnaStruct
    name::String
    twojmax::Int32
    ncoeff::Int32
    ncoeffall::Int32
    nperdim::Int32
    idxb_max::Int32
    idxu_max::Int32
    idxz_max::Int32
    idxcg_max::Int32
    ntypes::Int32
    nelements::Int32
    ndoubles::Int32
    ntriples::Int32
    beta_max::Int32
    bnormflag::Int32
    chemflag::Int32
    quadraticflag::Int32
    switchflag::Int32
    bzeroflag::Int32
    wselfallflag::Int32
    rcutfacflag::Int32
    twojmaxflag::Int32

    wself::Float64
    rmin0::Float64
    rfac0::Float64
    rcutfac::Float64
    rcutmax::Float64

    map::Vector{Int32}
    idx_max::Vector{Int32}
    idxz::Vector{Int32}
    idxz_block::Vector{Int32}
    idxb::Vector{Int32}
    idxb_block::Vector{Int32}
    idxu_block::Vector{Int32}
    idxcg_block::Vector{Int32}

    rcutsq::Vector{Float64}
    radelem::Vector{Float64}
    wjelem::Vector{Float64}
    bzero::Vector{Float64}
    coeffelem::Vector{Float64}
    rootpqarray::Vector{Float64}
    cglist::Vector{Float64}

    SnaStruct() = new();
end

function idxbcount(twojmax::Int32)
    count = 0
    for j1 = 0:twojmax
        for j2 = 0:j1
            for j = (j1-j2):2:min(twojmax, j1+j2)
                if (j >= j1)
                    count = count + 1
                end
            end
        end
    end
    return count
end

function idxucount(twojmax::Int32)
    count = 0
    for j = 0:twojmax
        for mb = 0:j
            for ma = 0:j
                count = count + 1
            end
        end
    end
    return count
end

function idxzcount(twojmax::Int32)
    count = 0
    for j1 = 0:twojmax
        for j2 = 0:j1
            for j = (j1-j2):2:min(twojmax, j1+j2)
                for mb = 0:j
                    if (2*mb <= j)
                        for ma = 0:j
                            count = count + 1
                        end
                    end
                end
            end
        end
    end
    return count
end

function idxcgcount(twojmax::Int32)
    count = 0
    for j1 = 0:twojmax
        for j2 = 0:j1
            for j = (j1-j2):2:min(twojmax, j1+j2)
                for mb = 0:j1
                    for ma = 0:j2
                        count = count + 1
                    end
                end
            end
        end
    end
    return count
end

function initializesna(ntypes::Int32, twojmax::Int32, chemflag::Int32)

jdim = twojmax + 1
jdimpq = twojmax + 2
if (chemflag == 0)
    nelements = 1;
else
    nelements = ntypes
end

sna = SnaStruct()

sna.name = "SNAP"
sna.idxb_max = idxbcount(twojmax)
sna.idxu_max = idxucount(twojmax)
sna.idxz_max = idxzcount(twojmax)
sna.idxcg_max = idxcgcount(twojmax)

sna.ntypes = ntypes
sna.twojmax = twojmax
sna.nelements = nelements;
sna.ndoubles = nelements*nelements;
sna.ntriples = nelements*nelements*nelements;
sna.ncoeff = sna.idxb_max*sna.ntriples
sna.ncoeffall = sna.ncoeff + 1

sna.map = zeros(Int32, ntypes+1)
sna.idxcg_block = zeros(Int32, jdim*jdim*jdim)
sna.idxz_block = zeros(Int32, jdim*jdim*jdim)
sna.idxb_block = zeros(Int32, jdim*jdim*jdim)
sna.idxu_block = zeros(Int32, jdim)
sna.idx_max = zeros(Int32, 5)
sna.idxz = zeros(Int32, sna.idxz_max*10)
sna.idxb = zeros(Int32, sna.idxb_max*3)

sna.rcutsq = zeros(Float64, (ntypes)*(ntypes))
sna.radelem = zeros(Float64, ntypes+1)
sna.wjelem = zeros(Float64, ntypes+1)
sna.bzero = zeros(Float64, jdim)
sna.rootpqarray = zeros(Float64, jdimpq*jdimpq)
sna.cglist = zeros(Float64, sna.idxcg_max)

return sna

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# bring in the previous definitions
include("rocketparts.jl")
# Space Vehicle Collection
struct NoPayload <: Payload
end

mass(payload::NoPayload) = 0.0

Rocket(tank::Tank, engine::Engine) = Rocket(NoPayload(), tank, engine)

const nopayload = NoPayload()
mutable struct SpaceVehicle
    active_stage::Payload
end

# Define empty case
SpaceVehicle() = SpaceVehicle(NoPayload())

# Cumbersome first attempt (code in book makes no sense)
# function stage_separate!(ship::SpaceVehicle)
#    stage = ship.active_stage
#    if stage isa Rocket
#       ship.active_stage = stage.payload
#       # these original two lines actually make no sense.
#       #stage.payload = NoPayload()
#       #payload = NoPayload()
#       # After the stage separate what should be returned is the payload.
#        # For the example below we return as a ship so we can iterate
#       SpaceVehicle(stage.payload)
#    else
#       nothing
#    end
# end

# Example
# uncomment to see how first version works
# merlin = Engine(845e3, 282, 470)
# tank = Tank(4e3, 111.5e3)
# r2 = Rocket(SpaceProbe(22e3), tank, merlin)
# r1 = Rocket(r2, tank, EngineCluster(merlin, 9))
# ship = SpaceVehicle(r1)

# This is the result you want to see after modification of Engheim's code
# julia> include("rocketcoll.jl")
# SpaceVehicle(Rocket(Rocket(SpaceProbe(22000.0), Tank(4000.0, 111500.0, 107500.0), SingleEngine(845000.0, 282.0, 470.0)), Tank(4000.0, 111500.0, 107500.0), EngineCluster(SingleEngine(845000.0, 282.0, 470.0), 9)))

# julia> post = stage_separate!(ship)
# SpaceVehicle(Rocket(SpaceProbe(22000.0), Tank(4000.0, 111500.0, 107500.0), SingleEngine(845000.0, 282.0, 470.0)))

# julia> postpost = stage_separate!(post)
# SpaceVehicle(SpaceProbe(22000.0))

# julia> postpostpost = stage_separate!(postpost)

# julia>

# Second attempt to work with push and pop
function stage_separate!(ship::SpaceVehicle)
   stage = ship.active_stage
   if stage isa Rocket
      ship.active_stage = stage.payload
      stage.payload
   else
      nothing
   end
end

import Base: popfirst!, pushfirst!

popfirst!(ship::SpaceVehicle) = stage_separate!(ship)

function pushfirst!(ship::SpaceVehicle, rockets...)
    for r in reverse(rockets)
      # the assumption is that only the last rocket can and must have a payload
      # so for all other rockets, the payload becomes the next stage rocket
      # you need the if condition so that the final payload isn't erased (as it is in the book's
      # code - see discussion in example below)
      if r.payload isa NoPayload
        r.payload = ship.active_stage
      end
      ship.active_stage = r
    end
    ship
end

function SpaceVehicle(rockets::Array{Rocket})
    ship = SpaceVehicle()
    pushfirst!(ship, rockets...)
    ship
end

# Example
# note the results in the book make no sense because you lose the SpaceProbe payload
# Notice result on p266 does not match result on p262.
# The inner rocket should have spaceprobe as load, but Engheim's code erases it!
# to fix this we add the if statement in the code above
merlin = Engine(845e3, 282, 470)
tank = Tank(4e3, 111.5e3)
r1 = Rocket(tank, EngineCluster(merlin, 9))
r2 = Rocket(SpaceProbe(22e3), tank, merlin)
ship = SpaceVehicle([r1,r2])

# import iterate interface and add required methods
import Base: iterate

# define single parameter methods
iterate(payload::Payload) = payload, nothing
iterate(r::Rocket) = r, r.payload
iterate(ship::SpaceVehicle) = iterate(ship.active_stage)
# define tuple method
function iterate(::Union{Rocket, SpaceVehicle}, state)
    if state === nothing || state isa NoPayload
        nothing
    else
        iterate(state)
    end
end

# test iteration
for r in ship
  println(mass(r))
end

# optional methods
import Base: IteratorSize

IteratorSize(::SpaceVehicle) = Base.SizeUnknown()
IteratorSize(::Rocket) = Base.SizeUnknown()
IteratorSize(::Payload) = Base.SizeUnknown()

# test map
map(mass, ship)

# set element type for Collection
# import is critical to get desired result
import Base: eltype

eltype(::Type{SpaceVehicle}) = Payload
eltype(::Type{Rocket}) = Payload
#eltype(::Type{Payload}) = Payload

collect(ship)

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
# Autogenerated wrapper script for primecount_jll for x86_64-apple-darwin
export libprimecount, libprimesieve, primecount

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("primecount")
JLLWrappers.@declare_library_product(libprimecount, "@rpath/libprimecount.6.dylib")
JLLWrappers.@declare_library_product(libprimesieve, "@rpath/libprimesieve.9.dylib")
JLLWrappers.@declare_executable_product(primecount)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libprimecount,
        "lib/libprimecount.6.5.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprimesieve,
        "lib/libprimesieve.9.6.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        primecount,
        "bin/primecount",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# ------------------------------------------------------------------------
# categorical

MMI.categorical(::FI, a...; kw...) = categorical(a...; kw...)

# ------------------------------------------------------------------------
# matrix

MMI.matrix(::FI, ::Val{:table}, X; kw...) = Tables.matrix(X; kw...)

# ------------------------------------------------------------------------
# int

MMI.int(::FI, x) = throw(
    DomainError(x, "Can only convert categorical elements to integers. "))

MMI.int(::FI, x::Missing)       = missing
MMI.int(::FI, x::AbstractArray) = int.(x)

# first line is no good because it promotes type to higher ineger type:
# MMI.int(::FI, x::CategoricalValue) = CategoricalArrays.levelcode(x)
MMI.int(::FI, x::CategoricalValue) = CategoricalArrays.level(x)

# ------------------------------------------------------------------------
# classes

MMI.classes(::FI, p::CategoricalPool) = [p[i] for i in 1:length(p)]
MMI.classes(::FI, x::CategoricalValue) = classes(CategoricalArrays.pool(x))
MMI.classes(::FI, v::CategoricalArray) = classes(CategoricalArrays.pool(v))

# ------------------------------------------------------------------------
# schema

MMI.schema(::FI, ::Val{:table}, X; kw...) = schema(X; kw...)

# ------------------------------------------------------------------------
# decoder

struct CategoricalDecoder{V,R}
    classes::CategoricalVector{V, R, V, CategoricalValue{V,R}, Union{}}
end

MMI.decoder(::FI, x) = CategoricalDecoder(classes(x))

(d::CategoricalDecoder{V,R})(i::Integer) where {V,R} =
    CategoricalValue{V,R}(d.classes[i])
(d::CategoricalDecoder)(a::AbstractArray{<:Integer}) = d.(a)

# ------------------------------------------------------------------------
# table

function MMI.table(::FI, cols::NamedTuple; prototype=NamedTuple())
    Tables.istable(prototype) || error("`prototype` is not a table. ")
    if !Tables.istable(cols)
        tuple_of_vectors = tuple((collect(v) for v in values(cols))...)
        names = keys(cols)
        cols = NamedTuple{names}(tuple_of_vectors)
    end
    return Tables.materializer(prototype)(cols)
end

function MMI.table(::FI, A::AbstractMatrix; names=nothing, prototype=nothing)
    if names === nothing
        _names = [Symbol(:x, j) for j in 1:size(A, 2)]
    else
        _names = collect(names)
    end
    matrix_table = Tables.table(A, header=_names)
    prototype === nothing && return matrix_table
    return Tables.materializer(prototype)(matrix_table)
end

# ------------------------------------------------------------------------
# nrows, selectrows, selectcols

function MMI.nrows(::FI, ::Val{:table}, X)
    Tables.columnaccess(X) || return length(collect(X))
    # if has columnaccess
    cols = Tables.columntable(X)
    !isempty(cols) || return 0
    return length(cols[1])
end

MMI.selectrows(::FI, ::Val{:table}, X, ::Colon) = X
MMI.selectcols(::FI, ::Val{:table}, X, ::Colon) = X

function MMI.selectrows(::FI, ::Val{:table}, X, r)
    r = r isa Integer ? (r:r) : r
    # next uncommented line is a hack; see
    # https://github.com/alan-turing-institute/MLJBase.jl/issues/151
    isdataframe(X) && return X[r, :]
    cols = Tables.columntable(X)
    new_cols = NamedTuple{keys(cols)}(tuple((c[r] for c in values(cols))...))
    return Tables.materializer(X)(new_cols)
end

function MMI.selectcols(::FI, ::Val{:table}, X, c::Union{Symbol,Integer})
    cols = Tables.columntable(X) # named tuple of vectors
    return cols[c]
end

function MMI.selectcols(::FI, ::Val{:table}, X, c::AbstractArray)
    cols = Tables.columntable(X) # named tuple of vectors
    newcols = project(cols, c)
    return Tables.materializer(X)(newcols)
end

# -------------------------------
# utils for select*

# project named tuple onto a tuple with only specified `labels` or indices:
project(t::NamedTuple, labels::AbstractArray{Symbol}) =
    NamedTuple{tuple(labels...)}(t)
project(t::NamedTuple, label::Colon) = t
project(t::NamedTuple, label::Symbol) = project(t, [label,])
project(t::NamedTuple, indices::AbstractArray{<:Integer}) =
    NamedTuple{tuple(keys(t)[indices]...)}(tuple([t[i] for i in indices]...))
project(t::NamedTuple, i::Integer) = project(t, [i,])

# utils for selectrows
typename(X)    = split(string(supertype(typeof(X)).name), '.')[end]
isdataframe(X) = typename(X) == "AbstractDataFrame"

# ----------------------------------------------------------------
# univariate finite

# see src/univariate_finite/


text:
@testset "Geometry" begin

@testset "Testing function definitions" begin
    @test isdefined(GRASS, :make_grid)
    @test isdefined(GRASS, :calc_r2)
    @test isdefined(GRASS, :calc_mu)
    @test isdefined(GRASS, :find_nearest_mu)
    @test isdefined(GRASS, :find_nearest_ax)
end

@testset "Testing back/forth conversions" begin
    t = (0.75, 0.5)
    @test (0.75^2 + 0.5^2) == GRASS.calc_r2(t)
    @test sqrt(1.0 - GRASS.calc_r2(t)) == GRASS.calc_mu(t)
    @test GRASS.make_grid(132) == GRASS.make_grid(N=132)
    @test GRASS.calc_mu(t) == GRASS.calc_mu(t[1], t[2])
    @test GRASS.calc_mu((0.0, 0.0)) == 1.0
    @test GRASS.calc_mu((1.0, 2.5)) == 0.0
end

@testset "Testing off-disk coordinate" begin
    t = (1.25, 1.25)
    @test iszero(GRASS.calc_mu(t))
    @test GRASS.calc_r2(t) > one(eltype(t))
end

@testset "Testing normalization" begin
    N = 132
    u1 = 0.0
    u2 = 0.0
    @test GRASS.calc_norm_term(0.0, 0.0, N, u1, u2) == pi ./ (2.0 * N^2)
end

end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
@doc raw"""
    AlternatingGradientProblem <: Problem

An alternating gradient problem consists of
* a `ProductManifold M` ``=\mathcal M = \mathcal M_1 × ⋯ × M_n``
* a cost function ``F(x)``
* a gradient ``\operatorname{grad}F`` that is either
  * given as one function ``\operatorname{grad}F`` returning a tangent vector `X` on `M` or
  * an array of gradient functions ``\operatorname{grad}F_i``, `ì=1,…,n` s each returning a component of the gradient
  which might be allocating or mutating variants, but not a mix of both.

!!! note

    This Problem requires the `ProductManifold` from `Manifolds.jl`, so `Manifolds.jl` to be loaded.

!!! note

    The input of each of the (component) gradients is still the whole vector `x`,
    just that all other then the `i`th input component are assumed to be fixed and just
    the `i`th components gradient is computed / returned.

# Constructors
    AlternatingGradientProblem(M::ProductManifold, F, gradF::Function;
        evaluation=AllocatingEvaluation()
    )
    AlternatingGradientProblem(M::ProductManifold, F, gradF::AbstractVector{<:Function};
        evaluation=AllocatingEvaluation()
    )

Create a alternating gradient problem with an optional `cost` and the gradient either as one
function (returning an array) or a vector of functions.
"""
struct AlternatingGradientProblem{T,MT<:ProductManifold,TCost,TGradient} <:
       AbstractGradientProblem{T}
    M::MT
    cost::TCost
    gradient!!::TGradient
end
function AlternatingGradientProblem(
    M::TM, F::TCost, gradF!!::G; evaluation::AbstractEvaluationType=AllocatingEvaluation()
) where {TM<:ProductManifold,G,TCost}
    return AlternatingGradientProblem{typeof(evaluation),TM,TCost,G}(M, F, gradF!!)
end
function AlternatingGradientProblem(
    M::TM,
    F::TCost,
    gradF!!::AbstractVector{<:TG};
    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
) where {TM<:ProductManifold,TCost,TG}
    return AlternatingGradientProblem{typeof(evaluation),TM,TCost,typeof(gradF!!)}(
        M, F, gradF!!
    )
end

@doc raw"""
    get_gradient(P::AlternatingGradientProblem, x)
    get_gradient!(P::AlternatingGradientProblem, Y, x)

Evaluate all summands gradients at a point `x` on the `ProductManifold M` (in place of `Y`)
"""
function get_gradient(
    p::AlternatingGradientProblem{AllocatingEvaluation,<:AbstractManifold,TC,<:Function}, x
) where {TC}
    return p.gradient!!(p.M, x)
end
function get_gradient(
    p::AlternatingGradientProblem{
        AllocatingEvaluation,<:AbstractManifold,TC,<:AbstractVector
    },
    x,
) where {TC}
    Y = ProductRepr([gi(p.M, x) for gi in p.gradient!!]...)
    return Y
end
function get_gradient!(
    p::AlternatingGradientProblem{AllocatingEvaluation,<:AbstractManifold,TC,<:Function},
    X,
    x,
) where {TC}
    copyto!(p.M, X, x, get_gradient(p, x))
    return X
end
function get_gradient!(
    p::AlternatingGradientProblem{
        AllocatingEvaluation,<:AbstractManifold,TC,<:AbstractVector
    },
    X,
    x,
) where {TC}
    copyto!(p.M, X, x, get_gradient(p, x))
    return X
end
function get_gradient(
    p::AlternatingGradientProblem{MutatingEvaluation,<:AbstractManifold,TC,<:Function}, x
) where {TC}
    Y = zero_vector(p.M, x)
    return p.gradient!!(p.M, Y, x)
end
function get_gradient(
    p::AlternatingGradientProblem{
        MutatingEvaluation,<:AbstractManifold,TC,<:AbstractVector
    },
    x,
) where {TC}
    Y = zero_vector(p.M, x)
    get_gradient!(p, Y, x)
    return Y
end
function get_gradient!(
    p::AlternatingGradientProblem{MutatingEvaluation,<:AbstractManifold,TC,<:Function}, X, x
) where {TC}
    return p.gradient!!(p.M, X, x)
end
function get_gradient!(
    p::AlternatingGradientProblem{
        MutatingEvaluation,<:AbstractManifold,TC,<:AbstractVector
    },
    X,
    x,
) where {TC}
    for (gi, Xi) in zip(p.gradient!!, submanifold_components(p.M, X))
        gi(p.M, Xi, x)
    end
    return X
end

@doc raw"""
    get_gradient(p::AlternatingGradientProblem, k, x)
    get_gradient!(p::AlternatingGradientProblem, Y, k, x)

Evaluate one of the component gradients ``\operatorname{grad}f_k``, ``k∈\{1,…,n\}``, at `x` (in place of `Y`).
"""
function get_gradient(
    p::AlternatingGradientProblem{AllocatingEvaluation,<:AbstractManifold,TC,<:Function},
    k,
    x,
) where {TC}
    return get_gradient(p, x)[p.M, k]
end
function get_gradient(
    p::AlternatingGradientProblem{
        AllocatingEvaluation,<:AbstractManifold,TC,<:AbstractVector
    },
    k,
    x,
) where {TC}
    return p.gradient!![k](p.M, x)
end
function get_gradient!(
    p::AlternatingGradientProblem{AllocatingEvaluation,<:AbstractManifold,TC,<:Function},
    X,
    k,
    x,
) where {TC}
    copyto!(p.M[k], X, p.gradient!!(p.M, x)[p.M, k])
    return X
end
function get_gradient!(
    p::AlternatingGradientProblem{
        AllocatingEvaluation,<:AbstractManifold,TC,<:AbstractVector
    },
    X,
    k,
    x,
) where {TC}
    copyto!(p.M[k], X, p.gradient!![k](p.M, x))
    return X
end
function get_gradient(
    p::AlternatingGradientProblem{MutatingEvaluation,<:AbstractManifold,TC}, k, x
) where {TC}
    X = zero_vector(p.M[k], x[p.M, k])
    get_gradient!(p, X, k, x)
    return X
end
function get_gradient!(
    p::AlternatingGradientProblem{MutatingEvaluation,<:AbstractManifold,TC,<:Function},
    X,
    k,
    x,
) where {TC}
    # this takes a lot more allocations than other methods, but the gradient can only be evaluated in full
    Xf = zero_vector(p.M, x)
    get_gradient!(p, Xf, x)
    copyto!(p.M[k], X, x[p.M, k], Xf[p.M, k])
    return X
end
function get_gradient!(
    p::AlternatingGradientProblem{
        MutatingEvaluation,<:AbstractManifold,TC,<:AbstractVector
    },
    X,
    k,
    x,
) where {TC}
    return p.gradient!![k](p.M, X, x)
end

"""
    AlternatingGradientDescentOptions <: AbstractGradientDescentOptions

Store the fields for an alternating gradient descent algorithm,
see also [`AlternatingGradientProblem`](@ref) and [`alternating_gradient_descent`](@ref).

# Fields
* `direction` (`AlternatingGradient(zero_vector(M, x))` a [`DirectionUpdateRule`](@ref)
* `evaluation_order` – (`:Linear`) – whether
* `inner_iterations`– (`5`) how many gradient steps to take in a component before alternating to the next
  to use a randomly permuted sequence (`:FixedRandom`), a per
  cycle newly permuted sequence (`:Random`) or the default `:Linear` evaluation order.
* `order` the current permutation
* `retraction_method` – (`default_retraction_method(M)`) a `retraction(M,x,ξ)` to use.
* `stepsize` ([`ConstantStepsize`](@ref)`(M)`) a [`Stepsize`](@ref)
* `stopping_criterion` ([`StopAfterIteration`](@ref)`(1000)`)– a [`StoppingCriterion`](@ref)
* `x` the current iterate
* `k`, ì` internal counters for the outer and inner iterations, respectively.

# Constructors

    AlternatingGradientDescentOptions(M, x; kwargs...)

Generate the options for point `x` and and where the keyword arguments
`inner_iterations`, `order_type`, `order`, `retraction_method`, `stopping_criterion`, and `stepsize``
are keyword arguments
"""
mutable struct AlternatingGradientDescentOptions{
    TX,
    TV,
    D<:DirectionUpdateRule,
    TStop<:StoppingCriterion,
    TStep<:Stepsize,
    RM<:AbstractRetractionMethod,
} <: AbstractGradientOptions
    x::TX
    gradient::TV
    direction::D
    stop::TStop
    stepsize::TStep
    order_type::Symbol
    order::Vector{<:Int}
    retraction_method::RM
    k::Int # current iterate
    i::Int # inner iterate
    inner_iterations::Int
end
@deprecate AlternatingGradientDescentOptions(
    x,
    X,
    direction::DirectionUpdateRule;
    retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
    kwargs...,
) AlternatingGradientDescentOptions(
    DefaultManifold(2), x, X, direction; retraction_method=retraction_method, kwargs...
)
function AlternatingGradientDescentOptions(
    M::AbstractManifold,
    x;
    inner_iterations::Int=5,
    order_type::Symbol=:Linear,
    order::Vector{<:Int}=Int[],
    retraction_method::AbstractRetractionMethod=default_retraction_method(M),
    stopping_criterion::StoppingCriterion=StopAfterIteration(1000),
    stepsize::Stepsize=ConstantStepsize(M),
)
    X = zero_vector(M, x)
    return AlternatingGradientDescentOptions{
        typeof(x),
        typeof(X),
        AlternatingGradient,
        typeof(stopping_criterion),
        typeof(stepsize),
        typeof(retraction_method),
    }(
        x,
        X,
        AlternatingGradient(zero_vector(M, x)),
        stopping_criterion,
        stepsize,
        order_type,
        order,
        retraction_method,
        0,
        0,
        inner_iterations,
    )
end
"""
    AlternatingGradient <: DirectionUpdateRule

The default gradient processor, which just evaluates the (alternating) gradient on one of
the components
"""
struct AlternatingGradient{T} <: AbstractStochasticGradientProcessor
    dir::T
end

function (s::AlternatingGradient)(
    p::AlternatingGradientProblem, o::AlternatingGradientDescentOptions, iter
)
    # at begin of inner iterations reset internal vector to zero
    (o.i == 1) && zero_vector!(p.M, s.dir, o.x)
    # update order(k)th component inplace
    get_gradient!(p, s.dir[p.M, o.order[o.k]], o.order[o.k], o.x)
    return o.stepsize(p, o, iter), s.dir # return urrent full gradient
end

# update Armijo to work on the kth gradient only.
function (a::ArmijoLinesearch)(
    p::AlternatingGradientProblem, o::AlternatingGradientDescentOptions, ::Int
)
    X = zero_vector(p.M, o.x)
    X[p.M, o.order[o.k]] .= get_gradient(p, o.order[o.k], o.x)
    a.last_stepsize = linesearch_backtrack(
        p.M,
        x -> p.cost(p.M, x),
        o.x,
        X,
        a.last_stepsize,
        a.sufficient_decrease,
        a.contraction_factor,
        a.retraction_method,
    )
    return a.last_stepsize
end


text:
using CSV
using DataFrames
using Random

# Load simulate_sample(N0, N1):
include("dgp-normal-normal1.jl");

# Save 100 simulated samples
begin
    N0, N1, Niter = 500, 0, 100
    Random.seed!(1)
    df = map(1:Niter) do iter
        _, _, y0, X0, _, _ = simulate_sample(N0, N1)
        df = DataFrame(X0, :auto)
        df[!, :iter] .= iter
        df[!, :y0] = y0
        df
    end;
    df = reduce(vcat, df)
    CSV.write("data/simulated-data-normal-normal.csv", df);
end


text:
module InteractiveFiction

using Printf, Parameters, Random, Pipe, DataStructures

include("types.jl")

export ObjectGraphNode, ObjectGraph, World, DynamicString, Text, Object, Optional
export @do_while, @register_attribute

function toggle_dbg!()
    global dbg_on
    dbg_on = !dbg_on
end
export toggle_dbg!

function Activity(name::Symbol, main::Function)
    a = Activity(name, create_blank_rulebook(:before), create_blank_rulebook(:carry_out), create_blank_rulebook(:after))
    add_rule!(a.carry_out_rules, :activity, main)
    a
end

function Activity(name::Symbol, before::Rulebook, carry::Rulebook, after::Rulebook)
    a = Activity(name, before, carry, after, Dict())
end

function Activity(name::Symbol, main::Function, before::Rulebook, after::Rulebook)
    a = Activity(name, before, create_blank_rulebook(:carry_out), after, Dict())
    add_rule!(a.carry_out_rules, :activity, main)
    a
end

function create_blank_rulebook(n::Symbol)
    Rulebook(n, nothing, [], [], [], Dict())
end

function create_blank_action(n::Symbol)::Action
    nstr = string(n)
    Action(n, [nstr], 1, create_blank_rulebook(:set_action_variables), create_blank_rulebook(:before_rules),
    create_blank_rulebook(:check_rules), create_blank_rulebook(:carry_out_rules),
    create_blank_rulebook(:report_rules), Dict())
end

const name(w::World, s::Symbol) = name(w.objects[s])
const name(o::Object)::String = o.name
const name(r::Rulebook)::String = string(r.name)
const name(r::Rule)::String = string(r.name)
const name(a::Action)::String = string(a.name)
const player(w::World)::Object = w.player
const type(o::Object)::Symbol = o.type

function get_object(w::World, o::Symbol)
    return w.objects[o]
end

function istype(o::Object, t::Symbol)
    w = o.world
    st = [type(o)]
    while length(st) > 0
        s = pop!(st)
        if s == t
            return true
        end
        map(y -> push!(st, y.type), w.object_hierarchy[s].parents)
    end
    return false
end

Base.isequal(o::Object, s::Symbol) = id(o) == s
Base.isequal(s::Symbol, o::Object) = id(o) == s
#
function istype(w::World, o::Symbol, t::Symbol)
    istype(w, w.objects[o], t)
end

function hasattr(o::Object, s::Symbol)
    haskey(o.attrs, s)
end

export name

const id(o::Object)::Symbol= o.id
export id

function make_blank_world(title::String)::World
    w = World(title = title, first_room=nothing, player=nothing)
    make_type!(w, :object, Vector{Symbol}(), make_object_type!)
    make_type!(w, :direction, :object, make_direction_type!)
    make_type!(w, :thing, :object, make_thing_type!)
    make_type!(w, :room, :object, make_room_type!)
    add_base_rulebooks!(w)
    create_directions!(w)
    add_base_actions!(w)
    add_base_activities!(w)
    player = add_thing!(w, "yourself"; id=:player, description = "As good-looking as ever")
    add_object!(w, :room, "Nowhere Room"; id=NOWHERE_ROOM, description = "The void. If you see this, you messed up.")
    w.first_room = w.objects[NOWHERE_ROOM]
    w.player = w.objects[player]
    location(w.player, w.first_room)
    w
end
export make_blank_world

function Base.show(io::IO, o::Object)
    print(io, @sprintf "%s (ID: %s, Type: %s)" o.name o.id string(o.type))
end

function Base.show(io::IO, ::MIME"text/plain", o::Object)
    print(io, @sprintf "%s (ID: %s, Type: %s)\n\nProperties:\n%s\n" o.name o.id string(o.type) join([string(k, " - ", v) for (k, v) in filter(x -> !(x[2] isa AliasOf) ,o.attrs)], "\n"))
end

function Base.show(io::IO, a::Action)
    print(io, a.name)
end

function indent!(w::World)
    w.message_indent += 1
end

function unindent!(w::World)
    w.message_indent -= 1
end

function indent_dbg!(w::World)
    w.message_indent += 1
end

function unindent_dbg!(w::World)
    w.message_indent -= 1
end
export indent!
export unindent!

buf = []

isnewline = true

function sayln!(w::World, str::Text, args...; kwargs...)
    w.print_counter +=1
    #print(repeat('\t', getindex(w.message_indent)))
    #printstyled(string(getindex(w.print_counter), " "); color = :yellow)
    println(str, args...; kwargs...)
    global buf
    global isnewline
    if isnewline
        push!(buf, str)
    else
        buf[end] = string(buf[end], str)
    end
    isnewline = true
    nothing
end

function say!(w::World, str::Text, args...; kwargs...)
    w.print_counter +=1
    #print(repeat('\t', getindex(w.message_indent)))
    #printstyled(string(getindex(w.print_counter), " "); color = :yellow)
    print(str, args...; kwargs...)
    global buf
    global isnewline
    if isnewline
        push!(buf, str)
    else
        buf[end] = string(buf[end], str)
    end
    isnewline = false
    nothing
end
export say!
export sayln!

dbg_on = false
function say_dbg!(w::World, str::Text, args...; kwargs...)
    if !dbg_on
        return
    end
    w.print_counter +=1
    print(repeat('\t', w.message_indent))
    printstyled(getindex(w.print_counter); color = :yellow)
    printstyled("MSG: "; color = :blue)
    print(str, args...; kwargs...)
end
export say_dbg!

function say_dbgln!(w::World, str::Text, args...; kwargs...)
    if !dbg_on
        return
    end
    w.print_counter +=1
    print(repeat('\t', getindex(w.message_indent)))
    printstyled(getindex(w.print_counter); color = :yellow)
    printstyled("MSG: "; color = :blue)
    println(str, args...; kwargs...)
end
export say_dbgln!

function err_dbg!(w::World, str::Text, args...; kwargs...)
    if !dbg_on
        return
    end
    w.print_counter +=1
    print(repeat('\t', getindex(w.message_indent)))
    printstyled(getindex(w.print_counter); color = :yellow)
    printstyled("ERR: "; color = :red)
    println(str, args...; kwargs...)
end
export err_dbg!

function make_object(w::World, t::Symbol, name::String, id=nothing)::Object
    iters = 0
    obj = Object(t, name, isnothing(id) ? Symbol(randstring(8)) : id, Dict(), w)
    stack = [w.object_hierarchy[t]]
    rev_stack = []
    @do_while length(stack) > 0 begin
        #update the type here
        iters += 1 #just a sanity check for loops
        current_type = popfirst!(stack)
        #current_type.build!(obj)
        push!(stack, current_type.parents...)
        push!(rev_stack, current_type)
        if iters > 1000
            println("LIKELY CYCLIC GRAPH PROBLEM")
        end
    end
    @do_while length(rev_stack) > 0 begin
        iters += 1 #just a sanity check for loops
        build_object!(w, obj, pop!(rev_stack))
        if iters > 1000
            println("LIKELY CYCLIC GRAPH PROBLEM")
        end
    end
    obj
end
export make_object

function build_object!(w::World, obj::Object, t::ObjectGraphNode)
    for x in t.attributes
        ops = w.attribute_store[x][t.type]
        #TODO: remember what this even is checking for
        if length(ops) == 3
            set_attr!(obj, x, ops[1])
        else
            set_attr!(obj, x, ops[2])
        end
        if ops[1] isa Vector
             #iterate each option and add an AliasOf
            if length(ops[1]) > 2
                println("cannot handle 3-option things rn")
            end
            set_attr!(obj, ops[1][1], AliasOf(x, ops[1][2]))
            set_attr!(obj, ops[1][2], AliasOf(x, ops[1][1]))
        end
    end
end
function make_type!(world::World,name::Symbol, parent::Symbol, build::Function)
    make_type!(world, name, [parent], build)
end

function make_type!(world::World, name::Symbol, parents::Vector{Symbol}, build!::Function)
    world.object_hierarchy[name] = ObjectGraphNode(name, [],
        map(x -> world.object_hierarchy[x], parents), Set())
    build!(world)
    for x in parents
        push!(world.object_hierarchy[x].children, world.object_hierarchy[name])
    end
end
export make_type!

function can_be!(t::Symbol, w::World, attr_name::Symbol, options::Symbol; usually=nothing::Optional{String})
    can_be!(t, w, attr_name, [options, Symbol(string("not_", string(options)))]; usually=usually)
end

function can_be!(t::Symbol, w::World, attr_name::Symbol, options::Vector; usually=nothing::Optional{Any},
    never=nothing::Optional{Any})
    if !haskey(w.attribute_store, attr_name)
        w.attribute_store[attr_name] = Dict{Symbol, Any}()
    end
    if !haskey(w.attribute_store, t)
        w.attribute_store[t] = Set()
    end
    is_dynamic = usually isa Function
    w.attribute_store[attr_name][t] = (options,
        isnothing(usually) ? options[1] : usually, never, is_dynamic)
    push!(w.object_hierarchy[t].attributes, attr_name)
    map(y -> push!(w.attribute_store[t], (y, attr_name)), options)
end
export can_be!

function has!(t::Symbol, w::World, attr_name::Symbol; usually=nothing, never=nothing)
    if !haskey(w.attribute_store, attr_name)
        w.attribute_store[attr_name] = Dict{Symbol, Any}()
    end
    if !haskey(w.attribute_store, t)
        w.attribute_store[t] = Set()
    end
    push!(w.object_hierarchy[t].attributes, attr_name)
    is_dynamic = usually isa Function
    w.attribute_store[attr_name][t] = (isnothing(usually) ? "" : usually, never, is_dynamic)
    push!(w.attribute_store[t], attr_name)
end
export has!

const NOTHING_ID = :NA
const NOWHERE_ROOM = :NOWHERE_ROOM

function get_attr(w::World, o::Symbol, a::Symbol)
    get_attr(w.objects[o], a)
end

function get_attr(o::Object, a::Symbol)
    if o.id == NOTHING_ID
        throw(ErrorException("can't do stuff with the nothing object"))
    end
    attr = o.attrs[a]
    #if it's dynamic, then we have a function we need to call instead.
    if attr[2]
        return attr[1](o)
    else
        return attr isa AliasOf ? get_attr(o, attr.attr) == a : attr[1]
    end
end
function make_object_type!(w::World)
    t = :object
    can_be!(t, w, :name_plurality, [:single_named, :plural_named]; usually=:single_named)
    can_be!(t, w, :proper_named, [:proper_named, :improper_named]; usually = :improper_named)
    has!(t, w, :indefinite_article; usually="the")
    has!(t, w, :understand_as)
end

function make_direction_type!(w::World)
    t = :direction
    has!(t, w, :opposite)
end

function make_room_type!(w::World)
    t = :room
    can_be!(t, w, :lighted_attr, [:lighted, :dark]; usually = :lighted)
    can_be!(t, w, :visited_attr, :visited; usually = :not_visited)
    has!(t, w, :located_at; usually=Set())
    has!(t, w, :description)
    has!(t, w, :map_connections; usually=Dict())
end

function make_thing_type!(w::World)
    thing = :thing
    can_be!(thing, w, :lit_attr, [:lit, :unlit]; usually = :unlit)
    can_be!(thing, w, :edible_attr, [:edible, :inedible]; usually = :inedible)
    can_be!(thing, w, :portable_attr, [:fixed_in_place, :portable]; usually = :portable)
    can_be!(thing, w, :wearable_attr, :wearable; usually = :not_wearable)
    has!(thing, w, :description; usually = "")
    has!(thing, w, :location; usually = NOWHERE_ROOM)
    can_be!(thing, w, :pushable_between, :pushable_between_rooms; usually = :pushable_between_rooms)
    can_be!(thing, w, :handled_attr, :handled)
    can_be!(thing, w, :described_attr, :described; usually = :described)
    can_be!(thing, w, :mentioned_attr, :mentioned; usually = :mentioned)
    can_be!(thing, w, :listing_marked, :marked_for_listing; usually = :not_marked_for_listing)
    has!(thing, w, :initial_appearance; usually = "")
end

function set_attr!(w::World, id::Symbol, name::Symbol, val::Any)
    set_attr!(w.objects[id], name, val)
end

function set_attr!(o::Object, name::Symbol, val::Any)
    if haskey(o.attrs, name) && o.attrs[name] isa AliasOf
        set_attr!(o, o.attrs[name].attr, val ? name : o.attrs[name].other)
    else
        o.attrs[name] = (val, val isa Function)
    end
end
export set_attr!

function add_object!(w::World, t::Symbol, n::String, args...; kwargs...)::Symbol
    obj = make_object(w, t, n, haskey(kwargs, :id) ? kwargs[:id] : nothing)
    for (k, v) in kwargs
        set_attr!(obj, k, v)
    end
    for k in args
        set_attr!(obj, k, true)
    end
    w.objects[obj.id] = obj
    if !haskey(w.objects_by_type, t)
        w.objects_by_type[t] = []
    end
    push!(w.objects_by_type[t],obj)
    return obj.id
end
export add_object!

function add_object!(w::World, o::Object)
    w.objects[id(o)] = o
    t = type(o)
    if !haskey(w.objects_by_type, t)
        w.objects_by_type[t] = []
    end
    push!(w.objects_by_type[t], o)
end

function add_thing!(w::World, n::String, args...; kwargs...)::Symbol
    add_object!(w, :thing, n, args...; kwargs...)
end
export add_thing!

function add_room!(w::World, n::String, args...; kwargs...)::Symbol
    s = add_object!(w, :room, n, args...; kwargs...)
    if id(w.first_room) == NOWHERE_ROOM
        w.first_room = w.objects[s]
    end
    s
end
export add_room!

function add_rulebook!(w::World, name::Symbol)::Rulebook
    w.rulebooks[name] = create_blank_rulebook(name)
    w.rulebooks[name]
end

function run(f, w::World)
    global buf
    for x in buf
        println(x)
    end
    buf = []
    w.message_indent = 0
    isnewline = true
    nothing
end

function follow_ruleset!(w::World, ruleset::Vector{Rule}, rb::Rulebook, args...)
    res = nothing
    for rule in ruleset
        say_dbgln!(w, @sprintf "Following the %s" name(rule))
        if haskey(rb.variables, :actor) && id(rb.variables[:actor]) != :player && !rule.any_actor
            println(rb.variables, rule.any_actor)
            continue
        end
        res = rule.rule(w, rb, args...)
        if res != nothing
            break
        end
    end
    res
end

function follow_rulebook!(world::World, rulebook::String, args...)
    follow_rulebook!(world, world.rulebooks[rulebook], args...)
end

function follow_rulebook!(w::World, rulebook::Rulebook, args...)
    say_dbgln!(w, @sprintf "Following the %s rulebook" name(rulebook))
    indent_dbg!(w)
    res = follow_ruleset!(w, rulebook.first_rules, rulebook, args...)
    if res == nothing
        res = follow_ruleset!(w, rulebook.rules, rulebook, args...)
        if res == nothing
            res = follow_ruleset!(w, rulebook.last_rules, rulebook, args...)
        end
    end
    unindent_dbg!(w)
    say_dbgln!(w, @sprintf("Finished following the %s rulebook with resolution: %s",
        name(rulebook), res == nothing ? "no outcome" : res))
    res
end
export follow_rulebook!

function do_activity!(w::World, n::Symbol, args...)
    act = w.activities[n]
    follow_rulebook!(w, act.before_rules, act, args...)
    follow_rulebook!(w, act.carry_out_rules, act, args...) #todo: make this "just the one most specific"
    follow_rulebook!(w, act.after_rules, act, args...)
end

function print_name!(w::World, n::Symbol, args...)
    do_activity!(w::World, :printing_name, n, args...)
end

const default_rule(w::World)::Nothing = nothing

const rulenotimpl(s::Symbol)::Function = (w, rb) -> err_dbg!(w, string("rule not implemented: ", string(s)))

export default_rule

function add_rule!(rulebook::Rulebook, name::Symbol, any_actor=true::Bool)
    add_rule!(rulebook, name, default_rule, any_actor)
end

function add_rule!(world::World, rulebook::Symbol, name::Symbol, rule::Function, any_actor=true::Bool)
    add_rule!(world.rulebooks[rulebook], name, rule, any_actor)
end

function add_rule!(rulebook::Rulebook, name::Symbol, rule::Function, any_actor=true::Bool)
    push!(rulebook.rules, Rule(name, rule, any_actor))
end
export add_rule!

function add_rule_first!(rulebook::Rulebook, name::Symbol, rule::Function, any_actor=true::Bool)
    push!(rulebook.first_rules, Rule(name, rule, any_actor))
end
export add_rule_first!

function add_rule_last!(rulebook::Rulebook, name::Symbol, any_actor=true::Bool)
    add_rule_last!(rulebook, name, default_rule, any_actor)
end

function add_rule_last!(rulebook::Rulebook, name::Symbol, rule::Function, any_actor=true::Bool)
    push!(rulebook.last_rules, Rule(name, rule, any_actor))
end
export add_rule_last!

function add_activity!(w::World, n::Symbol, main_rule::Function)
    w.activities[n] = Activity(n, main_rule)
end

function add_activity!(w::World, n::Symbol, before::Rulebook, carry::Rulebook, after::Rulebook)
    w.activities[n] = Activity(n, before, carry, after)
end

function default_rule(world) end

function when_play_begins!(world::World, rule_name::Symbol, rule = default_rule::Function, any_actor=true::Bool)
    add_rule!(world, :when_play_begins, rule_name, rule)
end
export when_play_begins!

function direction!(w::World, n::String, short::Symbol)::Symbol
    add_object!(w, :direction, n, short; understand_as=[string(short)])
end
export direction!

function room!(w::World, name::String, args...; kwargs...)
    if !(:description in keys(kwargs))
        kwargs = (kwargs..., description = @sprintf "It's the %s." name)
    end
    obj = add_object!(w, name, :room, args...; kwargs...)
    if w.first_room.id == "0xNOWHERE_ROOM"
        w.first_room = obj
    end
    obj
end
export room!

function move!(world::World, obj::Object, new_loc::Symbol)
    move!(world, obj, world.rooms[new_loc])
end

function move!(world::World, obj::Object, new_loc::Object)
    say_dbgln!(world, @sprintf "moved %s from %s to %s" name(obj) location(obj) name(new_loc))
    location(obj, id(new_loc))
end

function add_base_rulebooks!(w::World)

    function set_action_vars(world::World, r::Rulebook)
        action = r.variables[:action]
        action.action_variables[:actor] = r.variables[:actor]
        action.action_variables[:nouns] = r.variables[:nouns]
        follow_rulebook!(world, action.set_action_variables, action)
        for (k, v) in action.action_variables
            action.before_rules.variables[k] = v
            action.carry_out_rules.variables[k] = v
            action.report_rules.variables[k] = v
        end
    end

    function intro_text(world::World, rb::Rulebook)
        tot_len = length("-------") + length(world.title) + length("-------")
        say!(world, repeat("-", tot_len))
        say!(world, "\n")
        say!(world, "-------")
        say!(world, world.title)
        say!(world, "-------\n")
        say!(world, repeat("-", tot_len))
        say!(world, "\n\n\n")
    end

    #=
    def descend_processing(world, actor, action, nouns):
        return world.rulebooks["specific action processing rules"].follow(actor=actor, action=action, nouns=nobineuns)

    def work_out_details(world, actor, action, nouns):
        pass
    =#
    function position_player(world::World, r::Rulebook)::Nothing
        move!(world, world.player, world.first_room)
        nothing
    end

    function check_rules(world::World, r::Rulebook)
        #get the previous rulebook to get the variables there.
        ap = world.rulebooks[:action_processing]
        return follow_rulebook!(world, ap.variables[:action].check_rules)
    end

    function carry_out_rules(world::World, r::Rulebook)
        #get the previous rulebook to get the variables there.
        ap = world.rulebooks[:action_processing]
        return follow_rulebook!(world, ap.variables[:action].carry_out_rules)
    end

    function clean_actions(w::World, r::Rulebook)
        for k in keys(r.variables[:action].action_variables)
            delete!(r.variables[:action].action_variables, k)
        end
        for k in keys(r.variables)
            delete!(r.variables, k)
        end
    end
    begins = add_rulebook!(w, :when_play_begins)
    add_rule_first!(begins, :display_banner, intro_text)
    add_rule_first!(begins, :position_player_in_world, position_player)
    add_rule_first!(begins, :initial_room_description, (x, rb) -> begin try_action!(x, "looking"); nothing end)

    add_rulebook!(w, :before)
    add_rulebook!(w, :instead)

    ap = add_rulebook!(w, :action_processing)
    #announce multiple from list, set pronouns are skipped
    #what????
    add_rule_first!(ap, :set_action_variables, set_action_vars)
    add_rule_first!(ap, :before_stage, (x, rb) -> follow_rulebook!(x, x.rulebooks[:before]))
    add_rule!(ap, :carrying_requirements, rulenotimpl(:carrying_requirements))
    add_rule!(ap, :basic_visibility, rulenotimpl(:basic_visibility))
    add_rule!(ap, :basic_accessibility, rulenotimpl(:basic_accessibility))
    add_rule!(ap, :instead_stage, (x, rb) -> follow_rulebook!(x, x.rulebooks[:instead]))
    add_rule!(ap, :requested_actions_require_persuasion, rulenotimpl(:requested_actions_require_persuasion))
    add_rule!(ap, :carry_out_requested_actions, rulenotimpl(:carry_out_requested_actions))
    add_rule!(ap, :investigate_player_awareness, rulenotimpl(:investigate_player_awareness))
    add_rule!(ap, :check_stage, check_rules)
    add_rule!(ap, :carry_out_stage, carry_out_rules)
    add_rule!(ap, :after_stage, rulenotimpl(:after_stage))
    add_rule!(ap, :investigate_player_awareness_after, rulenotimpl(:investigate_player_awareness_after))
    #add_rule(sap, "report stage rule", lambda world, actor, action, nouns: action.report_rules.follow(action=action))
    add_rule_last!(ap, :clean_actions, clean_actions)
    add_rule_last!(ap, :end_action_processing, (x, rb) -> true)
end

function things_located_at(w::World, o::Object; include_indirect=false)
    things = []
    for x in located_at(o)
        push!(things, x)
    end
    things
end

function add_base_activities!(w::World)
    function printing_dark_room(w::World, rb::Rulebook)
        say!(w, "Darkness")
    end

    function printing_desc_dark_room(w::World, rb::Rulebook, ac)
        sayln!(w, "It is pitch dark, and you can't see a thing.")
    end

    function printing_name(w::World, rb::Rulebook, ac, n::Union{Symbol, Object})
        if n isa Symbol
            n = w.objects[n]
        end
        say!(w, name(n))
    end

    function describing_locale(w::World, rb::Rulebook, a::Activity, o::Object)
        pq = a.activity_variables[:locale_priorities]
        for (item, _) in pq
            do_activity!(w, :print_locale_paragraph_about, item)
        end
    end

    function also_see_locale(w::World, rb::Rulebook, a::Activity, o::Object)
        #rb.variables[:mentionable_count] = 0
        also_see = []
        pq = a.activity_variables[:locale_priorities]
        for (item, priority) in pq
            if priority < 0
                marked_for_listing!(w, item, a.activity_variables[:])
                push!(also_see, item)
            end
        end
        if length(also_see) > 0
            if istype(w, o, :supporter) || istype(w, o, :animal)
                say!(w, "On {the :supporter} you ", :supporter => o)
            elseif istype(w, o, :room) && id(o) == location(player)
                    say!(w, "You ")
            else
                say!(w, "In {the :obj} you ", :obj => o)
            end
            say!("can {if :pgraph_count}also {endif}see ", :pgraph_count => rb[:locale_paragraph_count] > 0)
        end
        common_holder = nothing
        contents_form_of_list = true
        for x in also_see
            if common_holder != held_by(x)
                if isnothing(common_holder)
                    common_holder = held_by(x)
                else
                    contents_form_of_list = false
                end
            end
        end
        if contents_form_of_list && !isnothing(common_holder)
            #list the contents of the common holder, as a sentence, including contents,
            #giving brief inventory information, tersely, not listing concealed items,
            #listing marked items only;
            sayln!(w, "print the list of things")
        else
            sayln!(w, "everything else here.")
        end
        if o == location(w.player)
            say!(w, " here")
        end
        sayln!(w, ".")
    end

    #TODO: need to feed activities the corresponding activity not the rulebook

    function init_locale(w::World, rb::Rulebook, a::Activity, o::Object)
        a.activity_variables[:locale_priorities] = PriorityQueue()
        a.activity_variables[:locale_paragraph_count] = 0
        if haskey(a.activity_variables, :mentionable_timestamp)
            a.activity_variables[:mention_timestamp] += 1
        else
            a.activity_variables[:mention_timestamp] = 1
        end
        for obj in things_located_at(w, o)
            a.activity_variables[:locale_priorities][id(obj)] = -5
        end
    end

    #function find_interesting_stuff(w::World, rb::Rulebook, o::Object)

    before_locale = create_blank_rulebook(:before_printing_locale)
    after_locale = create_blank_rulebook(:after_printing_locale)
    carry_out_locale = create_blank_rulebook(:carry_out_locale)

    add_rule_first!(before_locale, :initialise_locale, init_locale)
    add_rule!(carry_out_locale, :interesting_locale_paragraphs, describing_locale)
    add_rule!(carry_out_locale, :you_can_also_see, also_see_locale)

    add_activity!(w, :printing_name_of_a_dark_room, printing_dark_room)
    add_activity!(w, :printing_description_of_a_dark_room, printing_desc_dark_room)
    add_activity!(w, :printing_name, printing_name)
    add_activity!(w, :describing_locale,  before_locale, carry_out_locale,after_locale)
end
function add_base_actions!(w::World)
    create_looking!(w)
end

function create_directions!(world::World)
    n = direction!(world, "north", :n)
    ne = direction!(world, "northeast", :ne)
    nw = direction!(world, "northwest", :nw)
    s = direction!(world, "south", :s)
    se = direction!(world, "southeast", :se)
    sw = direction!(world, "southwest", :sw)
    e = direction!(world, "east", :e)
    w = direction!(world, "west", :w)
    u = direction!(world, "up", :u)
    d = direction!(world, "down", :d)
    i = direction!(world, "inside", :in)
    o = direction!(world, "outside", :out)
    opposite(world, n, s)
    opposite(world, n, s)
    opposite(world, s, n)
    opposite(world, ne, sw)
    opposite(world, nw, se)
    opposite(world, se, nw)
    opposite(world, sw, ne)
    opposite(world, e, w)
    opposite(world, w, e)
    opposite(world, i, o)
    opposite(world, o, i)
end

function run_world!(world::World)
    global buf
    global isnewline
    world.message_indent = 0
    buf = []
    isnewline = true
    follow_rulebook!(world, world.rulebooks[:when_play_begins])
    s = join(buf, "\n")
    println("----------")
    buf = []
    if dbg_on
        println(s)
    end
    isnewline = true
    s
end
export run_world!

function try_action!(w::World, a::Symbol)
    try_action!(w, w.actions[a], [], player(w))
end

function try_action!(w::World, a::String)
    try_action!(w, w.actions[Symbol(a)], [], player(w))
end

function try_action!(w::World, a::Action)
    try_action!(w, a, [], player(w))
end

function try_action!(world::World, action::Action, nouns::Vector)
    try_action!(world, action, nouns, player(world))
end

function try_action!(world::World, action::Action, nouns::Vector, actor::Object)
    n = (name(actor) == "yourself" ? "the player" : name(actor))
    say_dbgln!(world, @sprintf "%s is trying to do %s" n action.name)
    indent_dbg!(world)
    ap = world.rulebooks[:action_processing]
    ap.variables[:actor] = actor
    ap.variables[:action] = action
    ap.variables[:nouns] = nouns
    outcome = follow_rulebook!(world, ap)
    unindent_dbg!(world)
    if outcome != nothing
        say_dbgln!(world, @sprintf "outcome of %s trying to do %s is %s" n action outcome)
    end
    say_dbgln!(world, @sprintf "%s has finished trying to perform the action %s" n action)
    return outcome
end
export try_action!

const actor(r::Rulebook) = r.variables[:actor]
function create_looking!(w::World)
    look = create_blank_action(:looking)
    look.understand_as = ["look"]
    look.applies_to = 0
    w.actions[:looking] = look

    function find_visibility_holder(w::World, s::Symbol)
        find_visibility_holder(w, w.objects[s])
    end

    function find_visibility_holder(w::World, o::Object)
        if hasattr(o, :supported_by)
            return supported_by(o)
        elseif hasattr(o, :enclosed_by)
            return enclosed_by(o)
        else #TODO: add in opaque objects
            return location(o)
        end
    end
    #abbbrv form allowed is "are we using go"
    function desc_heading_rule(w::World, r::Rulebook)
        #we can safely asssume if we are here, the actor is the player.
        if r.variables[:visibility_level] == 0
            do_activity!(w, :printing_name_of_a_dark_room)
        elseif r.variables[:visibility_ceiling] == location(w.player)
            print_name!(w, location(w.player))
        else
            print_name!(w, r.variables[:visibility_ceiling])
            sayln!()
        end
        imd_lvl = find_visibility_holder(w, actor(r))
        for i = 2:r.variables[:visibility_level]
            if istype(w, imd_lvl, :supporter)
                say!(w, "(on ")
            else
                say!(w, "(in ")
            end
            print_name!(w, imd_lvl)
            say!(w, ")")
            imd_lvl = find_visibility_holder(w, imd_lvl)
        end
        sayln!(w, "\n") #TODO: "run paragraph on with special look spacing"?
    end

    function desc_body_rule(world::World, r::Rulebook)
        if r.variables[:visibility_level] == 0
            if world.room_descriptions == :abbreviated ||
                (world.room_descriptions == :sometimes_abbreviated && world.darkness_witnessed)
                return
            end
            do_activity!(w, :printing_description_of_a_dark_room)
        elseif r.variables[:visibility_ceiling] == location(world.player)
            if world.room_descriptions == :abbreviated ||
                (world.room_descriptions == :sometimes_abbreviated && r.variables[:room_describing_action] != :looking)
                return
            end
            say!(world, @pipe world.player |> location |> description(world, _))
        end
    end

    #I have 0 clue what this code does but it's what the spec says.
    function desc_obj_rule(world::World, r::Rulebook)
        if r.variables[:visibility_level] == 0
            return
        end
        imd_lvl = actor(r)
        ip_cnt = r.variables[:visibility_level]
        while ip_cnt > 0
            marked_for_listing(imd_lvl, true)
            imd_lvl = find_visibility_holder(w, imd_lvl)
            ip_cnt -= 1
        end
        ip_cnt2 = r.variables[:visibility_level]
        while ip_cnt2 > 0
            imd_lvl = actor(r)
            ip_cnt = 0
            while ip_cnt < ip_cnt2
                imd_lvl = find_visibility_holder(w, imd_lvl)
                ip_cnt += 1
                do_activity!(w, :describing_locale, get_object(w, imd_lvl))
                ip_cnt2 -= 1
            end
        end
    end

    function check_arrival_rule(w::World, r::Rulebook)
        if r.variables[:visibility_level] == 0
            w.darkness_witnessed = true
        else
            visited(w, location(w.player))
        end
    end

    function other_people_looking(world::World, r::Rulebook)
        if actor(r) != id(world.player)
            print_name!(w, actor(r))
            sayln!(w, " looks around.")
        end
    end

    function det_vis_ceil(w::World, rb::Rulebook, ac::Action)
        ac.action_variables[:visibility_ceiling] = location(w.player)
        ac.action_variables[:visibility_level] = 1
        ac.action_variables[:room_describing_action] = :looking
    end

    look.action_variables[:room_describing_action] = look
    look.action_variables[:visibility_level] = 0
    look.action_variables[:visibility_ceiling] = nothing
    add_rule!(look.set_action_variables, :determine_visibility_ceiling, det_vis_ceil)
    c = look.carry_out_rules
    add_rule!(c, :room_description_heading, desc_heading_rule, false)
    add_rule!(c, :room_description_body, desc_body_rule, false)
    add_rule!(c, :room_description_paragraphs_about_objects, desc_obj_rule, false)
    add_rule!(c, :check_new_arrivals, check_arrival_rule, false)
    add_rule!(look.report_rules, :other_people_looking, other_people_looking)

    look
end
export create_looking!

@register_attribute single_named
@register_attribute plural_named
@register_attribute proper_named
@register_attribute improper_named
@register_attribute name_plurality
@register_attribute understand_as
@register_attribute indefinite_article

@register_attribute opposite

@register_attribute lighted_attr
@register_attribute visited_attr
@register_attribute visited
@register_attribute unvisited
@register_attribute lighted
@register_attribute dark
@register_attribute located_at

@register_attribute lit_attr
@register_attribute lit
@register_attribute unlit
@register_attribute edible_attr
@register_attribute edible
@register_attribute inedible
@register_attribute fixed_in_place
@register_attribute portable_attr
@register_attribute portable
@register_attribute wearable
@register_attribute wearable_attr
@register_attribute not_wearable
@register_attribute description
@register_attribute location
@register_attribute pushable_between
@register_attribute pushable_between_rooms
@register_attribute not_pushable_between_rooms
@register_attribute handled_attr
@register_attribute handled
@register_attribute not_handled
@register_attribute described_attr
@register_attribute described
@register_attribute not_described
@register_attribute mentioned_attr
@register_attribute mentioned
@register_attribute not_mentioned
@register_attribute listing_marked
@register_attribute marked_for_listing
@register_attribute unmarked_for_listing
@register_attribute initial_appearance
end