text:
#f(x,y)=mean(KL(x,sigmoid(y))), where KL is the Kullback-Leibler divergence

FBinaryEntropyLossXsigmoidY(x,y)=([sum(x.*log(x)+(1.-x).*log(1.-x)-x.*y+log1pexp(y))/length(x)],nothing)
FBinaryEntropyLossXsigmoidY_inplace(value,aux,x,y)=copy!(value,sum(x.*log(x)+(1.-x).*log(1.-x)-x.*y+log1pexp(y))/length(x))

function DBinaryEntropyLossXsigmoidY(derivativeIDX,f_c,faux_c,grad_c,grad_n,x,y)
    if derivativeIDX==1
        axpy!(-grad_c[1]/length(x),y-log(x./(1-x)),grad_n)
    elseif derivativeIDX==2
        axpy!(grad_c[1]/length(x),sigmoid(y)-x,grad_n)
    end
end

if PROC=="GPU"
#    FBinaryEntropyLossXsigmoidY(x::CudaArray,y::CudaArray)=(binaryentropyXsigmoidY(x,y),nothing)
    FBinaryEntropyLossXsigmoidY_inplace(value,aux,x::CudaArray,y::CudaArray)=binaryentropyXsigmoidY!(x,y,value)

    function DBinaryEntropyLossXsigmoidY(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray,y::CudaArray)
        if derivativeIDX==1
            DXbinaryentropyXsigmoidY!(x,y,grad_c,grad_n)
        elseif derivativeIDX==2
            DYbinaryentropyXsigmoidY!(x,y,grad_c,grad_n)
        end
    end

#    FBinaryEntropyLossXsigmoidY_inplace(value,aux,x::CudaArray,y::CudaArray)=binaryentropyXsigmoidY!(x,y,value)

#    function DBinaryEntropyLossXsigmoidY(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray{Float64},y::CudaArray{Float64})
#        if derivativeIDX==1
#            DXbinaryentropyXsigmoidY!(x,y,grad_c,grad_n)
#        elseif derivativeIDX==2
#            DYbinaryentropyXsigmoidY!(x,y,grad_c,grad_n)
#        end
#    end

#    function DBinaryEntropyLossXsigmoidY(derivativeIDX,f_c,faux_c,grad_c::CudaArray{Float32},grad_n::CudaArray{Float32},x::CudaArray{Float32},y::CudaArray{Float32})
#        if derivativeIDX==1
#            DXbinaryentropyXsigmoidY_32!(x,y,grad_c,grad_n)
#        elseif derivativeIDX==2
#            DYbinaryentropyXsigmoidY_32!(x,y,grad_c,grad_n)
#        end
#    end

end

Derivative[FBinaryEntropyLossXsigmoidY]=DBinaryEntropyLossXsigmoidY
Inplace[FBinaryEntropyLossXsigmoidY]=FBinaryEntropyLossXsigmoidY_inplace

BinaryKullbackLeiblerLossXsigmoidY(nx,ny)=ADnode(FBinaryEntropyLossXsigmoidY,[nx ny]) # give it a better name
export BinaryKullbackLeiblerLossXsigmoidY

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
println([i for i in 1:5])

m = [3α + 2β for α in 1:3, β in 2:3]

println(m)

ext: jl

max_stars_repo_name: jameszhan/notes-ml

max_issues_repo_name: jameszhan/notes-ml

max_forks_repo_name: jameszhan/notes-ml


text:
@testset "MultiResolutionHybrid" begin
    img = generate_surface(SynthSphere(50), 0.5)
    p,q = MultiResolutionHybrid()(img)
    @test_broken 1 == 2
end


text:
@reexport module _Algebra

import Base: sqrt, rand, log, zero, one, convert, colon
import Nemo: root, roots

export sqrts, elements, multiple

using Brainstorm._Algorithm: @forcartesian
using Brainstorm._Math: factor
using AbstractAlgebra: GroupElem
using Nemo: characteristic, degree, fmpz, FinFieldElem, FqFiniteField, FqNmodFiniteField, gen, isone, iszero, order, remove

convert(::Type{Integer}, n::fmpz) = BigInt(n)

colon(a::Integer, b::fmpz) = colon(fmpz(a), b)

include("group.jl")
include("field.jl")

end


text:
using KitBase
import KitBase.BSON
cd(@__DIR__)

# save
set = Setup(space = "1d1f1v", boundary = "fix", maxTime = 0.1)
ps = PSpace1D(0, 1, 100, 1)
vs = VSpace1D(-5, 5, 28)
gas = Gas(Kn = 5e-4, K = 0.0, γ = 3.0)
ib = begin
    primL = [1.0, 0.0, 1 / 0.5]
    primR = [0.125, 0.0, 1 / 0.625]
    p = (γ = gas.γ, u = vs.u, primL = primL, primR = primR)
    _fw = function (x, p)
        prim = ifelse(x < 0.5, p.primL, p.primR)
        return prim_conserve(prim, p.γ)
    end
    _ff = function (x, p)
        prim = ifelse(x < 0.5, p.primL, p.primR)
        return maxwellian(p.u, prim)
    end
    _bc = function (x, p)
        prim = ifelse(x < 0.5, p.primL, p.primR)
        return prim
    end
    IB1F(_fw, _ff, _bc, p)
end
ks = SolverSet(set, ps, vs, gas, ib)
BSON.@save "ks.bson" ks

# load
using OffsetArrays
BSON.@load "ks.bson" ks

w0, f0, prim0 = ks.ib(1)


text:
# deprecated functions

mutable struct GenericWrapParam{T} <: FunctorInferenceType
  #TODO: <: FunctorIT cannot be right here -- Used in one of the unpacking converters
  usrfnc!::T
  params::Vector{Array{Float64,2}}
  varidx::Int
  particleidx::Int
  measurement::Tuple #Array{Float64,2}
  samplerfnc::Function # TODO -- remove, since no required. Direct multiple dispatch at solve
  specialzDim::Bool
  partial::Bool
  # hypoverts::Vector{Symbol}
  hypotheses::Union{Void, Distributions.Categorical}
  activehypo::Union{UnitRange{Int},Vector{Int}}
  factormetadata::FactorMetadata
  GenericWrapParam{T}() where {T} = new()
  GenericWrapParam{T}(fnc::T, t::Vector{Array{Float64,2}}) where {T} = new(fnc, t, 1,1, (zeros(0,1),) , +, false, false, nothing, 1:length(t), FactorMetadata())
  GenericWrapParam{T}(fnc::T, t::Vector{Array{Float64,2}}, varidx::Int, prtcl::Int) where {T} = new(fnc, t, varidx, prtcl, (zeros(0,1),) , +, false, false, nothing, 1:length(t), FactorMetadata())
  GenericWrapParam{T}(fnc::T, t::Vector{Array{Float64,2}}, i::Int, j::Int, meas::Tuple, smpl::Function) where {T} = new(fnc, t, i, j, meas, smpl, false, false, nothing, 1:length(t), FactorMetadata())
  GenericWrapParam{T}(fnc::T, t::Vector{Array{Float64,2}}, i::Int, j::Int, meas::Tuple, smpl::Function, szd::Bool) where {T} = new(fnc, t, i, j, meas, smpl, szd, false, nothing, 1:length(t), FactorMetadata())
  GenericWrapParam{T}(fnc::T, t::Vector{Array{Float64,2}}, i::Int, j::Int, meas::Tuple, smpl::Function, szd::Bool, partial::Bool) where {T} = new(fnc, t, i, j, meas, smpl, szd, partial, nothing, 1:length(t), FactorMetadata())
  GenericWrapParam{T}(fnc::T, t::Vector{Array{Float64,2}}, i::Int, j::Int, meas::Tuple, smpl::Function, szd::Bool, partial::Bool, mhcat::Union{Void,Categorical}) where {T} = new(fnc, t, i, j, meas, smpl, szd, partial, mhcat, 1:length(t), FactorMetadata())
  GenericWrapParam{T}(fnc::T, t::Vector{Array{Float64,2}}, i::Int, j::Int, meas::Tuple, smpl::Function, szd::Bool, partial::Bool, mhcat::Tuple) where {T} = new(fnc, t, i, j, meas, smpl, szd, partial, Categorical(mhcat), 1:length(t), FactorMetadata())
end

mutable struct FastRootGenericWrapParam{T} <: Function
  p::Vector{Int}
  perturb::Vector{Float64}
  X::Array{Float64,2}
  Y::Vector{Float64}
  xDim::Int
  zDim::Int
  gwp::GenericWrapParam{T}
  gg::Function
  res::Vector{Float64}
  FastRootGenericWrapParam{T}(xArr::Array{Float64,2}, zDim::Int, residfnc::GenericWrapParam{T}) where {T} =
      new(collect(1:size(xArr,1)), zeros(zDim), xArr, zeros(size(xArr,1)), size(xArr,1), zDim, residfnc, +, zeros(0))
end

function packmultihypo(fnc::GenericWrapParam{T}) where {T<:FunctorInferenceType}
  fnc.hypotheses != nothing ? string(fnc.hypotheses) : ""
end

function prepgenericwrapper(
            Xi::Vector{Graphs.ExVertex},
            usrfnc::UnionAll,
            samplefnc::Function;
            multihypo::Union{Void, Distributions.Categorical}=nothing )
      # multiverts::Vector{Symbol}=Symbol[]
  #
  error("prepgenericwrapper -- unknown type usrfnc=$(usrfnc), maybe the wrong usrfnc conversion was dispatched.  Place an error in your unpacking convert function to ensure that IncrementalInference.jl is calling the right unpacking conversion function.")
end

function prepgenericwrapper(
            Xi::Vector{Graphs.ExVertex},
            usrfnc::T,
            samplefnc::Function;
            multihypo::Union{Void, Distributions.Categorical}=nothing ) where {T <: FunctorInferenceType}
      # multiverts::Vector{Symbol}=Symbol[]
  #
  warn("GenericWrapParam deprecated, use CommonConvWrapper instead.")
  ARR = Array{Array{Float64,2},1}()
  maxlen, sfidx = prepareparamsarray!(ARR, Xi, 0, 0)
  # test if specific zDim or partial constraint used
  fldnms = fieldnames(usrfnc)
  # sum(fldnms .== :zDim) >= 1
  gwp = GenericWrapParam{T}(
            usrfnc,
            ARR,
            1,
            1,
            (zeros(0,1),),
            samplefnc,
            sum(fldnms .== :zDim) >= 1,
            sum(fldnms .== :partial) >= 1,
            multihypo
        )
    gwp.factormetadata.variableuserdata = []
    gwp.factormetadata.solvefor = :null
    for xi in Xi
      push!(gwp.factormetadata.variableuserdata, getData(xi).softtype)
    end
    return gwp
end

# will be deprecated
function convert(
            ::Type{IncrementalInference.GenericFunctionNodeData{IncrementalInference.GenericWrapParam{F},Symbol}},
            d::IncrementalInference.GenericFunctionNodeData{P,String} ) where {F <: FunctorInferenceType, P <: PackedInferenceType}
  #
  warn("Packing and unpacking of GenericWrapParam{T} will be deprecated, use CommonConvWrapper{T} instead.")
  usrfnc = convert(F, d.fnc)
  # @show d.multihypo
  mhcat = parsemultihypostr(d.multihypo)
  # @show typeof(mhcat)
  gwpf = prepgenericwrapper(Graphs.ExVertex[], usrfnc, getSample, multihypo=mhcat)

  # ccw = prepgenericconvolution(Graphs.ExVertex[], usrfnc, multihypo=mhcat)
  return FunctionNodeData{GenericWrapParam{typeof(usrfnc)}}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          Symbol(d.frommodule), gwpf)
end

function approxConvOnElements!(frl::FastRootGenericWrapParam{T},
                               elements::Union{Vector{Int}, UnitRange{Int}}  ) where {T <: FunctorPairwise}
  #

  for n in elements
    frl.gwp.particleidx = n
    numericRootGenericRandomizedFnc!( frl )
  end
  nothing
end

function approxConvOnElements!(frl::FastRootGenericWrapParam{T},
                               elements::Union{Vector{Int}, UnitRange{Int}}) where {T <: FunctorPairwiseMinimize}

  # TODO should not be claiming new memory every single time....
  # frl.res = zeros(frl.xDim)
  # frl.gg = (x) -> frl.gwp(frl.res, x) # TODO standardize this function into frl

  # TODO -- once Threads.@threads have been optmized JuliaLang/julia#19967, also see area4 branch
  for n in elements
    frl.gwp.particleidx = n
    numericRootGenericRandomizedFnc!( frl )
    # frl.res[1:frl.xDim] = 0.0
    # r = optimize( frl.gg, frl.X[1:frl.xDim, frl.gwp.particleidx] )
    # # TODO -- clearly lots of optmization to be done here
    # frl.Y[1:frl.xDim] = r.minimizer
    # frl.X[1:frl.xDim,frl.gwp.particleidx] = frl.Y
  end
  nothing
end

function prepareFastRootGWP(gwp::GenericWrapParam{T},
                            Xi::Vector{Graphs.ExVertex},
                            solvefor::Int,
                            N::Int  ) where {T <: FunctorInferenceType}
  #

  ARR = Array{Array{Float64,2},1}()
  maxlen, sfidx = prepareparamsarray!(ARR, Xi, N, solvefor)
  # should be selecting for the correct multihypothesis mode here with `gwp.params=ARR[??]`
  gwp.params = ARR
  gwp.varidx = sfidx
  gwp.measurement = gwp.samplerfnc(gwp.usrfnc!, maxlen)
  zDim = size(gwp.measurement[1],1) # TODO -- zDim aspect desperately needs to be redone
  if gwp.specialzDim
    zDim = gwp.usrfnc!.zDim[sfidx]
  end
  # Construct complete fr (with fr.gwp) object
  # TODO -- create FastRootGenericWrapParam at addFactor time only?
  fr = FastRootGenericWrapParam{T}(gwp.params[sfidx], zDim, gwp)
  fr.res = zeros(fr.xDim)
  fr.gg = (x) -> fr.gwp(fr.res, x)
  return fr, sfidx, maxlen
end

function computeAcrossHypothesis(frl::FastRootGenericWrapParam{T},
                                 allelements,
                                 activehypo,
                                 certainidx,
                                 sfidx) where {T <:Union{FunctorPairwise, FunctorPairwiseMinimize}}
  #

  count = 0
  for (mhidx, vars) in activehypo
    count += 1
    if sfidx in certainidx || mhidx in certainidx # certainidx[count] in vars
      # standard case mhidx, sfidx = $mhidx, $sfidx
      frl.gwp.activehypo = vars
      approxConvOnElements!(frl, allelements[count])
    elseif mhidx == sfidx
      # multihypo, do conv case, mhidx == sfidx
      frl.gwp.activehypo = sort(union([sfidx;], certainidx))
      approxConvOnElements!(frl, allelements[count])
    elseif mhidx != sfidx
      # multihypo, take other value case
      # sfidx=2, mhidx=3:  2 should take a value from 3
      # sfidx=3, mhidx=2:  3 should take a value from 2
      frl.gwp.params[sfidx][:,allelements[count]] = view(frl.gwp.params[mhidx],:,allelements[count])
      # frl.gwp.params[sfidx][:,allelements[count]] = frl.gwp.params[mhidx][:,allelements[count]]
    else
      error("computeAcrossHypothesis -- not dealing with multi-hypothesis case correctly")
    end
  end
  nothing
