text:
using OrdinaryDiffEq
using StaticArrays
using Plots

"""
    Duffing

A Duffing equation with periodic forcing of the form

```math
mx″ + cx′ + kx + μx³ = F₀cos(ωt)
```
"""
Base.@kwdef mutable struct Duffing
    m::Float64 = 1.0
    c::Float64 = 0.01
    k::Float64 = 2.0
    μ::Float64 = 1
    F₀::Float64 = 0.05
    ω::Float64 = 1.0
end

"""
    rhs_dt(u, p::Duffing, t)

Return the right-hand side of the Duffing.
"""
function rhs_dt(u, p::Duffing, t)
    x = u[1]
    x′ = u[2]
    f = p.F₀ * cos(p.ω * t)
    # Use a static vector (SVector) for speed; it's not necessary for any other reason
    return SVector(x′, (f - p.c * x′ - p.k * x - p.μ * x^3) / p.m)
end

"""
    integrate(u, p::Duffing, t)

Return a solution of the Duffing equation over the time interval `t` starting at the initial
condition `u` with the parameters of the Duffing equation taken from `p`.
"""
function integrate(u, p::Duffing, t)
    # Explicitly convert everything to Float64 in case integer values are used
    # Use SVectors for speed (not necessary)
    prob = ODEProblem(rhs_dt, SVector{2, Float64}(u), Float64.(t), p)
    return solve(prob, Tsit5())
end

"""
    example()

Run an example Duffing equation simulation and plot the results.
"""
function example()
    p = Duffing()
    sol = integrate([0, 0], p, (0, 100))
    plot(sol.t, sol[1,:])
end

# Code to benchmark the timings

# using BenchmarkTools
# p = Duffing()
# @benchmark integrate($([0, 0]), $p, (0, 100))

# Output on my laptop

# BenchmarkTools.Trial:
#   memory estimate:  60.34 KiB
#   allocs estimate:  462
#   --------------
#   minimum time:     73.100 μs (0.00% GC)
#   median time:      75.700 μs (0.00% GC)
#   mean time:        90.955 μs (3.73% GC)
#   maximum time:     3.807 ms (96.36% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1

max_stars_repo_name: dawbarton/julia-de-examples

max_issues_repo_name: dawbarton/julia-de-examples

max_forks_repo_name: dawbarton/julia-de-examples


text:
module APL

export @apl_str, parse_apl, eval_apl, init_repl

include("parser.jl")
include("eval.jl")
include("show.jl")
include("repl.jl")

macro apl_str(str)
    parse_apl(str) |> eval_apl |> esc
end

__init__() = isdefined(Base, :active_repl) ? init_repl() : nothing

end # module

ext: jl


text:
# =========================================================================== #

# =========================================================================== #

#Implementing an x0 viable solution huristic with a gluton method
#receiving an array of lenght()=m , a matrix of size m*n and thoose size
#returing a Bool array of lenght=n
#Asked question:
#Mettre en place une heuristique de construction d’une x0 réalisable.

function Glouton(cost, matrix, n, m)

    #Creating a set of lines that will be avaluated
    actif= zeros(Bool, m) #if actif[j]=0 the line will be avaluated
    stop1= ones(Bool,m)

    #Creation the an array of activated variables
    varactive = ones(Bool, n)
    stop2= zeros(Bool , n)
#    println("\n varactive: ", varactive, "\n stop2: ", stop2)

    #creating my utility array for my m variables
    util = zeros(Float64, n)

    #initialisation de la Solution
    SOL = zeros(Float64, n)

#    println( "\n Matrix: ", matrix, "\n Cost: ", cost,"\n n= ", n, "   m= ", m,"\n actif: ", actif, "\n stop: ", stop1, "\n util: ", util, "\n SOL: ", SOL, "\n \n")

    ite = 0
    #Création de la solution
    while actif!=stop1 && varactive!=stop2

#      println("\n ++++++++++++++++++++++++++++++++++++++++")
#      println("\n itération: ", ite)
      #println("\n actif : ", actif, "\n stop1 : " , "\n varactive : " , varactive , "\n stop2 : ")
#      println("\n ++++++++++++++++++++++++++++++++++++++++")

        #Foction d'utilité

        util = Utilite(cost, matrix, actif, n, m , varactive, stop1, stop2)

        #Choix parmis les candidats de util
#        PosCandidat::Int64
        PosCandidat = PosMax(util,n)

        #On ajoute le candidat à la solution
        SOL[PosCandidat] = 1

        #Desactive! le candidat selectionné
        Desactive!(PosCandidat, matrix, actif, m, varactive,n)
        ite=ite+1
    end
    #println(calculz(SOL,cost,n))
    #return(SOL)
    Z = calculz(SOL,cost,n)
    println("Solution avant améliration: " , Z)
#    return(SOL, actif, Z)
    return(SOL, Z)

end

# =========================================================================== #

#Detremine an utility  based on an active (actif) set of matrix's lines
function Utilite(cost, matrix, actif, n, m , varactive, stop1, stop2)
#    println("==========================")
    util = zeros(Float64, n) #réinitialisation du vecteur

    #On each column of matrix
    for j=1:n
            if varactive[j] == 1
                    #For each line
                    for i=1:m
                        #checking if the line is active
                        if actif[i]==0 #&& i <= m
                            K=matrix[i,j]
                            util[j]=util[j]+K
#                            println("====== \n", "i: ", i ,"    j: ", j , "   k: ", K )
#                            println("\n varactive: ", varactive)
#                            println("\n actif : ", actif)
#                            println("\n util : ", util)
                        end
                    end

                # dividing the number if iterations of the variable j by its cost (not a zero)

                if util[j]!=false #On évide de divisé par 0...
#                    println("\n util avant modif: ", util)
                    util[j] =  cost[j] / util[j]
#                    println("\n util après modif: ", util)

                else
                    varactive[j] = 0
                end
            end
    end
#    println("\n Util: ", util)
#    println("\n \n ===================Utilite fini =============")
    return util
end

# =========================================================================== #

function PosMax(util,n) #On choisie un candidat parmie les différents coûts, par ordre lexicographique
#println("%%%%%%%%%%%%%%%%%%%%")
#    Pos::Int64
    Pos=0
    ValCan=0.0

    for i=1:n

        if util[i] > ValCan

            ValCan = util[i]
            Pos = i
        end
    end
    #println("\n \n %%%%%% Pos max fini %%%%% \n pos: ", Pos)
    return Pos
end

# =========================================================================== #

function Desactive!(PosCandidat, matrix, actif, m , varactive , n) #On désactive les lignes où est le candidat
#println("\n======================= desactive ==============================")
#println("~~~~~~~~~~~~~~~~~~~~~~~~")
    for i=1:m
#        println("\n Boucle i: ", i)
        if matrix[i,PosCandidat] == 1
#            println("\n If matrix[i,PosCandidat] == 1: " , matrix[i,PosCandidat])
            actif[i] = 1
#            println("\n actif[i] = 1: ", actif[i])
            for j=1:n
#                println("\n for j=1:n ; j: ", j)
                if matrix[i,j] ==1
#                    println("\n if matrix[i,j] ==1 :", matrix[i,j])
                    varactive[j] = 0
#                    println("\n varactive[j] = 0: ", varactive)
                end
            end
        end
    end
#    println("\n \n ~~~~~~~~~Desactive fini! ~~~~~~~\n ")
#    println("varactive: ", varactive)
#    println("actif: ", actif)
end

#=========================================================================#

function calculz(x,costs,m)
    z = 0
    for i in 1:m
        z+= x[i]*costs[i]
    end
    return z
end

max_stars_repo_path: Archives/Truc qui fonctionne/DM1_1.jl

max_issues_repo_path: Archives/Truc qui fonctionne/DM1_1.jl

max_forks_repo_path: Archives/Truc qui fonctionne/DM1_1.jl


text:
# Tests conversion among basis-structure representations

@testset "Test Basis Structure Representations" begin
    mb = Basis(SplineParams(linspace(-3, 3, 9), 0, 1),
               ChebParams(7, 1e-4, 10.0),
               LinParams(linspace(-4, 2, 11), 0))

    # construct evaluation points
    X, x123 = nodes(mb)

    # construct expanded, direct, and tensor basis-structure representation
    Φ_expanded = BasisMatrix(mb, Expanded(), X)
    Φ_direct = BasisMatrix(mb, Direct(), X)
    Φ_tensor = BasisMatrix(mb, Tensor(), x123)

    # construct expanded, direct, and tensor basis-structure representation with 1D basis
    Φ_expanded_1d = BasisMatrix(mb[1], Expanded(), X[:,1])
    Φ_direct_1d = BasisMatrix(mb[1], Direct(), X[:,1])

    @testset "test standard Base methods" begin
        # test == and ndims, multiD
        for Φ in (Φ_expanded, Φ_direct, Φ_tensor)
            @test  Φ  ==  Φ
            @test ndims(Φ) == 3
        end

        @test Φ_expanded != Φ_direct
        @test Φ_expanded != Φ_tensor
        @test Φ_tensor != Φ_direct
        @test Φ_expanded_1d != Φ_expanded

        for Φ in (Φ_direct, Φ_tensor)
            @test ( Φ_expanded  ==  Φ ) == false
        end

        # test == and ndims, 1D
        for Φ in (Φ_expanded_1d, Φ_direct_1d)
            @test  Φ  ==  Φ
            @test ndims(Φ) == 1
        end

    end

    @testset "test convert methods" begin
        @test  Φ_direct  ==  convert(Direct, Φ_tensor)
        @test  Φ_expanded  ==  convert(Expanded, Φ_direct)
        @test  Φ_expanded  ==  convert(Expanded, Φ_tensor)
        @test ==(Φ_expanded,
                 convert(Expanded, convert(Direct, Φ_tensor)))

        # test the no-op covert method
        @test convert(Expanded, Φ_expanded) == Φ_expanded
        @test convert(Direct, Φ_direct) == Φ_direct
        @test convert(Tensor, Φ_tensor) == Φ_tensor

        # test that we can't do expanded -> (direct|tensor) or direct -> tensor
        @test_throws MethodError convert(Direct, Φ_expanded)
        @test_throws MethodError convert(Tensor, Φ_expanded)
        @test_throws MethodError convert(Tensor, Φ_direct)
    end

    @testset "test internal tools" begin
        ## test bmat_type
        for (TP, TM) in [(SplineParams{Vector{Float64}}, SparseMatrixCSC{Float64,Int}),
                         (SplineParams{Vector{Float32}}, SparseMatrixCSC{Float32,Int}),
                         (LinParams{AbstractVector{Float16}}, SparseMatrixCSC{Float16,Int}),
                         (ChebParams{Complex{Float64}}, Matrix{Complex{Float64}})]
            @test TM == @inferred BasisMatrices.bmat_type(TP, one(eltype(TP)))
        end

        @test BasisMatrices.bmat_type(SplineSparse, mb) == AbstractMatrix{Float64}
        @test BasisMatrices.bmat_type(mb) == AbstractMatrix{Float64}

        @test BasisMatrices.bmat_type(SplineSparse, mb[1]) == SplineSparse{Float64,Int}
        @test BasisMatrices.bmat_type(mb[1]) == SparseMatrixCSC{Float64,Int}

        @test BasisMatrices.bmat_type(SplineSparse, mb[2]) == Matrix{Float64}
        @test BasisMatrices.bmat_type(mb[2]) == Matrix{Float64}

        @test BasisMatrices.bmat_type(SplineSparse, mb[3]) == SplineSparse{Float64,Int}
        @test BasisMatrices.bmat_type(mb[3]) == SparseMatrixCSC{Float64,Int}

        ## test _checkx
        # create test data
        xm = rand(1, 2)
        xv = rand(2)
        xvv = [rand(2) for i=1:2]

        @test BasisMatrices._checkx(2, xm) == xm
        @test BasisMatrices._checkx(2, xv) == reshape(xv, 1, 2)
        @test BasisMatrices._checkx(2, xvv) == xvv
        @test BasisMatrices._checkx(1, xv) == xv

        @test_throws ErrorException BasisMatrices._checkx(2, rand(1, 3))
        @test_throws ErrorException BasisMatrices._checkx(2, rand(3))
        @test_throws ErrorException BasisMatrices._checkx(2, [rand(2) for i=1:3])

        ## test check_convert
        for Φ in (Φ_expanded, Φ_direct, Φ_tensor)
            @test BasisMatrices.check_convert(Φ, zeros(1, 3)) == (3, 1, 3)
            @test_throws ErrorException BasisMatrices.check_convert(Φ, zeros(1, 2))
            @test_throws ErrorException BasisMatrices.check_convert(Φ, fill(-1, 1, 3))
        end

        ## test check_basis_structure
        @testset "test check_basis_structure" begin
            order1 = zeros(Int, 1, 2)
            out1 = BasisMatrices.check_basis_structure(2, xm, order1)
            @test out1 == (1, order1, order1, [1 1], xm)

            # test when order::Int
            @test BasisMatrices.check_basis_structure(2, xm, 0) == out1

            # test the reshape(order, 1, N) branch (isa(order, Vector))
            @test BasisMatrices.check_basis_structure(2, xm, [0, 0]) == out1

            # check N=1 --> order = fill(order, 1, 1) branch
            out_1d = BasisMatrices.check_basis_structure(1, xv, 0)
            order_1d_out = fill(0, 1, 1)
            @test out_1d == (1, order_1d_out, order_1d_out, fill(1, 1, 1), xv)

            # check m > 1 branch
            order_m2 = [0 0; 1 1; 0 -1]
            out_m2 = BasisMatrices.check_basis_structure(2, xm, order_m2)
            @test out_m2 == (3, order_m2, [0 -1], [2 3], xm)

            order2 = zeros(Int, 1, 5)  # should throw error
            @test_throws ErrorException BasisMatrices.check_basis_structure(2, X,
                                                                       order2)
        end

    end

    @testset "constructors" begin
        # test constructor methods with ABSR types instead of instances
        @test BasisMatrix(mb, Expanded, X) == Φ_expanded
        @test BasisMatrix(mb, Direct, X) == Φ_direct
        @test BasisMatrix(mb, Tensor, x123) == Φ_tensor
        @test BasisMatrix(mb[1], Expanded, X[:,1]) == Φ_expanded_1d
        @test BasisMatrix(mb[1], Direct, X[:,1]) == Φ_direct_1d

        # test extra type methods
        @test BasisMatrix(Void, mb, Expanded, X) == Φ_expanded
        @test BasisMatrix(Void, mb, Direct, X) == Φ_direct
        @test BasisMatrix(Void, mb, Tensor, x123) == Φ_tensor
        @test BasisMatrix(Void, mb[1], Expanded, X[:,1]) == Φ_expanded_1d
        @test BasisMatrix(Void, mb[1], Direct, X[:,1]) == Φ_direct_1d

        ## Test SplineSparseMethods
        @test BasisMatrix(SplineSparse, mb, Expanded, X) == Φ_expanded

        # this one will be different -- need to break it out
        bmd = @inferred BasisMatrix(SplineSparse, mb, Direct, X)
        @test isa(bmd.vals[1], SplineSparse{Float64,Int})
        @test bmd.vals[1] == Φ_direct.vals[1]  # this works even though they have different types

        @test isa(bmd.vals[2], Matrix{Float64})
        @test bmd.vals[2] == Φ_direct.vals[2]

        @test isa(bmd.vals[3], SplineSparse{Float64,Int})
        @test bmd.vals[3] == Φ_direct.vals[3]

        bmt = @inferred BasisMatrix(SplineSparse, mb, Tensor, x123)
        @test isa(bmt.vals[1], SplineSparse{Float64,Int})
        @test bmt.vals[1] == Φ_tensor.vals[1]  # this works even though they have different types

        @test isa(bmt.vals[2], Matrix{Float64})
        @test bmt.vals[2] == Φ_tensor.vals[2]

        @test isa(bmt.vals[3], SplineSparse{Float64,Int})
        @test bmt.vals[3] == Φ_tensor.vals[3]

    end

    @testset "Test from PR #25" begin
        basisμ = Basis(Cheb, 20, 0.0, 1.0)
        basisσ = Basis(Cheb, 20, 0.0, 1.0)
        basis = Basis(basisμ, basisσ)
        S, (μs, σs) = nodes(basis)
        bs = BasisMatrix(basis, Expanded(), S, [0 2])
        @test isa(bs, BasisMatrix{Expanded}) == true
    end

    # call show (which calls writemime) so we can get 100% coverage
    @testset "Printing" begin
        iob = IOBuffer()
        show(iob, Φ_tensor)
    end

    @testset "non Float64 eltype" begin
        pc = ChebParams(10, Float32(-3.0), Float32(1.0))
        pl = LinParams(Float32[0.1, 0.2, 0.3, 0.4, 0.5])
        ps = SplineParams(Float32[0.1, 0.2, 0.3, 0.4, 0.5], 0, 3)
        b = Basis((pc, pl, ps))
        bs = @inferred BasisMatrix(b, Direct(), nodes(b)[1])
        @test eltype(bs.vals) == AbstractMatrix{Float32}

        b2 = Basis((pl, ps))
        bs2 = @inferred BasisMatrix(b2, Direct(), nodes(b2)[1])
        @test eltype(bs2.vals) == SparseMatrixCSC{Float32,Int}

        bs3 = @inferred BasisMatrix(SplineSparse, b2, Direct(), nodes(b2)[1])
        @test eltype(bs3.vals) == SplineSparse{Float32,Int}
    end

    @testset "no ABSR" begin
        bmd = BasisMatrix(mb, X)
        bmt = BasisMatrix(mb, x123)
        @test isa(bmd, BasisMatrix{Direct})
        @test isa(bmt, BasisMatrix{Tensor})

        bmd2 = BasisMatrix(SplineSparse, mb, X)
        bmt2 = BasisMatrix(SplineSparse, mb, x123)
        @test isa(bmd2.vals[1], SplineSparse)
        @test isa(bmd2.vals[2], Matrix)
        @test isa(bmd2.vals[3], SplineSparse)
        @test isa(bmt2.vals[1], SplineSparse)
        @test isa(bmt2.vals[2], Matrix)
        @test isa(bmt2.vals[3], SplineSparse)
    end

