text:
###################
# Helper functions
###################

determine_return_type(fun, S) = Base.promote_op(fun, S)
determine_return_type(fun, S::Type{Tuple{A}}) where {A} = Base.promote_op(fun, A)
determine_return_type(fun, S::Type{Tuple{A,B}}) where {A,B} = Base.promote_op(fun, A, B)
determine_return_type(fun, S::Type{Tuple{A,B,C}}) where {A,B,C} = Base.promote_op(fun, A, B, C)
determine_return_type(fun, S::Type{SVector{1,A}}) where {A} = Base.promote_op(fun, A)
determine_return_type(fun, S::Type{SVector{2,A}}) where {A} = Base.promote_op(fun, A, A)
determine_return_type(fun, S::Type{SVector{3,A}}) where {A} = Base.promote_op(fun, A, A, A)
determine_return_type(fun, S::Type{SVector{4,A}}) where {A} = Base.promote_op(fun, A, A, A, A)

promote_dictionary(dict, fun) = _promote_dictionary(dict, fun, determine_return_type(fun, domaintype(dict)))
_promote_dictionary(dict, fun, ::Union{}) = dict
_promote_dictionary(dict, fun, ::Type{T}) where {T} = BasisFunctions.ensure_coefficienttype(promote_type(codomaintype(dict),T), dict)

function directsolver(A::DiagonalOperator; verbose = false, options...)
    verbose && println("Using direct solver: inverse of diagonal matrix")
    inv(A)
end

function directsolver(A; directsolver = :svd, verbose = false, options...)
    if directsolver == :svd
        verbose && println("Using direct solver: SVD")
        SVD_solver(A)
    elseif directsolver == :qr
        verbose && println("Using direct solver: QR")
        QR_solver(A)
	elseif directsolver == :regsvd
        verbose && println("Using direct solver: regularized SVD")
        regularized_SVD_solver(A; verbose=verbose, options...)
    else
        DS = directsolver(A; verbose=verbose, options...)
        verbose && println("Using direct solver: $DS")
        DS
    end
end

function iterativesolver(A; iterativesolver = :lsqr, verbose = false, options...)
    if iterativesolver == :lsqr
        verbose && println("Using iterative solver: LSQR")
        LSQR_solver(A; verbose=verbose, options...)
    elseif iterativesolver == :lsmr
        verbose && println("Using iterative solver: LSMR")
        LSMR_solver(A; verbose=verbose, options...)
    else
        DS = iterativesolver(A)
        verbose && println("Using direct solver: $DS")
        DS
    end
end

####################
# The Fun interface
####################

# The `Fun` interface first turns the approximation problem into a subtype of
# `ApproximationProblem`. Next, it calls `approximate`. Finally, it discards all
# the operators that were computed and simply returns the function.
# Users wanting to access the operators can call `approximate` directly.

guess_coefficienttype(dict, fun) = promote_type(codomaintype(dict), determine_return_type(fun, domaintype(dict)))

function Fun(fun, dict::Dictionary, args...;
        coefficienttype = guess_coefficienttype(dict, fun),
        options...)
    ap = approximationproblem(coefficienttype, dict, args...)
    Fun(fun, ap; options...)
end

function Fun(fun, platform::Platform, args...; options...)
    ap = approximationproblem(platform, args...)
    Fun(fun, ap; options...)
end

function Fun(fun, ap::ApproximationProblem; verbose = false, options...)
    if verbose
        println("Fun: using the following dictionary:")
        show(dictionary(ap))
		println()
    end
    F, A, B, C, S = approximate(fun, ap; verbose=verbose, options...)
    F
end

function Fun(fun, ap::AdaptiveApproximation; verbose = false, options...)
    verbose && println("Fun: adaptive approximation using platform $(ap.platform)")
    approximate(fun, ap; verbose=verbose, options...)
end

function Fun(dict::Dictionary, coefficients::AbstractArray; verbose = false, options...)
	verbose && println("Fun: coefficients are given, no approximation problem is solved")
	Expansion(dict, coefficients)
end

############################
# The approximate function
############################

approximate(fun, dict::Dictionary, args...; options...) =
    approximate(fun, approximationproblem(dict, args...); options...)
approximate(fun, platform::Platform, args...; options...) =
    approximate(fun, approximationproblem(platform, args...); options...)

approximate(fun, ap::ApproximationProblem;
            samplingstyle = SamplingStyle(ap),
            solverstyle = SolverStyle(samplingstyle, ap),
            problemstyle = ProblemStyle(ap), options...) =
    approximate(problemstyle, promote(samplingstyle, solverstyle)..., fun, ap; options...)

# If the sampling has product structure and a single solverstyle is specified,
# we turn the solverstyle into a product style as well.
promote(samplingstyle::SamplingStyle, solverstyle::SolverStyle) = (samplingstyle,solverstyle)
promote(samplingstyle::ProductSamplingStyle, solverstyle::ProductSolverStyle) = (samplingstyle,solverstyle)
promote(samplingstyle::ProductSamplingStyle, solverstyle::SolverStyle) =
	(samplingstyle,ProductSolverStyle(map(x->solverstyle, samplingstyle.styles)))

# Construct the approximation problem and solve it
function approximate(::DictionaryOperatorStyle, samplingstyle::SamplingStyle, solverstyle::SolverStyle, fun, ap::ApproximationProblem;
            verbose = false, options...)
    dict = dictionary(ap)

    verbose && println("Approximate: using sampling style $samplingstyle")
	verbose && println("Approximate: using solver style $solverstyle")

    S = samplingoperator(samplingstyle, ap; verbose=verbose, options...)
	if verbose
		println()
		println("Approximate: sampling operator with size $(size(S)) is")
		println(S)
		println()
	end
    A, B = normalized_discretization(fun, samplingstyle, ap, S; verbose=verbose, options...)
    C = solve(solverstyle, ap, A, B; S=S, verbose=verbose, samplingstyle=samplingstyle, options...)
    F = Expansion(dictionary(ap), C)
    res = norm(A*C-B)
    verbose && println("Approximate: ended with residual $res\n")
    F, A, B, C, S, samplingparameter(samplingstyle, ap; verbose=verbose, options...)
end

# Construct the approximation problem and solve it
function approximate(pstyle::GenericOperatorStyle, samplingstyle::SamplingStyle, solverstyle::SolverStyle, fun, ap;
            verbose = false, options...)
    dict = dictionary(ap)

    verbose && println("Approximate: using problemstyle style $pstyle")
    verbose && println("Approximate: using sampling style $samplingstyle")
	verbose && println("Approximate: using solver style $solverstyle")

    S = samplingoperator(pstyle, ap; samplingstyle=samplingstyle, verbose=verbose, options...)
	if verbose
		println()
		println("Approximate: sampling operator with size $(size(S)) is")
		println(S)
		println()
	end

    A = AZ_A(pstyle, ap; samplingstyle=samplingstyle, options...)
    B = apply(S, fun; options...) # calculated just to keep the same interface.
    C = solve(solverstyle, ap, A, fun; verbose=verbose, problemstyle=pstyle, samplingstyle=samplingstyle, options...)
    F = Expansion(dictionary(ap), C)
    res = norm(apply(S,A; options...)*C-B)
    verbose && println("Approximate: ended with residual $res\n")
    F, A, B, C, S, samplingparameter(samplingstyle, ap; verbose=verbose, options...)
end

solve(style::SolverStyle, ap, A::DictionaryOperator, B; options...) =
    solver(style, ap, A; B=B, options...) * B
solve(style::SolverStyle, ap, A::AbstractOperator, fun; options...) =
    apply(solver(style, ap, A; options...), fun; options...)


text:
module NonlinearSchrodinger

using FFTW
using OrdinaryDiffEq, DiffEqOperators
using Memoization
using LaTeXStrings
using Elliptic
using Elliptic.Jacobi
using RecipesBase
using PolynomialRoots

export compute_IoM!, compute_spectrum!
export print
export solve!
export ψ₀_periodic, ψ₀_DT
export λ_maximal, λ_given_m, λ_given_f
export params
export PHF, find_peaks
#export solve_dt!

export Density_ψ

export Sim, Box, Calc
export T1A!, T1B!, T2A!, T2B!
export T4A_TJ!, T4B_TJ!, T6A_TJ!, T6B_TJ!, T8A_TJ!, T8B_TJ!
export T4A_SF!, T4B_SF!, T6A_SF!, T6B_SF!, T8A_SF!, T8B_SF!
export T4A_CMP!, T4B_CMP!, T6A_CMP!, T6B_CMP!, T8A_CMP!, T8B_CMP!
export T6A_Ss14!, T6B_Ss14!, T6A_Ys7!, T6B_Ys7!, T6A_KLs9!, T6B_KLs9!
export T8A_Ss15!, T8B_Ss15!, T8A_KLs17!, T8B_KLs17!
export T1A_H!, T2A_H!
export T1A_SS!

include("Types.jl")
include("Simulation.jl")
include("Utilities.jl")
include("Recipes.jl")
include("Darboux.jl")

end #module

ext: jl


text:
"""
	boost!(myscratch)

Function to perform a single boostingstep,
which is called by the wrapper-function 'boosting'.

# Arguments:
- `myscratch::Boostscratch -> Boostscratch-Data-Container for boosting`

# Returned objects
- `newselval::Float64`: Current score function value
- Update of myscratch (performed one boosting step)

# Examples
```julia-repl
julia> lab = ["X1","X2","X3","X4","X5","X6","X7","X8","X9","X10"]
julia> myscratch = Boostscratch(1,Array{Float64,1}(),Array{Float64,1}(),Unibeta(Array{Float64,1}(),Array{String,1}()),Covarmat(Array{Float64,2}(),Array{String,1}()),Array{Float64,1}(10),Array{String,1}(),Array{String,1}(),Array{String,1}(),0.1,10,1,Array{Float64,1}(),1,1,10,"Y",lab)

julia> myscratch.pooledunibeta = calc_unibeta(lab,"Y")
julia> myscratch.actualbeta = zeros(length(myscratch.pooledunibeta.unibeta))
julia> myscratch.actualnom = myscratch.pooledunibeta.unibeta

julia> myscratch.wantedlabels = selectionofcovs(myscratch.pooledunibeta, myscratch.a)
julia> myscratch.pooledcovarmat = calc_covarmat(myscratch.wantedlabels)
julia> myscratch.usedlabels = deepcopy(myscratch.wantedlabels)

julia> boost!(myscratch)
10.490144453484817
```
"""
function boost!(myscratch::Boostscratch)
	# Get score of all variables (should be devided by variance, but due to standardization this is 1)
	myscratch.actualscore = myscratch.actualnom.^2 # ./varx

	# Get the variable with the highest score - this is the selected one
	myscratch.actualsel = findmax(myscratch.actualscore)[2]

	# Save that this is the currently selected one in selections - vector
	myscratch.selections = vcat(myscratch.selections, myscratch.pooledunibeta.labels[myscratch.actualsel])

	# Calculate the corresponding update
	myscratch.actualupdates[myscratch.actualstepno] =
		myscratch.nu * myscratch.actualnom[myscratch.actualsel]

	# Calculate the update of the corresponding effect estimate
	myscratch.actualbeta[myscratch.actualsel] +=
		myscratch.actualupdates[myscratch.actualstepno]

	# Get covariances from all potential variables with selected one from covariance matrix
	covar = getcovarfromtriangular(myscratch.pooledcovarmat.covarmat,
				findfirst(isequal(myscratch.pooledunibeta.labels[myscratch.actualsel]),
					myscratch.pooledcovarmat.labels
				)
			)

	# Update corresponding score vector
	@simd for i = 1 : length(covar)
		@inbounds myscratch.actualnom[findfirst(isequal(myscratch.pooledcovarmat.labels[i]),
												myscratch.pooledunibeta.labels)] -=
			myscratch.actualupdates[myscratch.actualstepno] * covar[i]
	end
	# Calculate the current score of selected variable
	newselval = myscratch.actualnom[myscratch.actualsel]^2
	# Return current score of selected variable
	return newselval
end

"""
	reboost!(myscratch)

Function to repeat the boostingsteps (which were already performed)
with new data resulting from a new data call.

# Arguments
- `myscratch::Boostscratch`: Bootscratch-Data-Container for boosting

# Returned objects
- No objects returned.
- `myscratch.actualnom` is updated.

# Examples
```julia-repl
julia> lab = ["X1","X2","X3","X4","X5","X6","X7","X8","X9","X10"]
julia> myscratch = Boostscratch(1,Array{Float64,1}(),Array{Float64,1}(),Unibeta(Array{Float64,1}(),Array{String,1}()),Covarmat(Array{Float64}(undef,0,0),Array{String,1}()),Array{Float64,1}(10),Array{String,1}(),Array{String,1}(),Array{String,1}(),0.1,10,1,Array{Float64,1}(),2,1,10,"Y",lab)

julia> myscratch.pooledunibeta = calc_unibeta(lab,"Y")
julia> myscratch.actualbeta = zeros(length(myscratch.pooledunibeta.unibeta))
julia> myscratch.actualnom = myscratch.pooledunibeta.unibeta

julia> myscratch.wantedlabels = selectionofcovs(myscratch.pooledunibeta, myscratch.a)
julia> myscratch.pooledcovarmat = calc_covarmat(myscratch.wantedlabels)
julia> myscratch.usedlabels = deepcopy(myscratch.wantedlabels)

julia> boost!(myscratch)
julia> myscratch.wantedlabels = ["X9","X10"]
julia> calc_covarmat!(myscratch)
julia> reboost!(myscratch)
"""
function reboost!(myscratch::Boostscratch)
	# Repeat all previous boosting steps
	for i = 1 : (myscratch.actualstepno-1)
		# Get covariance of all variables with the one selected in boosting step i
		covar = getcovarfromtriangular(myscratch.pooledcovarmat.covarmat,
					findfirst(isequal(myscratch.selections[i]), myscratch.pooledcovarmat.labels))

		# Update the score functions one the newly called variables
		for j = 1 : length(myscratch.wantedlabels)
			myscratch.actualnom[findfirst(isequal(myscratch.wantedlabels[j]),myscratch.pooledunibeta.labels)] -=
				myscratch.actualupdates[i] * covar[findfirst(isequal(myscratch.wantedlabels[j]), myscratch.pooledcovarmat.labels)]
		end
	end
