text:
"""
    AbstractEdgeWeight

A supertype for the weights that are put on the edges of the clique graph.
"""
abstract type AbstractEdgeWeight end

"""
    AbstractComplexityWeight

A supertype for weights that are based on the complexity of clique-related operations in the
solver algorithm.

All subtypes have to define the function:
`compute_complexity_savings(n_1::Int, n_2::Int, n_m::Int, edge_weight::Subtype) --> Float64`

with
  - n_1: `|C_1|`
  - n_2: `|C_2|`
  - n_3: `|C_1 ∪ C_2|`
"""
abstract type AbstractComplexityWeight <: AbstractEdgeWeight end

"""
    ComplexityWeight

An edge weight that approximates the computational savings from merging two cliques
`Ci` and `Cj` by `e(Ci, Cj) = |Ci|^3 + |Cj|^3 - |Ci ∪ Cj|^3`.
"""
struct ComplexityWeight <: AbstractComplexityWeight end

"""
    NoMerge <: AbstractMergeStrategy

A strategy that does not merge cliques.
"""
struct NoMerge <: AbstractTreeBasedMerge end

"""
    CliqueGraphMerge(edge_weight::AbstractEdgeWeight = ComplexityWeight()) <: AbstractGraphBasedMerge

The (default) merge strategy based on the *reduced* clique graph ``\\mathcal{G}(\\mathcal{B}, \\xi)``, for a set of cliques ``\\mathcal{B} = \\{ \\mathcal{C}_1, \\dots, \\mathcal{C}_p\\}`` where the edge set ``\\xi`` is obtained by taking the edges of the union of clique trees.

Moreover, given an edge weighting function ``e(\\mathcal{C}_i,\\mathcal{C}_j) = w_{ij}``, we compute a weight for each edge that quantifies the computational savings of merging the two cliques.
After the initial weights are computed, we merge cliques in a loop:

**while** clique graph contains positive weights:
- select two permissible cliques with the highest weight ``w_{ij}``
- merge cliques ``\\rightarrow`` update clique graph
- recompute weights for updated clique graph

Custom edge weighting functions can be used by defining your own `CustomEdgeWeight <: AbstractEdgeWeight` and a corresponding `edge_metric` method. By default, the `ComplexityWeight <: AbstractEdgeWeight` is used which computes the weight based
on the cardinalities of the cliques: ``e(\\mathcal{C}_i,\\mathcal{C}_j)  = |\\mathcal{C}_i|^3 + |\\mathcal{C}_j|^3 - |\\mathcal{C}_i \\cup \\mathcal{C}_j|^3``.

See also: *Garstka, Cannon, Goulart - A clique graph based merging strategy for decomposable SDPs (2019)*
"""
mutable struct CliqueGraphMerge <: AbstractGraphBasedMerge
  stop::Bool                                  # a flag to indicate that merging should be stopped
  edges::SparseMatrixCSC{Float64, Int}      # the edges and weights of the reduced clique graph
  p::Array{Int, 1}                          # as a workspace variable to store the sorting of weights
  adjacency_table::Dict{Int, Set{Int}}    # a double structure of edges, to allow fast lookup of neighbors
  edge_weight::AbstractEdgeWeight             # used to dispatch onto the correct scoring function
  clique_tree_recomputed::Bool                # a flag to indicate whether a final clique tree has been recomputed from the clique graph
  function CliqueGraphMerge(; edge_weight = ComplexityWeight())
    new(false, spzeros(0, 0),  Int[], Dict{Int, Set{Int}}(), edge_weight, false)
  end
end

"""
    ParentChildMerge(t_fill = 8, t_size = 8) <: AbstractTreeBasedMerge

The merge strategy suggested in *Sun and Andersen - Decomposition in conic optimization with partially separable structure (2014)*.
The initial clique tree is traversed in topological order and a clique ``\\mathcal{C}_\\ell``  is greedily merged to its parent clique
``\\mathcal{C}_{par(\\ell)}`` if at least one of the two conditions are met

  - ``(| \\mathcal{C}_{par(\\ell)}| -| \\eta_\\ell|) (|\\mathcal{C}_\\ell| - |\\eta_\\ell|) \\leq t_{\\text{fill}}`` (fill-in condition)
  - ``\\max \\left\\{ |\\nu_{\\ell}|,  |\\nu_{par(\\ell)}|  \\right\\} \\leq t_{\\text{size}}`` (supernode size condition)
"""
mutable struct ParentChildMerge <: AbstractTreeBasedMerge
  stop::Bool
  clique_ind::Int
  t_fill::Int
  t_size::Int

  function ParentChildMerge(; t_fill = 8, t_size = 8)
    new(false, 2, t_fill, t_size)
  end
end

"Free unneccessary memory after merging completed."
function free_clique_graph!(strategy::AbstractGraphBasedMerge)
  strategy.edges = spzeros(0, 0)
  strategy.p = Int[]
  strategy.adjacency_table = Dict{Int, Set{Int}}()
  return nothing
end

# Main clique merging routine:
# 1. initialise!() - tree and strategy
# 2. traverse() - find the next merge candidates
# 3. evalute() - evaluate wether to merge them or not
# 4. log() - log the decision
# 5. update_strategy!() - update the tree/graph and strategy
function _merge_cliques!(t::SuperNodeTree, strategy::AbstractMergeStrategy)
  initialise!(t, strategy)

  while !strategy.stop
    # find merge candidates

    cand = traverse(t, strategy)
    # evaluate wether to merge the candidates
    do_merge = evaluate(t, strategy, cand)
    if do_merge
      merge_two_cliques!(t, cand, strategy)
    end
    log_merge!(t, do_merge, cand)

    # update strategy information after the merge
    update_strategy!(strategy, t, cand, do_merge)
    t.num == 1 && break
    strategy.stop && break

  end
  return nothing
end

"Given a supernode elimination tree `t` and a merge strategy `strategy`, merge the cliques in the tree."
merge_cliques!(t::SuperNodeTree) = merge_cliques!(t, t.strategy)
merge_cliques!(t::SuperNodeTree, strategy::NoMerge) = nothing

function merge_cliques!(t::SuperNodeTree, strategy::AbstractTreeBasedMerge)
  # call main merge routine
  _merge_cliques!(t, strategy)
  # do some merge strategy specific post-processing of the tree
  Nc = t.num
  # the merging creates empty supernodes and seperators, recalculate a post order for the supernodes
  snd_post = post_order(t.snd_par, t.snd_child, Nc)
  t.snd_post = snd_post

  return nothing
end

function merge_cliques!(t::SuperNodeTree, strategy::AbstractGraphBasedMerge)
  # call main merge routine
  _merge_cliques!(t, strategy)
  # since for now we have a graph, not a tree, a post ordering or a parent structure does not make sense. Therefore just number
  # the non-empty supernodes in t.snd
  t.snd_post = findall(x -> !isempty(x), t.snd)
  t.snd_par = -ones(Int, length(t.snd))

  # recomute a clique tree from the clique graph
  t.num > 1 && clique_tree_from_graph!(t)

  #turn the unordered sets for t.snd and t.sep into orderable arrays
  t.snd = sort.(collect.(t.snd))
  t.sep = sort.(collect.(t.sep))

  # free up memory
  free_clique_graph!(t.strategy)
  return nothing
end

# calculate block sizes (notice: the block sizes are stored in post order)
function calculate_block_dimensions!(t::SuperNodeTree)
  Nc = t.num
  t.nBlk = zeros(Nc)
  for iii = 1:Nc
    c = t.snd_post[iii]
    t.nBlk[iii] = length(t.sep[c]) + length(t.snd[c])
  end
end

" Merge two cliques `cand` of the tree `t` that are in a parent - child relationship."
function merge_child!(t::SuperNodeTree, cand::Array{Int, 1})
  # determine which clique is the parent
  p, ch = determine_parent(t, cand[1], cand[2])

  # merge child's vertex sets into parent's vertex set
  union!(t.snd[p], t.snd[ch])
  empty!(t.snd[ch])
  empty!(t.sep[ch])

  # update parent structure
  for grandch in t.snd_child[ch]
    t.snd_par[grandch] = p
  end
  t.snd_par[ch] = -1 #-1 instead of NaN, effectively remove that entry from the parent list

  # update children structure
  delete!(t.snd_child[p], ch)
  union!(t.snd_child[p], t.snd_child[ch])
  empty!(t.snd_child[ch])

  # decrement number of cliques in tree
  t.num -= 1
  return nothing
end

"Given the clique graph `t` merge the two cliques with indices in `cand`."
function merge_two_cliques!(t::SuperNodeTree, cand::Array{Int, 1}, strategy::AbstractGraphBasedMerge)
  c_1 = cand[1]
  c_2 = cand[2]
  # merge clique c_2 into c_1
  union!(t.snd[c_1], t.snd[c_2])
  empty!(t.snd[c_2])

  # decrement number of cliques in tree
  t.num -= 1

  return nothing
end

"Given the clique tree `t` merge the two cliques with indices in `cand` as parent-child."
merge_two_cliques!(t::SuperNodeTree, cand::Array{Int, 1}, strategy::ParentChildMerge) = merge_child!(t, cand)

# compute the edges and weights of the initial reduced clique graph
function initialise!(t::SuperNodeTree, strategy::CliqueGraphMerge)
  # save("rs1555_before_merging.jld2", "snd", t.snd, "sep", t.sep)
  # compute the edges and intersections of cliques in the reduced clique graph
  rows, cols = compute_reduced_clique_graph!(t.sep, t.snd)

  weights = compute_weights!(rows, cols, t.snd, strategy.edge_weight)

  strategy.edges = sparse(rows, cols, weights, t.num, t.num)
  strategy.p = zeros(Int, length(strategy.edges.nzval))
  strategy.adjacency_table = compute_adjacency_table(strategy.edges, t.num)
  return nothing
end

function initialise!(t::SuperNodeTree, strategy::ParentChildMerge)
  # start with node that has second highest order
  strategy.clique_ind = length(t.snd) - 1
end

"Compute the edge weight between all cliques specified by the edges (rows, cols)."
function compute_weights!(rows::Array{Int, 1}, cols::Array{Int, 1}, snd::Array{Set{Int}, 1}, edge_weight::AbstractComplexityWeight)
  weights = zeros(Float64, length(rows))
  for k = 1:length(rows)
    c_1 = snd[rows[k]]
    c_2 = snd[cols[k]]
    weights[k] = edge_metric(c_1, c_2, edge_weight)
  end
  return weights
end

"Find the next two cliques in the clique graph `t` to merge."
function traverse(t::SuperNodeTree, strategy::CliqueGraphMerge)

   p = strategy.p
   # find edge with highest weight, if permissible return cliques
   edge = max_elem(strategy.edges)
   ispermissible(edge, strategy.adjacency_table, t.snd) && return [edge[1]; edge[2]]
   # else: sort the weights in edges.nzval to find the permutation p
   sortperm!(view(p, 1:length(strategy.edges.nzval)), strategy.edges.nzval, alg = QuickSort, rev = true)

   # try edges with decreasing weight and check if the edge is permissible
  for k = 2:length(strategy.edges.nzval)
    edge = edge_from_index(strategy.edges, p[k])
    if ispermissible(edge, strategy.adjacency_table, t.snd)
      return [edge[1]; edge[2]]
    end
  end

end

" Traverse tree `t` in descending topological order and return parent and clique (root has highest order)."
function traverse(t::SuperNodeTree, strategy::ParentChildMerge)
  c = t.snd_post[strategy.clique_ind]
  return [t.snd_par[c]; c]
end

"Decide whether to merge the two cliques with clique indices `cand`."
function evaluate(t::SuperNodeTree, strategy::ParentChildMerge, cand::Array{Int, 1})
  strategy.stop && return false
  par = cand[1]
  c = cand[2]
  dim_par_snd, dim_par_sep = COSMO.clique_dim(t, par)
  dim_clique_snd, dim_clique_sep = COSMO.clique_dim(t, c)
  return fill_in(dim_clique_snd, dim_clique_sep, dim_par_snd, dim_par_sep) <= strategy.t_fill || max_snd_size(dim_clique_snd, dim_par_snd) <= strategy.t_size
end

function evaluate(t::SuperNodeTree, strategy::CliqueGraphMerge, cand::Array{Int, 1})
  do_merge = (strategy.edges[cand[1], cand[2]] >= 0)

  if !do_merge
    strategy.stop = true
  end
  return do_merge
end

" After a merge attempt, update the strategy information."
function update_strategy!(strategy::ParentChildMerge, t::SuperNodeTree, cand::Array{Int, 1}, do_merge::Bool)
  # try to merge last node of order 1, then stop
  if strategy.clique_ind == 1
    strategy.stop = true
    # otherwise decrement node index
  else
    strategy.clique_ind -= 1
  end
end

"After a merge happened, update the reduced clique graph."
function update_strategy!(strategy::CliqueGraphMerge, t::SuperNodeTree, cand::Array{Int, 1}, do_merge::Bool)

  # After a merge operation update the information of the strategy
  if do_merge

    c_1_ind = cand[1]
    c_removed = cand[2]
    edges = strategy.edges
    n = size(edges, 2)
    adjacency_table = strategy.adjacency_table

    c_1 = t.snd[c_1_ind]
    neighbors = adjacency_table[c_1_ind]
    # neighbors exclusive to the removed clique (and not c1)
    new_neighbors = setdiff(adjacency_table[c_removed], neighbors, c_1_ind)

    # recalculate edge values of all of c_1's neighbors
    for n_ind in neighbors
        if n_ind != c_removed
          neighbor = t.snd[n_ind]
          edges[max(c_1_ind, n_ind), min(c_1_ind, n_ind)] = COSMO.edge_metric(c_1, neighbor, strategy.edge_weight)
        end
    end

    # point edges exclusive to removed clique to "surviving" clique 1
    for n_ind in new_neighbors
        neighbor = t.snd[n_ind]
        edges[max(c_1_ind, n_ind), min(c_1_ind, n_ind)]  = COSMO.edge_metric(c_1, neighbor, strategy.edge_weight)
    end

    # overwrite the weight to any "deleted" edges that still contain a link to c_removed
    strategy.edges[c_removed+1:n, c_removed] .= 0
    strategy.edges[c_removed, 1:c_removed] .= 0
    dropzeros!(edges)

    # update adjacency table in a similar manner
    union!(adjacency_table[c_1_ind], new_neighbors)
    for new_neighbor in new_neighbors
      push!(adjacency_table[new_neighbor], c_1_ind)
    end
    delete!(adjacency_table, c_removed)
    for (key, set) in adjacency_table
      delete!(set, c_removed)
    end
  end

  return nothing
