text:
@doc doc"""
    ProductManifold{TM<:Tuple, TRanges<:Tuple, TSizes<:Tuple} <: Manifold

Product manifold $M_1 \times M_2 \times \dots \times M_n$ with product geometry.
`TRanges` and `TSizes` statically define the relationship between representation
of the product manifold and representations of point, tangent vectors
and cotangent vectors of respective manifolds.

# Constructor

    ProductManifold(M_1, M_2, ..., M_n)

generates the product manifold $M_1 \times M_2 \times \dots \times M_n$.
Alternatively, the same manifold can be contructed using the `×` operator:
`M_1 × M_2 × M_3`.
"""
struct ProductManifold{TM<:Tuple} <: Manifold
    manifolds::TM
end

ProductManifold(manifolds::Manifold...) = ProductManifold{typeof(manifolds)}(manifolds)

function cross(M1::Manifold, M2::Manifold)
    return ProductManifold(M1, M2)
end

function cross(M1::ProductManifold, M2::Manifold)
    return ProductManifold(M1.manifolds..., M2)
end

function cross(M1::Manifold, M2::ProductManifold)
    return ProductManifold(M1, M2.manifolds...)
end

function cross(M1::ProductManifold, M2::ProductManifold)
    return ProductManifold(M1.manifolds..., M2.manifolds...)
end

"""
    submanifold(M::ProductManifold, i::Integer)

Extract the `i`th factor of the product manifold `M`.
"""
function submanifold(M::ProductManifold, i::Integer)
    return M.manifolds[i]
end

"""
    submanifold(M::ProductManifold, i::Val)

Extract the factor of the product manifold `M` indicated by indices in `i`.
For example, for `i` equal to `Val((1, 3))` the product manifold constructed
from the first and the third factor is returned.
"""
function submanifold(M::ProductManifold, i::Val)
    return ProductManifold(select_from_tuple(M.manifolds, i))
end

"""
    submanifold(M::ProductManifold, i::AbstractVector)

Extract the factor of the product manifold `M` indicated by indices in `i`.
For example, for `i` equal to `[1, 3]` the product manifold constructed
from the first and the third factor is returned.

This function is not type-stable, for better preformance use
`submanifold(M::ProductManifold, i::Val)`.
"""
submanifold(M::ProductManifold, i::AbstractVector) = submanifold(M, Val(tuple(i...)))

function isapprox(M::ProductManifold, x, y; kwargs...)
    return all(t -> isapprox(t...; kwargs...), ziptuples(M.manifolds, x.parts, y.parts))
end

function isapprox(M::ProductManifold, x, v, w; kwargs...)
    return all(t -> isapprox(t...; kwargs...), ziptuples(M.manifolds, x.parts, v.parts, w.parts))
end

function representation_size(M::ProductManifold)
    return (mapreduce(m -> prod(representation_size(m)), +, M.manifolds),)
end

manifold_dimension(M::ProductManifold) = mapreduce(manifold_dimension, +, M.manifolds)

struct ProductMetric <: Metric end

function det_local_metric(M::MetricManifold{ProductManifold, ProductMetric}, x::ProductArray)
    dets = map(det_local_metric, M.manifolds, x.parts)
    return prod(dets)
end

function inner(M::ProductManifold, x, v, w)
    subproducts = map(inner, M.manifolds, x.parts, v.parts, w.parts)
    return sum(subproducts)
end

function norm(M::ProductManifold, x, v)
    norms_squared = map(norm, M.manifolds, x.parts, v.parts).^2
    return sqrt(sum(norms_squared))
end

function distance(M::ProductManifold, x, y)
    return sqrt(sum(map(distance, M.manifolds, x.parts, y.parts).^2))
end

function exp!(M::ProductManifold, y, x, v)
    map(exp!, M.manifolds, y.parts, x.parts, v.parts)
    return y
end

function exp(M::ProductManifold, x::ProductRepr, v::ProductRepr)
    return ProductRepr(map(exp, M.manifolds, x.parts, v.parts)...)
end

function log!(M::ProductManifold, v, x, y)
    map(log!, M.manifolds, v.parts, x.parts, y.parts)
    return v
end

function log(M::ProductManifold, x::ProductRepr, y::ProductRepr)
    return ProductRepr(map(log, M.manifolds, x.parts, y.parts)...)
end

function injectivity_radius(M::ProductManifold, x)
    return min(map(injectivity_radius, M.manifolds, x.parts)...)
end

function injectivity_radius(M::ProductManifold)
    return min(map(injectivity_radius, M.manifolds)...)
end

"""
    ProductRetraction(retractions::AbstractRetractionMethod...)

Product retraction of `retractions`. Works on [`ProductManifold`](@ref).
"""
struct ProductRetraction{TR<:Tuple} <: AbstractRetractionMethod
    retractions::TR
end

ProductRetraction(retractions::AbstractRetractionMethod...) = ProductRetraction{typeof(retractions)}(retractions)

function retract!(M::ProductManifold, y, x, v, method::ProductRetraction)
    map(retract!, M.manifolds, y.parts, x.parts, v.parts, method.retractions)
    return y
end

struct InverseProductRetraction{TR<:Tuple} <: AbstractInverseRetractionMethod
    inverse_retractions::TR
end

"""
    InverseProductRetraction(inverse_retractions::AbstractInverseRetractionMethod...)

Product inverse retraction of `inverse_retractions`.
Works on [`ProductManifold`](@ref).
"""
InverseProductRetraction(inverse_retractions::AbstractInverseRetractionMethod...) = InverseProductRetraction{typeof(inverse_retractions)}(inverse_retractions)

function inverse_retract!(M::ProductManifold, v, x, y, method::InverseProductRetraction)
    map(inverse_retract!, M.manifolds, v.parts, x.parts, y.parts, method.inverse_retractions)
    return v
end

function flat!(M::ProductManifold, v::FVector{CotangentSpaceType}, x, w::FVector{TangentSpaceType})
    vfs = map(u -> FVector(CotangentSpace, u), v.data.parts)
    wfs = map(u -> FVector(TangentSpace, u), w.data.parts)
    map(flat!, M.manifolds, vfs, x.parts, wfs)
    return v
end

function sharp!(M::ProductManifold, v::FVector{TangentSpaceType}, x, w::FVector{CotangentSpaceType})
    vfs = map(u -> FVector(TangentSpace, u), v.data.parts)
    wfs = map(u -> FVector(CotangentSpace, u), w.data.parts)
    map(sharp!, M.manifolds, vfs, x.parts, wfs)
    return v
end

"""
    check_manifold_point(M::ProductManifold, x; kwargs...)

Check whether `x` is a valid point on the [`ProductManifold`](@ref) `M`.

The tolerance for the last test can be set using the ´kwargs...`.
"""
function check_manifold_point(M::ProductManifold, x::ProductRepr; kwargs...)
    for t ∈ ziptuples(M.manifolds, x.parts)
        err = check_manifold_point(t...; kwargs...)
        err === nothing || return err
    end
    return nothing
end

function check_manifold_point(M::ProductManifold, x::ProductArray; kwargs...)
    for t ∈ ziptuples(M.manifolds, x.parts)
        err = check_manifold_point(t...; kwargs...)
        err === nothing || return err
    end
    return nothing
end

"""
    check_tangent_vector(M::ProductManifold, x, v; kwargs... )

Check whether `v` is a tangent vector to `x` on the [`ProductManifold`](@ref)
`M`, i.e. atfer [`check_manifold_point`](@ref)`(M, x)`, and all projections to
base manifolds must be respective tangent vectors.

The tolerance for the last test can be set using the ´kwargs...`.
"""
function check_tangent_vector(M::ProductManifold, x::ProductRepr, v::ProductRepr; kwargs...)
    perr = check_manifold_point(M, x)
    perr === nothing || return perr
    for t ∈ ziptuples(M.manifolds, x.parts, v.parts)
        err = check_tangent_vector(t...; kwargs...)
        err === nothing || return err
    end
    return nothing
end

function check_tangent_vector(M::ProductManifold, x::ProductArray, v::ProductArray; kwargs...)
    perr = check_manifold_point(M, x)
    perr === nothing || return perr
    for t ∈ ziptuples(M.manifolds, x.parts, v.parts)
        err = check_tangent_vector(t...; kwargs...)
        err === nothing || return err
    end
    return nothing
end

"""
    ProductPointDistribution(M::ProductManifold, distributions)

Product distribution on manifold `M`, combined from `distributions`.
"""
struct ProductPointDistribution{TM<:ProductManifold, TD<:(NTuple{N,Distribution} where N)} <: MPointDistribution{TM}
    manifold::TM
    distributions::TD
end

function ProductPointDistribution(M::ProductManifold, distributions::MPointDistribution...)
    return ProductPointDistribution{typeof(M), typeof(distributions)}(M, distributions)
end

function ProductPointDistribution(distributions::MPointDistribution...)
    M = ProductManifold(map(d -> support(d).manifold, distributions)...)
    return ProductPointDistribution(M, distributions...)
end

function support(d::ProductPointDistribution)
    return MPointSupport(d.manifold)
end

function rand(rng::AbstractRNG, d::ProductPointDistribution)
    return ProductRepr(map(d -> rand(rng, d), d.distributions)...)
end

function _rand!(rng::AbstractRNG, d::ProductPointDistribution, x::AbstractArray{<:Number})
    x .= rand(rng, d)
    return x
end

function _rand!(rng::AbstractRNG, d::ProductPointDistribution, x::ProductRepr)
    map(t -> _rand!(rng, t[1], t[2]), d.distributions, x.parts)
    return x
end

"""
    ProductFVectorDistribution([type::VectorBundleFibers], [x], distrs...)

Generates a random vector at point `x` from vector space (a fiber of a tangent
bundle) of type `type` using the product distribution of given distributions.

Vector space type and `x` can be automatically inferred from distributions `distrs`.
"""
struct ProductFVectorDistribution{TSpace<:VectorBundleFibers{<:VectorSpaceType, <:ProductManifold}, TD<:(NTuple{N,Distribution} where N), TX} <: FVectorDistribution{TSpace, TX}
    type::TSpace
    x::TX
    distributions::TD
