text:
# Utilities for creating and modifying vector spaces and TensorMaps.
# To be `included` in MERAKit.jl.

"""
An `AbstractTensorMap` from N indices to N indices.
"""
SquareTensorMap{N} = AbstractTensorMap{S1, N, N} where {S1}

"""
A `Union` type of the different `TensorKitManifolds` tangent types: `GrassmannTangent`,
`StiefelTangent`, and `UnitaryTangent`.
"""
Tangent = Union{Grassmann.GrassmannTangent, Stiefel.StiefelTangent, Unitary.UnitaryTangent}

"""
    disentangler_type(::Type{ST}, ::Type{ET}, Tan::Bool)

Given the `IndexSpace` type `ST` and storage type `ET`, return the concrete type of a 2-to-2
disentangler. `Tan` is a flag for whether we want the `TensorMap` type itself (`Tan =
false`) or the corresponding tangent type, i.e. a `StiefelTangent`.

See also: [`ternaryisometry_type`](@ref), [`binaryisometry_type`](@ref)
"""
function disentangler_type(::Type{ST}, ::Type{ET}, Tan::Bool) where {ST, ET}
    TensorType = tensormaptype(ST, 2, 2, ET)
    if Tan
        TA = tensormaptype(ST, 2, 2, ET)
        DisType = Stiefel.StiefelTangent{TensorType, TA}
    else
        DisType = TensorType
    end
    return DisType
end

"""
    ternaryisometry_type(::Type{ST}, ::Type{ET}, Tan::Bool)

Given the `IndexSpace` type `ST` and storage type `ET`, return the concrete type of a 3-to-1
isometry. `Tan` is a flag for whether we want the `TensorMap` type itself (`Tan =
false`) or the corresponding tangent type, i.e. a `GrassmannTangent`.

See also: [`disentangler_type`](@ref), [`binaryisometry_type`](@ref)
"""
function ternaryisometry_type(::Type{ST}, ::Type{ET}, Tan::Bool) where {ST, ET}
    TensorType = tensormaptype(ST, 3, 1, ET)
    if Tan
        TU = tensormaptype(ST, 3, 1, ET)
        SET = isreal(sectortype(ST)) ? real(ET) : ET
        TS = tensormaptype(ST, 1, 1, SET)
        TV = tensormaptype(ST, 1, 1, ET)
        IsoType = Grassmann.GrassmannTangent{TensorType, TU, TS, TV}
    else
        IsoType = TensorType
    end
    return IsoType
end

"""
    binaryisometry_type(::Type{ST}, ::Type{ET}, Tan::Bool)

Given the `IndexSpace` type `ST` and storage type `ET`, return the concrete type of a 2-to-1
isometry. `Tan` is a flag for whether we want the `TensorMap` type itself (`Tan =
false`) or the corresponding tangent type, i.e. a `GrassmannTangent`.

See also: [`disentangler_type`](@ref), [`ternaryinaryisometry_type`](@ref)
"""
function binaryisometry_type(::Type{ST}, ::Type{ET}, Tan::Bool) where {ST, ET}
    TensorType = tensormaptype(ST, 2, 1, ET)
    if Tan
        TU = tensormaptype(ST, 2, 1, ET)
        SET = isreal(sectortype(ST)) ? real(ET) : ET
        TS = tensormaptype(ST, 1, 1, SET)
        TV = tensormaptype(ST, 1, 1, ET)
        IsoType = Grassmann.GrassmannTangent{TensorType, TU, TS, TV}
    else
        IsoType = TensorType
    end
    return IsoType
end

"""
    randomisometry(T, Vout, Vin)

Given an element type `T` and two vector spaces, `Vin` for the domain and `Vout` for the
codomain, return a Haar random isometry.

The implementation uses a QR decomposition of a Gaussian random matrix.
"""
randomisometry(::Type{T}, Vout, Vin) where {T} = TensorMap(randhaar, T, Vout ← Vin)

"""
    initialize_disentangler(T, Vout, Vin, random::Bool)

Initialize a disentangler from `Vin ⊗ Vin` to `Vout ⊗ Vout`, of element type `T`.

The returned tensor is Haar random if `random = true`. If `random = false` and `Vin == Vout`
it is the identity. If `random = false` and `Vin != Vout` it is the tensor product of two
one-site Haar random unitaries.
"""
function initialize_disentangler(::Type{T}, Vout, Vin, random::Bool) where {T}
    if random
        u = randomisometry(T, Vout ⊗ Vout, Vin ⊗ Vin)
    else
        if Vin == Vout
            u = isomorphism(Matrix{T}, Vout ⊗ Vout, Vin ⊗ Vin)
        else
            uhalf = randomisometry(T, Vout, Vin)
            u = uhalf ⊗ uhalf
        end
    end
    return u
end

"""
    support(op)

Return the number of sites/indices `N` that the operator `op` is supported on, assuming it
is an operator from `N` sites to `N` sites.

See also: [`expand_support`](@ref)
"""
support(op::SquareTensorMap{N}) where {N} = N

"""
    expand_support(op, n::Integer)

Given an operator from a `N` indices to `N` indices, expand its support to a larger number
of indices `n` by tensoring with the identity. The different ways of doing the expansion,
e.g. I ⊗ op and op ⊗ I, are averaged over.

See also: [`support`](@ref)
"""
@inline expand_support(op::SquareTensorMap, n::Int) = _expand_support(op, Val(n))

@noinline function _expand_support(op::SquareTensorMap{N}, ::Val{n}) where {N, n}
    if n <= N
        return op
    else
        dom = ProductSpace(ntuple(i->space(op, 1), n)...)
        op2 = fill!(similar(op, dom, dom), 0)
        V = space(op, 1)
        eye = id(storagetype(op), V)
        for k = 0:n-N
            eyes1 = Base.fill_to_length((), eye, Val(k))
            eyes2 = Base.fill_to_length((), eye, Val(n-N-k))
            coeff = factorial(n-N)/(factorial(k)*factorial(n-N-k)) / 2^(n-N)
            axpy!(coeff, ⊗(eyes1..., op, eyes2...), op2)
            # TODO This would generate the uniform sum. See if it changes something.
            # axpy!(1/(n-N+1), ⊗(eyes1..., op, eyes2...), op2)
        end
        return op2
    end
end

"""
    remove_symmetry(V)

Strip a vector space of its symmetry structure, i.e. return the corresponding
`ℂ^n` or `ℝ^n`.
"""
remove_symmetry(V::ElementarySpace{ℝ}) = CartesianSpace(dim(V))
remove_symmetry(V::ElementarySpace{ℂ}) = ComplexSpace(dim(V), isdual(V))

"""
    remove_symmetry(t::AbstractTensorMap)

Strip an `AbstractTensorMap` of its internal symmetries, and return the corresponding
`TensorMap` that operates on `ComplexSpace` or `CartesianSpace`.
"""
function remove_symmetry(t::TensorMap)
    dom = domain(t)
    cod = codomain(t)
    S = field(spacetype(t)) === ℝ ? CartesianSpace : ComplexSpace
    domain_nosym = ProductSpace{S}(map(remove_symmetry, tuple(dom...)))
    codomain_nosym = ProductSpace{S}(map(remove_symmetry, tuple(cod...)))
    t_nosym = TensorMap(convert(Array, t), codomain_nosym ← domain_nosym)
    return t_nosym
end

"""
    expand_vectorspace(V, newdim)

Given a vector space `V` and a dictionary `newdim` of dimensions for the various irrep
sectors, return another vector space of the same type, but with these new dimension. If some
irrep sectors are not in the dictionary, the dimensions of the original space are used.
"""
function expand_vectorspace(V::CartesianSpace, newdim)
    d = length(newdim) > 0 ? first(values(newdim)) : dim(V)
    return typeof(V)(d)
end

function expand_vectorspace(V::ComplexSpace, newdim)
    d = length(newdim) > 0 ? first(values(newdim)) : dim(V)
    return typeof(V)(d, V.dual)
end

function expand_vectorspace(V::GeneralSpace, newdim)
    d = length(newdim) > 0 ? first(values(newdim)) : dim(V)
    return typeof(V)(d, V.dual, V.conj)
end

function expand_vectorspace(V::GradedSpace, newdims)
    olddims = Dict(s => dim(V, s) for s in sectors(V))
    sectordict = merge(olddims, newdims)
    return typeof(V)(sectordict; dual = V.dual)
end

# If the first argument given to depseudoserialize is a String, we assume its a
# representation of a an object that can `eval`uated. So we evaluate it and call
# depseudoserialize again.
depseudoserialize(str::String, args...) = depseudoserialize(eval(Meta.parse(str)), args...)

function pseudoserialize(t::T) where T <: TensorMap
    # We make use of the nice fact that many TensorKit objects return on repr strings that
    # are valid syntax to reconstruct these objects.
    domstr = repr(t.dom)
    codomstr = repr(t.codom)
    eltyp = eltype(t)
    if isa(t.data, AbstractArray)
        data = deepcopy(t.data)
    else
        data = Dict(repr(s) => deepcopy(d) for (s, d) in t.data)
    end
    return repr(T), domstr, codomstr, eltyp, data
end

function depseudoserialize(
    ::Type{T}, domstr, codomstr, eltyp, data
) where {T <: AbstractTensorMap}
    dom = eval(Meta.parse(domstr))
    codom = eval(Meta.parse(codomstr))
    t = TensorMap(zeros, eltyp, codom ← dom)
    if isa(t.data, AbstractArray)
        t.data[:] = data
    else
        for (irrepstr, irrepdata) in data
            irrep = eval(Meta.parse(irrepstr))
            t.data[irrep][:] = irrepdata
        end
    end
    return t
end

"""
    append_dummy_index(t::TensorMap)

Create a copy of `t`, but with an extra, trivial index as the last index of the domain.
"""
function append_dummy_index(t::TensorMap)
    return TensorMap(t.data, codomain(t), insertunit(domain(t)))
end

"""
    remove_dummy_index(t::TensorMap)

Assuming the last index of the domain of `t` is a trivial dummy index, return a `t` but with
that index removed.
"""
function remove_dummy_index(t::TensorMap{S}) where {S}
    @assert space(t, numind(t)) == oneunit(S)'
    dom = ProductSpace{S}(Base.front(domain(t).spaces))
    return TensorMap(t.data, codomain(t), dom)
end