end

function evalPotentialSpecific(Xi::Vector{Graphs.ExVertex},
                               gwp::GenericWrapParam{T},
                               solvefor::Int;
                               N::Int=100,
                               dbg::Bool=false ) where {T <: Union{FunctorPairwise, FunctorPairwiseMinimize}}
  #
  fnc = gwp.usrfnc!

  # Prep computation variables
  fr, sfidx, maxlen = prepareFastRootGWP(gwp, Xi, solvefor, N)
  certainidx, allelements, activehypo, mhidx = assembleHypothesesElements!(fr.gwp.hypotheses, maxlen, sfidx, length(Xi))

  # perform the numeric solutions on the indicated elements
  computeAcrossHypothesis(fr, allelements, activehypo, certainidx, sfidx)

  return fr.gwp.params[gwp.varidx]
end

function assembleNullHypothesis(fr::FastRootGenericWrapParam{T},
                                maxlen::Int,
                                spreadfactor::Float64 ) where {T}
  #

  nhc = rand(fr.gwp.usrfnc!.nullhypothesis, maxlen) - 1
  val = fr.gwp.params[fr.gwp.varidx]
  d = size(val,1)
  var = Base.var(val,2) + 1e-3
  ENT = Distributions.MvNormal(zeros(d), spreadfactor*diagm(var[:]))
  allelements = 1:maxlen
  return allelements, nhc, ENT
end

function computeAcrossNullHypothesis!(frl::FastRootGenericWrapParam{T},
                                      allelements,
                                      nhc,
                                      ENT  ) where {T <: FunctorPairwiseNH}
  #

  # TODO --  Threads.@threads see area4 branch
  for n in allelements
    # frl.gwp(x, res)
    if nhc[n] != 0
      frl.gwp.particleidx = n
      numericRootGenericRandomizedFnc!( frl )
    else
      frl.gwp.params[frl.gwp.varidx][:,n] += rand(ENT)
    end
  end
  nothing
end

function evalPotentialSpecific(Xi::Vector{Graphs.ExVertex},
                               generalwrapper::GenericWrapParam{T},
                               solvefor::Int;
                               N::Int=100,
                               spreadfactor::Float64=10.0,
                               dbg::Bool=false ) where {T <: FunctorSingletonNH}
  #

  fnc = generalwrapper.usrfnc!

  val = getVal(Xi[1])
  d = size(val,1)
  var = Base.var(val,2) + 1e-3

  # determine amount share of null hypothesis particles
  generalwrapper.measurement = generalwrapper.samplerfnc(generalwrapper.usrfnc!, N)
  # values of 0 imply null hypothesis
  # generalwrapper.usrfnc!.nullhypothesis::Distributions.Categorical
  nhc = rand(generalwrapper.usrfnc!.nullhypothesis, N) - 1

  # TODO -- not valid for manifold
  ENT = Distributions.MvNormal(zeros(d), spreadfactor*diagm(var[:]))

  for i in 1:N
    if nhc[i] == 0
      generalwrapper.measurement[1][:,i] = val[:,i] + rand(ENT)  # TODO use view and inplace add operation
    end
  end
  # TODO -- returning to memory location inside
  return generalwrapper.measurement[1]
end

function evalPotentialSpecific(Xi::Vector{Graphs.ExVertex},
                               generalwrapper::GenericWrapParam{T},
                               solvefor::Int;
                               N::Int=0,
                               dbg::Bool=false ) where {T <: FunctorSingleton}
  #

  fnc = generalwrapper.usrfnc!

  nn = (N <= 0 ? size(getVal(Xi[1]),2) : N)
  generalwrapper.measurement = generalwrapper.samplerfnc(generalwrapper.usrfnc!, nn)
  if !generalwrapper.partial
    return generalwrapper.measurement[1]
  else
    val = deepcopy(getVal(Xi[1]))
    i = 0
    for dimnum in fnc.partial
      i += 1
      val[dimnum,:] = generalwrapper.measurement[1][i,:]
    end
    return val
  end
end

function evalPotentialSpecific(Xi::Vector{Graphs.ExVertex},
                               gwp::GenericWrapParam{T},
                               solvefor::Int;
                               N::Int=100,
                               spreadfactor::Float64=10.0,
                               dbg::Bool=false ) where {T <: FunctorPairwiseNH}
  #

  # TODO -- could be constructed and maintained at addFactor! time
  fr, sfidx, maxlen = prepareFastRootGWP(gwp, Xi, solvefor, N)
  # prepare nullhypothesis
  allelements, nhc, ENT = assembleNullHypothesis(fr, maxlen, spreadfactor)

  # Compute across the true or null hypothesis
  computeAcrossNullHypothesis!(fr, allelements, nhc, ENT )

  return fr.gwp.params[gwp.varidx]
end

# Shuffle incoming X into random permutation in fr.Y
# shuffled fr.Y will be placed back into fr.X[:,fr.gwp.particleidx] upon fr.gwp.usrfnc(x, res)
function shuffleXAltD!(fr::FastRootGenericWrapParam, X::Vector{Float64})

  # populate defaults
  for i in 1:fr.xDim
    fr.Y[i] = fr.X[i, fr.gwp.particleidx]
  end
  # fr.Y[1:fr.xDim] = view(fr.X, 1:fr.xDim, fr.gwp.particleidx)
  # fr.Y[1:fr.xDim] = fr.X[1:fr.xDim,fr.gwp.particleidx]
  # copy!(fr.Y, fr.X[:,fr.gwp.particleidx])
  for i in 1:fr.zDim
    fr.Y[fr.p[i]] = X[i]
  end
  nothing
end

function (p::GenericWrapParam)(res::Vector{Float64}, x::Vector{Float64})

  # TODO -- move to inner lambda that is defined once against p.params...
  # approximates by not considering cross indices among parameters
  # @show length(p.params), p.varidx, p.particleidx, size(x), size(res), size(p.measurement)
  p.params[p.varidx][:, p.particleidx] = x
  # p.usrfnc!(res, p.particleidx, p.measurement, p.params...)
  # who are active hypotheses?  p.params[p.activehypo]...
  p.usrfnc!(res, p.factormetadata, p.particleidx, p.measurement, view(p.params,p.activehypo)...)
  # p.usrfnc!(res, p.factormetadata, p.particleidx, p.measurement, p.params[p.activehypo]...)
end

function (fr::FastRootGenericWrapParam)( res::Vector{Float64}, x::Vector{Float64} )

  shuffleXAltD!(fr, x)
  fr.gwp( res, fr.Y )
end

function numericRootGenericRandomizedFnc!(
            frl::FastRootGenericWrapParam{T};
            perturb::Float64=1e-10,
            testshuffle::Bool=false ) where {T <: FunctorPairwiseMinimize}
  #

  fill!(frl.res, 0.0) # 1:frl.xDim
  r = optimize( frl.gg, frl.X[:, frl.gwp.particleidx] ) # frl.gg
  # TODO -- clearly lots of optmization to be done here
  frl.Y[:] = r.minimizer
  frl.X[:,frl.gwp.particleidx] = frl.Y
  nothing
end

function numericRootGenericRandomizedFnc!(
            fr::FastRootGenericWrapParam{T};
            perturb::Float64=1e-10,
            testshuffle::Bool=false ) where {T <: FunctorPairwise}
  #
  # info("numericRootGenericRandomizedFnc! FastRootGenericWrapParam{$T}")
  # @show fr.zDim, fr.xDim, fr.gwp.partial, fr.gwp.particleidx

  if fr.zDim < fr.xDim && !fr.gwp.partial || testshuffle
    # less measurement dimensions than variable dimensions -- i.e. shuffle
    shuffle!(fr.p)
    for i in 1:fr.xDim
      fr.perturb[1:fr.zDim] = perturb*randn(fr.zDim)
      fr.X[fr.p[1:fr.zDim], fr.gwp.particleidx] += fr.perturb
      r = nlsolve(  fr,
                    fr.X[fr.p[1:fr.zDim], fr.gwp.particleidx] # this is x0
                 )
      if r.f_converged
        shuffleXAltD!( fr, r.zero )
        break;
      else
        # TODO -- report on this bottleneck, useful for optimization of code
        # @show i, fr.p, fr.xDim, fr.zDim
        temp = fr.p[end]
        fr.p[2:end] = fr.p[1:(end-1)]
        fr.p[1] = temp
        if i == fr.xDim
          error("numericRootGenericRandomizedFnc could not converge, i=$(i), fr.gwp.usrfnc!=$(typeof(fr.gwp.usrfnc!))")
        end
      end
    end
    #shuffleXAltD!( fr, r.zero ) # moved up
  elseif fr.zDim >= fr.xDim && !fr.gwp.partial
    # equal or more measurement dimensions than variable dimensions -- i.e. don't shuffle
    fr.perturb[1:fr.xDim] = perturb*randn(fr.xDim)
    fr.X[1:fr.xDim, fr.gwp.particleidx] += fr.perturb[1:fr.xDim] # moved up
    r = nlsolve( fr.gwp, fr.X[1:fr.xDim,fr.gwp.particleidx] )
    if sum(isnan.(( r ).zero)) == 0
      fr.Y[1:fr.xDim] = ( r ).zero
    else
      warn("got NaN, fr.gwp.particleidx = $(fr.gwp.particleidx), r=$(r)")
      @show fr.gwp.usrfnc!
      for thatlen in 1:length(fr.gwp.params)
        @show thatlen, fr.gwp.params[thatlen][:, fr.gwp.particleidx]
      end
    end
  elseif fr.gwp.partial
    # improve memory management in this function
    fr.p[1:length(fr.gwp.usrfnc!.partial)] = Int[fr.gwp.usrfnc!.partial...] # TODO -- move this line up and out of inner loop
    r = nlsolve(  fr,
                  fr.X[fr.p[1:fr.zDim], fr.gwp.particleidx] # this is x0
               )
    shuffleXAltD!( fr, r.zero )
  else
    error("Unresolved numeric solve case")
  end
  fr.X[:,fr.gwp.particleidx] = fr.Y
  nothing
end

# must set either dims or initval for proper initialization
# Add node to graph, given graph struct, labal, init values,
# std dev [TODO -- generalize], particle size and ready flag for concurrency
function addNode!(fg::FactorGraph,
      lbl::Symbol,
      initval::Array{Float64}=zeros(1,1),
      stdev::Array{Float64}=ones(1,1); # this is bad and should be removed TODO
      N::Int=100,
      ready::Int=1,
      labels::Vector{<: AbstractString}=String[],
      api::DataLayerAPI=dlapi,
      uid::Int=-1,
      dims::Int=-1  )
  #
  warn("this addNode! is deprecated, please use FactorGraph01.jl:addNode!(fg::FactorGraph, lbl::Symbol, softtype::Type{T}) instead.")
  currid = fg.id+1
  if uid==-1
    fg.id=currid
  else
    currid = uid
  end
  dims = dims != -1 ? dims : size(initval,1)

  lblstr = string(lbl)
  vert = ExVertex(currid,lblstr)
  addNewVarVertInGraph!(fg, vert, currid, lbl, ready, nothing)
  # dlapi.setupvertgraph!(fg, vert, currid, lbl) #fg.v[currid]
  dodims = fg.dimID+1
  # TODO -- vert should not loose information here
  setDefaultNodeData!(vert, initval, stdev, dodims, N, dims) #fg.v[currid]

  vnlbls = deepcopy(labels)
  push!(vnlbls, fg.sessionname)
  # addvert!(fg, vert, api=api)
  api.addvertex!(fg, vert, labels=vnlbls) #fg.g ##vertr =

  fg.dimID+=dims # rows indicate dimensions, move to last dimension
  push!(fg.nodeIDs, currid)

  return vert
end

function getVert(fgl::FactorGraph, lbl::T; api::DataLayerAPI=dlapi) where {T <: AbstractString}
  warn("IncrementalInference.getVert{T <: AbstractString}(fgl::FactorGraph, lbl::T; api::DataLayerAPI=dlapi) is deprecated, use lbl::Symbol instead")
  getVert(fgl, Symbol(lbl), api=api)
end

function getVal(fgl::FactorGraph, lbl::T; api::DataLayerAPI=dlapi) where {T <: AbstractString}
  warn("IncrementalInference.getVal{T <: AbstractString}(fgl::FactorGraph, lbl::T; api::DataLayerAPI=dlapi) is deprecated, use lbl::Symbol instead")
  getVal(fgl, Symbol(lbl),api=api)
end

function FNDencode{T <: FunctorInferenceType, P <: PackedInferenceType}(::Type{PackedFunctionNodeData{P}}, d::FunctionNodeData{T})
  warn("FNDencode deprecated, use the convert functions through dispatch instead, PackedFunctionNodeData{P=$(P)}.")
  return convert(PackedFunctionNodeData{P}, d) #PackedFunctionNodeData{P}
end
function FNDdecode{T <: FunctorInferenceType, P <: PackedInferenceType}(::Type{FunctionNodeData{T}}, d::PackedFunctionNodeData{P})
  warn("FNDdecode deprecated, use the convert functions through dispatch instead, FunctionNodeData{T=$(T)}.")
  return convert(FunctionNodeData{T}, d) #FunctionNodeData{T}
end

function FNDencode{T <: InferenceType, P <: PackedInferenceType}(::Type{PackedFunctionNodeData{P}}, d::FunctionNodeData{T})
  warn("FNDencode deprecated, use the convert functions through dispatch instead, PackedFunctionNodeData{P=$(P)}.")
  return convert(PackedFunctionNodeData{P}, d) #PackedFunctionNodeData{P}
end
function FNDdecode{T <: InferenceType, P <: PackedInferenceType}(::Type{FunctionNodeData{T}}, d::PackedFunctionNodeData{P})
  warn("FNDdecode deprecated, use the convert functions through dispatch instead, FunctionNodeData{T=$(T)}.")
  return convert(FunctionNodeData{T}, d) #FunctionNodeData{T}
end

# this will likely expand with more internal bells and whistles
# to perform in place memory operations for array values in

# see FastRootGenericWrapParam{T}