end  # testset


text:
###
# UCS Oracle Tests
###
facts("suppFcnTest UCS") do
	muhat  = [0.0, 1.2]
	Gamma1 = .002
	Gamma2 = .05
	covbar = [ 1.  0.2
	           0.2 1.3]

	covhat = covbar - Gamma2 * eye(2)  #covbar = covhat + Gamma2 I
	w = UCSOracle(muhat, covhat, Gamma1, Gamma2, .1)

	zstar, ustar = suppFcn([1, 1], w, :Max)
	@fact zstar --> roughly(6.132331444671241)
	@fact ustar[1] --> roughly(2.4661657223356204)
	@fact ustar[2] --> roughly(3.66616572233562)

	zstar, ustar = suppFcn([1, 1], w, :Min)
	@fact zstar --> roughly(-3.7323314446712406)
	@fact ustar[1] --> roughly(-2.4661657223356204)
	@fact ustar[2] --> roughly(-1.266165722335620)
end

#doesn't use any bounds
facts("portTest UCS No bounds") do
	srand(8675309); data = randn(500, 2)
	w = UCSOracle(data, .1, .1, .1)
	portTest(w, -2.4630938710200345, [0.5253456488147422, 0.4746543511852578],TOL=1e-2)
end

#VG Comment this back in when we have QP support in test suite
# facts("portTest2  UCS bounds") do
# 	srand(8675309); data = randn(500, 2)
# 	w = UCSOracle(data, .1, .1, .1)
# 	portTest(w, -2.4604223389522866, [0.530985977658841, 0.46901402234115896], TOL=1e-6,
# 			unc_lower=[-1e6, -1e6], unc_upper=[1e6, 1e6])
# end


text:
# Script to obtain the distribution of separations for random particles

push!(LOAD_PATH, "./models")
push!(LOAD_PATH, ".")
push!(LOAD_PATH, "./SparseInverseProblems/src")
using SparseInverseProblemsMod
using SuperResModels
using TestCases
using Utils

using PyCall
@pyimport numpy as np
@pyimport matplotlib.pyplot as plt

### Manual: This code is to generate files that contain the distribution of the
### distances of the generated particles. To have an efficient Sampling for plotting.
### To modify the interval of interest, change the lowLimit and upperLimit variables.

### Static parameters and container for use.
f_c = 20
x_max = 1.0
num_x = 10*f_c
filter_x =  ones(2*f_c+1);
model_static = SuperResModels.Fourier1d(f_c, x_max, filter_x, num_x)
### Dynamic parameters and container for use.
K = 2
tau = 1.0/(K*1.0)
# the maximum velocity on the domain is given by x_max/2 + K tau v_max = 1
v_max = 0.5
num_v = 10
model_dynamic = SuperResModels.DynamicFourier1d(model_static, v_max, tau, K, num_v*K)

# How particles are generated
minimumNumPartic = 4
maximumNumPartic = 10

test_case = () -> TestCases.cloud_1d_full(x_max, v_max, 1, 1, K, tau, rand(minimumNumPartic:maximumNumPartic))

### We generate particles and we evaluate their separation
N = 1000000
separations = zeros(N)
for i = 1:N
    if mod(i,100000)== 0
	println(i)
    end
    (thetas, weights) = test_case()
    separations[i] = Utils.separation_val(thetas,K,tau,x_max)
end

#plt.figure()
#plt.hist(separations,1000)
#plt.show()

