text:
@doc raw"""
    AsymptoticSKCETest(kernel::Kernel, predictions, targets)

Calibration hypothesis test based on the unbiased estimator of the squared kernel
calibration error (SKCE) with quadratic sample complexity.

# Details

Let ``\mathcal{D} = (P_{X_i}, Y_i)_{i=1,\ldots,n}`` be a data set of predictions and
corresponding targets. Denote the null hypothesis "the predictive probabilistic model is
calibrated" with ``H_0``.

The hypothesis test approximates the p-value ``ℙ(\mathrm{SKCE}_{uq} > c \,|\, H_0)``, where
``\mathrm{SKCE}_{uq}`` is the unbiased estimator of the SKCE, defined as
```math
\frac{2}{n(n-1)} \sum_{1 \leq i < j \leq n} h_k\big((P_{X_i}, Y_i), (P_{X_j}, Y_j)\big),
```
where
```math
\begin{aligned}
h_k\big((μ, y), (μ', y')\big) ={}&   k\big((μ, y), (μ', y')\big)
                                   - 𝔼_{Z ∼ μ} k\big((μ, Z), (μ', y')\big) \\
                                 & - 𝔼_{Z' ∼ μ'} k\big((μ, y), (μ', Z')\big)
                                   + 𝔼_{Z ∼ μ, Z' ∼ μ'} k\big((μ, Z), (μ', Z')\big).
\end{aligned}
```

The p-value is estimated based on the asymptotically valid approximation
```math
ℙ(n\mathrm{SKCE}_{uq} > c \,|\, H_0) \approx ℙ(T > c \,|\, \mathcal{D}),
```
where ``T`` is the bootstrap statistic
```math
T = \frac{2}{n} \sum_{1 \leq i < j \leq n} \bigg(h_k\big((P^*_{X_i}, Y^*_i), (P^*_{X_j}, Y^*_j)\big)
- \frac{1}{n} \sum_{r = 1}^n h_k\big((P^*_{X_i}, Y^*_i), (P_{X_r}, Y_r)\big)
- \frac{1}{n} \sum_{r = 1}^n h_k\big((P_{X_r}, Y_r), (P^*_{X_j}, Y^*_j)\big)
+ \frac{1}{n^2} \sum_{r, s = 1}^n h_k\big((P_{X_r}, Y_r), (P_{X_s}, Y_s)\big)\bigg)
```
for bootstrap samples ``(P^*_{X_i}, Y^*_i)_{i=1,\ldots,n}`` of ``\mathcal{D}``.
This can be reformulated to the approximation
```math
ℙ(n\mathrm{SKCE}_{uq}/(n - 1) - \mathrm{SKCE}_b > c \,|\, H_0) \approx ℙ(T' > c \,|\, \mathcal{D}),
```
where
```math
\mathrm{SKCE}_b = \frac{1}{n^2} \sum_{i, j = 1}^n h_k\big((P_{X_i}, Y_i), (P_{X_j}, Y_j)\big)
```
and
```math
T' = \frac{2}{n(n - 1)} \sum_{1 \leq i < j \leq n} h_k\big((P^*_{X_i}, Y^*_i), (P^*_{X_j}, Y^*_j)\big)
- \frac{2}{n^2} \sum_{i, r=1}^n h_k\big((P^*_{X_i}, Y^*_i), (P_{X_r}, Y_r)\big).
```

# References

Widmann, D., Lindsten, F., & Zachariah, D. (2019). [Calibration tests in multi-class
classification: A unifying framework](https://proceedings.neurips.cc/paper/2019/hash/1c336b8080f82bcc2cd2499b4c57261d-Abstract.html).
In: Advances in Neural Information Processing Systems (NeurIPS 2019) (pp. 12257–12267).

Widmann, D., Lindsten, F., & Zachariah, D. (2021). [Calibration tests beyond
classification](https://openreview.net/forum?id=-bxf89v3Nx).
"""
struct AsymptoticSKCETest{K<:Kernel,E,V,M} <: HypothesisTests.HypothesisTest
    """Kernel."""
    kernel::K
    """Calibration error estimate."""
    estimate::E
    """Test statistic."""
    statistic::V
    """Symmetric kernel matrix, consisting of pairwise evaluations of ``h_{ij}``."""
    kernelmatrix::M
end

function AsymptoticSKCETest(
    kernel::Kernel, predictions::AbstractVector, targets::AbstractVector
)
    # compute the calibration error estimate, the test statistic, and the kernel matrix
    estimate, statistic, kernelmatrix = estimate_statistic_kernelmatrix(
        kernel, predictions, targets
    )

    return AsymptoticSKCETest(kernel, estimate, statistic, kernelmatrix)
end

# HypothesisTests interface

HypothesisTests.default_tail(::AsymptoticSKCETest) = :right

function HypothesisTests.pvalue(test::AsymptoticSKCETest; kwargs...)
    return pvalue(Random.GLOBAL_RNG, test; kwargs...)
end
function HypothesisTests.pvalue(
    rng::Random.AbstractRNG, test::AsymptoticSKCETest; bootstrap_iters::Int=1_000
)
    return bootstrap_ccdf(rng, test.statistic, test.kernelmatrix, bootstrap_iters)
end

HypothesisTests.testname(::AsymptoticSKCETest) = "Asymptotic SKCE test"

# parameter of interest: name, value under H0, point estimate
function HypothesisTests.population_param_of_interest(test::AsymptoticSKCETest)
    return "SKCE", zero(test.estimate), test.estimate
end

function HypothesisTests.show_params(io::IO, test::AsymptoticSKCETest, ident="")
    return println(io, ident, "test statistic: $(test.statistic)")
end

@doc raw"""
    estimate_statistic_kernelmatrix(kernel, predictions, targets)

Compute the estimate of the SKCE, the test statistic, and the matrix of the evaluations of
the kernel function.

# Details

Let ``\mathcal{D} = (P_{X_i}, Y_i)_{i=1,\ldots,n}`` be a data set of predictions and
corresponding targets.

The unbiased estimator ``\mathrm{SKCE}_{uq}`` of the SKCE is defined as
```math
\frac{2}{n(n-1)} \sum_{1 \leq i < j \leq n} h_k\big((P_{X_i}, Y_i), (P_{X_j}, Y_j)\big),
```
where
```math
\begin{aligned}
h_k\big((μ, y), (μ', y')\big) ={}&   k\big((μ, y), (μ', y')\big)
                                   - 𝔼_{Z ∼ μ} k\big((μ, Z), (μ', y')\big) \\
                                 & - 𝔼_{Z' ∼ μ'} k\big((μ, y), (μ', Z')\big)
                                   + 𝔼_{Z ∼ μ, Z' ∼ μ'} k\big((μ, Z), (μ', Z')\big).
\end{aligned}
```

The test statistic is defined as
```math
\frac{n}{n-1} \mathrm{SKCE}_{uq} - \mathrm{SKCE}_b,
```
where
```math
\mathrm{SKCE}_b = \frac{1}{n^2} \sum_{i, j = 1}^n h_k\big((P_{X_i}, Y_i), (P_{X_j}, Y_j)\big)
```
(see [`AsymptoticSKCETest`](@ref)). This is equivalent to
```math
\frac{1}{n^2} \sum_{i, j = 1} h_k\big((P_{X_i}, Y_i), (P_{X_j}, Y_j)\big) \bigg(\frac{n^2}{(n - 1)^2} 1(i \neq j) - \bigg).
```

The kernelmatrix ``K \in \mathbb{R}^{n \times n}`` is defined as
```math
K_{ij} = h_k\big((P_{X_i}, Y_i), (P_{X_j}, Y_j)\big)
```
for ``i, j \in \{1, \ldots, n\}``.
"""
function estimate_statistic_kernelmatrix(kernel, predictions, targets)
    # obtain number of samples
    nsamples = length(predictions)
    nsamples > 1 || error("there must be at least two samples")

    # pre-computations
    α = (2 * nsamples - 1) / (nsamples - 1)^2

    @inbounds begin
        # evaluate the kernel function for the first pair of samples
        prediction = predictions[1]
        target = targets[1]

        # initialize the kernel matrix
        hij = unsafe_skce_eval(kernel, prediction, target, prediction, target)
        kernelmatrix = Matrix{typeof(hij)}(undef, nsamples, nsamples)
        kernelmatrix[1, 1] = hij

        # initialize the test statistic and the unbiased estimate of the SKCE
        statistic = -hij / 1
        estimate = zero(statistic)

        # add evaluations of all other pairs of samples
        nstatistic = 1
        nestimate = 0
        for i in 2:nsamples
            predictioni = predictions[i]
            targeti = targets[i]

            for j in 1:(i - 1)
                predictionj = predictions[j]
                targetj = targets[j]

                # evaluate the kernel function
                hij = unsafe_skce_eval(kernel, predictioni, targeti, predictionj, targetj)

                # update the kernel matrix
                kernelmatrix[j, i] = hij

                # update the estimate and the test statistic
                nstatistic += 2
                statistic += 2 * (α * hij - statistic) / nstatistic
                nestimate += 1
                estimate += (hij - estimate) / nestimate
            end

            # evaluate the kernel function for the `i`th sample
            hij = unsafe_skce_eval(kernel, predictioni, targeti, predictioni, targeti)

            # update the kernel matrix
            kernelmatrix[i, i] = hij

            # update the test statistic
            nstatistic += 1
            statistic -= (statistic + hij) / nstatistic
        end
    end

    # add lower triangle of the kernel matrix
    LinearAlgebra.copytri!(kernelmatrix, 'U')

    return estimate, statistic, kernelmatrix
end

@doc raw"""
    bootstrap_ccdf(rng::AbstractRNG, statistic, kernelmatrix, bootstrap_iters::Int)

Estimate the value of the inverse CDF of the test statistic under the calibration null
hypothesis by bootstrapping.

# Details

Let ``\mathcal{D} = (P_{X_i}, Y_i)_{i=1,\ldots,n}`` be a data set of predictions and
corresponding targets. Denote the null hypothesis "the predictive probabilistic model is
calibrated" with ``H_0``, and the test statistic with ``T``.

The value of the inverse CDF under the null hypothesis is estimated based on the
asymptotically valid approximation
```math
ℙ(T > c \,|\, H_0) \approx ℙ(T' > c \,|\, \mathcal{D}),
```
where the bootstrap statistic ``T'`` is defined as
```math
T' = \frac{2}{n(n - 1)} \sum_{1 \leq i < j \leq n} h_k\big((P^*_{X_i}, Y^*_i), (P^*_{X_j}, Y^*_j)\big)
- \frac{2}{n^2} \sum_{i, r=1}^n h_k\big((P^*_{X_i}, Y^*_i), (P_{X_r}, Y_r)\big)
```
for bootstrap samples ``(P^*_{X_i}, Y^*_i)_{i=1,\ldots,n}`` of ``\mathcal{D}``
(see [`AsymptoticSKCETest`](@ref)).

Let ``C_i`` be the number of times that data pair ``(P_{X_i}, Y_i)`` was resampled.
Then we obtain
```math
T' = \frac{1}{n^2} \sum_{i=1}^n C_i \sum_{j=1}^n \bigg(\frac{n}{n-1} (C_j - \delta_{i,j}) - 2\bigg) h_k\big((P_{X_i}, Y_i), (P_{X_j}, Y_j)\big).
```
"""
function bootstrap_ccdf(
    rng::Random.AbstractRNG, statistic, kernelmatrix, bootstrap_iters::Int
)
    # initialize array of counts of resampled indices
    nsamples = LinearAlgebra.checksquare(kernelmatrix)
    resampling_counts = Vector{Int}(undef, nsamples)

    # for each bootstrap sample
    α = nsamples / (nsamples - 1)
    extreme_count = 0
    sampler = Random.Sampler(rng, 1:nsamples)
    for _ in 1:bootstrap_iters
        # resample data set
        fill!(resampling_counts, 0)
        for _ in 1:nsamples
            idx = rand(rng, sampler)
            @inbounds resampling_counts[idx] += 1
        end

        # evaluate the bootstrap statistic
        z = zero(statistic)
        n = 0
        for i in 1:nsamples
            # check if the `i`th data pair was sampled
            ci = resampling_counts[i]
            iszero(ci) && continue

            zi = Statistics.mean(enumerate(resampling_counts)) do (j, cj)
                # obtain evaluation of the kernel function
                @inbounds hij = kernelmatrix[j, i]
                return ((cj - (i == j)) * α - 2) * hij
            end

            # update bootstrap statistic
            n += ci
            z += ci * (zi - z) / n
        end

        # check if the bootstrap statistic is ≥ the original statistic
        if z ≥ statistic
            extreme_count += 1
        end
    end

    return extreme_count / bootstrap_iters
