text:
module ParetoSmooth
    using Requires
    using DocStringExtensions

    function __init__()
        @require Turing =
            "fce5fe82-541a-59a6-adf8-730c64b5f9a0" include("TuringHelpers.jl")
        @require MCMCChains =
            "c7f686f2-ff18-58e9-bc7b-31028e88f75d" include("MCMCChainsHelpers.jl")
    end

    include("ESS.jl")
    include("GPD.jl")
    include("InternalHelpers.jl")
    include("ImportanceSampling.jl")
    include("LooStructs.jl")
    include("LeaveOneOut.jl")
    include("PublicHelpers.jl")
    include("LooCompare.jl")

end


text:
"""
This module calculates (generalized) hypergeometric functions:

    mFn(a;b;z) = ∑_{k=0}^∞ (a_1,k) ⋯ (a_m,k) / (b_1,k) ⋯ (b_n,k) zᵏ/k!
"""
module HypergeometricFunctions

using DualNumbers
import ApproxFun: @clenshaw, real, eps
import FastTransforms: pochhammer

export _₂F₁, _₃F₂, mFn

include("Gauss.jl")
include("specialfunctions.jl")

end #module

using .HypergeometricFunctions

export _₂F₁, _₃F₂, mFn

ext: jl


text:
using LinearAlgebra

"""
Compute the QAP objective function value.
"""
function qap_objective(A::AbstractArray{<:Real,2}, B::AbstractArray{<:Real,2}, P::Union{BitMatrix, Matrix{Bool}, Matrix{Int}})
	return tr(A * P * B' * P')
end

function qap_objective(A::AbstractArray{<:Real,2}, B::AbstractArray{<:Real,2}, perm::Array{Int,1})
	return sum(A .* B[perm,perm])
end

"""
Compute the match ratio of two permutations.
"""
function match_ratio(perm1::Array{Int,1}, perm2::Array{Int,1})
	@assert length(perm1) == length(perm2)
	@assert isperm(perm1) && isperm(perm2)
	N = length(perm1)
	mr = sum(perm1 .== perm2) / N
	return mr
end


text:
using BenchmarkTools
using Random
using Distances, LinearAlgebra

const SUITE = BenchmarkGroup()

Random.seed!(1234)

dim = 50
N1 = 1000;
N2 = 500;
alpha = 2.0

X = rand(Float64, N1, dim)
Y = rand(Float64, N2, dim)

KXY = rand(Float64, N1, N2)
KX = rand(Float64, N1, N1)
sKX = Symmetric(rand(Float64, N1, N1))
kX = rand(Float64, N1)

include("kernelmatrix.jl")
include("MLKernels.jl")


text:
abstract type abstractModulatedWaveFD <: AbstractPOFDProblem end
abstract type abstractModulatedWaveShooting <: AbstractShootingProblem end

"""
	pb = TWProblem(F, J, ∂::Tuple, u₀; DAE = 0)

This composite type implements a functional for freezing symmetries in order, for example, to compute travelling waves (TW). Note that you can freeze many symmetries, not just one, by passing many Lie generators. When you call `pb(x, par)`, it computes:

					┌                  ┐
					│ f(x, par) - s⋅∂x │
					│   <x - u₀, ∂u₀>  │
					└                  ┘

## Arguments
- `(x, p) -> F(x, p)` function with continuous symmetries
- `J` jacobian of `F`. Can be matrix based or matrix-free. The requirements are same as for [`newton`](@ref)
- `∂::Tuple = (T1, T2, ⋯)` tuple of Lie generators. In effect, each of these is an (differential) operator which can be specified as a (sparse) matrix or as an operator implementing `LinearAlgebra.mul!`.
- `u₀` reference solution

## Additional Constructor(s)

	pb = TWProblem(F, J, ∂, u₀; kw...)

This simplified call handles the case where a single symmetry needs to be frozen.

## Useful function

- `updateSection!(pb::TWProblem, u0)` updates the reference solution of the problem using `u0`.
- `nbConstraints(::TWProblem)` number of constraints (or Lie generators)

"""
@with_kw struct TWProblem{Tf, TJf, Tu0, TDu0, TD}
	F::Tf
	J::TJf = nothing
	∂::TD			# differential operator
	u₀::Tu0 		# reference solution, we only need one!
	∂u₀::TDu0 = (D * u₀,)
	DAE::Int = 0
	nc::Int = 1 	# number of constraints
	isinplace::Bool = false
	@assert 0 <= DAE <= 1
	@assert 0 < nc
end

function TWProblem(F, J, ∂::Tuple, u₀; DAE = 0, isinplace = false)
	# ∂u₀ = Tuple( apply(_D, u₀) for _D in ∂)
	∂u₀ = Tuple( mul!(zero(u₀), _D, u₀, 1, 0) for _D in ∂)
	return TWProblem(F = F, J = J, ∂ = ∂,
		u₀ = u₀,
		∂u₀ = ∂u₀,
		# u₀∂u₀ = Tuple( dot(u₀, u) for u in ∂u₀),
		DAE = DAE,
		isinplace = isinplace,
		nc = length(∂) )
end

# constructor
TWProblem(F, J, ∂, u₀; kw...) = TWProblem(F, J, (∂,), u₀; kw...)

@inline nbConstraints(pb::TWProblem) = pb.nc

# we put type information to ensure the user pass a correct u0
function updateSection!(pb::TWProblem{Tf, TJf, Tu0, TDu0, TD}, u₀::Tu0) where {Tf, TJf, Tu0, TDu0, TD}
	copyto!(pb.u₀, u₀)
	for (∂, ∂u₀) in zip(pb.∂, pb.∂u₀)
		# pb.u₀∂u₀ = Tuple( dot(u₀, u) for u in ∂u₀)
		copyto!(∂u₀, ∂ * u₀)
	end
end

"""
- `ss` tuple of speeds
- `D` tuple of Lie generators
"""
function applyD(pb::TWProblem, out, ss, u)
	for (D, s) in zip(pb.∂, ss)
		# out .-=  s .* (D * u)
		mul!(out, D, u, -s, 1)
	end
	out
end
applyD(pb::TWProblem, u) = applyD(pb, zero(u), 1, u)

# s is the speed.
# Return F(u, p) - s * D * u
@views function VFplusD(pb::TWProblem, u::AbstractVector, s::Tuple, pars)
	# apply the vector field
	out = pb.F(u, pars)
	# we add the freezing, it can be done now since out is filled by the previous call!!
	applyD(pb, out, s, u)
	return out
end

# function (u, p) -> F(u, p) - s * D * u to be used with shooting or Trapezoid
VFtw(pb::TWProblem, u::AbstractVector, parsFreez) = VFplusD(pb, u, parsFreez.s, parsFreez.user)

@views function (pb::TWProblem)(x::AbstractVector, pars)
	# number of constraints
	nc = pb.nc
	# number of unknowns
	N = length(x) - nc
	# array containing the result
	out = similar(x)
	u = x[1:N]
	outu = out[1:N]
	# get the speed
	s = Tuple(x[end-nc+1:end])
	# apply the vector field
	outu .= VFplusD(pb, u, s, pars)
	# we put the constraints
	for ii in 0:nc-1
		out[end-ii] = dot(u, pb.∂u₀[ii+1])
		if pb.DAE == 0
			out[end-ii] -= dot(pb.u₀, pb.∂u₀[ii+1])
		end
	end
	return out
end

@views function (pb::TWProblem)(x::AbstractVector, pars, dx::AbstractVector)
	# number of constraints
	nc = pb.nc
	# number of unknowns
	N = length(x) - nc
	# array containing the result
	out = similar(x)
	u = x[1:N]
	du = dx[1:N]
	outu = out[1:N]
	# get the speed
	s = Tuple(x[end-nc+1:end])
	ds = Tuple(dx[end-nc+1:end])
	# get the jacobian
	J = pb.J(u, pars)
	outu .= apply(J, du)
	applyD(pb, outu, s, du)
	applyD(pb, outu, ds, u)
	# we put the constraints
	for ii in 0:nc-1
		out[end-ii] = dot(du, pb.∂u₀[ii+1])
	end
	return out
end

