text:
"""
Baseclass for any kind of cost. Notable examples for
costs are `Loss` and `Penalty`.
"""
abstract type Cost end

"""
Baseclass for all losses. A loss is some (possibly simplified)
function `L(x, y, ŷ)`, of features `x`, targets `y` and outputs
`ŷ = f(x)` for some function `f`.
"""
abstract type Loss <: Cost end

"""
A loss is considered **supervised**, if all the information needed
to compute `L(x, y, ŷ)` are contained in `y` and `ŷ`, and thus allows
for the simplification `L(y, ŷ)`.
"""
abstract type SupervisedLoss <: Loss end

"""
A supervised loss that can be simplified to `L(y, ŷ) = L(ŷ - y)`
is considered **distance-based**.
"""
abstract type DistanceLoss <: SupervisedLoss end

"""
A supervised loss with targets `y ∈ {-1, 1}`, and which
can be simplified to `L(y, ŷ) = L(y⋅ŷ)` is considered
**margin-based**.
"""
abstract type MarginLoss <: SupervisedLoss end

"""
A loss is considered **unsupervised**, if all the information needed
to compute `L(x, y, ŷ)` are contained in `x` and `ŷ`, and thus allows
for the simplification `L(x, ŷ)`.
"""
abstract type UnsupervisedLoss <: Loss end

"""
Baseclass for all penalties.
"""
abstract type Penalty <: Cost end

"""
    value(loss, target, output) -> Number

Compute the (non-negative) numeric result for the `loss` function.
Note that `target` and `output` can be of different numeric type,
in which case promotion is performed in the manner appropriate for
the given loss.

    value(loss, targets, outputs) -> AbstractVector

Compute the result for each pair of values in `targets` and `outputs`.

    value(loss, targets, outputs, aggmode) -> Number

Compute the weighted or unweighted sum or mean (depending on
aggregation mode `aggmode`) of the individual values of the `loss`
function for each pair in `targets` and `outputs`. This method
will not allocate a temporary array.

In the case that the two parameters are arrays with a different
number of dimensions, broadcast will be performed. Note that the
given parameters are expected to have the same size in the
dimensions they share.

    value(loss, targets, outputs, aggmode, obsdim) -> AbstractVector

Compute the aggregated result along dimension with observations `obsdim`.
This method will not allocate a temporary array, but it will allocate
the resulting vector.

Both arrays have to be of the same shape and size. Furthermore they have
to have at least two dimensions (i.e. they must not be vectors).

## Notes

- New loss functions only need to implement the first method with single
  `target` and `output`. Fallback implementations are available for other
  methods in `LossFunctions.jl`.
"""
function value end

"""
    deriv(loss, target, output) -> Number

Compute the analytical derivative with respect to the `output` for the
`loss` function. Note that `target` and `output` can be of different
numeric type, in which case promotion is performed in the manner
appropriate for the given loss.

    deriv(loss, targets, outputs) -> AbstractVector

Compute the result for each pair of values in `targets` and `outputs`.

    deriv(loss, targets, outputs, aggmode) -> Number

Compute the weighted or unweighted sum or mean (depending on
aggregation mode `aggmode`) of the individual values of the `loss`
function for each pair in `targets` and `outputs`. This method
will not allocate a temporary array.

In the case that the two parameters are arrays with a different
number of dimensions, broadcast will be performed. Note that the
given parameters are expected to have the same size in the
dimensions they share.

    deriv(loss, targets, outputs, aggmode, obsdim) -> AbstractVector

Compute the aggregated result along dimension with observations `obsdim`.
This method will not allocate a temporary array, but it will allocate
the resulting vector.

Both arrays have to be of the same shape and size. Furthermore they have
to have at least two dimensions (i.e. they must not be vectors).

## Notes

- New loss functions only need to implement the first method with single
  `target` and `output`. Fallback implementations are available for other
  methods in `LossFunctions.jl`.
"""
function deriv end

"""
    deriv2(loss, target, output) -> Number

Compute the second derivative with respect to the `output` for the
`loss` function. Note that `target` and `output` can be of different
numeric type, in which case promotion is performed in the manner
appropriate for the given loss.

    deriv2(loss, targets, outputs) -> AbstractVector

Compute the result for each pair of values in `targets` and `outputs`.

    deriv2(loss, targets, outputs, aggmode) -> Number

Compute the weighted or unweighted sum or mean (depending on
aggregation mode `aggmode`) of the individual values of the `loss`
function for each pair in `targets` and `outputs`. This method
will not allocate a temporary array.

In the case that the two parameters are arrays with a different
number of dimensions, broadcast will be performed. Note that the
given parameters are expected to have the same size in the
dimensions they share.

    deriv2(loss, targets, outputs, aggmode, obsdim) -> AbstractVector

Compute the aggregated result along dimension with observations `obsdim`.
This method will not allocate a temporary array, but it will allocate
the resulting vector.

Both arrays have to be of the same shape and size. Furthermore they have
to have at least two dimensions (i.e. they must not be vectors).

## Notes

- New loss functions only need to implement the first method with single
  `target` and `output`. Fallback implementations are available for other
  methods in `LossFunctions.jl`.
"""
function deriv2 end

"""
    isconvex(loss) -> Bool

Return `true` if the given `loss` denotes a convex function.
A function `f: ℝⁿ → ℝ` is convex if its domain is a convex set
and if for all `x, y` in that domain, with `θ` such that for
`0 ≦ θ ≦ 1`, we have `f(θ x + (1 - θ) y) ≦ θ f(x) + (1 - θ) f(y)`.
"""
isconvex(loss::SupervisedLoss) = isstrictlyconvex(loss)

"""
    isstrictlyconvex(loss) -> Bool

Return `true` if the given `loss` denotes a strictly convex function.
A function `f : ℝⁿ → ℝ` is strictly convex if its domain is a convex
set and if for all `x, y` in that domain where `x ≠ y`, with `θ` such
that for `0 < θ < 1`, we have `f(θ x + (1 - θ) y) < θ f(x) + (1 - θ) f(y)`.
"""
isstrictlyconvex(loss::SupervisedLoss) = isstronglyconvex(loss)

"""
    isstronglyconvex(loss) -> Bool

Return `true` if the given `loss` denotes a strongly convex function.
A function `f : ℝⁿ → ℝ` is `m`-strongly convex if its domain is a convex
set, and if for all `x, y` in that domain where `x ≠ y`, and `θ` such that
for `0 ≤ θ ≤ 1`, we have
`f(θ x + (1 - θ)y) < θ f(x) + (1 - θ) f(y) - 0.5 m ⋅ θ (1 - θ) | x - y |₂²`

In a more familiar setting, if the loss function is differentiable we have
`(∇f(x) - ∇f(y))ᵀ (x - y) ≥ m | x - y |₂²`
"""
isstronglyconvex(loss::SupervisedLoss) = false

"""
    isdifferentiable(loss, [x]) -> Bool

Return `true` if the given `loss` is differentiable
(optionally limited to the given point `x` if specified).

A function `f : ℝⁿ → ℝᵐ` is differentiable at a point `x` in the interior
domain of `f` if there exists a matrix `Df(x) ∈ ℝ^(m × n)` such that
it satisfies:

`lim_{z ≠ x, z → x} (|f(z) - f(x) - Df(x)(z-x)|₂) / |z - x|₂ = 0`

A function is differentiable if its domain is open and it is
differentiable at every point `x`.
"""
isdifferentiable(loss::SupervisedLoss) = istwicedifferentiable(loss)
isdifferentiable(loss::SupervisedLoss, at) = isdifferentiable(loss)

"""
    istwicedifferentiable(loss, [x]) -> Bool

Return `true` if the given `loss` is differentiable
(optionally limited to the given point `x` if specified).

A function `f : ℝⁿ → ℝ` is said to be twice differentiable
at a point `x` in the interior domain of `f`, if the function
derivative for `∇f` exists at `x`: `∇²f(x) = D∇f(x)`.

A function is twice differentiable if its domain is open and it
is twice differentiable at every point `x`.
"""
istwicedifferentiable(loss::SupervisedLoss) = false
istwicedifferentiable(loss::SupervisedLoss, at) = istwicedifferentiable(loss)

"""
    islocallylipschitzcont(loss) -> Bool

Return `true` if the given `loss` function is locally-Lipschitz
continous.

A supervised loss `L : Y × ℝ → [0, ∞)` is called locally Lipschitz
continuous if for all `a ≥ 0` there exists a constant `cₐ ≥ 0`,
such that

`sup_{y ∈ Y} | L(y,t) − L(y,t′) | ≤ cₐ |t − t′|, t, t′ ∈ [−a,a]`

Every convex function is locally lipschitz continuous.
"""
islocallylipschitzcont(loss::SupervisedLoss) =
    isconvex(loss) || islipschitzcont(loss)

"""
    islipschitzcont(loss) -> Bool

Return `true` if the given `loss` function is Lipschitz continuous.

A supervised loss function `L : Y × ℝ → [0, ∞)` is Lipschitz continous,
if there exists a finite constant `M < ∞` such that
`|L(y, t) - L(y, t′)| ≤ M |t - t′|, ∀ (y, t) ∈ Y × ℝ`
"""
islipschitzcont(loss::SupervisedLoss) = false

"""
    isnemitski(loss) -> Bool

Return `true` if the given `loss` denotes a Nemitski loss function.

We call a supervised loss function `L : Y × ℝ → [0,∞)` a Nemitski
loss if there exist a measurable function `b : Y → [0, ∞)` and an
increasing function `h : [0, ∞) → [0, ∞)` such that
`L(y,ŷ) ≤ b(y) + h(|ŷ|), (y, ŷ) ∈ Y × ℝ`

If a loss if locally lipsschitz continuous then it is a Nemitski loss.
"""
isnemitski(loss::SupervisedLoss) = islocallylipschitzcont(loss)
isnemitski(loss::MarginLoss) = true

"""
    isunivfishercons(loss) -> Bool
"""
isunivfishercons(loss::Loss) = false

"""
    isfishercons(loss) -> Bool

Return `true` if the givel `loss` is Fisher consistent.

We call a supervised loss function `L : Y × ℝ → [0,∞)` a Fisher
consistent loss if the population minimizer of the risk `E[L(y,f(x))]`
for all measurable functions leads to the Bayes optimal decision rule.
"""
isfishercons(loss::Loss) = isunivfishercons(loss)

"""
    isclipable(loss) -> Bool

Return `true` if the given `loss` function is clipable. A
supervised loss `L : Y × ℝ → [0,∞)` can be clipped at `M > 0`
if, for all `(y,t) ∈ Y × ℝ`, `L(y, t̂) ≤ L(y, t)` where
`t̂` denotes the clipped value of `t` at `± M`.
That is
`t̂ = -M` if `t < -M`, `t̂ = t` if `t ∈ [-M, M]`, and `t = M` if `t > M`.
"""
isclipable(loss::SupervisedLoss) = false
isclipable(loss::DistanceLoss) = true # can someone please double check?

"""
    isdistancebased(loss) -> Bool

Return `true` if the given `loss` is a distance-based loss.

A supervised loss function `L : Y × ℝ → [0,∞)` is said to be
distance-based, if there exists a representing function `ψ : ℝ → [0,∞)`
satisfying `ψ(0) = 0` and `L(y, ŷ) = ψ (ŷ - y), (y, ŷ) ∈ Y × ℝ`.
"""
isdistancebased(loss::Loss) = false
isdistancebased(loss::DistanceLoss) = true

"""
    ismarginbased(loss) -> Bool

Return `true` if the given `loss` is a margin-based loss.

A supervised loss function `L : Y × ℝ → [0,∞)` is said to be
margin-based, if there exists a representing function `ψ : ℝ → [0,∞)`
satisfying `L(y, ŷ) = ψ(y⋅ŷ), (y, ŷ) ∈ Y × ℝ`.
"""
ismarginbased(loss::Loss) = false
ismarginbased(loss::MarginLoss) = true

"""
    isclasscalibrated(loss) -> Bool
"""
isclasscalibrated(loss::SupervisedLoss) = false
isclasscalibrated(loss::MarginLoss) =
    isconvex(loss) && isdifferentiable(loss, 0) && deriv(loss, 0) < 0

"""
    issymmetric(loss) -> Bool

Return `true` if the given loss is a symmetric loss.

A function `f : ℝ → [0,∞)` is said to be symmetric
about origin if we have `f(x) = f(-x), ∀ x ∈ ℝ`.

A distance-based loss is said to be symmetric if its
representing function is symmetric.
"""
issymmetric(loss::SupervisedLoss) = false

"""
    isminimizable(loss) -> Bool

Return `true` if the given `loss` is a minimizable loss.
"""
isminimizable(loss::SupervisedLoss) = isconvex(loss)


text:
import SimpleSMC: particleCopy!, _copyParticles!, _mutateParticles!,
  _logWeightParticles!, _intermediateOutput!, pickParticle!,
  _sampleCategoricalSorted!
import NonUniformRandomVariateGeneration: sampleCategorical, sampleBinomial

function _indexCoupledMutateParticles!(zetas1::Vector{Particle},
  zetas2::Vector{Particle}, M!::F, p::Int64, zetaAncs1::Vector{Particle},
  zetaAncs2::Vector{Particle}, pScratch::ParticleScratch, xref1::Particle,
  xref2::Particle) where {Particle, F<:Function, ParticleScratch}
  @inbounds particleCopy!(zetas1[1], xref1)
  @inbounds particleCopy!(zetas2[1], xref2)
  rng = getRNG()
  for j in 2:length(zetas1)
    @inbounds M!(zetas1[j], rng, p, zetaAncs1[j], pScratch)
    if p == 1 || zetaAncs1[j] == zetaAncs2[j]
      @inbounds particleCopy!(zetas2[j], zetas1[j])
    else
      @inbounds M!(zetas2[j], rng, p, zetaAncs2[j], pScratch)
    end
  end
end

function _rngCoupledMutateParticles!(zetas1::Vector{Particle},
  zetas2::Vector{Particle}, M!::F, p::Int64, zetaAncs1::Vector{Particle},
  zetaAncs2::Vector{Particle}, pScratch::ParticleScratch, xref1::Particle,
  xref2::Particle, extraRNG::RNG) where {Particle, F<:Function, ParticleScratch}
  @inbounds particleCopy!(zetas1[1], xref1)
  @inbounds particleCopy!(zetas2[1], xref2)
  rng = getRNG()
  copyto!(extraRNG, rng)
  for j in 2:length(zetas1)
    @inbounds M!(zetas1[j], rng, p, zetaAncs1[j], pScratch)
    @inbounds M!(zetas2[j], extraRNG, p, zetaAncs2[j], pScratch)
  end
end

function _computeWeights(ccsmcio::CCSMCIO, ws1::Vector{Float64},
  ws2::Vector{Float64})
  nws1 = ccsmcio.nws1
  nws2 = ccsmcio.nws2
  cws = ccsmcio.cws
  rws1 = ccsmcio.rws1
  rws2 = ccsmcio.rws2

  nws1 .= ws1
  nws1 ./= sum(nws1)

  nws2 .= ws2
  nws2 ./= sum(nws2)

  for i in 1:ccsmcio.N
    @inbounds cws[i] = min(nws1[i], nws2[i])
  end

  rws1 .= nws1 .- cws
  rws2 .= nws2 .- cws
end

function _coupledResample!(ccsmcio::CCSMCIO)
  _computeWeights(ccsmcio, ccsmcio.smcio1.ws, ccsmcio.smcio2.ws)
  rng = getRNG()
  cProb = sum(ccsmcio.cws)

  N::Int64 = ccsmcio.N
  nCommon::Int64 = sampleBinomial(N-1, cProb, rng)

  as1 = ccsmcio.smcio1.internal.as
  as2 = ccsmcio.smcio2.internal.as
  scratch1 = ccsmcio.smcio1.internal.scratch1
  scratch2 = ccsmcio.smcio1.internal.scratch2

  as1[1] = 1
  _sampleCategoricalSorted!(as1, ccsmcio.cws, nCommon, scratch1, scratch2, 1,
    rng)

  for i in 1:1+nCommon
    as2[i] = as1[i]
  end

  if nCommon < N - 1
    _sampleCategoricalSorted!(as1, ccsmcio.rws1, N-nCommon-1, scratch1,
      scratch2, 1+nCommon, rng)
    _sampleCategoricalSorted!(as2, ccsmcio.rws2, N-nCommon-1, scratch1,
      scratch2, 1+nCommon, rng)
  end

