text:
using luteos
using SymPy

# let # limit scope

Ps = [P1(), P2(), P3()]         # Range of polynomial order
Ns = [7, 13]#[9, 17, 33] # Range of grid size

dim = 3

# Material
mat = Material( ν=0.33, E=1.0, dim=3 )

## Set up problem
#   Get functions for exact solution
function ExactSol3D( Cstiff )

  dim = 3

  @syms x1 x2 x3

  u1 = 10 * (x2 - x2.^2) .* (x3 - x3.^2) .* sin.(pi*x1) .* (1 - x1) .* (1 - x2/2) .* (1 - x3/4)
  u2 =  2 * (x1 - x1.^2) .* (x3 - x3.^2) .* sin.(pi*x2) .* (1 - x2) .* (1 - x1/2) .* (1 - x3/4)
  u3 =  5 * (x1 - x1.^2) .* (x2 - x2.^2) .* cos.(pi*x3) .* (1 - x3/2) .* (1 - x1/2) .* (1 - x2/2)

  ∂u1∂x1 = diff( u1, x1 )
  ∂u1∂x2 = diff( u1, x2 )
  ∂u1∂x3 = diff( u1, x3 )
  ∂u2∂x1 = diff( u2, x1 )
  ∂u2∂x2 = diff( u2, x2 )
  ∂u2∂x3 = diff( u2, x3 )
  ∂u3∂x1 = diff( u3, x1 )
  ∂u3∂x2 = diff( u3, x2 )
  ∂u3∂x3 = diff( u3, x3 )

  ∇u = [ ∂u1∂x1 ∂u2∂x1 ∂u3∂x1; ∂u1∂x2 ∂u2∂x2 ∂u3∂x2; ∂u1∂x3 ∂u2∂x3 ∂u3∂x3 ]

  ϵ = 1/2 * ( ∇u + ∇u.' )

  σ = fill( 0*x1, dim, dim )

  for ii in 1:dim, jj in 1:dim, kk in 1:dim, ll in 1:dim
    σ[ii,jj] += (Cstiff[ii,jj,kk,ll] * ϵ[kk,ll])[1,1,1,1]
  end

  F  = fill( 0*x1, dim )

  ∇σ = fill( 0*x1, dim, dim, dim )

  for ii in 1:dim, jj in 1:dim
    ∇σ[ii,jj,1] = diff( σ[ii,jj], x1 )
    ∇σ[ii,jj,2] = diff( σ[ii,jj], x2 )
    ∇σ[ii,jj,3] = diff( σ[ii,jj], x3 )
  end

  for ii in 1:dim, jj in 1:dim
    F[ii] -= ∇σ[ii,jj,jj]
  end

  u1func_org = lambdify( u1, [x1, x2, x3] )
  u2func_org = lambdify( u2, [x1, x2, x3] )
  u3func_org = lambdify( u3, [x1, x2, x3] )

  F1func_org  = lambdify( F[1], [x1, x2, x3] )
  F2func_org  = lambdify( F[2], [x1, x2, x3] )
  F3func_org  = lambdify( F[3], [x1, x2, x3] )

  σ11func_org = lambdify( σ[1,1], [x1, x2, x3] )
  σ12func_org = lambdify( σ[1,2], [x1, x2, x3] )
  σ13func_org = lambdify( σ[1,3], [x1, x2, x3] )
  σ22func_org = lambdify( σ[2,2], [x1, x2, x3] )
  σ23func_org = lambdify( σ[2,3], [x1, x2, x3] )
  σ33func_org = lambdify( σ[3,3], [x1, x2, x3] )

  ϵ11func_org = lambdify( ϵ[1,1], [x1, x2, x3] )
  ϵ12func_org = lambdify( ϵ[1,2], [x1, x2, x3] )
  ϵ13func_org = lambdify( ϵ[1,3], [x1, x2, x3] )
  ϵ22func_org = lambdify( ϵ[2,2], [x1, x2, x3] )
  ϵ23func_org = lambdify( ϵ[2,3], [x1, x2, x3] )
  ϵ33func_org = lambdify( ϵ[3,3], [x1, x2, x3] )

  # make sure you can call it with matrices
  u1func = x -> u1func_org.(x[:,1], x[:,2], x[:,3])
  u2func = x -> u2func_org.(x[:,1], x[:,2], x[:,3])
  u3func = x -> u3func_org.(x[:,1], x[:,2], x[:,3])

  F1func = x -> F1func_org.(x[:,1], x[:,2], x[:,3])
  F2func = x -> F2func_org.(x[:,1], x[:,2], x[:,3])
  F3func = x -> F3func_org.(x[:,1], x[:,2], x[:,3])

  σ11func = x -> σ11func_org.(x[:,1], x[:,2], x[:,3])
  σ12func = x -> σ12func_org.(x[:,1], x[:,2], x[:,3])
  σ13func = x -> σ13func_org.(x[:,1], x[:,2], x[:,3])
  σ22func = x -> σ22func_org.(x[:,1], x[:,2], x[:,3])
  σ23func = x -> σ23func_org.(x[:,1], x[:,2], x[:,3])
  σ33func = x -> σ33func_org.(x[:,1], x[:,2], x[:,3])

  ϵ11func = x -> ϵ11func_org.(x[:,1], x[:,2], x[:,3])
  ϵ12func = x -> ϵ12func_org.(x[:,1], x[:,2], x[:,3])
  ϵ13func = x -> ϵ13func_org.(x[:,1], x[:,2], x[:,3])
  ϵ22func = x -> ϵ22func_org.(x[:,1], x[:,2], x[:,3])
  ϵ23func = x -> ϵ23func_org.(x[:,1], x[:,2], x[:,3])
  ϵ33func = x -> ϵ33func_org.(x[:,1], x[:,2], x[:,3])

  return( u1func, u2func, u3func,
          F1func, F2func, F3func,
          σ11func, σ12func, σ13func, σ22func, σ23func, σ33func,
          ϵ11func, ϵ12func, ϵ13func, ϵ22func, ϵ23func, ϵ33func )

end

( u1func, u2func, u3func,
  F1func, F2func, F3func,
  σ11func, σ12func, σ13func, σ22func, σ23func, σ33func,
  ϵ11func, ϵ12func, ϵ13func, ϵ22func, ϵ23func, ϵ33func ) = ExactSol3D( mat.Cstiff )

#   Setup boundary conditions
function funcB( p::Array{Float64} )
  return fill(0.0, size(p,1), 3)
end
bctype = [2,1,2,1,1,1]
funcTop = (p) -> [ σ13func(p)  σ23func(p)  σ33func(p)]
funcBot = (p) -> [-σ13func(p) -σ23func(p) -σ33func(p)]
#   Setup source function
source = (p) -> [F1func(p) F2func(p) F3func(p)]

## Compute error for different polynomial orders and grid sizes
#   Initialize arrays
Err_uh1 = fill( 0.0, length(Ps), length(Ns) )
Err_uh2 = fill( 0.0, length(Ps), length(Ns) )
Err_uh3 = fill( 0.0, length(Ps), length(Ns) )
Err_σh1 = fill( 0.0, length(Ps), length(Ns) )
Err_σh2 = fill( 0.0, length(Ps), length(Ns) )
Err_σh3 = fill( 0.0, length(Ps), length(Ns) )
Err_σh5 = fill( 0.0, length(Ps), length(Ns) )
Err_σh6 = fill( 0.0, length(Ps), length(Ns) )
Err_σh9 = fill( 0.0, length(Ps), length(Ns) )
# Err_ϵh1 = fill( 0.0, length(Ps), length(Ns) )
# Err_ϵh2 = fill( 0.0, length(Ps), length(Ns) )
# Err_ϵh3 = fill( 0.0, length(Ps), length(Ns) )
# Err_ϵh5 = fill( 0.0, length(Ps), length(Ns) )
# Err_ϵh6 = fill( 0.0, length(Ps), length(Ns) )
# Err_ϵh9 = fill( 0.0, length(Ps), length(Ns) )

#   Loop over polynomial order and grid size
for ii in 1:length(Ps), jj in 1:length(Ns)

  P = Ps[ii]; N = Ns[jj]

  @printf( " %6i %6i\n", P.p, N )

  mesh   = Mesh3D( "cube", P, N = N)
  master = Master3D( P )

  prob = Elas( @sprintf("Reg %i %i", P.p, N), mat, source, bctype, false, [funcBot, funcB, funcTop, funcB, funcB, funcB] )

  (uhath, uh, σh) = hdgSolve( master, mesh, prob )

  #   Initialize arrays
  err_uh = fill( 0.0, dim )
  err_σh = fill( 0.0, dim^2 )
  # err_ϵh = fill( 0.0, dim^2 )

  # preallocate
  jcwd = fill( 0.0, size(master.∇ϕ,2), size(master.∇ϕ,2) )
  ∂ξ∂x = fill( 0.0, size(master.∇ϕ,2), dim^2 )
  ∂x∂ξ = fill( 0.0, size(master.∇ϕ,2), dim^2 )

  for kk in 1:size(mesh.t,1)

    # Compute Jacobians
    luteos.compJacob!( master, mesh.nodes[:,:,kk], ∂ξ∂x, jcwd, ∂x∂ξ )

    # u
    Δuh1       = master.ϕ' * ( uh[:,1,kk] - u1func( mesh.nodes[:,:,kk] ) )
    Δuh2       = master.ϕ' * ( uh[:,2,kk] - u2func( mesh.nodes[:,:,kk] ) )
    Δuh3       = master.ϕ' * ( uh[:,3,kk] - u3func( mesh.nodes[:,:,kk] ) )
    err_uh[1] += Δuh1' * jcwd * Δuh1
    err_uh[2] += Δuh2' * jcwd * Δuh2
    err_uh[3] += Δuh3' * jcwd * Δuh3

    # σ
    Δσh1       = master.ϕ' * ( σh[:,1,kk] - σ11func( mesh.nodes[:,:,kk] ) )
    Δσh2       = master.ϕ' * ( σh[:,2,kk] - σ12func( mesh.nodes[:,:,kk] ) )
    Δσh3       = master.ϕ' * ( σh[:,3,kk] - σ13func( mesh.nodes[:,:,kk] ) )
    Δσh5       = master.ϕ' * ( σh[:,5,kk] - σ22func( mesh.nodes[:,:,kk] ) )
    Δσh6       = master.ϕ' * ( σh[:,6,kk] - σ23func( mesh.nodes[:,:,kk] ) )
    Δσh9       = master.ϕ' * ( σh[:,9,kk] - σ33func( mesh.nodes[:,:,kk] ) )
    err_σh[1] += Δσh1' * jcwd * Δσh1
    err_σh[2] += Δσh2' * jcwd * Δσh2
    err_σh[3] += Δσh3' * jcwd * Δσh3
    err_σh[5] += Δσh5' * jcwd * Δσh5
    err_σh[6] += Δσh6' * jcwd * Δσh6
    err_σh[9] += Δσh9' * jcwd * Δσh9

    # # ϵ
    # Δϵh1       = master.ϕ' * ( ϵh[:,1,kk] - ϵ11func( mesh.nodes[:,:,kk] ) )
    # Δϵh2       = master.ϕ' * ( ϵh[:,2,kk] - ϵ12func( mesh.nodes[:,:,kk] ) )
    # Δϵh3       = master.ϕ' * ( ϵh[:,3,kk] - ϵ13func( mesh.nodes[:,:,kk] ) )
    # Δϵh5       = master.ϕ' * ( ϵh[:,5,kk] - ϵ22func( mesh.nodes[:,:,kk] ) )
    # Δϵh6       = master.ϕ' * ( ϵh[:,6,kk] - ϵ23func( mesh.nodes[:,:,kk] ) )
    # Δϵh9       = master.ϕ' * ( ϵh[:,9,kk] - ϵ33func( mesh.nodes[:,:,kk] ) )
    # err_ϵh[1] += Δϵh1' * jcwd * Δϵh1
    # err_ϵh[2] += Δϵh2' * jcwd * Δϵh2
    # err_ϵh[3] += Δϵh3' * jcwd * Δϵh3
    # err_ϵh[5] += Δϵh5' * jcwd * Δϵh5
    # err_ϵh[6] += Δϵh6' * jcwd * Δϵh6
    # err_ϵh[9] += Δϵh9' * jcwd * Δϵh9

  end

  Err_uh1[ii,jj] = sqrt(err_uh[1])
  Err_uh2[ii,jj] = sqrt(err_uh[2])
  Err_uh3[ii,jj] = sqrt(err_uh[3])

  Err_σh1[ii,jj] = sqrt(err_σh[1])
  Err_σh2[ii,jj] = sqrt(err_σh[2])
  Err_σh3[ii,jj] = sqrt(err_σh[3])
  Err_σh5[ii,jj] = sqrt(err_σh[5])
  Err_σh6[ii,jj] = sqrt(err_σh[6])
  Err_σh9[ii,jj] = sqrt(err_σh[9])

  # Err_ϵh1[ii,jj] = sqrt(err_ϵh[1])
  # Err_ϵh2[ii,jj] = sqrt(err_ϵh[2])
  # Err_ϵh3[ii,jj] = sqrt(err_ϵh[3])
  # Err_ϵh5[ii,jj] = sqrt(err_ϵh[5])
  # Err_ϵh6[ii,jj] = sqrt(err_ϵh[6])
  # Err_ϵh9[ii,jj] = sqrt(err_ϵh[9])

end

# Compute convergence rates
h = 1 ./ ( Ns - 1 )
conv_uh1 = (log.( Err_uh1[:,end-1]) - log.( Err_uh1[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));
conv_uh2 = (log.( Err_uh2[:,end-1]) - log.( Err_uh2[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));
conv_uh3 = (log.( Err_uh3[:,end-1]) - log.( Err_uh3[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));

conv_σh1 = (log.( Err_σh1[:,end-1]) - log.( Err_σh1[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));
conv_σh2 = (log.( Err_σh2[:,end-1]) - log.( Err_σh2[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));
conv_σh3 = (log.( Err_σh3[:,end-1]) - log.( Err_σh3[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));
conv_σh5 = (log.( Err_σh5[:,end-1]) - log.( Err_σh5[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));
conv_σh6 = (log.( Err_σh6[:,end-1]) - log.( Err_σh6[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));
conv_σh9 = (log.( Err_σh9[:,end-1]) - log.( Err_σh9[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));

# conv_ϵh1 = (log.( Err_ϵh1[:,end-1]) - log.( Err_ϵh1[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));
# conv_ϵh2 = (log.( Err_ϵh2[:,end-1]) - log.( Err_ϵh2[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));
# conv_ϵh3 = (log.( Err_ϵh3[:,end-1]) - log.( Err_ϵh3[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));
# conv_ϵh5 = (log.( Err_ϵh5[:,end-1]) - log.( Err_ϵh5[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));
# conv_ϵh6 = (log.( Err_ϵh6[:,end-1]) - log.( Err_ϵh6[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));
# conv_ϵh9 = (log.( Err_ϵh9[:,end-1]) - log.( Err_ϵh9[:,end] ) ) / (log.( h[end-1]) - log.( h[end] ));

# Output to terminal
@printf("\n")
@printf("   Convergence rates for 3D Dirichlet + Neumann problem\n\n")
@printf("   ------------------------------------------\n\n")
@printf( "P    ")
for jj in 1:size(Ps,1)
  @printf( " %6i", Ps[jj].p )
end
@printf( "\n" )

#   u
@printf( "u₁   ")
for jj in 1:size(Ps,1)
  @printf( " %6.4f", conv_uh1[jj] )
end
@printf( "\n" )
@printf( "u₂   ")
for jj in 1:size(Ps,1)
  @printf( " %6.4f", conv_uh2[jj] )
end
@printf( "\n" )
@printf( "u₃   ")
for jj in 1:size(Ps,1)
  @printf( " %6.4f", conv_uh3[jj] )
end
@printf( "\n" )

#   σ
@printf( "σ₁₁  ")
for jj in 1:length(Ps)
  @printf( " %6.4f", conv_σh1[jj] )
end
@printf( "\n" )
@printf( "σ₁₂  ")
for jj in 1:length(Ps)
  @printf( " %6.4f", conv_σh2[jj] )
end
@printf( "\n" )
@printf( "σ₁₃  ")
for jj in 1:length(Ps)
  @printf( " %6.4f", conv_σh3[jj] )
end
@printf( "\n" )
@printf( "σ₂₂  ")
for jj in 1:length(Ps)
  @printf( " %6.4f", conv_σh5[jj] )
end
@printf( "\n" )
@printf( "σ₂₃  ")
for jj in 1:length(Ps)
  @printf( " %6.4f", conv_σh6[jj] )
end
@printf( "\n" )
@printf( "σ₃₃  ")
for jj in 1:length(Ps)
  @printf( " %6.4f", conv_σh9[jj] )
end
@printf( "\n" )

# #   ϵ
# @printf( "ϵ₁₁  ")
# for jj in 1:length(Ps)
#   @printf( " %6.4f", conv_ϵh1[jj] )
# end
# @printf( "\n" )
# @printf( "ϵ₁₂  ")
# for jj in 1:length(Ps)
#   @printf( " %6.4f", conv_ϵh2[jj] )
# end
# @printf( "\n" )
# @printf( "ϵ₁₃  ")
# for jj in 1:length(Ps)
#   @printf( " %6.4f", conv_ϵh3[jj] )
# end
# @printf( "\n" )
# @printf( "ϵ₂₂  ")
# for jj in 1:length(Ps)
#   @printf( " %6.4f", conv_ϵh5[jj] )
# end
# @printf( "\n" )
# @printf( "ϵ₂₃  ")
# for jj in 1:length(Ps)
#   @printf( " %6.4f", conv_ϵh6[jj] )
# end
# @printf( "\n" )
# @printf( "ϵ₃₃  ")
# for jj in 1:length(Ps)
#   @printf( " %6.4f", conv_ϵh9[jj] )
# end
# @printf( "\n" )

open("errors_Elas_Neumann3D.dat", "w") do f
  @printf(f, "P \t N \t E_uh1 \t E_uh2 \t E_uh3 \t E_σh1 \t E_σh2 \t E_σh3 \t E_σh5 \t E_σh6 \t E_σh9 \n")#\t E_ϵh1 \t E_ϵh2 \t E_ϵh3 \t E_ϵh5 \t E_ϵh6 \t E_ϵh9 \n")#\t E_J\n")
  for ii in 1:length(Ps), jj in 1:length(Ns)
    @printf(f, "%i \t %i \t %16.15e \t %16.15e \t %16.15e \t %16.15e \t %16.15e \t %16.15e \t %16.15e \t %16.15e \t %16.15e \n", #\t %16.15e \t %16.15e \t %16.15e \t %16.15e \t %16.15e \t %16.15e\n",
      Ps[ii].p, Ns[jj], Err_uh1[ii,jj], Err_uh2[ii,jj], Err_uh3[ii,jj],
      Err_σh1[ii,jj], Err_σh2[ii,jj], Err_σh3[ii,jj], Err_σh5[ii,jj], Err_σh6[ii,jj], Err_σh9[ii,jj] )#,
      # Err_ϵh1[ii,jj], Err_ϵh2[ii,jj], Err_ϵh3[ii,jj], Err_ϵh5[ii,jj], Err_ϵh6[ii,jj], Err_ϵh9[ii,jj] )
  end
end

# end # limit scope


text:
export BlockGP

###################### Implementation of AbstractGaussianProcess interface #################

"""
    BlockGP{Tfs<:AbstractVector{<:AbstractGP}} <: AbstractGaussianProcess

A collection of `AbstractGaussianProcess`es. Akin to the usual conception of a
Multi-Output Gaussian process.
"""
struct BlockGP{Tfs<:AV{<:AbstractGP}} <: AbstractGaussianProcess
    fs::Tfs
end
BlockGP(fs::Tuple) = BlockGP([fs...])

mean(f::BlockGP) = BlockMean(mean.(f.fs))
kernel(f::BlockGP) = BlockKernel(kernel.(f.fs), kernel.(f.fs, permutedims(f.fs)))
kernel(f1::BlockGP, f2::BlockGP) = BlockCrossKernel(kernel.(f1.fs, permutedims(f2.fs)))
kernel(f1::GP, f2::BlockGP) = BlockCrossKernel(kernel.(Ref(f1), permutedims(f2.fs)))
kernel(f1::BlockGP, f2::GP) = BlockCrossKernel(kernel.(f1.fs, Ref(f2)))

deconstruct(f::BlockGP) = (f.fs...,)

# Return a gpc object from one of the children.
function getproperty(f::BlockGP, d::Symbol)
    if d === :gpc
        return f.fs[1].gpc
    else
        return getfield(f, d)
    end
end

##################################### Syntactic Sugar ######################################

# Specialised implementation of `rand` for `BlockGP`s.
function rand(rng::AbstractRNG, f::BlockGP, N::Int)
    M = BlockArray(undef_blocks, AbstractMatrix{Float64}, length.(f.fs), [N])
    μ = mean_vec(f)
    for b in eachindex(f.fs)
        setblock!(M, getblock(μ, b) * ones(1, N), b, 1)
    end
    return M + chol(cov(f))' * BlockMatrix(randn.(Ref(rng), length.(f.fs), N))
end
rand(f::BlockGP, N::Int) = rand(Base.Random.GLOBAL_RNG, f, N)

function rand(rng::AbstractRNG, f::BlockGP)
    return mean_vec(f) + chol(cov(f))' * BlockVector(randn.(Ref(rng), length.(f.fs)))
end
rand(f::BlockGP) = rand(Base.Random.GLOBAL_RNG, f)

# Convenience methods for invoking `logpdf` and `rand` with multiple processes.
rand(rng::AbstractRNG, fs::AV{<:AbstractGP}) = vec.(rand(rng, fs, 1))
rand(rng::AbstractRNG, fs::AV{<:AbstractGP}, N::Int) = rand(rng, BlockGP(fs), N).blocks
rand(fs::AV{<:AbstractGP}, N::Int) = rand(Random.GLOBAL_RNG, fs, N)
rand(fs::AV{<:AbstractGP}) = vec.(rand(Random.GLOBAL_RNG, fs))
logpdf(fs::AV{<:AbstractGP}, ys::AV{<:AV{<:Real}}) = logpdf(BlockGP(fs), BlockVector(ys))

ext: jl


text:
# Double Pendulum
T = Float64
model = TrajectoryOptimization.Dynamics.doublependulum
model_d = rk3(model)
n = model.n; m = model.m

# costs
Qf = 100.0*Diagonal(I,n)
Q = (1.0e-2)*Diagonal(I,n)
R = (1.0e-2)*Diagonal(I,m)

x0 = [0.; 0.; 0.; 0.]
xf = [π; 0.; 0.; 0.];
N = 101
dt = 0.01
U0 = [ones(m) for k = 1:N-1]
obj = TrajectoryOptimization.LQRObjective(Q,R,Qf,xf,N)

doublependulum = TrajectoryOptimization.Problem(model_d, obj, x0=x0, xf=xf, N=N, dt=dt)
initial_controls!(doublependulum, U0)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
#=
--------------------------------------------------------------------------------------------------------------------------------
Base and generic math
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    isnothing(::Any)

True if the argument is ```nothing``` (false otherwise).
"""
isnothing(::Any) = false;
isnothing(::Nothing) = true;

"""
    verb_message(verb::Bool, message::String)

@info `message` if `verb` is true.
"""
verb_message(verb::Bool, message::String) = verb ? @info(message) : nothing;

"""
    check_bounds(X::Number, LB::Number, UB::Number)

Check whether `X` is larger or equal than `LB` and lower or equal than `UB`

    check_bounds(X::Number, LB::Number)

Check whether `X` is larger or equal than `LB`
"""
check_bounds(X::Number, LB::Number, UB::Number) = X < LB || X > UB ? throw(DomainError) : nothing
check_bounds(X::Number, LB::Number) = X < LB ? throw(DomainError) : nothing

"""
    error_info(err::Exception)
    error_info(err::RemoteException)

Return error main information
"""
error_info(err::Exception) = (err, err.msg, stacktrace(catch_backtrace()));
error_info(err::RemoteException) = (err.captured.ex, err.captured.ex.msg, [err.captured.processed_bt[i][1] for i=1:length(err.captured.processed_bt)]);

"""
    mean_skipmissing(X::AbstractArray{Float64,1})
    mean_skipmissing(X::AbstractArray{Union{Missing, Float64},1})

Compute the mean of the observed values in `X`.

    mean_skipmissing(X::AbstractArray{Float64})
    mean_skipmissing(X::AbstractArray{Union{Missing, Float64}})

Compute the mean of the observed values in `X` column wise.

# Examples
```jldoctest
julia> mean_skipmissing([1.0; missing; 3.0])
2.0

julia> mean_skipmissing([1.0 2.0; missing 3.0; 3.0 5.0])
3-element Array{Float64,1}:
 1.5
 3.0
 4.0
```
"""
mean_skipmissing(X::AbstractArray{Float64,1}) = mean(X);
mean_skipmissing(X::AbstractArray{Float64}) = mean(X, dims=2);
mean_skipmissing(X::AbstractArray{Union{Missing, Float64},1}) = mean(skipmissing(X));
mean_skipmissing(X::AbstractArray{Union{Missing, Float64}}) = vcat([mean_skipmissing(X[i,:]) for i=1:size(X,1)]...);

"""
    std_skipmissing(X::AbstractArray{Float64,1})
    std_skipmissing(X::AbstractArray{Union{Missing, Float64},1})

Compute the standard deviation of the observed values in `X`.

    std_skipmissing(X::AbstractArray{Float64})
    std_skipmissing(X::AbstractArray{Union{Missing, Float64}})

Compute the standard deviation of the observed values in `X` column wise.

# Examples
```jldoctest
julia> std_skipmissing([1.0; missing; 3.0])
1.4142135623730951

julia> std_skipmissing([1.0 2.0; missing 3.0; 3.0 5.0])
3-element Array{Float64,1}:
   0.7071067811865476
 NaN
   1.4142135623730951
```
"""
std_skipmissing(X::AbstractArray{Float64,1}) = std(X);
std_skipmissing(X::AbstractArray{Float64}) = std(X, dims=2);
std_skipmissing(X::AbstractArray{Union{Missing, Float64},1}) = std(skipmissing(X));
std_skipmissing(X::AbstractArray{Union{Missing, Float64}}) = vcat([std_skipmissing(X[i,:]) for i=1:size(X,1)]...);

"""
    is_vector_in_matrix(vect::AbstractVector, matr::AbstractMatrix)

Check whether the vector `vect` is included in the matrix `matr`.

# Examples
julia> is_vector_in_matrix([1;2], [1 2; 2 3])
true
"""
is_vector_in_matrix(vect::AbstractVector, matr::AbstractMatrix) = sum(sum(vect .== matr, dims=1) .== length(vect)) > 0;

"""
    isconverged(new::Float64, old::Float64, tol::Float64, ε::Float64, increasing::Bool)

Check whether `new` is close enough to `old`.

# Arguments
- `new`: new objective or loss
- `old`: old objective or loss
- `tol`: tolerance
- `ε`: small Float64
- `increasing`: true if `new` increases, at each iteration, with respect to `old`
"""
isconverged(new::Float64, old::Float64, tol::Float64, ε::Float64, increasing::Bool) = increasing ? (new-old)./(abs(old)+ε) <= tol : -(new-old)./(abs(old)+ε) <= tol;

"""
    soft_thresholding(z::Float64, ζ::Float64)

Soft thresholding operator.
"""
soft_thresholding(z::Float64, ζ::Float64) = sign(z)*max(abs(z)-ζ, 0);

#=
--------------------------------------------------------------------------------------------------------------------------------
Transformations
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    demean(X::Array{Float64,1})
    demean(X::JVector)

Demean data.

    demean(X::Array{Float64,2})
    demean(X::JArray)

Demean data.

# Examples
```jldoctest
julia> demean([1.0; 1.5; 2.0; 2.5; 3.0])
5-element Array{Float64,1}:
 -1.0
 -0.5
  0.0
  0.5
  1.0

julia> demean([1.0 3.5 1.5 4.0 2.0; 4.5 2.5 5.0 3.0 5.5])
2×5 Array{Float64,2}:
 -1.4   1.1  -0.9   1.6  -0.4
  0.4  -1.6   0.9  -1.1   1.4
```
"""
demean(X::Array{Float64,1}) = X .- mean(X);
demean(X::Array{Float64,2}) = X .- mean(X,dims=2);
demean(X::JVector) = X .- mean_skipmissing(X);
demean(X::JArray) = X .- mean_skipmissing(X);

#=
--------------------------------------------------------------------------------------------------------------------------------
Base: time series
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    interpolate(X::JArray{Float64}, n::Int64, T::Int64)

Interpolate each series in `X`, in turn, by replacing missing observations with the sample average of the non-missing values.

# Arguments
- `X`: observed measurements (`nxT`)
- `n` and `T` are the number of series and observations

    interpolate(X::Array{Float64}, n::Int64, T::Int64)

# Arguments
- `X`: observed measurements (`nxT`)
- `n` and `T` are the number of series and observations
"""
function interpolate(X::JArray{Float64}, n::Int64, T::Int64)
    data = copy(X);
    for i=1:n
        data[i, ismissing.(X[i, :])] .= mean_skipmissing(X[i, :]);
    end
    data = convert(Array{Float64}, data);
    return data;
end

interpolate(X::Array{Float64}, n::Int64, T::Int64) = X;

"""
    lag(X::Array, p::Int64)

Construct the data required to run a standard vector autoregression.

# Arguments
- `X`: observed measurements (`nxT`), where `n` and `T` are the number of series and observations.
- `p`: number of lags in the vector autoregression

# Output
- `X_{t}`
- `X_{t-1}`
"""
function lag(X::Array, p::Int64)

    # VAR(p) data
    X_t = X[:, 1+p:end];
    X_lagged = vcat([X[:, p-j+1:end-j] for j=1:p]...);

    # Return output
    return X_t, X_lagged;
end

"""
    companion_form(Ψ::Array{Float64,2}, Σ::SymMatrix)

Construct the companion form parameters of a VAR(p) with coefficients Ψ and var-cov matrix of the residuals Σ.
"""
function companion_form(Ψ::Array{Float64,2}, Σ::SymMatrix)

    # Dimensions
    n = size(Σ,2);
    p = Int64(size(Ψ,2)/n);
    np_1 = n*(p-1);

    # Companion form VAR(p)
    C = [Ψ; Matrix(I, np_1, np_1) zeros(np_1, n)];
    V = Symmetric([Σ zeros(n, np_1); zeros(np_1, n*p)])::SymMatrix;

    # Return output
    return C, V;
end

"""
    ext_companion_form(Ψ::Array{Float64,2}, Σ::SymMatrix)

Construct the companion form parameters of a VAR(p) with coefficients Ψ and var-cov matrix of the residuals Σ. The companion form is extend with additional n entries.
"""
function ext_companion_form(Ψ::Array{Float64,2}, Σ::SymMatrix)

    # Dimensions
    n = size(Σ,2);
    p = Int64(size(Ψ,2)/n);
    np = n*p;

    # Companion form VAR(p)
    C = [Ψ zeros(n, n); Matrix(I, np, np) zeros(np, n)];
    V = Symmetric([Σ zeros(n, np); zeros(np, np+n)])::SymMatrix;

    # Return output
    return C, V;
end


text:
#######
# useful concepts for dispatch
#######

"""
Advection terms

right now really only non-linear or ::Nothing
"""
abstract type AdvectionTerm end
struct NonLinearAdvectionTerm <: AdvectionTerm end

"""
Turbulence Closures

ways to handle drag and diffusion and such
"""
abstract type TurbulenceClosure end

struct LinearDrag{T} <: TurbulenceClosure
    λ::T
end

struct ConstantViscosity{T} <: TurbulenceClosure
    μ::T
    ν::T
    κ::T
    function ConstantViscosity{T}(;
        μ = T(1e-6),   # m²/s
        ν = T(1e-6),   # m²/s
        κ = T(1e-6),   # m²/s
    ) where {T <: AbstractFloat}
        return new{T}(μ, ν, κ)
    end
end

"""
Forcings

ways to add body terms and sources
"""
abstract type Forcing end
abstract type CoriolisForce <: Forcing end

struct fPlaneCoriolis{T} <: CoriolisForce
    fₒ::T
    β::T
    function fPlaneCoriolis{T}(;
        fₒ = T(1e-4), # Hz
        β = T(1e-11), # Hz/m
    ) where {T <: AbstractFloat}
        return new{T}(fₒ, β)
    end
end

struct SphereCoriolis{T} <: CoriolisForce
    Ω::T
    function SphereCoriolis{T}(;
        Ω = T(2π / 86400), # Hz
    ) where {T <: AbstractFloat}
        return new{T}(Ω)
    end
end
struct KinematicStress{T} <: Forcing
    τₒ::T
    function KinematicStress{T}(; τₒ = T(1e-4)) where {T <: AbstractFloat}
        return new{T}(τₒ)
    end
end

struct Buoyancy{T} <: Forcing
    α::T # 1/K
    g::T # m/s²
    function Buoyancy{T}(; α = T(2e-4), g = T(10)) where {T <: AbstractFloat}
        return new{T}(α, g)
    end
end

"""
Grouping structs
"""
abstract type AbstractModel end

Base.@kwdef struct SpatialModel{𝒜, ℬ, 𝒞, 𝒟, ℰ, ℱ} <: AbstractModel
    balance_law::𝒜
    physics::ℬ
    numerics::𝒞
    grid::𝒟
    boundary_conditions::ℰ
    parameters::ℱ
end

polynomialorders(model::SpatialModel) = convention(
    model.grid.resolution.polynomial_order,
    Val(ndims(model.grid.domain)),
)

abstract type ModelPhysics end

Base.@kwdef struct FluidPhysics{𝒪, 𝒜, 𝒟, 𝒞, ℬ} <: ModelPhysics
    orientation::𝒪 = ClimateMachine.Orientations.FlatOrientation()
    advection::𝒜 = NonLinearAdvectionTerm()
    dissipation::𝒟 = nothing
    coriolis::𝒞 = nothing
    buoyancy::ℬ = nothing
end

abstract type AbstractInitialValueProblem end

Base.@kwdef struct InitialValueProblem{𝒫, ℐ𝒱} <: AbstractInitialValueProblem
    params::𝒫 = nothing
    initial_conditions::ℐ𝒱 = nothing
end

abstract type AbstractSimulation end

struct Simulation{𝒜, ℬ, 𝒞, 𝒟, ℰ, ℱ} <: AbstractSimulation
    model::𝒜
    state::ℬ
    timestepper::𝒞
    initial_conditions::𝒟
    callbacks::ℰ
    time::ℱ
end

function Simulation(;
    model = nothing,
    state = nothing,
    timestepper = nothing,
    initial_conditions = nothing,
    callbacks = nothing,
    time = nothing,
)
    model = DGModel(model, initial_conditions = initial_conditions)

    FT = eltype(model.grid.vgeo)

    if state == nothing
        state = init_ode_state(model, FT(0); init_on_cpu = true)
    end
    # model = (discrete = dgmodel, spatial = model)
    return Simulation(
        model,
        state,
        timestepper,
        initial_conditions,
        callbacks,
        time,
    )
end

coordinates(simulation::Simulation) = coordinates(simulation.model.grid)
polynomialorders(simulation::Simulation) =
    polynomialorders(simulation.model.grid)

abstract type AbstractTimestepper end

Base.@kwdef struct TimeStepper{S, T} <: AbstractTimestepper
    method::S
    timestep::T
end

"""
calculate_dt(grid, wavespeed = nothing, diffusivity = nothing, viscocity = nothing, cfl = 0.1)
"""
function calculate_dt(
    grid;
    wavespeed = nothing,
    diffusivity = nothing,
    viscocity = nothing,
    cfl = 1.0,
)
    Δx = min_node_distance(grid)
    Δts = []
    if wavespeed != nothing
        push!(Δts, Δx / wavespeed)
    end
    if diffusivity != nothing
        push!(Δts, Δx^2 / diffusivity)
    end
    if viscocity != nothing
        push!(Δts, Δx^2 / viscocity)
    end
    if Δts == []
        @error("Please provide characteristic speed or diffusivities")
        return nothing
    end
    return cfl * minimum(Δts)
end

function calculate_dt(
    grid::DiscretizedDomain;
    wavespeed = nothing,
    diffusivity = nothing,
    viscocity = nothing,
    cfl = 1.0,
)
    return calculate_dt(
        grid.numerical;
        wavespeed = wavespeed,
        diffusivity = diffusivity,
        viscocity = viscocity,
        cfl = cfl,
    )
end


text:
ENV["GKSwstype"]="100" #src

# # GEMPIC
#
# ### Geometric ElectroMagnetic Particle-In-Cell Methods
#
# https://arxiv.org/abs/1609.03053
#
# Michael Kraus, Katharina Kormann, Philip J. Morrison, Eric Sonnendrücker
#
# Framework for Finite Element Particle-in-Cell methods based on
# the discretization of the underlying Hamiltonian structure of the
# Vlasov-Maxwell system.
#
# Install the GEMPIC package

using Pkg

Pkg.add("GEMPIC");

using ProgressMeter, Plots, GEMPIC

#md # ---

# # Strong Landau Damping
#
# The physical parameters

kx, α = 0.5, 0.5
xmin, xmax = 0, 2π/kx
domain = [xmin, xmax, xmax - xmin]

#md # --

# The numerical parameters

∆t = 0.05
nx = 32
n_particles = 100000
mesh = GEMPIC.Mesh( xmin, xmax, nx)
spline_degree = 3

#md # ---

# Initialize particles

mass, charge = 1.0, 1.0
particle_group = GEMPIC.ParticleGroup{1,2}( n_particles, mass, charge, 1)
sampler = LandauDamping( α, kx)

sample!(sampler, particle_group)

#md # --

# Particle-mesh coupling operators

kernel_smoother1 = ParticleMeshCoupling( domain, [nx], n_particles,
                                         spline_degree-1, :galerkin)

kernel_smoother0 = ParticleMeshCoupling( domain, [nx], n_particles,
                                         spline_degree, :galerkin)

# Allocate electrostatic fields and Maxwell solver

rho = zeros(Float64, nx)
efield_poisson = zeros(Float64, nx)

maxwell_solver = Maxwell1DFEM(domain, nx, spline_degree)

solve_poisson!( efield_poisson, particle_group,
                kernel_smoother0, maxwell_solver, rho)

#md # ---

# ### Charge density

xg = LinRange(xmin, xmax, nx)
sval = eval_uniform_periodic_spline_curve(spline_degree-1, rho)
plot( xg, sval, label="ρ")
#md savefig("rho.svg"); nothing # hide

#md # ![](rho.svg)

#md # ---

# ### Electric field

sval = eval_uniform_periodic_spline_curve(spline_degree-1, efield_poisson)
plot( xg, sval, label="efield")
#md savefig("ex.svg"); nothing # hide

#md # ![](ex.svg)

#md # ---

# Initialize the arrays for the spline coefficients of the fields
efield_dofs = [copy(efield_poisson), zeros(Float64, nx)]
bfield_dofs = zeros(Float64, nx)

propagator = HamiltonianSplitting( maxwell_solver,
                                   kernel_smoother0,
                                   kernel_smoother1,
                                   particle_group,
                                   efield_dofs,
                                   bfield_dofs,
                                   domain);

efield_dofs_n = propagator.e_dofs

thdiag = TimeHistoryDiagnostics( particle_group, maxwell_solver,
                        kernel_smoother0, kernel_smoother1 );

#md # ---

# ## Loop over time

steps, Δt = 500, 0.05

@showprogress 1 for j = 1:steps # loop over time

    ## Strang splitting
    strang_splitting!(propagator, Δt, 1)

    ## Diagnostics
    solve_poisson!( efield_poisson, particle_group,
                    kernel_smoother0, maxwell_solver, rho)

    write_step!(thdiag, j * Δt, spline_degree,
                    efield_dofs,  bfield_dofs,
                    efield_dofs_n, efield_poisson)

end

#md nothing # hide
#md # ---

# ## Diagnostics stored in a dataframe

using DataFrames
first(thdiag.data, 5)

#md # ---

plot(thdiag.data[!,:Time], log.(thdiag.data[!,:PotentialEnergyE1]))
#md savefig("mod_e.svg"); nothing # hide
#md # ![](mod_e.svg)

@test true #src

#md # Next  [Conclusion](/Numkin2019/06/build/index.html)


text:
"""
Input/output
"""
module Io

"""
    prettyarr(a::Array)

Get a string represantation of an array similar to the REPL output
"""
function prettyarr(a::Array)::String
    buf = IOBuffer()
    show(buf, MIME("text/plain"), a)
    read(seekstart(buf), String)
end

include("Io/parse.jl")
include("Io/dimacs.jl")
include("Io/dimacs10.jl")
include("Io/snap.jl")

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
""" [Nieuwstadt1993](@cite) """
Nieuwstadt_θ_liq_ice(::Type{FT}) where {FT} = z -> if z <= 1350.0
        FT(300)
    else
        FT(300) + 3 * (z - 1350) / 1000
    end

""" [Nieuwstadt1993](@cite) """
Nieuwstadt_u(::Type{FT}) where {FT} = z -> FT(0.01)

""" [Nieuwstadt1993](@cite) """
Nieuwstadt_tke(::Type{FT}) where {FT} = z -> if (z <= 1600.0)
        FT(0.1) * FT(1.46) * FT(1.46) * (1 - z / 1600)
    else
        FT(0)
    end


text:
@testset "reducedims" begin
    A = reshape(Vector(1:16), (4,4))
    Aindices = IArray(A, a = 2:5, b = [:one, :two, :three, :four])

    A1 = reduce(max, A, dims=1)
    A2 = reduce(max, A, dims=2)
    A3 = reduce(max, A, dims=:)

    Aind1 = reduce(max, Aindices, dims=:a)
    Aind2 = reduce(max, Aindices, dims=:b)
    Aind3 = reduce(max, Aindices, dims=:)

    @test A1 == Aind1
    @test A2 == Aind2
    @test A3 == Aind3
end

@testset "Base" begin
    A = [10 20; 31 40]
    Aindices = IArray(A, x = 2:3, y = [:one, :two])

    @testset "$f" for f in (sum, prod, maximum, minimum, extrema)
        @test f(A) == f(Aindices)
        @test f(Aindices; dims=:x) == f(Aindices; dims=1) == f(A; dims=1)

        #@test indexnames(f(Aindices; dims=:x)) == (:x, :y) == indexnames(f(Aindices; dims=1))
    end

    @testset "$f" for f in (cumsum, cumprod, sort)
        @test f(Aindices; dims=:x) == f(Aindices; dims=1) == f(A; dims=1)

        #@test indexnames(f(Aindices; dims=:x)) == (:x, :y) == indexnames(f(Aindices; dims=1))

        @test f([1, 4, 3]) == f(IArray([1, 4, 3], (:vec,)))
#        @test_throws UndefKeywordError f(nda)
#        @test_throws UndefKeywordError f(a)
    end

    #= TODO sort
    @testset "sort!" begin
        A = [1 9; 7 3]
        Anamed = IArray(A, (:x, :y))

        # Vector case
        veca = [1, 9, 7, 3]
        sort!(NamedDimsArray(veca, :vec); order=Base.Reverse)
        @test issorted(veca; order=Base.Reverse)

        # Higher-dim case: `dims` keyword in `sort!` requires Julia v1.1+
        if VERSION > v"1.1-"
            sort!(Anamed, dims=:y)
            @test issorted(a[2, :])
            @test_throws UndefKeywordError sort!(Anamed)

            sort!(Anamed; dims=:x, order=Base.Reverse)
            @test issorted(a[:, 1]; order=Base.Reverse)
        end
    end
    =#

    #= TODO eachslice
    @testset "eachslice" begin
        slices = [[111 121; 211 221], [112 122; 212 222]]
        A = cat(slices...; dims=3)
        Anamed = IArrayArray(A, (:a, :b, :c))

        @test (
            sum(eachslice(Anamed; dims=:c)) ==
            sum(eachslice(Anamed; dims=3)) ==
            sum(eachslice(A; dims=3)) ==
            slices[1] + slices[2]
        )
#        @test_throws ArgumentError eachslice(Anamed; dims=(1, 2))
#        @test_throws ArgumentError eachslice(a; dims=(1, 2))

#        @test_throws UndefKeywordError eachslice(Anamed)
#        @test_throws UndefKeywordError eachslice(A)

        @test (
            indexnames(first(eachslice(Anamed; dims=:b))) ==
            indexnames(first(eachslice(Anamed; dims=2))) ==
            (:a, :c)
        )
    end
    =#

    @testset "mapslices" begin
        A = [10 20; 31 40]
        Anamed = IArray(A, (:x, :y))

        @test (
            mapslices(join, Anamed; dims=:x) ==
            mapslices(join, Anamed; dims=1) ==
            ["1031" "2040"]
        )
        @test (
            mapslices(join, Anamed; dims=:y) ==
            mapslices(join, Anamed; dims=2) ==
            reshape(["1020", "3140"], Val(2))
        )
        @test (
            mapslices(join, Anamed; dims=(:x, :y)) ==
            mapslices(join, Anamed; dims=(1, 2)) ==
            reshape(["10312040"], (1, 1))
        )
#        @test_throws UndefKeywordError mapslices(join, nda)
#        @test_throws UndefKeywordError mapslices(join, a)

        #=@test (
            indexnames(mapslices(join, Anamed; dims=:y)) ==
            indexnames(mapslices(join, Anamed; dims=2)) ==
            (:x, :y)
        )
        =#
    end

    @testset "mapreduce" begin
        A = [10 20; 31 40]
        Aindices = IArray(A, (:x, :y))

        @test mapreduce(isodd, |, Aindices) == true == mapreduce(isodd, |, A)
        @test (mapreduce(isodd, |, Aindices; dims=:x) ==
               mapreduce(isodd, |, Aindices; dims=1) ==
               [true false])
        @test (mapreduce(isodd, |, Aindices; dims=:y) ==
               mapreduce(isodd, |, Aindices; dims=2) ==
               [false true]')
        #=@test (indexnames(mapreduce(isodd, |, Aindices; dims=:y)) ==
               indexnames(mapreduce(isodd, |, Aindices; dims=2)) ==
               (:x, :y))
        =#
    end

    @testset "zero" begin
        A = [10 20; 31 40]
        Aindices = IArray(A, (:x, :y))

        @test zero(Aindices) == [0 0; 0 0] == zero(A)
        #@test indexnames(zero(Aindices)) == (:x, :y)
    end

    @testset "count" begin
        A = [true false; true true]
        Aindices = IArray(A, (:x, :y))

        @test count(Aindices) == count(A) == 3
#        @test_throws ErrorException count(nda; dims=:x)
#        @test_throws ErrorException count(a; dims=1)
    end
end  # Base

@testset "Statistics" begin
    A = [10 20; 30 40]
    Aindices = IArray(A, (:x, :y))
    @testset "$f" for f in (mean, std, var, median)
        @test f(Aindices) == f(A)
        @test f(Aindices; dims=:x) == f(Aindices; dims=1) == f(A; dims=1)

        #@test indexnames(f(Aindices; dims=:x)) == (:x, :y) == indexnames(f(Aindices; dims=1))
    end
end


text:
function run_repeats()
    for seed in 1:100
        # create .sh file to submit jobs
        filename = "sim$seed.sh"
        open(filename, "w") do io
            println(io, "#!/bin/bash")
            println(io, "#")
            println(io, "#SBATCH --job-name=sim$seed")
            println(io, "#")
            println(io, "#SBATCH --time=1:00:00")
            println(io, "#SBATCH --cpus-per-task=8")
            println(io, "#SBATCH --mem-per-cpu=2G")
            println(io, "#SBATCH --partition=owners")
            println(io, "")
            println(io, "#save job info on joblog:")
            println(io, "echo \"Job \$JOB_ID started on:   \" `hostname -s`")
            println(io, "echo \"Job \$JOB_ID started on:   \" `date `")
            println(io, "")
            println(io, "# load the job environment:")
            println(io, "module load julia")
            println(io, "")
            println(io, "# run code")
            println(io, "export JULIA_NUM_THREADS=8")
            println(io, "echo 'julia extra_sim.jl $seed'")
            println(io, "julia extra_sim.jl $seed")
            println(io, "")
            println(io, "#echo job info on joblog:")
            println(io, "echo \"Job \$JOB_ID ended on:   \" `hostname -s`")
            println(io, "echo \"Job \$JOB_ID ended on:   \" `date `")
            println(io, "#echo \" \"")
        end

        # submit job
        run(`sbatch $filename`)
        println("submitted job $seed")
        rm(filename, force=true)
        # sleep(1)
    end
end
run_repeats()


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to string processing in HTML backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# This was adapted from Julia `escape_string` function. In case of HTML, we
# should replace some special characters like `<` and `>`.
#
# If `replace_newline` is `true`, then `\n` is replaced with `<br>`. Otherwise,
# it is escaped, leading to `\\n`.
function _str_html_escaped(
    io::IO,
    s::AbstractString,
    replace_newline::Bool = false,
    escape_html_chars::Bool = true,
)
    a = Iterators.Stateful(s)
    for c in a
        if isascii(c)
            c == '\n'          ? (replace_newline ? print(io, "<BR>") : print(io, "\\n")) :
            c == '&'           ? (escape_html_chars ? print(io, "&")  : print(io, c)) :
            c == '<'           ? (escape_html_chars ? print(io, "<")   : print(io, c)) :
            c == '>'           ? (escape_html_chars ? print(io, ">")   : print(io, c)) :
            c == '"'           ? (escape_html_chars ? print(io, """) : print(io, c)) :
            c == '\''          ? (escape_html_chars ? print(io, "'") : print(io, c)) :
            c == '\0'          ? print(io, escape_nul(peek(a))) :
            c == '\e'          ? print(io, "\\e") :
            c == '\\'          ? print(io, "\\\\") :
            '\a' <= c <= '\r'  ? print(io, '\\', "abtnvfr"[Int(c)-6]) :
            c == '%'           ? print(io, "\\%") :
            isprint(c)         ? print(io, c) :
                                 print(io, "\\x", string(UInt32(c), base = 16, pad = 2))
        elseif !isoverlong(c) && !ismalformed(c)
            isprint(c)         ? print(io, c) :
            c <= '\x7f'        ? print(io, "\\x", string(UInt32(c), base = 16, pad = 2)) :
            c <= '\uffff'      ? print(io, "\\u", string(UInt32(c), base = 16, pad = Base.need_full_hex(peek(a)) ? 4 : 2)) :
                                 print(io, "\\U", string(UInt32(c), base = 16, pad = Base.need_full_hex(peek(a)) ? 8 : 4))
        else # malformed or overlong
            u = bswap(reinterpret(UInt32, c))
            while true
                print(io, "\\x", string(u % UInt8, base = 16, pad = 2))
                (u >>= 8) == 0 && break
            end
        end
    end
end

function _str_html_escaped(
    s::AbstractString,
    replace_newline::Bool = false,
    escape_html_chars::Bool = true
)
    return sprint(
        _str_html_escaped,
        s,
        replace_newline,
        escape_html_chars;
        sizehint = lastindex(s)
    )
end


text:
const c = cosmology()

"""
    li_prior_wt(m1, m2, z)

Returns the LALInference prior in m1, m2, and z.

```math
p\\left( m_1, m_2, z \\right) \\propto \\left( 1 + z \\right)^2 d_L^2(z) \\frac{\\partial d_L}{\\partial z}
```
"""
function li_prior_wt(m1, m2, z)
    dL = ustrip(u"Gpc", luminosity_dist(c, z))
    dC = dL / (1+z)
    dH = ustrip(u"Gpc", hubble_dist(c, 0))
    ddLdz = (dC + (1+z)*dH/Cosmology.E(c, z))

    (1 + z)*(1 + z)*dL*dL*ddLdz
end

"""
    md_sfr(z[, lambda, z_p, kappa])

Returns the Madau-Dickinson SFR.

```math
\\frac{\\mathrm{d} N}{\\mathrm{d} V \\mathrm{d} t} \\propto \\frac{\\left( 1 + z \\right)^\\lambda}{1 + \\left( \\frac{1+z}{1+z_p}\\right)^\\kappa}
```

The default values are `(lambda, z_p, kappa) = (2.7, 1.9, 5.6)`.
"""
function md_sfr(z, lambda=2.7, z_p=1.9, kappa=5.6)
    zp1 = 1 + z
    zp1^lambda / (1 + (zp1 / (1 + z_p))^kappa)
end

"""
    md_sfr_zwt(z[, lambda, z_p, kappa])

Returns the weight in `z` corresponding to the Madau-Dickinson SFR.

```math
\\frac{\\mathrm{d} N}{\\mathrm{d} V \\mathrm{d} t} \\propto \\frac{\\left( 1 + z \\right)^\\lambda}{1 + \\left( \\frac{1+z}{1+z_p}\\right)^\\kappa}
```

The default values are `(lambda, z_p, kappa) = (2.7, 1.9, 5.6)`.
"""
function md_sfr_zwt(z, lambda=2.7, z_p=1.9, kappa=5.6)
    zp1 = 1 + z
    wt = zp1^lambda / (1 + (zp1 / (1 + z_p))^kappa)

    4*pi*md_sfr(z, lambda, z_p, kappa) * ustrip(u"Gpc^3", comoving_volume_element(c, z)) / zp1
end

"""
    read_selection(file)

Read the injection set from `file` and return `(m1s, m2s, zs, pdraw, Ndraw)` for sensitivity estimation.
"""
function read_selection(file)
    h5open(file, "r") do f
        i = f["injections"]
        m1s = convert(Vector{Float64}, read(i, "mass1_source"))
        m2s = convert(Vector{Float64}, read(i, "mass2_source"))
        zs = convert(Vector{Float64}, read(i, "redshift"))

        T = read(attributes(f), "analysis_time_s") / (356.25*24.0*3600.0) # yr

        pm1m2z = convert(Vector{Float64}, read(i, "mass1_source_mass2_source_sampling_pdf") .* read(i, "redshift_sampling_pdf"))
        pdraw = pm1m2z / T

        sel_flag = (read(i, "ifar_cwb") .> 1) .| (read(i, "ifar_gstlal") .> 1) .| (read(i, "ifar_mbta") .> 1) .| (read(i, "ifar_pycbc_bbh") .> 1) .| (read(i, "ifar_pycbc_hyperbank") .> 1)

        Ndraw = round(Int, read(attributes(i)["total_generated"]))

        m1s[sel_flag], m2s[sel_flag], zs[sel_flag], pdraw[sel_flag], Ndraw
    end
end

"""
    resample_selection(log_dN, m1s, m2s, pdraws, Ndraw)

Returns re-sampled `(m1s, m2s, pdraws, Ndraw)` as if injections had been drawn
from the given fiducial population.
"""
function resample_selection(log_dN, m1s, m2s, pdraws, Ndraw)
    pop_wts = exp.(log_dN.(m1s, m2s))
    wts = pop_wts ./ pdraws
    sum_wts = sum(wts)
    sum_wts2 = sum(wts.*wts)
    norm = sum_wts / Ndraw
    Neff_norm = round(Int, sum_wts*sum_wts / sum_wts2)

    inds = sample(1:length(m1s), weights(wts), Neff_norm)
    norm_pop_wts = pop_wts ./ norm
    m1s[inds], m2s[inds], norm_pop_wts[inds], Neff_norm
end


text:
using Docile, Lexicon, RCall

const api_directory = "api"
const modules = Docile.Collector.submodules(RCall)

cd(dirname(@__FILE__)) do

    # Run the doctests *before* we start to generate *any* documentation.
    for m in modules
        failures = failed(doctest(m))
        if !isempty(failures.results)
            println("\nDoctests failed, aborting commit.\n")
            display(failures)
            exit(1) # Bail when doctests fail.
        end
    end

    # Generate and save the contents of docstrings as markdown files.
    index  = Index()
    config = Config(md_subheader = :category, category_order = [:module,    :function, :method, :type,
                                                                :typealias, :macro,    :global])
    for mod in modules
        update!(index, save(joinpath(api_directory, "$(mod).md"), mod, config))
    end
    save(joinpath(api_directory, "index.md"), index, config)

end


text:
using Test, simgym
import simgym: reproduce!

testpos = [0.,0.]
testvel = [.0001,.0001]
testdeath = 2

cell = Cell[Cell(copy(testpos),copy(testvel),testdeath,30)]
s = Simulation(cell,len_jones)
s.steptotal .= 2
reproduce!(s,200)

@test length(s.cells) == 3

max_stars_repo_name: andi-bergeson/simgym

max_issues_repo_name: andi-bergeson/simgym

max_forks_repo_name: andi-bergeson/simgym


text:
@inline dot_product_score(q, k) = scaled_dot_product_score(q, k, true)

@inline function scaled_dot_product_score(q, k, s = sqrt(inv(size(k, 1))))
    return matmul(batched_transpose(k), q, s)
end

masked_score(mask::Union{AbstractAttenMaskOp, AbstractAttenMask, Nothing}) = masked_score $ mask
@inline masked_score(::Nothing, score, args...) = score(args...)
@inline masked_score(::AbstractAttenMaskOp, ::Nothing, score, args...) = score(args...)
@inline masked_score(mask::AbstractAttenMask, score, args...) = masked_score(MaskOp(mask), mask, score, args...)
@inline function masked_score(maskop::AbstractAttenMaskOp, mask::AbstractAttenMask, score, args...)
    return _collapsed_nonbatch_call(score(args...)) do sc
        apply_mask(maskop, mask, sc)
    end
end

normalized_score(norm) = normalized_score $ norm
@inline function normalized_score(norm, score, args...)
    return _collapsed_call(norm, score(args...))
end

@inline attention_score(f, args...) = f(args...)


text:
"""
    time-varying LQR
"""
function tvlqr(A, B, Q, R)
    T = length(Q)

    P = [zero(A[1]) for t = 1:T]
    K = [zero(B[1]') for t = 1:T-1]
    P[T] = Q[T]

    for t = T-1:-1:1
        K[t] = (R[t] + B[t]' * P[t+1] *  B[t]) \ (B[t]' * P[t+1] * A[t])
        P[t] = (Q[t] + K[t]' * R[t] * K[t]
                + (A[t] - B[t] * K[t])' * P[t+1] * (A[t] - B[t] * K[t]))
    end

    return K, P
end

function tvlqr(model::Model{<: Integration, FixedTime}, x, u, h, Q, R)
    A, B = jacobians(model, x, u, h)
    K, P = tvlqr(A, B, Q, R)
    return K, P
end

function tvlqr(model::Model{<: Integration, FreeTime}, x, u, h, Q, R)
    A, B = jacobians(model, x, u, h)
    K, P = tvlqr(A, B, Q, [r[1:end-1, 1:end-1] for r in R])
    return K, P
end

"""
    jacobians along trajectory
"""
function _jacobians(model, x, u, h)
    A = []
    B = []

    w = zeros(model.d)

    T = length(x)

    for t = 1:T-1
        xt = x[t]
        ut = u[t]
        xt⁺ = x[t+1]

        fx(z) = fd(model, xt⁺, z, ut, w, h, t)
        fu(z) = fd(model, xt⁺, xt, z, w, h, t)
        fx⁺(z) = fd(model, z, xt, ut, w, h, t)

        A⁺ = ForwardDiff.jacobian(fx⁺, xt⁺)
        push!(A, -1.0 * A⁺ \ ForwardDiff.jacobian(fx, xt))
        push!(B, -1.0 * A⁺ \ ForwardDiff.jacobian(fu, ut))
    end

    return A, B
end

jacobians(model::Model{<: Integration, FixedTime}, x, u, h) = _jacobians(model, x, u, h)

function jacobians(model::Model{<: Integration, FreeTime}, x, u, h)
    A, B = _jacobians(model, x, u, h)
    T = length(x)
    return A, [B[t][:, 1:end-1] for t = 1:T-1]
end

"""
    projection
        use time-varying LQR to get dynamically feasible solution
"""
function lqr_projection(model::Model{<: Integration, FixedTime}, x̄, ū, h̄, Q, R)
    K, P = tvlqr(model, x̄, ū, h̄, Q, R)

    x_proj = [copy(x̄[1])]
    u_proj = []

    T = length(x̄)

    for t = 1:T-1
        push!(u_proj, ū[t] - K[t] * (x_proj[end] - x̄[t]))
        push!(x_proj, propagate_dynamics(model, x_proj[end], u_proj[end],
            zeros(model.d), h̄[1], t, tol_r = 1.0e-12, tol_d = 1.0e-12))
    end

    return x_proj, u_proj
end

function lqr_projection(model::Model{<: Integration, FreeTime}, x̄, ū, h̄, Q, R)
    K, P = tvlqr(model, x̄, ū, h̄, Q, R)

    x_proj = [copy(x̄[1])]
    u_proj = []

    T = length(x̄)

    for t = 1:T-1
        push!(u_proj, [ū[t][1:end-1] - K[t] * (x_proj[end] - x̄[t]); h̄[1]])
        push!(x_proj, propagate_dynamics(model, x_proj[end], u_proj[end],
            zeros(model.d), h̄[1], t, tol_r = 1.0e-12, tol_d = 1.0e-12))
    end

    return x_proj, u_proj
end


text:
using ModelingToolkit, Symbolics, Test
using ModelingToolkit: get_namespace, has_var, inputs, outputs, is_bound, bound_inputs,
                       unbound_inputs, bound_outputs, unbound_outputs, isinput, isoutput

# Test input handling
@parameters tv
D = Differential(tv)
@variables x(tv) u(tv) [input = true]
@test isinput(u)

@named sys = ODESystem([D(x) ~ -x + u], tv) # both u and x are unbound
@named sys2 = ODESystem([D(x) ~ -sys.x], tv, systems = [sys]) # this binds sys.x in the context of sys2, sys2.x is still unbound
@named sys3 = ODESystem([D(x) ~ -sys.x + sys.u], tv, systems = [sys]) # This binds both sys.x and sys.u

@named sys4 = ODESystem([D(x) ~ -sys.x, u ~ sys.u], tv, systems = [sys]) # This binds both sys.x and sys3.u, this system is one layer deeper than the previous. u is directly forwarded to sys.u, and in this case sys.u is bound while u is not

@test has_var(x ~ 1, x)
@test has_var(1 ~ x, x)
@test has_var(x + x, x)
@test !has_var(2 ~ 1, x)

@test get_namespace(x) == ""
@test get_namespace(sys.x) == "sys"
@test get_namespace(sys2.x) == "sys2"
@test get_namespace(sys2.sys.x) == "sys2₊sys"

@test !is_bound(sys, u)
@test !is_bound(sys, x)
@test !is_bound(sys, sys.u)
@test is_bound(sys2, sys.x)
@test !is_bound(sys2, sys.u)
@test !is_bound(sys2, sys2.sys.u)

fsys2 = flatten(sys2)
@test is_bound(fsys2, sys.x)
@test !is_bound(fsys2, sys.u)
@test !is_bound(fsys2, sys2.sys.u)

@test is_bound(sys3, sys.u) # I would like to write sys3.sys.u here but that's not how the variable is stored in the equations
@test is_bound(sys3, sys.x)

@test is_bound(sys4, sys.u)
@test !is_bound(sys4, u)

fsys4 = flatten(sys4)
@test is_bound(fsys4, sys.u)
@test !is_bound(fsys4, u)

@test isequal(inputs(sys), [u])
@test isequal(inputs(sys2), [sys.u])

@test isempty(bound_inputs(sys))
@test isequal(unbound_inputs(sys), [u])

@test isempty(bound_inputs(sys2))
@test isempty(bound_inputs(fsys2))
@test isequal(unbound_inputs(sys2), [sys.u])
@test isequal(unbound_inputs(fsys2), [sys.u])

@test isequal(bound_inputs(sys3), [sys.u])
@test isempty(unbound_inputs(sys3))

# Test output handling
@parameters tv
D = Differential(tv)
@variables x(tv) y(tv) [output = true]
@test isoutput(y)
@named sys = ODESystem([D(x) ~ -x, y ~ x], tv) # both y and x are unbound
syss = structural_simplify(sys) # This makes y an observed variable

@named sys2 = ODESystem([D(x) ~ -sys.x, y ~ sys.y], tv, systems = [sys])

@test !is_bound(sys, y)
@test !is_bound(sys, x)
@test !is_bound(sys, sys.y)

@test !is_bound(syss, y)
@test !is_bound(syss, x)
@test !is_bound(syss, sys.y)

@test isequal(unbound_outputs(sys), [y])
@test isequal(unbound_outputs(syss), [y])

@test isequal(unbound_outputs(sys2), [y])
@test isequal(bound_outputs(sys2), [sys.y])

syss = structural_simplify(sys2)

@test !is_bound(syss, y)
@test !is_bound(syss, x)
@test is_bound(syss, sys.y)

@test isequal(unbound_outputs(syss), [y])
@test isequal(bound_outputs(syss), [sys.y])

## Code generation with unbound inputs

@variables t x(t)=0 u(t)=0 [input = true]
D = Differential(t)
eqs = [
    D(x) ~ -x + u,
]

@named sys = ODESystem(eqs)
f, dvs, ps = ModelingToolkit.generate_control_function(sys, expression = Val{false},
                                                       simplify = true)

@test isequal(dvs[], x)
@test isempty(ps)

p = []
x = [rand()]
u = [rand()]
@test f[1](x, u, p, 1) == -x + u

# more complicated system

@variables u(t) [input = true]

function Mass(; name, m = 1.0, p = 0, v = 0)
    @variables y(t) [output = true]
    ps = @parameters m = m
    sts = @variables pos(t)=p vel(t)=v
    eqs = [D(pos) ~ vel
           y ~ pos]
    ODESystem(eqs, t, [pos, vel], ps; name)
end

function Spring(; name, k = 1e4)
    ps = @parameters k = k
    @variables x(t) = 0 # Spring deflection
    ODESystem(Equation[], t, [x], ps; name)
end

function Damper(; name, c = 10)
    ps = @parameters c = c
    @variables vel(t) = 0
    ODESystem(Equation[], t, [vel], ps; name)
end

function SpringDamper(; name, k = false, c = false)
    spring = Spring(; name = :spring, k)
    damper = Damper(; name = :damper, c)
    compose(ODESystem(Equation[], t; name),
            spring, damper)
end

connect_sd(sd, m1, m2) = [sd.spring.x ~ m1.pos - m2.pos, sd.damper.vel ~ m1.vel - m2.vel]
sd_force(sd) = -sd.spring.k * sd.spring.x - sd.damper.c * sd.damper.vel

# Parameters
m1 = 1
m2 = 1
k = 1000
c = 10

@named mass1 = Mass(; m = m1)
@named mass2 = Mass(; m = m2)
@named sd = SpringDamper(; k, c)

eqs = [connect_sd(sd, mass1, mass2)
       D(mass1.vel) ~ (sd_force(sd) + u) / mass1.m
       D(mass2.vel) ~ (-sd_force(sd)) / mass2.m]
@named _model = ODESystem(eqs, t)
@named model = compose(_model, mass1, mass2, sd);

f, dvs, ps = ModelingToolkit.generate_control_function(model, expression = Val{false},
                                                       simplify = true)
@test length(dvs) == 4
@test length(ps) == length(parameters(model))
p = ModelingToolkit.varmap_to_vars(ModelingToolkit.defaults(model), ps)
x = ModelingToolkit.varmap_to_vars(ModelingToolkit.defaults(model), dvs)
u = [rand()]
@test f[1](x, u, p, 1) == [u; 0; 0; 0]

@parameters t
@variables x(t) u(t) [input = true]
eqs = [Differential(t)(x) ~ u]
@named sys = ODESystem(eqs, t)
structural_simplify(sys)


text:
using ArgParse
using BSON
using ConvHMM
using Distributed
using ExpFamilyDistributions
using Glob
using MarkovModels
using NaturalSort

# BUT specific package to work with the SGE cluster
# See https://github.com/BUTSpeechFIT/BUTSGEManager
using BUTSGEManager

function getargs()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--epochs", "-e"
            help = "number of epochs to train"
            arg_type = Int
            default = 1
        "--init", "-i"
            help = "initial model to compute the alignments"
            arg_type = String
            default = ""
        "--update-prec-rate", "-u"
            help = "update rate of the precision parameter"
            arg_type = Int
            default = 10
        "--jobs", "-j"
            help = "number of jobs to run in parallel"
            arg_type = Int
            default = 4
        "--jobs-args", "-a"
            help = "arguments to pass to the parallel engine"
            arg_type = String
            default = ""
        "uttids"
            help = "list of utterance ids to train on"
            arg_type = String
            required = true
        "feadir"
            help = "features directory"
            arg_type = String
            required = true
        "alidir"
            help = "alignment graphs directory"
            arg_type = String
            required = true
        "modeldir"
            help = "directory where was created the initial model"
            arg_type = String
            required = true
    end
    args = parse_args(s)
end

const args = getargs()
const initfile = args["init"]
const updaterate = args["update-prec-rate"]
const uttids = readlines(args["uttids"])
const feadir = args["feadir"]
const alidir = args["alidir"]
const outdir = args["modeldir"]
const epochs = args["epochs"]

const init = initfile == "" ? nothing : ConvHMM.load(initfile)

# Make sure the directory exists
run(`mkdir -p $outdir`)

########################################################################
# Setup the workers
@info "initializing the workers"
addprocs(SGEManager(args["jobs"]), args = args["jobs-args"],
         exeflags = "--project=$(Base.active_project())")

@everywhere using BSON
@everywhere using ExpFamilyDistributions
@everywhere using MarkovModels
@everywhere using ConvHMM
@everywhere const args = $args
@everywhere const feadir = $feadir
@everywhere const alidir = $alidir
@everywhere const init = $init
@everywhere const updaterate = $updaterate

########################################################################
# Functions to accumulate statistics and to update the parameters

# Sum dictionaries fo statistics
@everywhere function mergestats(
    s1::Dict{Int64, Dict{Int64, Vector{Float64}}},
    s2::Dict{Int64, Dict{Int64, Vector{Float64}}}
)
    retval = Dict{Int64, Dict{Int64, Vector{Float64}}}()
    for d in keys(s1)
        retval[d] = Dict{Int64, Vector{Float64}}()
        for s in keys(s1[d])
            retval[d][s] = s1[d][s] + s2[d][s]
        end
   end
   retval
end

@everywhere function accumulate_h(
    model::Vector{Vector{ARNormal1D{K}}},
    X::Matrix{T},
    γ::Matrix{T}
) where {K, T <: AbstractFloat}

    D = length(model)
    S = length(model[1])

    accstats = Dict{Int64, Dict{Int64, Vector{Float64}}}()
    for d in 1:D
        accstats[d] = Dict{Int64, Vector{Float64}}()
        for (s, stats) in zip(1:S, accstats_h(model[d], X[d, :], γ))
            accstats[d][s] = stats
        end
    end
    accstats
end

@everywhere function accumulate_λ(
    model::Vector{Vector{ARNormal1D{K}}},
    X::Matrix{T},
    γ::Matrix{T}
) where {K, T <: AbstractFloat}

    D = length(model)
    S = length(model[1])

    accstats = Dict{Int64, Dict{Int64, Vector{Float64}}}()
    for d in 1:D
        accstats[d] = Dict{Int64, Vector{Float64}}()
        for (s, stats) in zip(1:S, accstats_λ(model[d], X[d, :], γ))
            accstats[d][s] = stats
        end
    end
    accstats
end

function update_h!(
    model::Vector{Vector{ARNormal1D{K}}},
    stats::Dict{Int64, Dict{Int64, Vector{T}}},
) where {K, T <: AbstractFloat}

    D = length(model)
    S = length(model[1])
    for d in 1:D
        for s in 1:S
            m = model[d][s]
            s = stats[d][s]
            η₀ = naturalparam(m.hprior)
            update!(m.hposterior, η₀ .+ s)
        end
    end
end

function update_λ!(
    model::Vector{Vector{ARNormal1D{K}}},
    stats::Dict{Int64, Dict{Int64, Vector{T}}},
) where {K, T <: AbstractFloat}

    D = length(model)
    S = length(model[1])
    for d in 1:D
        for s in 1:S
            m = model[d][s]
            s = stats[d][s]
            η₀ = naturalparam(m.λprior)
            update!(m.λposterior, η₀ .+ s)
        end
    end
end

########################################################################
# Training
#
# NOTE: the training is an EM-like algorithm. For the M-step, we
# alternately update the filters (h) or the precision parameters (λ)

# KL-divergence post/prior of the parameters of the model
function kldiv_post(model)
    KL = 0.
    for d in 1:length(model)
        for m in model[d]
            KL += kldiv(m.hposterior, m.hprior)
            KL += kldiv(m.λposterior, m.λprior)
        end
    end
    KL
end

# Retrieve the last check-point
fnames = [basename(fname) for fname in glob(joinpath(outdir, "emissions_*bson"))]
fnames = sort(fnames, lt=natural)
epochno = replace(replace(fnames[end], ("emissions_" => "")), (".bson" => ""))
start = parse(Int64, epochno) + 1

if start > 1
    @info "found $(fnames[end]), starting training from epoch $start"
end

@everywhere reducer(a, b) = (mergestats(a[1], b[1]), a[2] + b[2], a[3] + b[3])
for e in start:epochs
    path = joinpath(outdir, "emissions_$(e-1).bson")
    @everywhere emissions = ConvHMM.load($path)

    @everywhere epoch = $e
    stats, totll, totN = @distributed reducer for uttid in uttids
        X = BSON.load(joinpath(feadir, uttid * ".bson"))[:data]
        ali = BSON.load(joinpath(alidir, uttid * ".bson"))[:ali]

        S = length(emissions[1])
        N = size(X, 2)

        # E-step

        # Special case if this is the first epoch and the user has
        # specified a model for initialization
        alimodel = epoch == 1 && ! isnothing(init) ? init : emissions

        lnαβ, totll =   αβrecursion(ali, alimodel(X))
        γ_sparse = resps(ali, lnαβ, dense = false)
        γ = zeros(S, N)
        for idx in keys(γ_sparse) γ[idx, :] = γ_sparse[idx] end

        if e % updaterate > 0
            accstats = accumulate_h(emissions, X, γ)
        else
            accstats = accumulate_λ(emissions, X, γ)
        end

        (accstats, totll, N)
    end

    kl = kldiv_post(emissions)
    𝓛 = (totll - kl) / totN
    @info "epoch $(e)/$(epochs) 𝓛 = $(round(𝓛, digits = 3)) llh = $(round(totll, digits = 3)) KL = $(round(kl, digits=3))"

    if e % updaterate > 0
        update_h!(emissions, stats)
    else
        update_λ!(emissions, stats)
    end

    save(joinpath(outdir, "emissions_$(e).bson"), emissions)
end

max_stars_repo_name: BUTSpeechFIT/ConvHMM

max_issues_repo_name: BUTSpeechFIT/ConvHMM

max_forks_repo_name: BUTSpeechFIT/ConvHMM


text:
function constraint(joint::Joint{T,Nλ,Nb,N,Nb½},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion,
        η, μ) where {T,Nλ,Nb,N,Nb½}

    e1 = joint_constraint(joint, xa, qa, xb, qb, η)
    e2 = minimal_coordinates(joint, xa, qa, xb, qb)

    s, γ = split_impulses(joint, η)

    return [
            s .* γ .- μ;
            s[SUnitRange(1,Nb½)] - (joint.joint_limits[2] .- e2);
            s[SUnitRange(Nb½+1,Nb)] - (e2 .- joint.joint_limits[1]);
            e1;
           ]
end

function constraint_jacobian_configuration(relative::Symbol, joint::Joint{T,Nλ,Nb},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion,
        η) where {T,Nλ,Nb}

    ∇comp = szeros(T,Nb,7)
    ∇mincoord = minimal_coordinates_jacobian_configuration(relative, joint, xa, qa, xb, qb, attjac=false)
    ∇unlim = joint_constraint_jacobian_configuration(relative, joint, xa, qa, xb, qb, η)

    return [∇comp; ∇mincoord; -∇mincoord; ∇unlim]
end

function add_limits(mech::Mechanism, joint::JointConstraint;
    # NOTE: this only works for joints between serial chains (ie, single child joints)
    tra_limits=joint.translational.joint_limits,
    rot_limits=joint.rotational.joint_limits)

    # update translational
    tra = joint.translational
    T = typeof(tra).parameters[1]
    Nλ = typeof(tra).parameters[2]
    Nb½ = length(tra_limits[1])
    Nb = 2Nb½
    N̄λ = 3 - Nλ
    N = Nλ + 2Nb
    tra_limit = (Translational{T,Nλ,Nb,N,Nb½,N̄λ}(tra.axis, tra.axis_mask1, tra.axis_mask2, tra.axis_mask3,
        tra.vertices, tra.spring, tra.damper, tra.spring_offset, tra_limits,
        tra.spring_type, tra.input), joint.parent_id, joint.child_id)

    # update rotational
    rot = joint.rotational
    T = typeof(rot).parameters[1]
    Nλ = typeof(rot).parameters[2]
    Nb½ = length(rot_limits[1])
    Nb = 2Nb½
    N̄λ = 3 - Nλ
    N = Nλ + 2Nb
    rot_limit = (Rotational{T,Nλ,Nb,N,Nb½,N̄λ}(rot.axis, rot.axis_mask1, rot.axis_mask2, rot.axis_mask3,
        rot.orientation_offset, rot.spring, rot.damper, rot.spring_offset, rot_limits,
        rot.spring_type, rot.input), joint.parent_id, joint.child_id)

    JointConstraint((tra_limit, rot_limit); name=joint.name)
end


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

name = "COINMetisBuilder"
version = v"1.3.5"

# Collection of sources required to build MetisBuilder
sources = [
    "https://github.com/coin-or-tools/ThirdParty-Metis/archive/releases/1.3.5.tar.gz" =>
    "98a6110d5d004a16ad42ee26cfac508477f44aa6fe296b90a6413fe0273ebe24",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd ThirdParty-Metis-releases-1.3.5/
./get.Metis
update_configure_scripts
for path in ${LD_LIBRARY_PATH//:/ }; do
    for file in $(ls $path/*.la); do
        echo "$file"
        baddir=$(sed -n "s|libdir=||p" $file)
        sed -i~ -e "s|$baddir|'$path'|g" $file
    done
done
mkdir build
cd build/
## STATIC BUILD START
if [ $target = "x86_64-apple-darwin14" ]; then
  export AR=/opt/x86_64-apple-darwin14/bin/x86_64-apple-darwin14-ar
fi

../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all
## STATIC BUILD END

## DYNAMIC BUILD START
#../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --enable-shared --disable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all
## DYNAMIC BUILD END

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

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

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

]

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


text:
const doublecolon = @eval (:(x::Int)).head

# -- ScannedConsts ------------------------------------------------------------

type ScannedConsts
    names::Vector{Symbol}  # const names in order
    declared::Set{Symbol}  # set of all declared (const) names

    ScannedConsts() = new(Symbol[], Set{Symbol}())
end

function addconst(c::ScannedConsts, name::Symbol)
    if !has(c.declared, name)
        add(c.declared, name)
        push(c.names, name)
    end
end

# -- scanconsts ---------------------------------------------------------------

function scanconsts_let(ex::Expr)
    @assert ex.head == :let
    c = ScannedConsts()
    scanconsts(c, ex.args[1]) # scan the body; header can't have consts
    c.names
end

# # todo: remove?
# function scanconsts(ex)
#     c = ScannedConsts()
#     scanconsts(c, ex)
#     c.names
# end

scanconsts(::ScannedConsts, ::Any) = nothing
function scanconsts(c::ScannedConsts, ex::Expr)
#     println("scanconsts: ex=$ex")

    # don't recurse into scope forming exprs
    if contains([:for, :while, :try, :let, :(->)], ex.head)
        return
    end
    if ex.head == :global; return; end

    # still have to watch out for e g f(x)=x^2
    if ex.head == :(=)
        scanconsts_lhs(c, ex.args[1])
        scanconsts(c, ex.args[2])
    elseif contains([:type, :abstract, :typealias, :const], ex.head)
        scanconsts_const(c, ex.args[1])
    elseif ex.head == :function
        scanconsts_const(c, ex.args[1].args[1]) # assume ex.args[1] is a :call
    elseif ex.head == :return
        error("@namespace: return out of namespace scope not suppported")
    else
        foreach(arg->scanconsts(c, arg), ex.args)
    end
end

scanconsts_lhs(c::ScannedConsts, lhs::Symbol) = nothing
function scanconsts_lhs(c::ScannedConsts, ex::Expr)
#     println("scanconsts_lhs: ex=$ex")

    if ex.head == :call
        scanconsts_const(c, ex.args[1])
        foreach(arg->scanconsts(c, arg), ex.args[2:end])
    else
        foreach(arg->scanconsts(c, arg), ex.args)
    end
end

scanconsts_const(c::ScannedConsts, sym::Symbol) = addconst(c, sym)
function scanconsts_const(c::ScannedConsts, ex::Expr)
#    println("scanconsts_local: ex = $ex")
    if ex.head == :(=)
        scanconsts_const(c, ex.args[1])
        scanconsts(c, ex.args[2])
    elseif ex.head == :comparison
        scanconsts_const(c, ex.args[1])
        scanconsts(c, ex.args[3])
    elseif ex.head == :curly
        scanconsts_const(c, ex.args[1])
        foreach(arg->scanconsts(c, arg), ex.args[2:end])
    elseif ex.head == :global
    elseif ex.head == :local
        foreach(arg->scanconsts_const(c, arg), ex.args)
    else
        error("scanconsts_const: unexpected ex = ", ex)
    end
end

# == @namespace ===============================================================

function declare_struct_type(structname, names, types)
    fields = {expr(doublecolon, name, t) for (name, t) in zip(names, types)}
    tdef = expr(:type, structname, expr(:block, fields))
    eval(tdef)
end

macro namespace(name::Symbol, body::Expr)
    namespace(name, body)
end

function namespace(typename::Symbol, body::Expr)
    if (body.head != :let)
        #error("@namespace: body must be a let block")
    end
    fieldnames = scanconsts_let(body)

    @gensym NamespaceStruct types x
    epilogue = quote

        $types = {typeof($x) for ($x) in {$fieldnames...}}
#        println($types)
        declare_struct_type(($expr(:quote, NamespaceStruct)),
                            {${expr(:quote, name) for name in fieldnames}...},
                            $types)
        ($NamespaceStruct)($fieldnames...)
    end

    augblock = expr(:block, body.args[1], epilogue)
    auglet   = expr(:let, augblock, body.args[2:end]...)

    quote
        $typename = $auglet
        nothing
    end
end

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
# Model not supporting Interval
MOIU.@model(SimpleModel,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone,
             MOI.RotatedSecondOrderCone, MOI.GeometricMeanCone,
             MOI.PositiveSemidefiniteConeTriangle, MOI.ExponentialCone),
            (MOI.PowerCone, MOI.DualPowerCone),
            (),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction))


text:
import UUIDs: UUID, uuid1
import .ExpressionExplorer: SymbolsState

"The building block of a `Notebook`. Contains code, output, reactivity data, mitochondria and ribosomes."
Base.@kwdef mutable struct Cell
    "Because Cells can be reordered, they get a UUID. The JavaScript frontend indexes cells using the UUID."
    cell_id::UUID=uuid1()

    code::String=""
    # code_last_updated::Float64=0
    # code_last_updated_by::Union{String,Nothing}=nothing
    code_folded::Bool=false

    output_repr::Union{Nothing,String,Vector{UInt8},Dict}=nothing
    repr_mime::MIME=MIME("text/plain")
    errored::Bool=false
    runtime::Union{Missing,UInt64}=missing
    queued::Bool=false
    running::Bool=false

    "Time that the last output was created, used only on the frontend to rerender the output"
    last_run_timestamp::Float64=0
    "Whether `this` inside `<script id=something>` should refer to the previously returned object in HTML output. This is used for fancy animations. true iff a cell runs as a reactive consequence."
    persist_js_state::Bool=false

    parsedcode::Union{Nothing,Expr}=nothing
    module_usings::Set{Expr}=Set{Expr}()
    rootassignee::Union{Nothing,Symbol}=nothing
    function_wrapped::Bool=false
end

Cell(cell_id, code) = Cell(cell_id=cell_id, code=code)
Cell(code) = Cell(uuid1(), code)

function Base.convert(::Type{Cell}, cell::Dict)
	Cell(
        cell_id=UUID(cell["cell_id"]),
        code=cell["code"],
        code_folded=cell["code_folded"],
    )
end
function Base.convert(::Type{UUID}, string::String)
    UUID(string)
end


text:
using DrWatson
@quickactivate "EPExperiment"

using EPExperiments
using AnnealedIS

function multiple_expectations()
    @expectation function expt_prog(y)
        x ~ Normal(0, 1)
        y ~ Normal(x, 1)
        return x, x^2, x^3
    end

    y_observed = 3
    expt_prog1, expr_prog2, expt_prog3 = expt_prog
    expct1 = expt_prog1(y_observed)

    expct1_estimate, diagnostics = estimate_expectation(
        expct1, TABI(TuringAlgorithm(AnIS(), 100))
    )
end

multiple_expectations()


text:
"Represents the order of "
struct Order
end

"reorder `ir` according to `order`"
function reorder(ir::IR, order::Order)
end


text:
# Algorithms for Maximum-a-posteriori inference

module MAP

using SPNetworks
import SPNetworks:
    Node, SumNode, ProductNode, LeafNode, CategoricalDistribution, IndicatorFunction, GaussianDistribution,
    isleaf, isprod, issum,
    logpdf!,
    vardims

include("MaxProduct.jl")
include("LocalSearch.jl")
include("BeliefPropagation.jl")

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Implementation of
#  "Table-driven Implementation of the Logarithm Function in IEEE Floating-point Arithmetic"
#  Tang, Ping-Tak Peter
#  ACM Trans. Math. Softw. (1990), 16(4):378--400
#  https://doi.org/10.1145/98267.98294

# Does not currently handle floating point flags (inexact, div-by-zero, etc).

import .Base.unsafe_trunc
import .Base.Math.@horner

# Float64 lookup table.
# to generate values:
  # N=39 # (can be up to N=42).
  # sN = 2.0^N
  # isN = 1.0/sN
  # s7 = 2.0^7
  # is7 = 1.0/s7
  # for j=0:128
  #   l_big = Base.log(big(1.0+j*is7))
  #   l_hi = isN*Float64(round(sN*l_big))
  #   l_lo = Float64(l_big-l_hi)
  #   j % 2 == 0 && print("\n    ")
  #   print("(",l_hi,",",l_lo,"),")
  # end

const t_log_Float64 = ((0.0,0.0),(0.007782140442941454,-8.865052917267247e-13),
    (0.015504186536418274,-4.530198941364935e-13),(0.0231670592820592,-5.248209479295644e-13),
    (0.03077165866670839,4.529814257790929e-14),(0.0383188643027097,-5.730994833076631e-13),
    (0.04580953603181115,-5.16945692881222e-13),(0.053244514518155484,6.567993368985218e-13),
    (0.06062462181580486,6.299848199383311e-13),(0.06795066190898069,-4.729424109166329e-13),
    (0.07522342123775161,-1.6408301585598662e-13),(0.08244366921098845,8.614512936087814e-14),
    (0.08961215869021544,-5.283050530808144e-13),(0.09672962645890948,-3.5836667430094137e-13),
    (0.10379679368088546,7.581073923016376e-13),(0.11081436634049169,-2.0157368416016215e-13),
    (0.11778303565552051,8.629474042969438e-13),(0.1247034785010328,-7.556920687451337e-14),
    (0.1315763577895268,-8.075373495358435e-13),(0.13840232285838283,7.363043577087051e-13),
    (0.14518200984457508,-7.718001336828099e-14),(0.15191604202664166,-7.996871607743758e-13),
    (0.15860503017574956,8.890223439724663e-13),(0.16524957289584563,-5.384682618788232e-13),
    (0.17185025692742784,-7.686134224018169e-13),(0.17840765747314435,-3.2605717931058157e-13),
    (0.18492233849428885,-2.7685884310448306e-13),(0.1913948530000198,-3.903387893794952e-13),
    (0.1978257433293038,6.160755775588723e-13),(0.20421554142922105,-5.30156516006026e-13),
    (0.21056476910780475,-4.55112422774782e-13),(0.21687393830143264,-8.182853292737783e-13),
    (0.22314355131493357,-7.238189921749681e-13),(0.22937410106533207,-4.86240001538379e-13),
    (0.23556607131286,-9.30945949519689e-14),(0.24171993688651128,6.338907368997553e-13),
    (0.24783616390413954,4.4171755371315547e-13),(0.25391520998164196,-6.785208495970588e-13),
    (0.25995752443668607,2.3999540484211735e-13),(0.2659635484978935,-7.555569400283742e-13),
    (0.27193371548310097,5.407904186145515e-13),(0.2778684510030871,3.692037508208009e-13),
    (0.28376817313073843,-9.3834172236637e-14),(0.28963329258294834,9.43339818951269e-14),
    (0.29546421289342106,4.148131870425857e-13),(0.3012613305781997,-3.7923164802093147e-14),
    (0.3070250352957373,-8.25463138725004e-13),(0.31275571000333,5.668653582900739e-13),
    (0.318453731119007,-4.723727821986367e-13),(0.32411946865431673,-1.0475750058776541e-13),
    (0.32975328637257917,-1.1118671389559323e-13),(0.33535554192167183,-5.339989292003297e-13),
    (0.3409265869704541,1.3912841212197566e-13),(0.3464667673470103,-8.017372713972018e-13),
    (0.35197642315688427,2.9391859187648e-13),(0.3574558889213222,4.815896111723205e-13),
    (0.3629054936900502,-6.817539406325327e-13),(0.36832556115950865,-8.009990055432491e-13),
    (0.3737164097929053,6.787566823158706e-13),(0.37907835293481185,1.5761203773969435e-13),
    (0.3844116989112081,-8.760375990774874e-13),(0.38971675114044046,-4.152515806343612e-13),
    (0.3949938082405424,3.2655698896907146e-13),(0.40024316412745975,-4.4704265010452445e-13),
    (0.4054651081078191,3.452764795203977e-13),(0.4106599249844294,8.390050778518307e-13),
    (0.4158278951435932,1.1776978751369214e-13),(0.4209692946442374,-1.0774341461609579e-13),
    (0.42608439531068143,2.186334329321591e-13),(0.43117346481813,2.413263949133313e-13),
    (0.4362367667745275,3.90574622098307e-13),(0.44127456080423144,6.437879097373207e-13),
    (0.44628710262804816,3.713514191959202e-13),(0.45127464413963025,-1.7166921336082432e-13),
    (0.4562374334818742,-2.8658285157914353e-13),(0.4611757151214988,6.713692791384601e-13),
    (0.46608972992544295,-8.437281040871276e-13),(0.4709797152190731,-2.821014384618127e-13),
    (0.4758459048698569,1.0701931762114255e-13),(0.4806885293455707,1.8119346366441111e-13),
    (0.4855078157816024,9.840465278232627e-14),(0.49030398804461583,5.780031989454028e-13),
    (0.49507726679803454,-1.8302857356041668e-13),(0.4998278695566114,-1.620740015674495e-13),
    (0.5045560107519123,4.83033149495532e-13),(0.5092619017905236,-7.156055317238212e-13),
    (0.5139457511013461,8.882123951857185e-13),(0.5186077642083546,-3.0900580513238243e-13),
    (0.5232481437651586,-6.10765519728515e-13),(0.5278670896204858,3.565996966334783e-13),
    (0.532464798869114,3.5782396591276384e-13),(0.5370414658973459,-4.622608700154458e-13),
    (0.5415972824321216,6.227976291722515e-13),(0.5461324375974073,7.283894727206574e-13),
    (0.5506471179523942,2.680964661521167e-13),(0.5551415075406112,-1.0960825046059278e-13),
    (0.5596157879353996,2.3119493838005378e-14),(0.5640701382853877,-5.846905800529924e-13),
    (0.5685047353526897,-2.1037482511444942e-14),(0.5729197535620187,-2.332318294558741e-13),
    (0.5773153650352469,-4.2333694288141915e-13),(0.5816917396350618,-4.3933937969737843e-13),
    (0.5860490450031648,4.1341647073835564e-13),(0.590387446602108,6.841763641591467e-14),
    (0.5947071077462169,4.758553400443064e-13),(0.5990081896452466,8.367967867475769e-13),
    (0.6032908514389419,-8.576373464665864e-13),(0.6075552502243227,2.1913281229340092e-13),
    (0.6118015411066153,-6.224284253643115e-13),(0.6160298772156239,-1.098359432543843e-13),
    (0.6202404097512044,6.531043137763365e-13),(0.6244332880123693,-4.758019902171077e-13),
    (0.6286086594227527,-3.785425126545704e-13),(0.6327666695706284,4.0939233218678666e-13),
    (0.636907462236195,8.742438391485829e-13),(0.6410311794206791,2.521818845684288e-13),
    (0.6451379613736208,-3.6081313604225574e-14),(0.649227946625615,-5.05185559242809e-13),
    (0.6533012720119586,7.869940332335532e-13),(0.6573580727090302,-6.702087696194906e-13),
    (0.6613984822452039,1.6108575753932459e-13),(0.6654226325445052,5.852718843625151e-13),
    (0.6694306539429817,-3.5246757297904794e-13),(0.6734226752123504,-1.8372084495629058e-13),
    (0.6773988235909201,8.860668981349492e-13),(0.6813592248072382,6.64862680714687e-13),
    (0.6853040030982811,6.383161517064652e-13),(0.6892332812385575,2.5144230728376075e-13),
    (0.6931471805601177,-1.7239444525614835e-13))

# Float32 lookup table
# to generate values:
  # N=16
  # sN = 2f0^N
  # isN = 1f0/sN
  # s7 = 2.0^7
  # is7 = 1.0/s7
  # for j=0:128
  #   j % 4 == 0 && print("\n    ")
  #   print(float64(Base.log(big(1.0+j*is7))),",")
  # end

const t_log_Float32 = (0.0,0.007782140442054949,0.015504186535965254,0.02316705928153438,
    0.030771658666753687,0.0383188643021366,0.0458095360312942,0.053244514518812285,
    0.06062462181643484,0.06795066190850775,0.07522342123758753,0.08244366921107459,
    0.08961215868968714,0.09672962645855111,0.10379679368164356,0.11081436634029011,
    0.11778303565638346,0.12470347850095724,0.13157635778871926,0.13840232285911913,
    0.1451820098444979,0.15191604202584197,0.15860503017663857,0.16524957289530717,
    0.17185025692665923,0.1784076574728183,0.184922338494012,0.19139485299962947,
    0.19782574332991987,0.2042155414286909,0.21056476910734964,0.21687393830061436,
    0.22314355131420976,0.22937410106484582,0.2355660713127669,0.24171993688714516,
    0.24783616390458127,0.25391520998096345,0.25995752443692605,0.26596354849713794,
    0.27193371548364176,0.2778684510034563,0.2837681731306446,0.28963329258304266,
    0.2954642128938359,0.3012613305781618,0.3070250352949119,0.3127557100038969,
    0.3184537311185346,0.324119468654212,0.329753286372468,0.3353555419211378,
    0.3409265869705932,0.34646676734620857,0.3519764231571782,0.3574558889218038,
    0.3629054936893685,0.3683255611587076,0.37371640979358406,0.37907835293496944,
    0.38441169891033206,0.3897167511400252,0.394993808240869,0.4002431641270127,
    0.4054651081081644,0.4106599249852684,0.415827895143711,0.42096929464412963,
    0.4260843953109001,0.4311734648183713,0.43623676677491807,0.4412745608048752,
    0.44628710262841953,0.45127464413945856,0.4562374334815876,0.46117571512217015,
    0.46608972992459924,0.470979715218791,0.4758459048699639,0.4806885293457519,
    0.4855078157817008,0.4903039880451938,0.4950772667978515,0.4998278695564493,
    0.5045560107523953,0.5092619017898079,0.5139457511022343,0.5186077642080457,
    0.5232481437645479,0.5278670896208424,0.5324647988694718,0.5370414658968836,
    0.5415972824327444,0.5461324375981357,0.5506471179526623,0.5551415075405016,
    0.5596157879354227,0.564070138284803,0.5685047353526688,0.5729197535617855,
    0.5773153650348236,0.5816917396346225,0.5860490450035782,0.5903874466021763,
    0.5947071077466928,0.5990081896460834,0.6032908514380843,0.6075552502245418,
    0.6118015411059929,0.616029877215514,0.6202404097518576,0.6244332880118935,
    0.6286086594223741,0.6327666695710378,0.6369074622370692,0.6410311794209312,
    0.6451379613735847,0.6492279466251099,0.6533012720127457,0.65735807270836,
    0.661398482245365,0.6654226325450905,0.6694306539426292,0.6734226752121667,
    0.6773988235918061,0.6813592248079031,0.6853040030989194,0.689233281238809,
    0.6931471805599453)

# truncate lower order bits (up to 26)
# ideally, this should be able to use ANDPD instructions, see #9868.
@inline function truncbits(x::Float64)
    reinterpret(Float64, reinterpret(UInt64,x) & 0xffff_ffff_f800_0000)
end

logb(::Type{Float32},::Val{2})  = 1.4426950408889634
logb(::Type{Float32},::Val{:ℯ}) = 1.0
logb(::Type{Float32},::Val{10}) = 0.4342944819032518
logbU(::Type{Float64},::Val{2})  = 1.4426950408889634
logbL(::Type{Float64},::Val{2})  = 2.0355273740931033e-17
logbU(::Type{Float64},::Val{:ℯ}) = 1.0
logbL(::Type{Float64},::Val{:ℯ}) = 0.0
logbU(::Type{Float64},::Val{10}) = 0.4342944819032518
logbL(::Type{Float64},::Val{10}) = 1.098319650216765e-17

# Procedure 1
# XXX we want to mark :consistent-cy here so that this function can be concrete-folded,
# because the effect analysis currently can't prove it in the presence of `@inbounds` or
# `:boundscheck`, but still the access to `t_log_Float64` is really safe here
Base.@assume_effects :consistent @inline function log_proc1(y::Float64,mf::Float64,F::Float64,f::Float64,base=Val(:ℯ))
    jp = unsafe_trunc(Int,128.0*F)-127

    ## Steps 1 and 2
    @inbounds hi,lo = t_log_Float64[jp]
    l_hi = mf* 0.6931471805601177 + hi
    l_lo = mf*-1.7239444525614835e-13 + lo

    ## Step 3
    # @inbounds u = f*c_invF[jp]
    # u = f/F
    # q = u*u*@horner(u,
    #                 -0x1.0_0000_0000_0001p-1,
    #                 +0x1.5_5555_5550_9ba5p-2,
    #                 -0x1.f_ffff_ffeb_6526p-3,
    #                 +0x1.9_99b4_dfed_6fe4p-3,
    #                 -0x1.5_5576_6647_2e04p-3)

    ## Step 3' (alternative)
    u = (2.0f)/(y+F)
    v = u*u
    q = u*v*@horner(v,
                    0.08333333333303913,
                    0.012500053168098584)

    ## Step 4
    m_hi = logbU(Float64, base)
    m_lo = logbL(Float64, base)
    return fma(m_hi, l_hi, fma(m_hi, (u + (q + l_lo)), m_lo*l_hi))
end

# Procedure 2
@inline function log_proc2(f::Float64,base=Val(:ℯ))
    ## Step 1
    g = 1.0/(2.0+f)
    u = 2.0*f*g
    v = u*u

    ## Step 2
    q = u*v*@horner(v,
                    0.08333333333333179,
                    0.012500000003771751,
                    0.0022321399879194482,
                    0.0004348877777076146)

    ## Step 3
    # based on:
    #   2(f-u) = 2(f(2+f)-2f)/(2+f) = 2f^2/(2+f) = fu
    #   2(f-u1-u2) - f*(u1+u2) = 0
    #   2(f-u1) - f*u1 = (2+f)u2
    #   u2 = (2(f-u1) - f*u1)/(2+f)

    m_hi = logbU(Float64, base)
    m_lo = logbL(Float64, base)
    return fma(m_hi, u, fma(m_lo, u, m_hi*fma(fma(-u,f,2(f-u)), g, q)))
end

# Procedure 1
# XXX we want to mark :consistent-cy here so that this function can be concrete-folded,
# because the effect analysis currently can't prove it in the presence of `@inbounds` or
# `:boundscheck`, but still the access to `t_log_Float32` is really safe here
Base.@assume_effects :consistent @inline function log_proc1(y::Float32,mf::Float32,F::Float32,f::Float32,base=Val(:ℯ))
    jp = unsafe_trunc(Int,128.0f0*F)-127

    ## Steps 1 and 2
    @inbounds hi = t_log_Float32[jp]
    l = mf*0.6931471805599453 + hi

    ## Step 3
    # @inbounds u = f*c_invF[jp]
    # q = u*u*@horner(u,
    #                 Float32(-0x1.00006p-1),
    #                 Float32(0x1.55546cp-2))

    ## Step 3' (alternative)
    u = (2f0f)/(y+F)
    v = u*u
    q = u*v*0.08333351f0

    ## Step 4
    Float32(logb(Float32, base)*(l + (u + q)))
end

# Procedure 2
@inline function log_proc2(f::Float32,base=Val(:ℯ))
    ## Step 1
    # compute in higher precision
    u64 = Float64(2f0*f)/(2.0+f)
    u = Float32(u64)
    v = u*u

    ## Step 2
    q = u*v*@horner(v,
                    0.08333332f0,
                    0.012512346f0)

    ## Step 3: not required

    ## Step 4
    Float32(logb(Float32, base)*(u64 + q))
end

log2(x::Float32)  = _log(x, Val(2),  :log2)
log(x::Float32)   = _log(x, Val(:ℯ), :log)
log10(x::Float32) = _log(x, Val(10), :log10)
log2(x::Float64)  = _log(x, Val(2),  :log2)
log(x::Float64)   = _log(x, Val(:ℯ), :log)
log10(x::Float64) = _log(x, Val(10), :log10)

function _log(x::Float64, base, func)
    if x > 0.0
        x == Inf && return x

        # Step 2
        if 0.9394130628134757 < x < 1.0644944589178595
            f = x-1.0
            return log_proc2(f, base)
        end

        # Step 3
        xu = reinterpret(UInt64,x)
        m = Int(xu >> 52) & 0x07ff
        if m == 0 # x is subnormal
            x *= 1.8014398509481984e16 # 0x1p54, normalise significand
            xu = reinterpret(UInt64,x)
            m = Int(xu >> 52) & 0x07ff - 54
        end
        m -= 1023
        y = reinterpret(Float64,(xu & 0x000f_ffff_ffff_ffff) | 0x3ff0_0000_0000_0000)

        mf = Float64(m)
        F = (y + 3.5184372088832e13) - 3.5184372088832e13 # 0x1p-7*round(0x1p7*y)
        f = y-F

        return log_proc1(y,mf,F,f,base)
    elseif x == 0.0
        -Inf
    elseif isnan(x)
        NaN
    else
        throw_complex_domainerror(func, x)
    end
end

function _log(x::Float32, base, func)
    if x > 0f0
        x == Inf32 && return x

        # Step 2
        if 0.939413f0 < x < 1.0644945f0
            f = x-1f0
            return log_proc2(f, base)
        end

        # Step 3
        xu = reinterpret(UInt32,x)
        m = Int(xu >> 23) & 0x00ff
        if m == 0 # x is subnormal
            x *= 3.3554432f7 # 0x1p25, normalise significand
            xu = reinterpret(UInt32,x)
            m = Int(xu >> 23) & 0x00ff - 25
        end
        m -= 127
        y = reinterpret(Float32,(xu & 0x007f_ffff) | 0x3f80_0000)

        mf = Float32(m)
        F = (y + 65536.0f0) - 65536.0f0 # 0x1p-7*round(0x1p7*y)
        f = y-F

        log_proc1(y,mf,F,f,base)
    elseif x == 0f0
        -Inf32
    elseif isnan(x)
        NaN32
    else
        throw_complex_domainerror(func, x)
    end
end

function log1p(x::Float64)
    if x > -1.0
        x == Inf && return x
        if -1.1102230246251565e-16 < x < 1.1102230246251565e-16
            return x # Inexact

        # Step 2
        elseif -0.06058693718652422 < x < 0.06449445891785943
            return log_proc2(x)
        end

        # Step 3
        z = 1.0 + x
        zu = reinterpret(UInt64,z)
        s = reinterpret(Float64,0x7fe0_0000_0000_0000 - (zu & 0xfff0_0000_0000_0000)) # 2^-m
        m = Int(zu >> 52) & 0x07ff - 1023 # z cannot be subnormal
        c = m > 0 ? 1.0-(z-x) : x-(z-1.0) # 1+x = z+c exactly
        y = reinterpret(Float64,(zu & 0x000f_ffff_ffff_ffff) | 0x3ff0_0000_0000_0000)

        mf = Float64(m)
        F = (y + 3.5184372088832e13) - 3.5184372088832e13 # 0x1p-7*round(0x1p7*y)
        f = (y - F) + c*s #2^m(F+f) = 1+x = z+c

        log_proc1(y,mf,F,f)
    elseif x == -1.0
        -Inf
    elseif isnan(x)
        NaN
    else
        throw_complex_domainerror(:log1p, x)
    end
end

function log1p(x::Float32)
    if x > -1f0
        x == Inf32 && return x
        if -5.9604645f-8 < x < 5.9604645f-8
            return x # Inexact
        # Step 2
        elseif -0.06058694f0 < x < 0.06449446f0
            return log_proc2(x)
        end

        # Step 3
        z = 1f0 + x
        zu = reinterpret(UInt32,z)
        s = reinterpret(Float32,0x7f000000 - (zu & 0xff80_0000)) # 2^-m
        m = Int(zu >> 23) & 0x00ff - 127 # z cannot be subnormal
        c = m > 0 ? 1f0-(z-x) : x-(z-1f0) # 1+x = z+c
        y = reinterpret(Float32,(zu & 0x007f_ffff) | 0x3f80_0000)

        mf = Float32(m)
        F = (y + 65536.0f0) - 65536.0f0 # 0x1p-7*round(0x1p7*y)
        f = (y - F) + s*c #2^m(F+f) = 1+x = z+c

        log_proc1(y,mf,F,f)
    elseif x == -1f0
        -Inf32
    elseif isnan(x)
        NaN32
    else
        throw_complex_domainerror(:log1p, x)
    end
end

#function make_compact_table(N)
#    table = Tuple{UInt64,Float64}[]
#    lo, hi = 0x1.69555p-1, 0x1.69555p0
#    for i in 0:N-1
#        # I am not fully sure why this is the right formula to use, but it apparently is
#        center = i/(2*N) + lo < 1 ? (i+.5)/(2*N) + lo : (i+.5)/N + hi -1
#        invc = Float64(center < 1 ? round(N/center)/N : round(2*N/center)/(N*2))
#        c = inv(big(invc))
#        logc = Float64(round(0x1p43*log(c))/0x1p43)
#        logctail = reinterpret(Float64, Float64(log(c) - logc))
#        p1 = (reinterpret(UInt64,invc) >> 45) % UInt8
#        push!(table, (p1|reinterpret(UInt64,logc),logctail))
#    end
#    return Tuple(table)
#end
#const t_log_table_compat = make_compact_table(128)
const t_log_table_compat = (
    (0xbfd62c82f2b9c8b5, 5.929407345889625e-15),
    (0xbfd5d1bdbf5808b4, -2.544157440035963e-14),
    (0xbfd57677174558b3, -3.443525940775045e-14),
    (0xbfd51aad872df8b2, -2.500123826022799e-15),
    (0xbfd4be5f957778b1, -8.929337133850617e-15),
    (0xbfd4618bc21c60b0, 1.7625431312172662e-14),
    (0xbfd404308686a8af, 1.5688303180062087e-15),
    (0xbfd3a64c556948ae, 2.9655274673691784e-14),
    (0xbfd347dd9a9880ad, 3.7923164802093147e-14),
    (0xbfd2e8e2bae120ac, 3.993416384387844e-14),
    (0xbfd2895a13de88ab, 1.9352855826489123e-14),
    (0xbfd2895a13de88ab, 1.9352855826489123e-14),
    (0xbfd22941fbcf78aa, -1.9852665484979036e-14),
    (0xbfd1c898c16998a9, -2.814323765595281e-14),
    (0xbfd1675cababa8a8, 2.7643769993528702e-14),
    (0xbfd1058bf9ae48a7, -4.025092402293806e-14),
    (0xbfd0a324e27390a6, -1.2621729398885316e-14),
    (0xbfd0402594b4d0a5, -3.600176732637335e-15),
    (0xbfd0402594b4d0a5, -3.600176732637335e-15),
    (0xbfcfb9186d5e40a4, 1.3029797173308663e-14),
    (0xbfcef0adcbdc60a3, 4.8230289429940886e-14),
    (0xbfce27076e2af0a2, -2.0592242769647135e-14),
    (0xbfcd5c216b4fc0a1, 3.149265065191484e-14),
    (0xbfcc8ff7c79aa0a0, 4.169796584527195e-14),
    (0xbfcc8ff7c79aa0a0, 4.169796584527195e-14),
    (0xbfcbc286742d909f, 2.2477465222466186e-14),
    (0xbfcaf3c94e80c09e, 3.6507188831790577e-16),
    (0xbfca23bc1fe2b09d, -3.827767260205414e-14),
    (0xbfca23bc1fe2b09d, -3.827767260205414e-14),
    (0xbfc9525a9cf4509c, -4.7641388950792196e-14),
    (0xbfc87fa06520d09b, 4.9278276214647115e-14),
    (0xbfc7ab890210e09a, 4.9485167661250996e-14),
    (0xbfc7ab890210e09a, 4.9485167661250996e-14),
    (0xbfc6d60fe719d099, -1.5003333854266542e-14),
    (0xbfc5ff3070a79098, -2.7194441649495324e-14),
    (0xbfc5ff3070a79098, -2.7194441649495324e-14),
    (0xbfc526e5e3a1b097, -2.99659267292569e-14),
    (0xbfc44d2b6ccb8096, 2.0472357800461955e-14),
    (0xbfc44d2b6ccb8096, 2.0472357800461955e-14),
    (0xbfc371fc201e9095, 3.879296723063646e-15),
    (0xbfc29552f81ff094, -3.6506824353335045e-14),
    (0xbfc1b72ad52f6093, -5.4183331379008994e-14),
    (0xbfc1b72ad52f6093, -5.4183331379008994e-14),
    (0xbfc0d77e7cd09092, 1.1729485484531301e-14),
    (0xbfc0d77e7cd09092, 1.1729485484531301e-14),
    (0xbfbfec9131dbe091, -3.811763084710266e-14),
    (0xbfbe27076e2b0090, 4.654729747598445e-14),
    (0xbfbe27076e2b0090, 4.654729747598445e-14),
    (0xbfbc5e548f5bc08f, -2.5799991283069902e-14),
    (0xbfba926d3a4ae08e, 3.7700471749674615e-14),
    (0xbfba926d3a4ae08e, 3.7700471749674615e-14),
    (0xbfb8c345d631a08d, 1.7306161136093256e-14),
    (0xbfb8c345d631a08d, 1.7306161136093256e-14),
    (0xbfb6f0d28ae5608c, -4.012913552726574e-14),
    (0xbfb51b073f06208b, 2.7541708360737882e-14),
    (0xbfb51b073f06208b, 2.7541708360737882e-14),
    (0xbfb341d7961be08a, 5.0396178134370583e-14),
    (0xbfb341d7961be08a, 5.0396178134370583e-14),
    (0xbfb16536eea38089, 1.8195060030168815e-14),
    (0xbfaf0a30c0118088, 5.213620639136504e-14),
    (0xbfaf0a30c0118088, 5.213620639136504e-14),
    (0xbfab42dd71198087, 2.532168943117445e-14),
    (0xbfab42dd71198087, 2.532168943117445e-14),
    (0xbfa77458f632c086, -5.148849572685811e-14),
    (0xbfa77458f632c086, -5.148849572685811e-14),
    (0xbfa39e87b9fec085, 4.6652946995830086e-15),
    (0xbfa39e87b9fec085, 4.6652946995830086e-15),
    (0xbf9f829b0e780084, -4.529814257790929e-14),
    (0xbf9f829b0e780084, -4.529814257790929e-14),
    (0xbf97b91b07d58083, -4.361324067851568e-14),
    (0xbf8fc0a8b0fc0082, -1.7274567499706107e-15),
    (0xbf8fc0a8b0fc0082, -1.7274567499706107e-15),
    (0xbf7fe02a6b100081, -2.298941004620351e-14),
    (0xbf7fe02a6b100081, -2.298941004620351e-14),
    (0x0000000000000080, 0.0),
    (0x0000000000000080, 0.0),
    (0x3f8010157589007e, -1.4902732911301337e-14),
    (0x3f9020565893807c, -3.527980389655325e-14),
    (0x3f98492528c9007a, -4.730054772033249e-14),
    (0x3fa0415d89e74078, 7.580310369375161e-15),
    (0x3fa466aed42e0076, -4.9893776716773285e-14),
    (0x3fa894aa149fc074, -2.262629393030674e-14),
    (0x3faccb73cdddc072, -2.345674491018699e-14),
    (0x3faeea31c006c071, -1.3352588834854848e-14),
    (0x3fb1973bd146606f, -3.765296820388875e-14),
    (0x3fb3bdf5a7d1e06d, 5.1128335719851986e-14),
    (0x3fb5e95a4d97a06b, -5.046674438470119e-14),
    (0x3fb700d30aeac06a, 3.1218748807418837e-15),
    (0x3fb9335e5d594068, 3.3871241029241416e-14),
    (0x3fbb6ac88dad6066, -1.7376727386423858e-14),
    (0x3fbc885801bc4065, 3.957125899799804e-14),
    (0x3fbec739830a2063, -5.2849453521890294e-14),
    (0x3fbfe89139dbe062, -3.767012502308738e-14),
    (0x3fc1178e8227e060, 3.1859736349078334e-14),
    (0x3fc1aa2b7e23f05f, 5.0900642926060466e-14),
    (0x3fc2d1610c86805d, 8.710783796122478e-15),
    (0x3fc365fcb015905c, 6.157896229122976e-16),
    (0x3fc4913d8333b05a, 3.821577743916796e-14),
    (0x3fc527e5e4a1b059, 3.9440046718453496e-14),
    (0x3fc6574ebe8c1057, 2.2924522154618074e-14),
    (0x3fc6f0128b757056, -3.742530094732263e-14),
    (0x3fc7898d85445055, -2.5223102140407338e-14),
    (0x3fc8beafeb390053, -1.0320443688698849e-14),
    (0x3fc95a5adcf70052, 1.0634128304268335e-14),
    (0x3fca93ed3c8ae050, -4.3425422595242564e-14),
    (0x3fcb31d8575bd04f, -1.2527395755711364e-14),
    (0x3fcbd087383be04e, -5.204008743405884e-14),
    (0x3fcc6ffbc6f0104d, -3.979844515951702e-15),
    (0x3fcdb13db0d4904b, -4.7955860343296286e-14),
    (0x3fce530effe7104a, 5.015686013791602e-16),
    (0x3fcef5ade4dd0049, -7.252318953240293e-16),
    (0x3fcf991c6cb3b048, 2.4688324156011588e-14),
    (0x3fd07138604d5846, 5.465121253624792e-15),
    (0x3fd0c42d67616045, 4.102651071698446e-14),
    (0x3fd1178e8227e844, -4.996736502345936e-14),
    (0x3fd16b5ccbacf843, 4.903580708156347e-14),
    (0x3fd1bf99635a6842, 5.089628039500759e-14),
    (0x3fd214456d0eb841, 1.1782016386565151e-14),
    (0x3fd2bef07cdc903f, 4.727452940514406e-14),
    (0x3fd314f1e1d3603e, -4.4204083338755686e-14),
    (0x3fd36b6776be103d, 1.548345993498083e-14),
    (0x3fd3c2527733303c, 2.1522127491642888e-14),
    (0x3fd419b423d5e83b, 1.1054030169005386e-14),
    (0x3fd4718dc271c83a, -5.534326352070679e-14),
    (0x3fd4c9e09e173039, -5.351646604259541e-14),
    (0x3fd522ae0738a038, 5.4612144489920215e-14),
    (0x3fd57bf753c8d037, 2.8136969901227338e-14),
    (0x3fd5d5bddf596036, -1.156568624616423e-14))

 @inline function log_tab_unpack(t::UInt64)
    invc = UInt64(t&UInt64(0xff)|0x1ff00)<<45
    logc = t&(~UInt64(0xff))
    return (reinterpret(Float64, invc), reinterpret(Float64, logc))
end

# Log implementation that returns 2 numbers which sum to give true value with about 68 bits of precision
# Since `log` only makes sense for positive exponents, we speed up the implimentation by stealing the sign bit
# of the input for an extra bit of the exponent which is used to normalize subnormal inputs.
# Does not normalize results.
# Adapted and modified from https://github.com/ARM-software/optimized-routines/blob/master/math/pow.c
# Copyright (c) 2018-2020, Arm Limited. (which is also MIT licensed)
# note that this isn't an exact translation as this version compacts the table to reduce cache pressure.
function _log_ext(xu)
    # x = 2^k z; where z is in range [0x1.69555p-1,0x1.69555p-0) and exact.
    # The range is split into N subintervals.
    # The ith subinterval contains z and c is near the center of the interval.
    tmp = reinterpret(Int64, xu - 0x3fe6955500000000) #0x1.69555p-1
    i = (tmp >> 45) & 127
    z = reinterpret(Float64, xu - (tmp & 0xfff0000000000000))
    k = Float64(tmp >> 52)
    # log(x) = k*Ln2 + log(c) + log1p(z/c-1).
    t, logctail = t_log_table_compat[i+1]
    invc, logc = log_tab_unpack(t)
    # Note: invc is j/N or j/N/2 where j is an integer in [N,2N) and
    # |z/c - 1| < 1/N, so r = z/c - 1 is exactly representible.
    r = fma(z, invc, -1.0)
    # k*Ln2 + log(c) + r.
    t1 = muladd(k, 0.6931471805598903, logc) #ln(2) hi part
    t2 = t1 + r
    lo1 = muladd(k, 5.497923018708371e-14, logctail) #ln(2) lo part
    lo2 = t1 - t2 + r
    ar = -0.5 * r
    ar2, lo3 = two_mul(r, ar)
    # k*Ln2 + log(c) + r + .5*r*r.
    hi = t2 + ar2
    lo4 = t2 - hi + ar2
    p = evalpoly(r, (-0x1.555555555556p-1, 0x1.0000000000006p-1, -0x1.999999959554ep-2, 0x1.555555529a47ap-2, -0x1.2495b9b4845e9p-2, 0x1.0002b8b263fc3p-2))
    lo = lo1 + lo2 + lo3 + muladd(r*ar2, p, lo4)
    return hi, lo
end

max_stars_repo_name: TechPenguineer/julia

max_issues_repo_name: TechPenguineer/julia

max_forks_repo_name: Nexuscompute/julia