# build the sparse jacobian of the freezed problem
function (pb::TWProblem)(::Val{:JacFullSparse}, ufreez::AbstractVector, par; δ = 1e-9)
	# number of constraints
	nc = nbConstraints(pb)
	# number of unknowns
	N = length(ufreez) - nc
	# get the speed
	s = Tuple(ufreez[end-nc+1:end])
	# get the state space vector
	u = ufreez[1:N]
	# the jacobian of the
	J1 = pb.J(u, par)
	# we add the Lie algebra generators
	rightpart = zeros(N, nc)
	for ii in 1:nc
		J1 = J1 - s[ii] * pb.∂[ii]
		mul!(view(rightpart, :, ii), pb.∂[ii], u, -1, 0)
	end
	J2 = hcat(J1, rightpart)
	for ii in 1:nc
		J2 = vcat(J2, vcat(pb.∂u₀[ii], zeros(nc))')
	end
	return J2
end

function newton(prob::TWProblem, orbitguess, par, optn::NewtonPar;
		jacobian = :MatrixFree, kwargs...)
	@assert jacobian in  (:MatrixFree, :MatrixFreeAD, :AutoDiff, :FullLU)
	if jacobian == :AutoDiff
		jac = (x, p) -> sparse(ForwardDiff.jacobian(z -> prob(z, p), x))
	elseif jacobian == :MatrixFreeAD
		jac = (x, p) -> (dx -> ForwardDiff.derivative(t -> prob(x .+ t .* dx, p), 0))
	elseif jacobian == :FullLU
		jac = (x, p) -> prob(Val(:JacFullSparse), x, p)
	# elseif jacobian == :FullSparseInplace
	elseif jacobian == :MatrixFree
		jac = (x, p) -> (dx ->  prob(x, p, dx))
	end
	return newton(prob, jac, orbitguess, par, optn; kwargs...,)
end

function continuation(prob::TWProblem,
		orbitguess, par, lens::Lens, contParams::ContinuationPar;
		jacobian = :MatrixFree, kwargs...)
	@assert jacobian in (:MatrixFree, :MatrixFreeAD, :AutoDiff, :FullLU)

	if jacobian == :AutoDiff
		jac = (x, p) -> sparse(ForwardDiff.jacobian(z -> prob(z, p), x))
	elseif jacobian == :MatrixFreeAD
		jac = (x, p) -> (dx -> ForwardDiff.derivative(t -> prob(x .+ t .* dx, p), 0))
	elseif jacobian == :FullLU
		jac = (x, p) -> prob(Val(:JacFullSparse), x, p)
	# elseif jacobian == :FullSparseInplace
	elseif jacobian == :MatrixFree
		jac = (x, p) -> (dx ->  prob(x, p, dx))
	end
	# define the mass matrix for the eigensolver
	N = length(orbitguess)
	B = Diagonal(vcat(ones(N-1),0))
	# convert eigsolver to generalised one
	old_eigsolver = contParams.newtonOptions.eigsolver
	contParamsWave = @set contParams.newtonOptions.eigsolver = convertToGEV(old_eigsolver, B)
	# call continuation
	branch, u, τ = continuation(prob, jac, orbitguess, par, lens, contParamsWave; kwargs...,)
	return setproperties(branch; type = :TravellingWave, functional = prob), u, τ
end


text:
# indicator of the zero cone

"""
  IndZero()

Returns the indicator function of the zero point, or "zero cone", i.e.,
  `g(x) = 0 if x = 0, +∞ otherwise`
"""

immutable IndZero <: IndicatorConvexCone end

function (f::IndZero){T <: RealOrComplex}(x::AbstractArray{T})
  for k in eachindex(x)
    if x[k] != zero(T)
      return Inf
    end
  end
  return 0.0
end

function prox!{T <: RealOrComplex}(f::IndZero, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
  for k in eachindex(x)
    y[k] = zero(T)
  end
  return 0.0
end

fun_name(f::IndZero) = "indicator of the zero cone"
fun_dom(f::IndZero) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndZero) = "x ↦ 0 if all(x = 0), +∞ otherwise"
fun_params(f::IndZero) = "none"

function prox_naive(f::IndZero, x::AbstractArray, gamma::Real=1.0)
  return zero(x), 0.0
end


text:
module YaoToEinsum

using Yao, OMEinsum
using LinearAlgebra

export yao2einsum, add_gate!

struct EinBuilder{T}
    slots::Vector{Int}
    labels::Vector{Vector{Int}}
    tensors::Vector{AbstractArray{T}}
    maxlabel::Base.RefValue{Int}
end

Yao.nqubits(eb::EinBuilder) = length(eb.slots)
function add_tensor!(eb::EinBuilder{T}, tensor::AbstractArray{T,N}, labels::Vector{Int}) where {N,T}
    @assert N == length(labels)
    push!(eb.tensors, tensor)
    push!(eb.labels, labels)
end

function EinBuilder(::Type{T}, n::Int) where T
    EinBuilder(collect(1:n), Vector{Int}[], AbstractArray{T}[], Ref(n))
end
newlabel!(eb::EinBuilder) = (eb.maxlabel[] += 1; eb.maxlabel[])

function add_gate!(eb::EinBuilder{T}, b::PutBlock{D,C}) where {T,D,C}
    return add_matrix!(eb, C, mat(T, b.content), collect(b.locs))
end
# general and diagonal gates
function add_matrix!(eb::EinBuilder{T}, k::Int, m::AbstractMatrix, locs::Vector) where T
    if !isdiag(m)
        nlabels = [newlabel!(eb) for _=1:k]
        add_tensor!(eb, reshape(Matrix{T}(m), fill(2, 2k)...), [nlabels..., eb.slots[locs]...])
        eb.slots[locs] .= nlabels
    else
        add_tensor!(eb, reshape(Vector{T}(diag(m)), fill(2, k)...), eb.slots[locs])
    end
    return eb
end
# swap gate
function add_gate!(eb::EinBuilder{T}, b::PutBlock{2,2,ConstGate.SWAPGate}) where {T,N}
    lj = eb.slots[b.locs[2]]
    eb.slots[b.locs[2]] = eb.slots[b.locs[1]]
    eb.slots[b.locs[1]] = lj
    return eb
end

# control gates
function add_gate!(eb::EinBuilder{T}, b::ControlBlock{BT,C,M}) where {T, BT,C,M}
    return add_controlled_matrix!(eb, M, mat(T, b.content), collect(b.locs), collect(b.ctrl_locs), collect(b.ctrl_config))
end
function add_controlled_matrix!(eb::EinBuilder{T}, k::Int, m::AbstractMatrix, locs::Vector, control_locs, control_vals) where T
    if length(control_locs) == 0
        return add_matrix!(eb, k, m, locs)
    end
    sig = eb.slots[control_locs[1]]
    val = control_vals[1]
    for i=1:length(control_locs)-1
        newsig = newlabel!(eb)
        add_tensor!(eb, and_gate(T, control_vals[i+1], val), [newsig,eb.slots[control_locs[i+1]],sig])
        sig = newsig
        val = 1
    end
    if !isdiag(m)
        t1 = reshape(Matrix{T}(m), fill(2, 2k)...)
        t2 = reshape(Matrix{T}(I, 1<<k, 1<<k), fill(2, 2k)...)
        if val == 1
            t1, t2 = t2, t1
        end
        nlabels = [newlabel!(eb) for _=1:k]
        add_tensor!(eb, cat(t1, t2; dims=2k+1), [nlabels..., eb.slots[locs]..., sig])
        eb.slots[locs] .= nlabels
    else
        t1 = reshape(Vector{T}(diag(m)), fill(2, k)...)
        t2 = reshape(ones(T, 1<<k), fill(2, k)...)
        if val == 1
            t1, t2 = t2, t1
        end
        add_tensor!(eb, cat(t1, t2; dims=k+1), [eb.slots[locs]..., sig])
    end
    return eb
end

function and_gate(::Type{T}, a::Int, b::Int) where T
    m = zeros(T, 2, 2, 2)
    for v1 in (0, 1)
        for v2 in (0, 1)
            # the first is output
            m[(v1==a && v2==b)+1, v1+1,v2+1] = 1
        end
    end
    return m
end

function add_gate!(eb::EinBuilder, b::ChainBlock)
    for ib in subblocks(b)
        add_gate!(eb, ib)
    end
    return eb
end

function add_gate!(eb::EinBuilder, b::AbstractBlock)
    B = Optimise.to_basictypes(b)
    if typeof(B) == typeof(b)
        throw("block of type `$(typeof(b))` can not be converted to tensor network representation!")
    else
        add_gate!(eb, B)
    end
    return eb
end

"""
    yao2einsum(circuit, initial_state::Dict, final_state::Dict)
    yao2einsum(circuit; initial_state=Dict(), final_state=Dict())

Transform a Yao `circuit` to a generalized tensor network (einsum) notation.
This function returns a 2-tuple of (einsum code, input tensors).
`initial_state` and `final_state` are dictionaries that specify the initial state and final state as product states,
e.g. a vector `Dict(1=>1, 2=>1, 3=>0, 4=>1)` specifies a product state `|1⟩⊗|1⟩⊗|0⟩⊗|1⟩`.
If an qubit in initial state or final state is not specified, it will be treated as an open edge.

```jldoctest
julia> using YaoToEinsum, Yao

julia> c = chain(3, put(3, 2=>X), put(3, 1=>Y), control(3, 1, 3=>Y))
nqubits: 3
chain
├─ put on (2)
│  └─ X
├─ put on (1)
│  └─ Y
└─ control(1)
   └─ (3,) Y

julia> yao2einsum(c; initial_state=Dict(1=>0, 2=>1), final_state=Dict(1=>ArrayReg([0.6, 0.8im]), 2=>1))
(1, 2, 4∘2, 5∘1, 6∘3∘5, 5, 4 -> 6∘3, AbstractArray{ComplexF64}[[1.0 + 0.0im, 0.0 + 0.0im], [0.0 + 0.0im, 1.0 + 0.0im], [0.0 + 0.0im 1.0 + 0.0im; 1.0 + 0.0im 0.0 + 0.0im], [0.0 + 0.0im 0.0 - 1.0im; 0.0 + 1.0im 0.0 + 0.0im], [1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0 + 0.0im;;; 0.0 + 0.0im 0.0 - 1.0im; 0.0 + 1.0im 0.0 + 0.0im], [0.6 + 0.0im, 0.0 + 0.8im], [0.0 + 0.0im, 1.0 + 0.0im]])
```
"""
function yao2einsum(circuit::AbstractBlock{D}; initial_state::Dict=Dict{Int,Int}(), final_state::Dict=Dict{Int,Int}()) where {D}
    T = promote_type(ComplexF64, dict_regtype(initial_state), dict_regtype(final_state), Yao.parameters_eltype(circuit))
    vec_initial_state = Dict{Int,ArrayReg{D,T}}([k=>render_single_qubit_state(T, v) for (k, v) in initial_state])
    vec_final_state = Dict{Int,ArrayReg{D,T}}([k=>render_single_qubit_state(T, v) for (k, v) in final_state])
    yao2einsum(circuit, vec_initial_state, vec_final_state)
end
dict_regtype(d::Dict) = promote_type(_regtype.(values(d))...)
_regtype(::ArrayReg{D,VT}) where {D,VT} = VT
_regtype(::Int) = ComplexF64
render_single_qubit_state(::Type{T}, x::Int) where T = x == 0 ? zero_state(T, 1) : product_state(T, bit"1")
render_single_qubit_state(::Type{T}, x::ArrayReg) where T = ArrayReg(collect(T, statevec(x)))

function yao2einsum(circuit::AbstractBlock{D}, initial_state::Dict{Int,<:ArrayReg{D,T}}, final_state::Dict{Int,<:ArrayReg{D,T}}) where {D,T}
    n = nqubits(circuit)
    eb = EinBuilder(T, n)
    openindices = Int[]
    for k=1:n
        if haskey(initial_state, k)
            add_tensor!(eb, statevec(initial_state[k]), [eb.slots[k]])
        else
            push!(openindices, eb.slots[k])
        end
    end
    add_gate!(eb, circuit)
    openindices2 = Int[]
    for k=1:n
        if haskey(final_state, k)
            add_tensor!(eb, statevec(final_state[k]), [eb.slots[k]])
        else
            push!(openindices2, eb.slots[k])
        end
    end
    return build_einsum(eb, vcat(openindices2, openindices))
end

function build_einsum(eb::EinBuilder, openindices)
    return EinCode(eb.labels, openindices), eb.tensors
end

end


text:
module Paddle

using PyCall

using ChainRulesCore
using DLPack
using Functors: @functor
using Adapt

using ..PaddleChainRules: PyAdaptor, fmap

const dlpack = PyNULL()
const paddle = PyNULL()
const ispysetup = Ref{Bool}(false)

pyto_dlpack(x) = @pycall dlpack.to_dlpack(x)::PyObject
pyfrom_dlpack(x) = @pycall dlpack.from_dlpack(x)::PyObject

include("Net.jl")

struct PaddleModuleWrapper
    NN::PaddleStatelessModule
    dtype::Type
    params::Vector
end

Base.show(io::IO, f::PaddleModuleWrapper) = print(io, f.NN, " ", f.dtype)

Base.length(f::PaddleModuleWrapper) = length(f.params)
Base.iterate(f::PaddleModuleWrapper) = iterate(f.params)
Base.iterate(f::PaddleModuleWrapper, state) = iterate(f.params, state)

function PaddleModuleWrapper(paddle_module)
    params = paddle_module.parameters()
    jlparams = map(x->DLPack.wrap(x, pyto_dlpack), params)
    dtype = eltype(jlparams[1])
    return PaddleModuleWrapper(PaddleStatelessModule(paddle_module), dtype, jlparams)
end

@functor PaddleModuleWrapper (params,)

function (wrap::PaddleModuleWrapper)(args...; kwargs...)
    out = wrap.NN(map(x -> DLPack.share(x, PyObject, pyfrom_dlpack), wrap.params), fmap(x -> DLPack.share(x, PyObject, pyfrom_dlpack), args)...)
    res = fmap(x->DLPack.wrap(x, pyto_dlpack), out)
    return res
end

function vjp(stateless_module::PaddleStatelessFCNet, pyparams::Vector, pyargs...; kwargs...) # grad wrt params and args
    res = stateless_module(pyparams, pyargs...; kwargs...)
    paramslen = length(pyparams)
    function vjp_func(Δ)
        # compute the grad of both params and args? will it cost more during training?
        grad = paddle.fluid.dygraph.grad(res, vcat(pyparams, pyargs...), Δ, retain_graph=true)
        return (grad[1:paramslen], grad[paramslen+1:end])
    end
    return res, vjp_func
end

function vjp(stateless_module::PaddleStatelessGeneralNet, pyparams::Vector, pyargs...; kwargs...) # grad wrt params and args
    res = stateless_module(pyparams, pyargs...; kwargs...)
    paramslen = length(pyparams)
    function vjp_func(Δ)
        # compute the grad of both params and args? will it cost more during training?
        grad = paddle.fluid.dygraph.grad(res, vcat(stateless_module.NN.parameters(), pyargs...), Δ, retain_graph=true)
        return (grad[1:paramslen], grad[paramslen+1:end])
    end
    return res, vjp_func
end

## grad wrt params only
# function vjp_wrt_params(stateless_module::PaddleStatelessModule, pyparams::Vector, pyargs...; kwargs...)
#     res = stateless_module(pyparams, pyargs...; kwargs...)
#     function vjp_func(Δ)
#         # compute the grad of params
#         grad = paddle.fluid.dygraph.grad(res, pyparams, Δ, retain_graph=true)
#         return grad
#     end
#     return res, vjp_func
# end

function to_paddle_tensor(x::AbstractArray, stop_grad::Bool)
    out = DLPack.share(x, PyObject, pyfrom_dlpack)
    out.stop_gradient = stop_grad
    return out
end

function ChainRulesCore.rrule(wrap::PaddleModuleWrapper, args...; kwargs...)
    T = typeof(first(wrap.params))
    pyparams = fmap(x -> to_paddle_tensor(x,false), wrap.params)
    pyargs = fmap(x -> to_paddle_tensor(x,false), args)
    paddle_primal, vjp_func = vjp(wrap.NN, pyparams, pyargs...; kwargs...)
    project = ProjectTo(args)
    function pullback(Δ)
        cΔ = fmap(x->Adapt.adapt(PyAdaptor{T}(), x), Δ)
        pycΔ = fmap(x->DLPack.share(x, PyObject, pyfrom_dlpack), cΔ)
        paddle_tangent_vals = vjp_func(pycΔ)
        jlparams_tangents = map(x -> DLPack.wrap(x, pyto_dlpack), paddle_tangent_vals[1])
        args_tangents = project(fmap(x -> DLPack.wrap(x, pyto_dlpack), paddle_tangent_vals[2]))
        return (Tangent{PaddleModuleWrapper}(; NN = NoTangent(), dtype = NoTangent(), params = jlparams_tangents), args_tangents...)
    end
    res = fmap(x->DLPack.wrap(x, pyto_dlpack), paddle_primal)
    return res, pullback
end

function __init__()
    try
        copy!(paddle, pyimport("paddle"))
        copy!(dlpack, pyimport("paddle.utils.dlpack"))
        ispysetup[] = true
    catch err
        @warn """PaddleChainRules.jl has failed to import paddle from Python.
                 Please make sure these are installed.
        """
        @debug err
        ispysetup[] = false
    end
end

end

max_stars_repo_name: PangHua/InfiniteScience

max_issues_repo_name: PangHua/InfiniteScience

max_forks_repo_name: PangHua/InfiniteScience


text:
using DataGraphs
using Documenter

DocMeta.setdocmeta!(DataGraphs, :DocTestSetup, :(using DataGraphs); recursive=true)

makedocs(;
  modules=[DataGraphs],
  authors="Matthew Fishman <mfishman@flatironinstitute.org> and contributors",
  repo="https://github.com/mtfishman/DataGraphs.jl/blob/{commit}{path}#{line}",
  sitename="DataGraphs.jl",
  format=Documenter.HTML(;
    prettyurls=get(ENV, "CI", "false") == "true",
    canonical="https://mtfishman.github.io/DataGraphs.jl",
    assets=String[],
  ),
  pages=["Home" => "index.md"],
)

deploydocs(; repo="github.com/mtfishman/DataGraphs.jl", devbranch="main")


text:
__precompile__()

module CheckImports

"""
`modulechildren(m::Module)`

Returns the modules in m
"""
modulechildren(m::Module) = filter(x->isa(x, Module), map(x->m.(x), names(m, true)))

"""
`isimported(pkg::AbstractString;m::Module=Main)`

Returns true if module "pkg" is defined in m, otherwise false. m defaults to Main
"""
isimported(pkg::AbstractString;m::Module=Main)= maximum(map(string,modulechildren(m)).==pkg)

export modulechildren, isimported

end # module

max_stars_repo_name: ChrisRackauckas/CheckImports

max_issues_repo_name: ChrisRackauckas/CheckImports

max_forks_repo_name: ChrisRackauckas/CheckImports


text:
using DiscreteFunctions, SimpleDrawing, Plots, SimpleGraphs

import SimpleDrawing.draw

function func2graph(f::DiscreteFunction)
    n = length(f)
    G = IntGraph(n)
    for a=1:n
        b = f(a)
        if a!=b
            add!(G,a,b)
        end
    end
    embed(G,:combined)
    return G
end

"""
`draw(f::DiscreteFunction)` draws a directed graph representation of `f`.
"""
function draw(f::DiscreteFunction)
    n = length(f)
    newdraw()

    G = func2graph(f)
    xy = G.cache[:GraphEmbedding].xy

    # xy = [ (n*sin(2*pi*t/n),n*cos(-2*pi*t/n)) for t=0:n-1 ]

    mu = 0.9
    for j=1:n
        fj = f(j)
        a,b = xy[j]
        c,d = xy[fj]
        if fj!=j
            draw_vector(mu*(c-a),mu*(d-b),a,b,color=:black)
        else
            draw_circle(a,b,0.1,color=:black)
        end
    end

    for j=1:n
        x,y = xy[j]
        # draw_point(x,y,color=:black,marker=1)
        annotate!(x,y,text(string(j),8,color=:red))
    end

    finish()
end


text:
#############################################################################
# CoordinatedPrecoding
# Description here
# See http://github.com/rasmusbrandt/CoordinatedPrecoding.jl
#############################################################################

VERSION >= v"0.4" && __precompile__()

module CoordinatedPrecoding

using Compat
import PyPlot, Lumberjack, ProgressMeter

# Operator extension
import Base: +, -

export

# Channels
    Channel,

    SinglecarrierChannel, MulticarrierChannel,

    get_channel_gains_dB,
    get_average_channel_gains_dB,

# Networks
    Network,

    InterferingBroadcastChannel, setup_interfering_broadcast_channel,
    IndoorsNetwork, setup_indoors_network,
    RandomLargeScaleNetwork, setup_random_large_scale_network,
    TriangularMacroNetwork, setup_triangularmacro_network,
    TriangularHetNetNetwork, setup_triangularhetnet_network,

    draw_user_drop!, draw_channel, plot_network_layout,

    get_num_MSs, get_num_BSs,

    get_assignment, set_assignment!,
    has_aux_network_param,
    get_aux_network_param, set_aux_network_param!,
    get_aux_network_params, set_aux_network_params!,
    has_aux_precoding_param,
    get_aux_precoding_param, set_aux_precoding_param!,
    get_aux_precoding_params, set_aux_precoding_params!,
    has_aux_assignment_param,
    get_aux_assignment_param, set_aux_assignment_param!,
    get_aux_assignment_params, set_aux_assignment_params!,

    get_num_MS_antennas, get_num_BS_antennas,
    get_transmit_powers, set_transmit_powers!,
    get_transmit_powers_dBm, set_transmit_powers_dBm!,
    get_receiver_noise_powers, set_receiver_noise_powers!,
    get_receiver_noise_powers_dBm, set_receiver_noise_powers_dBm!,

    get_num_streams, set_num_streams!,
    get_user_priorities, set_user_priorities!,

    get_average_SNRs, set_average_SNRs!,
    get_average_SNRs_dB, set_average_SNRs_dB!,

    require_equal_num_MS_antennas,
    require_equal_num_BS_antennas,
    require_single_antenna_MSs,
    require_single_antenna_BSs,
    require_equal_num_streams,
    require_single_stream,

# Assignment
    Assignment,

    IDCellAssignment!,
    LargeScaleFadingCellAssignment!,

    AssignmentResults,
    AuxAssignmentParams,

    active_BSs,
    num_served_MSs,
    serving_BS_id,
    served_MS_ids,
    served_MS_ids_except_me,

    coordinated_BS_ids,
    coordinated_MS_ids,

    require_equal_num_MSs_per_cell,

# Precoding
    Eigenprecoding,
    Gomadam2008_MaxSINR,
    Gomadam2008_MinWLI,
    Komulainen2013_WMMSE,
    Razaviyayn2013_MinMaxWMMSE,
    Shi2011_WMMSE,

    PrecodingResults,
    AuxPrecodingParams,

    initial_receivers,
    initial_MSE_weights,
    initial_precoders,

    calculate_logdet_rates,
    calculate_MMSE_rates,
    calculate_weighted_logdet_rates,
    calculate_weighted_MMSE_rates,
    calculate_allocated_power,

# Simulation
    simulate,
    postprocess,
    plot,

    simulate_assignment_convergence,
    postprocess_assignment_convergence,
    plot_assignment_convergence,

    simulate_precoding_convergence,
    postprocess_precoding_convergence,
    plot_precoding_convergence,

    timing,

# Utilities
    @defaultize_param!,
    clean_simulation_params_for_jld,
    Partition, Block, PartitionIterator, restricted_growth_string

##########################################################################
# General includes

# Networks
include("network/network.jl")

# Channels
include("channel.jl")

# Simulation
include("simulation/simulation.jl")

# Cell/cluster assignment
include("assignment.jl")

# Precoding algorithms
include("precoding/precoding.jl")

# Miscellaneous utilities
include("misc/partitions.jl")
include("misc/utils.jl")

##########################################################################
# Specific includes

# Networks
include("network/InterferingBroadcastChannel.jl")
include("network/IndoorsNetwork.jl")
include("network/RandomLargeScaleNetwork.jl")
include("network/TriangularMacroNetwork.jl")
include("network/TriangularHetNetNetwork.jl")

# Precoding methods
include("precoding/Eigenprecoding.jl")
include("precoding/Gomadam2008_MaxSINR.jl")
include("precoding/Gomadam2008_MinWLI.jl")
include("precoding/Komulainen2013_WMMSE.jl")
include("precoding/Razaviyayn2013_MinMaxWMMSE.jl")
include("precoding/Shi2011_WMMSE.jl")

# Simulation utilities
include("simulation/timing.jl")
include("simulation/visualization.jl")

end


text:
using LinearAlgebra
using PyPlot
using Statistics
using Distributions
using ProgressMeter

σ = 1
τ = sqrt(10)
μ = 5
n = 5

X = [9.37;10.18;9.16;11.60;10.33]
μn = mean(X)*n/σ^2/(n/σ^2 + 1/τ^2) + μ*(1/τ^2)/(n/σ^2 + 1/τ^2)
τn2 = 1/(n/σ^2+1/τ^2)

# Hyperparameters
# proposal step size
δ = 3/100
burnin = 2000
N = 10000
θ0 = 0.0

function logf(θ)
    -sum((X.-θ).^2)/2σ^2 - (θ-μ)^2/2τ^2
end

function proposal(x)
    x + (rand()-0.5)*2 * δ
end

sim = zeros(N)
sim[1] = θ0
@showprogress for i = 1:N-1
    x = sim[i]
    x_star = proposal(x)
    Δ =  logf(x_star) - logf(x)
    if log(rand())<Δ
        sim[i+1] = x_star
    else
        sim[i+1] = x
    end
end
L = logf.(sim)

# sim = sim[burnin+1:end]
figure(figsize=(12,4))
subplot(131)
plot(sim)
plot(1:length(sim), ones(length(sim))*μn, "--")
title("\$\\theta\$ Value")
xlabel("Iteration")
ylabel("\$\\theta\$")
subplot(132)
plot(L)
title("Log likelihood")
xlabel("Iteration")
ylabel("Log likelihood")
subplot(133)
hist(sim[burnin+1:end], density=true, bins=50)
y = pdf.(Normal(μn, sqrt(τn2)), LinRange(8.5,11.5,100))
plot(LinRange(8.5,11.5,100), y, label="Exact")
title("Distribution")
xlabel("\$\\theta\$")
ylabel("Density")
legend()
tight_layout()


text:
#
# Protein - TMAO
#

using MDDF

# Ploting
using Plots
ENV["GKSwstype"]="nul" # This supresses the need of a display while plotting

# Here we use the PDBTools package to read the pdb file and select the
# solute and solvent atoms (from http://github.com/m3g/PDBTools)
using PDBTools
atoms = PDBTools.readPDB("../NAMD/structure.pdb")

# The solute is a single protein molecule (infinte dilution case). In this case,
# use the option nmols=1
solute_indexes = PDBTools.selindex(atoms,"protein")
solute = MDDF.Selection( solute_indexes, nmols=1 )

# The solvent is TMAO, which has 14 atoms. Use the natomspermol to indicate how many
# atoms each molecule has, such that there is no ambiguity on how to split the coordinates
# of the selection into individual molecules.
solvent_indexes = PDBTools.selindex(atoms,"resname TMAO")
solvent = MDDF.Selection( solvent_indexes, natomspermol=14 )

# Input options for the calcualtion
options = MDDF.Options(binstep=0.2)

# Run MDDF calculation, and get the resutls in the R structure. Here we use the
# Chemfiles package to read the trajectory. This package reads the trajectories
# in many common simulation formats
trajectory = MDDF.ChemFile("../NAMD/trajectory.dcd",solute,solvent)
@time R = MDDF.mddf(trajectory,options)

plot(layout=(6,1))

sp=1
plot!(ylabel="MDDF",subplot=sp)
plot!(R.d,R.mddf,subplot=sp,label="mddf")
plot!(R.d,R.rdf,subplot=sp,label="rdf")
plot!(legend=:topright,subplot=sp)

sp=2
plot!(ylabel="KB",subplot=sp)
plot!(R.d,R.kb,subplot=sp,label="mddf")
plot!(R.d,R.kb_rdf,subplot=sp,label="rdf")
plot!(legend=:topright,subplot=sp)

sp=3
plot!(ylabel="Count",subplot=sp)
scatter!(R.d,R.md_count,subplot=sp,label="md")
scatter!(R.d,R.md_count_random,subplot=sp,label="rand")

sp=4
plot!(ylabel="Shell vol", subplot=sp)
plot!(R.d,R.volume.shell,subplot=sp,label="")

sp=5
plot!(ylabel="Sum MD", subplot=sp)
scatter!(R.d,R.sum_md_count,subplot=sp,label="md")

sp=6
plot!(ylabel="Sum RAND", subplot=sp)
scatter!(R.d,R.sum_md_count_random,subplot=sp,label="rand")
scatter!(R.d,R.sum_rdf_count,subplot=sp,label="rdf")
plot!(legend=:topleft,subplot=sp)

plot!(size=(800,1300))
savefig("./plots.pdf")


text:
"""
    grdimage(cmd0::String="", arg1=[], arg2=[], arg3=[]; kwargs...)

Produces a gray-shaded (or colored) map by plotting rectangles centered on each grid node and assigning them a gray-shade (or color) based on the z-value.

Full option list at [`grdimage`](http://gmt.soest.hawaii.edu/doc/latest/grdimage.html)

Parameters
----------

- **A** : **img_out** : **image_out** : -- Str --

    Save an image in a raster format instead of PostScript.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/grdimage.html#a)
- $(GMT.opt_J)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **D** : **img_in** : **image_in** : -- Str or [] --

    Specifies that the grid supplied is an image file to be read via GDAL.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/grdimage.html#d)
- **E** : **dpi** : -- Int or [] --

    Sets the resolution of the projected grid that will be created.
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/grdimage.html#e)
- **G** : -- Str or Int --

    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/grdimage.html#g)
- **I** : **shade** : **intensity** : **intensfile** : -- Str or GMTgrid --

    Gives the name of a grid file or GMTgrid with intensities in the (-1,+1) range,
    or a grdgradient shading flags.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/grdimage.html#i)
- **M** : **monochrome** : -- Bool or [] --

    Force conversion to monochrome image using the (television) YIQ transformation.
    [`-M`](http://gmt.soest.hawaii.edu/doc/latest/grdimage.html#m)
- **N** : **noclip** : -- Bool or [] --

    Do not clip the image at the map boundary.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/grdimage.html#n)
- $(GMT.opt_P)
- **Q** : **nan_t** : **nan_alpha** : -- Bool or [] --

    Make grid nodes with z = NaN transparent, using the colormasking feature in PostScript Level 3.
- $(GMT.opt_R)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_f)
- $(GMT.opt_n)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function grdimage(cmd0::String="", arg1=[], arg2=[], arg3=[], arg4=[]; K=false, O=false, first=true, kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("grdimage", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

	cmd, opt_B, = parse_BJR(d, "", "", O, " -JX12c/0")
	cmd = parse_UVXY(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_n(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd  = parse_params(cmd, d)

	cmd, K, O, = set_KO(cmd, opt_B, first, K, O)			# Set the K O dance

	cmd = add_opt(cmd, 'A', d, [:A :img_out :image_out])
	cmd = add_opt(cmd, 'D', d, [:D :img_in :image_in])
	cmd = add_opt(cmd, 'E', d, [:E :dpi])
	cmd = add_opt(cmd, 'G', d, [:G])
	cmd = add_opt(cmd, 'M', d, [:M :monochrome])
	cmd = add_opt(cmd, 'N', d, [:N :noclip])
	cmd = add_opt(cmd, 'Q', d, [:Q :nan_t :nan_alpha])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)		# Find how data was transmitted
	if (got_fname == 0 && isempty_(arg1))			# Than it must be using the three r,g,b grids
		cmd, got_fname, arg1, arg2, arg3 = find_data(d, cmd0, cmd, 3, arg1, arg2, arg3)
		if (got_fname == 0 && isempty_(arg1))
			error("No input data to use in grdimage.")
		end
	end

	if (isa(arg1, Array{<:Number}))
		arg1 = mat2grid(arg1)
		if (!isempty_(arg2) && isa(arg2, Array{<:Number}))  arg2 = mat2grid(arg2)  end
		if (!isempty_(arg3) && isa(arg3, Array{<:Number}))  arg3 = mat2grid(arg3)  end
	end

	N_used = got_fname == 0 ? 1 : 0		# To know whether a cpt will go to arg1 or arg2
	cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', N_used, arg1, arg2)
	if (!isempty_(arg3) && occursin("-C", cmd))		# This lieves out the case when the r,g,b were sent as a text.
		error("Cannot use the three R,G,B grids and a color table.")
	end

	for sym in [:I :shade :intensity :intensfile]
		if (haskey(d, sym))
			if (!isa(d[sym], GMTgrid))		# Uff, simple. Either a file name or a -A type modifier
				cmd = cmd * " -I" * arg2str(d[sym])
			else
				cmd, N = put_in_slot(cmd, d[sym], 'I', [arg1, arg2, arg3, arg4])
				if (N == 1)     arg1 = d[sym]
				elseif (N == 2) arg2 = d[sym]
				elseif (N == 3) arg3 = d[sym]
				elseif (N == 4) arg4 = d[sym]
				end
			end
			break
		end
	end

	cmd = finish_PS(d, cmd, output, K, O)
    return finish_PS_module(d, cmd, "", output, fname_ext, opt_T, K, "grdimage", arg1, arg2, arg3, arg4)
end

# ---------------------------------------------------------------------------------------------------
grdimage!(cmd0::String="", arg1=[], arg2=[], arg3=[], arg4=[]; K=true, O=true, first=false, kw...) =
	grdimage(cmd0, arg1, arg2, arg3, arg4; K=true, O=true, first=false, kw...)

grdimage(arg1, cmd0::String="", arg2=[], arg3=[], arg4=[]; K=false, O=false, first=true, kw...) =
	grdimage(cmd0, arg1, arg2, arg3, arg4; K=K, O=O, first=first, kw...)

grdimage!(arg1, cmd0::String="", arg2=[], arg3=[], arg4=[]; K=true, O=true, first=false, kw...) =
	grdimage(cmd0, arg1, arg2, arg3, arg4; K=K, O=O, first=first, kw...)


text:
using Waterfall
include(joinpath(WATERFALL_DIR,"bin","io.jl"))
include(joinpath(WATERFALL_DIR,"src","utils","correlate.jl"))

function cutoff(df::DataFrames.DataFrame; nstep=6, value=:Value)
    ii = [1; sortperm(abs.(df[2:end-1,value]); rev=true)[1:nstep].+1; size(df,1)]
    ii = 1:nstep
    return df[ii,:]
end

function plot(x::Waterfall.Cascade)
    A = copy.(x.correlation-I)
    df = convert(DataFrames.DataFrame, Waterfall.rowprod!(x))
    N = size(df,2)

    StatsPlots.gr(size=(1000,1000), markerstrokewidth=1)

    s = StatsPlots.@df df StatsPlots.corrplot(cols(1:N),
        grid = false,
        markeralpha = 0.4,
        tickfontsize = 8,
        xaxis = (tickfontrotation = 60.0),
    )

    [StatsPlots.annotate!(s[ii,jj],
        Statistics.mean(StatsPlots.xlims(s[ii,jj])),
        StatsPlots.ylims(s[ii,jj])[end],
        StatsPlots.text(Printf.@sprintf("%.2g", A[ii-1,jj-1]), :center, :middle, 8)
        )
        for ii in 3:N-1 for jj in 2:ii-1 if A[ii-1,jj-1]!==0.0
    ]

    StatsPlots.savefig("correlation.png")
end

nstep = 4
nsample = 1000
nperm = 10
ncor = 2

plot(Waterfall.Cascade(cutoff(df; nstep=nstep);
    interactivity=(0.5,0.75),
    nsample=nsample,
    ncor=ncor,
    kwargs...),
)

# plot_correlation(cascade)
# plot_cascade(cascade)


text:
"""
    Regions

Main module for Regions.jl - a set of types that model a discrete 2-dimensional region concept.

# Exports

* Run
* Region

# Dependencies

* Images.jl

"""
module Regions

using Images

export binarize, components

include("range.jl")
include("run.jl")
include("region.jl")
include("region_vector.jl")

"""
    binarize(image, predicate)

Binarize an image and return a region.

The predicate must be a funtion that takes a pixel value and returns a boolean result based
on the pixel value.

Here are some useful examples:

# the returned region will contain all pixels > 0.3
reg = binarize(img, x -> x > 0.3)

# the returned region will contain all pixels <= 0.5
reg = binarize(img, x -> x <= 0.5)

# the returned region will contain all pixels between 0.3 and 0.8
reg = binarize(img, x -> x > 0.3 && x < 0.8)

"""
function binarize(image, predicate)
    region = Region(Run[], false)

    rows, columns = size(image)

    for row in 1:rows
        inside_object = false
        start_column = 0
        for column in 1:columns
            if predicate(image[row, column])
                if false == inside_object
                    inside_object = true
                    start_column = column
                end
            else
                if true == inside_object
                    inside_object = false
                    push!(region.runs, Run(row, start_column:(column-1)))
                end
            end
        end
        # if still inside at the end of a line...
        if true == inside_object
            push!(region.runs, Run(row, start_column:columns))
        end
    end

    return region
end

"""
    components(region::Region, dx::Unsigned=1, dy::Unsigned=1)

This function splits a region into connected component objects that are returned as a vector of regions.

The dx and dy parameters specify the size of gaps that are still considered to be the same object.
"""
function components(region::Region, dx::Unsigned=unsigned(1), dy::Unsigned=unsigned(1))

    function uf_find_root(ufa, x)
        i = x
        while ufa[i] >= 0
            i = ufa[i]
        end
        return i
    end

    function uf_union(ufa, x, y)
        i = uf_find_root(ufa, x)
        j = uf_find_root(ufa, y)

        ti = x
        tmp = -1
        while ufa[ti] >= 0
            tmp = ti
            ti = ufa[ti]
            ufa[tmp] = i
        end

        ti = y
        tmp = -1
        while ufa[ti] >= 0
            tmp = ti
            ti = ufa[ti]
            ufa[tmp] = j
        end

        not_same_tree = i != j

        if not_same_tree
            ## this maintains the scan-order of the regions in the result
            if i < j
                ufa[j] = i
            else
                ufa[i] = j
            end
        end

        return not_same_tree

    end

    ## union_find: indexed by chord index in r.chords()
    ##
    ## value:
    ## <= -1:  this node is a root with abs(val) indicating tree depth of the associated tree
    ## >= 0 : index index (in union_find) of parent

    union_find = [-1 for _ = 1:length(region.runs)]

    run_index = 1

    ## see ngi::chord::are_chords_close; this is required for 4-connectivity
    dy = max(1, dy)
    for run in region.runs
        next_run_index = run_index + 1
        while next_run_index <= length(region.runs) &&
            region.runs[next_run_index].row <= run.row + dy
            if isclose(run, region.runs[next_run_index], dx, dy)
                uf_union(union_find, run_index, next_run_index)
            end
            next_run_index += 1
        end
        run_index += 1
    end

    ## Make contiguous labels for the roots: store a negative value in the union_find array
    ## at the root location:
    ##    -1 <-> region 0
    ##    -2 <-> region 1
    ## ....
    ## Also counts the # of roots (i.e. number of connected components).
    r_idx = -1

    for idx = 1:length(union_find)
        if union_find[idx] <= -1
            union_find[idx] = r_idx
            r_idx -= 1
        end
    end

    num_regions = r_idx * (-1) - 1;
    connected_objects = Region[Region(Run[]) for _ = 1:num_regions];

    for idx = 1:length(union_find)
        root = uf_find_root(union_find, idx)
        region_idx = union_find[root] * (-1)
        push!(connected_objects[region_idx].runs, region.runs[idx])
    end

    return connected_objects;
end

end # module


text:
"""
    μ = center!(X[, μ; obsdim, operate_on])

or

    μ = center!(D[, μ; operate_on])

where `X` is of type Matrix or Vector and `D` of type DataFrame.

Shift `X` along `obsdim` by `μ` according to X = X - μ

`μ`         :  Vector or value describing the translation.
               Defaults to mean(X, 2)

`obsdim`    :  Specify which axis corresponds to observations.
               Defaults to obsdim=2 (observations are columns of matrix)
               For DataFrames `obsdim` is obsolete and centering occurs
               column wise.

`operate_on`:  Specify the indices of columns or rows to be centered.
               Defaults to all columns/rows.
               For DataFrames this must be a vector of symbols, not indices
               E.g. `operate_on`=[1,3] will perform centering on columns
               with index 1 and 3 only (if obsdim=1, else rows 1 and 3)

Note on DataFrames:
Columns containing `NA` values are skipped.
Columns containing non numeric elements are skipped.

Examples:

    X = rand(4, 100)
    x = rand(10)
    D = DataFrame(A=rand(10), B=collect(1:10), C=[string(x) for x in 1:10])

    μ = center!(X, obsdim=2)
    μ = center!(X, ObsDim.First())
    μ = center!(X, obsdim=1, operate_on=[1,3]
    μ = center!(X, [7.0, 8.0], obsdim=1, operate_on=[1,3]
    μ = center!(D)
    μ = center!(D, operate_on=[:A, :B])
    μ = center!(D, [-1,-1], operate_on=[:A, :B])
"""

function center!(X, μ; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim)))
    center!(X, μ, convert(ObsDimension, obsdim), operate_on)
end

function center!(X; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim)))
    center!(X, convert(ObsDimension, obsdim), operate_on)
end

function center!(X::AbstractArray{T,N}, obsdim::ObsDim.Constant{M}; operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T,N,M}
    center!(X, ObsDim.Constant{M}(), operate_on)
end

function center!(X::AbstractArray{T,N}, obsdim::ObsDim.Last; operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T,N}
    center!(X, ObsDim.Constant{N}(), operate_on)
end

function center!(X::AbstractArray{T,N}, obsdim::ObsDim.Last, operate_on::AbstractVector) where {T,N}
    center!(X, ObsDim.Constant{N}(), operate_on)
end

function center!(X::AbstractArray{T,N}, obsdim::ObsDim.Constant{M}, operate_on::AbstractVector) where {T,N,M}
    μ = vec(mean(X, M))[operate_on]
    center!(X, μ, obsdim, operate_on)
end

function center!(X::AbstractArray{T,N}, μ::AbstractVector, obsdim::ObsDim.Constant{M}; operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T,N,M}
    center!(X, μ, ObsDim.Constant{M}(), operate_on)
end

function center!(X::AbstractArray{T,N}, μ::AbstractVector, obsdim::ObsDim.Last; operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T,N}
    center!(X, μ, ObsDim.Constant{N}(), operate_on)
end

function center!(X::AbstractMatrix, μ::AbstractVector, ::ObsDim.Constant{1}, operate_on)
    @assert length(μ) == length(operate_on)
    nObs, nVars = size(X)
    for (i, iVar) in enumerate(operate_on)
        @inbounds for iObs in 1:nObs
            X[iObs, iVar] = X[iObs, iVar] - μ[i]
        end
    end
    μ
end

function center!(X::AbstractMatrix, μ::AbstractVector, ::ObsDim.Constant{2}, operate_on)
    @assert length(μ) == length(operate_on)
    nVars, nObs = size(X)
    for iObs in 1:nObs
        @inbounds for (i, iVar) in enumerate(operate_on)
            X[iVar, iObs] = X[iVar, iObs] - μ[i]
        end
    end
    μ
end

function center!(x::AbstractVector; obsdim=LearnBase.default_obsdim(x), operate_on=default_scaleselection(x))
    center!(x, convert(ObsDimension, obsdim), operate_on)
end

function center!(x::AbstractVector, ::ObsDim.Constant, operate_on::AbstractVector)
    μ = mean(x)
    for iVar in operate_on
        x[iVar] = x[iVar] - μ
    end
    μ
end

function center!(x::AbstractVector, μ::AbstractVector, ::ObsDim.Constant{1}, operate_on::AbstractVector)
    @assert length(μ) == length(operate_on)
    @inbounds for (i, iVar) in enumerate(operate_on)
        x[iVar] = x[iVar] - μ[i]
    end
    μ
end

function center!(x::AbstractVector, μ::AbstractVector, ::ObsDim.Last, operate_on::AbstractVector)
    center!(x, μ, ObsDim.Constant{1}(), operate_on)
end

function center!(x::AbstractVector, μ::Real, ::ObsDim.Constant{1}, operate_on)
    @inbounds for i in operate_on
        x[i] = x[i] - μ
    end
    μ
end

function center!(x::AbstractVector, μ::Real, ::ObsDim.Last, operate_on)
    center!(x, μ, ObsDim.Constant{1}(), operate_on)
end

# --------------------------------------------------------------------

function center!(D::AbstractDataFrame; operate_on=default_scaleselection(D))
    center!(D, operate_on)
end

function center!(D::AbstractDataFrame, operate_on::AbstractVector{Symbol})
    μ_vec = Float64[]
    for colname in operate_on
        if eltype(D[colname]) <: Real
            μ = mean(D[colname])
            if isna(μ)
                warn("Skipping \"$colname\" because it contains NA values")
                continue
            end
            center!(D, μ, colname)
            push!(μ_vec, μ)
        else
            warn("Skipping \"$colname\" because data is not of type T <: Real.")
        end
    end
    μ_vec
end

function center!(D::AbstractDataFrame, μ::AbstractVector; operate_on=default_scaleselection(D))
    center!(D, μ, operate_on)
end

function center!(D::AbstractDataFrame, μ::AbstractVector, operate_on::AbstractVector{Symbol})
    for (icol, colname) in enumerate(operate_on)
        if eltype(D[colname]) <: Real
            center!(D, μ[icol], colname)
        else
            warn("Skipping \"$colname\" because data is not of type T <: Real.")
        end
    end
    μ
end

function center!(D::AbstractDataFrame, μ::Real, colname::Symbol)
    if sum([isna(value) for value in D[colname]]) > 0
        warn("Skipping \"$colname\" because it contains NA values")
    else
        newcol::Vector{Float64} = convert(Vector{Float64}, D[colname])
        center!(newcol, μ)
        D[colname] = newcol
    end
    μ
end


text:
module CmeSIM

export Track, SIMTrack, FluorophoreMarker, Fluorophore, Acquisition, SIMAcquisition, times

include("./marker/mod.jl")
include("./acquisition/mod.jl")
include("./tracks/mod.jl")

end


text:
export elements, benchmark_elements

function benchmark_elements(times::Int)
    r = [8.59072560e+02, -4.13720368e+03, 5.29556871e+03]
    v = [7.37289205e+00, 2.08223573e+00, 4.39999794e-01]
    mu = 3.986004418e5
    worst = -Inf
    best = Inf
    total = 0.0
    elements(r, v, mu)
    for i = 1:times
        gc_enable(false)
        t = @elapsed elements(r, v, mu)
        gc_enable(true)
        if t > worst
            worst = t
        end
        if t < best
            best = t
        end
        total += t
    end
    println("[$(total/times),$best,$worst]")
end

function elements(r::AbstractArray, v::AbstractArray, mu::Float64)
    rm = norm(r)
    vm = norm(v)
    h = cross(r, v)
    hm = norm(h)
    k = [0.0, 0.0, 1.0]
    n = cross(k, h)
    nm = norm(n)
    xi = vm^2/2 - mu/rm
    ec = ((vm^2 - mu/rm)*r - v*dot(r, v))/mu
    ecc = norm(ec)
    if ecc != 1
        sma = -mu/(2*xi)
    else
        sma = hm^2/mu
    end
    inc = acos(h[3]/hm)
    node = acos(n[1]/nm)
    peri = acos(dot(n, ec)/(ecc*nm))
    ano = acos(dot(ec, r)/(ecc*rm))
    if n[2] < 0
        node = 2*pi - node
    end
    if ec[3] < 0
        peri = 2*pi - peri
    end
    if dot(r, v) < 0
        ano = 2*pi - ano
    end
    return sma, ecc, inc, node, peri, ano
end

max_issues_repo_name: OpenAstrodynamics/benchmarks

max_forks_repo_name: OpenAstrodynamics/benchmarks


text:
module Day6
export day6

function readInput(inputfile="inputs/input06.txt")
    data = Any[]

    current = Any[]
    for line in eachline(inputfile)
        if line == "" && length(current) > 0
            push!(data, current)
            current = Any[]
        else
            push!(current, Set(line))
        end
    end
    if length(current) > 0
        push!(data, current)
    end
    data
end

function day6(inputfile="inputs/input06.txt")
    values = readInput(inputfile)
    (reduce(+, map(v -> length(union(v...)), values)),
     reduce(+, map(v -> length(intersect(v...)), values)))
end

end #module


text:
# ---------------------------------------------------------------------------- #
#
#   types_and_consts_egads.jl
#
#   Definition of types and constants in Julia from egads. This file is
#     automatically generated using the Clang.jl C-wrapper.
#
#   egads.jl
#   Spring 2018
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

# Automatically generated using Clang.jl wrap_c, version 0.0.0

using Compat

const EGADS_READERR = Int32(-32)
const EGADS_TESSTATE = Int32(-31)
const EGADS_EXISTS = Int32(-30)
const EGADS_ATTRERR = Int32(-29)
const EGADS_TOPOCNT = Int32(-28)
const EGADS_OCSEGFLT = Int32(-27)
const EGADS_BADSCALE = Int32(-26)
const EGADS_NOTORTHO = Int32(-25)
const EGADS_DEGEN = Int32(-24)
const EGADS_CONSTERR = Int32(-23)
const EGADS_TOPOERR = Int32(-22)
const EGADS_GEOMERR = Int32(-21)
const EGADS_NOTBODY = Int32(-20)
const EGADS_WRITERR = Int32(-19)
const EGADS_NOTMODEL = Int32(-18)
const EGADS_NOLOAD = Int32(-17)
const EGADS_RANGERR = Int32(-16)
const EGADS_NOTGEOM = Int32(-15)
const EGADS_NOTTESS = Int32(-14)
const EGADS_EMPTY = Int32(-13)
const EGADS_NOTTOPO = Int32(-12)
const EGADS_REFERCE = Int32(-11)
const EGADS_NOTXFORM = Int32(-10)
const EGADS_NOTCNTX = Int32(-9)
const EGADS_MIXCNTX = Int32(-8)
const EGADS_NODATA = Int32(-7)
const EGADS_NONAME = Int32(-6)
const EGADS_INDEXERR = Int32(-5)
const EGADS_MALLOC = Int32(-4)
const EGADS_NOTOBJ = Int32(-3)
const EGADS_NULLOBJ = Int32(-2)
const EGADS_NOTFOUND = Int32(-1)
const EGADS_SUCCESS = Int32(0)
const EGADS_OUTSIDE = Int32(1)
const EGADSMAJOR = Int32(1)
const EGADSMINOR = Int32(12)

# Skipping MacroDefinition: EGADSPROP EGADSprop : Revision 1.12 #

const MAGIC = Int32(98789)
const MTESSPARAM = Int32(2)
const CONTXT = Int32(0)
const TRANSFORM = Int32(1)
const TESSELLATION = Int32(2)
const NIL = Int32(3)
const EMPTY = Int32(4)
const REFERENCE = Int32(5)
const PCURVE = Int32(10)
const CURVE = Int32(11)
const SURFACE = Int32(12)
const NODE = Int32(20)
const EDGE = Int32(21)
const LOOP = Int32(22)
const FACE = Int32(23)
const SHELL = Int32(24)
const BODY = Int32(25)
const MODEL = Int32(26)
const LINE = Int32(1)
const CIRCLE = Int32(2)
const ELLIPSE = Int32(3)
const PARABOLA = Int32(4)
const HYPERBOLA = Int32(5)
const TRIMMED = Int32(6)
const BEZIER = Int32(7)
const BSPLINE = Int32(8)
const OFFSET = Int32(9)
const PLANE = Int32(1)
const SPHERICAL = Int32(2)
const CYLINDRICAL = Int32(3)
const REVOLUTION = Int32(4)
const TOROIDAL = Int32(5)
const CONICAL = Int32(10)
const EXTRUSION = Int32(11)
const SREVERSE = Int32(-1)
const NOMTYPE = Int32(0)
const SFORWARD = Int32(1)
const ONENODE = Int32(1)
const TWONODE = Int32(2)
const OPEN = Int32(3)
const CLOSED = Int32(4)
const DEGENERATE = Int32(5)
const WIREBODY = Int32(6)
const FACEBODY = Int32(7)
const SHEETBODY = Int32(8)
const SOLIDBODY = Int32(9)
const ATTRINT = Int32(1)
const ATTRREAL = Int32(2)
const ATTRSTRING = Int32(3)
const ATTRCSYS = Int32(12)
const ATTRPTR = Int32(13)
const SUBTRACTION = Int32(1)
const INTERSECTION = Int32(2)
const FUSION = Int32(3)
const BOX = Int32(1)
const SPHERE = Int32(2)
const CONE = Int32(3)
const CYLINDER = Int32(4)
const TORUS = Int32(5)
const UISO = Int32(0)
const VISO = Int32(1)
const NODEOFF = Int32(1)
const EDGEOFF = Int32(2)
const FACEDUP = Int32(3)
const FACECUT = Int32(4)
const FACEOFF = Int32(5)

mutable struct egAttr
    name::Cstring
    _type::Cint
    length::Cint
    vals::Void
end

mutable struct egAttrs
    nattrs::Cint
    attrs::Ptr{egAttr}
end

mutable struct egObject
    magicnumber::Cint
    oclass::Int16
    mtype::Int16
    attrs::Ptr{Void}
    blind::Ptr{Void}
    topObj::Ptr{egObject}
    tref::Ptr{egObject}
    prev::Ptr{egObject}
    next::Ptr{egObject}
end

const ego = Ptr{egObject}

const NULL_E = ego(0)

mutable struct egCntxt
    outLevel::Cint
    tess::NTuple{2, Cdouble}
    signature::Ptr{Cstring}
    usrPtr::Ptr{Void}
    pool::Ptr{egObject}
    last::Ptr{egObject}
end

mutable struct egFconn
    index::Cint
    nface::Cint
    faces::Ptr{Cint}
    tric::Ptr{Cint}
end

mutable struct egBary
    tri::Cint
    w::NTuple{2, Cdouble}
end

mutable struct egTess1D
    obj::Ptr{egObject}
    nodes::NTuple{2, Cint}
    faces::NTuple{2, egFconn}
    xyz::Ptr{Cdouble}
    t::Ptr{Cdouble}
    _global::Ptr{Cint}
    npts::Cint
end

mutable struct egPatch
    ipts::Ptr{Cint}
    bounds::Ptr{Cint}
    nu::Cint
    nv::Cint
end

mutable struct egTess2D
    mKnots::Ptr{egObject}
    xyz::Ptr{Cdouble}
    uv::Ptr{Cdouble}
    _global::Ptr{Cint}
    ptype::Ptr{Cint}
    pindex::Ptr{Cint}
    bary::Ptr{egBary}
    frame::Ptr{Cint}
    frlps::Ptr{Cint}
    tris::Ptr{Cint}
    tric::Ptr{Cint}
    patch::Ptr{egPatch}
    npts::Cint
    nframe::Cint
    nfrlps::Cint
    ntris::Cint
    npatch::Cint
    tfi::Cint
end

mutable struct egTessel
    src::Ptr{egObject}
    xyzs::Ptr{Cdouble}
    tess1d::Ptr{egTess1D}
    tess2d::Ptr{egTess2D}
    globals::Ptr{Cint}
    params::NTuple{6, Cdouble}
    tparam::NTuple{2, Cdouble}
    nGlobal::Cint
    nEdge::Cint
    nFace::Cint
    nu::Cint
    nv::Cint
    done::Cint
end

# Skipping MacroDefinition: CHECK_STATUS ( X ) if ( status < SUCCESS ) { printf ( "ERROR:: BAD STATUS = %d from %s (called from %s:%d)\n" , status , # X , routine , __LINE__ ) ; goto cleanup ; } #
# Skipping MacroDefinition: SET_STATUS ( STAT , X ) status = STAT ; printf ( "ERROR:: BAD STATUS = %d from %s (called from %s:%d)\n" , status , # X , routine , __LINE__ ) ; goto cleanup ; #
# Skipping MacroDefinition: MALLOC ( PTR , TYPE , SIZE ) DPRINT3 ( "mallocing %s in routine %s (size=%d)" , # PTR , routine , SIZE ) ; if ( PTR != NULL ) { printf ( "ERROR:: MALLOC overwrites for %s=%llx (called from %s:%d)\n" , # PTR , ( long long ) PTR , routine , __LINE__ ) ; status = BAD_MALLOC ; goto cleanup ; } PTR = ( TYPE * ) malloc ( ( SIZE ) * sizeof ( TYPE ) ) ; if ( PTR == NULL ) { printf ( "ERROR:: MALLOC PROBLEM for %s (called from %s:%d)\n" , # PTR , routine , __LINE__ ) ; status = BAD_MALLOC ; goto cleanup ; } #
# Skipping MacroDefinition: RALLOC ( PTR , TYPE , SIZE ) DPRINT3 ( "rallocing %s in routine %s (size=%d)" , # PTR , routine , SIZE ) ; if ( PTR == NULL ) { MALLOC ( PTR , TYPE , SIZE ) ; } else { realloc_temp = realloc ( PTR , ( SIZE ) * sizeof ( TYPE ) ) ; if ( PTR == NULL ) { printf ( "ERROR:: RALLOC PROBLEM for %s (called from %s:%d)\n" , # PTR , routine , __LINE__ ) ; status = BAD_MALLOC ; goto cleanup ; } else { PTR = ( TYPE * ) realloc_temp ; } } #
# Skipping MacroDefinition: FREE ( PTR ) if ( PTR != NULL ) { DPRINT2 ( "freeing %s in routine %s" , # PTR , routine ) ; free ( PTR ) ; } PTR = NULL ; /* macros for debug printing */
# Skipping MacroDefinition: DPRINT0 ( FORMAT ) #
# Skipping MacroDefinition: DPRINT0x ( FORMAT ) #
# Skipping MacroDefinition: DPRINT1 ( FORMAT , A ) #
# Skipping MacroDefinition: DPRINT1x ( FORMAT , A ) #
# Skipping MacroDefinition: DPRINT2 ( FORMAT , A , B ) #
# Skipping MacroDefinition: DPRINT2x ( FORMAT , A , B ) #
# Skipping MacroDefinition: DPRINT3 ( FORMAT , A , B , C ) #
# Skipping MacroDefinition: DPRINT3x ( FORMAT , A , B , C ) #
# Skipping MacroDefinition: DPRINT4 ( FORMAT , A , B , C , D ) #
# Skipping MacroDefinition: DPRINT4x ( FORMAT , A , B , C , D ) #
# Skipping MacroDefinition: DPRINT5 ( FORMAT , A , B , C , D , E ) #
# Skipping MacroDefinition: DPRINT6 ( FORMAT , A , B , C , D , E , F ) #
# Skipping MacroDefinition: DPRINT7 ( FORMAT , A , B , C , D , E , F , G ) #
# Skipping MacroDefinition: DPRINT7x ( FORMAT , A , B , C , D , E , F , G ) #
# Skipping MacroDefinition: DPRINT8 ( FORMAT , A , B , C , D , E , F , G , H ) #
# Skipping MacroDefinition: DPRINT9 ( FORMAT , A , B , C , D , E , F , G , H , I ) #
# Skipping MacroDefinition: DPRINT10 ( FORMAT , A , B , C , D , E , F , G , H , I , J ) #
# Skipping MacroDefinition: DPRINT11 ( FORMAT , A , B , C , D , E , F , G , H , I , J , K ) #
# Skipping MacroDefinition: DPRINT12 ( FORMAT , A , B , C , D , E , F , G , H , I , J , K , L ) #
# Skipping MacroDefinition: DPRINT13 ( FORMAT , A , B , C , D , E , F , G , H , I , J , K , L , M ) #
# Skipping MacroDefinition: DPRINT14 ( FORMAT , A , B , C , D , E , F , G , H , I , J , K , L , M , N ) #
# Skipping MacroDefinition: DPRINT15 ( FORMAT , A , B , C , D , E , F , G , H , I , J , K , L , M , N , O ) #
# Skipping MacroDefinition: ROUTINE ( NAME ) char routine [ ] = # NAME ; if ( routine [ 0 ] == '\0' ) printf ( "bad routine(%s)\n" , routine ) ; /* macros for status printing */
# Skipping MacroDefinition: SPRINT0 ( OL , FORMAT ) DPRINT0 ( FORMAT ) ; if ( outLevel >= OL ) { printf ( FORMAT ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT0x ( OL , FORMAT ) DPRINT0x ( FORMAT ) ; if ( outLevel >= OL ) { printf ( FORMAT ) ; fflush ( stdout ) ; } #
# Skipping MacroDefinition: SPRINT1 ( OL , FORMAT , A ) DPRINT1 ( FORMAT , A ) ; if ( outLevel >= OL ) { printf ( FORMAT , A ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT1x ( OL , FORMAT , A ) DPRINT1x ( FORMAT , A ) ; if ( outLevel >= OL ) { printf ( FORMAT , A ) ; fflush ( stdout ) ; } #
# Skipping MacroDefinition: SPRINT2 ( OL , FORMAT , A , B ) DPRINT2 ( FORMAT , A , B ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT2x ( OL , FORMAT , A , B ) DPRINT2x ( FORMAT , A , B ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B ) ; fflush ( stdout ) ; } #
# Skipping MacroDefinition: SPRINT3 ( OL , FORMAT , A , B , C ) DPRINT3 ( FORMAT , A , B , C ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT3x ( OL , FORMAT , A , B , C ) DPRINT3x ( FORMAT , A , B , C ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C ) ; fflush ( stdout ) ; } #
# Skipping MacroDefinition: SPRINT4 ( OL , FORMAT , A , B , C , D ) DPRINT4 ( FORMAT , A , B , C , D ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C , D ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT4x ( OL , FORMAT , A , B , C , D ) DPRINT4x ( FORMAT , A , B , C , D ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C , D ) ; fflush ( stdout ) ; } #
# Skipping MacroDefinition: SPRINT5 ( OL , FORMAT , A , B , C , D , E ) DPRINT5 ( FORMAT , A , B , C , D , E ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C , D , E ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT6 ( OL , FORMAT , A , B , C , D , E , F ) DPRINT6 ( FORMAT , A , B , C , D , E , F ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C , D , E , F ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT7 ( OL , FORMAT , A , B , C , D , E , F , G ) DPRINT7 ( FORMAT , A , B , C , D , E , F , G ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C , D , E , F , G ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT8 ( OL , FORMAT , A , B , C , D , E , F , G , H ) DPRINT8 ( FORMAT , A , B , C , D , E , F , G , H ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C , D , E , F , G , H ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT9 ( OL , FORMAT , A , B , C , D , E , F , G , H , I ) DPRINT9 ( FORMAT , A , B , C , D , E , F , G , H , I ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C , D , E , F , G , H , I ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT10 ( OL , FORMAT , A , B , C , D , E , F , G , H , I , J ) DPRINT10 ( FORMAT , A , B , C , D , E , F , G , H , I , J ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C , D , E , F , G , H , I , J ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT11 ( OL , FORMAT , A , B , C , D , E , F , G , H , I , J , K ) DPRINT11 ( FORMAT , A , B , C , D , E , F , G , H , I , J , K ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C , D , E , F , G , H , I , J , K ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT12 ( OL , FORMAT , A , B , C , D , E , F , G , H , I , J , K , L ) DPRINT12 ( FORMAT , A , B , C , D , E , F , G , H , I , J , K , L ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C , D , E , F , G , H , I , J , K , L ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT13 ( OL , FORMAT , A , B , C , D , E , F , G , H , I , J , K , L , M ) DPRINT13 ( FORMAT , A , B , C , D , E , F , G , H , I , J , K , L , M ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C , D , E , F , G , H , I , J , K , L , M ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT14 ( OL , FORMAT , A , B , C , D , E , F , G , H , I , J , K , L , M , N ) DPRINT14 ( FORMAT , A , B , C , D , E , F , G , H , I , J , K , L , M , N ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C , D , E , F , G , H , I , J , K , L , M , N ) ; printf ( "\n" ) ; } #
# Skipping MacroDefinition: SPRINT15 ( OL , FORMAT , A , B , C , D , E , F , G , H , I , J , K , L , M , N , O ) DPRINT15 ( FORMAT , A , B , C , D , E , F , G , H , I , J , K , L , M , N , O ) ; if ( outLevel >= OL ) { printf ( FORMAT , A , B , C , D , E , F , G , H , I , J , K , L , M , N , O ) ; printf ( "\n" ) ; } /* error codes */

const BAD_MALLOC = -900
const PI = 3.141592653589793
const TWOPI = 6.283185307179586
const PIo2 = 1.5707963267948966
const PIo4 = 0.7853981633974483
const PIo180 = 0.017453292519943295
const HUGEQ = 9.9999999e7
const HUGEI = 9999999
const EPS03 = 0.001
const EPS06 = 1.0e-6
const EPS12 = 1.0e-12
const EPS20 = 1.0e-20

# Skipping MacroDefinition: SQR ( A ) ( ( A ) * ( A ) ) #
# Skipping MacroDefinition: NINT ( A ) ( ( ( A ) < 0 ) ? ( int ) ( A - 0.5 ) : ( int ) ( A + 0.5 ) ) #
# Skipping MacroDefinition: MIN ( A , B ) ( ( ( A ) < ( B ) ) ? ( A ) : ( B ) ) #
# Skipping MacroDefinition: MAX ( A , B ) ( ( ( A ) < ( B ) ) ? ( B ) : ( A ) ) #
# Skipping MacroDefinition: MINMAX ( A , B , C ) MIN ( MAX ( ( A ) , ( B ) ) , ( C ) ) #
# Skipping MacroDefinition: SIGN ( A ) ( ( ( A ) < 0 ) ? - 1 : ( ( ( A ) > 0 ) ? + 1 : 0 ) ) #
# Skipping MacroDefinition: FSIGN ( A , B ) ( ( B ) >= 0 ? fabs ( A ) : - fabs ( A ) ) #

const CHUNK = 256
const MAXELEN = 2048
const DEGENUV = 1.0e-14

mutable struct KEY
    keys::NTuple{3, Cint}
end

mutable struct DATA
    close::Cint
    xyz::NTuple{3, Cdouble}
end

mutable struct ENTRY
    key::KEY
    data::DATA
end

mutable struct element
    item::ENTRY
    next::Ptr{element}
end

mutable struct ELEMENT
    item::ENTRY
    next::Ptr{element}
end

mutable struct triVert
    _type::Cint
    edge::Cint
    index::Cint
    xyz::NTuple{3, Cdouble}
    uv::NTuple{2, Cdouble}
end

mutable struct triTri
    indices::NTuple{3, Cint}
    neighbors::NTuple{3, Cint}
    mid::NTuple{3, Cdouble}
    area::Cdouble
    mark::Int16
    close::Int16
    hit::Int16
    count::Int16
end

mutable struct triSeg
    indices::NTuple{2, Cint}
    neighbor::Cint
    edge::Cint
    index::Cint
end

mutable struct triStruct
    face::Ptr{Cint}
    fIndex::Cint
    orUV::Cint
    planar::Cint
    phase::Cint
    VoverU::Cdouble
    maxlen::Cdouble
    chord::Cdouble
    dotnrm::Cdouble
    accum::Cdouble
    edist2::Cdouble
    eps2::Cdouble
    devia2::Cdouble
    minlen::Cdouble
    range::NTuple{4, Cdouble}
    qparm::NTuple{3, Cdouble}
    uvs::Ptr{Cdouble}
    orCnt::Cint
    maxPts::Cint
    mverts::Cint
    nverts::Cint
    verts::Ptr{triVert}
    mtris::Cint
    ntris::Cint
    tris::Ptr{triTri}
    msegs::Cint
    nsegs::Cint
    segs::Ptr{triSeg}
    nfrvrts::Cint
    mframe::Cint
    nframe::Cint
    frame::Ptr{Cint}
    mloop::Cint
    nloop::Cint
    loop::Ptr{Cint}
    lens::NTuple{5, Cint}
    numElem::Cint
    tfi::Cint
    hashTab::Ptr{Ptr{ELEMENT}}
end

mutable struct Front
    sleft::Cint
    i0::Cint
    i1::Cint
    sright::Cint
    snew::Int16
    mark::Int16
end

mutable struct fillArea
    mfront::Cint
    nfront::Cint
    mpts::Cint
    npts::Cint
    pts::Ptr{Cint}
    nsegs::Cint
    segs::Ptr{Cint}
    front::Ptr{Front}
end

mutable struct connect
    node1::Cint
    node2::Cint
    tri::Ptr{Cint}
    thread::Cint
end

mutable struct EMPtess
    mutex::Ptr{Void}
    master::Clong
    _end::Cint
    index::Cint
    ignore::Cint
    mark::Ptr{Cint}
    tess::Ptr{Cint}
    btess::Ptr{Cint}
    body::Ptr{Cint}
    faces::Ptr{Ptr{Cint}}
    edges::Ptr{Ptr{Cint}}
    params::Ptr{Cdouble}
    tparam::Ptr{Cdouble}
    qparam::NTuple{3, Cdouble}
    ptr::Ptr{Void}
end

const SUCCESS = 0
const MALLOC_ERROR = -901
const BAD_POINT_INDEX = -902
const BAD_TRIANGLE_INDEX = -903
const CANNOT_SWAP = -904
const COULD_NOT_RECOVER_BND = -905
const NUMBER_OF_POINT_MISMATCH = -906
const PRM_OK_ONEFACE = 1
const PRM_OK_SIMPLE = 2
const PRM_OK_JONES = 3
const PRM_OK_AXIAL = 4
const PRM_OK_PLANAR = 5
const PRM_OK_FLOATER = 6
const PRM_OK_UNROLLING = 7
const PRM_NOTCONVERGED = -401
const PRM_BADNUMVERTICES = -402
const PRM_ZEROPIVOT = -403
const PRM_NEGATIVEAREAS = -404
const PRM_TOLERANCEUNMET = -405
const PRM_NOGLOBALUV = -406
const PRM_BADPARAM = -407
const PRM_BADDIVISION = -408
const PRM_CANNOTFORMLOOP = -409
const PRM_WIGGLEDETECTED = -410
const PRM_INTERNAL = -499

mutable struct prmXYZ
    x::Cdouble
    y::Cdouble
    z::Cdouble
end

mutable struct prmUVF
    u0::Cdouble
    v0::Cdouble
    u1::Cdouble
    v1::Cdouble
    u2::Cdouble
    v2::Cdouble
end

mutable struct prmUV
    u::Cdouble
    v::Cdouble
end

mutable struct prmTri
    indices::NTuple{3, Cint}
    neigh::NTuple{3, Cint}
    own::Cint
end

mutable struct cfitCell
    wKnot::Cint
    eKnot::Cint
    wCell::Cint
    eCell::Cint
    child::Cint
    count::Cint
    umin::Cdouble
    umax::Cdouble
    rmserr::Cdouble
    maxerr::Cdouble
end

mutable struct cfitTree
    nvar::Cint
    nu::Cint
    ncel::Cint
    cell::Ptr{cfitCell}
    nknt::Cint
    knot::Ptr{Cdouble}
end

mutable struct gridCell
    swKnot::Cint
    seKnot::Cint
    nwKnot::Cint
    neKnot::Cint
    wCell::Cint
    eCell::Cint
    sCell::Cint
    nCell::Cint
    dtype::Cint
    child::Cint
    count::Cint
    umin::Cdouble
    umax::Cdouble
    vmin::Cdouble
    vmax::Cdouble
    rmserr::Cdouble
    maxerr::Cdouble
end

mutable struct gridTree
    nvar::Cint
    nu::Cint
    nv::Cint
    ncel::Cint
    cell::Ptr{gridCell}
    nknt::Cint
    knot::Ptr{Cdouble}
end


text:
using AutoPreallocation
using BenchmarkTools
using Test

# functions to test on
f_ones() = ones(64)
f_matmul() = ones(32,64) * ones(64, 2)

f_matmul_noprealloc() = ones(32,64) * @no_prealloc(ones(64, 2))

@testset "ones example" begin
    @assert (@ballocated f_ones()) === 624
    val, record = record_allocations(f_ones)
    @test (@ballocated avoid_allocations($record, f_ones)) <= 64
end

@testset "matmul example" begin
    @assert (@ballocated f_matmul()) === 18_304
    val, record = record_allocations(f_matmul)
    # NOTE: (@Roger-luo) not sure why this is 256 on my machine
    @test (@ballocated avoid_allocations($record, f_matmul)) <= 352

    _, pf = preallocate(f_matmul)
    @test (@ballocated $pf()) <= 352
end

@testset "noprealloc example" begin
    @assert (@ballocated f_matmul_noprealloc()) === 18_304
    val, record = record_allocations(f_matmul_noprealloc)
    @test length(record.allocations) == 2
    @test record.initial_sizes == [(32, 64), (32, 2)]
    @test (@ballocated avoid_allocations($record, f_matmul_noprealloc)) <= 1520
end

@testset "check thread-safe" begin
    f(x, y) = x * y
    n = Threads.nthreads()
    results = Vector{Any}(undef, n)
    As = [rand(4, 4) for _ in 1:n]
    Bs = [rand(4, 4) for _ in 1:n]
    _, pf = preallocate(f, As[1], Bs[1])
    Threads.@threads for k in 1:n
        results[k] = pf(As[k], Bs[k])
    end

    # if it's not thread-safe, the result
    # will be modified by other multiplications
    # since they will share the same memory
    for k in 1:n
        @test results[k] ≈ f(As[k], Bs[k])
    end
end


text:
using Documenter
using Pkg
Pkg.activate("..")
#push!(LOAD_PATH,"../src/")
using SymmetricPolynomials

makedocs(sitename="SymmetricPolynomials")

max_stars_repo_name: yushoteke/SymmetricPolynomials

max_issues_repo_name: yushoteke/SymmetricPolynomials

max_forks_repo_name: yushoteke/SymmetricPolynomials


text:
import Base.show

"""

# LooCompare

A struct containing the results of PsisLoo model comparisom.

$(FIELDS)

# Extended help

# Fields
  - `psis::Vector{PsisLoo}` : Vector of PsisLoo objects.
  - `table::KeyedArray` : Comparison table.

# Example of a comparison table
```
┌───────┬───────────┬─────────┬────────┐
│       │ elpd_diff │ se_diff │ weight │
├───────┼───────────┼─────────┼────────┤
│ m5_1t │      0.00 │    0.00 │   0.67 │
│ m5_3t │     -0.69 │    0.42 │   0.33 │
│ m5_2t │     -6.68 │    4.74 │   0.00 │
└───────┴───────────┴─────────┴────────┘
```

where:

1. `elpd_diff`  : Difference between total loo_est values between models.
2. `se_diff`    : Standard error of the difference in total l00_est values.
3. `weight`     : Relative support for each model.

In this example table the models have been sorted in ascending total loo_est values.
The PsisLoo objects in the field `psis` is sorted as listed in `table`.

See also: [`PsisLoo`](@ref).
"""
struct LooCompare
    psis::Vector{PsisLoo}
    table::KeyedArray
end

"""

# loo_compare

Construct a PsisLoo comparison table for loglikelihood matrices.

$(SIGNATURES)

# Extended help

### Required arguments
    - `loglikelihoods::Vector{Array{Float64, 3}}` : Vector of loglikelihood matrices

### Optional arguments
    - `model_names=nothing` : Optional specify models
    - `sort_models=true` : Sort models according to ascending elpd values

### Return values
    - `result::LooCompare` : LooCompare object

See also: [`LooCompare`](@ref).
"""
function loo_compare(
    loglikelihoods::Vector{Array{Float64, 3}};
    model_names=nothing,
    sort_models=true)

    nmodels = length(loglikelihoods)

    if isnothing(model_names)
        mnames = ["model_$i" for i in 1:nmodels]
    else
        mnames = model_names
    end

    psis = psis_loo.(loglikelihoods)

    psis_values = Vector{Float64}(undef, nmodels)
    se_values = Vector{Float64}(undef, nmodels)
    loos = Vector{Vector{Float64}}(undef, nmodels)

    for i in 1:nmodels
        psis_values[i] = psis[i].estimates(:loo_est, :total)
        se_values[i] = psis[i].estimates(:loo_est, :se_total)
        loos[i] = psis[i].pointwise(:loo_est)
    end

    if sort_models
        ind = sortperm([psis_values[i][1] for i in 1:nmodels]; rev=true)
        psis = psis[ind]
        psis_values = psis_values[ind]
        se_values = se_values[ind]
        loos = loos[ind]
        mnames = mnames[ind]
    end

    # Setup comparison vectors

    elpd_diff = zeros(nmodels)
    se_diff = zeros(nmodels)
    weight = ones(nmodels)

    # Compute comparison values

    for i in 2:nmodels
        elpd_diff[i] = psis_values[i] - psis_values[1]
        diff = loos[1] - loos[i]
        se_diff[i] = √(length(loos[i]) * var(diff; corrected=false))
    end
    data = elpd_diff
    data = hcat(data, se_diff)

    sumval = sum([exp(psis_values[i]) for i in 1:nmodels])
    @. weight = exp(psis_values) / sumval
    data = hcat(data, weight)

    # Create KeyedArray object

    table = KeyedArray(
        data,
        model = mnames,
        statistic = [:elpd_diff, :se_diff, :weight],
    )

    # Return LooCompare object

    LooCompare(psis, table)

end

function Base.show(io::IO, ::MIME"text/plain", loo_compare::LooCompare)
    table = loo_compare.table
    return pretty_table(
        table;
        compact_printing=false,
        header=table.statistic,
        row_names=table.model,
        formatters=ft_printf("%5.2f"),
        alignment=:r,
    )
end

export
    LooCompare,
    loo_compare


text:
using IterTools, Base.Test

import Base: IsInfinite, SizeUnknown, HasLength, iteratorsize, HasShape

import Base.Iterators: take, countfrom, drop

include("testing_macros.jl")

@testset "IterTools" begin
@testset "iterators" begin
    @testset "ncycle" begin
        ncy1 = ncycle(0:3,3)

        i = 0
        for j in ncy1
            @test j == i % 4
            i += 1
        end

        @test eltype(ncy1) == Int
        i = 0
        for j in collect(ncy1)
            @test j == i % 4
            i += 1
        end
    end

    @testset "repeatedly" begin
        i = 0
        for j = repeatedly(() -> 1, 10)
            @test j == 1
            i += 1
        end
        @test i == 10
        for j = repeatedly(() -> 1)
            @test j == 1
            i += 1
            i <= 10 || break
        end
    end

    @testset "chain" begin
        ch1 = chain(1:2:5, 0.2:0.1:1.6)

        @test eltype(ch1) == typejoin(Int, Float64)
        @test collect(ch1) == [1:2:5; 0.2:0.1:1.6]

        ch2 = chain(1:0, 1:2:5, 0.2:0.1:1.6)

        @test eltype(ch2) == typejoin(Int, Float64)
        @test collect(ch2) == [1:2:5; 0.2:0.1:1.6]
        @test length(ch2) == length(collect(ch2))
        @test iteratorsize(ch2) == HasLength()

        ch3 = chain(1:10, 1:10, 1:10)
        @test length(ch3) == 30
        @test iteratorsize(ch3) == HasLength()

        r = countfrom(1)
        ch4 = chain(1:10, countfrom(1))
        @test eltype(ch4) == Int
        @test_throws MethodError length(ch4)
        @assert iteratorsize(r) == IsInfinite()
        @test iteratorsize(ch4) == IsInfinite()

        ch5 = chain()
        @test length(ch5) == 0
        @test iteratorsize(ch5) == HasLength()

        c = chain(ch1, ch2, ch3)
        @test length(c) == length(ch1) + length(ch2) + length(ch3)
        @test collect(c) == [collect(ch1); collect(ch2); collect(ch3)]

        r = rand(2,2)
        c = chain(r, r)
        @test length(c) == 8
        @test collect(c) == [vec(r); vec(r)]
        @test iteratorsize(r) == HasShape()
        @test iteratorsize(c) == HasLength()

        r = distinct(collect(1:10))
        @test iteratorsize(r) == SizeUnknown() #lazy filtering
        c = chain(1:10, r)
        @test_throws MethodError length(c)
        @test length(collect(c)) == 20
        @test iteratorsize(c) == SizeUnknown()
    end

    @testset "product" begin
        x1 = 1:2:10
        x2 = 1:5

        p0 = product(x1, x2)

        @test eltype(p0) == Tuple{Int, Int}
        @test collect(p0) == vec([(y1, y2) for y1 in x1, y2 in x2])

        p1 = product()

        @test eltype(p1) == Tuple{}
        @test length(p1) == 1
        @test collect(p1) == [()]
    end

    @testset "distinct" begin
        x = [5, 2, 2, 1, 2, 1, 1, 2, 4, 2]
        di0 = distinct(x)
        @test eltype(di0) == Int
        @test collect(di0) == unique(x)
    end

    @testset "partition" begin
        pa0 = partition(take(countfrom(1), 6), 2)
        @test eltype(pa0) == Tuple{Int, Int}
        @test collect(pa0) == [(1,2), (3,4), (5,6)]

        pa1 = partition(take(countfrom(1), 4), 2, 1)
        @test eltype(pa1) == Tuple{Int, Int}
        @test collect(pa1) == [(1,2), (2,3), (3,4)]

        pa2 = partition(take(countfrom(1), 8), 2, 3)
        @test eltype(pa2) == Tuple{Int, Int}
        @test collect(pa2) == [(1,2), (4,5), (7,8)]

        pa3 = partition(take(countfrom(1), 0), 2, 1)
        @test eltype(pa3) == Tuple{Int, Int}
        @test collect(pa3) == []

        @test_throws ArgumentError partition(take(countfrom(1), 8), 2, 0)
    end

    @testset "imap" begin
        function test_imap(expected, input...)
            result = collect(imap(+, input...))
            @test result == expected
        end

        @testset "empty arrays" begin
            test_imap(
                Any[],
                []
            )

            test_imap(
                Any[],
                Union{}[]
            )
        end

        @testset "simple operation" begin
            test_imap(
                Any[1,2,3],
                [1,2,3]
            )
        end

        @testset "multiple arguments" begin
            test_imap(
                Any[5,7,9],
                [1,2,3],
                [4,5,6]
            )
        end

        @testset "different-length arguments" begin
            test_imap(
                Any[2,4,6],
                [1,2,3],
                countfrom(1)
            )
        end
    end

    @testset "groupby" begin
        function test_groupby(input, expected)
            result = collect(groupby(x -> x[1], input))
            @test result == expected
        end

        @testset "empty arrays" begin
            test_groupby(
                [],
                Any[]
            )

            test_groupby(
                Union{}[],
                Any[]
            )
        end

        @testset "singletons" begin
            test_groupby(
                ["xxx"],
                Any[["xxx"]]
            )
        end

        @testset "typical operation" begin
            test_groupby(
                ["face", "foo", "bar", "book", "baz"],
                Any[["face", "foo"], ["bar", "book", "baz"]]
            )
        end

        @testset "trailing singletons" begin
            test_groupby(
                ["face", "foo", "bar", "book", "baz", "xxx"],
                Any[["face", "foo"], ["bar", "book", "baz"], ["xxx"]]
            )
        end

        @testset "leading singletons" begin
            test_groupby(
                ["xxx", "face", "foo", "bar", "book", "baz"],
                Any[["xxx"], ["face", "foo"], ["bar", "book", "baz"]]
            )
        end

        @testset "middle singletons" begin
            test_groupby(
                ["face", "foo", "xxx", "bar", "book", "baz"],
                Any[["face", "foo"], ["xxx"], ["bar", "book", "baz"]]
            )
        end
    end

    @testset "subsets" begin
        @testset "all lengths" begin
            s0 = subsets(Any[])
            @test eltype(eltype(s0)) == Any
            @test collect(s0) == Vector{Any}[Any[]]

            s1 = subsets([:a])
            @test eltype(eltype(s1)) == Symbol
            @test collect(s1) == Vector{Symbol}[Symbol[], Symbol[:a]]

            s2 = subsets([:a, :b, :c])
            @test eltype(eltype(s2)) == Symbol
            @test collect(s2) == Vector{Symbol}[
                Symbol[], Symbol[:a], Symbol[:b], Symbol[:a, :b], Symbol[:c],
                Symbol[:a, :c], Symbol[:b, :c], Symbol[:a, :b, :c],
            ]
        end

        @testset "specific length" begin
            sk0 = subsets(Any[],0)
            @test eltype(eltype(sk0)) == Any
            @test collect(sk0) == Vector{Any}[Any[]]

            sk1 = subsets([:a, :b, :c], 1)
            @test eltype(eltype(sk1)) == Symbol
            @test collect(sk1) == Vector{Symbol}[Symbol[:a], Symbol[:b], Symbol[:c]]

            sk2 = subsets([:a, :b, :c], 2)
            @test eltype(eltype(sk2)) == Symbol
            @test collect(sk2) == Vector{Symbol}[Symbol[:a,:b], Symbol[:a,:c], Symbol[:b,:c]]

            sk3 = subsets([:a, :b, :c], 3)
            @test eltype(eltype(sk3)) == Symbol
            @test collect(sk3) == Vector{Symbol}[Symbol[:a,:b,:c]]

            sk4 = subsets([:a, :b, :c], 4)
            @test eltype(eltype(sk4)) == Symbol
            @test collect(sk4) == Vector{Symbol}[]

            @testset for i in 1:5
                sk5 = subsets(collect(1:4), i)
                @test eltype(eltype(sk5)) == Int
                @test length(collect(sk5)) == binomial(4, i)
            end

            @testset "implicit conversions" begin
                sk10 = subsets(1:4, 3)
                @test eltype(eltype(sk10)) == Int
                @test length(collect(sk10)) == binomial(4, 3)

                sk11 = subsets(1:3, Int32(2))
                @test eltype(eltype(sk11)) == Int
                @test length(collect(sk11)) == binomial(3, 2)
            end
        end
    end

    @testset "nth" begin
        @testset for xs in Any[[1, 2, 3], 1:3, reshape(1:3, 3, 1)]
            @test nth(xs, 3) == 3
            @test_throws BoundsError nth(xs, 0)
            @test_throws BoundsError nth(xs, 4)
        end

        @testset for xs in Any[take(1:3, 3), drop(-1:3, 2)]
            @test nth(xs, 3) == 3
            @test_throws BoundsError nth(xs, 0)
        end

        s = subsets([1, 2, 3])
        @test_throws BoundsError nth(s, 0)
        @test_throws BoundsError nth(s, length(s) + 1)

        # #100
        @test nth(drop(repeatedly(() -> 1), 1), 1) == 1
    end

    @testset "takenth" begin
        tn0 = takenth(Any[], 10)
        @test eltype(tn0) == Any
        @test collect(tn0) == Any[]

        tn1 = takenth(Int[], 10)
        @test eltype(tn1) == Int
        @test collect(tn1) == Int[]

        @test_throws ArgumentError takenth([], 0)

        tn2 = takenth(10:20, 3)
        @test eltype(tn2) == Int
        @test collect(tn2) == [12,15,18]

        tn3 = takenth(10:20, 1)
        @test eltype(tn3) == Int
        @test collect(tn3) == collect(10:20)
    end

    @testset "peekiter" begin
        pi0 = peekiter(1:10)
        @test eltype(pi0) == Int
        @test collect(pi0) == collect(1:10)

        pi1 = peekiter([])
        @test eltype(pi1) == eltype([])
        @test collect(pi1) == collect([])

        it = peekiter([:a, :b, :c])
        @test eltype(it) == Symbol
        s = start(it)
        @test get(peek(it, s)) == :a

        it = peekiter([])
        s = start(it)
        @test isnull(peek(it, s))

        it = peekiter(1:10)
        s = start(it)
        x, s = next(it, s)
        @test get(peek(it, s)) == 2
    end
end

@testset "@itr" begin
    @testset "@zip" begin
        @test_zip [1,2,3] [:a, :b, :c] ['x', 'y', 'z']
        @test_zip [1,2,3] [:a, :b] ['w', 'x', 'y', 'z']
        @test_zip [1,2,3] Any[] ['w', 'x', 'y', 'z']
        @test_zip [1,2,3] Union{}[] ['w', 'x', 'y', 'z']
    end

    @testset "@enumerate" begin
        @test_enumerate [:a, :b, :c]
        @test_enumerate Union{}[]
        @test_enumerate Any[]
    end

    @testset "@take" begin
        @test_take [:a, :b, :c] 2
        @test_take [:a, :b, :c] 5
        @test_take [:a, :b, :c] 0
        @test_take Any[] 2
        @test_take Union{}[] 2
        @test_take Any[] 0
        @test_take [(:a,1), (:b,2), (:c,3)] 2
    end

    @testset "@takestrict" begin
        @test_takestrict [:a, :b, :c] 2
        @test_takestrict [:a, :b, :c] 3
        @test_takestrict [:a, :b, :c] 5
        @test_takestrict [:a, :b, :c] 0
        @test_takestrict Any[] 2
        @test_takestrict Union{}[] 2
        @test_takestrict Any[] 0
        @test_takestrict [(:a,1), (:b,2), (:c,3)] 2
        @test_takestrict [(:a,1), (:b,2), (:c,3)] 3
        @test_takestrict [(:a,1), (:b,2), (:c,3)] 4
    end

    @testset "@drop" begin
        @test_drop [:a, :b, :c] 2
        @test_drop [:a, :b, :c] 5
        @test_drop [:a, :b, :c] 0
        @test_drop Any[] 2
        @test_drop Union{}[] 2
        @test_drop Any[] 0
        @test_drop [(:a,1), (:b,2), (:c,3)] 2
    end

    @testset "@chain" begin
        @test_chain [1,2,3] [:a, :b, :c] ['x', 'y', 'z']
        @test_chain [1,2,3] [:a, :b] ['w', 'x', 'y', 'z']
        @test_chain [1,2,3] Any[] ['w', 'x', 'y', 'z']
        @test_chain [1,2,3] Union{}[] ['w', 'x', 'y', 'z']
        @test_chain [1,2,3] 4 [('w',3), ('x',2), ('y',1), ('z',0)]
    end
end
end


text:
"""
$(SIGNATURES)

Transform local flow variables to global frame
"""
function global_frame(w::AV, cosa, sina)
    if eltype(w) <: Int
        G = similar(w, Float64)
    else
        G = similar(w)
    end

    if length(w) == 2
        G[1] = w[1] * cosa - w[2] * sina
        G[2] = w[1] * sina + w[2] * cosa
    elseif length(w) == 4
        G[1] = w[1]
        G[2] = w[2] * cosa - w[3] * sina
        G[3] = w[2] * sina + w[3] * cosa
        G[4] = w[4]
    else
        throw("local -> global: dimension dismatch")
    end

    return G
end

global_frame(w::AV, n::AV) = global_frame(w, n[1], n[2])

"""
$(SIGNATURES)
"""
function global_frame(w::AV, dirccos::AM)
    if eltype(w) <: Int
        G = similar(w, Float64)
    else
        G = similar(w)
    end

    if length(w) == 3
        G[1] = w[1] * dirccos[1, 1] + w[2] * dirccos[2, 1] + w[3] * dirccos[3, 1]
        G[2] = w[1] * dirccos[1, 2] + w[2] * dirccos[2, 2] + w[3] * dirccos[3, 2]
        G[3] = w[1] * dirccos[1, 3] + w[2] * dirccos[2, 3] + w[3] * dirccos[3, 3]
    elseif length(w) == 5
        G[1] = w[1]
        G[2] = w[2] * dirccos[1, 1] + w[3] * dirccos[2, 1] + w[4] * dirccos[3, 1]
        G[3] = w[2] * dirccos[1, 2] + w[3] * dirccos[2, 2] + w[4] * dirccos[3, 2]
        G[4] = w[2] * dirccos[1, 3] + w[3] * dirccos[2, 3] + w[4] * dirccos[3, 3]
        G[5] = w[5]
    else
        throw("local -> global: dimension dismatch")
    end

    return G
end

"""
$(SIGNATURES)

Transform global flow variables to local frame
"""
function local_frame(w::AV, cosa, sina)
    if eltype(w) <: Int
        L = similar(w, Float64)
    else
        L = similar(w)
    end

    if length(w) == 2
        L[1] = w[1] * cosa + w[2] * sina
        L[2] = w[2] * cosa - w[1] * sina
    elseif length(w) == 4
        L[1] = w[1]
        L[2] = w[2] * cosa + w[3] * sina
        L[3] = w[3] * cosa - w[2] * sina
        L[4] = w[4]
    else
        throw("global -> local: dimension dismatch")
    end

    return L
end

local_frame(w::AV, n) = local_frame(w::AV, n[1], n[2])

"""
$(SIGNATURES)
"""
function local_frame(w::AV, dirccos::AM)
    if eltype(w) <: Int
        L = similar(w, Float64)
    else
        L = similar(w)
    end

    if length(w) == 3
        L[1] = w[1] * dirccos[1, 1] + w[2] * dirccos[1, 2] + w[3] * dirccos[1, 3]
        L[2] = w[1] * dirccos[2, 1] + w[2] * dirccos[2, 2] + w[3] * dirccos[2, 3]
        L[3] = w[1] * dirccos[3, 1] + w[2] * dirccos[3, 2] + w[3] * dirccos[3, 3]
    elseif length(w) == 5
        L[1] = w[1]
        L[2] = w[2] * dirccos[1, 1] + w[3] * dirccos[1, 2] + w[4] * dirccos[1, 3]
        L[3] = w[2] * dirccos[2, 1] + w[3] * dirccos[2, 2] + w[4] * dirccos[2, 3]
        L[4] = w[2] * dirccos[3, 1] + w[3] * dirccos[3, 2] + w[4] * dirccos[3, 3]
        L[5] = w[5]
    else
        throw("global -> local: dimension dismatch")
    end

    return L
end

"""
$(SIGNATURES)
"""
function local_velocity(u, v, cosa, sina)
    vn = @. u * cosa + v * sina
    vt = @. v * cosa - u * sina

    return vn, vt
end

local_velocity(u, v, n) = local_velocity(u, v, n[1], n[2])

"""
$(SIGNATURES)

Calculate unit normal vector
"""
function unit_normal(p1::T, p2::T) where {T<:AV}
    Δ = p2 .- p1
    l = norm(Δ) + 1e-6

    return [-Δ[2], Δ[1]] ./ l
end

"""
$(SIGNATURES)
"""
function unit_normal(p1::T, p2::T, p3::T) where {T<:AV}
    v1 = p2 .- p1
    v2 = p3 .- p1

    n = cross(v1, v2)
    l = norm(n) + 1e-6

    return n ./ l
end

"""
$(SIGNATURES)

Calculate point-point/line/surface distance
"""
point_distance(p1::T, p2::T) where {T<:AV} = norm(p1 .- p2)

"""
$(SIGNATURES)
"""
function point_distance(p::T, p1::T, p2::T) where {T<:AV}
    x0, y0 = p
    x1, y1 = p1
    x2, y2 = p2

    return abs((x2 - x1) * (y1 - y0) - (x1 - x0) * (y2 - y1)) /
           sqrt((x2 - x1)^2 + (y2 - y1)^2 + 1e-6)
end


text:
using LicenseCheck
using Test

MIT = """
    MIT License Copyright (c) <year> <copyright holders>

    Permission is hereby granted, free of charge, to any person obtaining a copy
    of this software and associated documentation files (the "Software"), to deal
    in the Software without restriction, including without limitation the rights
    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
    copies of the Software, and to permit persons to whom the Software is furnished
    to do so, subject to the following conditions:

    The above copyright notice and this permission notice (including the next
    paragraph) shall be included in all copies or substantial portions of the
    Software.

    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
    FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS
    OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
    OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE."""

Latex2e = """
    Copyright (C) 2007, 2008, 2009, 2010 Karl Berry.

    Copyright (C) 1988, 1994, 2007 Stephen Gilmore.

    Copyright (C) 1994, 1995, 1996 Torsten Martinsen.

    Permission is granted to make and distribute verbatim copies of this manual
    provided the copyright notice and this permission notice are preserved on
    all copies.

    Permission is granted to copy and distribute modified versions of this manual
    under the conditions for verbatim copying, provided that the entire resulting
    derived work is distributed under the terms of a permission notice identical
    to this one.

    Permission is granted to copy and distribute translations of this manual into
    another language, under the above conditions for modified versions."""

dorian_gray = """
    The artist is the creator of beautiful things.  To reveal art and
    conceal the artist is art's aim.  The critic is he who can translate
    into another manner or a new material his impression of beautiful
    things.

    The highest as the lowest form of criticism is a mode of autobiography.
    Those who find ugly meanings in beautiful things are corrupt without
    being charming.  This is a fault.

    Those who find beautiful meanings in beautiful things are the
    cultivated.  For these there is hope.  They are the elect to whom
    beautiful things mean only beauty.

    There is no such thing as a moral or an immoral book.  Books are well
    written, or badly written.  That is all."""

@testset "LicenseCheck" begin
    @testset "`licensecheck`" begin
        result = licensecheck(MIT)
        @test result.licenses_found == ["MIT"]
        @test result.license_file_percent_covered ≈ 100.0 atol = 2

        result = licensecheck(MIT * "\n" * Latex2e)
        @test result.licenses_found == ["MIT", "Latex2e"]
        @test result.license_file_percent_covered ≈ 100.0 atol = 2

        result = licensecheck(MIT * "\n" * dorian_gray)
        @test result.licenses_found == ["MIT"]
        @test result.license_file_percent_covered ≈ 100 * length(MIT) / (length(dorian_gray) + length(MIT)) atol = 5

        result = licensecheck(MIT * "\n" * dorian_gray * "\n" * Latex2e)
        @test result.licenses_found == ["MIT", "Latex2e"]
        @test result.license_file_percent_covered ≈
            100 * (length(MIT) + length(Latex2e)) /
            (length(dorian_gray) + length(MIT) + length(Latex2e)) atol = 5
    end

    @testset "`is_osi_approved`" begin
        @test is_osi_approved("MIT") == true
        @test is_osi_approved("LGPL-3.0") == true
        @test is_osi_approved("ABC") == false

        @test is_osi_approved(find_license(pkgdir(LicenseCheck))) == true
        @test !is_osi_approved(nothing) # for if `find_license` returns `nothing`
        @test is_osi_approved((; licenses_found = ["MIT", "MIT"]))
        @test !is_osi_approved((; licenses_found = String[]))
    end

    @testset "`find_licenses_*`" begin
        fl = find_license(joinpath(@__DIR__, ".."))
        # check it found the right one
        @test fl.license_filename == "LICENSE"
        @test fl.licenses_found == ["MIT"]
        @test fl.license_file_percent_covered > 90

        for method in (find_licenses, dir -> find_licenses(dir; allow_brute=false), find_licenses_by_bruteforce, find_licenses_by_list_intersection)
            results = method(joinpath(@__DIR__, ".."))
            @test only(results) == fl
        end

        # this can return more than 1 result due to case-insenitivity issues
        results = find_licenses_by_list(joinpath(@__DIR__, ".."))
        @test fl ∈ results
    end
end


text:
module MochaClassification

using Compat
using Mocha
using MochaClassificationWrapper
using MLUtilities
using Utilities
using DataModel

export MochaClassificationModel

"""
Neural Network Classification using Mocha

fields:
    - backend: Mocha Backend (CPU or GPU)
    - mochaNet: Mocha wrapper
    - labels
    - ntrain: number of data items
    - batchsize
"""
type MochaClassificationModel <: AbstractDataModel
    backend::Backend
    mochaNet::MochaClassificationWrapper.MWrap
    labels::Array{Int64,1}
    ntrain::Float64 # training set size
    batchsize::Int64
end

function MochaClassificationModel(
    datax,datac,
    modelfactory::Function,
    backend::Backend;
    ntrain::Int = length(datac),
    batchsize::Int = 100,
    do_shuffle::Bool = true,
    do_accuracy::Bool = true,
    do_predprob::Bool = false
    )

    data_layer = MemoryDataLayer(name = "data",
                               data = Array[datax,datac],
                               batch_size = batchsize,
                               shuffle = do_shuffle)

    mochaNet = MochaClassificationWrapper.MWrap(data_layer,
                                 modelfactory,
                                 "MochaClassificationNet",
                                 do_accuracy,
                                 do_predprob,
                                 backend)
    MochaClassificationModel(backend,mochaNet,datac[:],ntrain,batchsize)
end

"""
`DataModel.getllik(dms::MochaClassificationModel)` returns a function that provides the full loglikelihood
"""
function DataModel.getllik(dms::MochaClassificationModel)
    function llik(x)
        (accuracy, loglikelihood) = MochaClassificationWrapper.evaluateTestNN(dms.mochaNet, x, dms.batchsize)
        return loglikelihood
    end
end

"""
`DataModel.getgrad(dms::MochaClassificationModel)` returns function returning the stochastic gradient function.
"""
function DataModel.getgrad(dms::MochaClassificationModel)
    function grad(x)
        (llik, grad) = MochaClassificationWrapper.forwardbackward(dms.mochaNet, x, regu_coef = 0.)
        return grad * dms.ntrain
    end
end

"""
`getparams(dms::MochaClassificationModel)` reads out the current weights.
"""
function getparams(dms::MochaClassificationModel)
    MochaClassificationWrapper.getparams(dms.mochaNet)
end

"""
`setparams!(dms::MochaClassificationModel,para::Array{Float64,1` sets the weights to `para`.
"""
function setparams!(dms::MochaClassificationModel,para::Array{Float64,1})
    MochaClassificationWrapper.setparams!(dms.mochaNet,para)
end
"""
`fetchnparams(dms::MochaClassificationModel)` returns the number of parameters of the model.
"""
function fetchnparams(dms::MochaClassificationModel)
    MochaClassificationWrapper.getnparams(dms.mochaNet)
end

#-------------------
# Initialisation functions
#-------------------
"""
Xavier Initialisation
"""
function init_xavier(dms::MochaClassificationModel)
    MochaClassificationWrapper.init_xavier(dms.mochaNet)
end
"""
Initialisation based on fan-in.
"""
function init_simple_fanin(dms::MochaClassificationModel)
    MochaClassificationWrapper.init_simple_fanin(dms.mochaNet)
end
"""
Gaussian Initialisation
"""
function init_gaussian(dms::MochaClassificationModel, initvar::Float64)
    MochaClassificationWrapper.init_gaussian(dms.mochaNet, initvar)
end
"""
Uniform Initialisation
"""
function init_uniform(dms::MochaClassificationModel, initvar::Float64)
    MochaClassificationWrapper.init_uniform(dms.mochaNet, initvar)
end

#-------------------
# evaluation functions
#-------------------
"""
Evaluate accuracy based on predictive probabilities
"""
function evaluate(dms::MochaClassificationModel,
    predprobs::Array{Float64,2})
    prediction =  (findmax(predprobs,1)[2]-1.) % size(predprobs)[1]
    accuracy = sum(dms.labels .== prediction[:])./(length(dms.labels)+0.0)
    loglikelihood = sum(log([predprobs[i,j] for (i,j) in zip(dms.labels+1,1:round(Int,dms.ntrain))]))
    @dict(accuracy, loglikelihood)
end

"""
Evaluate accuracy at weights `x`.
"""
function evaluate(dms::MochaClassificationModel,
    x::Vector{Float64})
    @assert length(x) == fetchnparams(dms)
    (accuracy, loglikelihood) = MochaClassificationWrapper.evaluateTestNN(dms.mochaNet, x, dms.batchsize)
    @dict(accuracy, loglikelihood)
end

"""
Get predictive probabilities given weights `x`
"""
function getpredprobs(dms::MochaClassificationModel,
    x::Vector{Float64})
    predictiveprobs = MochaClassificationWrapper.getpredprobs(dms.mochaNet, x, dms.batchsize)
    @dict(predictiveprobs)
end

end