text:
using Revise
using Test
using LinearAlgebra
using MFGnet
using Printf

d = 2
nex = 4
nTh = 3

Rs = [Float64;Float32]

for k=1:2
    @testset "derivative check $(Rs[k])" begin
    R = Rs[k]
    L = SingleLayer()
    N = ResNN(L,Vector(range(R(0.0),stop=R(1.0),length=5)))

    s = randn(R,d,nex)
    K = randn(R,d,d,nTh)
    b = randn(R,d,nTh)
    Θ = (K,b)
    w = randn(R,d)
    z = N(s,Θ)
    @test eltype(z) == R
    Φ = w'*z
    ds = randn(R,size(s))
    dΦ,d2Φ = getGradAndHessian(N,w,s,Θ)
    @test eltype(dΦ) == R
    @test eltype(d2Φ)==R
    dΦds = sum(ds.*dΦ,dims=1)

    curv = zeros(1,nex)
    for j=1:nex
        curv[1,j] = ds[:,j]'*d2Φ[:,:,j]*ds[:,j]
    end
    trH = MFGnet.getTraceHess(N,s,Θ)
    trHt, = MFGnet.getTraceHessAndGrad(N,s,Θ)
    @test norm(trH-trHt)/norm(trH) < sqrt(eps(R))
    for j=1:nex
        @test abs(trH[j] - tr(d2Φ[:,:,j]))/abs(tr(d2Φ[:,:,j])) < sqrt(eps(R))
    end

    for j=1:15
        h = R(2.0^(-j))
        Φt = w'*N(s + h .* ds,Θ)

        E0 = norm(Φ-Φt)
        E1 = norm(Φ+h*dΦds - Φt)
        E2 = norm(Φ+h*dΦds+0.5*h^2*curv - Φt)

        @printf("h=%1.3e\t\tE0=%1.3e\tE1=%1.3e\tE1=%1.3e\n",h,E0,E1,E2)
    end
    @test true
    end
end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
"""
    struct DEsopoPape <: ShortestPathAlgorithm

The structure used to configure and specify that [`shortest_paths`](@ref)
should use the [D'Esopo-Pape algorithm](http://web.mit.edu/dimitrib/www/SLF.pdf).
No fields are specified or required.

### Implementation Notes
`DEsopoPape` supports the following shortest-path functionality:
- non-negative distance matrices / weights
- all destinations
"""
struct DEsopoPape <: ShortestPathAlgorithm end
struct DEsopoPapeResult{T, U<:Integer} <: ShortestPathResult
    parents::Vector{U}
    dists::Vector{T}
end

function shortest_paths(g::AbstractGraph, src::Integer, distmx::AbstractMatrix, ::DEsopoPape)
    T = eltype(distmx)
    U = eltype(g)
    nvg = nv(g)
    (src in 1:nvg) || throw(DomainError(src, "src should be in between 1 and $nvg"))
    dists = fill(typemax(T), nvg)
    parents = zeros(U, nvg)
    state = Vector{Int8}()
    state = fill(Int8(2), nvg)
    q = U[src]
    @inbounds dists[src] = 0

    @inbounds while !isempty(q)
        u = popfirst!(q)
        state[u] = 0

        for v in outneighbors(g, u)
            alt = dists[u] + distmx[u, v]
            if (dists[v] > alt)
                dists[v] = alt
                parents[v] = u

                if state[v] == 2
                    state[v] = 1
                    push!(q, v)
                elseif state[v] == 0
                    state[v] = 1
                    pushfirst!(q, v)
                end
            end
        end
    end

    return DEsopoPapeResult{T, U}(parents, dists)
end

ext: jl


text:
module FDMTests

include("mpiexec.jl")
run_mpi_driver(procs=8,file="driver_fdm.jl")

end # module


text:
using Serialization
using Distributed
import Distributed.pmap

export runMC

@doc """
    runMC(param::Parameter)
    runMC(params::AbstractArray{Parameter}
          ;
          parallel::Bool=false,
          autoID::Bool=true)

Runs Monte Carlo simulation(s) and returns calculated observables.

# Restart

If a checkpoint file named `"\$(param["Checkpoint Filename Prefix"])_\$(param["ID"]).dat"` exists and
`param["Checkpoint Interval"] > 0.0`, `runMC` loads this file and restarts the pending simulation.
NOTE: Restart will fail if the version or the system image of julia change (see the doc of `Serialization.serialize` ).

# Keyward arguments
- `autoID`: If true, `"ID"`s will be set (overwritten) as `params[i]["ID"] = i`.
- `parallel`: If true, runs simulations in parallel (uses `pmap` instead of `map`).

# Required keys in `param`
- "Model"
- "Lattice"
- "Update Method"

# Optional keys in `param`
- "MCS": The number of Monte Carlo steps after thermalization
    - Default: `8192`
- "Thermalization": The number of Monte Carlo steps for thermalization
    - Default: `MCS>>3`
- "Seed": The initial seed of the random number generator, `MersenneTwister`
    - Default: determined randomly (see the doc of `Random.seed!`)
- "Checkpoint Filename Prefix": See the "Restart" section.
    - Default: `"cp"`
- "ID": See the "Restart" section.
    - Default: `0`
- "Checkpoint Interval": Time interval between writing checkpoint file in seconds.
    - Default: `0.0`, this means that NO checkpoint file will be loaded and saved.
"""
function runMC(params::AbstractArray{T}; parallel::Bool=false, autoID::Bool=true) where T<:Dict
    map_fn = ifelse(parallel, pmap, map)
    return map_fn(enumerate(params)) do (id,p)
        if autoID
            p["ID"] = id
        end
        return runMC(p)
    end
end

function runMC(param::Parameter)
    model = param["Model"](param)
    if "Seed" in keys(param)
        seed!(model, param["Seed"])
    end
    ret = runMC(model, param)
    return ret
end

function runMC(model, param::Parameter)
    MODEL = typeof(model)
    verbose = get(param, "Verbose", false) :: Bool
    if verbose
        println("Start: ", param)
    end
    cp_filename = @sprintf("%s_%d.dat",
                           get(param, "Checkpoint Filename Prefix", "cp")::String,
                           get(param, "ID", 0)::Int)
    cp_interval = get(param, "Checkpoint Interval", 0.0) :: Float64
    tm = time()

    MCS = get(param, "MCS", 8192) :: Int
    Therm = get(param, "Thermalization", MCS>>3) :: Int

    mcs = 0
    MCS += Therm
    obs = BinningObservableSet()
    makeMCObservable!(obs, "Time per MCS")
    makeMCObservable!(obs, "MCS per Second")

    if cp_interval > 0.0 && ispath(cp_filename)
        open(cp_filename) do io
            model = deserialize(io)
            obs = deserialize(io)
            mcs = deserialize(io)
        end
    end

    update! = param["Update Method"] :: Function
    if haskey(param, "Estimator")
        estimator = param["Estimator"] :: Function
    else
        estimator = default_estimator(model, update!)
    end
    pp = get(param, "Post Process", postproc) :: Function
    p = convert_parameter(model, param)

    while mcs < MCS
        if mcs < Therm
            update!(model, p...)
        else
            t = @elapsed begin
                st = update!(model,p...)
                localobs = estimator(model, p..., st)
            end
            obs["Time per MCS"] << t
            obs["MCS per Second"] << (1.0/t)
            accumulateObservables!(model, obs, localobs)
        end
        mcs += 1
        if cp_interval > 0.0 && time() - tm > cp_interval
            open(cp_filename, "w") do io
                serialize(io, model)
                serialize(io, obs)
                serialize(io, mcs)
            end
            tm += cp_interval
         end
    end

    if cp_interval > 0.0
        open(cp_filename, "w") do io
            serialize(io, model)
            serialize(io, obs)
            serialize(io, mcs)
        end
    end

    jk = pp(model, param, obs)

    if verbose
        println("Finish: ", param)
    end
    return jk
end

@doc """
    accumulateObservables!(model, obs::MCObservableSet, localobs::Dict)

Accumulates `localobs` into `obs`. For example, `obs["Energy"] << localobs["Energy"]`.
"""
function accumulateObservables!(::Model, obs::MCObservableSet, localobs::Measurement)
    if length(obs) < 3
        @inbounds for key in keys(localobs)
            makeMCObservable!(obs, key)
            obs[key] << localobs[key]
        end
    else
        @inbounds for key in keys(localobs)
            obs[key] << localobs[key]
        end
    end
    return obs
end


text:
using BNP_WMReg_Joint

using BSON
import Clustering.hclust, Clustering.cutree
import Distances.pairwise, Distances.Euclidean
import BayesInference.embed, BayesInference.etr, BayesInference.logsumexp
using PDMats
using StatsBase
using Statistics
using Random
using Dates

using LinearAlgebra
BLAS.set_num_threads(1)

const whichsim = ARGS[1]
const whichfit = ARGS[2]

BSON.@load "data/$(whichsim).bson" y X y_valid X_valid Y_valid ###### set directory correctly
BSON.@load "postsimB/postsim_$(whichsim)_$(whichfit).bson"  model sims mesg ###### set directory correctly

outfile = "KL/KL_BNP_WMAR_$(mesg).txt"

nprime = length(y_valid)
all( y[1:model.n] .== model.y ) || throw("Data mismatch.")

nrep = size(Y_valid, 1)
nsim = length(sims)

L = deepcopy(model.K)
H = deepcopy(model.H)

### Calculate llik ordinate for each Y_valid

using Distributions

function β_y_modify_γ(β_y::Array{T, 2}, γ::BitArray{1}) where T <: Real
    modify_indx = findall(.!(γ))
    βout = deepcopy(β_y)

    for k in modify_indx # this works even if no modifications are necessary
        βout[:,k] *= 0
    end
    return βout
end
function β_y_modify_γ(β_y::Array{T, 2}, γ::BitArray{2}) where T <: Real
    size(β_y) == size(γ) || throw("β and γ size mismatch.")
    return β_y .* γ
end

function ldnorm(y::Float64, mu::Float64, sig2::Float64)
    -0.5 * log(2.0π * sig2) - (0.5 * (y - mu)^2 / sig2)
end

ldw_mat = ldensweight_mat(X_valid[1:nprime, 1:L], sims[1:nsim], model.state.γδc)

ldens00 = Array{Float64,2}(undef, nrep, H)
ldens0 = Array{Float64,1}(undef, nprime)
fill!(ldens0, 0.0)

for j = 1:nprime

    for ii = 1:nsim
        βγ_y = β_y_modify_γ(sims[ii][:β_y], sims[ii][:γ])

        for h = 1:H
            Ey_h = sims[ii][:μ_y][h] - sum( βγ_y[h,:] .* (X_valid[j,1:L] - sims[ii][:μ_x][h,:]) )
            for jj = 1:nrep
                # ldens00[jj,h] = ldw_mat[ii, j, h] + logpdf(Normal(Ey_h, sqrt(sims[ii][:δ_y][h])), Y_valid[jj,j])
                ldens00[jj,h] = ldw_mat[ii, j, h] + ldnorm(Y_valid[jj,j], Ey_h, sims[ii][:δ_y][h])
            end
        end

        ldens0[j] += mean( logsumexp(ldens00, 2) ) # mean over replicates jj

    end
    if j % 50 == 0
        report_file = open(outfile, "a+")
        write(report_file, "Calculating density ordinates. Obs $(j) of $(nprime) at $(Dates.now()).\n")
        close(report_file)
    end
    # println("j = ", j, " of ", nprime)
