text:
using DiffEqFlux, OrdinaryDiffEq, Test
using DiffEqFlux: group_ranges

## Test group partitioning helper function
@test group_ranges(10, 4) == [1:4, 4:7, 7:10]
@test group_ranges(10, 5) == [1:5, 5:9, 9:10]
@test group_ranges(10, 10) == [1:10]
@test_throws DomainError group_ranges(10, 1)
@test_throws DomainError group_ranges(10, 11)

## Define initial conditions and time steps
datasize = 30
u0 = Float32[2.0, 0.0]
tspan = (0.0f0, 5.0f0)
tsteps = range(tspan[1], tspan[2], length = datasize)

# Get the data
function trueODEfunc(du, u, p, t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
prob_trueode = ODEProblem(trueODEfunc, u0, tspan)
ode_data = Array(solve(prob_trueode, Tsit5(), saveat = tsteps))

# Define the Neural Network
nn = FastChain((x, p) -> x.^3,
                FastDense(2, 16, tanh),
                FastDense(16, 2))
p_init = initial_params(nn)

neuralode = NeuralODE(nn, tspan, Tsit5(), saveat = tsteps)
prob_node = ODEProblem((u,p,t)->nn(u,p), u0, tspan, p_init)

function predict_single_shooting(p)
    return Array(neuralode(u0, p))
end

# Define loss function
function loss_function(data, pred)
	return sum(abs2, data - pred)
end

## Evaluate Single Shooting
function loss_single_shooting(p)
    pred = predict_single_shooting(p)
    l = loss_function(ode_data, pred)
    return l, pred
end

res_single_shooting = DiffEqFlux.sciml_train(loss_single_shooting, neuralode.p,
                                          ADAM(0.05),
										  maxiters = 300)

loss_ss, _ = loss_single_shooting(res_single_shooting.minimizer)
println("Single shooting loss: $(loss_ss)")

## Test Multiple Shooting
group_size = 3
continuity_term = 200

function loss_multiple_shooting(p)
    return multiple_shoot(p, ode_data, tsteps, prob_node, loss_function, Tsit5(),
                          group_size; continuity_term,
                          abstol=1e-8, reltol=1e-6) # test solver kwargs
end

res_ms = DiffEqFlux.sciml_train(loss_multiple_shooting, neuralode.p,
                                ADAM(0.05), maxiters = 300)

# Calculate single shooting loss with parameter from multiple_shoot training
loss_ms, _ = loss_single_shooting(res_ms.minimizer)
println("Multiple shooting loss: $(loss_ms)")
@test loss_ms < 10loss_ss

# Test with custom loss function
group_size = 4
continuity_term = 50

function continuity_loss_abs2(û_end, u_0)
    return sum(abs2, û_end - u_0) # using abs2 instead of default abs
end

function loss_multiple_shooting_abs2(p)
    return multiple_shoot(p, ode_data, tsteps, prob_node,
                          loss_function, continuity_loss_abs2, Tsit5(),
                          group_size; continuity_term)
end

res_ms_abs2 = DiffEqFlux.sciml_train(loss_multiple_shooting_abs2, neuralode.p,
                                     ADAM(0.05), maxiters = 300)

loss_ms_abs2, _ = loss_single_shooting(res_ms_abs2.minimizer)
println("Multiple shooting loss with abs2: $(loss_ms_abs2)")
@test loss_ms_abs2 < loss_ss

## Test different SensitivityAlgorithm (default is InterpolatingAdjoint)
function loss_multiple_shooting_fd(p)
    return multiple_shoot(p, ode_data, tsteps, prob_node,
                          loss_function, continuity_loss_abs2, Tsit5(),
                          group_size; continuity_term,
                          sensealg=ForwardDiffSensitivity())
end

res_ms_fd = DiffEqFlux.sciml_train(loss_multiple_shooting_fd, neuralode.p,
                                ADAM(0.05), maxiters = 300)

# Calculate single shooting loss with parameter from multiple_shoot training
loss_ms_fd, _ = loss_single_shooting(res_ms_fd.minimizer)
println("Multiple shooting loss with ForwardDiffSensitivity: $(loss_ms_fd)")
@test loss_ms_fd < 10loss_ss

# Integration return codes `!= :Success` should return infinite loss.
# In this case, we trigger `retcode = :MaxIters` by setting the solver option `maxiters=1`.
loss_fail, _ = multiple_shoot(p_init, ode_data, tsteps, prob_node, loss_function, Tsit5(),
                              datasize; maxiters=1, verbose=false)
@test loss_fail == Inf

## Test for DomainErrors
@test_throws DomainError multiple_shoot(p_init, ode_data, tsteps, prob_node,
                                        loss_function, Tsit5(), 1)
@test_throws DomainError multiple_shoot(p_init, ode_data, tsteps, prob_node,
                                        loss_function, Tsit5(), datasize + 1)

## Ensembles
u0s = [Float32[2.0, 0.0], Float32[3.0, 1.0]]
function prob_func(prob, i, repeat)
    remake(prob, u0 = u0s[i])
end
ensemble_prob = EnsembleProblem(prob_node, prob_func = prob_func)
ensemble_prob_trueODE = EnsembleProblem(prob_trueode, prob_func = prob_func)
ensemble_alg = EnsembleThreads()
trajectories = 2
ode_data_ensemble = Array(solve(ensemble_prob_trueODE, Tsit5(), ensemble_alg, trajectories = trajectories, saveat = tsteps))

group_size = 3
continuity_term = 200
function loss_multiple_shooting_ens(p)
    return multiple_shoot(p, ode_data_ensemble, tsteps, ensemble_prob, ensemble_alg,
                          loss_function, Tsit5(),
                          group_size; continuity_term,
                          trajectories,
                          abstol=1e-8, reltol=1e-6) # test solver kwargs
end

res_ms_ensembles = DiffEqFlux.sciml_train(loss_multiple_shooting_ens, neuralode.p,
                                ADAM(0.05), maxiters = 300)

loss_ms_ensembles, _ = loss_single_shooting(res_ms_ensembles.minimizer)

println("Multiple shooting loss with EnsembleProblem: $(loss_ms_ensembles)")

@test loss_ms_ensembles < 10loss_ss


text:
#

calc(::BM_TEST, ::FUNT, x, v, n, nel, ix, incx) = calc(BM_JBase1(), FUNT(), x, v, n, nel, ix, incx)
calc(::BM_JBase1, ::FUNT, x, v, n, nel, ix, incx) = scale!(view(x, ix:incx:n),v)
calc(::BM_JBase2, ::FUNT, x, v, n, nel, ix, incx) = scale!(v, view(x, ix:incx:n))
calc(::BM_BLAS, ::FUNT, x, v, n, nel, ix, incx) = BLAS.scal!(nel, v, pointer(x,ix), incx)
calc(::BM_Broadcast, ::FUNT, x, v, n, nel, ix, incx) = begin
    vx = view(x, ix:incx:n)
    broadcast!(*,vx,vx,v)
end
calc(::BM_Forloop, ::FUNT, x, v, n, nel, ix, incx) = scaleloopgen(x, v, n, ix, incx)
function scaleloopgen(x::Array, v, n, ix, incx)
    @inbounds for i=ix:incx:n
        x[i] *= v
    end
    return x
end
calc(::BM_FAO, ::FUNT, x, v, n, nel, ix, incx) = unsafe_fast_scale!(x, ix, incx, v, nel)

Base.start{T}(p::BenchCase{T,FUNT}, n::Int) =
    (rand(n), float64(1.000001), n, nmax2nel(2, 40, n), 2, 40)

function Base.run{Op,FUNT}(p::BenchCase{Op,FUNT}, n::Int, s)  # bench type, config, start value
    calc(Op(), FUNT(), s...)
    return nothing
end

ext: jl


text:
module TMLEEpistasis

using DataFrames
using CSV
using SnpArrays
using MLJ
using Distributions
using TMLE
using TOML
using BGEN
using HighlyAdaptiveLasso

import MLJ: fit, transform, target_scitype

###############################################################################
# GENERIC FUNCTIONS

is_binary(y) = sort(unique(y)) == [0, 1]

###############################################################################
# INCLUDES

include("ukbb.jl")
include("models.jl")
include("estimators.jl")
include("utils.jl")

###############################################################################
# EXPORTS

export UKBBVariantRun

end


text:
@testset "Read input" begin
    testinput = JuSwalbe.readinput("test.txt")
    @test isa(testinput, JuSwalbe.Inputconstants)
    @test testinput.lx == 10
    @test testinput.ly == 5
    @test testinput.maxruntime == 1000
    @test testinput.dumping == 100
    @test testinput.τ == 1.0
    @test testinput.gravity == 0.0
    @test testinput.γ == 0.01
    @test testinput.δ == 1.0
end

@testset "Find argument" begin
    arr = ["hello" 15; "my" 0.077; "wonderful" 10534; "world" 0.4231]
    @test findargument(arr, "wonderful") == 10534;
    @test findargument(arr, "hello") == 15;
    @test findargument(arr, "world") == 0.4231;
    @test findargument(arr, "my") == 0.077;
    # Returns nothing if the string can't be found
    @test findargument(arr, "lol") == nothing
    # Has to be a string
    @test_throws MethodError findargument(arr, 2)
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
import Base: pairs

mutable struct BinaryNode{K, V} <: AbstractNode{K, V}
    k::K
    v::V
    l::BinaryNode{K, V}
    r::BinaryNode{K, V}
    p::BinaryNode{K, V}
    function BinaryNode{K, V}() where {K, V}
        self = new{K, V}()
        self.l = self.r = self.p = n= self
        return self
    end
end

function divide_and_conquer(pred::Function, nil::BinaryNode{K, V}) where {K, V}
    predl, predr, k, v = pred()
    predl === nothing && predr === nothing && return (nil, 0)
    node = BinaryNode{K, V}()
    (node.k, node.v) = (k, v)
    node.l, cl = predl === nothing ? (nil, 0) : divide_and_conquer(predl, nil)
    node.r, cr = predr === nothing ? (nil, 0) : divide_and_conquer(predr, nil)
    node.l.p = node.r.p = node
    return node, cl + cr + 1
end

mutable struct BinaryTree{K, V} <: AbstractBinaryTree{K, V}
    nil::BinaryNode{K, V}
    root::BinaryNode{K, V}
    n::Int
    function BinaryTree{K, V}(pred::Function) where {K, V}
        nil = BinaryNode{K, V}()
        root, n = divide_and_conquer(pred, nil)
        return new(nil, root, n)
    end
end

Rectangle.isnil(t::BinaryTree, node::BinaryNode) = node === t.nil

function Base.pairs(t::BinaryTree{K, V}) where {K, V}
    ps = Vector{Tuple{K, V}}(undef, length(t))
    i = 0
    _inorder(t.root, t) do n
        ps[i+=1] = (n.k, n.v)
        return true
    end
    return ps
end

function Base.keys(t::BinaryTree{K, V}) where {K, V}
    ks = Vector{K}(undef, length(t))
    i = 0
    _inorder(t.root, t) do n
        ps[i+=1] = n.k
        return true
    end
    return ps
end

function Base.values(t::BinaryTree{K, V}) where {K, V}
    ks = Vector{V}(undef, length(t))
    i = 0
    _inorder(t.root, t) do n
        ps[i+=1] = n.v
        return true
    end
    return ps
end

struct XYData{T}
    dir::Int
    loc::T
    range::Tuple{T, T}
end

@inline function Line(xy::XYData{T}) where T
    dir, r1, r2, l = xy.dir, xy.range..., xy.loc
    return dir == 1 ? Line(r1, l, r2, l) : Line(l, r1, l, r2)
end

const XYNode{T, V} = BinaryNode{XYData{T}, V}
const XYTree{T, V} = BinaryTree{XYData{T}, V}

function get_values(t::XYTree{K, V}) where {K, V}
    lines, vs = Vector{Line{K}}(undef, length(t)), Vector{V}(undef, length(t))
    ps = pairs(t)
    for i = firstindex(ps):lastindex(ps)
        lines[i], vs[i] = Line(ps[i][1]), ps[i][2]
    end
    return lines, vs
end


text:
# Add functions for calculating the normalized compression distance (NCD).

# Compressors typically write a header to the compressed strings. In NCD calculations
# we do not want to include the length of this header. For each compressor type
# we try to find the length of the header by compressing some very short strings.
# We cache these min lengths in a hash for speedy access later.
const CompressorHeaderLengths = Dict{Any,Int}()

function header_length{C <: Compressor}(compressor::C)
	ct = typeof(compressor)
	if !haskey(CompressorHeaderLengths, ct)
		CompressorHeaderLengths[ct] = find_min_header_length(compressor)
	end
	return CompressorHeaderLengths[ct]
end

const ShortStrings = ["", "a", "b", "z", "A", "B", "Z", "0", "aa", "bb", "zz", "aaa", "zzz"]

function find_min_header_length{C <: Compressor}(compressor::C)
	lens = Int[]
	map(ShortStrings) do s
	  try
	  	push!(lens, clen(compressor, s))
	  catch ex
	  end
	end
	minimum(lens)
end

function ncdcalc(cx::Int, cy::Int, cxy::Int)
	(cxy - min(cx, cy)) / max(cx, cy)
end

function ncd{S <: AbstractString}(c::Compressor, x::S, y::S)
	hlen = header_length(c)
	cx = clen(c, x) - hlen
	cy = clen(c, y) - hlen
	cxy = clen(c, x, y) - hlen
	ncdcalc(cx, cy, cxy)
end

const DefaultCompressor = LibzCompressor()

ncd{S <: AbstractString}(x::S, y::S, c::Compressor = DefaultCompressor) = ncd(c, x, y)
ncd{S1 <: AbstractString, S2 <: Any}(x::S1, y::S2, c::Compressor = DefaultCompressor) = ncd(c, x, string(y))
ncd{S1 <: Any, S2 <: AbstractString}(x::S1, y::S2, c::Compressor = DefaultCompressor) = ncd(c, string(x), y)
ncd{S1 <: Any, S2 <: Any}(x::S1, y::S2, c::Compressor = DefaultCompressor) = ncd(c, string(x), string(y))

# NCD for lists is here called NCDL. It is an implementation of the "NCD for Multisets"
# as described in the paper:
#  A. Cohen & P. Vitanyi, "", 2012, http://arxiv.org/pdf/1212.5711.pdf
# This is a heuristic to calc the NCDL in O(n^2) operations since the def of NCDL implies
# a O(2^n) algorithm.
function ncdl{T <: Any}(Xs::Vector{T}, c::Compressor = LibzCompressor())
  ncdl(map(string, Xs), c)
end

function ncdl{S <: AbstractString}(Xstrs::Vector{S}, c::Compressor = LibzCompressor(), return_indexset = false)
  Gxs = map(x -> clen(c, x), Xstrs)
  n = length(Xstrs)

  maxncd1 = -10.0 # Start with value < 0.0 to ensure first one is larger (NCD values are between 0.0 and 1+epsilon).

  # We will deselect one element in each step but start with all selected:
  selected = collect(1:n) .> 0
  selectedindices = Set(collect(1:n))
  numselected = n
  best = largest = 0

  while numselected >= 2

    # Calc max G(X \ {x}) for the still selected ones
    maxGxminus1 = -Inf
    for i in selectedindices
      selected[i] = false
      newGxminus1 = clen(c, join(Xstrs[selected], ""))
      selected[i] = true
      if newGxminus1 > maxGxminus1
        maxGxminus1 = newGxminus1
        largest = i
      end
    end

    # NCD1 = (G(X) - minimum(G(X))) / maxGxminus1
    ncd1 = (clen(c, join(Xstrs[selected], "")) - minimum(Gxs[selected])) / maxGxminus1
    #println("$numselected: $ncd1")
    if ncd1 > maxncd1
      maxncd1 = ncd1
      best = copy(selectedindices)
    end
    selected[largest] = false
    selectedindices = setdiff(selectedindices, Set(largest))
    numselected -= 1

  end

  if return_indexset
    return (maxncd1, best)
  else
    return maxncd1
  end
end


text:
## Mesh interface functions that we think, which should be supported
# to plugin Mesh library to finch
# this is mostly influenced by the dendro loop strucutres. But We can modify this as we go along.

# dimention of the simplex
@enum SIMPLEX_DIM begin
    D0=0,
    D1=1,
    D2=2,
    D3=3,
    D4=4
end

# element (support upto 4D elements) loop types
@enum SIMPLEX4_LOOP_TYPE begin
    ALL         =0,
    LOCAL       =1,
    HALO        =2,
    INDEP       =3,
    DEP         =4,
    WRITABLE    =5
end

# simplex of dim 3
@enum SIMPLEX3_LOOP_TYPE begin
    ALL         =0,
    LOCAL       =1,
    HALO        =2,
    INDEP       =3,
    DEP         =4,
    WRITABLE    =5,
    BOUNDRY     =6
end

# simplex of dim 2
@enum SIMPLEX2_LOOP_TYPE begin
    ALL         =0,
    LOCAL       =1,
    HALO        =2,
    INDEP       =3,
    DEP         =4,
    WRITABLE    =5,
    BOUNDRY     =6
end

# simplex of dim 1
@enum SIMPLEX1_LOOP_TYPE begin
    ALL         =0,
    LOCAL       =1,
    HALO        =2,
    INDEP       =3,
    DEP         =4,
    WRITABLE    =5,
    BOUNDRY     =6
end

# simplex of dim 0 (vertices)
@enum SIMPLEX0_LOOP_TYPE begin
    ALL         =0,
    LOCAL       =1,
    HALO        =2,
    INDEP       =3,
    DEP         =4,
    WRITABLE    =5,
    BOUNDRY     =6
end

# init function to initialize the iterator loop.
function _init{LoopType}() where {LoopType}
end

# increment to the next iterator.
function _next{LoopType}() where {LoopType}
end

# returns the end iterator
function _end{LoopType}() where {LoopType}
end

# get the current iterator.
function _curr{LoopType}() where {LoopType}
end

# get the indices of the specifid iterator.
# curr: current simplex iterator.
# dim: get the indices of simplex dim for current simplex
# example: if you are iterating over D3 simplex,
# you can get the indices of D2(surface), D1(edges), D0(vertices)

#=
function get_indices{Iter}(Iter curr, SIMPLEX_DIM dim) where {Iter}
#end

# dim : dim of the geometric object corresponding to the vector (4d, 3d, 2d, 1d, 0d),
       #if it is 0 then it will be nodal vector, 1 then edge vector, 2 surface vector etc.
# with_halo : allocate mem. for halo regions as well.
# dof : degrees of freedoms
function alloc_vector{T}(SIMPLEX_DIM dim, with_halo=false, dof=1) where {T}
end

# deallocation for the vector.
function dealloc_vector{Vec}(Vec vec) where {Vec}
end

# perform the halo sync. between proc. for a given vector (async begin non blocking call).
function sync_vector_begin{Vec}(Vec vec, with_halo=false, dof=1) where {Vec}
end

# perform the halo sync. between proc. for a given vector (async end non blocking call).
function sync_vector_end{Vec}(Vec vec, with_halo=false, dof=1) where {Vec}
end

=#

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
"""
    GammaExponential(lnℓ::AbstractFloat, lnγ::AbstractFloat)
    <: RadialKernel{Euclidean}

The gamma exponential kernel with length scale
``\\exp(\\gamma\\ln\\ell) = \\ell^\\gamma > 0``.

```math
k_{\\gamma,\\ell}(x, y) =
k_{\\gamma,\\ell}(\\lVert x - y\\rVert_2) =
\\exp\\bigg\\lbrace -
\\frac{\\lVert x - y\\rVert_2^\\gamma}{\\ell^\\gamma}\\bigg\\rbrace
```
"""
mutable struct GammaExponential{F<:AbstractFloat} <: RadialKernel{Euclidean}
    dist::Euclidean
    lnℓ::F
    lnγ::F

    ℓᵞ::F
    γ::F
    function GammaExponential(
        dist::Euclidean,
        lnℓ::AbstractFloat,
        lnγ::AbstractFloat
    )
        lnℓ, lnγ = promote(lnℓ, lnγ)
        γ = exp(lnγ)
        new{typeof(lnℓ)}(dist, lnℓ, lnγ, exp(γ * lnℓ), γ)
    end
end

GammaExponential(lnℓ::AbstractFloat, lnγ::AbstractFloat) =
    GammaExponential(Euclidean(), lnℓ, lnγ)

@inline (k::GammaExponential)(τ::AbstractFloat) = exp(-τ^k.γ / k.ℓᵞ)

numparams(::GammaExponential) = (1, 1)
paramtypes(::GammaExponential{F}) where F = (F, F)
params(k::GammaExponential) = (lnℓ = k.lnℓ, lnγ = k.lnγ)
function setparams!(k::GammaExponential{F}, lnℓ::F, lnγ::F) where F
    lnℓ, lnγ = promote(lnℓ, lnγ)
    γ = exp(lnγ)
    k.lnℓ = lnℓ
    k.lnγ = lnγ
    k.ℓᵞ = exp(γ * lnℓ)
    k.γ = γ
end


text:
module bertoldi_compression_examples

using FinEtools
using FinEtools.MeshImportModule
using FinEtools.MeshImportModule: import_ABAQUS
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed2DStrain
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.MatDeforNeohookeanADModule: MatDeforNeohookeanAD
using FinEtoolsDeforNonlinear.MatDeforStVKModule: MatDeforStVK
using FinEtoolsDeforNonlinear.MatDeforStVKADModule: MatDeforStVKAD
using FinEtoolsDeforNonlinear.MatDeforI1RivlinADModule: MatDeforI1RivlinAD
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.AlgoDeforNonlinearModule: nonlinearstatics
using LinearAlgebra: norm
using Statistics: mean
using SparseArrays
using DelimitedFiles
using Interpolations
using UnicodePlots

filetag(v) = Float64(Int(round(v*100000))) / 100000

function neohookeanad_q4()
    mr = DeforModelRed2DStrain
    E, nu = 7.0*phun("MPa"), 0.4772727
    m = MatDeforNeohookeanAD(mr, E, nu)
    L= 101.0*phun("mm");
    W = 2/2*phun("mm");
    H = 2/2*phun("mm");
    umag = 10.1*phun("mm");# Magnitude of the displacement
    nincr = 48
    utol = 10e-7;
    graphics = ~true;
    maxdu_tol = W/1e7;
    maxiter = 15
    tolerance = 1.0 * phun("mm") / 10

    output = MeshImportModule.import_ABAQUS("./bertoldi-specimen-1-mesh-0.002.inp")
    fens, fes = output["fens"], output["fesets"][1]
    fens.xyz = fens.xyz[:, 1:2]
    box = boundingbox(fens.xyz)
    fens2, fes2 = mirrormesh(fens, fes, vec([-1.0 0.0]),
    	box[[1, 3]]; renumb = c -> c[end:-1:1])
    fens, fes1, fes2 = mergemeshes(fens, fes,  fens2, fes2, tolerance)
    fes = cat(fes1, fes2)
    fens2, fes2 = mirrormesh(fens, fes, vec([0.0 -1.0]),
    	box[[1, 3]]; renumb = c -> c[end:-1:1])
    fens, fes1, fes2 = mergemeshes(fens, fes,  fens2, fes2, tolerance)
    fes = cat(fes1, fes2)
    box = boundingbox(fens.xyz)

    l1  = selectnode(fens; box = [-Inf,Inf,box[3],box[3]], inflate  =  tolerance)
    e1 = FDataDict("node_list"=>l1, "component"=>1)
    e2 = FDataDict("node_list"=>l1, "component"=>2)
    l2  = selectnode(fens; box = [-Inf,Inf,box[4],box[4]], inflate  =  tolerance)
    e3 = FDataDict("node_list"=>l2, "component"=>1)
    table = LinearInterpolation([0, 1], umag*[0,-1])
    move(x, lambda) = table(lambda);
    movel1 = l2
    e4 = FDataDict("node_list"=>l2, "component"=>2, "displacement"=>move)
    l3  = selectnode(fens; box = [box[1],box[1],-Inf,Inf], inflate  =  tolerance)
    l4  = selectnode(fens; box = [box[2],box[2],-Inf,Inf], inflate  =  tolerance)

	xyz = deepcopy(fens.xyz)
    vtkexportmesh("bertoldi_compression-boundary.vtk", fens, meshboundary(fes))

# Perturbed the locations of the nodes in the interior
    fens.xyz .+= 0.1 .* phun("mm") .* 2.0 .* (rand(size(fens.xyz)...) .- 0.5)
    # Restore the coordinates of the points on the outer boundary
    fens.xyz[l1, :] = xyz[l1, :]
    fens.xyz[l2, :] = xyz[l2, :]
    fens.xyz[l3, :] = xyz[l3, :]
    fens.xyz[l4, :] = xyz[l4, :]

    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(2, 2), H), m)

    region1 = FDataDict("femm"=>femm)
    modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1],  "essential_bcs"=>[e1, e2, e3, e4])

    modeldata["load_multipliers"] = (1:nincr)./nincr*1.0;
    modeldata["maxdu_tol"] = maxdu_tol;
    modeldata["maxiter"] = maxiter;
    modeldata["line_search"]  = true;
    modeldata["iteration_observer"] = (lambda, iter, du, modeldata) -> begin
        @show iter, modeldata["maxdu"], modeldata["maxbal"]
    end
    Ux = FFlt[]; Rx = FFlt[]
    modeldata["increment_observer"] = (lambda, modeldata) -> begin
    	@show lambda
        push!(Ux, mean(modeldata["un1"].values[movel1,2]));
        push!(Rx, sum(modeldata["reactions"].values[movel1,2]));
        vtkexportmesh("bertoldi_compression-$(filetag(lambda)).vtk", fens, fes; vectors = [("u", modeldata["un1"].values)])
    end

    modeldata = nonlinearstatics(modeldata);
    # @show Ux, Rx
    @show (minimum(Ux), minimum(Rx), maximum(Ux), maximum(Rx))

    pl = lineplot(Ux ./ L, Rx, canvas = DotCanvas)
    display(pl)