end

"""
    isdisjoint(c_a, c_b, is_sorted = false)

Checks whether two sets `c_a`, `c_b` have no common elements. Assumes by default that the sets are sorted.
"""
function isdisjoint(c_a::AbstractVector, c_b::AbstractVector; is_sorted = false)
  if !is_sorted
    sort!(c_a)
    sort!(c_b)
  end
  m = length(c_a)
  n = length(c_b)

  i = 1
  j = 1
  while i <= m && j <= n
    if c_a[i] < c_b[j]
      i += 1
    elseif c_b[j] < c_a[i]
      j += 1
    else
      return false
    end
  end
  return true
end

"""
    edge_metric(c_a::AbstractVector, c_b::AbstractVector, edge_weight::AbstractComplexityWeight)

Given two cliques `c_a` and `c_b` return a value for their edge weight.
"""
function edge_metric(c_a::T, c_b::T, edge_weight::AbstractComplexityWeight) where {T <: Union{AbstractVector, AbstractSet}}
  n_1 = length(c_a)
  n_2 = length(c_b)

  # merged block size
  n_m = union_dim(c_a, c_b)
  return compute_complexity_savings(n_1, n_2, n_m, edge_weight)
end

# Assuming the complexity of the projection is roughly O(n^3), how many operations are saved by projection the merged cliques
# instead of the individual cliques
compute_complexity_savings(n_1::Int, n_2::Int, n_m::Int, edge_weight::ComplexityWeight) = n_1^3 + n_2^3 - n_m^3

"""
    max_elem(A::SparseMatrixCSC)

Find the matrix indices (i, j) of the first maximum element among the elements stored in A.nzval
"""
function max_elem(A::SparseMatrixCSC)
  length(A.nzval) == 0 && throw(DomainError("Sparse matrix A doesn't contain any entries"))
  n = size(A, 2)

  ~, ind = findmax(A.nzval)
  row = A.rowval[ind]

  col = 0
  for c = 1:n
    col_indices = A.colptr[c]:A.colptr[c+1]-1
    if in(ind, col_indices)
      col = c
      break;
    end
  end
  return (row, col)
end

function edge_from_index(A::SparseMatrixCSC{Float64, Int}, ind::Int)
  # find the edge for that value
  row = A.rowval[ind]
  n = size(A, 2)
  col = 0
  for c = 1:n
    col_indices = A.colptr[c]:A.colptr[c+1]-1
    if in(ind, col_indices)
      col = c
      break;
    end
  end
  return (row, col)
end

" Given two cliques `c1` and `c2` in the tree `t`, return the parent clique first."
function determine_parent(t::SuperNodeTree, c1::Int, c2::Int)
  if in(c2, t.snd_child[c1])
    return c1, c2
  else
    return c2, c1
  end
end

"""
find_neighbors(edges::SparseMatrixCSC, c::Int)

Find all the cliques connected to `c` which are given by the nonzeros in `(c, 1:c-1)` and `(c+1:n, c)`.
"""
function find_neighbors(edges::SparseMatrixCSC, c::Int)
  neighbors = zeros(Int, 0)
  m, n = size(edges)
  # find all nonzero columns in row c up to column c
  if c > 1
   neighbors = vcat(neighbors, findall(x -> x!= 0, edges[c, 1:c-1]))
 end
 # find all nonzero entries in column c below c
 if c < n
  rows = edges.rowval[edges.colptr[c]:edges.colptr[c+1]-1]
  if edges.colptr[c] <= edges.colptr[c+1] - 1
    neighbors = vcat(neighbors, rows)
  end
end
return neighbors
end

function clique_intersections!(E::SparseMatrixCSC, snd::Array{Set{Int}, 1})
  # iterate over the nonzeros of the connectivity matrix E which represents the clique graph and replace the value by
  # |C_i ∩ C_j|
  rows = rowvals(E)
  for col in 1:size(E, 2)
    for j in nzrange(E, col)
      row = rows[j]
      E[row, col] = intersect_dim(snd[row], snd[col])
    end
  end
  return nothing
end

"""
    kruskal!(E::SparseMatrixCSC, num_cliques::Int)

Kruskal's algorithm to find a maximum weight spanning tree from the clique intersection graph.

 `E[i,j]` holds the cardinalities of the intersection between two cliques (i, j). Changes the entries in the connectivity matrix `E` to a negative
 value if an edge between two cliques is included in the max spanning tree.

 This is a modified version of https://github.com/JuliaGraphs/LightGraphs.jl/blob/master/src/spanningtrees/kruskal.jl
 """
function kruskal!(E::SparseMatrixCSC, num_cliques::Int)
  num_initial_cliques = size(E, 2)
  connected_c = DataStructures.IntDisjointSets(num_initial_cliques)

  I, J, V = findnz(E)
  # sort the weights and edges from maximum to minimum value
  p = sortperm(V, rev = true)
  I = I[p]
  J = J[p]
  V = V[p]
  num_edges_found = 0
  # iterate through edges (I -- J) with decreasing weight
  for k = 1:length(V)
    row = I[k]
    col = J[k]
    if !in_same_set(connected_c, row, col)
      union!(connected_c, row, col)
      # we indicate an edge in the MST with a positive value in E (all other values are >= 0)
      E[row, col] = -1.0
      num_edges_found += 1
      # break when all cliques are connected in one tree
      num_edges_found >= num_cliques - 1 && break
    end
  end
  return nothing
end

function assign_children!(snd_par::Array{Int, 1}, snd_child::Array{Set{Int}, 1}, c::Int, edges::SparseMatrixCSC)
  # determine neighbors
  neighbors = find_neighbors(edges, c)
  for n in neighbors
    # conditions that there is a edge in the MST and that n is not the parent of c
    if edges[max(c, n), min(c, n)] == -1.0 && snd_par[c] != n
      snd_par[n] = c
      push!(snd_child[c], n)
      assign_children!(snd_par, snd_child, n, edges)
    end
  end
  return nothing
end

" Given the maximum weight spanning tree represented by `E`, determine a parent structure `snd_par` for the clique tree."
function determine_parent_cliques!(snd_par::Array{Int, 1}, snd_child::Array{Set{Int}, 1}, cliques::Array{Set{Int}, 1}, post::Array{Int, 1}, E::SparseMatrixCSC)
  # vertex with highest order
  v = post[end]
  c = 0
  # Find clique that contains that vertex
  for (k, clique) in enumerate(cliques)
    if v ∈ clique
      # set that clique to the root
      snd_par[k] = 0
      c = k
      break
    end
  end

  # recursively assign children to cliques along the MST defined by E
  assign_children!(snd_par, snd_child, c, E)
  return nothing
end

" Traverse the clique tree in descending topological order and split the clique sets into supernodes and separators."
function split_cliques!(snd::Array{Set{Int}, 1}, sep::Array{Set{Int}, 1}, snd_par::Array{Int,1}, snd_post::Array{Int}, num_cliques::Int)

  # travese in topological decending order through the clique tree and split the clique in snd and sep
  for j = 1:1:(num_cliques - 1)
    c_ind = snd_post[j]
    par_ind = snd_par[c_ind]

    # find intersection of clique with parent
    # FIXME: Can be made a lot faster by using non-standard functions
    sep[c_ind] = intersect(snd[c_ind], snd[par_ind])
    snd[c_ind] = filter!(x -> x ∉ sep[c_ind], snd[c_ind])
  end
  return nothing
end

"""
    clique_tree_from_graph!(clique_graph::SuperNodeTree)

Given the cliques and edges of a clique graph, this function computes a valid clique tree.

This is necessary to perform the psd completion step of the dual variable after solving the problem.
"""
function clique_tree_from_graph!(t::SuperNodeTree)
  # a clique tree is a maximum weight spanning tree of the clique graph where the edge weight is the cardinality of the intersection between two cliques
  # compute intersection value for each edge in clique graph
  # clique_intersections!(t.strategy.edges, t.snd)
  clique_intersections!(t.strategy.edges, t.snd)

  # find a maximum weight spanning tree of the clique graph using Kruskal's algorithm
  kruskal!(t.strategy.edges, t.num)

  # determine the root clique of the clique tree (it can be any clique, but we AbstractComplexityWeightose the clique that contains the vertex with the highest order)
  determine_parent_cliques!(t.snd_par, t.snd_child, t.snd, t.post, t.strategy.edges)

  # recompute a postorder for the supernodes
  t.snd_post = post_order(t.snd_par, t.snd_child, t.num)

  t.sep = [Set{Int}() for i = 1:length(t.snd)]
  # split clique sets back into seperators and supernodes
  split_cliques!(t.snd, t.sep, t.snd_par, t.snd_post, t.num)
  t.strategy.clique_tree_recomputed = true
  return nothing
end

"Return the number of elements in s ∩ s2."
function intersect_dim(s::Set, s2::Set)
  if length(s) < length(s2)
        sa = s
        sb = s2
    else
        sa = s2
        sb = s
    end
    dim = 0
    for e in sa
        e in sb && (dim += 1)
    end
    return dim
end

" Find the size of the set `A ∪ B` under the assumption that `A` and `B` only have unique elements."
function union_dim(A::T, B::T)  where {T <: Union{AbstractVector, AbstractSet}}
  dim = length(A)
  for elem in B
    !in(elem, A) && (dim += 1)
  end
  return dim
end

"Compute the amount of fill-in created by merging two cliques with the respective supernode and separator dimensions."
function fill_in(dim_clique_snd::Int, dim_clique_sep::Int, dim_par_snd::Int, dim_par_sep::Int)
  dim_par = dim_par_snd + dim_par_sep
  dim_clique = dim_clique_snd + dim_clique_sep
  return ((dim_par - dim_clique_sep) * (dim_clique - dim_clique_sep))::Int
end

max_snd_size(dim_clique_snd::Int, dim_par_snd::Int) = max(dim_clique_snd, dim_par_snd)
clique_dim(t, c_ind) = length(t.snd[c_ind]), length(t.sep[c_ind])