end

B_vec = ldens0 ./ float(nsim) # posterior mean

### Compute and aggregate KL
using Distributions

if whichsim == "Ricker_single"
    function ldens_true(y::Float64, x::Vector{Float64}, theta::Float64, sig::Float64)
            μ = x[2] * exp( theta - x[2] )
            return logpdf(Normal(μ, sig), y)
    end

elseif whichsim == "Ricker_single_lognormal"
    function ldens_true(y::Float64, x::Vector{Float64}, theta::Float64, sig::Float64)
            μ = log(x[2]) + theta - x[2]
            return logpdf(LogNormal(μ, sig), y)
    end

elseif whichsim == "Ricker_twolag_lognormal"
    function ldens_true(y::Float64, x::Vector{Float64}, theta::Float64, sig::Float64)
            μ = log(x[2]) + theta - x[2]
            σ = sig * x[1]
            return logpdf(LogNormal(μ, σ), y)
    end

end

ld_true = [ ldens_true(Y_valid[jj,j], X_valid[j,:], 2.6, 0.09) for jj = 1:nrep, j = 1:nprime ]

A_vec = mean(ld_true, dims=1)[1,:]
# plot(A_vec)

# using Plotly
# trace_true = scatter(Dict("x" => 1:nprime,
#     "y" => A_vec, "name" => "data",
#     "mode" => "markers", "type" => "scatter", "marker" => Dict("color" => "blue")))
# trace_mod = scatter(Dict("x" => 1:nprime,
#     "y" => B_vec, "name" => "validation",
#     "mode" => "markers", "type" => "scatter", "marker" => Dict("color" => "red")))
# plot([trace_true, trace_mod])

## sanity check
# trace = scatter(Dict("x" => X_valid[1:nprime,2],
#     "y" => mean(Y_valid[:,1:nprime], dims=1)[1,:], "name" => "data",
#     "mode" => "markers", "type" => "scatter", "marker" => Dict("color" => "blue")))
# trace = scatter(Dict("x" => X_valid[1:nprime,2],
#     "y" => Y_valid[1,1:nprime], "name" => "data",
#     "mode" => "markers", "type" => "scatter", "marker" => Dict("color" => "blue")))
# plot([trace])

# plot(A_vec - B_vec)

KL_valid = A_vec .- B_vec

KL_out = mean(KL_valid)

report_file = open(outfile, "a+")

### Break out over certain regions by indexing KL_valid
if whichsim == "Ricker_single_lognormal"

    x2_cutoff = 2.0

    low_x2_bool = X_valid[1:nprime, 2] .<= x2_cutoff
    low_x2fit_bool = X[1:model.n, 2] .<= x2_cutoff

    n_low2 = sum(low_x2_bool)
    n_high2 = sum(.!low_x2_bool)

    n_low2_fit = sum(low_x2fit_bool)
    n_high2_fit = sum(.!low_x2fit_bool)

    write(report_file, "\n\nX2 below $(x2_cutoff), n_fit: $(n_low2_fit), n_valid: $(n_low2), mean KL: $(mean(KL_valid[low_x2_bool]))\n")
    write(report_file, "\nX2 above $(x2_cutoff), n_fit: $(n_high2_fit), n_valid: $(n_high2), mean KL: $(mean(KL_valid[.!low_x2_bool]))\n")

elseif whichsim == "Ricker_twolag_lognormal"

    x1_cutoff = 2.0
    x2_cutoff = 3.0

    low_x1_bool = X_valid[1:nprime, 1] .<= x1_cutoff
    low_x2_bool = X_valid[1:nprime, 2] .<= x2_cutoff

    low_x1fit_bool = X[1:model.n, 1] .<= x1_cutoff
    low_x2fit_bool = X[1:model.n, 2] .<= x2_cutoff

    low1_low2_bool = [ low_x1_bool[j] && low_x2_bool[j] for j = 1:nprime ]
    low1_high2_bool = [ low_x1_bool[j] && !low_x2_bool[j] for j = 1:nprime ]
    high1_low2_bool = [ !low_x1_bool[j] && low_x2_bool[j] for j = 1:nprime ]
    high1_high2_bool = [ !low_x1_bool[j] && !low_x2_bool[j] for j = 1:nprime ]

    low1_low2_fit_bool = [ low_x1fit_bool[j] && low_x2fit_bool[j] for j = 1:model.n ]
    low1_high2_fit_bool = [ low_x1fit_bool[j] && !low_x2fit_bool[j] for j = 1:model.n ]
    high1_low2_fit_bool = [ !low_x1fit_bool[j] && low_x2fit_bool[j] for j = 1:model.n ]
    high1_high2_fit_bool = [ !low_x1fit_bool[j] && !low_x2fit_bool[j] for j = 1:model.n ]

    n_low1 = sum(low_x1_bool)
    n_high1 = sum(.!low_x1_bool)

    n_low2 = sum(low_x2_bool)
    n_high2 = sum(.!low_x2_bool)

    n_low1_low2 = sum( low1_low2_bool )
    n_low1_high2 = sum( low1_high2_bool )
    n_high1_low2 = sum( high1_low2_bool )
    n_high1_high2 = sum( high1_high2_bool )

    n_low1_fit = sum(low_x1fit_bool)
    n_high1_fit = sum(.!low_x1fit_bool)

    n_low2_fit = sum(low_x2fit_bool)
    n_high2_fit = sum(.!low_x2fit_bool)

    n_low1_low2_fit = sum( low1_low2_fit_bool )
    n_low1_high2_fit = sum( low1_high2_fit_bool )
    n_high1_low2_fit = sum( high1_low2_fit_bool )
    n_high1_high2_fit = sum( high1_high2_fit_bool )

    write(report_file, "\n\nX1 below $(x1_cutoff), n_fit: $(n_low1_fit), n_valid: $(n_low1), mean KL: $(mean(KL_valid[low_x1_bool]))\n")
    write(report_file, "\nX1 above $(x1_cutoff), n_fit: $(n_high1_fit), n_valid: $(n_high1), mean KL: $(mean(KL_valid[.!low_x1_bool]))\n")

    write(report_file, "\n\nX2 below $(x2_cutoff), n_fit: $(n_low2_fit), n_valid: $(n_low2), mean KL: $(mean(KL_valid[low_x2_bool]))\n")
    write(report_file, "\nX2 above $(x2_cutoff), n_fit: $(n_high2_fit), n_valid: $(n_high2), mean KL: $(mean(KL_valid[.!low_x2_bool]))\n")

    write(report_file, "\n\nX1 below $(x1_cutoff) and X2 below $(x2_cutoff), n_fit: $(n_low1_low2_fit), n_valid: $(n_low1_low2), mean KL: $(mean(KL_valid[low1_low2_bool]))\n")
    write(report_file, "\nX1 below $(x1_cutoff) and X2 above $(x2_cutoff), n_fit: $(n_low1_high2_fit), n_valid: $(n_low1_high2), mean KL: $(mean(KL_valid[low1_high2_bool]))\n")
    write(report_file, "\nX1 above $(x1_cutoff) and X2 below $(x2_cutoff), n_fit: $(n_high1_low2_fit), n_valid: $(n_high1_low2), mean KL: $(mean(KL_valid[high1_low2_bool]))\n")
    write(report_file, "\nX1 above $(x1_cutoff) and X2 above $(x2_cutoff), n_fit: $(n_high1_high2_fit), n_valid: $(n_high1_high2), mean KL: $(mean(KL_valid[high1_high2_bool]))\n")

end

write(report_file, "\n\nFinal KL: $(KL_out)\n")
close(report_file)


text:
using AbstractTrees
import AbstractTrees: children, printnode

struct File
    path::String
end

children(f::File) = ()

struct Directory
    path::String
end

function children(d::Directory)
    contents = readdir(d.path)
    children = Vector{Union{Directory,File}}(undef,length(contents))
    for (i,c) in enumerate(contents)
        path = joinpath(d.path,c)
        children[i] = isdir(path) ? Directory(path) : File(path)
    end
    return children
end

printnode(io::IO, d::Directory) = print(io, basename(d.path))
printnode(io::IO, f::File) = print(io, basename(f.path))

dirpath = realpath(joinpath(dirname(pathof(AbstractTrees)),".."))
d = Directory(dirpath)
print_tree(d)


text:
const IN_PLACE = MPI_IN_PLACE

"""
    Barrier(comm::Comm)

Blocks until `comm` is synchronized.

If `comm` is an intracommunicator, then it blocks until all members of the group have called it.

If `comm` is an intercommunicator, then it blocks until all members of the other group have called it.

# External links
$(_doc_external("MPI_Barrier"))
"""
function Barrier(comm::Comm)
    # int MPI_Barrier(MPI_Comm comm)
    @mpichk ccall((:MPI_Barrier, libmpi), Cint, (MPI_Comm,), comm)
    return nothing
end

"""
    Bcast!(buf[, count=length(buf)], root::Integer, comm::Comm)

Broadcast the first `count` elements of the buffer `buf` from `root` to all processes.

# External links
$(_doc_external("MPI_Bcast"))
"""
function Bcast!(buffer, count::Integer,
                root::Integer, comm::Comm)

    # int MPI_Bcast(void* buffer, int count, MPI_Datatype datatype, int root,
    #               MPI_Comm comm)
    @mpichk ccall((:MPI_Bcast, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, MPI_Comm),
                  buffer, count, Datatype(eltype(buffer)), root, comm)
    buffer
end

function Bcast!(buffer::Union{Ref, AbstractArray}, root::Integer, comm::Comm)
    Bcast!(buffer, length(buffer), root, comm)
end

function bcast(obj, root::Integer, comm::Comm)
    isroot = Comm_rank(comm) == root
    count = Ref{Cint}()
    if isroot
        buf = MPI.serialize(obj)
        count[] = length(buf)
    end
    Bcast!(count, root, comm)
    if !isroot
        buf = Array{UInt8}(undef, count[])
    end
    Bcast!(buf, root, comm)
    if !isroot
        obj = MPI.deserialize(buf)
    end
    obj
end