end

function ProductFVectorDistribution(type::VectorBundleFibers{<:VectorSpaceType, <:ProductManifold}, x::Union{AbstractArray, MPoint, ProductRepr}, distributions::FVectorDistribution...)
    return ProductFVectorDistribution{typeof(type), typeof(distributions), typeof(x)}(type, x, distributions)
end

function ProductFVectorDistribution(type::VectorBundleFibers{<:VectorSpaceType, <:ProductManifold}, distributions::FVectorDistribution...)
    x = ProductRepr(map(d -> support(d).x, distributions))
    return ProductFVectorDistribution(type, x, distributions...)
end

function ProductFVectorDistribution(distributions::FVectorDistribution...)
    M = ProductManifold(map(d -> support(d).space.M, distributions)...)
    VS = support(distributions[1]).space.VS
    all(d -> support(d).space.VS == VS, distributions) || error("Not all distributions have support in vector spaces of the same type, which is currently not supported")
    # Probably worth considering sum spaces in the future?
    x = ProductRepr(map(d -> support(d).x, distributions)...)
    return ProductFVectorDistribution(VectorBundleFibers(VS, M), x, distributions...)
end

function support(tvd::ProductFVectorDistribution)
    return FVectorSupport(tvd.type, ProductRepr(map(d -> support(d).x, tvd.distributions)...))
end

function rand(rng::AbstractRNG, d::ProductFVectorDistribution)
    return ProductRepr(map(d -> rand(rng, d), d.distributions)...)
end

function _rand!(rng::AbstractRNG, d::ProductFVectorDistribution, v::AbstractArray{<:Number})
    v .= rand(rng, d)
    return v
end

function _rand!(rng::AbstractRNG, d::ProductFVectorDistribution, v::ProductRepr)
    map(t -> _rand!(rng, t[1], t[2]), d.distributions, v.parts)
    return v
end


text:
# generative model interface functions

"""
    generate_s{S,A}(p::Union{POMDP{S,A},MDP{S,A}}, s::S, a::A, rng::AbstractRNG)

Return the next state given current state `s` and action taken `a`.
"""
function generate_s end

"""
    generate_o{S,A,O}(p::POMDP{S,A,O}, s::S, a::A, sp::S, rng::AbstractRNG)

Return the next observation given current state `s`, action taken `a` and next state `sp`.

Usually the observation would only depend on the next state `sp`.

    generate_o{S,A,O}(p::POMDP{S,A,O}, s::S, rng::AbstractRNG)

Return the observation from the current state. This should be used to generate initial observations.
"""
function generate_o end

"""
    generate_sr{S}(p::Union{POMDP{S},MDP{S}}, s, a, rng::AbstractRNG)

Return the next state `sp` and reward for taking action `a` in current state `s`.
"""
function generate_sr end

"""
    generate_so{S,A,O}(p::POMDP{S,A,O}, s::S, a::A, rng::AbstractRNG)

Return the next state `sp` and observation `o`.
"""
function generate_so end

"""
    generate_or{S,A,O}(p::POMDP{S,A,O}, s::S, a::A, sp::S, rng::AbstractRNG)

Return the observation `o` and reward for taking action `a` in current state `s` reaching state `sp`.
"""
function generate_or end

"""
    generate_sor{S,A,O}(p::POMDP{S,A,O}, s::S, a::A, rng::AbstractRNG)

Return the next state `sp`, observation `o` and reward for taking action `a` in current state `s`.
"""
function generate_sor end

"""
    initialstate{S}(p::Union{POMDP{S},MDP{S}}, rng::AbstractRNG)

Return the initial state for the problem `p`.

Usually the initial state is sampled from an initial state distribution.
"""
function initialstate end
@deprecate initial_state initialstate

ext: jl


text:
module CollisionDetection

using Compat
using LinearAlgebra

using StaticArrays

export Octree
export boxes, fitsinbox, boudingbox, boxesoverlap
export searchtree, find

include("octree.jl")
include("searcheq.jl")

end # module


text:
if VERSION < v"1.0.0"
    error("VS Code julia language server only works with julia 1.0.0+")
end

using Logging
global_logger(ConsoleLogger(stderr))

@info "Starting the Julia Language Server"

using InteractiveUtils, Sockets

include("../error_handler.jl")

struct LSPrecompileFailure <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::LSPrecompileFailure)
    print(io, ex.msg)
end

try
    if length(Base.ARGS) != 7
        error("Invalid number of arguments passed to julia language server.")
    end

    debug_mode = if Base.ARGS[2] == "--debug=yes"
        true
    elseif Base.ARGS[2] == "--debug=no"
        false
    else
        error("Invalid argument passed.")
    end

    detached_mode = if Base.ARGS[7] == "--detached=yes"
        true
    elseif Base.ARGS[7] == "--detached=no"
        false
    else
        error("Invalid argumentpassed.")
    end

    if debug_mode
        ENV["JULIA_DEBUG"] = "all"
    end

    if detached_mode
        serv = listen(7777)
        global conn_in = accept(serv)
        global conn_out = conn_in
    else
        global conn_in = stdin
        global conn_out = stdout
        (outRead, outWrite) = redirect_stdout()
    end

    try
        using LanguageServer, SymbolServer
    catch err
        if err isa ErrorException && startswith(err.msg, "Failed to precompile")
            throw(LSPrecompileFailure(err.msg))
        else
            rethrow(err)
        end
    end

    symserver_store_path = joinpath(ARGS[5], "symbolstorev5")

    if !ispath(symserver_store_path)
        mkpath(symserver_store_path)
    end

    @info "Symbol server store is at '$symserver_store_path'."

    server = LanguageServerInstance(
        conn_in,
        conn_out,
        Base.ARGS[1],
        Base.ARGS[4],
        (err, bt) -> global_err_handler(err, bt, Base.ARGS[3], "Language Server"),
        symserver_store_path,
        ARGS[6] == "download"
    )
    run(server)
catch err
    global_err_handler(err, catch_backtrace(), Base.ARGS[3], "Language Server")
end

max_stars_repo_name: novitk/julia-vscode

max_issues_repo_name: novitk/julia-vscode

max_forks_repo_name: novitk/julia-vscode


text:
import LinearAlgebra
using LinearAlgebra
using SparseArrays
using DelimitedFiles

using jInv
### matrix free S (slicing) matrix
export normalize
export getData
function getData(m::Array,pFor::PointCloudParam,doClear::Bool=false)
d = 0;
npc = pFor.npcAll;
Parray = pFor.P;
n_points = 0;
#for k=1:length(Parray)
#	n_points+= size(Parray[k],1);
#end
n_points = size(Parray,1);
d = zeros(Float32,n_points);
Jacobian = convert(SparseMatrixCSC{Float64,Int32},spzeros(length(d),length(m)));
Mesh = pFor.Mesh;
if pFor.method == MATFree
	if pFor.method == RBF10Based
		numParamOfRBF = 10;
	else
		numParamOfRBF = 5;
	end
	# global count = 1;
	# for i=1:npc
		# P = Parray[i];
		# dcurr = P*m;
		# d[count:(count + size(P,1) - 1)] = dcurr;
		# Jacobian[count:(count + size(P,1) - 1),:] = P;
		# count = count + size(P,1);
	# end
	# pFor.Jacobian = Jacobian;

	Jacobian = convert(SparseMatrixCSC{Float64,Int32},Matrix(1.0I,length(m),length(m)));
	pFor.Jacobian = Jacobian;
	println("size of u:",size(m)); println("size of jacobian:",size(Jacobian));
	d = m;
	return d,pFor;
end

if pFor.method == RBFBased || pFor.method == RBF10Based || pFor.method == RBF5Based
	Minv = Mesh;
	if pFor.method == RBF10Based
		numParamOfRBF = 10;
	else
		numParamOfRBF = 5;
	end
	nRBF 				= div(length(m)- 5*pFor.npcAll,numParamOfRBF);
	(m1,theta_phi,b) 	= splitRBFparamAndRotationsTranslations(m,nRBF,pFor.npcAll,numParamOfRBF);
	theta_phi 			= theta_phi[pFor.workerSubIdxs,:];
	b 					= b[pFor.workerSubIdxs,:];
	mrot,Jrot 			= rotateAndMoveRBF(m1,Minv,theta_phi,b;computeJacobian = 1,numParamOfRBF=numParamOfRBF);
	d,JacT 				= getPCDataRBF(pFor,mrot,theta_phi,b,numParamOfRBF);
	nRBF 				= div(length(m),numParamOfRBF) ;
	# multiply Jacobian with Jrot, and then take the transpose.

	JacobianT = convert(SparseMatrixCSC{Float64,Int32},spzeros(length(m),prod(size(d))));
	IpIdxs = getIpIdxs(pFor.workerSubIdxs,nRBF,pFor.npcAll,numParamOfRBF);

	JacobianT = Jrot'*JacT;
	#JacobianT[IpIdxs,:] = JacT;
	pFor.Jacobian = JacobianT';
end
return d,pFor;
end

function getPCDataRBF(pFor, m, theta_phi,b,numParamOfRBF = 5)
#Mesh = pFor.Mesh;
#n = pFor.Mesh.n;
Mesh = pFor.Mesh;
npc = pFor.npcAll;
#P is a sparse matrix of size k \times n, where k = #points on point cloud , n = volume size of the mesh
Jacobians = Array{SparseMatrixCSC{Float64,Int32}}(undef, npc);
count = 1;
Parray = pFor.P;
n_points = 0;
for k=1:length(Parray)
	n_points+= size(Parray[k],1);
end
d = zeros(Float32,n_points);
for i=1:npc
	P = Parray[i];
	sigmaH = getDefaultHeaviside();
	u,JBuilder = MeshFreeParamLevelSetModelFunc(Mesh,m[:,i],Xc = P,computeJacobian = 1,sigma = sigmaH,bf = 1,numParamOfRBF = numParamOfRBF);
	J1 = getSparseMatrixTransposed(JBuilder);
	#dcurr = P'*u;
	dcurr = u;
	d[count:(count + size(P,1) - 1)] = dcurr;
	count = count + size(P,1);
	Jacobians[i] = (J1);