end


text:
using LinearAlgebra
using Mueller
using RecipesBase: apply_recipe
using StableRNGs
using Test
using Unitful: °

rng = StableRNG(128584)

@testset "Mueller.jl" begin

    @testset "rotations" for EL in [linear_polarizer, waveplate, qwp, hwp, mirror]
        M = EL()
        for θ in 2 * π .* rand(rng, 1000)
            R = Mueller.rotation(θ)
            @test R' * M * R ≈ rotate(M, θ)
            @test rotate(rotate(M, θ), -θ) ≈ M
        end
    end

    @testset "linear polarizers" for T in [Float32, Float64, BigFloat]
        @testset "horizontal" begin
            horizontal = linear_polarizer(T)
            @test horizontal ≈ 0.5 * T[1 1 0 0;
                                       1 1 0 0
                                       0 0 0 0
                                       0 0 0 0]

            # parallal light perfect
            S = T[1, 1, 0, 0]
            Sp = horizontal * S
            @test Sp ≈ T[1, 1, 0, 0]

            # orthogonal light blocked
            S = T[1, -1, 0, 0]
            Sp = horizontal * S
            @test Sp ≈ zeros(T, 4)

            # mixed light attenuated
            S = T[1, 1/√2, 1/√2, 0]
            Sp = horizontal * S
            @test Sp ≈ 1/2 * T[1 + 1/√2, 1  + 1/√2, 0, 0]
        end
        @testset "vertical" begin
            vertical = linear_polarizer(T, 90°)
            @test vertical ≈ 0.5 * T[1 -1 0 0;
                                     -1  1 0 0
                                     0 0 0 0
                                     0 0 0 0]

            # parallal light perfect
            S = T[1, -1, 0, 0]
            Sp = vertical * S
            @test Sp ≈ T[1, -1, 0, 0]

            # orthogonal light blocked
            S = T[1, 1, 0, 0]
            Sp = vertical * S
            @test Sp ≈ zeros(T, 4)

            # mixed light attenuated
            S = T[1, -1/√2, 1/√2, 0]
            Sp = vertical * S
            @test Sp ≈ 1/2 * T[1 + 1/√2, -1 - 1/√2, 0, 0]
        end
        @testset "double polarization" begin
            M0 = linear_polarizer(T)
            M90 = linear_polarizer(T, 90°)
            M = M0 * M90
            @test all(≈(0), M)

            # more randomly
            for θ in 180 .* rand(rng, 1000)
                M1 = linear_polarizer(T, θ*°)
                M2 = linear_polarizer(T, (θ + 90)*°)
                M = M0 * M90
                @test all(≈(0), M)
            end
        end
    end

    @testset "waveplates" for T in [Float32, Float64, BigFloat]
        M = waveplate(T)
        @test eltype(M) == T

        @testset "hwp" begin
            M = hwp(T)
            @test eltype(M) == T
            @test M ≈ waveplate(T, π, 0)

            Sp = M * T[1, 1, 0, 0]
            @test Sp ≈ T[1, 1, 0, 0] atol=1e-10
            Sp = M * T[1, 0, 1, 0]
            @test Sp ≈ T[1, 0, -1, 0] atol=1e-10

            for angle in π/2 .* (rand(rng, 1000) .- 0.5)
                M = hwp(T, angle)
                Sp = M * T[1, 1, 0, 0]
                AoLP = 0.5 * atan(Sp[3], Sp[2])
                @test AoLP ≈ 2 * angle atol=1e-6
            end
        end

        @testset "qwp" begin
            M = qwp(T)
            @test eltype(M) == T
            @test M ≈ waveplate(T, π/2, 0)

            Sp = M * T[1, 1, 0, 0]
            @test Sp ≈ T[1, 1, 0, 0] atol=1e-10
            Sp = M * T[1, 0, 1, 0]
            @test Sp ≈ T[1, 0, 0, 1] atol=1e-10
        end

        @testset "generate circular polarization" begin
            M = qwp(T) * linear_polarizer(T, 45°)
            S = T[1, 0, 0, 0]
            Sp = M * S
            @test Sp ≈ T[0.5, 0, 0, 0.5] atol=1e-10
        end

    end

    @testset "plot recipe" begin
        S = [1, 0, 0, 1]
        obj = Mueller.PolEllipse([S])
        recipes = apply_recipe(Dict{Symbol,Any}(), obj)
        recipe = only(recipes)
        x, y = recipe.args[1:2]
        @test minimum(x) ≈ -1/sqrt(2) atol=1e-5
        @test maximum(x) ≈ 1/sqrt(2) atol=1e-5
        @test minimum(y) ≈ -1/sqrt(2) atol=1e-5
        @test maximum(y) ≈ 1/sqrt(2) atol=1e-5
    end
end


text:
module BoundaryGrids

using Gridap
using Gridap.Helpers

export BoundaryGrid
import Gridap: Triangulation
import Gridap: BoundaryTriangulation
import Gridap: points
import Gridap: cells
import Gridap: celltypes
import Gridap: cellorders

struct BoundaryGrid{D,Z} <: Grid{D,Z}
  grid::Grid{D,Z}
  descriptor::BoundaryDescriptor

  function BoundaryGrid(
    grid::Grid{D,Z}, descriptor::BoundaryDescriptor) where {D,Z}
    @assert D == Z + 1
    new{D,Z}(grid,descriptor)
  end

end

for op in (:points,:cells,:celltypes,:cellorders)
  @eval begin
    $op(g::BoundaryGrid) = $op(g.grid)
  end
end

function Triangulation(grid::BoundaryGrid)
  trian = Triangulation(grid.grid)
  BoundaryTriangulation(trian,grid.descriptor)
end

function BoundaryGrid(model::DiscreteModel)
  d = celldim(model)
  graph = GridGraph(model)
  facet_to_cells = connections(graph,d-1,d)
  nfacets = length(facet_to_cells)
  mask = fill(false,nfacets)
  _find_on_boundary!(mask,facet_to_cells)
  BoundaryGrid(model,mask)
end

function _find_on_boundary!(mask,facet_to_cells)
  for (facet,cells) in enumerate(facet_to_cells)
    if length(cells) == 1
      mask[facet] = true
    end
  end
end

function BoundaryGrid(model::DiscreteModel,tags,icell::Int=1)
  _tags = _setup_tags(FaceLabels(model),tags)
  BoundaryGrid(model,_tags,icell)
end

function BoundaryGrid(model::DiscreteModel,tags::Vector{Int},icell::Int=1)
  d = celldim(model)
  labels = FaceLabels(model)
  oldfacet_to_label = labels_on_dim(labels,d-1)
  tag_to_labels = labels.tag_to_labels
  nfacets = length(oldfacet_to_label)
  oldfacet_to_mask = fill(false,nfacets)
  _setup_mask!(oldfacet_to_mask,oldfacet_to_label,tag_to_labels,tags)
  BoundaryGrid(model,oldfacet_to_mask,icell)
end

function BoundaryGrid(model::DiscreteModel,mask::AbstractVector{Bool},icell::Int=1)
  d = celldim(model)
  oldfacet_to_mask = mask
  fgrid = Grid(model,d-1)
  facet_to_oldfacet = findall(oldfacet_to_mask)
  grid = GridPortion(fgrid,facet_to_oldfacet)
  graph = GridGraph(model)
  oldfacet_to_cells = connections(graph,d-1,d)
  cell_to_oldfacets = connections(graph,d,d-1)
  oldfacet_to_cell = get_local_item(oldfacet_to_cells,icell)
  oldfacet_to_lfacet = find_local_index(oldfacet_to_cell, cell_to_oldfacets)
  facet_to_cell = reindex(oldfacet_to_cell, facet_to_oldfacet)
  facet_to_lfacet = reindex(oldfacet_to_lfacet,facet_to_oldfacet)
  cellgrid = Grid(model,d)
  cell_to_extrussion = celltypes(cellgrid)
  cell_to_polytope = _cell_to_polytope(cell_to_extrussion)
  trian = Triangulation(cellgrid)
  phi = CellGeomap(trian)
  descriptor = BoundaryDescriptor(
    facet_to_cell,facet_to_lfacet,cell_to_polytope,phi)
  BoundaryGrid(grid,descriptor)
end

function BoundaryTriangulation(model::DiscreteModel,tags,icell::Int=1)
  grid = BoundaryGrid(model,tags,icell)
  Triangulation(grid)
end

function BoundaryTriangulation(model::DiscreteModel)
  grid = BoundaryGrid(model)
  Triangulation(grid)
end

function _setup_mask!(facet_to_mask,facet_to_label,tag_to_labels,tags)
  nfacets = length(facet_to_mask)
  for facet in 1:nfacets
    for tag in tags
      for label in tag_to_labels[tag]
        if facet_to_label[facet] == label
          facet_to_mask[facet] = true
        end
      end
    end
  end
end

function _cell_to_polytope(cell_to_extrussion)
  @notimplemented
end

function _cell_to_polytope(cell_to_extrussion::ConstantCellValue)
  extrussion = cell_to_extrussion.value
  l = cell_to_extrussion.length
  poly = Polytope(extrussion)
  ConstantCellValue(poly,l)
end

_setup_tags(facelabels,tags) = tags

_setup_tags(facelabels,tag::Integer) = [tag,]

function _setup_tags(facelabels,name::String)
  _setup_tags(facelabels,[name,])
end

function _setup_tags(facelabels,names::Vector{String})
  [ tag_from_name(facelabels,s) for s in names ]
end

end # module


text:
using BinaryBuilder

name = "Rmath"
version = v"0.2.2"

sources = [
    "https://github.com/JuliaStats/Rmath-julia/archive/v$(version).tar.gz" =>
        "6544f40e51999469873b0f28d4bdeecdc847d4b24250a65027ae07e7dccb9ccd",
]

script = raw"""
cd $WORKSPACE/srcdir/Rmath-julia-*
make
mkdir -p "${libdir}"
mv src/libRmath-julia.* "${libdir}"
"""

platforms = supported_platforms()

products = [
    LibraryProduct("libRmath-julia", :libRmath),
]

dependencies = []

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: tbeason/Yggdrasil

max_issues_repo_name: tbeason/Yggdrasil

max_forks_repo_name: tbeason/Yggdrasil


text:
function initialize!(p::GridVisualizer,::Type{VTKViewType})
    pctx=p.context
    VTKView=pctx[:Plotter]
    frame=VTKView.StaticFrame()
    pctx[:frame]=frame
    if pctx[:clear]
        VTKView.clear!(frame)
    end
    layout=pctx[:layout]
    tlayout=(layout[2],layout[1])
    for I in CartesianIndices(layout)
        ctx=p.subplots[I]
        ctx[:frame]=frame
    end
    VTKView.layout!(frame,tlayout...)
    VTKView.size!(frame,pctx[:resolution]...)
    pctx
end

function save(fname,p,::Type{VTKViewType})
    VTKView=p.context[:Plotter]
    base,ext=splitext(fname)
    if ext!=".png"
        error("VTKView can only save png files")
    end
    VTKView.writepng(p.context[:frame],fname)
end

function save(fname,scene,VTKView,::Type{VTKViewType})
    base,ext=splitext(fname)
    if ext!=".png"
        error("VTKView can only save png files")
    end
    VTKView.writepng(scene,fname)
end

function reveal(p::GridVisualizer,::Type{VTKViewType})
    VTKView=p.Plotter
    VTKView.display(p.context[:frame])
    p.context[:frame]
end

function reveal(ctx::SubVisualizer,TP::Type{VTKViewType})
    if ctx[:show]||ctx[:reveal]
        reveal(ctx[:GridVisualizer],TP)
    end
end

function gridplot!(ctx, TP::Type{VTKViewType},grid)
    VTKView=ctx[:Plotter]
    frame=ctx[:frame]
    if !haskey(ctx,:dataset)
        dataset=VTKView.DataSet()
        ctx[:dataset]=dataset
    end
    VTKView.simplexgrid!(ctx[:dataset],grid[Coordinates],grid[CellNodes])
    VTKView.boundarygrid!(ctx[:dataset],grid[BFaceNodes])
    VTKView.boundarymarker!(ctx[:dataset],grid[BFaceRegions])
    VTKView.cellmarker!(ctx[:dataset],grid[CellRegions])
    if !haskey(ctx,:gridview)
        ctx[:gridview]=VTKView.GridView()
        VTKView.data!(ctx[:gridview],ctx[:dataset])
        VTKView.addview!(frame,ctx[:gridview],ctx[:iplot]...)
    end
    reveal(ctx,TP)