"""
    Scatter!(sendbuf, recvbuf[, count=length(recvbuf)], root::Integer, comm::Comm)

Splits the buffer `sendbuf` in the `root` process into `Comm_size(comm)` chunks of length
`count`, and sends the `j`-th chunk to the process of rank `j` into the `recvbuf` buffer.

`sendbuf` on the root process should be a buffer of length `count*Comm_size(comm)`, and
on non-root processes it is ignored and can be `nothing`.

`recvbuf` can be `nothing` on the root process, in which case it is unmodified (this
corresponds the behaviour of `MPI_IN_PLACE` in `MPI_Scatter`). For example
```
if root == MPI.Comm_rank(comm)
    Scatter!(buf, nothing, count, root, comm)
else
    Scatter!(nothing, buf, count, root, comm)
end
```

`count` should be the same for all processes.

# See also
- [`Scatter`](@ref) to allocate the output buffer.
- [`Scatterv!`](@ref) if the number of elements varies between processes.

# External links
$(_doc_external("MPI_Scatter"))
"""
function Scatter!(sendbuf, recvbuf, count::Integer, root::Integer, comm::Comm)
    isroot = Comm_rank(comm) == root
    if isroot
        @assert sendbuf !== nothing
        @assert_minlength sendbuf count*Comm_size(comm)
        if recvbuf === nothing
            recvbuf = MPI_IN_PLACE
        end
    end
    @assert_minlength recvbuf count
    T = recvbuf isa SentinelPtr ? eltype(sendbuf) : eltype(recvbuf)
    # int MPI_Scatter(const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    #                 void* recvbuf, int recvcount, MPI_Datatype recvtype, int root,
    #                 MPI_Comm comm)
    @mpichk ccall((:MPI_Scatter, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, MPIPtr, Cint, MPI_Datatype, Cint, MPI_Comm),
                  sendbuf, count, Datatype(T), recvbuf, count, Datatype(T), root, comm)
    recvbuf
end

function Scatter!(sendbuf::AbstractArray, recvbuf::Union{Ref,AbstractArray}, root::Integer, comm::Comm)
    Scatter!(sendbuf, recvbuf, length(recvbuf), root, comm)
end

"""
    Scatter(sendbuf, count, root, comm)

Splits the buffer `sendbuf` in the `root` process into `Comm_size(comm)` chunks
and sends the `j`-th chunk to the process of rank `j`, allocating the output buffer.

# See also
- [`Scatter!`](@ref) for the mutating operation.
- [`Scatterv!`](@ref) if the number of elements varies between processes.

# External links
$(_doc_external("MPI_Scatter"))
"""
function Scatter(sendbuf, count::Integer, root::Integer, comm::Comm)
    Scatter!(sendbuf, similar(sendbuf, count), count, root, comm)
end

"""
    Scatterv!(sendbuf, recvbuf, counts, root, comm)

Splits the buffer `sendbuf` in the `root` process into `Comm_size(comm)` chunks
of length `counts[j]` and sends the j-th chunk to the process of rank j into the
`recvbuf` buffer, which must be of length at least `count`.

`recvbuf` can be `nothing` on the root process, in which case it is unmodified (this
corresponds the behaviour of `MPI_IN_PLACE` in `MPI_Scatterv`). For example
```
if root == MPI.Comm_rank(comm)
    Scatterv!(buf, nothing, counts, root, comm)
else
    Scatterv!(nothing, buf, counts, root, comm)
end
```

# See also
- [`Scatterv`](@ref) for the allocating operation
- [`Scatter!`](@ref)/[`Scatter`](@ref) if the counts are the same for all processes

# External links
$(_doc_external("MPI_Scatterv"))
"""
function Scatterv!(sendbuf, recvbuf, counts::Vector, root::Integer, comm::Comm)
    rank = Comm_rank(comm)
    isroot = rank == root
    if isroot
        @assert sendbuf !== nothing
        @assert_minlength sendbuf sum(counts)
    end
    if recvbuf === nothing
        recvbuf = MPI_IN_PLACE
    end
    @assert_minlength recvbuf counts[rank+1]
    T = recvbuf isa SentinelPtr ? eltype(sendbuf) : eltype(recvbuf)
    recvcnt = counts[Comm_rank(comm) + 1]
    disps = accumulate(+, counts) - counts
    # int MPI_Scatterv(const void* sendbuf, const int sendcounts[],
    #                  const int displs[], MPI_Datatype sendtype, void* recvbuf,
    #                  int recvcount, MPI_Datatype recvtype, int root, MPI_Comm comm)
    @mpichk ccall((:MPI_Scatterv, libmpi), Cint,
                  (MPIPtr, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, MPIPtr, Cint, MPI_Datatype, Cint, MPI_Comm),
                  sendbuf, counts, disps, Datatype(T), recvbuf, recvcnt, Datatype(T), root, comm)
    recvbuf
end

"""
    Scatterv(sendbuf, counts, root, comm)

Splits the buffer `sendbuf` in the `root` process into `Comm_size(comm)` chunks
of length `counts[j]` and sends the j-th chunk to the process of rank j, which
allocates the output buffer

# See also
- [`Scatterv!`](@ref) for the mutating operation
- [`Scatter!`](@ref)/[`Scatter`](@ref) if the counts are the same for all processes

# External links
$(_doc_external("MPI_Scatterv"))
"""
function Scatterv(sendbuf, counts::Vector, root::Integer, comm::Comm)
    recvbuf = similar(sendbuf, counts[Comm_rank(comm) + 1])
    Scatterv!(sendbuf, recvbuf, counts, root, comm)
end

"""
    Gather!(sendbuf, recvbuf[, count::Integer=length(sendbuf)], root::Integer, comm::Comm)

Each process sends the first `count` elements of the buffer `sendbuf` to the
`root` process. The `root` process stores elements in rank order in the buffer
buffer `recvbuf`.

`sendbuf` can be `nothing` on the root process, in which case the corresponding entries in
`recvbuf` are assumed to be already in place (this corresponds the behaviour of
`MPI_IN_PLACE` in `MPI_Gather`). For example
```
if root == MPI.Comm_rank(comm)
    Gather!(nothing, buf, count, root, comm)
else
    Gather!(buf, nothing, count, root, comm)
end
```

`recvbuf` on the root process should be a buffer of length `count*Comm_size(comm)`, and
on non-root processes it is ignored and can be `nothing`.

`count` should be the same for all processes.

# See also
- [`Gather`](@ref) for the allocating operation.
- [`Gatherv!`](@ref) if the number of elements varies between processes.
- [`Allgather!`](@ref) to send the result to all processes.

# External links
$(_doc_external("MPI_Gather"))
"""
function Gather!(sendbuf, recvbuf, count::Integer, root::Integer, comm::Comm)
    isroot = Comm_rank(comm) == root
    if isroot
        @assert recvbuf !== nothing
        @assert_minlength recvbuf count*Comm_size(comm)
        if sendbuf === nothing
            sendbuf = MPI_IN_PLACE
        end
    end
    @assert_minlength sendbuf count
    T = sendbuf isa SentinelPtr ? eltype(recvbuf) : eltype(sendbuf)
    # int MPI_Gather(const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    #                void* recvbuf, int recvcount, MPI_Datatype recvtype, int root,
    #                MPI_Comm comm)
    @mpichk ccall((:MPI_Gather, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, MPIPtr, Cint, MPI_Datatype, Cint, MPI_Comm),
                  sendbuf, count, Datatype(T), recvbuf, count, Datatype(T), root, comm)
    isroot ? recvbuf : nothing
end
function Gather!(sendbuf, recvbuf, root::Integer, comm::Comm)
    Gather!(sendbuf, recvbuf, length(sendbuf), root, comm)
end

"""
    Gather(sendbuf[, count=length(sendbuf)], root, comm)

Each process sends the first `count` elements of the buffer `sendbuf` to the `root`
process. The `root` allocates the output buffer and stores elements in rank order.

`count` should be the same for all processes.

# See also
- [`Gather!`](@ref) for the mutating operation.
- [`Gatherv!`](@ref)/[`Gatherv`](@ref) if the number of elements varies between processes.
- [`Allgather!`](@ref)/[`Allgather`](@ref) to send the result to all processes.

# External links
$(_doc_external("MPI_Gather"))
"""
function Gather(sendbuf, count::Integer, root::Integer, comm::Comm)
    Gather!(sendbuf, Comm_rank(comm) == root ? similar(sendbuf, Comm_size(comm) * count) : nothing, count, root, comm)
end
function Gather(sendbuf::AbstractArray, root::Integer, comm::Comm)
    Gather(sendbuf, length(sendbuf), root, comm)
end
function Gather(object::T, root::Integer, comm::Comm) where {T}
    Gather!(Ref(object), Comm_rank(comm) == root ? Vector{T}(undef, Comm_size(comm)) : nothing, 1, root, comm)
end

"""
    Allgather!(sendbuf, recvbuf[, count::Integer=length(sendbuf)], comm::Comm)
    Allgather!(sendrecvbuf, count::Integer, comm::Comm)

Each process sends the first `count` elements of `sendbuf` to the other processes, who
store the results in rank order into `recvbuf`.

`count` should be the same for all processes.

If only one buffer `sendrecvbuf` is provided, then each process send data is assumed to be
in the area where it would receive it's own contribution.

# See also
- [`Allgather`](@ref) for the allocating operation
- [`Allgatherv!`](@ref)/[`Allgatherv`](@ref) if the number of elements varies between processes.
- [`Gather!`](@ref) to send only to a single root process

# External links
$(_doc_external("MPI_Allgather"))
"""
function Allgather!(sendbuf, recvbuf, count::Integer, comm::Comm)
    @assert recvbuf !== nothing
    @assert_minlength recvbuf count*Comm_size(comm)
    @assert_minlength sendbuf count
    T = eltype(recvbuf)

    # int MPI_Allgather(const void* sendbuf, int sendcount,
    #                   MPI_Datatype sendtype, void* recvbuf, int recvcount,
    #                   MPI_Datatype recvtype, MPI_Comm comm)
    @mpichk ccall((:MPI_Allgather, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, MPIPtr, Cint, MPI_Datatype, MPI_Comm),
                  sendbuf, count, Datatype(T), recvbuf, count, Datatype(T), comm)
    recvbuf
end
function Allgather!(sendrecvbuf, count::Integer, comm::Comm)
    Allgather!(MPI.IN_PLACE, sendrecvbuf, count, comm)
end

"""
    Allgather(sendbuf[, count=length(sendbuf)], comm)

Each process sends the first `count` elements of `sendbuf` to the other processes, who
store the results in rank order allocating the output buffer.

`count` should be the same for all processes.

# See also
- [`Allgather!`](@ref) for the mutating operation
- [`Allgatherv!`](@ref)/[`Allgatherv`](@ref) if the number of elements varies between processes.
- [`Gather!`](@ref) to send only to a single root process

# External links
$(_doc_external("MPI_Allgather"))
"""
function Allgather(sendbuf, count::Integer, comm::Comm)
    Allgather!(sendbuf, similar(sendbuf, Comm_size(comm) * count), count, comm)
end
function Allgather(sendbuf::AbstractArray, comm::Comm)
    Allgather(sendbuf, length(sendbuf), comm)
end
function Allgather(object::T, comm::Comm) where T
    Allgather!(Ref(object), Vector{T}(undef, Comm_size(comm)), 1, comm)
end