end
JacT = blockdiag(Jacobians...);
return d,JacT
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
### TaskMCJob type is used for Monte Carlo jobs based on coroutines (tasks)

type TaskMCJob{M<:MCModel, S<:MCSampler, R<:MCRunner, T<:MCTuner} <: MCJob
  model::Vector{M}
  sampler::Vector{S}
  runner::Vector{R}
  tuner::Vector{T}
  heap::Vector{MCHeap}
  send::Function
  receive::Function
  reset::Function
  jobtype::Symbol
  dim::Int
  task::Vector{Task}

  function TaskMCJob(m::Vector{M}, s::Vector{S}, r::Vector{R}, t::Vector{T})
    job::TaskMCJob = new()

    job.dim = length(m)
    @assert job.dim == length(s) == length(r) == length(t) "Number of models, samplers, runners and tuners not equal."

    job.model, job.sampler, job.runner, job.tuner = m, s, r, t
    job.heap = MCHeap[initialize_heap(m[i], s[i], r[i], t[i]) for i = 1:job.dim]
    job.task = Task[Task(()->initialize_task!(job.heap[i], m[i], s[i], r[i], t[i])) for i = 1:job.dim]
    job.send = produce
    job.receive = (i::Int)->consume(job.task[i])
    job.reset = (i::Int, x::Vector{Float64})->reset(job.task[i], x)
    job.jobtype = :task

    job
  end
end

TaskMCJob{M<:MCModel, S<:MCSampler, R<:MCRunner, T<:MCTuner}(m::M, s::S, r::R, t::T) =
  TaskMCJob{M, S, R, T}([m], [s], [r], [t])


text:
using Agents
using BenchmarkTools
using Random
Random.seed!(10)

println("Forest fire model")
include("forestFire_defs.jl")
a = @benchmark data=run!(forest, dummystep, forest_step!, 100, adata=[(:status, length), (:status, count)], when=1:100) setup=(forest = model_initiation(f=0.05, d=0.8, p=0.01, griddims=(20, 20), seed=2))
display(a)

println("\n\nmove_agent!")
a = @benchmark move_agent!(model.agents[1], (3,4), model) setup=(model = model_initiation(f=0.1, d=0.8, p=0.1, griddims=(100, 50), seed=2))
display(a)

println("\n\ngetindex(model)")
a = @benchmark model[250] setup=(model = model_initiation(f=0.1, d=0.8, p=0.1, griddims=(100, 50), seed=2))
display(a)

println("\n\nkill_agent!")
a = @benchmarkable kill_agent!(model.agents[250], model) setup=(model = model_initiation(f=0.1, d=0.8, p=0.1, griddims=(100, 50), seed=2))
display(run(a))

println("\n\nSchelling model")
include("schelling_defs.jl")
a = @benchmark data=run!(model, agent_step!, 10, adata=[:pos, :mood, :group], when=1:10) setup=(model = instantiate_modelS(numagents=370, griddims=(20,20), min_to_be_happy=5))
display(a)

println("\n\nsample! 1 - no space - n == numagents")
include("sample_defs.jl")
a = @benchmark sample!(m, nagents(m)) setup=(m = sampleModelInitialize(n=100))
display(a)

println("\n\nsample! 2 - no space - n == numagents*2")
a = @benchmark sample!(m, nagents(m)*2) setup=(m = sampleModelInitialize(n=100))
display(a)

println("\n\nsample! 3 - no space - n == numagents/2")
a = @benchmark sample!(m, 50) setup=(m = sampleModelInitialize(n=100))
display(a)

println("\n\nsample! 4 - with space - n == numagents")
a = @benchmark sample!(m, nagents(m)) setup=(m = sampleModelInitializeS(n=100))
display(a)

println("\n\nsample! 5 - with space - n == numagents * 2")
a = @benchmark sample!(m, nagents(m)*2) setup=(m = sampleModelInitializeS(n=100))
display(a)

println("\n\nsample! 6 - with space - n == numagents/2")
a = @benchmark sample!(m, 50) setup=(m = sampleModelInitializeS(n=100))
display(a)


text:
const live_mode = true

import ..debug_header

@static if live_mode

    struct IODebug{T <: IO} <: IO
        io::T
    end

    logwrite(iod::IODebug, f, x) = show_io_debug(stdout, "➡ ", f, x)
    logread(iod::IODebug, f, x) = show_io_debug(stdout, "⬅ ", f, x)
    logunread(iod::IODebug, f, x) = show_io_debug(stdout, "♻ ", f, x)

else

    struct IODebug{T <: IO} <: IO
        io::T
        log::Vector{Tuple{String,String}}
    end

    IODebug(io::T) where T <: IO = IODebug{T}(io, [])

    logwrite(iod::IODebug, f, x) = push!(iod.log, ("➡ ", f, x))
    logread(iod::IODebug, f, x) = push!(iod.log, ("⬅ ", f, x))
    logunread(iod::IODebug, f, x) = push!(iod.log, ("♻ ", f, x))

end

Base.wait_close(iod::IODebug) = Base.wait_close(iod.io)

Base.write(iod::IODebug, a...) =
    (logwrite(iod, :write, join(a));
     write(iod.io, a...))

Base.write(iod::IODebug, a::Array) =
    (logwrite(iod, :write, join(a));
     write(iod.io, a))

Base.write(iod::IODebug, x::String) =
    (logwrite(iod, :write, x);
     write(iod.io, x))

Base.unsafe_write(iod::IODebug, x::Ptr{UInt8}, n::UInt) =
    (logwrite(iod, :unsafe_write, unsafe_string(x,n));
     unsafe_write(iod.io, x, n))

Base.unsafe_read(iod::IODebug, x::Ptr{UInt8}, n::UInt) =
    (r = unsafe_read(iod.io, x, n);
     logread(iod, :unsafe_read, unsafe_string(x,n)); r)

Base.read(iod::IODebug, n::Integer) =
    (r = read(iod.io, n);
     logread(iod, :read, String(copy(r))); r)

Base.readavailable(iod::IODebug) =
    (r = readavailable(iod.io);
     logread(iod, :readavailable, String(copy(r))); r)

IOExtras.unread!(iod::IODebug, bytes) =
    (logunread(iod, :unread!, String(copy(bytes)));
     unread!(iod.io, bytes))

Base.eof(iod::IODebug) = eof(iod.io)
Base.close(iod::IODebug) = close(iod.io)
Base.isopen(iod::IODebug) = isopen(iod.io)
Base.iswritable(iod::IODebug) = iswritable(iod.io)
Base.isreadable(iod::IODebug) = isreadable(iod.io)
IOExtras.startread(iod::IODebug) = startread(iod.io)
IOExtras.startwrite(iod::IODebug) = startwrite(iod.io)
IOExtras.closeread(iod::IODebug) = closeread(iod.io)
IOExtras.closewrite(iod::IODebug) = closewrite(iod.io)

Base.bytesavailable(iod::IODebug) = bytesavailable(iod.io)

Base.show(io::IO, iod::IODebug) = show(io, iod.io)

function show_log(io::IO, iod::IODebug)
    lock(io)
    println(io, "$(typeof(iod)):\nio:     $(iod.io)")
    prevop = ""
    for (operation, f, bytes) in iod.log
        if prevop != "" && prevop != operation
            println(io)
        end
        show_io_debug(io, operation, f, bytes)
        prevop = operation
    end
    println(io)
    unlock(io)
end

function show_io_debug(io::IO, operation, f, bytes)
    prefix = string(debug_header(), rpad(operation, 4))
    i = j = 1
    while i <= length(bytes)
        j = findnext(isequal('\n'), bytes, i)
        if j === nothing || j == 0
            j = prevind(bytes, length(bytes)+1)
        end
        println(io, prefix, "\"", escape_string(bytes[i:j]), "\"",
                i == 1 ? " ($f)" : "")
        if i == 1
            prefix = rpad("DEBUG:", length(prefix) - 1)
        end
        i = nextind(bytes, j)
    end
end


text:
using Dates
using DataFrames
using CategoricalArrays

"""
    strp_datetimes(datetimes::T, format::Union{String,DateFormat} = "y-m-d H:M:S") where T <: AbstractArray{<:AbstractString}

Convert an array of timestamps and to an array of `DateTime` objects.

Any of the strings it's unable to parse, will be replaced with `missing`.

# Examples

```julia-repl
julia> date_strings = [
    "2021-01-27 14:03:25",
    "1999-10-05 01:13:43",
    "abcdefg"
    ];

julia> strp_datetimes(date_strings)
3-element Array{Union{Missing, DateTime},1}:
 2021-01-27T14:03:25
 1999-10-05T01:13:43
 missing
```
"""
function strp_datetimes(datetimes::T, format::Union{String,DateFormat} = "y-m-d H:M:S") where T <: AbstractArray{<:AbstractString}
    if format isa String
        format = DateFormat(format)
    end
    broadcast(
        s::String -> begin
            try
                DateTime(s,format)
            catch e
                if e isa ArgumentError
                    missing
                else
                    rethrow(e)
                end
            end
        end,
        datetimes
    )
end

"""
    extract_datetime_features(datetimes::T) where T <: AbstractArray{<:DateTime}

Extract a `DataFrame` of features from an array of `DateTime` objects.
Features extracted:

* `year`: Year from `datetime`
* `month`: Month from `datetime`
* `dayofmonth`: Day of the month (0-31)
* `dayofweek`: Day of the week (ordered)
* `isweekend`: Is `datetime` a weekend?
* `quarter`: The quarter from datetimes
* `hour`: Hour of the day from `datetime`
* `minute`: Minute from `datetime`
* `second`: Second from `datetime`
* `isAM`: Is time AM (vs PM)?

The same as the following:
```julia-repl
julia>  hcat(
    extract_date_features(datetimes),
    extract_time_features(datetimes)
    )
```

# Examples

```julia-repl
julia> data = strp_datetimes([
    "2021-01-27 14:03:25",
    "1999-10-05 01:13:43",
    "2010-06-11 11:00:00"
]);

julia> extract_datetime_features(data)
3×10 DataFrame
 Row │ year   month    dayofmonth  dayofweek  isweekend  quarter  hour   minut ⋯
     │ Int64  Cat…     Int64       Cat…       Bool       Int64    Int64  Int64 ⋯
─────┼──────────────────────────────────────────────────────────────────────────
   1 │  2021  January          27  Wednesday      false        1     14        ⋯
   2 │  1999  October           5  Tuesday        false        4      1      1
   3 │  2010  June             11  Friday         false        2     11
                                                               3 columns omitted
```

See also: [`extract_date_features`](@ref), [`extract_time_features`](@ref)
"""
function extract_datetime_features(datetimes::T) where T <: AbstractArray{<:DateTime}
    hcat(
        extract_date_features(datetimes),
        extract_time_features(datetimes)
    )