"""
    pad_with_zeros_to(t::AbstractTensorMap, spacedict::Dict)

Transform `t` to change the vector spaces of its indices, by throwing elements away or
padding the tensor with zeros.

`spacedict` is a dictionary with index labels (1, 2, 3, ...) as keys, and `VectorSpace`s as
values. It tells us which indices should have their space changed, and to what. Instead of a
dictionary, a varargs of `Pair`s `index => vectorspace` also works.

For each index `i`, its current space `Vorig = space(t, i)` and new space
`Vnew = spacedict[i]` should be of the same type. If `Vnew` is strictly larger than `Vold`
then `t` is padded with zeros to fill in the new elements. Otherwise some elements of `t`
will be truncated away.
"""
function pad_with_zeros_to(t::AbstractTensorMap, spacedict::Dict)
    S = spacetype(t)
    newcodomainspaces = tuple(codomain(t)...)
    for i = 1:length(newcodomainspaces)
        if haskey(spacedict, i)
            @assert isdual(spacedict[i]) == isdual(space(t, i))
            # unpredictable behaviour otherwise
            newcodomainspaces = Base.setindex(newcodomainspaces, spacedict[i], i)
        end
    end
    newdomainspaces = tuple(domain(t)...)
    for i = 1:length(newdomainspaces)
        j = i + numout(t)
        if haskey(spacedict, j)
            @assert isdual(spacedict[j]) == isdual(space(t, j))
            # unpredictable behaviour otherwise
            newdomainspaces = Base.setindex(newdomainspaces, spacedict[j]', i)
        end
    end
    newcodomain = ProductSpace{S}(newcodomainspaces)
    newdomain = ProductSpace{S}(newdomainspaces)
    tnew = fill!(similar(t, newcodomain, newdomain), zero(eltype(t)))
    for (f1, f2) in fusiontrees(t)
        a = t[f1, f2]
        anew = tnew[f1, f2]
        axes = Base.OneTo.(min.(size(a), size(anew)))
        copyto!(view(anew, axes...), view(a, axes...))
    end
    return tnew
end

pad_with_zeros_to(t::AbstractTensorMap, spaces...) = pad_with_zeros_to(t, Dict(spaces))

"""
    gershgorin_bounds(t::AbstractTensorMap)

For a Hermitian square tensor (a square matrix after fusing the domain and codomain into
single indices), return a lower and upper bound between which all its eigenvalues lie.

This costs O(D^2) time, where D is the matrix dimension.

See also: [`gershgorin_discs`](@ref)
"""
function gershgorin_bounds(t::SquareTensorMap)
    v = real(first(last(first(blocks(t))))) # get first element from first block
    lb, ub = v, v
    largest_bound((lb1, ub1), (lb2, ub2)) = (min(lb1, lb2), max(ub1, ub2))
    for (c, b) in blocks(t)
        lb, ub = largest_bound((lb, ub), gershgorin_bounds(b))
    end
    return lb, ub
end

function gershgorin_bounds(a::Array{S, 2}) where {S}
    nonhermiticity = norm(a - a')/norm(a)
    if nonhermiticity > 1e-12
        msg = "Computing gershgorin_bounds, ignoring non-hermiticity of $(nonhermiticity)."
        @warn(msg)
    end
    discs = gershgorin_discs(a)
    lb = minimum(real(c) - r for (c, r) in discs)
    ub = maximum(real(c) + r for (c, r) in discs)
    return lb, ub
end

function gershgorin_discs(a::Array{S, 2}) where {S}
    d = size(a, 1)
    if size(a, 2) != d
        msg = "Can't apply Gershgorin disc theorem to a non-square matrix."
        throw(ArgumentError(msg))
    end
    centres = diag(a)
    abs_centres = abs.(centres)
    radii1 = dropdims(sum(abs.(a), dims = 1), dims = 1) .- abs_centres
    radii2 = dropdims(sum(abs.(a), dims = 2), dims = 2) .- abs_centres
    radii = min.(radii1, radii2)
    discs = collect(zip(centres, radii))
    return discs
end

"""
    densitymatrix_entropy(rho)

Compute the von Neumann entropy of a density matrix `rho`.
"""
function densitymatrix_entropy(rho)
    eigs = eigh(rho)[1]
    S = -dot(eigs, log(eigs))
    return S
end

"""
    precondition_tangent(X::Tangent, rho::AbstractTensorMap,
                         delta = precondition_regconst(X))

Precondition the tangent vector `X` with the metrix g(X, Y) =  Tr[X' Y rho], where `rho`
positive definite tensor.

This works for `X` being a `StiefelTangent`, `GrassmannTangent`, or `UnitaryTangent`.

Preconditioning requires inverting `rho`. `delta` is a threshold parameter for regularising
that inverse. The regularised inverse is S -> 1/sqrt(S^2 + delta^2).

See also: [`precondition_regconst`](@ref)
"""
function precondition_tangent(
    X::Stiefel.StiefelTangent, rho::AbstractTensorMap, delta = precondition_regconst(X)
)
    W, A, Z = X.W, X.A, X.Z
    E, U = eigh(rho)
    Einv = reginv_E(E, delta)
    rhoinv = U * Einv * U'
    Z_prec = projectcomplement!(Z * rhoinv, W)
    A_prec = projectantihermitian!(symmetric_sylvester(E, U, 2*A, delta))
    return Stiefel.StiefelTangent(W, A_prec, Z_prec)
end

function precondition_tangent(
    X::Grassmann.GrassmannTangent, rho::AbstractTensorMap, delta = precondition_regconst(X)
)
    W, Z = X.W, X.Z
    E, U = eigh(rho)
    Einv = reginv_E(E, delta)
    rhoinv = U * Einv * U'
    Z_prec = projectcomplement!(Z * rhoinv, W)
    return Grassmann.GrassmannTangent(W, Z_prec)
end

function precondition_tangent(
    X::Unitary.UnitaryTangent, rho::AbstractTensorMap, delta = precondition_regconst(X)
)
    W, A = X.W, X.A
    E, U = eigh(rho)
    A_prec = projectantihermitian!(symmetric_sylvester(E, U, 2*A, delta))
    return Unitary.UnitaryTangent(W, A_prec)
end

function reginv_E(E, delta)
    E_isreal = eltype(E) <: Real && isreal(sectortype(E))
    Ereg = sqrt(E^2 + delta^2*id(domain(E)))
    if E_isreal
        return inv(real(Ereg))
    else
        return inv(Ereg)
    end
end

"""
    precondition_regconst(X::Tangent)

The default regularisation constant to use when inverting the density matrix `rho` in
preconditioning.

See also: [`precondition_tangent`](@ref)
"""
function precondition_regconst(X::Tangent)
    delta = sqrt(eps(real(float(one(eltype(X.W))))))
    delta = max(delta, norm(X) / 100.0)
    return delta
end

"""
    symmetric_sylvester(E, U, C, delta)

Solve the Sylvester equation `A*X + X*A = C`, where we know that `A = A'` and the arguments
`E` and `U` are the eigenvalue decomposition of `A = U*E*U'`. This requires performing a
matrix inversion of `1 ⊗ A + A ⊗ 1`, and that inverse is regularised as
`X -> 1/sqrt(X^2 + delta^2)`.
"""
function symmetric_sylvester(
    E::AbstractTensorMap, U::AbstractTensorMap, C::AbstractTensorMap, delta
)
    cod = domain(C)
    dom = codomain(C)
    sylAB(c) = symmetric_sylvester(diag(block(E, c)), block(U, c), block(C, c), delta)
    data = TensorKit.SectorDict(c => sylAB(c) for c in blocksectors(cod ← dom))
    return TensorMap(data, cod ← dom)
end

function symmetric_sylvester(E::AbstractVector, U::AbstractMatrix, C::AbstractMatrix, delta)
    temp1 = typeof(C)(undef, size(C))
    temp2 = typeof(C)(undef, size(C))
    mul!(temp1, U', C)
    mul!(temp2, temp1, U)
    for i in CartesianIndices(temp2)
        temp2[i] /= sqrt((E[i[1]] + E[i[2]])^2 + delta^2)
    end
    mul!(temp1, U, temp2)
    mul!(temp2, temp1, U')
    return temp2
end


text:
function ln(x)
    sn = 0
    k = 0
    bk =  -(1 - x)
    while abs(bk) > eps()
        ak  = bk / (k+1)
        k += 1
        sn += ak
        bk *= (1-x)
    end

    return sn
end

function row2(x)
    sn = 1
    ak = -0.5 * x
    k = 1
    while abs(ak) > eps()
        sn += ak
        ak *= (x - k)*x
        k += 1
        ak = ak / k
    end

    return sn
end

function row3(x)
    uk = -x*x
    vk = 1
    wk = 1
    sn = 0
    k = 1
    while abs(uk*(vk + wk)) > eps()
        sn += uk*(vk + wk)
        uk = -uk * x * x
        vk = vk /(k + 1)
        wk = wk /2k/(2k + 1)
        k += 1
    end

    return sn
end

function test_rows()
    println(ln(0.5), " ", log(0.5))
    println(row2(0.5), " ", 1/sqrt(1.5))
    println(row3(0.5), " ", 0.5*sin(0.5) - exp(-0.25))
end

ext: jl


text:
using DataFrames, Query, Compat, GLM

function receiveLoop!(port,dataQueue,ctl)
  # Connect to server
  cli = connect(port)
  try
    while !ctl["shouldExit"]
      push!(dataQueue,readline(cli))
      sleep(0.1)
    end
  catch ex

  finally
    close(cli)
  end
end

function testServer(port,ctl)
  shouldExit = false
  # Open server for listening
  server = listen(port)
  try
    while !ctl["shouldExit"]
      sock = accept(server)
      while isopen(sock)
        println(sock,"hi there\n")
      end
      sleep(0.5)
    end
  catch ex

  finally
    close(server)
  end
end

function parseMsg(msg)
    println(msg)
end

###################################################

portnum = 5000
dataQueue = []
#threads = Dict("server"=>Dict{Any,Any}("ctl"=>Dict("shouldExit"=>false)),
#               "client"=>Dict{Any,Any}("ctl"=>Dict("shouldExit"=>false)))

#threads["server"]["thread"] = @async testServer(portnum,threads["server"]["ctl"])

threads = Dict("client"=>Dict{Any,Any}("ctl"=>Dict("shouldExit"=>false)))
threads["client"]["thread"] = @async receiveLoop!(portnum,
                                                  dataQueue,
                                                  threads["client"]["ctl"])

shouldExit = false
try
  while !shouldExit
    if !isempty(dataQueue)
      el = pop!(dataQueue)
      parseMsg(el)
    end
    sleep(0.1)
  end
catch ex
  #threads["server"]["ctl"]["shouldExit"] = true
  threads["client"]["ctl"]["shouldExit"] = true
end

max_stars_repo_name: umd-agrc/SimpleControlSim

max_issues_repo_name: umd-agrc/SimpleControlSim

max_forks_repo_name: umd-agrc/SimpleControlSim


text:
# Test basic methods
@testset "Basics" begin
    # test _macro_error
    @testset "_macro_error" begin
        source = LineNumberNode(42, :bob)
        @test_throws ErrorException InfiniteOpt._macro_error(:test, (:x,), source, "hello")
    end
    # test _esc_non_constant
    @testset "_esc_non_constant" begin
        @test InfiniteOpt._esc_non_constant(42) == 42
        @test InfiniteOpt._esc_non_constant(:(quote 42 end)) == :(quote 42 end)
        @test InfiniteOpt._esc_non_constant(:(x[i])) == esc(:(x[i]))
        @test InfiniteOpt._esc_non_constant(:x) == esc(:x)
    end
    # test _name_call
    @testset "_name_call" begin
        @test InfiniteOpt._name_call("bob", []) == "bob"
        @test InfiniteOpt._name_call("", []) == ""
        expected = :(string("x", "[", string($(Expr(:escape, "2"))), ",", string($(Expr(:escape, "1"))), "]"))
        @test InfiniteOpt._name_call("x", ["2", "1"]) == expected
    end
    # test _extract_kw_args
    @testset "_extract_kw_args" begin
        args = (:x, :(d=42), :(container=Cat))
        @test InfiniteOpt._extract_kw_args(args)[1] == [:x]
        @test InfiniteOpt._extract_kw_args(args)[2] == [:(d=42)]
        @test InfiniteOpt._extract_kw_args(args)[3] == :(Cat)
    end
end


text:
struct IntrinsicMetric{T, D, M<:Metric,
                       V<:SVector{D, T}, N<:NearestNeighbors.NNTree{V, M}} <: Metric
    points ::Vector{V}
    radius ::T
    metric ::M
    graph  ::SimpleWeightedGraph{Int, T}
    tree   ::N
end

adjgraph(im::IntrinsicMetric) = im.graph
points(im::IntrinsicMetric) = im.points
npoints(im::IntrinsicMetric) = length(im.points)
ambientdim(im::IntrinsicMetric{T, D}) where {T, D} = D
ambientmetric(im::IntrinsicMetric) = im.metric

function Base.show(io::IO, im::IntrinsicMetric{T, D, M, V, N}) where {T, D, M, V, N}
    println(io, "IntrinsicMetric{$T, $D, $M, $V, $N}")
    println(io, "  Number of points: $(length(im.points))")
    println(io, "  Dimensions: $D")
    println(io, "  Radius: $(im.radius)")
    println(io, "  Metric: $(im.metric)")
    print(io, "  Graph: $(im.graph)")
end

function IntrinsicMetric(data::AbstractMatrix{T}, r;
                         metric = Euclidean()) where T
    dim = size(data, 1)
    pts = reinterpret(SVector{dim, T}, data, (length(data) ÷ dim, ))
    IntrinsicMetric(pts, r, metric = metric)
end

function IntrinsicMetric{T, D, M<:Metric}(pts::AbstractVector{SVector{D, T}}, r;
                                          metric::M = Euclidean())
    tree = BallTree(pts, metric, reorder = false)

    neighs = inrange(tree, pts, r, false)
    is = Int[]; js = Int[]; ds = T[]
    for (i, ns) in enumerate(neighs)
        for j in ns
            i < j || continue
            append!(is, [i, j])
            append!(js, [j, i])
            d = evaluate(metric, pts[i], pts[j])
            append!(ds, [d, d])
        end
    end
    n = length(pts)
    g = SimpleWeightedGraph(sparse(is, js, ds, n, n))

    IntrinsicMetric{T, D, M, eltype(pts), typeof(tree)}(pts, T(r), metric, g, tree)
end

function Distances.evaluate(im::IntrinsicMetric, p1, p2)
    length(p1) == length(p2) == ambientdim(im) ||
        throw(DimensionMismatch("Incompatible dimensions!"))
    p1 == p2 && return 0

    nearest1 = knn(im.tree, p1, 1)
    nearest2 = knn(im.tree, p2, 1)

    total_dist = first(nearest1[2]) + first(nearest2[2])

    i1 = first(nearest1[1])
    i2 = first(nearest2[1])

    total_dist + dijkstra_shortest_paths(im.graph, i1).dists[i2]
end

# TODO: colwise: each col with each other col?

function Distances.pairwise(im::IntrinsicMetric,
                            a::AbstractMatrix, b::AbstractMatrix)
    res = zeros(promote_type(eltype(a), eltype(b)), size(a, 2), size(b, 2))
    pairwise!(res, im, a, b)
end

function Distances.pairwise!(res::AbstractMatrix, im::IntrinsicMetric,
                             a::AbstractMatrix, b::AbstractMatrix)
    n = size(a, 2)
    m = size(b, 2)
    size(a, 1) == size(b, 1) == ambientdim(im) ||
        throw(DimensionMismatch("Incompatible dimensions!"))
    size(res) == (n, m) ||
        throw(DimensionMismatch("Output matrix should be $n×$m, " *
                                "but is $(size(res, 1))×(size(res, 2))"))

    a_idx, a_dst = map.(first, knn(im.tree, a, 1))
    b_idx, b_dst = map.(first, knn(im.tree, b, 1))

    for i in 1:n
        dsts = dijkstra_shortest_paths(im.graph, a_idx[i]).dists[b_idx]
        for j in 1:m
            a[:, i] ≈ b[:, j] && continue
            res[i, j] = min(a_dst[i] + b_dst[j] + dsts[j])
        end
    end

    res
end

function Distances.pairwise(im::IntrinsicMetric, a::AbstractMatrix)
    res = zeros(eltype(a), size(a, 2), size(a, 2))
    pairwise!(res, im, a)
end

function Distances.pairwise!(res::AbstractMatrix, im::IntrinsicMetric,
                             a::AbstractMatrix)
    n = size(a, 2)
    size(a, 1) == ambientdim(im) ||
        throw(DimensionMismatch("Incompatible dimensions!"))
    size(res) == (n, n) ||
        throw(DimensionMismatch("Output matrix should be $n×$n, " *
                                "but is $(size(res, 1))×(size(res, 2))"))

    a_idx, a_dst = map.(first, knn(im.tree, a, 1))

    for i in 1:n
        dsts = dijkstra_shortest_paths(im.graph, a_idx[i]).dists[a_idx]
        for j in 1:i-1
            a[:, i] ≈ a[:, j] && continue
            d = min(a_dst[i] + a_dst[j] + dsts[j])
            res[i, j] = d
            res[j, i] = d
        end
    end

    res
end

function Distances.pairwise!(res::AbstractMatrix{T},
                             im::IntrinsicMetric{T, D}) where {T, D}
    pts = reinterpret(T, points(im), (D, npoints(im)))
    pairwise!(res, im, pts)
end

function Distances.pairwise(im::IntrinsicMetric{T, D}) where {T, D}
    pts = reinterpret(T, points(im), (D, npoints(im)))
    pairwise(im, pts)
end


text:
using DiffEqBase: AbstractODEProblem, DiscreteProblem
using Setfield: Lens, set, get, @set

const DEP{iip} = AbstractODEProblem{uType, tType, iip} where {uType, tType}

TimeKind(::Type{<: DiscreteProblem}) = Discrete()
TimeKind(::Type{<: AbstractODEProblem}) = Continuous()
StateKind(::Type{<: DEP{true}}) = MutableState()
StateKind(::Type{<: DEP{false}}) = ImmutableState()

struct DiffEqWrapper{P, L, D}
    de_prob::P
    param_axis::L
    param_diff::D
end

DiffEqWrapper(de_prob, param_axis) = DiffEqWrapper(de_prob, param_axis, nothing)

function diffeq_homotopy(H, x, p::DiffEqWrapper{<:DEP{true}}, t)
    q = set(p.de_prob.p, p.param_axis, t)
    p.de_prob.f(H, x, q, 0)
    maybe_subtract!(H, x, statekind(p.de_prob), timekind(p.de_prob))
end

function diffeq_homotopy(x, p::DiffEqWrapper{<:DEP{false}}, t)
    q = set(p.de_prob.p, p.param_axis, t)
    H = p.de_prob.f(x, q, 0)
    return maybe_subtract!(H, x, statekind(p.de_prob), timekind(p.de_prob))
end

maybe_subtract!(H, x, p) = maybe_subtract!(H, x, statekind(p), timekind(p))
maybe_subtract!(H, ::Any, ::StateKind, ::Continuous) = H
maybe_subtract!(H, x, ::MutableState, ::Discrete) = H .-= x
maybe_subtract!(H, x, ::ImmutableState, ::Discrete) = H .- x

struct MutableParamDiff{F, X, C}
    param_to_state::F
    x::X
    cfg::C
end

function MutableParamDiff(p::DiffEqWrapper)
    x = deepcopy(p.de_prob.u0)
    param_to_state = let x = x
        (H, t) -> diffeq_homotopy(H, x, p, t)
    end
    cfg = ForwardDiff.DerivativeConfig(param_to_state, x, 1.0)
    return MutableParamDiff(param_to_state, x, cfg)
end

function diffeq_homotopy_jacobian(H, J, x, p::DiffEqWrapper{<:DEP{true}}, t)
    p.param_diff :: MutableParamDiff
    p.param_diff.x .= x
    ForwardDiff.derivative!(
        (@view J[:, end]),
        p.param_diff.param_to_state,
        H,
        t,
        p.param_diff.cfg,
    )
    q = set(p.de_prob.p, p.param_axis, t)
    p.de_prob.f.jac((@view J[:, 1:end-1]), x, q, t)
    return (H, maybe_subtract!(J, Eye(size(J)...), p.de_prob))
end

"""
    BifurcationProblem(ode_or_map::AbstractODEProblem,
                       param_axis::Lens,
                       t_domain::Tuple;
                       <keyword arguments>)

# Arguments
- `ode_or_map`: An `ODEProblem` or `DiscreteProblem`.
- `param_axis :: Lens`: The lens to set/get a parameter of `ode_or_map`.
- `t_domain :: Tuple`: A pair of numbers specifying the lower and
  upper bound for `param_axis`.
"""
function BifurcationProblem(
        de_prob::DEP, param_axis::Lens, t_domain::Tuple;
        kwargs0...)
    de_prob = deepcopy(de_prob)
    u0 = de_prob.u0
    t0 = get(de_prob.p, param_axis)
    p0 = DiffEqWrapper(de_prob, param_axis)
    if de_prob.f.jac === nothing
        p = p0
        kwargs = kwargs0
    elseif de_prob isa DEP{true}
        p = @set p0.param_diff = MutableParamDiff(p0)
        kwargs = (; homotopy_jacobian=diffeq_homotopy_jacobian, kwargs0...)
    else
        # TODO: implement
        p = p0
        kwargs = kwargs0
    end
    return FixedPointBifurcationProblem(
        statekind(de_prob),
        timekind(de_prob),
        diffeq_homotopy, u0, t0,
        t_domain, p; kwargs...)
end


text:
@testset "Elliptic Curve Point Operations" begin
    @testset "Integer Type" begin
        @testset "Not Equal" begin
            a = ECC.Point(3, -7, 5, 7)
            b = ECC.Point(18, 77, 5, 7)
            @test a != b
            @test !(a != a)
        end
        @testset "On Curve?" begin
            @test_throws DomainError ECC.Point(-2, 4, 5, 7)
            @test typeof(ECC.Point(3, -7, 5, 7)) <: ECC.Point
            @test typeof(ECC.Point(18, 77, 5, 7)) <: ECC.Point
        end
        @testset "Addition" begin
            @testset "Base Case" begin
                a = ECC.Point(∞, ∞, 5, 7)
                b = ECC.Point(2, 5, 5, 7)
                c = ECC.Point(2, -5, 5, 7)
                @test a + b == b
                @test b + a == b
                @test b + c == a
            end

            @testset "Case 1" begin
                a = ECC.Point(3, 7, 5, 7)
                b = ECC.Point(-1, -1, 5, 7)
                @test a + b == ECC.Point(2, -5, 5, 7)
            end

            @testset "Case 2" begin
                a = ECC.Point(-1, 1, 5, 7)
                @test a + a == ECC.Point(18, -77, 5, 7)
            end
        end
    end;

    @testset "FiniteElement Type" begin
        @testset "On curve?" begin
            𝑝 = 223
            𝑎, 𝑏 = ECC.FieldElement(0, 𝑝), ECC.FieldElement(7, 𝑝)

            valid_points = ((192, 105), (17, 56), (1, 193))
            invalid_points = ((200, 119), (42, 99))

            for 𝑃 ∈ valid_points
                𝑥 = ECC.FieldElement(𝑃[1], 𝑝)
                𝑦 = ECC.FieldElement(𝑃[2], 𝑝)
                @test typeof(ECC.Point(𝑥, 𝑦, 𝑎, 𝑏)) <: ECC.Point
            end

            for 𝑃 ∈ invalid_points
                𝑥 = ECC.FieldElement(𝑃[1], 𝑝)
                𝑦 = ECC.FieldElement(𝑃[2], 𝑝)
                @test_throws DomainError ECC.Point(𝑥, 𝑦, 𝑎, 𝑏)
            end
        end
        @testset "Addition" begin
            𝑝 = 223
            𝑎 = ECC.FieldElement(0, 𝑝)
            𝑏 = ECC.FieldElement(7, 𝑝)

            additions = (
                (192, 105, 17, 56, 170, 142),
                (47, 71, 117, 141, 60, 139),
                (143, 98, 76, 66, 47, 71),
                )

            for 𝑛 ∈ additions
                𝑃 = ECC.Point(ECC.FieldElement(𝑛[1],𝑝),ECC.FieldElement(𝑛[2],𝑝),𝑎,𝑏)
                𝑄 = ECC.Point(ECC.FieldElement(𝑛[3],𝑝),ECC.FieldElement(𝑛[4],𝑝),𝑎,𝑏)
                𝑅 = ECC.Point(ECC.FieldElement(𝑛[5],𝑝),ECC.FieldElement(𝑛[6],𝑝),𝑎,𝑏)
                @test 𝑃 + 𝑄 == 𝑅
            end
        end
        @testset "Scalar Multiplication" begin
            𝑝 = 223
            𝑎 = ECC.FieldElement(0, 𝑝)
            𝑏 = ECC.FieldElement(7, 𝑝)

            multiplications = (
                (2, 192, 105, 49, 71),
                (2, 143, 98, 64, 168),
                (2, 47, 71, 36, 111),
                (4, 47, 71, 194, 51),
                (8, 47, 71, 116, 55),
                (21, 47, 71, ∞, ∞)
                )

            for 𝑛 ∈ multiplications
                λ = 𝑛[1]
                i = 2
                fieldelements = []
                while i < 6
                    if 𝑛[i] == ∞
                        push!(fieldelements, ∞)
                    else
                        push!(fieldelements, ECC.FieldElement(𝑛[i],𝑝))
                    end
                    i += 1
                end
                𝑃 = ECC.Point(fieldelements[1],fieldelements[2],𝑎,𝑏)
                𝑅 = ECC.Point(fieldelements[3],fieldelements[4],𝑎,𝑏)
                @test λ * 𝑃 == 𝑅
            end
        end
    end;
end


text:
# unary
@i @inline function NEG(a!::GVar)
    NEG(a!.x)
    NEG(a!.g)
end

function INV(x!::GVar{T}) where T
    x2 = x!.x ^ 2
    GVar(INV(x!.x), -x!.g * x2)
end

@i @inline function DEC(a!::GVar)
    DEC(a!.x)
end

# +-
@i @inline function :(-=)(identity)(a!::GVar, b::GVar)
    a!.x -= b.x
    b.g += a!.g
end

# inv
@eval @i @inline function :(-=)(inv)(out!::GVar{T}, y::GVar) where T
    out!.x -= inv(y.x)
    @routine @invcheckoff begin
        @zeros T a1
        a1 += y.x ^ 2
    end
    y.g -= out!.g / a1
    ~@routine
end

# +- (triple)
@i @inline function :(-=)(+)(out!::GVar, x::GVar, y::GVar)
    out!.x -= x.x + y.x
    x.g += out! |> grad
    y.g += out! |> grad
end

@i @inline function :(-=)(+)(out!::GVar, x::GVar, y::Real)
    out!.x -= (x |> value) + (y |> value)
    x.g += out! |> grad
end

@i @inline function :(-=)(+)(out!::GVar, x::Real, y::GVar)
    out!.x -= (x |> value) + (y |> value)
    y.g += out! |> grad
end

@i @inline function :(-=)(-)(out!::GVar, x::GVar, y::GVar)
    out!.x -= x.x - y.x
    x.g += out! |> grad
    y.g -= out! |> grad
end

@i @inline function :(-=)(-)(out!::GVar, x::Real, y::GVar)
    out!.x -= (x |> value) - y.x
    y.g -= out!.g
end

@i @inline function :(-=)(-)(out!::GVar, x::GVar, y::Real)
    out!.x -= x.x - (y |> value)
    x.g += out! |> grad
end

# NOTE: it will error on `SWAP(a!::GVar, b)` or `SWAP(a!, b:GVar)`
@i @inline function SWAP(a!::GVar, b!::GVar)
    SWAP(a! |> value,  b! |> value)
    SWAP(a!.g, b!.g)
end

# */
@i @inline function :(-=)(*)(out!::GVar, x::GVar, y::GVar)
    out!.x -= x.x * y.x
    x.g += out!.g * y.x
    y.g += x.x * out!.g
end

@i @inline function :(-=)(*)(out!::GVar, x::Real, y::GVar)
    out!.x -= (x |> value) * y.x
    y.g += (x |> value) * out!.g
end

@i @inline function :(-=)(*)(out!::GVar, x::GVar, y::Real)
    out!.x -= x.x * (y |> value)
    x.g += out!.g * (y |> value)
end

for DIV in [:/, :÷]
    @eval @i @inline function :(-=)($DIV)(out!::GVar{T}, x::GVar, y::GVar) where T
    out!.x -= $DIV(x.x, y.x)
    @routine @invcheckoff begin
        a1 ← zero(out! |> grad)
        a2 ← zero(out! |> grad)
        a1 += x.x * out!.g
        a2 += $DIV(a1, y.x)
    end
    x.g += $DIV(out!.g, y.x)
    y.g -= $DIV(a2, y.x)
    ~@routine
end

@eval @i @inline function :(-=)($DIV)(out!::GVar{T}, x::Real, y::GVar) where T
    out!.x -= $DIV(x, y.x)
    @routine @invcheckoff begin
        a1 ← zero(out!.g)
        a2 ← zero(out!.g)
        a1 += x * out!.g
        a2 += $DIV(a1, y.x)
    end
    y.g -= $DIV(a2, y.x)
    ~@routine
end

@eval @i @inline function :(-=)($DIV)(out!::GVar, x::GVar, y::Real)
    out!.x -= $DIV(x.x, y)
    x.g += $DIV(out!.g, y)
end
end

@i @inline function :(-=)(^)(out!::GVar{T}, x::GVar, n::GVar) where T
    # grad x
    @routine @invcheckoff begin
        @zeros T anc1 anc2 anc3 jac1 jac2 nx_1
        nx_1 += n.x - 1
        anc1 += x.x ^ nx_1
        jac1 += anc1 * n.x

        # get grad of n
        anc2 += log(x.x)
        anc3 += anc1 * x.x
        jac2 += anc3 * anc2
    end
    out!.x -= anc1 * x.x
    x.g += out!.g * jac1
    n.g += out!.g * jac2
    ~@routine
end

@i @inline function :(-=)(^)(out!::GVar{T}, x::GVar, n::Real) where T
    @routine @invcheckoff begin
        anc1 ← zero(x.x)
        jac ← zero(x.x)
	nx_1 ← zero(n)

	nx_1 += n - 1
        anc1 += x.x ^ nx_1
        jac += anc1 * n
    end
    out!.x -= anc1 * x.x
    x.g += out!.g * jac
    ~@routine
end

@i @inline function :(-=)(^)(out!::GVar{T}, x::Real, n::GVar) where T
    # get jac of n
    @routine @invcheckoff begin
        anc1 ← zero(x)
        anc2 ← zero(x)
        jac ← zero(x)

        anc1 += log(x)
        anc2 += x ^ n.x
        jac += anc1*anc2
    end
    out!.x -= anc2
    n.g += out!.g * jac
    ~@routine
end

for (OP, F) in [(:min, :<), (:max, :>)]
@eval @i @inline function :(-=)($OP)(out!::GVar{T}, x::GVar, y::GVar) where T
    if $F(x, y)
        out!.x -= x.x
        x.g += out!.g
    else
        out!.x -= y.x
        y.g += out!.g
    end
end

@eval @i @inline function :(-=)($OP)(out!::GVar{T}, x::GVar, y::Real) where T
    if $F(x, y)
        out!.x -= x.x
        x.g += out!.g
    else
        out!.x -= y.x
    end
end

@eval @i @inline function :(-=)($OP)(out!::GVar{T}, x::Real, y::GVar) where T
    if $F(x, y)
        out!.x -= x.x
    else
        out!.x -= y.x
        y.g += out!.g
    end
end
end

@i @inline function :(-=)(atan)(out!::GVar{T}, y::GVar, x::GVar) where T
    out!.x -= atan(y.x, x.x)
    @routine @invcheckoff begin
        @zeros T xy2 jac_x jac_y
        xy2 += abs2(x.x)
        xy2 += abs2(y.x)
        jac_y += x.x / xy2
        jac_x += (-y.x) / xy2
    end
    y.g += out!.g * jac_y
    x.g += out!.g * jac_x
    ~@routine
end
@i @inline function :(-=)(atan)(out!::GVar{T}, y::Real, x::GVar) where T
    out!.x -= atan(y, x.x)
    @routine @invcheckoff begin
        @zeros T xy2 jac_x
        xy2 += abs2(x.x)
        xy2 += abs2(y)
        jac_x += (-y) / xy2
    end
    x.g += out!.g * jac_x
    ~@routine
end
@i @inline function :(-=)(atan)(out!::GVar{T}, y::GVar, x::Real) where T
    out!.x -= atan(y.x, x)
    @routine @invcheckoff begin
        @zeros T xy2 jac_y
        xy2 += abs2(x)
        xy2 += abs2(y.x)
        jac_y += x / xy2
    end
    y.g += out!.g * jac_y
    ~@routine
end

@i @inline function :(-=)(atan)(out!::GVar{T}, x::GVar) where T
    out!.x -= atan(x.x)
    @routine @invcheckoff begin
        xy2 ← one(T)
        xy2 += abs2(x.x)
    end
    x.g += out!.g / xy2
    ~@routine
end

@i @inline function :(-=)(abs)(out!::GVar, x::GVar{T}) where T
    out!.x -= abs(x.x)
    if (x > 0, ~)
        x.g += out!.g
    else
        x.g -= out!.g
    end
end

@i @inline function :(-=)(abs2)(out!::GVar, x::GVar{T}) where T
    out!.x -= abs2(x.x)
    x.g += out!.g * x.x
    x.g += out!.g * x.x
end

for op in [:*, :/, :^, :+, :-, :atan, :max, :min]
    @eval @nograd :(-=)($op)(out!::GVar, x::Real, y::Real)
    @eval @nograd :(-=)($op)(out!::Real, x::Real, y::GVar)
    @eval @nograd :(-=)($op)(out!::Real, x::GVar, y::GVar)
    @eval @nograd :(-=)($op)(out!::Real, x::GVar, y::Real)
end

@i @inline function :(-=)(sqrt)(out!::GVar, x::GVar{T}) where T
    if x.x != 0
        @routine @invcheckoff begin
            @zeros T anc1 anc2
            anc1 += sqrt(x.x)
            anc2 += 2 * anc1
        end
        out!.x -= anc1
        x.g += out!.g / anc2
        ~@routine
    end
end

@i @inline function :(-=)(exp)(out!::GVar, x::GVar{T}) where T
    @routine @invcheckoff begin
        anc1 ← zero(T)
        anc1 += exp(x.x)
    end
    out!.x -= anc1
    x.g += out!.g * anc1
    ~@routine
end

@i @inline function :(-=)(log)(out!::GVar, x::GVar{T}) where T
    out!.x -= log(x.x)
    x.g += out!.g / x.x
end

@i @inline function :(-=)(sin)(out!::GVar, x::GVar{T}) where T
    @routine @invcheckoff begin
        @zeros T s c
        (s, c) += sincos(x.x)
    end
    out!.x -= s
    x.g += out!.g * c
    ~@routine
end

@i @inline function :(-=)(sinh)(out!::GVar, x::GVar{T}) where T
    out!.x -= sinh(x.x)
    @routine @invcheckoff begin
        anc1 ← zero(x.x)
        anc1 += cosh(x.x)
    end
    x.g += out!.g * anc1
    ~@routine
end

@i @inline function (:-=)(asin)(out!::GVar, x::GVar{T}) where T
    out!.x -= asin(x.x)
    @routine @invcheckoff begin
        @zeros T sqrt_1_x2 x2
        x2 += x.x ^ 2
        sqrt_1_x2 += sqrt(x2 |> NEG |> AddConst(1))
    end
    x.g += out!.g / sqrt_1_x2
    ~@routine
end

@i @inline function (:-=)(cos)(out!::GVar, x::GVar{T}) where T
    @routine @invcheckoff begin
        @zeros T s c
        (s, c) += sincos(x.x)
    end
    out!.x -= c
    x.g -= out!.g * s
    ~@routine
end

@i @inline function :(-=)(cosh)(out!::GVar, x::GVar{T}) where T
    out!.x -= cosh(x.x)
    @routine @invcheckoff begin
        anc1 ← zero(x.x)
        anc1 += sinh(x.x)
    end
    x.g += out!.g * anc1
    ~@routine
end

@i @inline function :(-=)(acos)(out!::GVar, x::GVar{T}) where T
    out!.x -= acos(x.x)
    @routine @invcheckoff begin
        @zeros T sqrt_1_x2 x2
        x2 += x.x ^ 2
        sqrt_1_x2 += sqrt(x2 |> NEG |> AddConst(1))
    end
    x.g -= out!.g / sqrt_1_x2
    ~@routine
end

@i @inline function :(-=)(tan)(out!::GVar, x::GVar{T}) where T
    @routine @invcheckoff begin
        anc1 ← zero(x.x)
        anc2 ← one(x.x)
        anc1 += tan(x.x)
        anc2 += anc1^2
    end
    out!.x -= anc1
    x.g += out!.g * anc2
    ~@routine
end

@i @inline function :(-=)(tanh)(out!::GVar, x::GVar{T}) where T
    @routine @invcheckoff begin
        anc1 ← zero(x.x)
        anc2 ← one(x.x)
        anc1 += tanh(x.x)
        anc2 -= anc1^2
    end
    out!.x -= anc1
    x.g += out!.g * anc2
    ~@routine
end

@i @inline function :(-=)(sincos)(out!::Tuple{T1,T1}, x::GVar{T}) where {T1<:GVar, T}
    @routine @invcheckoff begin
        s ← zero(T)
        c ← zero(T)
        (s, c) += sincos(x.x)
    end
    (out! .|> value) -= (s, c)
    x.g += (out!.:1 |> grad) * c
    x.g -= (out!.:2 |> grad) * s
    ~@routine
end

for op in [:sqrt, :exp, :log, :sin, :cos, :tanh, :abs, :abs2, :identity, :inv]
    @eval @nograd :(-=)($op)(out!::Real, x::GVar)
    @eval @nograd :(-=)($op)(out!::GVar, x::Real)
end

@nograd :(-=)(sincos)(out!::Tuple{<:Real,<:Real}, x::GVar)
@nograd :(-=)(sincos)(out!::Tuple{<:GVar,<:GVar}, x::Real)

@i @inline function IROT(a!::GVar, b!::GVar, θ::GVar)
    IROT(a!.x, b!.x, θ.x)
    NEG(θ |> value)
    θ.x -= π/2
    ROT(a!.g, b!.g, θ.x)
    θ.g += a!.x * a!.g
    θ.g += b!.x * b!.g
    θ.x += π/2
    NEG(θ |> value)
    ROT(a!.g, b!.g, π/2)
end

@i @inline function IROT(a!::GVar, b!::GVar, θ::Real)
    IROT(a!.x, b!.x, θ)
    NEG(θ)
    θ -= π/2
    ROT(a!.g, b!.g, θ)
    θ += π/2
    NEG(θ)
    ROT(a!.g, b!.g, π/2)
end

@nograd IROT(a!::Real, b!::Real, θ::GVar)

export primitive_grad
function primitive_grad end

@i @inline function (mf::MinusEq)(out!::GVar, args...; kwargs...)
    out!.x -= mf.f((args .|> value)...; kwargs...)
    (args .|> grad) .+= (@skip! ntuple(x->out!.g, length(args))) .* (@skip! primitive_grad(mf.f, (args .|> value)...; kwargs...))  # unsafe statement, error on recursive gradient
end

@i @inline function (mf::MinusEq)(out!::GVar, x::GVar; kwargs...)
    out!.x -= mf.f(x |> value; kwargs...)
    x.g += (@skip! out!.g) * (@skip! primitive_grad(mf.f, x.x; kwargs...))  # unsafe statement
end

@i @inline function :(-=)(convert)(out!::GVar{Tx, Tg}, y::GVar) where {Tx, Tg}
    out!.x -= convert(y.x)
    y.g += convert(out!.g)
end

@i @inline function HADAMARD(x::GVar, y::GVar)
    HADAMARD(x.x, y.x)
    HADAMARD(x.g, y.g)
end

@i @inline function (f::AddConst)(y::GVar)
    y.x += f.x
end

# more data views
for (DT, OP, NOP) in [(:AddConst, :+, :-, :add), (:SubConst, :-, :+)]
    @eval chfield(x::GVar, ac::$DT, xval::GVar) = GVar($NOP(xval.x, ac.x), xval.g)
end

#chfield(x::T, ::typeof(INV), xval::T) where T<:GVar = GVar(INV(xval.x), -xval.g*(xval.x^2))
#chfield(x::T, ::typeof(NEG), xval::T) where T<:GVar = GVar(-xval.x, -xval.g)

for F in [:INV, :NEG, :FLIP, :INC, :DEC]
    @eval NiLangCore.chfield(x::T, ::typeof($F), xval::T) where T<:GVar = (~$F)(xval)
end


text:
module LookupTables

mutable struct LookupTable2D
    xvar::StepRangeLen{Float64}
    yvar::StepRangeLen{Float64}

    xdict::Dict{Float64, Int64}
    ydict::Dict{Float64, Int64}

    vals::Array{Float64, 2}
end

function RangeToDict(z::StepRangeLen)
    zdict = Dict()
    for i in eachindex(z)
        zdict[z[i]] = i
    end
    return zdict
end

function LookupTable2D(xvar::StepRangeLen{Float64}, yvar::StepRangeLen{Float64}, vals::Array{Float64, 2})
    xdict = RangeToDict(xvar)
    ydict = RangeToDict(yvar)

    return LookupTable2D(xvar, yvar, xdict, ydict, vals)
end

function LookupTable2D(xn::Int64, yn::Int64)
    return LookupTable2D(range(0, stop=1, length=xn), range(0, stop=1, length=yn), zeros(xn, yn))
end

function lookup(ltable::LookupTable2D, x::Float64, y::Float64)
    xind = ltable.xdict[x]
    yind = ltable.ydict[y]

    return ltable.vals[xind, yind]
end

export LookupTable2D, lookup

end # module


text:
# This file is a part of ROOTFramework.jl, licensed under the MIT License (MIT).

using Cxx

cxxinclude("TBrowser.h")

export new_TBrowser

export TBrowserPtr
const TBrowserPtr = pcpp"TBrowser"

new_TBrowser() = icxx""" new TBrowser(); """


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ bb15be80-a6e8-11eb-0275-d1636e2a6445
begin
	using Plots
	using Printf
	using CSV
	using DataFrames
end

# ╔═╡ 4c484400-a6e8-11eb-343d-e50f20f9e747
md"""

## **_Blasius Flow_**

#### **Description:**

This notebook computes the Blasius Flow (two-dimensional flat plate flow)

		blasius(ηmax, N, itermax, ϵProfile, ϵBC)

If the arguments are missing, it will use the default values.

		blasius(ηmax=10, N=30, itermax=40, ϵProfile=1e-6, ϵBC=1e-6)

#### **Blasius Equation**
Flow over an infinite flat plate can be represented using the ordinary differential equation (ODE) below:

$$f''' + ½ f f'' = 0$$

Boundary Conditions:

$$f(0) = 0$$
$$f'(0) = 0$$
$$f'(∞) = 1$$

#### **Solution Method:**

Let's assume:

$$p = f(η)$$
$$h = f'(η)$$

Then, we obtain two ODEs equations and boundary conditions
ODEs:

$$h'' + p h' = 0$$
$$p' - h = 0$$

BCs:

$$p(0) = 0$$
$$h(0) = 0$$
$$h(∞) = 1$$

The equations given above are solved using Thomas algorithm.

Feel free to ask questions!

*Furkan Oz*
*Kursat Kara*

[foz@okstate.edu](foz@okstate.edu)
[kursat.kara@okstate.edu](kursat.kara@okstate.edu)

"""

# ╔═╡ c46503b0-a6e8-11eb-2f66-b7b20fe36bf8
"""
Computes the Blasius Flow (Two-dimensional flow over a flat plate)

    blasius(ηmax, N, itermax, ϵProfile, ϵBC)

If the arguments are missing, it will use the default values.

    blasius(ηmax=10, N=30, itermax=40, ϵProfile=1e-6, ϵBC=1e-6)

Furkan Oz,
[foz@okstate.edu](foz@okstate.edu),

"""
function blasius(ηmax=10, N=30, itermax=40, ϵProfile=1e-6, ϵBC=1e-6)

	Δη = ηmax/N
    Δη²= Δη^2;

    #itermax = 40
    #ϵProfile = 1e-6
    #ϵBC = 1e-6

    iter = 0
    errorProfile = 1.
    errorBC = 1.;

    # Initialization of the arrays
    A = zeros(N+1)
    B = zeros(N+1)
    C = zeros(N+1)
    D = zeros(N+1)

    G = zeros(N+1)
    H = zeros(N+1)

    p = zeros(N+1)
    h = zeros(N+1)
    η = zeros(N+1)

    η = [(i-1)*Δη for i=1:N+1];

    #BCs
    h[1]   = 0.0   # h(η=0) = 0
    h[N+1] = 1.0   # h(η=∞) = 1
    p[1]   = 0.0   # p(η=0) = 0

    G[1] = 0.0
    H[1] = 0.0;

    #Solution initialization

    # assume a linear profile for h(η)
    h = [(i-1)/N  for i=1:N+1]

    for i = 2:N+1
        p[i] = p[i-1] + (h[i] + h[i-1])*Δη/2
    end

    println("iter     error            convergence of h(η→∞)")
    println("-----------------------------------------------")

    #while ϵProfile<=errorProfile && ϵBC<=errorBC && iter<itermax
    while ϵProfile<=errorProfile && iter<itermax
        A = [ 2/Δη² + p[i]/(2*Δη) for i=1:N+1]
        B = [-4/Δη²  for i=1:N+1]
        C = [ 2/Δη² - p[i]/(2*Δη) for i=1:N+1]
        D = [0 for i=1:N+1]

        for i=2:N
            G[i] = - ( C[i]*G[i-1] + D[i] )/(B[i] + C[i] * H[i-1])
            H[i] = -                 A[i]  /(B[i] + C[i] * H[i-1])
        end

        hold = copy(h)

        for i=N:-1:2
            h[i] = G[i] + H[i] * h[i+1]
        end

        errorProfile = maximum(abs.(hold-h))

        errorBC = abs(h[N+1]-h[N])

        for i = 2:N
            p[i] = p[i-1] + (h[i] + h[i-1])*Δη/2
        end

        iter += 1

        #println("iter: $iter,  errorProfile: $errorProfile, errorBC: $errorBC")
        @printf("%4.4d %16.6e %16.6e \n", iter, errorProfile, errorBC)
    end

     if errorProfile<=ϵProfile
        println("")
        println("Solution converged!")
        println("The maximum change between consecutive profiles is less than the error criteria ϵProfile=$ϵProfile.")
     end

     if errorBC<=ϵBC
        println("")
        println("Solution for the boundary condition converged!")
        println("The difference between h(N) and h(N+1) is less than the error criteria ϵBC=$ϵBC.")
     end
    return η,h
end

# ╔═╡ 0a40028e-a6e9-11eb-0fa1-31a1e582b549
blasius();

# ╔═╡ 0d602e50-a6e9-11eb-3e38-9770e40c5747
ηtest, htest = blasius(10,40);

# ╔═╡ 18d364f2-a6e9-11eb-2177-0ba14d02a5b4
begin
	data = CSV.read("blasius-schlichting-dataset.csv",DataFrame)
	ηref = convert(Array,data[1:3:end,1])
	uref = convert(Array,data[1:3:end,2])
	plot_ref = plot([htest,uref],[ηtest*sqrt(1/2),ηref],
        title = "Blasius Similarity Profile",
        label = ["Thomas Algorithm" "Schlichting(2016)"],
        legend = :topleft,
        xlabel = "h = f '",
        ylabel = "\\eta",
		ylim = [0,5],
        linewidth = [2.0 0.000],
        linecolor = :auto,
        markershape = :auto,
		markeralpha = [1.0 0.6],
        markercolor = :auto,
	)
	#png(plot_ref, "blasiusprofile")
end

# ╔═╡ 38f42fd0-a6e9-11eb-26fc-8f7fbe77dbb0
begin
	η5, h5 = blasius(5);
	η10, h10 = blasius(10);
	η20, h20 = blasius(20);
end

# ╔═╡ 3ecfd20e-a6e9-11eb-02d9-9dbb4b44a2be
plot([h5, h10, h20], [η5, η10, η20],
        title = "Blasius Flow - Effect of Viscous Layer Height",
        label = ["\\etamax = 5" "\\etamax = 10" "\\etamax = 20"],
        legend = :topleft,
        xlabel = "h",
        ylabel = "\\eta",
        linewidth = 1,
        linecolor = :black,
        markershape = [:x :utriangle :dtriangle],
        markersize = [5 5 5],
        markeralpha = [0.9 0.6 0.6],
        markercolor = [:black :red :green]
    )

# ╔═╡ 46bd8b70-a6e9-11eb-226c-39e68021b054
begin
	ηN10, hN10 = blasius(10, 10);
	ηN20, hN20 = blasius(10, 20);
	ηN40, hN40 = blasius(10, 40);
	ηN80, hN80 = blasius(10, 80);
end

# ╔═╡ 4b765892-a6e9-11eb-0d51-f332409372a6
begin
	gr()
plot([hN10, hN20, hN40, hN80], [ηN10, ηN20, ηN40, ηN80],
        seriestype = :scatter,
        title = "Blasius Flow - Effect of Grid Resolution",
        label = ["N = 10" "N = 20" "N = 40" "N = 80"],
        legend = :topleft,
        xlabel = "h",
        ylabel = "\\eta",
        markershape = [:circle :utriangle :dtriangle :square],
        markersize = [10 8 6 3],
        markeralpha = [0.6 0.6 0.6 0.6],
        markercolor = [:purple :red :green :blue]
    )
end

# ╔═╡ Cell order:
# ╟─4c484400-a6e8-11eb-343d-e50f20f9e747
# ╠═bb15be80-a6e8-11eb-0275-d1636e2a6445
# ╠═c46503b0-a6e8-11eb-2f66-b7b20fe36bf8
# ╠═0a40028e-a6e9-11eb-0fa1-31a1e582b549
# ╠═0d602e50-a6e9-11eb-3e38-9770e40c5747
# ╠═18d364f2-a6e9-11eb-2177-0ba14d02a5b4
# ╠═38f42fd0-a6e9-11eb-26fc-8f7fbe77dbb0
# ╠═3ecfd20e-a6e9-11eb-02d9-9dbb4b44a2be
# ╠═46bd8b70-a6e9-11eb-226c-39e68021b054
# ╠═4b765892-a6e9-11eb-0d51-f332409372a6

max_stars_repo_name: frkanz/A-CFD-Tutorial-in-Julia

max_issues_repo_name: frkanz/A-CFD-Tutorial-in-Julia

max_forks_repo_name: frkanz/A-CFD-Tutorial-in-Julia


text:
const DATA_FILE = joinpath(@__DIR__, "dummy.csv")
const NUM_DIMENSIONS, NUM_OBSERVATIONS = size(load_data(DATA_FILE)[1])
const SOLVER = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
const INIT_STRAT = SimpleCombinedStrategy(RuleOfThumbScott(), BoundedTaxErrorEstimate(0.05, 0.02, 0.98))

experiment = Dict{Symbol, Any}(
    :hash => 1,
    :data_file => DATA_FILE,
    :data_set_name => "Dummy",
    :output_file => joinpath(@__DIR__, "example_query_synthesis.json"),
    :log_dir => joinpath(@__DIR__, "example_query_synthesis.log"),
    :split_strategy_name => "Sf",
    :initial_pool_strategy_name => "Pu",
    :model => Dict(:type => :SVDDneg,
                   :param => Dict{Symbol, Any}(),
                   :init_strategy => INIT_STRAT),
    :query_strategy => Dict(:type => :DecisionBoundaryQss,
                            :param => Dict{Symbol, Any}(
                                :optimizer => ParticleSwarmOptimization()
                            )),
    :split_strategy => OneClassActiveLearning.DataSplits(trues(NUM_OBSERVATIONS)),
    :oracle => QuerySynthesisOCCOracle(SVDDneg, INIT_STRAT, DATA_FILE, SOLVER),
    :param => Dict(:num_al_iterations => 10,
                   :solver => SOLVER,
                   :initial_pools => fill(:U, NUM_OBSERVATIONS),
                   :adjust_K => true,
                   :initial_pool_resample_version => 1))


text:
__precompile__()

module BioSequenceGraphs

using BioSequences

include("Nodes.jl")
include("Links.jl")
include("SequenceGraph.jl")
include("IO.jl")
include("GFA1/GFA1.jl")

end # module BioSequenceGraphs


text:
struct Circuit{T}
    netlist::Netlist
    conduction_matrix::AbstractMatrix{T}
    currents_matrix#::Union{AbstractVector{T},SparseVector{T,Int},SparseVector{Number,Int}}
    parameters
 end

 function Circuit(net::Netlist,s::Num)
    N = nnodes(net)
    M = spzeros(N,N)
    I = spzeros(N)
    Ps = [s]
    for component in net.components
        m,i,p = SymStamp(component,N,s)
        M += m
        I += i
        push!(Ps,p...)
    end

    Gm = M[2:end,2:end]
    Im = I[2:end]
    #println(Ps)
    Circuit(net,Gm,Im,Ps)
 end

 function Circuit(net::Netlist)
    N = nnodes(net)
    M = spzeros(N,N)
    I = spzeros(N)
    Ps = []
    for component in net.components
        m,i,p = SymStamp(component,N)
        M += m
        I += i
        push!(Ps,p...)
    end

    Gm = M[2:end,2:end]
    Im = I[2:end]
    #println(Ps)
    Circuit(net,Gm,Im,Ps)
end

function SymbolicMatrixes(ckt::Circuit)
    N = nnodes(ckt.netlist)-1
    Vvector = @variables V[1:N]
    Vvector  = Symbolics.scalarize(Vvector)
    M = copy(ckt.conduction_matrix)
    I = copy(ckt.currents_matrix)
    return Matrix(M),Vector(I),Vector(Vvector[1])
end

function SymbolicAuxEqs(ckt::Circuit,V)
    eq = []
    for comp in ckt.netlist.components
        dd = symSubs(comp,V)
        for d in dd
            push!(eq,first(d)~last(d))
        end
    end
    return eq
end

function CircuitFunction(ckt::Circuit)
    N = nnodes(ckt.netlist)-1 #exclude ground
    Vnames = map(x->x.name,ckt.netlist.nodes)[2:end]
    test = Symbolics.Symbol.(Vnames)
    Vvector = Vector(Symbolics.Variable.(test))
    M = copy(ckt.conduction_matrix)
    I = copy(ckt.currents_matrix)
    S = size(M)
    for comp in ckt.netlist.components
        for l in 1:S[1]
            I[l] = substitute(I[l],symSubs(comp,Vvector))
            for c in 1:S[2]
                M[l,c] = substitute(M[l,c],symSubs(comp,Vvector))
            end
        end
    end
    return Matrix(M),Vector(I),Vector(Vvector)
end

"""
 dc_op(ckt::Circuit; maxiter=1000,abs_tol=1e-6,rel_tol=1e-2)
 Calculates the DC Operating Point of a given circuit using NewtonRaphson method
"""
function dc_op(ckt::Circuit; maxiter=1000,abs_tol=1e-6,rel_tol=1e-2)
    Mm,Im,Vm = CircuitFunction(ckt)

    M_f = eval(build_function(Mm,Vm)[2])
    I_f = eval(build_function(Im,Vm)[2])

    v₀ = zeros(length(Vm))
    Mm = zeros(size(Mm))
    Im = zeros(size(Im))
    Base.invokelatest(M_f,Mm,v₀)
    Base.invokelatest(I_f,Im,v₀)
    v₁ = Mm\Im
    iter = 0
    abserror = 100
    relerror = 100
    while iter < maxiter && abserror > abs_tol && relerror > rel_tol
        v₀ = v₁
        Base.invokelatest(M_f,Mm,v₀)
        Base.invokelatest(I_f,Im,v₀)
        v₁ = Mm\Im
        iter += 1
        abserror = abs(maximum(v₁-v₀))
        relerror = abs(maximum((v₁-v₀)/v₁))
        if iter == maxiter
            @warn "Maximum iteration achieved"
        end
    end
    return v₁
end

struct CircuitProperties{T}
    ckt::Circuit
    Vi::Int
    Vo::Int
    Z::AbstractMatrix{T}
    Y::AbstractMatrix{T}
    T::AbstractMatrix{T}
    H::AbstractMatrix{T}
    VoltageGain::T
    TransconductanceGain::T
    Rin::T
    Rout::T
end

"""
CircuitProperties(ckt::Circuit,Vi::Int,Vo::Int)
Calculate The Circuit Parameters such as Voltage Gain, TransconductanceGain, R_in, R_out for a circuit with input node Vi and output node Vo
# Parameters
ckt::Circuit - Circuit to be anaylized
Vi::Int - Input Node number
Vo::Int - Output Node number
"""

function CircuitProperties(ckt::Circuit,Vi::Int,Vo::Int)
    Z = inv(Matrix(ckt.conduction_matrix))
    Z = Z[[Vi,Vo],[Vi,Vo]]
    Y = Z
    try
        Y = inv(Z)
    catch
        Y = Matrix(ckt.conduction_matrix)[[Vi,Vo],[Vi,Vo]]
    end
    # T Matrix
    y11 = Y[1,1]
    y12 = Y[1,2]
    y21 = Y[2,1]
    y22 = Y[2,2]
    A = -y22/y21
    B = (-1/y21)
    C = ((y12*y21-y11*y22)/y21)
    D = (-y11/y21)
    T = [A B;C D]
    #H Matrix

    A = 1/y11
    B = (-y12/y11)
    C = y21/y11
    D = ((y11*y22-y12*y21)/y11)
    H = [A B;C D]

    VoltageGain = simplify(1/T[1,1])
    TransconductanceGain = simplify(1/T[1,2])
    Rin = simplify(1/Y[1,1])
    Rout = simplify(Z[2,2])

    CircuitProperties(ckt,Vi,Vo,Z,Y,T,H,VoltageGain,TransconductanceGain,Rin,Rout)
end


text:
module TestTheories
using Test

using Catlab.Theories
using Catlab.Syntax

sexpr(expr::GATExpr) = sprint(show_sexpr, expr)
unicode(expr::GATExpr; all::Bool=false) =
  all ? sprint(show, MIME("text/plain"), expr) : sprint(show_unicode, expr)
latex(expr::GATExpr; all::Bool=false) =
  all ? sprint(show, MIME("text/latex"), expr) : sprint(show_latex, expr)

@testset "Categories" begin
  include("Category.jl")
end

@testset "MonoidalCategories" begin
  include("Monoidal.jl")
  include("MonoidalAdditive.jl")
end

@testset "HigherCategories" begin
  include("HigherCategory.jl")
end

@testset "Preorders" begin
  include("Preorders.jl")
end

@testset "Schemas" begin
  include("Schema.jl")
end

end


text:
haltonnumber = BlackBoxOptim.Utils.haltonnumber
haltonsequence = BlackBoxOptim.Utils.haltonsequence

# The Halton sequence with base 2 according to http://en.wikipedia.org/wiki/Halton_sequence
HaltonSequence2 = [
    1/2, 1/4, 3/4, 1/8, 5/8, 3/8, 7/8, 1/16, 9/16
]

# The Halton sequence with base 2 according to http://en.wikipedia.org/wiki/Halton_sequence
HaltonSequence3 = [
    1/3, 2/3, 1/9, 4/9, 7/9, 2/9, 5/9, 8/9, 1/27
]

@testset "Halton numbers and sequence" begin

    @testset "Halton numbers" begin

        map(1:length(HaltonSequence2)) do i
            @test isapprox(haltonnumber(2, i), HaltonSequence2[i])
        end

        map(1:length(HaltonSequence3)) do i
            @test isapprox(haltonnumber(3, i), HaltonSequence3[i])
        end

    end

    @testset "Halton sequences" begin

        hseq2 = haltonsequence(2, length(HaltonSequence2))
        @test length(hseq2) == length(HaltonSequence2)
        map(1:length(HaltonSequence2)) do i
            @test isapprox(hseq2[i], HaltonSequence2[i])
        end

        hseq3 = haltonsequence(3, length(HaltonSequence3))
        @test length(hseq3) == length(HaltonSequence3)
        map(1:length(HaltonSequence3)) do i
            @test isapprox(hseq3[i], HaltonSequence3[i])
        end

    end

end


text:
module NeuralNetwork

include("structs.jl")
include("regularization.jl")

using Flux
using Flux: onecold
using Flux: Optimiser
using Flux.Losses: logitcrossentropy
using LinearAlgebra
using MLDatasets
using ProgressMeter
using Statistics

export get_loss_and_accuracy, train

"""
	get_loss_and_accuracy(data_loader, model, spectral_decoupling; args...)

Calculates loss and accuracy of the given data for the input model.

# Examples
```julia-repl
julia> get_loss_and_accuracy(train_loader, model)
(0.1337, 1.0f0)
```
"""
function get_loss_and_accuracy(data_loader::Flux.Data.DataLoader, model, spectral_decoupling::Bool = false; args...)
	args = Args(; args...)
	@assert !(spectral_decoupling && isa(args.sd_λ, Float64)) "sd_λ must be specified if spectral decoupling is used"

	accuracy = 0.0f0
	loss = 0.0f0
	num = 0

	# iterate over batches
	for (x, y) in data_loader
		ŷ = model(y)
		loss += spectral_decoupling ?
			Regularization.spectral_decoupling(ŷ, y, args.sd_λ) :
			logitcrossentropy(ŷ, y, agg = sum)
		accuracy += sum(onecold(ŷ) .== onecold(y))
		num += size(x)[end]
	end

	# divide metrics cumulated over batches by number of total data points
	return spectral_decoupling ? loss : loss / num, accuracy / num
end

"""
	neural_network(500)

Creates a Flux neural network according to the topology used by Pezeshki et al.
That is a NN with 2 hidden layers, 500 units each and ReLU activation. The hidden
layer dimensions can be adjusted by the first and only parameter.
"""
function neural_network(hidden_dim::Int64 = 500)
	return Chain(
		Dense(2, hidden_dim, relu),
		Dense(hidden_dim, 2)
	)
end

"""
    train(train_loader, test_loader, optimizer, spectral_decoupling; args...)

Creates and trains a neural network as defined in [`neural_network`](@ref) using cross-entropy loss with
the specified optimizer and additionally spectral decoupling as a loss regularization.

# Examples
```julia-repl
julia> train(train_loader, test_loader, "SGD", false)
[ Info: opt = SGD, sd = false, lr = 0.01, bs = 50
Training... 100%|███████████████████████████████████████████████████████████████████████████| Time: 0:00:05
  epoch:     1000
  loss:      0.1337
  accuracy:  1.0
```
"""
function train(data_loader::Flux.Data.DataLoader, optimizer::String = "SGD", spectral_decoupling::Bool = false; args...)
	args = Args(; args...)

	@info "opt = $optimizer, sd = $spectral_decoupling, lr = $(args.learning_rate), bs = $(args.batchsize)" *
		(spectral_decoupling ? ", sd_λ = $(args.sd_λ)" : "") * (optimizer == "WD" ? ", wd_λ = $(args.wd_λ)" : "")

	model = neural_network()
	params = Flux.params(model)
	optimizer = get_optimizer(args.learning_rate, optimizer; args.wd_λ)

	# record loss and accuracy progression over epochs
	loss_prog, accuracy_prog = [], []
	# record weight evolution of learned features
	z1_prog, z2_prog = [], []

	# training
	prog = Progress(args.epochs, 0.25, "Training... ", 75)
	for epoch in 1:args.epochs
		for (x, y) in data_loader
			grads = Flux.gradient(params) do
				spectral_decoupling ?
					Regularization.spectral_decoupling(model(x), y, args.sd_λ) :
					logitcrossentropy(model(x), y)
			end

			Flux.Optimise.update!(optimizer, params, grads)
			# evaluate learned features
			# Y = Diagonal(y)
			# Φ = ? --> NTK regime, unable to obtain NTRF matrix to perform SVD
			# U = U of singular value decomposition

			# workaround
			# can only be done if hidden_dim == batchsize
			if (size(model.layers[1].W, 1) == size(x, 2))
				U = transpose(model.layers[1].W) .* model.layers[2].W
				zs = U .* model(x)
				push!(z1_prog, sum(abs.(zs[1,:])))
				push!(z2_prog, sum(abs.(zs[2,:])))
			end
		end

		# evaluate loss and accuracy
		loss, accuracy = get_loss_and_accuracy(data_loader, model)
		push!(loss_prog, loss)
		push!(accuracy_prog, accuracy)

		# show neat progress bar
		ProgressMeter.next!(prog; showvalues = [
			(:epoch, epoch),
			(:loss, loss),
			(:accuracy, accuracy)
			]
		)
	end

	return model, Dict("loss" => loss_prog, "accuracy" => accuracy_prog, "z1" => z1_prog, "z2" => z2_prog)
end

"""
	get_optimizer(learning_rate, optimizer, wd_λ)

Returns a Flux.Optimiser according to the one requested by acronym with fixed parameters except
for the given learning rate. Supported optimizers are:
	- Adaptive Moment Estimation (ADAM)
	- Gradient Descent (GD)
	- Stochastic Gradient Descent (SGD)
	- Weight Decay (WD)
"""
function get_optimizer(learning_rate::Float64, optimizer::String = "SGD"; args...)
	args = Args(; args...)
	optimizer = uppercase(optimizer)
	@assert optimizer in ["ADAM", "GD", "SGD", "WD"] "Requested optimizer '$optimizer' is not supported."
	if optimizer == "WD"
		@assert isa(args.wd_λ, Float64) "wd_λ must be specified if WD is used"
	end

	if optimizer == "ADAM"
		return ADAM(learning_rate)
	elseif optimizer == "GD"
		return Descent(learning_rate)
	elseif optimizer == "SGD"
		# stochastic gradient descent
		return Momentum(learning_rate, 0.9)
	elseif optimizer == "WD"
		return Optimiser(WeightDecay(args.wd_λ), Descent(learning_rate))
	end
end

end


text:
using CairoMakie
using CubedSphere

w(z) = (1 - cn(z, 1/2)) / sn(z, 1/2)

CubedSphere.visualize_conformal_mapping(w, x_min=-1.85, y_min=-1.85, x_max=1.85, y_max=1.85, n_lines=51, n_samples=101, filepath="square_to_disk_conformal_mapping.png")


text:
using Distributions
using JSON, ForwardDiff, Calculus, PDMats, Compat # test dependencies
using Base.Test

tests = [
    "types",
    "utils",
    "samplers",
    "categorical",
    "univariates",
    "continuous",
    "fit",
    "multinomial",
    "binomial",
    "poissonbinomial",
    "dirichlet",
    "dirichletmultinomial",
    "mvnormal",
    "mvlognormal",
    "mvtdist",
    "kolmogorov",
    "edgeworth",
    "matrix",
    "vonmisesfisher",
    "conversion",
    "mixture",
    "gradlogpdf",
    "truncate",
    "noncentralt",
    "locationscale",
    "quantile_newton",
]

print_with_color(:blue, "Running tests:\n")

if nworkers() > 1
    rmprocs(workers())
end

if Base.JLOptions().code_coverage == 1
    addprocs(Sys.CPU_CORES, exeflags = ["--code-coverage=user", "--inline=no", "--check-bounds=yes"])
else
    addprocs(Sys.CPU_CORES, exeflags = "--check-bounds=yes")
end

@everywhere using Distributions
@everywhere using JSON, ForwardDiff, Calculus, PDMats, Compat # test dependencies
@everywhere using Base.Test
@everywhere srand(345679)
res = pmap(tests) do t
    include(t*".jl")
    nothing
end

# print method ambiguities
println("Potentially stale exports: ")
display(Base.Test.detect_ambiguities(Distributions))
println()


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

@enum errorType begin
    unsupportedFMU
    unknownFMUType
    unknown
end

"""
Format the fmi2Status into a String
"""
function errorTypeString(type::errorType)
    fname = StackTraces.stacktrace()[3].func    # index 3 to step into calling function!

    if type == unsupportedFMU
        return "$fname() doesn't support FMUs with this version."
    elseif type == unknownFMUType
        return "Unknown FMU type in $fname(), is neigther CS nor ME."
    end

    "Unknwon Assertion in $fname()."
end

function assert(cond::Bool, type::errorType = unknwon)
    @assert cond [errorTypeString(type)]
end

function error(type::errorType = unknwon)
    @assert false [errorTypeString(type)]
end


text:
# This page was generated on DATEOFTODAY
using GeoPhyInv
using Statistics
#md using Plots; gr();
#md using ColorSchemes
using Test

# This tutorial demonstrates how an instance of `SeisForwExpt` can be used iteratively
# to perform forward modeling using `update!` for
# various bundles of medium parameters.

# ### Medium #1
medium = Medium(:elastic_homo2D, 5)

# ### Medium #2
medium_new = similar(medium)
update!(medium_new, [:vp, :rho, :vs], rectangle = [[-500, -500], [500, 500]], perc = 5.0) # perturb velocity

# ### AGeom
ageom = AGeom(medium.mgrid, :surf, SSrcs(3), Recs(100)); # surface seismic

# ### Plotting #1
#md p1=heatmap(medium, :rho, seriescolor=cgrad(:roma))
#md scatter!(ageom, SSrcs())
#md scatter!(ageom, Recs())
#md p2=heatmap(medium_new, :rho, seriescolor=cgrad(:roma))
#md scatter!(ageom, SSrcs())
#md scatter!(ageom, Recs())
#md plot(p1, p2, size=(800,300))

# ### SrcWav
tgrid = range(0.0, stop = 2.0, length = 2500) # generate a temporal grid
wav = ricker(15.0, tgrid, tpeak = 0.25); # Choose a source wavelet
srcwav = SrcWav(tgrid, ageom, [:vz]) # initialize
update!(srcwav, [:vz], wav) # distribute to all supersources

# ### SeisForwExpt
pa = SeisForwExpt(
    FdtdElastic(),
    medium = medium,
    ageom = ageom,
    srcwav = srcwav,
    tgrid = tgrid,
    rfields = [:vz],
    verbose = true,
);

# ### Modeling #1
@time update!(pa);
d1 = copy(pa[:data][:vz])
#md p1=heatmap(pa[:data], :vz, 99, 99, grid=true, legend=:none, seriescolor=cgrad(:seismic));

# ### Change medium in `pa` without memory allocation
update!(pa, medium_new)

# ### Modeling #2
@time update!(pa);
d2 = copy(pa[:data][:vz])
#md p2=heatmap(pa[:data], :vz, 99, 99, grid=true, legend=:none, seriescolor=cgrad(:seismic));

# Test
@test d1 ≠ d2

# ### Plotting #2
#md plot(p1,p2, size=(500, 300))


text:
# BSD 3-Clause License
#
# Copyright (c) 2017-2019, Battelle Memorial Institute; Lawrence Livermore National Security, LLC; Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

using BinaryBuilder

HELICS_VERSION = v"2.7.0"
HELICS_SHA = "ad005c0948ef4284417d429112772d0b63ebfbc62c9093c02ac10f4a333d70f4"

sources = [
    ArchiveSource("https://github.com/GMLC-TDC/HELICS/releases/download/v$HELICS_VERSION/Helics-v$HELICS_VERSION-source.tar.gz",
                  "$HELICS_SHA"),
]

script = raw"""
cd $WORKSPACE/srcdir

mkdir build
cd build
CMAKE_ARGS="-DCMAKE_FIND_ROOT_PATH="${prefix}" -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE="$CMAKE_TARGET_TOOLCHAIN" -DCMAKE_BUILD_TYPE=Release"
HELICS_ARGS="-DHELICS_BUILD_TESTS=OFF"
cmake ${CMAKE_ARGS} ${HELICS_ARGS} ..
make -j${nproc}
make install
if [[ "${target}" == *-mingw* ]]; then
    # Remove a broken link that we don't need anyway
    rm ${prefix}/bin/libzmq.dll.a
fi
"""

products = [
    LibraryProduct("libhelicsSharedLib", :libhelicsSharedLib),
]

platforms = expand_cxxstring_abis(supported_platforms(exclude = Sys.isfreebsd))

dependencies = [
    Dependency("ZeroMQ_jll"),
    BuildDependency("boost_jll"),
]

# Build 'em!
build_tarballs(
    ARGS,
    "HELICS",
    HELICS_VERSION,
    sources,
    script,
    platforms,
    products,
    dependencies,
    ; preferred_gcc_version=v"7",
)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
include(string(@__DIR__, "/PARAMETERS_arlett.jl"))
include(string(@__DIR__, "/locations_arlett.jl"))
include(string(@__DIR__, "/spawn_arlett.jl"))
include(string(@__DIR__, "/flow_arlett.jl"))
include(string(@__DIR__, "/arlett_sim.jl"))
include(string(@__DIR__, "/arlett_animation.jl"))
include(string(@__DIR__, "/arlett_noflow.jl"))

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
#############################################################################
#Copyright (c) 2016-2017 Charles Le Losq
#
#The MIT License (MIT)
#
#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 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.
#
# functions.jl contains several mathematic functions
#
#
#############################################################################

"""
	peakmeas(x::Array{Float64},y::Array{Float64};smoothing = "yes", filter = :SavitzkyGolay, M=5,N=2,ese_y=1.,y_smo_out=false)

The peakmeas function allows performing measurements of the position, width, intensity and centroïd of a dominant peak in a provided x-y signal.

It smooths the signal with a Savitzky-Golay filter prior to measuring the peak position, width and intensity. It is advised to check that the M and N values of the Savitzky-Golay filter are adequate for your problem before trusting the results from peakmeas. For that, just use the y_smo_out option.

half-width at half-maximum are calculated as the width of the peak at half its maximum intensity. This calculation is not affected by any asumption of peak symmetry (no fitting is done).

Centroïd is calculated as sum(y./sum(y).*x).

INPUTS:

	x: Array{Float64}, the x values;

	y: Array{Float64}, the y values.

OPTIONS:

	smoothing, String, triggers the smoothing of the spectrum if set to yes (default value);

	filter, Symbol, the filter that will be used. See the smooth function documentation;

	M=5, the M parameter for smoothing y with a Savitzky-Golay filter. See smooth function documentation;

	N=2, the M parameter for smoothing y with a Savitzky-Golay filter. See smooth function documentation;

	y_smo_out=false, the smoothed signal. Signal will be smoothed if set to true, using the SavitzkyGolayFilter function with the M and N values. y_smo output will also be provided.

OUTPUTS:

	intensity: Float64, the intensity of the peak;

	position: Float64, the position of the peak;

	hwhm: Float64, the half-width at half-maximum of the peak;

	centroïd: Float64, the centroïd of the peak;

	if y_smo_out is set to true, then another output is provided:

	y_smo: Array{Float64}, the smoothed y signal.
"""
function peakmeas(x::Array{Float64,1}, y::Array{Float64,1}; smoothing = "yes", method = "savgol", window_length=5, polyorder=2, ese_y=1., y_smo_out=false)
    ### PRELIMINARY CHECK: INCREASING SIGNAL
    if x[end,1] < x[1,1]
        x = flipdim(x,1)
        y = flipdim(y,1)
    end

	if smoothing == "yes"
    	y_smo = vec(smooth(x,y,method = method, window_length=window_length, polyorder=polyorder))
	else
		y_smo = collect(y)
	end

    x_maximum = x[y_smo .== maximum(y_smo)]

	if size(x_maximum,1) >= 2
		error("Something went wrong, is there two peaks with identical intensities in the signal? This function treats signal with one main peak...")
	end

    x_1 = x[x .<x_maximum];
	x_2 = x[x .>=x_maximum];
    y_first_portion = y_smo[x .<x_maximum];
	y_second_portion = y_smo[x .>=x_maximum];

    half_int = maximum(y_smo)/2
    idx_1 = findmin(abs.(y_first_portion.-half_int))
    idx_2 = findmin(abs.(y_second_portion.-half_int))

    hwhm = (x_2[idx_2[2]].-x_1[idx_1[2]])./2
    position = x_maximum[1]
    intensity = maximum(y_smo)
    centroid = sum(y_smo./sum(y_smo).*x)

    if y_smo_out == true
      return intensity, position, hwhm, centroid, y_smo
    elseif y_smo_out ==false
      return intensity, position, hwhm, centroid
    else
      error("Set y_smo_out to true or false.")
    end

end


text:
@testset "Function Call" begin
s = """
function f(x; y = 1, z::Int = 4)
    println(z)
    x + y
end

f(1; y = 3)
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f{T}(x::T, y::T)
    a = Array{T, 1}()
    append!(a, [x, y])
    a
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(::Type{Int}, x, y)
    x + y
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x, y, x)
    x + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E331

s = """
function f{Int64}(x::Int64, y::Int64)
    x + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E534
@test msgs[1].variable == "Int64"

s = """
function f{T<:Int64}(x::T, y::T)
    x + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E513
@test msgs[1].variable == "T <: Int64"

s = """
function f{Int<:Real}(x::Int, y::Int)
    x + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E536

s = """
function f(x, args...)
    x + length(args)
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x, args..., bogus...)
    x + length(args) + length(bogus)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E413

s = """
function f(x=1, y, args...)
    x + length(args) + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E411

s = """
function f(x, y; z, q=1)
    x + q + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E423

s = """
function f(x, y; args..., z=1)
    x + length(args) + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E412

s = """
function f(x, args...; kwargs...)
    x + length(args) + length(kwargs)
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x::Array{Number,1})
    length(x)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E533

s = """
function f(x::Dict{Symbol,Number})
    length(x)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E533

s = """
function f(x; y = 1, z::Int = 0.1)
    x + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E516

s = """
function f(x; y = 1, z::Int = error("You must provide z"))
    x + y + z
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
g(x) = -x
function f(arr::Array)
    map(g, arr)
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
g(; x=0, y=1) = x+y
function f()
    d = Dict{Symbol, Int}(:x=>6, :y=>4)
    g(; d...)
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f{T}(a::Vector{T})
    n = size(a)
    tmp  = Array(T, n)
    tmp2 = zeros(T, 1)
    tmp3 = zeros(1, 2, 3)
    @lintpragma("Info type a")
    @lintpragma("Info type n")
    @lintpragma("Info type tmp")
    @lintpragma("Info type T")
    @lintpragma("Info type tmp2")
    @lintpragma("Info type tmp3")
    tmp
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test msgs[2].code == :I271
if VERSION ≥ v"0.6-"
    @test msgs[3].code == :I271
end
@test msgs[4].code == :I271
if VERSION ≥ v"0.6-"
    @test msgs[5].code == :I271
end
@test msgs[6].code == :I271

s = """
function f(x)
    f = x
    f
end
"""
msgs = lintstr(s)
@test msgs[1].code == :W355
@test msgs[1].variable == "f"

s = """
function f(x)
    local f = x
    f
end
"""
msgs = lintstr(s)
@test_broken isempty(msgs)

s = """
function f(x)
    x.somefunc(1)
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x=1)
    @lintpragma("Info type x")
    return x
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271

s = """
function f(x::Int8=Int8(1))
    @lintpragma("Info type x")
    return x
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271

s = """
function f(c::Char)
    x = convert(Int, c)
    @lintpragma("Info type x")
    return x
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271

s = """
function f(args...; dict...)
    @lintpragma("Info type args")
    @lintpragma("Info type dict")
    length(args) + length(dict)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test msgs[2].code == :I271

s = """
function f(args::Float64...)
    @lintpragma("Info type args")
    length(args)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271

s = """
function f(args::Float64...)
    if args[1] == 1
        0
    else
        length(args)
    end
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x::Symbol)
    if x == "blah"
        0
    else
        1
    end
end
"""
msgs = lintstr(s)
@test msgs[1].code == :W542

s = """
function f(X::Int)
    Y = X + 1
    return Y
end
"""
msgs = lintstr(s)
@test isempty(msgs)

@test_broken isempty(lintstr("""
function f1(a::Float64, b::Float64)
  function f2(x::Float64, a=a)
    x + a^2
  end
  return f2(b)
end
"""))

s="""
function f1(a::Float64; b=string(a))
    string(a) * b
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
f(; x=1) = x
f(; 2)
"""
msgs = lintstr(s)
@test msgs[1].code == :E133

s = """
function () end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function foo end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
(==)(1, 1)
"""
msgs = lintstr(s)
@test isempty(msgs)

s="""
(==)(1)
"""
msgs = lintstr(s)
@test msgs[1].code == :E132

s = """
(==)((1, 1)...)
"""
msgs = lintstr(s)
@test isempty(msgs)
end


text:
module RocketFindIndexOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: find_index()" begin

    println("Testing: operator find_index()")

    run_testset([
        (
            source      = from(0:5) |> find_index(d -> d !== 0 && d % 2 == 0),
            values      = @ts([ 3, c ]),
            source_type = Int
        ),
        (
            source      = of(1) |> find_index(d -> d !== 0 && d % 2 == 0),
            values      = @ts(c),
            source_type = Int
        ),
        (
            source      = of("Something") |> find_index(d -> true),
            values      = @ts([ 1, c ]),
            source_type = Int
        ),
        (
            source      = completed() |> find_index(d -> d !== 0 && d % 2 == 0),
            values      = @ts(c),
            source_type = Int
        ),
        (
            source      = completed() |> find_index(d -> true),
            values      = @ts(c),
            source_type = Int
        ),
        (
            source      = never() |> find_index(d -> true),
            values      = @ts(),
            source_type = Int
        )
    ])

end

end


text:
using ExcelFiles
using ExcelReaders
using IteratorInterfaceExtensions
using TableTraits
using TableTraitsUtils
using Dates
using DataValues
using Test

@testset "ExcelFiles" begin

filename = normpath(dirname(pathof(ExcelReaders)), "..", "test", "TestData.xlsx")

efile = load(filename, "Sheet1")

@test isiterable(efile) == true

full_dfs = [create_columns_from_iterabletable(load(filename, "Sheet1!C3:O7")), create_columns_from_iterabletable(load(filename, "Sheet1"))]
for (df, names) in full_dfs
    @test length(df) == 13
    @test length(df[1]) == 4

    @test df[1] == [1., 1.5, 2., 2.5]
    @test df[2] == ["A", "BB", "CCC", "DDDD"]
    @test df[3] == [true, false, false, true]
    @test df[4] == [2, "EEEEE", false, 1.5]
    @test df[5] == [9., "III", NA, true]
    @test df[6] == [3., NA, 3.5, 4]
    @test df[7] == ["FF", NA, "GGG", "HHHH"]
    @test df[8] == [NA, true, NA, false]
    @test df[9] == [Date(2015,3,3), DateTime(2015,2,4,10,14), Date(1988,4,9), Dates.Time(15,2,0)]
    @test df[10] == [Date(1965,4,3), DateTime(1950,8,9,18,40), Dates.Time(19,0,0), NA]
    @test eltype(df[11]) == ExcelReaders.ExcelErrorCell
    @test df[12][1][] isa ExcelReaders.ExcelErrorCell
    @test df[12][2][] isa ExcelReaders.ExcelErrorCell
    @test df[12][3][] isa ExcelReaders.ExcelErrorCell
    @test df[12][4] == NA
    @test df[13] == [NA, 3.4, "HKEJW", NA]
end

df, names = create_columns_from_iterabletable(load(filename, "Sheet1!C4:O7", header=false))
@test names == [:x1,:x2,:x3,:x4,:x5,:x6,:x7,:x8,:x9,:x10,:x11,:x12,:x13]
@test length(df[1]) == 4
@test length(df) == 13
@test df[1] == [1., 1.5, 2., 2.5]
@test df[2] == ["A", "BB", "CCC", "DDDD"]
@test df[3] == [true, false, false, true]
@test df[4] == [2, "EEEEE", false, 1.5]
@test df[5] == [9., "III", NA, true]
@test df[6] == [3, NA, 3.5, 4]
@test df[7] == ["FF", NA, "GGG", "HHHH"]
@test df[8] == [NA, true, NA, false]
@test df[9] == [Date(2015, 3, 3), DateTime(2015, 2, 4, 10, 14), DateTime(1988, 4, 9), Dates.Time(15,2,0)]
@test df[10] == [Date(1965, 4, 3), DateTime(1950, 8, 9, 18, 40), Dates.Time(19,0,0), NA]
@test isa(df[11][1], ExcelReaders.ExcelErrorCell)
@test isa(df[11][2], ExcelReaders.ExcelErrorCell)
@test isa(df[11][3], ExcelReaders.ExcelErrorCell)
@test isa(df[11][4], ExcelReaders.ExcelErrorCell)
@test isa(df[12][1][], ExcelReaders.ExcelErrorCell)
@test isa(df[12][2][], ExcelReaders.ExcelErrorCell)
@test isa(df[12][3][], ExcelReaders.ExcelErrorCell)
@test DataValues.isna(df[12][4])
@test df[13] == [NA, 3.4, "HKEJW", NA]

good_colnames = [:c1, :c2, :c3, :c4, :c5, :c6, :c7, :c8, :c9, :c10, :c11, :c12, :c13]
df, names = create_columns_from_iterabletable(load(filename, "Sheet1!C4:O7", header=false, colnames=good_colnames))
@test names == good_colnames
@test length(df[1]) == 4
@test length(df) == 13
@test df[1] == [1., 1.5, 2., 2.5]
@test df[2] == ["A", "BB", "CCC", "DDDD"]
@test df[3] == [true, false, false, true]
@test df[4] == [2, "EEEEE", false, 1.5]
@test df[5] == [9., "III", NA, true]
@test df[6] == [3, NA, 3.5, 4]
@test df[7] == ["FF", NA, "GGG", "HHHH"]
@test df[8] == [NA, true, NA, false]
@test df[9] == [Date(2015, 3, 3), DateTime(2015, 2, 4, 10, 14), DateTime(1988, 4, 9), Dates.Time(15,2,0)]
@test df[10] == [Date(1965, 4, 3), DateTime(1950, 8, 9, 18, 40), Dates.Time(19,0,0), NA]
@test isa(df[11][1], ExcelReaders.ExcelErrorCell)
@test isa(df[11][2], ExcelReaders.ExcelErrorCell)
@test isa(df[11][3], ExcelReaders.ExcelErrorCell)
@test isa(df[11][4], ExcelReaders.ExcelErrorCell)
@test isa(df[12][1][], ExcelReaders.ExcelErrorCell)
@test isa(df[12][2][], ExcelReaders.ExcelErrorCell)
@test isa(df[12][3][], ExcelReaders.ExcelErrorCell)
@test DataValues.isna(df[12][4])
@test df[13] == [NA, 3.4, "HKEJW", NA]

# Too few colnames
@test_throws ErrorException create_columns_from_iterabletable(load(filename, "Sheet1!C4:O7", header=true, colnames=[:c1, :c2, :c3, :c4]))

# Test for constructing DataFrame with empty header cell
data, names = create_columns_from_iterabletable(load(filename, "Sheet2!C5:E7"))
@test names == [:Col1, :x1, :Col3]

end


text:
function mark_kross!(r::Robot)
    for side in (Nord,Sud,Ost,West)
        putmarkhole(r,side)
        backing(r, inverse(side))
    end
    putmarker!(r)
end

function putmarkhole(r::Robot,side::HorizonSide)
    while !isborder(r,side)
        move!(r,side)
        putmarker!(r)
    end
end

function backing(r::Robot,side::HorizonSide)
    while ismarker(r)
        move!(r,side)
    end
end

inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2, 4))

max_stars_repo_name: bubasooo/proggggggggg

max_issues_repo_name: bubasooo/proggggggggg

max_forks_repo_name: bubasooo/proggggggggg


text:
facts("DB") do
    const mock_db = DB(
        Dict{Char, Dict{AbstractString, Lemma}}(),
        Dict{Char, Dict{Int, Synset}}(),
		Dict{Tuple{Int,AbstractString}, AbstractString}()
    )

    context("path_to_data_file") do
        expected = joinpath("MockDB", "dict", "data.verb")
        @fact WordNet.path_to_data_file("MockDB", 'v') --> expected
    end

    context("path_to_index_file") do
        expected = joinpath("MockDB", "dict", "index.verb")
        @fact WordNet.path_to_index_file("MockDB", 'v') --> expected
    end

    context("show") do
        @fact io_to_string(io -> show(io, mock_db)) --> "WordNet.DB"
    end
end


text:
#
# runtests.jl --
#
# Exercises XPA communication in Julia.
#
#------------------------------------------------------------------------------
#
# This file is part of XPA.jl released under the MIT "expat" license.
# Copyright (C) 2016-2017, Éric Thiébaut (https://github.com/emmt/XPA.jl).
#
module XPATests

using XPA
import Base: RefValue

const VERBOSE = true

function sproc1(running::RefValue{Bool}, srv::XPA.Server, params::String,
                buf::XPA.SendBuffer)
    if running[]
        VERBOSE && println("send: $params")
        result = 42
        try
            XPA.store!(buf, result)
            return XPA.SUCCESS
        catch err
            error(srv, err)
            return XPA.FAILURE
        end
    end
end

function rproc1(running::RefValue{Bool}, srv::XPA.Server, params::String,
                buf::XPA.ReceiveBuffer)

    status = XPA.SUCCESS
    if running[]
        nbytes = sizeof(buf)
        VERBOSE && println("receive: $params [$nbytes byte(s)]")
        if params == "quit"
            running[] = false
        elseif params == "greetings"
            status = XPA.message(srv, "hello folks!")
        end
    end
    return status
end

function main1()
    running = Ref(true)
    server = XPA.Server("TEST", "test1", "help me!",
                        XPA.SendCallback(sproc1, running),
                        XPA.ReceiveCallback(rproc1, running))
    while running[]
        XPA.poll(-1, 1)
    end
    close(server)
end

function sproc2(::Nothing, srv::XPA.Server, params::String,
                buf::XPA.SendBuffer)
    VERBOSE && println("send: $params")
    result = 42
    try
        XPA.store!(buf, result)
        return XPA.SUCCESS
    catch err
        error(srv, err)
        return XPA.FAILURE
    end
end

function rproc2(::Nothing, srv::XPA.Server, params::String,
                buf::XPA.ReceiveBuffer)

    status = XPA.SUCCESS
    if params == "quit"
        close(srv)
    elseif params == "greetings"
        status = XPA.message(srv, "hello folks!")
    else
        nbytes = sizeof(buf)
        status = XPA.message(srv, "I received $nbytes from you, thanks!\n")
    end
    return status
end

function main2()
    server = XPA.Server("TEST", "test1", "help me!",
                        XPA.SendCallback(sproc2, nothing),
                        XPA.ReceiveCallback(rproc2, nothing))
    XPA.mainloop()
end

end