"Store information about the merge of the two merge candidates `cand`."
function log_merge!(t::SuperNodeTree, do_merge::Bool, cand::Array{Int, 1})
  t.merge_log.clique_pairs = vcat(t.merge_log.clique_pairs, cand')
  push!(t.merge_log.decisions, do_merge)
  do_merge && (t.merge_log.num += 1)
  return nothing
end

"""
    print_merge_logs(ws; print_max = 10)

Print the merge logs for each decomposed cone in the problem.
"""
function print_merge_logs(ws; print_max = 10)
  sp_arr = ws.ci.sp_arr
  # print the merge logs for each sparsity pattern
  println(">>> Merge Logs:")
  for (iii, sp) in enumerate(sp_arr)
    m_log = sp.sntree.merge_log
    println("Sparsity Pattern Nr. $(iii), Graph Size: $(length(sp.sntree.par))")
    println("\t Num merges: $(m_log.num)\n\t Num decisions: $(length(m_log.decisions))\n\t Pairs:")

    # only print a certain number
    print_max = min(print_max, length(m_log.decisions))
    for jjj = 1:print_max
      p = m_log.clique_pairs[jjj, :]
      println("\t \tC_A: $(p[1]), C_B: $(p[2]) ($(m_log.decisions[jjj]))")
    end

  end
end

"Return an array of `MergeLog`s corresponding to the decomposed cones."
get_merge_logs(ws) = map(x -> x.sntree.merge_log, ws.ci.sp_arr)


text:
using Dates
using Intervals
using Test
using YeSQL

queries = load_path(joinpath(@__DIR__, "sql", "queries.sql"))

@testset "Test query loading" begin
    @test !(:get_unique_instruments in keys(queries))
    # load all SQL files
    qpath = load_path(joinpath(@__DIR__, "sql"))
    @test :create_test_types in keys(qpath)
end

@testset "Test interval string representation" begin
    d1 = Date(2020,1,1)
    d2 = Date(2021,1,1)
    @test YeSQL._libpq_literal(d1..d2) == "[2020-01-01,2021-01-01]"
    @test YeSQL._libpq_literal([d1..d2, d1..nothing]) == ["[2020-01-01,2021-01-01]", "[2020-01-01,)"]
    @test YeSQL._libpq_literal([3..4, Interval{Closed,Open}(3,10)]) == ["[3,4]", "[3,10)"]
    @test YeSQL._libpq_literal(d1..nothing) == "[2020-01-01,)"
    @test YeSQL._libpq_literal(Interval{Date}(nothing, nothing)) == "(,)"
    @test YeSQL._libpq_literal(Interval{Date,Unbounded,Open}(nothing,d2)) == "(,2021-01-01)"
    @test YeSQL._libpq_literal(Interval{Date,Unbounded,Closed}(nothing,d2)) == "(,2021-01-01]"
end

@testset "Test parsing of query names" begin
    names = Dict(
        "-- name: normal_select" => (YeSQL.Select(), :normal_select),
        "-- name: select_value\$" => (YeSQL.SelectValue(), :select_value),
        "-- name: select_one^" => (YeSQL.SelectOne(), :select_one),
        "-- name: script#" => (YeSQL.Script(), :script),
        "-- name: insert!" => (YeSQL.InsertUpdateDelete(), :insert),
        "-- name: insert_many*!" => (YeSQL.InsertUpdateDeleteMany(), :insert_many),
        "-- name: insert_returning^!" => (YeSQL.InsertReturning(), :insert_returning),
        "-- name: insert_returning_value\$!" => (YeSQL.InsertReturningValue(), :insert_returning_value),
    )

    for (name_str, (target_op, target_name)) in names
        nm = match(YeSQL.NAME_REGEX, name_str)
        operation, query_name = YeSQL._interpret(nm[:query_name])
        @test operation == target_op
        @test query_name == target_name
    end
end

@testset "Test query extraction from string" begin
    s1 = "-- name: string_query1\n select count(*) from my_table where mycolumn = :myvalue"
    s2 = "-- name: string_query2\n select count(*) from my_table where mycolumn = :myvalue"
    queries_from_str = load_string(s1)
    load_string!(queries_from_str, s2)
    @test length(queries_from_str) == 2
    for q in values(queries_from_str)
        @test q.args == [:myvalue]
    end
end

@testset "Test that duplicate query warning is generated" begin
    s = "-- name: string_query1\n select count(*) from my_table where mycolumn = :myvalue"
    q = load_string(s)
    @test_logs (:warn, r"string_query1 already defined and will be overwritten") load_string!(q, s)
end

@testset "Test operation, argument, and docstring extraction" begin
    ops = Dict(
        :get_musicians => YeSQL.Select(),
        :get_musician_id => YeSQL.SelectValue(),
    )
    args = Dict(
        :get_musicians => (:employed_date,),
        :get_musician_id => (:first_name, :last_name),
    )

    docstrings = Dict(
        :get_musicians => "Return musicians active on :employed_date. A superfluous line to test multi-line docstrings.",
        :get_musician_id => "Return musician id",
    )

    for qn in keys(args)
        q = queries[qn]
        @test all(a in q.args for a in args[qn])
        @test q.docstring == docstrings[qn]
        expected_query = strip(replace(read(joinpath("./sql", "expected", string(qn) * ".sql"), String),
                                       r"\s+" => " "))
        @test expected_query == q.query
        @test q.operation == ops[qn]
    end
end

@testset "Test that error is thrown for invalid path" begin
    @test_throws ArgumentError load_path("garbage.sql")
end

include("libpq_tests.jl")

ext: jl


text:
# Test for GaussRadau()

n = 1
x, w = gaussradau(n)
@test x[1] ≈ -1.
@test w[1] ≈ 2.

n = 2
x, w = gaussradau(n)
@test x[1] ≈ -1.
@test x[2] ≈ 1/3
@test w[1] ≈ .5
@test w[2] ≈ 1.5

tol = 1e-14
n = 42
ntests = 4
x, w = gaussradau(n)
@test length(x) == n && length(w) == n
@test abs(x[37] - 0.908847278001044) < tol
@test abs(w[37] - 0.031190846817016) < tol
@test x[1] == -1

@test dot( w,exp.(x)) ≈ exp(1)-exp(-1)


text:
using Printf

# Source: https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/src/trials.jl
function prettytime(t)
    if t < 1e3
        value, units = t, "ns"
    elseif t < 1e6
        value, units = t / 1e3, "μs"
    elseif t < 1e9
        value, units = t / 1e6, "ms"
    else
        s = t / 1e9
        if s < 60
            value, units = s, "s"
        else
            value, units = (s / 60), "min"
        end
    end
    return string(@sprintf("%.3f", value), " ", units)
end


text:
# Modified from
# https://github.com/tpapp/AltDistributions.jl

# See also the Stan manual (the "Stanual", though nobody calls it that)
# https://mc-stan.org/docs/2_27/reference-manual/cholesky-factors-of-correlation-matrices-1.html

export LKJCholesky
using PositiveFactorizations

const CorrCholeskyUpper = TV.CorrCholeskyFactor

@doc """
    LKJCholesky(k=3, η=1.0)
    LKJCholesky(k=3, logη=0.0)

`LKJCholesky(k, ...)` gives the `k×k` LKJ distribution (Lewandowski et al 2009)
expressed as a Cholesky decomposition. As a special case, for
`C = rand(LKJCholesky(k=K, η=1.0))` (or equivalently
`C=rand(LKJCholesky{k}(k=K, logη=0.0))`), `C.L * C.U` is uniform over the set of
all `K×K` correlation matrices. Note, however, that in this case `C.L` and `C.U`
are **not** sampled uniformly (because the multiplication is nonlinear).

The `logdensity` method for this measure applies for `LowerTriangular`,
`UpperTriangular`, or `Diagonal` matrices, and will "do the right thing". The
`logdensity` **does not check if `L*U` yields a valid correlation matrix**.

Valid values are ``η > 0``. When ``η > 1``, the distribution is unimodal with a
peak at `I`, while ``0 < η < 1`` yields a trough. ``η = 2`` is recommended as a vague prior.

Adapted from
https://github.com/tpapp/AltDistributions.jl
""" LKJCholesky

@parameterized LKJCholesky(k,η)

@kwstruct LKJCholesky(k,η)
@kwstruct LKJCholesky(k,logη)

LKJCholesky(k::Integer) = LKJCholesky(k, 1.0)

asparams(::Type{<:LKJCholesky}, ::Val{:η}) = asℝ₊
asparams(::Type{<:LKJCholesky}, ::Val{:logη}) = asℝ

# Modified from
# https://github.com/tpapp/AltDistributions.jl

using LinearAlgebra

logdensity(d::LKJCholesky, C::Cholesky) = logdensity(d, C.UL)

function logdensity(d::LKJCholesky{(:k, :η,)}, L::Union{LinearAlgebra.AbstractTriangular, Diagonal})
    η = d.η
    # z = diag(L)
    # sum(log.(z) .* ((k:-1:1) .+ 2*(η-1)))

    # Note: https://github.com/cscherrer/MeasureTheory.jl/issues/100#issuecomment-852428192
    c = d.k + 2(η - 1)
    n = size(L,1)
    s = sum(1:n) do i
        (c - i) * @inbounds log(L[i,i])
    end
    return s
end

function logdensity(d::LKJCholesky{(:k, :logη)}, L::Union{LinearAlgebra.AbstractTriangular, Diagonal})
    c = d.k + 2 * expm1(d.logη)
    n = size(L,1)
    s = sum(1:n) do i
        (c - i) * @inbounds log(L[i,i])
    end
    return s
end

TV.as(d::LKJCholesky) = CorrCholesky(d.k)

function basemeasure(μ::LKJCholesky{(:k,:η)})
    t = as(μ)
    f(par) = Dists.lkj_logc0(par.k, par.η)
    base = Pushforward(t, Lebesgue(ℝ)^dimension(t), false)
    ParamWeightedMeasure(f, (k= μ.k, η= μ.η), base)
end

function basemeasure(μ::LKJCholesky{(:k,:logη)})
    t = as(μ)
    η = exp(μ.logη)
    f(par) = Dists.lkj_logc0(par.k, exp(par.logη))
    base = Pushforward(t, Lebesgue(ℝ)^dimension(t), false)
    ParamWeightedMeasure(f, (k= μ.k, logη= logη), base)
end

# Note (from @sethaxen): this can be done without the cholesky decomposition, by randomly
# sampling the canonical partial correlations, as in the LKJ paper, and then
# mapping them to the cholesky factor instead of the correlation matrix. Stan
# implements* this. But that can be for a future PR.
#
# * https://github.com/stan-dev/math/blob/develop/stan/math/prim/prob/lkj_corr_cholesky_rng.hpp
function Base.rand(rng::AbstractRNG, ::Type, d::LKJCholesky{(:k, :η,)})
    return cholesky(Positive, rand(rng, Dists.LKJ(d.k, d.η)))
end;

function Base.rand(rng::AbstractRNG, ::Type, d::LKJCholesky{(:k, :logη)})
    return cholesky(Positive, rand(rng, Dists.LKJ(d.k, exp(d.logη))))
end;

ConstructionBase.constructorof(::Type{L}) where {L<:LKJCholesky} = LKJCholesky


text:
using Gadfly
using LinearAlgebra
using VoronoiDelaunay
using GeometricalPredicates

function distmesh(fd, fh, h0, bbox, pfix) end

# Co jest ten broadcasting aka .+ w 12

function fh(x,y)
    return 0.2
end

function generate(fd, fh, bbox, h0, pfix)
    ttol = .1
    geps = .001 * h0
    h1 = h0 * (sqrt(3) / 2)
    v1 = bbox[1, 1]:h0:bbox[2, 1]
    v2 = bbox[1, 2]:h1:bbox[2, 2]
    x, y = meshgrid(v1, v2)
    x[2:2:end, :] = x[2:2:end, :] .+ h0 / 2
    p = [x[:] y[:]]
    z = [vcat(a) for a in eachrow(p) if fd(a) < geps]
    f = transpose(reshape(vcat(z...), 2, length(z)))
    r0 = 1 ./ fh(f[:,1], f[:,2]).^2
    r0_max = r0./max(r0)
    l = [vcat(b) for b in eachrow(f) if (rand(Float64, size(f)))[1] < r0_max ]
    g = transpose(reshape(vcat(z...), 2, length(z)))
    g
    N = size(g,1)
    N
    println(N)

    pold = Inf
    scale = abs(1/(bbox[1] - bbox[2]))
    transx = 1 + (0-bbox[1]*scale)
    transy = 1 + (0-bbox[1,2]*scale)

    maximum(sqrt.(sum((g.-pold).^2,dims = 2)/h0))>ttol
    tess = DelaunayTessellation()
    aa = map(el -> Point(el[1]*scale+transx, el[2]*scale+transy), eachrow(g))
    aa

    push!(tess, aa)
    zda = Array{Point2D,1}()
    for tri in tess
        push!(zda, centroid(Primitive(geta(tri), getb(tri), getc(tri))))
    end
    zza = filter(cent -> dc([getx(cent) gety(cent)]) > -geps, zda)
    zdda = Array{Float64,1}()
    barvec = Array{Array{Float64,1},1}()
    for edge in delaunayedges(tess)
        push!(zdda,sqrt(length2(Line(geta(edge), getb(edge))))/scale)
        push!(barvec, [(getx(getb(edge))-getx(geta(edge)))/scale, ((gety(getb(edge))-gety(geta(edge))))/scale])
    end
    hbars  = 0.2
    ss = 0.2*1.2*sqrt(sum(zdda.^2)/(0.2^2))
    ttt = Array{Float64,1}()
    for edge in delaunayedges(tess)
        push!(ttt,ss - sqrt(length2(Line(geta(edge), getb(edge))))/scale)
    end
    rttt = reshape(ttt,:,1)
    mzdda = reshape(zdda, :, 1)
    init_Fvec = rttt ./ mzdda
#    Fvec = Fvec.*barvec
    l
#    println(sort(with_f[:]))

    #mzzda
#    width = max_coord - min_coord
#    a = Point2D[Point(min_coord + rand() * width, min_coord + rand() * width) for i in 1:100]
#    push!(tess, a)
#    x, y = getplotxy(delaunayedges(tess))
#    Gadfly.plot(x=x, y=y, Geom.path)
end

function evaluate(fd, p, geps)
    for x in eachrow(p)
        if(fd(x)<geps)

        end
    end
end

function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where {T}
    m, n = length(vy), length(vx)

    gx = reshape(repeat(vx, inner = m, outer = 1), m, n)
    gy = reshape(repeat(vy, inner = 1, outer = n), m, n)

    return gx, gy
end

function feval(f, x...)
end

function generateNodes()
    x = [1, 2, 3, 4]
    y = [1, 1, 2, 3]
    plot(x, y)
end

function dc(p)
    return sqrt(sum(p.^2))-1
end

function dcircle(p, xc, yc, r)
    return sqrt.((p[:, 1] .- xc) .^ 2 .+ (p[:, 2] .- yc) .^ 2) .- r
end

max_stars_repo_name: jstarczewski/julia-fem

max_issues_repo_name: jstarczewski/julia-fem

max_forks_repo_name: jstarczewski/julia-fem


text:
# Import Labelled Intervals
function (load::CSVImporter{<:CSVSchema{<: GenericVendor, <: GenericProduct, <: IntervalLabels}})(path::Path)
    disable!(load)
    df = load_dataframe(path)
    _, ncol = size(df)
    if ncol != 6
        CImGui.OpenPopup("Have you opened the appropriate file?")
        return nothing
    else
        dict = Dict{String, LabelledInterval}()
        for row in eachrow(df)
            label = row[:label]
            start = Float64(row[:start])
            stop = Float64(row[:stop])
            i₀ = row[:first]
            s = row[:step]
            i₁ = row[:last]
            interval = i₀:s:i₁
            dict[label] = LabelledInterval(label, NestedInterval(start = start, stop = stop, interval = interval))
        end
        return LabelledIntervals("Conditions", dict)
    end
end


text:
# Bond Helpers
mutable struct FixedRateBondHelper{DC <: DayCount, P <: PricingEngine} <: BondHelper
  price::Quote
  bond::FixedRateBond{DC, P}
end
value(b::FixedRateBondHelper) = b.price.value
maturity_date(b::FixedRateBondHelper) = maturity_date(b.bond)

# bond helper functions
function implied_quote(bond_h::BondHelper, clean::Bool = true)
  bond = bond_h.bond
  recalculate!(bond)
  settlement_value = bond.settlementValue
  return clean ? clean_price(bond, settlement_value, settlement_date(bond)) : dirty_price(bond, settlement_value, settlement_date(bond))
end

function clone(fixedRateBondHelper::FixedRateBondHelper, ts::TermStructure)
  # first we have to clone the PE
  newPE = clone(fixedRateBondHelper.bond.pricingEngine, ts)
  # then we have to clone the bond
  newBond = clone(fixedRateBondHelper.bond, newPE)

  return FixedRateBondHelper{typeof(newBond.dc), typeof(newPE)}(fixedRateBondHelper.price, newBond)
end

update_termstructure(bondHelper::FixedRateBondHelper, ts::TermStructure) = clone(bondHelper, ts)


text:
"""
    @extern(fun, returntype, argtypes, args...)

Creates a call to an external function meant to be included at link time.
Use the same conventions as `ccall`.

This transforms into the following `ccall`:

    ccall("extern fun", llvmcall, returntype, argtypes, args...)
"""
macro extern(name, rettyp, argtyp, args...)
    externfun = string("extern ", name isa AbstractString || name isa Symbol ? name : name.value)
    Expr(:call, :ccall, externfun, esc(:llvmcall), esc(rettyp),
         Expr(:tuple, esc.(argtyp.args)...), esc.(args)...)
end


text:
## %PREDICTMD_GENERATED_BY%

import PredictMDExtra
PredictMDExtra.import_all()

import PredictMD
PredictMD.import_all()

# PREDICTMD IF INCLUDE TEST STATEMENTS
import CSVFiles
import DataFrames
import FileIO
import RDatasets
import Random
import Test
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

# PREDICTMD IF INCLUDE TEST STATEMENTS
logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
if isnothing(logger)
    logger_stream = devnull
else
    logger_stream = logger.stream
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### Begin project-specific settings

DIRECTORY_CONTAINING_THIS_FILE = @__DIR__
PROJECT_DIRECTORY = dirname(
    joinpath(splitpath(DIRECTORY_CONTAINING_THIS_FILE)...)
    )
PROJECT_OUTPUT_DIRECTORY = joinpath(
    PROJECT_DIRECTORY,
    "output",
    )
mkpath(PROJECT_OUTPUT_DIRECTORY)
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "data"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "models"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "plots"))