end


text:
"""
    struct CLagrangianFESpace{S} <: SingleFieldFESpace
      grid::Grid
      dof_to_node::Vector{Int}
      dof_to_comp::Vector{Int8}
      node_and_comp_to_dof::Vector{S}
      # + private fields
    end
"""
struct CLagrangianFESpace{S} <: SingleFieldFESpace
  space::UnconstrainedFESpace
  grid::Grid
  dof_to_node::Vector{Int}
  dof_to_comp::Vector{Int8}
  node_and_comp_to_dof::Vector{S}

  @doc """
      CLagrangianFESpace(::Type{T},grid::Grid) where T
  """
  function CLagrangianFESpace(::Type{T},grid::Grid) where T
    space, dof_to_node, dof_to_comp, node_and_comp_to_dof = _generate_clagrangian_fespace(zero(T),grid)
    S = eltype(node_and_comp_to_dof)
    new{S}(space,grid,dof_to_node,dof_to_comp,node_and_comp_to_dof)
  end
end

# FESpace interface

ConstraintStyle(::Type{<:CLagrangianFESpace}) = UnConstrained()

function get_free_dof_ids(f::CLagrangianFESpace)
  get_free_dof_ids(f.space)
end

function get_cell_shapefuns(f::CLagrangianFESpace)
  get_cell_shapefuns(f.space)
end

function get_cell_shapefuns_trial(f::CLagrangianFESpace)
  get_cell_shapefuns_trial(f.space)
end

get_vector_type(f::CLagrangianFESpace) = get_vector_type(f.space)

get_dof_value_type(f::CLagrangianFESpace) = get_dof_value_type(f.space)

get_triangulation(f::CLagrangianFESpace) = get_triangulation(f.space)

# SingleFieldFESpace interface

function get_cell_dof_ids(f::CLagrangianFESpace)
  get_cell_dof_ids(f.space)
end

function get_cell_dof_basis(f::CLagrangianFESpace)
  get_cell_dof_basis(f.space)
end

function get_dirichlet_dof_ids(f::CLagrangianFESpace)
  get_dirichlet_dof_ids(f.space)
end

function num_dirichlet_tags(f::CLagrangianFESpace)
  num_dirichlet_tags(f.space)
end

function get_dirichlet_dof_tag(f::CLagrangianFESpace)
  get_dirichlet_dof_tag(f.space)
end

function scatter_free_and_dirichlet_values(f::CLagrangianFESpace,free_values,dirichlet_values)
  scatter_free_and_dirichlet_values(f.space,free_values,dirichlet_values)
end

function gather_free_and_dirichlet_values!(free_values,dirichlet_values,f::CLagrangianFESpace,cell_vals)
  gather_free_and_dirichlet_values!(free_values,dirichlet_values,f.space,cell_vals)
end

function gather_dirichlet_values!(dirichlet_vals,f::CLagrangianFESpace,cell_vals)
    gather_dirichlet_values!(dirichlet_vals,f.space,cell_vals)
end

# Helpers

function _generate_clagrangian_fespace(z,grid)

  grid_reffes = get_reffes(grid)
  cell_to_ctype = get_cell_type(grid)
  function fun(reffe)
    p = get_polytope(reffe)
    orders = get_orders(reffe)
    LagrangianRefFE(typeof(z),p,orders)
  end
  # using map instead of lazy_map seems to kill the compiler.
  ctype_to_reffe = collect(lazy_map(fun,grid_reffes))
  cell_to_reffe = expand_cell_data(ctype_to_reffe,cell_to_ctype)

  nnodes = num_nodes(grid)
  dof_to_node, dof_to_comp, node_and_comp_to_dof = _generate_dof_layout_component_major(z,nnodes)

  cell_node_ids = get_cell_node_ids(grid)

  cell_dofs_ids = _generate_cell_dofs_clagrangian_fespace(
    z,
    cell_node_ids,
    ctype_to_reffe,
    cell_to_ctype,
    node_and_comp_to_dof)

  cell_map = get_cell_map(grid)
  cell_shapefuns, cell_dof_basis = compute_cell_space(cell_to_reffe, grid)

  nfree = length(dof_to_node)
  ndirichlet = 0

  dirichlet_dof_tag = Int8[]
  dirichlet_cells = Int32[]
  cell_is_dirichlet = Fill(false,num_cells(grid))
  ntags = 0
  vector_type = Vector{_dof_type(z)}

  space = UnconstrainedFESpace(
    vector_type,
    nfree,
    ndirichlet,
    cell_dofs_ids,
    cell_shapefuns,
    cell_dof_basis,
    cell_is_dirichlet,
    dirichlet_dof_tag,
    dirichlet_cells,
    ntags)

  space, dof_to_node, dof_to_comp, node_and_comp_to_dof

end

_dof_type(r::Number) = typeof(r)
_dof_type(r::MultiValue) = eltype(r)

function _generate_dof_layout_component_major(z,nnodes::Integer)
  ndofs = nnodes
  dof_to_comp = ones(Int8,ndofs)
  dof_to_node = collect(1:nnodes)
  node_and_comp_to_dof = dof_to_node
  (dof_to_node, dof_to_comp, node_and_comp_to_dof)
end

function _generate_dof_layout_component_major(z::MultiValue,nnodes::Integer)
  ncomps = num_components(z)
  V = change_eltype(z,Int)
  ndofs = ncomps*nnodes
  dof_to_comp = zeros(Int8,ndofs)
  dof_to_node = zeros(Int,ndofs)
  node_and_comp_to_dof = zeros(V,nnodes)
  m = zero(Mutable(V))
  for node in 1:nnodes
    o = ncomps*(node-1)
    for comp in 1:ncomps
      dof = comp+o
      dof_to_comp[dof] = comp
      dof_to_node[dof] = node
      m[comp] = dof
    end
    node_and_comp_to_dof[node] = m
  end
  (dof_to_node, dof_to_comp, node_and_comp_to_dof)
end

function _generate_cell_dofs_clagrangian_fespace(
  z,
  cell_node_ids,
  reffes,
  cell_to_ctype,
  node_and_comp_to_dof)

  cell_node_ids
end

function _generate_cell_dofs_clagrangian_fespace(
  z::MultiValue,
  cell_node_ids,
  reffes,
  cell_to_ctype,
  node_and_comp_to_dof)

  ncomps = num_components(z)

  ctype_to_lnode_to_comp_to_ldof = map(get_node_and_comp_to_dof,reffes)
  ctype_to_num_ldofs = map(num_dofs,reffes)

  cell_dofs = _allocate_cell_dofs_clagrangian_fespace(ctype_to_num_ldofs,cell_to_ctype)

  cache = array_cache(cell_node_ids)
  _fill_cell_dofs_clagrangian_fespace!(
    cell_dofs,
    cache,
    cell_node_ids,
    ctype_to_lnode_to_comp_to_ldof,
    cell_to_ctype,
    node_and_comp_to_dof,
    ncomps)

  cell_dofs
end

function _allocate_cell_dofs_clagrangian_fespace(ctype_to_num_ldofs,cell_to_ctype)
  ncells = length(cell_to_ctype)
  ptrs = zeros(Int32,ncells+1)
  for cell in 1:ncells
    ctype = cell_to_ctype[cell]
    nldofs = ctype_to_num_ldofs[ctype]
    ptrs[cell+1] = nldofs
  end
  length_to_ptrs!(ptrs)
  ndata = ptrs[end]-1
  data = zeros(Int32,ndata)
  Table(data,ptrs)
end

function _fill_cell_dofs_clagrangian_fespace!(
  cell_dofs,
  cache,
  cell_node_ids,
  ctype_to_lnode_and_comp_to_ldof,
  cell_to_ctype,
  node_and_comp_to_dof,
  ncomps)

  ncells = length(cell_to_ctype)

  for cell in 1:ncells
    nodes = getindex!(cache,cell_node_ids,cell)
    ctype = cell_to_ctype[cell]
    lnode_and_comp_to_ldof = ctype_to_lnode_and_comp_to_ldof[ctype]
    p = cell_dofs.ptrs[cell]-1
    for (lnode, node) in enumerate(nodes)
      for comp in 1:ncomps
        ldof = lnode_and_comp_to_ldof[lnode][comp]
        dof = node_and_comp_to_dof[node][comp]
        cell_dofs.data[p+ldof] = dof
      end
    end
  end

end

ext: jl


text:
# define a function that returns a Plots.Shape
function rectangle(w, h, x, y)

    return Shape(x .+ [0, w, w, 0], y .+ [0, 0, h, h])

end


text:
"""
TESTIFY
"""
module Client

using Dates
using LinearAlgebra # for dot etc...
using Distributed
using DataStructures

using cywwtaip  # HTW lib from AI professor

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

outputBacktrace(ex) = open(FILE_ERROR,"w+") do f Base.showerror(f, ex, catch_backtrace()) end

const ROOT = abspath(joinpath(@__DIR__,"../"))
const LOGS = abspath(joinpath(ROOT,"logs/"))
const FILE_ERROR = abspath(joinpath(LOGS, "error.log"))

name = nothing
winMsg  = nothing
server = nothing
client = nothing

latency = 0
playerNumber = 1 # 0 = nobody, 1 - 3 a player

mutable struct AIBot
  parent::Bot
  node::Union{Nothing, GraphNode}
  startnode::Union{Nothing, GraphNode}
  prevnode::Union{Nothing, GraphNode}
  target::Union{Nothing, GraphNode}
  targetEnergyField::Union{Nothing, GraphNode}
  movespeed::Float32

  AIBot(parent::Bot, node::Union{Nothing, GraphNode}) = new(parent,node,nothing,nothing,nothing,nothing,0)
end

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

getPlayerbotsTemplate() = Array{Union{Nothing,Array{Union{Nothing,Bot},1}},1}(undef,3)
getBotsTemplate() = Array{Union{Nothing,Bot},1}(undef,3)
getAIBotsTemplate() = Array{Union{Nothing,AIBot},1}(undef,3)

world = getWorldTemplate()
bots=getBotsTemplate()
playerbots=getPlayerbotsTemplate()
AIBots = getAIBotsTemplate()

worldlist=Dict{Symbol,Integer}()
listOfEnergyNodes = Array{GraphNode,1}()
listOfEnergyFields = Array{GraphNode,1}()
listOfWaypoints= Array{Array{GraphNode,1},1}()

ENERGYTIME = 0
ENERGYAREA = 0.94

lastNodesMax = 100
lastNodes = []

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

function resetWorld()
  open(FILE_ERROR,"w+") do f write(f,"") end

  global world, worldlist, playerbots, bots, AIBots
  world = getWorldTemplate()
  playerbots = getPlayerbotsTemplate()
  bots = getBotsTemplate()
  worldlist = Dict{Symbol,Integer}()
  AIBots = getAIBotsTemplate()

  fill!(world,nothing)
  fill!(playerbots,nothing)
  fill!(bots,nothing)
  fill!(AIBots,nothing)

  global listOfEnergyNodes = Array{GraphNode,1}()
  global listOfEnergyFields = Array{GraphNode,1}()
  global listOfWaypoints= [GraphNode[],GraphNode[],GraphNode[]]

  list1 = Array{Union{Nothing,GraphNode},1}(undef,  lastNodesMax)
  list2 = Array{Union{Nothing,GraphNode},1}(undef, lastNodesMax)
  list3 = Array{Union{Nothing,GraphNode},1}(undef, lastNodesMax)

  for i=1:lastNodesMax list1[i] = nothing end
  for i=1:lastNodesMax list2[i] = nothing end
  for i=1:lastNodesMax list3[i] = nothing end

  global lastNodes = [[0, list1],[0, list2],[0, list3]]
