text:
struct OrbitPlot
    sol
    body_names
    dim
end
@recipe function f(p::OrbitPlot)
    sol = p.sol
    body_names = p.body_names
    dim = p.dim
    @assert dim ∈ (2, 3)
    N = length(sol.u[1].x[1].x[1])
    body_names = body_names==nothing ? ["orbit $i" for i in 1:N] : body_names
    ind = i -> i:N:length(sol.u[1].x[1])
    for i in 1:N
        @series begin
            vars=(ind(i)...,)
            label --> body_names[i]
            vars --> vars
            sol
        end
    end
end

function orbitplot(sol::DESolution;body_names=nothing,dim=3,kwargs...)
    RecipesBase.plot(OrbitPlot(sol,body_names,dim);kwargs...)
end

export orbitplot

function plot_orbits(sol;body_names=nothing,dim=3,kwargs...)
    @assert dim ∈ (2, 3)
    N = length(sol.u[1].x[1].x[1])
    body_names = body_names==nothing ? ["orbit $i" for i in 1:N] : body_names
    ind = i -> i:N:length(sol.u[1].x[1])
    p = plot(sol, vars=(ind(1)...,), lab=body_names[1], kwargs...)
    for i in 2:N
        plot!(p, sol, vars=(ind(i)...,), lab=body_names[i], kwargs...)
    end
    p
end