end

"""
    extract_time_features(datetimes::T) where T <: AbstractArray{<:Union{Time,DateTime}}

Extract a `DataFrame` of features from an array of `DateTime` or `Time` objects.
Features extracted:

* `hour`: Hour of the day from `datetime`
* `minute`: Minute from `datetime`
* `second`: Second from `datetime`
* `isAM`: Is time AM (vs PM)?

# Examples

```julia-repl
julia> data = strp_datetimes([
    "2021-01-27 14:03:25",
    "1999-10-05 01:13:43",
    "2010-06-11 11:00:00"
]);

julia> extract_time_features(data)
3×4 DataFrame
 Row │ hour   minute  second   isAM
     │ Int64  Int64   Float64  Bool
─────┼───────────────────────────────
   1 │    14       3     25.0  false
   2 │     1      13     43.0   true
   3 │    11       0      0.0   true
```

See also: [`extract_datetime_features`](@ref), [`extract_date_features`](@ref)
"""
function extract_time_features(datetimes::T) where T <: AbstractArray{<:Union{Time,DateTime}}
    DataFrame(
        hour=hour.(datetimes),
        minute=minute.(datetimes),
        second=(second.(datetimes) .+ millisecond.(datetimes) ./ 1e3),
        isAM=hour.(datetimes) .< 12
    )
end

"""
    extract_date_features(datetimes::T) where T <: AbstractArray{<:Union{Date,DateTime}}

Extract a `DataFrame` of features from an array of `DateTime` or `Date` objects.
Features extracted:

* `year`: Year from `datetime`
* `month`: Month from `datetime`
* `dayofmonth`: Day of the month (1-31)
* `dayofweek`: Day of the week (ordered)
* `isweekend`: Is `datetime` a weekend?
* `quarter`: The quarter from datetimes

# Examples

```julia-repl
julia> data = strp_datetimes([
    "2021-01-27 14:03:25",
    "1999-10-05 01:13:43",
    "2010-06-11 11:00:00"
]);

julia> extract_date_features(data)
3×6 DataFrame
 Row │ year   month    dayofmonth  dayofweek  isweekend  quarter
     │ Int64  Cat…     Int64       Cat…       Bool       Int64
─────┼───────────────────────────────────────────────────────────
   1 │  2021  January          27  Wednesday      false        1
   2 │  1999  October           5  Tuesday        false        4
   3 │  2010  June             11  Friday         false        2
```

See also: [`extract_datetime_features`](@ref), [`extract_time_features`](@ref)
"""
function extract_date_features(datetimes::T) where T <: AbstractArray{<:Union{Date,DateTime}}
    DataFrame(
        year=year.(datetimes),
        month=get_month(datetimes),
        dayofmonth=day.(datetimes),
        dayofweek=get_weekday(datetimes),
        isweekend=broadcast(
            d -> d ∈ ("Saturday","Sunday"),
            get_weekday(datetimes)
        ),
        quarter=quarterofyear.(datetimes)
    )
end

"""
    get_month(datetimes::T) where T <: AbstractArray{<:Union{Date,DateTime}}

Return an ordered `CategoricalArray` of month names extracted from `datetimes`.

# Examples:

```julia-repl
julia> data = strp_datetimes([
    "2021-01-27 14:03:25",
    "1999-10-05 01:13:43",
    "2010-06-11 11:00:00"
]);

julia> get_month(data)
3-element CategoricalArrays.CategoricalArray{String,1,UInt32}:
 "January"
 "October"
 "June"
```

See also: [`extract_datetime_features`](@ref), [`extract_date_features`](@ref), [`get_weekday`](@ref)
"""
function get_month(datetimes::T) where T <: AbstractArray{<:Union{Date,DateTime}}
    categorical(
        monthname.(datetimes),
        levels=monthname.(1:12)
    )
end

"""
    get_weekday(datetimes::T) where T <: AbstractArray{<:Union{Date,DateTime}}

Return an ordered `CategoricalArray` of weekday names extracted from `datetimes`.

# Examples:

```julia-repl
julia> data = strp_datetimes([
    "2021-01-27 14:03:25",
    "1999-10-05 01:13:43",
    "2010-06-11 11:00:00"
]);

julia> get_weekday(data)
3-element CategoricalArrays.CategoricalArray{String,1,UInt32}:
 "Wednesday"
 "Tuesday"
 "Friday"
```

See also: [`extract_datetime_features`](@ref), [`extract_date_features`](@ref), [`get_weekday`](@ref)
"""
function get_weekday(datetimes::T) where T <: AbstractArray{<:Union{Date,DateTime}}
    categorical(
        dayname.(datetimes),
        levels=dayname.(1:7)
    )
end


text:
export SimRoundRobinF
"""
    SimRoundRobinF

 Simulates multiple round-robin competition between a list of players under some performance model.
 Differs from SimRoundRobin because it adds in a "home ground" advantage term, which is randomly
 distributed for each match up.

## Parameters
* `n::Int`: number of RR competitions to run
* `factor_scale::Float64`

## Limitations
* no ties
* just win/loss, so scores are trivial

```
"""
struct SimRoundRobinF <: SimulateRule
    n::Int
    factor_scale::Float64
    function SimRoundRobinF(n, factor_scale)
        @check_args(SimRoundRobin, n >= one(n))
        @check_args(SimRoundRobinF, factor_scale >= zero(factor_scale))
        new(n, factor_scale)
    end
end
default_θ = 400 / log(10) # this results in the standard values used for Elo in Chess and many other cases
SimRoundRobinF(; n::Int=1, factor_scale=default_θ) = SimRoundRobinF(n, factor_scale)

function simulate( r::RatingsList, model::SimRoundRobinF, perf_model::ContinuousUnivariateDistribution )
    m = length( r.players )
    L = Int( m*(m-1)/2 )
    df = DataFrame([String, String, Int, Int, Int, Union{Missing,Int}, Union{Missing,Int}], [PlayerA, PlayerB, Outcome, ScoreA, ScoreB, FactorA, FactorB], model.n*L)
    wins = Dict{String, Int}()

    # simulate all players playing each other once per round
    c = 1
    for i=1:model.n
        for j=1:m-1
            for k=j+1:m
                pA = r.players[j]
                pB = r.players[k]
                df[c, PlayerA] = pA
                df[c, PlayerB] = pB
                y = rand()
                if y <= 0.5
                    df[c, FactorA] = 1
                    df[c, FactorB] = 0
                else
                    df[c, FactorA] = 0
                    df[c, FactorB] = 1
                end
                o = outcome( r.ratings[pA], r.ratings[pB],  model.factor_scale*(df[c, FactorA]-df[c, FactorB]), perf_model )
                df[c, Outcome] = o
                if o == 1
                    df[c, ScoreA] = 1
                    df[c, ScoreB] = 0
                    increment!(wins, pA, 1)
                elseif o == -1
                    df[c, ScoreA] = 0
                    df[c, ScoreB] = 1
                    increment!(wins, pB, 1)
                end
                c += 1
            end
        end
    end

    # winners are players who wins most games
    #   if there is a tie, one is picked (by Julia not me)
    winner = findmax(wins)[2]
    return df, winner
end

# random outcome of a single match
#   should be able to push this into a single instantiation, or corresponding update rule?
function outcome( rA::Float64, rB::Float64, f::Float64, perf_model::ContinuousUnivariateDistribution)
    p = cdf( perf_model, rA - rB + f )
    x = rand()
    if x <= p
        o = 1
    else
        o = -1
    end
    return o
end


text:
# To convert video to png
# run(`ffmpeg -r 1 -i $folder/$file.avi -r 1 $folder/png/$file%03d.png`)

# using FileIO

"""
    Reads a video.
    Picks the red channels since all the channels were equal.
    Returns Array{Gray{T}, 3}
"""
function readvideo(filename)
    f = openvideo(filename)
    A1 = []
    while !eof(f)
        push!(A1, read(f))
    end
    # make A1 typed
    function checkgray(x)
        red(x) == green(x) == blue(x)
        Gray(red(x))
    end
    catfun([x for x in A1], checkgray)
end
function catfun(As::AbstractArray{S}, fun) where {TE, T <: Color{TE}, S <: AbstractArray{T}}
    AA = Array{Gray{TE}, 3}(undef, size(As[1])..., length(As))
    for (i, A) in enumerate(As)
        v = view(AA, :, :, i)
        @. v = fun(A)
    end
    AA
end


text:
export ghcubature, GaussHermiteCubature

import FastGaussQuadrature: gausshermite
import LinearAlgebra: mul!, axpy!

using Distributions

const product  = Iterators.product
const repeated = Iterators.repeated
const sqrtPI1  = sqrt(pi)

struct GaussHermiteCubature{PI, WI} <: AbstractApproximationMethod
    piter :: PI
    witer :: WI
end

GaussHermiteCubature(p::Int) = ghcubature(p)

approximation_name(approx::GaussHermiteCubature)       = "GaussHermite($(approx.p))"
approximation_short_name(approx::GaussHermiteCubature) = "GH$(approx.p)"

function ghcubature(p::Int)
    points, weights = gausshermite(p)
    return GaussHermiteCubature(points, weights)
end