# Spatial resolution of this approximation
M = 500
if x_max != 1.0
    warn(" Beware, the interval where we decided to consider the separations was tuned
    for the specific case of x_max equal to 1, since you modified this value, check if
        upper and lowerLimit is adequate for your use, just plot an histogram of the 		separations.")
end
upperLimit = 0.06
lowerLimit = 0.001
binss = linspace(lowerLimit,upperLimit,M)
vals = np.histogram(separations, bins = binss )
vals = vals[1]/(binss[2]-binss[1])/N

if minimum(vals) <= 0.00000001
error(" There are bin intervals without particles ")
end

plt.figure()
plt.plot(binss[1:end-1], vals)
plt.show()

np.save("data/1Dsimulations/separationDistribVal_06interval.npy", vals)
np.save("data/1Dsimulations/separationDistribBins_06interval.npy", binss)


text:
using ImageCore
using Test, Suppressor

refambs = @suppress_out detect_ambiguities(ImageCore, Base)
using SuperPixels
ambs = @suppress_out detect_ambiguities(ImageCore, Base, SuperPixels)

@testset "SuperPixels.jl" begin
    # check if SuperPixels.jl introduces new ambiguities
    @test isempty(setdiff(ambs, refambs))

    include("types.jl")
    include("synthesize.jl")
end

nothing


text:
#Test publish and subscribe ability
#works alongside echonode.py
#typegeneration.jl must be run first

using .geometry_msgs.msg

const Nmsgs = 10
const rate = 20. #Hz
const msgs = PoseStamped[]
const refs = Array{Vector3}(undef, Nmsgs)
const t0 = to_nsec(get_rostime())

for i=1:Nmsgs
    refs[i] = Vector3(rand(3)...)
end

const ros_pub = Publisher("vectors", Vector3, queue_size = 10)
rossleep(Duration(3.0))

function publish_messages(pubobj, msgs, rate_hz)
    r = Rate(rate_hz)
    for msg in msgs
        publish(pubobj, msg)
        rossleep(r)
    end
    rossleep(Duration(1.0))
end

function pose_cb(msg::PoseStamped, msgs::Vector{PoseStamped})
    mtime = to_nsec(msg.header.stamp) - t0
    mtime > 0 && println("Message received, time: ",mtime," nanoseconds")
    if msg.header.stamp.secs > 1.0
        push!(msgs, msg)
        println("Got message #",msg.header.seq)
    end
end
pose_cb(PoseStamped(), msgs) #warm up run

const ros_sub = Subscriber("poses", PoseStamped, pose_cb, (msgs,), queue_size = 10)

#First message doesn't go out for some reason
publish(ros_pub, Vector3(1.1,2.2,3.3))
rossleep(Duration(1.0))

#Test messages
publish_messages(ros_pub, refs, 20.0)
rossleep(Duration(1.0))

println("Received ",length(msgs)," / ",Nmsgs)

@test length(msgs) == Nmsgs
for i=1:Nmsgs
    @test msgs[i].pose.position.x ≈ refs[i].x
    @test msgs[i].pose.position.y ≈ refs[i].y
    @test msgs[i].pose.position.z ≈ refs[i].z
end
empty!(msgs)


text:
# SPDX-License-Identifier: MIT

"""
    abstract type AbstractParameter

Generic type for all parameters. Subtype of `AbstractParameter` should
implement the property `μ` which is the model-dependent
"canonical form" of the parameter.
"""
abstract type AbstractParameter <: BMObject end

"""
    ParameterList{N,T<:AbstractParam}(m1, m2, m3, ...)

Store a list of (parameters for an object's attribute. `N` is the
number of the parameters and `T` is the type of the parameters.
The list is immutable.

# Constructor

    ParameterList(m1, m2, m3, ...)

"""
const ParameterList{N,T<:AbstractParameter} = NTuple{N,T}
ParameterList(m...) = tuple(m...)

"""
    getparams(obj)

Returns a list of a all the parameters in `obj` and its
attributes. Note that array of parameters and array of sub-models
should be stored as [`BayesParameterList`](@ref) and [`ModelList`](@ref)
respectively.
"""
function getparams(obj)
    params = Set{AbstractParameter}()
    for name in fieldnames(typeof(obj))
        prop = getproperty(obj, name)
        if typeof(prop) <: AbstractParameter
            push!(params, prop)
            push!.([params], getparams(prop))
        elseif typeof(prop) <: BMObjectList
            for item in prop
                if typeof(item) <: AbstractParameter
                    push!(params, item)
                end
                push!.([params], getparams(item))
            end
        else
            push!.([params], getparams(prop))
        end
    end
    [params...]
end

"""
    statistics(param)

Returns the sufficient statistics of the `param` for the model it
belongs to.
"""
statistics

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
module EarthMoversDistance

import JuMP
import COSMO, MathOptInterface
import FluPredictibility.BioTools
using BioSequences: BioSequence

export EMD

"""
	EMD(f,P,Q)

Compute the earth mover's distance between distributions `P` and `Q`, using `f` as a metric.
Enumerating `P` (resp. `Q`) should yield pairs `(x,p)` (resp. `(y,q)`) representing state `x` and its weight `p`. Distance between elements is `f(x,y)`.
"""
function EMD(f::Function,P,Q)
	# Building the distance matrix
	d = Array{Float64, 2}(undef, length(P), length(Q))
	for (i,(x,p)) in enumerate(P)
		for (j,(y,q)) in enumerate(Q)
			d[i,j] = f(x,y)
		end
	end
	# Building the JuMP model
	model = JuMP.Model(COSMO.Optimizer)
	JuMP.set_silent(model)
	JuMP.@variable(model, w[i=1:length(P), j=1:length(Q)] >= 0)
	for (i,(x,p)) in enumerate(P)
		JuMP.@constraint(model, sum(w[i,j] for j in 1:length(Q)) ==  p)
	end
	for (j,(y,q)) in enumerate(Q)
		JuMP.@constraint(model, sum(w[i,j] for i in 1:length(P)) ==  q)
	end
	JuMP.@objective(model, Min, sum(sum(w.*d)))
	JuMP.optimize!(model)
	# Looking at the solution
	if JuMP.termination_status(model) != MathOptInterface.OPTIMAL
		@warn "Termination status of optimization: $(JuMP.termination_status(model))"
	end
	flow = JuMP.value.(w)
	emd = JuMP.objective_value(model)

	return (emd=emd, flow=flow)
end

##
## Specialized functions for BioSequences and BioTools
##
"""
	EMD(f::Function, X::Array{<:BioSequence,1},Y::Array{<:BioSequence,1})
	EMD(X::Array{<:BioSequence,1},Y::Array{<:BioSequence,1})
"""
function EMD(f::Function, X::Array{<:BioSequence,1},Y::Array{<:BioSequence,1})
	P = sequence_distribution(X)
	Q = sequence_distribution(Y)
	ret = EMD(f, P, Q)
	flow = Dict()
	for (i,(x,p)) in enumerate(P)
		for(j,(y,q)) in enumerate(Q)
			flow[x,y] = ret.flow[i,j]
		end
	end
	return (emd=ret.emd, flow=flow), P, Q
end
EMD(X::Array{<:BioSequence,1},Y::Array{<:BioSequence,1}) = EMD(BioTools.hamming, X, Y)
EMD(f::Function, X::Array{<:BioTools.AbstractStrain,1}, Y::Array{<:BioTools.AbstractStrain,1}) = EMD(f,
	[x.seq for x in X], [y.seq for y in Y])
EMD(X::Array{<:BioTools.AbstractStrain,1}, Y::Array{<:BioTools.AbstractStrain,1}) = EMD([x.seq for x in X], [y.seq for y in Y])
EMD(P::AbstractDict{<:BioSequence,Float64}, Q::AbstractDict{<:BioSequence,Float64}) = EMD(BioTools.hamming, P, Q)

"""
	sequence_distribution(X::Array{<:BioSequence})
	sequence_distribution(X::Array{<:AbstractStrain})

Convenience function. Transform an array `X` of sequences or strains into a distribution form usable by the EMD functions, *i.e.* a `Dict{BioSequence => weight::Float64}`.
"""
function sequence_distribution(X::Array{<:BioSequence})
	d = Dict{BioSequence, Float64}()
	for s in X
		d[s] = get(d, s, 0.) + 1. /length(X)
	end
	return d
end
function sequence_distribution(X::Array{<:BioTools.Strain})
	d = Dict{BioSequence, Float64}()
	for s in X
		d[s.seq] = get(d, s.seq, 0.) + 1. /length(X)
	end
	return d
end

end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
# Using Base modules.
using Random

# Load a plotting library.
using Plots

# Load the distributions library.
using Distributions

# Set the true probability of heads in a coin.
p_true = 0.5

# Iterate from having seen 0 observations to 100 observations.
Ns = 0:100;

# Draw data from a Bernoulli distribution, i.e. draw heads or tails.
Random.seed!(12)
data = rand(Bernoulli(p_true), last(Ns))

# Here's what the first five coin flips look like:
data[1:5]

# Our prior belief about the probability of heads in a coin toss.
prior_belief = Beta(1, 1);

# Import StatsPlots for animating purposes.
using StatsPlots

# Make an animation.
animation = @gif for (i, N) in enumerate(Ns)

    # Count the number of heads and tails.
    heads = sum(data[1:i-1])
    tails = N - heads

    # Update our prior belief in closed form (this is possible because we use a conjugate prior).
    updated_belief = Beta(prior_belief.α + heads, prior_belief.β + tails)

    # Plotting
    plot(updated_belief,
        size = (500, 250),
        title = "Updated belief after $N observations",
        xlabel = "probability of heads",
        ylabel = "",
        legend = nothing,
        xlim = (0,1),
        fill=0, α=0.3, w=3)
    vline!([p_true])
end

# Load Turing and MCMCChains.
using Turing, MCMCChains

# Load the distributions library.
using Distributions

# Load StatsPlots for density plots.
using StatsPlots

@model coinflip(y) = begin

    # Our prior belief about the probability of heads in a coin.
    p ~ Beta(1, 1)

    # The number of observations.
    N = length(y)
    for n in 1:N
        # Heads or tails of a coin are drawn from a Bernoulli distribution.
        y[n] ~ Bernoulli(p)
    end
end;

# Settings of the Hamiltonian Monte Carlo (HMC) sampler.
iterations = 1000
ϵ = 0.05
τ = 10

# Start sampling.
chain = sample(coinflip(data), HMC(ϵ, τ), iterations, progress=false);

# Construct summary of the sampling process for the parameter p, i.e. the probability of heads in a coin.
p_summary = chain[:p]
plot(p_summary, seriestype = :histogram)

# Compute the posterior distribution in closed-form.
N = length(data)
heads = sum(data)
updated_belief = Beta(prior_belief.α + heads, prior_belief.β + N - heads)

# Visualize a blue density plot of the approximate posterior distribution using HMC (see Chain 1 in the legend).
p = plot(p_summary, seriestype = :density, xlim = (0,1), legend = :best, w = 2, c = :blue)

# Visualize a green density plot of posterior distribution in closed-form.
plot!(p, range(0, stop = 1, length = 100), pdf.(Ref(updated_belief), range(0, stop = 1, length = 100)),
        xlabel = "probability of heads", ylabel = "", title = "", xlim = (0,1), label = "Closed-form",
        fill=0, α=0.3, w=3, c = :lightgreen)

# Visualize the true probability of heads in red.
vline!(p, [p_true], label = "True probability", c = :red)

isdefined(Main, :TuringTutorials) && Main.TuringTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: leachim/TuringTutorials

max_issues_repo_name: leachim/TuringTutorials

max_forks_repo_name: leachim/TuringTutorials


text:
module day06

include("inputs.jl")
ex1 = inputs.ex1_d6()
ex2 = inputs.ex2_d6()
input = inputs.input_d6()

# Part 1 & 2
function count_yesses(input)
    yesses = [0, 0] # [p1:everyone's, p2:anyone's]
    for group in input
        testset = Set.(group)
        yesses += [ length(union(testset...)),      # p1
                    length(intersect(testset...))]  # p2
    end
    return yesses
end

# A more concise version of the above:
yeps = [sum(length(union(group...)) for group in input),
        sum(length(intersect(group...)) for group in input)]

@assert count_yesses(ex1)[1] == 6
@assert count_yesses(ex2)[1] == 11
@assert count_yesses(ex2)[2] == 6
@info "Day 6, Part 1 answer: $(count_yesses(input)[1])"
@info "Day 6, Part 2 answer: $(count_yesses(input)[2])"

end


text:
using SurrealNumbers: NotImplementedError, @notimplemented_default

@testset "notimplemented_default" begin
	@test @notimplemented_default(1 + 2, 0) == 3
	@test @notimplemented_default(throw(NotImplementedError(nothing, ())), 0) == 0
	@test_throws DivideError @notimplemented_default(1 ÷ 0, nothing)
end


text:
# ### log log

#md # ![FILE_NAME.png](images/FILE_NAME.png)

f(x) = x^-1.7
plot(f, 1e-3, 3, scale=:log10, title="log-log")


text:
# Helper functions
function d_subproblem{T}(s_fft::Array{Complex{T}, 2}, mu_s::T, lambda_s_fft::Array{Complex{T}, 2},
                         ZtZ_fft::Array{Complex{T}, 3}, Ztx_fft::Array{Complex{T}, 2})
  convLen = size(s_fft, 2);
  dictSize = size(ZtZ_fft, 1);
  # (Eq. 14)
  d_fft = zeros(s_fft);
  A = zeros(Complex128, dictSize, dictSize); b = zeros(Complex128, dictSize, 1);
  for i=1:convLen;
    # d_fft[:, i] = ( ZtZ_fft[:, :, i] + mu_s.*I ) \ ( Ztx_fft[:, i] + mu_s.*s_fft[:, i] - lambda_s_fft[:, i] );
    copy!(A, ZtZ_fft[:, :, i]);
    for j=1:dictSize; A[j, j] += mu_s; end
    copy!(b, s_fft[:, i]);
    scale!(b, mu_s);
    b += Ztx_fft[:, i];
    b -= lambda_s_fft[:, i];
    d_fft[:, i] = A\b;
  end
  return d_fft::Array{Complex{T}, 2};
end

function s_subproblem{T}(d::Array{T, 3}, mu_s::T, lambda_s::Array{T, 3}, mask::Array{Int, 1})
  # (Eq. 20)
  # s = d + 1/mu_s.*lambda_s;
  s = copy(lambda_s);
  scale!(s, 1/mu_s);
  s += d;
  s[mask] = 0;
  # (Eq. 18)
  normalize!(s);
  return s::Array{T, 3};
end

function fcsc_learning_multi{T}(data::LearningData{T}, optParams::LearningOpts{T})
# Alg. parameters:
#   mu_max: 1e5
#   tau: [1.01 1.5] larger means faster convergence but sacrifices primal feasibility,
#        smaller means slower convergence

  const mu_max = 1e5;
  const tau = 1.01;

  const ZtZ_fft    = data.ZtZ_fft;
  const Ztx_fft    = data.Ztx_fft;
  const convH      = data.convH;
  const convW      = data.convW;
  const filterSize = data.filterSize;
  const d          = optParams.Dinit;
        s          = optParams.s;
        lambda_s   = optParams.lambda_s;
        mu_s::T    = optParams.mu_s;
  const max_iter   = optParams.max_iter;
  const tol        = optParams.tol;
  const DEBUG      = optParams.DEBUG;
        debugData  = optParams.debugData;

  dictSize = size(ZtZ_fft, 1);
  convLen = convH*convW;

  assert(size(ZtZ_fft, 3)==convLen, "ZtZ_fft does not have the correct depth.")
  assert(size(Ztx_fft, 1)==dictSize && size(Ztx_fft, 2)==convLen,
    "Ztx_fft does not have the correct dimensions.")
  assert(isempty(d) || size(d, 3)==dictSize,
    "Number of filters in Dinit does not match the number of channels in Z.")

  mask = ones(T, convH, convW, dictSize);
  mask[1:filterSize, 1:filterSize, :] = 0;
  mask = find(mask);

  # Initialization
  if isempty(lambda_s);
    lambda_s = zeros(T, convH, convW, dictSize);
  end
  lambda_s_fft = reshape(fft2(lambda_s), convLen, dictSize).';
  if isempty(d);
    d = zeros(T, convH, convW, dictSize);
  end
  d_fft = reshape(fft2(d), convLen, dictSize).';
  if isempty(s);
    s = s_subproblem(d, mu_s, lambda_s, mask);
  end
  s_fft = reshape(fft2(s), convLen, dictSize).';

  if DEBUG;
    debugData["d"] = d;
    debugData["s"] = s;
    debugData["lambda_s"] = lambda_s;
    (prev_obj, rerr, spar) = compute_lagrangian(debugData);
  end

  iter = 0;
  converged = false;
  while ~converged;
    iter += 1;

    if max_iter>1; @printf("\b\b\b%03d", iter); end

    # Solve for d* in the frequency domain
    d_fft = d_subproblem(s_fft, mu_s, lambda_s_fft, ZtZ_fft, Ztx_fft);
    d = real(ifft2(reshape(d_fft.', convH, convW, dictSize)));
    if DEBUG;
      debugData["d"] = d;
      (obj, rerr, spar) = compute_lagrangian(debugData);
      is_mono = (prev_obj>=obj) ? 'O' : 'X';
      @printf("          D: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
      prev_obj = obj;
    end

    # Solve for s* in the spatial domain
    s = s_subproblem(d, mu_s, lambda_s, mask);
    if DEBUG;
      debugData["s"] = s;
      (obj, rerr, spar) = compute_lagrangian(debugData);
      is_mono = (prev_obj>=obj) ? 'O' : 'X';
      @printf("          S: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
      prev_obj = obj;
    end

    if max_iter==1; return (d, s); end # Don't do unnecessary work
    s_fft = reshape(fft2(s), convLen, dictSize).';

    # Update the Lagrangian multipliers
    for i=1:convLen*dictSize; lambda_s[i] += mu_s*( d[i]-s[i] ); end
    lambda_s_fft = reshape(fft2(lambda_s), convLen, dictSize).';
    if DEBUG;
      debugData["lambda_s"] = lambda_s;
      (obj, rerr, spar) = compute_lagrangian(debugData);
      is_mono = (prev_obj>=obj) ? 'O' : 'X';
      @printf("          L: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
      prev_obj = obj;
    end

    # Update penalty term
    mu_s = min(tau*mu_s, mu_max);
    if DEBUG;
      debugData["mu_s"] = mu_s;
      (obj, rerr, spar) = compute_lagrangian(debugData);
      is_mono = (prev_obj>=obj) ? 'O' : 'X';
      @printf("          M: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
      prev_obj = obj;
    end

    if iter>=max_iter || mean((d-s).^2)<=tol;
      converged = true;
    end
  end

  return (d, s);

end

max_stars_repo_name: bkong/FCSC


text:
# Reference: https://github.com/mdavezac/spglib.jl/blob/master/src/spglib.jl#L70
# This is an internal function, do not export!
function cchars2string(itr)
    vec = collect(Char, Iterators.filter(!iszero, itr))
    return String(vec)
end

# See https://github.com/spglib/spglib/blob/444e061/python/spglib/spglib.py#L115-L165
"""
    get_symmetry(cell::Cell, symprec=1e-5)

Return the symmetry operations of a `cell`.
"""
function get_symmetry(cell::Cell, symprec = 1e-5)
    max_size = length(cell.types) * 48
    rotation = Array{Cint,3}(undef, 3, 3, max_size)
    translation = Array{Cdouble,2}(undef, 3, max_size)
    if iszero(cell.magmoms)
        return get_symmetry!(rotation, translation, cell, symprec)
    else
        equivalent_atoms = zeros(length(cell.magmoms))
        primitive_lattice = zeros(Cdouble, 3, 3)
        if ndims(cell.magmoms) == 1
            spin_flips = zeros(length(rotation))
        else
            spin_flips = nothing
        end
        # TODO: unfinished!
    end
end

function get_symmetry!(
    rotation::AbstractArray,
    translation::AbstractMatrix,
    cell::Cell,
    symprec = 1e-5,
)
    if size(rotation, 3) != size(translation, 2)
        throw(DimensionMismatch("`rotation` & `translation` have different max size!"))
    end
    if !(size(rotation, 1) == size(rotation, 2) == size(translation, 1) == 3)
        throw(DimensionMismatch("`rotation` & `translation` don't have the right size!"))
    end
    @unpack lattice, positions, types = _expand_cell(cell)
    rotation = Base.cconvert(Array{Cint,3}, rotation)
    translation = Base.cconvert(Matrix{Cdouble}, translation)
    max_size = Base.cconvert(Cint, size(rotation, 3))
    num_atom = Base.cconvert(Cint, length(types))
    num_sym = ccall(
        (:spg_get_symmetry, libsymspg),
        Cint,
        (
            Ptr{Cint},
            Ptr{Float64},
            Cint,
            Ptr{Float64},
            Ptr{Float64},
            Ptr{Cint},
            Cint,
            Float64,
        ),
        rotation,
        translation,
        max_size,
        lattice,
        positions,
        types,
        num_atom,
        symprec,
    )
    if num_sym == 0
        throw(SpglibError("Symmetry operation search failed!"))
    end
    return rotation[:, :, 1:num_sym], translation[:, 1:num_sym]
end

function get_symmetry_with_collinear_spin!(
    rotation::AbstractArray{T,3},
    translation::AbstractMatrix,
    equivalent_atoms::AbstractVector,
    max_size::Integer,
    cell::Cell,
    symprec = 1e-5,
) where {T}
    @unpack lattice, positions, types, magmoms = _expand_cell(cell)
    rotation = Base.cconvert(Array{Cint,3}, rotation)
    translation = Base.cconvert(Matrix{Cdouble}, translation)
    equivalent_atoms = Base.cconvert(Vector{Cint}, equivalent_atoms)
    max_size = Base.cconvert(Cint, max_size)
    num_atom = Base.cconvert(Cint, length(types))
    num_sym = ccall(
        (:spg_get_symmetry_with_collinear_spin, libsymspg),
        Cint,
        (
            Ptr{Cint},
            Ptr{Cdouble},
            Ptr{Cint},
            Cint,
            Ptr{Cdouble},
            Ptr{Cdouble},
            Ptr{Cint},
            Ptr{Cdouble},
            Cint,
            Cdouble,
        ),
        rotation,
        translation,
        equivalent_atoms,
        max_size,
        lattice,
        positions,
        types,
        magmoms,
        num_atom,
        symprec,
    )
    if num_sym == 0
        throw(SpglibError("Symmetry operation search failed!"))
    end
    return num_sym
end
function get_symmetry_with_collinear_spin(cell::Cell, symprec = 1e-5)
    num_atom = length(cell.types)
    max_size = num_atom * 48
    rotation = Array{Cint,3}(undef, 3, 3, max_size)
    translation = Matrix{Cdouble}(undef, 3, max_size)
    equivalent_atoms = Vector{Cint}(undef, num_atom)
    num_sym = get_symmetry_with_collinear_spin!(
        rotation,
        translation,
        equivalent_atoms,
        max_size,
        cell,
        symprec,
    )
    return rotation[:, :, 1:num_sym], translation[:, 1:num_sym], equivalent_atoms
end

"""
    get_symmetry_from_database(hall_number)

Return the symmetry operations given a `hall_number`.

This function allows to directly access to the space group operations in the
`spglib` database. To specify the space group type with a specific choice,
`hall_number` is used.
"""
function get_symmetry_from_database(hall_number)
    rotation = Array{Cint,3}(undef, 3, 3, 192)
    translation = Array{Cdouble,2}(undef, 3, 192)
    return get_symmetry_from_database!(rotation, translation, hall_number)
end

function get_symmetry_from_database!(
    rotation::AbstractArray,
    translation::AbstractMatrix,
    hall_number,
)
    if !(size(rotation, 3) == size(translation, 2) == 192)
        throw(
            DimensionMismatch(
                "`rotation` & `translation` should have space for 192 symmetry operations!",
            ),
        )
    end
    if !(size(rotation, 1) == size(rotation, 2) == size(translation, 1) == 3)
        throw(DimensionMismatch("`rotation` & `translation` don't have the right size!"))
    end
    rotation = Base.cconvert(Array{Cint,3}, rotation)
    translation = Base.cconvert(Matrix{Cdouble}, translation)
    hall_number = Base.cconvert(Cint, hall_number)
    num_sym = ccall(
        (:spg_get_symmetry_from_database, libsymspg),
        Cint,
        (Ptr{Cint}, Ptr{Float64}, Cint),
        rotation,
        translation,
        hall_number,
    )
    if num_sym == 0
        throw(SpglibError("Symmetry operation search failed!"))
    end
    return rotation[:, :, 1:num_sym], translation[:, 1:num_sym]
end

"""
    get_hall_number_from_symmetry(rotation::AbstractArray{T,3}, translation::AbstractMatrix, num_operations::Integer, symprec=1e-5) where {T}

Obtain `hall_number` from the set of symmetry operations.

This is expected to work well for the set of symmetry operations whose
distortion is small. The aim of making this feature is to find space-group-type
for the set of symmetry operations given by the other source than spglib. Note
that the definition of `symprec` is different from usual one, but is given in the
fractional coordinates and so it should be small like `1e-5`.
"""
function get_hall_number_from_symmetry(
    rotation::AbstractArray{T,3},
    translation::AbstractMatrix,
    num_operations::Integer,
    symprec = 1e-5,
) where {T}
    rotation = Base.cconvert(Array{Cint,3}, rotation)
    translation = Base.cconvert(Matrix{Cdouble}, translation)
    num_operations = Base.cconvert(Cint, num_operations)
    return ccall(
        (:spg_get_hall_number_from_symmetry, libsymspg),
        Cint,
        (Ptr{Cint}, Ptr{Float64}, Cint, Float64),
        rotation,
        translation,
        num_operations,
        symprec,
    )
end

"""
    get_multiplicity(cell::Cell, symprec=1e-5)

Return the exact number of symmetry operations. An error is thrown when it fails.
"""
function get_multiplicity(cell::Cell, symprec = 1e-5)
    @unpack lattice, positions, types = _expand_cell(cell)
    num_atom = Base.cconvert(Cint, length(types))
    num_sym = ccall(
        (:spg_get_multiplicity, libsymspg),
        Cint,
        (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint, Cdouble),
        lattice,
        positions,
        types,
        num_atom,
        symprec,
    )
    if num_sym == 0
        throw(SpglibError("Symmetry operation search failed!"))
    end
    return num_sym
end

"""
    get_dataset(cell::Cell, symprec=1e-5)

Search symmetry operations of an input unit cell structure.
"""
function get_dataset(cell::Cell, symprec = 1e-5)
    @unpack lattice, positions, types = _expand_cell(cell)
    num_atom = Base.cconvert(Cint, length(types))
    ptr = ccall(
        (:spg_get_dataset, libsymspg),
        Ptr{SpglibDataset},
        (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint, Cdouble),
        lattice,
        positions,
        types,
        num_atom,
        symprec,
    )
    if ptr == C_NULL
        return nothing
    else
        raw = unsafe_load(ptr)
        return convert(Dataset, raw)
    end
end

"""
    get_dataset_with_hall_number(cell::Cell, hall_number::Integer, symprec=1e-5)

Search symmetry operations of an input unit cell structure, using a given Hall number.
"""
function get_dataset_with_hall_number(cell::Cell, hall_number::Integer, symprec = 1e-5)
    @unpack lattice, positions, types = _expand_cell(cell)
    num_atom = Base.cconvert(Cint, length(types))
    hall_number = Base.cconvert(Cint, hall_number)
    ptr = ccall(
        (:spg_get_dataset_with_hall_number, libsymspg),
        Ptr{SpglibDataset},
        (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint, Cint, Cdouble),
        lattice,
        positions,
        types,
        num_atom,
        hall_number,
        symprec,
    )
    if ptr == C_NULL
        return nothing
    else
        raw = unsafe_load(ptr)
        return convert(Dataset, raw)
    end
end

"""
    get_spacegroup_type(hall_number::Integer)

Translate Hall number to space group type information.
"""
function get_spacegroup_type(hall_number::Integer)
    spgtype = ccall(
        (:spg_get_spacegroup_type, libsymspg),
        SpglibSpacegroupType,
        (Cint,),
        hall_number,
    )
    return convert(SpacegroupType, spgtype)
end

"""
    get_spacegroup_number(cell::Cell, symprec=1e-5)

Get the spacegroup number of a `cell`.
"""
function get_spacegroup_number(cell::Cell, symprec = 1e-5)
    dataset = get_dataset(cell, symprec)
    return dataset.spacegroup_number
end
"""
    get_spacegroup_type(cell::Cell, symprec=1e-5)

Get `SpacegroupType` from `cell`.
"""
function get_spacegroup_type(cell::Cell, symprec = 1e-5)  # See https://github.com/spglib/spglib/blob/444e061/python/spglib/spglib.py#L307-L324
    dataset = get_dataset(cell, symprec)
    return get_spacegroup_type(dataset.hall_number)
end

"""
    get_international(cell::Cell, symprec=1e-5)

Return the space group type in Hermann–Mauguin (international) notation.
"""
function get_international(cell::Cell, symprec = 1e-5)
    @unpack lattice, positions, types = _expand_cell(cell)
    symbol = Vector{Cchar}(undef, 11)
    exitcode = ccall(
        (:spg_get_international, libsymspg),
        Cint,
        (Ptr{Cchar}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint, Cdouble),
        symbol,
        lattice,
        positions,
        types,
        length(types),
        symprec,
    )
    if exitcode == 0
        throw(SpglibError("Spacegroup search failed!"))
    end
    return cchars2string(symbol)
end

"""
    get_schoenflies(cell::Cell, symprec=1e-5)

Return the space group type in Schoenflies notation.
"""
function get_schoenflies(cell::Cell, symprec = 1e-5)
    @unpack lattice, positions, types = _expand_cell(cell)
    symbol = Vector{Cchar}(undef, 7)
    exitcode = ccall(
        (:spg_get_schoenflies, libsymspg),
        Cint,
        (Ptr{Cchar}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint, Cdouble),
        symbol,
        lattice,
        positions,
        types,
        length(types),
        symprec,
    )
    if exitcode == 0
        throw(SpglibError("Spacegroup search failed!"))
    end
    return cchars2string(symbol)
end

function Base.convert(::Type{SpacegroupType}, spgtype::SpglibSpacegroupType)
    values = map(fieldnames(SpacegroupType)) do name
        value = getfield(spgtype, name)
        if value isa Cint
            value
        elseif value isa NTuple{N,Cchar} where {N}
            cchars2string(value)
        else  # This should never happen!
            @assert false "unexpected field type $(typeof(value))!"  # See https://discourse.julialang.org/t/please-stop-using-error-and-errorexception-in-packages-and-base/12096
        end
    end
    return SpacegroupType(values...)
end


text:
@fastmath function middleloop!( iz::Int, rb_tar_idx::Int, rb_src_idx::Int,
                                pssrb::PotentialSimulationSetupRB{T, Cartesian, 3, Array{T, 3}},
                                update_even_points::Val{even_points},
                                depletion_handling::Val{depletion_handling_enabled},
                                is_weighting_potential::Val{_is_weighting_potential},
                                only2d::Val{only_2d},
                                idx3iseven::Val{idx3_is_even})::Nothing where {T, even_points, depletion_handling_enabled, _is_weighting_potential, only_2d, idx3_is_even}
    @inbounds begin
        inz::Int = iz - 1

        pwwzr::T        = pssrb.geom_weights[3][1, inz]
        pwwzl::T        = pssrb.geom_weights[3][2, inz]
        pwΔmpz::T       = pssrb.geom_weights[3][3, inz]
        Δz_ext_inv_r::T = pssrb.geom_weights[3][4, inz + 1]
        Δz_ext_inv_l::T = pssrb.geom_weights[3][4, inz]

        line_weights::Array{T, 2} = Array{T, 2}(undef, size(pssrb.potential, 1) - 2, 6)
        # Even though this causes some allocations it is faster than using a predefined array, e.g. stored in pssrb
        # Especially when using multiple threads

        #=
            Splitting this into two loops over even and uneven iφ
            seems to increase the performance quite a bit
            by hard-coding oscillating type:
                zyi_is_even_t::Union{Val{true}, Val{false}} = Val(iseven(iz + iy) ? true : false)
        =#
        for iy in 2:2:(size(pssrb.potential, 2) - 1)
            iny::Int = iy - 1

            pwwyr::T  = pssrb.geom_weights[2][1, iny]
            pwwyl::T  = pssrb.geom_weights[2][2, iny]
            pwΔmpy::T = pssrb.geom_weights[2][3, iny]
            pwΔmpy_pwΔmpz::T = pwΔmpy * pwΔmpz
            Δy_ext_inv_r_pwΔmpz::T  = pssrb.geom_weights[2][4, iny + 1] * pwΔmpz
            Δy_ext_inv_l_pwΔmpz::T  = pssrb.geom_weights[2][4, iny]     * pwΔmpz
            Δz_ext_inv_r_pwΔmpy::T = Δz_ext_inv_r * pwΔmpy
            Δz_ext_inv_l_pwΔmpy::T = Δz_ext_inv_l * pwΔmpy

            pwwyr_pwwzr::T = pwwyr * pwwzr
            pwwyr_pwwzl::T = pwwyr * pwwzl
            pwwyl_pwwzr::T = pwwyl * pwwzr
            pwwyl_pwwzl::T = pwwyl * pwwzl

            load_weights_for_innerloop!(line_weights, pssrb, iy, iny, iz, inz,
                update_even_points, Val(idx3_is_even ? true : false),
                pwwzr, pwwzl, pwwyr, pwwyl,
                pwwyr_pwwzr, pwwyr_pwwzl, pwwyl_pwwzr, pwwyl_pwwzl,
                pwΔmpy_pwΔmpz,
                Δy_ext_inv_r_pwΔmpz, Δy_ext_inv_l_pwΔmpz,
                Δz_ext_inv_r_pwΔmpy, Δz_ext_inv_l_pwΔmpy,
            )

            innerloop!(line_weights, pssrb, iy, iny, iz, inz, rb_tar_idx, rb_src_idx,
                update_even_points, Val(idx3_is_even ? true : false),
                depletion_handling, is_weighting_potential, only2d)
        end
        for iy in 3:2:(size(pssrb.potential, 2) - 1)
            iny::Int = iy - 1

            pwwyr::T  = pssrb.geom_weights[2][1, iny]
            pwwyl::T  = pssrb.geom_weights[2][2, iny]
            pwΔmpy::T = pssrb.geom_weights[2][3, iny]
            pwΔmpy_pwΔmpz::T = pwΔmpy * pwΔmpz
            Δy_ext_inv_r_pwΔmpz::T  = pssrb.geom_weights[2][4, iny + 1] * pwΔmpz
            Δy_ext_inv_l_pwΔmpz::T  = pssrb.geom_weights[2][4, iny]     * pwΔmpz
            Δz_ext_inv_r_pwΔmpy::T = Δz_ext_inv_r * pwΔmpy
            Δz_ext_inv_l_pwΔmpy::T = Δz_ext_inv_l * pwΔmpy

            pwwyr_pwwzr::T = pwwyr * pwwzr
            pwwyr_pwwzl::T = pwwyr * pwwzl
            pwwyl_pwwzr::T = pwwyl * pwwzr
            pwwyl_pwwzl::T = pwwyl * pwwzl

            load_weights_for_innerloop!(line_weights, pssrb, iy, iny, iz, inz,
                update_even_points, Val(idx3_is_even ? false : true),
                pwwzr, pwwzl, pwwyr, pwwyl,
                pwwyr_pwwzr, pwwyr_pwwzl, pwwyl_pwwzr, pwwyl_pwwzl,
                pwΔmpy_pwΔmpz,
                Δy_ext_inv_r_pwΔmpz, Δy_ext_inv_l_pwΔmpz,
                Δz_ext_inv_r_pwΔmpy, Δz_ext_inv_l_pwΔmpy,
            )

            innerloop!(line_weights, pssrb, iy, iny, iz, inz, rb_tar_idx, rb_src_idx,
                update_even_points, Val(idx3_is_even ? false : true),
                depletion_handling, is_weighting_potential, only2d)
        end
    end
end

function load_weights_for_innerloop!(line_weights, pssrb::PotentialSimulationSetupRB{T, Cartesian, 3, Array{T, 3}},
        iy, iny, iz, inz,
        update_even_points, zyi_is_even_t::Val{zyi_is_even},
        pwwzr, pwwzl, pwwyr, pwwyl,
        pwwyr_pwwzr, pwwyr_pwwzl, pwwyl_pwwzr, pwwyl_pwwzl,
        pwΔmpy_pwΔmpz,
        Δy_ext_inv_r_pwΔmpz, Δy_ext_inv_l_pwΔmpz,
        Δz_ext_inv_r_pwΔmpy, Δz_ext_inv_l_pwΔmpy
        ) where {T, zyi_is_even}
    @fastmath @inbounds @simd ivdep for ix in 2:(size(pssrb.potential, 1) - 1)
        inx::Int = nidx(ix, update_even_points, zyi_is_even_t)::Int

        pwwxr::T        = pssrb.geom_weights[1][1, inx]
        pwwxl::T        = pssrb.geom_weights[1][2, inx]
        pwΔmpx::T       = pssrb.geom_weights[1][3, inx]
        Δx_ext_inv_r::T = pssrb.geom_weights[1][4, inx + 1]
        Δx_ext_inv_l::T = pssrb.geom_weights[1][4, inx]

        pwwxr_pwwzr::T = pwwxr * pwwzr
        pwwxl_pwwzr::T = pwwxl * pwwzr
        pwwxr_pwwzl::T = pwwxr * pwwzl
        pwwxl_pwwzl::T = pwwxl * pwwzl
        pwwxl_pwwyr::T = pwwxl * pwwyr
        pwwxr_pwwyr::T = pwwxr * pwwyr
        pwwxl_pwwyl::T = pwwxl * pwwyl
        pwwxr_pwwyl::T = pwwxr * pwwyl

        ϵ_rrr::T = pssrb.ϵ_r[ inx + 1,  iy, iz]
        ϵ_rlr::T = pssrb.ϵ_r[ inx + 1, iny, iz]
        ϵ_rrl::T = pssrb.ϵ_r[ inx + 1,  iy, inz ]
        ϵ_rll::T = pssrb.ϵ_r[ inx + 1, iny, inz ]
        ϵ_lrr::T = pssrb.ϵ_r[ inx,  iy,  iz ]
        ϵ_llr::T = pssrb.ϵ_r[ inx, iny,  iz ]
        ϵ_lrl::T = pssrb.ϵ_r[ inx,  iy, inz ]
        ϵ_lll::T = pssrb.ϵ_r[ inx, iny, inz ]

        # right weight in r: wrr
        wxr::T = ϵ_rrr * pwwyr_pwwzr
        wxr    = muladd(ϵ_rlr, pwwyl_pwwzr, wxr)
        wxr    = muladd(ϵ_rrl, pwwyr_pwwzl, wxr)
        wxr    = muladd(ϵ_rll, pwwyl_pwwzl, wxr)
        # left weight in r: wrr
        wxl::T = ϵ_lrr * pwwyr_pwwzr
        wxl    = muladd(ϵ_llr, pwwyl_pwwzr, wxl)
        wxl    = muladd(ϵ_lrl, pwwyr_pwwzl, wxl)
        wxl    = muladd(ϵ_lll, pwwyl_pwwzl, wxl)
        # right weight in φ: wφr
        wyr::T = ϵ_lrr * pwwxl_pwwzr
        wyr    = muladd(ϵ_rrr, pwwxr_pwwzr, wyr)
        wyr    = muladd(ϵ_lrl, pwwxl_pwwzl, wyr)
        wyr    = muladd(ϵ_rrl, pwwxr_pwwzl, wyr)
        # left weight in φ: wφl
        wyl::T = ϵ_llr * pwwxl_pwwzr
        wyl    = muladd(ϵ_rlr, pwwxr_pwwzr, wyl)
        wyl    = muladd(ϵ_lll, pwwxl_pwwzl, wyl)
        wyl    = muladd(ϵ_rll, pwwxr_pwwzl, wyl)
        # right weight in z: wzr
        wzr::T = ϵ_rrr * pwwxr_pwwyr
        wzr    = muladd(ϵ_rlr, pwwxr_pwwyl, wzr)
        wzr    = muladd(ϵ_lrr, pwwxl_pwwyr, wzr)
        wzr    = muladd(ϵ_llr, pwwxl_pwwyl, wzr)
        # left weight in z: wzr
        wzl::T = ϵ_rrl * pwwxr_pwwyr
        wzl    = muladd(ϵ_rll, pwwxr_pwwyl, wzl)
        wzl    = muladd(ϵ_lrl, pwwxl_pwwyr, wzl)
        wzl    = muladd(ϵ_lll, pwwxl_pwwyl, wzl)

        wxr *= pwΔmpy_pwΔmpz * Δx_ext_inv_r
        wxl *= pwΔmpy_pwΔmpz * Δx_ext_inv_l
        wyr *= Δy_ext_inv_r_pwΔmpz * pwΔmpx
        wyl *= Δy_ext_inv_l_pwΔmpz * pwΔmpx
        wzr *= Δz_ext_inv_r_pwΔmpy * pwΔmpx
        wzl *= Δz_ext_inv_l_pwΔmpy * pwΔmpx

        line_weights[ix-1, 1] = wxr
        line_weights[ix-1, 2] = wxl
        line_weights[ix-1, 3] = wyr
        line_weights[ix-1, 4] = wyl
        line_weights[ix-1, 5] = wzr
        line_weights[ix-1, 6] = wzl
    end
end

function innerloop!(line_weights, pssrb::PotentialSimulationSetupRB{T, Cartesian, 3, Array{T, 3}},
        iy, iny, iz, inz, rb_tar_idx, rb_src_idx,
        update_even_points, zyi_is_even_t::Val{zyi_is_even},
        depletion_handling::Val{depletion_handling_enabled},
        is_weighting_potential::Val{_is_weighting_potential},
        only2d::Val{only_2d}) where {T, depletion_handling_enabled, _is_weighting_potential, only_2d, zyi_is_even}
    @fastmath @inbounds @simd ivdep for ix in 2:(size(pssrb.potential, 1) - 1)
    # for ix in 2:(size(pssrb.potential, 1) - 1)
        ixr::Int = get_rbidx_right_neighbour(ix, update_even_points, zyi_is_even_t)

        wxr = line_weights[ix-1, 1]
        wxl = line_weights[ix-1, 2]
        wyr = line_weights[ix-1, 3]
        wyl = line_weights[ix-1, 4]
        wzr = line_weights[ix-1, 5]
        wzl = line_weights[ix-1, 6]

        vxr::T = pssrb.potential[    ixr,     iy,     iz, rb_src_idx]
        vxl::T = pssrb.potential[ixr - 1,     iy,     iz, rb_src_idx]
        vyr::T = pssrb.potential[     ix, iy + 1,     iz, rb_src_idx]
        vyl::T = pssrb.potential[     ix,    iny,     iz, rb_src_idx]
        vzr::T = pssrb.potential[     ix,     iy, iz + 1, rb_src_idx]
        vzl::T = pssrb.potential[     ix,     iy,    inz, rb_src_idx]

        new_potential::T = _is_weighting_potential ? zero(T) : (pssrb.q_eff_imp[ix, iy, iz, rb_tar_idx] + pssrb.q_eff_fix[ix, iy, iz, rb_tar_idx])
        new_potential = muladd( wxr, vxr, new_potential)
        new_potential = muladd( wxl, vxl, new_potential)
        new_potential = muladd( wyr, vyr, new_potential)
        new_potential = muladd( wyl, vyl, new_potential)
        new_potential = muladd( wzr, vzr, new_potential)
        new_potential = muladd( wzl, vzl, new_potential)

        new_potential *= pssrb.volume_weights[ix, iy, iz, rb_tar_idx]

        old_potential::T = pssrb.potential[ix, iy, iz, rb_tar_idx]

        new_potential -= old_potential
        new_potential = muladd(new_potential, pssrb.sor_const[1], old_potential)

        if depletion_handling_enabled
            vmin::T = min(vxr, vxl, vyr, vyl, vzr, vzl)
            vmax::T = max(vxr, vxl, vyr, vyl, vzr, vzl)

            new_point_type = pssrb.point_types[ix, iy, iz, rb_tar_idx]
            if new_potential < vmin || new_potential > vmax
                new_potential -= pssrb.q_eff_imp[ix, iy, iz, rb_tar_idx] * pssrb.volume_weights[ix, iy, iz, rb_tar_idx] * pssrb.sor_const[1]
                if (pssrb.point_types[ix, iy, iz, rb_tar_idx] & undepleted_bit == 0 &&
                    pssrb.point_types[ix, iy, iz, rb_tar_idx] & pn_junction_bit > 0)
                    new_point_type += undepleted_bit
                end # mark this point as undepleted
            elseif pssrb.point_types[ix, iy, iz, rb_tar_idx] & undepleted_bit > 0
                new_point_type -= undepleted_bit
            end
            pssrb.point_types[ix, iy, iz, rb_tar_idx] = new_point_type
        end

        pssrb.potential[ix, iy, iz, rb_tar_idx]::T = ifelse(pssrb.point_types[ix, iy, iz, rb_tar_idx] & update_bit > 0, new_potential, old_potential)
    end
end


text:
using Test

@timed_testset "utils/" begin
    include("utils/import_datasets.jl")
    include("utils/import_example_nets.jl")
end


text:
## indices
import ..ParallelStencil: INDICES, WITHIN_DOC
ix, iy, iz = INDICES[1], INDICES[2], INDICES[3]
ixi, iyi, izi = :($ix+1), :($iy+1), :($iz+1);
## 2nd-order accuracy for first-order derivatives
ix_1=:($ix+1);
iy_1=:($iy+1);
iz_1=:($iz+1);

macro   d_xi_1(A::Symbol)  esc(:( $A[$ix_1+1,$iyi ,$izi ] - $A[$ix_1  ,$iyi ,$izi ] )) end
macro   d_yi_1(A::Symbol)  esc(:( $A[$ixi ,$iy_1+1,$izi ] - $A[$ixi ,$iy_1  ,$izi ] )) end
macro   d_zi_1(A::Symbol)  esc(:( $A[$ixi ,$iyi ,$iz_1+1] - $A[$ixi ,$iyi ,$iz_1  ] )) end
## 8th-order accuracy for first-order derivatives
ix_8=:($ix+3);
iy_8=:($iy+3);
iz_8=:($iz+3);

macro   d_xi_8(A::Symbol)
    esc(:(1225/1024*($A[$ix_8+1,$iyi,$izi]-$A[$ix_8,$iyi,$izi])-
    245/3072*($A[$ix_8+2,$iyi,$izi]-$A[$ix_8-1,$iyi,$izi])+
    49/5120*($A[$ix_8+3,$iyi,$izi]-$A[$ix_8-2,$iyi,$izi])-
    5/7186*($A[$ix_8+4,$iyi,$izi]-$A[$ix_8-3,$iyi,$izi])));
end

macro   d_yi_8(A::Symbol)
    esc(:(1225/1024*($A[$ixi,$iy_8+1,$izi]-$A[$ixi,$iy_8,$izi])-
    245/3072*($A[$ixi,$iy_8+2,$izi]-$A[$ixi,$iy_8-1,$izi])+
    49/5120*($A[$ixi,$iy_8+3,$izi]-$A[$ixi,$iy_8-2,$izi])-
    5/7168*($A[$ixi,$iy_8+4,$izi]-$A[$ixi,$iy_8-3,$izi])));
end
macro   d_zi_8(A::Symbol)
    esc(:(1225/1024*($A[$ixi,$iyi,$iz_8+1]-$A[$ixi,$iyi,$iz_8])-
    245/3072*($A[$ixi,$iyi,$iz_8+2]-$A[$ixi,$iyi,$iz_8-1])+
    49/5120*($A[$ixi,$iyi,$iz_8+3]-$A[$ixi,$iyi,$iz_8-2])-
    5/7168*($A[$ixi,$iyi,$iz_8+4]-$A[$ixi,$iyi,$iz_8-3])));
end
## 12th-order accuracy for first-order derivatives
ix_12=:($ix+5);
iy_12=:($iy+5);
iz_12=:($iz+5);

macro   d_xi_12(A::Symbol)
    esc(:(160083/131072*($A[$ix_12+1,$iyi,$izi]-$A[$ix_12,$iyi,$izi])-
    12705/131072*($A[$ix_12+2,$iyi,$izi]-$A[$ix_12-1,$iyi,$izi])+
    22869/1310720*($A[$ix_12+3,$iyi,$izi]-$A[$ix_12-2,$iyi,$izi])-
    5445/1835008*($A[$ix_12+4,$iyi,$izi]-$A[$ix_12-3,$iyi,$izi])+
    847/2359296*($A[$ix_12+5,$iyi,$izi]-$A[$ix_12-4,$iyi,$izi])-
    63/2883584*($A[$ix_12+6,$iyi,$izi]-$A[$ix_12-5,$iyi,$izi])
    ));
end

macro   d_yi_12(A::Symbol)
    esc(:(160083/131072*($A[$ixi,$iy_12+1,$izi]-$A[$ixi,$iy_12,$izi])-
    12705/131072*($A[$ixi,$iy_12+2,$izi]-$A[$ixi,$iy_12-1,$izi])+
    22869/1310720*($A[$ixi,$iy_12+3,$izi]-$A[$ixi,$iy_12-2,$izi])-
    5445/1835008*($A[$ixi,$iy_12+4,$izi]-$A[$ixi,$iy_12-3,$izi])+
    847/2359296*($A[$ixi,$iy_12+5,$izi]-$A[$ixi,$iy_12-4,$izi])-
    63/2883584*($A[$ixi,$iy_12+6,$izi]-$A[$ixi,$iy_12-5,$izi])
    ));
end
macro   d_zi_12(A::Symbol)
    esc(:(160083/131072*($A[$ixi,$iyi,$iz_12+1]-$A[$ixi,$iyi,$iz_12])-
    12705/131072*($A[$ixi,$iyi,$iz_12+2]-$A[$ixi,$iyi,$iz_12-1])+
    22869/1310720*($A[$ixi,$iyi,$iz_12+3]-$A[$ixi,$iyi,$iz_12-2])-
    5445/1835008*($A[$ixi,$iyi,$iz_12+4]-$A[$ixi,$iyi,$iz_12-3])+
    847/2359296*($A[$ixi,$iyi,$iz_12+5]-$A[$ixi,$iyi,$iz_12-4])-
    63/2883584*($A[$ixi,$iyi,$iz_12+6]-$A[$ixi,$iyi,$iz_12-5])
    ));
end


text:
mutable struct BasicContUnvParameterNState{N<:Real} <: ParameterNState{Continuous, Univariate}
  value::Vector{N}
  loglikelihood::Vector{N}
  logprior::Vector{N}
  logtarget::Vector{N}
  gradloglikelihood::Vector{N}
  gradlogprior::Vector{N}
  gradlogtarget::Vector{N}
  tensorloglikelihood::Vector{N}
  tensorlogprior::Vector{N}
  tensorlogtarget::Vector{N}
  dtensorloglikelihood::Vector{N}
  dtensorlogprior::Vector{N}
  dtensorlogtarget::Vector{N}
  diagnosticvalues::Matrix
  monitor::Vector{Bool}
  n::Integer
  diagnostickeys::Vector{Symbol}

  function BasicContUnvParameterNState{N}(
    n::Integer,
    monitor::Vector{Bool}=[true; fill(false, 12)],
    diagnostickeys::Vector{Symbol}=Symbol[],
    ::Type{N}=Float64,
    diagnosticvalues::Matrix=Array{Any}(length(diagnostickeys), isempty(diagnostickeys) ? 0 : n)
  ) where N<:Real
    instance = new()

    l = Array{Integer}(13)
    for i in 1:13
      l[i] = (monitor[i] == false ? zero(Integer) : n)
    end

    fnames = fieldnames(BasicContUnvParameterNState)
    for i in 1:13
      setfield!(instance, fnames[i], Array{N}(l[i]))
    end

    instance.diagnosticvalues = diagnosticvalues
    instance.monitor = monitor
    instance.n = n
    instance.diagnostickeys = diagnostickeys

    instance
  end
end

BasicContUnvParameterNState(
  n::Integer,
  monitor::Vector{Bool}=[true; fill(false, 12)],
  diagnostickeys::Vector{Symbol}=Symbol[],
  ::Type{N}=Float64,
  diagnosticvalues::Matrix=Array{Any}(length(diagnostickeys), isempty(diagnostickeys) ? 0 : n)
) where N<:Real =
  BasicContUnvParameterNState{N}(n, monitor, diagnostickeys, N, diagnosticvalues)

function BasicContUnvParameterNState(
  n::Integer,
  monitor::Vector{Symbol},
  diagnostickeys::Vector{Symbol}=Symbol[],
  ::Type{N}=Float64,
  diagnosticvalues::Matrix=Array{Any}(length(diagnostickeys), isempty(diagnostickeys) ? 0 : n)
) where N<:Real
  fnames = fieldnames(BasicContUnvParameterNState)
  BasicContUnvParameterNState(n, [fnames[i] in monitor ? true : false for i in 1:13], diagnostickeys, N, diagnosticvalues)
end

const ContUnvMarkovChain = BasicContUnvParameterNState

function copy!(nstate::BasicContUnvParameterNState, state::BasicContUnvParameterState, i::Integer)
  fnames = fieldnames(BasicContUnvParameterNState)

  for j in 1:13
    if nstate.monitor[j]
      getfield(nstate, fnames[j])[i] = getfield(state, fnames[j])
    end
  end

  if !isempty(nstate.diagnosticvalues)
    nstate.diagnosticvalues[:, i] = state.diagnosticvalues
  end
end

eltype{N<:Real}(::Type{BasicContUnvParameterNState{N}}) = N
eltype{N<:Real}(::BasicContUnvParameterNState{N}) = N

==(z::S, w::S) where {S<:BasicContUnvParameterNState} = reduce(&, [getfield(z, n) == getfield(w, n) for n in fieldnames(S)[1:16]])

isequal(z::S, w::S) where {S<:BasicContUnvParameterNState} =
  reduce(&, [isequal(getfield(z, n), getfield(w, n)) for n in fieldnames(S)[1:16]])

function show(io::IO, nstate::BasicContUnvParameterNState{N}) where N<:Real
  fnames = fieldnames(BasicContUnvParameterNState)
  fbool = map(n -> !isempty(getfield(nstate, n)), fnames[1:13])
  indentation = "  "

  println(io, "BasicContUnvParameterNState:")

  println(io, indentation*"eltype: $(eltype(nstate))")
  println(io, indentation*"number of states = $(nstate.n)")

  print(io, indentation*"monitored components:")
  if !any(fbool)
    println(io, " none")
  else
    print(io, "\n")
    for i in 1:13
      if fbool[i]
        println(io, string(indentation^2, fnames[i]))
      end
    end
  end

  print(io, indentation*"diagnostics:")
  if isempty(nstate.diagnostickeys)
    print(io, " none")
  else
    for k in nstate.diagnostickeys
      print(io, "\n")
      print(io, string(indentation^2, k))
    end
  end
end


text:
__precompile__()

module IsoRank

using LinearMaps
using DataStructures

export isorank, greedyalign

"""
    kronlm([T], A, B)

Kronecker product of `A` and `B`, stored as a linear operator (from
[LinearMaps.jl](https://github.com/Jutho/LinearMaps.jl.git)) so that
you don't have to create the actual matrix, i.e. `kronlm(A,B)*x == kron(A,B)*x`.
This is much faster than
directly creating the matrix: `O(|V|^2+|V||E|)` instead of `O(|V|^2 + |E|^2)` for each
step of the power iteration where `|V|` and `|E|` are the number of nodes and edges
in the graphs.

# Arguments
- `A,B` : linear operators with multiply and transpose operations
- `T` : element type of the resulting linear operator
"""
function kronlm(::Type{T},A,B) where {T}
    f = (y,x) -> begin
        y[:] = vec(B * reshape(x, size(B,2), size(A,2)) * A')
        y end
    fc = (y,x) -> begin
        y[:] = vec(B' * reshape(x, size(B,1), size(A,1)) * A)
        y end
    m = size(A,1) * size(B,1)
    n = size(A,2) * size(B,2)
    LinearMap{T}(f, fc, m, n)
end
kronlm(A,B) = kronlm(promote_type(eltype(A),eltype(B)),A,B)

"""
    powermethod!(A, x, Ax=similar(x);
                 <keyword arguments>) -> radius, x, [log/history]

Performs power method in order to find the dominant eigenvector
of the linear operator A. Eigenvector is normalized w.r.t. L_1 norm.
Modifies initial eigenvector estimate x.

# Arguments
- `A` : linear operator
- `x` : initial estimate of the eigenvector, not necessarily normalized
- `Ax` : scratch space for the iteration process

# Keyword arguments
- `maxiter` : maximum # of iterations
- `tol=eps(Float64) * size(A,2)` : error tolerance in L_1 norm
- `log=true,verbose=true` : logging and printing
"""
function powermethod!(A, x, Ax=similar(x);
                      maxiter=15, tol=eps(Float64) * size(A,2),
                      log=true, verbose=true)
    T = eltype(A)
    x ./= norm(x,1)
    verbose && println("Running power method, maxiter = $maxiter, tol = $tol")
    history = Tuple{Int,T,T}[]
    iter = 0
    radius = zero(T)
    err = Inf
    while iter <= maxiter
        A_mul_B!(Ax, A, x)
        radius = norm(Ax,1)
        Ax ./= radius # want |x|_1 = 1
        err = norm(Ax-x,1)
        verbose && @show iter,err,radius
        log && push!(history,(iter,err,radius))
        copy!(x, Ax)
        err < tol && break
        iter += 1
    end
    isconverged = err < tol
    if log radius,x,history,isconverged else radius,x end
end

"""
    isorank(G1::SparseMatrixCSC, G2::SparseMatrixCSC,
            [alpha::Real=0.85, B::AbstractMatrix=ones(Float64,size(G1,1),size(G2,1))];
            <keyword arguments>) -> R [, res, L]

Creates the IsoRank matrix, which contains topological similarities of
all nodes across two networks (See Rohit Singh, Jinbo Xu, and Bonnie
Berger. (2008) Global alignment of multiple protein interaction
networks with application to functional orthology detection,
Proc. Natl. Acad. Sci. USA, 105:12763-12768.). That is, finds the
PageRank values of the modified adjacency matrix of the product graph
of G1 and G2.  `B`, containing prior node similarities, acts as the
personalization vector, allowing you to incorporate external
information. If you don't have node similarities, you can still create
a good IsoRank matrix by damping like PageRank does.

# Arguments
- `G1,G2` : two adjacency matrices of undirected networks
- `B` : matrix of node similarities between `G1` and `G2`, not necessarily normalized
- `alpha`: weight between edge and node conservation

# Keyword arguments
- `details=false` : If true, returns (R,res,L) where R is the IsoRank
    matrix of node similarities, res is the power method detailed
    results structure, L is the linear operator that the power method
    finds the eigenvector of; if false, returns R
- See [`powermethod!`](@ref) for other keyword arguments
"""
function isorank(G1::AbstractMatrix, G2::AbstractMatrix,
                 alpha::Real=0.85, B::AbstractMatrix=ones(Float64,size(G1,1),size(G2,1));
                 details=false, args...)
    A = kronlm(Float64,G2,G1)
    res = pagerank(A, alpha, vec(B); details=details, args...)
    if details
        R = reshape(res[1], size(G1,1), size(G2,1))
        R, res[2], res[3], A
    else
        R = reshape(res, size(G1,1), size(G2,1))
        R
    end
end

"""
     pagerank(A, alpha=0.85, p = fill(1/size(A,2),size(A,2)),
              x=copy(p), Ax=similar(x);
              <keyword args>) -> x [, res, L]

Creates PageRank vector.

# Arguments
- `A` : Adjacency matrix of the graph. A[u,v] = 1 if u -> v
- `alpha` : Damping.
- `p` : Initial probability vector, or personalization vector, not necessarily normalized.

# Keyword arguments
- `details=false` : If true, returns (x,res,L) where x is the PageRank
  vector, res is the power method detailed results structure, L is the linear
  operator that the power method finds the eigenvector of; if false,
  returns x.
- See [`powermethod!`](@ref) for other keyword arguments
"""
function pagerank(A, alpha=0.85, p = fill(1/size(A,2),size(A,2));
                  details=false, args...)
    S = 1.0 ./ (A * ones(Float64,size(A,2)))
    D = find(isinf,S) # S[i] = Inf if node i has no outlinks
    S[D] = 0.0
    p = p./norm(p,1)
    L = LinearMap{Float64}((y,x) -> begin
                           At_mul_B!(y, A, S .* x)
                           y .= alpha .* y .+ (alpha * sum(x[D]) + 1.0 - alpha) .* p
                           y
                           end, size(A,1), size(A,2))
    x = copy(p)
    res = powermethod!(L, x; args...)
    if details x, res, L else x end
end

"""
    greedyalign(R::AbstractMatrix,seeds=Vector{Tuple{Int,Int}}();
                maxiter=size(R,1)-length(seeds))

Given matrix `R`, find a alignment using the greedy method.

# Arguments
- `R` : m x n matrix
- `seeds` : Initial seed node pairs
- `maxiter` : Stop after aligning `maxiter` node pairs
"""
function greedyalign(R::AbstractMatrix,seeds=Vector{Tuple{Int,Int}}();
                     maxiter=size(R,1)-length(seeds))
    f = zeros(Int,size(R,1))
    n1 = size(R,1); n2 = size(R,2)
    kv = (((i,j),-R[i,j]) for i=1:n1, j=1:n2)
    Q = PriorityQueue{Tuple{Int,Int},eltype(R)}(kv)
    L1 = Set{Int}()
    L2 = Set{Int}() # nodes already aligned
    for k = 1:length(seeds)
        i,j = seeds[k]
        f[i] = j
        push!(L1,i)
        push!(L2,j)
        for ip = setdiff(1:n1,L1); dequeue!(Q,(ip,j)); end
        for jp = setdiff(1:n2,L2); dequeue!(Q,(i,jp)); end
    end
    println("Building priority queue")
    iter = 1
    while iter <= maxiter && !isempty(Q)
        i,j = dequeue!(Q)
        f[i] = j
        push!(L1,i)
        push!(L2,j)
        for ip = setdiff(1:n1,L1); dequeue!(Q,(ip,j)); end
        for jp = setdiff(1:n2,L2); dequeue!(Q,(i,jp)); end
        print("\rIteration $iter/$maxiter")
        iter += 1
    end
    println()
    f
end

end


text:
using PlanningDomains, Test
using PlanningDomains: clear_cache!, clear_all_caches!

@testset "JuliaPlanners repository" begin
    collection_names = list_collections(JuliaPlannersRepo)
    @test isempty(collection_names)
    domain_names = list_domains(JuliaPlannersRepo)
    @test !isempty(domain_names)
    for dname in domain_names
        problem_names = list_problems(JuliaPlannersRepo, dname)
        @test !isempty(problem_names)
        domain = load_domain(JuliaPlannersRepo, dname)
        problem = load_problem(JuliaPlannersRepo, dname, problem_names[1])
    end
    @test find_domains(JuliaPlannersRepo, "blocks") == ["blocksworld", "blocksworld-axioms"]
    @test find_problems(JuliaPlannersRepo, "blocksworld", "1") == ["problem-1"]
end

@testset "IPC Instances repository" begin
    collection_names = list_collections(IPCInstancesRepo)
    @test !isempty(collection_names)
    domain_names = list_domains(IPCInstancesRepo)
    @test !isempty(domain_names)
    domain_names = list_domains(IPCInstancesRepo, "ipc-2000")
    @test !isempty(domain_names)
    for dname in domain_names
        dname = "ipc-2000-$dname"
        problem_names = list_problems(IPCInstancesRepo, dname)
        @test !isempty(problem_names)
        domain = load_domain(IPCInstancesRepo, dname)
        problem = load_problem(IPCInstancesRepo, dname, problem_names[1])
    end
    @test find_domains(IPCInstancesRepo, "ipc-2000", "blocks") ==
        ["blocks-strips-typed", "blocks-strips-untyped"]
    @test find_problems(IPCInstancesRepo, "ipc-2000-blocks-strips-typed", "100") ==
        ["instance-100"]
    clear_cache!(IPCInstancesRepo())
end

@testset "Planning Domains repository" begin
    collection_names = list_collections(PlanningDomainsRepo)
    @test !isempty(collection_names)
    domain_names = list_domains(PlanningDomainsRepo)
    @test !isempty(domain_names)
    domain_names = list_domains(PlanningDomainsRepo, "IPC-2000")
    @test !isempty(domain_names)
    for dname in domain_names
        problem_names = list_problems(PlanningDomainsRepo, dname)
        @test !isempty(problem_names)
        domain = load_domain(PlanningDomainsRepo, dname)
        problem = load_problem(PlanningDomainsRepo, dname, problem_names[1])
    end
    @test find_domains(PlanningDomainsRepo, "IPC-2000", "blocks") ==
        ["112-blocks"]
    @test find_problems(PlanningDomainsRepo, "112-blocks", "probBLOCKS-4-0") ==
        ["3966-probBLOCKS-4-0.pddl"]
    clear_cache!(PlanningDomainsRepo())
end

@testset "Default repository" begin
    collection_names = list_collections()
    @test isempty(collection_names)
    domain_names = list_domains()
    @test !isempty(domain_names)
    for dname in domain_names
        dname = Symbol(dname)
        problem_names = list_problems(dname)
        @test !isempty(problem_names)
        domain = load_domain(dname)
        problem = load_problem(dname, problem_names[1])
    end
    @test find_domains("blocks") == ["blocksworld", "blocksworld-axioms"]
    @test find_problems("blocksworld", "1") == ["problem-1"]
end

clear_all_caches!()


text:
# Reduced mesh composed only of the elements that have shocks in them and their
# neigbors

#TODO: move this to Utils?

#------------------------------------------------------------------------------
# for DG-type shock capturing, need to know interfaces

"""
  Interface that also stores the index that interface corresponds to in the
  parent mesh.
"""
struct RelativeInterface
  iface::Interface  # needed for SBP
  idx_orig::Int
  #peeridx::Int32  # index of the peer on the *original* mesh
end

"""
  Boundary time that also contains the index of the boundary in the parent mesh.
"""
struct RelativeBoundary
  bndry::Boundary
  idx_orig::Int
  fac::Int8  # factor to multiply normal vector by, 1 or -1
end

function RelativeBoundary(bndry::Boundary, idx_orig::Integer)
  return RelativeBoundary(bndry, idx_orig, 1)
end

"""
  Stores a limited about of information about the mesh defined by the elements
  with shocks in them and their neighbors.

  Definitions:

    shocked
    neighbor

  The way this data structure is meant to be used is:

   1) at the beginning of a simulation, construct this object
   2) every time the residual/Jacobian is evaluated, call `push!` repeatedly
      to add the elements that have shocks in them
   3) call [`completeShockMesh`](@ref) to finish constructing the data
      structure
   4) evaluate the shock capturing scheme

  Note that steps 2-4 will happen repeatedly.  This data structure is optimized
  to be re-used several times, under the assumption that the shock will not
  change position too much.  The function [`reset`](@ref) should be called
  before step 2.

  Note that the sizes of the arrays are *minimum* sizes.  In practice they
  may be larger.  Do not call `size()` on these arrays, use the appropriate
  integer field.

  **Fields**

   * elnums_all: the element numbers on the original mesh of the elements on
                 the reduced mesh, length `shockmesh.numEl`  This maps element
                 numbers *from* the reduced mesh *to* the full mesh
   * elnum_mesh: the element numbers on the reduced mesh of the elements
                 on the reduced mesh, length `mesh.numEl`.  This maps element
                 numbers *from* the original mesh *to* the reduced mesh.
   * ifaces: array of [`RelativeInterface`](@ref) objects for the faces
            between elements in the reduced mesh (both elements with shocks in
            them and neighbors. length `shockmesh.numInterfaces`
   * bndryfaces: array of [`RelativeBoundary`](@ref) objects.  Contains only
                 boundary faces that are on the boundary of the original mesh
                 and the element has a shock in it (neighbor elements are
                 not included).
   * numShock: number of elements with shocks in them
   * numNeighbor: number of elements that do not have shocks in them, but
                  are adjacent to elements that do
   * numInterfaces: number of interfaces
   * numBoundaryFaces: number of boundary faces
   * numEl: total number of elements
"""
mutable struct ShockedElements{Tres}
  elnums_all::Vector{Int}  # union of elnums_shock and elnums_neighbor,
                           # and shared elements
  elnums_mesh::Vector{Int}  # temporary array, length mesh.numEl
                            # contains the indices of the selected elements
                            # in elnums_all
                            # contains the indices of elements in elnums_all
                            # TODO: can this be a BitArray?
  ifaces::Vector{RelativeInterface}
  bndryfaces::Vector{RelativeBoundary}

  shared_interfaces::Vector{Vector{RelativeInterface}}

  numShock::Int  # number of elements with shocks in them.  Because we try
                 # to pre-allocate elnums_shock, its length may be greater than
                 # the number of elements with shocks
  numNeighbor::Int  # number of elements that neighbor an element with a shock
                    # in it, but don't have a shock in it
  numShared::Vector{Int}  # number of elements that are owned by other processes
                          # but share a face with a shocked element
  numInterfaces::Int  # number of interfaces in ifaces
  numSharedInterfaces::Vector{Int}  # number of interfaces shared with each
                                    # peer process
  numBoundaryFaces::Int

  npeers::Int  # number of peers
  peer_indices::Vector{Int}  # array containing the *index* of the peer
                             # process in the original mesh list of peer
                             # processes

  numEl::Int  # total number of elements
  isNeumann::Bool  # if true, boundary is Neumann, else Dirichlet

  # useful ranges for iterating
  local_els::UnitRange{Int}
  neighbor_els::UnitRange{Int}
  shared_els::Vector{UnitRange{Int}}
  bndry_offsets::Vector{Int}

  # internal state used for push operations
  idx_all::Int  # current index in elnums_all
  sz_all::Int  # current size of elnums_all

  idx_if::Int  # current index in ifaces
  sz_if::Int  # current size of ifaces

  idx_sf::Vector{Int}  # current index in shared_interfaces
  sz_sf::Vector{Int}   # current size of shared_interfaces

  idx_b::Int  # current index in bndryfaces
  sz_b::Int  # current size of bndryfaces

  function ShockedElements{Tres}(mesh::AbstractMesh) where {Tres}

    # try to guess initial size
    size_guess = max(div(mesh.numEl, 10), 1)
    elnums_all = Array{Int}(size_guess)  # do this later
    elnums_mesh = zeros(mesh.numGlobalEl)
    ifaces = Array{RelativeInterface}(0)
    bndryfaces = Array{RelativeBoundary}(0)
    shared_interfaces = Vector{Vector{RelativeInterface}}(0)

    numShock = 0
    numNeighbor = 0
    numShared = Vector{Int}(0)
    numInterfaces = 0
    numSharedInterfaces = Vector{Int}(0)
    numBoundaryFaces = 0
    npeers = 0
    peer_indices = Vector{Int}(0)
    numEl = 0
    isNeumann = true  # was true
    local_els = 0:0
    neighbor_els = 0:0
    shared_els = Vector{UnitRange{Int}}(0)
    bndry_offsets = Vector{Int}(length(mesh.bndry_offsets))

    idx_all = 1
    sz_all = size_guess
    idx_if = 1
    sz_if = 0
    idx_sf = Vector{Int}(0)
    sz_sf = Vector{Int}(0)
    idx_b = 1
    sz_b = 0

    return new(elnums_all, elnums_mesh,
               ifaces, bndryfaces, shared_interfaces, numShock,
               numNeighbor, numShared,
               numInterfaces, numSharedInterfaces, numBoundaryFaces, npeers,
               peer_indices, numEl, isNeumann,
               local_els, neighbor_els, shared_els, bndry_offsets,
               idx_all, sz_all, idx_if, sz_if, idx_sf, sz_sf, idx_b, sz_b)
  end
end


text:
@testset "CauchyProblem" begin

    @testset "constant solution" begin
    prob = CauchyProblem(identity, b=1, d_dϕb=0)

    θ = solve(prob)

    ϕ = range(0, 20, length=100)

    @test all(@inferred θ.(ϕ) .≈ θ.i)
    @test all(@inferred θ.(ϕ) .≈ θ.b)
    @test all(@inferred d_dϕ.(θ, ϕ) .≈ 0)
    @test θ.iterations == 0
    @test isnan(@inferred θ(-1))
    end

    @testset "exact" begin
    # Reference: Philip (1960) Table 1, No. 13
    # https://doi.org/10.1071/PH600001
    D = θ -> 0.5*(1 - log(θ))

    prob1 = DirichletProblem(D, i=0, b=1)
    θ1 = solve(prob1)

    prob2 = CauchyProblem(D, b=θ1.b, d_dϕb=θ1.d_dϕb)
    θ2 = solve(prob2)

    @test θ2.i == θ1.i
    @test θ2.b == θ1.b
    @test θ2.d_dϕb == θ1.d_dϕb
    @test θ2.ϕi == θ1.ϕi
    @test θ2.ϕb == θ1.ϕb

    @test θ2.iterations == 0

    @test isnan(@inferred θ2(-1))

    end
end


text:
using StatsBase
using BenchmarkTools
using Random

# Types and basic methods
include("EventType.jl")
include("ObjectType.jl")
include("NeighbourObjType.jl")
include("UniverseType.jl")

# Functions
include("Event.jl")
include("Universe.jl")
include("NeighbourObj.jl")
include("BoundaryCondition.jl")
if boundaryCondition == "open"
    using .OpenBoundaryCondition
elseif boundaryCondition == "PPP"
    using .PPPBoundaryCondition
end
include("Object.jl")

function Evolve()
    universe = InitlUniverse()
    while true
        CustomEvents(universe)
        EssentialEvents(universe)
        StopEvents(universe)
        if universe.stop == true
            break
        end
    end
    universe
end

max_stars_repo_name: permissionx/mc-ke

max_issues_repo_name: permissionx/mc-ke

max_forks_repo_name: permissionx/mc-ke


text:
@testset "Save and Load" begin
    @testset "Basic ops" begin
        args = (rand(3, 4), rand(3, 4))
        ort_test(ONNX.add, args...)
    end
end


text:
using NBInclude, Test

@testset "Scoping Tests" begin
    VERSION >= v"1.0.0" ? eval(:(@test_throws LoadError @nbinclude(joinpath("includes", "scoping.ipynb")))) : true
    @nbinclude(joinpath("includes", "scoping.ipynb"); softscope = true)
    @test a == 11
end

@test include(joinpath("includes", "test1.jl")) == 314159
@test f(5) == 6
@test normpath(myfile) == abspath("test.ipynb") * ":In[6]"
@test normpath(myfile2) == abspath("test.ipynb") * ":In[+7]"

x=[]; @nbinclude("test2.ipynb")
@test x == [1, 2, 3, 4, 5, 6]

x=[]; @nbinclude("test2.ipynb"; counters = [1, 4, 5])
@test x == [1, 4, 5]

x=[]; @nbinclude("test2.ipynb"; regex=r"#.*executeme")
@test x == [2, 4]

x=[]; @nbinclude("test2.ipynb"; counters = [1, 4, 5], regex=r"#.*executeme")
@test x == [4]

z = 0; @nbinclude("test2.ipynb"; anshook = x -> (global z += 1))
@test z == 7

@test !in_nbinclude()

let desired_output = "#   This is an example notebook:\n\nfunction f(x)\n    return x + 1\nend\n\nf(3)\n\n# comment\n\nconst myfile = @__FILE__\n\nconst myfile2 = @__FILE__\nf(314158)\n# comment",
    test_ipynb = joinpath(@__DIR__, "test.ipynb"), outfilename = tempname()
    @test sprint(nbexport, test_ipynb) == desired_output
    try
        nbexport(outfilename, test_ipynb)
        @test read(outfilename, String) == desired_output
    finally
        rm(outfilename, force=true)
    end
end


text:
using Test

var"@testcontext" = var"@testset"


text:
using SVG
using Test

function wtest()
    z = IOBuffer()
    write(z, Svg([Circle(0,1,1), Line(0,1,1,2), Polyline([0,1,2], [3,4,5])]), 10, 10)
    seek(z, 0)
    read(z, String) == "<svg width=\"10\" height=\"10\">\n<circle cx=\"0\" cy=\"1\" r=\"1\" style=\"fill:none;stroke:black;stroke-width:1.0\" />\n<line x1=\"0\" y1=\"1\" x2=\"1\" y2=\"2\" style=\"fill:none;stroke:black;stroke-width:1.0\" />\n<polyline points=\"0,3 1,4 2,5 \" style=\"fill:none;stroke:black;stroke-width:1.0\" />\n</svg>\n"
end

@testset "SVG.jl" begin

   @test bounds(Circle(0,1,1)) == (xmin=-1, ymin=0, xmax=1, ymax=2)
   @test bounds(scaled(Circle(0,1,1), identity, y->10y)) == (xmin=-1, ymin=9, xmax=1, ymax=11)
   @test bounds(Line(0,1,1,2)) == (xmin=0, ymin=1, xmax=1, ymax=2)
   @test bounds(Polyline([0,1,2], [3,4,5])) == (xmin=0, ymin=3, xmax=2, ymax=5)
   @test bounds(Polyline([(0,3), (1,4), (2,5)])) == (xmin=0, ymin=3, xmax=2, ymax=5)
   @test bounds(Svg([Circle(0,1,1), Line(0,1,1,2), Polyline([0,1,2], [3,4,5])])) == (xmin=-1, ymin=0, xmax=2, ymax=5)
   @test bounds(scaled(Svg([Circle(0,1,1), Line(0,1,1,2), Polyline([0,1,2], [3,4,5])]), 99,200)) == (xmin=32, ymin=32, xmax=99, ymax=165)

   @test Polyline([0,1,2], [3,4,5]) == Polyline([(0,3), (1,4), (2,5)])
   @test Polyline([0,1,2], [3,4,5], identity, identity) == Polyline([(0,3), (1,4), (2,5)], identity)
   #@test wtest()

end


text:
module CMBLensing

using Adapt
using Base.Broadcast: AbstractArrayStyle, ArrayStyle, Broadcasted, broadcasted,
    DefaultArrayStyle, preprocess_args, Style
using Base.Iterators: repeated, product
using Base.Threads
using Base: @kwdef, @propagate_inbounds, Bottom, OneTo, showarg, show_datatype,
    show_default, show_vector, typed_vcat
using Combinatorics
using DataStructures
using DelimitedFiles
using Distributed: pmap, nworkers, myid, workers
using FFTW
using InteractiveUtils
using Interpolations
using KahanSummation
using Loess
using LinearAlgebra
using LinearAlgebra: diagzero
using MacroTools: @capture, combinedef, isdef, isexpr, postwalk, splitdef
using Match
using Markdown
using Measurements
using Memoization
using Parameters
using Printf
using ProgressMeter
using Random
using Random: seed!
using Roots
using Requires
using Setfield
using StaticArrays: @SMatrix, @SVector, SMatrix, StaticArray, StaticArrayStyle,
    StaticMatrix, StaticVector, SVector
using Statistics
using StatsBase
using Strided: capturestridedargs, make_capture, _mapreduce_fuse!, promoteshape, maybestrided, StridedView

import Adapt: adapt_structure
import Base: +, -, *, \, /, ^, ~, ≈,
    adjoint, axes, broadcast, broadcastable, BroadcastStyle, conj, convert,
    copy, copyto!, eltype, fill!, getindex, getproperty, hash, hcat, hvcat, inv,
    iterate, keys, lastindex, length, literal_pow, mapreduce, materialize!,
    materialize, one, permutedims, print_array, promote, promote_rule,
    promote_rule, promote_type, propertynames, real, setindex!, show,
    show_datatype, show_vector, similar, size, sqrt, string, summary, transpose,
    zero
import Base.Broadcast: instantiate, preprocess
import LinearAlgebra: diag, dot, isnan, ldiv!, logdet, mul!, pinv,
    StructuredMatrixStyle, structured_broadcast_alloc, tr
import Measurements: ±
import Statistics: std

export
    @namedtuple, @repeated, @unpack, animate, argmaxf_lnP, BandPassOp, cache,
    CachedLenseFlow, camb, cov_to_Cℓ, Cℓ_2D, Cℓ_to_Cov, DataSet, DerivBasis,
    diag, Diagonal, DiagOp, dot, EBFourier, EBMap, Field, FieldArray, fieldinfo,
    FieldMatrix, FieldOrOpArray, FieldOrOpMatrix, FieldOrOpRowVector,
    FieldOrOpVector, FieldRowVector, FieldTuple, FieldVector, FieldVector, Flat,
    FlatEB, FlatEBFourier, FlatEBMap, FlatFieldMap, FlatFieldFourier, FlatField,
    FlatFourier, FlatIEBCov, FlatIEBFourier, FlatIEBMap, FlatIQUFourier,
    FlatIQUMap, FlatMap, FlatQU, FlatQUFourier, FlatQUMap, FlatS0, FlatS02,
    FlatS2, FlatS2Fourier, FlatS2Map, Fourier, fourier∂, FuncOp, FΦTuple,
    get_Cℓ, get_Cℓ, get_Dℓ, get_αℓⁿCℓ, get_ρℓ, get_ℓ⁴Cℓ, gradhess, HighPass,
    IdentityOp, IEBFourier, IEBMap, InterpolatedCℓs, IQUFourier, IQUMap,
    LazyBinaryOp, LenseBasis, LenseFlow, LenseOp, LinOp, lnP, load_camb_Cℓs,
    load_sim_dataset, LowPass, make_mask, Map, MAP_joint, MAP_marg, map∂,
    MidPass, mix, nan2zero, noiseCℓs, NoLensing, OuterProdOp, ParamDependentOp,
    pixwin, PowerLens, QUFourier, QUMap, resimulate, RK4Solver, S0, S02, S2,
    sample_joint, shiftℓ, simulate, symplectic_integrate, Taylens, toCℓ, toDℓ,
    tuple_adjoint, ud_grade, unmix, Ð, Ł, δf̃ϕ_δfϕ, δfϕ_δf̃ϕ, ℓ², ℓ⁴, ∇, ∇², ∇¹,
    ∇ᵢ, ∇⁰, ∇ⁱ, ∇₀, ∇₁, ⋅, ⨳

# generic stuff
include("util.jl")
include("util_fft.jl")
include("numerical_algorithms.jl")
include("generic.jl")
include("cls.jl")
include("field_tuples.jl")
include("field_vectors.jl")
include("specialops.jl")

# lensing
include("lensing.jl")
include("lenseflow.jl")
include("powerlens.jl")

# flat-sky maps
include("flat_fftgrid.jl")
include("flat_s0.jl")
include("flat_s2.jl")
include("flat_s0s2.jl")
include("flat_generic.jl")
include("masking.jl")
include("taylens.jl")

# plotting
isjuno = false
function animate end
@init @require Juno="e5e0dc1b-0480-54bc-9374-aad01c23163d" isjuno=Juno.isactive()
@init @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" include("plotting.jl")

# sampling and maximizing the posteriors
include("dataset.jl")
include("posterior.jl")
include("maximization.jl")
include("sampling.jl")

# other estimates
include("quadratic_estimate.jl")

# curved-sky (not yet upgraded to new system)
# include("healpix.jl")

@init @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" include("zygote.jl")

# gpu
@init @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" include("gpu.jl")

end


text:
module fci

# Write your package code here.
#include("lanczos.jl")
include("olsen.jl")

end

max_stars_repo_name: nbraunsc/fci

max_issues_repo_name: nbraunsc/fci

max_forks_repo_name: nbraunsc/fci