end # function neohookeanad_q4

function neohookeanad_q8()
    mr = DeforModelRed2DStrain
    E, nu = 7.0*phun("MPa"), 0.4772727
    m = MatDeforNeohookeanAD(mr, E, nu)
    L= 101.0*phun("mm");
    W = 2/2*phun("mm");
    H = 2/2*phun("mm");
    umag = 10.1*phun("mm");# Magnitude of the displacement
    nincr = 48
    utol = 10e-7;
    graphics = ~true;
    maxdu_tol = W/1e7;
    maxiter = 15
    tolerance = 1.0 * phun("mm") / 10

    output = MeshImportModule.import_ABAQUS("./bertoldi-specimen-1-mesh-0.002.inp")
    fens, fes = output["fens"], output["fesets"][1]
    fens.xyz = fens.xyz[:, 1:2]
    box = boundingbox(fens.xyz)
    fens2, fes2 = mirrormesh(fens, fes, vec([-1.0 0.0]),
    	box[[1, 3]]; renumb = c -> c[end:-1:1])
    fens, fes1, fes2 = mergemeshes(fens, fes,  fens2, fes2, tolerance)
    fes = cat(fes1, fes2)
    fens2, fes2 = mirrormesh(fens, fes, vec([0.0 -1.0]),
    	box[[1, 3]]; renumb = c -> c[end:-1:1])
    fens, fes1, fes2 = mergemeshes(fens, fes,  fens2, fes2, tolerance)
    fes = cat(fes1, fes2)
    box = boundingbox(fens.xyz)

    fens, fes = Q4toQ8(fens, fes)

    l1  = selectnode(fens; box = [-Inf,Inf,box[3],box[3]], inflate  =  tolerance)
    e1 = FDataDict("node_list"=>l1, "component"=>1)
    e2 = FDataDict("node_list"=>l1, "component"=>2)
    l2  = selectnode(fens; box = [-Inf,Inf,box[4],box[4]], inflate  =  tolerance)
    e3 = FDataDict("node_list"=>l2, "component"=>1)
    table = LinearInterpolation([0, 1], umag*[0,-1])
    move(x, lambda) = table(lambda);
    movel1 = l2
    e4 = FDataDict("node_list"=>l2, "component"=>2, "displacement"=>move)
    l3  = selectnode(fens; box = [box[1],box[1],-Inf,Inf], inflate  =  tolerance)
    l4  = selectnode(fens; box = [box[2],box[2],-Inf,Inf], inflate  =  tolerance)

	xyz = deepcopy(fens.xyz)
    vtkexportmesh("bertoldi_compression-boundary.vtk", fens, meshboundary(fes))

    # Perturb the locations of the nodes in the interior
    fens.xyz .+= 0.05 .* phun("mm") .* 2.0 .* (rand(size(fens.xyz)...) .- 0.5)
    # Restore the coordinates of the points on the outer boundary
    fens.xyz[l1, :] = xyz[l1, :]
    fens.xyz[l2, :] = xyz[l2, :]
    fens.xyz[l3, :] = xyz[l3, :]
    fens.xyz[l4, :] = xyz[l4, :]
    # The mid-edge nodes should be at the midpoint
    for i in 1:count(fes)
    	k = fes.conn[i][5]
    	p = fes.conn[i][1]
    	q = fes.conn[i][2]
    	fens.xyz[k, :] .= (fens.xyz[p, :] + fens.xyz[q, :]) / 2
    	k = fes.conn[i][6]
    	p = fes.conn[i][2]
    	q = fes.conn[i][3]
    	fens.xyz[k, :] .= (fens.xyz[p, :] + fens.xyz[q, :]) / 2
    	k = fes.conn[i][7]
    	p = fes.conn[i][3]
    	q = fes.conn[i][4]
    	fens.xyz[k, :] .= (fens.xyz[p, :] + fens.xyz[q, :]) / 2
    	k = fes.conn[i][8]
    	p = fes.conn[i][4]
    	q = fes.conn[i][1]
    	fens.xyz[k, :] .= (fens.xyz[p, :] + fens.xyz[q, :]) / 2
    end

    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(2, 2), H), m)

    region1 = FDataDict("femm"=>femm)
    modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1],  "essential_bcs"=>[e1, e2, e3, e4])

    modeldata["load_multipliers"] = (1:nincr)./nincr*1.0;
    modeldata["maxdu_tol"] = maxdu_tol;
    modeldata["maxiter"] = maxiter;
    modeldata["line_search"]  = true;
    modeldata["iteration_observer"] = (lambda, iter, du, modeldata) -> begin
        @show iter, modeldata["maxdu"], modeldata["maxbal"]
    end
    Ux = FFlt[]; Rx = FFlt[]
    modeldata["increment_observer"] = (lambda, modeldata) -> begin
    	@show lambda
        push!(Ux, mean(modeldata["un1"].values[movel1,2]));
        push!(Rx, sum(modeldata["reactions"].values[movel1,2]));
        vtkexportmesh("bertoldi_compression-$(filetag(lambda)).vtk", fens, fes; vectors = [("u", modeldata["un1"].values)])
    end

    modeldata = nonlinearstatics(modeldata);
    # @show Ux, Rx
    @show (minimum(Ux), minimum(Rx), maximum(Ux), maximum(Rx))

    pl = lineplot(Ux ./ L, Rx, canvas = DotCanvas)
    display(pl)