end

function _coupledSampleIndices(ccsmcio, ws1, ws2)
  _computeWeights(ccsmcio, ws1, ws2)
  rng = getRNG()
  cProb = sum(ccsmcio.cws)
  if rand(rng) < cProb
    idx = sampleCategorical(ccsmcio.cws, rng)
    return (idx, idx)
  else
    idx1 = sampleCategorical(ccsmcio.rws1, rng)
    idx2 = sampleCategorical(ccsmcio.rws2, rng)
    return (idx1, idx2)
  end
end

function _trace!(path::Vector{Particle}, smcio::SMCIO{Particle}, k) where Particle
  n::Int64 = smcio.n
  allAs::Vector{Vector{Int64}} = smcio.allAs
  allZetas::Vector{Vector{Particle}} = smcio.allZetas
  @inbounds particleCopy!(path[n], allZetas[n][k])
  for p = n-1:-1:1
    @inbounds k = allAs[p][k]
    @inbounds particleCopy!(path[p], allZetas[p][k])
  end
end

function _pickParticles!(ccsmcio)
  ws1 = ccsmcio.smcio1.ws
  ws2 = ccsmcio.smcio2.ws
  k1, k2 = _coupledSampleIndices(ccsmcio, ws1, ws2)
  _trace!(ccsmcio.ref1, ccsmcio.smcio1, k1)
  _trace!(ccsmcio.ref2, ccsmcio.smcio2, k2)
end

function _pickParticlesBS!(ccsmcio::CCSMCIO{Particle}, lM::F) where
  {Particle, F<:Function}
  allZetas1::Vector{Vector{Particle}} = ccsmcio.smcio1.allZetas
  allZetas2::Vector{Vector{Particle}} = ccsmcio.smcio2.allZetas
  allWs1::Vector{Vector{Float64}} = ccsmcio.smcio1.allWs
  allWs2::Vector{Vector{Float64}} = ccsmcio.smcio2.allWs
  ref1 = ccsmcio.ref1
  ref2 = ccsmcio.ref2
  n = ccsmcio.n
  N = ccsmcio.N
  pScratch = ccsmcio.smcio1.internal.particleScratch

  ws1 = ccsmcio.smcio1.ws
  ws2 = ccsmcio.smcio2.ws
  k1, k2 = _coupledSampleIndices(ccsmcio, ws1, ws2)
  particleCopy!(ref1[n], allZetas1[n][k1])
  particleCopy!(ref2[n], allZetas2[n][k2])

  bws1 = ccsmcio.bws1
  bws2 = ccsmcio.bws2

  for p=n-1:-1:1
    bws1 .= log.(allWs1[p])
    bws2 .= log.(allWs2[p])
    for j in 1:N
      bws1[j] += lM(p+1, allZetas1[p][j], ref1[p+1], pScratch)
      bws2[j] += lM(p+1, allZetas2[p][j], ref2[p+1], pScratch)
    end
    m = maximum(bws1)
    bws1 .= exp.(bws1 .- m)
    m = maximum(bws2)
    bws2 .= exp.(bws2 .- m)
    k1, k2 = _coupledSampleIndices(ccsmcio, bws1, bws2)
    particleCopy!(ref1[p], allZetas1[p][k1])
    particleCopy!(ref2[p], allZetas2[p][k2])
  end
end

## index-coupled ancestor sampling
## this is slightly different to
## Jacob, P.E., Lindsten, F. and Schön, T.B., 2018. Smoothing with couplings of
## conditional particle filters. J. of the Amer. Statist. Assoc. 2018.
## who appear to use quantile coupling of the ancestor index in their code
function _ancestorSample!(ccsmcio, p, lM::F) where F<:Function
  zetas1 = ccsmcio.smcio1.zetas
  zetas2 = ccsmcio.smcio2.zetas
  zetaRef1 = ccsmcio.ref1[p+1]
  zetaRef2 = ccsmcio.ref2[p+1]
  bws1 = ccsmcio.bws1
  bws2 = ccsmcio.bws2
  ws1 = ccsmcio.smcio1.ws
  ws2 = ccsmcio.smcio2.ws
  pScratch = ccsmcio.smcio1.internal.particleScratch

  N = ccsmcio.N

  bws1 .= log.(ws1)
  bws2 .= log.(ws2)
  for j in 1:N
    bws1[j] += lM(p+1, zetas1[j], zetaRef1, pScratch)
    bws2[j] += lM(p+1, zetas2[j], zetaRef2, pScratch)
  end
  m = maximum(bws1)
  bws1 .= exp.(bws1 .- m)
  m = maximum(bws2)
  bws2 .= exp.(bws2 .- m)
  k1, k2 = _coupledSampleIndices(ccsmcio, bws1, bws2)

  ccsmcio.smcio1.internal.as[1] = k1
  ccsmcio.smcio2.internal.as[1] = k2
end

function checkEqual(v1::Vector{Particle}, v2::Vector{Particle}) where Particle
  for i = 1:length(v1)
    if v1[i] != v2[i]
      return false
    end
  end
  return true
end

ext: jl


text:
# syntax: proto3
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

type __enum_ExchangeType <: ProtoEnum
    APPENDENTRIES::Int32
    REQUESTVOTE::Int32
    INSTALLSNAPSHOT::Int32
    __enum_ExchangeType() = new(0,1,2)
end #type __enum_ExchangeType
const ExchangeType = __enum_ExchangeType()

type RPCRequest
    _type::Int32
    server_id::UInt64
    request::Array{UInt8,1}
    RPCRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type RPCRequest
hash(v::RPCRequest) = ProtoBuf.protohash(v)
isequal(v1::RPCRequest, v2::RPCRequest) = ProtoBuf.protoisequal(v1, v2)
==(v1::RPCRequest, v2::RPCRequest) = ProtoBuf.protoeq(v1, v2)