end

gridplot!(ctx, T::Type{VTKViewType}, ::Type{Val{2}},grid)=gridplot!(ctx, T,grid)
gridplot!(ctx, T::Type{VTKViewType}, ::Type{Val{3}},grid)=gridplot!(ctx, T,grid)

function scalarplot!(ctx, TP::Type{VTKViewType},grid,func)
    VTKView=ctx[:Plotter]
    frame=ctx[:frame]
    if !haskey(ctx,:dataset)
        ctx[:dataset]=VTKView.DataSet()
    end
    if !haskey(ctx,:grid)
        ctx[:grid]=grid
        VTKView.simplexgrid!(ctx[:dataset],grid[Coordinates],grid[CellNodes])
    end
    if !haskey(ctx,:scalarview)
        ctx[:scalarview]=VTKView.ScalarView()
        VTKView.addview!(frame,ctx[:scalarview],ctx[:iplot])
        VTKView.data!(ctx[:scalarview],ctx[:dataset],ctx[:label])
    end
    if !seemingly_equal(grid,ctx[:grid])
        VTKView.simplexgrid!(ctx[:dataset],grid[Coordinates],grid[CellNodes])
    end
    VTKView.pointscalar!(ctx[:dataset],func,ctx[:label])
    if dim_space(grid)==3
        VTKView.isolevels!(ctx[:scalarview],[ctx[:flevel]])
        VTKView.show_isosurfaces!(ctx[:scalarview],true)
    end
    reveal(ctx,TP)
end

scalarplot!(ctx, T::Type{VTKViewType}, ::Type{Val{2}},grid, func)=scalarplot!(ctx, T,grid,func)
scalarplot!(ctx, T::Type{VTKViewType}, ::Type{Val{3}},grid, func)=scalarplot!(ctx, T,grid,func)
gridplot!(ctx, T::Type{VTKViewType}, ::Type{Val{1}},grid)=nothing

function scalarplot!(ctx, TP::Type{VTKViewType}, ::Type{Val{1}},grid, func)
    VTKView=ctx[:Plotter]
    frame=ctx[:frame]
    if !haskey(ctx,:plot)
        ctx[:plot]=VTKView.XYPlot()
        VTKView.addview!(frame,ctx[:plot],ctx[:iplot])
        VTKView.xrange!(ctx[:plot],ctx[:xlimits]...)
        VTKView.yrange!(ctx[:plot],ctx[:limits]...)
        VTKView.linewidth!(ctx[:plot],1)
    end
    if ctx[:clear]
        VTKView.clear!(ctx[:plot])
    end
    plot=ctx[:plot]
    VTKView.plotcolor!(plot,rgbtuple(ctx[:color])...)

    VTKView.title!(plot,ctx[:title])

    if ctx[:label]!=""
        VTKView.plotlegend!(plot,ctx[:label])
    end

    VTKView.addplot!(plot,collect(grid[Coordinates][1,:]),collect(func))
    reveal(ctx,TP)
end