"""
    Gatherv!(sendbuf, recvbuf, counts, root, comm)

Each process sends the first `counts[rank]` elements of the buffer `sendbuf` to
the `root` process. The `root` stores elements in rank order in the buffer
`recvbuf`.

`sendbuf` can be `nothing` on the root process, in which case the corresponding entries in
`recvbuf` are assumed to be already in place (this corresponds the behaviour of
`MPI_IN_PLACE` in `MPI_Gatherv`). For example
```
if root == MPI.Comm_rank(comm)
    Gatherv!(nothing, buf, counts, root, comm)
else
    Gatherv!(buf, nothing, counts, root, comm)
end
```

# See also
- [`Gatherv`](@ref) for the allocating operation
- [`Gather!`](@ref)
- [`Allgatherv!`](@ref)/[`Allgatherv`](@ref) to send the result to all processes

# External links
$(_doc_external("MPI_Gatherv"))
"""
function Gatherv!(sendbuf, recvbuf, counts::Vector{Cint}, root::Integer, comm::Comm)
    isroot = Comm_rank(comm) == root
    displs = accumulate(+, counts) - counts
    sendcnt = counts[Comm_rank(comm) + 1]
    if isroot
        @assert recvbuf !== nothing
        @assert_minlength recvbuf sum(counts)
        if sendbuf === nothing
            sendbuf = MPI_IN_PLACE
        end
    end
    @assert_minlength sendbuf sendcnt
    T = sendbuf isa SentinelPtr ? eltype(recvbuf) : eltype(sendbuf)

    # int MPI_Gatherv(const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    #                 void* recvbuf, const int recvcounts[], const int displs[],
    #                 MPI_Datatype recvtype, int root, MPI_Comm comm)
    @mpichk ccall((:MPI_Gatherv, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, MPIPtr, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, Cint, MPI_Comm),
                  sendbuf, sendcnt, Datatype(T), recvbuf, counts, displs, Datatype(T), root, comm)
    isroot ? recvbuf : nothing
end

"""
    Gatherv(sendbuf, counts, root, comm)

Each process sends the first `counts[rank]` elements of the buffer `sendbuf` to
the `root` process. The `root` allocates the output buffer and stores elements
in rank order.

# See also
- [`Gatherv!`](@ref) for the mutating operation
- [`Gather!`](@ref)/[`Gather`](@ref)
- [`Allgatherv!`](@ref)/[`Allgatherv`](@ref) to send the result to all processes

# External links
$(_doc_external("MPI_Gatherv"))
"""
function Gatherv(sendbuf, counts::Vector{Cint}, root::Integer, comm::Comm)
    Gatherv!(sendbuf, Comm_rank(comm) == root ? similar(sendbuf, sum(counts)) : nothing, counts, root, comm)
end

"""
    Allgatherv!(sendbuf, recvbuf, counts, comm)
    Allgatherv!(sendrecvbuf, counts, comm)

Each process sends the first `counts[rank]` elements of the buffer `sendbuf` to
all other process. Each process stores the received data in rank order in the
buffer `recvbuf`.

If only one buffer `sendrecvbuf` is provided, then for each process, the data to be sent
is taken from the interval of `recvbuf` where it would store it's own data.

# See also
- [`Allgatherv`](@ref) for the allocating operation
- [`Gatherv!`](@ref)/[`Gatherv`](@ref) to send the result to a single process

# External links
$(_doc_external("MPI_Allgatherv"))
"""
function Allgatherv!(sendbuf, recvbuf, counts::Vector{Cint}, comm::Comm)
    displs = accumulate(+, counts) - counts
    sendcnt = counts[Comm_rank(comm) + 1]
    @assert recvbuf !== nothing
    @assert_minlength recvbuf sum(counts)
    @assert_minlength sendbuf sendcnt
    T = eltype(recvbuf)

    # int MPI_Allgatherv(const void* sendbuf, int sendcount,
    #                    MPI_Datatype sendtype, void* recvbuf, const int recvcounts[],
    #                    const int displs[], MPI_Datatype recvtype, MPI_Comm comm)
    @mpichk ccall((:MPI_Allgatherv, libmpi), Cint,
              (MPIPtr, Cint, MPI_Datatype, MPIPtr, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, MPI_Comm),
              sendbuf, sendcnt, Datatype(T), recvbuf, counts, displs, Datatype(T), comm)
    recvbuf
end
function Allgatherv!(sendrecvbuf, counts::Vector{Cint}, comm::Comm)
    Allgatherv!(MPI_IN_PLACE, sendrecvbuf, counts, comm)
end

"""
    Allgatherv(sendbuf, counts, comm)

Each process sends the first `counts[rank]` elements of the buffer `sendbuf` to
all other process. Each process allocates an output buffer and stores the
received data in rank order.

# See also
- [`Allgatherv!`](@ref) for the mutating operation.
- [`Gatherv!`](@ref)/[`Gatherv`](@ref) to send the result to a single process.

# External links
$(_doc_external("MPI_Allgatherv"))
"""
function Allgatherv(sendbuf, counts::Vector{Cint}, comm::Comm)
    recvbuf = similar(sendbuf, sum(counts))
    Allgatherv!(sendbuf, recvbuf, counts, comm)
end

"""
    Alltoall!(sendbuf, recvbuf, count::Integer, comm::Comm)
    Alltoall!(sendrecvbuf, count::Integer, comm::Comm)

Every process divides the buffer `sendbuf` into `Comm_size(comm)` chunks of
length `count`, sending the `j`-th chunk to the `j`-th process.
Every process stores the data received from the `j`-th process in the `j`-th
chunk of the buffer `recvbuf`.

```
rank    send buf                        recv buf
----    --------                        --------
 0      a,b,c,d,e,f       Alltoall      a,b,A,B,α,β
 1      A,B,C,D,E,F  ---------------->  c,d,C,D,γ,ψ
 2      α,β,γ,ψ,η,ν                     e,f,E,F,η,ν
```

If only one buffer `sendrecvbuf` then data is overwritten.

# See also
- [`Alltoall`](@ref) for the allocating operation

# External links
$(_doc_external("MPI_Alltoall"))
"""
function Alltoall!(sendbuf, recvbuf, count::Integer, comm::Comm)
    buflength = count * Comm_size(comm)
    @assert_minlength recvbuf buflength
    @assert_minlength sendbuf buflength
    sendbuf isa SentinelPtr || @assert eltype(sendbuf) == eltype(recvbuf)
    T = eltype(recvbuf)
    # int MPI_Alltoall(const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    #                  void* recvbuf, int recvcount, MPI_Datatype recvtype,
    #                  MPI_Comm comm)
    @mpichk ccall((:MPI_Alltoall, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, MPIPtr, Cint, MPI_Datatype, MPI_Comm),
                  sendbuf, count, Datatype(T), recvbuf, count, Datatype(T), comm)
    recvbuf
end
function Alltoall!(sendrecvbuf, count::Integer, comm::Comm)
    Alltoall!(MPI_IN_PLACE, recvbuf, count, comm)
end

"""
    Alltoall(sendbuf, count::Integer, comm::Comm)

Every process divides the buffer `sendbuf` into `Comm_size(comm)` chunks of
length `count`, sending the `j`-th chunk to the `j`-th process.
Every process allocates the output buffer and stores the data received from the
`j`-th process in the `j`-th chunk.

```
rank    send buf                        recv buf
----    --------                        --------
 0      a,b,c,d,e,f       Alltoall      a,b,A,B,α,β
 1      A,B,C,D,E,F  ---------------->  c,d,C,D,γ,ψ
 2      α,β,γ,ψ,η,ν                     e,f,E,F,η,ν
```

# See also
- [`Alltoall!`](@ref) for the mutating operation

# External links
$(_doc_external("MPI_Alltoall"))
"""
function Alltoall(sendbuf, count::Integer,  comm::Comm)
    recvbuf = similar(sendbuf, Comm_size(comm)*count)
    Alltoall!(sendbuf, recvbuf, count, comm)
end

"""
    Alltoallv!(sendbuf, recvbuf, scounts::Vector, rcounts::Vector, comm::Comm)

Similar to [`Alltoall!`](@ref), except with different size chunks per process.

# See also
- [`Alltoallv`](@ref) for the allocating operation

# External links
$(_doc_external("MPI_Alltoallv"))
"""
function Alltoallv!(sendbuf, recvbuf, scounts::Vector{Cint}, rcounts::Vector{Cint}, comm::Comm)
    @assert_minlength sendbuf sum(scounts)
    @assert_minlength recvbuf sum(rcounts)
    @assert eltype(sendbuf) == eltype(recvbuf)
    T = eltype(sendbuf)

    sdispls = accumulate(+, scounts) - scounts
    rdispls = accumulate(+, rcounts) - rcounts
    # int MPI_Alltoallv(const void* sendbuf, const int sendcounts[],
    #                   const int sdispls[], MPI_Datatype sendtype, void* recvbuf,
    #                   const int recvcounts[], const int rdispls[],
    #                   MPI_Datatype recvtype, MPI_Comm comm)
    @mpichk ccall((:MPI_Alltoallv, libmpi), Cint,
                  (MPIPtr, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, MPIPtr, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, MPI_Comm),
                  sendbuf, scounts, sdispls, Datatype(T), recvbuf, rcounts, rdispls, Datatype(T), comm)
    recvbuf
end

"""
    Alltoallv(sendbuf, recvbuf, scounts::Vector, rcounts::Vector, comm::Comm)

Similar to [`Alltoall`](@ref), except with different size chunks per process.

# See also
- [`Alltoallv!`](@ref) for the mutating operation

# External links
$(_doc_external("MPI_Alltoallv"))
"""
function Alltoallv(sendbuf, scounts::Vector{Cint},
                   rcounts::Vector{Cint}, comm::Comm)
    recvbuf = similar(sendbuf, sum(rcounts))
    Alltoallv!(sendbuf, recvbuf, scounts, rcounts, comm)
end

### Reduce/Scan

## Reduce

# mutating
"""
    Reduce!(sendbuf, recvbuf[, count::Integer=length(sendbuf)], op, root::Integer, comm::Comm)
    Reduce!(sendrecvbuf, op, root::Integer, comm::Comm)

Performs elementwise reduction using the operator `op` on the first `count` elements of
the buffer `sendbuf` and stores the result in `recvbuf` on the process of rank `root`.

On non-root processes `recvbuf` is ignored, and can be `nothing`.

To perform the reduction in place, provide a single buffer `sendrecvbuf`.

# See also
- [`Reduce`](@ref) to handle allocation of the output buffer.
- [`Allreduce!`](@ref)/[`Allreduce`](@ref) to send reduction to all ranks.
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Reduce"))
"""
function Reduce!(sendbuf, recvbuf, count::Integer, op::Union{Op,MPI_Op}, root::Integer, comm::Comm)
    isroot = Comm_rank(comm) == root
    @assert_minlength sendbuf count
    if isroot
        @assert recvbuf !== nothing
        @assert_minlength recvbuf count
    end
    T = sendbuf isa SentinelPtr ? eltype(recvbuf) : eltype(sendbuf)
    # int MPI_Reduce(const void* sendbuf, void* recvbuf, int count,
    #                MPI_Datatype datatype, MPI_Op op, int root, MPI_Comm comm)
    @mpichk ccall((:MPI_Reduce, libmpi), Cint,
                  (MPIPtr, MPIPtr, Cint, MPI_Datatype, MPI_Op, Cint, MPI_Comm),
                  sendbuf, recvbuf, count, Datatype(T), op, root, comm)
    recvbuf