text:
#-----------------------------------------------------------------------
#       BEGIN   KerrSchild.jl
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
#   KERR-SCHILD FUNCTIONS
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    rsq( X::RealVec , a::Real )
The `rsq` function takes a point `X` in Cartesian Kerr-Schild
coordinates and calculates the value of ``r^2`` at that point in a
Kerr spacetime with rotation parameter `a`
"""
function rsq( X::RealVec , a::Real )
    tpfl=typeof(X[1])
    x = X[2]
    y = X[3]
    z = X[4]

    return ( -a^2 + x^2 + y^2 + z^2 + sqrt(4*(a^2)*(z^2)
             + (-a^2 + x^2 + y^2 + z^2)^2))/2
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    fks( X::RealVec , a::Real , GM::Real=1 )
The `fks` function takes a point `X` in Cartesian Kerr-Schild
coordinates and calculates the value of function ``f`` at that point in
a Kerr spacetime with rotation parameter `a` and with the product of
gravitational constant and mass `GM`
"""
function fks( X::RealVec , a::Real , GM::Real=1 )
    tpfl=typeof(X[1])
    x = X[2]
    y = X[3]
    z = X[4]

    rs = rsq(X,a)
    r = sqrt(rs)
    return 2*GM*rs*r/(rs^2+(a^2)*(z^2))
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    kks( X::RealVec , a::Real )
The `kks` function takes a point `X` in Cartesian Kerr-Schild
coordinates and calculates the tensor product ``k_\\mu k_\\nu`` at that
point in a Kerr spacetime with rotation parameter `a`
"""
function kks( X::RealVec , a::Real )
    tpfl=typeof(X[1])
    k = zeros(tpfl,4)
    kk = zeros(tpfl,4,4)
    x = X[2]
    y = X[3]
    z = X[4]

    rs = rsq(X,a)
    r = sqrt(rs)

    k[1] = one(tpfl)
    k[2] = (r*x+a*y)/(rs+a^2)
    k[3] = (r*y-a*x)/(rs+a^2)
    k[4] = z/r

    for i=1:4
        for j=1:4
            kk[i,j] = k[i]*k[j]
        end
    end
    return kk
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   KERR-SCHILD METRIC COMPONENTS
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    gks( X::RealVec , a::Real=0 , GM::Real=1 )
The `gks` function takes a point `X` in Cartesian Kerr-Schild
coordinates and calculates the components of the Kerr-Schild metric at
that point in a Kerr spacetime with rotation parameter `a` and with the
product of gravitational constant and mass `GM`
"""
function gks( X::RealVec , a::Real=0 , GM::Real=1 )
    tpfl=typeof(X[1])

    return ημν(tpfl) + fks(X,a,GM)*kks(X,a)
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    ge( X::RealVec )
The `ge` function takes a point `X` in Cartesian Kerr-Schild
coordinates and calculates the components of the Kerr-Schild metric at
that point in a Kerr spacetime with the product of gravitational
constant and mass ``GM=1`` and with rotation parameter ``a=738``
(roughly corresponding to the Earth's angular momentum)
"""
function ge( X::RealVec )
    tpfl=typeof(X[1])

    a = tpfl(738)   #   This is the roughly the rotation parameter
                    #   for Earth in units where M_Earth = 1
                    #   For reference, the avg radius of Earth is
                    #   1.437e9

    return gks(X,a,1)
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#       END     KerrSchild.jl
#-----------------------------------------------------------------------

ext: jl


text:
"""
# BasicMD.jl

A collection of basic routines for Molecular Dynamics simulations written in
Julia.  These include:
* Euler–Maruyama (EM)
* Random Walk Metropolis (RWM)
* Metropolis Adjusted Langevin (MALA)
* Brünger, Brooks and Karplus (BBK)
* ABOBA, BAOAB
* Verlet
* Hamiltonian/Hybrid Monte Carlo (HMC)

## REPL help
`?` followed by an sampler name (`?RWM`)to obtain infromation about individual
samplers.  Use `?` on `sample_trajectory!` and `sample_trajectory` to obtain
information about running the samplers.  Finally, `?` on `MDOptions` will provide
information about how to set additional options, such as the number of
information.

## Performing Sampling

Having constructed a `sampler` structure, samplers are executed with the the
commands
```
sample_trajectory!(X, sampler, options=opts);
```
for an in place transformation of `X` or
```
Xvals = sample_trajectory(X₀, sampler);
```
to obtain a full trajectory.  Metropolis methods will also return the running
acceptance rate.

"""
module BasicMD

using LinearAlgebra

include("types.jl")
include("constraints.jl")
include("sample.jl")
include("sample_observables.jl")
include("sample_constraints.jl")
include("utils.jl")
export sample_trajectory, sample_trajectory!,
    sample_observables, MDOptions,
    Constraints, trivial_constraint!
# RWM methods
include("metropolis/zeroth_order/rwm.jl")
export RWM
# MALA methods
include("metropolis/first_order/mala.jl")
export MALA
# EM methods
include("nonmetropolis/first_order/em.jl")
export EM
# BBK methods
include("nonmetropolis/second_order/bbk.jl")
export BBK
# ABOBA methods
include("nonmetropolis/second_order/aboba.jl")
export ABOBA
# BAOAB methods
include("nonmetropolis/second_order/baoab.jl")
export BAOAB
# HMC methods
include("metropolis/second_order/hmc.jl")
export HMC
# G-JF methods
include("nonmetropolis/second_order/g_jf.jl")
export GJF
# Verlet methods
include("nonmetropolis/second_order/verlet.jl")
export Verlet
end


text:
@inline function Phi_se(CellData,s,z,Def)
    """
    Solid-Electrolyte Potential Transfer Function
    # Add Ins and Outs
        # Cell Data
        # Frequency Vector
        # Discretisation Locations
        # Electrode Definition
    """

 if Def == "Pos"
    Electrode = CellData.Pos #Electrode Length
 else
    Electrode = CellData.Neg #Electrode Length
 end

ζ = (1-CellData.Const.t_plus)/F    #Simplifying Variable
Rs = Electrode.Rs       # Particle radius [m]
Ds = Electrode.Ds       # Solid diffusivity [m^2/s]
CC_A = CellData.Const.CC_A   # Current-collector area [m^2]
as = 3*Electrode.ϵ_s/Rs # Specific interfacial surf. area
κ_eff = CellData.Const.κ*Electrode.ϵ_e^Electrode.κ_brug #Effective Electrolyte Conductivity
σ_eff = Electrode.σ*Electrode.ϵ_s^Electrode.σ_brug #Effective Electrode Conductivity

#Defining SOC
θ = CellData.Const.SOC * (Electrode.θ_100-Electrode.θ_0) + Electrode.θ_0

#Beta's
β = @. Rs*sqrt(s/Ds)

#Prepare for j0
ce0 = CellData.Const.ce0
cs_max = Electrode.cs_max
cs0 = cs_max * θ
α = Electrode.α

#Current Flux Density
if CellData.Const.CellTyp == "Doyle_94"
   κ = Electrode.k_norm/Electrode.cs_max/ce0^(1-α)
   j0 = κ*(ce0*(cs_max-cs0))^(1-α)*cs0^α
else
   j0 = Electrode.k_norm*(ce0*(cs0/cs_max*(1-cs0/cs_max)))^(1-α)
end

#Resistance
Rtot = R*CellData.Const.T/(j0*F^2) + Electrode.RFilm

∂Uocp_elc = CellData.Const.∂Uocp(Def,θ)/cs_max #Open Circuit Potential Partial
res0 = @. -3*∂Uocp_elc/(as*F*Electrode.L*CC_A*Rs) # residual for pole removal

ν = @. Electrode.L*sqrt((as/σ_eff+as/κ_eff)/(Rtot+∂Uocp_elc*(Rs/(F*Ds))*(tanh(β)/(tanh(β)-β)))) #Condensing Variable - eq. 4.13
ν_∞ = @. Electrode.L*sqrt(as*((1/κ_eff)+(1/σ_eff))/(Rtot))

ϕ_tf = @. Electrode.L/(CC_A*ν*sinh(ν))*((1/κ_eff)*cosh(ν*z)+(1/σ_eff)*cosh(ν*(z-1))) #Transfer Function - eq. 4.14
ϕ_tf = @. ϕ_tf - res0./s

zero_tf = @. (6*(5*Ds*F*Rtot-∂Uocp_elc*Rs)*σ_eff)/(30*CC_A*as*Ds*F*σ_eff*Electrode.L) + (5*as*Ds*F*Electrode.L^2*(σ_eff*(-1+3*z^2)+κ_eff*(2-6*z+3*z^2)))/(30*CC_A*as*Ds*F*σ_eff*κ_eff*Electrode.L)
D = @. Electrode.L/(CC_A*ν_∞*sinh(ν_∞))*((1/κ_eff)*cosh(ν_∞*z)+(1/σ_eff)*cosh(ν_∞*(z-1))) # Contribution to D as G->∞
D_term = "@. $(Electrode.L)/($CC_A*$ν_∞*sinh($ν_∞))*((1/$κ_eff)*cosh($ν_∞*$z)+(1/$σ_eff)*cosh($ν_∞*($z-1)))"
ϕ_tf[:,findall(s.==0)] .= zero_tf[:,findall(s.==0)]

if Def == "Pos" #Double check this implementation
   ϕ_tf = -ϕ_tf
   D = -D
   D_term = "@. -$(Electrode.L)/($CC_A*$ν_∞*sinh($ν_∞))*((1/$κ_eff)*cosh($ν_∞*$z)+(1/$σ_eff)*cosh($ν_∞*($z-1)))"
   if Debug == 1
      println("ϕ_tf:Phi_se:Pos",ϕ_tf[:,1])
      println("D:Phi_se:Pos",D)
      println("z:Phi_se:Pos",z)
      println("zero_tf:Phi_se:Pos",zero_tf)
      println("ν_∞:Phi_se:Pos",ν_∞)
      println("j0:Phi_se:Pos",j0)
      println("Rtot:Phi_se:Pos",Rtot)
      println("σ_eff:Phi_se:Pos",σ_eff)
      #println("res0:Phi_se:Pos",res0[:,1])
      println("∂Uocp_elc:Phi_se:Pos",∂Uocp_elc)
      println("L:Phi_se:Pos",Electrode.L)
  end
end
res0 = zeros(length(z))
return ϕ_tf, D, res0, D_term

end

max_stars_repo_name: BradyPlanden/RealisationAlgorithms

max_issues_repo_name: BradyPlanden/RealisationAlgorithms

max_forks_repo_name: BradyPlanden/RealisationAlgorithms


text:
struct BudgetedLongestPathInstance{T}
  graph::AbstractGraph{T}
  rewards::Dict{Edge{T}, Float64}
  weights::Dict{Edge{T}, Int}
  src::T
  dst::T

  budget::Int # weights * solution >= budget
  max_weight::Int

  function BudgetedLongestPathInstance(graph::AbstractGraph{T}, rewards::Dict{Edge{T}, Float64},
                                       weights::Dict{Edge{T}, Int}, src::T, dst::T;
                                       budget::Union{Nothing, Int}=nothing,
                                       max_weight::Union{Nothing, Int}=nothing) where T
    # Error checking.
    if src == dst
      error("Source node is the same as destination node.")
    end

    if ! isnothing(budget) && budget < 0
      error("Budget is present and negative.")
    end

    if any(collect(values(weights)) .< 0)
      error("At least a weight is negative.")
    end

    if length(rewards) != length(weights)
      error("Not the same number of values and weights; these two vectors must have the same size.")
    end

    # Complete the optional parameters.
    if isnothing(max_weight)
      max_weight = maximum(collect(values(weights)))
    end

    if isnothing(budget)
      d = length(values) # Dimension of the problem.
      budget = d * max_weight
    end

    # Return a new instance.
    new{T}(graph, rewards, weights, src, dst, budget, max_weight)
  end
end

graph(i::BudgetedLongestPathInstance{T}) where T = i.graph
rewards(i::BudgetedLongestPathInstance{T}) where T = i.rewards
weights(i::BudgetedLongestPathInstance{T}) where T = i.weights
src(i::BudgetedLongestPathInstance{T}) where T = i.src
dst(i::BudgetedLongestPathInstance{T}) where T = i.dst
budget(i::BudgetedLongestPathInstance{T}) where T = i.budget
max_weight(i::BudgetedLongestPathInstance{T}) where T = i.max_weight

dimension(i::BudgetedLongestPathInstance{T}) where T = ne(graph(i))
reward(i::BudgetedLongestPathInstance{T}, u::T, v::T) where T = rewards(i)[Edge(u, v)]
weight(i::BudgetedLongestPathInstance{T}, u::T, v::T) where T = weights(i)[Edge(u, v)]

struct BudgetedLongestPathSolution{T}
  instance::BudgetedLongestPathInstance{T}
  path::Vector{Edge{T}}
  states::Dict{Tuple{T, Int}, Float64}
  solutions::Dict{Tuple{T, Int}, Vector{Edge{T}}}
end

function paths_all_budgets(s::BudgetedLongestPathSolution{T}, max_budget::Int) where T
  return Dict{Int, Vector{Edge{T}}}(
    budget => s.solutions[s.instance.dst, budget] for budget in 0:max_budget)
end

function paths_all_budgets_as_tuples(s::BudgetedLongestPathSolution{T}, max_budget::Int) where T
  return Dict{Int, Vector{Tuple{T, T}}}(
    budget => [(src(e), dst(e)) for e in s.solutions[s.instance.dst, budget]]
    for budget in 0:max_budget)
end

function budgeted_lp_dp(i::BudgetedLongestPathInstance{T}) where T
  V = Dict{Tuple{T, Int}, Float64}()
  S = Dict{Tuple{T, Int}, Vector{Edge{T}}}()

  # Initialise. For β = 0, this is exactly Bellman-Ford algorithm with costs
  # (instead of rewards). Otherwise, use the same initialisation as
  # Bellman-Ford.
  β0 = lp_dp(ElementaryPathInstance(graph(i), rewards(i), src(i), dst(i)))
  for v in vertices(graph(i))
    S[v, 0] = β0.solutions[v]
    V[v, 0] = length(S[v, 0]) == 0 ? 0 : sum(rewards(i)[e] for e in S[v, 0])
  end

  for β in 1:budget(i)
    for v in vertices(graph(i))
      V[v, β] = -Inf
      S[v, β] = Edge{T}[]
    end

    V[src(i), β] = 0.0
  end

  # Dynamic part.
  for β in 1:budget(i)
    while true
      changed = false

      for v in vertices(graph(i))
        for w in inneighbors(graph(i), v)
          # Compute the remaining part of the budget still to use.
          remaining_budget = max(0, β - weight(i, w, v))

          # If the explored subproblem has no solution, skip it.
          if V[w, remaining_budget] == -Inf
            continue
          end

          # If using the solution to the currently explored subproblem would
          # lead to a cycle, skip it.
          if any(src(e) == v for e in S[w, remaining_budget])
            continue
          end

          # Compute the amount of budget already used by the solution to the currently explored subproblem.
          used_budget = 0
          if length(S[w, remaining_budget]) > 0
            used_budget = sum(weight(i, src(e), dst(e)) for e in S[w, remaining_budget])
          end
          if used_budget < remaining_budget
            continue
          end

          # Compute the maximum: is passing through w advantageous?
          if V[w, remaining_budget] + reward(i, w, v) > V[v, β] && used_budget >= remaining_budget
            changed = true

            V[v, β] = V[w, remaining_budget] + reward(i, w, v)
            S[v, β] = vcat(S[w, remaining_budget], Edge(w, v))
          end
        end
      end

      if ! changed
        break
      end
    end
  end

  return BudgetedLongestPathSolution(i, S[dst(i), budget(i)], V, S)
end


text:
module OperationsTest

import Random: randstring
import LibGit2
using Test
using UUIDs

using Pkg
using Pkg.Types

import Random: randstring
import LibGit2

include("utils.jl")

const TEST_PKG = (name = "Example", uuid = UUID("7876af07-990d-54b4-ab0e-23690620f79a"))

import Pkg.Types: semver_spec, VersionSpec
@testset "semver notation" begin
    @test semver_spec("^1.2.3") == VersionSpec("1.2.3-1")
    @test semver_spec("^1.2")   == VersionSpec("1.2.0-1")
    @test semver_spec("^1")     == VersionSpec("1.0.0-1")
    @test semver_spec("^0.2.3") == VersionSpec("0.2.3-0.2")
    @test semver_spec("^0.0.3") == VersionSpec("0.0.3-0.0.3")
    @test semver_spec("^0.0")   == VersionSpec("0.0.0-0.0")
    @test semver_spec("^0")     == VersionSpec("0.0.0-0")
    @test semver_spec("~1.2.3") == VersionSpec("1.2.3-1.2")
    @test semver_spec("~1.2")   == VersionSpec("1.2.0-1.2")
    @test semver_spec("~1")     == VersionSpec("1.0.0-1")
    @test semver_spec("1.2.3")  == semver_spec("^1.2.3")
    @test semver_spec("1.2")    == semver_spec("^1.2")
    @test semver_spec("1")      == semver_spec("^1")
    @test semver_spec("0.0.3")  == semver_spec("^0.0.3")
    @test semver_spec("0")      == semver_spec("^0")

    @test semver_spec("0.0.3, 1.2") == VersionSpec(["0.0.3-0.0.3", "1.2.0-1"])
    @test semver_spec("~1.2.3, ~v1") == VersionSpec(["1.2.3-1.2", "1.0.0-1"])

    @test   v"1.5.2"  in semver_spec("1.2.3")
    @test   v"1.2.3"  in semver_spec("1.2.3")
    @test !(v"2.0.0"  in semver_spec("1.2.3"))
    @test !(v"1.2.2"  in semver_spec("1.2.3"))
    @test   v"1.2.99" in semver_spec("~1.2.3")
    @test   v"1.2.3"  in semver_spec("~1.2.3")
    @test !(v"1.3"    in semver_spec("~1.2.3"))
    @test  v"1.2.0"   in semver_spec("1.2")
    @test  v"1.9.9"   in semver_spec("1.2")
    @test !(v"2.0.0"  in semver_spec("1.2"))
    @test !(v"1.1.9"  in semver_spec("1.2"))
    @test   v"0.2.3"  in semver_spec("0.2.3")
    @test !(v"0.3.0"  in semver_spec("0.2.3"))
    @test !(v"0.2.2"  in semver_spec("0.2.3"))
    @test   v"0.0.0"  in semver_spec("0")
    @test  v"0.99.0"  in semver_spec("0")
    @test !(v"1.0.0"  in semver_spec("0"))
    @test  v"0.0.0"   in semver_spec("0.0")
    @test  v"0.0.99"  in semver_spec("0.0")
    @test !(v"0.1.0"  in semver_spec("0.0"))

    @test semver_spec("<1.2.3") == VersionSpec("0.0.0 - 1.2.2")
    @test semver_spec("<1.2") == VersionSpec("0.0.0 - 1.1")
    @test semver_spec("<1") == VersionSpec("0.0.0 - 0")
    @test semver_spec("<2") == VersionSpec("0.0.0 - 1")
    @test semver_spec("<0.2.3") == VersionSpec("0.0.0 - 0.2.2")
    @test semver_spec("<2.0.3") == VersionSpec("0.0.0 - 2.0.2")
    @test   v"0.2.3" in semver_spec("<0.2.4")
    @test !(v"0.2.4" in semver_spec("<0.2.4"))

    @test semver_spec("=1.2.3") == VersionSpec("1.2.3")
    @test semver_spec("=1.2") == VersionSpec("1.2.0")
    @test semver_spec("  =1") == VersionSpec("1.0.0")
    @test   v"1.2.3" in semver_spec("=1.2.3")
    @test !(v"1.2.4" in semver_spec("=1.2.3"))
    @test !(v"1.2.2" in semver_spec("=1.2.3"))

    @test semver_spec("≥1.3.0") == semver_spec(">=1.3.0")

    @test semver_spec(">=   1.2.3") == VersionSpec("1.2.3-*")
    @test semver_spec(">=1.2  ") == VersionSpec("1.2.0-*")
    @test semver_spec("  >=  1") == VersionSpec("1.0.0-*")
    @test   v"1.0.0" in semver_spec(">=1")
    @test   v"0.0.1" in semver_spec(">=0")
    @test   v"1.2.3" in semver_spec(">=1.2.3")
    @test !(v"1.2.2" in semver_spec(">=1.2.3"))

    @test_throws ErrorException semver_spec("^^0.2.3")
    @test_throws ErrorException semver_spec("^^0.2.3.4")
    @test_throws ErrorException semver_spec("0.0.0")
end

# TODO: Should rewrite these tests not to rely on internals like field names
@testset "union, isjoinable" begin
    @test sprint(print, VersionRange("0-0.3.2")) == "0-0.3.2"
    # test missing paths on union! and isjoinable
    # there's no == for VersionBound or VersionRange
    unified_vr = union!([VersionRange("1.5-2.8"), VersionRange("2.5-3")])[1]
    @test unified_vr.lower.t == (UInt32(1), UInt32(5), UInt32(0))
    @test unified_vr.upper.t == (UInt32(3), UInt32(0), UInt32(0))
    unified_vr = union!([VersionRange("2.5-3"), VersionRange("1.5-2.8")])[1]
    @test unified_vr.lower.t == (UInt32(1), UInt32(5), UInt32(0))
    @test unified_vr.upper.t == (UInt32(3), UInt32(0), UInt32(0))
    unified_vr = union!([VersionRange("1.5-2.2"), VersionRange("2.5-3")])[1]
    @test unified_vr.lower.t == (UInt32(1), UInt32(5), UInt32(0))
    @test unified_vr.upper.t == (UInt32(2), UInt32(2), UInt32(0))
    unified_vr = union!([VersionRange("1.5-2.2"), VersionRange("2.5-3")])[2]
    @test unified_vr.lower.t == (UInt32(2), UInt32(5), UInt32(0))
    @test unified_vr.upper.t == (UInt32(3), UInt32(0), UInt32(0))
    unified_vb = Types.VersionBound(union!([v"1.5", v"1.6"])[1])
    @test unified_vb.t == (UInt32(1), UInt32(5), UInt32(0))
    unified_vb = Types.VersionBound(union!([v"1.5", v"1.6"])[2])
    @test unified_vb.t == (UInt32(1), UInt32(6), UInt32(0))
    unified_vb = Types.VersionBound(union!([v"1.5", v"1.5"])[1])
    @test unified_vb.t == (UInt32(1), UInt32(5), UInt32(0))
end

temp_pkg_dir() do project_path
    @testset "simple add and remove with preview" begin
        Pkg.init(project_path)
        Pkg.add(TEST_PKG.name; preview = true)
        @test !isinstalled(TEST_PKG)
        Pkg.add(TEST_PKG.name)
        @test isinstalled(TEST_PKG)
        @eval import $(Symbol(TEST_PKG.name))
        Pkg.rm(TEST_PKG.name; preview = true)
        @test isinstalled(TEST_PKG)
        Pkg.rm(TEST_PKG.name)
        @test !isinstalled(TEST_PKG)
    end

    @testset "package with wrong UUID" begin
        @test_throws ResolverError Pkg.add(PackageSpec(TEST_PKG.name, UUID(UInt128(1))))
    end

    @testset "adding and upgrading different versions" begin
        # VersionNumber
        Pkg.add(PackageSpec(TEST_PKG.name, v"0.3"))
        @test Pkg.installed()[TEST_PKG.name] == v"0.3"
        Pkg.add(PackageSpec(TEST_PKG.name, v"0.3.1"))
        @test Pkg.installed()[TEST_PKG.name] == v"0.3.1"
        Pkg.rm(TEST_PKG.name)

        # VersionRange
        Pkg.add(PackageSpec(TEST_PKG.name, VersionSpec(VersionRange("0.3.0-0.3.2"))))
        @test Pkg.installed()[TEST_PKG.name] == v"0.3.2"
        Pkg.up(; level = UpgradeLevel(:patch))
        @test Pkg.installed()[TEST_PKG.name] == v"0.3.3"
        Pkg.up(; level = UpgradeLevel(:minor))
        @test Pkg.installed()[TEST_PKG.name].minor != 3
        Pkg.rm(TEST_PKG.name)
    end

    @testset "testing" begin
        # TODO: Check that preview = true doesn't actually execute the test
        # TODO: Test-only dependencies
        Pkg.add(TEST_PKG.name)
        Pkg.test(TEST_PKG.name; coverage=true)
        pkgdir = Base.locate_package(Base.PkgId(TEST_PKG.uuid, TEST_PKG.name))
        # No coverage files being generated?
        @test_broken TEST_PKG.name * ".cov" in readdir(pkgdir)
        Pkg.rm(TEST_PKG.name)
    end

    @testset "pinning / freeing" begin
        Pkg.add(TEST_PKG.name)
        old_v = Pkg.installed()[TEST_PKG.name]
        Pkg.pin(PackageSpec(TEST_PKG.name, v"0.2"))
        @test Pkg.installed()[TEST_PKG.name].minor == 2
        Pkg.up(TEST_PKG.name)
        @test Pkg.installed()[TEST_PKG.name].minor == 2
        Pkg.free(TEST_PKG.name)
        Pkg.up()
        @test Pkg.installed()[TEST_PKG.name] == old_v
        Pkg.rm(TEST_PKG.name)
    end

    @testset "develop / freeing" begin
        Pkg.add(TEST_PKG.name)
        old_v = Pkg.installed()[TEST_PKG.name]
        Pkg.rm(TEST_PKG.name)
        mktempdir() do devdir
            withenv("JULIA_PKG_DEVDIR" => devdir) do
                Pkg.REPLMode.pkgstr("develop $(TEST_PKG.name)")
                @test isinstalled(TEST_PKG)
                @test Pkg.installed()[TEST_PKG.name] > old_v
                test_pkg_main_file = joinpath(devdir, TEST_PKG.name, "src", TEST_PKG.name * ".jl")
                @test isfile(test_pkg_main_file)
                # Pkg #152
                write(test_pkg_main_file,
                    """
                    module Example
                        export hello, domath
                        const example2path = joinpath(@__DIR__, "..", "deps", "deps.jl")
                        if !isfile(example2path)
                            error("Example is not installed correctly")
                        end
                        hello(who::String) = "Hello, \$who"
                        domath(x::Number) = x + 5
                    end
                    """)
                mkpath(joinpath(devdir, TEST_PKG.name, "deps"))
                write(joinpath(devdir, TEST_PKG.name, "deps", "build.jl"),
                    """
                    touch("deps.jl")
                    """
                )
                Pkg.build(TEST_PKG.name)
                @test isfile(joinpath(devdir, TEST_PKG.name, "deps", "deps.jl"))
                Pkg.test(TEST_PKG.name)
                Pkg.free(TEST_PKG.name)
                @test Pkg.installed()[TEST_PKG.name] == old_v
            end
        end
    end

    @testset "stdlibs as direct dependency" begin
        uuid_pkg = (name = "CRC32c", uuid = UUID("8bf52ea8-c179-5cab-976a-9e18b702a9bc"))
        Pkg.add("CRC32c")
        @test haskey(Pkg.installed(), uuid_pkg.name)
        Pkg.up()
        Pkg.test("CRC32c")
        Pkg.rm("CRC32c")
    end

    @testset "package name in resolver errors" begin
        try
            Pkg.add([PackageSpec(TEST_PKG.name, VersionSpec(v"55"))])
        catch e
            @test occursin(TEST_PKG.name, sprint(showerror, e))
        end
    end

    @testset "protocols" begin
        mktempdir() do devdir
            withenv("JULIA_PKG_DEVDIR" => devdir) do
                try
                    Pkg.setprotocol!("notarealprotocol")
                    # Pkg.develop is broken, update to use when fixed
                    @test_throws CommandError pkg"develop Example"
                    Pkg.setprotocol!()
                    pkg"develop Example"
                    @test isinstalled(TEST_PKG)
                finally
                    Pkg.setprotocol!()
                end
            end
        end
    end

    @testset "check logging" begin
        usage = Pkg.TOML.parse(String(read(joinpath(Pkg.logdir(), "manifest_usage.toml"))))
        @test any(x -> startswith(x, joinpath(project_path, "Manifest.toml")), keys(usage))
    end

    @testset "adding nonexisting packages" begin
        nonexisting_pkg = randstring(14)
        @test_throws CommandError Pkg.add(nonexisting_pkg)
        @test_throws CommandError Pkg.up(nonexisting_pkg)
    end

    Pkg.rm(TEST_PKG.name)

    @testset "legacy CI script" begin
        mktempdir() do dir
            LibGit2.with(LibGit2.clone("https://github.com/JuliaLang/Example.jl", joinpath(dir, "Example.jl"))) do r
                cd(joinpath(dir, "Example.jl")) do
                    let Pkg = Pkg
                        Pkg.clone(pwd())
                        Pkg.build("Example")
                        Pkg.test("Example"; coverage=true)
                        @test isfile(Pkg.dir("Example", "src", "Example.jl"))
                    end
                end
            end
        end
    end

    @testset "add julia" begin
        @test_throws CommandError Pkg.add("julia")
    end

    @testset "up in Project without manifest" begin
        mktempdir() do dir
            cp(joinpath(@__DIR__, "test_packages", "UnregisteredWithProject"), joinpath(dir, "UnregisteredWithProject"))
            cd(joinpath(dir, "UnregisteredWithProject")) do
               with_current_env() do
                    Pkg.up()
                    @test haskey(Pkg.installed(), "Example")
                end
            end
        end
    end

    @testset "failing building a package should throw" begin
        mktempdir() do path
            cd(path) do
                Pkg.generate("FailBuildPkg")
                cd("FailBuildPkg")
                with_current_env() do
                    write_build(pwd(), "error()")
                    @test_throws CommandError Pkg.build()
                end
            end
        end
    end
end

temp_pkg_dir() do project_path
    @testset "libgit2 downloads" begin
        Pkg.add(TEST_PKG.name; use_libgit2_for_all_downloads=true)
        @test haskey(Pkg.installed(), TEST_PKG.name)
        Pkg.rm(TEST_PKG.name)
    end
end

@testset "preview generate" begin
    mktempdir() do tmp
        cd(tmp) do
            Pkg.generate("Foo"; preview=true)
            @test !isdir(joinpath(tmp, "Foo"))
        end
    end
end

temp_pkg_dir() do project_path
    @testset "test should instantiate" begin
        mktempdir() do dir
            cp(joinpath(@__DIR__, "test_packages", "UnregisteredWithProject"), joinpath(dir, "UnregisteredWithProject"))
            cd(joinpath(dir, "UnregisteredWithProject")) do
                with_current_env() do
                    Pkg.add("Test") # test https://github.com/JuliaLang/Pkg.jl/issues/324
                    Pkg.test()
                end
            end
        end
    end
end

temp_pkg_dir() do project_path
    @testset "valid project file names" begin
        extract_uuid(toml_path) = begin
            uuid = ""
            for line in eachline(toml_path)
                m = match(r"uuid = \"(.+)\"", line)
                if m !== nothing
                    uuid = m.captures[1]
                    break
                end
            end
            return uuid
        end

        cd(project_path) do
            mktempdir() do tmp; cd(tmp) do
                pkg_name = "FooBar"
                # create a project and grab its uuid
                Pkg.generate(pkg_name)
                uuid = extract_uuid(joinpath(pkg_name, "Project.toml"))
                # activate project env
                Pkg.activate(abspath(pkg_name))
                # add an example project to populate manifest file
                Pkg.add("Example")
                Pkg.activate()
                # change away from default names
                mv(joinpath(pkg_name, "Project.toml"), joinpath(pkg_name, "JuliaProject.toml"))
                mv(joinpath(pkg_name, "Manifest.toml"), joinpath(pkg_name, "JuliaManifest.toml"))
                # make sure things still work
                Pkg.develop(abspath(pkg_name))
                @test isinstalled((name=pkg_name, uuid=UUID(uuid)))
                Pkg.rm(pkg_name)
                @test !isinstalled((name=pkg_name, uuid=UUID(uuid)))
            end end
        end # cd project_path
    end # @testset
end

include("repl.jl")

end # module


text:
"""
    rmvn(n[, μ], Σ)

Fast parrallel generation of multivariate normal samples.
"""
function rmvn end
for T in (Float64, Float32, Float16)
    @eval function rmvn(n::Int, μ::Vector{$T}, Σ::Matrix{$T})
        μ' .+ _rmvn(n, Σ)
    end
    @eval rmvn(n::Real, μ::Vector{$T}, Σ::Matrix{$T}) = rmvn(Int(n), μ, Σ)
    @eval function rmvn(n::Real, Σ::Matrix{$T})
        d = size(Σ, 2)
        rmvn(Int(n), zeros($T, d), Σ)
    end
end
function rmvn(n::Real, μ::Vector{<:Real}, Σ::Matrix{<:Real})
    rmvn(Int(n), Vector{Float64}(μ), Matrix{Float64}(Σ))
end


text:
"""
    gpsample([x, ]gp::FiniteGP; samples=1, orbit=0.0, kwargs...)

Plot sample(s) from the finite projection of the Gaussian process `gp` along `x`.

If `x` is not provided, it is set to `gp.x`.

The `orbit` keyword argument can be used to visualize a manifold of
similar samples.[^PH2013] Values are mapped to the interval ``[0, 1)`` which
correspond to angles in ``[0, 2\\pi)`` on the great circle of samples.

## Attributes

$(Makie.ATTRIBUTES)

[^PH2013]: Philipp Hennig (2013). [Animating Samples from Gaussian Distributions](http://mlss.tuebingen.mpg.de/2013/2013/Hennig_2013_Animating_Samples_from_Gaussian_Distributions.pdf). Technical Report No. 8 of the Max Planck Institute for Intelligent Systems.
"""
@recipe(GPSample, x, gp) do scene
    Attributes(;
        color=theme(scene, :linecolor),
        colormap=theme(scene, :colormap),
        colorrange=Makie.automatic,
        linestyle=theme(scene, :linestyle),
        linewidth=theme(scene, :linewidth),
        samples=Node(1),
        orbit=Node(0.0),
        cycle=[:color],
        inspectable=theme(scene, :inspectable),
    )
end

"""
    gpsample(x::AbstractVector, gp::AbstractGP; kwargs...)

Plot sample(s) from the finite projection `gp(x, 1e-9)` along `x`.
"""
gpsample(::AbstractVector, ::AbstractGP; kwargs...)

function Makie.plot!(plot::GPSample)
    @extract plot (x, gp)

    # Compute mean of the GP and the Cholesky decomposition of its covariance matrix
    mean_cov = lift(gp) do gp
        mean, Cmat = mean_and_cov(gp)
        C = cholesky(AbstractGPs._symmetric(Cmat))
        return mean, C
    end

    # Generate standard normally distributed samples and corresponding random tangents
    # We can then swipe a great circle of equivalent samples
    u_v_scale = lift(x, plot.samples) do x, nsamples
        # Generate samples and map them to the unit sphere
        u = randn(length(x), nsamples)
        scale = permutedims(map(norm, eachcol(u)))
        u ./= scale

        # Sample orthogonal tangent directions and map them to the unit sphere
        v = randn(size(u))
        for (ui, vi) in zip(eachcol(u), eachcol(v))
            vi .-= dot(ui, vi) .* ui
            vi ./= norm(vi)
        end

        return u, v, scale
    end

    orbit = plot.orbit
    points = lift(x, mean_cov, u_v_scale, orbit) do x, (mean, C), (u, v, scale), orbit
        # Ensure that `theta` is in [0, 2) (corresponding to angles in [0, 2π))
        theta = 2 * mod(orbit, 1)

        # Compute values on the orbit between `u` and `v`
        sin2pi_theta, cos2pi_theta = @static if VERSION < v"1.6.0-DEV.292"
            sinpi(theta), cospi(theta)
        else
            sincospi(theta)
        end
        R = scale .* (cos2pi_theta .* u .+ sin2pi_theta .* v)

        # Compute points to plot (NaN required to separate lines)
        d, n = size(R)
        z = fill(Point2f(NaN32), d + 1, n)
        zflat = reinterpret(Float32, z)

        # Fill x coordinates
        xcoords = view(zflat, 1:2:(2 * d), :)
        xcoords .= x

        # Compute y coordinates
        ycoords = view(zflat, 2:2:(2 * d), :)
        mul!(ycoords, C.U', R)
        ycoords .+= mean

        # Return vector of points
        return vec(z)
    end

    lines!(
        plot,
        points;
        color=plot.color,
        linestyle=plot.linestyle,
        linewidth=plot.linewidth,
        colormap=plot.colormap,
        colorrange=plot.colorrange,
        inspectable=plot.inspectable,
    )

    return plot
end


text:
function Base.promote_shape(a::Tuple{Vararg{Union{AbstractUnitRange,AbstractIndex}}}, b)
    return _promote_shape(a, b)
end
function Base.promote_shape(a, b::Tuple{Vararg{Union{AbstractUnitRange,AbstractIndex}}})
    return _promote_shape(a, b)
end
function Base.promote_shape(
    a::Tuple{Vararg{Union{AbstractUnitRange,AbstractIndex}}},
    b::Tuple{Vararg{Union{AbstractUnitRange,AbstractIndex}}}
   )
    return _promote_shape(a, b)
end

function _promote_shape(a, b)
    if length(a) < length(b)
        return promote_shape(b, a)
    end
    for i=1:length(b)
        if a[i] != b[i]
            throw(DimensionMismatch("dimensions must match"))
        end
    end
    for i=length(b)+1:length(a)
        if a[i] != 1:1
            throw(DimensionMismatch("dimensions must match"))
        end
    end
    return a
end


text:
using Bennu
using Test
using Pkg: Pkg

using Adapt
using Combinatorics
using CUDA
using CUDAKernels
using EzXML: EzXML
using FastGaussQuadrature: FastGaussQuadrature
using FillArrays
using KernelAbstractions
using LazyArrays
using LinearAlgebra
using Random
using SparseArrays
using StaticArrays
using StructArrays
using Tullio
using WriteVTK

CUDA.allowscalar(false)

include("arrays.jl")
include("cells.jl")
include("gridgenerators.jl")
include("grids.jl")
include("kernelabstractions.jl")
include("kroneckeroperators.jl")
include("metrics.jl")
include("operators.jl")
include("partitions.jl")
include("permutations.jl")
include("quadratures.jl")
include("sparsearrays.jl")
include("structarrays.jl")
include("tuples.jl")

@testset "examples" begin
    julia = Base.julia_cmd()
    base_dir = joinpath(@__DIR__, "..")

    for example_dir in readdir(joinpath(base_dir, "examples"), join=true)
        @testset "$example_dir" begin
            mktempdir() do tmp_dir
                # Change to temporary directory so that any files created by the
                # example get cleaned up after execution.
                cd(tmp_dir)
                example_project = Pkg.Types.projectfile_path(example_dir)
                tmp_project = Pkg.Types.projectfile_path(tmp_dir)
                cp(example_project, tmp_project)

                for script in filter!(s->endswith(s, ".jl"),
                                      readdir(example_dir, join=true))
                    cmd = `$julia --project=$tmp_project -e "import Pkg; Pkg.develop(path=raw\"$base_dir\"); Pkg.instantiate(); include(raw\"$script\")"`
                    @test success(pipeline(cmd, stderr=stderr, stdout=stdout))
                end
            end
        end
    end
end


text:
using FluxExtra, Documenter

makedocs(modules=[FluxExtra,FluxExtra.Normalizations],
    sitename = "FluxExtra.jl",
    pages = ["Home" => "index.md",
            "Layers" => "layers.md",
            "Functions" => "functions.md"],
    authors = "Open Machine Learning Association",
    format = Documenter.HTML(prettyurls = false)
)

deploydocs(
    repo = "github.com/OML-NPA/FluxExtra.jl.git",
    devbranch = "main"
)


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

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

"""

Create or update a dedicated host group. For details of Dedicated Host and Dedicated Host Groups please see [Dedicated Host Documentation] (https://go.microsoft.com/fwlink/?linkid=2082596)
Param: resourceGroupName::String (required)
Param: hostGroupName::String (required)
Param: parameters::DedicatedHostGroup (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DedicatedHostGroup
"""
function _swaggerinternal_dedicatedHostGroupsCreateOrUpdate(_api::DedicatedHostGroupsApi, resourceGroupName::String, hostGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", DedicatedHostGroup, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/hostGroups/{hostGroupName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "hostGroupName", hostGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function dedicatedHostGroupsCreateOrUpdate(_api::DedicatedHostGroupsApi, resourceGroupName::String, hostGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostGroupsCreateOrUpdate(_api, resourceGroupName, hostGroupName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function dedicatedHostGroupsCreateOrUpdate(_api::DedicatedHostGroupsApi, response_stream::Channel, resourceGroupName::String, hostGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostGroupsCreateOrUpdate(_api, resourceGroupName, hostGroupName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Delete a dedicated host group.
Param: resourceGroupName::String (required)
Param: hostGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_dedicatedHostGroupsDelete(_api::DedicatedHostGroupsApi, resourceGroupName::String, hostGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/hostGroups/{hostGroupName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "hostGroupName", hostGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function dedicatedHostGroupsDelete(_api::DedicatedHostGroupsApi, resourceGroupName::String, hostGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostGroupsDelete(_api, resourceGroupName, hostGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function dedicatedHostGroupsDelete(_api::DedicatedHostGroupsApi, response_stream::Channel, resourceGroupName::String, hostGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostGroupsDelete(_api, resourceGroupName, hostGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Retrieves information about a dedicated host group.
Param: resourceGroupName::String (required)
Param: hostGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: DedicatedHostGroup
"""
function _swaggerinternal_dedicatedHostGroupsGet(_api::DedicatedHostGroupsApi, resourceGroupName::String, hostGroupName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DedicatedHostGroup, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/hostGroups/{hostGroupName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "hostGroupName", hostGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function dedicatedHostGroupsGet(_api::DedicatedHostGroupsApi, resourceGroupName::String, hostGroupName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostGroupsGet(_api, resourceGroupName, hostGroupName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function dedicatedHostGroupsGet(_api::DedicatedHostGroupsApi, response_stream::Channel, resourceGroupName::String, hostGroupName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostGroupsGet(_api, resourceGroupName, hostGroupName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists all of the dedicated host groups in the specified resource group. Use the nextLink property in the response to get the next page of dedicated host groups.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DedicatedHostGroupListResult
"""
function _swaggerinternal_dedicatedHostGroupsListByResourceGroup(_api::DedicatedHostGroupsApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DedicatedHostGroupListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/hostGroups", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function dedicatedHostGroupsListByResourceGroup(_api::DedicatedHostGroupsApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostGroupsListByResourceGroup(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function dedicatedHostGroupsListByResourceGroup(_api::DedicatedHostGroupsApi, response_stream::Channel, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostGroupsListByResourceGroup(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists all of the dedicated host groups in the subscription. Use the nextLink property in the response to get the next page of dedicated host groups.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DedicatedHostGroupListResult
"""
function _swaggerinternal_dedicatedHostGroupsListBySubscription(_api::DedicatedHostGroupsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DedicatedHostGroupListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/hostGroups", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function dedicatedHostGroupsListBySubscription(_api::DedicatedHostGroupsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostGroupsListBySubscription(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function dedicatedHostGroupsListBySubscription(_api::DedicatedHostGroupsApi, response_stream::Channel, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostGroupsListBySubscription(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Update an dedicated host group.
Param: resourceGroupName::String (required)
Param: hostGroupName::String (required)
Param: parameters::DedicatedHostGroupUpdate (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DedicatedHostGroup
"""
function _swaggerinternal_dedicatedHostGroupsUpdate(_api::DedicatedHostGroupsApi, resourceGroupName::String, hostGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", DedicatedHostGroup, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/hostGroups/{hostGroupName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "hostGroupName", hostGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function dedicatedHostGroupsUpdate(_api::DedicatedHostGroupsApi, resourceGroupName::String, hostGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostGroupsUpdate(_api, resourceGroupName, hostGroupName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function dedicatedHostGroupsUpdate(_api::DedicatedHostGroupsApi, response_stream::Channel, resourceGroupName::String, hostGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostGroupsUpdate(_api, resourceGroupName, hostGroupName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export dedicatedHostGroupsCreateOrUpdate, dedicatedHostGroupsDelete, dedicatedHostGroupsGet, dedicatedHostGroupsListByResourceGroup, dedicatedHostGroupsListBySubscription, dedicatedHostGroupsUpdate


text:
len = length

using PyCall,DataFrames,RCall
unshift!(PyVector(pyimport("sys")["path"]), "")
@pyimport ncdata
filename = "dsmacc_out.nc"
filename = "../methane.nc"
data = ncdata.get(filename)

specs = names!(DataFrame(data["spec"]),[Symbol(i)for i in data["sc"]])
rates = names!(DataFrame(data["rate"]),[Symbol(i)for i in data["rc"]])

rates = rates[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["rc"])]]
specs = specs[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["sc"])]]

rates = rates[7:len(rates)]
rates = rates[:,[Symbol(i) for i in filter(x->sum(rates[x])>0,names(rates))]]

# reading in data and cleaning
joint = "\n"*(join(names(rates),"\n"))*"\n"

products= [split(i[4:len(i)],"+") for i in matchall(r"-->([A-z0-9+]*)",joint)]
reactants = [split(i[2:len(i)-1],"+") for i in matchall(r"\n([A-z0-9+]{1,60})([-->]{0,1})",joint)]

node_set = Set([ i[2:len(i)-1] for i in matchall(r"[\n+>-]{1}[A-z0-9]+[\n+>-]{1}",joint) ])
nodes = [i for i in node_set]
dict = Dict( [[i,x] for (x,i) in enumerate(nodes)])

#df = DataFrame(source = Char[], target = Char[], weight = Float16[])

flux = []
edges = []
for i in 1:length(reactants)
  rcol=[]
  #println(i)
  for j in reactants[i]
    coeff = match( r"([\d\s\.]*)(\D[\d\D]*)", j)
    dummy = specs[Symbol(coeff[2])]
    try
      push!(rcol,parse(Float32,coeff[1]))
    catch
      push!(rcol,1.0f0)
    end

    prod = 1
    for k in rcol
      prod*=k
    end

    push!(flux, prod*rates[i])
    # make edges array
    for l in products[i]
        push!(edges,[match( r"([\d\s\.]*)(\D[\d\D]*)", j)[2],match( r"([\d\s\.]*)(\D[\d\D]*)", l)[2],i])
    end

  end
end

normalise(x) = x/norm(x)

t = 180

links = filter(i -> flux[i][t]>0 , 1:len(flux))
tflux = [log10(flux[i][t]) for i in links]
weight = 1+normalise(tflux)

newflux = [flux[i][t] for i in 1:len(flux)]
counter = 0
for i in links
  newflux[i] = weight[counter+=1]
end

edge = filter(i -> newflux[i[3]]>0 , edges)
source = [i[1] for i in edge]
target=[i[2] for i in edge]
weighted = [newflux[i[3]]+0.0001 for i in edge]

@rput source
@rput target
@rput weighted

R"library(igraph)"

R"el <- structure(list(V1 = source, V2 = target, weight = weighted), .Names = c('V1',
'V2', 'weight'), class = 'data.frame', row.names = c(1:$(length(edge))
))"

R"g <- graph.data.frame(el)"

R"net = g"

# get node names !!!!
R"nodes = V(g)$name"
@rget nodes

trasivity_pernode = R"transitivity(net, type='local')"
triad_census =R"triad_census(net)"

degree = R"degree(net, mode='all')"
indegree = R"degree(net, mode='in')"
outdegree = R"degree(net, mode='out')"

closeness = R"closeness(net, mode='all')# ,weights=NA)"

ec = R"eigen_centrality(net, directed=TRUE, weights=E(net)$weight)$vector"
betweenness = R"betweenness(net, directed=T)"

hubs = R"hub_score(net)$vector"
authorities = R"authority_score(net, weights=NA)$vector"

pagerank= R"page_rank(g, algo = c('prpack', 'arpack', 'power'), vids = V(g),
     directed = TRUE, damping = 0.85, personalized = NULL, weights = NULL,
     options = NULL)$vector"

####clustering below:

f = open("centrality.out","w")

write(f,join(nodes,","))
write(f,"\n")

for t in 1:length(flux[1])

  links = filter(i -> flux[i][t]>0 , 1:len(flux))
  tflux = [log10(flux[i][t]) for i in links]
  weight = 1+normalise(tflux)

  newflux = [flux[i][t] for i in 1:len(flux)]
  counter = 0
  for i in links
    newflux[i] = weight[counter+=1]
  end

  edge = filter(i -> newflux[i[3]]>0 , edges)
  source = [i[1] for i in edge]
  target=[i[2] for i in edge]
  weighted = [newflux[i[3]]+0.0001 for i in edge]

    @rput source
    @rput target
    @rput weighted

    R"el <- structure(list(V1 = source, V2 = target, weight = weighted), .Names = c('V1',
    'V2', 'weight'), class = 'data.frame', row.names = c(1:$(length(edge))
    ))"
    R"g <- graph.data.frame(el)"
    R"net = g"

    #centrality = R"eigen_centrality(net, directed=TRUE, weights=E(net)$weight)$vector"
    centrality =R"betweenness(net, directed=T)"

write(f,join(centrality,","))
write(f,"\n")
end

close(f)

max_stars_repo_name: wolfiex/RopaJL

max_issues_repo_name: wolfiex/RopaJL

max_forks_repo_name: wolfiex/RopaJL


text:
#!/usr/bin/env julia

# This file tests generation of multiblock VTK files.
# It also tests the non-default values of the optional parameters of vtk_grid,
# i.e., compress=false and append=false.

using WriteVTK
import Compat: UTF8String, view
typealias FloatType Float32
const vtm_filename_noext = "multiblock"

# Like sub2ind, but specifically for 3D arrays.
mysub2ind(dims, i, j, k) = i + dims[1]*(j-1) + dims[1]*dims[2]*(k-1)

function first_block_data()
    const Ni, Nj, Nk = 10, 15, 20

    x = zeros(FloatType, Ni, Nj, Nk)
    y = copy(x)
    z = copy(x)

    # Just to test subarrays:
    qall = zeros(FloatType, 3Ni, 2Nj, 2Nk)
    q = view(qall, 1:Ni, 1:Nj, 1:Nk)

    for k = 1:Nk, j = 1:Nj, i = 1:Ni
        r::FloatType = 1 + (i - 1)/(Ni - 1)
        th::FloatType = 3*pi/4 * (j-1) / (Nj-1)
        x[i, j, k] = r * cos(th)
        y[i, j, k] = r * sin(th)
        z[i, j, k] = (k-1) / (Nk-1)
        q[i, j, k] = k * exp(j) * sqrt(i)
    end

    return x, y, z, q
end

function second_block_data()
    const Ni, Nj, Nk = 20, 16, 12

    x = zeros(FloatType, Ni)
    y = zeros(FloatType, Nj)
    z = zeros(FloatType, Nk)

    [x[i] =  1 +   (i-1) / (Ni-1) for i = 1:Ni]
    [y[j] = -2 + 2*(j-1) / (Nj-1) for j = 1:Nj]
    [z[k] =        (k-1) / (Nk-1) for k = 1:Nk]

    q = zeros(FloatType, Ni, Nj, Nk)

    for k = 1:Nk, j = 1:Nj, i = 1:Ni
        q[i, j, k] = k * exp(j) * sqrt(i)
    end

    return x, y, z, q
end

function third_block_data()
    # This is basically the same as in the unstructured example.
    const Ni, Nj, Nk = 40, 50, 20
    const dims = [Ni, Nj, Nk]

    # Create points and point data.
    pts = Array(FloatType, 3, Ni, Nj, Nk)
    pdata = Array(FloatType, Ni, Nj, Nk)

    for k = 1:Nk, j = 1:Nj, i = 1:Ni
        r = 1 + (i - 1)/(Ni - 1)
        th = 3*pi/4 * (j - 1)/(Nj - 1)
        pts[1, i, j, k] = r * cos(th)
        pts[2, i, j, k] = r * sin(th)
        pts[3, i, j, k] = -2 + (k - 1)/(Nk - 1)
        pdata[i, j, k] = i*i + k*sqrt(j)
    end

    # Create cells (all hexahedrons in this case) and cell data.
    const celltype = VTKCellTypes.VTK_HEXAHEDRON
    cells = MeshCell[]
    cdata = FloatType[]

    for k = 2:Nk, j = 2:Nj, i = 2:Ni
        # Define connectivity of cell.
        inds = Array(Int32, 8)
        inds[1] = mysub2ind(dims, i-1, j-1, k-1)
        inds[2] = mysub2ind(dims, i  , j-1, k-1)
        inds[3] = mysub2ind(dims, i  , j  , k-1)
        inds[4] = mysub2ind(dims, i-1, j  , k-1)
        inds[5] = mysub2ind(dims, i-1, j-1, k  )
        inds[6] = mysub2ind(dims, i  , j-1, k  )
        inds[7] = mysub2ind(dims, i  , j  , k  )
        inds[8] = mysub2ind(dims, i-1, j  , k  )

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells, c)
        push!(cdata, i*j*k)
    end

    return pts, cells, pdata, cdata
end

function main()
    # Initialise multiblock file.
    vtm = vtk_multiblock(vtm_filename_noext)

    # Create first block.
    x, y, z, q = first_block_data()
    vtk = vtk_grid(vtm, x, y, z; compress=false, append=false)
    vtk_point_data(vtk, q, "q_values")

    # Create second block.
    x, y, z, q = second_block_data()
    vtk = vtk_grid(vtm, x, y, z; compress=false)
    vtk_point_data(vtk, q, "q_values")

    # Create third block.
    points, cells, q, c = third_block_data()
    vtk = vtk_grid(vtm, points, cells; append=false)
    vtk_point_data(vtk, q, "q_values")
    vtk_cell_data(vtk, c, "c_values")

    # Saved multiblock file and included block files.
    outfiles = vtk_save(vtm)
    println("Saved:", [" "^3 * s for s in outfiles]...)

    return outfiles::Vector{UTF8String}
end

main()


text:
using CSV, DataFrames
using DataStructures

function _read_csv(filename)
    println("Reading CSV...")
    m = CSV.read(filename, DataFrame; header=false)
    return Matrix{Int}(m)
end

struct AirportNode
    matrix
    depth::Int
    terminals::Int
    g::Int
    h::Int
    f::Int

    function AirportNode(matrix, depth, terminals; g=0, h=0)
        return new(matrix, depth, terminals, g, h, g+h)
    end
end

function _get_neighbors(matrix, (row, col); diagonal=true)
    if diagonal
        neighbors = []
        for nrow=row-1:row+1
            for ncol=col-1:col+1
                if nrow == row && ncol == col
                    continue
                else
                    push!(neighbors, matrix[nrow, ncol])
                end
            end
        end
    else
        top = matrix[row-1, col]
        right = matrix[row, col+1]
        bottom = matrix[row+1, col]
        left = matrix[row, col-1]
        neighbors = [top, right, bottom, left]
    end

    return neighbors
end

"""
Rule 1 states that at any given index, a 0 can become a 2 if:
    - Bordering a 1 on side
    - Other 3 sides are 0
"""
function check_rule_1(node, (row, col))
    state = node.matrix

    if state[row, col] != 0
        return false
    end

    neighbors = _get_neighbors(node.matrix, (row, col); diagonal=false)

    # we're looking for [0, 0, 0, 1]
    return sort(neighbors) == [0, 0, 0, 1]
end

"""
Rule 2 states that at any given index, a 0 can become a 1 if:
    - Neighbors existing “1”, even if diagonally
    - Doesn’t cause a violation in Rule 1 (handled in transition function)
"""
function check_rule_2(node, (row, col))
    if node.matrix[row, col] != 0
        return false
    end

    neighbors = _get_neighbors(node.matrix, (row, col))

    return 1 in neighbors
end

"""
Transition function for adding new walkways.
"""
function transition_walkway(node, (row, col))
    state = node.matrix

    if state[row, col] != 0
        return
    end

    if !check_rule_2(node, (row, col))
        return
    end

    new_m = copy(node.matrix)
    new_m[row, col] = 1
    new_node = AirportNode(
        new_m, node.depth+1, node.terminals; g=node.g+1, h=node.h)

    # check the new_m for rule 1 violations
    # new node cannot be directly next to a 2
    neighbors = _get_neighbors(new_node.matrix, (row, col); diagonal=false)

    if 2 in neighbors
        return
    end

    return new_node
end

"""
Transition function for adding new terminals.
"""
function transition_terminal(node, (row, col))
    if !check_rule_1(node, (row, col))
        return
    end

    new_m = copy(node.matrix)
    new_m[row, col] = 2
    new_node = AirportNode(
        new_m, node.depth+1, node.terminals+1; g=node.g, h=node.h-1)

    return new_node
end

"""
Depth First Search Airport Terminal solver.
"""
function airport_dfs(terminals, seed_file, max_depth)
    s_fns = [
        transition_terminal,
        transition_walkway
    ]
    nodes = Stack{AirportNode}()
    initial_state = _read_csv(seed_file)
    println("Running Airport Terminal Solver...")
    seed = AirportNode(initial_state, 1, 0)
    push!(nodes, seed)
    last_node = seed
    i = 0

    while !isempty(nodes)
        node = pop!(nodes)
        last_node = node
        state = node.matrix

        i += 1
        if i % 100000 == 0
            println("\n\nIterations: $(i). Stack size: $(length(nodes))")
        end

        if node.depth > max_depth
            continue
        end

        _is_goal(node) = node.terminals == terminals

        if _is_goal(node)
            return node
        end

        # Run walkway and terminal transition functions on current matrix
        # and add valid nodes to the stack.
        rows, cols = size(state)
        for row=1:rows
            for col=1:cols
                for fn in s_fns
                    new_node = fn(node, (row, col))

                    if isnothing(new_node)
                        continue
                    end

                    push!(nodes, new_node)
                end
            end
        end
    end
    println("failed. final state:\n$(last_node.matrix)")
end

"""
A* Search Airport Terminal solver.
"""
function airport_a_star(terminals, seed_file, max_depth)
    s_fns = [
        transition_terminal,
        transition_walkway
    ]
    nodes = PriorityQueue{AirportNode, Int}()
    initial_state = _read_csv(seed_file)
    println("Running Airport Terminal Solver...")
    seed = AirportNode(initial_state, 1, 0; h=terminals, g=0)
    enqueue!(nodes, seed, seed.f)
    last_node = seed
    i = 0

    while !isempty(nodes)
        node = dequeue!(nodes)
        last_node = node
        state = node.matrix

        i += 1
        if i % 100000 == 0
            println("\n\nIterations: $(i). Stack size: $(length(nodes)). g: $(node.g), h: $(node.h)")
        end

        if node.depth > max_depth
            continue
        end

        _is_goal(node) = node.terminals == terminals

        if _is_goal(node)
            res = """
                Solution found! Iterations: $(i). g: $(node.g). h: $(node.h). \
                Depth: $(node.depth)
            """
            println(res)
            display(node.matrix)
            return node
        end

        # Run walkway and terminal transition functions on current matrix
        # and add valid nodes to the stack.
        rows, cols = size(state)
        for row=1:rows
            for col=1:cols
                for fn in s_fns
                    new_node = fn(node, (row, col))

                    if isnothing(new_node)
                        continue
                    end

                    enqueue!(nodes, new_node, new_node.f)
                end
            end
        end
    end
    println("failed. final state:\n$(last_node.matrix)")
end

max_stars_repo_name: camirmas/DesignAutomation

max_issues_repo_name: camirmas/DesignAutomation

max_forks_repo_name: camirmas/DesignAutomation


text:
# Sample file for module  MarkovChains_monotonicity
# We check monotonicites of matrices P1 ... P20 from article
# Lorek, P., Markowski, P. "Monotonicity requirements for efficient exact sampling with Markov chains"

# Remark
# Julia 0.5.0 does not search for modules
# in current path by default. You may be forced to use first the following:
push!(LOAD_PATH, pwd())

using MarkovChains_monotonicity;

# Check weak_up, weak_down, mobius_up, mobius_down and stochastic monotonicity for given C and P
function check_all_monotonicities(C::Array{Int64,2},P::Union{Array{Rational{Int64},2},Array{Float64,2}})
   println("\t Weak up: \t ", is_weak_up(C,P));
   println("\t Weak down: \t ", is_weak_down(C,P));
   println("\t Mobius up: \t ", is_mobius_up(C,P));
   println("\t Mobius down: \t ", is_mobius_down(C,P));
   println("\t Stochastic: \t ", is_stochastic(C,P));
end

# general partial orderings
C1=[1 1 1 1 1 1; 0 1 0 1 1 1; 0 0 1 1 1 1; 0 0 0 1 0 1; 0 0 0 0 1 1; 0 0 0 0 0 1]
C2=[1 1 1 1 1 1 1 1; 0 1 0 0 1 1 0 1; 0 0 1 0 1 0 1 1; 0 0 0 1 0 1 1 1; 0 0 0 0 1 0 0 1; 0 0 0 0 0 1 0 1; 0 0 0 0 0 0 1 1; 0 0 0 0 0 0 0 1]
C3=[1 1 1 1 1; 0 1 0 0 1; 0 0 1 0 1; 0 0 0 1 1; 0 0 0 0 1]

P1=[3//8 1//8 1//8 1//8 0 1//8 1//8 0; 1//4 1//4 0 1//8 1//8 0 1//8 1//8; 1//4 0 3//8 0 0 1//4 0 1//8; 1//4 0 1//8 0 1//8 1//4 1//8 1//8; 0 1//8 1//8 0 0 1//4 1//4 1//4; 0 0 1//4 1//4 0 1//8 0 3//8; 0 0 1//8 1//8 1//4 1//4 1//8 1//8; 0 0 1//8 0 1//8 1//8 0 5//8]
P2=[5//8 0 1//8 1//8 0 1//8 0 0; 1//8 1//8 1//4 1//4 1//8 1//8 0 0; 3//8 0 1//8 0 1//4 1//4 0 0; 1//4 1//4 1//4 0 0 1//8 1//8 0; 1//8 1//8 1//4 1//8 0 1//8 0 1//4; 1//8 0 1//4 0 0 3//8 0 1//4; 1//8 1//8 0 1//8 1//8 0 1//4 1//4; 0 1//8 1//8 0 1//8 1//8 1//8 3//8]
P3=[1//2 1//6 0 1//3 0 0; 1//3 1//6 1//6 1//3 0 0; 1//3 1//6 0 1//3 1//6 0; 1//6 1//6 1//6 0 1//6 1//3; 1//6 0 1//6 1//6 1//3 1//6; 0 1//6 1//3 0 1//6 1//3]
P4=[2//5 1//5 1//5 1//5 0; 2//5 1//5 1//5 1//5 0; 0 2//5 2//5 1//5 0; 0 2//5 1//5 2//5 0; 0 1//5 2//5 0 2//5]
P5=[2//5 0 2//5 1//5 0; 0 2//5 1//5 2//5 0; 0 1//5 2//5 2//5 0; 0 1//5 1//5 1//5 2//5; 0 1//5 1//5 1//5 2//5]
P6=[17//24 0 0 1//8 1//8 1//24; 1//8 5//16 5//16 1//12 1//12 1//12; 1//8 5//16 5//16 1//12 1//12 1//12; 1//12 1//12 1//12 5//16 5//16 1//8; 1//12 1//12 1//12 5//16 5//16 1//8; 1//24 1//8 1//8 0 0 17//24]
P7=[1//3 1//3 1//3 0 0 0; 1//3 1//3 0 1//3 0 0; 1//3 0 1//3 1//3 0 0; 0 1//3 1//3 1//3 0 0; 0 1//6 1//6 1//6 1//6 1//3; 0 1//6 1//6 1//6 1//6 1//3]
P8=[2//5 1//5 1//5 1//5 0; 2//5 1//5 1//5 1//5 0; 2//5 0 1//5 1//5 1//5; 1//5 1//5 2//5 1//5 0; 0 2//5 1//5 0 2//5]
P9=[2//5 0 1//5 2//5 0; 0 1//5 2//5 1//5 1//5; 1//5 1//5 1//5 0 2//5; 0 1//5 1//5 1//5 2//5; 0 1//5 1//5 1//5 2//5]
P11=[1//3 1//6 1//6 1//6 1//6 0; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 0 1//6 1//6 1//6 1//6 1//3]
P12=[1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 0 1//3 1//6 1//6 1//6 1//6]
P13=[1//6 1//6 1//6 1//6 1//3 0; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6]
P14=[1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6; 1//6 1//6 1//6 1//6 1//6 1//6]
P15=[17//24 0 0 1//8 1//8 1//24; 1//8 5//16 5//16 1//12 1//12 1//12; 1//8 5//16 5//16 1//12 1//12 1//12; 1//12 1//12 1//12 5//16 5//16 1//8; 1//12 1//12 1//12 5//16 5//16 1//8; 1//24 1//16 1//16 1//16 1//16 17//24]
P16=[17//24 1//16 1//16 1//16 1//16 1//24; 1//8 5//16 5//16 1//12 1//12 1//12; 1//8 5//16 5//16 1//12 1//12 1//12; 1//12 1//12 1//12 5//16 5//16 1//8; 1//12 1//12 1//12 5//16 5//16 1//8; 1//24 1//8 1//8 0 0 17//24]

# tree-like ordering

C4=[1 0 0 0 1 0 1; 0 1 0 0 1 0 1; 0 0 1 0 0 1 1; 0 0 0 1 0 1 1; 0 0 0 0 1 0 1; 0 0 0 0 0 1 1 ; 0 0 0 0 0 0 1]

P17=[1//2 1//2 0 0 0 0 0; 0 0 1//2 1//2 0 0 0; 1//2 1//2 0 0 0 0 0; 0 0 1//2 1//2 0 0 0;   0 0 0 0 1//2 1//2 0;  0 0 0 0 1//2 1//2 0;  0 0 0 0 1//2 1//2 0]
P18=[1//7 1//7 1//7 1//7 1//7 1//7 1//7 ;1//7 1//7 1//7 1//7 1//7 1//7 1//7 ;1//7 1//7 1//7 1//7 1//7 1//7 1//7 ;1//7 1//7 1//7 1//7 1//7 1//7 1//7 ;1//7 1//7 1//7 1//7 1//7 1//7 1//7 ;1//7 1//7 1//7 1//7 1//7 1//7 1//7 ;1//7 1//7 1//7 1//7 1//7 1//7 1//7 ]
P19=eye(7)
P20=[1//2 1//2 0 0 0 0 0; 0 0 1//2 1//2 0 0 0; 1//2 1//2 0 0 0 0 0; 0 0 1//2 1//2 0 0 0;   0 0 0 0 1//2 1//2 0;  0 0 0 0 1//2 1//2 0;  0 0 0 0 0 0 1]

println("\nChecking Example 1 (C2,P1) ..");
check_all_monotonicities(C2,P1)

println("\nChecking Example 2 (C2,P2) ..");
check_all_monotonicities(C2,P2)

println("\nChecking Example 3 (C1,P3) ..");
check_all_monotonicities(C1,P3)

println("\nChecking Example 4 (C3,P4) ..");
check_all_monotonicities(C3,P4)

println("\nChecking Example 5 (C3,P5) ..");
check_all_monotonicities(C3,P5)

println("\nChecking Example 6 (C1,P6) ..");
check_all_monotonicities(C1,P6)

println("\nChecking Example 7 (C1,P7) ..");
check_all_monotonicities(C1,P7)

println("\nChecking Example 8 (C3,P8) ..");
check_all_monotonicities(C3,P8)

println("\nChecking Example 9 (C3,P9) ..");
check_all_monotonicities(C3,P9)

println("\nChecking Example 11 (C1,P11) ..");
check_all_monotonicities(C1,P11)

println("\nChecking Example 12 (C1,P12) ..");
check_all_monotonicities(C1,P12)

println("\nChecking Example 13 (C1,P13) ..");
check_all_monotonicities(C1,P13)

println("\nChecking Example 14 (C1,P14) ..");
check_all_monotonicities(C1,P14)

println("\nChecking Example 15 (C1,P15) ..");
check_all_monotonicities(C1,P15)

println("\nChecking Example 16 (C1,P16) ..");
check_all_monotonicities(C1,P16)

println("\nChecking Example 17 (C4,P17) ..");
check_all_monotonicities(C4,P17)

println("\nChecking Example 18 (C4,P18) ..");
check_all_monotonicities(C4,P18)

println("\nChecking Example 19 (C4,P19) ..");
check_all_monotonicities(C4,P19)

println("\nChecking Example 20 (C4,P20) ..");
check_all_monotonicities(C4,P20)


text:
function fit(; y::Vector{M}, niters::Int, batchsize::Int, c::Constants,
             opt=ADAM(1e-2), print_freq::Int=10, nsave::Int=0,
             scale_logqy::Float64=1.0,
             init::Union{StateMP, Nothing}=nothing,
             seed=nothing, flushOutput::Bool=false,
             verbose::Int=1) where M <: Matrix

  if verbose >= 1
    println("opt: $opt")
    println("niters: $niters")
    println("batchsize: $batchsize")
    printConstants(c)
  end

  # Set a random seed for reproducibility if seed is provided
  if seed != nothing
    Random.seed!(seed)
  end

  # Set save frequency
  if 0 < nsave < niters
    save_freq = ceil(Int, niters / nsave)
  elseif nsave >= niters
    save_freq = 1
  else # nsave <= 1
    save_freq = niters
  end

  # Initialize if needed
  if init == nothing
    init = State(c)
  end

  # Create metrics dictionary
  metrics = Dict(m => Float64[] for m in (:ll, :lp, :lq, :elbo))

  # Make a deep copy of initial state
  state = deepcopy(init)

  # Collect all model parameters and variational parameters
  ps = ADVI.vparams(state)

  # Print message
  function printMsg(iter::Integer)
    if iter % print_freq == 0
      m = ["$(key): $(round(metrics[key][end] / sum(c.N), digits=3))"
           for key in keys(metrics)]
      println("$(ShowTime()) | $(iter)/$(niters) | $(join(m, " | "))")
    end
  end

  # Create loss function
  function loss(y_batch::T) where T
    out = -compute_elbo(state, y_batch, c, metrics, scale_logqy=scale_logqy) / sum(c.N)
    if isnan(out) || isinf(out)
      out = zero(out)
      printMsg(0)
      println("Skipping update -- Loss was NaN or Inf")
    end
    return out
  end

  # Create an object to store history of state
  state_hist = [deepcopy(init)]

  # Main loop
  println("Compiling model...")
  @time for iter in 1:niters
    # Sample minibatch
    y_mini = [if 0 < batchsize < c.N[i]
                idx = Distributions.sample(1:c.N[i], batchsize, replace=false)
                y[i][idx, :]
              else
                y[i]
              end for i in 1:c.I]

    # Compute and update gradients. Two different ways.
    # Flux.train!(loss, ps, [(y_mini, )], opt)
    gs = Tracker.gradient(() -> loss(y_mini), ps)
    Flux.Tracker.update!(opt, ps, gs)

    printMsg(iter)

    if iter % save_freq == 0
      append!(state_hist, [deepcopy(state)])
    end
    # println("eps: $(rsample(state)[2].eps)")

    if flushOutput
      flush(stdout)
    end
  end

  # Put things in output
  out = Dict(:metrics => metrics,
             :state_hist => state_hist,
             :state => state,
             :batchsize => batchsize,
             :niters => niters,
             :seed => seed,
             :c => c)

  return out
end


text:
function networkTiming(; iters = 100)
    for model in [vKITTI_Conv_D,vKITTI_Conv_RGB]
        x = KnetArray(zeros(Float32,224,224,3,1));
        w = Array[
            #Conv64
            xavier(Float32,3,3,3,64),        zeros(Float32,1,1,64,1),
            xavier(Float32,3,3,64,64),       zeros(Float32,1,1,64,1),
            #Conv128
            xavier(Float32,3,3,64,128),      zeros(Float32,1,1,128,1),
            xavier(Float32,3,3,128,128),     zeros(Float32,1,1,128,1),
            #Conv256
            xavier(Float32,3,3,128,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            #Conv512
            xavier(Float32,3,3,256,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            #Scores
            xavier(Float32,1,1,512,13),      zeros(Float32,1,1,13,1),
            #Deconv
            bilinear(Float32,32,32,13,13),   zeros(Float32,1,1,13,1),
        ]; w = map(KnetArray,w);

        println(model)
        times = []
        for i=1:iters
            push!(times,(@timed y = model(w,x))[2])
        end
        println("Minimum: $(minimum(times))s ($(1/minimum(times)) FPS) Maximum: $(maximum(times))s ($(1/maximum(times)) FPS) Average: $(mean(times))s ($(1/mean(times)) FPS)")
    end
    for model in [vKITTI_Conv_RGBD,]
        x = KnetArray(zeros(Float32,224,224,6,1));
        w = Array[
            #RGB
            #Conv64
            xavier(Float32,3,3,3,64),        zeros(Float32,1,1,64,1),
            xavier(Float32,3,3,64,64),       zeros(Float32,1,1,64,1),
            #Conv128
            xavier(Float32,3,3,64,128),      zeros(Float32,1,1,128,1),
            xavier(Float32,3,3,128,128),     zeros(Float32,1,1,128,1),
            #Conv256
            xavier(Float32,3,3,128,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            #Conv512
            xavier(Float32,3,3,256,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            #Colorized Depth
            #Conv64
            xavier(Float32,3,3,3,64),        zeros(Float32,1,1,64,1),
            xavier(Float32,3,3,64,64),       zeros(Float32,1,1,64,1),
            #Conv128
            xavier(Float32,3,3,64,128),      zeros(Float32,1,1,128,1),
            xavier(Float32,3,3,128,128),     zeros(Float32,1,1,128,1),
            #Conv256
            xavier(Float32,3,3,128,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            #Conv512
            xavier(Float32,3,3,256,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            #Scores
            xavier(Float32,1,1,512,13),      zeros(Float32,1,1,13,1),
            #Deconv
            bilinear(Float32,32,32,13,13),   zeros(Float32,1,1,13,1),
        ]; w = map(KnetArray,w);

        println(model)
        times = []
        for i=1:iters
            push!(times,(@timed y = model(w,x))[2])
        end
        println("Minimum: $(minimum(times))s ($(1/minimum(times)) FPS) Maximum: $(maximum(times))s ($(1/maximum(times)) FPS) Average: $(mean(times))s ($(1/mean(times)) FPS)")
    end
    println("***********************")
    h = KnetArray(zeros(Float32,28,28,512,1))
    c = KnetArray(zeros(Float32,28,28,512,1))
    for model in [vKITTI_LSTM_D,vKITTI_LSTM_RGB]
        x = KnetArray(zeros(Float32,224,224,3,1));
        w = Array[
            #Conv64
            xavier(Float32,3,3,3,64),        zeros(Float32,1,1,64,1),
            xavier(Float32,3,3,64,64),       zeros(Float32,1,1,64,1),
            #Conv128
            xavier(Float32,3,3,64,128),      zeros(Float32,1,1,128,1),
            xavier(Float32,3,3,128,128),     zeros(Float32,1,1,128,1),
            #Conv256
            xavier(Float32,3,3,128,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            #Conv512
            xavier(Float32,3,3,256,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            #LSTM
            xavier(Float32,3,3,512,512),     xavier(Float32,3,3,512,512),   zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     xavier(Float32,3,3,512,512),   zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     xavier(Float32,3,3,512,512),   zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     xavier(Float32,3,3,512,512),   zeros(Float32,1,1,512,1),
            #Scores
            xavier(Float32,1,1,512,13),      zeros(Float32,1,1,13,1),
            #Deconv
            bilinear(Float32,32,32,13,13),   zeros(Float32,1,1,13,1),
        ]; w = map(KnetArray,w);

        println(model)
        times = []
        for i=1:iters
            push!(times,(@timed y = model(w,h,c,x))[2])
        end
        println("Minimum: $(minimum(times))s ($(1/minimum(times)) FPS) Maximum: $(maximum(times))s ($(1/maximum(times)) FPS) Average: $(mean(times))s ($(1/mean(times)) FPS)")
    end
    for model in [vKITTI_LSTM_RGBD,]
        x = KnetArray(zeros(Float32,224,224,6,1));
        w = Array[
            #RGB
            #Conv64
            xavier(Float32,3,3,3,64),        zeros(Float32,1,1,64,1),
            xavier(Float32,3,3,64,64),       zeros(Float32,1,1,64,1),
            #Conv128
            xavier(Float32,3,3,64,128),      zeros(Float32,1,1,128,1),
            xavier(Float32,3,3,128,128),     zeros(Float32,1,1,128,1),
            #Conv256
            xavier(Float32,3,3,128,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            #Conv512
            xavier(Float32,3,3,256,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            #Colorized Depth
            #Conv64
            xavier(Float32,3,3,3,64),        zeros(Float32,1,1,64,1),
            xavier(Float32,3,3,64,64),       zeros(Float32,1,1,64,1),
            #Conv128
            xavier(Float32,3,3,64,128),      zeros(Float32,1,1,128,1),
            xavier(Float32,3,3,128,128),     zeros(Float32,1,1,128,1),
            #Conv256
            xavier(Float32,3,3,128,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
            #Conv512
            xavier(Float32,3,3,256,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
            #LSTM
            xavier(Float32,3,3,512,512),     xavier(Float32,3,3,512,512),   zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     xavier(Float32,3,3,512,512),   zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     xavier(Float32,3,3,512,512),   zeros(Float32,1,1,512,1),
            xavier(Float32,3,3,512,512),     xavier(Float32,3,3,512,512),   zeros(Float32,1,1,512,1),
            #Scores
            xavier(Float32,1,1,512,13),      zeros(Float32,1,1,13,1),
            #Deconv
            bilinear(Float32,32,32,13,13),   zeros(Float32,1,1,13,1),
        ]; w = map(KnetArray,w);

        println(model)
        times = []
        for i=1:iters
            push!(times,(@timed y = model(w,h,c,x))[2])
        end
        println("Minimum: $(minimum(times))s ($(1/minimum(times)) FPS) Maximum: $(maximum(times))s ($(1/maximum(times)) FPS) Average: $(mean(times))s ($(1/mean(times)) FPS)")
    end
end

max_stars_repo_name: eeyrdkl/MScThesis

max_issues_repo_name: eeyrdkl/MScThesis

max_forks_repo_name: eeyrdkl/MScThesis


text:
using VPTrees, DataStructures
using Base.Test

function euclidean(coords::Matrix, i::Int, j::Int)
    @inbounds x = coords[1, i] - coords[1, j]
    @inbounds y = coords[2, i] - coords[2, j]
    @inbounds z = coords[3, i] - coords[3, j]

    return sqrt(x^2 + y^2 + z^2)
end

function euclidean_periodic(coords::Matrix, i::Int, j::Int)
    xbw, ybw, zbw = 1.0, 1.0, 1.0

    xhbw = xbw / 2
    yhbw = ybw / 2
    zhbw = zbw / 2

    @inbounds x = coords[1, i] - coords[1, j]
    @inbounds y = coords[2, i] - coords[2, j]
    @inbounds z = coords[3, i] - coords[3, j]

    if x < -xhbw x += xbw elseif x > xhbw x -= xbw end
    if y < -yhbw y += ybw elseif y > yhbw y -= ybw end
    if z < -zhbw z += zbw elseif z > zhbw z -= zbw end

    return sqrt(x^2 + y^2 + z^2)
end

function is_consistent{T}(metric::Function, node::Node{T}, point_index::Int, τ::T, isless)
    node.index == 0 && return true

    μ = metric(node.index, point_index)
    #println("$(point_index) and $(node.index): $μ <= $τ is $(μ <= τ) and isless = $isless")
    ((μ <= τ) ⊻ isless) && return false

    node.isleaf && return true

    is_consistent(metric, node.inside, point_index, τ, isless) || return false
    is_consistent(metric, node.outside, point_index, τ, isless) || return false

    return true
end

function is_consistent{T}(metric::Function, node::Node{T})
    node.isleaf && return true

    is_consistent(metric, node.inside, node.index, node.radius, true) || return false
    is_consistent(metric, node.outside, node.index, node.radius, false) || return false
    #println()
    is_consistent(metric, node.inside) || return false
    is_consistent(metric, node.outside) || return false

    return true
end

is_consistent(tree::VPTree) = is_consistent(tree.metric, tree.root)

create_metric(distance::Function, points) = (i::Int, j::Int) -> distance(points, i, j)

function test_knn(distance::Function, num_points::Int, num_iters::Int)
    for s in 1:num_iters
        srand(s)
        points = rand(3, num_points)
        point_indices = collect(1:num_points)
        metric = create_metric(distance, points)

        num_neighbors = 5

        tree = VPTree(metric, point_indices)
        neighbors = NeighborList{Float64}(num_neighbors)

        for i in 1:num_points
            knn!(neighbors, tree, i)
            indices_test = indices(neighbors)

            pq = PriorityQueue(Int, Float64, Base.Order.Reverse)
            for j in 1:num_points
                i == j && continue

                r = metric(i, j)

                if length(pq) < num_neighbors
                    enqueue!(pq, j, r)
                elseif r < last(peek(pq))
                    dequeue!(pq)
                    enqueue!(pq, j, r)
                end
            end

            indices_actual = keys(pq)

            isempty(symdiff(indices_test, indices_actual)) || return false
        end
    end

    return true
end

function test_rangesearch(distance::Function, num_points::Int, num_iters::Int)
    for s in 1:num_iters
        srand(s)
        points = rand(3, num_points)
        point_indices = collect(1:num_points)
        metric = create_metric(distance, points)

        tree = VPTree(metric, point_indices)
        neighbors = NeighborList{Float64}(num_points)

        r = rand()
        for i in 1:num_points
            rangesearch!(neighbors, tree, i, r)
            indices_test = indices(neighbors)

            indices_actual = filter(j -> metric(i, j) <= r, 1:num_points)

            isempty(symdiff(indices_test, indices_actual)) || return false
        end
    end

    return true
end

function test_integrity(distance::Function, num_points::Int, num_iters::Int)
    for s in 1:num_iters
        srand(s)
        points = rand(3, num_points)
        point_indices = collect(1:num_points)
        metric = create_metric(distance, points)

        tree = VPTree(metric, point_indices)
        is_consistent(tree) || return false
    end

    return true
end

@testset "VP Tree tests" begin
    const num_points = 1000
    const num_iters = 5

    @testset "VP Tree integrity test" begin
        @test test_integrity(euclidean, num_points, num_iters)
        @test test_integrity(euclidean_periodic, num_points, num_iters)
    end

    @testset "k-nearest neighbors test" begin
        @test test_knn(euclidean, num_points, num_iters)
        @test test_knn(euclidean_periodic, num_points, num_iters)
    end

    @testset "Range search test" begin
        @test test_rangesearch(euclidean, num_points, num_iters)
        @test test_rangesearch(euclidean_periodic, num_points, num_iters)
    end
end


text:
"""
Get the latest VersionNumber for base_path/Versions.toml
"""
function _get_latest_version(base_path::AbstractString)
    versions_file_path = joinpath(base_path, "Versions.toml")

    if isfile(versions_file_path)
        versions_content = parsefile(versions_file_path)
        versions = [VersionNumber(v) for v in collect(keys(versions_content))]

        return first(findmax(versions))
    end
end

"""
    _find_latest_pkg_entry(pkg_name::Union{AbstractString, Missing}, pkg_uuid::Union{Missing, UUID}=missing; registries::Array{RegistryInstance}=reachable_registries())

Given a list of registries, find the entry corresponding to the given package (specified as a `name` and/or an `UUID`) with the latest version.

For example, if a package is registered in two registries with the same name and UUID, returns the `PkgEntry` corresponding to the registry which
has the higher version number of the package.
"""
function _find_latest_pkg_entry(pkg_name::Union{AbstractString, Missing}, pkg_uuid::Union{Missing, UUID}=missing; registries::Array{RegistryInstance}=reachable_registries())
    if ismissing(pkg_name) && ismissing(pkg_uuid)
        throw(ArgumentError("Must supply either a `pkg_name` or `pkg_uuid`!"))
    end

    entries = PkgEntry[]

    for rego in registries
        for (name, entry) in rego.pkgs
            if !ismissing(pkg_name)
                pkg_name == name || continue
            end
            if !ismissing(pkg_uuid)
                pkg_uuid == entry.uuid || continue
            end
            push!(entries, entry)
        end
    end

    if isempty(entries)
        throw(PackageNotInRegistry("No package found with supplied name and/or UUID. Got: `pkg_name`=$(pkg_name), `pkg_uuid`=$(pkg_uuid)."))
    elseif length(entries) == 1
        return only(entries)
    else
        base_path = pkg_entry -> joinpath(pkg_entry.registry_path, pkg_entry.path)

        # uses `Compat` for the 2-argument `argmax` introduced in Julia v1.7
        return argmax(_get_latest_version ∘ base_path, entries)
    end
end

"""
Use levenshtein distance to find packages closely named to pkg_to_compare

Rules:
- Ignore package name casings
- Only return matches that are between 1 and 2 distances
"""
function _find_alternative_packages(pkg_to_compare::String, packages::Array)
    MIN_LEVENSHTEIN = 1
    MAX_LEVENSHTEIN = 2

    pkg_to_compare = uppercase(pkg_to_compare)
    pkg_distances = [(REPL.levenshtein(uppercase(pkg), pkg_to_compare), pkg) for pkg in packages]

    return [x[2] for x in pkg_distances if MIN_LEVENSHTEIN <= x[1] <= MAX_LEVENSHTEIN]
end

"""
Get the package name from a UUID
"""
function _get_pkg_name(uuid::UUID, registries)
    for rego in registries
        for (pkg_name, pkg_entry) in rego.pkgs
            if pkg_entry.uuid == uuid
                return pkg_name
            end
        end
    end

    throw(NoUUIDMatch("No package found with the UUID $uuid"))
end

function _get_pkg_name(uuid::UUID; kwargs...)
    registries = reachable_registries(; kwargs...)
    return _get_pkg_name(uuid, registries)
end

"""
Get the UUID from a package name and the registry it is in.
Specify a registry name as well to avoid ambiguity with same package names in multiple registries.
"""
function _get_pkg_uuid(
    pkg_name::String, registry_name::String;
    depots::Union{String, Vector{String}}=Base.DEPOT_PATH,
)
    registry = only(reachable_registries(registry_name; depots=depots))
    return _get_pkg_uuid(pkg_name, registry)
end

function _get_pkg_uuid(pkg_name::String, registry::RegistryInstance)
    if haskey(registry.pkgs, pkg_name)
        return registry.pkgs[pkg_name].uuid
    else
        alt_packages = _find_alternative_packages(pkg_name, collect(keys(registry.pkgs)))

        warning = "The package $(pkg_name) was not found in the $(registry.name) registry."

        if !isempty(alt_packages)
            warning *= "\nPerhaps you meant: $(string(alt_packages))"
        end

        warning *= "\nOr you can search in another registry using `users(\"$(pkg_name)\", \"OtherRegistry\")`"
        throw(PackageNotInRegistry(warning))
    end
end


text:
export OneLine2Trajectory, oneLine2dNodes, oneLine2dDensity

"""
    OneLine2dTrajectory(numProfiles, numSamplingPerProfile
                  ; TE::Float64=0.0
                  , AQ::Float64=1.e-3
                  , angle::Float64=0.0
                  , kargs...)

returns a trajectory consisting of one arbitrarily rotated profile.

# Arguments
* `numProfiles::Int64`            - number of profiles
* `numSamplingPerProfile::Int64`  - number of sampling points per profile
* (`TE::Float64=0.0`)             - echo time in s
* (`AQ::Float64=1.e-3`)           - readout duration in s (per profile)
* (`angle::Float64=0.0`)          - angle of the profile (with respect to the x-axis) in radians
"""
function OneLine2dTrajectory(::Type{T}, numProfiles, numSamplingPerProfile
                  ; TE=0.0
                  , AQ=1.e-3
                  , angle=0.0
                  , kargs...) where T
  nodes = oneLine2dNodes(T, numProfiles, numSamplingPerProfile; angle=angle)
  times = readoutTimes(T, numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ)
  return  Trajectory("OneLine", nodes, times, TE, AQ, numProfiles, numSamplingPerProfile, 1, false, true)
end

function oneLine2dNodes(::Type{T}, numProfiles, numSamplingPerProfile; angle::Float64=0.0, kargs...) where T
  nodes = zeros(T, 2,numSamplingPerProfile, numProfiles)
  angles = zeros(T, numProfiles)
  angles[:] .= angle
  pos = collect((0:numSamplingPerProfile-1)/numSamplingPerProfile .- 0.5)
  for l = 1:numProfiles
    for k = 1:numSamplingPerProfile
      nodes[1,k,l] = (-1)^l * pos[k]*cos(angles[l])
      nodes[2,k,l] = (-1)^l * pos[k]*sin(angles[l])
    end
  end
  return reshape(nodes, 2, numSamplingPerProfile*numProfiles)
end

function oneLine2dDensity(::Type{T}, numSamplingPerProfile::Int64, numProfiles::Int64) where T
  density = zeros(T, tr.numSamplingPerProfile, numProfiles)
  pos = collect((0:numSamplingPerProfile-1)/numSamplingPerProfile .- 0.5)
  for l = 1:numProfiles
    for k = 1:numSamplingPerProfile
      density[k,l] = abs(pos[k])
    end
  end
  return vec(density)
end


text:
# Linear covariance function

@inline dotijp(X1::AbstractMatrix, X2::AbstractMatrix, i::Int, j::Int, p::Int) = X1[p,i]*X2[p,j]
@inline function dotij(X1::AbstractMatrix, X2::AbstractMatrix, i::Int, j::Int, dim::Int)
	s=zero(promote_type(eltype(X1), eltype(X2)))
	@inbounds @simd for p in 1:dim
		s+=dotijp(X1,X2,i,j,p)
	end
	return s
end
include("lin_iso.jl")
include("lin_ard.jl")

"""
    Lin(ll::Union{Real,Vector{Real}})

Create linear kernel with length scale `exp.(ll)`.

See also [`LinIso`](@ref) and [`LinArd`](@ref).
"""
Lin(ll::Real) = LinIso(ll)
Lin(ll::AbstractVector{<:Real}) = LinArd(ll)