text:
function  valueToString( v )
	#=
	if typeof(v) <: String
		vstr = [string("\"",v,"\"")];
	elseif typeof(v) <: Number
		if any(size(v)~=[1 1])
			vstr = '[';
			for i=1:size(v,1)
				for j=1:size(v,2)
					vstr = [vstr num2str(v(i,j)) ' '];
				end
				vstr = [vstr '; '];
			end
			vstr = [vstr ']'];
		else
			vstr = num2str(v);
		end
	elseif islogical(v)
		if v
			vstr = 'true';
		else
			vstr = 'false';
		end
	elseif isstruct(v)
		vstr = 'struct(';
		fields = fieldnames(v);
		for i = 1:numel(fields)
			vstr = [vstr '''' fields{i} ''', '];
			w = v.(fields{i});
			if iscell(w)
				error('A cell cannot have a cell as a value, since it messes up the whole struct.');
			end
			vstr = [vstr valueToString(w)];
			if i<numel(fields)
				vstr = [vstr ', '];
			end
		end
		vstr = [vstr ')'];
	elseif iscell(v)
		vstr = '{';
		for i = 1:length(v)
			vstr = [vstr valueToString(v{i})];
			if i<length(v)
				vstr = [vstr ', '];
			end
		end
		vstr = [vstr '}'];
	else
		error(['Add support for this type of object before passing it to valueToString:' class(v)]);
	end
	=#
	#Julia: simplified version, maybe need to specify
	vstr = string(v)
	return vstr
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
# prior factor that passes a density belief straight through to inference without resampling

export PartialPriorPassThrough, PackedPartialPriorPassThrough

struct PartialPriorPassThrough{B <: Union{<:HeatmapGridDensity,<:LevelSetGridNormal}, T <:Tuple} <: AbstractPrior
  Z::B
  partial::T
end

getManifold(pppt::PartialPriorPassThrough) = getManifold(pppt.Z)

# this step is skipped during main inference process
getSample(cf::CalcFactor{<:PartialPriorPassThrough}) = sampleTangent(getManifold(cf.factor), cf.factor.Z)

## ====================================================================================================
## Serialize PartialPriorPassThrough
## ====================================================================================================

"""
    $TYPEDEF

Required internal density to store its type
"""
mutable struct PackedPartialPriorPassThrough <: AbstractPackedFactor
  Z::String # PackedHeatmapGridDensity
  partial::Vector{Int}
end

function convert( ::Union{Type{<:AbstractPackedFactor}, Type{<:PackedPartialPriorPassThrough}},
                  obj::PartialPriorPassThrough )
  #

  po = convert(PackedSamplableBelief, obj.Z)
  str = JSON2.write(po)
  PackedPartialPriorPassThrough(str, Int[obj.partial...])
end

function convert( ::Union{Type{<:AbstractFactor}, Type{<:PartialPriorPassThrough}},
                  obj::PackedPartialPriorPassThrough )
  #

  # get as bland obj to extract type
  _up = JSON2.read(obj.Z)
  _typ = DFG.getTypeFromSerializationModule(_up._type)
  # now redo with type
  pdens = JSON2.read(obj.Z, _typ)
  dens = convert(SamplableBelief, pdens)
  PartialPriorPassThrough(dens, tuple(obj.partial...))
end

#


text:
module Errors

# use Cenum for error codes because it allows overriding show, unlike Base.Enum
using CEnum

using ..LibPQ: Connection, Result, libpq_c, error_message, error_field

"Base abstract type for all custom exceptions thrown by LibPQ.jl"
abstract type LibPQException <: Exception end

"An exception with an error message generated by PostgreSQL"
abstract type PostgreSQLException <: LibPQException end

# PostgreSQL errors have trailing newlines
# https://www.postgresql.org/docs/10/libpq-status.html#LIBPQ-PQERRORMESSAGE
Base.showerror(io::IO, err::PostgreSQLException) = print(io, chomp(err.msg))

"An exception generated by LibPQ.jl"
abstract type JLClientException <: LibPQException end

"An error regarding a connection reported by PostgreSQL"
struct PQConnectionError <: PostgreSQLException
    msg::String
end

function PQConnectionError(jl_conn::Connection)
    return PQConnectionError(error_message(jl_conn))
end

"An error from parsing connection parameter strings reported by PostgreSQL"
struct ConninfoParseError <: PostgreSQLException
    msg::String
end

"An error regarding a connection generated by LibPQ.jl"
struct JLConnectionError <: JLClientException
    msg::String
end

"An error regarding a query result generated by LibPQ.jl"
struct JLResultError <: JLClientException
    msg::String
end

"""
An error regarding a query result generated by PostgreSQL

The `Code` parameter represents the PostgreSQL error code as defined in
[Appendix A. PostgreSQL Error Codes](https://www.postgresql.org/docs/devel/errcodes-appendix.html).
The `Class` parameter is the first two characters of that code, also listed on that page.

For a list of all error aliases, see `src/error_codes.jl`, which was generated using the
PostgreSQL documentation linked above.

```jldoctest
julia> try execute(conn, "SELORCT NUUL;") catch err println(err) end
LibPQ.Errors.SyntaxError("ERROR:  syntax error at or near \\"SELORCT\\"\\nLINE 1: SELORCT NUUL;\\n        ^\\n")

julia> LibPQ.Errors.SyntaxError
LibPQ.Errors.PQResultError{LibPQ.Errors.C42,LibPQ.Errors.E42601}
```
"""
struct PQResultError{Class, Code} <: PostgreSQLException
    msg::String
    verbose_msg::Union{String, Nothing}

    function PQResultError{Class_, Code_}(msg, verbose_msg) where {Class_, Code_}
        return new{Class_::Errors.Class, Code_::Errors.ErrorCode}(
            convert(String, msg),
            convert(Union{String, Nothing}, verbose_msg),
        )
    end
end

include("error_codes.jl")

# avoid exposing the meaningless integer value of the enum
function Base.show(io::IO, ::MIME"text/plain", class::Class)
    print(io, class, "::", typeof(class))
end
function Base.show(io::IO, ::MIME"text/plain", code::ErrorCode)
    print(io, code, "::", typeof(code))
end

Base.parse(::Type{Class}, str::AbstractString) = getfield(Errors, Symbol("C", str))
Base.parse(::Type{ErrorCode}, str::AbstractString) = getfield(Errors, Symbol("E", str))

function PQResultError{Class, Code}(msg::String) where {Class, Code}
    return PQResultError{Class, Code}(msg, nothing)
end

function PQResultError(result::Result; verbose=false)
    msg = error_message(result; verbose=false)
    verbose_msg = verbose ? error_message(result; verbose=true) : nothing
    code_str = something(error_field(result, libpq_c.PG_DIAG_SQLSTATE), "UNOWN")
    class = parse(Class, code_str[1:2])
    code = parse(ErrorCode, code_str)

    return PQResultError{class, code}(msg, verbose_msg)
end

error_class(err::PQResultError{Class_}) where {Class_} = Class_::Class
error_code(err::PQResultError{Class_, Code_}) where {Class_, Code_} = Code_::ErrorCode

function Base.showerror(io::IO, err::T) where T <: PQResultError
    msg = err.verbose_msg === nothing ? err.msg : err.verbose_msg

    print(io, ERROR_NAMES[T], ": ", chomp(msg))
end

function Base.show(io::IO, err::T) where T <: PQResultError
    print(io, "LibPQ.Errors.", ERROR_NAMES[T], '(', repr(err.msg))

    if err.verbose_msg !== nothing
        print(io, ", ", repr(err.verbose_msg))
    end

    print(io, ')')
end

end


text:
###############################################################################
# General info, including the Git version
function getInfo(features, pkg_path)
  cd(pkg_path)
  gitsha = ""
  try
    gitsha = readall(`git log -1 --format="%H"`)
  catch
    # NOP
  end
  features[:GITSHA] = gitsha
end

###############################################################################
# REQUIRE file (the one in the root directory of the package itself)
function checkREQUIRE(features, pkg_path)
  REQUIRE_path = joinpath(pkg_path, "REQUIRE")
  if isfile(REQUIRE_path)
    # Exists
    features[:REQUIRE_EXISTS] = true
    # Does it specifiy a Julia dependency?
    grep_result = readall(REQUIRE_path)
    features[:REQUIRE_VERSION] = ismatch(r"julia", grep_result)
  else
    # Does not exist
    features[:REQUIRE_EXISTS] = false
    features[:REQUIRE_VERSION] = false
  end
end

###############################################################################
# License file
function checkLicense(features, pkg_path)

  features[:LICENSE_EXISTS] = false
  features[:LICENSE] = "Unknown"
  features[:LICENSE_FILE] = ""

  # Test for some sort of license file first
  possible_files = ["LICENSE",
                    "LICENSE.md",
                    "License.md",
                    "LICENSE.txt",
                    "README",
                    "README.md",
                    "README.txt",
                    "COPYING",
                    "COPYING.md",
                    "COPYING.txt"]
  for filename in possible_files
    fullfilename = joinpath(pkg_path, filename)
    if isfile(fullfilename)
      if guessLicense(features, fullfilename)
        # Stop once we identify a license
        features[:LICENSE_FILE] = filename
        return
      end
    end
  end
end

function guessLicense(features, filename)
  # Be greedy, read the whole file
  text = lowercase(readall(filename))

  # MIT License
  if ismatch(r"mit license", text) ||
     ismatch(r"mit expat license", text) ||
     ismatch(r"mit \"expat\" license", text) ||
     ismatch(r"permission is hereby granted, free of charge, to any person", text)
    features[:LICENSE_EXISTS] = true
    features[:LICENSE] = "MIT"
  # GPL v2
  elseif ismatch(r"gpl version 2", text) ||
         ismatch(r"gnu general public license\s+version 2", text) ||
         ismatch(r"gnu general public license, version 2", text)
    features[:LICENSE_EXISTS] = true
    features[:LICENSE] = "GPL v2"
  # GPL v3
  elseif ismatch(r"gpl version 3", text) ||
         ismatch(r"http://www.gnu.org/licenses/gpl-3.0.txt", text) ||
         ismatch(r"gnu general public license\s+version 3", text) ||
         ismatch(r"gpl v3", text)
    features[:LICENSE_EXISTS] = true
    features[:LICENSE] = "GPL v3"
  # LGPL v2.1
  elseif ismatch(r"lgpl version 2.1", text) ||
         ismatch(r"gnu lesser general public license\s+version 2\.1", text)
    features[:LICENSE_EXISTS] = true
    features[:LICENSE] = "LGPL v2.1"
  # LGPL v3.0
  elseif ismatch(r"lgpl-3.0", text)
    features[:LICENSE_EXISTS] = true
    features[:LICENSE] = "LGPL v3.0"
  # BSD
  elseif ismatch(r"bsd", text)
    features[:LICENSE_EXISTS] = true
    features[:LICENSE] = "BSD"
  # GNU Affero
  elseif ismatch(r"gnu affero general public license", text)
    features[:LICENSE_EXISTS] = true
    features[:LICENSE] = "GNU Affero"
  # Romantic WTF public license (!!) (see TOML.jl)
  elseif ismatch(r"romantic wtf public license", text)
    features[:LICENSE_EXISTS] = true
    features[:LICENSE] = "Romantic WTF"
  # No license identified
  else
    return false
  end

  return true
end

###############################################################################
# Testing folder/files
function checkTesting(features, pkg_path, pkg_name)
  features[:TEST_EXIST] = false

  # First of all, check for the existence of a test folder with Julia
  # files in it
  possible_folders = [joinpath(pkg_path, "test"),
                      joinpath(pkg_path, "tests")]
  for folder in possible_folders
    try
      # ls the folder
      file_list = split(readall(`ls $folder`), "\n")
      # if that didn't die, look through the files
      for file in file_list
        # Any Julia files in there?
        if file[(end-2):(end)] == ".jl"
          features[:TEST_EXIST] = true
          break
        end
      end
    catch
      # No test folder at all
    end
  end

  # Look for a master test file
  possible_files = [
    joinpath(pkg_path, "test",  "runtests.jl"),
    joinpath(pkg_path, "tests", "runtests.jl"),
    joinpath(pkg_path,          "runtests.jl"),
    joinpath(pkg_path, "test",  "run_tests.jl"),
    joinpath(pkg_path, "tests", "run_tests.jl"),
    joinpath(pkg_path,          "run_tests.jl"),
    joinpath(pkg_path, "test",  "tests.jl"),
    joinpath(pkg_path, "tests", "tests.jl"),
    joinpath(pkg_path,          "tests.jl"),
    joinpath(pkg_path, "test",  "test.jl" ),
    joinpath(pkg_path, "tests", "test.jl" ),
    joinpath(pkg_path,          "test.jl" )]
  features[:TEST_MASTERFILE] = ""
  for filename in possible_files
    if isfile(filename)
      features[:TEST_MASTERFILE] = filename
      features[:TEST_EXIST] = true
      break
    end
  end

  # Do they have a .travis.yml file?
  travis_file = joinpath(pkg_path, ".travis.yml")
  features[:TEST_TRAVIS] = isfile(travis_file)

  features[:TEST_STATUS] = ""

  # Are tests even meaningful?
  features[:TEST_POSSIBLE] = true
  features[:TEST_POSSIBLE] &= !(pkg_name == "Arduino")      # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "Clang")        # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "CPLEX")        # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "CLFFT")        # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "CUDA")         # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "FITSIO")       # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "Gaston")       # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "Gtk")          # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "Gurobi")       # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "ImageView")    # Reason: Tk
  features[:TEST_POSSIBLE] &= !(pkg_name == "MATLAB")       # Reason: MATLAB
  features[:TEST_POSSIBLE] &= !(pkg_name == "MATLABCluster") # Reason: MATLAB
  features[:TEST_POSSIBLE] &= !(pkg_name == "Mongo")        # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "Mongrel2")     # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "Mosek")        # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "OpenCL")       # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "OpenGL")       # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "Pandas")       # Reason: python
  features[:TEST_POSSIBLE] &= !(pkg_name == "ProfileView")  # Reason: Tk
  features[:TEST_POSSIBLE] &= !(pkg_name == "PyCall")       # Reason: python
  features[:TEST_POSSIBLE] &= !(pkg_name == "PyLexYacc")    # Reason: python
  features[:TEST_POSSIBLE] &= !(pkg_name == "PyPlot")       # Reason: python
  features[:TEST_POSSIBLE] &= !(pkg_name == "PySide")       # Reason: python
  features[:TEST_POSSIBLE] &= !(pkg_name == "SemidefiniteProgramming") # Reason: binaries
  features[:TEST_POSSIBLE] &= !(pkg_name == "Tk")           # Reason: Tk
  features[:TEST_POSSIBLE] &= !(pkg_name == "Winston")      # Reason: Tk
  features[:TEST_POSSIBLE] &= !(pkg_name == "Vega")         # Reason: weird build

  if !features[:TEST_POSSIBLE]
    features[:TEST_STATUS] = "not_possible"
    return
  end

  if features[:TEST_MASTERFILE] == ""
    # Couldn't find a master test file
    # See if "using" works
    testoutput = ""
    try
      fp = open("testusing.jl","w")
      write(fp, "using $pkg_name\n")
      close(fp)
      testoutput = readall(`julia testusing.jl`)
      features[:TEST_STATUS] = "using_pass"
    catch
      # Didn't run without errors
      features[:TEST_STATUS] = "using_fail"
    end
  else
    # Found a master test file
    # Run it to see if it works
    # Run the tests to see how they go
    testoutput = ""
    try
      # Move into the package directory in case tests rely on that
      # First saw in ASCIIPlots
      curdir = strip(readall(`pwd`))
      cd(Pkg.dir(pkg_name))
      testoutput = readall(`julia $(features[:TEST_MASTERFILE])`)
      cd(curdir)
      features[:TEST_STATUS] = "full_pass"
    catch
      # Has tests, and they failed
      features[:TEST_STATUS] = "full_fail"
    end
  end

end


text:
using UnPack: @unpack
using MuladdMacro: @muladd
using StaticArrays: SMatrix
using OrdinaryDiffEq: OrdinaryDiffEq
using LinearAlgebra: Hermitian, tr, Eigen, dot

using NQCDynamics.DynamicsMethods: MappingVariableMethods
using NQCModels: nstates, NQCModels

"""
    MInt <: OrdinaryDiffEq.OrdinaryDiffEqAlgorithm

Second order symplectic momentum integral algorithm.

# Reference

[J. Chem. Phys. 148, 102326 (2018)](https://doi.org/10.1063/1.5005557)
"""
struct MInt <: OrdinaryDiffEq.OrdinaryDiffEqAlgorithm end

OrdinaryDiffEq.isfsal(::MInt) = false

mutable struct MIntCache{uType,T} <: OrdinaryDiffEq.OrdinaryDiffEqMutableCache
    u::uType
    uprev::uType
    tmp::uType
    C::Matrix{T}
    D::Matrix{T}
    Γ::Matrix{T}
    Ξ::Matrix{T}
    tmp_mat::Matrix{T}
    tmp_vec1::Vector{T}
    tmp_vec2::Vector{T}
end

function OrdinaryDiffEq.alg_cache(::MInt,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
    tmp = zero(u)
    n = NQCModels.nstates(p.calculator.model)
    C = zeros(n,n)
    D = zeros(n,n)
    Γ = zeros(n,n)
    Ξ = zeros(n,n)
    tmp_mat = zeros(n,n)
    tmp_vec1 = zeros(n)
    tmp_vec2 = zeros(n)
    MIntCache(u, uprev, tmp, C, D, Γ, Ξ, tmp_mat, tmp_vec1, tmp_vec2)
end

function OrdinaryDiffEq.initialize!(_, ::MIntCache) end

@muladd function OrdinaryDiffEq.perform_step!(integrator, cache::MIntCache, repeat_step=false)
    @unpack dt,uprev,u,p = integrator
    @unpack tmp, Γ, Ξ = cache

    calc = p.calculator

    rtmp = DynamicsUtils.get_positions(tmp)
    vprev = DynamicsUtils.get_velocities(uprev)
    rprev = DynamicsUtils.get_positions(uprev)
    v = DynamicsUtils.get_velocities(u)
    r = DynamicsUtils.get_positions(u)
    X = MappingVariableMethods.get_mapping_positions(u)
    P = MappingVariableMethods.get_mapping_momenta(u)

    copyto!(v, vprev)
    copyto!(r, rprev)

    step_A!(rtmp, rprev, dt/2, v)

    propagate_mapping_variables!(cache, calc, X, P, rtmp, dt)

    eigs = Calculators.get_eigen(calc, rtmp)
    adiabatic_derivative = Calculators.get_adiabatic_derivative(calc, rtmp)
    state_independent_derivative = zero(rtmp)
    NQCModels.state_independent_derivative!(calc.model, state_independent_derivative, rtmp)
    X = MappingVariableMethods.get_mapping_positions(u)
    P = MappingVariableMethods.get_mapping_momenta(u)
    ∂V = Calculators.get_derivative(calc, rtmp)
    for i=1:natoms(p)
        for j=1:ndofs(p)

            set_gamma!(Γ, adiabatic_derivative[j,i], eigs.values, dt)
            transform_matrix!(Γ, eigs.vectors, cache.tmp_mat)
            E = Γ

            set_xi!(Ξ, adiabatic_derivative[j,i], eigs.values, dt)
            transform_matrix!(Ξ, eigs.vectors, cache.tmp_mat)
            F = Ξ

            force = get_mapping_nuclear_force(X, P, E, F, cache.tmp_vec1)
            v[j,i] -= force / p.atoms.masses[i]
            v[j,i] -= state_independent_derivative[j,i] / p.atoms.masses[i] * dt
            v[j,i] += tr(∂V[j,i]) * p.method.γ / p.atoms.masses[i] * dt / 2
        end
    end

    step_A!(r, rtmp, dt/2, v)

    return nothing
end

function propagate_mapping_variables!(cache, calc, X, P, rtmp, dt)

    @unpack C, D, tmp_vec1, tmp_vec2 = cache

    eigen = Calculators.get_eigen(calc, rtmp)

    set_C_propagator!(C, cache, eigen, dt)
    set_D_propagator!(D, cache, eigen, dt)

    # tmp_vec1 = C*X - D*P
    mul!(tmp_vec1, C, X)
    mul!(tmp_vec1, D, P, -1.0, 1.0)

    # tmp_vec2 = C*P + D*X
    mul!(tmp_vec2, C, P)
    mul!(tmp_vec2, D, X, 1.0, 1.0)

    copyto!(X, tmp_vec1)
    copyto!(P, tmp_vec2)
end

"Get the `C` propagator for the mapping variables."
function set_C_propagator!(C, cache, eigen::Eigen, dt::Real)
    fill!(C, zero(eltype(C)))
    for i in axes(C,1)
        C[i,i] = cos(eigen.values[i] * dt)
    end
    transform_matrix!(C, eigen.vectors, cache.tmp_mat)
end

"Get the `D` propagator for the mapping variables."
function set_D_propagator!(D, cache, eigen::Eigen, dt::Real)
    fill!(D, zero(eltype(D)))
    for i in axes(D,1)
        D[i,i] = sin(-eigen.values[i] * dt)
    end
    transform_matrix!(D, eigen.vectors, cache.tmp_mat)
end

function transform_matrix!(M, transform, tmp_mat)
    mul!(tmp_mat, M, transform')
    mul!(M, transform, tmp_mat)
    return nothing
end

"Get the `Γ` variable used to calculate the nuclear propagators."
function set_gamma!(Γ::AbstractMatrix, W::AbstractMatrix, Λ::AbstractVector, dt::Real)
    for i in axes(Γ, 2)
        Γ[i,i] = W[i,i] * dt
        for j in axes(Γ, 1)
            if j != i
                Γ[j,i] = sin((Λ[i] - Λ[j])*dt) * W[j,i] / (Λ[i] - Λ[j])
            end
        end
    end
    return Γ
end

"Get the `Ξ` variable used to calculate the nuclear propagators."
function set_xi!(Ξ::AbstractMatrix, W::AbstractMatrix, Λ::AbstractVector, dt::Real)
    for i in axes(Ξ, 2)
        for j in axes(Ξ, 1)
            if j != i
                Ξ[j,i] = (1 - cos((Λ[i] - Λ[j])*dt)) * W[j,i] / (Λ[i] - Λ[j])
            end
        end
    end
    return Ξ
end

"""
Get the force due to the mapping variables.

Equivalent to this but doesn't allocate:
    return 0.5 * (q'*E*q + p'*E*p) - q'*F*p
"""
function get_mapping_nuclear_force(q::AbstractVector, p::AbstractVector,
                           E::AbstractMatrix, F::AbstractMatrix, tmp_vec)
    force = zero(eltype(q))
    mul!(tmp_vec, E, q)
    force += dot(q, tmp_vec)
    mul!(tmp_vec, E, p)
    force += dot(p, tmp_vec)
    mul!(tmp_vec, F, p)
    force -= 2dot(q, tmp_vec)
    return force / 2
end

DynamicsMethods.select_algorithm(::Simulation{<:MappingVariableMethods.SpinMappingW}) = MInt()


text:
using Test
using CompScienceMeshes
using StaticArrays

import BEAST
BE = BEAST

Γ = readmesh(joinpath(dirname(@__FILE__),"assets","sphere2.in"))
nc = numcells(Γ)
t = chart(Γ, Γ.faces[1])
s = chart(Γ, Γ.faces[nc])

X = BEAST.raviartthomas(Γ)
x = BEAST.refspace(X)

κ = 1.0
op = BEAST.MWSingleLayer3D(κ)

n = BE.numfunctions(x)
z1 = zeros(ComplexF64, n, n)
z2 = zeros(z1)

tqd = BE.quadpoints(x, [t], (12,13))
bqd = BE.quadpoints(x, [t], (13,))

SE_strategy = BE.WiltonSEStrategy(
  tqd[2,1],
  BE.DoubleQuadStrategy(
	tqd[1,1],
	bqd[1,1],
  ),
)
BEAST.momintegrals!(op, x, x, t, t, z1, SE_strategy)

EE_strategy = BEAST.BogaertSelfPatchStrategy(20)
BEAST.momintegrals!(op, x, x, t, t, z2, EE_strategy)

Γ = meshrectangle(1.0, 1.0, 1.0)
t = chart(Γ, Γ.faces[1])
s = chart(Γ, Γ.faces[2])

z3 = zeros(ComplexF64, n, n)
EE_strategy = BEAST.BogaertEdgePatchStrategy(13,30)
BEAST.momintegrals!(op, x, x, t, s, z3, EE_strategy)

z4 = zeros(ComplexF64, n, n)
tqd = BE.quadpoints(x, [t], (12,13))
bqd = BE.quadpoints(x, [s], (13,))
SE_strategy = BE.WiltonSEStrategy(
  tqd[2,1],
  BE.DoubleQuadStrategy(
	tqd[1,1],
	bqd[1,1],
  ),
)
BEAST.momintegrals!(op, x, x, t, s, z4, SE_strategy)

p = [
    point(0.0, 0.0, 0.0),
    point(0.0, 1.0, 0.0),
    point(1.0, 1.0, 0.0),
    point(0.0, -1.0, 0.0),
    point(1.0, -1.0, 0.0)]

t = simplex(p[1], p[2], p[3])
s = simplex(p[1], p[4], p[5])

z5 = zeros(ComplexF64, n, n)
EE_strategy = BEAST.BogaertPointPatchStrategy(9,10)
BEAST.momintegrals!(op, x, x, t, s, z5, EE_strategy)

z6 = zeros(ComplexF64, n, n)
tqd = BE.quadpoints(x, [t], (12,13))
bqd = BE.quadpoints(x, [s], (13,))
SE_strategy = BE.WiltonSEStrategy(
  tqd[2,1],
  BE.DoubleQuadStrategy(
	tqd[1,1],
	bqd[1,1],
  ),
)
BEAST.momintegrals!(op, x, x, t, s, z6, SE_strategy)

@test norm(z1-z2)/norm(z1) < 1.0e-6
@test norm(z3-z4)/norm(z3) < 3.0e-6
@test norm(z5-z6)/norm(z5) < 1.0e-8

@show norm(z1-z2)/norm(z1)
@show norm(z3-z4)/norm(z3)
@show norm(z5-z6)/norm(z5)

# repeat the test for the MFIE interactions
κ = 1.0
op = BEAST.MWDoubleLayer3D(0.0)

p = [
    point(0.0, 0.0, 0.0),
    point(1.0, 0.0, 0.0),
    point(1.0, 1.0, 0.0),
    point(0.0, 1.0, 1.0),
    point(1.0,-1.0, 0.0)
]

t1 = simplex(p[1], p[2], p[3])
t2 = simplex(p[1], p[3], p[4])
t3 = simplex(p[1], p[2], p[5])

# test the edge patch integral
z1 = zeros(ComplexF64, BE.numfunctions(x), BE.numfunctions(x))
z2 = zeros(ComplexF64, BE.numfunctions(x), BE.numfunctions(x))
z3 = zeros(ComplexF64, BE.numfunctions(x), BE.numfunctions(x))

s1 = BEAST.BogaertEdgePatchStrategy(13, 30)
tqd = BE.quadpoints(x, [t1], (12,13))
bqd = BE.quadpoints(x, [t2], (13,))
s2 = BE.WiltonSEStrategy(
  tqd[2,1],
  BE.DoubleQuadStrategy(
	tqd[1,1],
	bqd[1,1],
  ),
)
s3 = BE.DoubleQuadStrategy(
  tqd[1,1],
  tqd[1,1],
)

BEAST.momintegrals!(op, x, x, t1, t2, z1, s1)
BEAST.momintegrals!(op, x, x, t1, t2, z2, s2)
BEAST.momintegrals!(op, x, x, t1, t2, z3, s3)

@test norm(z2-z1)/norm(z2) < 1.0e-3

s1 = BEAST.BogaertPointPatchStrategy(9,10)
tqd = BE.quadpoints(x, [t2], (12,13))
bqd = BE.quadpoints(x, [t3], (13,))
s2 = BE.WiltonSEStrategy(
  tqd[2,1],
  BE.DoubleQuadStrategy(
	tqd[1,1],
	bqd[1,1],
  ),
)
s3 = BE.DoubleQuadStrategy(
  tqd[1,1],
  tqd[1,1],
)

fill!(z1, 0); BEAST.momintegrals!(op, x, x, t2, t3, z1, s1)
fill!(z2, 0); BEAST.momintegrals!(op, x, x, t2, t3, z2, s2)
fill!(z3, 0); BEAST.momintegrals!(op, x, x, t2, t3, z3, s3)

@test norm(z2-z1)/norm(z2) < 2.0e-6

s1 = BEAST.BogaertSelfPatchStrategy(20)
tqd = BE.quadpoints(x, [t1], (12,13))
bqd = BE.quadpoints(x, [t1], (13,))
s2 = BE.WiltonSEStrategy(
  tqd[2,1],
  BE.DoubleQuadStrategy(
	tqd[1,1],
	bqd[1,1],
  ),
)
s3 = BE.DoubleQuadStrategy(
  tqd[1,1],
  tqd[1,1],
)

fill!(z1, 0); BEAST.momintegrals!(op, x, x, t1, t1, z1, s1)
fill!(z2, 0); BEAST.momintegrals!(op, x, x, t1, t1, z2, s2)
fill!(z3, 0); BEAST.momintegrals!(op, x, x, t1, t1, z3, s3)

# Test if IBPP gives the correct value for ∫∫∇G
t2 = simplex(p[3], p[4], p[1])
t3 = simplex(p[2], p[5], p[1])

s1 = BEAST.BogaertPointPatchStrategy(9, 10)
_, I = BEAST.GetIntegrals(t2[1],t2[2],t2[3],t3[1],t3[2], 0.0, s1)
#I = I.GG

Rpnt = eltype(p)
Cplx = typeof(complex(zero(eltype(p[1]))))
Cpnt = similar_type(Rpnt, Cplx, Size(length(Rpnt),))

u, w = BEAST.trgauss(10)
J = zeros(Cpnt, 3, 3)
for i in 1:length(w)
    p = neighborhood(t2, u[:,i])
    x = cartesian(p)
    λ = barycentric(p)
    dx = w[i] * jacobian(p)
    for j in 1:length(w)
        q = neighborhood(t3, u[:,j])
        y = cartesian(q)
        μ = barycentric(q)
        dy = w[j] * jacobian(q)
        r = norm(x-y)
        for i in 1:3
            for j in 1:3
                J[i,j] += dx * dy / r^3 * (x-y) / 4π * λ[i] * μ[j]
            end
        end
    end
end

@test maximum([real(norm(x-y)) for (x,y) in zip(I,J)]) < 2.0e-6


text:
using JSON

clear_gpu = clear_gpu_memory(conn)
@test typeof(clear_gpu) == Nothing

v = String(JSON.read("test/data/omnisci_vega.json"))
veg = render_vega(conn, v, 1)
@test typeof(veg) == OmniSci.TRenderResult
# using ImageMagick
# ImageMagick.readblob(veg.image)

#Can't test this on Travis, use OmniSci Jenkins?
#sql_execute_gdf


text:
# This file defines the EarthEngine types

# list of main earthengine types
maintypes = [
    :Number,
    :String,
    :Image,
    :Feature,
    :ImageCollection,
    :FeatureCollection,
    :Date,
    :List,
    :Filter,
    :Geometry,
    :Dictionary,
]

# types that get defined on-the-fly during Initialization
psuedotypes = [
    :Array,
    :Blob,
    :DateRange,
    :Classifier,
    :Clusterer,
    :ConfusionMatrix,
    :ErrorMargin,
    :Join,
    :Kernel,
    :Model,
    :PixelType,
    :Projection,
    :Reducer,
]

# primitive types that other types inherit from
basetypes = [:Collection, :Element, :ComputedObject, :Function]

# concat types together in one vector
# order matters!!! basetypes have to be after everything else
types = vcat(psuedotypes, maintypes, basetypes)

# create the Julia types from Python objects
for type in types
    @eval begin
        @pytype $(type) () -> ee.$(type)
    end
end

# define helper constructor for empty types that wrap more cryptic constructors
# this allows users to simply call EE.Type() to get a type value without having to provide inputs
# empty constructors for base types
EarthEngine.ComputedObject() = EarthEngine.ComputedObject("", "")
EarthEngine.Element() = EarthEngine.Element("", "")
EarthEngine.Collection() = EarthEngine.Collection("", "")

# empty constructors for psuedo types
EarthEngine.Reducer() = EarthEngine.Reducer(EarthEngine.ComputedObject())
EarthEngine.Array() = EarthEngine.Array([])
EarthEngine.Blob() = EarthEngine.Blob("")
EarthEngine.DateRange() = EarthEngine.DateRange("", "")
EarthEngine.Classifier() = EarthEngine.Classifier(EarthEngine.ComputedObject())
EarthEngine.Clusterer() = EarthEngine.Clusterer(EarthEngine.ComputedObject())
EarthEngine.ConfusionMatrix() = EarthEngine.ConfusionMatrix("")
EarthEngine.ErrorMargin() = EarthEngine.ErrorMargin(NaN)
EarthEngine.Join() = EarthEngine.Join(EarthEngine.ComputedObject())
EarthEngine.Kernel() = EarthEngine.Kernel(EarthEngine.ComputedObject())
EarthEngine.Model() = EarthEngine.Model(EarthEngine.ComputedObject())
EarthEngine.PixelType() = EarthEngine.PixelType(EarthEngine.ComputedObject())
EarthEngine.Projection() = EarthEngine.Projection(EarthEngine.ComputedObject())

# empty constructors for main types
EarthEngine.Date() = EarthEngine.Date("")
EarthEngine.List() = EarthEngine.List([])
EarthEngine.Geometry() = EarthEngine.Geometry(Point(NaN, NaN))
EarthEngine.Feature() = EarthEngine.Feature(EarthEngine.Geometry())
EarthEngine.FeatureCollection() = EE.FeatureCollection([])
EarthEngine.ImageCollection() = EarthEngine.ImageCollection("")
EarthEngine.String() = EarthEngine.String("")
EarthEngine.Number() = EarthEngine.Number(NaN)


text:
"""
The `Header` type is a struct representing the essential metadata for a single
record in a tar file with this definition:

    struct Header
        path :: String # path relative to the root
        type :: Symbol # type indicator (see below)
        mode :: UInt16 # mode/permissions (best viewed in octal)
        size :: Int64  # size of record data in bytes
        link :: String # target path of a symlink
    end

Types are represented with the following symbols: `file`, `hardlink`, `symlink`,
`chardev`, `blockdev`, `directory`, `fifo`, or for unknown types, the typeflag
character as a symbol. Note that [`extract`](@ref) refuses to extract records
types other than `file`, `symlink` and `directory`; [`list`](@ref) will only
list other kinds of records if called with `strict=false`.

The tar format includes various other metadata about records, including user and
group IDs, user and group names, and timestamps. The `Tar` package, by design,
completely ignores these. When creating tar files, these fields are always set
to zero/empty. When reading tar files, these fields are ignored aside from
verifying header checksums for each header record for all fields.
"""
struct Header
    path::String
    type::Symbol
    mode::UInt16
    size::Int64
    link::String
end

function Base.show(io::IO, hdr::Header)
    show(io, Header)
    print(io, "(")
    show(io, hdr.path)
    print(io, ", ")
    show(io, hdr.type)
    print(io, ", 0o", string(hdr.mode, base=8, pad=3), ", ")
    show(io, hdr.size)
    print(io, ", ")
    show(io, hdr.link)
    print(io, ")")
end

const TYPE_SYMBOLS = (
    '0'  => :file,
    '\0' => :file, # legacy encoding
    '1'  => :hardlink,
    '2'  => :symlink,
    '3'  => :chardev,
    '4'  => :blockdev,
    '5'  => :directory,
    '6'  => :fifo,
)

function to_symbolic_type(type::Char)
    for (t, s) in TYPE_SYMBOLS
        type == t && return s
    end
    isascii(type) ||
        error("invalid type flag: $(repr(type))")
    return Symbol(type)
end

function from_symbolic_type(sym::Symbol)
    for (t, s) in TYPE_SYMBOLS
        sym == s && return t
    end
    str = String(sym)
    ncodeunits(str) == 1 && isascii(str[1]) ||
        throw(ArgumentError("invalid type symbol: $(repr(sym))"))
    return str[1]
end

function check_header(hdr::Header)
    errors = String[]
    err(e::String) = push!(errors, e)

    # error checks
    isempty(hdr.path) &&
        err("path is empty")
    0x0 in codeunits(hdr.path) &&
        err("path contains NUL bytes")
    0x0 in codeunits(hdr.link) &&
        err("link contains NUL bytes")
    !isempty(hdr.path) && hdr.path[1] == '/' &&
        err("path is absolute")
    occursin(r"(^|/)\.\.(/|$)", hdr.path) &&
        err("path contains '..' component")
    hdr.type in (:file, :symlink, :directory) ||
        err("unsupported entry type")
    hdr.type ∉ (:hardlink, :symlink) && !isempty(hdr.link) &&
        err("non-link with link path")
    hdr.type == :symlink && hdr.size != 0 &&
        err("symlink with non-zero size")
    hdr.type == :directory && hdr.size != 0 &&
        err("directory with non-zero size")
    hdr.type != :directory && endswith(hdr.path, "/") &&
       err("non-directory path ending with '/'")
    hdr.type != :directory && (hdr.path == "." || endswith(hdr.path, "/.")) &&
       err("non-directory path ending with '.' component")
    hdr.size < 0 &&
       err("negative file size")
    isempty(errors) && return

    # contruct error message
    if length(errors) == 1
        msg = errors[1] * "\n"
    else
        msg = "tar header with multiple errors:\n"
        for e in sort!(errors)
            msg *= " * $e\n"
        end
    end
    msg *= repr(hdr)
    error(msg)
end


text:
function activate(url, sha, location = ".projects")
    # directory where a package is going to be checkedout
    env = joinpath(@__DIR__, location, sha)

    # if the directory exists, assume it's been cloned, instantiated, and precompiled
    if (isdir(env))
        Pkg.activate(env)
    else
        # if not checkout, activate, instantiate, and precompile it
        mkpath(env) # make a dir
        LibGit2.clone(url, env)
        LibGit2.checkout!(LibGit2.GitRepo(env), sha, force = true)
        Pkg.activate(env)
        pkg"instantiate"
        pkg"precompile"
    end
end


text:
@testset "kernelsum" begin
    k1 = LinearKernel()
    k2 = SqExponentialKernel()
    k = KernelSum(k1, k2)
    @test k == KernelSum([k1, k2]) == KernelSum((k1, k2))
    @test length(k) == 2
    @test string(k) == (
        "Sum of 2 kernels:\n" *
        "\tLinear Kernel (c = 0.0)\n" *
        "\tSquared Exponential Kernel (metric = Euclidean(0.0))"
    )

    # Standardised tests.
    TestUtils.test_interface(k, Float64)
    test_ADs(x -> KernelSum(SqExponentialKernel(), LinearKernel(; c=exp(x[1]))), rand(1))

    test_params(k1 + k2, (k1, k2))
end


text:
"""
Broadcast a message to all connected web pages to be interpetted by WebSocket listener. For example, in JavaScript:

var sock = new WebSocket('ws://'+window.location.host);
sock.onmessage = function( message ){
    var msg = JSON.parse(message.data);
    console.log(msg);
}
"""
function broadcast(route,args::Dict)
    if haskey(pages,route)
        for (cid,client) in pages[route].sessions
            if isopen(client)
                write(client, json(args))
            end
        end
    end
end
broadcast(r,t,d) = broadcast(r,Dict("type"=>t,"data"=>d))

function broadcast(args::Dict)
    for key in keys(pages)
        for (cid,client) in pages[key].sessions
            if isopen(client)
                write(client, json(args))
            end
        end
    end
    for (cid,client) in external
        if isopen(client)
            write(client, json(args))
        end
    end
end
broadcast(t,d) = broadcast(Dict("type"=>t,"data"=>d))
"""
Send a message to the specified connection to be interpetted by WebSocket listener. For example, in JavaScript:

var sock = new WebSocket('ws://'+window.location.host);
sock.onmessage = function( message ){
    var msg = JSON.parse(message.data);
    console.log(msg);
}
"""
function message(route,args::Dict)
    id = pop!(args,"id");
    if haskey(pages,route)
        sessions = pages[route].sessions
        if haskey(sessions,id)
            client = sessions[id]
            if isopen(client)
                write(client,json(args))
            end
        end
    end
end

"""
Block Julia control flow until until callback["notify"](name) is called.
"""
function block(f::Function,name)
    conditions[name] = Condition()
    f()
    wait(conditions[name])
    delete!(conditions,name)
    return nothing
end

"""Add a JS library to the current page from a url."""
function add_library(url)
    name = basename(url)
    block(name) do
        Pages.broadcast("script","""
            var script = document.createElement("script");
            script.charset = "utf-8";
            script.type = "text/javascript";
            script.src = "$(url)";
            script.onload = function() {
                Pages.notify("$(name)");
            };
            document.head.appendChild(script);
        """)
    end
end

mutable struct Element
    id::String
    tag::String
    name::String
    attr::Dict{String,String}
    style::Dict{String,String}
    html::String
    text::String
    parent::String

    function Element(;id = "",tag = "div",name = "element",attr = Dict{String,String}(),style = Dict{String,String}(),html = "",text = "",parent = "body")
        new(id,tag,name,attr,style,html,text,parent)
    end
end

function assign(io::IO,element::Element)
    # ==========================================================================
    # Add attributes to element
    for prop in ["attr","style"]
        field = getfield(element,Symbol(prop))
        for key in keys(field)
            print(io,"""
                $(element.name).$(prop)("$(key)","$(field[key])");
            """)
        end
    end
    # ==========================================================================
    # Add html & text to element
    for prop in ["html","text"]
        field = getfield(element,Symbol(prop))
        if !isempty(field)
            print(io,"""
                $(element.name).$(prop)("$(field)");
            """)
        end
    end
    io
end

function add(io::IO,element::Element)
    # ==========================================================================
    # Get or create element
    print(io,"""
        var parent = d3.select("$(element.parent)")
        var $(element.name) = null;
    """)
    if !isempty(element.id)
        print(io,"""
            var check = document.getElementById("$(element.id)");
            if (check === null) {
                $(element.name) = parent.append("$(element.tag)").attr("id","$(element.id)");
            } else {
                $(element.name) = d3.select("#$(element.id)");
            };
        """)
    else
        print(io,"""
            $(element.name) = d3.select("$(element.parent)").append("$(element.tag)");
        """)
    end
    assign(io,element)
end
function add(element::Element)
    Pages.broadcast("script",String(take!(add(IOBuffer(),element))))
end

function append(io::IO,element::Element;parent = """d3.select("body")""")
    if isempty(element.name)
        print(io,"""
            $(parent).append("$(element.tag)");
        """)
    else
        print(io,"""
            $(element.name) = $(parent).append("$(element.tag)");
        """)
    end
    assign(io,element)
end

function remove(io::IO,tag;parent = """d3.select("body")""")
    print(io,"""
        $(parent).selectAll("$(tag)").remove();
    """)
    io
end
function remove(tag;parent = """d3.select("body")""")
    Pages.broadcast("script",String(take!(remove(IOBuffer(),tag,parent=parent))))
end

function add_select(io::IO,options,element::Element)
    element.tag == "select" || return warn("Element must have tag = select.")
    add(io,element)
    remove(io,"option",parent=element.name)
    # print(io,"""
    #     $(element.name).selectAll("option").remove();
    # """)
    for key in keys(options)
        print(io,"""
            $(element.name).append("option").attr("value","$(key)").text("$(options[key])");
        """)
    end
    io
end
function add_select(options,element::Element)
    Pages.broadcast("script",String(take!(add_select(IOBuffer(),options,element))))
end

# function add_table(io::IO,df::DataFrame;table = Element(tag="table",name="table"),tr = Element(tag="tr",name="row"),th = Element(tag="th",name="header"),td = Element(tag="td",name="cell"))
#     table.tag == "table" || return warn("Element must have tag = table.")
#     tr.tag == "tr" || return warn("Element must have tag = tr.")
#     th.tag == "th" || return warn("Element must have tag = th.")
#     td.tag == "td" || return warn("Element must have tag = td.")
#     add(io,table)
#     remove(io,"tr",parent=table.name)
#     # ==========================================================================
#     # Add header
#     print(io,"""
#         var $(tr.name) = null;
#     """)
#     append(io,tr,parent=table.name)
#     print(io,"""
#         var $(th.name) = null;
#     """)
#     for name in names(df)
#         th.html = string(name)
#         append(io,th,parent=tr.name)
#     end
#     # ==========================================================================
#     # Add data
#     print(io,"""
#         var $(td.name) = null;
#     """)
#     for irow in 1:size(df,1)
#         row = df[irow,:]
#         append(io,tr,parent=table.name)
#         for name in names(df)
#             td.html = string(row[name][1])
#             append(io,td,parent=tr.name)
#         end
#     end
#     io
# end
# function add_table(df::DataFrame;table = Element(tag="table",name="table"),tr = Element(tag="tr",name="row"),th = Element(tag="th",name="header"),td = Element(tag="td",name="cell"))
#     Pages.broadcast("script",String(take!(add_table(IOBuffer(),df,table=table,tr=tr,th=th,td=td))))
# end


text:
function getdocs(ctx::Context, cursor::CXCursor)
	getblock(ctx).flags.nodocs && return nothing

	comment = clang_Cursor_getParsedComment(cursor)
	kind = clang_Comment_getKind(comment)
	if kind == CXComment_FullComment
		return Markdown.MD(ctx, cursor, comment)
	elseif kind == CXComment_Null
		return Markdown.MD(ctx, cursor)
	else
		@warn "Expected a full comment, but got $(kind)"
		return nothing
	end
end

function Markdown.MD(ctx::Context, cursor::CXCursor)
	code = getcode(ctx, cursor)
	code = Markdown.Code(String(language(ctx)), code)

	loc = getlocation(cursor)
	if isnothing(loc)
		loc = "Defined by the parser/compiler"
	else
		loc = first(loc)
		loc = loc.file == header(ctx) ? getblock(ctx, loc).loc : loc
		loc = Markdown.Link("Defined at $(basename(String(loc.file))):$(loc.line)", "file://$(loc.file)")
	end

	return Markdown.MD(code, Markdown.Paragraph(loc))
end

function Markdown.MD(ctx::Context, cursor::CXCursor, comment::CXComment)
	hasDetails = false
	hasParams  = false
	hasReturns = false

	contents = []
	for ind in 1:clang_Comment_getNumChildren(comment)
		child = clang_Comment_getChild(comment, ind-1)
		kind = clang_Comment_getKind(child)

		if !hasDetails && !isempty(contents)
			hasDetails = true
			push!(contents, Markdown.Header("Details", 2))
		end

		if kind == CXComment_Paragraph
			Bool(clang_Comment_isWhitespace(child)) || push!(contents, Markdown.Paragraph(ctx, child))
		elseif kind == CXComment_BlockCommand
			para = clang_BlockCommandComment_getParagraph(child)
			para = Markdown.Paragraph(ctx, para)

			num = clang_BlockCommandComment_getNumArgs(child)
			cmd = _string(clang_BlockCommandComment_getCommandName, child)
			if cmd == "brief" || cmd == "par" || cmd == "paragraph"
				push!(contents, para)
			elseif cmd == "bug" || cmd == "note" || cmd == "warning" || cmd == "deprecated" || cmd == "attention"
				category =
					cmd == "bug" ? "danger" :
					cmd == "note" ? "info" :
					cmd == "warning" ? "warning" :
					cmd == "deprecated" ? "warning" :
					cmd == "attention" ? "danger" : "info"
				para = Markdown.Admonition("danger", titlecase(cmd), [para])
				push!(contents, para)
			elseif cmd == "sa" || cmd == "see"
				para = isempty(para.content) ? para : Markdown.Paragraph("See also: [`$(strip(first(para.content)))`](@ref)")
				push!(contents, para)
			elseif cmd == "returns" || cmd == "return"
				if !hasReturns
					hasReturns = true

					hdr = Markdown.Header("Returns", 2)
					if !isempty(contents) && contents[end] isa Markdown.Header
						contents[end] = hdr
					else
						push!(contents, hdr)
					end
				end

				push!(contents, para)
			elseif cmd == "li"
				if !isempty(contents) && contents[end] isa Markdown.List
					push!(contents[end].items, para)
				else
					push!(contents, Markdown.List(para))
				end
			else
				getblock(ctx).flags.notify && @warn "Unhandled comment block-command: $(cmd)"
			end
		elseif kind == CXComment_ParamCommand
			if !hasParams
				hasParams = true

				hdr = Markdown.Header("Parameters", 2)
				if !isempty(contents) && contents[end] isa Markdown.Header
					contents[end] = hdr
				else
					push!(contents, hdr)
				end
			end

			addParameter(ctx, contents, child)
		elseif kind == CXComment_VerbatimBlockCommand
			push!(contents, Markdown.Code(ctx, child))
		elseif kind == CXComment_VerbatimLine
			text = _string(clang_VerbatimLineComment_getText, child)
			push!(contents, Markdown.Paragraph(text))
		else
			getblock(ctx).flags.notify && @warn "Unhandled comment document child: $(kind)"
		end
	end

	if length(contents) > 1 && contents[end] isa Markdown.Header
		contents = contents[1:end-1]
	end

	return Markdown.MD([Markdown.MD(ctx, cursor).content..., contents...])
end

function Markdown.Paragraph(ctx::Context, cxcomment::CXComment)
	contents = []

	for ind in 1:clang_Comment_getNumChildren(cxcomment)
		child = clang_Comment_getChild(cxcomment, ind-1)
		kind = clang_Comment_getKind(child)

		if kind == CXComment_Text
			text = _string(clang_TextComment_getText, child)
			push!(contents, text)
		elseif kind == CXComment_InlineCommand
			addInline(ctx, contents, child)
		elseif kind == CXComment_HTMLStartTag || kind == CXComment_HTMLEndTag
			# TODO: handle HTML stuff...
		else
			getblock(ctx).flags.notify && @warn "Unhandled comment paragraph child: $(kind)"
		end
	end

	return Markdown.Paragraph(contents)
end

function Markdown.Code(ctx::Context, cxcomment::CXComment)
	lines = []

	for ind in 1:clang_Comment_getNumChildren(cxcomment)
		child = clang_Comment_getChild(cxcomment, ind-1)
		kind = clang_Comment_getKind(child)

		if kind == CXComment_VerbatimBlockLine
			line = _string(clang_VerbatimBlockLineComment_getText, child)

			# a line with newlines is probably incorrectly parsed, so only keep before a newline
			push!(lines, first(split(line, '\n', limit=2)))
		else
			getblock(ctx).flags.notify && @warn "Unhandled comment code-block child: $(kind)"
		end
	end

	return Markdown.Code(join(lines, '\n'))
end

function addInline(ctx::Context, contents, cxcomment)
	num = clang_InlineCommandComment_getNumArgs(cxcomment)
	cmd = _string(clang_InlineCommandComment_getCommandName, cxcomment)

	text = _string(clang_InlineCommandComment_getArgText, cxcomment, 0)
	kind = clang_InlineCommandComment_getRenderKind(cxcomment)
	if kind == CXCommentInlineCommandRenderKind_Bold
		push!(contents, Markdown.Bold(text))
	elseif kind == CXCommentInlineCommandRenderKind_Emphasized
		push!(contents, Markdown.Italic(text))
	elseif kind == CXCommentInlineCommandRenderKind_Monospaced
		push!(contents, Markdown.Code(text))  # TODO: maybe need Code(language, text)
	else
		push!(contents, text)
	end
end

function addParameter(ctx::Context, contents, cxcomment)
	num = clang_Comment_getNumChildren(cxcomment)
	num == 1 || error("Incorrect number of parameter-command comment children")

	para = clang_Comment_getChild(cxcomment, 0)
	clang_Comment_getKind(para) == CXComment_Paragraph || error("Expected a parameter-command comment paragraph")

	param = _string(clang_ParamCommandComment_getParamName, cxcomment)
	param = Markdown.Paragraph([Markdown.Code(param), ":", Markdown.Paragraph(ctx, para).content...,])

	if !isempty(contents) && contents[end] isa Markdown.List
		push!(contents[end].items, param)
	else
		push!(contents, Markdown.List(param))
	end
end


text:
"""
```
setuptables()
```

Initialize a database with the tables needed to run ICDataServer:
  - `users`
  - `servers`
  - `instrumentkinds`
  - `instruments`
  - `jobs`
  - `notes`

The creation of the tables can proceed without additional input, but in some
cases it may make sense for the tables to be automatically filled with data.

If new instrument kinds are added to `deps/instrumentkinds.json`, they will be
added to `instrumentkinds` table. If kinds are encountered in the database that
are not in the file, they remain in the database and are not deleted.

If new users are added to `deps/users.json`, they will be "upserted" into the
`users` table (users are created if necessary; if an existing user name is
attempted to be inserted into the table, then that user has their info updated
from the json file). Existing users not found in users.json remain in the table.
"""
function setuptables(dsn, path="")
    if path == ""
        # request path from user
        println("Path to directory with setup files? [defaults to ICDataServer/deps]:")
        path = chomp(readline(STDIN))

        # default
        if path == ""
            path = joinpath(Pkg.dir("ICDataServer"), "deps")
        end
    end
    !ispath(path) && error("could not find path.")
    println(path)

    userstable(dsn, path; filename = "users.json")
    servertable(dsn, path; filename = "servers.json")
    restrictiontable(dsn, path; filename = "instrumentkinds.json",
        tablename = "instrumentkinds", colname = "kind")
    # restrictiontable(dsn, path; filename = "instrumentmakes.json",
    #     tablename = "instrumentmakes", colname = "make")
    # restrictiontable(dsn, path; filename = "instrumentmodels.json",
    #     tablename = "instrumentmodels", colname = "model")
    instrtable(dsn)
    jobstable(dsn)
    notestable(dsn)
end

function jobstable(dsn)
    ODBC.execute!(dsn,
    """
    CREATE TABLE IF NOT EXISTS jobs(
        job_id      bigserial PRIMARY KEY,
        cryostat    character varying,
        username    character varying REFERENCES users(username),
        device      character varying,
        nmeas       integer,
        jobsubmit   timestamp NOT NULL,
        jobstart    timestamp,
        jobstop     timestamp,
        jobstatus   smallint,
        dataserver  character varying NOT NULL REFERENCES servers(alias)
    );
    """)
end

function instrtable(dsn)
    ODBC.execute!(dsn,
    """
    CREATE TABLE IF NOT EXISTS instruments(
        alias	    character varying PRIMARY KEY,
        make        character varying,
        model	    character varying,
        serialnum	character varying,
        kind	    character varying REFERENCES instrumentkinds(kind),
        hw_server	character varying REFERENCES servers(alias),
        protocol	character varying,
        address     character varying
    );
    """)
end

function restrictiontable(dsn, path; filename="", tablename="", colname="")
    ODBC.execute!(dsn,
    """
    CREATE TABLE IF NOT EXISTS $tablename(
        $colname    character varying PRIMARY KEY
    );
    """)

    ipath = joinpath(path, filename)
    if isfile(ipath)
        stmt = ODBC.prepare(dsn,
            "INSERT INTO $tablename VALUES (?) ON CONFLICT DO NOTHING;")
        for t in JSON.parsefile(ipath)[colname]
            ODBC.execute!(stmt, (t,))
        end
    end
end

function userstable(dsn, path; filename="")

    ODBC.execute!(dsn,
    """
    CREATE TABLE IF NOT EXISTS users(
        username  character varying PRIMARY KEY,
        first     character varying,
        middle    character varying,
        last      character varying,
        email     character varying,
        phone     character varying,
        office    character varying
    );
    """)

    # add "default" user in case it is absent from the config file.
    ODBC.execute!(dsn, """
    INSERT INTO users (username, first, last) VALUES ('default', 'No', 'Name')
    ON CONFLICT (username) DO NOTHING;
    """)

    ipath = joinpath(path, filename)
    if isfile(ipath)
        for d in JSON.parsefile(ipath)["users"]
            k = keys(d)
            kstr = reduce((a,b)->a*","*b, "$ki" for ki in k)
            vstr = reduce((a,b)->a*","*b, "'$(d[ki])'" for ki in k)
            pstr = reduce((a,b)->a*","*b, "$ki = '$(d[ki])'" for ki in k)
            query = """
            INSERT INTO users ($kstr) VALUES ($vstr) ON CONFLICT (username)
            DO UPDATE SET $pstr ;
            """
            ODBC.execute!(dsn, query)
        end
    end
end

function servertable(dsn, path; filename="")
    ODBC.execute!(dsn,
    """
    CREATE TABLE IF NOT EXISTS servers(
        alias       character varying PRIMARY KEY,
        address     character varying NOT NULL,
        port        integer
    );
    """)

    ipath = joinpath(path, filename)
    if isfile(ipath)
        for d in JSON.parsefile(ipath)["servers"]
            k = keys(d)
            kstr = reduce((a,b)->a*","*b, "$ki" for ki in k)
            vstr = reduce((a,b)->a*","*b, "'$(d[ki])'" for ki in k)
            pstr = reduce((a,b)->a*","*b, "$ki = '$(d[ki])'" for ki in k)
            query = """
            INSERT INTO servers ($kstr) VALUES ($vstr) ON CONFLICT (alias)
            DO UPDATE SET $pstr ;
            """
            ODBC.execute!(dsn, query)
        end
    end
end

function notestable(dsn)
    ODBC.execute!(dsn,
    """
    CREATE TABLE IF NOT EXISTS notes(
        note_id  serial PRIMARY KEY,
        job_id   bigint REFERENCES jobs(job_id) NOT NULL,
        username character varying REFERENCES users(username) NOT NULL,
        dt       timestamp NOT NULL,
        note     character varying NOT NULL
    );
    """)
end


text:
#HW2 Problem 2
#Minhao Yan my497, Wentong Chen wc422
#1. Generate array of random numbers named x
x = rand(Float64, (20000,20000))

#2. Column by Column
exp_column = zeros((20000,20000))
function exp_cols(x)
    for i in 1:20000
        exp_column[:,i] = exp.(x[:,i])
    end
    return exp_column
end

#3. Row by Row
exp_row = zeros((20000,20000))
function exp_rows(x)
    for i in 1:20000
        exp_row[i,:] = exp.(x[i,:])
    end
    return exp_row
end

#4.Calculate time
@time exp_cols(x)
#23.460586 seconds (523.19 k allocations: 5.984 GiB, 4.02% gc time)
@time exp_rows(x)
#320.930191 seconds (273.41 k allocations: 5.970 GiB, 0.56% gc time)

#5. The columnn by column one is much faster.

#6. Exponentiate the full array
@time exp_full = exp.(x)
#55.045882 seconds (314.79 k allocations: 2.996 GiB, 0.01% gc time)
#Overall, the column by column method is the fastest, then the overall one,
#the slowest is the row by row method.


text:
include("./imports.jl")

@testset ExtendedTestSet "generate_spatial_bounding_box" begin
    @testset ExtendedTestSet "generate_spatial_bounding_box" begin
    x = [0 1 1 0; 0 1 1 0]
    roi_start, roi_end = generate_spatial_bounding_box(x, x -> x > 0)
    @test roi_start == [1, 2] && roi_end == [2, 3]
    end
end


text:
## ---
	using Flux: onehot
	using Transformers: Vocabulary, decode
	using Random
	using Flux: throttle, @epochs, ADAMW, ADAM, RADAM, Momentum, Optimiser, logitcrossentropy, params
	using Flux: softmax, onecold, WeightDecay
	using DataStructures: Queue
	Random.seed!(42)

## ---
module Datasets include("src/datasets.jl") end
dataSetlength, vocab, sample = Datasets.openDataset("input.txt")
batch(eachSampleSize,batchSize ) = hcat((sample(eachSampleSize) for x in 1:batchSize)...)

## ---
module Model include("src/model.jl") end

lengthOfSentences=280
numberOfSentencesInBatch=15
model=Model.minGPT(maxBlockSize=lengthOfSentences, vocabSize=length(vocab))
@show NumberOfParameters=sum(length(p) for p in params(model))

## ---

# for a batch of sentences, extract 1st to one before last character as input
# and 2nd to last char of each sentence as prediction.  Shape the predictred
# values as a long one dim series of values, one hot encode them using the vocab

trainingDataPair(b)=(b[1:end-1, :], onehot(vocab, reshape(b[2:end,:], :)))
trainingData(eachSampleSize=lengthOfSentences,batchSize=numberOfSentencesInBatch)= [trainingDataPair(batch(eachSampleSize,batchSize))
		for i in 1:dataSetlength/batchSize]  # count each sentence, no matter its length as one sample

function loss(x,y)
	yhat=model(x)
	logitcrossentropy(reshape(yhat,size(yhat)[1], :), y)
end

function validate()
	s=collect("Wher")
	for i in 1:20
		out=model(vocab(s))
		push!(s, decode(vocab,onecold(softmax(out)[:,end])))
	end
	join(s)
end
opt = ADAM(0.00058, (0.9, 0.95))

## ---
module Train include("src/train.jl") end

## ---

decay = 0.1
optWithWeightDecay=Optimiser(opt,
			Train.SelectiveWeightDecay(WeightDecay(decay), Model.noDecayList(model)))

Train.train!(loss,
 				params(model),
				trainingData(),  # use the defaults
				optWithWeightDecay,
				validate, repeat=2)

## ---

# pass a string to the model and see the top prediction for the next char in sequence
out=model(vocab(collect("Where ")))
@show sort(collect(zip(softmax(out[:,end]), decode(vocab, collect(1:size(out)[1])))), rev=true)
@show decode(vocab,onecold(softmax(out)[:,end]))

## ---
print(validate())

## ---

#trying to track allocated
d = [trainingData(10, 2)[1]]   # use small sentence size, 10, and 2 sentences in batch to get idea of timing
@timed Train.train!(loss,
 				params(model),
				d,
				opt)
## ---
using Profile
numberOfSentencesInBatch=10
lengthOfSentences=
d = [trainingData(512, 10)[1]]  # 10 sentnces of 512 char

Profile.clear()

@profile Train.train!(loss,
 				params(model),
				d,
				opt)

# (Since Julia's garbage collector is written in C, such events can be detected using the C=true
# output mode described below, or by using ProfileView.jl.)
#https://docs.julialang.org/en/v1/manual/profile/#Memory-allocation-analysis-1

# format – Introduced above, determines whether backtraces are printed with (default, :tree) or without (:flat) indentation indicating tree structure.
# C – If true, backtraces from C and Fortran code are shown (normally they are excluded). Try running the introductory example with Profile.print(C = true). This can be extremely helpful in deciding whether it's Julia code or C code that is causing a bottleneck; setting C = true also improves the interpretability of the nesting, at the cost of longer profile dumps.
# combine – Some lines of code contain multiple operations; for example, s += A[i] contains both an array reference (A[i]) and a sum operation. These correspond to different lines in the generated machine code, and hence there may be two or more different addresses captured during backtraces on this line. combine = true lumps them together, and is probably what you typically want, but you can generate an output separately for each unique instruction pointer with combine = false.
# maxdepth – Limits frames at a depth higher than maxdepth in the :tree format.
# sortedby – Controls the order in :flat format. :filefuncline (default) sorts by the source line, whereas :count sorts in order of number of collected samples.
# noisefloor – Limits frames that are below the heuristic noise floor of the sample (only applies to format :tree). A suggested value to try for this is 2.0 (the default is 0). This parameter hides samples for which n <= noisefloor * √N, where n is the number of samples on this line, and N is the number of samples for the callee.
# mincount – Limits frames with less than mincount occurrences.
#

Juno.profiler(;C = false)

####### TODO:  Try https://github.com/astrieanna/TypeCheck.jl   to make sure types are stable in the network

## ---

function sampleFromModel(s="Where ", minprob=.05)
	function explore(x, minprob)
		x[1] < minprob && return currentRes
		out=model(vocab(collect(x[2])))
		choices=sort(collect(zip(softmax(out[:,end]), decode(vocab, collect(1:size(out)[1])))), rev=true)
		options = [(x[1]*np, x[2]*nc) for (np,nc) in choices if (x[1]*np) > minprob ]
	end
	toexplore=Queue{Tuple{Float32,String}}()
	enqueue!(toexplore,(1, s))  #start with the initial string and prob 1
	results=[]
	for i in toexplore
		nxt=dequeue!(toexplore)
		res=explore(nxt, minprob)
		if isempty(res)  # go as long as you can't add any more character and satisfy the prob
			push!(results, nxt)
			continue #continue emptying the queue
		end
		for r in res
			enqueue!(toexplore, r)  ## add elements for further exploration
		end
	end
    #elements of the results are (prob, string), sort from high to low of the
	#longest match with highest prob
	sort(results, by=x->(length(x[2]), x[1]), rev=true)
end

@show res = sampleFromModel("VINCENTIO: Oh God, My lord!", 0.01)

## ---
# Try a random model with the data and see what the losses would look like
using Statistics: norm, mean, std
function justRandomData()
	function readData(f)
		s = open(f) do file
			read(file, String)
		end
	end
	function generateRandomY()

		randomExpectedString = rand(s, lengthOfSentences-1, numberOfSentencesInBatch)
		randomExpectedy= reshape(onehot(vocab, vocab.(randomExpectedString)), length(vocab), :)
	end

	# let model output  be a complete random numbers, => uniform distribution
	#randModel(x,y)=rand(Float32, length(vocab), lengthOfSentences-1, numberOfSentencesInBatch), y

	# shuffle the model output, so it is not uniform, but  has all the distributions of the model, but
	# result is all in wrong place (this would be a case where model outout is not properly wired!)
	#randModel(x,y)=shuffle(model(x)), y

	# let the output be a random sample of the text.  So we respect the distribution of the input
	# so the output is wrong but it is wrong with the same distribution as the data
	# result is uniform+entropy of the system
	s=readData("input.txt")
	randModel(x,y)=model(x), generateRandomY()
	losses=[]
	ctr=0
	for (x,y) in trainingData(lengthOfSentences, numberOfSentencesInBatch)
		yhat, y =randModel(x,y)
		@show l=logitcrossentropy(reshape(yhat,size(yhat)[1], :), y)
		push!(losses, l)
		ctr +=1
		if  ctr > 10 break end
	end
	@show mean(losses), std(losses)
end

justRandomData()

## ---
using StatsBase
using Plots
using LinearAlgebra
function getOutputResults(m, data, limit=5)
	output=[]
	ctr=0
	for (x,y) in data
		out=softmax(m(x))
		out=reshape(out,size(out)[1], :)
		push!(output, onecold(out)...)  ## get the encoded vocab using straight onecold, not onecold(vocab, ..) version

		ctr +=1
		if  ctr > limit break end
	end
	convert(Array{Int}, output)
end
function OutputHistogram(m)
	outChars=getOutputResults(m, d)
	fit(Histogram, outChars, 0:1:length(vocab))
end

#Try a untrained model
raw=OutputHistogram(Model.minGPT(maxBlockSize=lengthOfSentences, vocabSize=length(vocab)))

#... and trained model
trained=OutputHistogram(model)

inputData=fit(Histogram, vocab(collect(s)), 0:1:length(vocab))
plot(normalize(raw, mode=:pdf))
plot!(normalize(trained, mode=:pdf))
plot!(normalize(inputData, mode=:pdf))

max_stars_repo_name: daryoush/minGPT

max_issues_repo_name: daryoush/minGPT

max_forks_repo_name: daryoush/minGPT


text:
export Controller

struct Controller
    cell_division::Function # cell params -> bool
    new_cell_params::Function # cell params -> child params
    cell_death::Function # cell params -> bool
    cell_state_update::Function # cell params, state, input -> state, output
    cell_param_update::Function # cell params, state -> params

    connect::Function # cell 1 params, cc 2 params -> bool
    new_conn_params::Function # cell 1 params, cc 2 params -> conn params
    disconnect::Function # c1 params, c2 params, conn params, state -> bool
    conn_state_update::Function # conn params, state, cell input, conn input -> state, output
    conn_param_update::Function # c1 params, c2 params, conn params, state -> params
end

function Controller(cfg::Dict)
    nouts = output_lengths(cfg)
    cell_division(x::Array{Float64}) = false
    new_cell_params(x::Array{Float64}) = zeros(nouts[2])
    cell_death(x::Array{Float64}) = false
    cell_state_update(x::Array{Float64}) = zeros(nouts[4])
    cell_param_update(x::Array{Float64}) = zeros(nouts[5])
    connect(x::Array{Float64}) = false
    new_conn_params(x::Array{Float64}) = zeros(nouts[7])
    disconnect(x::Array{Float64}) = false
    conn_state_update(x::Array{Float64}) = zeros(nouts[9])
    conn_param_update(x::Array{Float64}) = zeros(nouts[10])

    Controller(cell_division, new_cell_params, cell_death,
               cell_state_update, cell_param_update,
               connect, new_conn_params, disconnect,
               conn_state_update, conn_param_update)
end


text:
for n in [2, 20, 200, 2000]
    x = randn(Float32, n, n)
    model = Dense(n, n)
    println("CPU n=$n")
    run_benchmark(model, x, cuda=false)
    println("CUDA n=$n")
    run_benchmark(model, x, cuda=true)
end


text:
"""
Please add the latest Kinetic.jl package as it fixes the issue with the regime classifier.
Backend requirement:
- KitBase: v0.7.6
- KitML: v0.4.4
"""

using Kinetic, NPZ, Plots
using KitBase.ProgressMeter: @showprogress
cd(@__DIR__)

set = Setup(case = "sod", space = "1d1f1v", maxTime = 0.15)
ps = PSpace1D(0.0, 1.0, 100, 1)
vs = VSpace1D(-5.0, 5.0, 64)
gas = Gas(Kn = 1e-3, K = 0, γ = 3)
ib = IB1F(ib_sod(set, ps, vs, gas)...)
ks = SolverSet(set, ps, vs, gas, ib)
ctr, face = init_fvm(ks)

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(ks.set.maxTime ÷ dt) + 1
res = zero(ctr[1].w)

X, Y, Z, I = begin
    prim = [1.0, 0.0, 1.0]
    w = prim_conserve(prim, ks.gas.γ)
    sw = zeros(3)
    M = maxwellian(ks.vs.u, prim)
    [[1e-4, 0.0, 1e-4]; zeros(3); 1.0], 0.0, M, 0
end

@showprogress for iter = 1:nt
    evolve!(ks, ctr, face, dt)
    update!(ks, ctr, face, dt, res)

    if iter % 3 == 0
        for i = 1:ks.ps.nx
            sw = (ctr[i+1].w .- ctr[i-1].w) / ks.ps.dx[i] / 2.0
            x, y = regime_data(ks, ctr[i].w, sw, ctr[i].f)
            global X = hcat(X, x)
            global Y = hcat(Y, y)
            global Z = hcat(Z, ctr[i].f)
            global I = hcat(I, iter)
        end
    end
end

#plot(ks, ctr)

npzwrite("data.npz", Dict("I" => I, "X" => X, "Y" => Y, "Z" => Z))

max_stars_repo_path: src/1d/sod/let's analyze sampler/sod.jl

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_path: src/1d/sod/let's analyze sampler/sod.jl

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_path: src/1d/sod/let's analyze sampler/sod.jl

max_forks_repo_name: CSMMLab/Flowmachine


text:
using Test

include("accumulate.jl")

@testset "multiply" begin
    @test accumulate(collect(1:10), *) == [1, 4, 9, 16, 25, 36, 49, 64, 81, 100]
end

@testset "divide" begin
    @test accumulate(collect(1:10), /) == [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
end

@testset "subtract" begin
    @test accumulate(collect(1:10), -) == [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
end

@testset "add" begin
    @test accumulate(collect(1:10), +) == [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
end

@testset "mod" begin
    @test accumulate(collect(1:10), %) == [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
end

@testset "square" begin
    @test accumulate(collect(1:10), ^) == [1, 4, 27, 256, 3125, 46656, 823543, 16777216, 387420489, 10000000000]
end


text:
# Taken from https://github.com/JuliaData/DataFrames.jl/blob/main/test/runtests.jl
using SatelliteImages
using Test

fatalerrors = length(ARGS) > 0 && ARGS[1] == "-f"
quiet = length(ARGS) > 0 && ARGS[1] == "-q"
anyerrors = false

my_tests = ["aggregate.jl", "coordinate_system.jl", "data_io.jl", "f_apply.jl", "mask_area.jl", "polygons.jl", "sparse_cube.jl"]

println("Running tests:")

for my_test in my_tests
    try
        include(my_test)
        println("\t\033[1m\033[32mPASSED\033[0m: $(my_test)")
    catch e
        global anyerrors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(my_test)")
        if fatalerrors
            rethrow(e)
        elseif !quiet
            showerror(stdout, e, backtrace())
            println()
        end
    end
end

if anyerrors
    throw("Tests failed")
end


text:
using Metal

@show devices()
dev = MtlDevice(1)

#lib = MTL.LibraryWithFile(d, "default.metallib")
src = read(dirname(pathof(Metal))*"/Metal/kernels/add.metal", String)

bufferSize = 128
bufferA = MtlArray{Float32,1}(undef, tuple(bufferSize), storage=Shared)
bufferB = MtlArray{Float32,1}(undef, tuple(bufferSize), storage=Shared)
bufferC = MtlArray{Float32,1}(undef, tuple(bufferSize), storage=Shared)

vecA = unsafe_wrap(Vector{Float32}, bufferA.buffer, tuple(bufferSize))
vecB = unsafe_wrap(Vector{Float32}, bufferB.buffer, tuple(bufferSize))
vecC = unsafe_wrap(Vector{Float32}, bufferC.buffer, tuple(bufferSize))

using Random
rand!.([vecA, vecB])

## Setup
opts = MtlCompileOptions()
lib = MtlLibrary(dev, src, opts)

fun = MtlFunction(lib, "add_arrays")
pip_addfun = MtlComputePipelineState(dev, fun)
queue = global_queue(dev) #MtlCommandQueue(dev)

##
vecA .= 0.0; vecB .= 0.0; vecC .= 0.0;
cmd = MTL.commit!(queue) do cmdbuf
    MtlComputeCommandEncoder(cmdbuf) do enc
        MTL.set_function!(enc, pip_addfun)
        MTL.set_buffer!(enc, bufferA.buffer, 0, 1)
        MTL.set_buffer!(enc, bufferB.buffer, 0, 2)
        MTL.set_buffer!(enc, bufferC.buffer, 0, 3)
        #MTL.set_buffers!(enc,
        #                        [bufferA.buffer, bufferB.buffer, bufferC.buffer],
        #                        [0,0,0], 1:3)
        gridSize = MtSize(length(vecA), 1, 1)
        threadGroupSize = min(length(vecA), pip_addfun.maxTotalThreadsPerThreadgroup)
        threadGroupSize = MTL.MtSize(threadGroupSize, 1, 1)
        @info threadGroupSize
        MTL.append_current_function!(enc, gridSize, threadGroupSize)
    end
end

# Execute
wait(cmd)

@show vecC


text:
module lab7

using Distributions
using QuadGK

include("model_spectrum.jl")

end # module