text:
## Environment and packages
using Statistics
using LinearAlgebra
using Random
using OrdinaryDiffEq
using DiffEqFlux
using Flux
using Tullio
using RecursiveArrayTools

const t₁ = 5                 # number of simulation years
const ρ = 900                     # Ice density [kg / m^3]
const g = 9.81                    # Gravitational acceleration [m / s^2]
const n = 3                      # Glen's flow law exponent
const maxA = 8e-16
const minA = 3e-17
const maxT = 1
const minT = -25

@views avg(A) = 0.25 .* ( A[1:end-1,1:end-1] .+ A[2:end,1:end-1] .+ A[1:end-1,2:end] .+ A[2:end,2:end] )
@views avg_x(A) = 0.5 .* ( A[1:end-1,:] .+ A[2:end,:] )
@views avg_y(A) = 0.5 .* ( A[:,1:end-1] .+ A[:,2:end] )
@views diff_x(A) = (A[begin + 1:end, :] .- A[1:end - 1, :])
@views diff_y(A) = (A[:, begin + 1:end] .- A[:, 1:end - 1])
@views inn(A) = A[2:end-1,2:end-1]

function loss_iceflow(θ, context, UA)
    H = context[2]
    tspan = (0.0,t₁)

    iceflow_UDE!(dH, H, θ, t) = iceflow_NN!(dH, H, θ, t, context, temp_series[5], UA) # closure
    iceflow_prob = ODEProblem(iceflow_UDE!,H,tspan,θ)
    H_pred = solve(iceflow_prob, BS3(), u0=H, p=θ, reltol=1e-6,
                    sensealg = InterpolatingAdjoint(autojacvec=ZygoteVJP()), save_everystep=false,
                    progress=true, progress_steps = 10)

    # Compute loss function for the full batch
    l_H = sum(H_pred.u[end])

    return l_H
end

function iceflow_NN!(dH, H, θ, t, context, temps, UA)
    # ArrayPartition(B, H, current_year, temp_series, batch_idx)

    year = floor(Int, t) + 1
    if year <= t₁
        temp = temps[year]
    else
        temp = temps[year-1]
    end

    A = predict_A̅(UA, θ, [temp]) # FastChain prediction requires explicit parameters
    # Compute the Shallow Ice Approximation in a staggered grid
    dH .= SIA(dH, H, A, context)
end

# Function without mutation for Zygote, with context as an ArrayPartition
function SIA(dH, H, A, context)
    # Retrieve parameters
    B = context[1]
    # Update glacier surface altimetry
    S = B .+ H

    # All grid variables computed in a staggered grid
    # Compute surface gradients on edges
    dSdx = diff_x(S) / Δx
    dSdy = diff_y(S) / Δy
    ∇S = (avg_y(dSdx).^2 .+ avg_x(dSdy).^2).^((n - 1)/2)

    Γ = 2 * A * (ρ * g)^n / (n+2) # 1 / m^3 s
    D = Γ .* avg(H).^(n + 2) .* ∇S

    # Compute flux components
    dSdx_edges = diff_x(S[:,2:end - 1]) / Δx
    dSdy_edges = diff_y(S[2:end - 1,:]) / Δy
    Fx = .-avg_y(D) .* dSdx_edges
    Fy = .-avg_x(D) .* dSdy_edges

    #  Flux divergence
    @tullio dH[i,j] := -(diff_x(Fx)[pad(i-1,1,1),pad(j-1,1,1)] / Δx + diff_y(Fy)[pad(i-1,1,1),pad(j-1,1,1)] / Δy) # MB to be added here

    return dH
end

predict_A̅(UA, θ, temp) = UA(temp, θ) .* 1e-16

function fake_temp_series(t, means=Array{Float64}([0,-2.0,-3.0,-5.0,-10.0,-12.0,-14.0,-15.0,-20.0]))
    temps, norm_temps, norm_temps_flat = [],[],[]
    for mean in means
       push!(temps, mean .+ rand(t).*1e-1) # static
    end

    return temps
end

##################################################
#### Generate reference dataset ####
##################################################

nx = ny = 100
const B = zeros(Float64, (nx, ny))
const σ = 1000
H₀ = Matrix{Float64}([ 250 * exp( - ( (i - nx/2)^2 + (j - ny/2)^2 ) / σ ) for i in 1:nx, j in 1:ny ])
Δx = Δy = 50 #m
ensemble = EnsembleSerial()
const minA_out = 0.3
const maxA_out = 8
sigmoid_A(x) = minA_out + (maxA_out - minA_out) / ( 1 + exp(-x) )

const temp_series = fake_temp_series(t₁)

# Train UDE
UA = FastChain(
        FastDense(1,3, x->tanh.(x)),
        FastDense(3,10, x->tanh.(x)),
        FastDense(10,3, x->tanh.(x)),
        FastDense(3,1, sigmoid_A)
    )
θ = initial_params(UA)

H = deepcopy(H₀)
current_year = 0f0
# Tuple with all the temp series and H_refs
context = (B, H, current_year, temp_series)
loss(θ) = loss_iceflow(θ, context, UA) # closure

# Debugging
println("Gradients: ", gradient(loss, θ))


text:
"""
This module implements Value Iteration for large/continuous state spaces by solving for some
subspace of the state space and interpolating the value function over the rest of the state space
"""
module LocalApproximationValueIteration

using Random
using Printf
using POMDPs
using POMDPModelTools
using LocalFunctionApproximation
using POMDPLinter: @POMDP_require, @req, @subreq, @warn_requirements

import POMDPs: Solver, solve, Policy, action, value

# Exports related to solver
export
    LocalApproximationValueIterationPolicy,
    LocalApproximationValueIterationSolver

include("local_approximation_vi.jl")

end # module

ext: jl


text:
#This file contains a wrapper object and some helper functions
#relating to endowment files
#creates a table of summary values
function summarizeNCCSFile(data::NCCSFile)::DataFrame

  dataDF::DataFrame = data.df #for convenience since we are just reading, not writing
  summaryDF::DataFrame = metaTable(data.meta) #get the metadata in tabular format
  fields::Vector{Symbol} = summaryDF[:fieldNames]
  #the below is an UGLY hack for CSV compat. Fix it later!
  summaryDF[:longDescriptions] .= (s::String->strip("""γ,γ$s""")).(summaryDF[:longDescriptions])

  NData::Int = size(dataDF,1)
  NFields::Int = length(fields)
  #allocate the new columns
  fileRows::Vector{MInt} = fill(NData, NFields)
  N::Vector{MInt} = missings(NFields)
  pMissing::Vector{MFloat64} = missings(NFields)
  Unique::Vector{MInt} = missings(NFields)
  Mean::Vector{MFloat64} = missings(NFields)
  MeanSq::Vector{MFloat64} = missings(NFields)
  Std::Vector{MFloat64} = missings(NFields)
  Skew::Vector{MFloat64} = missings(NFields)
  Kurtosis::Vector{MFloat64} = missings(NFields)
  Min::Vector{MFloat64} = missings(NFields)
  p05::Vector{MFloat64} = missings(NFields)
  p25::Vector{MFloat64} = missings(NFields)
  p50::Vector{MFloat64} = missings(NFields)
  p75::Vector{MFloat64} = missings(NFields)
  p95::Vector{MFloat64} = missings(NFields)
  Max::Vector{MFloat64} = missings(NFields)

  for i::Int ∈ 1:NFields
    #first get some type info
    T::Type = data.meta.dictionary[fields[i]]
    TSub::Type = typeintersect(T, Union{Float64, Bool, Int, String, CategoricalValue}) #get the nonmissable type
    TSub = CategoricalValue <: TSub || TSub <: CategoricalValue ? Symbol : TSub

    #compute if the type is numeric, then get the column type
    numeric::Bool = (T <: MFloat64) || (T <: MInt)
    VecT = (T <: CategoricalValue) || (CategoricalValue <: T) ? Vector{Union{Symbol,Missing}} : Vector{T}

    #form the rows and do the applicable computations
    col::VecT = VecT(undef, NData)
    col .= VecT(dataDF[fields[i]])
    sansMissing::Vector{TSub} = collect(skipmissing(col))
    if length(sansMissing) == 0
      N[i] = 0.0
      pMissing[i] = 1.0
    elseif numeric
      N[i] = length(sansMissing)
      pMissing[i] = 1.0 - length(sansMissing) / length(col)
      Unique[i] = length(unique(col))
      Mean[i] = mean(sansMissing)
      MeanSq[i] = mean(sansMissing .* sansMissing)
      Std[i] = std(sansMissing)
      Skew[i] = skewness(sansMissing)
      Kurtosis[i] = kurtosis(sansMissing)
      Min[i] = minimum(sansMissing)
      p05[i] = percentile(sansMissing,0.05)
      p25[i] = percentile(sansMissing,0.25)
      p50[i] = percentile(sansMissing,0.50)
      p75[i] = percentile(sansMissing,0.75)
      p95[i] = percentile(sansMissing,0.95)
      Max[i] = maximum(sansMissing)
    else
      N[i] = length(sansMissing)
      pMissing[i] = 1.0 - length(sansMissing) / length(col)
      Unique[i] = length(unique(col))
    end
  end

  #finish building the dataframe
  summaryDF = hcat(summaryDF, DataFrame(fileRows=fileRows, N=N, pMissing=pMissing, Unique=Unique, Mean=Mean,
    MeanSq=MeanSq, Std=Std, Skew=Skew, Kurtosis=Kurtosis, Min=Min, p05=p05, p25=p25,
    p50=p50, p75=p75, p95=p95, Max=Max))

  return summaryDF
end

#convenience method since loading many files is memory intensive
summarizeNCCSFile(meta::MetaDictionary) =
  summarizeNCCSFile(loadNCCSFile(meta, reserializeData=false, returnNCCSFile=true))::DataFrame

#vectorized method for creating summary tables
function summarizeNCCSFiles(metadata::Vector{MetaDictionary}; mode=:serial)::DataFrame
  #get the summary tables

  if mode == :serial #terminal branch, summarize the data
    NData::Int = length(metadata)
    #println("typeof: $(summarizeNCCSFile(metadata[1]))")
    summaryTable::DataFrame = vcat((summarizeNCCSFile).(metadata)...)
    return summaryTable

  elseif mode == :parallel #relaunch the dedicated parallel function
    numFiles::Int = length(metadata)
    pids = workers()
    np::Int = length(pids)

    #divide up the work
    assignments::Vector{Vector{Int}} = ((i::Int)->Vector{Int}()).(1:np)
    for i::Int ∈ 1:numFiles
      push!(assignments[i%np+1], i) #assigns a file to a pid
    end

    #set up a vector for pushing and receiving assignments
    futures::Vector{Future} = Vector{Future}()
    sizehint!(futures, np)

    for i::Int ∈ 1:np
      if length(assignments[i]) > 0
        push!(futures,
          (@spawnat pids[i] summarizeNCCSFiles(metadata[assignments[i]], mode=:serial)))
      end
    end
    outDFs::Vector{DataFrame} = (fetch).(futures)
    outDF::DataFrame = vcat(outDFs...) #concatenate the results
    return outDF
  else
    error("Incorrect mode for function summarizeNCCSFiles")
  end

end