end

# Convert user-provided functions to MPI.Op
function Reduce!(sendbuf, recvbuf, count::Integer, opfunc, root::Integer, comm::Comm)
    T = sendbuf isa SentinelPtr ? eltype(recvbuf) : eltype(sendbuf)
    Reduce!(sendbuf, recvbuf, count, Op(opfunc, T), root, comm)
end
function Reduce!(sendbuf::AbstractArray, recvbuf::Union{Nothing,AbstractArray},
                 op, root::Integer, comm::Comm)
    Reduce!(sendbuf, recvbuf, length(sendbuf), op, root, comm)
end

# inplace
function Reduce!(buf, op, root::Integer, comm::Comm)
    if Comm_rank(comm) == root
        Reduce!(MPI_IN_PLACE, buf, length(buf), op, root, comm)
    else
        Reduce!(buf, nothing, length(buf), op, root, comm)
    end
end

# allocating
"""
    recvbuf = Reduce(sendbuf, op, root::Integer, comm::Comm)

Performs elementwise reduction using the operator `op` on the buffer `sendbuf`, returning
the result `recvbuf` on the process of rank `root`, and `nothing` on non-root processes.

`sendbuf` can also be a scalar, in which case `recvbuf` will be a value of the same type.

# See also
- [`Reduce!`](@ref) for mutating and in-place operations
- [`Allreduce!`](@ref)/[`Allreduce`](@ref) to send reduction to all ranks.
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Reduce"))
"""
function Reduce(sendbuf::AbstractArray, op, root::Integer, comm::Comm)
    Reduce!(sendbuf, Comm_rank(comm) == root ? similar(sendbuf) : nothing, length(sendbuf), op, root, comm)
end
function Reduce(object::T, op, root::Integer, comm::Comm) where {T}
    if Comm_rank(comm) == root
        Reduce!(Ref(object), Ref{T}(), 1, op, root, comm)[]
    else
        Reduce!(Ref(object), nothing, 1, op, root, comm)
    end
end

## Allreduce

# mutating
"""
    Allreduce!(sendbuf, recvbuf[, count=length(sendbuf)], op, comm)
    Allreduce!(sendrecvbuf, op, comm)

Performs elementwise reduction using the operator `op` on the first `count` elements of
the buffer `sendbuf`, storing the result in the `recvbuf` of all processes in the group.

`Allreduce!` is equivalent to a [`Reduce!`](@ref) operation followed by
a [`Bcast!`](@ref), but can lead to better performance.

If only one `sendrecvbuf` buffer is provided, then the operation is performed in-place.

# See also
- [`Allreduce`](@ref), to handle allocation of the output buffer.
- [`Reduce!`](@ref)/[`Reduce`](@ref) to send reduction to a single rank.
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Allreduce"))
"""
function Allreduce!(sendbuf, recvbuf, count::Integer, op::Union{Op,MPI_Op}, comm::Comm)
    @assert_minlength sendbuf count
    @assert_minlength recvbuf count
    sendbuf isa SentinelPtr || @assert eltype(sendbuf) == eltype(recvbuf)
    T = eltype(recvbuf)
    # int MPI_Allreduce(const void* sendbuf, void* recvbuf, int count,
    #                   MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
    @mpichk ccall((:MPI_Allreduce, libmpi), Cint,
                  (MPIPtr, MPIPtr, Cint, MPI_Datatype, MPI_Op, MPI_Comm),
                  sendbuf, recvbuf, count, Datatype(T), op, comm)
    recvbuf
end
function Allreduce!(sendbuf, recvbuf, count::Integer, opfunc, comm::Comm)
    Allreduce!(sendbuf, recvbuf, count, Op(opfunc, eltype(recvbuf)), comm)
end

function Allreduce!(sendbuf::AbstractArray, recvbuf::AbstractArray, op, comm::Comm)
    Allreduce!(sendbuf, recvbuf, length(recvbuf), op, comm)
end

# inplace
function Allreduce!(buf, op, comm::Comm)
    Allreduce!(MPI_IN_PLACE, buf, length(buf), op, comm)
end

# allocating
"""
    recvbuf = Allreduce(sendbuf, op, comm)

Performs elementwise reduction using the operator `op` on the buffer `sendbuf`, returning
the result in the `recvbuf` of all processes in the group.

`sendbuf` can also be a scalar, in which case `recvbuf` will be a value of the same type.

# See also
- [`Allreduce!`](@ref) for mutating or in-place operations.
- [`Reduce!`](@ref)/[`Reduce`](@ref) to send reduction to a single rank.
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Allreduce"))
"""
function Allreduce(sendbuf::AbstractArray, op, comm::Comm)
    Allreduce!(sendbuf, similar(sendbuf), length(sendbuf), op, comm)
end
function Allreduce(obj::T, op, comm::Comm) where {T}
    Allreduce!(Ref(obj), Ref{T}(), 1, op, comm)[]
end

## Scan

# mutating
"""
    Scan!(sendbuf, recvbuf[, count::Integer], op, comm::Comm)
    Scan!(buf[, count::Integer], op, comm::Comm)

Inclusive prefix reduction (analagous to `accumulate` in Julia): `recvbuf` on rank `i`
will contain the the result of reducing `sendbuf` by `op` from ranks `0:i`.

If only a single buffer is provided, then operations will be performed in-place in `buf`.

# See also
- [`Scan`](@ref) to handle allocation of the output buffer
- [`Exscan!`](@ref)/[`Exscan`](@ref) for exclusive scan
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Scan"))
"""
function Scan!(sendbuf, recvbuf, count::Integer,
               op::Union{Op,MPI_Op}, comm::Comm)
    T = eltype(recvbuf)
    # int MPI_Scan(const void* sendbuf, void* recvbuf, int count,
    #              MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
    @mpichk ccall((:MPI_Scan, libmpi), Cint,
                  (MPIPtr, MPIPtr, Cint, MPI_Datatype, MPI_Op, MPI_Comm),
                  sendbuf, recvbuf, count, Datatype(T), op, comm)
    recvbuf
end
function Scan!(sendbuf, recvbuf, count::Integer, opfunc, comm::Comm)
    Scan!(sendbuf, recvbuf, count, Op(opfunc, eltype(recvbuf)), comm)
end
function Scan!(sendbuf::AbstractArray, recvbuf, op, comm::Comm)
    Scan!(sendbuf, recvbuf, length(sendbuf), op, comm)
end

# inplace
function Scan!(buf, count::Integer, opfunc, comm::Comm)
    Scan!(MPI_IN_PLACE, buf, count, Op(opfunc, eltype(sendbuf)), comm)
end
function Scan!(buf, opfunc, comm::Comm)
    Scan!(buf, length(buf), Op(opfunc, eltype(sendbuf)), comm)
end

# allocating
"""
    recvbuf = Scan(sendbuf, op, comm::Comm)

Inclusive prefix reduction (analagous to `accumulate` in Julia): `recvbuf` on rank `i`
will contain the the result of reducing `sendbuf` by `op` from ranks `0:i`.

`sendbuf` can also be a scalar, in which case `recvbuf` will also be a scalar of the same
type.

# See also
- [`Scan!`](@ref) for mutating or in-place operations
- [`Exscan!`](@ref)/[`Exscan`](@ref) for exclusive scan
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Scan"))
"""
function Scan(sendbuf::AbstractArray, op, comm::Comm)
    Scan!(sendbuf, similar(sendbuf), op, comm)
end
function Scan(object::T, op, comm::Comm) where {T}
    Scan!(Ref(object), Ref{T}(), 1, op, comm)[]
end

## Exscan

# mutating

"""
    Exscan!(sendbuf, recvbuf[, count::Integer], op, comm::Comm)
    Exscan!(buf[, count::Integer], op, comm::Comm)

Exclusive prefix reduction (analagous to `accumulate` in Julia): `recvbuf` on rank `i`
will contain the the result of reducing `sendbuf` by `op` from ranks `0:i-1`. The
`recvbuf` on rank `0` is ignored, and the `recvbuf` on rank `1` will contain the contents
of `sendbuf` on rank `0`.

If only a single `buf` is provided, then operations are performed in-place, and `buf` on
rank 0 will remain unchanged.

# See also
- [`Exscan`](@ref) to handle allocation of the output buffer
- [`Scan!`](@ref)/[`Scan`](@ref) for inclusive scan
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Exscan"))
"""
function Exscan!(sendbuf, recvbuf, count::Integer,
                op::Union{Op,MPI_Op}, comm::Comm)
    T = eltype(recvbuf)
    # int MPI_Exscan(const void* sendbuf, void* recvbuf, int count,
    #                MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
    @mpichk ccall((:MPI_Exscan, libmpi), Cint,
          (MPIPtr, MPIPtr, Cint, MPI_Datatype, MPI_Op, MPI_Comm),
          sendbuf, recvbuf, count, Datatype(T), op, comm)
    recvbuf
end
function Exscan!(sendbuf, recvbuf, count::Integer, opfunc, comm::Comm)
    Exscan!(sendbuf, recvbuf, count, Op(opfunc, eltype(recvbuf)), comm)
end
function Exscan!(sendbuf::AbstractArray, recvbuf::AbstractArray, op, comm::Comm)
    Exscan!(sendbuf, recvbuf, length(sendbuf), op, comm)
end

# inplace
function Exscan!(buf, count::Integer, opfunc, comm::Comm)
    Exscan!(MPI_IN_PLACE, buf, count, Op(opfunc, eltype(buf)), comm)
end
function Exscan!(buf, op, comm::Comm)
    Exscan!(buf, length(buf), op, comm)
end

# allocating
"""
    recvbuf = Exscan(sendbuf, op, comm::Comm)

Exclusive prefix reduction (analagous to `accumulate` in Julia): `recvbuf` on rank `i`
will contain the the result of reducing `sendbuf` by `op` from ranks `0:i-1`. The
`recvbuf` on rank `0` is undefined, and the `recvbuf` on rank `1` will contain the contents
of `sendbuf` on rank `0`.

# See also
- [`Exscan!`](@ref) for mutating and in-place operations
- [`Scan!`](@ref)/[`Scan`](@ref) for inclusive scan
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Scan"))
"""
function Exscan(sendbuf::AbstractArray, op, comm::Comm)
    Exscan!(sendbuf, similar(sendbuf), op, comm)
end
function Exscan(object::T, op, comm::Comm) where {T}
    Exscan!(Ref(object), Ref{T}(), 1, op, comm)[]
end


text:
# Print out functions defined in a .jl file
# Usage: julia functions.jl source.jl