function getweights(gh::GaussHermiteCubature, mean::T, variance::T) where { T <: Real }
    return Base.Generator(gh.witer) do weight
        return weight / sqrtPI1
    end
end

function getweights(gh::GaussHermiteCubature, mean::AbstractVector{T}, covariance::AbstractMatrix{T}) where { T <: Real }
    sqrtpi = (pi ^ (length(mean) / 2))
    return Base.Generator(product(repeated(gh.witer, length(mean))...)) do pweight
        return prod(pweight) / sqrtpi
    end
end

function getpoints(gh::GaussHermiteCubature, mean::T, variance::T) where { T <: Real }
    sqrt2V = sqrt(2 * variance)
    return Base.Generator(gh.piter) do point
        return mean + sqrt2V * point
    end
end

function getpoints(cubature::GaussHermiteCubature, mean::AbstractVector{T}, covariance::AbstractMatrix{T}) where { T <: Real }
    sqrtP = cholsqrt(covariance)
    sqrt2 = sqrt(2)

    tbuffer = similar(mean)
    pbuffer = similar(mean)
    return Base.Generator(product(repeated(cubature.piter, length(mean))...)) do ptuple
        copyto!(pbuffer, ptuple)
        copyto!(tbuffer, mean)
        return mul!(tbuffer, sqrtP, pbuffer, sqrt2, 1.0) # point = m + sqrt2 * sqrtP * p
    end
end


text:
module DESPOTAlpha

using POMDPs
using BeliefUpdaters
using Parameters
using CPUTime
using ParticleFilters
using D3Trees
using Random
using Printf
using POMDPModelTools
using POMDPSimulators
using LinearAlgebra
using Statistics
using Distributions
using Plots
using Plots.PlotMeasures

using MCTS
using BasicPOMCP
import BasicPOMCP: SolvedFORollout, SolvedFOValue, convert_estimator

export
    DESPOTAlphaSolver,
    DESPOTAlphaPlanner,
    DESPOTAlphaTree,

    default_action,

    IndependentBounds,
    bounds,
    init_bounds,
    bound,
    init_bound,

    FOValue,
    FORollout,

    info_analysis,
    hist_analysis

"""
    DESPOTAlphaSolver(<keyword arguments>)

Each field may be set via keyword argument. The fields that correspond to algorithm
parameters match the definitions in the paper exactly.

# Fields
- `epsilon_0`
- `xi`
- `K`
- `D`
- `lambda`
- `T_max`
- `max_trials`
- `bounds`
- `default_action`
- `rng`
- `random_source`
- `tree_in_info`

Further information can be found in the field docstrings (e.g.
`?DESPOTAlphaSolver.xi`)
"""
@with_kw struct DESPOTAlphaSolver{R<:AbstractRNG} <: Solver
    "The target gap between the upper and the lower bound at the root of the DESPOTAlpha tree."
    epsilon_0::Float64                      = 1e-3

    "The minimum relative gap required for a branch to be expanded."
    xi::Float64                             = 0.95

    "The number of particles used for approximating beliefs"
    K::Int                                  = 100

    "Maximum number of observations per action sequences"
    C::Int                                  = 20

    "If true, pp_alpha_upper will be updated during backup, which is not necessary when the MDP upper bound is used."
    pp_update::Bool                         = false

    "The maximum depth of the belief tree."
    max_depth::Int                          = 90

    "The maximum online planning time per step."
    T_max::Float64                          = 1.0

    "The maximum number of trials of the planner."
    max_trials::Int                         = typemax(Int)

    "A representation for the upper and lower bound on the discounted value (e.g. `IndependentBounds`)."
    bounds::Any                             = IndependentBounds(-1e6, 1e6)

    """A default action to use if algorithm fails to provide an action because of an error.

    This can either be an action object, i.e. `default_action=1` if `actiontype(pomdp)==Int` or a function `f(pomdp, b, ex)` where b is the belief and ex is the exception that caused the planner to fail.
    """
    default_action::Any                     = ExceptionRethrow()

    "A random number generator for the internal sampling processes."
    rng::R                                  = MersenneTwister(rand(UInt32))

    "If true, a reprenstation of the constructed DESPOT is returned by POMDPModelTools.action_info."
    tree_in_info::Bool                      = false

    "Issue an warning when the planning time surpass the time limit by `timeout_warning_threshold` times"
    timeout_warning_threshold::Float64     = T_max * 2.0

    "Number of pre-allocated belief nodes"
    num_b::Int                              = 50_000
end

mutable struct DESPOTAlphaTree{S,A,O}
    # belief nodes
    weights::Vector{Vector{Float64}} # stores weights for *belief node*
    alpha::Vector{Vector{Float64}} # stores alpha vectors for each belief node
    children::Vector{UnitRange{Int}} # to children *ba nodes*
    parent::Vector{Int} # maps to the parent *ba node*
    Delta::Vector{Int}
    u::Vector{Float64}
    l::Vector{Float64}
    obs_prob::Vector{Float64}

    # action nodes
    ba_children::Vector{UnitRange{Int}}
    ba_parent::Vector{Int} # maps to parent *belief node*
    ba_u::Vector{Float64}
    ba_l::Vector{Float64}
    ba_action::Vector{A}

    # action sequences
    as_map::Dict{Vector{A}, Int}
    particles::Vector{Vector{S}}
    obs::Vector{Vector{O}}
    r::Vector{Vector{Float64}}
    L::Vector{Matrix{Float64}}
    def_alpha::Vector{Vector{Float64}}
    def_alpha_upper::Vector{Vector{Float64}}
    pp_alpha_upper::Vector{Vector{Float64}}

    root_particles::Vector{S}
    b::Int
    ba::Int
    as::Int
end

mutable struct DESPOTAlphaPlanner{S, A, O, P<:POMDP{S,A,O}, B, OD, RNG<:AbstractRNG} <: Policy
    sol::DESPOTAlphaSolver{RNG}
    pomdp::P
    bounds::B
    discounts::Vector{Float64}
    rng::RNG
    # The following attributes are used to avoid reallocating memory
    obs_ind_dict::Dict{O, Int}
    obs_dists::Vector{OD}
    tree::Union{Nothing, DESPOTAlphaTree{S,A,O}}
end

function DESPOTAlphaPlanner(sol::DESPOTAlphaSolver, pomdp::POMDP{S,A,O}) where {S,A,O}
    rng = deepcopy(sol.rng)
    bounds = init_bounds(sol.bounds, pomdp, sol, rng)
    discounts = discount(pomdp) .^[0:(sol.max_depth+1);]

    obs_dists = Vector{Union{typeof(observation(pomdp, first(actions(pomdp)), rand(initialstate(pomdp)))), Nothing}}(undef, sol.K)
    return DESPOTAlphaPlanner(deepcopy(sol), pomdp, bounds, discounts, rng,
                            Dict{O, Int}(), obs_dists, nothing)
end

solver(p::DESPOTAlphaPlanner) = p.sol

include("bounds.jl")
include("tree.jl")
include("planner.jl")
include("pomdps_glue.jl")
include("visualization.jl")
include("analysis.jl")

end # module


text:
abstract type AbstractDist{E} end

np_float(::AbstractDist{E}) where {E} = E
np_fixed(::AbstractDist{E}) where {E} = !E

include("dist_tanh.jl")

struct EqualSpacing{E} <: AbstractDist{E}
    np::Int
    dsmax::Float64
end

function EqualSpacing(np::Int=0; dsmax::Real=-1)
    E = np < 2
    if E
        dsmax <= 0 && "`dsmax`(> 0) must be specified if `np` < 2 (auto determin)"
    else
        dsmax <= 0 || @warn "`dsmax` is ignored since `np` has been specified"
    end
    return EqualSpacing{E}(np, dsmax)
end

isvalid(d::EqualSpacing{true}) = d.dsmax > 0
isvalid(d::EqualSpacing{false}) = d.np > 1

function nodes(d::EqualSpacing{false}; len::L=1.0) where {L<:Real}
    return collect(LinRange(0, len, d.np))
end

function nodes(::EqualSpacing, np::Int; len::L=1.0) where {L<:Real}
    return collect(LinRange(0, len, np))
end

function nodes(d::EqualSpacing{true}; len::L=1.0) where {L<:Real}
    np = ceil(Int, len / d.dsmax) + 1
    return nodes(d, np; len=len)
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

function get_flattened_index(node)
    @assert isa(node, Expr) && node.head == :ref
    idx = []
    if isexpr(node.args[1]) && node.args[1].head == :ref
        append!(idx, get_flattened_index(node.args[1]))
    end
    append!(idx, node.args[2:end])
    idx
    # filter((i) -> i != 1, idx)
end

function search_for_symbol(node, cbdata, index, top_level, read)
    if isa(node, Symbol) && node == cbdata[:expected]
        cbdata[:found] = true
    end
    ASTWALK_RECURSE
end

function contains_index(actual_idxs, expected_idxs)
    function check_contains_index(actual, expected)
        cbdata = Dict(:expected => expected, :found => false)
        AstWalk(actual, search_for_symbol, cbdata)
        cbdata[:found]
    end
    all(map(check_contains_index, actual_idxs, expected_idxs))
end

function get_target(A, offsets)
    @assert isa(A, Expr) && A.head == :ref
    idx = A.args[2:end]
    offset = 0
    if length(idx) > length(offsets)
        for i in length(idx):-1:length(offsets)+1
            if offset == 0
                offset = :(($(idx[i]) - 1) * size($(A.args[1]), $(i - 1)))
            else
                offset = :(($offset + ($(idx[i]) - 1)) * size($(A.args[1]), $(i - 1)))
            end
        end
    end
    for i in length(offsets):-1:1
        if offsets[i] != 1
            inner_size = 1
            for d in 1:i - 1
                inner_size = :(size($(A.args[1]), $d) * $inner_size)
            end
            offset = :(($offset + ($(offsets[i]) - 1)))
        end
        if i > 1 && offset != 0
            offset = :($offset * size($(A.args[1]), $(i - 1)))
        end
    end
    return :(pointer($(A.args[1]), $offset))