#NOTE this ugly hack works around an issue in CSV.write
#if deleting this, be sure to get the flag where it is inserted elsehwere in the code
function csvWorkaroundWrite(filePath::String, df::DataFrame)

    b = IOBuffer()
    CSV.write(b, df)
    #b = Stream(format"CSV", IOBuffer())
    #save(b, df)
    outString::String = String(take!(copy(b)))
    outString = replace(outString, "\r\n"=>"InnerNewLineρη")
    outString = replace(outString, "\n"=>"\r\n")
    outString = replace(outString, "InnerNewLineρη"=>"\n")
    outString = replace(outString, "γ,γ"=>"")
    write(filePath, outString)
    #end ugly hack

    return nothing
end

#this function constructs and write a summary data file
function writeNCCSFileSummary(metadata::Vector{MetaDictionary};
  documentationPath::String =DOCUMENTATION_PATH,  summaryName = SUMMARY_NAME)::DataFrame

  summaryDF::DataFrame = summarizeNCCSFiles(metadata, mode=:serial)
  csvWorkaroundWrite("$documentationPath\\$summaryName.csv", summaryDF) #workaround on CSV

  #also make a serialized version
  oStream::IOStream = open("$documentationPath\\$summaryName.jls", "w")
  serialize(oStream, summaryDF)
  close(oStream)

  println("Summary table written.")

  return summaryDF
end

#consolidates the summary table
function consolidateSummary(metadata::Vector{MetaDictionary};
  refreshSummary::Bool = true, documentationPath::String = DOCUMENTATION_PATH,
  summaryName::String = SUMMARY_NAME, consolidatedName::String = CONSOLIDATED_NAME)

  #load and/or generate the table
  if !isfile("$documentationPath\\$summaryName.jls") || refreshSummary
    summaryDF::DataFrame = writeNCCSFileSummary(metadata)
  else
    iStream::IOStream = open("$documentationPath\\$summaryName.jls")
    summaryDF = deserialize(iStream)
    close(iStream)
  end

  summaryDF = summaryDF[(!).(summaryDF[:fullFiles]),:]
  NSummary::Int = size(summaryDF,1)

  #copy and do some initial filtering
  uniqueDF::DataFrame = deepcopy(summaryDF)
  deletecols!(uniqueDF, [:Unique, :Skew, :Kurtosis, :p05, :p25, :p50, :p75, :p95])
  uniqueDF[[:fileRows, :N, :pMissing, :Mean, :MeanSq, :Std, :Min, :Max]] =  missing
  NUnique::Int = size(uniqueDF,1)

  #create a smaller version with only unique columns
  uniqueDF[:keep] = trues(NUnique)
  uniqueDFBy::GroupedDataFrame = groupby(uniqueDF, [:fieldNames, :nonprofitTypes])
  for df::AbstractDataFrame in uniqueDFBy
    maxYear::Int = maximum(df[:fileYears])

    df[:keep] = (maxYear .== df[:fileYears])
    if sum(df[:keep]) > 1
      @warn "Multiple rows kept for field $(df[1, :fieldNames]), type $(df[1, :nonprofitTypes]),
        fullFile $(df[1, :fullFiles]), num kept $(sum(df[:keep])) in file $(df[1, :fileNames])"
    end
  end

  #make a dictionary to get the address of unique values

  uniqueDF = uniqueDF[uniqueDF[:keep], :]
  NUnique = size(uniqueDF, 1)
  uniqueIndex::Dict = Dict((uniqueDF[i, :fieldNames], uniqueDF[i, :nonprofitTypes])=>i for i ∈ 1:NUnique)

  uniqueDF[:fileYearsCovered] = Vector{MInt}(missing,NUnique)
  uniqueDF[:fileYearsMin] = Vector{MInt}(missing,NUnique)
  uniqueDF[:fileYearsMax] = Vector{MInt}(missing,NUnique)
  uniqueDF[:fileYearsGap] = Vector{MInt}(missing,NUnique)
  uniqueDF[:UniqueMin] = Vector{MInt}(missing,NUnique)
  uniqueDF[:UniqueMax] = Vector{MInt}(missing,NUnique)
  uniqueDF[:inconsistentType] = Vector{MBool}(missing,NUnique)
  uniqueDF[:inconsistentNumeric] = Vector{MBool}(missing,NUnique)

  #now get the aggregate summary stats
  #first make sure we only select columns for which we have information
  #validSummaryCols::Vector{Bool} =  ((f::Symbol, t::Symbol) -> haskey(uniqueIndex, (f,t))).(
  #  summaryDF[:fieldNames], summaryDF[:nonprofitTypes])
  summaryDFBy = groupby(summaryDF, [:fieldNames, :nonprofitTypes])
  for df::SubDataFrame in summaryDFBy
    #by(summaryDF, [:fieldNames, :nonprofitTypes]) do df::AbstractDataFrame

    row::Int = uniqueIndex[(df[1, :fieldNames], df[1, :nonprofitTypes])]

    #now aggregate the summary info
    uniqueDF[row,:fileYearsCovered] = size(df,1)
    uniqueDF[row,:fileYearsMin] = minimum(skipmissing(df[:fileYears]))
    uniqueDF[row,:fileYearsMax] = maximum(skipmissing(df[:fileYears]))
    uniqueDF[row, :fileYearsGap] =  #count the missing years
      (uniqueDF[row,:fileYearsMax]-uniqueDF[row,:fileYearsMin] + 1) - uniqueDF[row,:fileYearsCovered]

    #check if the metadata is at least somewhat consistent

    #make sure we don't reduce over 0D
    if length(df[(!ismissing).(df[:Unique]),:Unique]) > 0
      uniqueDF[row, :UniqueMin] = maximum(skipmissing(df[:Unique]))
      uniqueDF[row, :UniqueMax] = minimum(skipmissing(df[:Unique]))
    end
    uniqueDF[row,:inconsistentType] = length(unique(df[:nonprofitTypes])) > 1
    uniqueDF[row,:inconsistentNumeric] = length(unique(df[:numeric])) > 1

    #Get information on the quantity of missing values
    uniqueDF[row, :fileRows] = sum(skipmissing(df[:fileRows]))
    uniqueDF[row, :N] = sum(skipmissing(df[:N]))
    uniqueDF[row, :pMissing] =
      (uniqueDF[row, :N] > 0) ? (uniqueDF[row, :fileRows] - uniqueDF[row, :N]
      ) / uniqueDF[row, :fileRows] : 0.0

    #do the numeric processing
    if uniqueDF[row, :numeric]
      uniqueDF[row, :Mean] = sum(skipmissing(df[:Mean] .* df[:N])) / uniqueDF[row, :N]
      uniqueDF[row, :MeanSq] = sum(skipmissing(df[:MeanSq] .* df[:N])) / uniqueDF[row, :N]
      uniqueDF[row, :Std] = (uniqueDF[row, :MeanSq] - uniqueDF[row, :Mean]^2)^0.5
      uniqueDF[row, :Min] = minimum(skipmissing(df[:Min]))
      uniqueDF[row, :Max] = maximum(skipmissing(df[:Max]))
    end
  end

  csvWorkaroundWrite("$documentationPath\\$consolidatedName.csv", uniqueDF)

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
import Base.cat
using DelimitedFiles
using WriteVTK
using SparseArrays
using SymRCM
using MeshCore: AbsShapeDesc, P1, L2, T3, Q4, T4, H8
using MeshCore: VecAttrib
using MeshCore: shapedesc, nshapes,  IncRel
using MeshCore: ShapeColl
using LinearAlgebra: norm

# using UnicodePlots # 4 debugging only

"""
    transform(ir, T = x -> x)

Change the locations of the vertices through the transformation `T`.
"""
function transform(ir, T = x -> x)
    locs = ir.right.attributes["geom"]
    N = length(locs[1])
    TC = eltype(locs[1])
    nlocs =  VecAttrib([SVector{N, TC}(T(locs[i])) for i in 1:length(locs)])
    ir.right.attributes["geom"] = nlocs
end

"""
    vconnected(ir)

Find whether or not the vertices are connected to any shape on the left of the
incidence relation.

- `isconnected` = vector is returned which is for the node `k` either true
  (vertex `k` is connected), or false (vertex `k` is not connected).
"""
function vconnected(ir)
    isconnected = falses(nshapes(ir.right));
    for i in 1:nrelations(ir)
        for j in 1:nentities(ir, i)
            isconnected[ir[i, j]] = true
        end
    end
    return isconnected
end

"""
    vnewnumbering(ir, isconnected)

Compute the new numbering obtained by deleting unconnected vertices.
"""
function vnewnumbering(ir, isconnected)
    @_check length(isconnected) == nshapes(ir.right)
    p = fill(zero(indextype(ir)), nshapes(ir.right));
    id = 1;
    for i in 1:length(isconnected)
        if (isconnected[i])
            p[i] = id;
            id = id+1;
        end
    end
    return p
end

"""
    compactify(ir, new_numbering)

Compact the finite element node set by deleting unconnected nodes.

`fens` = array of finite element nodes
`connected` = The array element `connected[j]` is either 0 (when `j` is an
  unconnected node), or a positive number (when node `j` is connected to
  other nodes by at least one finite element)

# Output
`fens` = new set of finite element nodes
`new_numbering`= array which tells where in the new `fens` array the
     connected nodes are (or 0 when the node was unconnected). For instance,
     node 5 was connected, and in the new array it is the third node: then
     `new_numbering[5]` is 3.

# Examples
Let us say there are nodes not connected to any finite element that you
would like to remove from the mesh: here is how that would be
accomplished.
```
connected = findunconnnodes(fens, fes);
fens, new_numbering = compactnodes(fens, connected);
fes = renumberconn!(fes, new_numbering);
```
Finally, check that the mesh is valid:
```
validate_mesh(fens, fes);
```
"""
function compactify(ir, new_numbering)
    @_check length(new_numbering) == nshapes(ir.right)
    locs = ir.right.attributes["geom"] # this is the old geometry attribute
    N = length(locs[1]) # number of spatial dimensions
    TC = eltype(locs[1])
    # Now generate a vector of the locations for connected vertices
    v = [SVector{N, TC}(locs[i]) for i in 1:length(locs) if new_numbering[i] != 0]
    locs =  VecAttrib(v) # this is the new geometry attribute
    vertices = ShapeColl(P1, length(locs), "vertices")
    vertices.attributes["geom"] = locs
    # Elements
    C = [SVector{3}(ir[idx]) for idx in 1:nrelations(ir)]
    elements = ShapeColl(shapedesc(ir.left), length(C), "elements")
    return IncRel(elements, vertices, C)
end

# This code was taken over from FinEtools: refer to fusenodes
function _fusen(xyz1, id1, xyz2, id2, tolerance)
    dim = size(xyz1,2);
    nn1 = size(xyz1, 1)
    nn2 = size(xyz2, 1)
    # Decide which nodes should be checked for proximity
    ib = intersectboxes(inflatebox!(boundingbox(xyz1), tolerance), inflatebox!(boundingbox(xyz2), tolerance))
    node1in = fill(false, nn1);
    node2in = fill(false, nn2);
    if length(ib) > 0
        for i=1:nn1
            node1in[i] = inbox(ib, @view xyz1[i, :])
        end
        for i=1:nn2
            node2in[i] = inbox(ib, @view xyz2[i, :])
        end
    end
    # Mark nodes from the first array that are duplicated in the second
    if (tolerance > 0.0) # should we attempt to merge nodes?
        for i=1:nn1
            if node1in[i]
                breakoff = false
                for rx=1:nn2
                    if node2in[rx]
                        distance = 0.0
                        for cx=1:dim
                            distance = distance + abs(xyz2[rx,cx]-xyz1[i,cx]);
                            if (distance >= tolerance) # shortcut: if the distance is already too large, stop checking
                                break
                            end
                        end
                        if (distance < tolerance)
                            id1[i] = -rx; breakoff = true;
                        end
                    end
                    if breakoff
                        break
                    end
                end
            end
        end
    end
    # Generate  fused arrays of the nodes. First copy in the nodes from the second set...
    xyzm = zeros(eltype(xyz1),nn1+nn2,dim);
    for rx = 1:nn2
        for cx = 1:dim
            xyzm[rx,cx] = xyz2[rx,cx];
        end
    end
    # idm = zeros(eltype(id1),nn1+nn2);
    # for rx = 1:nn2
    #     idm[rx] = rx;
    # end
    mid=nn2+1;
    # ...and then we add in only non-duplicated nodes from the first set
    for i=1:nn1
        if id1[i]>0
            id1[i] = mid;
            # idm[mid] = mid;
            for cx = 1:dim
                xyzm[mid,cx] = xyz1[i,cx];
            end
            mid = mid+1;
        else
            id1[i] = id2[-id1[i]];
        end
    end
    nnodes = mid-1;
    # The set 1 is described by these locations. The new numbering applies also
    # to set 1.
    xyzm = xyzm[1:nnodes,:];
    new_indexes_of_set1_nodes = deepcopy(id1);
    # The node set 2 numbering stays the same
    return xyzm, new_indexes_of_set1_nodes