function functions(ex::Expr,s::Set)
    if ((ex.head == :function || ex.head == :(=)) && isa(ex.args[1],Expr) && ex.args[1].head == :call)
        fn = ex.args[1].args[1]
        isa(fn, Expr) && fn.head == :curly && (fn = fn.args[1])
        if !in(fn,s)
            push!(s,fn)
            if !isa(fn,Symbol)
                println("# Not a symbol: $fn")
            elseif !isdefined(Main,fn)
                println("# Not exported: $fn")
            else
                println("# $fn")
            end
        end
    else
        for a in ex.args
            if isa(a, Expr)
                functions(a,s)
            end
        end
    end
end

ex = Meta.parse("module Foo\n"*read(ARGS[1],String)*"\nend\n")
functions(ex,Set())


text:
# we will want to overload these operators and functions for our custom types
import Base: ==, ≈, *, +, -, isless, length, adjoint, print, zero, one
export comm

==(A::scal,B::scal) = A.v == B.v
==(A::OpProd,B::OpProd) = A.A == B.A && A.B == B.B
==(A::OpSum, B::OpSum)  = A.A == B.A && A.B == B.B
==(A::OpSumAnalytic, B::OpSumAnalytic) = A.A == B.A && A.ind == B.ind
==(A::ExpVal, B::ExpVal) = A.A == B.A
==(A::Corr, B::Corr) = A.A == B.A
==(A::param, B::param) = (A.name,A.inds,A.state) == (B.name,B.inds,B.state)

≈(A::Operator,B::Operator) = A == B
≈(A::scal,B::scal) = A.v ≈ B.v
≈(A::OpProd,B::OpProd) = A.A ≈ B.A && A.B ≈ B.B
≈(A::OpSum, B::OpSum)  = A.A ≈ B.A && A.B ≈ B.B
≈(A::OpSumAnalytic, B::OpSumAnalytic) = A.A ≈ B.A && A.ind == B.ind
≈(A::ExpVal, B::ExpVal) = A.A ≈ B.A
≈(A::Corr, B::Corr) = A.A ≈ B.A

# we sometimes want to sort integers and symbols together, with integers coming first
sortsentinel(x::Integer) = (1,x)
sortsentinel(x::Symbol) = (2,x)

OpOrder = (scal,δ,param,ExpVal,Corr,BosonCreate,BosonDestroy,FermionCreate,FermionDestroy,σplus,σminus,σ,OpProd,OpSumAnalytic,OpSum)
for (ii,op1) in enumerate(OpOrder)
    for op2 in OpOrder[ii+1:end]
        @eval isless(::$op1,::$op2) = true
        @eval isless(::$op2,::$op1) = false
    end
end
isless(A::scal,B::scal) = (real(A.v),imag(A.v)) < (real(B.v),imag(B.v))

for op in (BosonCreate,BosonDestroy,FermionCreate,FermionDestroy)
    @eval isless(A::$op,B::$op) = (A.name,sortsentinel.(A.inds)) < (B.name,sortsentinel.(B.inds))
end
# use conjugation state as last order parameter
isless(A::param,B::param) = (A.name,sortsentinel.(A.inds),A.state) < (B.name,sortsentinel.(B.inds),B.state)
isless(A::σ,B::σ) = (sortsentinel.(A.inds),A.a) < (sortsentinel.(B.inds),B.a)
for op in (δ,σminus,σplus)
    @eval isless(A::$op,B::$op) = sortsentinel.(A.inds) < sortsentinel.(B.inds)
end
for op in (ExpVal,Corr,OpSumAnalytic)
    @eval isless(A::$op,B::$op) = A.A < B.A
end

# return -1 if A<B, 0 if A==B, 1 if B<A
function iterlesseq(A,B)::Int
    resa = iterate(A)
    resb = iterate(B)
    while resa !== nothing && resb !== nothing
        vala, statea = resa
        valb, stateb = resb
        vala == valb || return isless(vala,valb) ? -1 : 1
        resa = iterate(A,statea)
        resb = iterate(B,stateb)
    end
    resa === nothing && resb === nothing && return 0
    # the longer iterator is larger
    resa === nothing ? 1 : -1
end

function isless(A::OpProd,B::OpProd)
    # only evaluate each part that we need for each step of the comparison to avoid unnecessary work
    # order operator products first by number of operators (also within expectation values)
    lA,lB = length(A), length(B)
    lA == lB || return lA < lB
    # then by operators
    c1 = iterlesseq(opiter(A,true),opiter(B,true))
    c1 == 0 || return c1 < 0
    # then by expectation values
    c2 = iterlesseq(expiter(A,true),expiter(B,true))
    c2 == 0 || return c2 < 0
    # then by reversed prefactors (to order params, not scalars)
    rpA = reverse(preftuple(A))
    rpB = reverse(preftuple(B))
    rpA == rpB || return rpA < rpB
    # finally without swallowing opsumanalytic to make sure we always have definite ordering
    iterlesseq(proditer(A,false),proditer(B,false)) < 0
end

# prefactors do not count for length calculation
length(::Union{scal,δ,param}) = 0
length(::BaseOperator) = 1
length(A::Union{ExpVal,Corr,OpSumAnalytic}) = length(A.A)
length(A::OpProd) = length(A.A) + length(A.B)