end
function collapse_blocks(node)
    body = []
    for arg in node.args
        if isa(arg, Expr) && arg.head == :block
            append!(body, collapse_blocks(arg))
        elseif arg != :nothing
            push!(body, arg)
        end
    end
    body
end

function isa_gemm_statement(statement)
    return isa(statement, Expr) &&
           statement.head == :(+=) &&
           isa(statement.args[2], Expr) &&
           statement.args[2].head == :call &&
           statement.args[2].args[1] == :(*) &&
           length(statement.args[2].args) == 3
end

function replace_index_with_expr(expr, idx, replace_with)
    function walker(node, cbdata, index, top_level, read)
        if isa(node, Symbol) && haskey(cbdata, node)
            return cbdata[node]
        end
        ASTWALK_RECURSE
    end
    cbdata = Dict(idx => replace_with)
    AstWalk(expr, walker, cbdata)
end

simple_ijk_orders = [(:i, :j, :k), (:j, :k, :i), (:j, :i, :k), (:k, :j, :i),
                     (:i, :k, :j), (:k, :i, :j)]
for (var1, var2, var3) in simple_ijk_orders
    fn = symbol(:pattern_match_gemm, var1, var2, var3)
    @eval function $fn(node, cbdata, index, top_level, read)
        @recurse_unless isa(node, Expr)
        @recurse_unless node.head == :for
        $(symbol(var1,:_loop)) = node
        $(symbol(var1,:_loopvar)) = get_loopvar(node)

        $(symbol(var2,:_loop)) = collapse_blocks($(symbol(var1,:_loop)).args[2])[1]
        @recurse_unless $(symbol(var2,:_loop)).head == :for
        $(symbol(var2,:_loopvar)) = get_loopvar($(symbol(var2,:_loop)))

        $(symbol(var3,:_loop)) = collapse_blocks($(symbol(var2,:_loop)).args[2])[1]
        @recurse_unless $(symbol(var3,:_loop)).head == :for
        $(symbol(var3,:_loopvar)) = get_loopvar($(symbol(var3,:_loop)))
        N = get_loop_length(i_loop)
        M = get_loop_length(j_loop)
        K = get_loop_length(k_loop)

        @recurse_unless length($(symbol(var3,:_loop)).args[2].args) == 1
        accum_stmt = $(symbol(var3,:_loop)).args[2].args[1]

        @recurse_unless isa_gemm_statement(accum_stmt)

        accum_target = accum_stmt.args[1]
        C = accum_target
        @recurse_unless length(get_flattened_index(accum_target)) >= 2
        C_idx = get_flattened_index(accum_target)
        @recurse_unless contains_index(C_idx[1:2], [j_loopvar, i_loopvar])

        i_start = get_loop_start(i_loop)
        j_start = get_loop_start(j_loop)
        k_start = get_loop_start(k_loop)

        (A, B) = accum_stmt.args[2].args[2:3]
        @recurse_unless length(get_flattened_index(A)) >= 2
        @recurse_unless length(get_flattened_index(B)) >= 2
        A_idx = get_flattened_index(A)[1:2]
        B_idx = get_flattened_index(B)[1:2]
        if contains_index(B_idx, [k_loopvar, i_loopvar])
            transB = 'N'
            B_offset = map(replace_index_with_expr, B_idx,[k_loopvar, i_loopvar], [k_start, i_start])
        elseif contains_index(B_idx, [i_loopvar, k_loopvar])
            transB = 'T'
            B_offset = map(replace_index_with_expr, B_idx,[i_loopvar, k_loopvar], [i_start, k_start])
        else
            return ASTWALK_RECURSE
        end
        if contains_index(A_idx, [j_loopvar, k_loopvar])
            transA = 'N'
            A_offset = map(replace_index_with_expr, A_idx, [j_loopvar, k_loopvar], [j_start, k_start])
        elseif contains_index(A_idx, [k_loopvar, j_loopvar])
            transA = 'T'
            A_offset = map(replace_index_with_expr, A_idx, [k_loopvar, j_loopvar], [k_start, j_start])
        else
            return ASTWALK_RECURSE
        end
        C_offset = map(replace_index_with_expr, C_idx[1:2], [j_loopvar, i_loopvar], [j_start, i_start])
        if length(C_idx) > 2 && C_idx[3] == k_loopvar &&
                contains(string(C.args[1]), "∇") &&
                length(split(string(C.args[1]), "∇")[2]) > 0
            push!(C_offset, 1)
        end
        return :(gemm!($transA, $transB, $M, $N, $K, 1.0f0,
                       $(get_target(A, A_offset)), size($(A.args[1]), 1),
                       $(get_target(B, B_offset)), size($(B.args[1]), 1), 1.0f0,
                       $(get_target(C, C_offset)), size($(C.args[1]), 1)))
    end
end

"""
The case of the "flattened" j loop
"""
function pattern_match_gemm4(node, cbdata, index, top_level, read)
    @recurse_unless isa(node, Expr)
    @recurse_unless node.head == :for
    i_loopvar = get_loopvar(node)
    i_loop = node

    j1_loop = collapse_blocks(node.args[2])[1]
    @recurse_unless j1_loop.head == :for
    j1_loopvar = get_loopvar(j1_loop)

    j2_loop = collapse_blocks(j1_loop.args[2])[1]
    @recurse_unless j2_loop.head == :for
    j2_loopvar = get_loopvar(j2_loop)

    k_loop = collapse_blocks(j2_loop.args[2])[1]
    @recurse_unless k_loop.head == :for
    k_loopvar = get_loopvar(k_loop)

    N = get_loop_length(i_loop)
    M1 = get_loop_length(j1_loop)
    M2 = get_loop_length(j2_loop)
    M = :($M1 * $M2)
    K = get_loop_length(k_loop)

    @recurse_unless length(k_loop.args[2].args) == 1
    accum_stmt = k_loop.args[2].args[1]

    @recurse_unless isa_gemm_statement(accum_stmt)

    accum_target = accum_stmt.args[1]
    C = accum_target
    @recurse_unless length(get_flattened_index(accum_target)) >= 3
    (j2, j1, i) = get_flattened_index(accum_target)[1:3]
    C_idx = [j2, j1, i]
    @recurse_unless contains_index(C_idx, [j2_loopvar, j1_loopvar, i_loopvar])

    i_start = get_loop_start(i_loop)
    j1_start = get_loop_start(j1_loop)
    j2_start = get_loop_start(j2_loop)
    k_start = get_loop_start(k_loop)

    (B, A) = accum_stmt.args[2].args[2:3]
    @recurse_unless length(get_flattened_index(A)) >= 3
    A_idx = get_flattened_index(A)[1:3]
    B_idx = get_flattened_index(B)[1:2]
    if contains_index(B_idx, [k_loopvar, i_loopvar])
        transB = 'N'
        B_offset = map(replace_index_with_expr, B_idx,[k_loopvar, i_loopvar], [k_start, i_start])
    elseif contains_index(B_idx, [i_loopvar, k_loopvar])
        transB = 'T'
        B_offset = map(replace_index_with_expr, B_idx,[i_loopvar, k_loopvar], [i_start, k_start])
    else
        return ASTWALK_RECURSE
    end
    if contains_index(A_idx, [j2_loopvar, j1_loopvar, k_loopvar])
        transA = 'N'
        lda = :(size($(A.args[1]), 1) * size($(A.args[1]), 2))
        A_offset = map(replace_index_with_expr, A_idx, [j2_loopvar, j1_loopvar, k_loopvar], [j2_start, j2_start, k_start])
    elseif contains_index(A_idx, [k_loopvar, j2_loopvar, j1_loopvar])
        transA = 'T'
        lda = :(size($(A.args[1]), 1))
        A_offset = map(replace_index_with_expr, A_idx, [k_loopvar, j2_loopvar, j1_loopvar], [k_start, j2_start, j2_start])
    else
        return ASTWALK_RECURSE
    end
    C_offset = map(replace_index_with_expr, C_idx, [j2_loopvar, j1_loopvar, i_loopvar], [j2_start, j1_start, i_start])
    return :(gemm!($transA, $transB, $M, $N, $K, 1.0f0,
             $(get_target(A, A_offset)), $lda,
             $(get_target(B, B_offset)), size($(B.args[1]), 1), 1.0f0,
             $(get_target(C, C_offset)), size($(C.args[1]), 1) * size($(C.args[1]), 2)))
end