end

"""
    fusevertices(locs1, locs2, tolerance)

Fuse together vertices from two vortex sets.

Fuse two vertex sets. If necessary, by gluing together vertices located within
`tolerance` of each other. The two vertex sets are fused
together by merging the vertices that fall within a box of size `tolerance`. The
merged vertex set, `fens`, and the new  indexes of the nodes in the set `fens1`
are returned.

The set `fens2` will be included unchanged, in the same order,
in the node set `fens`.
The indexes of the node set `fens1` will have changed.

# Example
After the call to this function we have
`k=new_indexes_of_fens1_nodes[j]` is the node in the node set `fens` which
used to be node `j` in node set `fens1`.
The finite element set connectivity that used to refer to `fens1`
needs to be updated to refer to the same nodes in  the set `fens` as
     `updateconn!(fes, new_indexes_of_fens1_nodes);`
"""
function fusevertices(locs1, locs2, tolerance)
    @_check length(locs1[1]) == length(locs2[1])
    dim = length(locs1[1]);
    nn1 = length(locs1)
    nn2 = length(locs2)
    TC = eltype(locs1[1])
    xyz1 = zeros(TC,nn1,dim);
    for i in 1:length(locs1)
        xyz1[i, :] = locs1[i]
    end
    xyz2 = zeros(TC,nn2,dim);
    for i in 1:length(locs2)
        xyz2[i, :] = locs2[i]
    end
    id1 = collect(1:nn1);
    id2 = collect(1:nn2);
    xyzm, new_indexes_of_set1_nodes = _fusen(xyz1, id1, xyz2, id2, tolerance)
    N = size(xyzm, 2)
    nlocs1 = VecAttrib([SVector{N, TC}(xyzm[i, :]) for i in 1:size(xyzm, 1)])
    return nlocs1, new_indexes_of_set1_nodes
end

"""
    withvertices(ir, locs)

Create a new incidence relation referring to a different set of vertices.

Presumably the set of vertices is simply broadened from the one used by the
incidence relation `ir`.
"""
function withvertices(ir, locs)
    vertices = ShapeColl(P1, length(locs), "vertices")
    vertices.attributes["geom"] = locs
    elements = ShapeColl(shapedesc(ir.left), nrelations(ir), "elements")
    return IncRel(elements, vertices, [ir[idx] for idx in 1:nrelations(ir)])
end

"""
    renumbered(ir, p)

Renumber the connectivity of the shapes based on a new numbering for the
vertices.

- `p` = new serial numbers for the vertices.  The connectivity
          should be changed as `conn[j]` --> `p[conn[j]]`

Returns new incidence relation with renumbered connectivity.

# Example

To do: Revise this example.
Let us say there are nodes not connected to any finite element that you would
like to remove from the mesh: here is how that would be accomplished.
```
connected = findunconnnodes(fens, fes);
fens, p = compactfens(fens, connected);
fes = renumberconn!(fes, p);
```
Finally, check that the mesh is valid:
```julia
validate_mesh(fens, fes);
```
"""
function renumbered(ir, p)
    N = nentities(ir, 1)
    C = SVector{N, indextype(ir)}[]
    for i in 1:nrelations(ir)
        c = ir[i]
        push!(C, SVector{N}(p[c]))
    end
    locs = ir.right.attributes["geom"]
    vertices = ShapeColl(P1, length(locs), "vertices")
    vertices.attributes["geom"] = locs
    elements = ShapeColl(shapedesc(ir.left), length(C), "elements")
    return IncRel(elements, vertices, C)
end

"""
    cat(ir1::T, ir2::T) where {T<:IncRel}

Concatenate the connectivities of two shape collections.

The two shape collections must be of the same shape descriptor, and they must
refer to the same attribute describing the locations of the vertices.
"""
function cat(ir1::T, ir2::T) where {T<:IncRel}
    # First, the left shapes must be of the same description
    @_check shapedesc(ir1.left) == shapedesc(ir2.left)
    # The shape collections on the right must be vertices
    @_check shapedesc(ir1.right) == shapedesc(ir2.right) == P1
    # The shape collections must refer to the same set of vertices. The length
    # is a cheap proxy for checking that condition.
    @_check length(ir1.right.attributes["geom"]) == length(ir2.right.attributes["geom"])
    # The shape collections must have fixed cardinality.
    N = nentities(ir1, 1)
    C = SVector{N, indextype(ir1)}[]
    for i in 1:nrelations(ir1)
        push!(C, ir1[i])
    end
    for i in 1:nrelations(ir2)
        push!(C, ir2[i])
    end
    locs = ir1.right.attributes["geom"]
    vertices = ShapeColl(P1, length(locs), "vertices")
    vertices.attributes["geom"] = locs
    elements = ShapeColl(shapedesc(ir1.left), length(C), "elements")
    return IncRel(elements, vertices, C)
end

"""
    mergeirs(conn1, conn2, tolerance = eps())

Merge two incidence relations, fuse vertices closer then tolerance.
"""
function mergeirs(conn1, conn2, tolerance = eps())
    locs1 = conn1.right.attributes["geom"]
    locs2 = conn2.right.attributes["geom"]
    nlocs1, ni1 = fusevertices(locs1, locs2, tolerance)
    conn1 = withvertices(conn1, nlocs1)
    conn2 = withvertices(conn2, nlocs1)
    conn1 = renumbered(conn1, ni1)
    return cat(conn1, conn2)
end

"""
    minimize_profile(conn)

Re-number the vertices so that the profile can be minimized.

# Output

Renumbered incidence relation.
"""
function minimize_profile(conn)
    I = fill(zero(Int), 0)
    J = fill(zero(Int), 0)
    sizehint!(I, nshapes(conn.right))
    for k in 1:nrelations(conn)
        ne = nentities(conn, k)
        for i in 1:ne
            append!(I, conn[k])
            for m in 1:ne
                push!(J, conn[k, i])
            end
        end
    end
    V = fill(1.0, length(I))
    S = sparse(I, J, V, nshapes(conn.right), nshapes(conn.right))
    # display(spy(S))
    # find the new numbering (permutation)
    p = symrcm(S)
    # display(spy(S[p, p]))
    # number the vertices of the shapes on the left using the new permutation
    conn = renumbered(conn, p)
    # reorder the vertices attribute: the order of the vertices changed
    ip = similar(p) # inverse permutation
    ip[p] = 1:length(p)
    locs = conn.right.attributes["geom"]
    newlocs = VecAttrib([locs[k] for k in ip])
    conn.right.attributes["geom"] = newlocs
    return conn
end


text: println("RaawwWWWWWRRRR!!")

max_stars_repo_name: bneron/apps

max_issues_repo_name: bneron/apps


text:
using DemoActions
using Test

@testset "DemoActions.jl" begin
    @test 2*5 == 10
end

# """
# The `mysolvers_inplace` test set verifies that the higher order solvers
# have lower error when comparing against the reference DifferentialEquations.jl solve.

# The following algorithms will be tested:

#     * eulers_inplace_f!
#     * midpoint_inplace_f!
#     * rk4_inplace_f!
#     * dp5! - this should error currently

# On the following models (found in src/models.jl):
#     * lotka!

#todo test that the inplace vs allocing return same answers
# """
@testset "mysolvers" begin

    test_algs = [eulers_inplace_f!, midpoint_inplace_f!, rk4_inplace_f!, dp5_inplace!]
    test_models = [lotka!]
    u0 = [1., 1.]
    p = [1.5, 1.0, 3.0, 1.0]
    tspan = (0., 10.)
    dt = 0.25

    for alg in test_algs
        @testset "$(alg)_inplace" begin
            for f in test_models
                @testset "$f" begin
                    sol = mysolve_inplace(alg, f, u0, tspan, dt; p=p)
                    @test !isnothing(sol)
                end
            end
        end
    end
end


text:
info("Testing ", MOD.AdditiveKernel)
for kernel_obj in additive_kernels
    info("Testing ", kernel_obj)

    # Test constructors
    for T in FloatingPointTypes
        default_floats, default_others = all_default_args[kernel_obj]
        default_args = (T[default_floats...]..., default_others...)
        fields = fieldnames(kernel_obj)
        k = length(fields) == 0 ? (kernel_obj){T}() : (kernel_obj)(default_args...)

        @test eltype(k) == T

        if length(fields) == 0 && T == Float64
            @test eltype((kernel_obj)()) == T
        end

        for i in eachindex(fields)
            @test getfield(k, fields[i]).value == default_args[i]
        end

        for (x,y) in ((zero(T),zero(T)), (zero(T),one(T)), (one(T),zero(T)), (one(T),one(T)))
            f = all_phifunctions[kernel_obj]
            @test_approx_eq MOD.phi(k, x, y) f(default_args..., x, y)
        end
    end

    # Test conversions
    for T in FloatingPointTypes
        k = (kernel_obj)()
        for U in FloatingPointTypes
            @test U == eltype(convert(Kernel{U}, k))
        end
    end

    # Test Properties
    properties = all_kernelproperties[kernel_obj]
    k = (kernel_obj)()

    @test MOD.ismercer(k)        == properties[1]
    @test MOD.isnegdef(k)        == properties[2]
    @test MOD.attainsnegative(k) == properties[3]
    @test MOD.attainszero(k)     == properties[4]
    @test MOD.attainspositive(k) == properties[5]

    # Test that output does not create error
    show(DevNull, k)

end


text:
export
    Path,
    Encounter,
    update!

#=
This file defines the encounter type, which is used to
track information about an encounter in a simulation.
All information that should be tracked, but that is not
necessary for the uavs to make decisions or for the dynamics
to step the system forward should be stored in an encounter.
=#

#=
This type is used to track the trajectory of a uav.
=#
type Path
    x_coordinates::Vector{Float64}
    y_coordinates::Vector{Float64}
end

#=
Encounter type tracks information about an encounter.
=#
type Encounter
    paths::Vector{Path}
    actions::Matrix{Float64}
    signals::Array{Symbol}
    min_horizontal_dist::Float64
    utilities::Array{Array{Float64}}
    belief_states::Array{SparseMatrixCSC}
    states::Array{State}
    min_horizontal_dists::Array{Float64}
    nmac::Bool
    function Encounter()
        paths = [Path([], []), Path([], [])]
        actions = Array(Float64, 2, 0)
        signals = [:no_signal, :no_signal]
        return new(paths, actions, signals, Inf, [], [], [], [], false)
    end