# TODO -- part of speed and flexibility refactoring exercise
# mutable struct FastGenericRoot{T} <: Function
#   p::Vector{Int}
#   perturb::Vector{Float64}
#   X::Vector{Float64}
#   Y::Vector{Float64}
#   xDim::Int
#   zDim::Int
#   usrfnc::T
#   FastGenericRoot{T}(xDim::Int, zDim::Int, residfnc::T) where {T} =
#       new(collect(1:xDim), zeros(zDim), zeros(xDim), zeros(xDim), xDim, zDim, residfnc)
# end
# function shuffleXAltD!(fr::FastGenericRoot, X::Vector{Float64})
#   copy!(fr.Y, fr.X)
#   for i in 1:fr.zDim
#     fr.Y[fr.p[i]] = X[i]
#   end
#   nothing
# end
# function (fr::FastGenericRoot)( res::Vector{Float64}, x::Vector{Float64} )
#   shuffleXAltD!(fr, x)  #(fr.Y, x, fr.x0, fr.zDim, fr.p)
#   fr.usrfnc( res, fr.Y )
# end
#
# """
#     $(SIGNATURES)
#
# DEPRECATED!
# Solve free variable x by root finding residual function fgr.usrfnc(x, res)
# randomly shuffle x dimensions if underconstrained by measurement z dimensions
# small random perturbation used to prevent trivial solver cases, div by 0 etc.
# result stored in fgr.Y
# """
# function numericRootGenericRandomizedFnc!(
#       fgr::FastGenericRoot{T};
#       perturb::Float64=1e-5,
#       testshuffle::Bool=false ) where {T}
#   #
#   warn("numericRootGenericRandomizedFnc!(fgr::FastGenericRoot{T}...) deprecated, use numericRootGenericRandomizedFnc!(fgr::FastRootGenericWrapRoot{T}...) instead.")
#   fgr.perturb[1:fgr.zDim] = perturb*randn(fgr.zDim)
#   if fgr.zDim < fgr.xDim || testshuffle
#     shuffle!(fgr.p)
#     r = nlsolve(  fgr,
#                 fgr.X[fgr.p[1:fgr.zDim]] + fgr.perturb # this is x0
#              )
#     # copy!(fgr.X, x0) # should need this line?
#     shuffleXAltD!( fgr, r.zero )
#     # copy!(fgr.X, r.zero)
#   else
#     # @show "direct solve"
#     fgr.Y = ( nlsolve(  fgr.usrfnc, fgr.X + fgr.perturb ) ).zero
#     # copy!(fgr.X, y)
#   end
#   nothing
# end
#
#
# function numericRootGenericRandomizedFnc(
#       residFnc!::Function,
#       zDim::Int,
#       xDim::Int,
#       x0::Vector{Float64};
#       perturb::Float64=1e-5,
#       testshuffle::Bool=false   )
#   #
#   # TODO -- this only start of refactoring for inplace, more to come
#   # xDim = length(x0)
#   fgr = FastGenericRoot{typeof(residFnc!)}(xDim, zDim, residFnc!)
#   shuffle!(fgr.p);
#   fgr.perturb[1:fgr.zDim] = perturb*randn(fgr.zDim)
#   copy!(fgr.X, x0)
#
#   numericRootGenericRandomizedFnc!( fgr, perturb=perturb, testshuffle=testshuffle )
#   fgr.Y
# end

# function evalPotentialSpecific(
#       fnc::T,
#       Xi::Vector{Graphs.ExVertex},
#       gwp::GenericWrapParam{T},
#       solvefor::Int;
#       N::Int=100  ) where {T <: FunctorPairwiseMinimize}
#   #
#   # TODO -- could be constructed and maintained at addFactor! time
#   fr, sfidx, maxlen = prepareFastRootGWP(gwp, Xi, solvefor, N)
#   certainidx, allelements, activehypo, mhidx = assembleHypothesesElements!(fr.gwp.hypotheses, maxlen, sfidx, length(Xi))
#
#   # perform the numeric solutions on the indicated elements
#   computeAcrossHypothesis(T, fr, allelements, activehypo, certainidx, sfidx)
#
#   return fr.gwp.params[gwp.varidx]
# end

# function computeAcrossHypothesis(T::Type{<:FunctorPairwiseMinimize}, fr, allelements, activehypo, certainidx, sfidx)
#   count = 0
#   for (mhidx, vars) in activehypo
#     count += 1
#     # if length(allelements[count]) > 0
#     #   fr.gwp.activehypo = vars
#     #   approxConvMinimizeOnElements!(fr, allelements[count])
#     # end
#     if sfidx in certainidx || mhidx in certainidx # certainidx[count] in vars
#       # standard case mhidx, sfidx = $mhidx, $sfidx
#       fr.gwp.activehypo = vars
#       approxConvOnElements!(fr, allelements[count])
#     elseif mhidx == sfidx
#       # multihypo, do conv case, mhidx == sfidx
#       fr.gwp.activehypo = sort(union([sfidx;], certainidx))
#       approxConvOnElements!(fr, allelements[count])
#     elseif mhidx != sfidx
#       # multihypo, take other value case
#       # sfidx=2, mhidx=3:  2 should take a value from 3
#       # sfidx=3, mhidx=2:  3 should take a value from 2
#       fr.gwp.params[sfidx][:,allelements[count]] = fr.gwp.params[mhidx][:,allelements[count]]
#     else
#       error("computeAcrossHypothesis -- not dealing with multi-hypothesis case correctly")
#     end
#   end
#   nothing
# end


text:
#=
This is an experiment to see if one player starts out holding just four aces
and the other player has the rest of the deck, how do they fare?
=#

using PlayingCards52, War, ProgressMeter

