text:
# distance functions
@test edit_distance("AAAATTTTCCCCGGGG", "AAAATTATCCCCGGGG") == 1
@test edit_distance("AAAATTTTCCCCGGGG", "AAAATTTCCCCGGGG") == 1
@test edit_distance("AAAATTTTCCCCGGGG", "AAAATATCCCCGGGG") == 2
@test hamming_distance("AAAATTTTCCCCGGGG", "AAAATTATCCCCGGGG") == 1
@test hamming_distance("AAAATTTTCCCCGGGG", "AAAATTTCCCCGGGG") == 2

function test_merge()
    r1 = dna("CCACCCAAAGCCGCGGGCGCTGATGATGGGTGAGGAGGGGGCGGCAAGATTTCGGGCGCCCCTGCCCTGAACGCCCTCAGCTGCTGCCGCCGGGGCCGCTCCGGAGGCGGGAGCCGGTAGCCGAGCCGGGCGACCTAGAGAACGAGCG")
    r2 = dna("CACTTCATTCAGCGGACCGACAGAGTGGCCGACGCTGAGCCTGACCCGCTCGTTCTCTAGGTCGCCCGGCTCGGCTACCGGCTCCCGCCTCCGGAGCGGCCCCGGCGGCAGCAGCTGAGGGCGTTCAGGGCAGGGGCGCCCGAAATC")
    merged_truth = dna("CCACCCAAAGCCGCGGGCGCTGATGATGGGTGAGGAGGGGGCGGCAAGATTTCGGGCGCCCCTGCCCTGAACGCCCTCAGCTGCTGCCGCCGGGGCCGCTCCGGAGGCGGGAGCCGGTAGCCGAGCCGGGCGACCTAGAGAACGAGCGGGTCAGGCTCAGCGTCGGCCACTCTGTCGGTCCGCTGAATGAAGTG")
    return try_merge(r1, r2) == merged_truth
end

function test_overlap()
    r1 = dna("CCACCCAAAGCCGCGGGCGCTGATGATGGGTGAGGAGGGGGCGGCAAGATTTCGGGCGCCCCTGCCCTGAACGCCCTCAGCTGCTGCCGCCGGGGCCGCTCCGGAGGCGGGAGCCGGTAGCCGAGCCGGGCGACCTAGAGAACGAGCGCCC")
    r2 = dna("CGCTCGTTCTCTAGGTCGCCCGGCTCGGCTACCGGCTCCCGCCTCCGGAGCGGCCCCGGCGGCAGCAGCTGAGGGCGTTCAGGGCAGGGGCGCCCGAAATCTTGCCGCCCCCTCCTCACCCATCATCAGCGCCCGCGGCTTTGGGTGGTTT")
    offset, overlap_len, distance = overlap(r1, r2)
    return offset == -3 && distance == 0
end

@test test_merge()
@test test_overlap()


text:
module ProxyInterfaces
using Compat
export @forward

"""
    @forward T.x functions...
Define methods for `functions` on type `T`, which call the relevant function
on the field `x`.
# Example
```julia
struct Wrapper
    x
end
@forward Wrapper.x  Base.sqrt                                  # now sqrt(Wrapper(4.0)) == 2.0
@forward Wrapper.x  Base.length, Base.getindex, Base.iterate   # several forwarded functions are put in a tuple
@forward Wrapper.x (Base.length, Base.getindex, Base.iterate)  # equivalent to above
```
"""
macro forward(ex, fs)
  check = ex isa Expr && ex.head == :. && ex.args[1] isa Symbol && ex.args[2] isa QuoteNode
  check || error("Syntax: @forward T.x f, g, h")
  T = esc(ex.args[1])
  field = ex.args[2].value  # no escape needed for . syntax
  fs = (fs isa Expr && fs.head == :tuple) ? map(esc, fs.args) : [esc(fs)]
  :($([:($f(x::$T, args...; kwargs...) = (Base.@_inline_meta; $f(x.$field, args...; kwargs...)))
       for f in fs]...);
    nothing)
end

include("iterator.jl")
include("indexing.jl")
include("array.jl")  # depends on @iterator
include("dict.jl")

end # module

ext: jl


text:
# foveation example

using Caesar
# using Neo4j
# using CloudGraphs

addrdict = nothing
include(joinpath(dirname(@__FILE__),"..","database","blandauthremote.jl"))
# prepare the factor graph with just one node
# (will prompt on stdin for db credentials)
cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)

neoids, syms = foveateQueryToPoint(cloudGraph,
        ["SESSTURT21";"SESSTURT38";"SESSTURT45"],
        point=[-9.0;9.0], fovrad=0.5)

# @async cloudimshow(cloudGraph, neoid=neoids[1])
for neoid in neoids
  cloudimshow(cloudGraph, neoid=neoid)
end

#


text:
# An example of the Burgers code for the
# propagation of Burgers equation in time
# studies the pointwise convergence order of a solution
# to Burgers equation as a function of resolution

using Burgers
using Printf
using PyPlot
include("Utils.jl")
using .Utils

function convergence_ratio(it::I, sc::SlopeCalc) where {T, I<:Int}
	# calculates the pointwise convergence ratio function
	# it     : the index for the time at which to slice u(x, t)
	#          must satisfy 0 <= t0 <= Ttot
	# sc     : the type of slope calculation (linear, minmod, constant)

	Ttot     = 2.0   # total duration
	Nt_b     = 512   # number of time divisions for lowest resolution
	Ncells_b = 128   # number of position divisions for lowest resolution
	Nghost   = 2     # padding on position divisions for boundary conditions

	if !(1 <= it <= Nt_b)
		error("The input time index must be between 1 and $Nt_base")
	end

	# The scaling factors
	scale1 = 8
	scale2 = 2*scale1
	scale3 = 2*scale2

	# Creating the Resolution instances
	res1 = Resolution{Float64, Int64}(Ttot, scale1*Nt_b, scale1*Ncells_b, Nghost)
	res2 = Resolution{Float64, Int64}(Ttot, scale2*Nt_b, scale2*Ncells_b, Nghost)
	res3 = Resolution{Float64, Int64}(Ttot, scale3*Nt_b, scale3*Ncells_b, Nghost)

	# Evaluating the solution to Burgers equation at t = t0
	ufun1 = get_evolved_slice(get_sinewave, res1, scale1*it, sc)
	ufun2 = get_evolved_slice(get_sinewave, res2, scale2*it, sc)
	ufun3 = get_evolved_slice(get_sinewave, res3, scale3*it, sc)

	# Calculating the absolute difference between solutions
	# of successive resolutions
	udiff12 = absolute_difference(ufun1, ufun2, ufun1.xdat)
	udiff23 = absolute_difference(ufun2, ufun3, ufun1.xdat)

	# stabilize the solutions, prevents the denominator from being too small
	eps = 1e-6
	uratio = (eps + udiff12)/(eps + udiff23)

	return uratio
end

function plot_convergence(cfun::FunArr{T, U}, imagename::AbstractString) where {T, U}
	# plots and saves the binary logarithm of the pointwise convergence function
	PyPlot.clf()

	xlabel("x")
	xlim(0, 1)
	xticks(range(0, 1, length = 10 + 1))

	ylabel("C(x)")
	ylim(0, 6)
	yticks(range(0, 6, length = 6 + 1))

	plot(cfun.xdat, log.(2, cfun.udat))
	savefig("images/" * imagename * ".png", dpi = 300)
end

# collect the command line arguments
# calculate and plot the convergence ratio
it, sc, imagename = commandlineparser(true)
convergence       = convergence_ratio(it, sc)
plot_convergence(convergence, imagename)


text:
#=
shepplogan.jl
=#

#using MIRTjim: jim, prompt
using ImagePhantoms # many
using Test: @test, @testset, @test_throws, @inferred

@testset "shepp" begin
    for case in (SheppLogan, SheppLoganToft, SheppLoganEmis, SheppLoganBrainWeb)
        @test ImagePhantoms.shepp_logan_values(case()) isa Vector
        @test ellipse_parameters(case()) isa Matrix
    end
    @test ellipse_parameters(SouthPark()) isa Matrix

    image0 = @NOTinferred shepp_logan(256, SheppLoganEmis())
    @test image0 isa Matrix
    image1 = @NOTinferred shepp_logan(256, SheppLoganEmis(); oversample=1)
    image2 = @NOTinferred shepp_logan(256, SheppLoganEmis(); oversample=3)
    @test image0 == image2