"""
The case of the "flattened" k loop
"""
function pattern_match_gemm5(node, cbdata, index, top_level, read)
    @recurse_unless isa(node, Expr)
    @recurse_unless node.head == :for
    i_loopvar = get_loopvar(node)
    i_loop = node

    k1_loop = collapse_blocks(node.args[2])[1]
    @recurse_unless k1_loop.head == :for
    k1_loopvar = get_loopvar(k1_loop)

    k2_loop = collapse_blocks(k1_loop.args[2])[1]
    @recurse_unless k2_loop.head == :for
    k2_loopvar = get_loopvar(k2_loop)

    j_loop = collapse_blocks(k2_loop.args[2])[1]
    @recurse_unless j_loop.head == :for
    j_loopvar = get_loopvar(j_loop)

    N = get_loop_length(i_loop)
    M = get_loop_length(j_loop)
    K1 = get_loop_length(k1_loop)
    K2 = get_loop_length(k2_loop)
    K = :($K1 * $K2)

    @recurse_unless length(j_loop.args[2].args) == 1
    accum_stmt = j_loop.args[2].args[1]

    @recurse_unless isa_gemm_statement(accum_stmt)

    accum_target = accum_stmt.args[1]
    C = accum_target
    @recurse_unless length(get_flattened_index(accum_target)) >= 2
    C_idx = get_flattened_index(accum_target)[1:2]
    @recurse_unless contains_index(C_idx, [j_loopvar, i_loopvar])
    (j, i) = C_idx

    i_start = get_loop_start(i_loop)
    j_start = get_loop_start(j_loop)
    k1_start = get_loop_start(k1_loop)
    k2_start = get_loop_start(k2_loop)

    (A, B) = accum_stmt.args[2].args[2:3]
    @recurse_unless length(get_flattened_index(A)) >= 3
    @recurse_unless length(get_flattened_index(B)) >= 3
    A_idx = get_flattened_index(A)[1:3]
    B_idx = get_flattened_index(B)[1:3]
    if contains_index(B_idx, [k2_loopvar, k1_loopvar, i_loopvar])
        transB = 'N'
        ldb = :(size($(B.args[1]), 1) * size($(B.args[1]), 2))
        B_offset = map(replace_index_with_expr, B_idx, [k2_loopvar, k1_loopvar, i_loopvar], [k2_start, k1_start, i_start])
    elseif contains_index(B_idx, [i_loopvar, k2_loopvar, k1_loopvar])
        transB = 'T'
        ldb = :(size($(B.args[1]), 1))
        B_offset = map(replace_index_with_expr, B_idx, [i_loopvar, k2_loopvar, k1_loopvar], [i_start, k2_start, k1_start])
    else
        return ASTWALK_RECURSE
    end
    if contains_index(A_idx, [j_loopvar, k2_loopvar, k1_loopvar])
        transA = 'N'
        lda = :(size($(A.args[1]), 1))
        A_offset = map(replace_index_with_expr, A_idx, [j_loopvar, k2_loopvar, k1_loopvar], [j_start, k2_start, k2_start])
    elseif contains_index(A_idx, [k2_loopvar, k1_loopvar, j_loopvar])
        transA = 'T'
        lda = :(size($(A.args[1]), 1) * size($(A.args[1]), 2))
        A_offset = map(replace_index_with_expr, A_idx, [k2_loopvar, k1_loopvar, j_loopvar], [k2_start, k2_start, j_start])
    else
        return ASTWALK_RECURSE
    end
    C_offset = map(replace_index_with_expr, C_idx, [j_loopvar, i_loopvar], [j_start, i_start])
    return :(gemm!($transA, $transB, $M, $N, $K, 1.0f0,
             $(get_target(A, A_offset)), $lda,
             $(get_target(B, B_offset)), $ldb, 1.0f0,
             $(get_target(C, C_offset)), size($(C.args[1]), 1)))
end

function pattern_match_gemm6(node, cbdata, index, top_level, read)
    @recurse_unless isa(node, Expr)
    @recurse_unless node.head == :for
    k_loopvar = get_loopvar(node)
    k_loop = node

    i1_loop = collapse_blocks(node.args[2])[1]
    @recurse_unless i1_loop.head == :for
    i1_loopvar = get_loopvar(i1_loop)

    i2_loop = collapse_blocks(i1_loop.args[2])[1]
    @recurse_unless i2_loop.head == :for
    i2_loopvar = get_loopvar(i2_loop)

    j_loop = collapse_blocks(i2_loop.args[2])[1]
    @recurse_unless j_loop.head == :for
    j_loopvar = get_loopvar(j_loop)

    N1 = get_loop_length(i1_loop)
    N2 = get_loop_length(i2_loop)
    N = :($N1 * $N2)
    M = get_loop_length(j_loop)
    K = get_loop_length(k_loop)

    @recurse_unless length(j_loop.args[2].args) == 1
    accum_stmt = j_loop.args[2].args[1]

    @recurse_unless isa_gemm_statement(accum_stmt)

    accum_target = accum_stmt.args[1]
    C = accum_target
    @recurse_unless length(get_flattened_index(accum_target)) >= 3
    C_idx = get_flattened_index(accum_target)[1:3]
    @recurse_unless contains_index(C_idx, [j_loopvar, i2_loopvar, i1_loopvar])

    i1_start = get_loop_start(i1_loop)
    i2_start = get_loop_start(i2_loop)
    j_start = get_loop_start(j_loop)
    k_start = get_loop_start(k_loop)

    (A, B) = accum_stmt.args[2].args[2:3]
    @recurse_unless length(get_flattened_index(B)) >= 3
    A_idx = get_flattened_index(A)[1:2]
    B_idx = get_flattened_index(B)[1:3]
    if contains_index(B_idx, [k_loopvar, i2_loopvar, i1_loopvar])
        transB = 'N'
        ldb = :(size($(B.args[1]), 1))
        B_offset = map(replace_index_with_expr, B_idx, [k_loopvar, i2_loopvar, i1_loopvar], [k_start, i2_start, i1_start])
    elseif contains_index(B_idx, [i2_loopvar, i1_loopvar, k_loopvar])
        transB = 'T'
        ldb = :(size($(B.args[1]), 1) * size($(B.args[1]), 2))
        B_offset = map(replace_index_with_expr, B_idx, [i2_loopvar, i1_loopvar, k_loopvar], [i2_start, i1_start, k_start])
    else
        return ASTWALK_RECURSE
    end
    if contains_index(A_idx, [j_loopvar, k_loopvar])
        transA = 'N'
        A_offset = map(replace_index_with_expr, A_idx, [j_loopvar, k_loopvar], [j_start, k_start])
    elseif contains_index(A_idx, [k_loopvar, j_loopvar])
        transA = 'T'
        A_offset = [k_start, j_start]
        A_offset = map(replace_index_with_expr, A_idx, [k_loopvar, j_loopvar], [k_start, j_start])
    else
        return ASTWALK_RECURSE
    end
    C_offset = map(replace_index_with_expr, C_idx, [j_loopvar, i2_loopvar, i1_loopvar], [j_start, i2_start, i1_start])
    return :(gemm!($transA, $transB, $M, $N, $K, 1.0f0,
             $(get_target(A, A_offset)), size($(A.args[1]), 1),
             $(get_target(B, B_offset)), $ldb, 1.0f0,
             $(get_target(C, C_offset)), size($(C.args[1]), 1)))
end


text:
# ---
# title: 1631. Path With Minimum Effort
# id: problem1631
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Binary Search, Depth-first Search, Union Find, Graph
# link: <https://leetcode.com/problems/path-with-minimum-effort/description/>
# hidden: true
# ---
#
# You are a hiker preparing for an upcoming hike. You are given `heights`, a 2D
# array of size `rows x columns`, where `heights[row][col]` represents the
# height of cell `(row, col)`. You are situated in the top-left cell, `(0, 0)`,
# and you hope to travel to the bottom-right cell, `(rows-1, columns-1)` (i.e.,
# **0-indexed** ). You can move **up** , **down** , **left** , or **right** ,
# and you wish to find a route that requires the minimum **effort**.
#
# A route's **effort** is the **maximum absolute difference** **** in heights
# between two consecutive cells of the route.
#
# Return _the minimum **effort** required to travel from the top-left cell to
# the bottom-right cell._
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/04/ex1.png)
#
#
#
#     Input: heights = [[1,2,2],[3,8,2],[5,3,5]]
#     Output: 2
#     Explanation: The route of [1,3,5,3,5] has a maximum absolute difference of 2 in consecutive cells.
#     This is better than the route of [1,2,2,2,5], where the maximum absolute difference is 3.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/04/ex2.png)
#
#
#
#     Input: heights = [[1,2,3],[3,8,4],[5,3,5]]
#     Output: 1
#     Explanation: The route of [1,2,3,4,5] has a maximum absolute difference of 1 in consecutive cells, which is better than route [1,3,5,3,5].
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/04/ex3.png)
#
#
#
#     Input: heights = [[1,2,1,1,1],[1,2,1,2,1],[1,2,1,2,1],[1,2,1,2,1],[1,1,1,2,1]]
#     Output: 0
#     Explanation: This route does not require any effort.
#
#
#
#
# **Constraints:**
#
#   * `rows == heights.length`
#   * `columns == heights[i].length`
#   * `1 <= rows, columns <= 100`
#   * `1 <= heights[i][j] <= 106`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using GeoInterface
using Test

@testset "Comparison operators" begin
    # Points are now compared by value:
    pt1, pt2, pt3 = Point([0.0, 0.0]), Point([0.0, 0.0]), Point([0.0, 1.0])
    @test pt1 == pt1
    @test pt1 == pt2
    @test pt1 != pt3

    # Also, the hash is based on the coordinates
    @test hash(pt1) == hash(pt1)
    @test hash(pt1) == hash(pt2)
    @test hash(pt1) != hash(pt3)

    # Implicitly, this should also work for `isequal`:
    @test isequal(pt1, pt1)
    @test isequal(pt1, pt2)
    @test !isequal(pt1, pt3)

    # Can also do approximate comparisons
    pt4 = Point([0.0, 1.001])
    @test pt3 != pt4
    @test pt3 ≉ pt4  atol=0.0001
    @test pt3 ≉ pt4  rtol=0.0001
    @test pt3 ≈ pt4  atol=0.001
    @test pt3 ≈ pt4  rtol=0.001

    # The same is not true for other geometry types: The representation is not
    # unique, so comparing the coordinates directly might be misleading.
    pg1 = Polygon([[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]])
    pg2 = Polygon([[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]])
    pg3 = Polygon([[[1.0, 0.0], [0.0, 1.0], [0.0, 0.0]]])
    @test pg1 == pg1   # same objects
    @test pg1 != pg2   # same values, but different objects
    @test pg1 != pg3   # equivalent, but not same values
end


text:
export Vandermonde
"""
[`Vandermonde` matrix](http://en.wikipedia.org/wiki/Vandermonde_matrix)

```julia
julia> a = collect(1.0:5.0)
julia> A = Vandermonde(a)
5×5 Vandermonde{Float64}:
 1.0  1.0   1.0    1.0    1.0
 1.0  2.0   4.0    8.0   16.0
 1.0  3.0   9.0   27.0   81.0
 1.0  4.0  16.0   64.0  256.0
 1.0  5.0  25.0  125.0  625.0
```

Adjoint Vandermonde:
```julia
julia> A'
5×5 LinearAlgebra.Adjoint{Float64,Vandermonde{Float64}}:
 1.0   1.0   1.0    1.0    1.0
 1.0   2.0   3.0    4.0    5.0
 1.0   4.0   9.0   16.0   25.0
 1.0   8.0  27.0   64.0  125.0
 1.0  16.0  81.0  256.0  625.0
```

The backslash operator \\ is overloaded to solve Vandermonde and adjoint
Vandermonde systems in `O(n^2)` time using the algorithm of
Björck & Pereyra (1970), https://doi.org/10.2307/2004623.

```julia
julia> A\a
5-element Array{Float64,1}:
 0.0
 1.0
 0.0
 0.0
 0.0

julia> r2 = A[2,:]
julia> A'\r2
5-element Array{Float64,1}:
 0.0
 1.0
 0.0
 0.0
 0.0
```
"""
struct Vandermonde{T} <: AbstractArray{T, 2}
    c :: Vector{T}