type RPCReply
    _type::Int32
    server_id::UInt64
    reply::Array{UInt8,1}
    RPCReply(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type RPCReply
hash(v::RPCReply) = ProtoBuf.protohash(v)
isequal(v1::RPCReply, v2::RPCReply) = ProtoBuf.protoisequal(v1, v2)
==(v1::RPCReply, v2::RPCReply) = ProtoBuf.protoeq(v1, v2)

type LogEntry
    index::UInt64
    term::UInt64
    key::AbstractString
    value::Array{UInt8,1}
    LogEntry(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type LogEntry
hash(v::LogEntry) = ProtoBuf.protohash(v)
isequal(v1::LogEntry, v2::LogEntry) = ProtoBuf.protoisequal(v1, v2)
==(v1::LogEntry, v2::LogEntry) = ProtoBuf.protoeq(v1, v2)

type AppendEntriesRequest
    term::UInt64
    leaderId::UInt64
    prevLogIndex::UInt64
    prevLogTerm::UInt64
    entry::Array{LogEntry,1}
    leaderCommit::UInt64
    AppendEntriesRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type AppendEntriesRequest
hash(v::AppendEntriesRequest) = ProtoBuf.protohash(v)
isequal(v1::AppendEntriesRequest, v2::AppendEntriesRequest) = ProtoBuf.protoisequal(v1, v2)
==(v1::AppendEntriesRequest, v2::AppendEntriesRequest) = ProtoBuf.protoeq(v1, v2)

type AppendEntriesReply
    term::UInt64
    success::Bool
    AppendEntriesReply(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type AppendEntriesReply
hash(v::AppendEntriesReply) = ProtoBuf.protohash(v)
isequal(v1::AppendEntriesReply, v2::AppendEntriesReply) = ProtoBuf.protoisequal(v1, v2)
==(v1::AppendEntriesReply, v2::AppendEntriesReply) = ProtoBuf.protoeq(v1, v2)

type RequestVoteRequest
    term::UInt64
    candidateId::UInt64
    lastLogIndex::UInt64
    lastLogTerm::UInt64
    RequestVoteRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type RequestVoteRequest
hash(v::RequestVoteRequest) = ProtoBuf.protohash(v)
isequal(v1::RequestVoteRequest, v2::RequestVoteRequest) = ProtoBuf.protoisequal(v1, v2)
==(v1::RequestVoteRequest, v2::RequestVoteRequest) = ProtoBuf.protoeq(v1, v2)

type RequestVoteReply
    term::UInt64
    voteGranted::Bool
    RequestVoteReply(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type RequestVoteReply
hash(v::RequestVoteReply) = ProtoBuf.protohash(v)
isequal(v1::RequestVoteReply, v2::RequestVoteReply) = ProtoBuf.protoisequal(v1, v2)
==(v1::RequestVoteReply, v2::RequestVoteReply) = ProtoBuf.protoeq(v1, v2)

type InstallSnapshotRequest
    term::UInt64
    leaderId::UInt64
    lastIncludedIndex::UInt64
    lastIncludedTerm::UInt64
    offset::UInt64
    data::Array{UInt8,1}
    done::Bool
    InstallSnapshotRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type InstallSnapshotRequest
hash(v::InstallSnapshotRequest) = ProtoBuf.protohash(v)
isequal(v1::InstallSnapshotRequest, v2::InstallSnapshotRequest) = ProtoBuf.protoisequal(v1, v2)
==(v1::InstallSnapshotRequest, v2::InstallSnapshotRequest) = ProtoBuf.protoeq(v1, v2)

type InstallSnapshotReply
    term::UInt64
    InstallSnapshotReply(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type InstallSnapshotReply
hash(v::InstallSnapshotReply) = ProtoBuf.protohash(v)
isequal(v1::InstallSnapshotReply, v2::InstallSnapshotReply) = ProtoBuf.protoisequal(v1, v2)
==(v1::InstallSnapshotReply, v2::InstallSnapshotReply) = ProtoBuf.protoeq(v1, v2)

export ExchangeType, RPCRequest, RPCReply, LogEntry, AppendEntriesRequest, AppendEntriesReply, RequestVoteRequest, RequestVoteReply, InstallSnapshotRequest, InstallSnapshotReply


text:
"""
    cudnnRNNForward(w, x; hiddenSize, o...)
    cudnnRNNForward!(y, w, x; hiddenSize, o...)
    cudnnRNNForward(w, x, d::cudnnRNNDescriptor; o...)
    cudnnRNNForward!(y, w, x, d::cudnnRNNDescriptor; o...)

Apply the RNN specified with weights `w` and configuration given by `d` or keyword options
to input `x`.

Keyword arguments for hidden input/output:
* `hx=nothing`: initialize the hidden vector if specified (by default initialized to 0).
* `cx=nothing`: initialize the cell vector (only in LSTMs) if specified (by default initialized to 0).
* `hy=nothing`: return the final hidden vector in `hy` if set to `Ref{Any}()`.
* `cy=nothing`: return the final cell vector in `cy` (only in LSTMs) if set to `Ref{Any}()`.

Keyword arguments specifying the RNN when `d::cudnnRNNDescriptor` is not given:
* `hiddenSize::Integer`: hidden vector size, which must be supplied when `d` is not given
* `algo::cudnnRNNAlgo_t = CUDNN_RNN_ALGO_STANDARD`: RNN algo (`CUDNN_RNN_ALGO_STANDARD`, `CUDNN_RNN_ALGO_PERSIST_STATIC`, or `CUDNN_RNN_ALGO_PERSIST_DYNAMIC`).
* `cellMode::cudnnRNNMode_t = CUDNN_LSTM`: Specifies the RNN cell type in the entire model (`CUDNN_RNN_RELU`, `CUDNN_RNN_TANH`, `CUDNN_LSTM`, `CUDNN_GRU`).
* `biasMode::cudnnRNNBiasMode_t = CUDNN_RNN_DOUBLE_BIAS`: Sets the number of bias vectors (`CUDNN_RNN_NO_BIAS`, `CUDNN_RNN_SINGLE_INP_BIAS`, `CUDNN_RNN_SINGLE_REC_BIAS`, `CUDNN_RNN_DOUBLE_BIAS`). The two single bias settings are functionally the same for `RELU`, `TANH` and `LSTM` cell types. For differences in `GRU` cells, see the description of `CUDNN_GRU` in cudnn docs.
* `dirMode::cudnnDirectionMode_t = CUDNN_UNIDIRECTIONAL`: Specifies the recurrence pattern: `CUDNN_UNIDIRECTIONAL` or `CUDNN_BIDIRECTIONAL`. In bidirectional RNNs, the hidden states passed between physical layers are concatenations of forward and backward hidden states.
* `inputMode::cudnnRNNInputMode_t = CUDNN_LINEAR_INPUT`: Specifies how the input to the RNN model is processed by the first layer. When `inputMode` is `CUDNN_LINEAR_INPUT`, original input vectors of size `inputSize` are multiplied by the weight matrix to obtain vectors of `hiddenSize`. When `inputMode` is `CUDNN_SKIP_INPUT`, the original input vectors to the first layer are used as is without multiplying them by the weight matrix.
* `mathPrec::DataType = eltype(x)`: This parameter is used to control the compute math precision in the RNN model. For `Float16` input/output can be `Float16` or `Float32`, for `Float32` or `Float64` input/output, must match the input/output type.
* `mathType::cudnnMathType_t = math_mode()`: Sets the preferred option to use NVIDIA Tensor Cores accelerators on Volta (SM 7.0) or higher GPUs. When `dataType` is `CUDNN_DATA_HALF`, the `mathType` parameter can be `CUDNN_DEFAULT_MATH` or `CUDNN_TENSOR_OP_MATH`. The `ALLOW_CONVERSION` setting is treated the same `CUDNN_TENSOR_OP_MATH` for this data type. When `dataType` is `CUDNN_DATA_FLOAT`, the `mathType` parameter can be `CUDNN_DEFAULT_MATH` or `CUDNN_TENSOR_OP_MATH_ALLOW_CONVERSION`. When the latter settings are used, original weights and intermediate results will be down-converted to `CUDNN_DATA_HALF` before they are used in another recursive iteration. When dataType is `CUDNN_DATA_DOUBLE`, the `mathType` parameter can be `CUDNN_DEFAULT_MATH`.
* `inputSize::Integer = size(x,1)`: Size of the input vector in the RNN model. When the `inputMode=CUDNN_SKIP_INPUT`, the `inputSize` should match the `hiddenSize` value.
* `projSize::Integer = hiddenSize`: The size of the LSTM cell output after the recurrent projection. This value should not be larger than `hiddenSize`. It is legal to set `projSize` equal to `hiddenSize`, however, in this case, the recurrent projection feature is disabled. The recurrent projection is an additional matrix multiplication in the LSTM cell to project hidden state vectors `ht` into smaller vectors `rt = Wr * ht`, where `Wr` is a rectangular matrix with `projSize` rows and hiddenSize columns. When the recurrent projection is enabled, the output of the LSTM cell (both to the next layer and unrolled in-time) is `rt` instead of `ht`. The recurrent projection can be enabled for LSTM cells and `CUDNN_RNN_ALGO_STANDARD` only.
* `numLayers::Integer = 1`: Number of stacked, physical layers in the deep RNN model. When `dirMode= CUDNN_BIDIRECTIONAL`, the physical layer consists of two pseudo-layers corresponding to forward and backward directions.
* `dropout::Real = 0`: When non-zero, dropout operation will be applied between physical layers. A single layer network will have no dropout applied. Dropout is used in the training mode only.
* `auxFlags::Integer = CUDNN_RNN_PADDED_IO_ENABLED`: Miscellaneous switches that do not require additional numerical values to configure the corresponding feature. In future cuDNN releases, this parameter will be used to extend the RNN functionality without adding new API functions (applicable options should be bitwise OR-ed). Currently, this parameter is used to enable or disable padded input/output (`CUDNN_RNN_PADDED_IO_DISABLED`, `CUDNN_RNN_PADDED_IO_ENABLED`). When the padded I/O is enabled, layouts `CUDNN_RNN_DATA_LAYOUT_SEQ_MAJOR_UNPACKED` and `CUDNN_RNN_DATA_LAYOUT_BATCH_MAJOR_UNPACKED` are permitted in RNN data descriptors.

Other keyword arguments:
* `layout::cudnnRNNDataLayout_t = CUDNN_RNN_DATA_LAYOUT_SEQ_MAJOR_UNPACKED`: The memory layout of the RNN data tensor. Options are `CUDNN_RNN_DATA_LAYOUT_SEQ_MAJOR_UNPACKED`: Data layout is padded, with outer stride from one time-step to the next; `CUDNN_RNN_DATA_LAYOUT_SEQ_MAJOR_PACKED`: The sequence length is sorted and packed as in the basic RNN API; `CUDNN_RNN_DATA_LAYOUT_BATCH_MAJOR_UNPACKED`: Data layout is padded, with outer stride from one batch to the next.
* `seqLengthArray::Vector{Cint} = nothing`: An integer array with `batchSize` number of elements. Describes the length (number of time-steps) of each sequence. Each element in `seqLengthArray` must be greater than or equal to 0 but less than or equal to `maxSeqLength`. In the packed layout, the elements should be sorted in descending order, similar to the layout required by the non-extended RNN compute functions. The default value `nothing` assumes uniform `seqLength`s, no padding.
* `devSeqLengths::CuVector{Cint} = nothing`: Device copy of `seqLengthArray`
* `fwdMode::cudnnForwardMode_t = CUDNN_FWD_MODE_INFERENCE`: set to `CUDNN_FWD_MODE_TRAINING` when training
"""
cudnnRNNForward, cudnnRNNForward!

# Public methods
cudnnRNNForward(w, x; hiddenSize, o...)     = cudnnRNNForwardWithDefaults(w, x; hiddenSize, o...)
cudnnRNNForward!(y, w, x; hiddenSize, o...) = cudnnRNNForwardWithDefaults(w, x; y, hiddenSize, o...)
cudnnRNNForward(w, x, d::cudnnRNNDescriptor; o...)     = cudnnRNNForwardWithDefaults(w, x; rnnDesc=d, o...)
cudnnRNNForward!(y, w, x, d::cudnnRNNDescriptor; o...) = cudnnRNNForwardWithDefaults(w, x; y, rnnDesc=d, o...)

# Private method
function cudnnRNNForwardWithDefaults(
    w, x;

    # input hidden vectors
    hx = nothing,
    cx = nothing,

    # output buffers
    y = nothing,
    hy = nothing,
    cy = nothing,

    # rnnDescriptor parameters
    # TODO: look into GetClip, SetClip
    algo::cudnnRNNAlgo_t = CUDNN_RNN_ALGO_STANDARD,
    cellMode::cudnnRNNMode_t = CUDNN_LSTM,
    biasMode::cudnnRNNBiasMode_t = CUDNN_RNN_DOUBLE_BIAS,
    dirMode::cudnnDirectionMode_t = CUDNN_UNIDIRECTIONAL,
    inputMode::cudnnRNNInputMode_t = CUDNN_LINEAR_INPUT,
    dataType::DataType = eltype(x),
    mathPrec::DataType = dataType, # has to match dataType with one extra possibility dt=Float16 => mp=Float16|Float32
    mathType::cudnnMathType_t = math_mode(),
    inputSize::Integer = size(x,1),
    hiddenSize::Integer = 0,
    projSize::Integer = hiddenSize,
    numLayers::Integer = 1,
    dropout::Real = 0,
    auxFlags::Integer = CUDNN_RNN_PADDED_IO_ENABLED, # When the padded I/O is enabled, layouts CUDNN_RNN_DATA_LAYOUT_SEQ_MAJOR_UNPACKED and CUDNN_RNN_DATA_LAYOUT_BATCH_MAJOR_UNPACKED are permitted in RNN data descriptors.

    # rnnDescriptor
    rnnDesc::cudnnRNNDescriptor = cudnnRNNDescriptor(algo, cellMode, biasMode, dirMode, inputMode, cudnnDataType(dataType), cudnnDataType(mathPrec), mathType, Int32(inputSize), checkHidden(hiddenSize), Int32(projSize), Int32(numLayers), cudnnDropoutDescriptor(Cfloat(dropout)), UInt32(auxFlags)),

    # rnnData parameters:
    layout::cudnnRNNDataLayout_t = CUDNN_RNN_DATA_LAYOUT_SEQ_MAJOR_UNPACKED, # padded [X,B,T] array
    seqLengthArray::Union{Nothing,Vector{Cint}} = nothing,         # assume no padding by default
    paddingFill::Ptr{Cvoid} = C_NULL,

    # rnnForward parameters
    fwdMode::cudnnForwardMode_t = CUDNN_FWD_MODE_INFERENCE, # set to CUDNN_FWD_MODE_TRAINING when training
    devSeqLengths::Union{Nothing,CuArray{Cint,1}} = nothing,
    reserveSpace::Union{CuArray,Nothing} = nothing,
    workspace::Union{CuArray,Nothing} = nothing,

    # gradient buffers: layer designers may want to preallocate, so leave them as kwargs
    dw = Ref{Any}(nothing),
    dx = Ref{Any}(nothing),
    dhx = Ref{Any}(nothing),
    dcx = Ref{Any}(nothing),
)
    # Verify all inputs: they should be compatible with rnnDesc (in case it is supplied), not necessarily with kwargs:
    rd = cudnnGetRNNDescriptor_v8(rnnDesc)
    @assert rd.hiddenSize > 0  "hiddenSize > 0 must be provided"
    @assert cudnnDataType(eltype(x)) == rd.dataType "Input x type not compatible with RNN"
    @assert size(x,1) == rd.inputSize               "Input x size not compatible with RNN"
    ydims = (rd.projSize << (rd.dirMode === CUDNN_BIDIRECTIONAL), size(x)[2:end]...)
    if y !== nothing
        @assert cudnnDataType(eltype(y)) == rd.dataType  "Output y type not compatible with RNN"
        @assert size(y) == ydims                         "Output y size not compatible with RNN or input x"
    else
        y = similar(x, ydims)
    end
    if layout === CUDNN_RNN_DATA_LAYOUT_SEQ_MAJOR_UNPACKED
        if seqLengthArray === nothing; seqLengthArray = fill(Cint(size(x,3)), size(x,2)); end
        @assert size(x,2) == length(seqLengthArray)  "Input x batchsize not compatible with seqLengthArray"
        @assert size(x,3) >= maximum(seqLengthArray) "Input x seqLength not compatible with seqLengthArray"
        xDesc = cudnnRNNDataDescriptor(rd.dataType, layout, Cint(size(x,3)), Cint(size(x,2)), Cint(size(x,1)), seqLengthArray, paddingFill)
        yDesc = cudnnRNNDataDescriptor(rd.dataType, layout, Cint(size(y,3)), Cint(size(y,2)), Cint(size(y,1)), seqLengthArray, paddingFill)
    elseif layout === CUDNN_RNN_DATA_LAYOUT_BATCH_MAJOR_UNPACKED
        if seqLengthArray === nothing; seqLengthArray = fill(Cint(size(x,2)), size(x,3)); end
        @assert size(x,3) == length(seqLengthArray)  "Input x batchsize not compatible with seqLengthArray"
        @assert size(x,2) >= maximum(seqLengthArray) "Input x seqLength not compatible with seqLengthArray"
        xDesc = cudnnRNNDataDescriptor(rd.dataType, layout, Cint(size(x,2)), Cint(size(x,3)), Cint(size(x,1)), seqLengthArray, paddingFill)
        yDesc = cudnnRNNDataDescriptor(rd.dataType, layout, Cint(size(y,2)), Cint(size(y,3)), Cint(size(y,1)), seqLengthArray, paddingFill)
    elseif layout === CUDNN_RNN_DATA_LAYOUT_SEQ_MAJOR_PACKED
        if seqLengthArray === nothing; seqLengthArray = fill(Cint(size(x,3)), size(x,2)); end
        @assert length(x)÷size(x,1) == sum(seqLengthArray) "Input x size not compatible with seqLengthArray"
        xDesc = cudnnRNNDataDescriptor(rd.dataType, layout, maximum(seqLengthArray), Cint(length(seqLengthArray)), Cint(size(x,1)), seqLengthArray, paddingFill)
        yDesc = cudnnRNNDataDescriptor(rd.dataType, layout, maximum(seqLengthArray), Cint(length(seqLengthArray)), Cint(size(y,1)), seqLengthArray, paddingFill)
    else
        error("Unknown layout $layout")
    end
    if devSeqLengths === nothing; devSeqLengths = CuArray(seqLengthArray); end

    hdims = (rd.projSize, length(seqLengthArray), rd.numLayers << (rd.dirMode === CUDNN_BIDIRECTIONAL))
    cdims = (rd.hiddenSize, length(seqLengthArray), rd.numLayers << (rd.dirMode === CUDNN_BIDIRECTIONAL))
    if hx !== nothing
        @assert cudnnDataType(eltype(hx)) == rd.dataType "Hidden hx type not compatible with RNN"
        @assert size(hx) == hdims                        "Hidden hx size not compatible with RNN"
    end
    if hy !== nothing
        @assert hy isa Ref{Any}
        if isassigned(hy) && hy[] !== nothing
            @assert cudnnDataType(eltype(hy[])) == rd.dataType "Hidden hy type not compatible with RNN"
            @assert size(hy[]) == hdims                        "Hidden hy size not compatible with RNN"
        else
            hy[] = similar(y, hdims)
        end
    end
    if rd.cellMode === CUDNN_LSTM
        if cx !== nothing
            @assert cudnnDataType(eltype(cx)) == rd.dataType "Hidden cx type not compatible with RNN"
            @assert size(cx) == cdims                        "Hidden cx size not compatible with RNN"
        end
        if cy !== nothing
            @assert cy isa Ref{Any}
            if isassigned(cy) && cy[] !== nothing
                @assert cudnnDataType(eltype(cy[])) == rd.dataType "Hidden cy type not compatible with RNN"
                @assert size(cy[]) == cdims                        "Hidden cy size not compatible with RNN"
            else
                cy[] = similar(y, cdims)
            end
        end
    end
    hDesc = cudnnTensorDescriptor(CUDNN_TENSOR_NCHW, rd.dataType, Cint(3), Cint[reverse(hdims)...])
    cDesc = cudnnTensorDescriptor(CUDNN_TENSOR_NCHW, rd.dataType, Cint(3), Cint[reverse(cdims)...])

    weightSpaceSize = cudnnRNNWeightSpaceSize(rnnDesc)
    @assert sizeof(w) >= weightSpaceSize "RNN weights should be at least $weightSpaceSize bytes."

    # Backward called separately on each variable. We will calculate all gradients on first call. Use `dready` to avoid subsequent calls.
    dready = Ref{Bool}(false)   # this will be turned to `true` by the first backward call.

    y_h_c = cudnnRNNForwardAD(w, x, hx, cx; rnnDesc, fwdMode, devSeqLengths, xDesc, yDesc, y, hDesc, hy=(hy isa Ref ? hy[] : hy), cDesc, cy=(cy isa Ref ? cy[] : cy), workspace, reserveSpace, dw, dx, dhx, dcx, dready)
    if hy isa Ref; hy[] = y_h_c[2]; end
    if cy isa Ref && rd.cellMode === CUDNN_LSTM; cy[] = y_h_c[3]; end
    return y_h_c[1]             # only return y; hy and cy can be accessed through keyword arguments. They still need to be in AutoGrad return value to be included in gradient calc.
end

# AD method

function cudnnRNNForwardAD(w, x, hx, cx; rnnDesc, fwdMode, devSeqLengths, xDesc, yDesc, y, hDesc, hy, cDesc, cy, workspace, reserveSpace, dw, dx, dhx, dcx, dready)
    (workspaceSize, reserveSpaceSize) = cudnnRNNTempSpaceSizes(rnnDesc, fwdMode, xDesc)
    if reserveSpaceSize > 0 && reserveSpace === nothing; reserveSpace = cudnnTempSpace(reserveSpaceSize); end
    @assert sizeof(reserveSpace) >= reserveSpaceSize  "reserveSpace should be at least $reserveSpaceSize bytes"
    # Cannot use @workspace here because it is shared between forw and back calls
    if workspaceSize > 0 && workspace === nothing; workspace = cudnnTempSpace(workspaceSize); end
    @assert sizeof(workspace) >= workspaceSize  "workspace should be at least $workspaceSize bytes"
    cudnnRNNForward(handle(), rnnDesc, fwdMode, devSeqLengths, xDesc, x, yDesc, y, hDesc, something(hx, CU_NULL), something(hy, CU_NULL), cDesc, something(cx, CU_NULL), something(cy, CU_NULL), sizeof(w), w, sizeof(workspace), something(workspace, CU_NULL), sizeof(reserveSpace), something(reserveSpace, CU_NULL))
    return (y, hy, cy)
end

# Helper methods

function cudnnRNNWeightSpaceSize(rnnDesc::cudnnRNNDescriptor)
    ws = Csize_t[0]
    cudnnGetRNNWeightSpaceSize(handle(), rnnDesc, ws)
    ws[1]
end

function cudnnRNNTempSpaceSizes(rnnDesc::cudnnRNNDescriptor, fwdMode::cudnnForwardMode_t, xDesc::cudnnRNNDataDescriptor)
    ws = Csize_t[0]; rs = Csize_t[0]
    cudnnGetRNNTempSpaceSizes(handle(), rnnDesc, fwdMode, xDesc, ws, rs)
    ws[1], rs[1]
end

function cudnnGetRNNDescriptor_v8(rnnDesc::cudnnRNNDescriptor)
    (algo, cellMode, biasMode, dirMode, inputMode, dataType, mathPrec, mathType, inputSize, hiddenSize, projSize, numLayers, dropout, auxFlags) = (Ref{cudnnRNNAlgo_t}(), Ref{cudnnRNNMode_t}(), Ref{cudnnRNNBiasMode_t}(), Ref{cudnnDirectionMode_t}(), Ref{cudnnRNNInputMode_t}(), Ref{cudnnDataType_t}(), Ref{cudnnDataType_t}(), Ref{cudnnMathType_t}(), Ref{Int32}(), Ref{Int32}(), Ref{Int32}(), Ref{Int32}(), Ref{Ptr{Nothing}}(), Ref{UInt32}())
    cudnnGetRNNDescriptor_v8(rnnDesc, algo, cellMode, biasMode, dirMode, inputMode, dataType, mathPrec, mathType, inputSize, hiddenSize, projSize, numLayers, dropout, auxFlags)
    (algo, cellMode, biasMode, dirMode, inputMode, dataType, mathPrec, mathType, inputSize, hiddenSize, projSize, numLayers, dropout, auxFlags) = (algo[], cellMode[], biasMode[], dirMode[], inputMode[], dataType[], mathPrec[], mathType[], inputSize[], hiddenSize[], projSize[], numLayers[], dropout[], auxFlags[])
    (; rnnDesc, algo, cellMode, biasMode, dirMode, inputMode, dataType, mathPrec, mathType, inputSize, hiddenSize, projSize, numLayers, dropout, auxFlags)
end

checkHidden(h) = (h > 0 ? Int32(h) : error("hiddenSize > 0 is required"))

"""
    cudnnGetRNNWeightParams(w, d::cudnnRNNDescriptor)
    cudnnGetRNNWeightParams(w; hiddenSize, o...)

Return an array of weight matrices and bias vectors of an RNN specified by `d` or keyword
options as views into `w`. The keyword arguments and defaults in the second form are the
same as those in cudnnRNNForward specifying the RNN.

In the returned array `a[1,l,p]` and `a[2,l,p]` give the weight matrix and bias vector for
the l'th layer and p'th parameter or `nothing` if the specified matrix/vector does not
exist. Note that the matrices should be transposed for left multiplication, e.g. `a[1,l,p]'
* x`

The `l` index refers to the pseudo-layer number. In uni-directional RNNs, a pseudo-layer is
the same as a physical layer (pseudoLayer=1 is the RNN input layer, pseudoLayer=2 is the
first hidden layer). In bi-directional RNNs, there are twice as many pseudo-layers in
comparison to physical layers:

    pseudoLayer=1 refers to the forward direction sub-layer of the physical input layer
    pseudoLayer=2 refers to the backward direction sub-layer of the physical input layer
    pseudoLayer=3 is the forward direction sub-layer of the first hidden layer, and so on

The `p` index refers to the weight matrix or bias vector linear ID index.

If cellMode in rnnDesc was set to CUDNN_RNN_RELU or CUDNN_RNN_TANH:

    Value 1 references the weight matrix or bias vector used in conjunction with the input from the previous layer or input to the RNN model.
    Value 2 references the weight matrix or bias vector used in conjunction with the hidden state from the previous time step or the initial hidden state.

If cellMode in rnnDesc was set to CUDNN_LSTM:

    Values 1, 2, 3 and 4 reference weight matrices or bias vectors used in conjunction with the input from the previous layer or input to the RNN model.
    Values 5, 6, 7 and 8 reference weight matrices or bias vectors used in conjunction with the hidden state from the previous time step or the initial hidden state.
    Value 9 corresponds to the projection matrix, if enabled (there is no bias in this operation).

Values and their LSTM gates:

    Values 1 and 5 correspond to the input gate.
    Values 2 and 6 correspond to the forget gate.
    Values 3 and 7 correspond to the new cell state calculations with hyperbolic tangent.
    Values 4 and 8 correspond to the output gate.

If cellMode in rnnDesc was set to CUDNN_GRU:

    Values 1, 2 and 3 reference weight matrices or bias vectors used in conjunction with the input from the previous layer or input to the RNN model.
    Values 4, 5 and 6 reference weight matrices or bias vectors used in conjunction with the hidden state from the previous time step or the initial hidden state.

Values and their GRU gates:

    Values 1 and 4 correspond to the reset gate.
    Values 2 and 5 reference to the update gate.
    Values 3 and 6 correspond to the new hidden state calculations with hyperbolic tangent.

"""
function cudnnGetRNNWeightParams(
    w;
    hiddenSize::Integer,
    inputSize::Integer = hiddenSize,
    projSize::Integer = hiddenSize,
    algo::cudnnRNNAlgo_t = CUDNN_RNN_ALGO_STANDARD,
    cellMode::cudnnRNNMode_t = CUDNN_LSTM,
    biasMode::cudnnRNNBiasMode_t = CUDNN_RNN_DOUBLE_BIAS,
    dirMode::cudnnDirectionMode_t = CUDNN_UNIDIRECTIONAL,
    inputMode::cudnnRNNInputMode_t = CUDNN_LINEAR_INPUT,
    dataType::DataType = Float32,
    mathPrec::DataType = dataType,
    mathType::cudnnMathType_t = math_mode(),
    numLayers::Integer = 1,
    dropout::Real = 0,
    auxFlags::Integer = CUDNN_RNN_PADDED_IO_ENABLED,
)
    cudnnGetRNNWeightParams(w, cudnnRNNDescriptor(algo, cellMode, biasMode, dirMode, inputMode, cudnnDataType(dataType), cudnnDataType(mathPrec), mathType, Int32(inputSize), checkHidden(hiddenSize), Int32(projSize), Int32(numLayers), cudnnDropoutDescriptor(Cfloat(dropout)), UInt32(auxFlags)))
end

function cudnnGetRNNWeightParams(w, rnnDesc::cudnnRNNDescriptor)
    d = cudnnGetRNNDescriptor_v8(rnnDesc)
    T = juliaDataType(d.dataType)
    weightSpace = reinterpret(T, w)
    nlayers = d.numLayers << (d.dirMode === CUDNN_BIDIRECTIONAL)
    nparams = (d.cellMode === CUDNN_RNN_RELU || d.cellMode === CUDNN_RNN_TANH ? 2 :
               d.cellMode === CUDNN_LSTM ? 9 : d.cellMode === CUDNN_GRU  ? 6 :
               error("Unknown cellMode $(d.cellMode)"))
    a = Array{Any}(undef, 2, nlayers, nparams)
    p = Ref{Ptr{Cvoid}}(0)
    cudnnCreateTensorDescriptor(p); mDesc = cudnnTensorDescriptor(p[])
    cudnnCreateTensorDescriptor(p); bDesc = cudnnTensorDescriptor(p[])
    mAddr = Ref{CuPtr{Cvoid}}(0)
    bAddr = Ref{CuPtr{Cvoid}}(0)
    for l in 1:nlayers, p in 1:nparams
        cudnnGetRNNWeightParams(handle(), rnnDesc, l-1, sizeof(weightSpace), weightSpace, p-1, mDesc, mAddr, bDesc, bAddr)
        mT,mD,mS = cudnnGetTensorDescriptor(mDesc)
        bT,bD,bS = cudnnGetTensorDescriptor(bDesc)
        @assert mT === bT === T
        if mAddr[] === CU_NULL
            a[1,l,p] = nothing
        else
            m0 = (mAddr[] - pointer(weightSpace)) ÷ sizeof(T) |> Int
            a[1,l,p] = reshape(view(weightSpace, (m0+1):(m0+prod(mD))), (mD[1],mD[2]))
        end
        if bAddr[] === CU_NULL
            a[2,l,p] = nothing
        else
            b0 = (bAddr[] - pointer(weightSpace)) ÷ sizeof(T) |> Int
            a[2,l,p] = view(weightSpace, (b0+1):(b0+prod(bD)))
        end
    end
    cudnnDestroyTensorDescriptor.((mDesc,bDesc))
    return a
end


text:
function velocity_index(mechanism::Mechanism{T,Nn,Ne}) where {T,Nn,Ne}
    ind = []
    off = 0
	for id in mechanism.root_to_leaves
        (id > Ne) && continue # only treat joints
        joint = mechanism.joints[id]
        nu = input_dimension(joint)
        push!(ind, Vector(off + nu .+ (1:nu)))
        off += 2nu
    end
    return vcat(ind...)
end

# find all the joints parents of a body
function parent_joints(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}, body::Body) where {T,Nn,Ne,Nb,Ni}
	ids = parents(mechanism.system, body.id)
	ids = intersect(ids, 1:Ne) # filter out the bodies
	return [get_node(mechanism, id) for id in ids]
end

# dimensions
"""
    minimal_dimension(mechanism)

    dimension of a mechanism's minimal representation

    mechanism: Mechanism
"""
function minimal_dimension(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}) where {T,Nn,Ne,Nb,Ni}
    nx = 0
    free_rot_base = false # we are going to check if the link attached to the base has free orientation
    nx = 2 * input_dimension(mechanism,
        ignore_floating_base=false)
    free_rot_base && (nx += 1)
    return nx
end

"""
    maximal_dimension(mechanism)

    dimension of a mechanism's maximal representation

    mechanism: Mechanism
"""
maximal_dimension(mechanism::Mechanism{T,Nn,Ne,Nb}; attjac::Bool=false) where {T,Nn,Ne,Nb} = attjac ? 12Nb : 13Nb

"""
    input_dimension(mechanism)

    return the number of inputs for mechanism

    mechanism: Mechanism
"""
function input_dimension(mechanism::Mechanism{T,Nn,Ne,Nb,Ni};
    ignore_floating_base::Bool=false) where {T,Nn,Ne,Nb,Ni}
    nu = 0
    for joint in mechanism.joints
        nu += input_dimension(joint, ignore_floating_base=ignore_floating_base)
    end
    return nu
end

"""
    set_floating_base(mechanism, name)

    returns a mechanism with modified joints having body identified with 'name' as the floating base

    mechanism: Mechanism
    name: Symbol, identifier for floating-base body
"""
function set_floating_base(mechanism::Mechanism, name::Symbol)
    # copy
    mechanism = deepcopy(mechanism)

    # original floating base
    body1 = mechanism.origin

    # new floating base
    body2 = get_body(mechanism, name)

    # find joint for new floating base
    body = body2
    joint = mechanism.joints[findall(j -> j.child_id == body.id, mechanism.joints)[1]]
    joint_list = JointConstraint[]
    push!(joint_list, joint)

    while joint_list[end].parent_id != body1.id
        push!(joint_list, mechanism.joints[findall(j -> j.child_id == joint_list[end].parent_id, mechanism.joints)[1]])
    end

    # reverse JointConstraint (p_old -> c_old) to (c_old -> p_old)
    new_joint_list = JointConstraint[]
    for j in joint_list
        j = deepcopy(j)
        if j.parent_id == body1.id
            j.child_id = body2.id
            push!(new_joint_list, j)
            continue
        else
            parent_id = j.parent_id
            child_id = j.child_id
            pbody = get_body(mechanism, parent_id)
            cbody = get_body(mechanism, child_id)

            # translational conversion
            t = j.translational
            t_reverse = Translational{typeof(t).parameters[1], typeof(t).parameters[2]}(cbody, pbody;
                    parent_vertex=t.vertices[2],
                    child_vertex=t.vertices[1],
                    axis=-t.axis,
                    spring=t.spring,
                    damper=t.damper,
                    spring_offset=t.spring_offset,
                    joint_limits=t.joint_limits,
                    spring_type=t.spring_type)[1]

            r = j.rotational
            r_reverse = Rotational{typeof(r).parameters[1], typeof(r).parameters[2]}(cbody, pbody;
                    axis=-r.axis,
                    orientation_offset=inv(r.orientation_offset),
                    spring=r.spring,
                    damper=r.damper,
                    spring_offset=r.spring_offset,
                    joint_limits=r.joint_limits,
                    spring_type=r.spring_type)[1]

            j.parent_id = child_id
            j.child_id = parent_id
            j.translational = t_reverse
            j.rotational = r_reverse
        end
        push!(new_joint_list, j)
    end

    # set joints
    for j in new_joint_list
        mechanism.joints[j.id] = j
    end

    return Mechanism(mechanism.origin, mechanism.bodies, mechanism.joints, mechanism.contacts,
        gravity=mechanism.gravity,
        timestep=mechanism.timestep)
end

# function reduce_fixed_joints(origin, bodies, joints)
#     remaining_bodies = ones(Bool,length(bodies))
#     remaining_joints = ones(Bool,length(joints))

#     for (j,joint) in enumerate(joints)
#         if typeof(joint) <: JointConstraint{T,6} where T # i.e., Fixed joint
#             parent_body = get_origin_or_body_from_id(joint.parent_id, origin, bodies)
#             child_body = get_origin_or_body_from_id(joint.child_id, origin, bodies)

#             v1, v2 = joint.translational.vertices
#             q_offset = joint.rotational.orientation_offset

#             child_body_com = v1 + vector_rotate(v1,q_offset) # in parent_body's frame
#             if parent_body == origin
#                 new_body_com = zeros(3)
#             else
#                 parent_m, child_m = parent_body.mass, child_body.mass
#                 parent_J, child_J = parent_body.inertia, child_body.inertia

#                 new_body_com = child_body_com*child_m/(parent_m+child_m) # in parent_body's frame

#                 parent_body.mass = parent_m + child_m
#                 new_body_J1 = parent_J + parent_m*skew(new_body_com)'*skew(new_body_com) # in parent_body's frame
#                 new_body_J2 = matrix_rotate(child_J,q_offset) + child_m*skew(-new_body_com)'*skew(-new_body_com) # in parent_body's frame
#                 parent_body.inertia = new_body_J1 + new_body_J2 # in parent_body's frame
#             end

#             parent_body.name = Symbol(String(parent_body.name)*"_merged_with_"*String(child_body.name))

#             for joint2 in joints
#                 joint2 == joint && continue

#                 v12, v22 = joint2.translational.vertices
#                 q_offset2 = joint2.rotational.orientation_offset

#                 if joint2.parent_id == parent_body.id
#                     joint2.translational.vertices = (v12-new_body_com, v22)
#                 elseif joint2.child_id == parent_body.id
#                     joint2.translational.vertices = (v12,v22-new_body_com)
#                 elseif joint2.parent_id == child_body.id
#                     joint2.parent_id = parent_body.id
#                     joint2.translational.vertices = (vector_rotate(v12,q_offset)+child_body_com-new_body_com, v22)
#                     joint2.rotational.orientation_offset = q_offset*q_offset2 # correct?
#                 elseif joint2.child_id == child_body.id
#                     joint2.child_id = parent_body.id
#                     joint2.translational.vertices = (v12,vector_rotate(v22,q_offset)+child_body_com-new_body_com)
#                     joint2.rotational.orientation_offset = q_offset*q_offset2 # correct?
#                 end
#             end
#             remaining_joints[j] = false
#             remaining_bodies[findfirst(x->x==child_body,bodies)] = false
#         end
#     end

#     return origin, bodies[remaining_bodies], joints[remaining_joints]
# end

# function get_origin_or_body_from_id(id, origin, bodies)
#     origin.id == id && (return origin)
#     for body in bodies
#         body.id == id && (return body)
#     end
# end


text:
# Autogenerated wrapper script for SassC_jll for x86_64-unknown-freebsd
export sassc

using libsass_jll
JLLWrappers.@generate_wrapper_header("SassC")
JLLWrappers.@declare_executable_product(sassc)
function __init__()
    JLLWrappers.@generate_init_header(libsass_jll)
    JLLWrappers.@init_executable_product(
        sassc,
        "bin/sassc",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Foldable

using Funs

export freduce

freduce{Z,FT}(f::Callable, zero::Z, xs::FT, yss...; R::Type=eltype(f)) =
    error("freduce not specialized for types ($Z,$FT)")

# Array
@generated function freduce{Z,T,D}(f::Callable, zero::Z, xs::Array{T,D},
        yss::Array...; R::Type=eltype(f))
    quote
        tuple($([:(@assert length(yss[$n]) == length(xs))
            for n in 1:length(yss)]...))
        r = zero
        @inbounds #=@simd=# for i in eachindex(xs)
            r = f(r, xs[i], $([:(yss[$n][i]) for n in 1:length(yss)]...))
        end
        r::R
    end
end

# Set
function freduce{Z,T}(f::Callable, zero::Z, xs::Set{T}; R::Type=eltype(f))
    r = zero
    for x in xs
        r = f(r, x)
    end
    r::R
end

# Tuple
@generated function freduce{Z}(f::Callable, zero::Z, xs::Tuple, yss::Tuple...;
        R::Type=eltype(f))
    quote
        tuple($([:(@assert length(yss[$n]) == length(xs))
            for n in 1:length(yss)]...))
        r = zero
        @inbounds #=@simd=# for i in eachindex(xs)
            r = f(r, xs[i], $([:(yss[$n][i]) for n in 1:length(yss)]...))
        end
        r::R
    end
end

end


text:
using Diversity.API

import Base.IteratorSize, Base.length, Base.IteratorEltype, Base.eltype

abstract type AbstractIterator{M <: AbstractMetacommunity} end

struct TypeIterator{M} <: AbstractIterator{M}
    metacommunity::M
    viewfn::Function
end

function TypeIterator(fn::Function, meta::M) where M <: AbstractMetacommunity
    n = ndims(fn(meta))
    if n == 2
        viewfn = (meta, i) -> view(fn(meta), i, :)
    elseif n == 1
        viewfn = (meta, i) -> view(fn(meta), i)
    else
        error("Can't iterate over types for the function of the metacommunity - $fn")
    end
    return TypeIterator(meta, viewfn)
end

function TypeIterator(meta::M) where M <: AbstractMetacommunity
    return TypeIterator(getabundance, meta)
end

import Base.iterate
function iterate(ti::TypeIterator, state = 1)
    if state > counttypes(ti.metacommunity)
        return nothing
    end
    return ti.viewfn(ti.metacommunity, state), state + 1
end

function IteratorSize(ti::Type{TypeIterator})
    return HasLength()
end

function length(ti::TypeIterator)
    return counttypes(ti.metacommunity)
end

function IteratorEltype(ti::Type{TypeIterator})
    return HasEltype()
end

function eltype(ti::TypeIterator)
    return eltype(ti.viewfn(ti.metacommunity, 1))
end

struct SubcommunityIterator{M} <: AbstractIterator{M}
    metacommunity::M
    viewfn::Function
end

function SubcommunityIterator(fn::Function,
                              meta::M) where M <: AbstractMetacommunity
    n = ndims(fn(meta))
    if n == 2
        viewfn = (meta, i) -> view(fn(meta), :, i)
    elseif n == 1
        viewfn = (meta, i) -> view(fn(meta), i)
    else
        error("Can't iterate over subcommunities for the " *
              "function of the metacommunity - $fn")
    end
    return SubcommunityIterator(meta, viewfn)
end

function SubcommunityIterator(meta::M) where M <: AbstractMetacommunity
    return SubcommunityIterator(getabundance, meta)
end

function iterate(si::SubcommunityIterator, state = 1)
    if state > countsubcommunities(si.metacommunity)
        return nothing
    end
    return si.viewfn(si.metacommunity, state), state + 1
end

function IteratorSize(si::Type{SubcommunityIterator})
    return HasLength()
end

function length(si::SubcommunityIterator)
    return countsubcommunities(si.metacommunity)
end

function IteratorEltype(si::Type{SubcommunityIterator})
    return HasEltype()
end

function eltype(si::SubcommunityIterator)
    return eltype(si.viewfn(si.metacommunity, 1))
end


text:
struct CSPData <: BaseData
    name::String
    capacity::Int64
    weights::Vector{Int64}
    demands::Vector{Int64}

    lb::Int64
    ub::Int64
end

function Base.show(io::IO, data::CSPData)
    print(io, "CSP Data $(data.name)")
    print(io, " ($(length(data.weights)) elements,")
    print(io, " $(sum(data.demands)) demands,")
    print(io, " capacity = $(data.capacity))")
    if data.ub != typemax(Int64)
        print(io, " [$(data.lb), $(data.ub)]")
    else
        print(io, " [$(data.lb), Inf]")
    end
end


text:
using Aqua
using PersistenceDiagramsBase

Aqua.test_all(PersistenceDiagramsBase)


text:
module SpringCollab2020SwapDashTrigger

using ..Ahorn, Maple

@mapdef Trigger "SpringCollab2020/SwapDashTrigger" SwapDashTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight)

const placements = Ahorn.PlacementDict(
    "Swap Dash (Spring Collab 2020)" => Ahorn.EntityPlacement(
        SwapDashTrigger,
        "rectangle"
    )
)

end


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

name = "Giflib"
version = v"5.1.4"

# Collection of sources required to build Giflib
sources = [
    ArchiveSource("https://downloads.sourceforge.net/project/giflib/giflib-$(version).tar.bz2",
                  "df27ec3ff24671f80b29e6ab1c4971059c14ac3db95406884fc26574631ba8d5"),

    # Apply security patch
    ArchiveSource("https://deb.debian.org/debian/pool/main/g/giflib/giflib_5.1.4-3.debian.tar.xz",
                  "767ea03c1948fa203626107ead3d8b08687a3478d6fbe4690986d545fb1d60bf"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/giflib-*/
atomic_patch -p1 ../debian/patches/CVE-2016-3977.patch

# We need to massage configure script to convince it to build the shared library
# for PowerPC.
if [[ "${target}" == powerpc64le-* ]]; then
    autoreconf -vi
fi

update_configure_scripts
./configure --prefix=${prefix} --host=${target}
make -j${nproc}
make install
install_license COPYING
"""

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

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

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

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

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# Wisdom & Hernandez version of Kepler solver, with Rein & Tamayo convergence test.
# Now using \gamma = \sqrt{\abs{\beta}}s rather than s now to solve Kepler's equation.

using ForwardDiff, DiffResults

include("g3.jl")

function cubic1(a::T, b::T, c::T) where {T <: Real}
a3 = a*third
Q = a3^2 - b*third
R = a3^3 + 0.5*(-a3*b + c)
if R^2 < Q^3
# println("Error in cubic solver ",R^2," ",Q^3)
 return -c/b
else
 A = -sign(R)*cbrt(abs(R) + sqrt(R*R - Q*Q*Q))
 if A == 0.0
   B = 0.0
  else
   B = Q/A
 end
 x1 = A + B - a3
 return x1
end
return
end

function jac_delxv_gamma!(x0::Array{T,1},v0::Array{T,1},k::T,h::T,drift_first::Bool;grad::Bool=false,auto::Bool=false,dlnq::T=convert(T,0.0),debug=false) where {T <: Real}
# Using autodiff, computes Jacobian of delx & delv with respect to x0, v0, k & h.

# Autodiff requires a single-vector input, so create an array to hold the independent variables:
  input = zeros(T,8)
  input[1:3]=x0; input[4:6]=v0; input[7]=k; input[8]=h
  if grad
    if debug
      # Also output gamma, r, fm1, dfdt, gmh, dgdtm1, and for debugging:
      delxv_jac = zeros(T,12,8)
    else
      # Output \delta x & \delta v only:
      delxv_jac = zeros(T,6,8)
    end
  end

# Create a closure so that the function knows value of drift_first:

  function delx_delv(input::Array{T2,1}) where {T2 <: Real} # input = x0,v0,k,h,drift_first
  # Compute delx and delv from h, s, k, beta0, x0 and v0:
  x0 = input[1:3]; v0 = input[4:6]; k = input[7]; h = input[8]
  # Compute r0:
  drift_first ?  r0 = norm(x0-h*v0) : r0 = norm(x0)
  # And its inverse:
  r0inv = inv(r0)
  # Compute beta_0:
  beta0 = 2k*r0inv-dot(v0,v0)
  beta0inv = inv(beta0)
  signb = sign(beta0)
  sqb = sqrt(signb*beta0)
  zeta = k-r0*beta0
  gamma_guess = zero(T2)
  # Compute \eta_0 = x_0 . v_0:
  drift_first ?  eta = dot(x0-h*v0,v0) : eta = dot(x0,v0)
  if zeta != zero(T2)
    # Make sure we have a cubic in gamma (and don't divide by zero):
    gamma_guess = cubic1(3eta*sqb/zeta,6r0*signb*beta0/zeta,-6h*signb*beta0*sqb/zeta)
  else
    # Check that we have a quadratic in gamma (and don't divide by zero):
    if eta != zero(T2)
      reta = r0/eta
      disc = reta^2+2h/eta
      disc > zero(T2) ?  gamma_guess = sqb*(-reta+sqrt(disc)) : gamma_guess = h*r0inv*sqb
    else
      gamma_guess = h*r0inv*sqb
    end
  end
  gamma  = copy(gamma_guess)
  # Make sure prior two steps differ:
  gamma1 = 2*copy(gamma)
  gamma2 = 3*copy(gamma)
  iter = 0
  ITMAX = 20
  # Compute coefficients: (8/28/19 notes)
#  c1 = k; c2 = -zeta; c3 = -eta*sqb; c4 = sqb*(eta-h*beta0); c5 = eta*signb*sqb
  c1 = k; c2 = -2zeta; c3 = 2*eta*signb*sqb; c4 = -sqb*h*beta0; c5 = 2eta*signb*sqb
  # Solve for gamma:
  while true
    gamma2 = gamma1
    gamma1 = gamma
    xx = 0.5*gamma
#    xx = gamma
    if beta0 > 0
      sx = sin(xx); cx = cos(xx)
    else
      sx = sinh(xx); cx = exp(-xx)+sx
    end
#    gamma -= (c1*gamma+c2*sx+c3*cx+c4)/(c2*cx+c5*sx+c1)
    gamma -= (k*gamma+c2*sx*cx+c3*sx^2+c4)/(2signb*zeta*sx^2+c5*sx*cx+r0*beta0)
    iter +=1
    if iter >= ITMAX || gamma == gamma2 || gamma == gamma1
      break
    end
  end
#  if typeof(gamma) ==  Float64
#    println("s: ",gamma/sqb)
#  end
  # Set up a single output array for delx and delv:
  if debug
    delxv = zeros(T2,12)
  else
    delxv = zeros(T2,6)
  end
  # Since we updated gamma, need to recompute:
  xx = 0.5*gamma
  if beta0 > 0
    sx = sin(xx); cx = cos(xx)
  else
    sx = sinh(xx); cx = exp(-xx)+sx
  end
  # Now, compute final values.  Compute Wisdom/Hernandez G_i^\beta(s) functions:
  g1bs = 2sx*cx/sqb
  g2bs = 2signb*sx^2*beta0inv
  g0bs = one(T2)-beta0*g2bs
  g3bs = G3(gamma,beta0)
  h1 = zero(T2); h2 = zero(T2)
#  if typeof(g1bs) == Float64
#    println("g1: ",g1bs," g2: ",g2bs," g3: ",g3bs)
#  end
  # Compute r from equation (35):
  r = r0*g0bs+eta*g1bs+k*g2bs
#  if typeof(r) == Float64
#    println("r: ",r)
#  end
  rinv = inv(r)
  dfdt = -k*g1bs*rinv*r0inv # [x]
  if drift_first
    # Drift backwards before Kepler step: (1/22/2018)
    fm1 = -k*r0inv*g2bs # [x]
    # This is given in 2/7/2018 notes: g-h*f
#    gmh = k*r0inv*(r0*(g1bs*g2bs-g3bs)+eta*g2bs^2+k*g3bs*g2bs)  # [x]
#    println("Old gmh: ",gmh," new gmh: ",k*r0inv*(h*g2bs-r0*g3bs))  # [x]
    gmh = k*r0inv*(h*g2bs-r0*g3bs)  # [x]
  else
    # Drift backwards after Kepler step: (1/24/2018)
    # The following line is f-1-h fdot:
    h1= H1(gamma,beta0); h2= H2(gamma,beta0)
#    fm1 =  k*rinv*(g2bs-k*r0inv*H1(gamma,beta0))  # [x]
    fm1 =  k*rinv*(g2bs-k*r0inv*h1)  # [x]
    # This is g-h*dgdt
#    gmh = k*rinv*(r0*H2(gamma,beta0)+eta*H1(gamma,beta0)) # [x]
    gmh = k*rinv*(r0*h2+eta*h1) # [x]
  end
  # Compute velocity component functions:
  if drift_first
    # This is gdot - h fdot - 1:
#    dgdtm1 = k*r0inv*rinv*(r0*g0bs*g2bs+eta*g1bs*g2bs+k*g1bs*g3bs) # [x]
#    println("gdot-h fdot-1: ",dgdtm1," alternative expression: ",k*r0inv*rinv*(h*g1bs-r0*g2bs))
    dgdtm1 = k*r0inv*rinv*(h*g1bs-r0*g2bs)
  else
    # This is gdot - 1:
    dgdtm1 = -k*rinv*g2bs # [x]
  end
#  if typeof(fm1) == Float64
#    println("fm1: ",fm1," dfdt: ",dfdt," gmh: ",gmh," dgdt-1: ",dgdtm1)
#  end
  @inbounds for j=1:3
  # Compute difference vectors (finish - start) of step:
    delxv[  j] = fm1*x0[j]+gmh*v0[j]        # position x_ij(t+h)-x_ij(t) - h*v_ij(t) or -h*v_ij(t+h)
  end
  @inbounds for j=1:3
    delxv[3+j] = dfdt*x0[j]+dgdtm1*v0[j]    # velocity v_ij(t+h)-v_ij(t)
  end
  if debug
    delxv[7] = gamma
    delxv[8] = r
    delxv[9] = fm1
    delxv[10] = dfdt
    delxv[11] = gmh
    delxv[12] = dgdtm1
  end
  if grad == true && auto == false && dlnq == 0.0
    # Compute gradient analytically:
    jac_mass = zeros(T,6)
    compute_jacobian_gamma!(gamma,g0bs,g1bs,g2bs,g3bs,h1,h2,dfdt,fm1,gmh,dgdtm1,r0,r,r0inv,rinv,k,h,beta0,beta0inv,eta,sqb,zeta,x0,v0,delxv_jac,jac_mass,drift_first,debug)
  end
  return delxv
  end

# Use autodiff to compute Jacobian:
if grad
  if auto
#    delxv_jac = ForwardDiff.jacobian(delx_delv,input)
    if debug
      delxv = zeros(T,12)
    else
      delxv = zeros(T,6)
    end
    out = DiffResults.JacobianResult(delxv,input)
    ForwardDiff.jacobian!(out,delx_delv,input)
    delxv_jac = DiffResults.jacobian(out)
    delxv = DiffResults.value(out)
  elseif dlnq != 0.0
# Use finite differences to compute Jacobian:
    if debug
      delxv_jac = zeros(T,12,8)
    else
      delxv_jac = zeros(T,6,8)
    end
    delxv = delx_delv(input)
    @inbounds for j=1:8
      # Difference the jth component:
      inputp = copy(input); dp = dlnq*inputp[j]; inputp[j] += dp
      delxvp = delx_delv(inputp)
      inputm = copy(input); inputm[j] -= dp
      delxvm = delx_delv(inputm)
      delxv_jac[:,j] = (delxvp-delxvm)/(2*dp)
    end
  else
# If grad = true and dlnq = 0.0, then the above routine will compute Jacobian analytically.
    delxv = delx_delv(input)
  end
# Return Jacobian:
  return  delxv::Array{T,1},delxv_jac::Array{T,2}
else
  return delx_delv(input)::Array{T,1}
end
end

function jac_delxv_gamma!(x0::Array{T,1},v0::Array{T,1},k::T,h::T,drift_first::Bool,delxv::Array{T,1},delxv_jac::Array{T,2},jac_mass::Array{T,1},debug::Bool) where {T <: Real}
# Analytically computes Jacobian of delx & delv with respect to x0, v0, k & h.

  # Compute r0:
  drift_first ?  r0 = norm(x0-h*v0) : r0 = norm(x0)
  # And its inverse:
  r0inv = inv(r0)
  # Compute beta_0:
  beta0 = 2k*r0inv-dot(v0,v0)
  beta0inv = inv(beta0)
  signb = sign(beta0)
  sqb = sqrt(signb*beta0)
  zeta = k-r0*beta0
  gamma_guess = zero(T)
  # Compute \eta_0 = x_0 . v_0:
  drift_first ?  eta = dot(x0-h*v0,v0) : eta = dot(x0,v0)
  if zeta != zero(T)
    # Make sure we have a cubic in gamma (and don't divide by zero):
    gamma_guess = cubic1(3eta*sqb/zeta,6r0*signb*beta0/zeta,-6h*signb*beta0*sqb/zeta)
  else
    # Check that we have a quadratic in gamma (and don't divide by zero):
    if eta != zero(T)
      reta = r0/eta
      disc = reta^2+2h/eta
      disc > zero(T) ?  gamma_guess = sqb*(-reta+sqrt(disc)) : gamma_guess = h*r0inv*sqb
    else
      gamma_guess = h*r0inv*sqb
    end
  end
  gamma  = copy(gamma_guess)
  # Make sure prior two steps differ:
  gamma1 = 2*copy(gamma)
  gamma2 = 3*copy(gamma)
  iter = 0
  ITMAX = 20
  # Compute coefficients: (8/28/19 notes)
#  c1 = k; c2 = -zeta; c3 = -eta*sqb; c4 = sqb*(eta-h*beta0); c5 = eta*signb*sqb
  c1 = k; c2 = -2zeta; c3 = 2*eta*signb*sqb; c4 = -sqb*h*beta0; c5 = 2eta*signb*sqb
  # Solve for gamma:
  while true
    gamma2 = gamma1
    gamma1 = gamma
    xx = 0.5*gamma
#    xx = gamma
    if beta0 > 0
      sx = sin(xx); cx = cos(xx)
    else
      sx = sinh(xx); cx = exp(-xx)+sx
    end
#    gamma -= (c1*gamma+c2*sx+c3*cx+c4)/(c2*cx+c5*sx+c1)
    gamma -= (k*gamma+c2*sx*cx+c3*sx^2+c4)/(2signb*zeta*sx^2+c5*sx*cx+r0*beta0)
    iter +=1
    if iter >= ITMAX || gamma == gamma2 || gamma == gamma1
      break
    end
  end
  # Set up a single output array for delx and delv:
  fill!(delxv,zero(T))
  # Since we updated gamma, need to recompute:
  xx = 0.5*gamma
  if beta0 > 0
    sx = sin(xx); cx = cos(xx)
  else
    sx = sinh(xx); cx = exp(-xx)+sx
  end
  # Now, compute final values.  Compute Wisdom/Hernandez G_i^\beta(s) functions:
  g1bs = 2sx*cx/sqb
  g2bs = 2signb*sx^2*beta0inv
  g0bs = one(T)-beta0*g2bs
  g3bs = G3(gamma,beta0)
  h1 = zero(T); h2 = zero(T)
  # Compute r from equation (35):
  r = r0*g0bs+eta*g1bs+k*g2bs
  rinv = inv(r)
  dfdt = -k*g1bs*rinv*r0inv # [x]
  if drift_first
    # Drift backwards before Kepler step: (1/22/2018)
    fm1 = -k*r0inv*g2bs # [x]
    # This is given in 2/7/2018 notes: g-h*f
#    gmh = k*r0inv*(r0*(g1bs*g2bs-g3bs)+eta*g2bs^2+k*g3bs*g2bs)  # [x]
    gmh = k*r0inv*(h*g2bs-r0*g3bs)  # [x]
  else
    # Drift backwards after Kepler step: (1/24/2018)
    # The following line is f-1-h fdot:
    h1= H1(gamma,beta0); h2= H2(gamma,beta0)
#    fm1 =  k*rinv*(g2bs-k*r0inv*H1(gamma,beta0))  # [x]
    fm1 =  k*rinv*(g2bs-k*r0inv*h1)  # [x]
    # This is g-h*dgdt
#    gmh = k*rinv*(r0*H2(gamma,beta0)+eta*H1(gamma,beta0)) # [x]
    gmh = k*rinv*(r0*h2+eta*h1) # [x]
  end
  # Compute velocity component functions:
  if drift_first
    # This is gdot - h fdot - 1:
#    dgdtm1 = k*r0inv*rinv*(r0*g0bs*g2bs+eta*g1bs*g2bs+k*g1bs*g3bs) # [x]
    dgdtm1 = k*r0inv*rinv*(h*g1bs-r0*g2bs)
  else
    # This is gdot - 1:
    dgdtm1 = -k*rinv*g2bs # [x]
  end
  @inbounds for j=1:3
  # Compute difference vectors (finish - start) of step:
    delxv[  j] = fm1*x0[j]+gmh*v0[j]        # position x_ij(t+h)-x_ij(t) - h*v_ij(t) or -h*v_ij(t+h)
  end
  @inbounds for j=1:3
    delxv[3+j] = dfdt*x0[j]+dgdtm1*v0[j]    # velocity v_ij(t+h)-v_ij(t)
  end
  if debug
    delxv[7] = gamma
    delxv[8] = r
    delxv[9] = fm1
    delxv[10] = dfdt
    delxv[11] = gmh
    delxv[12] = dgdtm1
  end
  # Compute gradient analytically:
  compute_jacobian_gamma!(gamma,g0bs,g1bs,g2bs,g3bs,h1,h2,dfdt,fm1,gmh,dgdtm1,r0,r,r0inv,rinv,k,h,beta0,beta0inv,eta,sqb,zeta,x0,v0,delxv_jac,jac_mass,drift_first,debug)
end

function compute_jacobian_gamma!(gamma::T,g0::T,g1::T,g2::T,g3::T,h1::T,h2::T,dfdt::T,fm1::T,gmh::T,dgdtm1::T,r0::T,r::T,r0inv::T,rinv::T,k::T,h::T,beta::T,betainv::T,eta::T,
  sqb::T,zeta::T,x0::Array{T,1},v0::Array{T,1},delxv_jac::Array{T,2},jac_mass::Array{T,1},drift_first::Bool,debug::Bool) where {T <: Real}
# Computes Jacobian:
if drift_first
  # First, x0 derivatives:
  #  delxv[  j] = fm1*x0[j]+gmh*v0[j]        # position x_ij(t+h)-x_ij(t) - h*v_ij(t) or -h*v_ij(t+h)
  #  delxv[3+j] = dfdt*x0[j]+dgdtm1*v0[j]    # velocity v_ij(t+h)-v_ij(t)
  # First, the diagonal terms:
  # Now the off-diagonal terms:
  d   = (h + eta*g2 + 2*k*g3)*betainv
  c1 = d-r0*g3
  c2 = eta*g0+g1*zeta
  c3  = d*k+g1*r0^2
  c4 = eta*g1+2*g0*r0
  c13 = g1*h-g2*r0
#  c9  = g2*r-h1*k
#  c10 = c3*c9*rinv+g2*r0*h-k*(2*g2*h+3*g3*r0)*betainv
  c9 = 2*g2*h-3*g3*r0
  c10 = k*r0inv^4*(-g2*r0*h+k*c9*betainv-c3*c13*rinv)
  c24 = r0inv^3*(r0*(2*k*r0inv-beta)*betainv-g1*c3*rinv/g2)
  h6 = H6(gamma,beta)
  # Derivatives of \delta x with respect to x0, v0, k & h:
  dfm1dxx = fm1*c24
  dfm1dxv = -fm1*(g1*rinv+h*c24)
  dfm1dvx = dfm1dxv
#  dfm1dvv = fm1*(2*betainv-g1*rinv*(d/g2-2*h)+h^2*c24)
  dfm1dvv = fm1*rinv*(-r0*g2 + k*h6*betainv/g2 + h*(2*g1+h*r*c24))
  dfm1dh  = fm1*(g1*rinv*(1/g2+2*k*r0inv-beta)-eta*c24)
  dfm1dk  = fm1*(1/k+g1*c1*rinv*r0inv/g2-2*betainv*r0inv)
#  dfm1dk2  = 2g2*betainv-c1*g1*rinv
  h4 = -H1(gamma,beta)*beta
#  h5 = g1*g2-g3*(2+g0)
  h5 = H5(gamma,beta)
#  println("H5 : ",g1*g2-g3*(2+g0)," ",H2(gamma,beta)-2*G3(gamma,beta)," ",h5)
#  h6 = 2*g2^2-3*g1*g3
#  println("H6: ",2*g2^2-3*g1*g3," ",h6)
#  dfm1dk2  = (r0*h4+k*h6)*betainv*rinv
  dfm1dk2  = (r0*h4+k*h6)
#  println("dfm1dk2: ",dfm1dk2," ", (r0*h4+k*h6)*betainv*rinv)
  dgmhdxx = c10
  dgmhdxv =  -g2*k*c13*rinv*r0inv-h*c10
  dgmhdvx =  dgmhdxv
  # dgmhdvv =  -d*k*c13*rinv*r0inv+2*g2*h*k*c13*rinv*r0inv+k*c9*betainv*r0inv+h^2*c10
  h8 = H8(gamma,beta)
  dgmhdvv =  2*g2*h*k*c13*rinv*r0inv+h^2*c10+
              k*betainv*rinv*r0inv*(r0^2*h8-beta*h*r0*g2^2 + (h*k+eta*r0)*h6)
  dgmhdh  =  g2*k*r0inv+k*c13*rinv*r0inv+g2*k*(2*k*r0inv-beta)*c13*rinv*r0inv-eta*c10
  dgmhdk  =  r0inv*(k*c1*c13*rinv*r0inv+g2*h-g3*r0-k*c9*betainv*r0inv)
#  dgmhdk2  =  c1*c13*rinv-c9*betainv
#  dgmhdk2 = -betainv*rinv*(h6*g3*k^2+eta*r0*(h6+g2*h4)+r0^2*g0*h5+k*eta*g2*h6+(g1*h6+g3*h4)*k*r0)
  dgmhdk2 = -(h6*g3*k^2+eta*r0*(h6+g2*h4)+r0^2*g0*h5+k*eta*g2*h6+(g1*h6+g3*h4)*k*r0)
#  println("dgmhdk2: ",dgmhdk2," ",-betainv*rinv*(h6*g3*k^2+eta*r0*(h6+g2*h4)+r0^2*g0*h5+k*eta*g2*h6+(g1*h6+g3*h4)*k*r0))
  @inbounds for j=1:3
    # First, compute the \delta x-derivatives:
    delxv_jac[  j,  j] = fm1
    delxv_jac[  j,3+j] = gmh
    @inbounds for i=1:3
      delxv_jac[j,  i] += (dfm1dxx*x0[i]+dfm1dxv*v0[i])*x0[j] + (dgmhdxx*x0[i]+dgmhdxv*v0[i])*v0[j]
      delxv_jac[j,3+i] += (dfm1dvx*x0[i]+dfm1dvv*v0[i])*x0[j] + (dgmhdvx*x0[i]+dgmhdvv*v0[i])*v0[j]
    end
    delxv_jac[  j,  7] = dfm1dk*x0[j] + dgmhdk*v0[j]
    delxv_jac[  j,  8] = dfm1dh*x0[j] + dgmhdh*v0[j]
    # Compute the mass jacobian separately since otherwise cancellations happen in kepler_driftij_gamma:
    jac_mass[  j] = (GNEWT*r0inv)^2*betainv*rinv*(dfm1dk2*x0[j]+dgmhdk2*v0[j])
  end
  # Derivatives of \delta v with respect to x0, v0, k & h:
  c5 = (r0-k*g2)*rinv/g1
  c6 = (r0*g0-k*g2)*betainv
  c7 = g2*(1/g1+c2*rinv)
  c8 = (k*c6+r*r0+c3*c5)*r0inv^3
  c12 = g0*h-g1*r0
  c17 = r0-r-g2*k
  c18 = eta*g1+2*g2*k
  c20 = k*(g2*k+r)-g0*r0*zeta
  c21 = (g2*k-r0)*(beta*c3-k*g1*r)*betainv*rinv^2*r0inv^3/g1+eta*g1*rinv*r0inv^2-2r0inv^2
  c22 = rinv*(-g1-g0*g2/g1+g2*c2*rinv)
  c25 = k*rinv*r0inv^2*(-g2+k*(c13-g2*r0)*betainv*r0inv^2-c13*r0inv-c12*c3*rinv*r0inv^2+
                      c13*c2*c3*rinv^2*r0inv^2-c13*(k*(g2*k+r)-g0*r0*zeta)*betainv*rinv*r0inv^2)
  c26 = k*rinv^2*r0inv*(-g2*c12-g1*c13+g2*c13*c2*rinv)
  ddfdtdxx = dfdt*c21
  ddfdtdxv = dfdt*(c22-h*c21)
  ddfdtdvx = ddfdtdxv
#  ddfdtdvv = dfdt*(betainv*(1-c18*rinv)+d*(g1*c2-g0*r)*rinv^2/g1-h*(2c22-h*c21))
  c34 = (-beta*eta^2*g2^2-eta*k*h8-h6*k^2-2beta*eta*r0*g1*g2+(g2^2-3*g1*g3)*beta*k*r0
        - beta*g1^2*r0^2)*betainv*rinv^2+(eta*g2^2)*rinv/g1 + (k*h8)*betainv*rinv/g1
  ddfdtdvv = dfdt*(c34 - 2*h*c22 +h^2*c21)
  ddfdtdk  = dfdt*(1/k-betainv*r0inv-c17*betainv*rinv*r0inv-c1*(g1*c2-g0*r)*rinv^2*r0inv/g1)
#  ddfdtdk2  = -g1*(-betainv*r0inv-c17*betainv*rinv*r0inv-c1*(g1*c2-g0*r)*rinv^2*r0inv/g1)
  #ddfdtdk2  = -(g2*k-r0)*(g1*r-beta*c1)*betainv*rinv^2*r0inv
  ddfdtdk2 = -(g2*k-r0)*(beta*r0*(g3-g1*g2)-beta*eta*g2^2+k*H3(gamma,beta))*betainv*rinv^2*r0inv
  ddfdtdh  = dfdt*(g0*rinv/g1-c2*rinv^2-(2*k*r0inv-beta)*c22-eta*c21)
  dgdtmhdfdtm1dxx = c25
  dgdtmhdfdtm1dxv = c26-h*c25
  dgdtmhdfdtm1dvx = c26-h*c25
  h2 = H2(gamma,beta)
#  dgdtmhdfdtm1dvv = d*k*rinv^3*r0inv*(c13*c2-r*c12)+k*(c13*(r0*g0-k*g2)-g2*r*r0)*betainv*rinv^2*r0inv-2*h*c26+h^2*c25
#  dgdtmhdfdtm1dvv = d*k*rinv^3*r0inv*k*(r0*(g1*g2-g3)+eta*g2^2+k*g2*g3)+k*(c13*(r0*g0-k*g2)-g2*r*r0)*betainv*rinv^2*r0inv-2*h*c26+h^2*c25
#  println("\dot g - h \dot f -1, dv0 terms: ",d*k*rinv^3*r0inv*c13*c2," ",-d*k*rinv^3*r0inv*r*c12," ",k*c13*r0*g0*betainv*rinv^2*r0inv," ",-k*c13*k*g2*betainv*rinv^2*r0inv," ",-k*g2*r*r0*betainv*rinv^2*r0inv," ",-2*h*c26," ",h^2*c25)
#  println("second version of dv0 terms:     ",d*k*rinv^3*r0inv*k*r0*g1*g2," ",-d*k*rinv^3*r0inv*k*r0*g3," ",d*k*rinv^3*r0inv*k*eta*g2^2," ",d*k*rinv^3*r0inv*k*k*g2*g3," ",-k*eta*k*g1*g2^2*betainv*rinv^2*r0inv," ",-k*g1*g2*g3*k^2*betainv*rinv^2*r0inv," ",-k*r0*eta*beta*g1*g2^2*betainv*rinv^2*r0inv," ",-r0*k*k*g1*h2*betainv*rinv^2*r0inv," ",-beta*k*g2^2*g0*r0^2*betainv*rinv^2*r0inv," ",-2*h*c26," ",h^2*c25)
  c33 = d*k*rinv^3*r0inv*k*(h*g2- r0*g3)+k*(-eta*k*g1*g2^2-g1*g2*g3*k^2-r0*eta*beta*g1*g2^2-r0*k*g1*h2 - beta*g2^2*g0*r0^2)*betainv*rinv^2*r0inv
  dgdtmhdfdtm1dvv = c33-2*h*c26+h^2*c25
  dgdtmhdfdtm1dk = rinv*r0inv*(-k*(c13-g2*r0)*betainv*r0inv+c13-k*c13*c17*betainv*rinv*r0inv+k*c1*c12*rinv*r0inv-k*c1*c2*c13*rinv^2*r0inv)
#  dgdtmhdfdtm1dk2 = -(c13-g2*r0)*betainv*r0inv-c13*c17*betainv*rinv*r0inv+c1*c12*rinv*r0inv-c1*c2*c13*rinv^2*r0inv
  #dgdtmhdfdtm1dk2 = g2*betainv+rinv*r0inv*(c1*c12+c13*((k*g2-r0)*betainv-c1*c2*rinv))
  h3 = H3(gamma,beta)
  dgdtmhdfdtm1dk2 = k*betainv*rinv^2*r0inv*(-beta*eta^2*g2^4+eta*g2*(g1*g2^2+g1^2*g3-5*g2*g3)*k+g2*g3*h3*k^2+
        2eta*r0*beta*g2^2*(g3-g1*g2)+(4g3-g0*g3-g1*g2)*(g3-g1*g2)*r0*k+beta*(2g1*g3*g2-g1^2*g2^2-g3^2)*r0^2)
  dgdtmhdfdtm1dh = g1*k*rinv*r0inv+k*c12*rinv^2*r0inv-k*c2*c13*rinv^3*r0inv-(2*k*r0inv-beta)*c26-eta*c25
  @inbounds for j=1:3
    # Next, compute the \delta v-derivatives:
    delxv_jac[3+j,  j] = dfdt
    delxv_jac[3+j,3+j] = dgdtm1
    @inbounds for i=1:3
      delxv_jac[3+j,  i] += (ddfdtdxx*x0[i]+ddfdtdxv*v0[i])*x0[j] + (dgdtmhdfdtm1dxx*x0[i]+dgdtmhdfdtm1dxv*v0[i])*v0[j]
      delxv_jac[3+j,3+i] += (ddfdtdvx*x0[i]+ddfdtdvv*v0[i])*x0[j] + (dgdtmhdfdtm1dvx*x0[i]+dgdtmhdfdtm1dvv*v0[i])*v0[j]
    end
    delxv_jac[3+j,  7] = ddfdtdk*x0[j] + dgdtmhdfdtm1dk*v0[j]
    delxv_jac[3+j,  8] = ddfdtdh*x0[j] + dgdtmhdfdtm1dh*v0[j]
    # Compute the mass jacobian separately since otherwise cancellations happen in kepler_driftij_gamma:
    jac_mass[3+j] = GNEWT^2*r0inv*rinv*(ddfdtdk2*x0[j]+dgdtmhdfdtm1dk2*v0[j])
  end
  if debug
    # Now include derivatives of gamma, r, fm1, dfdt, gmh, and dgdtmhdfdtm1:
    @inbounds for i=1:3
      delxv_jac[ 7,i] = -sqb*rinv*((g2-h*c3*r0inv^3)*v0[i]+c3*x0[i]*r0inv^3); delxv_jac[7,3+i] = sqb*rinv*((-d+2*g2*h-h^2*c3*r0inv^3)*v0[i]+(-g2+h*c3*r0inv^3)*x0[i])
      delxv_jac[ 8,i] = (c20*betainv-c2*c3*rinv)*r0inv^3*x0[i]+((eta*g2+g1*r0)*rinv+h*r0inv^3*(c2*c3*rinv-c20*betainv))*v0[i]
#      delxv_jac[8,3+i] = (g1-g2*c2*rinv+h*r0inv^3*(c2*c3*rinv-c20*betainv))*x0[i]+(-2g1*h+c18*betainv+(2g2*h-d)*c2*rinv+h^2*r0inv^3*(c20*betainv-c2*c3*rinv))*v0[i]
#      delxv_jac[8,3+i] = ((g1*r0+eta*g2)*rinv+h*r0inv^3*(c2*c3*rinv-c20*betainv))*x0[i]+(-2g1*h+c18*betainv+(2g2*h-d)*c2*rinv+h^2*r0inv^3*(c20*betainv-c2*c3*rinv))*v0[i]
      drdv0x0 = (beta*g1*g2+((eta*g2+k*g3)*eta*g0*c3)*rinv*r0inv^3 + (g1*g0*(2k*eta^2*g2+3eta*k^2*g3))*betainv*rinv*r0inv^2-
                          k*betainv*r0inv^3*(eta*g1*(eta*g2+k*g3)+g3*g0*r0^2*beta+2h*g2*k)+(g1*zeta)*rinv*((h*c3)*r0inv^3 - g2) -
                         (eta*(beta*g2*g0*r0+k*g1^2)*(eta*g1+k*g2))*betainv*rinv*r0inv^2)
#      delxv_jac[8,3+i] = drdv0x0*x0[i]+ (k*betainv*rinv*(eta*(2g0*g3-g1*g2+g3) - h6*k^2 + (g2^2 - 2*g1*g3)*beta*k*r0) + h*drdv0x0)*v0[i]
      delxv_jac[8,3+i] = drdv0x0*x0[i] - (k*betainv*rinv*(eta*(beta*g2*g3-h8) - h6*k + (g2^2 - 2*g1*g3)*beta*r0) + h*drdv0x0)*v0[i]
#                         +(-2g1*h+c18*betainv+(2g2*h-d)*c2*rinv+h^2*r0inv^3*(c20*betainv-c2*c3*rinv))*v0[i]
#      delxv_jac[8,3+i] = ((eta*g2+g1*r0)*rinv+h*betainv*rinv*r0inv^3*((3g1*g3 - 2g2^2)*k^3 - k^2*eta*h8 +
#            beta*k^2*(g2^2 - 3*g1*g3)*r0 - beta*r0^3 - k*g2*beta*(eta^2*g2 + 2eta*g1*r0 + g0*r0^2)))*x0[i]+(-2g1*h+c18*betainv+((g1*r0-g3*k-2*g0*h)*c2)*betainv*rinv +
#         h^2*((2*g2^2 - 3*g1*g3)*k^3 + beta*r0^3 + k^2*(eta*(g1*g2 - 3*g0*g3) + beta*(3*g1*g3 - g2^2)*r0) +
#         k*beta*g2*(eta*(eta*g2 + g1*r0) + r0*(eta*g1 + g0*r0)))*betainv*rinv*r0inv^3)*v0[i]
      delxv_jac[ 9,i] = dfm1dxx*x0[i]+dfm1dxv*v0[i]; delxv_jac[ 9,3+i]=dfm1dvx*x0[i]+dfm1dvv*v0[i]
      delxv_jac[10,i] = ddfdtdxx*x0[i]+ddfdtdxv*v0[i]; delxv_jac[10,3+i]=ddfdtdvx*x0[i]+ddfdtdvv*v0[i]
      delxv_jac[11,i] = dgmhdxx*x0[i]+dgmhdxv*v0[i]; delxv_jac[11,3+i]=dgmhdvx*x0[i]+dgmhdvv*v0[i]
      delxv_jac[12,i] = dgdtmhdfdtm1dxx*x0[i]+dgdtmhdfdtm1dxv*v0[i]; delxv_jac[12,3+i]=dgdtmhdfdtm1dvx*x0[i]+dgdtmhdfdtm1dvv*v0[i]
    end
    delxv_jac[ 7,7] = sqb*c1*r0inv*rinv; delxv_jac[7,8] = sqb*rinv*(1+eta*c3*r0inv^3+g2*(2k*r0inv-beta))
#    delxv_jac[ 8,7] = (c17*betainv+c1*c2*rinv)*r0inv
    delxv_jac[ 8,7] = betainv*r0inv*rinv*(-g2*r0^2*beta-eta*g1*g2*(k+beta*r0)+eta*g0*g3*(2*k+zeta)-
           g2^2*(beta*eta^2+2*k*zeta)+g1*g3*zeta*(3*k-beta*r0)); delxv_jac[8,8] = c2*rinv
    delxv_jac[ 8,8] = ((r0*g1+eta*g2)*rinv)*(beta-2*k*r0inv)+c2*rinv+eta*r0inv^3*(c2*c3*rinv-c20*betainv)
    delxv_jac[ 9,7] = dfm1dk; delxv_jac[ 9,8] = dfm1dh
    delxv_jac[10,7] = ddfdtdk; delxv_jac[10,8] = ddfdtdh
    delxv_jac[11,7] = dgmhdk; delxv_jac[11,8] = dgmhdh
    delxv_jac[12,7] = dgdtmhdfdtm1dk; delxv_jac[12,8] = dgdtmhdfdtm1dh
  end
else
  # Now compute the Kepler-Drift Jacobian terms:
  # First, x0 derivatives:
  #  delxv[  j] = fm1*x0[j]+gmh*v0[j]        # position x_ij(t+h)-x_ij(t) - h*v_ij(t) or -h*v_ij(t+h)
  #  delxv[3+j] = dfdt*x0[j]+dgdtm1*v0[j]    # velocity v_ij(t+h)-v_ij(t)
  # First, the diagonal terms:
  # Now the off-diagonal terms:
  d   = (h + eta*g2 + 2*k*g3)*betainv
  c1 = d-r0*g3
  c2 = eta*g0+g1*zeta
  c3  = d*k+g1*r0^2
  c4 = eta*g1+2*g0*r0
  c6 = (r0*g0-k*g2)*betainv
  c9  = g2*r-h1*k
  c14 = r0*g2-k*h1
  c15 = eta*h1+h2*r0
  c16 = eta*h2+g1*gamma*r0/sqb
  c17 = r0-r-g2*k
  c18 = eta*g1+2*g2*k
  c19 = 4*eta*h1+3*h2*r0
  c23 = h2*k-r0*g1
  h6 = H6(gamma,beta)
  h8 = H8(gamma,beta)
  # Derivatives of \delta x with respect to x0, v0, k & h:
  dfm1dxx = k*rinv^3*betainv*r0inv^4*(k*h1*r^2*r0*(beta-2*k*r0inv)+beta*c3*(r*c23+c14*c2)+c14*r*(k*(r-g2*k)+g0*r0*zeta))
  dfm1dxv = k*rinv^2*r0inv*(k*(g2*h2+g1*h1)-2g1*g2*r0+g2*c14*c2*rinv)
  dfm1dvx = dfm1dxv
#  dfm1dvv = k*r0inv*rinv^2*betainv*(r*(2*g2*r0-4*h1*k)+d*beta*c23-c18*c14+d*beta*c14*c2*rinv)
  dfm1dvv = k*r0inv*rinv^2*betainv*(2eta*k*(g2*g3-g1*h1)+(3g3*h2-4h1*g2)*k^2 +
    beta*g2*r0*(3h1*k-g2*r0)+c14*rinv*(-beta*g2^2*eta^2+eta*k*(2g0*g3-h2)-
     h6*k^2+(-2eta*g1*g2+k*(h1-2g1*g3))*beta*r0-beta*g1^2*r0^2))
  dfm1dh  = (g1*k-h2*k^2*r0inv-k*c14*c2*rinv*r0inv)*rinv^2
  dfm1dk  = rinv*r0inv*(4*h1*k^2*betainv*r0inv-k*h1-2*g2*k*betainv+c14-k*c14*c17*betainv*rinv*r0inv+
        k*(g1*r0-k*h2)*c1*rinv*r0inv-k*c14*c1*c2*rinv^2*r0inv)
#  dfm1dk2_old  = 4*h1*k*betainv*r0inv-h1-2*g2*betainv-c14*c17*betainv*rinv*r0inv+ (g1*r0-k*h2)*c1*rinv*r0inv-c14*c1*c2*rinv^2*r0inv
  # New expression for d(f-1-h \dot f)/dk with cancellations of higher order terms in gamma is:
  dfm1dk2  = betainv*r0inv*rinv^2*(r*(2eta*k*(g1*h1-g3*g2)+(4g2*h1-3g3*h2)*k^2-eta*r0*beta*g1*h1 + (g3*h2-4g2*h1)*beta*k*r0 + g2*h1*beta^2*r0^2) -
  # In the following line I need to replace 3g0*g3-g1*g2 by -H8:
          c14*(-eta^2*beta*g2^2 - k*eta*h8 - k^2*h6 - eta*r0*beta*(g1*g2 + g0*g3) + 2*(h1 - g1*g3)*beta*k*r0 - (g2 - beta*g1*g3)*beta*r0^2))
#  println("dfm1dk2: old ",dfm1dk2_old," new: ",dfm1dk2)
  dgmhdxx = k*rinv*r0inv*(h2+k*c19*betainv*r0inv^2-c16*c3*rinv*r0inv^2+c2*c3*c15*(rinv*r0inv)^2-c15*(k*(g2*k+r)-g0*r0*zeta)*betainv*rinv*r0inv^2)
  dgmhdxv = k*rinv^2*(h1*r-g2*c16-g1*c15+g2*c2*c15*rinv)
  dgmhdvx = dgmhdxv
#  dgmhdvv = k*rinv^2*(-d*c16-c15*c18*betainv+r*c19*betainv+d*c2*c15*rinv)
  dgmhdvv = k*betainv*rinv^2*(2*eta^2*(g1*h1-g2*g3)+eta*k*(4g2*h1-3h2*g3)+r0*eta*(4g0*h1-2g1*g3)+
  # In the following lines I need to replace g1*g2-3g0*g3-g1*g2 by H8:
           3r0*k*((g1+beta*g3)*h1-g3*g2)+(g0*h8-beta*g1*(g2^2+g1*g3))*r0^2 -
           c15*rinv*(beta*g2^2*eta^2+eta*k*h8+h6*k^2+(2eta*g1*g2-k*(g2^2-3g1*g3))*beta*r0+beta*g1^2*r0^2))
  dgmhdk  = rinv*(k*c1*c16*rinv*r0inv+c15-k*c15*c17*betainv*rinv*r0inv-k*c19*betainv*r0inv-k*c1*c2*c15*rinv^2*r0inv)
#  dgmhdk2_old  = c1*c16*rinv-c15*c17*betainv*rinv-c19*betainv-c1*c2*c15*rinv^2
  h7 = H7(gamma,beta)
  dgmhdk2 =  betainv*rinv^2*(r*(2eta^2*(g3*g2-g1*h1) + eta*k*(3g3*h2 - 4g2*h1) +
    r0*eta*(beta*g3*(g1*g2 + g0*g3) - 2g0*h6) + (-h6*(g1 + beta*g3) + g2*(2g3 - h2))*r0*k +
    (h7 - beta^2*g1*g3^2)*r0^2)- c15*(-beta*eta^2*g2^2 + eta*k*(-h2 + 2g0*g3) - h6*k^2 -
    r0*eta*beta*(h2 + 2g0*g3) + 2beta*(2*h1 - g2^2)*r0*k + beta*(beta*g1*g3 - g2)*r0^2))
#  println("gmhgdot: old ",dgmhdk2_old," new: ",dgmhdk2)
  dgmhdh  = k*rinv^3*(r*c16-c2*c15)
  @inbounds for j=1:3
    # First, compute the \delta x-derivatives:
    delxv_jac[  j,  j] = fm1
    delxv_jac[  j,3+j] = gmh
    @inbounds for i=1:3
      delxv_jac[j,  i] += (dfm1dxx*x0[i]+dfm1dxv*v0[i])*x0[j] + (dgmhdxx*x0[i]+dgmhdxv*v0[i])*v0[j]
      delxv_jac[j,3+i] += (dfm1dvx*x0[i]+dfm1dvv*v0[i])*x0[j] + (dgmhdvx*x0[i]+dgmhdvv*v0[i])*v0[j]
    end
    delxv_jac[  j,  7] = dfm1dk*x0[j] + dgmhdk*v0[j]
    delxv_jac[  j,  8] = dfm1dh*x0[j] + dgmhdh*v0[j]
    jac_mass[  j] = GNEWT^2*rinv*r0inv*(dfm1dk2*x0[j]+dgmhdk2*v0[j])
  end
  # Derivatives of \delta v with respect to x0, v0, k & h:
  c5 = (r0-k*g2)*rinv/g1
  c7 = g2*(1/g1+c2*rinv)
  c8 = (k*c6+r*r0+c3*c5)*r0inv^3
  c12 = g0*h-g1*r0
  c20 = k*(g2*k+r)-g0*r0*zeta
  ddfdtdxx = dfdt*(eta*g1*rinv-2-g0*c3*rinv*r0inv/g1+c2*c3*r0inv*rinv^2-k*(k*g2-r0)*betainv*rinv*r0inv)*r0inv^2
  ddfdtdxv = -dfdt*(g0*g2/g1+(r0*g1+eta*g2)*rinv)*rinv
  ddfdtdvx = ddfdtdxv
#  ddfdtdvv = dfdt*(betainv-d*g0*rinv/g1-c18*betainv*rinv+d*c2*rinv^2)
  ddfdtdvv = -k*rinv^3*r0inv*betainv*((beta*eta*g2^2+k*h8)*(r0*g0+k*g2)+
       g1*(- h6*k^2 + (-2eta*g1*g2+(h1-2g1*g3)*k)*beta*r0 - beta*g1^2*r0^2))
  ddfdtdk  = dfdt*(1/k+c1*(r0-g2*k)*r0inv*rinv^2/g1-betainv*r0inv*(1+c17*rinv))
#  ddfdtdk2  = -g1*(c1*(r0-g2*k)*r0inv*rinv^2/g1-betainv*r0inv*(1+c17*rinv))
#  ddfdtdk2  = r0inv*(g1*c17*betainv*rinv+g1*betainv-g1*c1*c2*rinv^2-c1*g0*rinv)
  h3 = H3(gamma,beta)
  ddfdtdk2  = (r0-g2*k)*betainv*r0inv*rinv^2*(-eta*beta*g2^2+h3*k+(g3-g1*g2)*beta*r0)
  ddfdtdh  = dfdt*(r0-g2*k)*rinv^2/g1
  dgdotm1dxx = rinv^2*r0inv^3*((eta*g2+g1*r0)*k*c3*rinv+g2*k*(k*(g2*k-r)-g0*r0*zeta)*betainv)
  dgdotm1dxv = k*g2*rinv^3*(r*g1+r0*g1+eta*g2)
  dgdotm1dvx = dgdotm1dxv
#  dgdotm1dvv = k*rinv^2*(d*g1+g2*c18*betainv-2*r*g2*betainv-d*g2*c2*rinv)
  dgdotm1dvv = k*betainv*rinv^3*(eta^2*beta*g2^3-eta*k*g2*h3+3r0*eta*beta*g1*g2^2 +
         r0*k*(-g0*h6+3beta*g1*g2*g3)+beta*g2*(g0*g2+g1^2)*r0^2)
  dgdotm1dk = rinv*r0inv*(-r0*g2+g2*k*(r+r0-g2*k)*betainv*rinv-k*g1*c1*rinv+k*g2*c1*c2*rinv^2)
#  dgdotm1dk2 = rinv*(g2*(r+r0-g2*k)*betainv-g1*c1+g2*c1*c2*rinv)
  dgdotm1dk2 = betainv*rinv^2*(-beta*eta^2*g2^3+eta*k*g2*h3+eta*r0*beta*g2*(g3-2g1*g2)+
       (h6-beta*g2^3)*r0*k + beta*g1*(g3 - g1*g2)*r0^2)
#       (g2^2*(1+g0)-3*g1*g3)*r0*k + beta*g1*(g3 - g1*g2)*r0^2)
  dgdotm1dh = k*rinv^3*(g2*c2-r*g1)
  @inbounds for j=1:3
    # Next, compute the \delta v-derivatives:
    delxv_jac[3+j,  j] = dfdt
    delxv_jac[3+j,3+j] = dgdtm1
    @inbounds for i=1:3
      delxv_jac[3+j,  i] += (ddfdtdxx*x0[i]+ddfdtdxv*v0[i])*x0[j] + (dgdotm1dxx*x0[i]+dgdotm1dxv*v0[i])*v0[j]
      delxv_jac[3+j,3+i] += (ddfdtdvx*x0[i]+ddfdtdvv*v0[i])*x0[j] + (dgdotm1dvx*x0[i]+dgdotm1dvv*v0[i])*v0[j]
    end
    delxv_jac[3+j,  7] = ddfdtdk*x0[j] + dgdotm1dk*v0[j]
    delxv_jac[3+j,  8] = ddfdtdh*x0[j] + dgdotm1dh*v0[j]
    jac_mass[3+j] = GNEWT^2*rinv*r0inv*(ddfdtdk2*x0[j]+dgdotm1dk2*v0[j])
  end
  if debug
    # Now include derivatives of gamma, r, fm1, gmh, dfdt, and dgdtmhdfdtm1:
    @inbounds for i=1:3
      delxv_jac[ 7,i] = -sqb*rinv*(g2*v0[i]+c3*x0[i]*r0inv^3); delxv_jac[7,3+i] = -sqb*rinv*(d*v0[i]+g2*x0[i])
      delxv_jac[ 8,i] = (c20*betainv-c2*c3*rinv)*r0inv^3*x0[i]+((r0*g1+eta*g2)*rinv)*v0[i]
      delxv_jac[8,3+i] = (c18*betainv-d*c2*rinv)*v0[i]+((r0*g1+eta*g2)*rinv)*x0[i]
#      delxv_jac[8,3+i] = (g2*(2k*(r-r0)+beta*(r0^2+eta^2*g2)-zeta*eta*g1) + c2*g3*(2k+zeta))*betainv*rinv*r0inv*v0[i]+((r0*g1+eta*g2)*rinv)*x0[i]
      delxv_jac[ 9,i] = dfm1dxx*x0[i]+dfm1dxv*v0[i]; delxv_jac[ 9,3+i]=dfm1dvx*x0[i]+dfm1dvv*v0[i]
      delxv_jac[10,i] = ddfdtdxx*x0[i]+ddfdtdxv*v0[i]; delxv_jac[10,3+i]=ddfdtdvx*x0[i]+ddfdtdvv*v0[i]
      delxv_jac[11,i] = dgmhdxx*x0[i]+dgmhdxv*v0[i]; delxv_jac[11,3+i]=dgmhdvx*x0[i]+dgmhdvv*v0[i]
      delxv_jac[12,i] = dgdotm1dxx*x0[i]+dgdotm1dxv*v0[i]; delxv_jac[12,3+i]=dgdotm1dvx*x0[i]+dgdotm1dvv*v0[i]
    end
    delxv_jac[ 7,7] = sqb*c1*r0inv*rinv; delxv_jac[7,8] = sqb*rinv
#    delxv_jac[ 8,7] = (c17*betainv+c1*c2*rinv)*r0inv; delxv_jac[8,8] = c2*rinv
    delxv_jac[ 8,7] = betainv*r0inv*rinv*(-g2*r0^2*beta-eta*g1*g2*(k+beta*r0)+eta*g0*g3*(2*k+zeta)-
           g2^2*(beta*eta^2+2*k*zeta)+g1*g3*zeta*(3*k-beta*r0)); delxv_jac[8,8] = c2*rinv
    delxv_jac[ 9,7] = dfm1dk; delxv_jac[ 9,8] = dfm1dh
    delxv_jac[10,7] = ddfdtdk; delxv_jac[10,8] = ddfdtdh
    delxv_jac[11,7] = dgmhdk; delxv_jac[11,8] = dgmhdh
    delxv_jac[12,7] = dgdotm1dk; delxv_jac[12,8] = dgdotm1dh
  end
end
#return delxv_jac::Array{T,2}
return
end


text:
module TestSndfile

include("testhelpers.jl")

using AudioIO
using FactCheck
import AudioIO: DeviceInfo, render, AudioSample, AudioBuf

facts("WAV file write/read") do
    fname = Pkg.dir("AudioIO", "test", "sinwave.wav")

    srate = 44100
    freq = 440
    t = [0 : 2 * srate - 1] / srate
    phase = 2 * pi * freq * t
    reference = int16((2 ^ 15 - 1) * sin(phase))

    AudioIO.open(fname, "w") do f
        write(f, reference)
    end

    # test basic reading
    AudioIO.open(fname) do f
        @fact f.sfinfo.channels => 1
        @fact f.sfinfo.frames => 2 * srate
        actual = read(f)
        @fact length(reference) => length(actual)
        @fact reference => actual[:, 1]
        @fact samplerate(f) => srate
    end

    # test seeking

    # test rendering as an AudioNode
    AudioIO.open(fname) do f
        # pretend we have a stream at the same rate as the file
        bufsize = 1024
        input = zeros(AudioSample, bufsize)
        test_info = DeviceInfo(srate, bufsize)
        node = FilePlayer(f)
        # convert to floating point because that's what AudioIO uses natively
        expected = convert(AudioBuf, reference ./ (2^15))
        buf = render(node, input, test_info)
        @fact expected[1:bufsize] => buf[1:bufsize]
        buf = render(node, input, test_info)
        @fact expected[bufsize+1:2*bufsize] => buf[1:bufsize]
    end
end

facts("Stereo file reading") do
    fname = Pkg.dir("AudioIO", "test", "440left_880right.wav")
    srate = 44100
    t = [0 : 2 * srate - 1] / srate
    expected = int16((2^15-1) * hcat(sin(2pi*t*440), sin(2pi*t*880)))

    AudioIO.open(fname) do f
        buf = read(f)
        @fact buf => mse(expected, 5)
    end
end

# note - currently AudioIO just mixes down to Mono. soon we'll support this
# new-fangled stereo sound stuff
facts("Stereo file rendering") do
    fname = Pkg.dir("AudioIO", "test", "440left_880right.wav")
    srate = 44100
    bufsize = 1024
    input = zeros(AudioSample, bufsize)
    test_info = DeviceInfo(srate, bufsize)
    t = [0 : 2 * srate - 1] / srate
    expected = convert(AudioBuf, 0.5 * (sin(2pi*t*440) + sin(2pi*t*880)))

    AudioIO.open(fname) do f
        node = FilePlayer(f)
        buf = render(node, input, test_info)
        @fact buf[1:bufsize] => mse(expected[1:bufsize])
        buf = render(node, input, test_info)
        @fact buf[1:bufsize] => mse(expected[bufsize+1:2*bufsize])
    end
end

end # module TestSndfile


text:
immutable Location
  file::String
  line::Int
end

const nofile = Location("",-1)

type Branch
  loc::Location
end

immutable History
  data::Vector{Tuple{Branch,Bool}}
end

History() = History([])

Base.push!(h::History, x::Tuple{Branch,Bool}) = push!(h.data, x)

branch!(h::History, b::Branch, x::Bool) = (push!(h, (b, x)); x)

@inline trace(f, hist, xs...) = f(xs...)

function tracemethod!(l::LambdaInfo)
  F = l.specTypes.parameters[1]
  switchfunc!(l, trace)
  insertarg!(l, 2, (symbol("#unused#"),F))
  insertarg!(l, 3, (:history,History))
  loc = nofile
  map!(l.code) do ex
    isexpr(ex, :line) && (loc = Location(string(ex.args[2]), ex.args[1]))
    isa(ex, LineNumberNode) && (loc = Location(loc.file, ex.line))
    isexpr(ex, :gotoifnot) || return ex
    Expr(:gotoifnot,
         Expr(:call, GlobalRef(HotFuzz, :branch!), SlotNumber(3), Branch(loc), ex.args[1]),
         ex.args[2])
  end
  return l
end

function trace!(meth::Method)
  linfo = copyil(meth.lambda_template)
  tracemethod!(linfo)
  jl_method_def(linfo)
  return
end

function runtrace(f, args...)
  h = History()
  trace(f, h, args...), h.data
end

# trace!(@which(gcd(1, 1)))

# runtrace(gcd, rand(0:100), rand(0:100))


text: abstract type CollisionOperator{M,N,ℳ,𝒩,RT,CT} end


text:
# time julia -e 'using TensorCast; TensorCast._macro(:(  Z[i,k][j] := fun(A[i,:], B[j])[k] + C[k]^2  ))'
# 10.1 sec without this, 7.2 sec with (on 2nd run)

_macro(:(  Z[i,j] := A[i] + B[j]  ))
_macro(:(  Z[i,_,k'] := A[i] + B[k'] / log(2) ))
_macro(:(  Z[(i,j)] := A[j,-i] ))
_macro(:(  Z[i,k][j] := fun(A[i,:], B[j])[k] + C[k]^2 ))
_macro(:(  S[i,j] = sum(k)  ),:(  A[i,j] + B[j,k]  ), call=CallInfo(:reduce))
_macro(:(  S[i,j] := sum(k)  ),:(  (B[j]+ C[k])[i]  ), call=CallInfo(:reduce))

pretty(:(  (B[j]+ C[k])[i] ))
pretty(@macroexpand @cast A[i,j] := B[j,i] + 1 )


text:
"""
    Rf, m = calib_force(POSES, F, m0::Real = 0.3; offset=true)
    Rf, m = calib_force(POSES, F, g::Vector; offset=true)

`POSES` ∈ ℜ(4,4,N) are the positions of the tool frame from the robot forward kinematics, 4x4 transformation matrices
`F` ∈ ℜ(N,3) vector of forces (accepts ℜ(Nx6) matrix with torques also)
# Usage
```
Rf*force[i,1:3] + forceoffs = POSES[1:3,1:3,i]'*[0, 0, mf*-9.82]
```
If `m0` is supplied, it's assumed that the gravity vector is [0,0,-g], or in words, that the gravity is acting along the negative z axis. A custom gravity vector `g` can also be supplied.

This function can also be used to calibrate an accelerometer.

[Bagge Carlson, F.](https://www.control.lth.se/staff/fredrik-bagge-carlson/), ["Machine Learning and System Identification for Estimation in Physical Systems"](https://lup.lub.lu.se/search/publication/ffb8dc85-ce12-4f75-8f2b-0881e492f6c0) (PhD Thesis 2018).

See also `calib_force_iterative` and `calib_force_eigen`.
"""
function calib_force(POSES,F,m0::Real=0.3; kwargs...)
    g0 = [0,0,-m0*9.82]
    calib_force(POSES,F,g0; kwargs...)
end
function calib_force(POSES,F,g::AbstractVector; offset=true, verbose=true)

    N = size(POSES,3)

    local forceoffs
    I = Robotlib.I3
    A  = Array{eltype(F)}(undef, 3N, offset ? 12 : 9)
    B  = Array{eltype(F)}(undef, 3N)
    A2 = Array{eltype(F)}(undef, 3N, offset ? 4 : 1)
    B2 = Array{eltype(F)}(undef, 3N)
    At = offset ? (F,i) -> [F[i,1]*I   F[i,2]*I    F[i,3]*I   -I] : (F,i) -> [F[i,1]*I   F[i,2]*I    F[i,3]*I]

    for i = 1:N
        RA               = POSES[1:3,1:3,i]'
        b                = RA*g
        A[3(i-1)+1:3i,:] = At(F,i)
        B[3(i-1)+1:3i]   = b
    end
    verbose && @info("Condition number of Gram matrix: ", cond(A'A))
    w  = tls(A,B)
    Rf = reshape(w[1:9],3,3)
    verbose && @info("Determinant of Rf before orthonormalization: ", det(Rf), opnorm(Rf))
    if det(Rf) < 0
        verbose && @warn("det(Rf) < 0, left handed coordinate system? Trying to solve the problem with reversed gravity (-g)")
        return calib_force(POSES,F,-g; offset=offset, verbose=verbose)
    end
    orthonormal!(Rf)
    At = offset ? RA -> [RA*g/norm(g) -I] : RA -> RA*g/norm(g)

    for i = 1:N
        RA                = POSES[1:3,1:3,i]'
        b                 = Rf*F[i,1:3]
        A2[3(i-1)+1:3i,:] = At(RA)
        B2[3(i-1)+1:3i]   = b
    end
    w         = A2\B2
    m         = w[1]/9.82
    offset && (forceoffs = w[2:4])

    m < 0 && @error("Estimated mass is negative, left handed coordinate system for sensor?")
    if offset
        return Rf,m,forceoffs
    else
        return Rf, m
    end
end

"""
    Rf, g, m = calib_force_eigen(POSES, F, g)

`POSES` ∈ ℜ(4,4,N) or ℜ(3,3,N) are the orientations of the tool frame from the robot forward kinematics, 4×4 transformation matrices or 3×3 rotation matrices.
`F` ∈ ℜ(N,3) vector of forces (accepts ℜ(N×6) matrix with torques also)
´g´ is an initial guess for the gravity vector.
# Usage
```
Rf*force[i,1:3] = POSES[1:3,1:3,i]'*g
```
See also `calib_force_eigen`.
"""
function calib_force_iterative(POSES, F, g; trace = false)
    N = size(POSES, 3)
    I = Robotlib.I3
    A = Array{eltype(F)}(undef, 3N, 3)
    B = F'[:]
    local m, Rf
    trace && (Rg = [])
    for iter = 1:6
        Rf, m = calib_force(
            POSES,
            F,
            g;
            offset = false,
            verbose = false,
        )
        for i = 1:N
            A[3(i-1)+1:3i, :] = Rf'POSES[:, :, i]'
        end
        g = A \ B
        trace && push!(Rg, (Rf, g))
    end
    trace && (return Rf, g, m, Rg)
    Rf, g, m
end

function eigenR(K)
    v = real(eigen(K).vectors[:, end])
    R = reshape(v, 3, 3)
    det(R) < 0 && (R .*= -1)
    orthonormal(R)
end

"""
    calib_force_eigen(POSES, F)

`POSES` ∈ ℜ(4,4,N) or ℜ(3,3,N) are the orientations of the tool frame from the robot forward kinematics, 4×4 transformation matrices or 3×3 rotation matrices.
`F` ∈ ℜ(N,3) vector of forces (accepts ℜ(N×6) matrix with torques also)
# Usage
```
Rf*force[i,1:3] + forceoffs = POSES[1:3,1:3,i]'*[0, 0, mf*-9.82]
```
"""
function calib_force_eigen(POSES, forces)
    N = size(POSES, 3)
    D = reshape(POSES[1:3,1:3,:], 3, 3N)'
    F = kron(forces, I(3))
    K = F \ D * (D \ F)
    R̂ = eigenR(K)
    ĝ = D \ F * vec(R̂)
    R̂, ĝ
end

import Robotlib: skew

"""
    This function uses Cayley transform to solve for R. It requires a known mass so it is recommended to use `calib_force` instead.
"""
function calib_force2(POSES,F,m,g  = -9.82)
    N  = size(POSES,3)

    mg = m*g

    I  = Robotlib.I3
    A  = Array{eltype(F)}(undef, 3N, 3)
    B  = Array{eltype(F)}(undef, 3N)

    for i = 1:N
        RA               = POSES[3,1:3,i]
        A[3(i-1)+1:3i,:] = skew(F[i,1:3][:] + RA*mg)
        B[3(i-1)+1:3i]   = RA*mg - F[i,1:3]
    end
    @info("Condition number of Gram matrix: ", cond(A'A))
    w  = tls(A,B)
    S  = skew(w)
    Rf = (I+S)\(I-S)
end

max_forks_repo_name: baggepinnen/Robotlib


text:
"""
This funtion does nothing if the product is not an operating market product.
"""
function get_maxperc(product::T, scenario_name::String, num_years::Int64, num_hours::Int64) where T <: Product
    return
end

"""
This funtion returns the array of capacity factors if product is energy.
"""
function get_maxperc(product::Energy, scenario_name::String, num_years::Int64, num_hours::Int64)
    return get_capacity_factors(product)[scenario_name]
end

"""
This funtion returns the array of maximum percentage limit if product is operating reserve.
"""
function get_maxperc(product::OperatingReserve, scenario_name::String, num_years::Int64, num_hours::Int64)
    return ones(num_years, num_hours) * get_max_limit(product)
end

"""
This function does nothing if the product is not an operating market product.
"""
function push_operating!(array::Array{OperatingProduct}, product::T) where T <: Product
    return
end

"""
This function pushes the product to the array if product is an operating market product.
"""
function push_operating!(array::Array{OperatingProduct}, product::OperatingProduct)
    push!(array, product)
    return
end

function push_energy_product!(array::Array{Energy}, product::OperatingProduct)
    return
end

function push_energy_product!(array::Array{Energy}, product::Energy)
    push!(array, product)
    return
end

function find_energy_product(products::Array{OperatingProduct})
    array = Energy[]

    for product in products
        push_energy_product!(array, product)
    end

    return array
end

"""
This function returns an array of only operating market products from an array of all products.
"""
function find_operating_products(products::Array{Product})
    operating_products = OperatingProduct[]

    for i = 1:length(products)
        push_operating!(operating_products, products[i])
    end

return  operating_products
end