adjoint(A::scal) = scal(conj(A.v))
adjoint(A::param) = A.state=='r' ? A : param(A.name,A.state=='n' ? 'c' : 'n',A.inds)
adjoint(A::δ) = A
adjoint(A::BosonDestroy) = BosonCreate(A.name,A.inds)
adjoint(A::BosonCreate) = BosonDestroy(A.name,A.inds)
adjoint(A::FermionDestroy) = FermionCreate(A.name,A.inds)
adjoint(A::FermionCreate) = FermionDestroy(A.name,A.inds)
adjoint(A::σminus) = σplus(A.inds)
adjoint(A::σplus) = σminus(A.inds)
adjoint(A::σ) = A
adjoint(A::OpSum) = A.A' + A.B'
adjoint(A::OpProd) = A.B' * A.A'
adjoint(A::ExpVal) = ExpVal(A.A')
adjoint(A::Corr) = Corr(A.A')
adjoint(A::OpSumAnalytic) = OpSumAnalytic(A.ind,A.A')

zero(::Type{<:Operator}) = scal(0)
zero(::Operator) = scal(0)
one(::Type{<:Operator}) = scal(1)
one(::Operator) = scal(1)

*(A::Operator) = A # needed to make prod((A,)) work
*(A::OpProd,B::Operator) = A.A*(A.B*B)
*(A::OpSum,B::Operator)      = A.A*B + A.B*B
*(A::OpSum,B::OpSum)         = A.A*B + A.B*B # resolve ambiguity
*(A::OpSum,B::OpSumAnalytic) = A.A*B + A.B*B # resolve ambiguity
*(A::Operator,     B::OpSum) = A*B.A + A*B.B
*(A::OpProd,       B::OpSum) = A*B.A + A*B.B # resolve ambiguity
*(A::OpSumAnalytic,B::OpSum) = A*B.A + A*B.B # resolve ambiguity

# allow addition, substraction, and multiplication with a number x by promoting it to scal(x) operator
*(x::Number,A::Operator) = scal(x)*A
*(A::Operator,x::Number) = scal(x)*A
+(x::Number,A::Operator) = scal(x)+A
+(A::Operator,x::Number) = scal(x)+A
-(x::Number,A::Operator) = scal(x)-A
-(A::Operator,x::Number) = A-scal(x)

function *(A::Operator,B::Operator)::Operator
    if A isa scal && A.v==0
        A
    elseif A isa scal && A.v==1
        B
    elseif A isa scal && B isa scal
        scal(A.v*B.v)
    elseif A isa scal && B isa OpProd && B.A isa scal
        scal(A.v*B.A.v)*B.B
    elseif B isa OpProd && (A*B.A) != OpProd(A,B.A)
        # if A*B.A is not ordered as we want, evaluate A*B.A first
        (A*B.A)*B.B
    elseif A isa δ && A.inds[2] in indextuple(B)
        # if second index of δ_iA,iB shows up on RHS, replace by iA
        # (relies on indices in δ being sorted)
        A * replace_index(B,A.inds[2],A.inds[1])
    elseif A isa σ && B isa σ && A.inds == B.inds
        # σa σb = δab + i ϵabc σc = δab + 1/2 [σa,σb]
        A.a==B.a ? scal(1) : scal(1//2)*comm(A,B)
    elseif A isa Union{σplus,σminus,FermionDestroy,FermionCreate} && A == B
        scal(0)
    elseif A>B
        B*A + comm(A,B)
    else
        OpProd(A,B)
    end
end

# these are the only overrides for minus, everything else falls back to +
-(A::Operator) = scal(-1)*A
-(A::Operator,B::Operator) = A + scal(-1)*B

separate_prefac(A::Operator) = (1,A)
separate_prefac(A::OpProd) = A.A isa scal ? (A.A.v,A.B) : (1,A)

+(A::OpSum,B::Operator) = A.A + (A.B + B)
function +(A::Operator,B::Operator)::Operator
    A isa scal && A.v == 0 && return B
    A isa scal && B isa scal && return scal(A.v+B.v)
    A isa scal && B isa OpSum && B.A isa scal && return scal(A.v+B.A.v) + B.B

    # ignore scalars in ordering
    sA, oA = separate_prefac(A)
    if B isa OpSum
        sB, oB = separate_prefac(B.A)
        if oA > oB
            return B.A + (A + B.B)
        elseif oA == oB
            return scal(sA+sB)*oA + B.B
        end
    else
        sB, oB = separate_prefac(B)
        if oA > oB
            return B + A
        elseif oA == oB
            return scal(sA+sB)*oA
        end
    end
    oA < oB || error("invariant oA < oB violated! oA: $oA, oB: $oB, oA<oB: $(oA<oB), oA==oB: $(oA==oB), oA>oB: $(oA>oB)")

    return OpSum(A,B)
end

comm(A::Scalar,B::Operator) = scal(0)
comm(A::Scalar,B::OpProd)   = scal(0)
comm(A::Scalar,B::OpSum)    = scal(0)
comm(A::Scalar,B::OpSumAnalytic) = scal(0)
comm(A::Operator,B::Scalar) = scal(0)
comm(A::OpSum,   B::Scalar) = scal(0)
comm(A::OpProd,  B::Scalar) = scal(0)
comm(A::OpSumAnalytic,B::Scalar) = scal(0)
comm(A::Scalar,B::Scalar)   = scal(0)
comm(A::OpProd,B::Operator) = comm(A.A,B)*A.B + A.A*comm(A.B,B)
comm(A::OpProd,B::OpSumAnalytic) = comm(A.A,B)*A.B + A.A*comm(A.B,B)
comm(A::OpProd,B::OpProd)   = comm(A.A,B)*A.B + A.A*comm(A.B,B)
comm(A::Operator,B::OpProd) = comm(A,B.A)*B.B + B.A*comm(A,B.B)
comm(A::OpSumAnalytic, B::OpProd) = comm(A,B.A)*B.B + B.A*comm(A,B.B)
comm(A::OpSum,B::Operator)  = comm(A.A,B) + comm(A.B,B)
comm(A::OpSum,B::OpSum)     = comm(A.A,B) + comm(A.B,B)
comm(A::OpSum,B::OpProd)    = comm(A.A,B) + comm(A.B,B)
comm(A::Operator,B::OpSum)  = comm(A,B.A) + comm(A,B.B)
comm(A::OpProd,  B::OpSum)  = comm(A,B.A) + comm(A,B.B)
comm(A::OpSumAnalytic, B::OpSum) = comm(A,B.A) + comm(A,B.B)
# different types of operators commute
commgroups = (Union{BosonDestroy,BosonCreate},Union{FermionDestroy,FermionCreate},Union{σ,σminus,σplus})
for (ii,op) in enumerate(commgroups)
    for op2 in commgroups[ii+1:end]
        @eval comm(A::$op,B::$op2) = scal(0)
        @eval comm(A::$op2,B::$op) = scal(0)
    end
end
# these operators commute with themselves
for op in (BosonDestroy,BosonCreate,σplus,σminus)
    @eval comm(A::$op,B::$op) = scal(0)
end
comm(A::BosonDestroy,B::BosonCreate) = A.name == B.name ? δ(A.inds,B.inds) : scal(0)
comm(A::BosonCreate,B::BosonDestroy) = -comm(B,A)
comm(A::σplus,B::σminus) = δ(A.inds,B.inds)*σz(A.inds)
comm(A::σminus,B::σplus) = -comm(B,A)

# {f_i, fdag_j} = f_i fdag_j + fdag_j f_i = δ_{i,j}
# f_i fdag_j = δ_{i,j} - fdag_j f_i
# [f_i, fdag_j] = f_i fdag_j - fdag_j f_i = δ_{i,j} - fdag_j f_i - fdag_j f_i = δ_{i,j} - 2 fdag_j f_i
# if they correspond to operators on different systems (i.e., with different names), also Fermion operators commute
comm(A::FermionDestroy,B::FermionCreate) = A.name == B.name ? δ(A.inds,B.inds) - 2*B*A : scal(0)
comm(A::FermionCreate,B::FermionDestroy) = -comm(B,A)
# {f_i,f_j} = f_i f_j + f_j f_i = 0
# we assume that if we need the commutator, it's because of exchanging order
# [f_i,f_j] = f_i f_j - f_j f_i = -f_j f_i - f_j f_i
comm(A::FermionDestroy,B::FermionDestroy) = A.name == B.name ? -2*B*A : scal(0)
comm(A::FermionCreate,B::FermionCreate) = A.name == B.name ? -2*B*A : scal(0)

# levicivita_lut[a,b] contains the Levi-Cevita symbol ϵ_abc
# for c=6-a-b, i.e, when a,b,c is a permutation of 1,2,3
const levicivita_lut = [0 1 -1; -1 0 1; 1 -1 0]

function comm(A::σ,B::σ)
    if A.a == B.a
        scal(0)
    else
        a = Int(A.a)
        b = Int(B.a)
        # a+b+c == 6 (since a,b,c is a permutation of 1,2,3)
        c = 6 - a - b
        s = levicivita_lut[a,b]
        δ(A.inds,B.inds)*scal(2im*s)*σ(c,A.inds)
    end
end

# when multiplying (or commuting) with an operator with an index, take into account that the term in the sum with equal index has to be treated specially
*(A::Union{param,ExpVal,Corr,BaseOperator},B::OpSumAnalytic) = (B = ensure_compatible_sumind(B,A); OpSumAnalytic(B.ind,A*B.A))
# no need to check indices here since we just dispatch to another routine
*(A::OpSumAnalytic,B::OpProd) = (A*B.A)*B.B
*(A::OpSumAnalytic,B::Operator) = (A = ensure_compatible_sumind(A,B); OpSumAnalytic(A.ind,A.A*B))

comm(A::Union{param,ExpVal,Corr,BaseOperator},B::OpSumAnalytic) = (B = ensure_compatible_sumind(B,A); OpSumAnalytic(B.ind,comm(A,B.A)))
comm(A::OpSumAnalytic,B::OpSumAnalytic) = (B = ensure_compatible_sumind(B,A); OpSumAnalytic(B.ind,comm(A,B.A)))
comm(A::OpSumAnalytic,B::Operator) = -comm(B,A)


text:
@testset "Valuation" begin
    @testset "Correctness (zero valuation)" begin
        @polyvar x y
        f = [x^2 - 2, x + y - 1]

        tracker, starts = coretracker_startsolutions(
            f;
            homotopy = (g, f) -> StraightLineHomotopy(g, f; gamma = 1.0 + im),
            seed = 12345,
            min_step_size = eps()^2,
            log_homotopy = true,
        )

        S = collect(starts)
        state = tracker.state
        val = HC.Valuation(S[1])
        init!(tracker, S[1], 0.0, 28.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s)
            end
        end
        @test val.ν ≈ [0.0, 0.0] atol = 1e-8
        @test norm(val.ν̇) < 1e-2
        @test norm(val.ν̈) < 1e-2

        init!(val)
        init!(tracker, S[1], 0.0, 28.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s, tracker.predictor)
            end
        end
        @test val.ν ≈ [0.0, 0.0] atol = 1e-8
        @test norm(val.ν̇) < 1e-7
        @test norm(val.ν̈) < 1e-6

        init!(val)
        init!(tracker, S[2], 0.0, 28.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s, tracker.predictor)
            end
        end
        @test val.ν ≈ [0.0, 0.0] atol = 1e-8
        @test norm(val.ν̇) < 1e-7
        @test norm(val.ν̈) < 1e-6

        ## Multi-projective tracking

        # affine valuation
        @polyvar x y v w

        tracker, starts = coretracker_startsolutions(
            [x * y - 6 * v * w, x^2 - 5 * v^2],
            variable_groups = [(x, v), (y, w)],
            log_homotopy = true,
        )
        S = collect(starts)
        state = tracker.state
        val = HC.Valuation(S[1]; affine = true)
        init!(val)
        init!(tracker, S[1], 0.0, 30.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s, tracker.predictor)
            end
        end
        @test is_success(status(tracker))
        @test val.ν ≈ [0.0, 0.0] atol = 1e-6

        val = HC.Valuation(S[1]; affine = false)
        init!(val)
        init!(tracker, S[1], 0.0, 30.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s, tracker.predictor)
            end
        end
        @test is_success(status(tracker))
        @test val.ν ≈ [0.0, 0.0, 0.0, 0.0] atol = 1e-6
    end

    @testset "Correctness (non-zero valuation)" begin
        f = equations(griewank_osborne())
        tracker, starts = coretracker_startsolutions(
            f;
            homotopy = (g, f) -> StraightLineHomotopy(g, f; gamma = 1.0 + im),
            seed = 12345,
            min_step_size = eps()^2,
            log_homotopy = true,
        )

        S = collect(starts)
        state = tracker.state
        val = HC.Valuation(S[1])
        init!(tracker, S[1], 0.0, 20.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s)
            end
        end
        @test val.ν ≈ [-0.5, -1] atol = 1e-4
        @test norm(val.ν̇) < 1e-4
        @test norm(val.ν̈) < 1e-4

        # Use analytic estimates for ν̇ and ν̈
        init!(val)
        @test all(isnan, val.ν)
        init!(tracker, S[1], 0.0, 20.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s, tracker.predictor)
            end
        end
        @test val.ν ≈ [-0.5, -1] atol = 1e-4
        @test norm(val.ν̇) < 1e-4
        @test norm(val.ν̈) < 1e-4
        @test HC.judge(val; tol = 1e-4, tol_at_infinity = 1e-4) == HC.VAL_AT_INFINITY

        init!(val)
        init!(tracker, S[2], 0.0, 25.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s, tracker.predictor)
            end
        end
        @test val.ν ≈ [1 // 3, 2 // 3] atol = 1e-3
        @test norm(val.ν̇) < 1e-3
        @test norm(val.ν̈) < 1e-3
        @test HC.judge(val; tol = 1e-3) == HC.VAL_FINITE

        init!(val)
        init!(tracker, S[2], 0.0, 25.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s)
            end
        end
        @test val.ν ≈ [1 // 3, 2 // 3] atol = 1e-3
        @test norm(val.ν̇) < 1e-3
        @test norm(val.ν̈) < 1e-3
        @test HC.judge(val; tol = 1e-3) == HC.VAL_FINITE
        @test HC.judge(val; tol = 1e-10) == HC.VAL_INDECISIVE

        init!(val)
        init!(tracker, S[3], 0.0, 25.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s)
            end
        end
        @test val.ν ≈ [2, -1] atol = 1e-5
        @test norm(val.ν̇) < 1e-5
        @test norm(val.ν̈) < 1e-5
        @test HC.judge(val; tol_at_infinity = 1e-4) == HC.VAL_AT_INFINITY
    end

    @testset "Correctness projective (non-zero valuation)" begin
        f = homogenize(equations(griewank_osborne()))
        tracker, starts = coretracker_startsolutions(
            f;
            homotopy = (g, f) -> StraightLineHomotopy(g, f; gamma = 1.0 + im),
            min_step_size = eps()^2,
            log_homotopy = true,
        )

        S = collect(starts)
        state = tracker.state
        val = HC.Valuation(state.x; affine = true)
        init!(tracker, S[1], 0.0, 25.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s)
            end
        end
        @test val.ν ≈ [-0.5, -1] atol = 1e-4
        @test norm(val.ν̇) < 1e-4
        @test norm(val.ν̈) < 1e-4
        @test HC.judge(val; tol_at_infinity = 1e-4) == HC.VAL_AT_INFINITY

        # Use analytic estimates for ν̇ and ν̈
        init!(val)
        @test all(isnan, val.ν)
        init!(tracker, S[1], 0.0, 20.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s, tracker.predictor)
            end
        end
        @test val.ν ≈ [-0.5, -1] atol = 1e-4
        @test norm(val.ν̇) < 1e-4
        @test norm(val.ν̈) < 1e-4
        @test HC.judge(val; tol_at_infinity = 1e-4) == HC.VAL_AT_INFINITY
        @test HC.judge(val; tol_at_infinity = 1e-10) == HC.VAL_INDECISIVE

        # PROJECTIVE VALUATIONS

        val = HC.Valuation(state.x; affine = false)
        init!(tracker, S[1], 0.0, 25.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s)
            end
        end
        @test val.ν ≈ [0.5, 0, 1] atol = 1e-3
        @test norm(val.ν̇) < 1e-3
        @test norm(val.ν̈) < 1e-3
        @test HC.judge(val; tol = 1e-2) == HC.VAL_FINITE

        # Use analytic estimates for ν̇ and ν̈
        init!(val)
        @test all(isnan, val.ν)
        init!(tracker, S[1], 0.0, 25.0)
        for _ in tracker
            if !state.last_step_failed
                HC.update!(val, state.x, state.ẋ, state.s, tracker.predictor)
            end
        end
        @test val.ν ≈ [0.5, 0, 1] atol = 1e-3
        @test norm(val.ν̇) < 1e-3
        @test norm(val.ν̈) < 1e-3
        @test HC.judge(val; tol = 1e-2) == HC.VAL_FINITE
    end
end


text:
using LinearAlgebra

# k_local(Elastic Modulus,thickness,height,length,number of nodes,row/column one, row/column two)

#creates local stiffness matrix and maps to global
function k_local(k,l,n,rc_1,rc_2)
    k_global = zeros(n,n)
    k_global[rc_1,rc_1] = k
    k_global[rc_1,rc_2] = -k
    k_global[rc_2,rc_1] = -k
    k_global[rc_2,rc_2] = k
    return k_global
end

# constants from problem statement
E = 29*10^6
l = 2
k1 = 25
k2 = 35
k3 = 45
k4 = 30
# combining all globally mapped matrices
kG = k_local(k1,l,5,1,2)+
     k_local(k2,l,5,2,3)+
     k_local(k3,l,5,3,4)+
     k_local(k4,l,5,3,5)
kG[4,:] = Float64[0 0 0 1 0] #adjustment for Reaction from fixed end
kG[5,:] = Float64[0 0 0 0 1] #adjustment for Reaction from fixed end
kG[:,4] = transpose(Float64[0 0 0 1 0]) #adjustment for Reaction from fixed end
kG[:,5] = transpose(Float64[0 0 0 0 1]) #adjustment for Reaction from fixed end

# Force Vector
f = Float64[75 0 0 0 0]
F = transpose(f)

# Displacement vector
u = kG\F
println("Nodal displacement vector (inches):")
display(u)

# Strain Function
ϵ(x2,x1,l) = (x2 - x1)/l

#Strain Vector
ϵM = [ϵ(u[2],u[1],l)
      ϵ(u[3],u[2],l)
      ϵ(u[4],u[3],l)
      ϵ(u[5],u[3],l)]
println("Element strain vector (inch/inch):")
display(ϵM)

# Stress Function
σ(E,ϵ) = E*ϵ

# Stress Vector
σM = [σ(E,ϵM[1])
      σ(E,ϵM[2])
      σ(E,ϵM[3])
      σ(E,ϵM[4])]
println("Element Stress vector (psi):")
display(σM)

#Reaction Forces vector
Fapp(k,x2,x) = k*(x2-x)
FappM = [Fapp(k12,u[2],u[1])
         Fapp(k12,u[3],u[2])
         Fapp(k3,u[4],u[3])
         Fapp(k4,u[5],u[3])]
println("Reaction Forces vector (lbs):")
display(FappM)


text:
# Promotion with PolyVar and Monomial
Base.promote_rule(::Type{PolyVar{C}}, ::Type{PolyVar{C}}) where {C} = PolyVar{C}
Base.promote_rule(::Type{<:DMonomialLike{C}}, ::Type{<:DMonomialLike{C}}) where {C} = Monomial{C}
#promote_rule{S<:Union{Monomial, PolyVar}, T<:Union{Monomial, PolyVar}}(::Type{S}, ::Type{T}) = Monomial{iscomm{S}}

Base.promote_rule(::Type{Polynomial{C, S}}, ::Type{Polynomial{C, T}}) where {S, C, T} = Polynomial{C, promote_type(S, T)}
Base.promote_rule(::Type{<:AbstractPolynomialLike{S}}, ::Type{Polynomial{C, T}}) where {S, C, T} = Polynomial{C, promote_type(S, T)}
Base.promote_rule(::Type{Polynomial{C, T}}, ::Type{<:AbstractPolynomialLike{S}}) where {S, C, T} = Polynomial{C, promote_type(S, T)}
Base.promote_rule(::Type{<:DMonomialLike{S}}, ::Type{Polynomial{C, T}}) where {S, C, T} = Polynomial{C, promote_type(T, Int)}
Base.promote_rule(::Type{Polynomial{C, T}}, ::Type{<:DMonomialLike{S}}) where {S, C, T} = Polynomial{C, promote_type(T, Int)}

MP.promote_rule_constant(::Type{T}, ::Type{<:DMonomialLike{C}}) where {C, T} = Term{C, promote_type(T, Int)}
MP.promote_rule_constant(::Type{S}, ::Type{Term{C, T}}) where {S, C, T} = Term{C, promote_type(S, T)}
MP.promote_rule_constant(::Type{S}, ::Type{<:TermPoly{C, T}}) where {S, C, T} = Polynomial{C, promote_type(S, T)}

# Promotion with Term
Base.promote_rule(::Type{Term{C, S}}, ::Type{Term{C, T}}) where {C,S,T} = Term{C, promote_type(S, T)}
Base.promote_rule(::Type{<:DMonomialLike{C}}, ::Type{Term{C, T}}) where {C,T} = Term{C, promote_type(T, Int)}
Base.promote_rule(::Type{Term{C, T}}, ::Type{<:DMonomialLike{C}}) where {C,T} = Term{C, promote_type(T, Int)}

# Promotion with Polynomial
#Base.promote_rule(::Type{Polynomial{C, S}}, ::Type{Term{C, T}}) where {C, S, T} = Polynomial{C, promote_type(S, T)}
#Base.promote_rule(::Type{Term{C, T}}, ::Type{Polynomial{C, S}}) where {C, S, T} = Polynomial{C, promote_type(S, T)}

#Base.promote_rule(::Type{<:TermPoly{C, S}}, ::Type{<:TermPoly{C, T}}) where {C, S, T} = Polynomial{C, promote_type(S, T)}
#Base.promote_rule(::Type{<:TermPoly{C, T}}, ::Type{<:DMonomialLike{C}}) where {C, T} = Polynomial{C, promote_type(T, Int)}
#Base.promote_rule(::Type{<:DMonomialLike{C}}, ::Type{<:TermPoly{C, T}}) where {C, T} = Polynomial{C, promote_type(T, Int)}


text:
using Test: @test, @testset

using PartiallyObservedInverseGames.TrajectoryVisualization: visualize_trajectory
using PartiallyObservedInverseGames.ForwardGame: IBRGameSolver, KKTGameSolver, solve_game
using PartiallyObservedInverseGames.InverseGames:
    InverseKKTConstraintSolver, InverseKKTResidualSolver, solve_inverse_game

unique!(push!(LOAD_PATH, joinpath(@__DIR__, "utils")))
import TestUtils
import TestDynamics
import CollisionAvoidanceGame

control_system =
    TestDynamics.ProductSystem([TestDynamics.Unicycle(0.25), TestDynamics.Unicycle(0.25)])

x0 = vcat([-1, 0, 0.1, 0 + deg2rad(10)], [0, -1, 0.1, pi / 2 + deg2rad(10)])
T = 25
player_cost_models = let
    cost_model_p1 = CollisionAvoidanceGame.generate_player_cost_model(;
        player_idx = 1,
        control_system,
        T,
        goal_position = [1, 0],
    )
    cost_model_p2 = CollisionAvoidanceGame.generate_player_cost_model(;
        player_idx = 2,
        control_system,
        T,
        goal_position = [0, 1],
    )

    (cost_model_p1, cost_model_p2)
end

@testset "Product Dynamics" begin
    @test control_system.n_states == 8
    @test control_system.n_controls == 4

    @test all(i in TestDynamics.state_indices(control_system, 1) for i in 1:4)
    @test all(i in TestDynamics.state_indices(control_system, 2) for i in 5:8)
    @test all(i in TestDynamics.input_indices(control_system, 1) for i in 1:2)
    @test all(i in TestDynamics.input_indices(control_system, 2) for i in 3:4)
end

@testset "Forward Game" begin
    @testset "IBR" begin
        ibr_converged, ibr_solution, ibr_models =
            solve_game(IBRGameSolver(), control_system, player_cost_models, x0, T)
        global ibr_solution

        @test ibr_converged
    end

    @testset "KKT" begin
        kkt_converged, kkt_solution, kkt_model = solve_game(
            KKTGameSolver(),
            control_system,
            player_cost_models,
            x0,
            T;
            init = ibr_solution,
        )
        global kkt_solution

        @test kkt_converged
    end
end

@testset "Inverse Game" begin
    @testset "Perfect Observation" begin
        observation_model = (; σ = 0, expected_observation = identity)

        @testset "Inverse KKT Constraints" begin
            converged, inverse_kkt_solution, inverse_kkt_model = solve_inverse_game(
                InverseKKTConstraintSolver(),
                kkt_solution.x;
                control_system,
                observation_model,
                player_cost_models,
            )

            @test converged

            for (cost_model, weights) in
                zip(player_cost_models, inverse_kkt_solution.player_weights)
                TestUtils.test_inverse_solution(weights, cost_model.weights)
            end

            TestUtils.test_data_fidelity(
                inverse_kkt_model,
                observation_model,
                kkt_solution.x,
                kkt_solution.x,
            )
        end

        @testset "Invsere KKT Residuals" begin
            converged, inverse_kkt_solution, inverse_kkt_model = solve_inverse_game(
                InverseKKTResidualSolver(),
                kkt_solution.x,
                kkt_solution.u;
                control_system,
                player_cost_models,
            )

            @test converged

            for (cost_model, weights) in
                zip(player_cost_models, inverse_kkt_solution.player_weights)
                TestUtils.test_inverse_solution(weights, cost_model.weights)
            end

            TestUtils.test_data_fidelity(
                inverse_kkt_model,
                observation_model,
                kkt_solution.x,
                kkt_solution.x,
            )
        end
    end

    @testset "Noisy Full Observation" begin
        observation_model = (; σ = 0.01, expected_observation = identity)
        y_obs = TestUtils.noisy_observation(observation_model, kkt_solution.x)

        converged, inverse_kkt_solution, inverse_kkt_model = solve_inverse_game(
            InverseKKTConstraintSolver(),
            y_obs;
            control_system,
            observation_model,
            player_cost_models,
        )

        @test converged

        TestUtils.test_data_fidelity(inverse_kkt_model, observation_model, kkt_solution.x, y_obs)
    end

    @testset "Noise-Free Partial Observation" begin
        observation_model = (; σ = 0.0, expected_observation = x -> x[[1, 2, 4, 5, 6, 8], :])
        y_obs = TestUtils.noisy_observation(observation_model, kkt_solution.x)

        converged, inverse_kkt_solution, inverse_kkt_model = solve_inverse_game(
            InverseKKTConstraintSolver(),
            y_obs;
            control_system,
            observation_model,
            player_cost_models,
            init_with_observation = true,
            max_observation_error = 0.1,
        )

        @test converged

        TestUtils.test_data_fidelity(inverse_kkt_model, observation_model, kkt_solution.x, y_obs)
    end

    @testset "Noisy Partial Observation" begin
        observation_model = (; σ = 0.01, expected_observation = x -> x[[1, 2, 4, 5, 6, 8], :])
        y_obs = TestUtils.noisy_observation(observation_model, kkt_solution.x)

        converged, inverse_kkt_solution, inverse_kkt_model = solve_inverse_game(
            InverseKKTConstraintSolver(),
            y_obs;
            control_system,
            observation_model,
            player_cost_models,
            init_with_observation = true,
            max_observation_error = 0.1,
        )

        @test converged

        TestUtils.test_data_fidelity(inverse_kkt_model, observation_model, kkt_solution.x, y_obs)
    end
end