end # function neohookeanad_q8

function allrun()
    println("#####################################################")
    println("# neohookeanad_q4 ")
    neohookeanad_q4()
    println("#####################################################")
    println("# neohookeanad_q8 ")
    neohookeanad_q8()
    return true
end # function allrun

end # module tension_compression_examples


text:
"""
    inflation_factor(m::LinearMixedModel, blups=ranef(m), resids=residuals(m)))

Compute how much the standard deviation of the BLUPs/conditional modes and residuals
needs to be inflated in order to match the (restricted) maximum-likelihood estimate.

Due to the shrinkage associated with the random effects in a mixed model, the
variance of the conditional modes / BLUPs / random intercepts and slopes is less
than the variance estimated by the model and displayed in the model summary or
via `MixedModels.VarCorr`. This shrinkage also impacts the observational level
residuals. To compensate for this, the resampled residuals and groups need to be
scale-inflated so that their standard deviation matches that of the estimates in
the original model.

The factor for scale inflation is returned as a `Vector` of inflation factors
for each of the random-effect strata (blocking variables) and the observation-level
variability. The factor is on the standard deviation scale (lower Cholesky factor
in the case of vector-valued random effects).
"""
function inflation_factor(m::LinearMixedModel, blups=ranef(m), resids=residuals(m))
# FIXME I'm not sure this is correct
#       the nonparametric bootstrap underestimates variance components
#       compared to the parametricbootstrap

    σ = sdest(m)
    σres = std(resids; corrected=false)
    inflation = map(zip(m.reterms, blups)) do (trm, re)
        # inflation
        λmle =  trm.λ * σ                               # L_R in CGR
        λemp = cholesky(cov(re'; corrected=false)).L    # L_S in CGR
        # no transpose because the RE are transposed relativ to CGR
        λmle \ λemp
    end

    return [inflation; σ / σres]
end


text:
include("/Users/ericagol/Computer/Julia/regress.jl")

n = 8
include("ttv.jl")
t0 = 7257.93115525
#h  = 0.12
h  = 0.05
#h  = 0.025
tmax = 600.0
#tmax = 80.0

# Read in initial conditions:
elements = readdlm("../test/elements.txt",',')

# Make an array, tt,  to hold transit times:
# First, though, make sure it is large enough:
ntt = zeros(Int64,n)
for i=2:n
  ntt[i] = ceil(Int64,tmax/elements[i,2])+3
end
println("ntt: ",ntt)
tt  = zeros(n,maximum(ntt))
tt1 = zeros(n,maximum(ntt))
tt2 = zeros(n,maximum(ntt))
tt3 = zeros(n,maximum(ntt))
tt4 = zeros(n,maximum(ntt))
# Save a counter for the actual number of transit times of each planet:
count = zeros(Int64,n)
count1 = zeros(Int64,n)
# Call the ttv function:
rstar = 1e12
dq = ttv_elements!(n,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar)
@time dq = ttv_elements!(n,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar)
# Now, try calling with kicks between planets rather than -drift+Kepler:
pair_input = ones(Bool,n,n)
# We want Keplerian between star & planets, and impulses between
# planets.  Impulse is indicated with 'true', -drift+Kepler with 'false':
for i=2:n
  pair_input[1,i] = false
  # We don't need to define this, but let's anyways:
  pair_input[i,1] = false
end
# Now, only include Kepler solver for adjacent planets:
for i=2:n-1
  pair_input[i,i+1] = false
  pair_input[i+1,i] = false
end

# Now call with smaller timestep:
count2 = zeros(Int64,n)
count3 = zeros(Int64,n)
dq = ttv_elements!(n,t0,h/4.,tmax,elements,tt2,count2,0.0,0,0,rstar)
for i=2:n
  println("Timing error -drift+Kepler: ",i-1," ",maximum(abs.(tt2[i,:]-tt1[i,:]))*24.*3600.)
end
@time dq = ttv_elements!(n,t0,h,tmax,elements,tt3,count1,0.0,0,0,rstar;pair=pair_input)
dq = ttv_elements!(n,t0,h/4.,tmax,elements,tt4,count2,0.0,0,0,rstar;pair=pair_input)
for i=2:n
  println("Timing error kickfast:           ",i-1," ",maximum(abs.(tt3[i,:]-tt4[i,:]))*24.*3600.)
end
for i=2:n
  println("Timing kickfast-(-drift+Kepler): ",i-1," ",maximum(abs.(tt1[i,:]-tt3[i,:]))*24.*3600.)
end

# Now, compute derivatives (with respect to initial cartesian positions/masses):
dtdq1 = zeros(n,maximum(ntt),7,n)
dtdq2 = zeros(n,maximum(ntt),7,n)
dtdq_kick1 = zeros(n,maximum(ntt),7,n)
dtdq_kick2 = zeros(n,maximum(ntt),7,n)
dtdelements1 = zeros(n,maximum(ntt),7,n)
dtdelements2 = zeros(n,maximum(ntt),7,n)
dtdelements_kick1 = zeros(n,maximum(ntt),7,n)
dtdelements_kick2 = zeros(n,maximum(ntt),7,n)

dtdelements1 = ttv_elements!(n,t0,h,tmax,elements,tt,count,dtdq1,rstar)
@time dtdelements1 = ttv_elements!(n,t0,h,tmax,elements,tt,count,dtdq1,rstar)
dtdelements2 = ttv_elements!(n,t0,h/2,tmax,elements,tt,count,dtdq2,rstar)
@time dtdelements_kick1 = ttv_elements!(n,t0,h,tmax,elements,tt,count,dtdq_kick1,rstar;pair=pair_input)
@time dtdelements_kick2 = ttv_elements!(n,t0,h/2,tmax,elements,tt,count,dtdq_kick2,rstar;pair=pair_input)
println(maximum(abs.(asinh.(dtdelements1)-asinh.(dtdelements2))))
println(maximum(abs.(asinh.(dtdelements_kick1)-asinh.(dtdelements_kick2))))
println(maximum(abs.(asinh.(dtdelements2)-asinh.(dtdelements_kick2))))

Profile.clear()
Profile.init(10^7,0.01)
#@profile dtdelements0 = ttv_elements!(n,t0,h,tmax,elements,tt,count,dtdq0,rstar);
@profile dtdelements0 = ttv_elements!(n,t0,h,tmax,elements,tt,count,dtdq1,rstar);
Profile.print()


text:
using Random
using MLJBase
using Test

mutable struct Foo <: MLJType
    rng::AbstractRNG
    x::Int
    y::Int
end

mutable struct Bar <: MLJType
    rng::AbstractRNG
    x::Int
    y::Int
end

f1 = Foo(MersenneTwister(7), 1, 2)
f2 = Foo(MersenneTwister(8), 1, 2)
@test f1.rng != f2.rng
@test f1 == f2
f1.x = 10
@test f1 != f2
b = Bar(MersenneTwister(7), 1, 2)
@test f2 != b

@test is_same_except(f1, f2, :x)
f1.y = 20
@test f1 != f2
@test is_same_except(f1, f2, :x, :y)

# test for nested fields

mutable struct Super <: MLJType
    sub::Foo
    z::Int
end

f1 = Foo(MersenneTwister(7), 1, 2)
f2 = Foo(MersenneTwister(8), 1, 2)
s1 = Super(f1, 20)
s2 = Super(f2, 20)
@test s1 == s2
s2.sub.x = 10
@test f1 != f2

true


text:
#=
# TODO: Infer dimensionality of user function from context type
struct RelaxMeta{N,T}
    p::SVector{N,T}
    pL::SVector{N,T}
    pU::SVector{N,T}
    v::VariableValues{Float64}
    s::Vector{Int}
    ϵ::Float64
    subgradient_refinement::Bool
end

@context RelaxCtx

function ctx_subgradient_refine(c::RelaxCtx{RelaxMeta{N,T},Nothing}, z::MC{N,T}) where {N,T}
    m = ctx.metadata
    v = ctx.metadata.v
    s = ctx.metadata.s
    ϵ = ctx.metadata.ϵ
    return m.subgradient_refinement ? set_value_post(z, m.v, m.s, m.ϵ), z
end

for k in UNIVARIATE_ATOM_TYPES
    f = UNIVARIATE_ATOM_DICT[k]
    @eval function Cassette.overdub(ctx::RelaxCtx{RelaxMeta{N,T},Nothing}, ::typeof($f), x::MC{N,T}) where {N,T}
        return ctx_subgradient_refine(ctx, ($f)(x))
    end
end

for k in BIVARIATE_ATOM_TYPES
    f = BIVARIATE_ATOM_DICT[k]
    @eval function Cassette.overdub(ctx::RelaxCtx{RelaxMeta{N,T},Nothing}, ::typeof($f), x::MC{N,T}, y::NumberNotRelax) where {N,T}
        return ctx_subgradient_refine(ctx, ($f)(x, y))
    end
    @eval function Cassette.overdub(ctx::RelaxCtx{RelaxMeta{N,T},Nothing}, ::typeof($f), y::NumberNotRelax, x::MC{N,T}) where {N,T}
        return ctx_subgradient_refine(ctx, ($f)(y, x))
    end
end

#=
for k in NARITY_ATOM_TYPES
    f = NARITY_ATOM_DICT[k]
    @eval function Cassette.overdub(ctx::RelaxCtx{RelaxMeta{N,T},Nothing}, ::typeof($f), ...) where {N,T}
    end
end
=#

# commonly used storage objects automatically get promoted to use the set-valued storage
Cassette.overdub(::F_RELAX_CTX{T}, ::typeof(zeros), n::Int) where {T<:Number} = zeros(T, n)
Cassette.overdub(::F_RELAX_CTX{T}, ::typeof(zeros), dims...) where {T<:Number} = zeros(T, dims...)
Cassette.overdub(::F_RELAX_CTX{T}, ::typeof(fill), v::Any, dims::Int...) where {T<:Number} = fill(convert(T,v), dims)
function Cassette.overdub(::F_RELAX_CTX{T}, ::typeof(hcat), A::AbstractArray, Bs::AbstractArray...) where {T<:Number}
    vA = hcat(A,Bs...)
    vR = zeros(T, size(vA)...)
    return copyto!(vR,vA)
end
function Cassette.overdub(::F_RELAX_CTX{T}, ::typeof(vcat), A::AbstractArray, Bs::AbstractArray...) where {T<:Number}
    vA = vcat(A,Bs...)
    vR = zeros(T, size(vA)...)
    return copyto!(vR,vA)
end
=#
#=
# solution object
function Cassette.overdub(::F_RELAX_CTX{T}, ::typeof(ODEProblem), f, u0, tspan) where {T<:Number}
end
function Cassette.overdub(::F_RELAX_CTX{T}, ::typeof(solve), prob, alg, kwags...) where {T<:Number}
end
=#

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
# test to check if preGen helps reduce the solution time
using Distributed;
addprocs(20);
@everywhere include("loadMod.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();

N = 5;
dataList = Dict();

caseList = [13,33,123];
T = 96;
τ = Int64(1/6*T);
Δt = 0.25;
iterMax = 20;
pathTrain = load("pathHist_600.jld");

for ci in 1:length(caseList)
    pathDict = pathTrain["pathDict"][ci][T];
    for j in procs()
        remotecall_fetch(readInData,j,ci,caseList,T,τ);
    end

    cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(T, Δt, 30, false,false, 2, 2);

    startPGT = time();
    cutDictPG = preGen(T, Δt, N, iterMax);
    preGenT = time() - startPGT;

    startT = time();
    # cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, true, false, 10, 10, cutDictPG,
    #     false, 200, [], 0, pathDict);
    cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(T, Δt, N, false, false,
        max(Int64(round(500/N)),20), max(Int64(round(500/N)),20), cutDictPG, false, 200, [], 0, pathDict);
    elapsedT = time() - startT;
    dataList["preGen"] = [LBHist,UBHist,UBuHist,UBlHist,timeHist,elapsedT,preGenT];

    save("pgResults_$(ci).jld","data",dataList);
end


text:
module TracerAdvectionDiffusion

export
   Problem,
   set_c!,
   updatevars!

using
   DocStringExtensions,
   Reexport

@reexport using FourierFlows

import LinearAlgebra: mul!, ldiv!

# --
# Problems
# --

"""
    Problem(; parameters...)

Construct a constant diffusivity problem with steady or time-varying flow.
"""
noflow(args...) = 0.0 # used as defaults for u, v functions in Problem()

function Problem(dev;
          nx = 128,
          Lx = 2π,
          ny = nx,
          Ly = Lx,
           κ = 0.1,
           η = κ,
           u = noflow,
           v = noflow,
          dt = 0.01,
     stepper = "RK4",
  steadyflow = false,
           T = Float64
  )

  grid = TwoDGrid(dev, nx, Lx, ny, Ly; T=T)
  params = steadyflow==true ? ConstDiffSteadyFlowParams(η, κ, u, v, grid) : ConstDiffParams(η, κ, u, v)
  vars = Vars(dev, grid)
  equation = Equation(params, grid)

  return FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev)
end

# --
# Params
# --

abstract type AbstractTracerParams <: AbstractParams end
abstract type AbstractConstDiffParams <: AbstractParams end
abstract type AbstractSteadyFlowParams <: AbstractParams end

"""
    struct ConstDiffParams{T} <: AbstractConstDiffParams

A struct containing the parameters for a constant diffusivity problem with time-varying flow.
Included are:

$(TYPEDFIELDS)
"""
struct ConstDiffParams{T} <: AbstractConstDiffParams
    "isotropic horizontal diffusivity coefficient"
       η :: T
    "isotropic vertical diffusivity coefficient"
       κ :: T
    "isotropic hyperdiffusivity coefficient"
      κh :: T
    "isotropic hyperdiffusivity order"
     nκh :: Int
    "function returning the x-component of advecting flow"
       u :: Function
    "function returning the y-component of advecting flow"
       v :: Function
end

"""
    ConstDiffParams(η, κ, u, v)

The constructor for the `params` struct for constant diffusivity problem and time-varying flow.
"""
ConstDiffParams(η, κ, u, v) = ConstDiffParams(η, κ, 0η, 0, u, v)

"""
    struct ConstDiffParams{T} <: AbstractConstDiffParams

A struct containing the parameters for a constant diffusivity problem with steady flow.
Included are:

$(TYPEDFIELDS)
"""
struct ConstDiffSteadyFlowParams{T,A} <: AbstractSteadyFlowParams
  "isotropic horizontal diffusivity coefficient"
     η :: T
  "isotropic vertical diffusivity coefficient"
     κ :: T
  "isotropic hyperdiffusivity coefficient"
    κh :: T
  "isotropic hyperdiffusivity order"
   nκh :: Int
   "x-component of advecting flow"
      u :: A
   "y-component of advecting flow"
      v :: A
end

"""
    ConstDiffSteadyFlowParams(η, κ, κh, nκh, u::Function, v::Function, grid)
    ConstDiffSteadyFlowParams(η, κ, u, v, grid)

The constructor for the `params` struct for constant diffusivity problem and steady flow.
"""
function ConstDiffSteadyFlowParams(η, κ, κh, nκh, u::Function, v::Function, grid)
   x, y = gridpoints(grid)
  ugrid = u.(x, y)
  vgrid = v.(x, y)

  return ConstDiffSteadyFlowParams(η, κ, κh, nκh, ugrid, vgrid)
end

ConstDiffSteadyFlowParams(η, κ, u, v, grid) = ConstDiffSteadyFlowParams(η, κ, 0η, 0, u, v, grid)

# --
# Equations
# --

"""
    Equation(params, grid)

Returns the equation for constant diffusivity problem with params p and grid g.
"""
function Equation(params::ConstDiffParams, grid)
  L = @. -params.η * grid.kr^2 - params.κ * grid.l^2 - params.κh * grid.Krsq^params.nκh

  return FourierFlows.Equation(L, calcN!, grid)
end

function Equation(params::ConstDiffSteadyFlowParams, grid)
  L = @. -params.η * grid.kr^2 - params.κ * grid.l^2 - params.κh * grid.Krsq^params.nκh

  return FourierFlows.Equation(L, calcN_steadyflow!, grid)
end

# --
# Vars
# --

"""
    struct Vars{Aphys, Atrans} <: AbstractVars

The variables for TracerAdvectionDiffussion problems.

$(FIELDS)
"""
struct Vars{Aphys, Atrans} <: AbstractVars
    "tracer concentration"
       c :: Aphys
    "tracer concentration x-derivative, ∂c/∂x"
      cx :: Aphys
    "tracer concentration y-derivative, ∂c/∂y"
      cy :: Aphys
    "Fourier transform of tracer concentration"
      ch :: Atrans
    "Fourier transform of tracer concentration x-derivative, ∂c/∂x"
     cxh :: Atrans
    "Fourier transform of tracer concentration y-derivative, ∂c/∂y"
     cyh :: Atrans
end

"""
    Vars(dev, grid)

Returns the variables `vars` for a constant diffusivity problem on `grid` and device `dev`.
"""
function Vars(::Dev, grid::AbstractGrid{T}) where {Dev, T}
  @devzeros Dev T (grid.nx, grid.ny) c cx cy
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) ch cxh cyh

  return Vars(c, cx, cy, ch, cxh, cyh)
end

# --
# Solvers
# --

"""
    calcN!(N, sol, t, clock, vars, params, grid)

Calculate the advective terms for a tracer equation with constant diffusivity and time-varying flow.
"""
function calcN!(N, sol, t, clock, vars, params::AbstractConstDiffParams, grid)
  @. vars.cxh = im * grid.kr * sol
  @. vars.cyh = im * grid.l  * sol

  ldiv!(vars.cx, grid.rfftplan, vars.cxh) # destroys vars.cxh when using fftw
  ldiv!(vars.cy, grid.rfftplan, vars.cyh) # destroys vars.cyh when using fftw

  x, y = gridpoints(grid)
  @. vars.cx = -params.u(x, y, clock.t) * vars.cx - params.v(x, y, clock.t) * vars.cy # copies over vars.cx so vars.cx = N in physical space
  mul!(N, grid.rfftplan, vars.cx)

  return nothing
end

"""
    calcN_steadyflow!(N, sol, t, clock, vars, params, grid)

Calculate the advective terms for a tracer equation with constant diffusivity and time-constant flow.
"""
function calcN_steadyflow!(N, sol, t, clock, vars, params::AbstractSteadyFlowParams, grid)
  @. vars.cxh = im * grid.kr * sol
  @. vars.cyh = im * grid.l  * sol

  ldiv!(vars.cx, grid.rfftplan, vars.cxh) # destroys vars.cxh when using fftw
  ldiv!(vars.cy, grid.rfftplan, vars.cyh) # destroys vars.cyh when using fftw

  @. vars.cx = -params.u * vars.cx - params.v * vars.cy # copies over vars.cx so vars.cx = N in physical space
  mul!(N, grid.rfftplan, vars.cx)

  return nothing
end

# --
# Helper functions
# --

"""
    updatevars!(prob)

Update the vars in v on the grid g with the solution in sol.
"""
function updatevars!(prob)
  vars, grid, sol = prob.vars, prob.grid, prob.sol

  @. vars.ch = sol

  ldiv!(vars.c, grid.rfftplan, deepcopy(vars.ch))

  return nothing
end

"""
    set_c!(prob, c)

Set the solution sol as the transform of c and update variables v
on the grid g.
"""
function set_c!(prob, c)
  sol, vars, grid = prob.sol, prob.vars, prob.grid

  mul!(sol, grid.rfftplan, c)

  updatevars!(prob)

  return nothing
end

end # module

max_stars_repo_name: FourierFlows/PassiveTracerFlows

max_issues_repo_name: FourierFlows/PassiveTracerFlows

max_forks_repo_name: FourierFlows/PassiveTracerFlows


text:
using Pkg; Pkg.activate(@__DIR__)
using Comrade
using PyCall
using Plots

using Zygote
using Distributions, DistributionsAD
using ComradeAHMC
using ComradeGalactic
using GalacticBBO
using ForwardDiff
using GalacticOptimJL

load_ehtim()
@pyimport eht_dmc as ed

# To download the data visit https://doi.org/10.25739/g85n-f134
obs = ehtim.obsdata.load_uvfits(joinpath(@__DIR__, "SR1_M87_2017_096_hi_hops_netcal_StokesI.uvfits"))
obs.add_scans()

obsavg = obs.flag_uvdist(uv_min=0.1e9).avg_coherent(0.0, scan_avg=true).add_fractional_noise(0.02)

obslist = obsavg.split_obs()
obss = ehtim.obsdata.merge_obs(obslist)
damp = extract_amp(obss)
dcp = extract_cphase(obss)
st = scantable(damp)

gcache = GainCache(st)

struct Model{G,C}
    gcache::G
    cache::C
    fovx::Float64
    fovy::Float64
    nx::Int
    ny::Int
end

function (model::Model)(θ)
    (;c, gamp) = θ
    g = exp.(gamp)
    cmat = reshape(c, model.ny, model.nx)
    img = IntensityMap(cmat, model.fovx, model.fovy, BSplinePulse{3}())

    m = modelimage(img, model.cache)
    return GainModel(model.gcache, g, m)
end

# define the priors
distamp = (AA = Normal(0.0, 0.1),
         AP = Normal(0.0, 0.1),
         LM = Normal(0.0, 1.0),
         AZ = Normal(0.0, 0.1),
         JC = Normal(0.0, 0.1),
         PV = Normal(0.0, 0.1),
         SM = Normal(0.0, 0.1)
         )
distphase = (AA = Normal(0.0, 0.0001),
             AP = Normal(0.0, 1π),
             LM = Normal(0.0, 1π),
             AZ = Normal(0.0, 1π),
             JC = Normal(0.0, 1π),
             PV = Normal(0.0, 1π),
             SM = Normal(0.0, 1π)
             )
# Now form the posterior
npix = 12
prior = (
          c = MvLogNormal(fill(-log(4*npix*npix), npix*npix), fill(2.0, npix*npix)),
          #c = DistributionsAD.TuringDirichlet(npix*npix, 1.0),
          #f = Uniform(0.0, 1.0),
          gamp = Comrade.GainPrior(distamp, st),
        )

# Now form the posterior

fovxy = μas2rad(70.0)
cache = create_cache(Comrade.DFTAlg(damp), IntensityMap(zeros(npix,npix), fovxy, fovxy, BSplinePulse{3}()))
mms = Model(gcache, cache, fovxy, fovxy, npix, npix)

#mms = NoGModel(cache, fovxy, fovxy, npix, npix)
lklhd = RadioLikelihood(damp, dcp)
post = Posterior(lklhd, prior, mms)
tpost = asflat(post)
lp = logdensityof(tpost)

# We will use HMC to sample the posterior.
# First to reduce burn in we use pathfinder
ndim = dimension(tpost)
f = OptimizationFunction(tpost, GalacticOptim.AutoZygote())
x0 = Comrade.HypercubeTransform.inverse(tpost, prior_sample(post))

#prob = OptimizationProblem(f, zeros(ndim), nothing, lb=fill(-5.0, ndim), ub=fill(5.0, ndim))
#sol = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited(); maxiters=1_000_000)

prob = OptimizationProblem(f, 0.5*randn(ndim), nothing)
prob = OptimizationProblem(f, sol.u, nothing)
function cb(x, args...)
    @info "lpost: $(lp(x))"
    return false
end
sol = solve(prob, LBFGS(); g_tol=1e-2, maxiters=2_000, callback=cb)

xopt = transform(tpost, sol)

residual(mms(xopt), dcp)
residual(mms(xopt), damp)

plot(mms(xopt), xlims=(-40.0,40.0), ylims=(-40.0,40.0))

# now we sample using hm
metric = DenseEuclideanMetric(ndim)
chain, stats = sample(post, AHMC(;metric), 6000; nadapts=4000, init_params=xopt)


text:
using ESDL
using OnlineStats
using WeightedOnlineStats
using Test
using DataFrames
using CSV
using StatsBase
using RDatasets
using Statistics

iris = dataset("datasets", "iris")

@testset "TableAggregator functions" begin
    fitMean = fittable(DataFrames.eachrow(iris), Mean, :SepalWidth)
    @test value(fitMean) ≈ mean(iris[!,:SepalWidth])

    minSepalLength, maxSepalLength = extrema(iris[!,:SepalLength])
    weightVector = map(x -> (x - minSepalLength) /
                    (maxSepalLength - minSepalLength), iris[!,:SepalLength]
                    )
    fitMeanWeight = fittable(DataFrames.eachrow(iris), WeightedMean, :SepalWidth,
                    weight=(x -> (x.SepalLength - minSepalLength) /
                        (maxSepalLength - minSepalLength)
                        )
                    )
    @test value(fitMeanWeight) ≈ mean(iris[!,:SepalWidth], weights(weightVector))

    # meanBy = combine(groupby(iris, :Species), SepalWidth => mean)
    # cn = names(meanBy)[2]
    # fitMeanBy = fittable(DataFrames.eachrow(iris), Mean, :SepalWidth, by=(:Species,))
    # for (key, value) in value(fitMeanBy)
    #     @test value ≈ meanBy[meanBy[:Species] .== key, cn][1]
    # end

    # meanByWeight = by(iris, :Species,
    #                     SepalWidth_mean = [:SepalWidth, :SepalLength] =>
    #                         (x -> mean(x.SepalWidth,
    #                             weights(map(x -> (x - minSepalLength) /
    #                                         (maxSepalLength - minSepalLength),
    #                                         x.SepalLength)
    #                                     )
    #                             )
    #                         )
    #                 )
    # fitMeanByWeight = fittable(DataFrames.eachrow(iris), WeightedMean,
    #                             :SepalWidth, by=(:Species,),
    #                             weight=(x -> (x.SepalLength - minSepalLength) /
    #                                         (maxSepalLength - minSepalLength)
    #                                     )
    #                             )
    # for (key, value) in value(fitMeanByWeight)
    #     @test value ≈ meanByWeight[meanByWeight[!,:Species] .== key, :SepalWidth_mean][1]
    # end
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.7.0
#     language: julia
#     name: julia-1.7
# ---

# %% [markdown]
# https://twitter.com/mkashi/status/1469153893179596801

# %%
using Plots

f(x) = x^3 - 3x^2 + 3x - 3
df(x) = 3x^2 - 6x + 3

x = 2.1
dx = @. 2.0 ^ (-60:4)

finitediff = @. (f(x + dx) - f(x)) / dx
abserr_finirediff = @. abs(finitediff - df(x))

plot(dx, abserr_finirediff; xscale=:log10, yscale=:log10, label="")
plot!(; xtick=@.(1e1^(-20:2:2)), ytick=@.(1e1^(-20:2:2)))
plot!(; xlabel="dx", ylabel="absolute error")

# %%
df(x)

# %%
using ForwardDiff
ForwardDiff.derivative(f, x)

# %%
@show finitediff;

# %%
using FiniteDifferences
@show [central_fdm(k, 1)(f, x) for k in 2:10];

# %%

max_stars_repo_path: 0025/finite difference and autimatic differentiation.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
module PowerSystemDataStructure

#using #dependencies

#import #methods to overload

export Bus, ConductorModel, ThermalModel, Line, Gen,
        Load, Wind, Storage, Params, Control,
        OptParams, CostFunction, Optimize, Time, System,
        mat2sys

# BY JONAS KERSULIS, Jan. 22 2015
# This module defines the types used to characterize a transmission system.

# FROM MADS'S WORK
# See runLevel3vFull.m.
# The majority of System is created in setup files which then call
# runLevel3vFull.m (e.g. testLevel3 which calls systemSetup_RTS96 and
# then calls runLevel3vFull).
# Fields which are marked with an asterisk '*' are initialized or added to
# System in runLevel3vFull.  Unmarked fields should be initialized prior to
# runLevel3vFull being called.

type Bus # (number of buses in network)x1
    num::Integer        # bus number
    ty::Integer         # bus type (1 = PQ, 2 = PV, 3 = slack)
    vBase               # base voltage in kV
    vMag                # voltage magnitude in puu
    vAng                # voltage angle in rad
    g                   # shunt conductance in pu
    b                   # shunt susceptance in pu
end

type ConductorModel # 1x1 structure containing...
    R_unit              # unit resistance of conductor [ohms/m]
    I_max               # current limit of network connection [A] (note that this is related to sLim field above. sLim = 0 or very large for voltage base will cause I_max to have a VERY large current rating (~50K) in order to approximate infinity.
    D                   # diameter of complete conductor [mm]
    Al_m                # unit mass of aluminum portion of ACSR [kg/m]
    St_m                # unit mass of steel portion of ACSR [kg/m]
    bundle              # number of conductors used in the network connection

    ConductorModel() = new()
end

type ThermalModel # 1x1 structure containing...
    T_amb               # ambient temperature [C]
    T_max               # ampacity-limited steady-state temperature [C] (related to sLim and I_max)
    η_c                 # conductive heat factor (varies w/ temp) [W/m-C]
    η_r                 # radiative heat factor [W/m-C^4]
    η_s                 # linear solar heat gain rate [W/m]
    mCp                 # line heat capacity per unit length (J/m-C)
    τ                   # (unitless) coefficient for temperature dynamics
    ρ                   # (C-m/W) coefficient for temperature dynamics
    γ                   # (unitless) coefficient for temperature dynamics
    fLossLim            # per-phase active power loss linearization point [W/m]

    ThermalModel() = new()
end

type Line # (number of lines in network)x1 structure
    # Unified branch model is chosen as:
    #         k  tmk:1    p     Zkm      q               m
    #         |----8 8------O-----ZZZ------O-------------|
    #   Tap   +             |              |             +   Fixed
    #   bus   Vk     ykm^sh Z              Z ykm^sh      Vm   bus
    #         -             |              |             -
    #        ---------------------------------------------
    # with tmk = tap*exp(j*phase) is pu (complex) tap-ratio
    # o when tmk = 1          ===> transmission line
    # o when tmk = real       ===> TCT ===> tap-ratio (~1) [assumes ykm^sh=0]
    # o when tmk = complex    ===> PST ===> angle (rads)   [assumes ykm^sh=0]

    from::Integer                   # from bus number
    to::Integer                     # to bus number
    r                               # resistance (pu)
    x                               # reactance (pu)
    b                               # shunt susceptance (pu)
    sLim                            # 3-phase power rating of line (MVA)
    length                          # length of line (miles) (Level 3 can estimate lengths if data not available)
    flow                            # active power flow on line (pu)
    fLoss                           # active power losses on line (pu)
    q                               # reactive power flow on line (pu)
    qLoss                           # reactive power losses on line (pu)
    delT                            # degrees above thermal line limit (C)
    isXfrmr::Bool                   # logical specifying if line is a transformer
    isTCXfrmr::Bool                 # logical specifying if transformer is tap-changing
    isPSXfrmr::Bool                 # logical specifying if transformer is phase-shifting
    tap                             # tap ratio (for Xfrmrs)
    phase                           # phase ratio (for Xfrmrs) (rad)
    maxPhase                        # maximum phase shift (rad)
    conductorModel::ConductorModel  # 1x1 structure containing...
    thermalModel::ThermalModel      # 1x1 structure containing...
    Δθ                              # width of PWL segments approximating abs(theta_{ij})
end

type Gen # (number of generators in network)x1 structure
    bus::Integer            # bus number of generator
    ID::Integer             # generator identifier (not needed)
    status::Integer         # 0 = out of service, 1 = in service (not needed: Level 3 assumes all gens are available)
    Pinj                    # active power injection (pu)
    Qinj                    # reactive power injection (pu)
    distFactor              # distribution factor (droop response)
    Pramp                   # active power ramping rate (pu/hr)
    vSet                    # voltage setpoint of generator (pu)
    Pmin                    # min active power injection (pu)
    Pmax                    # max active power injection (pu)
    Qmin                    # min reactive power injection (pu) (should be <= 0)
    Qmax                    # max reactive power injection (pu) (should be >= 0)
    pUpLim                  # active power ramp rate limit (pu/hr) (for increasing output)
    pDownLim                # active power ramp rate limit (pu/hr) (for decreasing output)
    qUpLim                  # reactive power ramp rate limit (pu/hr) (for increasing output)
    qDownLim                # reactive power ramp rate limit (pu/hr) (for decreasing output)
end

 type Load # (number of loads in network)x1 structure
    bus::Integer            # bus number of load
    Pload                   # active power demand forecast (pu) 1x(number of time intervals) vector
    Qload                   # reactive power demand (pu) 1x(number of time intervals) vector
    Pdemand                 # present active power load (pu)
    Qdemand                 # present reactive power load (pu)
    Pshed                   # active power shed from load through demand response (pu)
    Qshed                   # reactive power shed from load through demand response (pu)
end

type Wind # (number of wind farms in network)x1 structure
    bus         # bus number of wind generator
    Pforecast   # forecasted active power production (pu) 1x(number of time intervals) vector
    Pinj        # actual active power production (pu)
    Qinj
    Pshed       # active power shed from wind forecast through curtailment (pu)

    Wind() = new()
end

type Storage # (number of storage devices in network)x1 structure
    bus             # bus number of storage
    Pdemand         # active power being consumed by device [measured at bus] (pu)
    powerRate       # power rating (pu)
    energyRate      # energy rating (pu*hr)
    energyLvl       # energy currently stored in device (pu*hr)
    chargeEff       # charging efficiency [between 0 and 1]
    dischargeEff    # discharging efficiency [between 0 and 1]

    Storage() = new()
end

type Params # 1x1 structure
    sBase                   # 3-phase power system base (MVA)
    numOfBus::Integer       # No. of buses in system
    numOfLine::Integer      # No. of lines in system
    numOfGen::Integer       # No. of generators in system
    numOfLoad::Integer      # No. of loads in system
    numOfWind::Integer      # No. of wind generators in system
    numOfStorage::Integer   # No. of storage devices in system
    numOfTimeInt::Integer   # No. of time intervals (for load and wind forecasts)
    hasLosses::Integer      # integer indicating how many segments to use in piecewise-linear model of system losses. (0 indicates that PWL model is not used)
    phaseMax                # maximum phase angle difference allowed across a single line in the network [rad]
    NAI                     # * node arc incidence matrix (nbus x nline, 1 = from -1 = to)
    NGI                     # * node generator incidence matrix (nbus x ngen, 1 = incidence)
    NLI                     # * node load incidence matrix (nbus x nload, 1 = incidence)
    NWI                     # * node wind incidence matrix (nbus x nwind, 1 = incidence)
    NSI                     # * node storage incidence matrix (nbus x nstorage, 1 = incidence)
    phPW_b                  # * PWL approximation of abs(theta_{ij}) for each line

    # Allow partial initialization:
    Params() = new()
end

type Control # 1x1 structure
    T_s             # discrete sampling time (for MPC) [sec] (Must set to 60 for now)
    MPChorizon      # receeding prediction horizon for MPC
    MPCsafeTime     # number of MPC iterations needed without line trip for system to be considered safe
    MPCmaxTime      # maximum number of MPC iterations allowed
    lineOut         # nline x 1 logical specifying which lines are out of service
    storageTol      # scalar specifying tolerance around 0 for which storage devices are allowed to simulcharge
    overloadTol     # acceptable line flow overload allowed at safe MPC termination [%]

    Control() = new()
end

type OptParams # 1x1 structure containing GUROBI parameters
    outputflag      # decide whether GUROBI steps should display on console window (1=yes,0=no)
    timelimit       # max time GUROBI is allowed to operate (sec)
    aggregate       # determine aggregator setting (set to 1 to aggregate)

    OptParams() = new()
end

type CostFunction # 1x1 structure containing...
    delThat         # penalty for temperature overload [$/C^2]
    Elvl            # (terminal) cost of deviating from slow-timescale (interpolated) energy storage levels [$/pu^2]
    fG              # (terminal) penalizing generator levels from interpolated references [$/pu^2]
    dfG             # generator ramping cost factor that adds to quadratic gen cost term [$/pu^2]
    windSpill       # penalty for spilling wind from wind farms [$/pu^2]
    loadShed        # penalty for shedding load through demand response [$/pu^2]
    fQ              # storage utility cost [$/pu^2] (i.e. keeping storage level constant is NOT penalized)
    PSTshift        # phase-shift from PST wrt Level 1/2 reference [$/rad^2]

    CostFunction() = new()
end

type Optimize # 1x1 structure
    QPsolver        # optimization solver 2 == GUROBI (1 == CPLEX, currently unavailable)
    params::OptParams
    costFunction::CostFunction

    Optimize() = new()
end

type Time # 1x1 structure
    day     # scalar current day of year
    hour    # scalar current hour of day
    minute  # scalar current minute of hour

    Time() = new()
end

type System
    buses::Array{Bus}
    lines::Array{Line}
    gens::Array{Gen}
    loads::Array{Load}
    wind::Array{Wind}
    storage::Array{Storage}
    params::Params
    control::Control
    optimize::Optimize
    time::Time

    System() = new()
end

function mat2sys(mpc)
    """ Convert from caseformat to System structure.
    """
    baseMVA = mpc["baseMVA"]

    # System.bus
    bus = mpc["bus"]
    BUS = Array(Bus,0)

    for i in 1:size(bus,1)
        push!(BUS,Bus(
            bus[i,1],
            bus[i,2],
            bus[i,10],
            bus[i,8],
            bus[i,9]*π/180.,
            bus[i,5],
            bus[i,6]
            ))
    end

    # System.line
    branch = mpc["branch"]
    LINE = Array(Line,0)

    tap = branch[:,9]
    tap[tap.==0] = 1

    for i in 1:size(branch,1)
        push!(LINE,Line(
            branch[i,1],
            branch[i,2],
            branch[i,3],
            branch[i,4],
            branch[i,5],
            branch[i,6],
            NaN,
            0,
            0,
            0,
            0,
            -1,
            branch[i,9]!=0,
            branch[i,9]!=0,
            false,
            tap[i],
            branch[i,10]*pi/180,
            0,
            ConductorModel(),
            ThermalModel(),
            NaN
            ))
    end

    # System.gen
    gen = mpc["gen"]
    Qmin = gen[:,5]
    Qmax = gen[:,4]
    GEN = Array(Gen,0)

    for i in 1:size(gen,1)
        push!(GEN,Gen(
            gen[i,1],
            1,
            gen[i,8],
            gen[i,2]./baseMVA,
            gen[i,3]./baseMVA,
            1/size(gen,1),
            0,
            gen[i,6],
            gen[i,10]/baseMVA,
            gen[i,9]/baseMVA,
            Qmin,#-9999*ones(length(gen[:,1])),
            Qmax,#9999*ones(length(gen[:,1])),
            gen[i,9]/(2*baseMVA),
            -gen[i,9]/(2*baseMVA),
            gen[i,4]/(2*baseMVA),
            gen[i,5]/(2*baseMVA)
            ))
    end

    # System.load
    LOAD = Array(Load,0)
    for i in 1:size(bus,1)
        push!(LOAD,Load(
        bus[i,1],
        bus[i,3]/100*ones(1,24),
        bus[i,4]/100*ones(1,24),
        bus[i,3]/baseMVA,
        bus[i,4]/baseMVA,
        0,
        0
        ))
    end

    # System.params
    PARAMS = Params()
    PARAMS.sBase = baseMVA
    PARAMS.numOfBus = length(bus[:,1])
    PARAMS.numOfLine = length(branch[:,1])
    PARAMS.numOfGen = length(gen[:,1])
    PARAMS.numOfLoad = length(bus[:,1])

    # System
    SYSTEM = System()
    SYSTEM.buses = BUS
    SYSTEM.lines = LINE
    SYSTEM.gens = GEN
    SYSTEM.loads = LOAD
    SYSTEM.params = PARAMS

    return SYSTEM
end

end


text:
using Jive

a = 1

@__REPL__

@info :a a


text:
# cuda21: Kernels for Array->Vector reduction.

using Knet: reduction_ops

function cuda21src(f, j, op, f1, v0; BLK=64, THR=64) # BLK not used, determined by ny
    sprint() do s
        for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
            print(s,
"""
__device__ void _$(F)_21_0(volatile $T *x, int i) {
//for optimizing warps, volatile must be used as register optimization will lead to wrong answers
  $T ai, xi;
  ai=x[i]; xi=x[i+32]; x[i]=$op;
  ai=x[i]; xi=x[i+16]; x[i]=$op;
  ai=x[i]; xi=x[i+ 8]; x[i]=$op;
  ai=x[i]; xi=x[i+ 4]; x[i]=$op;
  ai=x[i]; xi=x[i+ 2]; x[i]=$op;
  ai=x[i]; xi=x[i+ 1]; x[i]=$op;
}

__global__ void _$(F)_21(int nx, $T *x, int sy, int ny, $T *y) {
  // x[i] goes into y[(i/sy)%ny]
  __shared__ $T buffer[$THR];
  int t = threadIdx.x;
  int b = blockIdx.x;
  $T ai, xi;

  // sum the elements assigned to this thread
  ai = $v0;
  if (sy == 1) {
     int istep = $THR*ny;
     for (int i=b+t*ny; i<nx; i+=istep) {
        xi=x[i]; xi=$f1; ai=$op;
     }
  } else {
    int jstep = sy*ny;
    for (int j=0; j<nx; j+=jstep) {
      int i0 = j+b*sy;
      int i1 = i0+sy;
      for (int i=i0+t; i<i1; i+=$THR) {
        xi=x[i]; xi=$f1; ai=$op;
      }
    }
  }
  buffer[t] = ai;
  __syncthreads();

  // help sum the entries in the block
  for(int stride=$THR/2; stride>32; stride>>=1) {
    if(t < stride) {
      ai=buffer[t]; xi=buffer[stride+t]; buffer[t]=$op;
    }
    __syncthreads();   // Q: can this be outside the for loop?
  }

  if(t<32) {
    _$(F)_21_0(buffer,t);  // This reuses warpSum from 20 scalar reduction.
  }
  __syncthreads();

  if(t==0) {  // the first thread in the block writes the block result to y
    y[blockIdx.x]=buffer[0];
  }
}

#ifdef __cplusplus
extern "C" {
#endif
  void $(F)_21(int nx, $T *x, int sy, int ny, $T *y) {
  // x[i] goes into y[(i/sy)%ny]
  //  _$(F)_21<<<$BLK,$THR>>>(nx,x,sy,ny,y);
  _$(F)_21<<<ny,$THR>>>(nx,x,sy,ny,y);
}
#ifdef __cplusplus
}
#endif
""")
        end
    end
end

for a in reduction_ops
    if !isa(a,Tuple); a=(a,); end
    print(cuda21src(a...))
end


text:
__precompile__()

module ProbNumODE

using Reexport
@reexport using DiffEqBase
import DiffEqBase: check_error!

import Base: copy, copy!

using LinearAlgebra
using Measurements
using ForwardDiff
@reexport using StructArrays
using UnPack
using StaticArrays
using RecipesBase
using Distributions
@reexport using GaussianDistributions
copy(P::Gaussian) = Gaussian(copy(P.μ), copy(P.Σ))
copy!(dst::Gaussian, src::Gaussian) = (copy!(dst.μ, src.μ); copy!(dst.Σ, src.Σ); nothing)
using ModelingToolkit
using DiffEqDevTools
using Optim
using PDMats

using UUIDs, ProgressLogging

import Base: copy

stack(x) = copy(reduce(hcat, x)')

include("filtering.jl")
include("steprules.jl")
include("priors.jl")
# include("measurement_model.jl")
include("sigmas.jl")

include("integrator_type.jl")
export EKF0, EKF1, ODEFilter
include("integrator_interface.jl")
include("integrator_utils.jl")
include("perform_step.jl")
include("solve.jl")
include("preconditioning.jl")
include("postprocessing.jl")
include("solution.jl")

include("dev/problems.jl")
export exponential_decay, logistic_equation, brusselator, fitzhugh_nagumo, lotka_volterra, van_der_pol, fitzhugh_nagumo_iip
include("dev/visualization.jl")
export hairer_plot
include("dev/evaluation.jl")
export MyWorkPrecision, MyWorkPrecisionSet, plot_wps

include("utils/progressbar.jl")
include("utils/rhs_derivatives.jl")

end

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
module Space

import Base:occursin, size, ==
export BoxSpace, DiscreteSpace, MultiBinarySpace, MultiDiscreteSpace, gymspace2jlspace, sample

include("abstractspace.jl")
include("boxspace.jl")
include("discretespace.jl")
include("multibinaryspace.jl")
include("multidiscretespace.jl")

# Tuple Support
sample(s::Tuple{Vararg{<:AbstractSpace}}) = map(sample, s)
occursin(a::Tuple, b::Tuple{Vararg{<:AbstractSpace}}) = length(a) == length(b) &&
    all(map((x, y) -> occursin(x, y), a, b))

# Dict Support
sample(s::Dict{String}) = Dict(map((k, v) -> (k, sample(v)), s))
occursin(a::Dict{String}, b::Dict{String}) = length(a) == length(b) &&
    all(p -> haskey(a, p.first) ?
            occursin(a[p.first], p.second) :
            false,
        b)

include("utils.jl")
end


text:
#! /usr/bin/julia

# Rosetta Code, CSV data manipulation

ifn = "csv_data_manipulation_in.dat"
ofn = "csv_data_manipulation_out.dat"

ifile = open(ifn, "r")
(a, h) = readcsv(ifile, Int, header=true)
close(ifile)

a = hcat(a, sum(a, 2))
h = hcat(h, "SUM")
a = vcat(h, a)

ofile = open(ofn, "w")
writecsv(ofile, a)
close(ofile)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
import Base: iterate, length, show

struct EdgeIter{T}
    g::T
end

length(eit::EdgeIter) = ne(eit.g)

show(io::IO, eit::EdgeIter) = write(io, "EdgeIter $(length(eit))")

function iterate(eit::EdgeIter, state=(1,1))
    fadjlist = fadj(eit.g)
    n = nv(eit.g)
    (i,j) = state
    n == 0 && return nothing

    while i <= n
        if j > length(fadjlist[i])
            i += 1
            j = 1
            continue
        end
        e = Edge(i, fadjlist[i][j])
        state = (i,j+1)
        return (e, state)
    end

    return nothing
end


text:
using Turing, Random, Test
using Turing: Sampler, NUTS

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@testset "hmc.jl" begin
    @numerical_testset "constrained bounded" begin
        obs = [0,1,0,1,1,1,1,1,1,1]

        @model constrained_test(obs) = begin
            p ~ Beta(2,2)
            for i = 1:length(obs)
                obs[i] ~ Bernoulli(p)
            end
            p
        end

        chain = sample(
            constrained_test(obs),
            HMC(1000, 1.5, 3)) # using a large step size (1.5)

        check_numerical(chain, [:p], [10/14], eps=0.1)
    end
    @numerical_testset "contrained simplex" begin
        obs12 = [1,2,1,2,2,2,2,2,2,2]

        @model constrained_simplex_test(obs12) = begin
            ps ~ Dirichlet(2, 3)
            for i = 1:length(obs12)
                obs12[i] ~ Categorical(ps)
            end
            return ps
        end

        chain = sample(
            constrained_simplex_test(obs12),
            HMC(1000, 0.75, 2))

        check_numerical(chain, ["ps[1]", "ps[2]"], [5/16, 11/16], eps=0.015)
    end
    @numerical_testset "hmc reverse diff" begin
        alg = HMC(4000, 0.15, 10)
        res = sample(gdemo_default, alg)
        check_gdemo(res, eps=0.1)
    end
    @turing_testset "matrix support" begin
        @model hmcmatrixsup() = begin
            v ~ Wishart(7, [1 0.5; 0.5 1])
            v
        end

        model_f = hmcmatrixsup()
        vs = []
        chain = nothing
        τ = 3000
        for _ in 1:5
            chain = sample(model_f, HMC(τ, 0.1, 3))
            r = reshape(chain[:v].value, τ, 2, 2)
            push!(vs, reshape(mean(r, dims = [1]), 2, 2))
        end

        @test maximum(abs, mean(vs) - (7 * [1 0.5; 0.5 1])) <= 0.5
    end
    @turing_testset "multivariate support" begin
        function sigmoid(t)
            return 1 / (1 + exp.(-t))
        end

        # Define NN flow
        function nn(x, b1, w11, w12, w13, bo, wo)
            h = tanh.([w11 w12 w13]' * x + b1)
            return sigmoid((wo' * h)[1] + bo)
        end

        # Generating training data
        N = 20
        M = round(Int64, N / 4)
        x1s = rand(M) * 5
        x2s = rand(M) * 5
        xt1s = Array([[x1s[i]; x2s[i]] for i = 1:M])
        append!(xt1s, Array([[x1s[i] - 6; x2s[i] - 6] for i = 1:M]))
        xt0s = Array([[x1s[i]; x2s[i] - 6] for i = 1:M])
        append!(xt0s, Array([[x1s[i] - 6; x2s[i]] for i = 1:M]))

        xs = [xt1s; xt0s]
        ts = [ones(M); ones(M); zeros(M); zeros(M)]

        # Define model

        alpha = 0.16            # regularizatin term
        var = sqrt(1.0 / alpha) # variance of the Gaussian prior

        @model bnn(ts) = begin
            b1 ~ MvNormal([0 ;0; 0],
                [var 0 0; 0 var 0; 0 0 var])
            w11 ~ MvNormal([0; 0], [var 0; 0 var])
            w12 ~ MvNormal([0; 0], [var 0; 0 var])
            w13 ~ MvNormal([0; 0], [var 0; 0 var])
            bo ~ Normal(0, var)

            wo ~ MvNormal([0; 0; 0],
                [var 0 0; 0 var 0; 0 0 var])
            for i = rand(1:N, 10)
                y = nn(xs[i], b1, w11, w12, w13, bo, wo)
                ts[i] ~ Bernoulli(y)
            end
            b1, w11, w12, w13, bo, wo
        end

        # Sampling
        chain = sample(bnn(ts), HMC(10, 0.1, 5))
    end
    Random.seed!(123)
    @numerical_testset "hmcda inference" begin
        alg1 = HMCDA(3000, 1000, 0.8, 0.015)
        # alg2 = Gibbs(3000, HMCDA(1, 200, 0.8, 0.35, :m), HMC(1, 0.25, 3, :s))
        alg3 = Gibbs(500,
            PG(30, 10, :s),
            HMCDA(1, 500, 0.8, 0.005, :m))
        # alg3 = Gibbs(2000, HMC(1, 0.25, 3, :m), PG(30, 3, :s))
        # alg3 = PG(50, 2000)

        res1 = sample(gdemo_default, alg1)
        check_gdemo(res1)

        # res2 = sample(gdemo([1.5, 2.0]), alg2)
        #
        # @test mean(res2[:s]) ≈ 49/24 atol=0.2
        # @test mean(res2[:m]) ≈ 7/6 atol=0.2

        res3 = sample(gdemo_default, alg3)
        check_gdemo(res3)
    end
    @turing_testset "hmcda constructor" begin
        alg = HMCDA(1000, 0.8, 0.75)
        println(alg)
        sampler = Sampler(alg)

        alg = HMCDA(1000, 200, 0.8, 0.75)
        println(alg)
        sampler = Sampler(alg)

        alg = HMCDA(1000, 200, 0.8, 0.75, :s)
        println(alg)
        sampler = Sampler(alg)

        @test isa(alg, HMCDA)
        @test isa(sampler, Sampler{<:Turing.Hamiltonian})
    end
    @numerical_testset "nuts inference" begin
        alg = NUTS(6000, 1000, 0.8)
        res = sample(gdemo_default, alg)
        check_gdemo(res)
    end
    @turing_testset "nuts constructor" begin
        alg = NUTS(1000, 200, 0.65)
        sampler = Sampler(alg)

        alg = NUTS(1000, 0.65)
        sampler = Sampler(alg)

        alg = NUTS(1000, 200, 0.65, :m)
        sampler = Sampler(alg)
    end
    @turing_testset "check discard" begin
        alg = NUTS(500, 100, 0.8)

        c1 = sample(gdemo_default, alg, discard_adapt = true)
        c2 = sample(gdemo_default, alg, discard_adapt = false)

        @test size(c1, 1) == 400
        @test size(c2, 1) == 500
    end
end


text:
export CompletedObservable, completed

import Base: ==
import Base: show

"""
    CompletedObservable{D, H}(scheduler::H)

Observable that emits no items to the Actor and just sends a complete notification on subscription.

# Constructor arguments
- `scheduler`: Scheduler-like object

See also: [`Subscribable`](@ref), [`completed`](@ref)
"""
@subscribable struct CompletedObservable{D, H} <: ScheduledSubscribable{D}
    scheduler :: H
end

getscheduler(observable::CompletedObservable) = observable.scheduler

function on_subscribe!(observable::CompletedObservable, actor, scheduler)
    complete!(actor, scheduler)
    return voidTeardown
end

"""
    completed(::Type{T} = Any; scheduler::H = AsapScheduler()) where { T, H <: AbstractScheduler }

Creation operator for the `CompletedObservable` that emits no items to the Actor and immediately sends a complete notification on subscription.

# Arguments
- `T`: type of output data source, optional, `Any` is the default
- `scheduler`: optional, scheduler-like object

# Examples

```jldoctest
using Rocket

source = completed(Int)
subscribe!(source, logger())
;

# output

[LogActor] Completed

```

See also: [`CompletedObservable`](@ref), [`subscribe!`](@ref), [`logger`](@ref)
"""
completed(::Type{T} = Any; scheduler::H = AsapScheduler()) where { T, H <: AbstractScheduler } = CompletedObservable{T, H}(scheduler)

Base.:(==)(::CompletedObservable{T, H}, ::CompletedObservable{T, H}) where { T, H } = true
Base.:(==)(::CompletedObservable,       ::CompletedObservable)                      = false

Base.show(io::IO, ::CompletedObservable{T, H}) where { T, H } = print(io, "CompletedObservable($T, $H)")


text:
@info "Running functionality tests..."
@testset "Stochastic Programs: Functionality" begin
    tol = 1e-2
    for (model, _scenarios, res, name) in problems
        sp = instantiate(model, _scenarios, optimizer = GLPK.Optimizer)
        @testset "SP Constructs: $name" begin
            optimize!(sp, cache = true)
            @test termination_status(sp) == MOI.OPTIMAL
            @test isapprox(optimal_decision(sp), res.x̄, rtol = tol)
            for i in 1:num_scenarios(sp)
                @test isapprox(optimal_recourse_decision(sp, i), res.ȳ[i], rtol = tol)
            end
            @test isapprox(objective_value(sp), res.VRP, rtol = tol)
            @test isapprox(EWS(sp), res.EWS, rtol = tol)
            @test isapprox(EVPI(sp), res.EVPI, rtol = tol)
            @test isapprox(VSS(sp), res.VSS, rtol = tol)
            @test isapprox(EV(sp), res.EV, rtol = tol)
            @test isapprox(EEV(sp), res.EEV, rtol = tol)
        end
        @testset "Inequalities: $name" begin
            @test EWS(sp) <= VRP(sp)
            @test VRP(sp) <= EEV(sp)
            @test VSS(sp) >= 0
            @test EVPI(sp) >= 0
            @test VSS(sp) <= EEV(sp) - EV(sp)
            @test EVPI(sp) <= EEV(sp) - EV(sp)
        end
        @testset "Copying: $name" begin
            sp_copy = copy(sp, optimizer = GLPK.Optimizer)
            add_scenarios!(sp_copy, scenarios(sp))
            @test num_scenarios(sp_copy) == num_scenarios(sp)
            generate!(sp_copy)
            @test num_subproblems(sp_copy) == num_subproblems(sp)
            optimize!(sp)
            optimize!(sp_copy)
            @test termination_status(sp_copy) == MOI.OPTIMAL
            @test isapprox(optimal_decision(sp_copy), optimal_decision(sp), rtol = tol)
            for i in 1:num_scenarios(sp)
                @test isapprox(optimal_recourse_decision(sp_copy, i), optimal_recourse_decision(sp, i), rtol = tol)
            end
            @test isapprox(objective_value(sp_copy), objective_value(sp), rtol = tol)
            @test isapprox(EWS(sp_copy), EWS(sp), rtol = tol)
            @test isapprox(EVPI(sp_copy), EVPI(sp), rtol = tol)
            @test isapprox(VSS(sp_copy), VSS(sp), rtol = tol)
            @test isapprox(EV(sp_copy), EV(sp), rtol = tol)
            @test isapprox(EEV(sp_copy), EEV(sp), rtol = tol)
        end
    end
    @testset "Sampling" begin
        sampled_sp = instantiate(simple, sampler, 100, optimizer = GLPK.Optimizer)
        generate!(sampled_sp)
        @test num_scenarios(sampled_sp) == 100
        @test isapprox(stage_probability(sampled_sp), 1.0)
        StochasticPrograms.sample!(sampled_sp, sampler, 100)
        generate!(sampled_sp)
        @test num_scenarios(sampled_sp) == 200
        @test isapprox(stage_probability(sampled_sp), 1.0)
    end
    @testset "Instant" begin
        optimize!(simple_sp)
        @test termination_status(simple_sp) == MOI.OPTIMAL
        @test isapprox(optimal_decision(simple_sp), simple_res.x̄, rtol = tol)
        for i in 1:num_scenarios(simple_sp)
            @test isapprox(optimal_recourse_decision(simple_sp, i), simple_res.ȳ[i], rtol = tol)
        end
        @test isapprox(objective_value(simple_sp), simple_res.VRP, rtol = tol)
        @test isapprox(EWS(simple_sp), simple_res.EWS, rtol = tol)
        @test isapprox(EVPI(simple_sp), simple_res.EVPI, rtol = tol)
        @test isapprox(VSS(simple_sp), simple_res.VSS, rtol = tol)
        @test isapprox(EV(simple_sp), simple_res.EV, rtol = tol)
        @test isapprox(EEV(simple_sp), simple_res.EEV, rtol = tol)
    end
    @testset "SMPS" begin
        simple_smps = read("io/smps/simple.smps", StochasticProgram, optimizer = GLPK.Optimizer)
        optimize!(simple_smps)
        @test termination_status(simple_smps) == MOI.OPTIMAL
        @test isapprox(optimal_decision(simple_smps), simple_res.x̄, rtol = tol)
        for i in 1:num_scenarios(simple_smps)
            @test isapprox(optimal_recourse_decision(simple_smps, i), simple_res.ȳ[i], rtol = tol)
        end
        @test isapprox(objective_value(simple_smps), simple_res.VRP, rtol = tol)
        @test isapprox(EWS(simple_smps), simple_res.EWS, rtol = tol)
        @test isapprox(EVPI(simple_smps), simple_res.EVPI, rtol = tol)
        @test isapprox(VSS(simple_smps), simple_res.VSS, rtol = tol)
        @test isapprox(EV(simple_smps), simple_res.EV, rtol = tol)
        @test isapprox(EEV(simple_smps), simple_res.EEV, rtol = tol)
    end
end


text:
module NewOnly

import MLJModelInterface

f(x) = x + x

end # module


text:
using Test
using CubedSphere

B_Rancic_correct = [
    0.00000000000000,
    0.67698819751739,
    0.11847293456554,
    0.05317178134668,
    0.02965810434052,
    0.01912447304028,
    0.01342565621117,
    0.00998873323180,
    0.00774868996406,
    0.00620346979888,
    0.00509010874883,
    0.00425981184328,
    0.00362308956077,
    0.00312341468940,
    0.00272360948942,
    0.00239838086555,
    0.00213001905118,
    0.00190581316131,
    0.00171644156404,
    0.00155493768255,
    0.00141600715207,
    0.00129556597754,
    0.00119042140226,
    0.00109804711790,
    0.00101642216628,
    0.00094391366522,
    0.00087919021224,
    0.00082115710311,
    0.00076890728775,
    0.00072168382969,
    0.00067885087750
]

@testset "CubedSphere" begin
    @testset "Rančić et al. (1996) Taylor coefficients" begin
        for k in 1:length(B_Rancic_correct)
            @test CubedSphere.B_Rancic[k] ≈ B_Rancic_correct[k]
        end
    end
    @testset "Rančić et al. (1996) inverse mapping" begin
        xv=collect(0:0.1:1)
        yv=collect(0:0.1:1)
        errmax=0;
        tol = 1.e-11
        for x in xv, y in yv
          err= (x,y) .- conformal_cubed_sphere_inverse_mapping(conformal_cubed_sphere_mapping(x,y)...)
          errmax=max(errmax,sum(abs.(err)))
        end
        @test errmax < tol
    end

end


text:
arg_dict = Dict{Any,Any}(
"run_type" => 5,
"var2" => "a",
"var4" => [1,2,3],
"Flux_name" => "IRFlux",
"BC1" => [0],
"solve" => false,
"delta_t" => 0.005,
"BC1_name" => "FreeStreamBC",
"physics" => "Euler",
"jac_type" => 1,
"volume_integral_type" => 2,
"face_integral_type" => 2,
"Volume_flux_name" => "IRFlux",
"FaceElementIntegral_name" => "ESLFFaceIntegral",
"use_DG" => true,
"operator_type" => "SBPOmega",
"operator_type2" => "SBPOmega", # make both operators the same for testing
"res_abstol" => 1.0e-12,
"order" => 2,
"order2" => 2,
"res_reltol" => -1.0,
"var1" => 1,
"itermax" => 20,
"smb_name" => "SRCMESHES/tri2l.smb",
"jac_method" => 2,
"var3" => 7,
"numBC" => 1,
"step_tol" => 1.0e-9,
"IC_name" => "ICFreeStream",
"t_max" => 500.0,
"write_face_vertnums" => true,
)


text:
export DenseConvDims

"""
    DenseConvDims

Concrete subclass of `ConvDims` for a normal, dense, conv2d/conv3d.
"""
struct DenseConvDims{N,K,C_in,C_out,S,P,D,F} <: ConvDims{N,S,P,D,F}
    I::NTuple{N,Int}
end

# Getters for the fields
input_size(c::DenseConvDims) = c.I
kernel_size(c::DenseConvDims{N,K,C_in,C_out}) where {N,K,C_in,C_out} = K
channels_in(c::DenseConvDims{N,K,C_in,C_out}) where {N,K,C_in,C_out} = C_in::Int
channels_out(c::DenseConvDims{N,K,C_in,C_out}) where {N,K,C_in,C_out} = C_out::Int

# Convenience wrapper to create DenseConvDims objects
function DenseConvDims(x_size::NTuple{M}, w_size::NTuple{M};
                       stride=1, padding=0, dilation=1, flipkernel::Bool=false) where M
    # Do common parameter validation
    stride, padding, dilation = check_spdf(x_size, w_size, stride, padding, dilation)

    # Ensure channels are equal
    if x_size[end-1] != w_size[end-1]
        xs = x_size[end-1]
        ws = w_size[end-1]
        throw(DimensionMismatch("Input channels must match! ($xs vs. $ws)"))
    end

    # The type parameters are what
    return DenseConvDims{
        M - 2,
        w_size[1:end-2],
        x_size[end-1],
        w_size[end],
        stride,
        padding,
        dilation,
        flipkernel
    }(
        # Input spatial size
        x_size[1:end-2],
    )
end

# Auto-extract sizes and sub out to big brother above
function DenseConvDims(x::AbstractArray, w::AbstractArray; kwargs...)
    if ndims(x) != ndims(w)
        throw(DimensionMismatch("Rank of x and w must match! ($(ndims(x)) vs. $(ndims(w)))"))
    end
    return DenseConvDims(size(x), size(w); kwargs...)
end

# Useful for constructing a new DenseConvDims that has only a few elements different
# from the original progenitor object that it inherits shapes from.
function DenseConvDims(c::ConvDims; N=spatial_dims(c), I=input_size(c), K=kernel_size(c),
                       C_in=channels_in(c), C_out=channels_out(c), S=stride(c),
                       P=padding(c), D=dilation(c), F=flipkernel(c))
    return DenseConvDims{N, K, C_in, C_out, S, P, D, F}(I)
end

function check_dims(x::NTuple{M}, w::NTuple{M}, y::NTuple{M}, cdims::DenseConvDims) where {M}
    # First, check that channel counts are all correct:
    @assert x[M-1] == channels_in(cdims) DimensionMismatch("Data input channel count ($(x[M-1]) vs. $(channels_in(cdims)))")
    @assert y[M-1] == channels_out(cdims) DimensionMismatch("Data output channel count ($(y[M-1]) vs. $(channels_out(cdims)))")
    @assert w[M-1] == channels_in(cdims) DimensionMismatch("Kernel input channel count ($(w[M-1]) vs. $(channels_in(cdims)))")
    @assert w[M] == channels_out(cdims) DimensionMismatch("Kernel output channel count ($(w[M]) vs. $(channels_out(cdims)))")

    # Next, check that the spatial dimensions match up
    @assert x[1:M-2] == input_size(cdims) DimensionMismatch("Data input spatial size ($(x[1:M-2]) vs. $(input_size(cdims)))")
    @assert y[1:M-2] == output_size(cdims) DimensionMismatch("Data output spatial size ($(y[1:M-2]) vs. $(output_size(cdims)))")
    @assert w[1:M-2] == kernel_size(cdims) DimensionMismatch("Kernel spatial size ($(w[1:M-2]) vs. $(kernel_size(cdims)))")

    # Finally, check that the batch size matches
    @assert x[M] == y[M] DimensionMismatch("Batch size ($(x[M]) vs. $(y[M]))")
end


text:
function Div(Px, Py, bound = "sym", order = 1)
    """
        div - divergence operator

        fd = div(Px,Py, options);
        fd = div(P, options);

          options.bound = 'per' or 'sym'
          options.order = 1 (backward differences)
                        = 2 (centered differences)

          Note that the -div and grad operator are adjoint
          of each other such that
              <grad(f),g>=<f,-div(g)>

          See also: grad.

        Copyright (c) 2007 Gabriel Peyre
    """

    # retrieve number of dimensions
    nbdims = ndims(Px)

    if nbdims >= 3
        if nbdims == 3
            Py = P[:, :, 2]
            Px = P[:, :, 1]
            nbdims = 2
        else
            Pz = P[:, :, :, 3]
            Py = P[:, :, :, 2]
            Px = P[:, :, :, 1]
            nbdims = 3
        end
    end

    if bound == "sym"
        nx = size(Px)[1]
        if order == 1
            fx = Px .- Px[vcat(1, collect(1, nx - 1)), :]
            fx[1, :] = Px[1, :]                        # boundary
            fx[nx, :] = - Px[nx - 1, :]

            if nbdims >= 2
                ny = size(Py)[2]
                fy = Py - Py[:, vcat(1, collect(1, ny - 1))]
                fy[:, 1] = Py[:, 1]                    # boundary
                fy[:, ny] = - Py[:, ny - 1]
            end

            if nbdims >= 3
                nz = size(Pz)[3]
                fz = Pz - Pz[:, :, vcat(1, collect(1 : nz - 1))]
                fz[:, :, 1] = Pz[:, :, 1]                # boundary
                fz[:, :, nz] = - Pz[:, :, nz - 1]
            end
        else
            fx = (Px[vcat(collect(2 : nx), nx), :] - Px[vcat([1], collect(1, nx-1)), :])./2.
            fx[1, :] = + Px[1, :]./2. + Px[1, :]           # boundary
            fx[1, :] = + Px[2, :]./2. - Px[1, :]
            fx[nx, :] = - Px[nx, :] - Px[nx - 1, :]./2.
            fx[nx - 1, :] = + Px[nx, :] - Px[nx - 2, :]./2.

            if nbdims >= 2
                ny = size(Py)[2]
                fy = (Py[:, vcat(collect(2 : ny), ny)] - Py[:, vcat(1, collect(1 : ny-1))])./2.
                fy[:, 1] = + Py[:, 2]./2. + Py[:, 1]       # boundary
                fy[:, 2] = + Py[:, 3]./2. - Py[:,1]
                fy[:, ny] = - Py[:, ny] - Py[:, ny-1]./2.
                fy[:, ny - 1] = + Py[:, ny] - Py[:, ny - 2]./2.
            end

            if nbdims >= 3
                nz = size(Pz)[3]
                fz = (Pz[:, :, vcat(collect(2 : nz), nz)] - Pz[:, :, vcat(1, collect(1 : nz - 1))])./2.
                fz[:, :, 1] = + Pz[:, :, 2]./2. + Pz[:, :, 1] # boundary
                fz[:, :, 2] = + Pz[:, :, 3]./2. - Pz[:, :, 1]
                fz[:, :, ny] = - Pz[:, :, nz] - Pz[:, :, nz - 1]./2.
                fz[:, :, ny - 1] = + Pz[:, :, nz] - Pz[:, :, nz - 2]./2.
            end
        end

    else
        if order == 1
            nx = size(Px)[1]
            fx = Px - Px[vcat(nx, collect(1 : nx-1)), :]

            if nbdims >= 2
                ny = size(Py)[2]
                fy = Py - Py[:, vcat(ny, collect(1 : ny-1))]
            end

            if nbdims>=3
                nz = size(Pz)[3]
                fz = Pz - Pz[:, :, vcat(nz, collect(1 : nz-1))]
            end

        else
            nx = size(Px)[1]
            fx = (Px[vcat(collect(2 : nx), 1), :]) - (Px[vcat(nx, collect(1, nx - 1)), :])

            if nbdims >= 2
                ny = size(Py)[2]
                fy = (Py[:, vcat(collect(2 : ny), 1)]) - (Py[:, vcat(ny, collect(1, ny - 1))])
            end

            if nbdims >= 3
                nz = size(Pz)[3]
                fz = (Pz[:, :, vcat(collect(2, nz), 1)]) - (Pz[:, :, vcat(nz, collect(1 : nz - 1))])
            end
        end
    end

    # gather result
    if nbdims == 3
        fd = fx + fy + fz

    elseif nbdims == 2
        fd = fx + fy

    else
        fd = fx
    end

    return fd
end


text:
using Test
using Distances: Euclidean, CosineDist
using Random
using SparseArrays
using LinearAlgebra
using ManifoldLearning
using SimilaritySearch
using SimSearchManifoldLearning

include("manifoldlearning_test.jl")

using SimSearchManifoldLearning: initialize_embedding, fuzzy_simplicial_set, compute_membership_strengths, smooth_knn_dist_opt_binsearch, spectral_layout, optimize_embedding, combine_fuzzy_sets, fit_ab, SMOOTH_K_TOLERANCE

include("utils_tests.jl")
include("umap_tests.jl")