"""
    ace_race(reps::Int = 10_000)

returns the frequency that four aces vs the rest of the deck wins at war.
Optionally, `ace_race(reps, skip_rank::Int)` to skip a different rank, e.g.,
`ace_race(skip_rank=13)` for kings.
"""
function ace_race(reps::Int = 10_000; skip_rank = 1)

    d = deck(false)
    d2 = [c for c in d if rank(c) != skip_rank]
    d1 = [Card(s, skip_rank) for s in PlayingCards52.suit_list]

    a_wins = 0
    PM = Progress(reps)

    for k = 1:reps
        shuffle!(d2)
        r = play_war(d1, d2)
        if r[end] > 0
            a_wins += 1
        end
        next!(PM)
    end

    a_pct = a_wins / reps
    return a_pct

    # a_len = a_steps / reps
    # b_len = b_steps / reps

    # println("Four aces wins $a_pct% of the time")
    # println("Average length of a winning game: $a_len")
    # println("Average length of a losing game:  $b_len")

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct VirtualMachineRunCommandsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets specific run command for a subscription in a location.
Param: location::String (required)
Param: commandId::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: RunCommandDocument
"""
function virtualMachineRunCommandsGet(_api::VirtualMachineRunCommandsApi, location::String, commandId::String, api_version::String, subscriptionId::String; _mediaType=nothing)

    _ctx = Swagger.Ctx(_api.client, "GET", RunCommandDocument, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/runCommands/{commandId}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "commandId", commandId)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Lists all available run commands for a subscription in a location.
Param: location::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: RunCommandListResult
"""
function virtualMachineRunCommandsList(_api::VirtualMachineRunCommandsApi, location::String, api_version::String, subscriptionId::String; _mediaType=nothing)

    _ctx = Swagger.Ctx(_api.client, "GET", RunCommandListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/runCommands", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export virtualMachineRunCommandsGet, virtualMachineRunCommandsList


text:
#
# Calculates numerically:∫ N(t⋅x | μ, Σ) dt
# (see also function analyticalintegraloft)
#
function numericalintegraloft(x, μ, Σ)

    g = MvNormal(μ, Σ)

    hquadrature(t -> pdf(g, t*x), -150, 150)[1]

end

#
# Calculates the integral: ∫ N(t⋅x | μ, Σ) dt
# (see also function loganalyticalintegraloft)
#
function analyticalintegraloft(x, μ, Σ)

    # see https://en.wikipedia.org/wiki/Gaussian_integral#Generalizations

    D = length(μ) ; @assert(size(Σ, 1) == size(Σ, 2) == D == length(x))

    a =  0.5*x'*(Σ\x)

    b =  x'*(Σ\μ)

    c = -0.5*μ'*(Σ\μ)

    return 1/sqrt((2π)^D*det(Σ)) * sqrt(π/a) * exp(b^2/(4*a)+c)

end

#
# Calculates: log ∫ N(t⋅x | μ, Σ) dt
# (see also function analyticalintegraloft)
#
function loganalyticalintegraloft(x, μ, Σ)

    # see https://en.wikipedia.org/wiki/Gaussian_integral#Generalizations

    D = length(μ) ; @assert(size(Σ, 1) == size(Σ, 2) == D == length(x))

    a =  0.5*x'*(Σ\x)

    b =  x'*(Σ\μ)

    c = -0.5*μ'*(Σ\μ)

    return - 0.5*D*log(2π) - 0.5*logdet(Σ) + log(sqrt(π/a)) + (b^2/(4*a)+c)

end

function difficultintegral4expectation(post, x, S)

    rg = MersenneTwister(1)

    bcandidates = [rand(rg, post)[4:6] for _ in 1:S]

    μ_a_given_b, Σ_a_given_b = gaussconditional(mean(post), cov(post), [1;2;3], [4;5;6])

    sub(b) = analyticalintegraloft(x-b, μ_a_given_b(b), Σ_a_given_b)

    return mean(map(sub, bcandidates))

end

function difficultintegral5(post, numsamples, Trange, Srange)

    # Fix random seed
    rg = MersenneTwister(1)

    # Draw samples for b
    bcandidates = [rand(rg, post)[4:6] for _ in 1:numsamples]

    # Calculate conditional Gaussian
    μ_a_given_b, Σ_a_given_b = gaussconditional(mean(post), cov(post), [1;2;3], [4;5;6])

    aux(x, b) = analyticalintegraloft(x-b, μ_a_given_b(b), Σ_a_given_b)

    # Store here results
    TS = zeros(length(Trange), length(Srange))

    # Report progress
    pr = Progress(length(Trange))

    # Loop over temperatures
    for i in 1:length(Trange)

        # generate galaxy vector
        obsgal = observedunitgalaxyvector(Trange[i])

        for j in 1:length(Srange)

            # scale galaxy vector
            x = obsgal * Srange[j]

            # calculate density
            TS[i, j] = mean(map(b -> aux(x, b), bcandidates))

        end

        # Update progress bar

        bestTsofar = Trange[argmax(vec(sum(TS,dims=2)))]

        ProgressMeter.next!(pr; showvalues = [(:Current_Temperature, Trange[i]), (:Best_so_far, bestTsofar)])

    end

    return TS

end

function difficultintegral6(post, numsamples, Trange)

    # Fix random seed
    rg = MersenneTwister(1)

    # Draw samples for b
    bcandidates = [rand(rg, post)[4:6] for _ in 1:numsamples]

    # Calculate conditional Gaussian
    μ_a_given_b, Σ_a_given_b = gaussconditional(mean(post), cov(post), [1;2;3], [4;5;6])

    aux(x, b) = loganalyticalintegraloft(x-b, μ_a_given_b(b), Σ_a_given_b)

    # Store here results
    TS = zeros(length(Trange))

    # Report progress
    pr = Progress(length(Trange))

    # Loop over temperatures
    Threads.@threads for i in 1:length(Trange)

        # generate galaxy vector
        obsgal = observedunitgalaxyvector(Trange[i])

        # define objective, remember we are minimising
        f = s -> -1.0 * logsumexp(map(b -> aux(obsgal * s, b), bcandidates))

        # approximate integral with value at mode
        TS[i] = exp(-1.0 * optimize(f, 0.001, 200).minimum)

        # Update progress bar
        bestTsofar = Trange[argmax(vec(sum(TS, dims=2)))]

        ProgressMeter.next!(pr; showvalues = [(:Best_so_far, bestTsofar)])

    end

    return TS

end

function difficultintegralVI(x, posterior, S, maxiter)

    function logp(a, b)
         logpdf(posterior, [a;b])
     end

    aux(param) = logp(param[1:3], param[4:end])

    μVI, ΣVI, negVI = VI(aux, randn(6), optimiser=LBFGS(), S=S, maxiter=maxiter)

    return μVI, ΣVI, negVI

end


text:
@doc raw"""
`ParametersDGVI`: Parameters for right-hand side function of Discontinuous Galerkin Variational Integrator.

### Parameters

* `DT`: data type
* `TT`: parameter type
* `D`: dimension of the system
* `S`: number of basis nodes
* `R`: number of quadrature nodes

### Fields

* `Θ`:  function of the noncanonical one-form (∂L/∂v)
* `f`:  function of the force (∂L/∂q)
* `g`:  function of the projection ∇ϑ(q)⋅v
* `Δt`: time step
* `b`:  quadrature weights
* `c`:  quadrature nodes
* `m`:  mass matrix
* `a`:  derivative matrix
* `r⁻`: reconstruction coefficients, jump lhs value
* `r⁺`: reconstruction coefficients, jump rhs value
* `t`:  initial time
* `q`:  solution of q  at time t
* `q⁻`: solution of q⁻ at time t
* `q⁺`: solution of q⁺ at time t
"""
mutable struct ParametersDGVI{DT, TT, D, S, R, ET <: NamedTuple} <: Parameters{DT,TT}
    equs::ET
    Δt::TT

    b::Vector{TT}
    c::Vector{TT}
    m::Matrix{TT}
    a::Matrix{TT}
    r⁻::Vector{TT}
    r⁺::Vector{TT}

    t::TT
    q::Vector{DT}
    p::Vector{DT}

    q⁻::Vector{DT}
    q⁺::Vector{DT}

    function ParametersDGVI{DT,D}(equs::ET, Δt::TT, b, c, m, a, r⁻, r⁺) where {DT, TT, D, ET <: NamedTuple}
        @assert length(b)  == length(c)
        @assert length(r⁻) == length(r⁺)
        new{DT,TT,D,length(r⁻),length(c),ET}(equs, Δt, b, c, m, a, r⁻, r⁺, zero(TT), zeros(DT,D), zeros(DT,D), zeros(DT,D), zeros(DT,D))
    end

    function ParametersDGVI{DT,D}(equs::NamedTuple, Δt::TT, basis::Basis{TT}, quadrature::QuadratureRule{TT}) where {DT,TT,D}
        # compute coefficients
        m = zeros(TT, nnodes(quadrature), nbasis(basis))
        a = zeros(TT, nnodes(quadrature), nbasis(basis))
        r⁻= zeros(TT, nbasis(basis))
        r⁺= zeros(TT, nbasis(basis))

        for i in eachindex(basis)
            for j in eachindex(quadrature)
                m[j,i] = basis[nodes(quadrature)[j], i]
                a[j,i] = basis'[nodes(quadrature)[j], i]
            end
            r⁻[i] = basis[one(TT), i]
            r⁺[i] = basis[zero(TT), i]
        end

        if get_config(:verbosity) > 1
            println()
            println("  Discontinuous Galerkin Variational Integrator")
            println("  =============================================")
            println()
            println("    b = ", b)
            println("    c = ", c)
            println("    m = ", m)
            println("    a = ", a)
            println("    r⁻= ", r⁻)
            println("    r⁺= ", r⁺)
            println()
        end

        ParametersDGVI{DT,D}(equs, Δt, weights(quadrature), nodes(quadrature), m, a, r⁻, r⁺)
    end
end

function update_params!(params::ParametersDGVI, sol::AtomicSolutionPODE)
    # set time for nonlinear solver and copy previous solution
    params.t  = sol.t
    params.q .= sol.q
    params.p .= sol.p
end

@doc raw"""
Nonlinear function cache for Discontinuous Galerkin Variational Integrator.

### Parameters

* `ST`: data type
* `D`: number of dimensions
* `S`: number of degrees of freedom
* `R`: number of nodes of quadrature formula

### Fields

* `X`: degrees of freedom
* `Q`: solution at quadrature nodes
* `V`: velocity at quadrature nodes
* `P`: one-form at quadrature nodes
* `F`: forces at quadrature nodes
* `q`:  current solution of ``q_{n}``
* `q⁻`: current solution of ``q_{n}^{-}``
* `q⁺`: current solution of ``q_{n}^{+}``
* `q̄`:  current solution of ``q_{n+1}``
* `q̄⁻`: current solution of ``q_{n+1}^{-}``
* `q̄⁺`: current solution of ``q_{n+1}^{+}``
* `ϕ`:  average of the solution at ``t_{n}``
* `ϕ̅`:  average of the solution at ``t_{n+1}``
* `λ`:  jump of the solution at ``t_{n}``
* `λ̄`:  jump of the solution at ``t_{n+1}``
* `θ`:  one-form evaluated across at ``t_{n}``
* `Θ̅`:  one-form evaluated across at ``t_{n+1}``
* `g`:  projection evaluated across at ``t_{n}``
* `ḡ`:  projection evaluated across at ``t_{n+1}``
"""
struct IntegratorCacheDGVI{ST,D,S,R} <: IODEIntegratorCache{ST,D}
    X::Vector{Vector{ST}}
    Q::Vector{Vector{ST}}
    V::Vector{Vector{ST}}
    P::Vector{Vector{ST}}
    F::Vector{Vector{ST}}

    q̃::Vector{ST}
    p̃::Vector{ST}
    ṽ::Vector{ST}
    f̃::Vector{ST}
    s̃::Vector{ST}

    q::Vector{ST}
    q⁻::Vector{ST}
    q⁺::Vector{ST}
    q̄::Vector{ST}
    q̄⁻::Vector{ST}
    q̄⁺::Vector{ST}

    p::Vector{ST}
    p̄::Vector{ST}

    ϕ::Vector{ST}
    ϕ⁻::Vector{ST}
    ϕ⁺::Vector{ST}
    ϕ̅::Vector{ST}
    ϕ̅⁻::Vector{ST}
    ϕ̅⁺::Vector{ST}

    λ::Vector{ST}
    λ⁻::Vector{ST}
    λ⁺::Vector{ST}
    λ̄::Vector{ST}
    λ̄⁻::Vector{ST}
    λ̄⁺::Vector{ST}

    θ::Vector{ST}
    θ⁻::Vector{ST}
    θ⁺::Vector{ST}
    Θ̅::Vector{ST}
    Θ̅⁻::Vector{ST}
    Θ̅⁺::Vector{ST}

    g::Vector{ST}
    g⁻::Vector{ST}
    g⁺::Vector{ST}
    ḡ::Vector{ST}
    ḡ⁻::Vector{ST}
    ḡ⁺::Vector{ST}

    h⁺::Vector{ST}
    h̅⁻::Vector{ST}

    function IntegratorCacheDGVI{ST,D,S,R}() where {ST,D,S,R}
        # create internal stage vectors
        X = create_internal_stage_vector(ST,D,S)
        Q = create_internal_stage_vector(ST,D,R)
        V = create_internal_stage_vector(ST,D,R)
        P = create_internal_stage_vector(ST,D,R)
        F = create_internal_stage_vector(ST,D,R)

        # create temporary vectors
        q̃ = zeros(ST,D)
        p̃ = zeros(ST,D)
        ṽ = zeros(ST,D)
        f̃ = zeros(ST,D)
        s̃ = zeros(ST,D)

        # create solution vectors
        q  = zeros(ST,D)
        q⁻ = zeros(ST,D)
        q⁺ = zeros(ST,D)
        q̄  = zeros(ST,D)
        q̄⁻ = zeros(ST,D)
        q̄⁺ = zeros(ST,D)

        p = zeros(ST,D)
        p̄ = zeros(ST,D)

        # create jump vectors
        ϕ  = zeros(ST,D)
        ϕ⁻ = zeros(ST,D)
        ϕ⁺ = zeros(ST,D)
        ϕ̅  = zeros(ST,D)
        ϕ̅⁻ = zeros(ST,D)
        ϕ̅⁺ = zeros(ST,D)

        λ  = zeros(ST,D)
        λ⁻ = zeros(ST,D)
        λ⁺ = zeros(ST,D)
        λ̄ = zeros(ST,D)
        λ̄⁻ = zeros(ST,D)
        λ̄⁺ = zeros(ST,D)

        θ  = zeros(ST,D)
        θ⁻ = zeros(ST,D)
        θ⁺ = zeros(ST,D)
        Θ̅  = zeros(ST,D)
        Θ̅⁻ = zeros(ST,D)
        Θ̅⁺ = zeros(ST,D)

        g  = zeros(ST,D)
        g⁻ = zeros(ST,D)
        g⁺ = zeros(ST,D)
        ḡ  = zeros(ST,D)
        ḡ⁻ = zeros(ST,D)
        ḡ⁺ = zeros(ST,D)

        h⁺ = zeros(ST,D)
        h̅⁻ = zeros(ST,D)

        new(X, Q, V, P, F,
            q̃, p̃, ṽ, f̃, s̃,
            q, q⁻, q⁺, q̄, q̄⁻, q̄⁺,
            p, p̄,
            ϕ, ϕ⁻, ϕ⁺, ϕ̅, ϕ̅⁻, ϕ̅⁺,
            λ, λ⁻, λ⁺, λ̄, λ̄⁻, λ̄⁺,
            θ, θ⁻, θ⁺, Θ̅, Θ̅⁻, Θ̅⁺,
            g, g⁻, g⁺, ḡ, ḡ⁻, ḡ⁺,
            h⁺, h̅⁻)
    end
end

function IntegratorCache{ST}(params::ParametersDGVI{DT,TT,D,S,R}; kwargs...) where {ST,DT,TT,D,S,R}
    IntegratorCacheDGVI{ST,D,S,R}(; kwargs...)
end

@inline CacheType(ST, params::ParametersDGVI{DT,TT,D,S,R}) where {DT,TT,D,S,R} = IntegratorCacheDGVI{ST,D,S,R}

@doc raw"""
`IntegratorDGVI`: Discontinuous Galerkin Variational Integrator.

The DGVI integrators arise from the discretization of the action integral
```math
\mathcal{A} [q] = \int \limits_{0}^{T} L(q(t), \dot{q}(t)) \, dt ,
```
with ``L`` a fully degenerate Lagrangian of the form
```math
L(q, \dot{q}) = \vartheta (q) \cdot \dot{q} - H(q) ,
```
where ``\vartheta (q)`` denotes the Cartan one-form and ``H(q)`` the Hamiltonian,
which is usually given by the total energy of the system.

### Discretization

Within each interval ``(t_{n}, t_{n+1})`` a piecewise-polynomial approximation ``q_h``
of the trajectory ``q`` is constructed using ``S`` basis functions ``\varphi_{i}``,
```math
q_h(t) \vert_{(t_{n}, t_{n+1})} = \sum \limits_{i=1}^{S} x_{n,i} \, \bar{\varphi}_{n,i} (t) ,
```
where ``\bar{\varphi}_{n,i} (t)`` is a rescaled basis function, defined by
```math
\bar{\varphi}_{n,i} (t) = \varphi_{i} \bigg( \frac{t - t_{n}}{t_{n+1} - t_{n}} \bigg) ,
```
and it is assumed that ``\varphi_{i}`` is compactly supported on ``[0,1]``.
These approximations ``q_h(t)`` are not assumed to be continuous across interval
boundaries ``t_{n}`` but usuaslly have jumps.

The integral over ``(t_{n}, t_{n+1})`` is approximated by a quadrature rule with
``R`` nodes ``c_i`` and weights ``b_i``.
Denote by ``m`` and ``a`` mass and derivative matrices, respectively, whose elements
 are given by
```math
m_{ij} = \varphi_j (c_i) ,
\qquad
a_{ij} = \varphi_j' (c_i) ,
\qquad
i = 1, ..., R ,
\;
j = 1, ..., S .
```
With that, the solution and its time derivative at the quadrature points can be written as
```math
Q_{n,i} \equiv q_h(t_n + c_i h) = m_{ij} x_{n,j} ,
\qquad
V_{n,i} \equiv \dot{q}_h (t_n + c_i h) = a_{ij} x_{n,j} ,
```
where
```math
x_{n} = ( x_{n,1}, ..., x_{n,S} )^T
```
is the vector containing the degrees of freedom of ``q_h \vert_{[t_{n}, t_{n+1}]}``.
The limits of ``q_h(t)`` at ``t_{n}`` and ``t_{n+1}`` are given by
```math
q_{n}^{+} = \lim \limits_{t \downarrow t_{n}} q_h(t) = \sum \limits_{j=1}^{S} r^{+}_{j} \delta x_{n,j} ,
\hspace{3em}
q_{n+1}^{-} = \lim \limits_{t \uparrow t_{n+1}} q_h(t) = \sum \limits_{j=1}^{S} r^{-}_{j} \delta x_{n,j} .
```

The discrete action reads
```math
\mathcal{A}_d [x_d] = h \sum \limits_{n=0}^{N-1} \bigg[
     \sum \limits_{i=1}^{R} b_i \big[ \vartheta (Q_{n,i}) \cdot V_{n,i} - H(Q_{n,i}) \big] \\
     + \frac{\vartheta (q_n) + \vartheta (q_n^+)}{2} \cdot (q_n^+ - q_n)
     + \frac{\vartheta (q_{n+1}^-) + \vartheta (q_{n+1})}{2} \cdot (q_{n+1} - q_{n+1}^-)
\bigg] ,
```
so that using the relations
```math
\delta Q_{n,i} = m_{ij} \delta x_{n,j} ,
\qquad
\delta V_{n,i} = \frac{a_{ij}}{h} \delta x_{n,j} ,
\qquad
\delta q_{n}^- = r^{-}_{j} \delta x_{n-1,j} ,
\qquad
\delta q_{n}^+ = r^{+}_{j} \delta x_{n,j} ,
```
the discrete action principle leads to the discrete equations of motion,
```math
0 = \sum \limits_{i=1}^{R} b_i \big[ h m_{ij} \nabla \vartheta (Q_{n,i}) \cdot V_{n,i} + a_{ij} \vartheta (Q_{n,i}) - h m_{ij} \nabla H(Q_{n,i}) \big] \\
+ r^{+}_{j} \, \frac{ \vartheta ( q_{n  }   ) + \vartheta( q_{n  }^+ ) }{2}
- r^{-}_{j} \, \frac{ \vartheta ( q_{n+1}^- ) + \vartheta( q_{n+1}   ) }{2}
   \big] \\
+ h r^{+}_{j} \, \nabla \vartheta (q_{n  }^+) \cdot (q_{n  }^+ - q_{n  }  )
+ h r^{-}_{j} \, \nabla \vartheta (q_{n+1}^-) \cdot (q_{n+1}   - q_{n+1}^-) ,
```
and
```math
\vartheta(q_{n}^+) = \vartheta (q_{n}^-) + \nabla \vartheta (q_{n}) \cdot (q_{n}^+ - q_{n}^-) ,
```
for all ``n`` and all ``j``.
Let us introduce the variable ``p_n`` as
```math
p_{n} = \vartheta (q_{n}^-) + \nabla \vartheta (q_{n}) \cdot (q_{n} - q_{n}^-) ,
```
so that
```math
\vartheta(q_{n}^+) = p_{n} + \nabla \vartheta (q_{n}) \cdot (q_{n}^+ - q_{n}) .
```
Then the above equations provide a map ``(q_{n}, p_{n}) \mapsto (q_{n+1}, p_{n+1})``.
In order to solve these equations, initial conditions ``q_{0}`` and
``p_{0} = \vartheta(q_{0})`` have to be prescribed.

### Fields

* `equation`: Implicit Ordinary Differential Equation
* `basis`: piecewise polynomial basis
* `quadrature`: numerical quadrature rule
* `Δt`: time step
* `params`: ParametersDGVI
* `solver`: nonlinear solver
* `iguess`: initial guess
* `q`: current solution vector for trajectory
* `p`: current solution vector for one-form
* `cache`: temporary variables for nonlinear solver
"""
struct IntegratorDGVI{DT, TT, D, S, R,
                BT <: Basis,
                PT <: ParametersDGVI{DT,TT,D,S},
                ST <: NonlinearSolver{DT},
                IT <: InitialGuessODE{TT}} <: IODEIntegrator{DT,TT}
    basis::BT
    quadrature::QuadratureRule{TT,R}

    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorDGVI(basis::BT, quadrature::QuadratureRule{TT,R}, params::ParametersDGVI{DT,TT,D,S},
                    solver::ST, iguess::IT, caches) where {DT,TT,D,S,R,BT,ST,IT}
        new{DT, TT, D, S, R, BT, typeof(params), ST, IT}(basis, quadrature, params, solver, iguess, caches)
    end

    function IntegratorDGVI{DT,D}(equations::NamedTuple, basis::Basis{TT}, quadrature::QuadratureRule{TT,R}, Δt::TT;
                                  interpolation=HermiteExtrapolation{DT}) where {DT,TT,D,R}

        # get number of stages
        S = nbasis(basis)

        # create params
        params = ParametersDGVI{DT,D}(equations, Δt, basis, quadrature)

        # create cache dict
        caches = CacheDict(params)

        # create nonlinear solver
        solver = create_nonlinear_solver(DT, D*(S+1), params, caches)

        # create initial guess
        iguess = InitialGuessODE(get_config(:ig_extrapolation), equations[:v̄], Δt)

        # create integrator
        IntegratorDGVI(basis, quadrature, params, solver, iguess, caches)
    end

    function IntegratorDGVI(equation::Union{IODE{DT}, LODE{DT}}, basis::Basis, quadrature::QuadratureRule, Δt; kwargs...) where {DT}
        IntegratorDGVI{DT, ndims(equation)}(get_functions(equation), basis, quadrature, Δt; kwargs...)
    end
end

@inline equation(integrator::IntegratorDGVI, i::Symbol) = integrator.params.equs[i]
@inline equations(integrator::IntegratorDGVI) = integrator.params.equs
@inline timestep(integrator::IntegratorDGVI) = integrator.params.Δt
@inline Base.ndims(::IntegratorDGVI{DT,TT,D}) where {DT,TT,D} = D

function IntegratorCache(int::IntegratorDGVI{DT,TT}) where {DT,TT}
    IntegratorCacheDGVI{DT, TT, ndims(int), nbasis(int.basis), nnodes(int.quadrature)}()
end

function initialize!(int::IntegratorDGVI, sol::AtomicSolutionPODE)
    sol.t̄ = sol.t - timestep(int)

    equation(int, :v̄)(sol.t, sol.q, sol.v)

    initialize!(int.iguess, sol.t, sol.q, sol.v,
                            sol.t̄, sol.q̄, sol.v̄)
end

function initial_guess!(int::IntegratorDGVI{DT,TT,D,S,R}, sol::AtomicSolutionPODE{DT,TT},
                        cache::IntegratorCacheDGVI{DT}=int.caches[DT]) where {DT,TT,D,S,R}
    if nbasis(int.basis) > 0
        for i in eachindex(int.basis)
            evaluate!(int.iguess, sol.q̄, sol.v̄,
                                  sol.q, sol.v,
                                  cache.q̃,
                                  grid(int.basis)[i])

            for k in 1:D
                int.solver.x[D*(i-1)+k] = cache.q̃[k]
            end
        end
    else
        for i in 1:S
            for k in 1:D
                int.solver.x[D*(i-1)+k] = 0
            end
        end
    end

    evaluate!(int.iguess, sol.q̄, sol.v̄,
                          sol.q, sol.v,
                          cache.q̃,
                          one(TT))

    for k in 1:D
        int.solver.x[D*S+k] = cache.q̃[k]
    end
end

"Compute stages of variational partitioned Runge-Kutta methods."
function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersDGVI{DT,TT,D,S,R},
                caches::CacheDict) where {ST,DT,TT,D,S,R}
    @assert length(x) == length(b)

    # get cache for internal stages
    cache = caches[ST]

    # compute stages from nonlinear solver solution x
    compute_stages!(x, cache, params)

    # compute rhs b of nonlinear solver
    compute_rhs!(b, cache, params)