end

getindex(V::Vandermonde, i::Int, j::Int) = V.c[i]^(j-1)
isassigned(V::Vandermonde, i::Int, j::Int) = isassigned(V.c, i)

size(V::Vandermonde, r::Int) = (r==1 || r==2) ? length(V.c) :
    throw(ArgumentError("Invalid dimension $r"))
size(V::Vandermonde) = length(V.c), length(V.c)

function Matrix(V::Vandermonde{T}) where T
    n=size(V, 1)
    M=Array{T}(undef, n, n)
    M[:,1] .= 1
    for j=2:n
        for i=1:n
	    M[i,j] = M[i,j-1]*V.c[i]
        end
    end
    M
end

function Matrix(V::Adjoint{T,Vandermonde{T}}) where T
    n=size(V, 1)
    M=Array{T}(undef, n, n)
    M[1,:] .= 1
    for j=1:n
        for i=2:n
	    M[i,j] = M[i-1,j]*adjoint(V.parent.c[j])
        end
    end
    M
end

function Matrix(V::Transpose{T,Vandermonde{T}}) where T
    n=size(V, 1)
    M=Array{T}(undef, n, n)
    M[1,:] .= 1
    for j=1:n
        for i=2:n
	    M[i,j] = M[i-1,j]*V.parent.c[j]
        end
    end
    M
end

function \(V::Adjoint{T1,Vandermonde{T1}}, y::AbstractVecOrMat{T2}) where T1 where T2
    T = vandtype(T1,T2)
    x = Array{T}(undef, size(y))
    copyto!(x, y)
    pvand!(adjoint(V.parent.c), x)
    return x
end

function \(V::Transpose{T1,Vandermonde{T1}}, y::AbstractVecOrMat{T2}) where T1 where T2
    T = vandtype(T1,T2)
    x = Array{T}(undef, size(y))
    copyto!(x, y)
    pvand!(V.parent.c, x)
    return x
end

function \(V::Vandermonde{T1}, y::AbstractVecOrMat{T2}) where T1 where T2
    T = vandtype(T1,T2)
    x = Array{T}(undef, size(y))
    copyto!(x, y)
    dvand!(V.c, x)
    return x
end

function vandtype(T1::Type, T2::Type)
    # Figure out the return type of Vandermonde{T1} \ Vector{T2}
    T = promote_type(T1, T2)
    S = typeof(oneunit(T)/oneunit(T1))
    return S
end

"""
    pvand!(a, b) -> b

Solves system ``A^T*x = b`` in-place.

``A^T`` is transpose of Vandermonde matrix ``A_{ij} = a_i^{j-1}``.

Algorithm by Bjorck & Pereyra,
Mathematics of Computation, Vol. 24, No. 112 (1970), pp. 893-903,
https://doi.org/10.2307/2004623
"""
function pvand!(alpha, B)
    n = length(alpha);
    if n != size(B,1)
        throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
    end
    nrhs = size(B,2)
    @inbounds begin
        for j=1:nrhs
            for k=1:n-1
                for i=n:-1:k+1
                    B[i,j] = B[i,j]-alpha[k]*B[i-1,j]
                end
            end
            for k=n-1:-1:1
                for i=k+1:n
                    B[i,j] = B[i,j]/(alpha[i]-alpha[i-k])
                end
                for i=k:n-1
                    B[i,j] = B[i,j]-B[i+1,j]
                end
            end
        end
    end
end
"""
    dvand!(a, b) -> b

Solves system ``A*x = b`` in-place.

``A`` is Vandermonde matrix ``A_{ij} = a_i^{j-1}``.

Algorithm by Bjorck & Pereyra,
Mathematics of Computation, Vol. 24, No. 112 (1970), pp. 893-903,
https://doi.org/10.2307/2004623
"""
function dvand!(alpha, B)
    n = length(alpha)
    if n != size(B,1)
        throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
    end
    nrhs = size(B,2)
    @inbounds begin
        for j=1:nrhs
            for k=1:n-1
                for i=n:-1:k+1
                    B[i,j] = (B[i,j]-B[i-1,j])/(alpha[i]-alpha[i-k])
                end
            end
            for k=n-1:-1:1
                for i=k:n-1
                    B[i,j] = B[i,j]-alpha[k]*B[i+1,j]
                end
            end
        end
    end
end


text:
"""
TableData(c; half=false, level=5)

Makes LaTeX tables for the paper
"""
function TableData()
cvals=[.5, .99, 1.0]
for ic=1:3
    MakeTable(cvals(ic))
end
end

function MakeTable(c,tex=true)
fname=Fname(c,Float16)
c==1 ? (maxit=30) : (maxit=10)
DataC=zeros(maxit+1,5,2)
tcap=string("Half Precision Computed Convergence Rates: c=",string(c))
readmpdata(fname,DataC)
for it=1:2
for iz=1:5
DataC[:,iz,it] = DataC[:,iz,it]./DataC[1,iz,it];
end
end
qresid=DataC[:,:,1]
qrate=qresid[2:maxit+1,:]./qresid[1:maxit,:]
itnum=collect(Int,1:maxit);
Rates=[itnum qrate]
if tex
headers = ["n", "1024", "2048", "4096", "8192", "16384"]
formats = "%d & %9.5e & %9.5e & %9.5e & %9.5e & %9.5e"
fprintTeX(headers, formats, Rates)
else
headers = ["1024", "2048", "4096", "8192", "16384"]
printhist(qrate,headers)
end
end

"""
fprintTeX(headers,formats,data)

Print a LaTeX table from a Julia array.

Inputs:
headers: the titles for the columns
          example: headers=["foo", "bar"]
formats: c-style formatting for the columns.
          fprintTeX will add the carriage returns for you.
          example: formats="%d & %7.2e";
"""
function fprintTeX(headers, formats, data)
    (mr, mc) = size(data)
    @printf("\\begin{tabular}{")
    for i = 1:mc
        @printf("l")
    end
    @printf("} \n")
    for i = 1:mc-1
        @printf("%9s &", headers[i])
    end
    @printf("%9s \\\\ \n" , headers[mc])
    @printf("\\hline \n")
    #
    # I am not sure why @printf needs this, but it does.
    # See https://github.com/JuliaLang/julia/issues/4248
    #
    printf(fmt::String, args...) = @eval @printf($fmt, $(args...))
    #
    bigform = string(formats, "   \\\\ \n")
    for i = 1:mr
        printf(bigform, data[i, :]...)
    end
    @printf("\\hline \n")
    @printf("\\end{tabular} \n")
end

"""
printhist(tablein,headers;TeX=false,figures=5)

You are not supposed to look at this code. It is mostly repulsive
bookkeeping.

Inputs:
tablein: columns of equal length with the data
         Pad the short ones with NaNs. My codes which all this
         do that. The formatting turns the NaNs into spaces of
         the correct length.
headers: the titles for the columns
         example: headers=["foo", "bar"]
         Do not add the header for the iteration counter. I do that.
TeX:     Format the tabel in LaTeX or not?
figures: The printf format will be (figures+7).figures e
         So if figures = 5 the data will be formatted in %12e5.
         This makes everything line up.

Compare iteration histories from SIAMFANLEquations family of solvers.
This is usually called by something else.
Nothing to see here. Move along.
"""
function printhist(tablein, headers; TeX = false, figures = 5)
    ntab = length(tablein[1, :])
    bighead = Array{String,2}(undef, 1, ntab + 1)
    bighead[1] = "n"
    for ih = 2:ntab+1
        bighead[ih] = headers[ih-1]
    end
    if ntab > 5
        error("Too many columns for the table. Use fewer.")
    end
    fmtout = buildformat(ntab, TeX, figures)
    tabfmt = fmtout.format
    nanspace = fmtout.nanspace
    headerfmt = fmtout.headerfmt
    printf(fmt::String, args...) = @eval @printf($fmt, $(args...))
    sprintf(fmt::String, args...) = @eval @sprintf($fmt, $(args...))
    itmax = length(tablein[:, 1])
    itc = 0:itmax-1
    if TeX
        @printf("\\begin{tabular}{")
        for i = 1:ntab+1
            @printf("l")
        end
        @printf("} \n")
        printf(headerfmt, bighead...)
    else
        printf(headerfmt, bighead...)
    end
    for it = 1:itmax
        st = sprintf(tabfmt, it, tablein[it, :]...)
        snan = findfirst(isequal('N'), st)
        lt = length(st)
        while typeof(snan) != Nothing
            st = string(st[1:snan-1], nanspace, st[snan+3:lt])
            snan = findfirst(isequal('N'), st)
        end
        printf("%s", st)
    end
    if TeX
        @printf("\\hline \n")
        @printf("\\end{tabular} \n")
    end
end

function buildformat(ntab, TeX, figures)
    format = "%2d "
    nanspace = "   "
    fwide = figures + 7
    fm1 = string(fwide)
    fm2 = string(figures)
    basefmt = string(" %", fm1, ".", fm2, "e ")
    hfmt = string(" %", fm1, "s")
    headerfmt = " %s"
    if TeX
        for i = 1:ntab
            format = string(format, "&", basefmt)
            headerfmt = string(headerfmt, "&", hfmt)
        end
        format = string(format, " \\\\ \n")
        headerfmt = string(headerfmt, " \\\\ \\hline \n")
    else
        for i = 1:ntab
            format = string(format, basefmt)
            headerfmt = string(headerfmt, hfmt)
        end
        format = string(format, " \n")
        headerfmt = string(headerfmt, " \n \n")
    end
    return (format = format, nanspace = nanspace, headerfmt = headerfmt)
end

max_stars_repo_name: ctkelley/MPResults