# PREDICTMD IF INCLUDE TEST STATEMENTS
@debug("PROJECT_OUTPUT_DIRECTORY: ", PROJECT_OUTPUT_DIRECTORY,)
if PredictMD.is_travis_ci()
    PredictMD.cache_to_path!(
        ;
        from = ["cpu_examples", "breast_cancer_biopsy", "output",],
        to = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### End project-specific settings

### Begin Knet neural network classifier code

Random.seed!(999)

trainingandtuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_features_df.csv",
    )
trainingandtuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_labels_df.csv",
    )
testing_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_features_df.csv",
    )
testing_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_labels_df.csv",
    )
training_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_features_df.csv",
    )
training_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_labels_df.csv",
    )
tuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_features_df.csv",
    )
tuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_labels_df.csv",
    )
trainingandtuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
trainingandtuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )
testing_features_df = DataFrames.DataFrame(
    FileIO.load(
        testing_features_df_filename;
        type_detect_rows = 100,
        )
    )
testing_labels_df = DataFrames.DataFrame(
    FileIO.load(
        testing_labels_df_filename;
        type_detect_rows = 100,
        )
    )
training_features_df = DataFrames.DataFrame(
    FileIO.load(
        training_features_df_filename;
        type_detect_rows = 100,
        )
    )
training_labels_df = DataFrames.DataFrame(
    FileIO.load(
        training_labels_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )

smoted_training_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "smoted_training_features_df.csv",
    )
smoted_training_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "smoted_training_labels_df.csv",
    )
smoted_training_features_df = DataFrames.DataFrame(
    FileIO.load(
        smoted_training_features_df_filename;
        type_detect_rows = 100,
        )
    )
smoted_training_labels_df = DataFrames.DataFrame(
    FileIO.load(
        smoted_training_labels_df_filename;
        type_detect_rows = 100,
        )
    )

categorical_feature_names_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "categorical_feature_names.jld2",
    )
continuous_feature_names_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "continuous_feature_names.jld2",
    )
categorical_feature_names = FileIO.load(
    categorical_feature_names_filename,
    "categorical_feature_names",
    )
continuous_feature_names = FileIO.load(
    continuous_feature_names_filename,
    "continuous_feature_names",
    )
feature_names = vcat(categorical_feature_names, continuous_feature_names)

single_label_name = :Class
negative_class = "benign"
positive_class = "malignant"
single_label_levels = [negative_class, positive_class]

categorical_label_names = Symbol[single_label_name]
continuous_label_names = Symbol[]
label_names = vcat(categorical_label_names, continuous_label_names)

knet_mlp_predict_function_source = """
function knetmlp_predict(
        w,
        x0::AbstractArray;
        probabilities::Bool = true,
        )
    x1 = Knet.relu.( w[1]*x0 .+ w[2] )
    x2 = Knet.relu.( w[3]*x1 .+ w[4] )
    x3 = w[5]*x2 .+ w[6]
    unnormalizedlogprobs = x3
    if probabilities
        normalizedlogprobs = Knet.logp(unnormalizedlogprobs; dims = 1)
        normalizedprobs = exp.(normalizedlogprobs)
        return normalizedprobs
    else
        return unnormalizedlogprobs
    end
end
"""

knet_mlp_loss_function_source = """
function knetmlp_loss(
        predict::Function,
        modelweights,
        x::AbstractArray,
        ytrue::AbstractArray;
        L1::Real = Float64(0),
        L2::Real = Float64(0),
        )
    loss = Knet.nll(
        predict(
            modelweights,
            x;
            probabilities = false,
            ),
        ytrue;
        dims = 1,
        )
    if L1 != 0
        loss += L1 * sum(sum(abs, w_i) for w_i in modelweights[1:2:end])
    end
    if L2 != 0
        loss += L2 * sum(sum(abs2, w_i) for w_i in modelweights[1:2:end])
    end
    return loss
end
"""

feature_contrasts = PredictMD.generate_feature_contrasts(
    smoted_training_features_df,
    feature_names,
    )

knetmlp_modelweights = Any[
    Float64.(
        0.1f0*randn(Float64,64,feature_contrasts.num_array_columns_without_intercept)
        ),
    Float64.(
        fill(Float64(0),64,1)
        ),
    Float64.(
        0.1f0*randn(Float64,32,64)
        ),
    Float64.(
        fill(Float64(0),32,1)
        ),
    Float64.(
        0.1f0*randn(Float64,2,32)
        ),
    Float64.(
        fill(Float64(0),2,1)
        ),
    ]

knetmlp_losshyperparameters = Dict()
knetmlp_losshyperparameters[:L1] = Float64(0.0)
knetmlp_losshyperparameters[:L2] = Float64(0.0)

knetmlp_optimizationalgorithm = :Momentum
knetmlp_optimizerhyperparameters = Dict()
knetmlp_minibatchsize = 48

knet_mlp_classifier =
    PredictMD.single_labelmulticlassdataframeknetclassifier(
        feature_names,
        single_label_name,
        single_label_levels;
        package = :Knet,
        name = "Knet MLP",
        predict_function_source = knet_mlp_predict_function_source,
        loss_function_source = knet_mlp_loss_function_source,
        losshyperparameters = knetmlp_losshyperparameters,
        optimizationalgorithm = knetmlp_optimizationalgorithm,
        optimizerhyperparameters = knetmlp_optimizerhyperparameters,
        minibatchsize = knetmlp_minibatchsize,
        modelweights = knetmlp_modelweights,
        printlosseverynepochs = 1,
        maxepochs = 50,
        feature_contrasts = feature_contrasts,
        )

PredictMD.parse_functions!(knet_mlp_classifier)

PredictMD.fit!(
    knet_mlp_classifier,
    smoted_training_features_df,
    smoted_training_labels_df,
    tuning_features_df,
    tuning_labels_df,
    )

PredictMD.set_max_epochs!(knet_mlp_classifier, 100)

PredictMD.fit!(
    knet_mlp_classifier,
    smoted_training_features_df,
    smoted_training_labels_df,
    tuning_features_df,
    tuning_labels_df,
    )

knet_learningcurve_lossvsepoch = PredictMD.plotlearningcurve(
    knet_mlp_classifier,
    :loss_vs_epoch;
    );

# PREDICTMD IF INCLUDE TEST STATEMENTS
filename = string(
    tempname(),
    "_",
    "knet_learningcurve_lossvsepoch",
    ".pdf",
    )
rm(filename; force = true, recursive = true,)
@debug("Attempting to test that the file does not exist...", filename,)
Test.@test(!isfile(filename))
@debug("The file does not exist.", filename, isfile(filename),)
PredictMD.save_plot(filename, knet_learningcurve_lossvsepoch)
if PredictMD.is_force_test_plots()
    @debug("Attempting to test that the file exists...", filename,)
    Test.@test(isfile(filename))
    @debug("The file does exist.", filename, isfile(filename),)
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

display(knet_learningcurve_lossvsepoch)
PredictMD.save_plot(
    joinpath(
        PROJECT_OUTPUT_DIRECTORY,
        "plots",
        "knet_learningcurve_lossvsepoch.pdf",
        ),
    knet_learningcurve_lossvsepoch,
    )

knet_learningcurve_lossvsepoch_skip10epochs = PredictMD.plotlearningcurve(
    knet_mlp_classifier,
    :loss_vs_epoch;
    startat = 10,
    endat = :end,
    );

# PREDICTMD IF INCLUDE TEST STATEMENTS
filename = string(
    tempname(),
    "_",
    "knet_learningcurve_lossvsepoch_skip10epochs",
    ".pdf",
    )
rm(filename; force = true, recursive = true,)
@debug("Attempting to test that the file does not exist...", filename,)
Test.@test(!isfile(filename))
@debug("The file does not exist.", filename, isfile(filename),)
PredictMD.save_plot(filename, knet_learningcurve_lossvsepoch_skip10epochs)
if PredictMD.is_force_test_plots()
    @debug("Attempting to test that the file exists...", filename,)
    Test.@test(isfile(filename))
    @debug("The file does exist.", filename, isfile(filename),)
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

display(knet_learningcurve_lossvsepoch_skip10epochs)
PredictMD.save_plot(
    joinpath(
        PROJECT_OUTPUT_DIRECTORY,
        "plots",
        "knet_learningcurve_lossvsepoch_skip10epochs.pdf",
        ),
    knet_learningcurve_lossvsepoch_skip10epochs,
    )

knet_learningcurve_lossvsiteration = PredictMD.plotlearningcurve(
    knet_mlp_classifier,
    :loss_vs_iteration;
    window = 50,
    sampleevery = 10,
    );

# PREDICTMD IF INCLUDE TEST STATEMENTS
filename = string(
    tempname(),
    "_",
    "knet_learningcurve_lossvsiteration",
    ".pdf",
    )
rm(filename; force = true, recursive = true,)
@debug("Attempting to test that the file does not exist...", filename,)
Test.@test(!isfile(filename))
@debug("The file does not exist.", filename, isfile(filename),)
PredictMD.save_plot(filename, knet_learningcurve_lossvsiteration)
if PredictMD.is_force_test_plots()
    @debug("Attempting to test that the file exists...", filename,)
    Test.@test(isfile(filename))
    @debug("The file does exist.", filename, isfile(filename),)
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

display(knet_learningcurve_lossvsiteration)
PredictMD.save_plot(
    joinpath(
        PROJECT_OUTPUT_DIRECTORY,
        "plots",
        "knet_learningcurve_lossvsiteration.pdf",
        ),
    knet_learningcurve_lossvsiteration,
    )

knet_learningcurve_lossvsiteration_skip100iterations =
    PredictMD.plotlearningcurve(
        knet_mlp_classifier,
        :loss_vs_iteration;
        window = 50,
        sampleevery = 10,
        startat = 100,
        endat = :end,
        );

# PREDICTMD IF INCLUDE TEST STATEMENTS
filename = string(
    tempname(),
    "_",
    "knet_learningcurve_lossvsiteration_skip100iterations",
    ".pdf",
    )
rm(filename; force = true, recursive = true,)
@debug("Attempting to test that the file does not exist...", filename,)
Test.@test(!isfile(filename))
@debug("The file does not exist.", filename, isfile(filename),)
PredictMD.save_plot(filename, knet_learningcurve_lossvsiteration_skip100iterations)
if PredictMD.is_force_test_plots()
    @debug("Attempting to test that the file exists...", filename,)
    Test.@test(isfile(filename))
    @debug("The file does exist.", filename, isfile(filename),)
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

display(knet_learningcurve_lossvsiteration_skip100iterations)
PredictMD.save_plot(
    joinpath(
        PROJECT_OUTPUT_DIRECTORY,
        "plots",
        "knet_learningcurve_lossvsiteration_skip100iterations.pdf",
        ),
    knet_learningcurve_lossvsiteration_skip100iterations,
    )

knet_mlp_classifier_hist_training =
    PredictMD.plotsinglelabelbinaryclassifierhistogram(
        knet_mlp_classifier,
        smoted_training_features_df,
        smoted_training_labels_df,
        single_label_name,
        single_label_levels,
        );

# PREDICTMD IF INCLUDE TEST STATEMENTS
filename = string(
    tempname(),
    "_",
    "knet_mlp_classifier_hist_training",
    ".pdf",
    )
rm(filename; force = true, recursive = true,)
@debug("Attempting to test that the file does not exist...", filename,)
Test.@test(!isfile(filename))
@debug("The file does not exist.", filename, isfile(filename),)
PredictMD.save_plot(filename, knet_mlp_classifier_hist_training)
if PredictMD.is_force_test_plots()
    @debug("Attempting to test that the file exists...", filename,)
    Test.@test(isfile(filename))
    @debug("The file does exist.", filename, isfile(filename),)
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

display(knet_mlp_classifier_hist_training)
PredictMD.save_plot(
    joinpath(
        PROJECT_OUTPUT_DIRECTORY,
        "plots",
        "knet_mlp_classifier_hist_training.pdf",
        ),
    knet_mlp_classifier_hist_training,
    )

knet_mlp_classifier_hist_testing =
        PredictMD.plotsinglelabelbinaryclassifierhistogram(
        knet_mlp_classifier,
        testing_features_df,
        testing_labels_df,
        single_label_name,
        single_label_levels,
        );

# PREDICTMD IF INCLUDE TEST STATEMENTS
filename = string(
    tempname(),
    "_",
    "knet_mlp_classifier_hist_testing",
    ".pdf",
    )
rm(filename; force = true, recursive = true,)
@debug("Attempting to test that the file does not exist...", filename,)
Test.@test(!isfile(filename))
@debug("The file does not exist.", filename, isfile(filename),)
PredictMD.save_plot(filename, knet_mlp_classifier_hist_testing)
if PredictMD.is_force_test_plots()
    @debug("Attempting to test that the file exists...", filename,)
    Test.@test(isfile(filename))
    @debug("The file does exist.", filename, isfile(filename),)
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

display(knet_mlp_classifier_hist_testing)
PredictMD.save_plot(
    joinpath(
        PROJECT_OUTPUT_DIRECTORY,
        "plots",
        "knet_mlp_classifier_hist_testing.pdf",
        ),
    knet_mlp_classifier_hist_testing,
    )