#   jim(jim(image1), jim(image2))

    ob = shepp_logan(SheppLoganEmis())
    x = LinRange(-1,1,201) * 0.5
    y = LinRange(-1,1,200) * 0.5
    image = phantom(x, y, ob)
    @test image isa Matrix
    image = phantom(ob).(x,y')
    @test image isa Matrix
#   jim(x, y, image)

    r = LinRange(-0.5,0.5,101)
    ϕ = deg2rad.(0:180)
    sino = radon(ob).(r,ϕ')
    @test sino isa Matrix
    sino = radon(r, ϕ, ob)
    @test sino isa Matrix
#   jim(r, ϕ, sino; aspect_ratio=:none, yflip=false)

    kx = LinRange(-1,1,100) * 9
    ky = LinRange(-1,1,101) * 9
    kspace = spectrum(ob).(kx, ky')
    @test kspace isa Matrix
    kspace = spectrum(kx, ky, ob)
    @test kspace isa Matrix
#   jim(kx, ky, kspace)

    image = shepp_logan(80, 100, SouthPark(), fovs=(1,1))
#   jim(image)
end


text:
struct Population
  risks::AbstractDataFrame
  distances::Union{Nothing, AbstractArray}
end

Population(risks) = Population(risks, nothing)

function individuals(x::Population)
  return size(x.risks, 1)
end

function Base.show(io::IO, x::Population)
  return print(io, "Population object (n=$(individuals(x)))")
end


text:
struct TRSolver{T} <: NearlyExactTRSP
    abstol::T
    maxiter::Integer
end
function (ms::TRSolver)(∇f, H, Δ, p)
    T = eltype(p)
    x, info = trs(H, ∇f, Δ)
    p .= x[:, 1]

    m = dot(∇f, p) + dot(p, H * p) / 2
    interior = norm(p, 2) ≤ Δ
    return (
        p = p,
        mz = m,
        interior = interior,
        λ = info.λ,
        hard_case = info.hard_case,
        solved = true,
    )
end
"""
    initial_safeguards(B, g, Δ)

Returns a tuple of initial safeguarding values for λ. Newton's method might not
work well without these safeguards when the Hessian is not positive definite.
"""
function initial_safeguards(B, g, Δ)
    # equations are on p. 560 of [MORESORENSEN]
    T = eltype(g)
    λS = maximum(-diag(B))

    # they state on the first page that ||⋅|| is the Euclidean norm
    gnorm = norm(g)
    Bnorm = opnorm(B, 1)
    λL = max(T(0), λS, gnorm / Δ - Bnorm)
    λU = gnorm / Δ + Bnorm
    (L = λL, U = λU, S = λS)
end
function safeguard_λ(λ::T, λsg) where {T}
    # p. 558
    λ = min(max(λ, λsg.L), λsg.U)
    if λ ≤ λsg.S
        λ = max(T(1) / 1000 * λsg.U, sqrt(λsg.L * λsg.U))
    end
    λ
end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
# Load Julia packages (libraries) needed  for the snippets in chapter 0

using StatisticalRethinking
import LogDensityProblems: logdensity, logdensity_and_gradient
import StatisticalRethinking: HMC2, generate_n_samples

# CmdStan uses a tmp directory to store the output of cmdstan

ProjDir = @__DIR__
cd(ProjDir)

# ### snippet 9.3

# Construct the logdensity problem

struct clip_9_3_model{TY <: AbstractVector, TX <: AbstractVector}
    "Observations."
    y::TY
    "Covariate"
    x::TX
end

# Make the type callable with the parameters *as a single argument*.

function (problem:: clip_9_3_model)(θ)
    @unpack y, x, = problem    # extract the data
    @unpack muy, mux = θ     # works on the named tuple too
    ll = 0.0
    ll += loglikelihood(Normal(mux, 1), x)
    ll += loglikelihood(Normal(muy, 1), y)
    ll += logpdf(Normal(0, 1), mux)
    ll += logpdf(Normal(0, 1), muy)
    ll
end

# Instantiate the model with data and inits.

Random.seed!(1234591)

N = 100
x = rand(Normal(0, 1), N)
y = rand(Normal(0, 1), N)

p = clip_9_3_model(y, x)
θ = (muy = 0.0, mux=0.0)
p(θ)

# Write a function to return properly dimensioned transformation.

problem_transformation(p::clip_9_3_model) =
    as((muy = asℝ, mux = asℝ))

# Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(problem_transformation(p), p)
∇P = ADgradient(:ForwardDiff, P);

# Tune and sample.

chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000);

# We use the transformation to obtain the posterior from the chain.

posterior = TransformVariables.transform.(Ref(problem_transformation(p)), get_position.(chain));

# Extract the posterior means,

[mean(first, posterior), mean(last, posterior)]

# Draw 200 samples:

function draw_n_samples(model, grad;
  epsilon = 0.03, # Step size
  L = 11, # No of leapfrog steps
  n_samples = 1000, # No of samples
  q = [-0.1, 0.2]) # Initial position

  samples = zeros(n_samples, 2)
  for i in 1:n_samples
    q, ptraj, qtraj, accept, dH = HMC(model, grad, 0.03, 11, q)
    samples[i, :] = q
  end

  samples

end

samples = draw_n_samples(p, ∇P; n_samples=200);
mean(samples, dims=1)

# End of `09/clip-03.jl`


text:
# This file was generated, do not modify it. # hide
mpg = auto.MPG
mpg = auto[:, 1]
mpg = auto[:, :MPG]
mpg |> mean


text:
import InteratomicPotentials: NonTrainablePotential, energy_and_force

# Integration to use DFTK as an interatomic potential
# This is useful for AIMD simulations

Base.@kwdef struct DFTKPotential <: NonTrainablePotential
    functionals::Vector{Symbol}    = [:gga_x_pbe, :gga_c_pbe]  # default to model_PBE
    model_kwargs::Dict{Symbol,Any} = Dict{Symbol,Any}()
    basis_kwargs::Dict{Symbol,Any} = Dict{Symbol,Any}()
    scf_kwargs::Dict{Symbol,Any}   = Dict{Symbol,Any}()
end
function DFTKPotential(Ecut, kgrid; kwargs...)
    p = DFTKPotential(; kwargs...)
    p.basis_kwargs[:Ecut]  = Ecut
    p.basis_kwargs[:kgrid] = kgrid
    p
end

function energy_and_force(system::AbstractSystem, potential::DFTKPotential)
    model = model_DFT(system, potential.functionals; potential.model_kwargs...)
    basis = PlaneWaveBasis(model; potential.basis_kwargs...)

    scfres = self_consistent_field(basis; potential.scf_kwargs...)
    # cache ψ and ρ as starting point for next calculation
    potential.scf_kwargs[:ψ] = scfres.ψ
    potential.scf_kwargs[:ρ] = scfres.ρ

    (; e=scfres.energies.total * u"hartree", f=compute_forces_cart(scfres) * u"hartree/bohr")
end


text:
@testset "79.word-search.jl" begin
    @test is_word_exist(hcat([['A','B','C','E'],['S','F','C','S'],['A','D','E','E']]...), "ABCCED")
    @test is_word_exist(hcat([['A','B','C','E'],['S','F','C','S'],['A','D','E','E']]...), "SEE")

    @test !is_word_exist(hcat([['A','B','C','E'],['S','F','C','S'],['A','D','E','E']]...), "ABCB")
end


text:
# Test evaluation performance of learned SPNs

@testset "Evaluation of large tree-shaped SPNs" begin

    using DataFrames, CSV

    SPN = SumProductNetwork(normpath("$(@__DIR__)/../assets/nltcs.spn"))

    testset = DataFrame!(CSV.File(normpath("$(@__DIR__)/../assets/nltcs.test.csv"),
                             header=collect(1:16) # columns names
                             ));
    # println("Testset has $(size(testset,1)) instances, $(size(testset,2)) columns.")
    test = convert(Matrix, testset) .+ 1
    @test logpdf(SPN,test) ≈ -19582.020235794218 #-21281.999990461
    @test plogpdf(SPN,test) ≈ -19582.020235794218 #-21281.999990461
    # @btime logpdf($SPN,$test)

    # SPN = SumProductNetwork(normpath("$(@__DIR__)/../assets/nltcs.learnspn.em.spn"))
    # @test logpdf(SPN,test) ≈ -22284.141587948394
    # @btime logpdf($SPN,$test)
end


text:
# POSCAR_volumes.jl
# Demonstration / scratchpad for using JuliaPhonons module

push!(LOAD_PATH,"../src") # Temporary versions of modules in PWD
using JuliaPhonons

# OK, a pretty weird thing to do!
# But the POSCAR reader, part of JuliaPhonons, can also do basic stats on the unit cell.

for filename in ARGS
    P=read_POSCAR(open(filename))
    @printf("%s Vol: %f A^3 CubicVector: %f A\n",filename,P.volume,P.volume^(1/3))
end

max_stars_repo_name: jarvist/Julia-Phonons

max_issues_repo_name: lucydot/Julia-Phonons

max_forks_repo_name: lucydot/Julia-Phonons


text:
@doc raw"""
First and second Poincaré invariant for a guiding center particle in an
axisymmetric magnetic field of the form ``B(x,y,z) = B_0 (1 + x^2 + y^2) e_z``.

The loop for the first invariant is initialized by
```math
q (\tau) = \begin{pmatrix}
r_x \cos (2\pi \tau) \\
r_y \sin (2\pi \tau) \\
z_0 + z_1 \sin (2\pi \tau) \\
u_0 + u_1 \cos (2\pi \tau) \\
\end{pmatrix}
```
with parameters
```math
B_0 = 1, \quad
r_x = 0.5, \quad
r_y = 0.3, \quad
z_0 = 0.0, \quad
z_1 = 0.1, \quad
u_0 = 0.5, \quad
u_1 = 0.05, \quad
\mu = 0.01 .
```

The surface for the second invariant is initialized by
```math
q (\tau) = \begin{pmatrix}
r_0 (\sigma - 0.5) \\
r_0 (\tau   - 0.5) \\
z_0 + z_1 \cos (2\pi \sigma) \cos (2\pi \tau) \\
u_0 + u_1 \sin (2\pi \sigma) \sin (2\pi \tau) \\
\end{pmatrix}
```
with parameters
```math
B_0 = 1, \quad
r_0 = 0.5, \quad
z_0 = 0.0, \quad
z_1 = 0.1, \quad
u_0 = 0.5, \quad
u_1 = 0.01, \quad
\mu = 0.01 .
```
"""
module SymmetricField

    import ElectromagneticFields.SymmetricQuadratic

    SymmetricQuadratic.@code() # inject magnetic field code

    μ_loop() = 1E-2
    μ_surface() = 1E-2

    function f_loop(s)
        rx = 0.5
        ry = 0.3
        z0 = 0.0
        z1 = 0.1
        u0 = 5E-1
        u1 = 5E-2

        xs = rx*cos(2π*s)
        ys = ry*sin(2π*s)
        zs = z0 + z1 * sin(2π*s)
        us = u0 + u1 * cos(2π*s)

        qs = [xs, ys, zs, us]

        return qs
    end

    function f_surface(s,t)
        r0 = 0.5
        z0 = 0.0
        z1 = 0.1
        u0 = 5E-1
        u1 = 1E-2

        x  = r0*(s-0.5)
        y  = r0*(t-0.5)
        z  = z0 + z1 * cos(2π*s) * cos(2π*t)
        u  = u0 + u1 * sin(2π*s) * sin(2π*t)

        q  = [x, y, z, u]

        return q
    end

    include("guiding_center_4d_common.jl")
    include("guiding_center_4d_equations.jl")
    include("guiding_center_4d_loop.jl")
    include("guiding_center_4d_surface.jl")
    include("guiding_center_4d_diagnostics.jl")

end


text:
abstract type AbstractPermutation{Ti} end

struct DomainPermutation{Ti} <: AbstractPermutation{Ti}
    prm::Vector{Ti}
end

struct MapPermutation{Ti} <: AbstractPermutation{Ti}
    asg::Vector{Ti}
end

function Base.:(==)(σ::MapPermutation, σ′::MapPermutation)
    return σ.asg == σ′.asg
end

function Base.:(==)(σ::DomainPermutation, σ′::DomainPermutation)
    return σ.prm == σ′.prm
end

function Base.convert(T::Type{<:DomainPermutation}, σ::MapPermutation{Ti}) where {Ti}
    return T(invperm(σ.asg))
end

function Base.convert(T::Type{<:MapPermutation}, σ::DomainPermutation{Ti}) where {Ti}
    return T(invperm(σ.prm))
end

function perm(σ::AbstractPermutation)
    return convert(DomainPermutation, σ).prm
end

struct IdentityPermuter end

function permute_stripe(A::AbstractMatrix, ::IdentityPermuter; kwargs...)
    (m, n) = size(A)
    return DomainPermutation(collect(1:n))
end

function permute_plaid(A::AbstractMatrix, ::IdentityPermuter; kwargs...)
    (m, n) = size(A)
    return (DomainPermutation(collect(1:m)), DomainPermutation(collect(1:n)))
end

function Base.getindex(σ::AbstractPermutation{Ti}, Π::SplitPartition) where {Ti}
    @inbounds return DomainPartition(Π.K, (1:(Π.spl[end] - 1))[perm(σ)], Π.spl)
end

function Base.getindex(σ::AbstractPermutation, Π::DomainPartition)
    @inbounds return DomainPartition(Π.K, Π.prm[perm(σ)], Π.spl)
end

function Base.getindex(σ::AbstractPermutation, Π::MapPartition)
    @inbounds return MapPartition(Π.K, Π.asg[convert(MapPermutation, σ).asg])
end


text:
using Literate

for tut1 in readdir(".")
    if occursin(r".*-tutorial.jl", tut1)
        # println("\nExample $ex1 in $(pwd())\n")
        Literate.markdown(tut1, "."; documenter=false);
    end
end


text:
# This file is part of Kpax3. License is MIT.

function logdPriorRow(p::Vector{Int},
                      ep::EwensPitmanPAUT)
  n = length(p)
  k = 0

  m = zeros(Int, n)

  for a in 1:n
    m[p[a]] += 1

    if m[p[a]] == 1
      k += 1
    end
  end

  logp = (k - 1) * log(ep.α) +
         SpecialFunctions.loggamma(ep.θ / ep.α + k) -
         SpecialFunctions.loggamma(ep.θ / ep.α + 1) +
         SpecialFunctions.loggamma(ep.θ + 1) -
         SpecialFunctions.loggamma(ep.θ + n) -
         k * SpecialFunctions.loggamma(1 - ep.α)

  for a in 1:n
    if m[a] > 0
      logp += SpecialFunctions.loggamma(m[a] - ep.α)
    end
  end

  logp
end

function logdPriorRow(n::Int,
                      k::Int,
                      m::Vector{Int},
                      ep::EwensPitmanPAUT)
  logp = (k - 1) * log(ep.α) +
         SpecialFunctions.loggamma(ep.θ / ep.α + k) -
         SpecialFunctions.loggamma(ep.θ / ep.α + 1) +
         SpecialFunctions.loggamma(ep.θ + 1) -
         SpecialFunctions.loggamma(ep.θ + n) -
         k * SpecialFunctions.loggamma(1 - ep.α)

  for a in 1:length(m)
    if m[a] > 0
      logp += SpecialFunctions.loggamma(m[a] - ep.α)
    end
  end

  logp
end

function logdPriorRow(p::Vector{Int},
                      ep::EwensPitmanPAZT)
  n = length(p)
  k = 0

  m = zeros(Int, n)

  for a in 1:n
    m[p[a]] += 1

    if m[p[a]] == 1
      k += 1
    end
  end

  logp = (k - 1) * log(ep.α) +
         SpecialFunctions.loggamma(k) -
         SpecialFunctions.loggamma(n) -
         k * SpecialFunctions.loggamma(1 - ep.α)

  for a in 1:n
    if m[a] > 0
      logp += SpecialFunctions.loggamma(m[a] - ep.α)
    end
  end

  logp
end

function logdPriorRow(n::Int,
                      k::Int,
                      m::Vector{Int},
                      ep::EwensPitmanPAZT)
  logp = (k - 1) * log(ep.α) +
         SpecialFunctions.loggamma(k) -
         SpecialFunctions.loggamma(n) -
         k * SpecialFunctions.loggamma(1 - ep.α)

  for a in 1:length(m)
    if m[a] > 0
      logp += SpecialFunctions.loggamma(m[a] - ep.α)
    end
  end

  logp
end

function logdPriorRow(p::Vector{Int},
                      ep::EwensPitmanZAPT)
  n = length(p)
  k = 0

  m = zeros(Int, n)

  for a in 1:n
    m[p[a]] += 1

    if m[p[a]] == 1
      k += 1
    end
  end

  logp = k * log(ep.θ) +
         SpecialFunctions.loggamma(ep.θ) -
         SpecialFunctions.loggamma(ep.θ + n)

  for a in 1:n
    if m[a] > 0
      logp += SpecialFunctions.loggamma(m[a])
    end
  end

  logp
end

function logdPriorRow(n::Int,
                      k::Int,
                      m::Vector{Int},
                      ep::EwensPitmanZAPT)
  logp = k * log(ep.θ) +
         SpecialFunctions.loggamma(ep.θ) -
         SpecialFunctions.loggamma(ep.θ + n)

  for a in 1:length(m)
    if m[a] > 0
      logp += SpecialFunctions.loggamma(m[a])
    end
  end

  logp
end

function logdPriorRow(p::Vector{Int},
                      ep::EwensPitmanNAPT)
  n = length(p)
  k = 0

  m = zeros(Int, n)

  for a in 1:n
    m[p[a]] += 1

    if m[p[a]] == 1
      k += 1
    end
  end

  log(prod((1:(k - 1)) .- ep.L) * ep.α^(k - 1) *
      exp(sum(SpecialFunctions.loggamma.(m[m .> 0] .- ep.α)) -
          k * SpecialFunctions.loggamma(1 - ep.α)) /
      prod((1:(n - 1)) .- ep.α * ep.L))
end

function logdPriorRow(n::Int,
                      k::Int,
                      m::Vector{Int},
                      ep::EwensPitmanNAPT)
  log(prod((1:(k - 1)) .- ep.L) * ep.α^(k - 1) *
      exp(sum(SpecialFunctions.loggamma.(m[m .> 0] .- ep.α)) -
          k * SpecialFunctions.loggamma(1 - ep.α)) /
      prod((1:(n - 1)) .- ep.α * ep.L))
end

"""
# Density of the Ewens-Pitman distribution

## Description

Probability of a partition according to the Ewens-Pitman distribution.

## Usage

dPriorRow(ep, p)
dPriorRow(ep, n, k, m)

## Arguments

* `ep` Object of (super)type EwensPitman
* `p` Vector of integers representing a partition
* `n` Set size (Integer)
* `k` Number of blocks (Integer)
* `m` Vector of integers representing block sizes

## Details

## Examples

"""
function dPriorRow(p::Vector{Int},
                   ep::EwensPitman)
  exp(logdPriorRow(p, ep))
end

function dPriorRow(n::Int,
                   k::Int,
                   m::Vector{Int},
                   ep::EwensPitman)
  exp(logdPriorRow(n, k, m, ep))
end


text:
using ..SearchProblem: AbstractSearchProblem

abstract type AbstractTreeSearchAlgorithm{S,A} end

init!(search_algo::AbstractTreeSearchAlgorithm, search_prob::AbstractSearchProblem) = error("not implemented")
stop_condition(search_algo::AbstractTreeSearchAlgorithm, search_prob::AbstractSearchProblem) = error("not implemented")
node_expansion_policy(search_algo::AbstractTreeSearchAlgorithm, search_prob::AbstractSearchProblem) = error("not implemented")
before_node_expansion!(search_algo::AbstractTreeSearchAlgorithm, node::TreeSearchNode,search_prob::AbstractSearchProblem) = error("not implemented")
on_node_generation!(search_algo::AbstractTreeSearchAlgorithm, child_node::TreeSearchNode,search_prob::AbstractSearchProblem) = error("not implemented")
on_node_expansion_finish!(search_algo::AbstractTreeSearchAlgorithm, node::TreeSearchNode,search_prob::AbstractSearchProblem) = error("not implemented")
on_stop!(search_algo::AbstractTreeSearchAlgorithm, search_prob::AbstractSearchProblem) = error("not implemented")

function graph_search!(search_algo::AbstractTreeSearchAlgorithm{S,A}, search_prob::AbstractSearchProblem{S,A}) where {S,A}
    init!(search_algo, search_prob)
    while !stop_condition(search_algo, search_prob)
        node::TreeSearchNode{S,A} = node_expansion_policy(search_algo, search_prob)
        before_node_expansion!(search_algo, node, search_prob)
        for child_node::TreeSearchNode{S,A} in get_children_nodes(search_prob, node)
            on_node_generation!(search_algo, child_node, search_prob)
        end
        on_node_expansion_finish!(search_algo, node, search_prob)
    end
    on_stop!(search_algo, search_prob)
    return search_algo
end


text:
"""
    Haversine(radius=6_371_000)

The haversine distance between two locations on a sphere of given `radius`, whose
default value is 6,371,000, i.e., the Earth's (volumetric) mean radius in meters; cf.
[NASA's Earth Fact Sheet](https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html).

Locations are described with longitude and latitude in degrees.
The computed distance has the unit of the radius.
"""
struct Haversine{T<:Number} <: Metric
    radius::T
end
Haversine{T}() where {T<:Number} = Haversine(T(6_371_000))
Haversine() = Haversine{Int}()

function (dist::Haversine)(x, y)
    length(x) == length(y) == 2 || haversine_error(dist)

    λ₁, φ₁ = x
    λ₂, φ₂ = y

    Δλ = λ₂ - λ₁  # longitudes
    Δφ = φ₂ - φ₁  # latitudes

    # haversine formula
    a = sind(Δφ/2)^2 + cosd(φ₁)*cosd(φ₂)*sind(Δλ/2)^2

    # distance on the sphere
    2 * (dist.radius * asin( min(√a, one(a)) )) # take care of floating point errors
end

haversine(x, y, radius::Number=6_371_000) = Haversine(radius)(x, y)

@noinline haversine_error(dist) = throw(ArgumentError("expected both inputs to have length 2 in $dist distance"))

result_type(::Haversine{T1}, ::Type{T2}, ::Type{T3}) where {T1<:Number,T2<:Number,T3<:Number} =
    float(promote_type(T1, T2, T3))

"""
    SphericalAngle()

The spherical angle distance between two locations on a sphere.

Locations are described with two angles, longitude and latitude, in radians.
The distance is computed with the haversine formula and also has units of radians.
"""
struct SphericalAngle <: Metric end

function (dist::SphericalAngle)(x, y)
    length(x) == length(y) == 2 || haversine_error(dist)

    λ₁, φ₁ = x
    λ₂, φ₂ = y

    Δλ = λ₂ - λ₁  # longitudes
    Δφ = φ₂ - φ₁  # latitudes

    # haversine formula
    a = sin(Δφ/2)^2 + cos(φ₁)*cos(φ₂)*sin(Δλ/2)^2

    # distance on the sphere
    2 * asin( min(√a, one(a)) ) # take care of floating point errors
end

const spherical_angle = SphericalAngle()

result_type(::SphericalAngle, ::Type{T1}, ::Type{T2}) where {T1<:Number,T2<:Number} = float(promote_type(T1, T2))


text:
module NMF_Exact

include("../Misc.jl")
using .Misc

using Distributions, SpecialFunctions, Einsum, Combinatorics
import Base.Iterators: product

export generate, log_marginal

function generate(I::Ƶ, J::Ƶ, K::Ƶ; μ::ℜ=1.0, γ::ℜ=0.1) where {ℜ<:Real, Ƶ<:Int}
    a::ℜ, b::ℜ = I*J*γ, γ/μ
    α_K, α_JK, α_IK = fill(a/K,K), fill(a/(K*J),J,K), fill(a/(K*I),I,K)

    λ = rand(Gamma(a,1.0/b))
    θ_K = rand(Dirichlet(α_K))
    θ_JK = reshape([θ_jk for k=1:K for θ_jk=rand(Dirichlet(α_JK[:,k]))],J,K)
    θ_IK = reshape([θ_ik for k=1:K for θ_ik=rand(Dirichlet(α_IK[:,k]))],I,K)

    θ_IJ = Array{ℜ}(undef,I,J)
    @einsum θ_IJ[i,j] = θ_K[k]*θ_JK[j,k]*θ_IK[i,k]
    Λ_IJ = λ.*θ_IJ

    X::Array{ℜ} = map.(λ_ij -> rand(Poisson(λ_ij)), Λ_IJ)
    return X, λ, (θ_K, θ_IK, θ_JK)
end

function generate(I::Ƶ, J::Ƶ, K::Ƶ, T::Ƶ; γ::ℜ=0.1) where {ℜ<:Real, Ƶ<:Int}
    a::ℜ = I*J*γ
    α_K, α_JK, α_IK = fill(a/K,K), fill(a/(K*J),J,K), fill(a/(K*I),I,K)

    θ_K = rand(Dirichlet(α_K))
    θ_JK = reshape([θ_jk for k=1:K for θ_jk=rand(Dirichlet(α_JK[:,k]))],J,K)
    θ_IK = reshape([θ_ik for k=1:K for θ_ik=rand(Dirichlet(α_IK[:,k]))],I,K)

    θ_IJ = Array{ℜ}(undef,I,J)
    @einsum θ_IJ[i,j] = θ_K[k]*θ_JK[j,k]*θ_IK[i,k]
    θ_flat = reshape(θ_IJ, I*J)

    X::Array{ℜ} = reshape(rand(Multinomial(T,θ_flat)),I,J)
    return X, (θ_K, θ_IK, θ_JK)
end

function log_marginal(X::Array{ℜ,2}, K::Ƶ; μ::ℜ=nanmean(X), γ::ℜ=0.1) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ, J::Ƶ = size(X)
    X_flat = reshape(Int.(X),I*J)

    a::ℜ, b::ℜ = I*J*γ, γ/μ
    α_K, α_JK, α_IK = fill(a/K,K), fill(a/(K*J),J,K), fill(a/(K*I),I,K)

    s, S_K, S_JK, S_IK, S⁺ = zeros(ℜ,K), zeros(ℜ,K), zeros(ℜ,J,K), zeros(ℜ,I,K), sum(X)
    log_PX, log_PS = -Inf, -Inf
    i::Ƶ, j::Ƶ = 0, 0

    log_C = a*log(b) - (a + S⁺)*log(b + 1.0) #+ lgamma(a + S⁺) - lgamma(a)
    log_C += sum(lgamma.(α_K)) - sum(lgamma.(α_IK)) - sum(lgamma.(α_JK))

    if K == 1
        S_K .= S⁺
        S_JK[:,1] = sum(X,dims=1)
        S_IK[:,1] = sum(X,dims=2)
        log_PS = log_C
        log_PS += sum(lgamma.(α_IK .+ S_IK)) + sum(lgamma.(α_JK .+ S_JK)) - sum(lgamma.(α_K .+ S_K))
        log_PS -= sum(lgamma.(X .+ 1.0))
        return log_PS
    end

    for S_div = product(map(X_ij -> combinations(1:(X_ij+K-1),K-1), X_flat)...)
        S_K .= α_K
        S_JK .= α_JK
        S_IK .= α_IK
        log_PS = 0.0
        for (ij, part) = enumerate(S_div)
            i, j = mod(ij-1,I)+1, div(ij-1,I)+1
            s[1] = part[1]-1
            for k=2:K-1
                s[k] = part[k] - part[k-1] - 1
            end
            s[K] = X[i,j] + K - 1 - part[K-1]

            S_K .+= s
            S_JK[j,:] .+= s
            S_IK[i,:] .+= s
            log_PS -= sum(lgamma.(s .+ 1.0))
        end
        log_PS += sum(lgamma.(S_IK)) + sum(lgamma.(S_JK)) - sum(lgamma.(S_K)) + log_C
        log_PX = logsumexp([log_PX,log_PS])
    end

    return log_PX
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
# include standard julia modules, and outside packages -
using JSON
using Dates
using DelimitedFiles
using MAT
using Logging
using LightXML
using TOML

# constants -
const path_to_package = dirname(pathof(@__MODULE__))

# include my files -
include("Types.jl")
include("Extensions.jl")
include("Report.jl")
include("Parser.jl")
include("Checks.jl")
include("Utility.jl")
include("Sequence.jl")
include("Metabolism.jl")
include("Cobra.jl")

# language specific generation -
include("MakeJulia.jl")
include("MakeMatlab.jl")
include("MakeOctave.jl")
include("MakePython.jl")
include("MakeVFF.jl")
include("MakeSBML.jl")
include("MakeDefaultProject.jl")

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule gsasl_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("gsasl")
JLLWrappers.@generate_main_file("gsasl", UUID("2e8251e1-4612-5cea-b4a3-ba3e3b422b6d"))
end  # module gsasl_jll


text:
@testset "matrix/cube" begin

    @testset "flatten" begin
        @test flatten(ones(3, 4, 4)) == ones(3, 16)
        @test flatten(ones(3, 4)) == ones(3, 4)
    end

    @testset "expand" begin
        @test expand(ones(3, 16)) == ones(3, 4, 4)
        @test expand(ones(3, 4, 4)) == ones(3, 4, 4)

        @test_throws ErrorException expand(ones(3, 15))
    end

    # Simple regression
    X = rand(rng, 10, 512, 512)
    @test expand(flatten(X)) == X
end

@testset "derotate" begin
    X = rand(rng, 10, 512, 512)

    # test no-op
    @test derotate(X, zeros(10)) === X

    # Simple regression (but be careful to only test some inner part that won't be clipped)
    X = rand(rng, 10, 512, 512)
    θ = 90 .* rand(rng, 10)
    B = derotate!(derotate(X, θ), -θ)
    # take some inner indices not wiped out due to cropping
    idxs = Colon(), 250:262, 250:262
    @test_broken X[idxs...] ≈ B[idxs...] rtol=1e-3

    # make sure it mutates
    targ = deepcopy(X)
    derotate!(X, θ)
    @test X !== targ

    # trivial
    X = zeros(1, 3, 3)
    X[1, 1, 2] = 1
    # expect the 1 to rotate 90° ccw physical, cw image
    @test derotate(X, [90])[1, 2, 3] ≈ 1 rtol=1e-3

end

@testset "collapse" begin
    X = rand(rng, 10, 512, 512)
    angles = sort!(90 .* rand(rng, 10))

    @test collapse(X) == collapse(X; method = median)
    @test collapse(X; method = mean) == mean(X, dims = 1)[1, :, :]

    @test !(collapse(X, angles) ≈ collapse(derotate(X, angles)))
    @test collapse(X, angles, method=median) ≈ collapse(derotate(X, angles))
    @test any(collapse(X, angles, fill=NaN) .=== NaN)

    @test collapse(ones(10, 512, 512), angles) == mean(derotate!(ones(10, 512, 512), angles), dims=1)[1, :, :]

    X_copy = deepcopy(X)
    @test collapse(X, angles) == collapse!(X_copy, angles)
    @test X_copy != X

end

@testset "shift frame" begin
    X = zeros(3, 3)
    X[2, 2] = 1
    for dx in -1:1, dy in -1:1
        row = 2 + dy
        col = 2 + dx
        @test shift_frame(X, dx, dy)[row, col] == 1
        @test shift_frame(X, (dx, dy)) == shift_frame(X, dx, dy)
    end

    @test shift_frame(X, 1, 1, fill=NaN)[3, 1] === NaN

    # cube
    cube = zeros(10, 3, 3)
    cube[:, 2, 2] .= 1
    for dx in -1:1, dy in -1:1
        row = 2 + dy
        col = 2 + dx
        @test shift_frame(cube, dx, dy)[:, row, col] == ones(10)
        shift_frame(cube, (dx, dy)) == shift_frame(cube, dx, dy)
        shift_frame(cube, fill(dx, 10), fill(dy, 10)) == shift_frame(cube, dx, dy)
        shift_frame(cube, fill((dx, dy), 10)) == shift_frame(cube, dx, dy)
    end
    @test all(shift_frame(cube, 1, 1, fill=NaN)[:, 3, 1] .=== NaN)
end

@testset "crop & cropview" begin
    cube = zeros(10, 101, 101)

    @test crop(cube, 51) == crop(cube, (51, 51)) == cube[:, 26:76, 26:76]
    @test crop(cube, 51, center=(52, 50)) == cube[:, 27:77, 25:75]

    @test size(crop(cube, 51)) == (10, 51, 51)
    @test_logs (:info, "adjusted size to (51, 51) to avoid uneven (odd) cropping") begin
        @test size(crop(cube, 50)) == (10, 51, 51)
    end
    @test size(crop(cube, 50; force=true)) == (10, 50, 50)

    cube_crop_view = cropview(cube, 51)
    @test parent(cube_crop_view) === cube
    @test collect(cube_crop_view) == Array(cube_crop_view) == crop(cube, 51)

    ### frames

    frame = zeros(101, 101)

    @test crop(frame, 51) == crop(frame, (51, 51)) == frame[26:76, 26:76]
    @test crop(frame, 51, center=(52, 50)) == frame[27:77, 25:75]

    @test size(crop(frame, 51)) == (51, 51)
    @test_logs (:info, "adjusted size to (51, 51) to avoid uneven (odd) cropping") begin
        @test size(crop(frame, 50)) == (51, 51)
    end
    @test size(crop(frame, 50; force=true)) == (50, 50)

    frame_crop_view = cropview(frame, 51)
    @test parent(frame_crop_view) === frame
    @test collect(frame_crop_view) == Array(frame_crop_view) == crop(frame, 51)
end


text:
"""

The exported names are
$(EXPORTS)
"""
module FESpaces

using DocStringExtensions
using Test
using FillArrays
using SparseArrays
using LinearAlgebra
using BlockArrays

using Gridap.Inference
using Gridap.Helpers
using Gridap.Arrays
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Fields
using Gridap.Integration
using Gridap.Algebra
using Gridap.Polynomials
using Gridap.TensorValues
using Gridap.CellData

using Gridap.ReferenceFEs: evaluate_dof_array

using Gridap.Arrays: _split
using Gridap.Arrays: Reindexed
using Gridap.Arrays: IdentityVector

import Gridap.Arrays: get_array
import Gridap.Arrays: array_cache
import Gridap.Arrays: getindex!
import Gridap.Arrays: kernel_cache
import Gridap.Arrays: apply_kernel!
import Gridap.Arrays: kernel_return_type
import Gridap.Arrays: kernel_testitem!
import Gridap.Arrays: reindex
import Gridap.Arrays: apply
import Gridap.CellData: get_cell_map
import Gridap.Geometry: get_reffes
import Gridap.Geometry: get_cell_type
import Gridap.CellData: RefStyle
import Gridap.CellData: get_cell_axes
import Gridap.CellData: change_ref_style
import Gridap.Helpers: operate
import Gridap.Geometry: similar_object
import Gridap.Geometry: restrict
import Gridap.Geometry: get_cell_id
import Gridap.Fields: integrate
import Gridap.Fields: evaluate
import Gridap.Fields: gradient
import Gridap.Fields: grad2curl
import Gridap.Fields: evaluate_field_array
import Gridap.ReferenceFEs: evaluate_dof_array
import Gridap.Fields: field_cache
import Gridap.Fields: evaluate_field!
import Gridap.Fields: field_gradient

import Gridap.Algebra: allocate_residual
import Gridap.Algebra: allocate_jacobian
import Gridap.Algebra: residual!
import Gridap.Algebra: jacobian!
import Gridap.Algebra: residual
import Gridap.Algebra: jacobian
import Gridap.Algebra: residual_and_jacobian!
import Gridap.Algebra: residual_and_jacobian
import Gridap.Algebra: zero_initial_guess
import Gridap.Algebra: get_matrix
import Gridap.Algebra: get_vector
import Gridap.Algebra: solve!
import Gridap.Algebra: solve
import Gridap.Algebra: allocate_vector
import Gridap.Algebra: allocate_matrix
import Gridap.Algebra: allocate_matrix_and_vector

export FEFunctionStyle
export is_a_fe_function
export get_free_values
export get_fe_space
export get_cell_values
export test_fe_function

export FESpace
export FEFunction
export EvaluationFunction
export num_free_dofs
export get_cell_basis
export get_cell_dofs
export zero_free_values
export constraint_style
export has_constraints
export get_cell_isconstrained
export get_cell_constraints
export get_cell_axes
export get_cell_axes_with_constraints
export test_fe_space

export Assembler
export AssemblyStrategy
export row_map
export col_map
export row_mask
export col_mask
export DefaultAssemblyStrategy
export get_test
export get_trial
export assemble_matrix!
export assemble_matrix_add!
export assemble_matrix
export assemble_vector!
export assemble_vector_add!
export assemble_vector
export assemble_matrix_and_vector!
export assemble_matrix_and_vector_add!
export assemble_matrix_and_vector
export allocate_vector
export allocate_matrix
export allocate_matrix_and_vector
export test_assembler
export get_matrix_type
export get_vector_type
export count_matrix_nnz_coo
export count_matrix_and_vector_nnz_coo
export fill_matrix_coo_symbolic!
export fill_matrix_and_vector_coo_symbolic!
export fill_matrix_coo_numeric!
export fill_matrix_and_vector_coo_numeric!
export test_sparse_matrix_assembler

export SingleFieldFESpace
export num_dirichlet_dofs
export zero_dirichlet_values
export gather_free_and_dirichlet_values
export gather_free_and_dirichlet_values!
export scatter_free_and_dirichlet_values
export get_dirichlet_values
export gather_dirichlet_values
export gather_dirichlet_values!
export num_dirichlet_tags
export gather_free_values
export gather_free_values!
export get_dirichlet_dof_tag
export compute_free_and_dirichlet_values
export compute_dirichlet_values
export compute_free_values
export compute_dirichlet_values_for_tags
export compute_dirichlet_values_for_tags!
export test_single_field_fe_space
export interpolate
export interpolate!
export interpolate_everywhere
export interpolate_everywhere!
export interpolate_dirichlet
export interpolate_dirichlet!
export get_cell_dof_basis

export SingleFieldFEFunction

export UnconstrainedFESpace
export GradConformingFESpace
export DiscontinuousFESpace

export FECellBasisStyle
export is_a_fe_cell_basis

export TrialFESpace
export TrialFESpace!
export HomogeneousTrialFESpace
export HomogeneousTrialFESpace!
export TestFESpace
export compute_conforming_cell_dofs
export SparseMatrixAssembler
export GenericSparseMatrixAssembler

export FEOperator
export test_fe_operator
export AffineFEOperator
export FEOperatorFromTerms
export get_algebraic_operator

export FESolver
export LinearFESolver
export NonlinearFESolver
export test_fe_solver

export FETerm
export AffineFETerm
export LinearFETerm
export FESource
export AffineFETermFromCellMatVec
export FETermFromCellJacRes
export get_cell_matrix
export get_cell_vector
export get_cell_jacobian
export get_cell_jacobian_and_residual
export get_cell_residual
export collect_cell_matrix
export collect_cell_vector
export collect_cell_matrix_and_vector
export collect_cell_jacobian
export collect_cell_jacobian_and_residual
export collect_cell_residual

export FESpaceWithConstantFixed
export FESpaceWithLastDofRemoved

export ZeroMeanFESpace
export CLagrangianFESpace
export ConformingFESpace
export DirichletFESpace
export FESpaceWithLinearConstraints
export ExtendedFESpace

export autodiff_cell_residual_from_energy
export autodiff_cell_jacobian_from_energy
export autodiff_cell_jacobian_from_residual

export FEEnergy

include("FEFunctions.jl")

include("FECellBases.jl")

include("FESpacesInterfaces.jl")

include("SingleFieldFESpaces.jl")

include("SingleFieldFEFunctions.jl")

include("FESpaceFactories.jl")

include("UnconstrainedFESpaces.jl")

include("ConformingFESpaces.jl")

include("TrialFESpaces.jl")

include("DiscontinuousFESpaces.jl")

include("Assemblers.jl")

include("SparseMatrixAssemblers.jl")

include("FETerms.jl")

include("FEOperators.jl")

include("AffineFEOperators.jl")

include("FEOperatorsFromTerms.jl")

include("FESolvers.jl")

include("FESpacesWithConstantFixed.jl")

include("ZeroMeanFESpaces.jl")

include("CLagrangianFESpaces.jl")

include("DirichletFESpaces.jl")

include("ExtendedFESpaces.jl")

include("FESpacesWithLinearConstraints.jl")

include("FEAutodiff.jl")

include("FETermsWithAutodiff.jl")

end # module


text:
module Knet

"Construct a path relative to Knet root, e.g. Knet.dir(\"examples\") => \"~/.julia/dev/Knet/examples\""
dir(path...)=joinpath(dirname(@__DIR__),path...)

"Default array and element type used by Knet, override by setting Knet.atype() or Knet.array_type[]"
atype() = array_type[]
atype(x) = convert(atype(),x)
const array_type = Ref{Type}(Array{Float32})

include("libknet8/LibKnet8.jl")
include("knetarrays/KnetArrays.jl")
include("cuarrays/CuArrays.jl")
include("autograd_gpu/AutoGrad_gpu.jl")
include("ops20/Ops20.jl")
include("ops20_gpu/Ops20_gpu.jl")
include("ops21/Ops21.jl")
include("ops21_gpu/Ops21_gpu.jl")
include("fileio_gpu/FileIO_gpu.jl")
include("train20/Train20.jl")
# include("layers21/Layers21.jl")

# See if we have a gpu at initialization:
import AutoGrad, CUDA
function __init__()
    if CUDA.functional()
        if isempty(Knet.LibKnet8.libknet8)
            @warn "libknet8 library not found, some GPU functionality may not be available, try reinstalling Knet."
        end
        Knet.array_type[] = Knet.KnetArrays.KnetArray{Float32}
        AutoGrad.set_gc_function(Knet.KnetArrays.cuallocator[] ? Knet.AutoGrad_gpu.gcnode : Knet.AutoGrad_gpu.knetgcnode)
        if CUDA.has_nvml() # Pick the device with highest memory
            mem(d) = CUDA.NVML.memory_info(CUDA.NVML.Device(CUDA.uuid(d))).free
            CUDA.device!(argmax(Dict(d=>mem(d) for d in CUDA.devices())))
        end
    end
end

# Match export list with v1.3.9 for backward compatibility
using AutoGrad #: @diff, AutoGrad, Param, cat1d, grad, gradloss, params, value
using Knet.LibKnet8 #: libknet8, @knet8, @knet8r, gpu
using Knet.KnetArrays #: KnetArray, gc, knetgc, ka, setseed, seed!
using Knet.FileIO_gpu #: cpucopy, gpucopy
using Knet.Ops20 #: RNN, accuracy, batchnorm, bce, bmm, bnmoments, bnparams, conv4, deconv4, dropout, elu, invx, logistic, logp, logsoftmax, logsumexp, mat, nll, pool, relu, rnnforw, rnninit, rnnparam, rnnparams, selu, sigm, softmax, unpool, zeroone
using Knet.Train20 #: Adadelta, Adagrad, Adam, Momentum, Nesterov, Rmsprop, SGD, Sgd, adadelta, adadelta!, adagrad, adagrad!, adam, adam!, atype, bilinear, converge, converge!, gaussian, goldensection, hyperband, minibatch, momentum, momentum!, nesterov, nesterov!, optimizers, param, param0, progress, progress!, rmsprop, rmsprop!, sgd, sgd!, train!, training, update!, xavier, xavier_normal, xavier_uniform

export @diff, Adadelta, Adagrad, Adam, AutoGrad, Knet, KnetArray, Momentum, Nesterov, Param, RNN, Rmsprop, SGD, Sgd, accuracy, adadelta, adadelta!, adagrad, adagrad!, adam, adam!, batchnorm, bce, bilinear, bmm, bnmoments, bnparams, cat1d, conv4, converge, converge!, cpucopy, deconv4, dropout, elu, gaussian, goldensection, gpu, gpucopy, grad, gradloss, hyperband, invx, ka, knetgc, logistic, logp, logsoftmax, logsumexp, mat, minibatch, momentum, momentum!, nesterov, nesterov!, nll, optimizers, param, param0, params, pool, progress, progress!, relu, rmsprop, rmsprop!, rnninit, rnnparam, rnnparams, selu, setseed, sgd, sgd!, sigm, softmax, train!, training, unpool, update!, value, xavier, xavier_normal, xavier_uniform, zeroone

# This is assumed by some old scripts:
export rnnforw

end # module


text:
struct Settings
    horizon::Base.RefValue{Int}
    time_series_cache_size::Int
    warm_start::Base.RefValue{Bool}
    initial_time::Base.RefValue{Dates.DateTime}
    optimizer::Union{Nothing, MOI.OptimizerWithAttributes}
    direct_mode_optimizer::Bool
    optimizer_solve_log_print::Bool
    detailed_optimizer_stats::Bool
    calculate_conflict::Bool
    system_to_file::Bool
    initialize_model::Bool
    initialization_file::String
    deserialize_initial_conditions::Bool
    export_pwl_vars::Bool
    allow_fails::Bool
    check_numerical_bounds::Bool
    ext::Dict{String, Any}
end

function Settings(
    sys;
    initial_time::Dates.DateTime=UNSET_INI_TIME,
    time_series_cache_size::Int=IS.TIME_SERIES_CACHE_SIZE_BYTES,
    warm_start::Bool=true,
    horizon::Int=UNSET_HORIZON,
    optimizer=nothing,
    direct_mode_optimizer::Bool=false,
    optimizer_solve_log_print::Bool=false,
    detailed_optimizer_stats::Bool=false,
    calculate_conflict::Bool=false,
    system_to_file::Bool=true,
    initialize_model::Bool=true,
    initialization_file="",
    deserialize_initial_conditions::Bool=false,
    export_pwl_vars::Bool=false,
    allow_fails::Bool=false,
    check_numerical_bounds=true,
    ext=Dict{String, Any}(),
)
    if time_series_cache_size > 0 &&
       sys.data.time_series_storage isa IS.InMemoryTimeSeriesStorage
        @info "Overriding time_series_cache_size because time series is stored in memory"
        time_series_cache_size = 0
    end

    if isa(optimizer, MOI.OptimizerWithAttributes) || optimizer === nothing
        optimizer_ = optimizer
    elseif isa(optimizer, DataType)
        optimizer_ = MOI.OptimizerWithAttributes(optimizer)
    else
        error(
            "The provided input for optimizer is invalid. Provide a JuMP.OptimizerWithAttributes object or a valid Optimizer constructor (e.g. GLPK.Optimizer).",
        )
    end

    return Settings(
        Ref(horizon),
        time_series_cache_size,
        Ref(warm_start),
        Ref(initial_time),
        optimizer_,
        direct_mode_optimizer,
        optimizer_solve_log_print,
        detailed_optimizer_stats,
        calculate_conflict,
        system_to_file,
        initialize_model,
        initialization_file,
        deserialize_initial_conditions,
        export_pwl_vars,
        allow_fails,
        check_numerical_bounds,
        ext,
    )
end

function log_values(settings::Settings)
    text = Vector{String}()
    for (name, type) in zip(fieldnames(Settings), fieldtypes(Settings))
        val = getfield(settings, name)
        if type <: Base.RefValue
            val = val[]
        end
        push!(text, "$name = $val")
    end

    @debug "Settings: $(join(text, ", "))" _group = LOG_GROUP_OPTIMIZATION_CONTAINER
end

function copy_for_serialization(settings::Settings)
    vals = []
    for name in fieldnames(Settings)
        if name == :optimizer
            # Cannot guarantee that the optimizer can be serialized.
            val = nothing
        else
            val = getfield(settings, name)
        end

        push!(vals, val)
    end

    return deepcopy(Settings(vals...))
end

function restore_from_copy(
    settings::Settings;
    optimizer::Union{Nothing, MOI.OptimizerWithAttributes},
)
    vals = Dict{Symbol, Any}()
    for name in fieldnames(Settings)
        if name == :optimizer
            vals[name] = optimizer
        elseif name == :ext
            continue
        else
            val = getfield(settings, name)
            vals[name] = isa(val, Base.RefValue) ? val[] : val
        end
    end

    return vals
end

get_horizon(settings::Settings) = settings.horizon[]
get_initial_time(settings::Settings)::Dates.DateTime = settings.initial_time[]
get_optimizer(settings::Settings) = settings.optimizer
get_ext(settings::Settings) = settings.ext
get_warm_start(settings::Settings) = settings.warm_start[]
get_system_to_file(settings::Settings) = settings.system_to_file
get_initialize_model(settings::Settings) = settings.initialize_model
get_initialization_file(settings::Settings) = settings.initialization_file
get_deserialize_initial_conditions(settings::Settings) =
    settings.deserialize_initial_conditions
get_export_pwl_vars(settings::Settings) = settings.export_pwl_vars
get_check_numerical_bounds(settings::Settings) = settings.check_numerical_bounds
get_allow_fails(settings::Settings) = settings.allow_fails
get_optimizer_solve_log_print(settings::Settings) = settings.optimizer_solve_log_print
get_calculate_conflict(settings::Settings) = settings.calculate_conflict
get_detailed_optimizer_stats(settings::Settings) = settings.detailed_optimizer_stats
get_direct_mode_optimizer(settings::Settings) = settings.direct_mode_optimizer
use_time_series_cache(settings::Settings) = settings.time_series_cache_size > 0

function set_horizon!(settings::Settings, horizon::Int)
    settings.horizon[] = horizon
    return
end

function set_initial_time!(settings::Settings, initial_time::Dates.DateTime)
    settings.initial_time[] = initial_time
    return
end

function set_warm_start!(settings::Settings, warm_start::Bool)
    settings.warm_start[] = warm_start
    return
end


text:
#####
##### ordering of faces in terms of vertices
#####
function face_vertices(elem::Line)
    return 1,2
end

function find_face_nodes(elem::Tri,r,s,tol=50*eps())
    e1 = findall(@. abs(s+1)<tol)
    e2 = findall(@. abs(r+s)<tol)
    e3 = findall(@. abs(r+1)<tol)
    return e1,reverse(e2),reverse(e3)
end

function find_face_nodes(elem::Quad,r,s,tol=50*eps())
    e1 = findall(@. abs(s+1)<tol)
    e2 = findall(@. abs(r-1)<tol)
    e3 = findall(@. abs(s-1)<tol)
    e4 = findall(@. abs(r+1)<tol)
    return e1,e2,reverse(e3),reverse(e4)
end

function find_face_nodes(elem::Hex,r,s,t,tol=50*eps())
    fv1 = findall(@. abs(r+1) < tol)
    fv2 = findall(@. abs(r-1) < tol)
    fv3 = findall(@. abs(s+1) < tol)
    fv4 = findall(@. abs(s-1) < tol)
    fv5 = findall(@. abs(t+1) < tol)
    fv6 = findall(@. abs(t-1) < tol)
    return fv1,fv2,fv3,fv4,fv5,fv6
end

function find_face_nodes(elem::Tet,r,s,t,tol=50*eps())
    fv1 = findall(@. abs(s+1) < tol)
    fv2 = findall(@. abs(r+s+t+1) < tol)
    fv3 = findall(@. abs(r+1) < tol)
    fv4 = findall(@. abs(t+1) < tol)
    return fv1,fv2,fv3,fv4
end

# face vertices = face nodes of degree 1
face_vertices(elem) = find_face_nodes(elem,nodes(elem,1)...)

#####
##### face data for diff elements
#####

function map_face_nodes(elem::Tri, face_nodes)
    r1D = face_nodes
    e = ones(size(r1D)) # vector of all ones
    rf = [r1D; -r1D; -e];
    sf = [-e; r1D; -r1D];
    return rf,sf
end
function map_face_nodes(elem::Quad, face_nodes)
    r1D = face_nodes
    e = ones(size(r1D))
    rf = [r1D; e; -r1D; -e]
    sf = [-e; r1D; e; -r1D]
    return rf,sf
end
function map_face_nodes(elem::Hex, face_nodes...)
    r,s = face_nodes
    e = ones(size(r))
    rf = [-e; e; r; r; r; r]
    sf = [r; r; -e; e; s; s]
    tf = [s; s; s; s; -e; e]
    return rf,sf,tf
end
function map_face_nodes(elem::Tet, face_nodes...)
    r,s = face_nodes
    e = ones(size(r))
    rf = [r; r; -e; r]
    sf = [-e; s; r; s]
    tf = [s; -(e + r + s); s; -e]
    return rf,sf,tf
end

"""
    function inverse_trace_constant(rd::RefElemData)

Returns the degree-dependent constant in the inverse trace equality over the reference element (as
reported in ["GPU-accelerated dG methods on hybrid meshes"](https://doi.org/10.1016/j.jcp.2016.04.003)
by Chan, Wang, Modave, Remacle, Warburton 2016).

Can be used to estimate dependence of maximum stable timestep on degree of approximation.
"""
inverse_trace_constant(rd::RefElemData{1}) = (rd.N+1) * (rd.N+2) / 2
inverse_trace_constant(rd::RefElemData{2, Quad}) = (rd.N+1) * (rd.N + 2)
inverse_trace_constant(rd::RefElemData{3, Hex}) = 3 * (rd.N + 1) * (rd.N + 2) / 2
inverse_trace_constant(rd::RefElemData{1, Line, SBP{TensorProductLobatto}}) = rd.N * (rd.N + 1) / 2
inverse_trace_constant(rd::RefElemData{2, Quad, SBP{TensorProductLobatto}}) = rd.N * (rd.N + 1)
inverse_trace_constant(rd::RefElemData{3, Hex, SBP{TensorProductLobatto}}) = 3 * rd.N * (rd.N + 1) / 2

# precomputed
_inverse_trace_constants(rd::RefElemData{2, Tri, Polynomial}) = (6.0, 10.898979485566365, 16.292060161853993, 23.999999999999808, 31.884512140579055, 42.42373503225737, 52.88579066878113, 66.25284319164409, 79.3535377715693, 95.53911875636945)
_inverse_trace_constants(rd::RefElemData{3, Tet, Polynomial}) = (10., 16.892024376045097, 23.58210016200093, 33.828424659883034, 43.40423356477473, 56.98869932201791, 69.68035962892684)
inverse_trace_constant(rd::RefElemData{2, Tri, Polynomial}) where {Dim} = _inverse_trace_constants(rd)[rd.N]
inverse_trace_constant(rd::RefElemData{3, Tet, Polynomial}) where {Dim} = _inverse_trace_constants(rd)[rd.N]

# generic fallback
function inverse_trace_constant(rd::RefElemData)
    return maximum(eigvals(Matrix(rd.Vf' * diagm(rd.wf) * rd.Vf), Matrix(rd.Vq' * diagm(rd.wq) * rd.Vq)))
end


text:
"""
    unwrap_laplacian(
        phas::AbstractArray{<:AbstractFloat, N ∈ (3, 4)},
        mask::AbstractArray{Bool, 3},
        vsz::NTuple{3, Real};
        solver::Symbol = :mgpcg
    ) -> typeof(similar(phas))

Laplacian phase unwrapping [1].

The Laplacian is computed using second order central finite differences
on the complex phase.
The resulting Poisson's equation can be solved with homogeneous Dirichlet
(`solver = :mgpcg`), Neumann (`:dct`), or periodic (`:fft`) boundary
conditions (BCs). Neumann and periodic BCs are imposed on the array while
Dirichlet BCs are imposed on a region-of-interest (ROI) (`mask`). The boundary
of the ROI is set such that values outside of it (`mask = 0`) are taken as
boundary points and values inside of it (`mask = 1`) as interior points, ie.
BC: `uphas[!mask] = 0`. This method combines phase unwrapping [1] and
harmonic background field removing [2].

### Arguments
- `phas::AbstractArray{<:AbstractFloat, N ∈ (3, 4)}`: wrapped (multi-echo) phase
- `mask::AbstractArray{Bool, 3}`: binary mask of region of interest
- `vsz::NTuple{3, Real}`: voxel size

### Keywords
- `solver::Symbol = :mgpcg`: solver for Poisson equation
    - `:dct`: homogeneous Neumann boundary condition on array
    - `:fft`: periodic boundary condition on array
    - `:mgpcg`: homogeneous Dirichlet boundary condition on `mask`
        (multigrid-preconditioned conjugate gradient method)

### Returns
- `typeof(similar(phas))`: unwrapped phase

### References
[1] Schofield MA, Zhu Y. Fast phase unwrapping algorithm for interferometric
    applications. Optics letters. 2003 Jul 15;28(14):1194-6.

[2] Zhou D, Liu T, Spincemaille P, Wang Y. Background field removal by solving
    the Laplacian boundary value problem. NMR in Biomedicine. 2014 Mar;27(3):312-9.
"""
function unwrap_laplacian(
    phas::AbstractArray{T, N},
    mask::AbstractArray{Bool, 3},
    vsz::NTuple{3, Real};
    solver::Symbol = :mgpcg
) where {T<:AbstractFloat, N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    size(mask) == size(phas)[1:3] || throw(DimensionMismatch())

    solver ∈ (:dct, :fft, :mgpcg) ||
        throw(ArgumentError("solver must be one of :dct, :fft, :mgpcg"))

    d2uphas = wrapped_laplacian(phas, vsz)

    if solver == :dct
        d2uphas = wrapped_laplacian_boundary_neumann!(d2uphas, phas, vsz)
        uphas = solve_poisson_dct(d2uphas, vsz)

    elseif solver == :fft
        d2uphas = wrapped_laplacian_boundary_periodic!(d2uphas, phas, vsz)
        uphas = solve_poisson_fft(d2uphas, vsz)

    elseif solver == :mgpcg
        # stc/utils/poisson_solver/mgpcg.jl
        nlevels = __DEFAULT_MGPCG_DEPTH(mask)

        opts = (
            # stc/bgremove/lbv.jl
            presmoother  = __DEFAULT_LBV_PRE(nlevels),
            coarsesolver = __DEFAULT_LBV_COARSE(),
            postsmoother = __DEFAULT_LBV_POST(nlevels),
            atol = sqrt(eps(Float64)),
            rtol = sqrt(eps(Float64)),
            maxit = maximum(size(mask)),
        )

        # set boundaries
        @inbounds for t in axes(d2uphas, 4)
            d2ut = @view(d2uphas[:,:,:,t])
            @batch for I in eachindex(d2ut, mask)
                d2ut[I] *= mask[I]
            end
        end

        uphas = solve_poisson_mgpcg(d2uphas, mask, vsz; opts...)
    end

    return uphas
end

function wrapped_laplacian(
    u::AbstractArray{<:AbstractFloat, N},
    dx::NTuple{3, Real}
) where {N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    return wrapped_laplacian!(tzero(u), u, dx)
end

function wrapped_laplacian!(
    d2u::AbstractArray{<:AbstractFloat, 3},
    u::AbstractArray{T, 3},
    dx::NTuple{3, Real}
) where {T<:AbstractFloat}
    size(d2u) == size(u) || throw(DimensionMismatch())

    τ = convert(T, 2π)

    nx, ny, nz = size(u)
    dx2 = convert(T, inv(dx[1]*dx[1]))
    dy2 = convert(T, inv(dx[2]*dx[2]))
    dz2 = convert(T, inv(dx[3]*dx[3]))

    tsz = padded_tilesize(T, (2, 2, 2), 1)
    R = vec(collect(TileIterator((2:nx-1, 2:ny-1, 2:nz-1), tsz)))

    @inbounds @batch per=thread for (I, J, K) in R
        for k in K
            for j in J
                for i in I
                    du = u[i,j,k] - u[i,j,k-1]
                    du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                    Δ = -dz2*du

                    du = u[i,j,k] - u[i,j-1,k]
                    du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                    Δ = muladd(-dy2, du, Δ)

                    du = u[i,j,k] - u[i-1,j,k]
                    du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                    Δ = muladd(-dx2, du, Δ)

                    du = u[i+1,j,k] - u[i,j,k]
                    du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                    Δ = muladd(dx2, du, Δ)

                    du = u[i,j+1,k] - u[i,j,k]
                    du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                    Δ = muladd(dy2, du, Δ)

                    du = u[i,j,k+1] - u[i,j,k]
                    du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                    Δ = muladd(dz2, du, Δ)

                    d2u[i,j,k] = Δ
                end
            end
        end
    end

    return d2u
end

function wrapped_laplacian!(
    d2u::AbstractArray{<:AbstractFloat, 4},
    u::AbstractArray{T, 4},
    dx::NTuple{3, Real}
) where {T<:AbstractFloat}
    size(d2u) == size(u) || throw(DimensionMismatch())

    τ = convert(T, 2π)

    nx, ny, nz, _ = size(u)
    dx2 = convert(T, inv(dx[1]*dx[1]))
    dy2 = convert(T, inv(dx[2]*dx[2]))
    dz2 = convert(T, inv(dx[3]*dx[3]))

    tsz = padded_tilesize(T, (2, 2, 2), 1)
    R = vec(collect(TileIterator((2:nx-1, 2:ny-1, 2:nz-1), tsz)))

    @inbounds for t in axes(u, 4)
        _u = @view(u[:,:,:,t])
        _d2u = @view(d2u[:,:,:,t])
        @batch per=thread for (I, J, K) in R
            for k in K
                for j in J
                    for i in I
                        du = _u[i,j,k] - _u[i,j,k-1]
                        du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                        Δ = -dz2*du

                        du = _u[i,j,k] - _u[i,j-1,k]
                        du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                        Δ = muladd(-dy2, du, Δ)

                        du = _u[i,j,k] - _u[i-1,j,k]
                        du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                        Δ = muladd(-dx2, du, Δ)

                        du = _u[i+1,j,k] - _u[i,j,k]
                        du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                        Δ = muladd(dx2, du, Δ)

                        du = _u[i,j+1,k] - _u[i,j,k]
                        du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                        Δ = muladd(dy2, du, Δ)

                        du = _u[i,j,k+1] - _u[i,j,k]
                        du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                        Δ = muladd(dz2, du, Δ)

                        _d2u[i,j,k] = Δ
                    end
                end
            end
        end
    end

    return d2u
end

@inline __wrap(x::T) where {T} = x < T(-π) || x > T(π)

#####
##### Boundary treatment
#####

# unroll? probably too aggressive
function wrapped_laplacian_boundary_neumann!(
    d2u::AbstractArray{<:AbstractFloat, N},
    u::AbstractArray{T, N},
    dx::NTuple{3, Real}
) where {T<:AbstractFloat, N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    size(d2u) == size(u) || throw(DimensionMismatch())

    τ = convert(T, 2π)

    sz = size(u)
    nx, ny, nz = sz[1:3]

    dx2 = convert(T, inv(dx[1]*dx[1]))
    dy2 = convert(T, inv(dx[2]*dx[2]))
    dz2 = convert(T, inv(dx[3]*dx[3]))

    _zero = zero(T)

    outer = CartesianIndices(ntuple(n -> 1:sz[n], Val(3)))
    inner = CartesianIndices(ntuple(n -> 2:sz[n]-1, Val(3)))
    E = EdgeIterator(outer, inner)

    R = Vector{NTuple{3, Int}}(undef, length(E))
    @inbounds for (i, I) in enumerate(E)
        R[i] = Tuple(I)
    end

    @inbounds for t in axes(u, 4)
        _u = @view(u[:,:,:,t])
        _d2u = @view(d2u[:,:,:,t])
        @batch per=thread for (i, j, k) in R
            Δ = _zero

            if k > 1
                du = _u[i,j,k] - _u[i,j,k-1]
                du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                Δ = muladd(-dz2, du, Δ)
            end

            if j > 1
                du = _u[i,j,k] - _u[i,j-1,k]
                du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                Δ = muladd(-dy2, du, Δ)
            end

            if i > 1
                du = _u[i,j,k] - _u[i-1,j,k]
                du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                Δ = muladd(-dx2, du, Δ)
            end

            if i < nx
                du = _u[i+1,j,k] - _u[i,j,k]
                du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                Δ = muladd(dx2, du, Δ)
            end

            if j < ny
                du = _u[i,j+1,k] - _u[i,j,k]
                du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                Δ = muladd(dy2, du, Δ)
            end

            if k < nz
                du = _u[i,j,k+1] - _u[i,j,k]
                du = __wrap(du) ? rem(du, τ, RoundNearest) : du
                Δ = muladd(dz2, du, Δ)
            end

            _d2u[i,j,k] = Δ
        end
    end

    return d2u
end

function wrapped_laplacian_boundary_periodic!(
    d2u::AbstractArray{<:AbstractFloat, N},
    u::AbstractArray{T, N},
    dx::NTuple{3, Real}
) where {T<:AbstractFloat, N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    size(d2u) == size(u) || throw(DimensionMismatch())

    τ = convert(T, 2π)

    sz = size(u)
    nx, ny, nz = sz[1:3]

    dx2 = convert(T, inv(dx[1]*dx[1]))
    dy2 = convert(T, inv(dx[2]*dx[2]))
    dz2 = convert(T, inv(dx[3]*dx[3]))

    outer = CartesianIndices(ntuple(n -> 1:sz[n], Val(3)))
    inner = CartesianIndices(ntuple(n -> 2:sz[n]-1, Val(3)))
    E = EdgeIterator(outer, inner)

    R = Vector{NTuple{3, Int}}(undef, length(E))
    @inbounds for (i, I) in enumerate(E)
        R[i] = Tuple(I)
    end

    @inbounds for t in axes(u, 4)
        _u = @view(u[:,:,:,t])
        _d2u = @view(d2u[:,:,:,t])
        @batch per=thread for (i, j, k) in R
            du = k == 1 ? _u[i,j,k] - _u[i,j,end] : _u[i,j,k] - _u[i,j,k-1]
            du = __wrap(du) ? rem(du, τ, RoundNearest) : du
            Δ = -dz2 * du

            du = j == 1 ? _u[i,j,k] - _u[i,end,k] : _u[i,j,k] - _u[i,j-1,k]
            du = __wrap(du) ? rem(du, τ, RoundNearest) : du
            Δ = muladd(-dy2, du, Δ)

            du = i == 1 ? _u[i,j,k] - _u[end,j,k] : _u[i,j,k] - _u[i-1,j,k]
            du = __wrap(du) ? rem(du, τ, RoundNearest) : du
            Δ = muladd(-dx2, du, Δ)

            du = i == nx ? _u[1,j,k] - _u[i,j,k] : _u[i+1,j,k] - _u[i,j,k]
            du = __wrap(du) ? rem(du, τ, RoundNearest) : du
            Δ = muladd(dx2, du, Δ)

            du = j == ny ? _u[i,1,k] - _u[i,j,k] : _u[i,j+1,k] - _u[i,j,k]
            du = __wrap(du) ? rem(du, τ, RoundNearest) : du
            Δ = muladd(dy2, du, Δ)

            du = k == nz ? _u[i,j,1] - _u[i,j,k] : _u[i,j,k+1] - _u[i,j,k]
            du = __wrap(du) ? rem(du, τ, RoundNearest) : du
            Δ = muladd(dz2, du, Δ)

            _d2u[i,j,k] = Δ
        end
    end

    return d2u
end


text:
export LSFManager, addprocs_lsf

struct LSFManager <: ClusterManager
    np::Integer
    bsub_flags::Cmd
    ssh_cmd::Cmd
    retry_delays
    throttle::Integer
end

struct LSFException <: Exception
    msg
end

function bpeek(manager, jobid, iarray)
    old_stderr = stderr
    rd,_ = redirect_stderr()
    try
        io = open(`$(manager.ssh_cmd) bpeek $(jobid)\[$iarray\]`)
        success(io) || throw(LSFException(String(readavailable(rd))))
        return io
    finally
        redirect_stderr(old_stderr)
    end
end

function _launch(manager, launched, c, jobid, iarray)
    config = WorkerConfig()

    io = retry(()->bpeek(manager, jobid, iarray),
               delays=manager.retry_delays,
               check=(s,e)->occursin("Not yet started", e.msg))()
    port_host_regex = r"julia_worker:([0-9]+)#([0-9.]+)"
    for line in eachline(io)
        mm = match(port_host_regex, line)
        isnothing(mm) && continue
        config.host = mm.captures[2]
        config.port = parse(Int, mm.captures[1])
        break
    end
    config.userdata = `$jobid\[$iarray\]`

    push!(launched, config)
    notify(c)
end

function launch(manager::LSFManager, params::Dict, launched::Array, c::Condition)
    try
        dir = params[:dir]
        exename = params[:exename]
        exeflags = params[:exeflags]

        np = manager.np

        jobname = `julia-$(getpid())`

        cmd = `$exename $exeflags $(worker_arg())`
        bsub_cmd = `$(manager.ssh_cmd) bsub $(manager.bsub_flags) -cwd $dir -J $(jobname)\[1-$np\] "$cmd"`

        line = open(readline, bsub_cmd)
        m = match(r"Job <([0-9]+)> is submitted", line)
        jobid = m.captures[1]

        asyncmap((i)->_launch(manager, launched, c, jobid, i),
                 1:np;
                 ntasks=manager.throttle)

    catch e
        println("Error launching workers")
        println(e)
    end
end

manage(manager::LSFManager, id::Int64, config::WorkerConfig, op::Symbol) = nothing

kill(manager::LSFManager, id::Int64, config::WorkerConfig) = remote_do(exit, id)

"""
    addprocs_lsf(np::Integer;
                 bsub_flags::Cmd=``,
                 ssh_cmd::Cmd=``,
                 retry_delays=ExponentialBackOff(n=10,
                                                 first_delay=1, max_delay=512,
                                                 factor=2),
                 throttle::Integer=np,
                 params...) =

Launch `np` workers on a cluster managed by IBM's Platform Load Sharing
Facility.  `bsub_flags` can be used to pass flags to `bsub` that are specific
to your cluster or workflow needs.  `ssh_cmd` can be used to launch workers
from other than the cluster head node (e.g. your personal workstation).
`retry_delays` is a vector of numbers specifying in seconds how long to
repeatedly wait for a worker to start.  `throttle` specifies how many workers
to launch at once.

# Examples

```
addprocs_lsf(1000; ssh_cmd=`ssh login`, throttle=10)
```
"""
addprocs_lsf(np::Integer;
             bsub_flags::Cmd=``,
             ssh_cmd::Cmd=``,
             retry_delays=ExponentialBackOff(n=10,
                                             first_delay=1, max_delay=512,
                                             factor=2),
             throttle::Integer=np,
             params...) =
        addprocs(LSFManager(np, bsub_flags, ssh_cmd, retry_delays, throttle); params...)


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module Vis

# If using precompiled system image (which we always are) you have to run AbstractPlotting.__init__() after loading Makie
# during precompilation, the display stack gets shuffled around such that the Makie display does not take priority.
# See https://discourse.julialang.org/t/makie-doesnt-display-plot-when-using-a-custom-julia-sysimage/38515.
__init__() = AbstractPlotting.__init__()

include("Visualization.jl")
include("Emitters.jl")

end # module Vis
export Vis


text:
"""
Subroutine dprsrch

This subroutine uses a projected search to compute a step
that satisfies a sufficient decrease condition for the quadratic

    q(s) = 0.5*s'*A*s + g'*s,

where A is a symmetric matrix in compressed column storage,
and g is a vector. Given the parameter alpha, the step is

    s[alpha] = P[x + alpha*w] - x,

where w is the search direction and P the projection onto the
n-dimensional interval [xl,xu]. The final step s = s[alpha]
satisfies the sufficient decrease condition

    q(s) <= mu_0*(g'*s),

where mu_0 is a constant in (0,1).

The search direction w must be a descent direction for the
quadratic q at x such that the quadratic is decreasing
in the ray x + alpha*w for 0 <= alpha <= 1.

MINPACK-2 Project. March 1999.
Argonne National Laboratory.
Chih-Jen Lin and Jorge J. More'.
"""
function dprsrch(n,x,xl,xu,A,g,w,wa1,wa2)
    one = 1.0
    p5 = 0.5

    # Constant that defines sufficient decrease.
    mu0 = 0.01

    # Interpolation factor.

    interpf = 0.5

    # Set the initial alpha = 1 because the quadratic function is
    # decreasing in the ray x + alpha*w for 0 <= alpha <= 1.

    alpha = one
    nsteps = 0

    # Find the smallest break-point on the ray x + alpha*w.
    nbrpt,brptmin,brptmax = dbreakpt(n,x,xl,xu,w)

    search = true
    while (search && alpha > brptmin)

        # Calculate P[x + alpha*w] - x and check the sufficient
        # decrease condition.

        nsteps = nsteps + 1
        dgpstep(n,x,xl,xu,alpha,w,wa1)
        dssyax(n,A,wa1,wa2)
        gts = ddot(n,g,1,wa1,1)
        q = p5*ddot(n,wa1,1,wa2,1) + gts
        if q <= mu0*gts
            search = false
        else

            # This is a crude interpolation procedure that
            # will be replaced in future versions of the code.

            alpha = interpf*alpha
        end
    end

    # Force at least one more constraint to be added to the active
    # set if alpha < brptmin and the full step is not successful.
    # There is sufficient decrease because the quadratic function
    # is decreasing in the ray x + alpha*w for 0 <= alpha <= 1.

    if (alpha < one && alpha < brptmin)
        alpha = brptmin
    end

    # Compute the final iterate and step.

    dgpstep(n,x,xl,xu,alpha,w,wa1)
    daxpy(n,alpha,w,1,x,1)
    dmid(n,x,xl,xu)
    dcopy(n,wa1,1,w,1)

    return
end


text:
using PhylogeneticTrees, JuMP, Gurobi
using CSV, DataFrames

paramsdf = try CSV.read(ARGS[1]) catch; CSV.read("params-miqo.csv") end
values = [try parse(Int, paramsdf[i,:value]) catch; paramsdf[i,:value] end for i in 1:nrow(paramsdf)]
params = Dict(zip(paramsdf[:param], values))

env = Gurobi.Env()
admix = (params["granularity"] > 1) && (params["admixture_events"] > 0)

# load data
pd = PhylogeneticTrees.PopulationData(params["mean_file"], params["cov_file"])

bt = PhylogeneticTrees.BinaryTree(params["depth"])
tp = PhylogeneticTrees.TreeProblem(pd, bt,
    solver = GurobiSolver(
        env,
        LogFile = params["log_file"],
        TimeLimit = params["time_limit"]),
    nlevels = params["granularity"],
    nmixtures = pd.npop + params["admixture_events"])
if admix
    if Bool(params["warm_start"])
        PhylogeneticTrees.warmstartunmixed(tp,
            solver = GurobiSolver(
                env,
                LogFile = params["log_file"],
                TimeLimit = params["warm_start_time_limit"]))
    end
    if !ismissing(params["unmixed_pops"])
        unmixpops = split(params["unmixed_pops"])
        for pop in unmixpops
            PhylogeneticTrees.unmix(tp, String(pop))
        end
    end
end
solve(tp.model)
PhylogeneticTrees.printnodes(tp, outfile = params["output_file"])


text:
using DashBootstrapComponents, DashHtmlComponents

list_group = html_div([
    dbc_listgroup([
        dbc_listgroupitem("Internal link", href="/l/components/list_group"),
        dbc_listgroupitem("External link", href="https://google.com"),
        dbc_listgroupitem("Disabled link", href="https://google.com", disabled=true),
        dbc_listgroupitem("Button", id="button-item", n_clicks=0, action=true),
    ]),
    html_p(id="counter"),
]);

callback!(app, Output("counter", "children"), Input("button-item", "n_clicks")) do n
    return "Button clicked $n times"
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
#!/usr/bin/env sh
#SBATCH -J bwbench
#SBATCH -N 1
#SBATCH -n 128
#SBATCH -A pc2-mitarbeiter
#SBATCH -p ALL
#SBATCH -t 01:30:00
#SBATCH --hint=nomultithread
#=
ml load lang
ml load Julia/1.7.2-linux-x86_64
export JULIA_NUM_THREADS=128
julia --project $(scontrol show job $SLURM_JOBID | awk -F= '/Command=/{print $2}')
exit
# =#

println("Node: ", gethostname(), "\n")

using ThreadPinning
using BandwidthBenchmark
using UnicodePlots

# Compact Pinning
println("----------------- Compact Pinning -----------------")
pinthreads(:compact)
threadinfo(; color=false)
println("Check compact pinning: ", getcpuids(), "\n")

# bwbench(; verbose=true);
println("bwbench():")
bwbench(; verbose=true, write_allocate=true);

println("bwscaling():")
data = bwscaling()
p = lineplot(data[:, 1], data[:, 2], title="Bandwidth Scaling, Compact Pinning", xlabel="# cores", ylabel="MB/s", border=:ascii, canvas=AsciiCanvas)
println(p)

println("flopsscaling():")
data = BandwidthBenchmark.flopsscaling()
p = lineplot(data[:, 1], data[:, 2], title="Flops Scaling, Compact Pinning", xlabel="# cores", ylabel="MFlops/s", border=:ascii, canvas=AsciiCanvas)
println(p)

# Scatter Pinning
println("----------------- Scatter Pinning -----------------")
pinthreads(:scatter)
threadinfo(; color=false)
println("Check scatter pinning: ", getcpuids(), "\n")

# bwbench(; verbose=true);
println("bwbench():")
bwbench(; verbose=true, write_allocate=true);

println("bwscaling():")
data = bwscaling()
p = lineplot(data[:, 1], data[:, 2], title="Bandwidth Scaling, Scattered Pinning", xlabel="# cores", ylabel="MB/s", border=:ascii, canvas=AsciiCanvas)
println(p)

println("flopsscaling():")
data = BandwidthBenchmark.flopsscaling()
p = lineplot(data[:, 1], data[:, 2], title="Flops Scaling, Scattered Pinning", xlabel="# cores", ylabel="MFlops/s", border=:ascii, canvas=AsciiCanvas)
println(p)

max_stars_repo_name: carstenbauer/BandwidthBenchmark

max_issues_repo_name: carstenbauer/BandwidthBenchmark

max_forks_repo_name: carstenbauer/BandwidthBenchmark


text:
using Printf
using Coverage

const default_pid = 52343

function analyze_dir(dir; pid=default_pid)
  pid_string = string(pid)
  filter(x->occursin(pid_string, x.filename), Coverage.analyze_malloc(dir))
end

# Projects containing code
dirs = [homedir() * "/projects/distributed_sensor_coverage",
        homedir() * "/projects/center_of_mass",
        homedir() * "/.julia/dev/MCTS"]

info = vcat(map(analyze_dir, dirs)...)

# sort by decreasing number of bytes
sort!(info, by=x->x.bytes, rev=true)

for line in info[1:30]
  @printf("%12d", line.bytes)
  println("  ", line.filename, "(", line.linenumber, ")")
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
## Put the file in a nicer form

L = readlines("../original-data/author-dblp/all_author_data.txt")

Author2ID = Dict()
Author2papers = Vector{Vector{Int64}}()
Author2Gender = Vector{Int64}()         # 1 is female, 0 is male

## Process line by line
nextID = 1  # current node id
removelist = Set{Int64}()   # set of papers to remove because we don't have gender data for an author
known = Set{Int64}()        # set of authors for which we have gender data
numpapers = 148251          # number of papers, known from processing data
GenderCounts = zeros(numpapers,3)
total_n = 0
total_w = 0
m = 148251
for j = 2:length(L)
    global nextID
    Lpiece = split(L[j],",")
    author_name = String(strip(Lpiece[4]))
    paperID = parse(Int64,Lpiece[1])
    gender = Lpiece[6]
    if occursin('M',gender)
        gid = 0
    elseif occursin('F',gender)
        gid = 1
    elseif occursin('-',gender)
        push!(removelist,paperID)
        gid = -1
    else
        @show gender
        println("oops")
    end

    id = get(Author2ID,author_name,0)
    if id == 0
        Author2ID[author_name] = nextID
        id = nextID
        nextID += 1

        # if new author, start new paper list
        paperlist = [paperID]
        push!(Author2papers,paperlist)

        # and add their gender
        push!(Author2Gender,gid)

        if gid == 0 || gid == 1
            push!(known,id)
        end
    else
        # add to existing paper list for this person
        push!(Author2papers[id], paperID)
    end
end

## Turn into hypergraph
using SparseArrays
n = nextID-1
I = Vector{Int64}()
J = Vector{Int64}()

for i = 1:n
    papers = Author2papers[i]

    for p in papers
        # put an edge in the bipartite representation for this author-paper
        push!(I,i)
        push!(J,p)
    end

end

B = sparse(J,I,ones(length(I)),m,n)

# Some authors show up in the author list twice, leading to a 2 in an
# entry of B. On two occassions, this is because two authors of the
# same name were on the same paper (often this is avoided by using numbers next
# to author, but evidently there's still some noise). In all 5 other cases,
# I checked the paper, and there's just a mistake in the data. So turning the
# 2 into a 1 is the right thing to do.

## Get rid of papers that have an unknown author on them, and get rid of
# authors of unknown gender. Save author names.
names = Vector{String}()
for j = 1:n
    push!(names,"")
end

for k in keys(Author2ID)
    id = Author2ID[k]
    names[id] = k
    # @show id,k
end

## save the data as a hypergraph

keepedges = setdiff(1:m,removelist)
H = B[keepedges,collect(known)]
m,n = size(H)
I,J,V = findnz(H)
H = sparse(I,J,ones(length(I)),m,n)
name = names[collect(known)]
class = Author2Gender[collect(known)]

matwrite("DBLP-Author-Gender_H.mat", Dict("H"=>H, "name"=>name,"class"=>class))

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily


text:
module extensive_form_game

# package code goes here
filename = "/Users/christiankroer/Documents/research/equilibrium_computation/efg_solve/games/kuhn.txt"

abstract Game
type ExtensiveFormGame
    domain::Tuple{Treeplex,Treeplex}
    tree::GameTree
    name::String

end

end # module


text:
struct Contact
    mechanism::Mechanism
    body::RigidBody
    point::Point3D
    obstacle::Obstacle
end

function contact_basis(normal::FreeVector3D, μ::Real, motion_type::Symbol)
    if motion_type == :xz
        return hcat(Rotations.RotY(π/2) * normal, Rotations.RotY(-π/2) * normal)
    elseif motion_type == :xy
        return hcat(Rotations.RotZ(π/2) * normal, Rotations.RotZ(-π/2) * normal)
    elseif motion_type == :yz
        return hcat(Rotations.RotX(π/2) * normal, Rotations.RotX(-π/2) * normal)
    elseif motion_type == :xyz
        R = Rotations.rotation_between([0., 0., 1.], normal.v)
        return hcat([R * RotZ(θ) * [1., 0., 0.] for θ in 0:π/2:3π/2]...)
    else
        throw(ArgumentError("Unrecognized motion type: $motion_type (should be :xy, :xz, :yz, or :xyz)"))
    end
end

mutable struct ContactJacobian{T}
    contact::Contact
    F::Matrix{Float64}
    ϕ::T
    N::Matrix{T}
    contact_rot::Matrix{T}
    contact_trans::Vector{T}
    geo_jacobian::GeometricJacobian{Matrix{T}}
    geo_jacobian_surface::Union{Nothing,GeometricJacobian{Matrix{T}}}
    P::Matrix{T}
    J::Matrix{T}

    function ContactJacobian{T}(contact::Contact) where T
        total_weight = mass(contact.mechanism) * norm(contact.mechanism.gravitational_acceleration)
        world = root_body(contact.mechanism)
        world_frame = default_frame(world)
        nv = num_velocities(contact.mechanism)
        β_dim = size(contact.obstacle.basis, 2)

        F = total_weight * hcat(contact.obstacle.normal.v, contact.obstacle.basis)

        ϕ = zero(T)
        N = Matrix{T}(undef, 1, nv)
        contact_rot = Matrix{T}(undef, 3, 3)
        contact_trans = Vector{T}(undef, 3)
        geo_jacobian = GeometricJacobian(world_frame, default_frame(contact.body), root_frame(contact.mechanism), Matrix{T}(undef, 3, nv), Matrix{T}(undef, 3, nv))
        geo_jacobian_surface = GeometricJacobian(world_frame, default_frame(contact.obstacle.body), root_frame(contact.mechanism), Matrix{T}(undef, 3, nv), Matrix{T}(undef, 3, nv))
        P = zeros(T, 3, 3)
        J = Matrix{T}(undef, nv, 1 + β_dim)

        new{T}(contact,F,ϕ,N,contact_rot,contact_trans,geo_jacobian,geo_jacobian_surface,P,J)
    end
end

function contact_jacobian!(cj::ContactJacobian, state::MechanismState)
    contact = cj.contact
    world = root_body(contact.mechanism)
    world_frame = default_frame(world)

    cp = transform(state, contact.point, contact.obstacle.normal.frame)

    cj.ϕ = separation(contact.obstacle, cp)
    cj.N .= contact.obstacle.normal.v' * point_jacobian(state, path(contact.mechanism, contact.obstacle.body, contact.body), cp).J

    cj.contact_rot = relative_transform(state, contact.obstacle.normal.frame, world_frame).mat[1:3,1:3]
    cj.contact_trans = relative_transform(state, contact.point.frame, world_frame).mat[1:3,4]

    p = contact.point.v + cj.contact_trans
    cj.P[1,2] = -p[3]
    cj.P[1,3] = p[2]
    cj.P[2,1] = p[3]
    cj.P[2,3] = -p[1]
    cj.P[3,1] = -p[2]
    cj.P[3,2] = p[1]

    geometric_jacobian!(cj.geo_jacobian, state, path(contact.mechanism, contact.body, world))
    cj.J .= (cj.geo_jacobian.linear' + cj.geo_jacobian.angular' * cj.P) * cj.contact_rot * cj.F
    if contact.obstacle.is_floating
        geometric_jacobian!(cj.geo_jacobian_surface, state, path(contact.mechanism, contact.obstacle.body, world))
        cj.J .+= -(cj.geo_jacobian_surface.linear' + cj.geo_jacobian_surface.angular' * cj.P) * cj.contact_rot * cj.F
    end
end

function contact_τ(cj::ContactJacobian, c_n, β)
    cj.J * vcat(c_n, β)
end


text:
@testset "ModelKit - Symbolic" begin
    @testset "SymEngine" begin
        @test Expression(MathConstants.catalan) isa Expression
        @test Expression(MathConstants.e) isa Expression
        @test Expression(MathConstants.pi) isa Expression
        @test Expression(MathConstants.γ) isa Expression
        @test Expression(Float16(2.1)) isa Expression
        @test Expression(Float32(2.1)) isa Expression
        @test Expression(big(2.1)) isa Expression
        @test Expression(Int32(2)) isa Expression
        @test Expression(Int16(2)) isa Expression
        @test Expression(Int128(2)) isa Expression
        @test Expression(2 // 3) isa Expression
        @test convert(BigFloat, Expression(big(2.1))) == big(2.1)
        @test ModelKit.is_number(Expression(2))
        @test ModelKit.to_number(Expression(big(2.1))) isa BigFloat
        @test ModelKit.to_number(Expression(2 // 3)) isa Rational{Int32}
        @test ModelKit.to_number(Expression(2) // Expression(3)) isa Rational{Int32}
        @test ModelKit.to_number(Expression(big(2)^129)) isa BigInt

        @test BigFloat(Expression(big(2.1))) == big(2.1)
        @test Float64(Expression(2)) == 2.0
        @test Rational{Int32}(Expression(2 // 3)) == 2 // 3
        @test UInt(Expression(2)) == UInt(2)

        @test complex(Expression(1), Expression(2)) == Expression(1) + im * 2
        @test complex(Expression(1), 2) == Expression(1) + im * 2
        @test complex(1, Expression(2)) == Expression(1) + im * 2
        @test +Expression(2) == 2
        @test Expression(HC.DoubleDouble.DoubleF64(2)) == big(2.0)
    end

    @testset "Variables" begin
        @var a b x[1:2] y[1:2, 1:3]

        @test a isa Variable
        @test b isa Variable
        @test x isa Vector{Variable}
        @test length(x) == 2
        @test y isa Matrix{Variable}
        @test size(y) == (2, 3)
        @test join(x, " ") == "x₁ x₂"

        @test reshape(sort(vec(y)), 2, 3) == y
        @var x[9:11]
        @test sort(x) == x

        @var c, d
        @test c isa Variable
        @test d isa Variable
        @test variables(c) == [c]
        @test c === c
        @test c !== copy(c)

        let
            c2 = @unique_var c
            @test c != c2
        end

        f = a + b
        @test variables(f) == [a, b]
        @test nvariables(f) == 2
        g = x[1]
        @test variables([f, g]) == [a, b, x[1]]
        @test nvariables([f, g]) == 3
    end

    @testset "Subs" begin
        @var x y w z u
        f = x^2 * (x + y * w)

        @test subs(f, x => z) == z^2 * (z + w * y)
        @test subs([f], x => z) == [z^2 * (z + w * y)]
        @test subs(f, [x, y] => [z^2, z + 2]) == z^4 * (w * (2 + z) + z^2)
        @test subs(f, [x, y] => [z^2, z + 2], w => u) == z^4 * (u * (2 + z) + z^2)
        @test subs(f, x => z^2, y => 3, w => u) == z^4 * (3 * u + z^2)
        @test subs(f, Dict(x => z^2, y => 3, w => u)) == z^4 * (3 * u + z^2)
    end

    @testset "Evaluation" begin
        @var x y w
        f = x^2 * (x + y * w)
        @test evaluate([f, f], [x, y, w] => [2, 3, -5]) == [-52, -52]
        @test [f, 2f]([x, y, w] => [2, 3, -5]) == [-52, -104]
    end

    @testset "Linear Algebra" begin
        @var x[1:2, 1:2]
        @test det(x[1:1, 1:1]) == x[1, 1]
        @test det(x) == -x[2, 1] * x[1, 2] + x[2, 2] * x[1, 1]
        @test det(x') == -x[2, 1] * x[1, 2] + x[2, 2] * x[1, 1]
        @test det(transpose(x)) == -x[2, 1] * x[1, 2] + x[2, 2] * x[1, 1]
    end

    @testset "Differentiation" begin
        @var x y

        f = x^2 + y^2
        g = x^3 + y^3

        @test differentiate(f, x) == 2x
        @test differentiate(f, [x, y]) == [2x, 2y]

        @test differentiate([f, g], x) == [2x, 3 * x^2]
        @test differentiate([f, g], [x, y]) == [2x 2y; 3*x^2 3*y^2]
    end

    @testset "Expand" begin
        @var x y
        @test expand((x + y)^2) == 2 * x * y + x^2 + y^2
    end

    @testset "to_dict" begin
        @var x y a
        @test ModelKit.to_dict(x^2 + y * x * a + y * x, [x, y]) ==
              Dict([2, 0] => Expression(1), [1, 1] => a + 1)
        @test ModelKit.degree(x * y^5 + y^2, [x, y]) == 6
    end

    @testset "Modeling" begin
        @testset "Bottleneck" begin
            @var x y z
            f = [
                (0.3 * x^2 + 0.5z + 0.3x + 1.2 * y^2 - 1.1)^2 +
                (0.7 * (y - 0.5x)^2 + y + 1.2 * z^2 - 1)^2 - 0.3,
            ]

            I = let x = [x, y, z]
                n, m = length(x), length(f)
                @unique_var y[1:n] v[1:m] w[1:m]
                J = differentiate(f, x)
                f′ = f(x => y)
                J′ = J(x => y)
                Nx = (x - y) - J' * v
                Ny = (x - y) - J′' * w
                System([f; f′; Nx; Ny], [x; y; v; w])
            end
            @test I isa System
            @test size(I) == (8, 8)
        end

        @testset "Steiner" begin
            @var x[1:2] a[1:5] c[1:6] y[1:2, 1:5]

            #tangential conics
            f = sum([a; 1] .* monomials([x; 1], 2; affine = false))
            ∇ = differentiate(f, x)
            #5 conics
            g = sum(c .* monomials([x; 1], 2; affine = false))
            ∇_2 = differentiate(g, x)
            #the general system
            #f_a_0 is tangent to g_b₀ at x₀
            function Incidence(f, a₀, g, b₀, x₀)
                fᵢ = f(x => x₀, a => a₀)
                ∇ᵢ = ∇(x => x₀, a => a₀)
                Cᵢ = g(x => x₀, c => b₀)
                ∇_Cᵢ = ∇_2(x => x₀, c => b₀)

                [fᵢ; Cᵢ; det([∇ᵢ ∇_Cᵢ])]
            end
            @var v[1:6, 1:5]
            I = vcat(map(i -> Incidence(f, a, g, v[:, i], y[:, i]), 1:5)...)
            F = System(I, [a; vec(y)], vec(v))
            @test size(F) == (15, 15)
        end

        @testset "Reach plane curve" begin
            @var x y
            f = (x^3 - x * y^2 + y + 1)^2 * (x^2 + y^2 - 1) + y^2 - 5
            ∇ = differentiate(f, [x; y]) # the gradient
            H = differentiate(∇, [x; y]) # the Hessian

            g = ∇ ⋅ ∇
            v = [-∇[2]; ∇[1]]
            h = v' * H * v
            dg = differentiate(g, [x; y])
            dh = differentiate(h, [x; y])

            ∇σ = g .* dh - ((3 / 2) * h) .* dg

            F = System([v ⋅ ∇σ; f], [x, y])
            @test size(F) == (2, 2)
        end
    end

    @testset "Horner" begin
        # use fano scheme as test
        n = 4
        s = 1
        @var x[1:n]
        ν = monomials(x, 5; affine = false)
        N = length(ν)
        q₀ = randn(ComplexF64, N - 1)
        @var q[1:N-1]
        F = sum(q[i] * ν[i] for i = 1:N-1) + 1

        @var a[1:n-1] b[1:n-1] t
        L = [a; 1] .* t + [b; 0]

        G = subs(F, x => L)

        FcapL = collect(values(ModelKit.to_dict(expand(G), [t])))

        @test all(FcapL) do f
            g = horner(f)
            expand(g) == f
        end
        @test all(FcapL) do f
            g = horner(f, [a; b])
            expand(g) == f
        end
    end

    @testset "Rand / dense poly" begin
        @var x y
        f, c = dense_poly([x, y], 3)
        @test length(c) == 10
        g = rand_poly(Float64, [x, y], 3)
        @test subs(f, c => coefficients(g, [x, y])) == g
        _, coeffs = exponents_coefficients(g, [x, y]; expanded = true)
        @test subs(f, c => coeffs) == g

        @var x[1:3]
        f, c = dense_poly(x, 3, coeff_name = :c)
        g = x[1]^3 + x[2]^3 + x[3]^3 - 1
        gc = coeffs_as_dense_poly(g, x, 3)
        @test subs(f, c => gc) == g
    end

    @testset "System" begin
        @var x y a b
        f = [(x + y)^3 + x^2 + x + 5y + 3a, 2 * x^2 + b]
        F = System(f, [x, y], [b, a])
        @test nvariables(F) == 2
        @test nparameters(F) == 2
        @test variables(F) == [x, y]
        @test parameters(F) == [b, a]
        show_F = """
        System of length 2
         2 variables: x, y
         2 parameters: b, a

         3*a + x + 5*y + x^2 + (x + y)^3
         b + 2*x^2"""
        @test sprint(show, F) == show_F
        @test degrees(F) == [3, 2]

        T = CompiledSystem(F; optimizations = false)
        F2 = System(T)
        @test F == F2
        @test sprint(show, T) == "Compiled: $show_F"
        @test size(T) == size(F) == (2, 2)

        T2 = CompiledSystem(F; optimizations = true)
        F3 = System(T2)
        @test expand.(F.expressions) == expand.(F3.expressions)

        F4 = System(f, parameters = [b, a])
        @test variables(F4) == [x, y]

        @var x y a b
        f = Any[(x+y)^3+x^2+x+5y+3a, 2*x^2+b]
        F = System(f, [x, y], [b, a])
        @test F isa System
    end

    @testset "System variables groups + homogeneous" begin
        @var x y z
        f = System([
            (x^2 + y^2 + x * y - 3 * z^2) * (x + 3z),
            (x^2 + y^2 + x * y - 3 * z^2) * (y - x + 2z),
            2x + 5y - 3z,
        ])
        @test is_homogeneous(f)
        multi_degrees(f) == [1 1; 2 0]

        @var x y z v w
        g = System([x * y - 2v * w, x^2 - 4 * v^2], variable_groups = [[x, v], [y, w]])
        @test is_homogeneous(g)
        multi_degrees(g) == [1 1; 2 0]
    end

    @testset "Homotopy" begin
        @var x y z t

        h = [x^2 + y + z + 2t, 4 * x^2 * z^2 * y + 4z - 6x * y * z^2]
        H = Homotopy(h, [x, y, z], t)

        @test nvariables(H) == 3
        @test nparameters(H) == 0
        @test variables(H) == [x, y, z]
        @test parameters(H) == []
        show_H = """
        Homotopy in t of length 2
         3 variables: x, y, z

         2*t + y + z + x^2
         4*z - 6*x*y*z^2 + 4*x^2*y*z^2"""

        @test sprint(show, H) == show_H

        H_any = Homotopy(convert(Vector{Any}, h), [x, y, z], t)
        @test H_any isa Homotopy

        T = CompiledHomotopy(H)
        H2 = ModelKit.interpret(T)
        @test H == H2
        @test sprint(show, T) == "Compiled: $show_H"
        @test size(T) == size(H) == (2, 3)
    end

    @testset "Optimizations for rational system" begin
        @var y[1:6] q[1:8]
        y₁, y₂, y₃, y₄, y₅, y₆ = y
        q₁, q₂, q₃, q₄, q₅, q₆, q₇, q₈ = q

        f =
            q₁ / y₁ - q₂ / (-y₁ + y₂) +
            q₅ * y₄ / (y₁ * y₄ - y₃ * y₂) +
            q₈ * y₆ / (y₁ * y₆ - y₂ * y₅)
        F = ModelKit.System([f], y, q)
        I = ModelKit.interpreter(ModelKit.optimize(F))
        u = Expression[0]
        ModelKit.execute!(u, I, y, q, ModelKit.InterpreterCache(Expression, I))

        @test expand(u[1] - f) == 0
    end
end


text:
# Place here all functions that are listed for deprication

function add_field(grid::Grid, field_name::String, field_type::String,
                                                                field_data)

  @warn("`add_field` function with no entry type argument is depricated."*
          " Entry type will default to `node`")

  add_field(grid, field_name, field_type, field_data, "node")
end

function calculate_field(grid::Grid, f, field_name::String, field_type::String)
  @warn("`calculate_field` function with no entry type argument is depricated."*
          " Entry type will default to `node`")

  calculate_field(grid, f, field_name, field_type, "node")
end