end

toKey(list::AbstractArray) = Symbol(join(list, ','))

getNode(id::Integer) = id>0 && id<=length(world) ? world[id] : nothing

function getNode(position::Array{Float32,1})
  key = toKey(position)
  haskey(worldlist, key) ? getNode(worldlist[key]) : nothing
end

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

distance(p1::Position,p2::Position) = √((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2 + (p2[3] - p1[3])^2)
distance(p1::GraphNode,p2::GraphNode) = distance(p1.position,p2.position)
distance(p1::GraphNode,p2::Position) = distance(p1.position,p2)
distance(p1::Position,p2::GraphNode) = distance(p1,p2.position)

# https://stackoverflow.com/questions/14066933/direct-way-of-computing-clockwise-angle-between-2-vectors
# https://stackoverflow.com/questions/5188561/signed-angle-between-two-3d-vectors-with-same-origin-within-the-same-plane
function getAngle(Va::Position,Vb::Position,Vn=[0,1,0]; mul=1)::Float32
  #angle = atan(Vn ⋅ (Va × Vb), Va ⋅ Vb)
  #angle = atan(Vn ⋅ (Va × Vb) / (Va ⋅ Vb))
  #return angle

  vcross = Va × Vb
  vdot = (Va ⋅ Vb) * (1.0/mul)
  sign = Vn ⋅ vcross

  normed = norm(Va) * norm(Vb)
  sina = norm(vcross) / normed
  cosa = vdot / normed
  angle = atan(sina, cosa)
  if sign<0 angle=-angle end

  angle
end

# http://geomalgorithms.com/a05-_intersect-1.html
# https://stackoverflow.com/questions/5666222/3d-line-plane-intersection
function getIntersectionPlanePoint(pos::Position,target::Position,up::Position)::Array{Float32,1}
  n = up
  V0 = pos
  P0 = target
  P1 = (target+n)
  w = P0 - V0 # ok
  u = P1 - P0 # ok
  s = -(n ⋅ w) / (n ⋅ u) # ok
  p = P0 + s*u # ok
end

function localAngleToTarget(pos::Position,target::Position,forward::Position=[0,0,1]; mul=1)::Float32
  global_up = up = normalize(pos)
  #right = forward × global_up
  #up = right × forward #
  #right = (forward × up) #* -1
  intersect = getIntersectionPlanePoint(pos,target,up)
  getAngle(pos,intersect,forward*-1; mul=mul)
end

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

function moveAI(client::NetworkClient, aibot::AIBot)
  target = aibot.target == nothing ? aibot.startnode : aibot.target # aibot.targetEnergyField
  #dist = distance(aibot.parent.position,target)
  #if dist < 0.051 println("$(aibot.parent.id): $(dist)") end
  # 0.024, 0.016, 0.01
  mul = 1 #((1000 * aibot.parent.speed)/dist)
  if mul < 1 mul=1 end

  position = getBotPosition(client, playerNumber, aibot.parent.id)
  direction = getBotDirection(client, aibot.parent.id)
  angle = localAngleToTarget(position, target.position, direction; mul=mul)
  #println("$(aibot.parent.id): "*string(angle))
  #println("$(aibot.parent.id): $(aibot.parent.speed): $(aibot.movespeed)")
  changeMoveDirection(client, aibot.parent.id, angle)
end

function moveAI(client::NetworkClient, aibot::AIBot, target::Position)
  mul = 100
  if mul < 1 mul=1 end
  position = getBotPosition(client, playerNumber, aibot.parent.id)
  direction = getBotDirection(client, aibot.parent.id)
  angle = localAngleToTarget(position, target, direction; mul=mul)
  changeMoveDirection(client, aibot.parent.id, angle)
end

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

upgradeToAIBots(bots::Array{Bot,1},nodes::Array{GraphNode,1}) = global AIBots = [AIBot(bot,node) for (bot,node) in zip(bots,nodes)]

function updateBots(client::NetworkClient)
  #cywwtaip.updateBots(client)
  #cywwtaip.updateBots(client, bots)

  for aibot in AIBots
    prevposition = aibot.parent.position
    start=time()
    cywwtaip.updateBot(client, aibot.parent)
    latency=time()-start
    aibot.movespeed = min(norm(prevposition - aibot.parent.position),1)
    aibot.parent.position += aibot.parent.direction*latency # adjust to future position because we have latency

    aibot.prevnode = nothing
    aibot.node = best = searchCurrentNode(aibot)
    if aibot.prevnode == nothing aibot.prevnode = best end
    #aibot.node=getNode(aibot.parent.position)
    if aibot.node == nothing @warn "bot $(aibot.parent.id): no start node found!" end
    if aibot.startnode == nothing aibot.startnode = aibot.node end
    searchNext(aibot)
    moveAI(client, aibot)
  end
end

function saveClosestNode(bots, nodes, current)
  for bot in bots
    prev = nodes[bot.id]
    if prev == nothing || (prev != nothing && distance(bot.position,current) < distance(bot.position, prev))
      nodes[bot.id] = current
    end
  end
end

function createWorld(client::NetworkClient)
  global world, worldlist, bots, playerbots

  bots = cywwtaip.getBots(client, playerNumber)
  nodes = Union{Nothing,GraphNode}[nothing,nothing,nothing]

  world = getGraph(client)
  if length(world) <= 0 return end

  f=world[1]
  bX=f; sX=f; bY=f; sY=f; bZ=f; sZ=f

  # save to worldlist
  i=0; for node in world
    i+=1
    node.id=i
    node.position
    worldlist[toKey(node.position)] = i
    (x,y,z) = node.position
    #println("owner($i):"*string(node.owner))

    # energy fields with abs max x,y or z
    if bX.position[1] < x bX=node end
    if sX.position[1] > x sX=node end
    if bY.position[2] < y bY=node end
    if sY.position[2] > y sY=node end
    if bZ.position[3] < z bZ=node end
    if sZ.position[3] > z sZ=node end

      # all energy field nodes
    if ((x>y&&x>z) && x>ENERGYAREA) || ((y>x&&y>z) && y>ENERGYAREA) || ((z>y&&z>x) && z>ENERGYAREA)
      push!(listOfEnergyNodes,node)
    end

    #node.position = cywwtaip.getPosition(node.handle)
    saveClosestNode(bots, nodes, node)
  end

  # save energy fields with abs max x,y or z
  push!(listOfEnergyFields,bX); push!(listOfEnergyFields,sX)
  push!(listOfEnergyFields,bY); push!(listOfEnergyFields,sY)
  push!(listOfEnergyFields,bZ); push!(listOfEnergyFields,sZ)

  # replace neighbors with GraphNodes
  for parent in world
    j=0; for node in parent.neighbors
      j+=1
      if isa(node,cywwtaip.JGraphNode) # only replace if type is wrong
        node = getNode(getPosition(node))
        node.parent = parent
        parent.neighbors[j] = node
      end
    end
  end

  for node in nodes
    if node == nothing error("Start Nodes not complete!") end
  end

  #playerbots = getBots(client)
  upgradeToAIBots(bots,Array{GraphNode,1}(nodes))
end

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

function searchCurrentNode(aibot::AIBot, position::Position)
  (old,best,bestdist) = (nothing,aibot.node, distance(aibot.node,position))

  for node in world
    dist = distance(node,position)
    if dist < bestdist; (old,best,bestdist) = (best,node,dist); end
  end

  best
end

function searchCurrentNode(aibot::AIBot)
  (old,best,bestdist) = (nothing,aibot.node, distance(aibot.node,aibot.parent.position))

  for node in world
    dist = distance(node,aibot.parent.position)
    if dist < bestdist; (old,best,bestdist) = (best,node,dist); end
  end

  best
end

function reconstruct_path_old!(cameFrom::AbstractDict, current::GraphNode)
  total_path = []
  for (_,node) in cameFrom
    push!(total_path,node)
  end
  push!(total_path,current)
  total_path
end

function reconstruct_path!(cameFrom::AbstractDict, current::GraphNode)
  total_path = [current]
  ks = keys(cameFrom)
  key = last(collect(ks))
  id = cameFrom[key]
  node = getNode(id)
  prepend!(total_path,[node])

  for i=1:length(cameFrom)
    if !haskey(cameFrom, id)
      @warn "Node not found in path?"
      break
    end

    old=id
    id = cameFrom[id]
    node = getNode(id)
    prepend!(total_path,[node])

    if old == id break end # found start
  end

  total_path
end

function getValue(set::SortedSet,key)
  for (fvalue,fkey) in set
    if fkey == key return fvalue end
  end
  nothing
end

function setValue!(set::SortedSet, key, value)
  rm=nothing
  for t in set
    (fvalue,fkey)=t
    if fkey == key
      if fvalue != value
        rm=t
      else
        return false # same entry?
      end
      break
    end
  end
  if rm != nothing delete!(set,rm) end
  push!(set,(value,key))
  true
end

# still looking for a perfect solution here...
function clustering()
 # energy field are already defined in the code
 # i calculate the distance at runtime, but i see now a pre-calculation with nodes might be better here
 # i have already a search node function with gets me the closest node for a bot
 # i need to connect clustering and aStarSearch somehow, to make it faster
 # position of a bot is already in relation with latency
end

# filter not owned fields and non blocked fields
filterFreeFields() = filter((x)->!x.owner.blocked && x.owner!=playerNumber,world)

#getNode(first(filter(x->in(x[2],openSet),collect(fScore)))[2])
function getNodeByLowestScore_old(openSet, fScore) #the node in openSet having the lowest fScore[] value
  for (_,node_id) in fScore
    if in(node_id,openSet) return getNode(node_id) end
  end
  @warn "No fScore node in openSet!"
  nothing #length(openSet)>0 ? getNode(first(openSet)) : nothing
end

function getLowestScore(openSet, fScore) #the node in openSet having the lowest fScore[] value
  for t in fScore
   if !in(t[2],openSet) continue end
   return t
  end
  #@warn "No fScore node in openSet!"
  nothing
end

#getNode(world, first(filter(x->in(x[2],openSet),collect(fScore)))[2])
function getNodeByLowestScore(openSet, fScore) #the node in openSet having the lowest fScore[] value
  t = getLowestScore(openSet, fScore)
  if t !=nothing return getNode(t[2]) end
  nothing #length(openSet)>0 ? getNode(world, first(openSet)) : nothing
end

function aStarSearch(start::GraphNode, goal::GraphNode)
  # The set of nodes already evaluated
  closedSet = Int32[]
  id=start.id

  # The set of currently discovered nodes that are not evaluated yet.
  # Initially, only the start node is known.
  openSet = Int32[id]

  # For each node, which node it can most efficiently be reached from.
  # If a node can be reached from many nodes, cameFrom will eventually contain the
  # most efficient previous step.
  cameFrom = OrderedDict{Int32,Int32}() # an empty map
  cameFrom[id] = id

  # For each node, the cost of getting from the start node to that node.
  gScore = SortedSet{Tuple{Float32,Int32}}() # map with default value of Infinity

  # The cost of going from start to start is zero.
  setValue!(gScore, id, 0)

  # For each node, the total cost of getting from the start node to the goal
  # by passing by that node. That value is partly known, partly heuristic.
  fScore = SortedSet{Tuple{Float32,Int32}}() # map with default value of Infinity

  # For the first node, that value is completely heuristic.
  setValue!(fScore, id, distance(start, goal)) # heuristic_cost_estimate = distance

  found = false
  stuck = false
  list = [start.id,goal.id]

  for i=1:length(world) #while length(openSet)>0 # avoid endless loop
    #sort!(s,alg=InsertionSort,lt=(x,y)->x[2]<y[2])

    current = getNodeByLowestScore(openSet,fScore) #the node in openSet having the lowest fScore[] value

    if current == nothing
      @warn "Cannot find a way to goal!"
      stuck = true
      break
    end

    if current == goal
      list = reconstruct_path!(cameFrom, current)
      found = true
      break
    end

    cid=current.id

    openSet = filter!(e->e ≠ cid, openSet) # delete
    #delete!(openSet, current.id)
    push!(closedSet, cid)

    for neighbor in current.neighbors # || neighbor.owner != playerNumber
      nid=neighbor.id
      if in(nid, closedSet) || neighbor.blocked continue end # Ignore the neighbor which is already evaluated or blocked.

      tscore = getValue(gScore,cid)
      if tscore == nothing
        tscore=0
        @warn "Error: Could not find Neighbor???"
        break
      end

      # The distance from start to a neighbor
      tentative_gScore = tscore + distance(current, neighbor) # dist_between = distance
      nscore = getValue(gScore,nid)

      if nscore != nothing && tentative_gScore >= nscore continue end
      if nscore == nothing nscore=0 end

      if !in(nid, openSet) push!(openSet, nid) end #neighbor not in openSet	# Discover a new node

    #  if !in(nid, openSet) #neighbor not in openSet	# Discover a new node
    #    push!(openSet, nid)
    #  elseif nscore != nothing && tentative_gScore >= nscore
    #    continue
    #  end

      # This path is the best until now. Record it!
      cameFrom[nid] = cid
      setValue!(gScore, nid, tentative_gScore)
      setValue!(fScore, nid, tentative_gScore*0 + distance(neighbor, goal)) # heuristic_cost_estimate = distance
    end
  end

  if !found && !stuck
    list = reconstruct_path!(cameFrom, goal)
    @warn "A*-Search Incomplete: $(length(list)) -> $(in(goal,list))"
  end

  #list = reconstruct_path!(cameFrom, goal)
  #@warn "a*-search result $(length(list)) -> $(in(goal,list))"
  list
end

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

function rechargeWhenLow(aibot)
  if (ENERGYTIME-time()) >= 9 return false end

  if aibot.targetEnergyField != nothing && distance(aibot.node, aibot.targetEnergyField) < 0.3
    aibot.targetEnergyField=nothing
    global ENERGYTIME = time()+10 # get next 10 sec
    return false
  end

  if aibot.targetEnergyField != nothing return true end

  best = (999,nothing)
  for node in listOfEnergyFields
    #node = getNode(index)
    dist = distance(aibot.node, node)
    if best[1] > dist best=(dist,node) end
  end
  #println(best[2].position)
  aibot.targetEnergyField = best[2] #rand(listOfEnergyFields)
  #aibot.target = aibot.targetEnergyField
  true
end

#bot 1 search and draw over enemy colors
#bot 2 when low: search energy felds suchen, else: draw non owned fields
#bot 3 draw all non owned fields

function searchForNotOwnedFields(aibot::AIBot, enemyOnly=false)
  current = aibot.node
  #println("$(aibot.parent.id): p:$playerNumber")
  for i=1:length(world) # do not go to infinity
    for node in current.neighbors
      if node.blocked || node.owner == playerNumber continue end
      if enemyOnly && node.owner == 0 continue end
      #println("$(aibot.parent.id): owner:$(node.owner)")
      return node
    end
    current=current.parent
  end
  nothing
end

function searchNext(aibot::AIBot)
  # run towards nearest energy target when low
  if aibot.parent.id == 2 && rechargeWhenLow(aibot) return end

  waypoints = listOfWaypoints[aibot.parent.id]

  if length(waypoints)>0
    global lastNodes, lastNodesMax
    counter = lastNodes[aibot.parent.id][1]
    list = lastNodes[aibot.parent.id][2]
    foundnode = searchCurrentNode(aibot, getBotPosition(client,playerNumber,aibot.parent.id))
    if foundnode == aibot.node
      counter+=1
      if counter >= lastNodesMax
        counter=0
        listOfWaypoints[aibot.parent.id]=[]
        #println("$(aibot.parent.id): got stuck!")
      end
    else counter=0
    end
    lastNodes[aibot.parent.id][1] = counter
  end

  if length(waypoints)==0
    colorTarget = searchForNotOwnedFields(aibot, aibot.parent.id == 1)
    if colorTarget==nothing colorTarget = first(listOfEnergyFields) end # default target
    waypoints = listOfWaypoints[aibot.parent.id] = aStarSearch(aibot.node,colorTarget)
    #println("$(aibot.parent.id): $(length(waypoints)) waitpoints!")
  end

  foundnode = searchCurrentNode(aibot, getBotPosition(client,playerNumber,aibot.parent.id))

  best = (999,0,nothing)
  for (index,waypoint) in enumerate(waypoints)
    index+=1
    dist = distance(foundnode, waypoint)
    if best[1]>dist best=(dist,index,waypoint) end
    if waypoint == foundnode best=(dist,index,waypoint); break end
  end

  if best[3] != nothing
    waypoints=listOfWaypoints[aibot.parent.id]=waypoints[best[2]:end]
    #println("$(aibot.parent.id): $(length(waypoints)) W")
  end

  if length(waypoints)==0
    waypoint = aibot.targetEnergyField
  else
    waypoint = first(waypoints)
  end

  aibot.target = waypoint
end

function updateWorld(client::NetworkClient)
  global world, bsave
  n=1+1
  updateGraph(client, world)
  updateBots(client)
end

function printWorld(printer::Union{Nothing,Function})
  if printer == nothing return end

str="""
# $(replace(string(now()),"T"=>" "))
-----------------------------------------
Position
 1:[$(string(join(bots[1].position,", ")))]
 2:[$(string(join(bots[2].position,", ")))]
 3:[$(string(join(bots[3].position,", ")))]

Direction
 1:[$(string(join(bots[1].direction,", ")))]
 2:[$(string(join(bots[2].direction,", ")))]
 3:[$(string(join(bots[3].direction,", ")))]

Speed
 1:[$(string(bots[1].speed))]
 2:[$(string(bots[2].speed))]
 3:[$(string(bots[3].speed))]
-----------------------------------------
"""
  @sync @distributed for node in world str*=toInfoString(node)*"\n" end
  str*="_________________________________________\n"

  printer(str)
end

####################################################################################################
using Revise
function main(args::Array{String,1}; printer::Union{Nothing,Function}=nothing)
  println("--------------------------------------------------------------------------------")
  @info "Start Game..."
  global playerNumber, latency

  len = length(args)
  global server = len>1 ? args[1] : "localhost" # 141.45.214.125
  global name = len>2 ? args[2] : "MyTeam"
  global winMsg = len>3 ? args[3] : "MyText"
  global client = nothing

  # Time limits
  serverTimeout = 1
  timeLimit = 60 # 60s
  timeLimitBeforeStart = 10 # 60s
  latency = 0

  timeLimitBroke = false
  prev_position = nothing
  DUMMY_CLIENTS = true
  startTime = 0
  tickTime = 0

  lstart=0
  lend=0

  resetWorld()

  try
    cywwtaip.init()
    cywwtaip.reset()

    @info "Connect to server '$server' (Timeout in $serverTimeout seconds)..."
    if DUMMY_CLIENTS
      try
        println("1. Dummy")
        client = NetworkClient(server, "Dummy1", "", 1)
      catch ex
      end
      try
        println("2. Dummy")
        client = NetworkClient(server, "Dummy2", "", 1)
      catch ex
      end
      println("Final Client")
    end
    client = NetworkClient(server, name, winMsg, serverTimeout)

    playerNumber = getMyPlayerNumber(client)
    createWorld(client) # create world first

    @info "Wait for start..."
    startTime = time()
    while isAlive(client)
      currentTime = (time() - startTime)

      lstart = time()
      position = getBotPositions(client, playerNumber)
      lend = time()

      if prev_position != nothing && prev_position != position break end
      prev_position = position

      if (time() - tickTime) > 1
        println("Seconds: $(timeLimitBeforeStart - round(currentTime))")
        tickTime=time()
      end

      if currentTime > timeLimitBeforeStart timeLimitBroke=true; break end # auto close
    end

    latency = lend-lstart
    startTime = lstart-1
    global ENERGYTIME = startTime + 10
    tickTime=0

    if !timeLimitBroke
      @info "Main Loop..."
      println("\n[START]\n------------------------------")

      while isAlive(client)
        currentTime = (time() - startTime)

        if (time() - tickTime) > 1
          #println("Seconds: $(timeLimit - round(currentTime)), latency: $latency\n")
          tickTime=time()
          #printWorld(printer) # 0.3
        end

        if currentTime > timeLimit timeLimitBroke=true; break end #after Xs break this loop

        updateWorld(client)
      end
    end

  catch ex
    outputBacktrace(ex)
    error("NetworkClient error! For more details look in $FILE_ERROR")
  end

  print("Game ended.")
  #cywwtaip.destroy() # cleanup cywwtaip
end

end # Client

max_stars_repo_name: Gilga/AI-Cywwtaip

max_issues_repo_name: Gilga/AI-Cywwtaip

max_forks_repo_name: Gilga/AI-Cywwtaip


text:
mutable struct RunningAvgPolicy{P<:AbstractPolicy} <: AbstractPolicy
    policy::P
    running_avg::Vector{Float32}
    smoothing::Float32
end

RunningAvgPolicy(; policy::AbstractPolicy, smoothing::Float32, na::Int) =
    RunningAvgPolicy(policy, zeros(Float32, na), smoothing)

function (p::RunningAvgPolicy)(env)
    a = p.policy(env)
    p.running_avg .= p.smoothing .* p.running_avg .+ (1 .- p.smoothing) .* a
end

function RLBase.update!(
    p::RunningAvgPolicy,
    traj::CircularArraySARTTrajectory,
    env::AbstractEnv,
    stage::PreActStage,
)
    update!(p.policy, traj, env, stage)
end

log_agent(p::RunningAvgPolicy) = log_agent(p.policy)

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
"""
    Hyperbolic{T <: Number} <: Construct{T}

A hyperbolic Cayley-Dickson construct as a pair.
"""
struct Hyperbolic{T <: Number} <: Construct{T}
    l::T
    r::T
end

function (==)(x::Hyperbolic, y::Hyperbolic)
    x.l == y.l && x.r == y.r
end

function real(z::Hyperbolic)
    real(z.l)
end

"""
    quadrance(z::CayleyDickson.Hyperbolic)

For hyperbolic Cayley-Dickson constructs, return the difference of the squared
absolute value of the left and right components of `z`.
Note that this can be positive, negative, or zero.
"""
function quadrance(z::Hyperbolic)
    quadrance(z.l) - quadrance(z.r)
end

function iszerodivisor(z::Hyperbolic)
    quadrance(z.l) == quadrance(z.r)
end

function zero(::Type{Hyperbolic{T}}) where T <: Number
    Hyperbolic{T}(zero(T), zero(T))
end

function zero(::Hyperbolic{T}) where T <: Number
    zero(Hyperbolic{T})
end

function one(::Type{Hyperbolic{T}}) where T <: Number
    Hyperbolic{T}(one(T), zero(T))
end

function one(::Hyperbolic{T}) where T <: Number
    one(Hyperbolic{T})
end

function (-)(z::Hyperbolic)
    Hyperbolic(-(z.l), -(z.r))
end

function conj(z::Hyperbolic)
    Hyperbolic(conj(z.l), -(z.r))
end

function (+)(x::Hyperbolic{T}, y::Hyperbolic{T}) where T <: Number
    Hyperbolic{T}(x.l + y.l, x.r + y.r)
end

function (+)(a::T, z::Hyperbolic{T}) where T <: Number
    Hyperbolic{T}(a + z.l, z.r)
end

function (+)(z::Hyperbolic{T}, a::T) where T <: Number
    Hyperbolic{T}(z.l + a, z.r)
end

function (+)(a::Real, z::Hyperbolic)
    l = a + z.l
    Hyperbolic(l, oftype(l, z.r))
end

function (+)(z::Hyperbolic, a::Real)
    l = z.l + a
    Hyperbolic(l, oftype(l, z.r))
end

function (-)(x::Hyperbolic{T}, y::Hyperbolic{T}) where T <: Number
    Hyperbolic{T}(x.l - y.l, x.r - y.r)
end

function (-)(a::T, z::Hyperbolic{T}) where T <: Number
    Hyperbolic{T}(a - z.l, -(z.r))
end

function (-)(z::Hyperbolic{T}, a::T) where T <: Number
    Hyperbolic{T}(z.l - a, z.r)
end

function (-)(a::Real, z::Hyperbolic)
    l = a - z.l
    Hyperbolic(l, -(oftype(l, z.r)))
end

function (-)(z::Hyperbolic, a::Real)
    l = z.l - a
    Hyperbolic(l, oftype(l, z.r))
end

function (*)(x::Hyperbolic{T}, y::Hyperbolic{T}) where T <: Number
    Hyperbolic{T}(
        (x.l * y.l) + (conj(y.r) * x.r),
        (y.r * x.l) + (x.r * conj(y.l))
    )
end

function (*)(a::T, z::Hyperbolic{T}) where T <: Number
    Hyperbolic{T}(a * z.l, z.r * a)
end

function (*)(z::Hyperbolic{T}, a::T) where T <: Number
    Hyperbolic{T}(z.l * a, z.r * conj(a))
end

function (*)(a::Real, z::Hyperbolic)
    Hyperbolic(a * z.l, a * z.r)
end

function (*)(z::Hyperbolic, a::Real)
    Hyperbolic(z.l * a, z.r * a)
end

function (/)(z::Hyperbolic, a::Real)
    if iszero(a)
        error(ZeroDenominator)
    end

    Hyperbolic(z.l / a, z.r / a)
end

function (\)(a::Real, z::Hyperbolic)
    if iszero(a)
        error(ZeroDenominator)
    end

    Hyperbolic(a \ z.l, a \ z.r)
end

function widen(::Type{Hyperbolic{T}}) where T
    Hyperbolic{widen(T)}
end

function float(::Type{Hyperbolic{T}}) where T <: AbstractFloat
    Hyperbolic{T}
end

function float(::Type{Hyperbolic{T}}) where T
    Hyperbolic{float(T)}
end

"""
    cloak(a::Hyperbolic)

For hyperbolic Cayley-Dickson constructs, equivalent to the minus operation.
"""
function cloak(a::Hyperbolic)
    -(a)
end

"""
    dagger(a::Hyperbolic)

For hyperbolic Cayley-Dickson constructs, equivalent to the identity operation.
"""
function dagger(a::Hyperbolic)
    a
end

"""
    star(a::Hyperbolic)

For hyperbolic Cayley-Dickson constructs, the identity operation.
"""
function star(a::Hyperbolic)
    a
end


text:
using SMLMData
using Optim
using Statistics

"""
    k_on, k_off, k_bleach, p_miss, nemitters =
        estimateparams(smld::SMLMData.SMLD2D, clusterdata::Vector{Matrix{Float32}})

Estimate rate parameters from the clusters in `smld` and `clusterdata`.

# Description
The rate parameters `k_on`, `k_off`, and `k_bleach`, the probability of missing
a localization for a visible emitter `p_miss`, and the underlying number of
emitters `nemitters` are estimated from the pre-clusters in
`smld`/`clusterdata`.  This is done by assuming that pre-clusters are, on
average, representative of a single blinking event of a single emitter.

# Citation
Schodt David J., Lidke Keith A., "Spatiotemporal Clustering of Repeated
Super-Resolution Localizations via Linear Assignment Problem",
Front. Bioinform., 20 October 2021
https://doi.org/10.3389/fbinf.2021.724325
"""
function estimateparams(smld::SMLMData.SMLD2D,
    clusterdata::Vector{Matrix{Float32}})

    # Compute some quantities needed from each of the preclusters.
    nclusters = size(clusterdata)[1]
    clusterdurations, nobservations = computeclusterinfo(clusterdata)

    # Estimate the sum of the off rate and the bleaching rate, as well as the
    # probability of missing a localization.
    k_offpbleach = -log(1 - 1 / Statistics.mean(clusterdurations))
    k_offpbleach = isinf(k_offpbleach) ? 1.0 : k_offpbleach
    p_miss = 1 - Statistics.mean(nobservations ./ clusterdurations)

    # Compute other parameters by fitting the observed number of localizations
    # over time (see citation for description of this math).
    nlocperframe = counts(smld.framenum)
    nlocperframe = nlocperframe[nlocperframe.!=0]
    frames = convert(Array{Float32}, unique(smld.framenum))
    nloccumulative = cumsum(nlocperframe)
    k_off = nclusters / nloccumulative[end]
    k_bleach = maximum([10^-5 k_offpbleach - k_off])
    k(k_on) = k_on + k_offpbleach
    l1(k_on) = k_on * k_bleach / k(k_on)
    l2(k_on) = k(k_on) - l1(k_on)
    model(x) = ceil(x[1]) * (1.0 .- p_miss) * (x[2] / k(x[2])) *
               ((1 / l1(x[2])) * (1.0 .- exp.(-l1(x[2]) * (frames .- 1.0))) .-
                (1 / l2(x[2])) * (1.0 .- exp.(-l2(x[2]) * (frames .- 1.0))))
    costfunction(x) = Statistics.mean((nloccumulative .- model(x)) .^ 2)
    lowerbound = [maximum(nlocperframe) 10^-5]
    upperbound = [nclusters nloccumulative[end] / frames[end]]
    x1_guess = ceil(nclusters * k_bleach)
    x1_guess = (x1_guess>lowerbound[1]) & (x1_guess<upperbound[1]) ?
               x1_guess : (upperbound[1]-lowerbound[1]) / 2.0
    initialguess = [x1_guess 1.0 / frames[end]]
    inner_optimizer = GradientDescent()
    optimresults = optimize(costfunction, lowerbound, upperbound, initialguess,
        Fminbox(inner_optimizer))
    xhat = optimresults.minimizer
    nemitters = xhat[1]
    k_on = xhat[2]

    return k_on, k_off, k_bleach, p_miss, nemitters
end


text:
module TestPresentation

using Test
using Catlab, Catlab.Doctrines

# Presentation
##############

A, B, C = Ob(FreeCategory, :A, :B, :C)
f = Hom(:f, A, B)
g = Hom(:g, B, C)

# Generators
pres = Presentation()
@test !has_generator(pres, :A)
add_generator!(pres, A)
@test generators(pres) == [ A ]
@test generator(pres, :A) == A
@test has_generator(pres, :A)
add_generator!(pres, B)
@test generators(pres) == [ A, B ]
@test_throws Exception add_generator!(pres, A)

add_generators!(pres, (f,g))
@test generators(pres) == [ A, B, f, g ]
@test generators(pres, FreeCategory.Ob) == [ A, B ]
@test generators(pres, FreeCategory.Hom) == [ f, g ]
@test generators(pres, (:A, :B)) == [ A, B ]
@test generators(pres, [:f, :g]) == [ f, g ]

# Merge presentations
pres = Presentation()
add_generators!(pres, (A, B, f))
other = Presentation()
add_generators!(other, (B, C, g))
merge_presentation!(pres, other)
@test generator(pres, :f) == f
@test generator(pres, :g) == g
@test generators(pres, FreeCategory.Ob) == [ A, B, C ]
@test generators(pres, FreeCategory.Hom) == [ f, g ]

# Presentation macro
####################

@present Company(FreeCategory) begin
  # Primitive concepts.
  Employee::Ob
  Department::Ob
  Str::Ob

  first_name::Hom(Employee, Str)
  last_name::Hom(Employee, Str)
  manager::Hom(Employee, Employee)
  works_in::Hom(Employee, Department)
  secretary::Hom(Department, Employee)

  # Defined concepts.
  second_level_manager := compose(manager, manager)
  third_level_manager := compose(manager, manager, manager)

  # Abbreviations (no syntactic term for LHS).
  boss = manager

  # Managers work in the same department as their employees.
  compose(boss, works_in) == works_in
  # The secretary of a department works in that department.
  compose(secretary, works_in) == id(Department)
end

# Check generators.
Employee, Department, Str = Ob(FreeCategory, :Employee, :Department, :Str)
@test collect(generators(Company)) == [
  Employee,
  Department,
  Str,
  Hom(:first_name, Employee, Str),
  Hom(:last_name, Employee, Str),
  Hom(:manager, Employee, Employee),
  Hom(:works_in, Employee, Department),
  Hom(:secretary, Department, Employee),
  Hom(:second_level_manager, Employee, Employee),
  Hom(:third_level_manager, Employee, Employee),
]

# Check equations.
manager = Hom(:manager, Employee, Employee)
works_in = Hom(:works_in, Employee, Department)
secretary = Hom(:secretary, Department, Employee)
@test collect(equations(Company)) == Equation[
  Hom(:second_level_manager, Employee, Employee) => compose(manager, manager),
  Hom(:third_level_manager, Employee, Employee) => compose(manager, manager, manager),
  compose(manager, works_in) => works_in,
  compose(secretary, works_in) => id(Department),
]

# Serialization
###############

to_json(expr) = to_json_sexpr(Company, expr)
from_json(sexpr) = parse_json_sexpr(Company, FreeCategory, sexpr)

# To JSON
to_json(generator(Company, :Employee)) == "Employee"
to_json(generator(Company, :manager)) == "manager"
to_json(compose(generator(Company, :manager), generator(Company, :manager))) ==
  ["compose", "manager", "manager"]

# From JSON
@test from_json("Employee") == generator(Company, :Employee)
@test from_json("manager") == generator(Company, :manager)
@test from_json(["compose", "manager", "manager"]) ==
  compose(generator(Company, :manager), generator(Company, :manager))

end


text:
using CairoMakie, Random, Dates, TimeSeries, MarketData, StatsBase
tempo = string.(timestamp(ohlc))
lentime = length(tempo)
slice_dates = range(1, lentime, step = 3)
colors = values(ohlc.Close) .> values(ohlc.Open)
lowV = values(ohlc.Low)
highV = values(ohlc.High)
linesegs = []
for i in 1:lentime
    push!(linesegs, Point2f0(i, lowV[i]))
    push!(linesegs, Point2f0(i, highV[i]))
end
linesegs = Point2f0.(linesegs)
cmap = ((:orangered, 0.85), (:dodgerblue, 0.85))

set_theme!(theme_dark())
fig = Figure(figure_padding = (25,15,60,15), resolution = (600, 400))
ax = Axis(fig[1,1]; xgridstyle=:dash, ygridstyle=:dash)
barplot!(ax, 1:lentime, values(ohlc.Open), fillto = values(ohlc.Close),
    width = 0.65, color = colors, strokewidth = 0.5, strokecolor = colors,
    colormap = cmap)
linesegments!(ax, linesegs; linewidth = 3, color = colors, colormap = cmap)
lines!(5:500, values(moving(mean, ohlc.Close, 5)); color = :yellow)
xlims!(ax, 194.5, 225.5)
ylims!(ax, 17,24)
ax.xticks = (slice_dates, tempo[slice_dates])
ax.xticklabelrotation = π/4
save("condleStick.png", fig)
fig


text:
function _load_config()
    if _is_calibration_loaded() # Not sure this is needed. Is this function ever called multiple times?
        @info "Calibration has already been loaded"
    else
        function_call_root_dir = @__DIR__

        config_file_relative_path = "calibration/config.jl"
        config_file_paths = String[]
        push!(config_file_paths, joinpath(pwd(), config_file_relative_path));
        push!(config_file_paths, joinpath(function_call_root_dir, config_file_relative_path))
        push!(config_file_paths, joinpath(function_call_root_dir, "..", config_file_relative_path))

        for file_path in config_file_paths
            if isfile(file_path)
                @info "Calibration directory found at " * file_path
                include(file_path)
                _set_calibration_loaded()
                return nothing # only load from the first found calibration directory
            end
        end
    end
    return nothing
end

function _is_calibration_loaded()::Bool
    global CALIBRATION_LOADED
    is_loaded = CALIBRATION_LOADED[]::Bool

    return is_loaded
end

function _set_calibration_loaded()
    global CALIBRATION_LOADED
    CALIBRATION_LOADED[] = true
    return nothing
end

function get_calibration(hydro_id, preamp_id = :empty)
    call_calibration(s) = Base.invokelatest(calibration, Val(s))
    return call_calibration(hydro_id), call_calibration(preamp_id)
end

macro calibration(device, calibration_data)
    @eval begin
        """
        `$($device)`
        """
        calibration(::Val{$device}) = $calibration_data
    end
end

calibration(::Val{:empty}) = []


text:
using Manopt
using Test
@testset "Manifold Tests     " begin
    include("manifolds/testGr.jl")
    include("manifolds/testGraph.jl")
    include("manifolds/testHn.jl")
    include("manifolds/testRn.jl")
    include("manifolds/testS1.jl")
    include("manifolds/testSn.jl")
    include("manifolds/testSO.jl")
    include("manifolds/testSPD.jl")
    include("manifolds/testStiefel.jl")
    include("manifolds/testSym.jl")
    include("manifolds/testGraphConstruction.jl")
    #
    include("manifolds/testManifold.jl")
    include("manifolds/testExtended.jl")
    include("manifolds/testCombined.jl")
end
@testset "Plan Tests         " begin
    include("plans/testDebug.jl")
    include("plans/testRecord.jl")
    include("plans/testGradientPlan.jl")
    include("plans/testSubGradientPlan.jl")
end
@testset "Algorithm Tests    " begin
    include("algorithms/testBasicsAlgs.jl")
end
@testset "Function Tests     " begin
    include("functions/testAdjointDifferentials.jl")
    include("functions/testDifferentials.jl")
    include("functions/testCosts.jl")
    include("functions/testGradients.jl")
    include("functions/testProximalMaps.jl")
end
@testset "Helper & Data Tests" begin
    include("helpers/testErrorMeasures.jl")
    include("helpers/testData.jl")
end
@testset "Solver Tests       " begin
    include("solvers/testDR.jl")
    include("solvers/testCPP.jl")
    include("solvers/testGradDesc.jl")
end


text:
using Libtask
using Turing, DynamicPPL, AdvancedPS
using BenchmarkTools

@model gdemo(x, y) = begin
    # Assumptions
    σ ~ InverseGamma(2,3)
    μ ~ Normal(0,sqrt(σ))
    # Observations
    x ~ Normal(μ, sqrt(σ))
    y ~ Normal(μ, sqrt(σ))
end

# Case 1: Sample from the prior.
m = Turing.Core.TracedModel(gdemo(1.5, 2.), SampleFromPrior(), VarInfo())
f = m.evaluator[1];
args = m.evaluator[2:end];

println("Directly call...")
@btime f(args...)
# (2.0, VarInfo (2 variables (μ, σ), dimension 2; logp: -6.162))
println("TapedTask construction...")
t = @btime TapedTask(f, args...)
println("Run a tape...")
@btime t.tf(args...)

# Case 2: SMC sampler
m = Turing.Core.TracedModel(gdemo(1.5, 2.), Sampler(SMC(50)), VarInfo());
f = m.evaluator[1];
args = m.evaluator[2:end];

println("Directly call...")
@btime f(args...)
println("TapedTask construction...")
t = @btime TapedTask(f, args...)
println("Run a tape...")
@btime t.tf(args...)

max_stars_repo_name: TuringLang/Libtask

max_issues_repo_name: TuringLang/Libtask

max_forks_repo_name: TuringLang/Libtask


text:
precompile(DirectImage, (Matrix{Float32},))
precompile(DirectImage, (Matrix{Float64},))
precompile(DirectImage, (Matrix{ComplexF32},))
precompile(DirectImage, (Matrix{ComplexF64},))

max_stars_repo_name: sefffal/SpaceImages


text:
#TODO: convert two function calls into Union{COOTen,ThirdOrderSymTensor}
#TODO: update new DistributedTensorConstructio types in docs
struct TAME_Return{T} <: returnType
	matchScore::Int
	motifCounts::Tuple{Int,Int}
	matching::Union{Dict{Int,Int},Vector{Int}}
	embedding::Matrix{T}
	profile::Union{Nothing,Vector{Tuple{String,Dict{String,Union{Array{Float64,1},Array{Array{Float64,1},1}}}}}}
end

#=------------------------------------------------------------------------------
              Routines for searching over alpha/beta parameters
------------------------------------------------------------------------------=#

#add in SparseSymmetricTensors.jl function definitions
function TAME_param_search(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor;
                           iter::Int = 15,tol::Float64=1e-6, alphas::Array{F,1}=[.5,1.0],
						   betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
						   kwargs...) where {F <: AbstractFloat}
	A_motifs = size(A.indices,2)
	B_motifs = size(B.indices,2)
	max_triangle_match = min(A_motifs,B_motifs)
    total_triangles = size(A.indices,2) + size(B.indices,2)
	best_TAME_PP_tris::Int = -1
	best_matching = Dict{Int,Int}()

	m = A.n
	n = B.n

	best_TAME_PP_x = Array{Float64,2}(undef,m,n)

    for α in alphas
        for β in betas

			x, triangle_count, matching = TAME(A,B,β,iter,tol,α;W = ones(A.n,B.n),kwargs...)

			if triangle_count > best_TAME_PP_tris
				best_matching = matching
                best_TAME_PP_tris = triangle_count
				best_TAME_PP_x = copy(x)
            end
			println("α:$(α) -- β:$β finished -- tri_match:$(triangle_count) -- max_tris $(max_triangle_match) -- best tri_match:$(best_TAME_PP_tris)")
        end

    end

	return TAME_Return(best_TAME_PP_tris,(A_motifs,B_motifs),best_matching,best_TAME_PP_x,nothing)
end

function TAME_param_search_profiled(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor;
                                    iter::Int = 15,tol::Float64=1e-6, alphas::Array{F,1}=[.5,1.0],
						            betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
						            profile::Bool=false,profile_aggregation="all",
						            kwargs...) where {F <: AbstractFloat}
	A_motifs = size(A.indices,2)
	B_motifs = size(B.indices,2)
	max_triangle_match = min(A_motifs,B_motifs)
    total_triangles = size(A.indices,2) + size(B.indices,2)
	best_TAME_PP_tris = -1
	best_matching = Dict{Int,Int}()

	m = A.n
	n = B.n

	best_TAME_PP_x = Array{Float64,2}(undef,m,n)
	experiment_profiles = Array{Tuple{String,Dict{String,Union{Array{Float64,1},Array{Array{Float64,1},1}}}},1}(undef,0)

    for α in alphas
        for β in betas

			x, triangle_count, matching, experiment_profile = TAME_profiled(A,B,β,iter,tol,α;W = ones(m,n),kwargs...)
			push!(experiment_profiles,("α:$(α)_β:$(β)",experiment_profile))

			if triangle_count > best_TAME_PP_tris
				best_matching = matching
                best_TAME_PP_tris = triangle_count
				best_TAME_PP_x = copy(x)
            end
			println("α:$(α) -- β:$β finished -- tri_match:$(best_TAME_PP_tris) -- max_tris $(max_triangle_match)")
        end

    end

	return TAME_Return(best_TAME_PP_tris,(A_motifs,B_motifs),best_matching,best_TAME_PP_x,experiment_profiles)

end

function TAME_param_search(A::SymTensorUnweighted{M},B::SymTensorUnweighted{M};
                           iter::Int = 15,tol::Float64=1e-6, alphas::Array{F,1}=[.5,1.0],
						   betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
						   kwargs...) where {F <: AbstractFloat,M <:Motif}
	A_motifs = size(A.indices,2)
	B_motifs = size(B.indices,2)
	max_motif_match = min(A_motifs,B_motifs)
	best_TAME_PP_motif::Int = -1
	best_matching = Dict{Int,Int}()

	m = A.n
	n = B.n

	best_TAME_PP_x = Array{Float64,2}(undef,m,n)

    for α in alphas
        for β in betas

			x, motif_count, matching = TAME(A,B,β,iter,tol,α;W = ones(A.n,B.n),kwargs...)

			if motif_count > best_TAME_PP_motif
				best_matching = matching
                best_TAME_PP_motif = motif_count
				best_TAME_PP_x = copy(x)
            end
			println("α:$(α) -- β:$β finished -- tri_match:$(motif_count) -- max_tris $(max_motif_match) -- best tri_match:$(best_TAME_PP_motif)")
        end

    end

	return TAME_Return(best_TAME_PP_motif,(A_motifs,B_motifs),best_matching,best_TAME_PP_x,nothing)
end

function TAME_param_search_profiled(A::SymTensorUnweighted{M},B::SymTensorUnweighted{M};
                                    iter::Int = 15,tol::Float64=1e-6, alphas::Array{F,1}=[.5,1.0],
						            betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
						            profile::Bool=false,profile_aggregation="all",
						            kwargs...) where {F <: AbstractFloat, M <: Motif}
	A_motifs = size(A.indices,2)
	B_motifs = size(B.indices,2)
	max_motif_match = min(A_motifs,B_motifs)
    total_motif = size(A.indices,2) + size(B.indices,2)
	best_TAME_PP_motif = -1
	best_matching = Dict{Int,Int}()

	m = A.n
	n = B.n

	best_TAME_PP_x = Array{Float64,2}(undef,m,n)
	experiment_profiles = Array{Tuple{String,Dict{String,Union{Array{Float64,1},Array{Array{Float64,1},1}}}},1}(undef,0)

    for α in alphas
        for β in betas

			x, motif_count, matching, experiment_profile = TAME_profiled(A,B,β,iter,tol,α;W = ones(m,n),kwargs...)
			push!(experiment_profiles,("α:$(α)_β:$(β)",experiment_profile))

			if motif_count > best_TAME_PP_motif
				best_matching = matching
                best_TAME_PP_motif = motif_count
				best_TAME_PP_x = copy(x)
            end
			println("α:$(α) -- β:$β finished -- tri_match:$(best_TAME_PP_motif) -- max_motifs $(max_motif_match)")
        end

    end

	return TAME_Return(best_TAME_PP_motif,(A_motifs,B_motifs),best_matching,best_TAME_PP_x,experiment_profiles)

end

#=------------------------------------------------------------------------------
             		    Spectral Relaxation Routines
------------------------------------------------------------------------------=#

function setup_tame_data(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor)
	return _index_triangles_nodesym(A.n,A.indices), _index_triangles_nodesym(B.n,B.indices)
end

function setup_tame_data(A::SymTensorUnweighted{M},B::SymTensorUnweighted{M}) where {M <: Motif}
	return _index_motifs_nodesym(A.n,A.indices), _index_motifs_nodesym(B.n,B.indices)
end

function _index_triangles_nodesym(n,Tris::Array{Int,2})

	Ti = [ Vector{Tuple{Int,Int}}(undef, 0) for i in 1:n ]
	for (ti,tj,tk) in eachcol(Tris)
		push!(Ti[ti], (tj,tk))
		push!(Ti[tj], (ti,tk))
		push!(Ti[tk], (ti,tj))
	end
	sort!.(Ti)
	return Ti
end

function _index_motifs_nodesym(n,Motifs::Array{Int,2})

	free_modes = collect(size(Motifs,1):-1:1)
	contracted_modes = collect(combinations(1:size(Motifs,1),size(Motifs,1)-1))

	Mi = [Vector{Vector{Int}}(undef, 0) for i in 1:n ]
	for edge in eachcol(Motifs)
		for (free_mode,contracted_mode) in zip(free_modes,contracted_modes)
			push!(Mi[edge[free_mode]], edge[contracted_mode])
		end

	end
	sort!.(Mi)
	return Mi
end

"""------------------------------------------------------------------------------
  An implementation of TAME, for our experiments we use the original C++ code,
  but we include these routines for any additional experimentation desired.

  Inputs
  ------
  * A,B - (ThirdOrderSymTensor):
	The tensors to align against one another.
  * W - (Array{Float,2}):
	The starting iteration, Iterate is normalized before the iterations begin.
  * β - (Float):
	The shift to use on the iterations.
  * α -(float):
	The mixing parameter for combining in the starting iterates.
  * 'max_iter' - (Int):
	The maximum number of iterations to run.
  * tol - (Float):
	The tolerence to solve the algorithm to. Computes the tolerance by measuring
	the absolute value of the difference between the computed eigenvalues.
  * 'update_user' - (Int):
	Specifies how frequently output messages should be printed to the user.
	Default is -1 which means no output messages are printed, else if the value
	is k, then the kth iterate will print out an update.
  * 'no_matching' - (Bool):
	Specifies whether or not to run the matching and scoring portions of the
	algorithm. Useful if only the iterates are desired.

  Output
  ------
	* 'best_x'- (Array{Float,2})
	  Returns the components to the iteration which matched the most triangles.
	  Reshapes the iterate x into a matrix.
	* 'best_motif_count' - (Int)
	  The maximum number of triangles matched.
	* 'best_matching' - (Dict{Int,Int})
	  The matching computed between the two graphs, maps the vertices of A to the
	  vertices of B.
------------------------------------------------------------------------------"""
function TAME(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,β::F, max_iter::Int,
			  tol::F,α::F;update_user::Int=-1,W::Array{F,2}=ones(A.n,B.n),
			  no_matching=false,kwargs...) where {F <:AbstractFloat}

    dimension = minimum((A.n,B.n))

	A_Ti, B_Ti = setup_tame_data(A,B)

	best_x = Array{Float64,2}(undef,A.n,B.n)
    best_motif_count = -1
	best_matching = Dict{Int,Int}()

    x0 = reshape(W,A.n*B.n)
	x0 ./=norm(x0)
    x_k = copy(x0)

    i = 1
    lambda = Inf

    while true

	    x_k_1 = impTTVnodesym(A.n, B.n, x_k, A_Ti, B_Ti)
		#x_k_1 = implicit_contraction(A,B,x_k)

		#println("norm diff is:",norm(x_k_1_test - x_k_1)/norm(x_k_1_test))
        new_lambda = dot(x_k_1,x_k)

        if β != 0.0
            x_k_1 .+= β * x_k
        end

        if α != 1.0
            x_k_1 = α * x_k_1 + (1 - α) * x0
        end

		x_k_1 ./= norm(x_k_1)

		if !no_matching

			X = reshape(x_k_1,A.n,B.n)
			triangles, missed_motifs, matching =  TAME_score(A,B,X;kwargs...)

			if update_user != -1 && i % update_user == 0
				println("finished iterate $(i):tris:$triangles -- gaped_t:$missed_motifs")
			end

			if triangles > best_motif_count
				best_matching = matching
				best_x = copy(x_k_1)
				best_motif_count = triangles
				best_iterate = i
			end

		end

		if update_user != -1 && i % update_user == 0
			println("λ: $(new_lambda)")
		end

        if abs(new_lambda - lambda) < tol || i >= max_iter
  			return reshape(best_x,A.n,B.n), best_motif_count, best_matching
        else
            x_k = copy(x_k_1)
            lambda = new_lambda
            i += 1
        end

    end

end

function TAME(A::SymTensorUnweighted{M}, B::SymTensorUnweighted{M}, β::F, max_iter::Int,
	                   tol::F,α::F;update_user::Int=-1,W::Array{F,2} = ones(A.n,B.n),
					   no_matching::Bool=false,kwargs...) where {F <:AbstractFloat, M <: Motif}

	A_Mi, B_Mi = setup_tame_data(A,B)

	best_X = Array{Float64,2}(undef,A.n,B.n)
    best_motif_count::Int = -1
	best_matching = Dict{Int,Int}()

    X_k = copy(W)
	X_k ./=norm(X_k)

	if α != 1.0
		X_0_with_alpha = copy(X_k)*(1-α)
	end

    i = 1
    lambda = Inf

    while true

		X_k_1 = impTTVnodesym(A.n, B.n, size(A.indices,1), X_k, A_Mi, B_Mi)

        new_lambda = dot(X_k_1,X_k)

        if β != 0.0
            X_k_1 .+= β * X_k
        end

        if α != 1.0
            X_k_1 .*=  α
			X_k_1 .+= X_0_with_alpha
        end

        X_k_1 ./= norm(X_k_1)

		if !no_matching

			matched_motifs, missed_motifs, matching = TAME_score(A,B,X_k_1;kwargs...)

			if update_user != -1 && i % update_user == 0
				println("finished iterate $(i): matched motifs:$(matched_motifs) -- missed motifs:$(missed_motifs)")
			end

			if matched_motifs > best_motif_count
				best_matching = matching
				best_X = copy(X_k_1)
				best_motif_count = matched_motifs
			end
		end

		if update_user != -1 && i % update_user == 0
			println("λ_$i: $(new_lambda)")
		end

        if abs(new_lambda - lambda) < tol || i >= max_iter
   			return best_X, best_motif_count, best_matching
        else
            X_k = copy(X_k_1)
            lambda = new_lambda
            i += 1
        end

    end

end

"""------------------------------------------------------------------------------
  The low rank implementation of TAME, computes the terms use the mixed property
  generalized to the tensor case.

  Inputs
  ------
  * A,B - (ThirdOrderSymTensor):
	The tensors to align against one another.
  * W - (Array{Float,2}):
	The starting iteration, Iterate is normalized before the iterations begin.
  * β - (Float):
	The shift to use on the iterations.
  * α -(float):
	The mixing parameter for combining in the starting iterates.
  * 'max_iter' - (Int):
	The maximum number of iterations to run.
  * tol - (Float):
	The tolerence to solve the algorithm to. Computes the tolerance by measuring
	the absolute value of the difference between the computed eigenvalues.
  * 'update_user' - (Int):
	Specifies how frequently output messages should be printed to the user.
	Default is -1 which means no output messages are printed, else if the value
	is k, then the kth iterate will print out an update.
  * 'no_matching' - (Bool):
	Specifies whether or not to run the matching and scoring portions of the
	algorithm. Useful if only the iterates or profiling is desired.

  Output
  ------
    * 'best_x'- (Array{Float,2})
      Returns the components to the iteration which matched the most triangles.
      Reshapes the iterate x into a matrix.
	* 'best_motif_count' - (Int)
	  The maximum number of triangles matched.
	* 'best_matching' - (Dict{Int,Int})
	  The matching computed between the two graphs, maps the vertices of A to the
	  vertices of B.
	* 'experiment_profile' - (Dict{String,Union{Array{F,1},Array{Array{F,1}}}}):
	  The experiment profile computed, keys for the experiment data collected are
	  as follows.
	  +'contraction_timings' - Time taken to compute each contraction.
	  +'matched_tris' - The number of triangles matched by each iterate.
	  +'gaped triangles' - The number of unmatched triangles, equal to the
	                       maximum - matched triangles.
	  +'sing_vals' - The singular values of each iterate X_k.
	  +'ranks' - The ranks of each iterate X_k.
	  +'matching_timings' - Time taken to solve the matchings.
	  +'scoring_timings' - Time taken to score each matching.

------------------------------------------------------------------------------"""
function TAME_profiled(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor, β::F, max_iter::Int,
	                   tol::F,α::F;update_user::Int=-1,W::Array{F,2} = ones(A.n,B.n),
					   no_matching::Bool=false,kwargs...) where {F <:AbstractFloat}

    dimension = minimum((A.n,B.n))

    experiment_profile = Dict{String,Union{Array{F,1},Array{Array{F,1}}}}(
		"contraction_timings"=>Array{F,1}(undef,0),
		"matching_timings"=>Array{F,1}(undef,0),
		"scoring_timings"=>Array{F,1}(undef,0),
		"matched triangles"=>Array{F,1}(undef,0),
		"gaped triangles"=>Array{F,1}(undef,0),
		"sing_vals"=>Array{Array{F,1},1}(undef,0),
		"ranks"=>Array{F,1}(undef,0)
	)

	A_Ti, B_Ti = setup_tame_data(A,B)

	best_x = Array{Float64,2}(undef,A.n,B.n)
    best_motif_count::Int = -1
	best_index = -1
	best_matching = Dict{Int,Int}()

    x0 = reshape(W,A.n*B.n)
	x0 ./=norm(x0)
    x_k = copy(x0)

    i = 1
    lambda = Inf

    while true

		x_k_1,t = @timed impTTVnodesym(A.n, B.n, x_k, A_Ti, B_Ti)
		push!(experiment_profile["contraction_timings"],t)

        new_lambda = dot(x_k_1,x_k)

        if β != 0.0
            x_k_1 .+= β * x_k
        end

        if α != 1.0
            x_k_1 = α * x_k_1 + (1 - α) * x0
        end

        x_k_1 ./= norm(x_k_1)

		S = svdvals(reshape(x_k_1,(A.n,B.n)))

		rank = 0.0
		for i in 1:length(S)
			if S[i] > S[1]*eps(Float64)*dimension
				rank = rank + 1
			end
		end

		push!(experiment_profile["sing_vals"],S)
		push!(experiment_profile["ranks"],rank)

		if !no_matching

			triangles, missed_motifs, matching, matching_time, scoring_time = TAME_score(A,B,reshape(x_k_1,A.n,B.n);return_timings=returnTimings(),kwargs...)

			push!(experiment_profile["matching_timings"],matching_time)
			push!(experiment_profile["scoring_timings"],scoring_time)
			push!(experiment_profile["matched triangles"],float(triangles))
			push!(experiment_profile["gaped triangles"],float(missed_motifs))

			if update_user != -1 && i % update_user == 0
				println("finished iterate $(i):tris:$(triangles) -- gaped_t:$(missed_motifs)")
			end

			if triangles > best_motif_count
				best_matching = matching
				best_x = copy(x_k_1)
				best_motif_count = triangles
				best_iterate = i
			end
		end

		if update_user != -1 && i % update_user == 0
			println("λ_$i: $(new_lambda)")
		end

        if abs(new_lambda - lambda) < tol || i >= max_iter
   			return reshape(best_x,A.n,B.n), best_motif_count, best_matching, experiment_profile
        else
            x_k = copy(x_k_1)
            lambda = new_lambda
            i += 1
        end

    end

end

# currently using types on A and B to alert user against passing Vector{SymTensorUnweighted{M}} in.
function TAME_profiled(A::SymTensorUnweighted{M}, B::SymTensorUnweighted{M}, β::F, max_iter::Int,
	                   tol::F,α::F;update_user::Int=-1,W::Array{F,2} = ones(A.n,B.n),
					   no_matching::Bool=false,kwargs...) where {F <:AbstractFloat, M <: Motif}

    dimension = minimum((A.n,B.n))

    experiment_profile = Dict{String,Union{Array{F,1},Array{Array{F,1}}}}(
		"contraction_timings"=>Array{F,1}(undef,0),
		"matching_timings"=>Array{F,1}(undef,0),
		"scoring_timings"=>Array{F,1}(undef,0),
		"svd_timings"=>Array{F,1}(undef,0),
		"matched motifs"=>Array{F,1}(undef,0),
		"missed motifs"=>Array{F,1}(undef,0),
		"sing_vals"=>Array{Array{F,1},1}(undef,0),
		"ranks"=>Array{F,1}(undef,0)
	)

	A_Mi, B_Mi = setup_tame_data(A,B)

	best_X = Array{Float64,2}(undef,A.n,B.n)
    best_motif_count::Int = -1
	best_index = -1
	best_matching = Dict{Int,Int}()

    X_k = copy(W)
	X_k ./=norm(X_k)

	if α != 1.0
		X_0_with_alpha = copy(X_k)*(1-α)
	end

    i = 1
    lambda = Inf

    while true

		X_k_1,t = @timed impTTVnodesym(A.n, B.n, size(A.indices,1), X_k, A_Mi, B_Mi)
		push!(experiment_profile["contraction_timings"],t)

        new_lambda = dot(X_k_1,X_k)

        if β != 0.0
            X_k_1 .+= β * X_k
        end

        if α != 1.0
            X_k_1 .*=  α
			X_k_1 += X_0_with_alpha
        end

        X_k_1 ./= norm(X_k_1)

		S,t = @timed svdvals(X_k_1)
		push!(experiment_profile["svd_timings"],t)

		rank = 0.0
		for i in 1:length(S)
			if S[i] > S[1]*eps(Float64)*dimension
				rank = rank + 1
			end
		end

		push!(experiment_profile["sing_vals"],S)
		push!(experiment_profile["ranks"],rank)

		if !no_matching

			matched_motifs, missed_motifs, matching, matching_time, scoring_time = TAME_score(A,B,X_k_1;return_timings=returnTimings(),kwargs...)

			push!(experiment_profile["matching_timings"],matching_time)
			push!(experiment_profile["scoring_timings"],scoring_time)
			push!(experiment_profile["matched motifs"],float(matched_motifs))
			push!(experiment_profile["missed motifs"],float(missed_motifs))

			if update_user != -1 && i % update_user == 0
				println("finished iterate $(i):tris:$(matched_motifs) -- gaped_t:$(missed_motifs)")
			end

			if matched_motifs > best_motif_count
				best_matching = matching
				best_X = copy(X_k_1)
				best_motif_count = matched_motifs
			end
		end

		if update_user != -1 && i % update_user == 0
			println("λ_$i: $(new_lambda)")
		end

        if abs(new_lambda - lambda) < tol || i >= max_iter
   			return best_X, best_motif_count, best_matching, experiment_profile
        else
            X_k = copy(X_k_1)
            lambda = new_lambda
            i += 1
        end

    end

end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
# zero!(tmpv)

@testset "Exploratory EV gradient" begin
    T = Float64
    t = tmpv
    p = prim
    σ = σv
    roy = 0.2
    geoid = 2
    itype = (geoid, roy,)
    θ1 = similar(θt)
    θ2 = similar(θt)
    EV1 = similar(evs.EV)
    EV2 = similar(evs.EV)
    fdEV = similar(evs.dEV)
    nsexp = ShaleDrillingModel._nSexp(wp)

    zero!(evs.EV)
    zero!(evs.dEV)
    zero!(evs.dEVσ)
    zero!(tmpv)
    zero!(EV1)
    zero!(EV2)
    zero!(fdEV )

    # run through θt
    for k = 1:length(θt)
        h = peturb(θt[k])
        θ1 .= θt
        θ2 .= θt
        θ1[k] -= h
        θ2[k] += h
        hh = θ2[k] - θ1[k]

        solve_vf_terminal!(evs, prim)
        solve_vf_infill!( evs, tmpv, prim, θ1, σv, false, itype)
        learningUpdate!(  evs, tmpv, prim,     σv, false)
        solve_vf_explore!(evs, tmpv, prim, θ1, σv, false, itype)
        fdEV[:,:,k,:] .= -evs.EV

        solve_vf_terminal!(evs, prim)
        solve_vf_infill!( evs, tmpv, prim, θ2, σv, false, itype)
        learningUpdate!(  evs, tmpv, prim,     σv, false)
        solve_vf_explore!(evs, tmpv, prim, θ2, σv, false, itype)

        fdEV[:,:,k,:] .+= evs.EV
        fdEV[:,:,k,:] ./= hh
    end

    # ----------------- analytic -----------------

    solve_vf_terminal!(evs, prim)
    solve_vf_infill!( evs, tmpv, prim, θt, σv, true, itype)
    learningUpdate!(  evs, tmpv, prim,     σv, true)
    solve_vf_explore!(evs, tmpv, prim, θt, σv, true, itype)

    # check infill + learning portion gradient
    @views maxv, idx = findmax(abs.(fdEV[:,:,:,nsexp+1:end].-evs.dEV[:,:,:,nsexp+1:end]))
    println("worst value is $maxv at $(CartesianIndices(fdEV)[idx]) for infill")

    # check exploration portion of gradient
    @views maxv, idx = findmax(abs.(fdEV[:,:,:,1:nsexp].-evs.dEV[:,:,:,1:nsexp]))
    println("worst value is $maxv at $(CartesianIndices(fdEV)[idx]) for exploratory")

    @test 0.0 < maxv < 0.1
    @test all(isfinite.(evs.dEV))
    @test all(isfinite.(fdEV))

    @test fdEV ≈ evs.dEV || maxv < 1e-5
    println("dEV/dθ looks ok! :)")
end

@testset "dEV/dσ" begin
    let T = Float64,
        EV1 = zeros(T, size(evs.EV)),
        EV2 = zeros(T, size(evs.EV)),
        t = tmpv,
        p = prim,
        nsexp1 = ShaleDrillingModel._nSexp(wp),
        σ = σv,
        roy = 0.2,
        geoid = 2,
        itype = (geoid, roy,),
        θ1 = similar(θt),
        θ2 = similar(θt),
        fdEVσ = zeros(T, size(evs.EV)),
        fdubv = zeros(T,size(tmpv.ubVfull)),
        fdEVσvw = @view(fdEVσ[:,:,1:nsexp1])

        println("testing dEV/dσ")

        # now do σ
        h = peturb(σv)
        σ1 = σv - h
        σ2 = σv + h
        hh = σ2 - σ1

        dograd = false

        zero!(evs)
        zero!(tmpv)
        solve_vf_all!(evs, tmpv, prim, θt, σ1, itype, dograd)
        # solve_vf_terminal!(evs, prim)
        # solve_vf_infill!( evs, tmpv, prim, θt, σ1, dograd, itype)
        # learningUpdate!(  evs, tmpv, prim,     σ1, dograd)
        # solve_vf_explore!(evs, tmpv, prim, θt, σ1, dograd, itype)

        fdEVσ .= -evs.EV

        zero!(evs)
        zero!(tmpv)
        solve_vf_all!(evs, tmpv, prim, θt, σ2, itype, dograd)
        # solve_vf_terminal!(evs, prim)
        # solve_vf_infill!( evs, tmpv, prim, θt, σ2, dograd, itype)
        # learningUpdate!(  evs, tmpv, prim,     σ2, dograd)
        # solve_vf_explore!(evs, tmpv, prim, θt, σ2, dograd, itype)

        fdEVσ .+= evs.EV
        fdEVσ ./= hh

        # ----------------- analytic -----------------

        dograd = true

        zero!(evs)
        zero!(tmpv)
        solve_vf_all!(evs, tmpv, prim, θt, σv, itype, dograd)
        # solve_vf_terminal!(evs, prim)
        # solve_vf_infill!(  evs, tmpv, prim, θt, σv, dograd, itype)
        # learningUpdate!(   evs, tmpv, prim,     σv, dograd)
        # solve_vf_explore!( evs, tmpv, prim, θt, σv, dograd, itype)

        @test size(fdEVσvw) == size(evs.dEVσ)
        @test all(isfinite.(evs.dEVσ))
        @test all(isfinite.(fdEVσvw))
        @show extrema(fdEVσvw)
        @show extrema(evs.dEVσ)

        println("Testing evs.dEVσ[:,:,2:end] .- fdEVσvw[:,:,2:end]")
        @show maximum(abs.((evs.dEVσ .- fdEVσvw)[:,:,2:end]))
        @test @views fdEVσvw[:,:,2:end] ≈ evs.dEVσ[:,:,2:end]

        maxv, idx = findmax(abs.(fdEVσvw.-evs.dEVσ))
        println("worst value is $maxv at $(CartesianIndices(fdEVσvw)[idx]) for dσ")
        @test 0.0 < maxv < 0.1
        @test @views fdEVσvw[:,:,1] ≈ evs.dEVσ[:,:,1]

        @test fdEVσvw ≈ evs.dEVσ # || maxv < 1.5e-6
        println("dEV/dσ looks ok! :)")
    end
end


text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
  x1'(t) = a1 * (x2(t) - x1(t)) - (ka * n * x1(t)) / (kc * ka + kc * x3(t) + ka * x1(t)),
  x2'(t) = a1 * (x1(t) - x2(t)),
  x3'(t) = b1 * (x4(t) - x3(t)) - (kc * n * x3(t)) / (kc * ka + kc * x3(t) + ka * x1(t)),
  x4'(t) = b2 * (x3(t) - x4(t)),
  y1(t) = x1(t)
)

res = identifiability_ode(ode, get_parameters(ode); p=0.99, p_mod=2^29 - 3, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
"""
    MicrophysicsParameters

Module containing 1-moment bulk microphysics parameters.
"""
module MicrophysicsParameters
using CLIMA.ParametersType

@exportparameter MP_n_0          8e6 * 2         "Marshall-Palmer distribution n_0 coeff [1/m4]"
@exportparameter C_drag          0.55            "drag coefficient for rain drops [-]"

@exportparameter τ_cond_evap     10              "condensation/evaporation timescale [s]"

@exportparameter q_liq_threshold 5e-4            "autoconversion threshold [-]  ∈(0.5, 1) * 1e-3 "
@exportparameter τ_acnv          1e3             "autoconversion timescale [s]  ∈(1e3, 1e4) "

@exportparameter E_col           0.8             "collision efficiency [-]"

@exportparameter a_vent          1.5             "ventilation factor coefficient [-]"
@exportparameter b_vent          0.53            "ventilation factor coefficient [-]"
@exportparameter K_therm         2.4e-2          "thermal conductivity of air [J/m/s/K] "
@exportparameter D_vapor         2.26e-5         "diffusivity of water vapor [m2/s]"
@exportparameter ν_air           1.6e-5          "kinematic viscosity of air [m2/s]"
@exportparameter N_Sc            ν_air/D_vapor   "Schmidt number [-]"
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
struct htest
    title::String
    d::OrderedDict
end

import Base.show

Base.show(io::IO, ::MIME"text/plain", h::htest) = show(io,h)

function show(io::IO, h::htest)
    print(io,"\n",h.title,"\n\n")

    for k in keys(h.d)
        if k == "note"
            continue
        end
        print(io,lpad(k,13)," = ",h.d[k],"\n")
    end

    if haskey(h.d,"note") && h.d["note"] != ""
        print(io,"\nNOTE: ", h.d["note"],"\n")
    end
end

function check_args(; kwargs...)
    for (key,val) in kwargs
        if key in (:k, :n, :n1, :n2, :w, :f, :f2, :r) && val == 0
            error("`",string(key),"` must be specified and greater than zero")
        elseif key in (:u, :v, :df) && val < 1
            error("Degree of freedom `",key,"` must be at least 1")
        elseif key == :w && val <= 0
            error("`w` must be positive")
        elseif key in (:alpha,:power) && (val <= 0.0 || val >= 1.0)
            error("`",key,"` must be a number in [0,1]")
        elseif key == :N && val < 1
            error("Number of observations `N` must be at least 1")
        end
    end
end

max_stars_repo_name: mwsohn/pwr

max_issues_repo_name: mwsohn/pwr

max_forks_repo_name: mwsohn/pwr


text:
using PlotuDiff
using PlotuDiff: Insertion, Deletion, Equality, Myers
using PlotuDiff: stringdiff_without_viewer as sdiff
using Test

@testset "Myers algorithm" begin
	a = "I LIKE PIZZA"
	b = "YOU LIKE MEZZE"

	@testset "Two lines" begin

		diffs = sdiff(a, b; alg=Myers())

		@test diffs == [
			Deletion('I'),
			[Insertion(c) for c in "YOU"]...,
			[Equality(c) for c in " LIKE "]...,
			Deletion('P'),
			Deletion('I'),
			Insertion('M'),
			Insertion('E'),
			Equality('Z'),
			Equality('Z'),
			Deletion('A'),
			Insertion('E'),
		]
	end

	@testset "Merge diffs" begin
		diffs = sdiff(a, b; alg=Myers()) |> PlotuDiff.merge_diffs

		@test diffs == [
			Deletion('I'),
			Insertion("YOU"),
			Equality(" LIKE "),
			Deletion("PI"),
			Insertion("ME"),
			Equality("ZZ"),
			Deletion('A'),
			Insertion('E'),
		]
	end
end