end

#=
Description:
Updates the encounter object to include the actions

Parameters:
- encounter: Encounter type to add actions to
- actions: actions to add

Side Effects:
- updates encounter to include the actions
=#
function update!(encounter::Encounter, actions::Vector{Float64})
    encounter.actions = hcat(encounter.actions, actions)
end

#=
Description:
Updates the encounter object to include the actions

Parameters:
- encounter: Encounter type to add actions to
- actions: actions to add
- utilities: the utilities for both uavs
- belief_state: the belief_state for both uavs

Side Effects:
- updates encounter to include the actions, beliefs, utilities
=#
function update!(encounter::Encounter, uav_actions::UAVActions,
            utilities::Array{Float64}, belief_state::SparseMatrixCSC)
    encounter.actions = hcat(encounter.actions, uav_actions.actions)
    encounter.signals = hcat(encounter.signals, uav_actions.signals)
    push!(encounter.utilities, utilities)
    push!(encounter.belief_states, belief_state)
end

#=
Description:
Updates the encounter object to include the uav x and y coordinates
from the state.

Parameters:
- encounter: Encounter type to add actions to
- state: state containing x and y coordinates to add
- min_horizontal_dist: minimum horz dist to add to list

Side Effects:
- updates encounter to include the x and y coordinates and horz dist
- also to include state, min_horizontal_dist
=#
function update!(encounter::Encounter, state::State, min_horizontal_dist::Float64 = -1.)
    # update trajectory information
    push!(encounter.states, state)
    for (uidx, uav_state) in enumerate(state.uav_states)
        push!(encounter.paths[uidx].x_coordinates, uav_state.x)
        push!(encounter.paths[uidx].y_coordinates, uav_state.y)
    end

    # update horizontal distance array
    if min_horizontal_dist == -1.
        min_horizontal_dist = get_horizontal_distance(state)
    end
    push!(encounter.min_horizontal_dists, min_horizontal_dist)
end

function save(encounters::Vector{Encounter}, output_filepath::String)
    JLD.save(output_filepath, "encounters", encounters)
end