show(
    logger_stream, PredictMD.singlelabelbinaryclassificationmetrics(
        knet_mlp_classifier,
        smoted_training_features_df,
        smoted_training_labels_df,
        single_label_name,
        positive_class;
        sensitivity = 0.95,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

show(
    logger_stream, PredictMD.singlelabelbinaryclassificationmetrics(
        knet_mlp_classifier,
        testing_features_df,
        testing_labels_df,
        single_label_name,
        positive_class;
        sensitivity = 0.95,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

knet_mlp_classifier_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "knet_mlp_classifier.jld2",
    )

PredictMD.save_model(knet_mlp_classifier_filename, knet_mlp_classifier)

# PREDICTMD IF INCLUDE TEST STATEMENTS
knet_mlp_classifier = nothing
Test.@test isnothing(knet_mlp_classifier)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### End Knet neural network classifier code

# PREDICTMD IF INCLUDE TEST STATEMENTS
if PredictMD.is_travis_ci()
    PredictMD.path_to_cache!(
        ;
        to = ["cpu_examples", "breast_cancer_biopsy", "output",],
        from = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

## %PREDICTMD_GENERATED_BY%


text:
module Allnix

function loop(count::Int64 = 1000)
    n = 1000001

    x = Array{Float64,1}(n)
    y = Array{Float64,1}(n)
    a = 0.5::Float64
    for i in 1:n
        x[i] = i-1
        y[i] = 1.0
    end

    #tic()

    @timev for j in 1:count
    #    y = a*x + y # this is slower
         @inbounds for i in 1:n
            y[i] = a*x[i] + y[i]
        end
    end
    #toc()

    println(y[n])
end

# package code goes here
include("mod1ch2.jl")
include("julia_doc.jl")
include("blas.jl")
include("Test.jl")

#@code_llvm loop(1000)

end # module


text:
@testset "Continuous space" begin

  # Basic model initialization
  space1 = ContinuousSpace(2; periodic = true, extend = (1, 1))
  space2 = ContinuousSpace(2; periodic = false, extend = (1, 1))
  space3 = ContinuousSpace(2; periodic = true, extend = (2, 1))
  space4 = ContinuousSpace(2; periodic = true, extend = (1, 2))
  space5 = ContinuousSpace(2; periodic = true, extend = (2, 2))
  space6 = ContinuousSpace(1; periodic = true, extend = (1,))
  space7 = ContinuousSpace(3; periodic = true, extend = (1,1,1))
  space8 = ContinuousSpace(3; periodic = false, extend = (1,1,1))

  @test space1.D == 2
  @test space2.D == 2
  @test space3.D == 2
  @test space4.D == 2
  @test space5.D == 2
  @test space6.D == 1
  @test space7.D == 3
  @test space8.D == 3

  model1 = ABM(Agent6, space1)
  model2 = ABM(Agent6, space2)
  model3 = ABM(Agent6, space3)
  model4 = ABM(Agent6, space4)
  model5 = ABM(Agent6, space5)
  model6 = ABM(Agent6, space6)

  @test nagents(model1) == 0
  @test Agents.collect_ids(DBInterface.execute(model1.space.db, "select id from tab")) == []

  # add_agent! with no existing agent (the agent is created)
  pos = (0.5, 0.5)
  vel = (0.2, 0.1)
  dia = 0.01
  add_agent!(pos, model1, vel, dia)
  @test Agents.collect_ids(DBInterface.execute(model1.space.db, "select id from tab")) == [1]
  dbrow = DBInterface.execute(model1.space.db, "select * from tab") |> DataFrame;
  @test dbrow[1, :a] == 0.5
  @test dbrow[1, :b] == 0.5

  # move_agent! without provided update_vel! function
  move_agent!(model1.agents[1], model1)
  dbrow = DBInterface.execute(model1.space.db, "select * from tab") |> DataFrame;
  @test dbrow[1, :a] == 0.7
  @test dbrow[1, :b] == 0.6
  @test dbrow[1, :a] == model1.agents[1].pos[1]
  @test dbrow[1, :b] == model1.agents[1].pos[2]

  kill_agent!(model1.agents[1], model1)
  dbrow = DBInterface.execute(model1.space.db, "select * from tab") |> DataFrame
  @test size(dbrow) == (0,0)

  # add_agent! with an existing agent
  agent = Agent6(2, pos, vel, dia)
  add_agent!(agent, model1)
  @test Agents.defvel(agent, model) == nothing
  @test Agents.collect_ids(DBInterface.execute(model1.space.db, "select id from tab")) == [2]

  # agents within some range are found correctly
  agent2 = model1.agents[2]
  agent3 = Agent6(3, agent2.pos .+ 0.005, vel, dia)
  add_agent_pos!(agent3, model1)
  n_ids = space_neighbors(agent2, model1, agent2.weight)
  @test length(n_ids) == 1
  @test n_ids[1] == 3
  n_ids = space_neighbors(agent2, model1, agent2.weight/10)
  @test length(n_ids) == 0

  # test that it finds both
  n_ids = space_neighbors(agent2.pos, model1, agent2.weight)
  @test sort!(n_ids) == [2, 3]

  # test various metrics
  c = ABM(Agent6, ContinuousSpace(2; extend = (2.0, 2.0), metric = :cityblock))
  e = ABM(Agent6, ContinuousSpace(2; extend = (2.0, 2.0), metric = :euclidean))
  r = sqrt(2) - 0.2
  for (i, φ) in enumerate(range(0; stop = 2π, length = 10))
    a = Agent6(i, (1, 1) .+ r .* sincos(φ), (0.0, 0.0), 0.0)
    add_agent_pos!(a, c)
    add_agent_pos!(a, e)
  end
  @test length(space_neighbors((1, 1), c, r)) == 10
  @test 4 < length(space_neighbors((1, 1), e, r)) < 10
end

@testset "Interacting pairs" begin
  space = ContinuousSpace(2, extend = (10, 10), periodic = false, metric = :euclidean)
  model = ABM(Agent6, space; scheduler = fastest)
  model = ABM(Agent6, space; scheduler = model -> sort!(collect(keys(model.agents));rev=true))
  pos = [
    (7.074386436066224, 4.963014649338054)
    (5.831962448496828, 4.926297135685473)
    (5.122087781793935, 5.300031210394806)
    (3.9715633336430156, 4.8106570045816675)
  ]
  for i in 1:4
    add_agent_pos!(Agent6(i+2, pos[i], (0.0, 0.0), 0), model)
  end
  pairs = interacting_pairs(model, 2.0, :scheduler).pairs
  @test length(pairs) == 2
  fi = [p[1] for p in pairs]
  se = [p[2] for p in pairs]
  @test fi == unique(fi)
  @test se == unique(se)
  for id in fi
    @test id ∉ se
  end
  pairs = interacting_pairs(model, 2.0, :all).pairs
  @test length(pairs) == 5
  @test (3, 6) ∉ pairs

  space2 = ContinuousSpace(2, extend = (10, 10), periodic = false, metric = :euclidean)
  model2 = ABM(Agent6, space2; scheduler = by_id)
  for i in 1:4
    add_agent_pos!(Agent6(i, pos[i], (0.0, 0.0), 0), model2)
  end
  # Note that length here is not the same as the test above with the same function
  # call. This is due to the `scheduler` order of operation
  pairs = interacting_pairs(model2, 2.0, :scheduler).pairs
  @test length(pairs) == 3
  # A more expensive search (in memory, not time), but guarantees true nearest neighbors
  pairs = interacting_pairs(model2, 2.0, :nearest).pairs
  @test length(pairs) == 2
  pairs = interacting_pairs(model2, 2.0, :all).pairs
  @test length(pairs) == 5
  @test (1, 4) ∉ pairs
end


text:
#This file implements a number of popular covariance functions
using SpecialFunctions
using IntelVectorMath
function matern( r::Float64, l::Float64, nu::Float64 )
  if r == 0
    return 1.0
  else
    return 2.0^(1.0 - nu) / gamma(nu) * (sqrt( 2.0 * nu ) * r / l )^nu * besselk( nu, sqrt( 2.0 * nu ) * r / l )
  end
end

function matern12( r, l )
  exp( - r / l )
end

function matern12!( r, l )
  r ./= -l
  IntelVectorMath.exp!(r)
end

function matern32( r, l )
  ( 1 + sqrt(3) * r / l ) * exp( - sqrt(3) * r / l )
end

function matern32!( r, l )
  r .= - sqrt(3) .* r ./ l
  rtemp = 1 .- r
  IntelVectorMath.exp!(r)
  r .*= rtemp
end

function matern52( r, l )
  ( 1 + sqrt(5) * r / l + 5 * r^2 / 3 / l^2 ) * exp( - sqrt(5) * r / l )
end

function matern52!( r, l )
  rtemp = ( 1 .+ sqrt(5) * r / l .+ 5 * r.^2 / 3 / l^2 )
  r .= -sqrt(5) .* r ./ l
  IntelVectorMath.exp!(r)
  r .*= rtemp
end

function gaussian( r::Float64, l::Float64 )
  return exp( - ( r/ l )^2 )
end

function exponential( r::Float64, l::Float64 )
  return exp( -r/l )
end

function exponential!(r, l)
  r ./= -l
  IntelVectorMath.exp!(r)
end

function invMultiquadratic( r::Float64, l::Float64, c::Float64 = 1. )
  return 1/sqrt( c + ( r / l )^2 )
end

function ratQuadratic( r::Float64, l::Float64, c::Float64 = 1. )
  return 1 - (r/l)^2/( (r/l)^2 + c )
end

function cauchy( r::Float64, l::Float64, alpha::Float64, beta::Float64 )
  return ( 1 + (r/l)^alpha )^(-beta/alpha)
end


text:
# Discrete Laplacian

function spdiagm_nonsquare(m, n, args...)
    I, J, V = SparseArrays.spdiagm_internal(args...)
    return sparse(I, J, V, m, n)
end

# returns -∇² (discrete Laplacian, real-symmetric positive-definite) on n₁×n₂ grid
function ∇²(n₁,n₂)
    o₁ = ones(n₁)
    ∂₁ = spdiagm_nonsquare(n₁+1,n₁,-1=>-o₁,0=>o₁)
    o₂ = ones(n₂)
    ∂₂ = spdiagm_nonsquare(n₂+1,n₂,-1=>-o₂,0=>o₂)
    return kron(sparse(I,n₂,n₂), ∂₁'*∂₁) + kron(∂₂'*∂₂, sparse(I,n₁,n₁))
end

max_stars_repo_path: Code/Numerical Analysis & Machine Learning/DiscreteLaplacian.jl

max_stars_repo_name: BambooFlower/Math-Scripts

max_issues_repo_name: BambooFlower/Math-Scripts

max_forks_repo_name: BambooFlower/Math-Scripts


text:
using ArgParse
using YAML
using Compat
import Compat: UTF8String, ASCIIString, view, readstring, foreach

include("parsing.jl")
include("commands.jl")

function main()
    parsed_args = parse_args(build_arg_table())

    command = parsed_args["%COMMAND%"]

    if command == "ls"
        foreach(x -> println(x), ls())
    elseif command == "exp"
        bootstrap_exp(
            parsed_args[command]["analysis_file"],
            parsed_args[command]["output_file"],
            parsed_args[command]["addprocs"],
            parsed_args[command]["debug"]
        )
    else
        bootstrap_config(
            parsed_args[command]["config_file"],
            parsed_args[command]["output_dir"],
            parsed_args[command]["no-graph"]
        )
    end
end

# fire up simulation if run using command line
if !isinteractive()
    main()
end


text:
abstract type AbstractCSGNode end

struct CSGNode <: AbstractCSGNode
    data::Union{AbstractImplicitSurface, Function}
    children::Array{CSGNode,1}
end

function evaluate(tree::CSGNode, coords)
    # tree node doesn't have children
    isempty(tree.children) && return evaluate(tree.data, coords)

    # node.data isa Function -> it has children
    setop = tree.data::Function
    return setop(tree.children..., coords)
end

AbstractTrees.children(tree::AbstractCSGNode) = tree.children
AbstractTrees.printnode(io::IO, tree::AbstractCSGNode) = print(io, tree.data)

function normal(tree::CSGNode, coords)
    return normal(evaluate(tree, coords), coords)
end

function valueandnormal(tree::CSGNode, coords)
    val = evaluate(tree, coords)
    n = normal(val, coords)
    return (value(val), n)
end

# compatibility for dual contouring
function isosurface(tree::CSGNode, isolevel, bounding_box, resolution; df=nothing)
    isosurface(p->value(evaluate(tree, p)), isolevel, bounding_box, resolution; df=df)
end

function depth(tree::AbstractCSGNode)
    childs = collect(children(tree))
    i = 0
    while true
        if childs == fill([], size(childs))
            return i
        end
        i += 1

        newchilds = Array{typeof(tree),1}(undef,0)
        for a in childs
            append!(newchilds, children(a))
        end
        childs = newchilds
    end
end

function treesize(tree::AbstractCSGNode)
    s = 0
    for _ in PreOrderDFS(tree)
        s+=1
    end
    return s
end

function numberofnodes(it::AbstractTrees.TreeIterator)
    n = 0
    for i in it
        isempty(i.children) && continue
        n += 1
    end
    return n
end

function selectfirstchildnode(it::AbstractTrees.TreeIterator, ind)
    n = 0
    for i in it
        isempty(i.children) && continue
        n += 1
        if n == ind
            return i.children[1]
        end
    end
    throw(BoundsError(it, ind))
end

function make(surf, n)
    if n == 0 || rand() > 0.7
        return rand(surf)
    else
        # more recursive calls
        op = rand(CSGOperations)
        if op == complement
            return CSGNode(op, [make(surf, n-1)])
        else
            return CSGNode(op, [make(surf, n-1), make(surf, n-1)])
        end
    end
end

function randomtree(surfaces, maxdepth::Int)
    @assert maxdepth > 0 "Maximum depth should be at least 1!"
    return make(surfaces, maxdepth)
end

"""
    evaldistance(tree::CSGNode, points)

Compute the sum of squared distances to the surface of the `tree` from every point of `points`.
"""
function evaldistance(tree::CSGNode, points)
    sum = 0.0
    for p in points
        d = value(evaluate(tree, p))
        sum+=d^2
    end
    return sum
end


text:
using Makie, GLMakie, Random
Random.seed!(3)
# patterns
# `'/'`, `'\\'`, `'-'`, `'|'`, `'x'`, and `'+'`
directions = [Vec2f0(1), Vec2f0(1, -1), Vec2f0(1, 0), Vec2f0(0, 1),
    [Vec2f0(1), Vec2f0(1, -1)], [Vec2f0(1, 0), Vec2f0(0, 1)]]
colors = [:white, :red, (:green, 0.5), :white, (:navy, 0.85),:black]
patterns = [Makie.LinePattern(direction= hatch; width = 5, tilesize=(20,20),
    linecolor = colors[indx], background_color = colors[end-indx+1])
    for (indx, hatch) in enumerate(directions)]
fig = Figure(resolution = (1200,800), fontsize = 32)
ax = Axis(fig[1,1])
for (idx, pattern) in enumerate(patterns)
    barplot!(ax, [idx], [idx*(2rand()+1)], color = pattern, strokewidth = 2)
end
ax.xticks = (1:6, ["/", "\\", "-", "|", "x", "+"])
save("strippedPatterns.png", fig)
fig


text:
using .DimensionalData

import StanSample: convert_a3d

@dim iteration XDim "iterations"
@dim chain YDim "chains"
@dim param ZDim "parameters"

function convert_a3d(a3d_array, cnames, ::Val{:dimarrays})
    psymbols= Symbol.(cnames)
    pa = permutedims(a3d_array, [1, 3, 2])

    DimArray(pa, (iteration, chain, param(psymbols)); name=:draws)
end

function convert_a3d(a3d_array, cnames, ::Val{:dimarray})
    psymbols= Symbol.(cnames)

    # Permute [draws, params, chains] to [draws, chains, params]
    a3dp = permutedims(a3d_array, [1, 3, 2])

    # Append all chains
    iters, chains, pars = size(a3dp)
    a3dpa = reshape(a3dp, iters*chains, pars)

    # Create the DimArray
    DimArray(a3dpa, (iteration, param(psymbols)); name=:draws)
end

function matrix(da::DimArray, sym::Union{Symbol, String})
    n = string.(dims(da, :param).val)
    syms = string(sym)
    sel = String[]
    for (i, s) in enumerate(n)
        if length(s) > length(syms) && syms == n[i][1:length(syms)] &&
            n[i][length(syms)+1] in ['[', '.', '_']
            append!(sel, [n[i]])
        end
    end
    length(sel) == 0 && error("$syms not in $n")
    da[param=At(Symbol.(sel))]
end


text:
module HydrostaticFreeSurfaceModels

export
    HydrostaticFreeSurfaceModel, VectorInvariant,
    ExplicitFreeSurface, ImplicitFreeSurface,
    PrescribedVelocityFields

using KernelAbstractions: @index, @kernel, Event, MultiEvent
using KernelAbstractions.Extras.LoopInfo: @unroll

using Oceananigans.Utils: launch!

import Oceananigans: fields, prognostic_fields

# This is only used by the cubed sphere for now.
fill_horizontal_velocity_halos!(args...) = nothing

#####
##### HydrostaticFreeSurfaceModel definition
#####

FreeSurfaceDisplacementField(velocities, arch, grid) = ReducedField(Center, Center, Nothing, arch, grid; dims=3)

include("compute_w_from_continuity.jl")

include("rigid_lid.jl")
include("explicit_free_surface.jl")

# Implicit solver functionality
include("implicit_free_surface.jl")
include("compute_vertically_integrated_lateral_face_areas.jl")
include("compute_vertically_integrated_volume_flux.jl")
include("implicit_free_surface_solver.jl")

include("hydrostatic_free_surface_velocity_fields.jl")
include("hydrostatic_free_surface_tendency_fields.jl")
include("hydrostatic_free_surface_model.jl")
include("show_hydrostatic_free_surface_model.jl")
include("set_hydrostatic_free_surface_model.jl")

#####
##### Time-stepping HydrostaticFreeSurfaceModels
#####

"""
    fields(model::HydrostaticFreeSurfaceModel)

Returns a flattened `NamedTuple` of the fields in `model.velocities` and `model.tracers`.
"""
fields(model::HydrostaticFreeSurfaceModel) = merge(prognostic_fields(model), model.auxiliary_fields)

"""
    prognostic_fields(model::HydrostaticFreeSurfaceModel)

Returns a flattened `NamedTuple` of the prognostic fields associated with `HydrostaticFreeSurfaceModel`.
"""
prognostic_fields(model::HydrostaticFreeSurfaceModel) =
    hydrostatic_prognostic_fields(model.velocities, model.free_surface, model.tracers)

hydrostatic_prognostic_fields(velocities, free_surface, tracers) = merge((u = velocities.u,
                                                                          v = velocities.v,
                                                                          η = free_surface.η),
                                                                          tracers)

displacement(free_surface) = free_surface.η
displacement(::Nothing) = nothing

include("barotropic_pressure_correction.jl")
include("hydrostatic_free_surface_advection.jl")
include("hydrostatic_free_surface_tendency_kernel_functions.jl")
include("calculate_hydrostatic_free_surface_tendencies.jl")
include("update_hydrostatic_free_surface_model_state.jl")
include("hydrostatic_free_surface_ab2_step.jl")
include("prescribed_hydrostatic_velocity_fields.jl")

#####
##### Some diagnostics
#####

include("vertical_vorticity_field.jl")

end # module


text:
abstract type AbstractPolicy end
abstract type AbstractPolicyParameters end
abstract type ExecutionPolicy end

# function initialize!(policy::AbstractPolicy,x₀::Vector{Float64})
#     error("No initialize! function for policy ", typeof(policy))
# end

function initialize!() end

params(policy::AbstractPolicy) = nothing

function init_wrapper(xbegin::Ptr{Cdouble},xend::Ptr{Cdouble},policy_data::Ptr{Cvoid})::Nothing
    policy = unsafe_pointer_to_objref(policy_data)::AbstractPolicy
    ptrdiff = Int(xend - xbegin)
    N = divrem(ptrdiff, sizeof(Cdouble))[1]
    x₀ = unsafe_wrap(Array, xbegin, N)
    initialize!(policy, x₀)
    nothing
end


text:
@doc """
### EulerEquationMod.calcIsentropicVortex

  This function calculates the isentropic vortex solution to the Euler
  equations at a node.  This function uses an inner radius of 1, an inner
  radius density of 2, an inner radius Mach number of 0.95, and an inner radius
  pressure of 1/gamma.  The denstiy as a function of radius r can be found in,
  for example,
  "Output Error Estimates for Summation-by-parts Finite-difference Schemes",
  JE Hicken.

  This function uses conservative variables regardless of the static parameter
  of params.

  Inputs:
    coords: a vector of length 2 containing the x and y coordinates of the point
    params: the params object.

  Inputs/Outputs:
    sol: vector of length 4 to be populated with the solution

  Aliasing restrictions: none

"""->
function calcIsentropicVortex(params::ParamType2,
                  coords::AbstractArray{Tmsh},
                  sol::AbstractVector{Tsol}) where {Tmsh, Tsol}
  # calculates the solution at a point of the isentropic vortex
  # 2D only

  # unpack arguments
  x = coords[1]
  y = coords[2]

  # get some values out of eqn
  cv = params.cv
  R = params.R
  gamma = params.gamma

  # the (hard coded) parameters are
  r_in = 1  # inner radius of sector of circle
  rho_in = 2 # density at inner radius
  M_in = 0.95  # Mach number
  p_in =  1/gamma

  # calculate r, theta coordinates from x,y
  r = sqrt(x*x + y*y)
  theta = atan2(y, x)  # angle in radians

  # calculate values at r radius
  tmp1 = ((gamma-1)/2)*M_in*M_in
  rho_r = rho_in*(1 + tmp1*(1- (r_in*r_in)/(r*r)))^(1/(gamma-1))

  p_r = p_in*(rho_r/rho_in)^gamma

  a_r = sqrt( gamma*p_r/rho_r )

  M_r = sqrt( (2/(gamma-1))*((rho_in/rho_r)^(gamma-1))*(1 + tmp1) - 2/(gamma-1) )
  U_r = M_r*a_r  # velocity magnitude

  u_r = U_r*sin(theta)
  v_r = -U_r*cos(theta)
  e_r = cv*p_r/(rho_r*R)
  E_r = rho_r*e_r + 0.5*rho_r*U_r*U_r

  # save solution to sol
  sol[1] = rho_r
  sol[2] = rho_r*u_r
  sol[3] = rho_r*v_r
  sol[4] = E_r

  return nothing

end

"""
  Reverse mode
"""
function calcIsentropicVortex_rev(params::ParamType2,
                  coords::AbstractArray{Tmsh}, coords_bar::AbstractVector{Tres},
                  sol_bar::AbstractVector{Tsol}) where {Tmsh, Tsol, Tres}
  # calculates the solution at a point of the isentropic vortex
  # 2D only

  # unpack arguments
  x = coords[1]
  y = coords[2]

  # get some values out of eqn
  cv = params.cv
  R = params.R
  gamma = params.gamma
  gamma_1 = params.gamma_1

  # the (hard coded) parameters are
  r_in = 1  # inner radius of sector of circle
  rho_in = 2 # density at inner radius
  M_in = 0.95  # Mach number
  p_in =  1/gamma

  # calculate r, theta coordinates from x,y
  r = sqrt(x*x + y*y)
  theta = atan2(y, x)  # angle in radians

  # calculate values at r radius
  tmp1 = ((gamma-1)/2)*M_in*M_in
  rho_r = rho_in*(1 + tmp1*(1- (r_in*r_in)/(r*r)))^(1/(gamma-1))

  p_r = p_in*(rho_r/rho_in)^gamma

  a_r = sqrt( gamma*p_r/rho_r )

  M_r = sqrt( (2/(gamma-1))*((rho_in/rho_r)^(gamma-1))*(1 + tmp1) - 2/(gamma-1) )
  U_r = M_r*a_r  # velocity magnitude

  u_r = U_r*sin(theta)
  v_r = -U_r*cos(theta)
  e_r = cv*p_r/(rho_r*R)
  E_r = rho_r*e_r + 0.5*rho_r*U_r*U_r

  # save solution to sol
  #sol[1] = rho_r
  #sol[2] = rho_r*u_r
  #sol[3] = rho_r*v_r
  #sol[4] = E_r

  #----------------------------------------------------------------------------
  # reverse sweep

  E_r_bar = zero(Tres)
  rho_r_bar = zero(Tres)
  v_r_bar = zero(Tres)
  u_r_bar = zero(Tres)
  U_r_bar = zero(Tres)
  e_r_bar = zero(Tres)
  p_r_bar = zero(Tres)
  a_r_bar = zero(Tres)
  M_r_bar = zero(Tres)
  r_bar   = zero(Tres)
  x_bar   = zero(Tres)
  y_bar   = zero(Tres)
  theta_bar = zero(Tres)

  E_r_bar   += sol_bar[4]

  rho_r_bar += v_r*sol_bar[3]
  v_r_bar   += rho_r*sol_bar[3]

  rho_r_bar += u_r*sol_bar[2]
  u_r_bar   += rho_r*sol_bar[2]

  rho_r_bar += sol_bar[1]

  # E_r through u_r
  rho_r_bar += (e_r + 0.5*U_r*U_r)*E_r_bar
  e_r_bar   += rho_r*E_r_bar
  U_r_bar   += rho_r*U_r*E_r_bar

  p_r_bar   += cv/(rho_r*R)*e_r_bar
  rho_r_bar += -e_r/rho_r*e_r_bar

  U_r_bar +=  sin(theta)*u_r_bar
  U_r_bar +=  -cos(theta)*v_r_bar
  theta_bar += U_r*sin(theta)*v_r_bar + U_r*cos(theta)*u_r_bar

  # calculating values at r radius
  a_r_bar += M_r*U_r_bar
  M_r_bar += a_r*U_r_bar

  rho_r_bar += (M_r_bar/M_r)*((rho_in/rho_r)^(gamma - 2))*(-rho_in/(rho_r*rho_r))*(1 + tmp1)
#  tmp1_bar   = M_r_bar/(M_r*gamma_1)*(rho_in/rho_r)^(gamma_1)

  p_r_bar   += (0.5/a_r)*gamma*a_r_bar/rho_r
  rho_r_bar += (0.5/a_r)*(-gamma*p_r/(rho_r*rho_r))*a_r_bar

  rho_r_bar += p_in*gamma*((rho_r/rho_in)^gamma_1)*p_r_bar/rho_in

  t2 = (rho_in/gamma_1)*(1 + tmp1*(1 - (r_in*r_in)/(r*r)))^(1/gamma_1 - 1)
#  tmp1_bar  += t2*(1 - (r_in*r_in)/(r*r))
  r_bar += rho_r_bar*t2*2*tmp1*(r_in*r_in)/(r*r*r)

  # r, theta
  theta, y_bar, x_bar = atan2_rev(y, x, theta_bar)
  x_bar += (0.5/r)*2*x*r_bar
  y_bar += (0.5/r)*2*y*r_bar

  coords_bar[1] += x_bar
  coords_bar[2] += y_bar

  return nothing
end

function calcIsentropicVortex(params::ParamType3,
                  coords::AbstractArray{Tmsh},
                  sol::AbstractVector{Tsol}) where {Tmsh, Tsol}
# calculates the solution at a point of the isentropic vortex

  # unpack arguments
  x = coords[1]; x_orig = x
  y = coords[2]; y_orig = y
  z = coords[3]; z_orig = z;

  #println("x = ", x, ", y = ", y, ", z = ", z)

  phi_z = pi/4  # angle of axis relative to z axis
  sgn = 1.0

  # calculate new x coordinate in rotated coordinate system
  theta1 = atan2(z, x)
  #println("theta1 = ", theta1)
  phi2 = 0.5*pi - theta1
  #println("phi2 = ", phi2)
  r_xz = sqrt(x*x + z*z)  # distance from origin to point in xz plane
  #println("r_xz = ", r_xz)
  x = r_xz*sin(phi_z + phi2)  # x coordinate in rotated system
  #println("x = ", x)
  #println("diff = ", abs(x - z_orig))

  theta3 = theta1 + phi_z + phi2 - 0.5*pi

  # get some values out of params
  cv = params.cv
  R = params.R
  gamma = params.gamma

  # the (hard coded) parameters are
  r_in = 1  # inner radius of sector of circle
  rho_in = 2 # density at inner radius
  M_in = 0.95  # Mach number
  p_in =  1/gamma

  # calculate r, theta coordinates from x,y
  r = sqrt(x*x + y*y)
  if phi_z > 0
    theta = atan2(x, y)
  else
    theta = atan2(y, x)
  end

  # calculate values at r radius
  tmp1 = ((gamma-1)/2)*M_in*M_in
  rho_r = rho_in*(1 + tmp1*(1- (r_in*r_in)/(r*r)))^(1/(gamma-1))
  p_r = p_in*(rho_r/rho_in)^gamma
  a_r = sqrt( gamma*p_r/rho_r )
  M_r = sqrt( (2/(gamma-1))*((rho_in/rho_r)^(gamma-1))*(1 + tmp1) - 2/(gamma-1) )
  U_r = M_r*a_r  # velocity magnitude

  if phi_z > 0
    v_r = U_r*sin(theta)
    u_r = -U_r*cos(theta)
  else
    u_r = U_r*sin(theta)
    v_r = -U_r*cos(theta)
  end
  e_r = cv*p_r/(rho_r*R)
  E_r = rho_r*e_r + 0.5*rho_r*U_r*U_r

  # convert velocities back to base frame of reference
  w_r = u_r*sin(theta3)
  u_r = u_r*cos(theta3)

  # save solution to sol
  sol[1] = rho_r
  sol[2] = rho_r*u_r
  sol[3] = rho_r*v_r
  sol[4] = rho_r*w_r
  sol[5] = E_r

  return nothing

end

@doc """
### EulerEquationMod.calcFreeStream

  This function calculates the free stream solution for an airfoil problem
  based on the angle of attack and Mach number in nondimensionalized variables.

  Density and energy are set to params.rho_free (usually 1.0) and params.E_free,
  (usually 1/(gamma*gamma_1) + 0.5*Ma*Ma), and the x and y momenta as

  rho*Ma*cos(angle of attack)  and rho*Ma*sin(angle of attack).

  The angle of attack must be in radians.

  This function uses conservative variables regardless of the static parameter
  of params.

  Inputs:
    coords: a vector of length 2 containing the x and y coordinates of the point
    params: the params object.

  Inputs/Outputs:
    sol: vector of length 4 to be populated with the solution

  Aliasing restrictions: none

"""->
function calcFreeStream(params::ParamType2,
            coords::AbstractArray{Tmsh, 1},
            sol::AbstractArray{Tsol, 1}) where {Tmsh, Tsol}
# calculate the free stream conditions using the fields of params

  rho = sol[1] = params.rho_free
  E = sol[4] = params.E_free

  Ma = params.Ma

  sol[2] = rho*Ma*cos(params.aoa)
  sol[3] = rho*Ma*sin(params.aoa)

  return nothing
end

function calcFreeStream(params::ParamType3,
            coords::AbstractArray{Tmsh, 1},
            sol::AbstractArray{Tsol, 1}) where {Tmsh, Tsol}
# calculate the free stream conditions using the fields of params

  # calculate the free stream conditions using the fields of params

  rho = sol[1] = params.rho_free
  E = sol[5] = params.E_free

  Ma = params.Ma

  sol[2] = rho*Ma*cos(params.aoa)
  sol[3] = 0.0
  sol[4] = -rho*Ma*sin(params.aoa)

  return nothing
end

"""
  Like [`calcFreeStream`](@ref), but assumes zero angle of attack, regardless
  of what the options dictionary says.
"""
function calcFreeStream0(params::ParamType2,
            coords::AbstractArray{Tmsh, 1},
            sol::AbstractArray{Tsol, 1}) where {Tmsh, Tsol}
# calculate the free stream conditions using the fields of params

  rho = sol[1] = params.rho_free
  E = sol[4] = params.E_free

  Ma = params.Ma

  sol[2] = rho*Ma
  sol[3] = 0

  return nothing
end

@doc """
### EulerEquationMod.calcFreeStream_daoa

  This function calculates the free stream solution for an airfoil problem
  based on the angle of attack and Mach number in nondimensionalized variables.

  Density and energy are set to params.rho_free (usually 1.0) and params.E_free,
  (usually 1/(gamma*gamma_1) + 0.5*Ma*Ma), and the x and y momenta as

  rho*Ma*cos(angle of attack)  and rho*Ma*sin(angle of attack).

  The angle of attack must be in radians.

  This function uses conservative variables regardless of the static parameter
  of params.

  Inputs:
    coords: a vector of length 2 containing the x and y coordinates of the point
    params: the params object.

  Inputs/Outputs:
    sol: vector of length 4 to be populated with the solution

  Aliasing restrictions: none

"""->

function calcFreeStream_dAlpha(params::ParamType2,
                   coords::AbstractArray{Tmsh, 1},
                   sol::AbstractArray{Tsol, 1}) where {Tmsh, Tsol}
# calculate the free stream conditions using the fields of params

  rho = params.rho_free
  Ma = params.Ma

  sol[2] = -rho*Ma*sin(params.aoa)
  sol[3] = -rho*Ma*cos(params.aoa)

  return nothing
end

function calcFreeStream_daoa(params::ParamType3,
                 coords::AbstractArray{Tmsh, 1},
                 sol::AbstractArray{Tsol, 1}) where {Tmsh, Tsol}

  # calculate the free stream conditions using the fields of params

  rho = params.rho_free
  Ma = params.Ma

  sol[2] = -rho*Ma*sin(params.aoa)
  sol[3] = 0.0
  sol[4] = -rho*Ma*cos(params.aoa)

  return nothing
end

@doc """
### EulerEquationMod.calcUnsteadyVortex

  This function calculates the unsteady vortex solution to the Euler equations
  at time params.t, where the vortex was centered at x = params.vortex_x0 at
  time t=0.  The analytical solution can be found in, for example,
  K. Mattsson et al./ Computers & Fluxs 36 (2007) 636-649

  This function uses conservative variables regardless of the static parameter
  of params.

  Inputs:
    coords: a vector of length 2 containing the x and y coordinates of the point
    params: the params object.

  Inputs/Outputs:
    sol: vector of length 4 to be populated with the solution

  Aliasing restrictions: none

"""->
function calcUnsteadyVortex(params::ParamType{Tdim},
          coords::AbstractArray{Tmsh, 1},
          sol::AbstractArray{Tsol, 1}) where {Tmsh, Tsol, Tdim}

  fval = vortex_f(coords, params)
  t = params.t
  epsilon = params.vortex_strength
  Ma = params.Ma
  gamma = params.gamma
  gamma_1 = params.gamma_1
  x = coords[1]
  y = coords[2]
  x0 = params.vortex_x0

  coeff1 = epsilon*epsilon*gamma_1*Ma*Ma/(8*pi*pi)
  rho = (1 -coeff1* e^(fval))^(1/gamma_1)

  coeff2 = epsilon*y/(2*pi)
  u = 1 - coeff2*e^(fval/2)

  coeff3 = epsilon*(x - x0 - t)/(2*pi)
  v = coeff3*e^(fval/2)

  q2 = rho*u
  q3 = rho*v

  # calculate E
  # the paper gives pressure, so calculate E from that and the
  # pressure expression in terms of E
  term1 = (rho^gamma)/(gamma_1*gamma*Ma*Ma)
  term2 = 0.5*(q2*q2 + q3*q3)/rho
  E = term1 + term2

  sol[1] = rho
  sol[2] = q2
  sol[3] = q3

  if Tdim == 2
    sol[4] = E
  else
    sol[4] = 0
    sol[5] = E
  end

  return nothing
end

function vortex_f(coords, params)
  t = params.t
  x = coords[1]
  y = coords[2]
  x0 = params.vortex_x0
  return 1 - ( (  (x-x0) - t)^2 + y*y)
end

"""
  Unsteady vortex for Carpenters paper "Entropy Stable Staggered Grid Spectral
  Collocation for hte Burgers and Compressible Navier-Stokes Equations

  It is similar to the other unsteady vortex, except the vortex travels
  at an angle

  Does *not* use params.Ma or params.vortex_strength, or vortex.x0

"""
function calcUnsteadyVortex2(params::ParamType{Tdim},
          coords::AbstractArray{Tmsh, 1},
          sol::AbstractArray{Tsol, 1}) where {Tmsh, Tsol, Tdim}
  t = params.t
  gamma = params.gamma
  gamma_1 = params.gamma_1
  x = coords[1]
  y = coords[2]
  x0 = 0.0
  y0 = 0
  epsilon = 5.0
  Ma = 0.5
#  alpha = 45.0  # degrees
  sind_alpha = 0.7071067811865476
  cosd_alpha = 0.7071067811865476
#  cinf = 1.0  # assumption, I think this is a free parameter
#  Uinf = Ma*cinf

  Uinf = 1.0
  cinf = Ma/Uinf

  ycoeff = y - y0 - Uinf*sind_alpha*t
  xcoeff = x - x0 - Uinf*cosd_alpha*t
  f = 1 - ( xcoeff^2 + ycoeff^2 )

  coeff1 = epsilon*epsilon*Ma*Ma*gamma_1/(8*pi*pi)
  T = 1 - coeff1*exp(f)
#  T = 1 - coeff1*e^f
  rho = T^(1/gamma_1)

#  u1 = Uinf - (epsilon*y/(2*pi))*e^(f/2)
  u1 = Uinf*cosd_alpha - (epsilon*ycoeff/(2*pi))*exp(f/2)
#  u2 = (epsilon*xcoeff/(2*pi))*e^(f/2)
  u2 = Uinf*sind_alpha + ((epsilon*xcoeff)/(2*pi))*exp(f/2)

  q2 = rho*u1
  q3 = rho*u2

#  p = rho*params.R*T
#  p = (rho^gamma)/(gamma*Ma*ma)
#  term1 = (p)/(gamma_1*gamma*Ma*Ma)
  term1 = (rho^gamma)/(gamma_1*gamma*Ma*Ma)
  term2 = 0.5*(q2*q2 + q3*q3)/rho
  E = term1 + term2

  sol[1] = rho
  sol[2] = q2
  sol[3] = q3
  sol[4] = E

  return nothing
end

@doc """
### EulerEquationMod.calcRho1Energy2

  This function sets the density at a node to 1, energy to 2 and the momenta to
  zero.

  This function uses conservative variables regardless of the static parameter
  of params.

  Inputs:
    coords: a vector of length 2 containing the x and y coordinates of the point
    params: the params object.

  Inputs/Outputs:
    sol: vector of length 4 to be populated with the solution

  Aliasing restrictions: none

"""->
function calcRho1Energy2(params::ParamType2,
             coords::AbstractArray{Tmsh, 1},
             sol::AbstractArray{Tsol,1}) where {Tmsh, Tsol}
  # for square test case with rho = 1, everything else  = 0

  sol[1] = 1.0
  sol[2] = 0.0
  sol[3] = 0.0
  sol[4] = 2.0

  return nothing
end
function calcRho1Energy2(params::ParamType3,
             coords::AbstractArray{Tmsh, 1},
             sol::AbstractArray{Tsol,1}) where {Tmsh, Tsol}
  # for square test case with rho = 1, everything else  = 0

  sol[1] = 1.0
  sol[2] = 0.0
  sol[3] = 0.0
  sol[4] = 0.0
  sol[5] = 2.0

  return nothing
end

@doc """
### EulerEquationMod.calcOnes

  This function sets all the solution variables at a node to 1.0

  This function uses conservative variables regardless of the static parameter
  of params.

  Inputs:
    coords: a vector of length 2 containing the x and y coordinates of the point
    params: the params object.

  Inputs/Outputs:
    sol: vector of length 4 to be populated with the solution

  Aliasing restrictions: none

"""->
function calcOnes(params::ParamType2,
      coords::AbstractArray{Tmsh, 1},
      sol::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

  fill!(sol, 1.0)

  return nothing
end  # end function calcOnes

@doc """
### EulerEquationMod.calcZeros

  This function sets all the solution variables at a node to 0.0

  This function uses conservative variables regardless of the static parameter
  of params.

  Inputs:
    coords: a vector of length 2 containing the x and y coordinates of the point
    params: the params object.

  Inputs/Outputs:
    sol: vector of length 4 to be populated with the solution

  Aliasing restrictions: none

"""->

function calcZeros(params::ParamType2,
       coords::AbstractArray{Tmsh, 1},
       sol::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

  fill!(sol, 0.0)

  return nothing
end  # end function calcZeros

@doc """
### EulerEquationMod.calcRho1Energy2U1VW0

  Sets the density values 1.0, x momentum to 1.0,
  v & w momenta to 0.0, and energy to 2.0 at a node.

  It should work for 2D and 3D meshes.

  This function uses conservative variables regardless of the static parameter
  of params.

  Inputs:
    coords: a vector of length 2 containing the x and y coordinates of the point
    params: the params object.

  Inputs/Outputs:
    sol: vector of length 4 to be populated with the solution

  Aliasing restrictions: none

"""->
function calcRho1Energy2U1VW0(params::ParamType2,
                           coords::AbstractArray{Tmsh},
                           sol::AbstractArray{Tsol, 1}) where {Tmsh, Tsol}
  # for square test case with rho = 1, digonal momentum, energy

  sol[1] = 1.0
  sol[2] = 1.0
  sol[3] = 0.0
  sol[4] = 2.0

  return nothing
end

function calcRho1Energy2U1VW0(params::ParamType3,
                           coords::AbstractArray{Tmsh},
                           sol::AbstractArray{Tsol, 1}) where {Tmsh, Tsol}
  # for square test case with rho = 1, digonal momentum, energy

  sol[1] = 1.0
  sol[2] = 1.0
  sol[3] = 0.0
  sol[4] = 0.0
  sol[5] = 2.0

  return nothing
end

@doc """
### EulerEquationMod.calcRho1Energy2U3

  Sets the density values 1.0, momenta to 0.35355, and
  energy to 2.0 at a node.

  This function uses conservative variables regardless of the static parameter
  of params.

  Inputs:
    coords: a vector of length 2 containing the x and y coordinates of the point
    params: the params object.

  Inputs/Outputs:
    sol: vector of length 4 to be populated with the solution

  Aliasing restrictions: none

"""->
function calcRho1Energy2U3(params::ParamType2,
               coords::AbstractArray{Tmsh},
               sol::AbstractArray{Tsol, 1}) where {Tmsh, Tsol}
  # for square test case with rho = 1, digonal momentum, energy

  sol[1] = 1.0
  sol[2] = 0.35355
  sol[3] = 0.35355
  sol[4] = 2.0

  return nothing
end

function calcRho1Energy2U3(params::ParamType3,
               coords::AbstractArray{Tmsh},
               sol::AbstractArray{Tsol, 1}) where {Tmsh, Tsol}
  # for square test case with rho = 1, digonal momentum, energy

  sol[1] = 1.0
  sol[2] = 0.35355
  sol[3] = 0.35355
  sol[4] = 0.35355
  sol[5] = 2.0

  return nothing
end

"""
  Alias for [`calcRho1Energy2U3`](@ref) used by initial conditions.
"""
const calcRho1E2U3 = calcRho1Energy2U3

@doc """
### EulerEquationMod.calcVortex

  Sets the density 1.0, energy to 2.0 at a node.  The momenta are calculated
  according to solid body rotation with an angular velocity of 0.5 centered
  at x = 0.

  This function uses conservative variables regardless of the static parameter
  of params.

  Inputs:
    coords: a vector of length 2 containing the x and y coordinates of the point
    params: the params object.

  Inputs/Outputs:
    sol: vector of length 4 to be populated with the solution

  Aliasing restrictions: none

"""->
function calcVortex(params::ParamType2,
        coords::AbstractArray{Tmsh,1},
        sol::AbstractArray{Tsol,1}) where {Tmsh, Tsol}
# solid body rotation
  x = coords[1]
  y = coords[2]

  r = sqrt(x*x + y*y)
  theta = atan2(y, x)

  omega = 0.5

  u_norm = omega*r

  u = -u_norm*sin(theta)
  v = u_norm*cos(theta)

  rho = 1.0
  E = 2.0

  sol[1] = rho
  sol[2] = rho*u
  sol[3] = rho*v
  sol[4] = E

  return nothing
end

"""
  Calculates a manufactured solution based on exponentials
"""
function calcExp(params::ParamType2, coords::AbstractArray{Tmsh,1},
                 q::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

  x = coords[1]
  y = coords[2]
  # a and b are parameters determine the scale and offset of the solution
  # the must be the same here and in the source term
  af = 1/5  # a = 1/af
  b = 0.01
  gamma_1 = params.gamma_1
  q[1] = exp(af*x*y + b)
  q[2] = exp(af*2*x*y + b)
  q[3] = exp(af*3*x*y + b)
  q[4] = (1/gamma_1 + 0.5)*exp(af*5*x*y + b) + 0.5*exp(af*3*x*y + b)

  return nothing
end

function calcExp_rev(params::ParamType2, coords::AbstractArray{Tmsh, 1},
                 coords_bar::AbstractArray{Tmsh, 1},
                 q_bar::AbstractArray{Tsol, 1}) where {Tmsh, Tsol}

  af = 1/5  # a = 1/af
  b = 0.01
  gamma_1 = params.gamma_1

  x = coords[1]
  y = coords[2]

  x_bar = zero(Tmsh)
  y_bar = zero(Tmsh)

  # some reused quantities
  t1 = (1/gamma_1 + 0.5)
  e1 = exp(af*x*y + b)
  e2 = exp(af*2*x*y + b)
  e3 = exp(af*3*x*y + b)
  e5 = exp(af*5*x*y + b)

  x_bar += (t1*af*5*y*e5 + 0.5*af*3*y*e3)*q_bar[4]
  y_bar += (t1*af*5*x*e5 + 0.5*af*3*x*e3)*q_bar[4]

  x_bar += af*3*y*e3*q_bar[3]
  y_bar += af*3*x*e3*q_bar[3]

  x_bar += af*2*y*e2*q_bar[2]
  y_bar += af*2*x*e2*q_bar[2]

  x_bar += af*y*e1*q_bar[1]
  y_bar += af*x*e1*q_bar[1]

  coords_bar[1] += x_bar
  coords_bar[2] += y_bar

  return nothing
end

function calcExp(params::ParamType3, coords::AbstractArray{Tmsh,1},
                 q::AbstractArray{Tsol,1}) where {Tmsh, Tsol}
  x = coords[1]
  y = coords[2]
  z = coords[3]

  # constant parameters
  a = MMSExp_a
  b = MMSExp_b
  c1 = MMSExp_c1
  c2 = MMSExp_c2
  c3 = MMSExp_c3
  c4 = MMSExp_c4
  c5 = MMSExp_c5
  d1 = MMSExp_d1
  d2 = MMSExp_d2
  d3 = MMSExp_d3
  d4 = MMSExp_d4
  d5 = MMSExp_d5

  # af = 1/50, b = 0.01, c = 8, d = 0.5, f = 10 works
  gamma_1 = params.gamma_1

  t2 = exp(b);
  t3 = a*c1*x*y*z;
  q[1] = d1*t2*exp(t3);
  q[2] = d2*t2*exp(a*c2*x*y*z);
  q[3] = d3*t2*exp(a*c3*x*y*z);
  q[4] = d4*t2*exp(a*c4*x*y*z);
  q[5] = (t2*exp(-t3)*((d2*d2)*exp(a*c2*x*y*z*2.0)+(d3*d3)*exp(a*c3*x*y*z*2.0)+(d4*d4)*exp(a*c4*x*y*z*2.0))*(1.0/2.0))/d1+(d5*t2*exp(a*c5*x*y*z))/gamma_1;

  return nothing
end

"""
  Calculates a manufactured solution from Gassner, Winters, Kopriva: Split
  Form Nodal DG Schemes with SBP Propertiy for the Compressible Euler Equations.
  This is typically used with a mesh that spans [-1, 1] in all directions

"""
function calcPeriodicMMS(params::ParamType2,
             coords::AbstractArray{Tmsh,1},
             q::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

  x = coords[1]
  y = coords[2]
  t = params.t

  t7 = t*2.0;
  t2 = -t7+x+y;
  t3 = 3.141592653589793*t2;
  t4 = sin(t3);
  t5 = t4*(1.0/1.0E1);
  t6 = t5+2.0;
  q[1] = t6;
  q[2] = t6;
  q[3] = t6;
  q[4] = t6*t6;

  return nothing
end

function calcPeriodicMMS(params::ParamType3,
             coords::AbstractArray{Tmsh,1},
             q::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

  x = coords[1]
  y = coords[2]
  z = coords[3]
  t = params.t
#=
  rho = 2 + 0.1*sin(pi*(x + y + z))
#  rho = 2 + 0.5*(x*x + y*y + z*z)
  q[1] = rho
  q[2] = 0
  q[3] = 0
  q[4] = 0
  q[5] = 10
=#
  t7 = t*2.0;
  t2 = -t7+x+y+z;
  t3 = 3.141592653589793*t2;
  t4 = sin(t3);
  t5 = t4*(1.0/1.0E1);
  t6 = t5+2.0;

  q[1] = t6
  q[2] = t6
  q[3] = t6
  q[4] = t6
  q[5] = t6*t6

  return nothing
end

"""
  Manufactured condition for a channel with y in [0, 1].
  No source term.  2D only
"""
function calcChannelMMS(params::ParamType2,
            coords::AbstractArray{Tmsh,1},
            q::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

  R = params.R
  gamma_1 = params.gamma_1

  rho_inf = 1
  u_inf = 1
  offset = 0.1
  T_inf = 1  # make this smaller to reduce Mach number

  y = coords[2]

  q[1] = rho_inf
  q[2] = rho_inf*u_inf*( y*(1-y) + offset)
  q[3] = 0
  q[4] = 0.5*rho_inf*u_inf*u_inf*( offset - y*(y-1))^2 + (R*T_inf*rho_inf)/gamma_1

  return nothing
end

"""
  One dimensional square wave, 2D only

  The base flow is uniform, and if  -0.05 < x < 0.05, then it is perturbed
  by a small amount.
"""
function calcSquare1D(params::ParamType2,
            coords::AbstractArray{Tmsh,1},
            q::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

  # apply base state
  q[1] = 1.0
  q[2] = 0.3
  q[3] = 0.3
  q[4] = 5.0

  x = coords[1]
  if x > -0.05 && x < 0.05
    # apply pertubration
    q[1] += 0.1
    q[2] += 0.1
    q[3] += 0.1
    q[4] += 0.1
  end

  return nothing
end

"""
  Two dimensional square wave, 2D only

  The base flow is uniform, and if  -0.05 < x, y < 0.05, then it is perturbed
  by a small amount.
"""
function calcSquare2D(params::ParamType2,
            coords::AbstractArray{Tmsh,1},
            q::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

  # apply base state
  q[1] = 1.0
  q[2] = 0.3
  q[3] = 0.3
  q[4] = 5.0

  x = coords[1]
  y = coords[2]
  if x > 7.5 && x < 12.5 && y > -1.25 && y < 1.25
    # apply pertubration
    q[1] += 0.1
    q[2] += 0.1
    q[3] += 0.1
    q[4] += 0.1
  end

  return nothing
end

"""
  Sedov explosion, 2D only

  Uniform fluid at rest with energetic region in the circle of radius
  0.05 centered at the origin

"""
function calcSedovExplosion(params::ParamType2,
                        coords::AbstractArray{Tmsh,1},
                        q::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

  # parameters
  E = 1.0 # total amount of energy contained in the blast (not the usual E)
  nu = 2  # cylindrical blast wave
  dr = 0.05  # radius of initial blast

  q[1] = 1.0
  q[2] = 0.0
  q[3] = 0.0
  q[4] = 100*(1e-5)/params.gamma_1

  x = coords[1]
  y = coords[2]
  if x*x + y*y < dr*dr
    q[4] = 100*3*params.gamma_1*E/( (nu+1)*Float64(pi)*(dr^nu) )
  end

  return nothing
end

"""
  Initial condition for SU2 inviscid channel
"""
function calcInvChannelIC(params::ParamType2,
             coords::AbstractArray{Tmsh,1},
             q::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

  rho = 1.22531
  rhou = rho*170.104
  rhov = 0
  E = 101300/params.gamma_1 + 0.5*rho*(rhou*rhou + rhov*rhov)

  q[1] = rho
  q[2] = rhou
  q[3] = rhov
  q[4] = E

  return nothing
end

"""
  Free stream conditions for SU2 inviscid channel
"""
function calcInvChannelFreeStream(params::ParamType2,
                         coords::AbstractArray{Tmsh,1},
                         q::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

  rho = 1.22531
  rhou = rho*170.104
  rhov = 0
  E = 101100/params.gamma_1 + 0.5*rho*(rhou*rhou + rhov*rhov)

  q[1] = rho
  q[2] = rhou
  q[3] = rhov
  q[4] = E

  return nothing
end

"""
  Used for testing shock capturing scheme
"""
function calcLaplaceSolution(params::ParamType,
                         coords::AbstractArray{Tmsh,1},
                         q::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

  x = coords[1]; y = coords[2]
  val = x*x*x + y*y*y
  val = sin(x)*sin(2*y)
  for i=1:length(q)
    q[i] = val + i
  end

  return nothing
end

global const wedge20_tanbeta = tand(53.422940527228654)
global const tan20 = tand(20)
global const cos20 = cosd(20)
global const sin20 = sind(20)

"""
  Compute the exact solution for a supersonic wedge where the semi-angle is
  20 degrees.  The vertex of the wedge must be at x, y = (-1, 0), and the
  free stream Mach number must be 2.0
"""
function calcWedge20(params::ParamType2,
                     coords::AbstractArray{Tmsh,1},
                     q::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

#  @assert params.Ma == 2.0
#  @assert params.aoa == 0.0

  # figure out if the given point is post-shock or pre-shock

  x = coords[1]; y = coords[2]; tanbeta = wedge20_tanbeta

  if (y > 0 && y > tanbeta*x + tanbeta) || (y <= 0 && y < -tanbeta*x - tanbeta)
    # pre-shock
    calcFreeStream(params, coords, q)
  else
    # computed from the oblique shock relations in Anderson's Aerodynamics
    M2 = 1.2102184008268027
    a2 = 1.1799115865336027

    q[1] = 2.0420057206059625
    q[2] = q[1]*M2*a2*cos20
    q[3] = q[1]*M2*a2*sin20
    q[4] = 7.1584095248438615

    if y < 0
      q[3] = -q[3]
    end
  end

  return nothing
end

"""
  Computes the potential flow solution for flow over a 20 degree wedge.
  Potential flow assumes incompressible flow, so the Mach number should be
  small when running this case
"""
function calcWedge20Potential(params::ParamType2,
                     coords::AbstractArray{Tmsh,1},
                     q::AbstractArray{Tsol,1}) where {Tmsh, Tsol}

  n = 9/8  # this corresponds to a 20 degree semi-angle wedge
  rprime = 2*sqrt(1 + 1.5*1.5)  # define pressure = free stream pressure at this
                                # point
  Vprime = n*params.Ma*rprime^(n-1)
  Vprime2 = Vprime*Vprime

  # Bernoulli's equation constant at rprime.  Use this for *all* streamlines
  c1 = params.p_free/params.rho_free + 0.5*Vprime2

  # the potential flow solution has the wedge geometry as:
  #  \                                    /
  #   \             rather than          /
  #    \                                /
  #     \________                ______/
  # so we transform the coordinates at the beginning of the computation
  # and then transform the velocity at the end.
  # Also, the potential flow solution as the origin at the corner of the wedge,
  # not the lower right corner

  x = -(coords[1] + 1.0); y = coords[2]

  # The potential flow solution si f(z) = z^n, where z = x + i*y, but computing
  # it this way would not be complex-stepable,  Instead apply Eulers identity

  # compute r-theta coordinates
  r = sqrt(x*x + y*y); nrn = n*params.Ma*r^(n-1)
  theta = atan2(y, x); nm1_theta = (n-1)*theta

  # compute velocity
  u = nrn*cos(nm1_theta); v = -nrn*sin(nm1_theta)

  # transform velocity from the potential flow geometry to the real geometry
  v = -v

  # compute pressure using Bernoulli's equation, then energy from it (using
  # ideal gas law)
  v_squared = u*u + v*v
  p = params.rho_free*(c1 - 0.5*v_squared)
  E = p/params.gamma_1 + 0.5*params.rho_free*v_squared

  q[1] = params.rho_free
  q[2] = params.rho_free*u
  q[3] = params.rho_free*v
  q[4] = E

  return nothing
end