end

function compute_stages!(x, cache::IntegratorCacheDGVI{ST,D,S}, params::ParametersDGVI{DT,TT,D,S}) where {ST,DT,TT,D,S}
    # copy x to X
    for i in 1:S
        for k in 1:D
            cache.X[i][k] = x[D*(i-1)+k]
        end
    end

    # copy x to q̄=qₙ+₁
    for k in 1:D
        cache.q̄[k] = x[D*S+k]
    end

    # compute Q, qₙ⁺ and qₙ₊₁⁻
    compute_stages_q!(cache, params)

    # compute V
    compute_stages_v!(cache, params)

    # compute P and F
    compute_stages_p!(cache, params)

    # compute jump
    compute_stages_λ!(cache, params)
end

"Compute solution at quadrature nodes and across jump."
function compute_stages_q!(cache::IntegratorCacheDGVI{ST,D,S,R},
                           params::ParametersDGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    local q::ST
    local q⁺::ST
    local q̄⁻::ST

    local X = cache.X
    local Q = cache.Q

    # copy q and q⁻
    cache.q  .= params.q
    cache.p  .= params.p
    cache.q⁻ .= params.q⁻

    # compute Q
    for i in 1:R
        for k in 1:D
            q = 0
            for j in 1:S
                q += params.m[i,j] * X[j][k]
            end
            Q[i][k] = q
        end
    end

    # compute qₙ⁺ and qₙ₊₁⁻
    for k in 1:D
        q⁺ = 0
        q̄⁻ = 0
        for i in 1:S
            q⁺ += params.r⁺[i] * X[i][k]
            q̄⁻ += params.r⁻[i] * X[i][k]
        end
        cache.q⁺[k] = q⁺
        cache.q̄⁻[k] = q̄⁻
    end
end

"Compute velocities at quadrature nodes."
function compute_stages_v!(cache::IntegratorCacheDGVI{ST,D,S,R},
                           params::ParametersDGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local v::ST

    for i in 1:R
        for k in 1:D
            v = 0
            for j in 1:S
                v += params.a[i,j] * cache.X[j][k]
            end
            cache.V[i][k] = v / params.Δt
        end
    end
end

"Compute one-form and forces at quadrature nodes."
function compute_stages_p!(cache::IntegratorCacheDGVI{ST,D,S,R},
                           params::ParametersDGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    local tᵢ::TT

    # compute P=ϑ(Q) and F=f(Q)
    for i in 1:R
        tᵢ = params.t + params.Δt * params.c[i]
        params.equs[:ϑ](tᵢ, cache.Q[i], cache.V[i], cache.P[i])
        params.equs[:f](tᵢ, cache.Q[i], cache.V[i], cache.F[i])
    end
end

function compute_stages_λ!(cache::IntegratorCacheDGVI{ST,D,S,R},
                           params::ParametersDGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    local t₀::TT = params.t
    local t₁::TT = params.t + params.Δt

    # compute ϕ and ϕ̅
    cache.ϕ  .= 0.5 * (cache.q⁺ .+ cache.q⁻)
    cache.ϕ⁻ .= 0.5 * (cache.q⁻ .+ cache.q )
    cache.ϕ⁺ .= 0.5 * (cache.q  .+ cache.q⁺)
    cache.ϕ̅⁻ .= 0.5 * (cache.q̄⁻ .+ cache.q̄ )

    # compute λ and λ̄
    cache.λ  .= cache.q⁺ .- cache.q⁻
    cache.λ⁻ .= cache.q  .- cache.q⁻
    cache.λ⁺ .= cache.q⁺ .- cache.q
    cache.λ̄⁻ .= cache.q̄  .- cache.q̄⁻

    # compute ϑ
    params.equs[:ϑ](t₀, cache.q,  cache.q,  cache.θ)
    params.equs[:ϑ](t₀, cache.q⁻, cache.q⁻, cache.θ⁻)
    params.equs[:ϑ](t₀, cache.q⁺, cache.q⁺, cache.θ⁺)

    params.equs[:ϑ](t₁, cache.q̄,  cache.q̄,  cache.Θ̅)
    params.equs[:ϑ](t₁, cache.q̄⁻, cache.q̄⁻, cache.Θ̅⁻)

    # compute projection
    params.equs[:g](t₀, cache.q,  cache.λ,  cache.g)
    params.equs[:g](t₀, cache.q⁻, cache.λ⁻, cache.g⁻)
    params.equs[:g](t₀, cache.q⁺, cache.λ⁺, cache.g⁺)
    params.equs[:g](t₁, cache.q̄⁻, cache.λ̄⁻, cache.ḡ⁻)
    params.equs[:g](t₀, cache.q,  cache.λ⁺, cache.h⁺)
    params.equs[:g](t₁, cache.q̄,  cache.λ̄⁻, cache.h̅⁻)

    # # compute ϑ
    # params.equs[:ϑ](t₀, cache.ϕ⁻, cache.ϕ⁻, cache.θ⁻)
    # params.equs[:ϑ](t₀, cache.ϕ⁺, cache.ϕ⁺, cache.θ⁺)
    # params.equs[:ϑ](t₁, cache.ϕ̅⁻, cache.ϕ̅⁻, cache.Θ̅⁻)
    #
    # # compute projection
    # params.equs[:g](t₀, cache.ϕ⁻, cache.λ⁻, cache.g⁻)
    # params.equs[:g](t₀, cache.ϕ⁺, cache.λ⁺, cache.g⁺)
    # params.equs[:g](t₁, cache.ϕ̅⁻, cache.λ̄⁻, cache.ḡ⁻)

    # compute pₙ₊₁ = ϑ(qₙ₊₁⁻) + ∇ϑ(qₙ₊₁)⋅(qₙ₊₁-qₙ₊₁⁻)
    #    i.e. p̄ = ϑ(q̄⁻) + ∇ϑ(q̄)⋅(q̄-q̄⁻)
    cache.p̄ .= cache.Θ̅⁻ .+ cache.h̅⁻
end

function compute_rhs!(b::Vector{ST}, cache::IntegratorCacheDGVI{ST,D,S,R},
                params::ParametersDGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    local z::ST

    # compute b = - [(P-AF)]
    for i in 1:S
        for k in 1:D
            z = 0
            for j in 1:R
                z += params.b[j] * params.m[j,i] * cache.F[j][k] * params.Δt
                z += params.b[j] * params.a[j,i] * cache.P[j][k]
            end

            z += params.r⁺[i] * 0.5 * ( cache.θ[k] + cache.θ⁺[k] )
            z -= params.r⁻[i] * 0.5 * ( cache.Θ̅[k] + cache.Θ̅⁻[k] )

            z += params.r⁺[i] * 0.5 * cache.g⁺[k]
            z += params.r⁻[i] * 0.5 * cache.ḡ⁻[k]

            # z += params.r⁺[i] * cache.θ⁺[k]
            # z -= params.r⁻[i] * cache.Θ̅⁻[k]
            #
            # z += params.r⁺[i] * cache.g⁺[k]
            # z += params.r⁻[i] * cache.ḡ⁻[k]

            b[D*(i-1)+k] = z
        end
    end

    # compute b = ϑ(qₙ⁺) - pₙ - ∇ϑ(qₙ)⋅(qₙ⁺-qₙ)
    for k in 1:D
        b[D*S+k] = cache.θ⁺[k] - cache.p[k] - cache.h⁺[k]
    end
end

function update_solution!(sol::AtomicSolutionPODE, cache::IntegratorCacheDGVI)
    sol.q .= cache.q̄
    sol.p .= cache.p̄
end

function integrate_step!(int::IntegratorDGVI{DT,TT}, sol::AtomicSolutionPODE{DT,TT},
                         cache::IntegratorCacheDGVI{DT}=int.caches[DT]) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, sol)

    # compute initial guess
    initial_guess!(int, sol, cache)

    # reset cache
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, cache, int.params)

    # copy solution from cache to integrator
    update_solution!(sol, cache)

    # copy solution to initial guess
    update_vector_fields!(int.iguess, sol.t, sol.q, sol.v)
end


text:
global const APPDIR = dirname(dirname(realpath(@__FILE__())))
@show APPDIR

# @ms() macro
include("$APPDIR/utils/time.jl")

# import and execute
@ms include("$APPDIR/src/main.jl")
main()

max_stars_repo_name: antcrepin/flight-route-optimizer

max_issues_repo_name: antcrepin/flight-route-optimizer

max_forks_repo_name: antcrepin/flight-route-optimizer


text:
"""
    VectorSpaceType

Abstract type for tangent spaces, cotangent spaces, their tensor products,
exterior products, etc.

Every vector space `fiber` is supposed to provide:
* a method of constructing vectors,
* basic operations: addition, subtraction, multiplication by a scalar
  and negation (unary minus),
* [`zero_vector!(fiber, X, p)`](@ref) to construct zero vectors at point `p`,
* `allocate(X)` and `allocate(X, T)` for vector `X` and type `T`,
* `copyto!(X, Y)` for vectors `X` and `Y`,
* `number_eltype(v)` for vector `v`,
* [`vector_space_dimension(::VectorBundleFibers{<:typeof(fiber)}) where fiber`](@ref).

Optionally:
* inner product via `inner` (used to provide Riemannian metric on vector
  bundles),
* [`flat`](@ref) and [`sharp`](@ref),
* `norm` (by default uses `inner`),
* [`project_vector`](@ref) (for embedded vector spaces),
* [`representation_size`](@ref) (if support for [`ProductArray`](@ref) is desired),
* broadcasting for basic operations.
"""
abstract type VectorSpaceType end

struct TangentSpaceType <: VectorSpaceType end

struct CotangentSpaceType <: VectorSpaceType end

TCoTSpaceType = Union{TangentSpaceType,CotangentSpaceType}

const TangentSpace = TangentSpaceType()
const CotangentSpace = CotangentSpaceType()

"""
    TensorProductType(spaces::VectorSpaceType...)

Vector space type corresponding to the tensor product of given vector space
types.
"""
struct TensorProductType{TS<:Tuple} <: VectorSpaceType
    spaces::TS
end

TensorProductType(spaces::VectorSpaceType...) = TensorProductType{typeof(spaces)}(spaces)

"""
    VectorBundleFibers(fiber::VectorSpaceType, M::Manifold)

Type representing a family of vector spaces (fibers) of a vector bundle over `M`
with vector spaces of type `fiber`. In contrast with `VectorBundle`, operations
on `VectorBundleFibers` expect point-like and vector-like parts to be
passed separately instead of being bundled together. It can be thought of
as a representation of vector spaces from a vector bundle but without
storing the point at which a vector space is attached (which is specified
separately in various functions).
"""
struct VectorBundleFibers{TVS<:VectorSpaceType,TM<:Manifold}
    fiber::TVS
    manifold::TM
end

const TangentBundleFibers{M} = VectorBundleFibers{TangentSpaceType,M} where {M<:Manifold}

TangentBundleFibers(M::Manifold) = VectorBundleFibers(TangentSpace, M)

const CotangentBundleFibers{M} =
    VectorBundleFibers{CotangentSpaceType,M} where {M<:Manifold}

CotangentBundleFibers(M::Manifold) = VectorBundleFibers(CotangentSpace, M)

"""
    VectorSpaceAtPoint(fiber::VectorBundleFibers, p)

A vector space (fiber type `fiber` of a vector bundle) at point `p` from
the manifold `fiber.manifold`.
"""
struct VectorSpaceAtPoint{TFiber<:VectorBundleFibers,TX}
    fiber::TFiber
    point::TX
end

"""
    TangentSpaceAtPoint(M::Manifold, p)

Return an object of type [`VectorSpaceAtPoint`](@ref) representing tangent
space at `p`.
"""
TangentSpaceAtPoint(M::Manifold, p) = VectorSpaceAtPoint(TangentBundleFibers(M), p)

"""
    CotangentSpaceAtPoint(M::Manifold, p)

Return an object of type [`VectorSpaceAtPoint`](@ref) representing cotangent
space at `p`.
"""
CotangentSpaceAtPoint(M::Manifold, p) = VectorSpaceAtPoint(CotangentBundleFibers(M), p)

"""
    VectorBundle(M::Manifold, type::VectorSpaceType)

Vector bundle on manifold `M` of type `type`.
"""
struct VectorBundle{TVS<:VectorSpaceType,TM<:Manifold} <: Manifold
    type::TVS
    manifold::TM
    fiber::VectorBundleFibers{TVS,TM}
end

function VectorBundle(fiber::TVS, M::TM) where {TVS<:VectorSpaceType,TM<:Manifold}
    return VectorBundle{TVS,TM}(fiber, M, VectorBundleFibers(fiber, M))
end

const TangentBundle{M} = VectorBundle{TangentSpaceType,M} where {M<:Manifold}

TangentBundle(M::Manifold) = VectorBundle(TangentSpace, M)

const CotangentBundle{M} = VectorBundle{CotangentSpaceType,M} where {M<:Manifold}

CotangentBundle(M::Manifold) = VectorBundle(CotangentSpace, M)

"""
    FVector(type::VectorSpaceType, data)

Decorator indicating that the vector `data` is from a fiber of a vector bundle
of type `type`.
"""
struct FVector{TType<:VectorSpaceType,TData}
    type::TType
    data::TData
end

const TFVector = FVector{TangentSpaceType}
const CoTFVector = FVector{CotangentSpaceType}

struct VectorBundleBasisData{BBasis<:CachedBasis, TBasis<:CachedBasis}
    base_basis::BBasis
    vec_basis::TBasis
end

(+)(X::FVector, Y::FVector) = FVector(X.type, X.data + Y.data)

(-)(X::FVector, Y::FVector) = FVector(X.type, X.data - Y.data)
(-)(X::FVector) = FVector(X.type, -X.data)

(*)(a::Number, X::FVector) = FVector(X.type, a * X.data)

function copyto!(X::FVector, Y::FVector)
    copyto!(X.data, Y.data)
    return X
end

base_manifold(B::VectorBundleFibers) = base_manifold(B.manifold)
base_manifold(B::VectorSpaceAtPoint) = base_manifold(B.fiber)
base_manifold(B::VectorBundle) = base_manifold(B.manifold)

"""
    bundle_projection(B::VectorBundle, x::ProductRepr)

Projection of point `p` from the bundle `M` to the base manifold.
Returns the point on the base manifold `B.manifold` at which the vector part
of `p` is attached.
"""
bundle_projection(B::VectorBundle, p) = submanifold_component(B.manifold, p, Val(1))

"""
    distance(B::VectorBundleFibers, p, X, Y)

Distance between vectors `X` and `Y` from the vector space at point `p`
from the manifold `B.manifold`, that is the base manifold of `M`.
"""
distance(B::VectorBundleFibers, p, X, Y) = norm(B, p, X - Y)
@doc raw"""
    distance(B::VectorBundle, p, q)

Distance between points $x$ and $y$ from the
vector bundle `B` over manifold `B.fiber` (denoted $\mathcal M$).

Notation:
  * The point $p = (x_p, V_p)$ where $x_p  ∈ \mathcal M$ and $V_p$ belongs to the
    fiber $F=π^{-1}(\{x_p\})$ of the vector bundle $B$ where $π$ is the
    canonical projection of that vector bundle $B$.
    Similarly, $q = (x_q, V_q)$.

The distance is calculated as

$d_B(x, y) = \sqrt{d_M(x_p, x_q)^2 + d_F(V_p, V_{q←p})^2}$

where $d_\mathcal M$ is the distance on manifold $\mathcal M$, $d_F$ is the distance
between two vectors from the fiber $F$ and $V_{q←p}$ is the result
of parallel transport of vector $V_q$ to point $x_p$. The default
behavior of [`vector_transport_to`](@ref) is used to compute the vector
transport.
"""
function distance(B::VectorBundle, p, q)
    xp, Vp = submanifold_components(B.manifold, p)
    xq, Vq = submanifold_components(B.manifold, q)
    dist_man = distance(B.manifold, xp, xq)
    vy_x = vector_transport_to(B.manifold, xq, Vq, xp)
    dist_vec = distance(B.fiber, xp, Vp, vy_x)
    return sqrt(dist_man^2 + dist_vec^2)
end

function number_eltype(::Type{FVector{TType,TData}}) where {TType<:VectorSpaceType,TData}
    return number_eltype(TData)
end
number_eltype(v::FVector) = number_eltype(v.data)

@doc raw"""
    exp(B::VectorBundle, p, X)

Exponential map of tangent vector $X$ at point $p$ from
vector bundle `B` over manifold `B.fiber` (denoted $\mathcal M$).

Notation:
  * The point $p = (x_p, V_p)$ where $x_p ∈ \mathcal M$ and $V_p$ belongs to the
    fiber $F=π^{-1}(\{x_p\})$ of the vector bundle $B$ where $π$ is the
    canonical projection of that vector bundle $B$.
  * The tangent vector $X = (V_{X,M}, V_{X,F}) ∈ T_pB$ where
    $V_{X,M}$ is a tangent vector from the tangent space $T_{x_p}\mathcal M$ and
    $V_{X,F}$ is a tangent vector from the tangent space $T_{V_p}F$ (isomorphic to $F$).

The exponential map is calculated as

$\exp_p(X) = (\exp_{x_p}(V_{X,M}), V_{\exp})$

where $V_{\exp}$ is the result of vector transport of $V_p + V_{X,F}$
to the point $\exp_{x_p}(V_{X,M})$.
The sum $V_p + V_{X,F}$ corresponds to the exponential map in the vector space $F$.
"""
exp(::VectorBundle, ::Any)

function exp!(B::VectorBundle, q, p, X)
    xp, Xp = submanifold_components(B.manifold, p)
    xq, Xq = submanifold_components(B.manifold, q)
    VXM, VXF = submanifold_components(B.manifold, X)
    exp!(B.manifold, xq, xp, VXM)
    vector_transport_to!(B.manifold, Xq, xp, Xp + VXF, xq)
    return q
end

@doc raw"""
    flat(M::Manifold, p, X::FVector)

Compute the flat isomorphism (one of the musical isomorphisms) of tangent vector `X`
from the vector space of type `M` at point `p` from the underlying [`Manifold`](@ref).

The function can be used for example to transform vectors
from the tangent bundle to vectors from the cotangent bundle
$♭ : T\mathcal M → T^{*}\mathcal M$
"""
function flat(M::Manifold, p, X::FVector)
    ξ = allocate_result(M, flat, X, p)
    return flat!(M, ξ, p, X)
end

function flat!(M::Manifold, ξ::FVector, p, X::FVector)
    error(
        "flat! not implemented for vector bundle fibers space " *
        "of type $(typeof(M)), vector of type $(typeof(ξ)), point of " *
        "type $(typeof(p)) and vector of type $(typeof(X)).",
    )
end

@decorator_transparent_signature flat!(
    M::AbstractDecoratorManifold,
    ξ::CoTFVector,
    p,
    X::TFVector,
)

function get_basis(M::VectorBundle, p, B::AbstractBasis)
    xp1 = submanifold_component(p, Val(1))
    base_basis = get_basis(M.manifold, xp1, B)
    vec_basis = get_basis(M.fiber, xp1, B)
    return CachedBasis(B, VectorBundleBasisData(base_basis, vec_basis))
end
function get_basis(M::VectorBundle, p, B::CachedBasis)
    return invoke(get_basis, Tuple{Manifold, Any, CachedBasis}, M, p, B)
end

function get_basis(M::VectorBundle, p, B::DiagonalizingOrthonormalBasis)
    xp1 = submanifold_component(p, Val(1))
    bv1 = DiagonalizingOrthonormalBasis(submanifold_component(B.frame_direction, Val(1)))
    b1 = get_basis(M.manifold, xp1, bv1)
    bv2 = DiagonalizingOrthonormalBasis(submanifold_component(B.frame_direction, Val(2)))
    b2 = get_basis(M.fiber, xp1, bv2)
    return CachedBasis(B, VectorBundleBasisData(b1, b2))
end

for BT in [
    DefaultOrthonormalBasis,
    ProjectedOrthonormalBasis{:gram_schmidt,ℝ},
    ProjectedOrthonormalBasis{:svd,ℝ},
]
    eval(quote
        @invoke_maker 3 AbstractBasis get_basis(M::VectorBundle, p, B::$BT)
    end)
end
function get_basis(M::TangentBundleFibers, p, B::AbstractBasis)
    return get_basis(M.manifold, p, B)
end

function get_coordinates!(M::VectorBundle, Y, p, X, B::AbstractBasis)
    px, Vx = submanifold_components(M.manifold, p)
    VXM, VXF = submanifold_components(M.manifold, X)
    n = manifold_dimension(M.manifold)
    get_coordinates!(M.manifold, view(Y, 1:n), px, VXM, B)
    get_coordinates!(M.fiber, view(Y, n+1:length(Y)), px, VXF, B)
    return Y
end
function get_coordinates!(
    M::VectorBundle,
    Y,
    p,
    X,
    B::CachedBasis{<:AbstractBasis,<:VectorBundleBasisData},
) where {N}
     px, Vx = submanifold_components(M.manifold, p)
     VXM, VXF = submanifold_components(M.manifold, X)
     n = manifold_dimension(M.manifold)
     get_coordinates!(M.manifold, view(Y, 1:n), px, VXM, B.data.base_basis)
     get_coordinates!(M.fiber, view(Y, n+1:length(Y)), px, VXF, B.data.vec_basis)
     return Y
end
for BT in [
    DefaultBasis,
    DefaultOrthogonalBasis,
    DefaultOrthonormalBasis,
    ProjectedOrthonormalBasis{:gram_schmidt,ℝ},
    ProjectedOrthonormalBasis{:svd,ℝ},
    VeeOrthogonalBasis,
]
    eval(quote
        @invoke_maker 5 AbstractBasis get_coordinates!(M::VectorBundle, Y, p, X, B::$BT)
    end)
end
for BT in [
    CachedBasis{<:AbstractBasis{ℝ},<:VectorBundleBasisData},
    CachedBasis{<:ManifoldsBase.AbstractOrthogonalBasis{ℝ},<:VectorBundleBasisData},
    CachedBasis{<:ManifoldsBase.AbstractOrthonormalBasis{ℝ},<:VectorBundleBasisData},
    CachedBasis{<:AbstractBasis{ℂ},<:VectorBundleBasisData},
]
    eval(quote
        @invoke_maker 5 CachedBasis{<:AbstractBasis,<:VectorBundleBasisData} get_coordinates!(
                M::VectorBundle,
                Y,
                p,
                X,
                B::$BT,
            )
    end)
end

function get_coordinates!(M::TangentBundleFibers, Y, p, X, B::ManifoldsBase.all_uncached_bases) where {N}
    return get_coordinates!(M.manifold, Y, p, X, B)
end

function get_vector!(M::VectorBundle, Y, p, X, B::DefaultOrthonormalBasis)
    n = manifold_dimension(M.manifold)
    xp1 = submanifold_component(p, Val(1))
    get_vector!(M.manifold, submanifold_component(Y, Val(1)), xp1, X[1:n], B)
    get_vector!(M.fiber, submanifold_component(Y, Val(2)), xp1, X[n+1:end], B)
    return Y
end
function get_vector!(
    M::VectorBundle,
    Y,
    p,
    X,
    B::CachedBasis{<:AbstractBasis,<:VectorBundleBasisData},
) where {N}
    n = manifold_dimension(M.manifold)
    xp1 = submanifold_component(p, Val(1))
    get_vector!(
        M.manifold,
        submanifold_component(Y, Val(1)),
        xp1,
        X[1:n],
        B.data.base_basis,
    )
    get_vector!(
        M.fiber,
        submanifold_component(Y, Val(2)),
        xp1,
        X[n+1:end],
        B.data.vec_basis,
    )
    return Y
end
function get_vector!(M::TangentBundleFibers, Y, p, X, B::ManifoldsBase.all_uncached_bases) where {N}
    return get_vector!(M.manifold, Y, p, X, B)
end

function get_vectors(
    M::VectorBundle,
    p,
    B::CachedBasis{<:AbstractBasis,<:VectorBundleBasisData},
)
    xp1 = submanifold_component(p, Val(1))
    zero_m = zero_tangent_vector(M.manifold, xp1)
    zero_f = zero_vector(M.fiber, xp1)
    vs = typeof(ProductRepr(zero_m, zero_f))[]
    for bv in get_vectors(M.manifold, xp1, B.data.base_basis)
        push!(vs, ProductRepr(bv, zero_f))
    end
    for bv in get_vectors(M.fiber, xp1, B.data.vec_basis)
        push!(vs, ProductRepr(zero_m, bv))
    end
    return vs
end

function get_vectors(M::VectorBundleFibers, p, B::CachedBasis)
    return get_vectors(M.manifold, p, B)
end
Base.@propagate_inbounds getindex(x::FVector, i) = getindex(x.data, i)

"""
    inner(B::VectorBundleFibers, p, X, Y)

Inner product of vectors `X` and `Y` from the vector space of type `B.fiber`
at point `p` from manifold `B.manifold`.
"""
function inner(B::VectorBundleFibers, p, X, Y)
    error(
        "inner not defined for vector space family of type $(typeof(B)), " *
        "point of type $(typeof(p)) and " *
        "vectors of types $(typeof(X)) and $(typeof(Y)).",
    )
end
inner(B::VectorBundleFibers{<:TangentSpaceType}, p, X, Y) = inner(B.manifold, p, X, Y)
function inner(B::VectorBundleFibers{<:CotangentSpaceType}, p, X, Y)
    return inner(
        B.manifold,
        p,
        sharp(B.manifold, p, FVector(CotangentSpace, X)).data,
        sharp(B.manifold, p, FVector(CotangentSpace, Y)).data,
    )
end
@doc raw"""
    inner(B::VectorBundle, p, X, Y)

Inner product of tangent vectors `X` and `Y` at point `p` from the
vector bundle `B` over manifold `B.fiber` (denoted $\mathcal M$).

Notation:
  * The point $p = (x_p, V_p)$ where $x_p ∈ \mathcal M$ and $V_p$ belongs to the
    fiber $F=π^{-1}(\{x_p\})$ of the vector bundle $B$ where $π$ is the
    canonical projection of that vector bundle $B$.
  * The tangent vector $v = (V_{X,M}, V_{X,F}) ∈ T_{x}B$ where
    $V_{X,M}$ is a tangent vector from the tangent space $T_{x_p}\mathcal M$ and
    $V_{X,F}$ is a tangent vector from the tangent space $T_{V_p}F$ (isomorphic to $F$).
    Similarly for the other tangent vector $w = (V_{Y,M}, V_{Y,F}) ∈ T_{x}B$.

The inner product is calculated as

$⟨X, Y⟩_p = ⟨V_{X,M}, V_{Y,M}⟩_{x_p} + ⟨V_{X,F}, V_{Y,F}⟩_{V_p}.$
"""
function inner(B::VectorBundle, p, X, Y)
    px, Vx = submanifold_components(B.manifold, p)
    VXM, VXF = submanifold_components(B.manifold, X)
    VYM, VYF = submanifold_components(B.manifold, Y)
    return inner(B.manifold, px, VXM, VYM) + inner(B.fiber, Vx, VXF, VYF)
end

function isapprox(B::VectorBundle, p, q; kwargs...)
    xp, Vp = submanifold_components(B.manifold, p)
    xq, Vq = submanifold_components(B.manifold, q)
    return isapprox(B.manifold, xp, xq; kwargs...) && isapprox(Vp, Vq; kwargs...)
end
function isapprox(B::VectorBundle, p, X, Y; kwargs...)
    px, Vx = submanifold_components(B.manifold, p)
    VXM, VXF = submanifold_components(B.manifold, X)
    VYM, VYF = submanifold_components(B.manifold, Y)
    return isapprox(B.manifold, VXM, VYM; kwargs...) &&
           isapprox(B.manifold, px, VXF, VYF; kwargs...)
end

@doc raw"""
    log(B::VectorBundle, p, q)

Logarithmic map of the point `y` at point `p` from
vector bundle `B` over manifold `B.fiber` (denoted $\mathcal M$).

Notation:
  * The point $p = (x_p, V_p)$ where $x_p ∈ \mathcal M$ and $V_p$ belongs to the
    fiber $F=π^{-1}(\{x_p\})$ of the vector bundle $B$ where $π$ is the
    canonical projection of that vector bundle $B$.
    Similarly, $q = (x_q, V_q)$.

The logarithmic map is calculated as

$\log_p q = (\log_{x_p}(x_q), V_{\log} - V_p)$

where $V_{\log}$ is the result of vector transport of $V_q$ to the point $x_p$.
The difference $V_{\log} - V_p$ corresponds to the logarithmic map in the vector space $F$.
"""
log(::VectorBundle, ::Any...)

function log!(B::VectorBundle, X, p, q)
    px, Vx = submanifold_components(B.manifold, p)
    py, Vy = submanifold_components(B.manifold, q)
    VXM, VXF = submanifold_components(B.manifold, X)
    log!(B.manifold, VXM, px, py)
    vector_transport_to!(B.manifold, VXF, py, Vy, px)
    copyto!(VXF, VXF - Vx)
    return X
end

function manifold_dimension(B::VectorBundle)
    return manifold_dimension(B.manifold) + vector_space_dimension(B.fiber)
end

"""
    norm(B::VectorBundleFibers, p, q)

Norm of the vector `X` from the vector space of type `B.fiber`
at point `p` from manifold `B.manifold`.
"""
norm(B::VectorBundleFibers, p, X) = sqrt(inner(B, p, X, X))
norm(B::VectorBundleFibers{<:TangentSpaceType}, p, X) = norm(B.manifold, p, X)

@doc raw"""
    project(B::VectorBundle, p)

Project the point `p` from the ambient space of the vector bundle `B`
over manifold `B.fiber` (denoted $\mathcal M$) to the vector bundle.

Notation:
  * The point $p = (x_p, V_p)$ where $x_p$ belongs to the ambient space of $\mathcal M$
    and $V_p$ belongs to the ambient space of the
    fiber $F=π^{-1}(\{x_p\})$ of the vector bundle $B$ where $π$ is the
    canonical projection of that vector bundle $B$.

The projection is calculated by projecting the point $x_p$ to the manifold $\mathcal M$
and then projecting the vector $V_p$ to the tangent space $T_{x_p}\mathcal M$.
"""
project(::VectorBundle, ::Any)

function project!(B::VectorBundle, q, p)
    px, pVx = submanifold_components(B.manifold, p)
    qx, qVx = submanifold_components(B.manifold, q)
    project!(B.manifold, qx, px)
    project!(B.manifold, qVx, qx, pVx)
    return q
end

@doc raw"""
    project(B::VectorBundle, p, X)

Project the element `X` of the ambient space of the tangent space $T_p B$
to the tangent space $T_p B$.

Notation:
  * The point $p = (x_p, V_p)$ where $x_p ∈ \mathcal M$ and $V_p$ belongs to the
    fiber $F=π^{-1}(\{x_p\})$ of the vector bundle $B$ where $π$ is the
    canonical projection of that vector bundle $B$.
  * The vector $x = (V_{X,M}, V_{X,F})$ where $x_p$ belongs to the ambient space of $T_{x_p}\mathcal M$
    and $V_{X,F}$ belongs to the ambient space of the
    fiber $F=π^{-1}(\{x_p\})$ of the vector bundle $B$ where $π$ is the
    canonical projection of that vector bundle $B$.

The projection is calculated by projecting $V_{X,M}$ to tangent space $T_{x_p}\mathcal M$
and then projecting the vector $V_{X,F}$ to the fiber $F$.
"""
project(::VectorBundle, ::Any, ::Any)

function project!(B::VectorBundle, Y, p, X)
    px, Vx = submanifold_components(B.manifold, p)
    VXM, VXF = submanifold_components(B.manifold, X)
    VYM, VYF = submanifold_components(B.manifold, Y)
    project!(B.manifold, VYM, px, VXM)
    project!(B.manifold, VYF, px, VXF)
    return Y
end

"""
    project_vector(B::VectorBundleFibers, p, X)

Project vector `X` from the vector space of type `B.fiber` at point `p`.
"""
function project_vector(B::VectorBundleFibers, p, X)
    Y = allocate_result(B, project_vector, p, X)
    return project_vector!(B, Y, p, X)
end

function project_vector!(B::VectorBundleFibers{<:TangentSpaceType}, Y, p, X)
    return project!(B.manifold, Y, p, X)
end
function project_vector!(B::VectorBundleFibers, Y, p, X)
    error("project_vector! not implemented for vector space family of type $(typeof(B)), output vector of type $(typeof(Y)) and input vector at point $(typeof(p)) with type of w $(typeof(X)).")
end

Base.@propagate_inbounds setindex!(x::FVector, val, i) = setindex!(x.data, val, i)

function representation_size(B::VectorBundleFibers{<:TCoTSpaceType})
    return representation_size(B.manifold)
end
function representation_size(B::VectorBundle)
    len_manifold = prod(representation_size(B.manifold))
    len_vs = prod(representation_size(B.fiber))
    return (len_manifold + len_vs,)
end

@doc raw"""
    sharp(M::Manifold, p, ξ::FVector)

Compute the sharp isomorphism (one of the musical isomorphisms) of vector `ξ`
from the vector space `M` at point `p` from the underlying [`Manifold`](@ref).

The function can be used for example to transform vectors
from the cotangent bundle to vectors from the tangent bundle
$♯ : T^{*}\mathcal M → T\mathcal M$
"""
function sharp(M::Manifold, p, ξ::FVector)
    X = allocate_result(M, sharp, ξ, p)
    return sharp!(M, X, p, ξ)
end

function sharp!(M::Manifold, X::FVector, p, ξ::FVector)
    error(
        "sharp! not implemented for vector bundle fibers space " *
        "of type $(typeof(M)), vector of type $(typeof(X)), point of " *
        "type $(typeof(p)) and vector of type $(typeof(ξ)).",
    )
end

@decorator_transparent_signature sharp!(
    M::AbstractDecoratorManifold,
    X::TFVector,
    p,
    ξ::CoTFVector,
)

show(io::IO, ::TangentSpaceType) = print(io, "TangentSpace")
show(io::IO, ::CotangentSpaceType) = print(io, "CotangentSpace")
function show(io::IO, tpt::TensorProductType)
    print(io, "TensorProductType(", join(tpt.spaces, ", "), ")")
end
function show(io::IO, fiber::VectorBundleFibers)
    print(io, "VectorBundleFibers($(fiber.fiber), $(fiber.manifold))")
end
function show(io::IO, mime::MIME"text/plain", vs::VectorSpaceAtPoint)
    summary(io, vs)
    println(io, "\nFiber:")
    pre = " "
    sf = sprint(show, "text/plain", vs.fiber; context = io, sizehint = 0)
    sf = replace(sf, '\n' => "\n$(pre)")
    println(io, pre, sf)
    println(io, "Base point:")
    sp = sprint(show, "text/plain", vs.point; context = io, sizehint = 0)
    sp = replace(sp, '\n' => "\n$(pre)")
    print(io, pre, sp)
end
show(io::IO, vb::VectorBundle) = print(io, "VectorBundle($(vb.type), $(vb.manifold))")
show(io::IO, vb::TangentBundle) = print(io, "TangentBundle($(vb.manifold))")
show(io::IO, vb::CotangentBundle) = print(io, "CotangentBundle($(vb.manifold))")

allocate(x::FVector) = FVector(x.type, allocate(x.data))
allocate(x::FVector, ::Type{T}) where {T} = FVector(x.type, allocate(x.data, T))

"""
    allocate_result(B::VectorBundleFibers, f, x...)

Allocates an array for the result of function `f` that is
an element of the vector space of type `B.fiber` on manifold `B.manifold`
and arguments `x...` for implementing the non-modifying operation
using the modifying operation.
"""
function allocate_result(B::VectorBundleFibers, f, x...)
    T = allocate_result_type(B, f, x)
    return allocate(x[1], T)
end
function allocate_result(M::Manifold, ::typeof(flat), w::TFVector, x)
    return FVector(CotangentSpace, allocate(w.data))
end
function allocate_result(M::Manifold, ::typeof(sharp), w::CoTFVector, x)
    return FVector(TangentSpace, allocate(w.data))
end

"""
    allocate_result_type(B::VectorBundleFibers, f, args::NTuple{N,Any}) where N

Returns type of element of the array that will represent the result of
function `f` for representing an operation with result in the vector space `fiber`
for manifold `M` on given arguments (passed at a tuple).
"""
function allocate_result_type(B::VectorBundleFibers, f, args::NTuple{N,Any}) where {N}
    T = typeof(reduce(+, one(number_eltype(eti)) for eti ∈ args))
    return T
end

size(x::FVector) = size(x.data)

function submanifold_component(M::Manifold, x::FVector, i::Val)
    return submanifold_component(M, x.data, i)
end
submanifold_component(x::FVector, i::Val) = submanifold_component(x.data, i)

submanifold_components(M::Manifold, x::FVector) = submanifold_components(M, x.data)
submanifold_components(x::FVector) = submanifold_components(x.data)

"""
    vector_space_dimension(B::VectorBundleFibers)

Dimension of the vector space of type `B`.
"""
function vector_space_dimension(B::VectorBundleFibers)
    error("vector_space_dimension not implemented for vector space family $(typeof(B)).")
end
function vector_space_dimension(B::VectorBundleFibers{<:TCoTSpaceType})
    return manifold_dimension(B.manifold)
end
function vector_space_dimension(B::VectorBundleFibers{<:TensorProductType})
    dim = 1
    for space in B.fiber.spaces
        dim *= vector_space_dimension(VectorBundleFibers(space, B.manifold))
    end
    return dim
end

"""
    zero_vector(B::VectorBundleFibers, p)

Compute the zero vector from the vector space of type `B.fiber` at point `p`
from manifold `B.manifold`.
"""
function zero_vector(B::VectorBundleFibers, p)
    X = allocate_result(B, zero_vector, p)
    return zero_vector!(B, X, p)
end

"""
    zero_vector!(B::VectorBundleFibers, X, p)

Save the zero vector from the vector space of type `B.fiber` at point `p`
from manifold `B.manifold` to `X`.
"""
function zero_vector!(B::VectorBundleFibers, X, p)
    error("zero_vector! not implemented for vector space family of type $(typeof(B)).")
end
function zero_vector!(B::VectorBundleFibers{<:TangentSpaceType}, X, p)
    return zero_tangent_vector!(B.manifold, X, p)
end

@doc raw"""
    zero_tangent_vector(B::VectorBundle, p)

Zero tangent vector at point `p` from the vector bundle `B`
over manifold `B.fiber` (denoted $\mathcal M$). The zero vector belongs to the space $T_{p}B$

Notation:
  * The point $p = (x_p, V_p)$ where $x_p ∈ \mathcal M$ and $V_p$ belongs to the
    fiber $F=π^{-1}(\{x_p\})$ of the vector bundle $B$ where $π$ is the
    canonical projection of that vector bundle $B$.

The zero vector is calculated as

$\mathbf{0}_{p} = (\mathbf{0}_{x_p}, \mathbf{0}_F)$

where $\mathbf{0}_{x_p}$ is the zero tangent vector from $T_{x_p}\mathcal M$ and
$\mathbf{0}_F$ is the zero element of the vector space $F$.
"""
zero_tangent_vector(::VectorBundle, ::Any...)

function zero_tangent_vector!(B::VectorBundle, X, p)
    xp, Vp = submanifold_components(B.manifold, p)
    VXM, VXF = submanifold_components(B.manifold, X)
    zero_tangent_vector!(B.manifold, VXM, xp)
    zero_vector!(B.fiber, VXF, Vp)
    return X
end


text:
######################################################################################
## Tcmd matrix
######################################################################################
using CpelNano

# Parameter vectors
a1 = -10.0; b1 = -100.0; c1 = 20.0; ϕ1 = 0.5 * [a1,b1,c1];
a2 = 10.0; b2 = 100.0; c2 = 20.0; ϕ2 = 0.5 * [a2,b2,c2];

# Get group data
s1 = 4; s2 = 4; nse_sc = 1.0; seed = 123;
ms_g1, ms_g2 = CpelNano.gen_grp_comp_data(s1, s2, ϕ1, ϕ2, nse_sc, seed);

# Compute matrix of CMDs
cmd_tbl = CpelNano.get_unmat_cmd_tbl(ms_g1, ms_g2)

######################################################################################
## All unmatched tests
######################################################################################
using CpelNano

# Parameter vectors
a1 = -10.0; b1 = -100.0; c1 = 20.0; ϕ1 = 0.5 * [a1,b1,c1];
a2 = 10.0; b2 = 100.0; c2 = 20.0; ϕ2 = 0.5 * [a2,b2,c2];

# Get group data
s1 = 4; s2 = 4; nse_sc = 0.01; seed = 123;
ms_g1, ms_g2 = CpelNano.gen_grp_comp_data(s1, s2, ϕ1, ϕ2, nse_sc, seed);

# Compute statistic
test_cont_1 = CpelNano.unmat_est_reg_test(ms_g1, ms_g1)
test_cont_2 = CpelNano.unmat_est_reg_test(ms_g2, ms_g2)
test_treat = CpelNano.unmat_est_reg_test(ms_g1, ms_g2)

######################################################################################
## All matched tests
######################################################################################
using CpelNano

# Parameter vectors
a1 = -10.0; b1 = -100.0; c1 = 20.0; ϕ1 = 0.05 * [a1,b1,c1];
a2 = 10.0; b2 = 100.0; c2 = 20.0; ϕ2 = 0.05 * [a2,b2,c2];

# Get group data
s1 = 4; s2 = 4; nse_sc = 1.0; seed = 123;
ms_g1, ms_g2 = CpelNano.gen_grp_comp_data(s1, s2, ϕ1, ϕ2, nse_sc, seed);

# Perform testing
test_cont_1 = CpelNano.mat_est_reg_test(ms_g1, ms_g1)
test_cont_2 = CpelNano.mat_est_reg_test(ms_g2, ms_g2)
test_treat = CpelNano.mat_est_reg_test(ms_g1, ms_g2)

######################################################################################
## All two-sample tests
######################################################################################
using CpelNano

L = 20; N = L; M = 25; pobs = 1.0; μ_m = 2.0; σ_m = 1.0; μ_u = -2.0; σ_u = 1.0;
Nl = fill(1.0, L); ρl = vcat(fill(1 / 1000, 8), fill(1 / 10, 4), fill(1 / 1000, 8)...);
dl = fill(150.0, L - 1); chrst = 1; chrend = 3000; cpg_pos = collect(100:150:chrend);

# Config
config = CpelNano.CpelNanoConfig()
config.max_em_iters = 50
config.max_em_init = 8
config.verbose = true

# Parameter vector
a = 0.0; b = 100.0; c = 0.0; ϕ = [a,b,c];

# Other parameters
x = CpelNano.RegStruct(); x.N = N; x.L = L; x.Nl = Nl; x.ρl = ρl; x.dl = dl;
x.chrst = chrst; x.chrend = chrend; x.cpg_pos = cpg_pos; x.ϕhat = ϕ;
α, β = CpelNano.get_αβ_from_ϕ(ϕ, x);
CpelNano.get_stat_sums!(x, config);

# Simulate data
rs = CpelNano.cpel_samp_ont(M, α, β, pobs, μ_m, σ_m, μ_u, σ_u);
rs.L = x.L; rs.Nl = x.Nl; rs.ρl = x.ρl; rs.dl = x.dl; rs.N = x.L; rs.ρn = x.ρl; rs.dn = x.dl;
rs.chrst = chrst; rs.chrend = chrend; rs.cpg_pos = cpg_pos;

# Run EM algorithm
CpelNano.get_ϕhat!(rs,config); rs.ϕhat;
CpelNano.get_stat_sums!(rs, config);
rs_1 = rs;

# Parameter vector
a = 10.0; b = 150.0; c = 50.0; ϕ = [a,b,c];

# Other parameters
x = CpelNano.RegStruct(); x.N = N; x.L = L; x.Nl = Nl; x.ρl = ρl; x.dl = dl;
x.chrst = chrst; x.chrend = chrend; x.cpg_pos = cpg_pos; x.ϕhat = ϕ;
α, β = CpelNano.get_αβ_from_ϕ(ϕ, x);
CpelNano.get_stat_sums!(x, config);

# Simulate data
rs = CpelNano.cpel_samp_ont(M, α, β, pobs, μ_m, σ_m, μ_u, σ_u);
rs.L = x.L; rs.Nl = x.Nl; rs.ρl = x.ρl; rs.dl = x.dl; rs.N = x.L; rs.ρn = x.ρl; rs.dn = x.dl;
rs.chrst = chrst; rs.chrend = chrend; rs.cpg_pos = cpg_pos;

# Run EM algorithm
CpelNano.get_ϕhat!(rs,config); rs.ϕhat;
CpelNano.get_stat_sums!(rs, config);
rs_2 = rs;

# Config
config.LMAX_TWO_SAMP = 50
config.max_em_iters = 50
config.max_em_init = 4
config.verbose = false

# Perform testing
test_treat = CpelNano.pmap_diff_two_samp_comp_test(rs_1, rs_2, config)

######################################################################################
## User interface
######################################################################################
using CpelNano

## Unmatched

# Fasta
fasta = "examples/full_example/reference/hg38_chr17_43023997_43145780.fa"

# Model files
mod_fls_g1 = ["examples/grp_comp_example/mod_files/g1_rep$(i).txt" for i = 1:6]
mod_fls_g2 = ["examples/grp_comp_example/mod_files/g2_rep$(i).txt" for i = 1:6]

# Config
config = CpelNano.CpelNanoConfig()
config.out_dir = "examples/grp_comp_example/tests/"
config.out_prefix = "cpelnano_grp_comp_unmat"
config.matched = false
config.verbose = true

# Differential analysis
CpelNano.diff_grp_comp(mod_fls_g1,mod_fls_g2,fasta,config)

## Matched

# Fasta
fasta = "examples/full_example/reference/hg38_chr17_43023997_43145780.fa"

# Model files
mod_fls_g1 = ["examples/grp_comp_example/mod_files/g1_rep$(i).txt" for i = 1:6]
mod_fls_g2 = ["examples/grp_comp_example/mod_files/g2_rep$(i).txt" for i = 1:6]

# Config
config = CpelNano.CpelNanoConfig()
config.out_dir = "examples/grp_comp_example/tests/"
config.out_prefix = "cpelnano_grp_comp_mat"
config.matched = true
config.verbose = true

# Differential analysis
CpelNano.diff_grp_comp(mod_fls_g1,mod_fls_g2,fasta,config)

#################################################################################################
## Two sample test
#################################################################################################
using Distributed
addprocs(8)
@everywhere using Pkg
@everywhere Pkg.activate("./")
@everywhere using CpelNano
using CpelNano

# Fasta
fasta = "examples/full_example/reference/hg38_chr17_43023997_43145780.fa"

# Model files
nano_s1 = "examples/full_example/nanopolish/full_example_noise2.0_methylation_calls.sorted.tsv"
nano_s2 = "examples/full_example/nanopolish/full_example_noise2.0_methylation_calls.sorted.tsv"
mod_path_s1 = "examples/two_samp_example/cpelnano/full_example_mod_nse_true_theta.txt"
mod_path_s2 = "examples/two_samp_example/cpelnano/full_example_mod_nse_true_theta.txt"

# Config
config = CpelNano.CpelNanoConfig()
# config.out_dir = "examples/two_samp_example/tests/"
config.out_dir = "/Users/jordiabante/Desktop/two_samp_example/"
config.out_prefix = "cpelnano_two_samp"
config.pval_comp = false
config.max_em_iters = 25
config.max_em_init = 4
config.verbose = true

# Differential analysis
CpelNano.diff_two_samp_comp(mod_path_s1,mod_path_s2,nano_s1,nano_s2,fasta,config)


text:
# install NodeJS modules
using NodeJS

run(Cmd(`$(npm_cmd()) install --scripts-prepend-node-path=true --production --no-package-lock --no-optional mathjax-node-cli`, dir=@__DIR__))

max_stars_repo_name: BenLauwens/NativeSVG

max_issues_repo_name: BenLauwens/NativeSVG

max_forks_repo_name: BenLauwens/NativeSVG


text:
using MNIST

### TYPE DEFINITION ###

@with_kw type MasterState
	master_mailbox = RemoteChannel(() -> PABASTO.Mailbox(), myid())
	shared_pserver_mailbox = RemoteChannel(() -> PABASTO.IntervalMailbox(), myid())
	workers::Array = Tuple{Int, Any, Any}[]
	paramservers::Array = Tuple{Int, Any, Any}[]
	num_train_examples
	num_processed_examples=0
	num_epoch=1
	max_num_epochs
	time_var=now()
	compute_loss_timeout=5
	compute_loss_last=now()
	# parameters for adaptive control policy
	tau

	num_workers
	num_paramservers
	# number of examples the master sends to worker in response to ExamplesRequestMessage
	examples_batch_size
	num_live_workers=0
	num_live_pservers=0
	final_params=nothing
	starting_params::Parameter
	num_processed_params=0
end

### INITIALIZATION ###

function initialize_nodes(state::MasterState)
	add_paramservers(state, state.num_paramservers)
	add_workers(state, state.num_workers)
end

function add_paramservers(state::MasterState, count)
	ids = add_procs(count)
	for id in ids
		index=length(state.paramservers)+1
		pserver_mailbox = RemoteChannel(() -> PABASTO.Mailbox(), id)
		#remotecall_fetch(paramserver_setup, id, state.master_mailbox, pserver_mailbox,index)
		ref = @spawnat id PABASTO.paramserver(state.master_mailbox, state.shared_pserver_mailbox,pserver_mailbox,index,state.starting_params)
		@async wait(ref)
		push!(state.paramservers, (id,ref,pserver_mailbox))
	end
	state.num_live_pservers += count
end

function remove_paramserver(state::MasterState)
	id, ref, pserver_mailbox = pop!(state.paramservers)
	put!(pserver_mailbox, CeaseOperationMessage())
end

function add_workers(state::MasterState, count)
	ids = add_procs(count)
	for id in ids
		worker_mailbox = RemoteChannel(() -> PABASTO.Mailbox(), id)
		ref = @spawnat id PABASTO.worker(id, state.master_mailbox, worker_mailbox,state.starting_params)
		@async wait(ref)
		push!(state.workers, (id, ref, worker_mailbox))
	end
	state.num_live_workers += count
end

function remove_worker(state::MasterState)
	id, ref, worker_mailbox = pop!(state.workers)
	put!(worker_mailbox, FinishOperationMessage())
end

### REQUEST HANDLING ###

# Handle request from workers for more examples
function handle(state::MasterState,request::ExamplesRequestMessage)

	if state.num_processed_examples >= state.num_train_examples && state.num_epoch < state.max_num_epochs
		state.num_processed_examples = 0;
		state.num_epoch = state.num_epoch + 1;
	elseif state.num_processed_examples >= state.num_train_examples
		# notify all workers to stop
		for (id, ref, worker_mailbox) in state.workers
			put!(worker_mailbox, FinishOperationMessage())
		end
		return
	end

	id = request.id
	worker_mailbox = request.worker_mailbox

	count = 0;
	examples = []
	while count < state.examples_batch_size && state.num_processed_examples < state.num_train_examples
		example_id = state.num_processed_examples + 1
		push!(examples, example_id)
		count = count + 1
		state.num_processed_examples = state.num_processed_examples + 1
	end

	put!(worker_mailbox, ExampleIndicesMessage(examples))
	println("[MASTER] Assigned examples $(examples[1])-$(examples[end]) to worker $(id)")
end

function handle(state::MasterState,msg::FinishedOperationMessage)
	state.num_live_workers -= 1
end

function handle(state::MasterState,msg::CeasedOperationMessage)
	state.num_live_pservers -= 1
end

function handle(state::MasterState,msg::GradientUpdateMessage)
	println("[MASTER] Dispatching GradientUpdateMessage")
	#any paramserver between 1 and 10000 can handle the gradinet update message
	put!(state.shared_pserver_mailbox,msg,1,10000)
end

function handle(state::MasterState,msg::ParameterUpdateRequestMessage)
	println("[MASTER] Dispatching ParameterUpdateRequestMessage")
	put!(state.shared_pserver_mailbox,msg,1,10000)
end

function handle(state::MasterState,msg::ParameterUpdateMessage)
	state.final_params=msg.parameters
end

function handle(state::MasterState,msg::TestLossMessage)
	println("[MASTER] Test loss is $(msg.loss)")
end

function handle(state::MasterState,msg::Void)
	println("[MASTER] Spinning")
	sleep(1)
end

function master(;starting_params=SimpleParameter(Any[PABASTO.dummy_weights1,PABASTO.dummy_biases1]))
	train_examples, train_labels = traindata()

	state = MasterState(
	num_train_examples=10000,#size(train_examples,2),
	max_num_epochs=1000,
	tau=20.0,
	num_workers=8,
	num_paramservers=3,
	examples_batch_size=100,
	starting_params=starting_params
	)
	initialize_nodes(state)

	while state.num_live_workers > 0
		handle(state,take!(state.master_mailbox))
		if Int(now()-state.compute_loss_last)>state.compute_loss_timeout*1000
			(id,ref,mailbox)=rand(state.workers)
			put!(mailbox,TestExampleIndicesMessage(1:100))
			state.compute_loss_last=now()
		end
	end

	state.num_paramservers = state.num_live_pservers

	# query the lead parameter server for
	put!(state.paramservers[1][3], ParameterUpdateRequestMessage(state.master_mailbox))

	while state.final_params == nothing
		handle(state,take!(state.master_mailbox))
	end

	# We got parameters so shut down all the paramservers
	for (id, ref, pserver_mailbox) in state.paramservers
		put!(pserver_mailbox, CeaseOperationMessage())
	end

	# write params to disk
	f = open("params.jls", "w")
	serialize(f, state.final_params)
	close(f)

	for (id,ref,mailbox) in cat(1,state.paramservers,state.workers)
		wait(ref)
	end
end


text:
@testset "sparse terms" begin
    t = PauliTerm("IXIYIZ", 2.0 + im)
    @test t == sparse_op(t)
    @test dense_op(sparse_op(t)) == t
    @test typeof(dense_op(sparse_op(t))) == typeof(t)
    @test typeof(sparse_op(t)) != typeof(t)

    t1 = PauliTerm("IXIYIZ", 2.0 + im)
    t2 = PauliTerm("XIYYZY", 1.0 + 0im)
    f1 = FermiTerm("+I+I-I-", 2.0 + im)
    f2 = FermiTerm("--I-+++", 2.0 + im)
    (st1, st2, sf1, sf2) = sparse_op.((t1, t2, f1, f2))
    @test t1 * t2 == st1 * st2
    ## Make sure that sparse and dense multiplication use different methods
    @test which(*, (typeof(t1), typeof(t2))) !=  which(*, (typeof(st1), typeof(st2)))
    @test f1 * f2 == sf1 * sf2
    @test which(*, (typeof(f1), typeof(f2))) !=  which(*, (typeof(sf1), typeof(sf2)))

    @test isone(t2 * t2)

    ps = rand_op_sum(Pauli, 100, 100)
    sps = sparse_op(ps)
    @test sort!(sps) == sparse_op(ps)
end


text:
# Give access to BSR lib code
pkgdir = joinpath(Pkg.dir(), "BSR")
bsrlibdir = joinpath(pkgdir, "deps/downloads/BSR/grouping/lib")
@mput bsrlibdir
@matlab addpath(bsrlibdir)

##############################################################################
##
## gpb
##
##############################################################################

function gpb{T<:Real}(img::Array{T})
    # img can either be grayscale (two dimensions) or color (three dimensions)
    @mput img
    @mquiet gpb_orient = globalPb_im(img)
    @mget gpb_orient
end

##############################################################################
##
## gpb-owt-ucm
##
##############################################################################

function gpb_owt_ucm_labels_ucm2{T<:Real}(img::Array{T}; scale::T=0.0)
    gpb_res = gpb(img)
    @mput gpb_res
    eval_string("ucm2 = contours2ucm(gpb_res, 'doubleSize');")
    eval_string("labels2 = bwlabel(ucm2 <= $scale);")
    eval_string("labels = labels2(2:2:end, 2:2:end);")

    labels = convert(Matrix{Int}, @mget labels)
    ucm2 = @mget ucm2

    labels, ucm2
end

function gpb_owt_ucm{T<:Real}(img::Array{T}; scale::T=0.0)
    gpb_owt_ucm_labels_ucm2(img, scale=scale)[1]
end

function gpb_owt_ucm_strength{T<:Real}(img::Array{T}; scale::T=0.0)
    # return label along with corresponding strengths between
    # each region
    labels, ucm2 = gpb_owt_ucm_labels_ucm2(img, scale=scale)

    hstrength = ucm2[2:2:end, 3:2:end]
    vstrength = ucm2[3:2:end, 2:2:end]

    nsp = maximum(labels)
    strengthdict = Dict{NTuple{2,Int},Float64}()
    nrows, ncols = size(labels)
    for j = 1:ncols
        for i = 1:nrows
            # left to right
            if j <= ncols - 1
                left = labels[i, j]
                right = labels[i, j+1]
                if left != right
                    strengthdict[(left, right)] =
                    strengthdict[(right, left)] =
                        hstrength[i, j]
                end
            end

            # top to bottom
            if i <= nrows - 1
                top = labels[i, j]
                bottom = labels[i+1, j]
                if top != bottom
                    strengthdict[(top, bottom)] =
                    strengthdict[(bottom, top)] =
                        vstrength[i, j]
                end
            end
        end
    end

    labels, strengthdict
end