text:
using GLMakie
function peaks(; n = 49)
    x = LinRange(-3, 3, n)
    y = LinRange(-3, 3, n)
    a = 3 * (1 .- x').^2 .* exp.(-(x'.^2) .- (y .+ 1).^2)
    b = 10 * (x' / 5 .- x'.^3 .- y.^5) .* exp.(-x'.^2 .- y.^2)
    c = 1 / 3 * exp.(-(x' .+ 1).^2 .- y.^2)
    return (x, y, a .- b .- c)
end
x, y, z = peaks(; n=30)
set_theme!(theme_dark())
fig =  Figure(resolution = (1200,800))
axs = [Axis3(fig[1,i]; aspect = :data) for i in 1:3]
contour3d!(axs[1], x, y, z; levels = 20, transparency = true)
contour!(axs[1], x, y, z; levels = 20, transformation = (:xy, minimum(z)),
    transparency = true)
lines!(axs[2], cat(x, NaN, dims=1), y, vcat(z, fill(NaN,30)'),
    color = maximum(z, dims=1)[1,:], transparency = false)
lines!(axs[2], cat(x, NaN, dims=1), y,
  	vcat(z, fill(NaN,30)')*0 .+ minimum(z), color = maximum(z, dims=1)[1,:],
  	transparency = true)
wireframe!(axs[3], x, y, z; color = :grey90, transparency = true)
wireframe!(axs[3], x, y, z*0 .+ minimum(z); color = :grey90,
	transparency = true)
[hidedecorations!(ax; grid = false) for ax in axs]
save("topoLines.png", fig)
fig


text:
export getSensMatVecEik, getSensTMatVecEik

function getOp(opCode::Int8,dim::Int64)
if dim==1
	op = mod(opCode,5)-2;
elseif dim==2
	op = mod(opCode,25);
	op = div(op,5)-2;
elseif dim==3
	op = div(opCode,25)-2;
end
end

function getOpCode(op1,op2,op3)
opCode = (op1+2) + 5*(op2+2) + 25*(op3+2);
return opCode;
end

function getSensMatVecEik(b::Array{Float64,1},res::Array{Float64,1},pEik::EikonalParam)
# This is basically a lower triangular solution.
	order = pEik.ordering;
	n = pEik.Mesh.n.+1;
	if pEik.Mesh.dim==2
		rowLength = 5;
		if pEik.HO
			getRowOfSensitivity = getRowOfSensitivity2D;
		else
			getRowOfSensitivity = getRowOfSensitivity2DFO;
		end
		srcCS = loc2cs(pEik.src[1],pEik.src[2],n);
		offset = [1;n[1]];
	else
		rowLength = 7;
		getRowOfSensitivity = getRowOfSensitivity3D;
		srcCS = loc2cs3D(pEik.src[1],pEik.src[2],pEik.src[3],n);
		offset = [1;n[1];n[1]*n[2]];
	end
	G0 = zeros(pEik.Mesh.dim);
	Idxs = zeros(Int64,pEik.Mesh.dim);
	colIdxs = zeros(Int64,rowLength);
	vals 	= zeros(Float64,rowLength);
	############ Invert a lower triangular Matrix ##########################
	res[:] .= 0.0;
	for j=order
		r = b[j];
		rowLength_j = getRowOfSensitivity(pEik,j,vals,colIdxs,srcCS,offset,G0,Idxs,n);
		for i = 1:rowLength_j
			r -= vals[i]*res[colIdxs[i]];
		end
		res[j] = r/vals[1];
	end
	if pEik.Mesh.dim==2
		selfMultiplyWithAnalyticSolution2D(n, pEik.Mesh.h,pEik.src,res);
	else
		selfMultiplyWithAnalyticSolution3D(n, pEik.Mesh.h,pEik.src,res);
	end

	return res;
end

function getSensTMatVecEik(b::Array{Float64,1},res::Array{Float64,1},pEik::EikonalParam)
# This is basically an upper triangular solution. This function ruins b.
	order = pEik.ordering;
	n = pEik.Mesh.n.+1;
	if pEik.Mesh.dim==2
		rowLength = 5;
		if pEik.HO
			getRowOfSensitivity = getRowOfSensitivity2D;
		else
			getRowOfSensitivity = getRowOfSensitivity2DFO;
		end
		srcCS = loc2cs(pEik.src[1],pEik.src[2],n);
		offset = [1;n[1]];
	else
		rowLength = 7;
		getRowOfSensitivity = getRowOfSensitivity3D;
		srcCS = loc2cs3D(pEik.src[1],pEik.src[2],pEik.src[3],n);
		offset = [1;n[1];n[1]*n[2]];
	end
	G0 = zeros(pEik.Mesh.dim);
	Idxs = zeros(Int64,pEik.Mesh.dim);
	colIdxs = zeros(Int64,rowLength);
	vals 	= zeros(Float64,rowLength);

	############ Invert an upper triangular Matrix ##########################
	N = length(b)
	res[:] = b[:]; # up to index i (backwards, starting from n) res holds the result. From that to 1 it holds the residual.
	if pEik.Mesh.dim==2
		selfMultiplyWithAnalyticSolution2D(n, pEik.Mesh.h,pEik.src,res);
	else
		selfMultiplyWithAnalyticSolution3D(n, pEik.Mesh.h,pEik.src,res);
	end

	for k=0:(N-1)
		j = order[N-k];
		rowLength_j = getRowOfSensitivity(pEik,j,vals,colIdxs,srcCS,offset,G0,Idxs,n);
		d = res[colIdxs[1]]/vals[1];
		res[colIdxs[1]] = d;
		for i = 2:rowLength_j
			res[colIdxs[i]] -= vals[i]*d;
		end
	end
	return res;
end

function getRowOfSensitivity3D(pEik::EikonalParam,rowIdx::EikIdxType,vals::Array{Float64},colIdxs::Array{Int64},srcCS::Int64,offsetVec::Array{Int64,1},
								G0::Array{Float64,1},Idxs::Array{Int64,1},n::Array{Int64,1});
# the ordering is : [diag,1st cell dx,2nd cell dx,1st cell dy,2nd cell dy]
# in 2D the length of the rows should be of length 5 at most.

OP = pEik.OP;

h = pEik.Mesh.h;
src = pEik.src;
rowIdx = convert(Int64,rowIdx);
colIdxs[1] = rowIdx;
T1 = pEik.T1;
if srcCS==rowIdx
	vals[1] = 2.0*T1[rowIdx];
	return 1;
end

vals[1] = 0.0;
cs2loc3D(Idxs,rowIdx,n);
@inbounds T0loc = analyticLocal3D(h,src,Idxs[1],Idxs[2],Idxs[3],G0);

k = 1;

for dim = 1:3
	hinvij  = T0loc./h[dim];
	hinv2ij = 0.5*hinvij;
	op = getOp(OP[rowIdx],dim);
	offset = offsetVec[dim];
	if op == 1 # forward 1
		diagDhat  = -hinvij + G0[dim];
		diagSensValue = 2.0*(diagDhat*T1[rowIdx] + hinvij*T1[rowIdx+offset]);
		vals[1]                 += diagDhat*diagSensValue;
		k+=1;
		vals[k]     = hinvij*diagSensValue;
		colIdxs[k]  = rowIdx + offset;
	elseif op == -1 # backward 1
		diagDhat  = hinvij + G0[dim];
		diagSensValue = 2.0*(diagDhat*T1[rowIdx] - hinvij*T1[rowIdx - offset]);
		vals[1]                 += diagDhat*diagSensValue;
		k+=1;
		vals[k]     = -hinvij*diagSensValue;
		colIdxs[k]  = rowIdx - offset;
	elseif op == -2 # backward 2
		diagDhat  = hinv2ij*3.0 + G0[dim];
		diagSensValue = 2.0*(diagDhat*T1[rowIdx] + hinv2ij*(-4.0*T1[rowIdx-offset] + T1[rowIdx-2*offset]));
		vals[1]                 += diagDhat*diagSensValue;
		k+=1;
		vals[k]     = -hinv2ij*4.0*diagSensValue;
		colIdxs[k]  = rowIdx - offset;
		k+=1;
		vals[k]    = hinv2ij*diagSensValue;
		colIdxs[k] = rowIdx - 2*offset;
	elseif op == 2 # forward 2
		diagDhat  = -hinv2ij*3.0 + G0[dim];
		diagSensValue = 2.0*(diagDhat*T1[rowIdx] + hinv2ij*(4.0*T1[rowIdx+offset] - T1[rowIdx+2*offset]));
		vals[1]                 += diagDhat*diagSensValue;
		k+=1;
		vals[k]     = hinv2ij*4.0*diagSensValue;
		colIdxs[k]  = rowIdx + offset;
		k+=1;
		vals[k]    = -hinv2ij*diagSensValue;
		colIdxs[k] = rowIdx + 2*offset;
	end
end
return k;
end

function getRowOfSensitivity2D(pEik::EikonalParam,rowIdx::EikIdxType,vals::Array{Float64},colIdxs::Array{Int64},srcCS::Int64,offsetVec::Array{Int64,1},
								G0::Array{Float64,1},Idxs::Array{Int64,1},n::Array{Int64,1});
# the ordering is : [diag,1st cell dx,2nd cell dx,1st cell dy,2nd cell dy]
# in 2D the length of the rows should be of length 5 at most.

OP = pEik.OP;

h = pEik.Mesh.h;
src = pEik.src;

rowIdx = convert(Int64,rowIdx);
colIdxs[1] = rowIdx;
T1 = pEik.T1;
if srcCS==rowIdx
	vals[1]    = 2.0*T1[rowIdx];
	return 1;
end
vals[1] = 0.0;

cs2loc(Idxs,rowIdx,n);
T0loc = analyticLocal(h,src,Idxs[1],Idxs[2],G0);

k = 1;

for dim = 1:2
	hinvij  = T0loc./h[dim];
	hinv2ij = 0.5*hinvij;
	op = getOp(OP[rowIdx],dim);
	offset = offsetVec[dim];
	if op == 1 # forward 1
		diagDhat  = -hinvij + G0[dim];
		diagSensValue = 2.0*(diagDhat*T1[rowIdx] + hinvij*T1[rowIdx+offset]);
		vals[1]                 += diagDhat*diagSensValue;
		k+=1;
		vals[k]     = hinvij*diagSensValue;
		colIdxs[k]  = rowIdx + offset;
	elseif op == -1 # backward 1
		diagDhat  = hinvij + G0[dim];
		diagSensValue = 2.0*(diagDhat*T1[rowIdx] - hinvij*T1[rowIdx - offset]);
		vals[1]                 += diagDhat*diagSensValue;
		k+=1;
		vals[k]     = -hinvij*diagSensValue;
		colIdxs[k]  = rowIdx - offset;
	elseif op == -2 # backward 2
		diagDhat  = hinv2ij*3.0 + G0[dim];
		diagSensValue = 2.0*(diagDhat*T1[rowIdx] + hinv2ij*(-4.0*T1[rowIdx-offset] + T1[rowIdx-2*offset]));
		vals[1]                 += diagDhat*diagSensValue;
		k+=1;
		vals[k]     = -hinv2ij*4.0*diagSensValue;
		colIdxs[k]  = rowIdx - offset;
		k+=1;
		vals[k]    = hinv2ij*diagSensValue;
		colIdxs[k] = rowIdx - 2*offset;
	elseif op == 2 # forward 2
		diagDhat  = -hinv2ij*3.0 + G0[dim];
		diagSensValue = 2.0*(diagDhat*T1[rowIdx] + hinv2ij*(4.0*T1[rowIdx+offset] - T1[rowIdx+2*offset]));
		vals[1]                 += diagDhat*diagSensValue;
		k+=1;
		vals[k]     = hinv2ij*4.0*diagSensValue;
		colIdxs[k]  = rowIdx + offset;
		k+=1;
		vals[k]    = -hinv2ij*diagSensValue;
		colIdxs[k] = rowIdx + 2*offset;
	end
end
return k;
end

function getRowOfSensitivity2DFO(pEik::EikonalParam,rowIdx::EikIdxType,vals::Array{Float64},colIdxs::Array{Int64},srcCS::Int64,offsetVec::Array{Int64,1},
								G0::Array{Float64,1},Idxs::Array{Int64,1},n::Array{Int64,1});
# the ordering is : [diag,1st cell dx,2nd cell dx,1st cell dy,2nd cell dy]
# in 2D the length of the rows should be of length 5 at most.

OP  = pEik.OP;
h   = pEik.Mesh.h;
src = pEik.src;
T1  = pEik.T1;
rowIdx = convert(Int64,rowIdx);

@inbounds colIdxs[1] = rowIdx;

if srcCS==rowIdx
	@inbounds vals[1]    = 2.0*T1[rowIdx];
	return 1;
end
@inbounds vals[1] = 0.0;

cs2loc(Idxs,rowIdx,n);

@inbounds T0loc = analyticLocal(h,src,Idxs[1],Idxs[2],G0);

k = 1;

for dim = 1:2
	@inbounds op = getOp(OP[rowIdx],dim);
	if op != 0
		@inbounds firstNeighbor = rowIdx + op*offsetVec[dim];
		@inbounds hinvij  = op*(T0loc./h[dim]);

		@inbounds diagDhat  = G0[dim]-hinvij;
		@inbounds diagSensValue = 2.0*(diagDhat*T1[rowIdx] + hinvij*T1[firstNeighbor]);
		@inbounds vals[1]     += diagDhat*diagSensValue;
		k+=1;
		@inbounds vals[k]     = hinvij*diagSensValue;
		@inbounds colIdxs[k]  = firstNeighbor;
	end
end
return k;
end


text:
using NBody

function main()
    while !eof(STDIN)
        str = readline(STDIN)
        nb = NBody.read_nbody_json(str)

        dump_nb(nb)
    end
end

function dump_nb(nb)
    for b in nb.bodies
        map(x-> @printf("%24.16f",x), b.pos)
        println()
    end
end

main()


text:
# ## Benchmark

# This tutorial follows closely the tutorial [OptimizationProblems.jl](https://juliasmoothoptimizers.github.io/OptimizationProblems.jl/dev/benchmark/).
# We show how to use [SolverBenchmark.jl](https://juliasmoothoptimizers.github.io/SolverBenchmark.jl) to benchmark solvers that take a NLPStopping as input.

using LinearAlgebra, NLPModels, Stopping
using JSOSolvers, StoppingInterface
using DataFrames, Printf, SolverBenchmark
using ADNLPModels, OptimizationProblems

# We select the problems from [OptimizationProblems.jl](https://juliasmoothoptimizers.github.io/OptimizationProblems.jl) that are unconstrained and scalable.

df = OptimizationProblems.meta
names_pb_vars = df[(df.variable_nvar .== true) .& (df.ncon .== 0), :name]

ad_problems = (
  OptimizationProblems.ADNLPProblems.eval(Symbol(problem))(n = 31) for problem ∈ names_pb_vars
)

# Then, we prepare the solvers we will benchmark.
# Here, we use [JSOSolvers.jl](https://github.com/JuliaSmoothOptimizers/JSOSolvers.jl) that are made Stopping-compatible using [StoppingInterface.jl](https://github.com/SolverStoppingJulia/StoppingInterface.jl).

solvers = Dict(
  :lbfgs => model -> stopping_to_stats(StoppingInterface.lbfgs(NLPStopping(model); mem=5, atol=1e-5, rtol=0.0, max_time = 5.)),
  :trunk => model -> stopping_to_stats(StoppingInterface.trunk(NLPStopping(model); atol=1e-5, rtol=0.0, max_time = 5.)),
)

# The main function used from [SolverBenchmark.jl](https://juliasmoothoptimizers.github.io/SolverBenchmark.jl).

stats = bmark_solvers(
  solvers, ad_problems
)

# The output of `bmark_solvers` can then be analyzed for the results as a table

cols = [:id, :name, :nvar, :objective, :dual_feas, :neval_obj, :neval_grad, :neval_hess, :iter, :elapsed_time, :status]
header = Dict(
  :nvar => "n",
  :objective => "f(x)",
  :dual_feas => "‖∇f(x)‖",
  :neval_obj => "# f",
  :neval_grad => "# ∇f",
  :neval_hess => "# ∇²f",
  :elapsed_time => "t",
)

for solver ∈ keys(solvers)
  pretty_stats(stats[solver][!, cols], hdr_override=header)
end

first_order(df) = df.status .== :first_order
unbounded(df) = df.status .== :unbounded
solved(df) = first_order(df) .| unbounded(df)
costnames = ["time", "obj + grad + hess"]
costs = [
  df -> .!solved(df) .* Inf .+ df.elapsed_time,
  df -> .!solved(df) .* Inf .+ df.neval_obj .+ df.neval_grad .+ df.neval_hess,
]

# or as a performance profile

using Plots
gr()

profile_solvers(stats, costs, costnames)


text:
#=
Settings for results in the paper "QuantumCircuitOpt: An Open-source Frameworkfor Provably Optimal
Quantum Circuit Design"

Version of QCOpt: v0.3.0

Last updated: Oct 12, 2021
=#

function decompose_controlled_Z()

    println(">>>>> Controlled-Z Gate <<<<<")

    return Dict{String, Any}(
    "num_qubits" => 2,
    "maximum_depth" => 4,
    "elementary_gates" => ["U3_1", "U3_2", "CNot_1_2", "Identity"],
    "target_gate" => QCOpt.CZGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact",

    "U3_θ_discretization" => -π:π/2:π,
    "U3_ϕ_discretization" => -π:π/2:π,
    "U3_λ_discretization" => -π:π/2:π,)
end

function decompose_controlled_V()

    println(">>>>> Controlled-V Gate <<<<<")

    return Dict{String, Any}(
    "num_qubits" => 2,
    "maximum_depth" => 7,
    "elementary_gates" => ["H_1", "H_2", "T_1", "T_2", "Tdagger_1", "Tdagger_2", "CNot_1_2", "CNot_2_1", "Identity"],
    "target_gate" => QCOpt.CVGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact")
end

function decompose_controlled_H()

    println(">>>>> Controlled-H Gate <<<<<")

    return Dict{String, Any}(
    "num_qubits" => 2,
    "maximum_depth" => 5,
    "elementary_gates" => ["U3_1", "U3_2", "CNot_1_2", "CNot_2_1", "Identity"],
    "target_gate" => QCOpt.CHGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact",

    "U3_θ_discretization" => -2*π:π/4:2*π,
    "U3_ϕ_discretization" => [0],
    "U3_λ_discretization" => [0],
    )
end

function decompose_magic_using_U3_CNot_2_1()

    println(">>>>> Magic basis using U3 and CNot_2_1 <<<<<")

    return Dict{String, Any}(
        "num_qubits" => 2,
        "maximum_depth" => 5,
        "elementary_gates" => ["U3_1", "U3_2", "CNot_2_1", "Identity"],
        "target_gate" => QCOpt.MGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",

        "U3_θ_discretization" => -π:π/2:π,
        "U3_ϕ_discretization" => -π:π/2:π,
        "U3_λ_discretization" => -π:π/2:π,
    )
end

function decompose_iSwapGate()

    println(">>>>> iSwap Gate <<<<<")

    return Dict{String, Any}(
        "num_qubits" => 2,
        "maximum_depth" => 10,
        "elementary_gates" => ["T_1", "T_2", "Tdagger_1", "Tdagger_2", "H_1", "H_2", "CNot_1_2", "CNot_2_1", "Identity"],
        "target_gate" => QCOpt.iSwapGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact")

end

function decompose_GroverDiffusion_using_Clifford()

    println(">>>>> Grover's Diffusion Operator <<<<<")

    return Dict{String, Any}(
        "num_qubits" => 2,
        "maximum_depth" => 6,
        "elementary_gates" => ["X_1", "X_2", "H_1", "H_2", "S_1", "S_2", "T_1", "T_2", "Y_1", "Y_2", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.GroverDiffusionGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact")
end

function decompose_GroverDiffusion_using_U3()

    println(">>>>> Grover's Diffusion Operator using U3 gate <<<<<")

    return Dict{String, Any}(
        "num_qubits" => 2,
        "maximum_depth" => 10,
        "elementary_gates" => ["U3_1", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.GroverDiffusionGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",

        "U3_θ_discretization" => -π:π/2:π,
        "U3_ϕ_discretization" => -π:π/2:π,
        "U3_λ_discretization" => [0])
end

function decompose_magic_using_SH_CNot_1_2()

    println(">>>>> M gate using S, H and CNOT_1_2 Gate <<<<<")

    return Dict{String, Any}(
        "num_qubits" => 2,
        "maximum_depth" => 10,
        "elementary_gates" => ["S_1", "S_2", "H_1", "H_2", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.MGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact")
end

function decompose_magic_using_SH_CNot_2_1()

    println(">>>>> M gate using S, H and CNOT_2_1 Gate <<<<<")

    return Dict{String, Any}(
        "num_qubits" => 2,
        "maximum_depth" => 10,
        "elementary_gates" => ["S_1", "S_2", "H_1", "H_2", "CNot_2_1", "Identity"],
        "target_gate" => QCOpt.MGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact")
end

function decompose_magic_using_U3_CNot_1_2()

    println(">>>>> Magic basis using U3 and CNot_1_2 <<<<<")

    return Dict{String, Any}(
        "num_qubits" => 2,
        "maximum_depth" => 5,
        "elementary_gates" => ["U3_1", "U3_2", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.MGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",

        "U3_θ_discretization" => -π:π/2:π,
        "U3_ϕ_discretization" => -π:π:π,
        "U3_λ_discretization" => -π:π/2:π,)
end

function decompose_toffoli_with_controlled_gates()

    # Reference: https://doi.org/10.1109/TCAD.2005.858352
    println(">>>>> Toffoli gate using controlled gates <<<<<")

    return Dict{String, Any}(
    "num_qubits" => 3,
    "maximum_depth" => 5,
    "elementary_gates" => ["CV_1_3", "CV_2_3", "CV_1_2", "CVdagger_1_3", "CVdagger_2_3", "CVdagger_1_2", "CNot_2_1", "CNot_1_2", "Identity"],
    "target_gate" => QCOpt.ToffoliGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact")
end

function decompose_Fredkin()

    println(">>>>> Fredkin gate using controlled gates <<<<<")

    return Dict{String, Any}(
    "num_qubits" => 3,
    "maximum_depth" => 7,
    # Reference: https://doi.org/10.1103/PhysRevA.53.2855
    "elementary_gates" => ["CV_1_2", "CV_2_3", "CV_1_3", "CVdagger_1_2", "CVdagger_2_3", "CVdagger_1_3", "CNot_1_2", "CNot_3_2", "CNot_2_3", "CNot_1_3", "Identity"],
    "target_gate" => QCOpt.CSwapGate(), #also Fredkin
    "objective" => "minimize_depth",
    "decomposition_type" => "exact")
end

function decompose_double_toffoli()

    println(">>>>> Double Toffoli using controlled gates <<<<<")

    # Reference: https://doi.org/10.1109/TCAD.2005.858352

    num_qubits = 4

    function target_gate()
        CV_3_4 = QCOpt.get_full_sized_gate("CV_3_4", num_qubits);
        CV_2_4 = QCOpt.get_full_sized_gate("CV_2_4", num_qubits);
        CVdagger_2_4 = QCOpt.get_full_sized_gate("CVdagger_2_4", num_qubits);
        CNot_1_3 = QCOpt.get_full_sized_gate("CNot_1_3", num_qubits);
        CNot_3_2 = QCOpt.get_full_sized_gate("CNot_3_2", num_qubits);

        return CV_2_4 * CNot_1_3 * CNot_3_2 * CV_3_4 * CVdagger_2_4 * CNot_3_2 * CNot_1_3
    end

    return Dict{String, Any}(
    "num_qubits" => 4,
    "maximum_depth" => 7,
    "elementary_gates" => ["CV_1_2", "CV_2_4", "CV_3_4", "CVdagger_1_2", "CVdagger_2_4", "CVdagger_3_4", "CNot_1_3", "CNot_3_2", "CNot_2_3", "Identity"],
    "target_gate" => target_gate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact")
end

function decompose_quantum_fulladder()
    #Reference-1: https://doi.org/10.1109/DATE.2005.249
    #Reference-2: https://doi.org/10.1109/TCAD.2005.858352

    println(">>>>> Quantum Full adder using controlled gates <<<<<")

    num_qubits = 4

    function target_gate()

        CV_1_2 = QCOpt.get_full_sized_gate("CV_1_2", num_qubits);
        CV_4_2 = QCOpt.get_full_sized_gate("CV_4_2", num_qubits);
        CVdagger_1_2 = QCOpt.get_full_sized_gate("CVdagger_1_2", num_qubits);
        CVdagger_3_2 = QCOpt.get_full_sized_gate("CVdagger_3_2", num_qubits);
        CNot_3_1 = QCOpt.get_full_sized_gate("CNot_3_1", num_qubits);
        CNot_4_3 = QCOpt.get_full_sized_gate("CNot_4_3", num_qubits);
        CNot_2_4 = QCOpt.get_full_sized_gate("CNot_2_4", num_qubits);

        return CNot_2_4 * CV_4_2 * CVdagger_1_2 * CVdagger_3_2 * CNot_4_3 * CNot_3_1 * CV_1_2
    end

    return Dict{String, Any}(
        "num_qubits" => num_qubits,
        "maximum_depth" => 7,
        "elementary_gates" => ["CV_1_2", "CV_4_2", "CV_3_2", "CVdagger_1_2", "CVdagger_4_2", "CVdagger_3_2", "CNot_3_1", "CNot_4_3", "CNot_2_4", "CNot_4_1", "Identity"],
        "target_gate" => target_gate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact")
end


text:
if !isfile(joinpath(dirname(@__FILE__), "Defaults.jl"))
    info("Default units, dimensions, and logic are set in ",
        "$(escape_string(joinpath(dirname(@__FILE__), "Defaults.jl")))")
    open(joinpath(dirname(@__FILE__), "Defaults.jl"), "w") do f
        print(f, """
        # Default dimensions and their abbreviations.
        # The dimension symbols are generated by tab completion: \\mbfL is 𝐋, etc.
        # This is kind of a compromise on polluting the namespace.
        @dimension 𝐋 "𝐋" Length
        @dimension 𝐌 "𝐌" Mass
        @dimension 𝐓 "𝐓" Time
        @dimension 𝐈 "𝐈" Current
        @dimension 𝚯 "𝚯" Temperature    # This one is \\mbfTheta
        @dimension 𝐉 "𝐉" Luminosity
        @dimension 𝐍 "𝐍" Amount

        # Define derived dimensions.
        @derived_dimension Area             𝐋^2
        @derived_dimension Volume           𝐋^3
        @derived_dimension Frequency        inv(𝐓)
        @derived_dimension Force            𝐌*𝐋/𝐓^2
        @derived_dimension Pressure         𝐌*𝐋^-1*𝐓^-2
        @derived_dimension Energy           𝐌*𝐋^2/𝐓^2
        @derived_dimension Momentum         𝐌*𝐋/𝐓
        @derived_dimension Power            𝐋^2*𝐌*𝐓^-3
        @derived_dimension Charge           𝐈*𝐓
        @derived_dimension Voltage          𝐈^-1*𝐋^2*𝐌*𝐓^-3
        @derived_dimension Resistance       𝐈^-2*𝐋^2*𝐌*𝐓^-3
        @derived_dimension Capacitance      𝐈^2*𝐋^-2*𝐌^-1*𝐓^4
        @derived_dimension Inductance       𝐈^-2*𝐋^2*𝐌*𝐓^-2
        @derived_dimension MagneticFlux     𝐈^-1*𝐋^2*𝐌*𝐓^-2
        @derived_dimension HField           𝐈/𝐋
        @derived_dimension BField           𝐈^-1*𝐌*𝐓^-2

        # Define base units. This is not to imply g is the base SI unit instead of kg.
        # See the documentation for further details.
        # #key:   Symbol  Display  Name      Dimension   Prefixes?
        @refunit  m       "m"      Meter     𝐋           true
        @refunit  s       "s"      Second    𝐓           true
        @refunit  A       "A"      Ampere    𝐈            true
        @refunit  K       "K"      Kelvin    𝚯           true
        @refunit  cd      "cd"     Candela   𝐉            true
        @refunit  g       "g"      Gram      𝐌           true
        @refunit  mol     "mol"    Mole      𝐍           true

        # Specify preferred unit for promotion.
        # This is separate from the @refunit macro for flexibility; consider that
        # the SI unit of mass is not g but instead kg, and yet some people use cgs units.
        # This macro should only be used with units having "pure" dimensions like 𝐋, 𝐓, 𝐈, etc.
        @preferunit m
        @preferunit s
        @preferunit A
        @preferunit K
        @preferunit cd
        @preferunit kg
        @preferunit mol

        # These lines allow for μ to be typed with option-m on a Mac.
        # The character encodings are different here so this is less crazy than it looks
        const µm = μm

        # Length
        #key: Symbol Display    Name        Equivalent to           10^n prefixes?
        @unit mi     "mi"       Mile        (201168//125)*m         false
        @unit yd     "yd"       Yard        (9144//10000)*m         false
        @unit ft     "ft"       Foot        (3048//10000)*m         false
        @unit inch   "in"       Inch        (254//10000)*m          false

        # Area
        @unit ac     "ac"       Acre        (316160658//78125)*m^2  false
        @unit a      "a"        Are         100m^2                  false

        # The hectare is used more frequently than any other power-of-ten of an are.
        const ha = Unitful.Units{(Unitful.Unit{:Are}(2,1//1),), typeof(𝐋)}()

        # Time
        @unit minute "min"      Minute      60s                     false
        @unit hr     "hr"       Hour        3600s                   false
        @unit d      "dy"       Day         86400s                  false
        @unit wk     "wk"       Week        604800s                 false

        # Angle
        @unit rad     "rad"     Radian      1                       false
        @unit °       "°"       Degree      pi/180                  false
        import Base: sin, cos, tan, cot, sec, csc
        for _y in [:sin, :cos, :tan, :cot, :sec, :csc]
            @eval (\$_y)(x::DimensionlessQuantity) = (\$_y)(uconvert(NoUnits, x))
        end

        # Temperature
        @unit °Ra    "°Ra"      Rankine     (5//9)*K                false
        @unit °C     "°C"       Celsius     1K                      false
        Unitful.offsettemp(::Unitful.Unit{:Celsius}) = 27315//100
        @unit °F     "°F"       Fahrenheit  (5//9)*K                false
        Unitful.offsettemp(::Unitful.Unit{:Fahrenheit}) = 45967//100

        # SI and related units
        @unit N      "N"        Newton      1kg*m/s^2               true
        @unit Pa     "Pa"       Pascal      1N/m^2                  true
        @unit J      "J"        Joule       1N*m                    true
        @unit W      "W"        Watt        1J/s                    true
        @unit C      "C"        Coulomb     1A*s                    true
        @unit V      "V"        Volt        1W/A                    true
        @unit Hz     "Hz"       Hertz       1/s                     true
        const q = 1.6021766208e-19*C        # CODATA 2014; `e` means 2.718...
        @unit eV     "eV"       eV          q*V                     true
        @unit Ω      "Ω"        Ohm         1V/A                    true
        @unit S      "S"        Siemens     1/Ω                     true
        @unit F      "F"        Farad       1s^4*A^2/(kg*m^2)       true
        @unit H      "H"        Henry       1J/(A^2)                true
        @unit T      "T"        Tesla       1kg/(A*s^2)             true
        @unit Wb     "Wb"       Weber       1kg*m^2/(A*s^2)         true

        # Constants (2014 CODATA values)    (uncertainties in final digits)
        const c0 = 299_792_458*m/s          # exact
        const c  = c0
        const μ0 = 4π*(1//10)^7*H/m         # exact (but gets promoted to Float64...)
        const µ0 = μ0                       # magnetic constant
        const ɛ0 = 1/(μ0*c^2)               # exact, electric constant
        const ϵ0 = ɛ0
        const Z0 = μ0*c                     # exact, impedance of free space
        const G  = 6.674_08e-11*m^3/kg/s^2  # (31) gravitational constant
        const h  = 6.626_070_040e-34*J*s    # (81) Planck constant
        const ħ  = h/2π                     # hbar
        const q  = 1.602_176_620_8e-19*C    # (98)  `e` means 2.718... so we use q
                                            # elementary charge
        const Φ0 = h/(2q)                   # Superconducting magnetic flux quantum
        const me = 9.109_383_56e-31*kg      # (11) electron rest mass
        const mn = 1.674_927_471e-27*kg     # (21) neutron rest mass
        const mp = 1.672_621_898e-27*kg     # (21) proton rest mass
        const μB = e*ħ/(2*me)               # Bohr magneton
        const Na = 6.022_140_857e23/mol       # (74) Avogadro constant
        const R  = 8.314_459_8*J/(mol*K)    # (48) molar gass constant
        const k  = 1.380_648_52e-23*(J/K)   # (79) Boltzmann constant
        const σ  = π^2*k^4/(60*ħ^3*c^2)     # Stefan-Boltzmann constant

        # Promotion rules

        # By default, pick the units specified by the @preferunit macro.
        # Our use of promote_rule here is only via promote_type;
        # We will never be promoting unit objects themselves.
        function promote_rule{S<:DimensionedUnits,T<:DimensionedUnits}(::Type{S}, ::Type{T})
            dS = dimension(S())
            dT = dimension(T())
            dS != dT && error("Dimensions are unequal in call to `promote_rule`.")
            typeof(dim2refunits(dS))
        end

        # You could also add rules like the following, which will not interfere with
        # the generic behavior for other dimensions:
        promote_rule{S<:EnergyUnit, T<:EnergyUnit}(::Type{S}, ::Type{T}) = typeof(J)
        promote_rule{S<:ForceUnit, T<:ForceUnit}(::Type{S}, ::Type{T}) = typeof(N)
        promote_rule{S<:PowerUnit, T<:PowerUnit}(::Type{S}, ::Type{T}) = typeof(W)
        promote_rule{S<:PressureUnit, T<:PressureUnit}(::Type{S}, ::Type{T}) = typeof(Pa)
        promote_rule{S<:ChargeUnit, T<:ChargeUnit}(::Type{S}, ::Type{T}) = typeof(C)
        promote_rule{S<:VoltageUnit, T<:VoltageUnit}(::Type{S}, ::Type{T}) = typeof(V)
        promote_rule{S<:ResistanceUnit, T<:ResistanceUnit}(::Type{S}, ::Type{T}) = typeof(Ω)
        promote_rule{S<:CapacitanceUnit, T<:CapacitanceUnit}(::Type{S}, ::Type{T}) = typeof(F)
        promote_rule{S<:InductanceUnit, T<:InductanceUnit}(::Type{S}, ::Type{T}) = typeof(H)
        promote_rule{S<:MagneticFluxUnit, T<:MagneticFluxUnit}(::Type{S}, ::Type{T}) = typeof(Wb)
        promote_rule{S<:BFieldUnit, T<:BFieldUnit}(::Type{S}, ::Type{T}) = typeof(T)
        """)
    end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions that creates a view of the string on the IO.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _view!(pagerd::Pager)

Write the view of pager `pagerd` to the view buffer.

"""
function _view!(pagerd::Pager)
    # Get the available display size.
    rows, cols = _get_pager_display_size(pagerd)

    # Get the necessary variables.
    start_row              = pagerd.start_row
    start_col              = pagerd.start_col
    lines                  = pagerd.lines
    num_lines              = pagerd.num_lines
    active_search_match_id = pagerd.active_search_match_id
    search_matches         = pagerd.search_matches
    buf                    = pagerd.buf
    freeze_columns         = pagerd.freeze_columns
    freeze_rows            = pagerd.freeze_rows
    title_rows             = pagerd.title_rows
    draw_ruler             = pagerd.draw_ruler

    # Make sure that the argument values are correct.
    start_row < 1 && (start_row = 1)
    start_col < 1 && (start_col = 1)

    # Printed lines.
    num_printed_lines = 0

    # Regex that matches ANSI escape characters.
    regex_ansi = r"\x1B(?:[@-Z\\-_]|\[[0-?]*[ -/]*[@-~])"

    # Indicate whether a line has been cropped.
    lines_cropped = 0

    # Indicate if all columns were displayed.
    columns_cropped = 0

    # Get the active search match.
    active_search_match = if active_search_match_id == 0
        (0, 0, 0)
    else
        search_matches[active_search_match_id]
    end

    # Assemble the vector with the lines to be printed.
    if freeze_rows > 0
        freeze_rows > rows && (freeze_rows = rows)
        start_row ≤ freeze_rows && (start_row = freeze_rows + 1)

        lines_indices = vcat(1:freeze_rows, start_row:num_lines)
    else
        lines_indices = collect(start_row:num_lines)
    end

    # Clamp title rows to the allowed interval. The title rows are select within
    # the frozen rows.
    title_rows = clamp(title_rows, 0, freeze_rows)

    # Store the last decoration applied to a line. It is required to draw the
    # ruler without interfering with the line decoration.
    last_decoration = ""

    for i ∈ lines_indices
        # If `i` is larger than `num_lines`, then it means that the user
        # requested to freeze more lines than we currently have.
        i > num_lines && break

        line = lines[i]

        # Check if we need to draw the ruler.
        if draw_ruler
            write(buf, _reset_crayon, _ruler_crayon)
            _draw_vertical_ruler!(buf, i, num_lines)
            write(buf, _reset_crayon, last_decoration)
        end

        # Get all search matches in this line.
        matches_i = filter(x -> x[1] == i, search_matches)

        # Split the lines into escape sequence and text.
        line_tokens, decoration, cropped_chars_i = _printing_recipe(
            line,
            i ≤ title_rows ? 1 : start_col,
            cols,
            matches_i,
            active_search_match,
            freeze_columns
        )

        columns_cropped = max(columns_cropped, cropped_chars_i)

        printable_chars = 0

        # Print the line.
        for j = 1:length(line_tokens)
            last_decoration = string(decoration[j])
            write(buf, last_decoration)
            write(buf, line_tokens[j])
        end

        # Check if we have a last decoration to apply.
        if length(line_tokens) < length(decoration)
            last_decoration = string(decoration[end])
            write(buf, last_decoration)
        end

        write(buf, '\n')
        num_printed_lines += 1

        if num_printed_lines ≥ rows
            lines_cropped = num_lines - (rows - freeze_rows) - (start_row - 1)
            break
        end
    end

    # Write the information to the sturcture.
    pagerd.lines_cropped   = lines_cropped
    pagerd.columns_cropped = columns_cropped

    # Since we modified the `buf`, we need to request redraw.
    _request_redraw!(pagerd)

    return nothing
end


text:
#=
This file handles the import/export of Variant Call Format files
=#

export bcf, vcf

infile = "/home/pdimens/PopGen.jl/data/source/filtered_oyster.vcf"

"""
    bcf(infile::String; silent::Bool = false)
Load a BCF file into memory as a PopObj object. Population and [optional]
location information need to be provided separately. Use `silent=true` to supress
printing during file loading.
- `infile` : path to VCF file
"""
function bcf(infile::String; silent::Bool = false)
    vcf_file = BCF.Reader(open(infile, "r"))

    # get sample names from header
    sample_names = header(vcf_file).sampleID

    # fill in pop/lat/long with missing
    population = fill("missing", length(sample_names))
    loc_xy = Vector{Union{Missing,Float32}}(undef, length(sample_names))

    # get loci names
    locinames = Vector{String}()

    ## genotype dataframe
    geno_df = DataFrame()

    # get genotypes
    for record in vcf_file
        # fix locus names to be syntax-safe
        chr_safe = replace(BCF.chrom(record), r"\.|\-|\=|\/" => "_")
        chr_safer = replace(chr_safe, "|" => "_")
        pos = VCF.pos(record) |> string
        push!(locinames, chr_safer*"_"*pos)

        # get the genotype information
        geno_raw = [split(i, ('/', '|')) for i in BCF.genotype(record, :, "GT")] |> sort

        # change missing data "." to "-1"
        geno_corr_miss = map(i -> replace(i, "." => "-1"), geno_raw)

        # convert everything to an integer
        geno_int = map(i -> parse.(Int8, i), geno_corr_miss)

        # add 1 to shift genos so 0 is 1 and -1 is 0 etc.
        geno_shift = map(i -> i .+ Int8(1), geno_int)
        geno_out = Vector{Union{Missing, Genotype}}()
        @inbounds for i in geno_shift
            if all(iszero.(i)) == true
                push!(geno_out, missing)
            else
                push!(geno_out, Tuple(i))
            end
        end
        insertcols!(geno_df, Symbol(chr_safer*"_"*pos) => geno_out)
    end
    close(vcf_file)
    if !silent
        @info "\n$(abspath(infile))\n$(length(sample_names)) samples detected\npopulation info must be added <---\n$(length(locinames)) loci detected"
    end

    # create loci dataframe
    insertcols!(geno_df, 1, :name => sample_names, :population => population)
    geno_parse = DataFrames.stack(geno_df, DataFrames.Not(1:2))
    rename!(geno_parse, [:name, :population, :locus, :genotype])
    categorical!(geno_parse, [:name, :population, :locus], compress = true)

    # make sure levels are sorted by order of appearance
    levels!(geno_parse.locus, unique(geno_parse.locus))
    levels!(geno_parse.name, unique(geno_parse.name))
    ploidy = DataFrames.combine(
        groupby(geno_parse, :name),
        :genotype => (i -> find_ploidy(i[i .!== missing])) => :ploidy
    ).ploidy

    # create samples df
    samples_df = DataFrame(
        name = sample_names,
        population = population,
        ploidy = ploidy,
        latitude = loc_xy,
        longitude = loc_xy
    )
    PopData(samples_df, geno_parse)
end

### VCF parsing ###

"""
    vcf(infile::String; silent::Bool = false)
Load a VCF file into memory as a PopObj object. Population and [optional]
location information need to be provided separately. Use `silent=true` to supress
printing during file loading.
- `infile` : path to VCF file
"""
function vcf(infile::String, silent::Bool = false)
    vcf_file = VCF.Reader(open(infile, "r"))

    # get sample names from header
    sample_names = header(vcf_file).sampleID

    # fill in pop/lat/long with missing
    population = fill("missing", length(sample_names))
    loc_xy = Vector{Union{Missing,Float32}}(undef, length(sample_names))

    # get loci names
    locinames = Vector{String}()

    ## genotype dataframe
    geno_df = DataFrame()

    # get genotypes
    for record in vcf_file
        # fix locus names to be syntax-safe
        chr_safe = replace(VCF.chrom(record), r"\.|\-|\=|\/" => "_")
        chr_safer = replace(chr_safe, "|" => "_")
        pos = VCF.pos(record) |> string
        push!(locinames, chr_safer*"_"*pos)

        # get the genotype information
        geno_raw = [split(i, ('/', '|')) for i in VCF.genotype(record, :, "GT")] |> sort

        # change missing data "." to "-1"
        geno_corr_miss = map(i -> replace(i, "." => "-1"), geno_raw)

        # convert everything to an integer
        geno_int = map(i -> parse.(Int8, i), geno_corr_miss)

        # add 1 to shift genos so 0 is 1 and -1 is 0 etc.
        geno_shift = map(i -> i .+ Int8(1), geno_int)
        geno_out = Vector{Union{Missing, Genotype}}()
        @inbounds for i in geno_shift
            if all(iszero.(i)) == true
                push!(geno_out, missing)
            else
                push!(geno_out, Tuple(i))
            end
        end
        insertcols!(geno_df, Symbol(chr_safer*"_"*pos) => geno_out)
    end
    close(vcf_file)
    if !silent
        @info "\n$(abspath(infile))\n$(length(sample_names)) samples detected\npopulation info must be added <---\n$(length(locinames)) loci detected"
    end

    # create loci dataframe
    insertcols!(geno_df, 1, :name => sample_names, :population => population)
    geno_parse = DataFrames.stack(geno_df, DataFrames.Not(1:2))
    rename!(geno_parse, [:name, :population, :locus, :genotype])
    categorical!(geno_parse, [:name, :population, :locus], compress = true)

    # make sure levels are sorted by order of appearance
    levels!(geno_parse.locus, unique(geno_parse.locus))
    levels!(geno_parse.name, unique(geno_parse.name))
    ploidy = DataFrames.combine(
        groupby(geno_parse, :name),
        :genotype => (i -> find_ploidy(i[i .!== missing])) => :ploidy
    ).ploidy

    # create samples df
    samples_df = DataFrame(
        name = sample_names,
        population = population,
        ploidy = ploidy,
        latitude = loc_xy,
        longitude = loc_xy
    )
    PopData(samples_df, geno_parse)
end


text:
module EquivariantOperators

include("./neural_networks.jl")
export Grid, put_point_source!, field, field_rank, field_norm
export FieldProd, field_prod, FieldConv, field_conv
export EquivLayer, LinearOperator
#
end # module

max_stars_repo_name: aced-differentiate/equiv-docs

max_issues_repo_name: aced-differentiate/equiv-docs

max_forks_repo_name: aced-differentiate/equiv-docs


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_script/patterns.jl
# Patterns to use in transormations that flatten expressions.
#############################################################################

# (1) register log(a^x) = x*log(a) DONE
src_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :log),
                                   2 => Template_Node(:op, :^),
                                   3 => Template_Node(:num, :a; check = a -> a >= 0.0),
                                   4 => Template_Node(:expr, :x))
src_dag = [4 => 2, 3 => 2, 2 => 1]
src = Template_Graph(src_nds, src_dag)
dest_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :*),
                                    2 => Template_Node(:expr, :x),
                                    3 => Template_Node(:op, :log),
                                    4 => Template_Node(:num, :a))
dest_dag = [4 => 3, 3 => 1, 2 => 1]
dest = Template_Graph(dest_nds, dest_dag)
register_substitution!(src, dest)

# (2) register exp(x)*exp(y) -> exp(x+y)  DONE
src_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :*),
                                   2 => Template_Node(:op, :exp),
                                   3 => Template_Node(:op, :exp),
                                   4 => Template_Node(:expr, :x),
                                   5 => Template_Node(:expr, :y))
src_dag = [5 => 3, 4 => 2, 3 => 1, 2 => 1]
src = Template_Graph(src_nds, src_dag)
dest_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :exp),
                                    2 => Template_Node(:op, :+),
                                    3 => Template_Node(:expr, :x),
                                    4 => Template_Node(:expr, :y))
dest_dag = [4 => 2, 3 => 2, 2 => 1]
dest = Template_Graph(dest_nds, dest_dag)
register_substitution!(src, dest)

# (3) register (a^{x})^{b} = (a^{b})^{x} # DONE
src_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :^),
                                   2 => Template_Node(:op, :^),
                                   3 => Template_Node(:num, :b),
                                   4 => Template_Node(:num, :a),
                                   5 => Template_Node(:expr, :x))
src_dag = [5 => 2, 4 => 2, 3 => 1, 2 => 1]
src = Template_Graph(src_nds, src_dag)
dest_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :^),
                                    2 => Template_Node(:op, :^),
                                    3 => Template_Node(:expr, :x),
                                    4 => Template_Node(:num, :a),
                                    5 => Template_Node(:num, :b))
dest_dag = [5 => 2, 4 => 2, 3 => 1, 2 => 1]
dest = Template_Graph(dest_nds, dest_dag)
register_substitution!(src, dest)

# (4) register (x^{a})^{b} = x^{(ab)}
#=
src_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :^),
                                   2 => Template_Node(:op, :^),
                                   3 => Template_Node(:num, :a),
                                   4 => Template_Node(:expr, :x),
                                   5 => Template_Node(:num, :b))
src_dag = [5 => 2, 4 => 2, 3 => 1, 2 => 1]
src = Template_Graph(src_nds, src_dag)
dest_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :^),
                                    2 => Template_Node(:expr, :x),
                                    3 => Template_Node(:op, :*),
                                    4 => Template_Node(:num, :a),
                                    5 => Template_Node(:num, :b))
dest_dag = [5 => 3, 4 => 3, 3 => 1, 2 => 1]
dest = Template_Graph(dest_nds, dest_dag)
register_substitution!(src, dest)
=#

# (5) register a^{\log(x)} = x^{\log(a)} DONE
src_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :^),
                                   2 => Template_Node(:num, :a),
                                   3 => Template_Node(:op, :log),
                                   4 => Template_Node(:expr, :x))
src_dag = [4 => 3, 3 => 1, 2 => 1]
src = Template_Graph(src_nds, src_dag)
dest_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :^),
                                    2 => Template_Node(:expr, :x),
                                    3 => Template_Node(:op, :log),
                                    4 => Template_Node(:num, :a))
dest_dag = [4 => 3, 3 => 1, 2 => 1]
dest = Template_Graph(dest_nds, dest_dag)
register_substitution!(src, dest)

# (6) register \log(xy) = \log(x) + \log(y) DONE
src_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :log),
                                   2 => Template_Node(:op, :*),
                                   3 => Template_Node(:expr, :x),
                                   4 => Template_Node(:expr, :y))
src_dag = [4 => 2, 3 => 2, 2 => 1]
src = Template_Graph(src_nds, src_dag)
dest_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :+),
                                    2 => Template_Node(:op, :log),
                                    3 => Template_Node(:op, :log),
                                    4 => Template_Node(:expr, :x),
                                    5 => Template_Node(:expr, :y))
dest_dag = [5 => 3, 4 => 2, 3 => 1, 2 => 1]
dest = Template_Graph(dest_nds, dest_dag)
register_substitution!(src, dest)

# (7) register \log(x/y) = \log(x) - \log(y)  DONE
src_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :log),
                                   2 => Template_Node(:op, :/),
                                   3 => Template_Node(:expr, :x),
                                   4 => Template_Node(:expr, :y))
src_dag = [4 => 2, 3 => 2, 2 => 1]
src = Template_Graph(src_nds, src_dag)
dest_nds = Dict{Int, Template_Node}(1 => Template_Node(:op, :-),
                                    2 => Template_Node(:op, :log),
                                    3 => Template_Node(:op, :log),
                                    4 => Template_Node(:expr, :x),
                                    5 => Template_Node(:expr, :y))
dest_dag = [5 => 3, 4 => 2, 3 => 1, 2 => 1]
dest = Template_Graph(dest_nds, dest_dag)
register_substitution!(src, dest)


text:
#=
Queue
=#

using SimpleDataStructures

sq = SimpleQueue{Int}()

enqueue!(sq, 1)
enqueue!(sq, 2)
dequeue!(sq)
enqueue!(sq, 3)
enqueue!(sq, 4)
dequeue!(sq)
dequeue!(sq)
dequeue!(sq)
dequeue!(sq)


text:
# Internal function to run perft
function runperft(board::Board, depth::Int, ply::Int, movestack::Vector{MoveStack})
    if isone(depth)
        n = count_moves(board)
        return n
    else
        movelist = movestack[ply + 1]
        gen_moves!(movelist, board)
        num_moves = 0
        @inbounds for i in eachindex(movelist)
            u = apply_move!(board, movelist[i])
            num_moves += runperft(board, depth - 1, ply + 1, movestack)
            undo_move!(board, movelist[i], u)
        end
        clear!(movelist)
        return num_moves
    end
end

"""
    perft(board::Board, depth::Int)

Run the perft routine to a given `depth`.
"""
function perft(board::Board, depth::Int)
    movestack = [MoveStack(MAX_MOVES) for d in 1:depth]
    runperft(board, depth, 0, movestack)
end
# Perft also generates zobrist keys, and PSQT + PVAL evaluations alnog the way.

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
module SfWeaponry

using Unitful
using ..SfElectro

export cratering_depth, cratering_strength, cratering_volume, jet_penetrator, newtonian_penetrator, coilgun_velocity, coilgun_field_strength,
	coilgun_length

function cratering_strength(yield_strength::Unitful.Pressure)
    3u"J/m^3" * Unitful.ustrip(u"Pa", yield_strength)
end

function cratering_volume(e::Unitful.Energy, yield_strength::Unitful.Pressure)
    e / cratering_strength(yield_strength)
end

function cratering_depth(e::Unitful.Energy, yield_strength::Unitful.Pressure)
    r = cratering_volume(e, yield_strength)
    # hemisphere volume = (2/3)πr^3

    cbrt(3r/2π)
end

"""
    jet_penetrator(l::Unitful.Length, ρj::Unitful.Density, ρt::Unitful.Density)

Approximation for penetration depth of hypervelocity jet of density `ρj` and length `l` into solid material of density `ρt`.
"""
jet_penetrator(l::Unitful.Length, ρj::Unitful.Density, ρt::Unitful.Density) = l * sqrt(ρj / ρt) |> unit(l)

"""
    newtonian_penetrator(l::Unitful.Length, ρp::Unitful.Density, ρt::Unitful.Density)

Approximation for penetration depth of high velocity impactor of density `ρj` and length `l` into solid material of density `ρt`.
"""
newtonian_penetrator(l::Unitful.Length, ρp::Unitful.Density, ρt::Unitful.Density) = l * (ρp / ρt) |> unit(l)

"""
    coilgun_velocity(m::Unitful.Mass, r::Unitful.Length, l::Unitful.Length, b::Tesla)

Use Luke Campbell's coilgun approximation to compute muzzle velocity of a cylindrical projectile of mass `m` and radius `r` firing down a barrel of length `l` in a unitform magnetic field of strength `b`.
"""
function coilgun_velocity(m::Unitful.Mass, r::Unitful.Length, l::Unitful.Length, b::Unitful.BField)
	ed = energy_density(b)
	v = π * r^2 * l
	e = ed * v

	return sqrt(2e / m) |> u"km/s"
end

"""
    coilgun_field_strength(m::Unitful.Mass, r::Unitful.Length, l::Unitful.Length, v::Unitful.Velocity)

Use Luke Campbell's coilgun approximation to compute magnetic field strength of a coilgun of length `l` firing a cylindrical projectile of mass `m` and radius `r` at velocity `v`.
"""
function coilgun_field_strength(m::Unitful.Mass, r::Unitful.Length, l::Unitful.Length, v::Unitful.Velocity)
	ke = 0.5m*v^2
	v = π * r^2 * l

	return field_strength(ke/v)
end

"""
    coilgun_length(m::Unitful.Mass, r::Unitful.Length, b::Unitful.BField, v::Unitful.Velocity)

Use Luke Campbell's coilgun approximation to compute length of a coilgun with a field strength of `b`, firing a cylindrical projectile of mass `m` and radius `r` at muzzle velocity `v`.
"""
function coilgun_length(m::Unitful.Mass, r::Unitful.Length, b::Unitful.BField, v::Unitful.Velocity)
	ed = energy_density(b)
	ke = 0.5m*v^2
	a = π * r^2

	return ke / (ed*a)|>u"m"
end

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics