text:
struct LineInfoNode
    mod::Module
    method::Symbol
    file::Symbol
    line::Int
    inlined_at::Int
end
const NullLineInfo = LineInfoNode(@__MODULE__, Symbol(""), Symbol(""), 0, 0)

include("compiler/ssair/ir.jl")
include("compiler/ssair/domtree.jl")
include("compiler/ssair/slot2ssa.jl")
include("compiler/ssair/queries.jl")
include("compiler/ssair/passes.jl")
include("compiler/ssair/verify.jl")
include("compiler/ssair/legacy.jl")

macro show(s)
    # return :(println($(QuoteNode(s)), " = ", $(esc(s))))
end

function normalize(@nospecialize(stmt), meta::Vector{Any}, table::Vector{LineInfoNode}, loc::RefValue{Int})
    if isa(stmt, Expr)
        if stmt.head == :meta
            args = stmt.args
            if length(args) > 0
                a1 = args[1]
                if a1 === :push_loc
                    let
                        current = loc[]
                        filename = args[2]::Symbol
                        methodname = NullLineInfo.method
                        mod = table[current].mod
                        line = 0
                        for i = 3:length(args)
                            ai = args[i]
                            if ai isa Symbol
                                methodname = ai
                            elseif ai isa Int32
                                line = Int(ai)
                            elseif ai isa Int64
                                line = Int(ai)
                            elseif ai isa Module
                                mod = ai
                            end
                        end
                        push!(table, LineInfoNode(mod, methodname, filename, line, current))
                        loc[] = length(table)
                    end
                elseif a1 === :pop_loc
                    n = (length(args) > 1) ? args[2]::Int : 1
                    for i in 1:n
                        current = loc[]
                        current = table[current].inlined_at
                        current == 0 && break
                        loc[] = current
                    end
                else
                    push!(meta, stmt)
                end
            end
            return nothing
        elseif stmt.head === :line
            return nothing # deprecated - we shouldn't encounter this
        elseif stmt.head === :gotoifnot
            return GotoIfNot(stmt.args...)
        elseif stmt.head === :return
            return ReturnNode{Any}(stmt.args...)
        end
    elseif isa(stmt, LabelNode)
        return nothing
    elseif isa(stmt, LineNumberNode)
        let # need to expand this node so that it is source-location independent
            current = loc[]
            info = table[current]
            methodname = info.method
            mod = info.mod
            file = stmt.file
            file isa Symbol || (file = info.file)
            line = stmt.line
            push!(table, LineInfoNode(mod, methodname, file, line, info.inlined_at))
            loc[] = length(table)
        end
        return nothing
    end
    return stmt
end

function run_passes(ci::CodeInfo, nargs::Int, linetable::Vector{LineInfoNode})
    mod = linetable[1].mod
    ci.code = copy(ci.code)
    meta = Any[]
    lines = fill(0, length(ci.code))
    let loc = RefValue(1)
        for i = 1:length(ci.code)
            stmt = ci.code[i]
            stmt = normalize(stmt, meta, linetable, loc)
            ci.code[i] = stmt
            if !(stmt === nothing)
                lines[i] = loc[]
            end
        end
    end
    ci.code = strip_trailing_junk!(ci.code, lines)
    cfg = compute_basic_blocks(ci.code)
    defuse_insts = scan_slot_def_use(nargs, ci)
    domtree = construct_domtree(cfg)
    ir = let code = Any[nothing for _ = 1:length(ci.code)]
             argtypes = ci.slottypes[1:(nargs+1)]
            IRCode(code, lines, cfg, argtypes, mod, meta)
        end
    ir = construct_ssa!(ci, ir, domtree, defuse_insts, nargs)
    ir = compact!(ir)
    verify_ir(ir)
    ir = type_lift_pass!(ir)
    ir = compact!(ir)
    verify_ir(ir)
    return ir
end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
#------------------------------------------------------------------------------#
#                                    Module                                    #
#------------------------------------------------------------------------------#

"""
# Description

    module IOT

High-level interface module to the `IsoOrthoTensor` package.

# Usage

This module is automatically `include`d by  the  `IsoOrthoTensor`  package,  and
provides the high-level interface to  the  Isotropic  and  Orthogonality  tensor
submodules.

```julia-repl
julia> using IsoOrthoTensor

julia> typeof(IOT)
Module

```
"""
module IOT

#------------------------------------------------------------------------------#
#                                   Imports                                    #
#------------------------------------------------------------------------------#

using Reexport

#------------------------------------------------------------------------------#
#                                   Includes                                   #
#------------------------------------------------------------------------------#

include("basicTensor.jl")
@reexport using .basicTensor

include("iso.jl")
@reexport using .iso

include("orθ.jl")   # θ: U+3b8
@reexport using .orθ

#------------------------------------------------------------------------------#
#                                   Exports                                    #
#------------------------------------------------------------------------------#

# NOTE: This module is @reexport'ed by the main package module

#------------------------------------------------------------------------------#
#                             Auxiliary Functions                              #
#------------------------------------------------------------------------------#

#------------------------------------------------------------------------------#
#                             Interface Functions                              #
#------------------------------------------------------------------------------#

end # module

ext: jl


text:
module GenSymx
    export gensymx, @gensymx

    numsubscript(s) = begin
        subs = Dict("0"=>"₀", "1"=>"₁", "2"=>"₂", "3"=>"₃", "4"=>"₄", "5"=>"₅", "6"=>"₆", "7"=>"₇", "8"=>"₈", "9"=>"₉");
        replace(s, r"0|1|2|3|4|5|6|7|8|9" => x->subs[x])
    end

    """
    gensymx([tag])
    Generates a symbol which will not conflict with other variable names.
    """
    gensymx(s) = begin
        sym = string(gensym(s))
        i = collect(findlast("#", sym))[1]

        nums = sym[i+1:end]
        sym = "ₓₓ" * sym[3:i-1] * "ₓ" * numsubscript(nums)

        Symbol(sym)
    end

    gensymx() = begin
        sym = string(gensym())
        Symbol("ₓₓ" * numsubscript(sym[3:end]))
    end

    """
        @gensymx
    Generates a gensymx symbol for a variable. For example, `@gensymx x y`
    is transformed into `x = gensymx("x"); y = gensymx("y")`.
    """
    macro gensymx(names...)
    Expr(:block,
       [:($(esc(i)) = gensymx($(string(i)))) for i in names]...,
       nothing
    )
    end

end


text:
function genser_type_for_mime(::Type{M}) where M <: MIME
    throw(ArgumentError("no genser type registered for $M"))
end

@inline function genser_type_for_mime(m::MIME{s}) where s
    genser_type_for_mime(typeof(m))
end

@inline function genser_type_for_mime(m::M) where M <: AbstractString
    genser_type_for_mime(MIME(m))
end

function genser_converter_for_encoding(::Type{E}) where E <: Encoding
    throw(ArgumentError("no genser encoder registered for $E"))
end

@inline function genser_converter_for_encoding(e::Encoding{s}) where s
    genser_converter_for_encoding(typeof(e))
end

@inline function genser_converter_for_encoding(e::E) where E <: AbstractString
    genser_converter_for_encoding(Encoding(e))
end


text:
@inline UNITLESS_ABS2(x::Number) = abs2(x)
@inline UNITLESS_ABS2(x::AbstractArray) = sum(UNITLESS_ABS2, x)
@inline UNITLESS_ABS2(x::ArrayPartition) = sum(UNITLESS_ABS2, x.x)

@inline recursive_length(u::AbstractArray{<:Number}) = length(u)
@inline recursive_length(u::AbstractArray{<:AbstractArray}) = sum(recursive_length, u)
@inline recursive_length(u::ArrayPartition) = sum(recursive_length, u.x)
@inline recursive_length(u::VectorOfArray) = sum(recursive_length, u.u)

@inline ODE_DEFAULT_NORM(u::Union{AbstractFloat,Complex},t) = @fastmath abs(u)
@inline ODE_DEFAULT_NORM(u::Array{T},t) where T<:Union{AbstractFloat,Complex} =
                                         sqrt(real(sum(abs2,u)) / length(u))
@inline ODE_DEFAULT_NORM(u::StaticArray{T},t) where T<:Union{AbstractFloat,Complex} =
                                            sqrt(real(sum(abs2,u)) / length(u))

@inline ODE_DEFAULT_NORM(u::AbstractArray,t) = sqrt(UNITLESS_ABS2(u)/recursive_length(u))
@inline ODE_DEFAULT_NORM(u,t) = norm(u)

@inline ODE_DEFAULT_ISOUTOFDOMAIN(u,p,t) = false
@inline ODE_DEFAULT_PROG_MESSAGE(dt,u,p,t) =
           "dt="*string(dt)*"\nt="*string(t)*"\nmax u="*string(maximum(abs.(u)))
@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u,p,t) = false
@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u::AbstractFloat,p,t) = isnan(u)
@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u::Float64,p,t) =
                                                any(x->(isnan(x) || x>1e50),u)
@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u::AbstractArray{T},p,t) where
                                    {T<:AbstractFloat} = any(isnan,u)
@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u::ArrayPartition,p,t) =
                                                 any(any(isnan,x) for x in u.x)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkInterfaceAssociation <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    securityRules::Any # spec type: Union{ Nothing, Vector{SecurityRule} } # spec name: securityRules

    function NetworkInterfaceAssociation(;id=nothing, securityRules=nothing)
        o = new()
        validate_property(NetworkInterfaceAssociation, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(NetworkInterfaceAssociation, Symbol("securityRules"), securityRules)
        setfield!(o, Symbol("securityRules"), securityRules)
        o
    end
end # type NetworkInterfaceAssociation

const _property_map_NetworkInterfaceAssociation = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("securityRules")=>Symbol("securityRules"))
const _property_types_NetworkInterfaceAssociation = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("securityRules")=>"Vector{SecurityRule}")
Base.propertynames(::Type{ NetworkInterfaceAssociation }) = collect(keys(_property_map_NetworkInterfaceAssociation))
Swagger.property_type(::Type{ NetworkInterfaceAssociation }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_NetworkInterfaceAssociation[name]))}
Swagger.field_name(::Type{ NetworkInterfaceAssociation }, property_name::Symbol) =  _property_map_NetworkInterfaceAssociation[property_name]

function check_required(o::NetworkInterfaceAssociation)
    true
end

function validate_property(::Type{ NetworkInterfaceAssociation }, name::Symbol, val)
end


text:
using lolPackage
using Test

    #TEST 1 WITH DERIVATIVES
    f(x)=x^2
    fprime(x)=2*x
    @test Newtonsroot1(f,fprime,xiv=0.1)[1]≈0.0 atol=0.00001

    f(x)=x^2-16
    fprime(x)=2*x
    @test Newtonsroot1(f,fprime,xiv=3.1)[1]≈4.0 atol=0.00001

    f(x)=(x-2)^2
    fprime(x)=2*(x-2)
    @test Newtonsroot1(f,fprime,xiv=1.0)[1]≈2.0 atol=0.00001

    #TEST 2 WITHOUT DERIVATIVES
    f(x)=x^2

    @test Newtonsroot1(f,xiv=1.0)[1]≈0.0 atol=0.00001

    f(x)=x^2-16
    @test Newtonsroot1(f,xiv=1.0)[1]≈4.0 atol=0.00001

    f(x)=(x-2)^2
    @test Newtonsroot1(f,xiv=1.0)[1]≈2.0 atol=0.00001

    #TEST 3 BIGFLOAT
    f(x)=(x-2)^2
    a=BigFloat(2.0)
    @testset "BigFloat" begin
        @test    Newtonsroot1(f,xiv=1.0)[1]≈2.0 atol=0.00001
        @test    Newtonsroot1(f,xiv=1.0)[1]≈a atol=0.00001
    end

    #TEST 4 tolerance (Accuracy dependent on tolerance)
    f(x)=4x^3-16x+10
    a=Newtonsroot1(f,xiv=1)[1]
    b=Newtonsroot1(f,xiv=1, tol=0.0001)[1]
    c=Newtonsroot1(f,xiv=1, tol=0.01)[1]
    @test f(a)<f(b)<f(c)

    #TEST 5 Non-Convergence
    #Test non-convergence (return nothing)
    f(x)=2+x^2
    @test Newtonsroot1(f,xiv=0.2)==nothing

    #TEST 6 Maxiter
    f(x)=log(x)-20
    a=Newtonsroot1(f,xiv=0.2)[1] #Algorithm needs 17 iterations in this case
    b=Newtonsroot1(f,xiv=0.2,maxiter=5)
    @testset "maxiter" begin
    @test a≈4.851651954097909e8 atol=0.000001
    @test b==nothing
    end;

max_stars_repo_name: anubhavpcjha/lolPackage

max_issues_repo_name: anubhavpcjha/lolPackage

max_forks_repo_name: anubhavpcjha/lolPackage


text:
module TestExamples2DEulerAcoustics

using Test
using Trixi

include("test_trixi.jl")

EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem")

@testset "Acoustic perturbation coupled with compressible Euler" begin
  @trixi_testset "elixir_euleracoustics_co-rotating_vortex_pair.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euleracoustics_co-rotating_vortex_pair.jl"),
                        initial_refinement_level=5,
                        tspan1=(0.0, 1.0), tspan_averaging=(0.0, 1.0), tspan=(0.0, 1.0),
      l2 = [0.00013268029905807722, 0.0001335062197031223, 0.00021776333678401362, 13.000001753042364, 26.00000080243847, 38.00000884725549, 51.000000003859995],
      linf = [0.22312716933051027, 0.1579924424942319, 0.25194831158255576, 13.468872744263273, 26.54666679978679, 38.139032147739684, 51.378134660241294]
    )
  end
end

end # module


text:
#=
"""
    solve(prob::FODESystem, h, PECE())

Use the Adams-Bashforth-Moulton method to solve the system of FODEs.

### References

```tex
@inproceedings{Garrappa2018NumericalSO,
  title={Numerical Solution of Fractional Differential Equations: A Survey and a Software Tutorial},
  author={Roberto Garrappa},
  year={2018}
}
```
"""
=#
#struct PECE <: FractionalDiffEqAlgorithm end

mutable struct M
    an
    bn
    a0
    halpha1
    halpha2
    mu
    mu_tol
    r
    index_fft
    an_fft
    bn_fft
end

function solve(prob::FODESystem, h, ::PECE)
    @unpack f, α, u0, tspan = prob
    t0 = tspan[1]; T = tspan[2]
    METH = M(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)# Initialization
    mu_tol = 1.0e-6
    mu = 1
    α = α[:]

    # Check compatibility size of the problem with number of fractional orders
    alpha_length = length(α)
    problem_size = size(u0, 1)

    m_alpha = ceil.(Int, α)
    m_alpha_factorial = zeros(alpha_length, maximum(m_alpha))
    for i = 1 : alpha_length
        for j = 0 : m_alpha[i]-1
            m_alpha_factorial[i, j+1] = factorial(j)
        end
    end

    f_temp = zeros(size(u0[:, 1]))
    #f_temp = sysf_vectorfield(t0, u0[:, 1], f)
    f(f_temp, u0[:, 1], nothing, t0)

    r::Int = 16
    N::Int = ceil(Int64, (T-t0)/h)
    Nr::Int = ceil(Int64, (N+1)/r)*r
    Qr::Int = ceil(Int64, log2(Nr/r)) - 1
    NNr::Int = 2^(Qr+1)*r

    # Preallocation of some variables
    y = zeros(problem_size, N+1)
    fy = zeros(problem_size, N+1)
    zn_pred = zeros(problem_size, NNr+1)
    mu > 0 ?  (zn_corr = zeros(problem_size, NNr+1)) : (zn_corr = 0)

    # Evaluation of coefficients of the PECE method
    nvett = 0:NNr+1 |> collect
    bn = zeros(alpha_length, NNr+1); an = copy(bn); a0 = copy(bn)
    for i_alpha = 1:alpha_length
        find_alpha = Float64[]
        if α[i_alpha] == α[1:i_alpha-1]
            push!(find_alpha, i_alpha)
        end

        if isempty(find_alpha) == false
            bn[i_alpha, :] = bn[find_alpha[1], :]
            an[i_alpha, :] = an[find_alpha[1], :]
            a0[i_alpha, :] = a0[find_alpha[1], :]
        else
            nalpha = nvett.^α[i_alpha]
            nalpha1 = nalpha.*nvett
            bn[i_alpha, :] = nalpha[2:end] - nalpha[1:end-1]
            an[i_alpha, :] = [1; (nalpha1[1:end-2] - 2*nalpha1[2:end-1] + nalpha1[3:end]) ]
            a0[i_alpha, :] = [0; nalpha1[1:end-2]-nalpha[2:end-1].*(nvett[2:end-1].-α[i_alpha].-1)]
        end
    end
    METH.bn = bn; METH.an = an; METH.a0 = a0
    METH.halpha1 = h.^α./gamma.(α.+1)
    METH.halpha2 = h.^α./gamma.(α.+2)
    METH.mu = mu; METH.mu_tol = mu_tol

    # Evaluation of FFT of coefficients of the PECE method
    METH.r = r
    if Qr >= 0
        index_fft = zeros(Int, 2, Qr+1)
        for l = 1 : Qr+1
            if l == 1
                index_fft[1, l] = 1; index_fft[2, l] = r*2
            else
                index_fft[1, l] = index_fft[2, l-1]+1; index_fft[2, l] = index_fft[2, l-1]+2^l*r
            end
        end

        bn_fft = zeros(Complex, alpha_length, index_fft[2, Qr+1]); an_fft = copy(bn_fft)
        for l = 1 : Qr+1
            coef_end = 2^l*r
            for i_alpha = 1 : alpha_length
                find_alpha = Float64[]
                if α[i_alpha] == α[1:i_alpha-1]
                    push!(find_alpha, i_alpha)
                end
                if isempty(find_alpha) == false
                    bn_fft[i_alpha, index_fft[1, l]:index_fft[2, l]] = bn_fft(find_alpha(1),index_fft(1,l):index_fft(2,l)) ;
                    an_fft[i_alpha, index_fft[1, l]:index_fft[2, l]] = an_fft(find_alpha(1),index_fft(1,l):index_fft(2,l)) ;
                else
                    bn_fft[i_alpha, index_fft[1, l]:index_fft[2, l]] = ourfft(METH.bn[i_alpha, 1:coef_end], coef_end)
                    an_fft[i_alpha, index_fft[1, l]:index_fft[2, l]] = ourfft(METH.an[i_alpha, 1:coef_end], coef_end)
                end
            end
        end
        METH.bn_fft = bn_fft ; METH.an_fft = an_fft ; METH.index_fft = index_fft ;
    end

    # Initializing solution and proces of computation
    t = t0 .+ collect(0:N)*h
    y[:, 1] = u0[:, 1]
    fy[:, 1] = f_temp
    (y, fy) = ABMTriangolo(1, r-1, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f) ;

    # Main process of computation by means of the FFT algorithm
    ff = zeros(1, 2^(Qr+2)); ff[1:2] = [0; 2] ; card_ff = 2
    nx0::Int = 0; ny0::Int = 0
    for qr = 0 : Qr
        L = 2^qr
        (y, fy) = DisegnaBlocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f)
        ff[1:2*card_ff] = [ff[1:card_ff]; ff[1:card_ff]]
        card_ff = 2*card_ff
        ff[card_ff] = 4*L
    end

    # Evaluation solution in T when T is not in the mesh
    if T < t[N+1]
        c = (T - t[N])/h
        t[N+1] = T
        y[:, N+1] = (1-c)*y[:, N] + c*y[:, N+1]
    end

    t = t[1:N+1]; y = y[:, 1:N+1]
    return t, y

end

function DisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f)

    nxi::Int = nx0; nxf::Int = nx0 + L*r - 1
    nyi::Int = ny0; nyf::Int = ny0 + L*r - 1
    is::Int = 1
    s_nxi = zeros(N)
    s_nxf = zeros(N)
    s_nyi = zeros(N)
    s_nyf = zeros(N)
    s_nxi[is] = nxi; s_nxf[is] = nxf; s_nyi[is] = nyi; s_nyf[is] = nyf

    i_triangolo = 0; stop = false
    while stop == false

        stop = (nxi+r-1 == nx0+L*r-1) || (nxi+r-1>=Nr-1)

        (zn_pred, zn_corr) = ABMQuadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length)

        (y, fy) = ABMTriangolo(nxi, nxi+r-1, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f)
        i_triangolo = i_triangolo + 1

        if stop == false
            if nxi+r-1 == nxf
                i_Delta = ff[i_triangolo]
                Delta = i_Delta*r
                nxi = s_nxf[is]+1; nxf = s_nxf[is]  + Delta
                nyi = s_nxf[is] - Delta +1; nyf = s_nxf[is]
                s_nxi[is] = nxi; s_nxf[is] = nxf; s_nyi[is] = nyi; s_nyf[is] = nyf ;
            else # Il triangolo finisce prima del quadrato -> si fa un quadrato accanto
                nxi = nxi + r ; nxf = nxi + r - 1 ; nyi = nyf + 1 ; nyf = nyf + r  ;
                is = is + 1 ;
                s_nxi[is] = nxi ; s_nxf[is] = nxf ; s_nyi[is] = nyi ; s_nyf[is] = nyf ;
            end
        end

    end
    return y, fy
end

function ABMQuadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length)
    coef_end::Int = nxf-nyi+1
    i_fft::Int = log2(coef_end/METH.r)
    funz_beg::Int = nyi+1; funz_end::Int = nyf+1
    Nnxf::Int = min(N, nxf)

    # Evaluation convolution segment for the predictor
    vett_funz = fy[:, funz_beg:funz_end]
    vett_funz_fft = rowfft(vett_funz, coef_end)
    zzn_pred = zeros(problem_size, coef_end)
    for i = 1 : problem_size
        i_alpha = min(alpha_length,i)
        Z = METH.bn_fft[i_alpha, METH.index_fft[1, i_fft]:METH.index_fft[2, i_fft]].*vett_funz_fft[i, :]
        zzn_pred[i, :] = real.(ourifft(Z, coef_end))
    end
    zzn_pred = zzn_pred[:, nxf-nyf:end-1]
    zn_pred[:, nxi+1:Nnxf+1] = zn_pred[:, nxi+1:Nnxf+1] + zzn_pred[:, 1:Nnxf-nxi+1]

    # Evaluation convolution segment for the corrector
    if METH.mu > 0
        if nyi == 0 # Evaluation of the lowest square
            vett_funz = [zeros(problem_size, 1) fy[:, funz_beg+1:funz_end]]
            vett_funz_fft = rowfft(vett_funz, coef_end)
        end
        zzn_corr = zeros(problem_size, coef_end)
        for i = 1 : problem_size
            i_alpha = min(alpha_length,i)
            Z = METH.an_fft[i_alpha,METH.index_fft[1, i_fft]:METH.index_fft[2, i_fft]].*vett_funz_fft[i, :]
            zzn_corr[i, :] = real.(ourifft(Z, coef_end))
        end
        zzn_corr = zzn_corr[:, nxf-nyf+1:end]
        zn_corr[:, nxi+1:Nnxf+1] = zn_corr[:, nxi+1:Nnxf+1] + zzn_corr[:, 1:Nnxf-nxi+1]
    else
        zn_corr = 0
    end
    return zn_pred, zn_corr
end

function ABMTriangolo(nxi, nxf, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f)

    for n = nxi:min(N, nxf)

        # Evaluation of the predictor
        Phi = zeros(problem_size, 1)
        if nxi == 1 # Case of the first triangle
            j_beg::Int = 0
        else # Case of any triangle but not the first
            j_beg = nxi
        end
        for j = j_beg:n-1
            Phi = Phi + METH.bn[1:alpha_length,n-j].*fy[:, j+1]
        end
        St = StartingTerm(t[n+1], u0, m_alpha, t0, m_alpha_factorial)
        y_pred = St + METH.halpha1.*(zn_pred[:, n+1] + Phi)
        f_pred = zeros(length(y_pred))
        #f_pred = sysf_vectorfield(t[n+1], y_pred, f)
        f(f_pred, y_pred, nothing, t[n+1])

        # Evaluation of the corrector
        if METH.mu == 0
            y[:, n+1] = y_pred
            fy[:, n+1] = f_pred
        else
            j_beg = nxi
            Phi = zeros(problem_size, 1)
            for j = j_beg : n-1
                Phi = Phi + METH.an[1:alpha_length, n-j+1].*fy[:, j+1]
            end
            Phi_n = St + METH.halpha2.*(METH.a0[1:alpha_length, n+1].*fy[:, 1] + zn_corr[:, n+1] + Phi)
            yn0 = y_pred; fn0 = f_pred
            stop = false; mu_it = 0
            while stop == false
                global yn1 = Phi_n + METH.halpha2.*fn0
                mu_it = mu_it + 1
                if METH.mu == Inf
                    stop = norm(yn1-yn0, Inf) < METH.mu_tol
                    if mu_it > 100 && ~stop
                        stop = 1
                    end
                else
                    stop = (mu_it == METH.mu)
                end
                global fn1 = zeros(length(yn1))
                #sysf_vectorfield(t[n+1], yn1, f)
                f(fn1, yn1, nothing, t[n+1])
                yn0 = yn1; fn0 = fn1
            end
            y[:, n+1] = yn1
            fy[:, n+1] = fn1
        end
    end
    return y, fy
end

sysf_vectorfield(t, y, f_fun) = f_fun(t, y)

function  StartingTerm(t, u0, m_alpha, t0, m_alpha_factorial)
    ys = zeros(size(u0, 1), 1)
    for k = 1 : maximum(m_alpha)
        if length(m_alpha) == 1
            ys = ys + (t-t0)^(k-1)/m_alpha_factorial[k]*u0[:, k]
        else
            i_alpha = findall(x -> x>=k, m_alpha)
            ys[i_alpha, 1] = ys[i_alpha, 1] + (t-t0)^(k-1)*u0[i_alpha, k]./m_alpha_factorial[i_alpha, k]
        end
    end
    return ys
end


text:
"""
温度流的Bubble，ultra low temperature
"""

"""
在极低温的情况下，将只有k和k′都在费米面上面的时候才有贡献
"""

"""
获取所有的Bubble(pp, fs, ex)
qpp = k1+k2; qfs = k3-k2; qex = k1-k3
"""
function all_bubble_tf_ult_mt(Γ4::Gamma4{T, P}, cert::Float64) where {T, P}
    #
    brlu_area = area(Γ4.model.brillouin)
    #因为现在所有的能带都必须有同一个lpats，所以只要每一个patch中有哪些tri
    tris_pat = Γ4.ltris_pat
    #
    bubbval_pp = Array{Float64, 7}(
        undef,
        #alpha, beta, b1, b2
        Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum,
        #nidx(alpha), k1(b1), k2(b2)
        Γ4.patchnum, Γ4.patchnum, Γ4.patchnum
    )
    Threads.@threads for idxs in CartesianIndices(bubbval_pp)
        alpha, beta, b1, b2, i_n, i1, i2 = Tuple(idxs)
        k1, k2 = Γ4.patches[b1][i1], Γ4.patches[b2][i2]
        q_pp = Γ4.model.kadd(k1, k2)
        bubbres = pi_αβ_minus_tf_ult(
            tris_pat[i_n], brlu_area, cert,
            q_pp, Γ4.model.dispersion[alpha], Γ4.model.dispersion[beta]
        )
        bubbval_pp[alpha, beta, b1, b2, i_n, i1, i2] = bubbres
    end
    bubb_qpp = Bubble{:tf, :minus}(Inf, bubbval_pp)
    #
    bubbval_fs = Array{Float64, 7}(
        undef,
        #alpha, beta, b2, b3
        Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum,
        #nidx(alpha), k2(b2), k3(b3)
        Γ4.patchnum, Γ4.patchnum, Γ4.patchnum
    )
    Threads.@threads for idxs in CartesianIndices(bubbval_fs)
        alpha, beta, b2, b3, i_n, i2, i3 = Tuple(idxs)
        k2, k3 = Γ4.patches[b2][i2], Γ4.patches[b3][i3]
        q_fs = Γ4.model.kadd(k3, -k2)
        bubbres = pi_αβ_plus_tf_ult(
            tris_pat[i_n], brlu_area, cert,
            q_fs, Γ4.model.dispersion[alpha],
            Γ4.model.dispersion[beta]
        )
        bubbval_fs[alpha, beta, b2, b3, i_n, i2, i3] = bubbres
    end
    bubb_qfs = Bubble{:tf, :plus}(Inf, bubbval_fs)
    #
    bubbval_ex = Array{Float64, 7}(
        undef,
        #alpha, beta, b1, b3
        Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum,
        #nidx(alpha), k1(b3), k3(b3)
        Γ4.patchnum, Γ4.patchnum, Γ4.patchnum
    )
    Threads.@threads for idxs in CartesianIndices(bubbval_ex)
        alpha, beta, b1, b3, i_n, i1, i3 = Tuple(idxs)
        k1, k3 = Γ4.patches[b1][i1], Γ4.patches[b3][i3]
        q_ex = Γ4.model.kadd(k1, -k3)
        bubbres = pi_αβ_plus_tf_ult(
            tris_pat[i_n], brlu_area, cert,
            q_ex, Γ4.model.dispersion[alpha],
            Γ4.model.dispersion[beta]
        )
        bubbval_ex[alpha, beta, b1, b3, i_n, i1, i3] = bubbres
    end
    bubb_qex = Bubble{:tf, :plus}(Inf, bubbval_ex)
    return bubb_qpp, bubb_qfs, bubb_qex
end

"""
温度流的Bubble, 这里必须要保证ltris在费米面上面
"""
function pi_αβ_plus_tf_ult(
    ltris::Vector{T},
    area::Float64, cert::Float64,
    qval::Point2D,
    dispα::Function, dispβ::Function) where T <: Basics.AbstractTriangle
    """温度流的+
    这里的lamb就是T，ltris中的所有三角都应该要在同一个patch中,
    tarea是每个小三角形的面积(已经不用了)，
    dispa是和k相关的那个能带，dispb是k-q相关的
    """
    nega_q = -qval
    result = 0.
    for tri in ltris
        #这个小三角形的k值
        kval = tri.center
        #k-q
        kprim = kval + nega_q
        #kprim = kadd(kval, nega_q)
        #epsilon_k
        eps_k = dispα(kval.x, kval.y)
        if !isapprox(eps_k, 0., atol=cert)
            continue
        end
        #epsilon_{k-q}
        eps_kp = dispβ(kprim.x, kprim.y)
        if !isapprox(eps_kp, 0., atol=cert)
            continue
        end
        #如果有一个严格等于0,就换个位置算
        vidx = 1
        while isapprox(eps_kp, 0.) || isapprox(eps_k, 0.)
            kval = tri.vertex[vidx]
            kprim = kval + nega_q
            #kprim = kadd(kval, nega_q)
            eps_k = dispα(kval.x, kval.y)
            eps_kp = dispβ(kprim.x, kprim.y)
            vidx += 1
        end
        #必须反号
        signed_k = eps_k / eps_kp
        if signed_k > 0
            continue
        end
        #这个小区域的贡献
        result += 2 * tri.edges[1].length / (abs(signed_k) + 1)
    end#对ltris的循环
    result = result / area
    return result
end

"""
温度流的Bubble, 这里必须要保证ltris在费米面上面
"""
function pi_αβ_minus_tf_ult(
    ltris::Vector{T},
    area::Float64, cert::Float64,
    qval::Point2D,
    dispα::Function, dispβ::Function) where T <: Basics.AbstractTriangle
    result = 0.
    for tri in ltris
        #这个小三角形的k值
        kval = tri.center
        nega_k = -kval
        #-k+q
        kprim = nega_k + qval
        #kprim = kadd(nega_k, qval)
        #epsilon_k
        eps_k = dispα(kval.x, kval.y)
        if !isapprox(eps_k, 0., atol=cert)
            continue
        end
        #-epsilon_{-k+q}
        neps_kp = -dispβ(kprim.x, kprim.y)
        if !isapprox(neps_kp, 0., atol=cert)
            continue
        end
        #如果其中一个完全等于0,就换个位置算
        vidx = 1
        while isapprox(neps_kp, 0.) || isapprox(eps_k, 0.)
            kval = tri.vertex[vidx]
            nega_k = -kval
            kprim = nega_k + qval
            #kprim = kadd(nega_k, qval)
            eps_k = dispα(kval.x, kval.y)
            neps_kp = -dispβ(kprim.x, kprim.y)
            vidx += 1
        end
        #必须反号，注意这里neps_kp是已经带负号的
        signed_k = eps_k / neps_kp
        if signed_k > 0
            continue
        end
        #计算这个小区域的贡献
        result += 2 * tri.edges[1].length / (abs(signed_k) + 1)
    end#对ltris的循环
    result = result / area
    return result
end


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ 775d21d8-2609-11ec-16c3-f9468cf21068
if isdefined(@__MODULE__, :PlutoRunner)
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.develop(; url = dirname(@__DIR__))
	Pkg.add("PlutoTest")
	Pkg.add("Revise")
	using Revise
	using PlutoSerialization
	using PlutoTest
end

# ╔═╡ 4df95dc3-6106-4150-a1c6-2735d1d38106
module M

struct T
	x
end

end

# ╔═╡ e55bcebd-6efa-4c34-8019-b51c2bc20302
# serialize("M.T.jls", M.T([]))

# ╔═╡ 264a1a7b-1296-4eb9-9fee-38f7ec383287
@test deserialize("M.T.jls") isa M.T

# ╔═╡ 63db7529-7c44-4f75-a0fb-acb6753177d9
@test deserialize("M.T.jls").x isa Vector{Any}

# ╔═╡ 07a5518e-235e-4cdc-a7cb-cfb402837a73
func = x -> x + 1

# ╔═╡ 12442b83-5459-4f54-9ee3-9a2e87665bd3
# serialize("func.jls", M.T(func))

# ╔═╡ 1e2eba1f-b5fc-4f9b-9cc7-2088965373b9
@test deserialize("func.jls") isa M.T

# ╔═╡ fca9ab30-0819-47d4-a74d-72ac6bc2b58b
@test deserialize("func.jls").x isa Function

# ╔═╡ 0af7c4d0-08e8-421e-b158-e98190cefcdf
@test deserialize("func.jls").x(1) == 2

# ╔═╡ 4a64a0e8-6254-44fa-a50d-0976e1d8e8e9
function closure(x)
	function (y)
		x + y
	end
end

# ╔═╡ dfd9093d-5e7e-4db1-9d9e-dc1880105510
# serialize("closure.jls", M.T(closure(big"1")))

# ╔═╡ 062caef8-9831-4b8e-b05b-1eb7162bd38c
@test deserialize("closure.jls") isa M.T

# ╔═╡ a9ace2cc-8814-4341-abc4-5c3be976b34f
@test deserialize("closure.jls").x isa Function

# ╔═╡ 9fffa093-ed3a-4647-986f-addb484b1fb8
@test deserialize("closure.jls").x(big"1") isa BigInt

# ╔═╡ 89e7273f-8d43-46b8-a44e-281dd45698f7
@test deserialize("closure.jls").x(1) == big"2"

# ╔═╡ Cell order:
# ╠═775d21d8-2609-11ec-16c3-f9468cf21068
# ╠═4df95dc3-6106-4150-a1c6-2735d1d38106
# ╠═e55bcebd-6efa-4c34-8019-b51c2bc20302
# ╠═264a1a7b-1296-4eb9-9fee-38f7ec383287
# ╠═63db7529-7c44-4f75-a0fb-acb6753177d9
# ╠═07a5518e-235e-4cdc-a7cb-cfb402837a73
# ╠═12442b83-5459-4f54-9ee3-9a2e87665bd3
# ╠═1e2eba1f-b5fc-4f9b-9cc7-2088965373b9
# ╠═fca9ab30-0819-47d4-a74d-72ac6bc2b58b
# ╠═0af7c4d0-08e8-421e-b158-e98190cefcdf
# ╠═4a64a0e8-6254-44fa-a50d-0976e1d8e8e9
# ╠═dfd9093d-5e7e-4db1-9d9e-dc1880105510
# ╠═062caef8-9831-4b8e-b05b-1eb7162bd38c
# ╠═a9ace2cc-8814-4341-abc4-5c3be976b34f
# ╠═9fffa093-ed3a-4647-986f-addb484b1fb8
# ╠═89e7273f-8d43-46b8-a44e-281dd45698f7


text:
@auto_hash_equals struct Location
    file::String
    line::Int
end
file(loc::Location) = loc.file
line(loc::Location) = loc.line
file(x) = file(location(x))
line(x) = line(location(x))
function Base.show(io::IO, loc::Location)
    if loc == UNKNOWN_LOCATION
        print(io, "unknown")
    else
        print(io, file(loc), ":", line(loc))
    end
end

const UNKNOWN_LOCATION = Location("unknown", -1)


text:
module Coord

using Gadfly
using Compose
using DataArrays
import Gadfly.Maybe
import Iterators: cycle

export cartesian

# Cartesian coordinates with position given by the x and y (and similar)
# aesthetics.
#
# Args:
#   xvars: Aesthetics to consider when choosing x bounds.
#   yvars: Aesthetics to consider when choosing y bounds.
#   xmin, xmax, ymin, ymax:
#     Force a particular x or y bound, rather than trying to choose it
#     from the data.
#
#
immutable Cartesian <: Gadfly.CoordinateElement
    xvars::Vector{Symbol}
    yvars::Vector{Symbol}
    xmin
    xmax
    ymin
    ymax
    xflip
    yflip

    function Cartesian(
            xvars=[:x, :xviewmin, :xviewmax, :xmin, :xmax],
            yvars=[:y, :yviewmin, :yviewmax, :ymin, :ymax, :middle,
                   :lower_hinge, :upper_hinge, :lower_fence, :upper_fence, :outliers];
            xflip::Bool=false, yflip::Bool=false,
            xmin=nothing, xmax=nothing, ymin=nothing, ymax=nothing)
        new(xvars, yvars, xmin, xmax, ymin, ymax, xflip, yflip)
    end
end

const cartesian = Cartesian

# Produce a canvas with suitable cartesian coordinates.
#
# Args:
#   coord: cartesian coordinate instance.
#
# Returns:
#   A compose Canvas.
#
function apply_coordinate(coord::Cartesian, aess::Gadfly.Aesthetics...)
    xmin = nothing
    xmax = nothing
    for var in coord.xvars
        for aes in aess
            if getfield(aes, var) === nothing
                continue
            end

            vals = getfield(aes, var)
            if !isa(vals, AbstractArray)
                vals = {vals}
            end

            for val in vals
                if !Gadfly.isconcrete(val)
                    continue
                end

                if xmin === nothing || val < xmin
                    xmin = val
                end

                if xmax === nothing || val > xmax
                    xmax = val
                end
            end
        end
    end

    ymin = nothing
    ymax = nothing
    for var in coord.yvars
        for aes in aess
            if getfield(aes, var) === nothing
                continue
            end

            # Outliers is an odd aesthetic that needs special treatment.
            if var == :outliers
                for vals in aes.outliers
                    for val in vals
                        if !Gadfly.isconcrete(val)
                            continue
                        end

                        if ymin === nothing || val < ymin
                            ymin = val
                        end

                        if ymax === nothing || val > ymax
                            ymax = val
                        end
                    end
                end

                continue
            end

            vals = getfield(aes, var)
            if !isa(vals, AbstractArray)
                vals = {vals}
            end

            for val in vals
                if !Gadfly.isconcrete(val)
                    continue
                end

                if ymin === nothing || val < ymin
                    ymin = val
                end

                if ymax === nothing || val > ymax
                    ymax = val
                end
            end
        end
    end

    xmin = coord.xmin === nothing ? xmin : coord.xmin
    xmax = coord.xmax === nothing ? xmax : coord.xmax
    ymin = coord.ymin === nothing ? ymin : coord.ymin
    ymax = coord.ymax === nothing ? ymax : coord.ymax

    if xmin === nothing || !isfinite(xmin)
        xmin = 0.0
        xmax = 1.0
    end

    if ymin === nothing || !isfinite(ymin)
        ymin = 0.0
        ymax = 1.0
    end

    # A bit of kludge. We need to extend a bit to be able to fit bars when
    # using discrete scales. TODO: Think more carefully about this. Is there a
    # way for the geometry to let the coordinates know that a little extra room
    # is needed to draw everything?
    if all([aes.x === nothing || typeof(aes.x) <: PooledDataArray for aes in aess])
        xmin -= 0.5
        xmax += 0.5
        xpadding = 0.0
    else
        xpadding = 0.03 * (xmax - xmin)
    end

    if all([aes.y === nothing || typeof(aes.y) <: PooledDataArray for aes in aess])
        ymin -= 0.5
        ymax += 0.5
        ypadding = 0.0
    else
        ypadding = 0.03 * (ymax - ymin)
    end

    width  = xmax - xmin + 2.0 * xpadding
    height = ymax - ymin + 2.0 * ypadding

    compose(
        canvas(unit_box=UnitBox(
            coord.xflip ? xmax + xpadding : xmin - xpadding,
            coord.yflip ? ymin - ypadding : ymax + ypadding,
            coord.xflip ? -width : width,
            coord.yflip ? height : -height)),
        svgclass("plotpanel"))
end

end # module Coord


text:
# ================================================
# Internal forward network functions
# ================================================

# output of neural network for a single input
function forward(nnet::Network, x0::Vector{<:Number})
    x = x0
    @inbounds for layer in nnet.layers
        W = layer.weights
        b = layer.bias
        x = layer.activation(W * x + b)
    end
    return x
end

# ================================================
# Composite methods to compute the network output
# ================================================

@with_kw struct SplitSolver{S<:Solver, FS, FM} <: Solver
    solver::S
    split_fun::FS
    merge_fun::FM
end

function SplitSolver(solver)
    # default: box approximation and split in two sets per dimension
    split_fun = X -> split(box_approximation(X), 2 * ones(Int, dim(X)))
    # default: box approximation of the union
    merge_fun = X -> box_approximation(X)
    return SplitSolver(solver, split_fun, merge_fun)
end

function NeuralVerification.forward_network(solver::SplitSolver, nnet::Network, X0)
    X0_split = solver.split_fun(X0)
    Y_union = UnionSetArray()
    for X in X0_split
        Y = forward_network(solver.solver, nnet, X)
        push!(array(Y_union), Y)
    end
    Y_merged = solver.merge_fun(Y_union)
    return Y_merged
end

# ============================================================================
# Methods to approximate the network output (without mathematical guarantees)
# ============================================================================

# solver using the CH of the sampled outputs as an inner approx of the real output
@with_kw struct SampledApprox <: Solver
    nsamples::Int = 10000
    include_vertices::Bool = true
    directions = OctDirections
end

function NeuralVerification.forward_network(solver::SampledApprox, nnet, input)
    samples = sample(input, solver.nsamples;
                     include_vertices=solver.include_vertices)

    m = output_dim(nnet)
    if m == 1
        MIN = Inf
        MAX = -Inf
        for sample in samples
            output = first(NV.compute_output(nnet, sample))
            MIN = min(MIN, output)
            MAX = max(MAX, output)
        end
        return Interval(MIN, MAX)
    else
        vlist = Vector{Vector{eltype(samples[1])}}(undef, length(samples))
        @inbounds for (i, sample) in enumerate(samples)
            vlist[i] = NV.compute_output(nnet, sample)
        end
        convex_hull!(vlist)
        P = VPolytope(vlist)
        Z = overapproximate(P, Zonotope, solver.directions)
        return Z
    end
end

# ==========================================================
# Methods to handle networks with ReLU activation functions
# ==========================================================

# solver that computes the box approximation in each layer
# it exploits that box(relu(X)) == relu(box(X))
struct BoxSolver <: Solver end

function NeuralVerification.forward_network(solver::BoxSolver, nnet::Network, X0)
    X = X0
    for layer in nnet.layers
        # affine map and box approximation
        W = layer.weights
        b = layer.bias
        X_am = AffineMap(W, X, b)

        # activation function
        if layer.activation isa Id
            X = X_am
            continue
        end

        @assert layer.activation isa ReLU "unsupported activation function"
        X = rectify(box_approximation(X_am))
    end
    return X
end

@with_kw struct ConcreteReLU <: Solver
    concrete_intersection::Bool = false
    convexify::Bool = false
end

function NeuralVerification.forward_network(solver::ConcreteReLU, nnet::Network, X0)
    X = [X0]
    for layer in nnet.layers
        if typeof(X[1]) <: UnionSetArray
            X = [x.array for x in X]
            X = reduce(vcat, X)
        end
        X = affine_map.(Ref(layer.weights), X, Ref(layer.bias))

        # activation function
        if layer.activation isa Id
            continue
        end
        @assert layer.activation isa ReLU "unsupported activation function"
        X = rectify.(X, solver.concrete_intersection)
    end
    return solver.convexify ? ConvexHullArray(X) : X
end

# solver that propagates the vertices, computes their convex hull, and applies
# some postprocessing to the result
@with_kw struct VertexSolver{T} <: Solver
    postprocessing::T = x -> x  # default: identity (= no postprocessing)
    apply_convex_hull::Bool = false
end

function NeuralVerification.forward_network(solver::VertexSolver, nnet::Network, X0)
    N = eltype(X0)
    P = X0

    for layer in nnet.layers
        # apply affine map
        W = layer.weights
        b = layer.bias
        P = convert(VPolytope, P)
        Q_am = affine_map(W, P, b)

        # activation function
        if layer.activation isa Id
            P = Q_am
            continue
        end

        @assert layer.activation isa ReLU "unsupported activation function"

        # compute Q_chull = convex_hull(Q_am, rectify(Q_am))
        vlist = Vector{Vector{N}}()
        vlist_rect = Vector{Vector{N}}()
        for v in vertices(Q_am)
            push!(vlist, v)
            v_rect = rectify(v)
            push!(vlist_rect, v_rect)
            if v != v_rect
                push!(vlist, v_rect)
            end
        end
        if solver.apply_convex_hull || true
            convex_hull!(vlist)
        end
        Q_chull = VPolytope(vlist)

        # filter out negative part
#         Q_pos = box_approximation(VPolytope(vlist_rect))  # alternative
        n = dim(Q_am)
        Q_pos = HPolyhedron(
            [HalfSpace(SingleEntryVector(i, n, -one(N)), zero(N)) for i in 1:n])
        P = intersection(Q_chull, Q_pos)
    end
    Q = solver.postprocessing(P)
    return Q
end

# ==============================================================================
# Methods to handle networks with sigmoid activation functions from [VER19]
#
# [VER19]: Ivanov, Radoslav, et al. "Verisig: verifying safety properties of
# hybrid systems with neural network controllers." Proceedings of the 22nd ACM
# International Conference on Hybrid Systems: Computation and Control. 2019.
# ==============================================================================

# ref. Eq (6) in [VER19]
# d(σ(x))/dx = σ(x)*(1-σ(x))
# g(t, x) = σ(tx) = 1 / (1 + exp(-tx))
# dg(t, x)/dt = g'(t, x) = x * g(t, x) * (1 - g(t, x))
@taylorize function sigmoid!(dx, x, p, t)
    xᴶ, xᴾ = x
    dx[1] = zero(xᴶ)
    dx[2] = xᴶ *(xᴾ - xᴾ^2)
end

# footnote (3) in [VER19]
# d(tanh(x))/dx = 1 - tanh(x)^2
# g(t, x) = tanh(tx)
# dg(t, x)/dt = g'(t, x) = x * (1 - g(t, x)^2)
@taylorize function tanh!(dx, x, p, t)
    xᴶ, xᴾ = x
    dx[1] = zero(xᴶ)
    dx[2] = xᴶ *(1 - xᴾ^2)
end

const HALFINT = IA.Interval(0.5, 0.5)
const ZEROINT = IA.Interval(0.0, 0.0)
const ACTFUN = Dict(Tanh() => (tanh!, ZEROINT),
                    Sigmoid() => (sigmoid!, HALFINT))

# Method: Cartesian decomposition (intervals for each one-dimensional subspace)
# Only Tanh, Sigmoid and Id functions are supported
function forward(nnet::Network, X0::LazySet;
                 alg=TMJets(abstol=1e-14, orderQ=2, orderT=6))

    # initial states
    xᴾ₀ = _decompose_1D(X0)
    xᴾ₀ = LazySets.array(xᴾ₀)  # see https://github.com/JuliaReach/ReachabilityAnalysis.jl/issues/254
    xᴾ₀ = [x.dat for x in xᴾ₀] # use concrete inteval matrix-vector operations

    for layer in nnet.layers  # loop over layers
        W = layer.weights
        m, n = size(W)
        b = layer.bias
        act = layer.activation

        xᴶ′ = W * xᴾ₀ + b  # (scalar matrix) * (interval vector) + (scalar vector)

        if act == Id()
            xᴾ₀ = copy(xᴶ′)
            continue
        end

        activation!, ival = ACTFUN[act]
        xᴾ′ = fill(ival, m)

        for i = 1:m  # loop over coordinates
            X0i = xᴶ′[i] × xᴾ′[i]
            ivp = @ivp(x' = activation!(x), dim=2, x(0) ∈ X0i)
            sol = RA.solve(ivp, tspan=(0., 1.), alg=alg)

            # interval overapproximation of the final reach-set along
            # dimension 2, which corresponds to xᴾ
            xᴾ_end = sol.F.ext[:xv][end][2]
            xᴾ′[i] = xᴾ_end
        end
        xᴾ₀ = copy(xᴾ′)
    end
    return CartesianProductArray([Interval(x) for x in xᴾ₀])
end

# =====================================
# Methods to handle arbitrary networks
# =====================================

struct BlackBoxSolver <: Solver end

struct BlackBoxController{FT} <: AbstractNetwork
    f::FT
end

function NeuralVerification.forward_network(solver::BlackBoxSolver,
                                            bbc::BlackBoxController, X0)
    return bbc.f(X0)
end

function forward(bbc::BlackBoxController, X0)
    return bbc.f(X0)
end

# =============================
# Handling of singleton inputs
# =============================

function forward(nnet::Network, X0::AbstractSingleton)
    x0 = element(X0)
    x1 = forward(nnet, x0)
    return Singleton(x1)
end

for SOLVER in LazySets.subtypes(Solver, true)
    @eval function NeuralVerification.forward_network(solver::$SOLVER,
                                                      nnet::Network,
                                                      X0::AbstractSingleton)
              return forward(nnet, X0)
          end
end


text:
import Base: +
struct Vertex2D{T<:Real}
    x::T
    y::T
end
Base.show(io::IO, v::Vertex2D) = print(io, "(x = $(v.x),y = $(v.y))")
Base.isless(a::Vertex2D, b::Vertex2D) = Base.isless(a.x, b.x) && Base.isless(a.y, b.y)
Base.isequal(a::Vertex2D, b::Vertex2D) = Base.isequal(a.x, b.x) && Base.isequal(a.y, b.y)
Base.promote_rule(::Type{Vertex2D{T}}, ::Type{Vertex2D{S}}) where {T<:Real,S<:Real} =
    Vertex2D{promote_type(T, S)}
Vertex2D(x::T, y::T) where {T<:Integer} = Vertex2D(convert(Float64, x), convert(Float64, y))
+(a::Vertex2D, b::Vertex2D) = Vertex2D(a.x + b.x, a.y + b.y)

Vertex2D(v::Vector{T}) where {T} =  length(v) == 2 ? Vertex2D(v[1],v[2]) : throw(0Base.DimensionMismatch("only 2-d vectors allowed"))
Base.vec(v::Vertex2D{T} where {T}) = Base.vec([v.x, v.y])


text:
module BackedUpImmutable

import Base.getindex, Base.setindex!, Base.get!, Base.get, Base.empty!, Base.pop!
import Base.haskey, Base.delete!, Base.iterate, Base.length

export StaticDict, BackedUpImmutableDict, restore!, restoreall!

""" Another name for ImmutableDict, but here as with an extra constructor. """
StaticDict = Base.ImmutableDict

""" Constructor for StaticDict / ImmutableDict to take an array of key value pairs. """
function Base.ImmutableDict(pairs::Vector{Pair{K,V}}) where V where K
    id = Base.ImmutableDict(pairs[1][1] => pairs[1][2])
    for p in pairs[2:end]
        id = StaticDict(id, p[1] => p[2])
    end
    id
end

""" Constructor for StaticDict to take a series of pairs, varargs style """
function Base.ImmutableDict(pairs...)
    pairvect = [pairs...]
    id = Base.ImmutableDict(pairvect[1][1] => pairvect[1][2])
    for p in pairvect[2:end]
        id = StaticDict(id, p[1] => p[2])
    end
    id
end

"""
    # BackedUpImmutableDict{K, V}
    * Combines a key, not value, immutable hash dictionary with a backup of the original value defaults.
    * For configuration data storage, with a simple restore to default
"""
mutable struct BackedUpImmutableDict{K, V} <: AbstractDict{K,V}
    d::StaticDict
    defaults::Dict{K, V}
end

"""
    Makes a BackedUpImmutableDict from a vector of key, value pairs
"""
BackedUpImmutableDict{K,V}(pairs::Vector{Pair{K,V}}) where V where K =
    BackedUpImmutableDict(StaticDict(pairs), Dict{K,V}(pairs...))

"""
    Makes a BackedUpImmutableDict from a tuple of key, value pairs (varargs style)
"""
BackedUpImmutableDict{K,V}(pairs...) where V where K = BackedUpImmutableDict{K,V}([pairs...])

Base.haskey(dic::BackedUpImmutableDict, k) = haskey(dic.d[k])
Base.getindex(dic::BackedUpImmutableDict, k) = getindex(dic.d, k)
Base.get(dic::BackedUpImmutableDict) = get(dic.d, k)
Base.delete!(dic::BackedUpImmutableDict, k) = throw("Cannot delete from an ImmutableDict")
Base.empty!(dic::BackedUpImmutableDict) = throw("Cannot empty! an ImmutableDict")
Base.pop!(dic::BackedUpImmutableDict) = throw("Cannot pop! from an ImmutableDict")
Base.length(dic::BackedUpImmutableDict) = length(dic.d)
Base.iterate(dic::BackedUpImmutableDict) = iterate(dic.d)
Base.iterate(dic::BackedUpImmutableDict, s) = iterate(dic.d, s)

function Base.setindex!(dic::BackedUpImmutableDict, v, k)
    if haskey(dic.d, k)
        id = Base.ImmutableDict(dic.d, k => v)
        dic.d = id
    else
        throw("Cannot add key $k to ImmutableDict")
    end
end

function Base.get!(dic::BackedUpImmutableDict, k, v)
    if haskey(dic.d, k)
        return get(dic.c, k)
    else
        throw("Cannot add key $k to ImmutableDict")
    end
end

"""
   # Restore a key's backed up default value.
"""
function restore!(dic, k)
    if haskey(dic.defaults, k)
        dic[k] = (v = dic.defaults[k])
        return v
    end
end

"""
   # Restore all values back to defaults
"""
function restoreall!(dic)
    dic.d = StaticDict(collect(dic.defaults))
end

end # module


text:
# generic solve_hjb( ) Implicit method.
# TODO: ns states/nc controls
# TODO: add uncertainty, nx shocks...

using LinearAlgebra, SparseArrays, Plots

# NGM ###############################
if 1==1
    σ= 2.0; ρ = 0.05; δ = 0.05; A = 1.0; α= 0.3;

    # Closed form. Doesn't work.
    # σ= 4.0;  δ = 0.05; A = 1.0; α= 0.3;
    # ρ = δ*(α*σ - 1) # restriction for closed form
    # FOC(v;σ=σ) = abs(v)^(-1.0/σ)                  # FOC, u'^-1()

    r(c; σ=σ)  = (c^(1-σ))/(1-σ)             # return fcn.      make it r(s,c)
    μ(s,c; α=α,δ=δ,A=A) = A*(s^α) - δ*s -c;  # Transition fcn.  ṡ=μ(s,c)

    dr(c;σ=σ) = c^(-σ)                       # derivative of return fcn.
    FOC(v;σ=σ) = v^(-1.0/σ)                  # FOC, u'^-1()

    μ_inv(s,ṡ; α=α,δ=δ,A=A) = A*(s^α) - δ*s -ṡ; # Inv Trans c=μ(s,ṡ), solve ṡ=μ(s,c) for c
    # con0(s;α=α,δ=δ,A=A)     = A*(s^α) - δ*s      # solve: ṡ=μ(s,c) for c(s) if ṡ=0.
    # μ_inv(s,0) = con0(s)

    s_ss = (α*A/(ρ+δ))^(1/(1-α))
    s_min = 0.001*s_ss
    s_max = 2.000*s_ss
    H = 10_000;
    s = LinRange(s_min, s_max, H)
    s = convert(Array, s)
    ds = (s_max-s_min)/(H-1)
end

# con-sav ###############################
if 1==1
    a_bar = -0.02;
    #a_bar = 0.1;
    σ= 2.0; ρ = 0.05; A=0.045; w = 0.1;
    r(c;σ=σ)  = (c^(1-σ))/(1-σ)            # return fcn.
    μ(s,c; w=w,A=A) = w + A*s - c;         # Transition fcn.  ṡ=μ(s,c)

    dr(c;σ=σ) = c^(-σ)                     # deriv return fcn.
    FOC(v;σ=σ) = v^(-1.0/σ)                # FOC, u'^-1()

    μ_inv(s,ṡ; w=w,A=A) = w + A*s -ṡ;  # Inv Trans c=μ(s,ṡ), solve ṡ=μ(s,c) for c
    # μ_inv(s,0) = con0(s)
    #con0(s;α=α,δ=δ,A=A) = w + A*s          # solve: ṡ=μ(s,c) for c(s) if ṡ=0.

    s_min = a_bar
    μ_inv(s_min,0) > 0

    s_max = 1.0
    H = 10_000;
    s = LinRange(s_min, s_max, H)
    s = convert(Array, s)
    ds = (s_max-s_min)/(H-1)
end

##################################
Δ = 1_000
maxit = 10_000
ε = 10e-6
dVf, dVb     = [zeros(H,1) for i in 1:2]
dV_Upwind, c = [zeros(H,1) for i in 1:2]
v0 = @. r(μ_inv(s,0))/ρ #initial guess for V
v = v0
dist=[]

for n=1:maxit
    #println(n)
	V=v
    dV = (V[2:H]-V[1:H-1])/ds

    # forward difference
	dVf[1:H-1] = dV
	dVf[H]     = dr(μ_inv(s_max,0))      # u'(c(s_max)) state constraint, for stability
    #dVf[H]= 0;

	# backward difference
	dVb[2:H] = dV
	dVb[1]   = dr(μ_inv(s_min,0))        # state constraint, for stability
    #dVb[1] = (A*s_min)^(-σ)

	I_concave = dVb .> dVf
    # scatter(I_concave) #1 everywhere EXCEPT @ last point H.

    # consumption and savings with forward difference
    cf  = FOC.(dVf)
    μ_f = μ.(s, cf)

    # consumption and savings with backward difference
    cb  = FOC.(dVb)
    μ_b = μ.(s, cb)

    # c if k̇=0 & V = u'(c)
	c0  = μ_inv.(s,0)   # con0.(s)
    dV0 = dr.(c0)

    # Now to make a choice between forward and backward difference
    If = μ_f .> 0
    Ib = μ_b .< 0
    I0 = (1.0 .- If - Ib)

    dV_Upwind= dVf.*If + dVb.*Ib + dV0.*I0
    c = FOC.(dV_Upwind)
    u = r.(c)

    # plot(dVf, lab="dVf")
    # plot!(dVb, lab="dVb")
    # plot!(dV0, lab="dV0")
    # plot!(dV_Upwind, lab="dV_Upwind")

    # plot(cf, lab="dVf")
    # plot!(cb, lab="dVb")
    # plot!(c0, lab="dV0")
    # plot!(c, lab="dV_Upwind")

    # plot(μ_f, lab="dVf")
    # plot!(μ_b, lab="dVb")
    # plot!(zero(μ_f), lab="dV0")

    # create the transition matrix
    X = -min.(μ_b,0)/ds
    Y = -max.(μ_f,0)/ds + min.(μ_b,0)/ds
    Z = max.(μ_f,0)/ds

    a1 = sparse(Diagonal((Y[:])))
    a2 = [zeros(1,H); sparse(Diagonal((X[2:H]))) zeros(H-1,1)]
    a3 = [zeros(H-1,1) sparse(Diagonal((Z[1:H-1]))); zeros(1,H)]
    AA = a1 + a2 + a3

    B = (ρ + 1/Δ)*sparse(I,H,H) - AA
    b = u + V./Δ

    # Solve: B V = b
    # => (ρ + 1/Δ)*sparse(I,H,H)*V - AA*V = u + V./Δ
    # => (ρ + 1/Δ)*V = u + V./Δ +  AA*V
    # => (ρ)*V + V./Δ = u + V./Δ +  AA*V
    # => (ρ)*V = u  +  AA*V

    V = B \ b
    V_change = V-v
    v = V

	push!(dist,findmax(abs.(V_change))[1])
    println(n, " ", dist[n])
	if dist[n] .< ε
		println("Value Function Converged Iteration=")
		println(n)
		break
	end
end

s_dot = μ.(s,c)
v_err = r.(c) + dV_Upwind.*s_dot - ρ.*v # approx @ borrowing constraint
# TODO: SIMULATE!

plot(dist,
		xlabel="Iteration", ylabel="||V^{n+1} - V^n||",
		#ylims=(-0.001,0.030),
		legend=false, title="")
#png("Convergence")

plot(s, v_err,
		xlabel="s", ylabel="Error in the HJB equation",
		xlims=(s_min,s_max),
		legend=false, title="")
#png("HJB_error")

plot(s, v,
		xlabel="k", ylabel="V(k)",
		xlims=(s_min,s_max),
		legend=false, title="")
#png("Value_function_vs_k")

plot(s, c,
		xlabel="k", ylabel="c(k)",
		xlims=(s_min,s_max),
		legend=false, title="")
#png("c(k)_vs_k")

plot(s, s_dot,
		xlabel="k", ylabel="s(k)",
		xlims=(s_min,s_max), title="", label="s(k)", legend=false)
plot!(s, zeros(H,1), label="", line=:dash)
#png("stateconstraint")

# Simulate.
# s_0 = 0.10;
# c_0 # = interpolate(s_0, (s,c(s)))
# sp = μ(s_0, c)

# Special case of NGM has closed-form sol. Compare!
σ= 4.0;  δ = 0.05; A = 1.0; α= 0.3;
ρ = δ*(α*σ - 1) # restriction for closed form

# Policy: c = (1-s)(k^α)
cp(s) = (1-((α*δ)/(ρ+δ)))*s^α
s_d(s) = (s^α) - δ*s - cp(s)
plot(s, cp)
plot!(s, s_d)

#
s_ss = (α*A/(ρ+δ))^(1/(1-α))
z_ss = (s_ss)^(1-α)
c_ss = (s_ss)^(α)  - δ*s_ss

s_0 = 5;
z_0 = (s_0)^(1-α)

λ=(1-α)*δ
z(t) = z_ss + exp(-λ*t)*(z_0 - z_ss)
kk(t) = z(t)^(1/(1-α))
cc(t) = (1-((α*δ)/(ρ+δ)))*kk(t)^α

t_sim = 0.01:0.01:100
plot(legend=:topleft)
plot!(t_sim, z, lab="z")
plot!(t_sim, kk, lab="k")
plot!(t_sim, cc, lab="c")
plot!([z_ss s_ss c_ss],  seriestype = :hline, lab="", color="grey")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/BoundingSphere.jl/blob/master/LICENSE

"""
    Ritter()

## Pros
* extremly fast
* simple

## Cons
* Very inaccurate.
"""
struct Ritter <: BoundingSphereAlg end

function max_distance_point(pts, pt1)
    pt_best = first(pts)
    d2_best = sqdist(pt_best, pt1)
    for pt in pts
        d2 = sqdist(pt, pt1)
        if d2 > d2_best
            pt_best = pt
            d2_best = d2
        end
    end
    pt_best
end

function sqsphere_two_points(pt1,pt2)
    c = map(middle, pt1, pt2)
    r2 = sqdist(pt1, pt2) / 4
    c, r2
end

@noinline function ritter(pts)
    pt1 = first(pts)
    pt2 = max_distance_point(pts, pt1)
    c, r2 = sqsphere_two_points(pt1, pt2)
    for pt in pts
        if sqdist(pt, c) > r2
            direction = (c - pt) / norm(c - pt)
            r = sqrt(r2)
            pt_op = c + direction * r
            c,r2 = sqsphere_two_points(pt, pt_op)
        end
    end
    r = sqrt(r2)
    c, r
end

boundingsphere(pts, alg::Ritter) = ritter(pts)


text:
include("types.jl")
include("shared.jl")
include("util.jl")
include("adv_cg.jl")

## set

solver = :mosek
# solver = :gurobi

log = 0
psdtol = 1e-6
perturb = 1e-12
obj_reltol_cv = 0.0
obj_reltol_test = 0.0

verbose = false
### prepare data

dname = "glass"
D_all = readcsv("data-example/" * dname * ".csv")
id_train = readcsv("data-example/" * dname * ".train")
id_test = readcsv("data-example/" * dname * ".test")

id_train = round(Int64, id_train)
id_test = round(Int64, id_test)

println(dname)

### Cross Validation, using first split
## First stage

id_tr = vec(id_train[1,:])
id_ts = vec(id_test[1,:])
X_train = D_all[id_tr,1:end-1]
y_train = round(Int, D_all[id_tr, end])

X_test = D_all[id_ts,1:end-1]
y_test = round(Int, D_all[id_ts, end])

X_train, mean_vector, std_vector = standardize(X_train)
X_test = standardize(X_test, mean_vector, std_vector)

Cs =  [2.0^i for i=0:3:12]
ncs = length(Cs)

# fold
n_train = size(X_train, 1)
n_test = size(X_test, 1)
kf = 5

# k folds
folds = k_fold(n_train, kf)

loss_list = zeros(ncs)
loss01_list = zeros(ncs)

# The first stage of CV
idx = randperm(n_train)
X_train = X_train[idx,:]
y_train = y_train[idx]

for i = 1:ncs

  println(i, " | Adversarial | C = ", string(Cs[i]))

  losses = zeros(n_train)
  losses1 = zeros(n_train)
  # k fold
  for j = 1:kf
    # prepare training and validation
    id_tr = vcat(folds[[1:j-1; j+1:end]]...)
    id_val = folds[j]

    X_tr = X_train[id_tr, :];  y_tr = y_train[id_tr]
    X_val = X_train[id_val, :];  y_val = y_train[id_val]

    print("    ",j, "-th fold : ")
    @time model = train_adv_cg(X_tr, y_tr, Cs[i], perturb=perturb, obj_reltol=obj_reltol_cv, solver=solver, log=log, psdtol=psdtol, verbose=verbose)

    _, ls, _, ls1, _, _ = test_adv(model, X_val, y_val)

    losses[id_val] = ls
    losses1[id_val] = ls1

  end

  loss_list[i] = mean(losses)
  loss01_list[i] = mean(losses1)
  # println("loss : ", string(mean(losses)))
  println("  => loss01 : ", string(mean(losses1)))
  println()

end

ind_max= indmin(loss01_list)
C0 = Cs[ind_max]
Cs =  [C0*2.0^(i-3) for i=1:5]
ncs = length(Cs)

## Second stage
idx = randperm(n_train)
X_train = X_train[idx,:]
y_train = y_train[idx]

for i = 1:ncs

  println(i, " | Adversarial | C = ", string(Cs[i]))

  losses = zeros(n_train)
  losses1 = zeros(n_train)
  # k fold
  for j = 1:kf
    # prepare training and validation
    id_tr = vcat(folds[[1:j-1; j+1:end]]...)
    id_val = folds[j]

    X_tr = X_train[id_tr, :];  y_tr = y_train[id_tr]
    X_val = X_train[id_val, :];  y_val = y_train[id_val]

    print("    ",j, "-th fold : ")
    @time model = train_adv_cg(X_tr, y_tr, Cs[i], perturb=perturb, obj_reltol=obj_reltol_cv, solver=solver, log=log, psdtol=psdtol, verbose=verbose)

    _, ls, _, ls1, _, _ = test_adv(model, X_val, y_val)
    losses[id_val] = ls
    losses1[id_val] = ls1
    #println("loss : ", string(ls))
    #println("loss01 : ", string(ls))

  end

  loss_list[i] = mean(losses)
  loss01_list[i] = mean(losses1)
  # println("loss : ", string(mean(losses)))
  println("  => loss01 : ", string(mean(losses1)))
  println()

end

ind_max= indmin(loss01_list)
C_best = Cs[ind_max]

### Evaluation

n_split = size(id_train, 1)

v_model = Vector{ClassificationModel}()
v_result = Vector{Tuple}()
v_acc = zeros(n_split)
v_cs_result = zeros(n_split, 5)

for i = 1:n_split
  # standardize
  id_tr = vec(id_train[i,:])
  id_ts = vec(id_test[i,:])
  X_train = D_all[id_tr,1:end-1]
  y_train = round(Int, D_all[id_tr, end])

  X_test = D_all[id_ts,1:end-1]
  y_test = round(Int, D_all[id_ts, end])

  X_train, mean_vector, std_vector = standardize(X_train)
  X_test = standardize(X_test, mean_vector, std_vector)

  #train and test
  @time model = train_adv_cg(X_train, y_train, C_best, perturb=perturb, obj_reltol=obj_reltol_test, solver=solver, log=log, psdtol=psdtol, verbose=verbose)

  result = test_adv(model, X_test, y_test)
  loss01 = result[3]
  acc = 1.0 - loss01
  cs_result = count_constraints(model, y_train)

  println("accuracy : ", acc)

  push!(v_model, model)
  push!(v_result, result)
  v_acc[i] = acc
  v_cs_result[i, :] = collect(cs_result)
end

println(dname)
println("mean accuracy : ", mean(v_acc))
println("std accuracy : ", std(v_acc))

max_stars_repo_name: Koprvhdix/adversarial-multiclass

max_issues_repo_name: Koprvhdix/adversarial-multiclass

max_forks_repo_name: Koprvhdix/adversarial-multiclass


text:
using Wasmtime
using Test

@testset "Call foreign" begin
    engine = WasmEngine()
    store = Wasmtime.WasmtimeStore(engine)
    code = wat"""
    (module
        (func $add (param $lhs i32) (param $rhs i32) (result i32)
            local.get $lhs
            local.get $rhs
            i32.add)
        (export "add" (func $add))
    )
    """
    modu = Wasmtime.WasmtimeModule(engine, code)

    instance = Wasmtime.WasmtimeInstance(store, modu)
    add = exports(instance).add

    res = add(Int32(1), Int32(2))

    @test res == Int32(3)
end

# include("./table.jl")
include("./import_export.jl")
include("./wat2wasm.jl")
include("./wasi.jl")


text:
include("ModuleA.jl")

import MyModule

mytype = MyModule.MyType(1)

#@show foo(mytype)
#@show MyType
#@show bar

@show MyModule.foo(mytype)
@show MyModule.MyType
@show MyModule.bar(2)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using SDFResults
using Test
using Unitful
using PICDataStructures
using RecursiveArrayTools: recursive_bottom_eltype

@testset "SDFResults.jl" begin
    dir = "gauss"
    sim = read_simulation(dir)
    @test sim isa EPOCHSimulation
    file = sim[1]
    @test file isa SDFFile

    # test for different code paths in expensive grid detection
    Ex, Ey = sim[1][:ex, :ey]
    Ez = sim[1][:ez]
    @test Ex isa ScalarField{3}
    @test Ey isa ScalarField{3}
    @test Ez isa ScalarField{3}
    @test unit(eltype(Ex)) == u"V/m"
    @test unit(recursive_bottom_eltype(getdomain(Ex))) == u"m"

    vars = (:grid, Symbol("py/electron"), Symbol("pz/electron"))
    @task all(vars .∈ keys(file))
    (x,y,z), py, pz = read(file, vars...)
    @test all(unit.(x) .== u"m")
    @test all(unit.(py) .== u"kg*m/s")

    t = get_time(file)
    @test (t |> u"fs") ≈ 10u"fs" atol = 0.1u"fs"

    nx, ny, nz = get_parameter.((file,), (:nx, :ny, :nz))
    @test nx == ny == nz == 10

    λ = get_parameter(file, :laser, :lambda)
    @test (λ |> u"nm") ≈ 800u"nm"

    @test ndims(file) == 3
    @test cell_volume(file) ≠ 0
end


text:
using PlotlyJS, JLD2, ORCA

function plot_error(finame::String, type::String; show::Bool=true, foname::String="")
    f = jldopen(finame)
    dvVlist = f["info/dvV"]
    dvVlist = dvVlist[1:5:end]
    noiselist = f["info/noise"]

    errors = zeros(length(dvVlist), length(noiselist))
    # assemble matrix for heatmap
    for (i,dvV) in enumerate(dvVlist)
        for (j,noise) in enumerate(noiselist)
            errors[i,j] = f["$type/$dvV.$noise"]
        end
    end

    # set up plot
    layout = Layout(width=1200, height=800,
                    xaxis_title="dv/v",
                    yaxis_title="Noise level",
                    font=attr(size=18),
                    showlegend=true,
                    title="Stretching Error")
    p = PlotlyJS.plot([NaN], layout)
    trace = PlotlyJS.heatmap(x=dvVlist .* 100, y=noiselist .* 100, z=errors)
    addtraces!(p, trace)
    deletetraces!(p, 1)

    if show==true
        display(p)
        readline()
    end
    if foname !== ""
        savefig(p, foname)
    end
end

plot_error("verificationDataError_DTW.jld2", "realData", foname="DTWError_stretchXcorr.pdf")

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
using Pkg
using RegistryTools.Compress: compress_versions, load, compress

@testset "compress_versions()" begin
    # Test exact version matching
    vs = [v"1.1.0", v"1.1.1", v"1.1.2"]
    @test compress_versions(vs, [vs[2]]) == Pkg.Types.VersionSpec("1.1.1")

    # Test holes
    vs = [v"1.1.0", v"1.1.1", v"1.1.4"]
    @test compress_versions(vs, [vs[2]]) == Pkg.Types.VersionSpec("1.1.1")

    # Test patch variation with length(subset) > 1
    vs = [v"1.1.0", v"1.1.1", v"1.1.2", v"1.1.3", v"1.2.0"]
    @test compress_versions(vs, [vs[2], vs[3]]) == Pkg.Types.VersionSpec("1.1.1-1.1.2")

    # Test minor variation
    vs = [v"1.1.0", v"1.1.1", v"1.2.0"]
    @test compress_versions(vs, [vs[2]]) == Pkg.Types.VersionSpec("1.1.1-1.1")

    # Test major variation
    vs = [v"1.1.0", v"1.1.1", v"1.2.0", v"2.0.0"]
    @test compress_versions(vs, [vs[2], vs[3]]) == Pkg.Types.VersionSpec("1.1.1-1")

    # Test build numbers and prerelease values are ignored
    vs = [v"1.1.0-alpha", v"1.1.0+0", v"1.1.0+1"]
    @test compress_versions(vs, [vs[2]]) == Pkg.Types.VersionSpec("1")
end

@testset "Compress.load" begin
    mktempdir(@__DIR__) do temp_dir
        versions = [v"1.0.0", v"1.1.0", v"1.2.0", v"1.3.0"]
        compat_file = joinpath(temp_dir, "Compat.toml")
        good_compat = """
            ["1-1.1"]
            julia = "1"

            ["1.2-1"]
            julia = "1.3.0-1"
            """
        write(compat_file, good_compat)
        compat = load(compat_file, versions)
        @test compat[v"1.0.0"] == Dict("julia" => "1")
        @test compat[v"1.1.0"] == Dict("julia" => "1")
        @test compat[v"1.2.0"] == Dict("julia" => "1.3.0-1")
        @test compat[v"1.3.0"] == Dict("julia" => "1.3.0-1")

        # Overlapping ranges.
        bad_compat = """
            ["1-1.1.2"]
            julia = "1"

            ["1.1.0-1"]
            julia = "1.3.0-1"
            """
        write(compat_file, bad_compat)
        @test_throws ErrorException load(compat_file, versions)
    end
end

@testset "issue#46: compress with mixed input data" begin
    pkglibdl_str = Dict("Pkg" => "44cfe95a-1eb2-52ea-b672-e2afdf69b78f",
                        "Libdl" => "8f399da3-3557-5675-b5ff-fb832c97cbdb")
    pkglibdl_uuid = Dict{String, Base.UUID}(k => Base.UUID(v) for (k, v) in pkglibdl_str)
    uncompressed = Dict{VersionNumber, Dict{Any, Any}}(
        v"0.21.1+0" => pkglibdl_str,
        v"0.22.0+0" => pkglibdl_str,
        v"0.22.1+0" => pkglibdl_uuid
        )
    versions = VersionNumber[ v"0.21.1+0", v"0.22.0+0"]
    compressed = compress("/dev/null", uncompressed, versions)
    @test compressed == Dict{String,Any}("0" => pkglibdl_str)
end

import RegistryTools.Compress
import Pkg.TOML

@testset "Compress roundtripping of a registry" begin
    git = Sys.which("git")
    git === nothing && return
    mktempdir() do dir
        run(`$git clone https://github.com/JuliaRegistries/General $dir`)
        reg = TOML.parsefile(abspath(dir, "Registry.toml"))
        for (uuid, data) in reg["packages"]
            for f in abspath.("General", data["path"], ("Deps.toml", "Compat.toml"))
                if isfile(f)
                    Compress.save(f, Compress.load(f))
                end
            end
        end
        @test success(`$git -C $dir diff-index --quiet HEAD --`)
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "FINUFFTBuilder"
version = v"0.1.4"

# Collection of sources required to build FINUFFTBuilder
sources = [
    "https://github.com/flatironinstitute/finufft/archive/1.1.1.zip" =>
    "a51e8d56544378ace1e5c4177717dfe0943ea3aab08018094f2a9c4375576ed6",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd finufft-1.1.1/
# Scrub out OMP suffix, since FFTW compiled with --with-combined-threads
sed -i s/_\$\(FFTWOMPSUFFIX\)// makefile
# -march=native appears to be Intel-only
if [[ $target == aarch64-* ]] || [[ $target == arm-* ]] || [[ $target == powerpc* ]]; then
   sed -i s/-march=native// makefile
fi
# Static libgcc linking for MacOS
#if [[ ${target} == x86_64-apple* ]]; then
#    static="-static-libgcc -static-libstdc++"
#fi
# Compile
make lib/libfinufft.so LIBRARY_PATH=$prefix/lib/ CPATH=$prefix/include/ LIBS="$static -lm -L$prefix/lib $LDFLAGS"
# Copy in lib with right suffix
if [[ ${target} == x86_64-apple* ]]; then
   cp `g++ --print-file-name=libgomp.dylib` $prefix/lib/libgomp.1.dylib
   cp lib/libfinufft.so $prefix/lib/libfinufft.dylib
else
   cp lib/libfinufft.so $prefix/lib/
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libfinufft", :libfinufft)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://raw.githubusercontent.com/JuliaMath/FFTW.jl/v0.2.4/deps/build_fftw.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: ludvigak/FINUFFTBuilder

max_issues_repo_name: ludvigak/FINUFFTBuilder

max_forks_repo_name: ludvigak/FINUFFTBuilder


text:
# NOTE: We put ContextType first, above, since that's the paramater that users are more
# likely to tweak, whereas the Storage type is more of an internal detail.

# Users call this `Runtime()` function, and we generate them a `_TopLevelRuntime()`.
function Runtime{CT,ST}(ctx::CT = CT(), st::ST = ST()) where {CT,ST<:AbstractSalsaStorage}
    return _TopLevelRuntime{CT,ST}(ctx, st)
end

# By default, Runtime() use the DefaultStorage provided by Salsa in default_storage.jl, and
# does not specify any custom context.
Runtime{CT}(ctx::CT = CT()) where {CT} = Runtime{CT,DefaultStorage}(ctx, DefaultStorage())
Runtime(st) = Runtime{EmptyContext,DefaultStorage}(EmptyContext(), st)
Runtime() = Runtime(DefaultStorage())  # Equivalent to DefaultRuntime()

# This is the top-level Runtime object which users will instantiate for using Salsa. We mark
# the _TopLevelRuntime as mutable, so that we can take its pointer_from_objref(). It will
# very likely be not isbits anyway, since the users' storage objects are unlikely to be
# isbits. It's a bit annoying that this is mutable, because its contents really musn't ever
# change, but here we are.
mutable struct _TopLevelRuntime{CT,ST<:AbstractSalsaStorage} <: Runtime{CT,ST}
    # Store the user-provided context object, where they can store per-runtime meta-data
    # that does not affect Salsa invalidation / re-evaluation.
    context::CT

    # The storage is where all the tracking of state and invalidation happens.
    storage::ST
end

########## Implementation of Runtime API

context(rt::_TopLevelRuntime) = rt.context

storage(rt::_TopLevelRuntime) = rt.storage

# Top-level runtimes do not support tracing. Explicit overloads for safety here.
trace(::_TopLevelRuntime) = error("Attempted to call `trace` on a top-level runtime.")
collect_call_stack(::_TopLevelRuntime) = error("Attempted to call `collect_call_stack` on a top-level runtime.")
collect_trace(::_TopLevelRuntime) = error("Attempted to call `collect_trace` on a top-level runtime.")
destruct_trace!(::_TopLevelRuntime) = error("Attempted to call `destruct_trace!` on a top-level runtime.")

function unmemoized_input_lookup(rt::_TopLevelRuntime, key)
    return _unwrap_salsa_value(rt, _unmemoized_input_lookup_internal(rt, key))
end

# Safety overload.
function previous_output(::_TopLevelRuntime)
    error("`previous_output(rt)` may only be called from inside a derived function.")
end


text:
module DeterministicPolicyGradient
using Distributions
export DPGopts, DPGfuns, DPGstate, dpg, meshgrid, meshgrid2, get_centers_multi, quadform

include("utilities.jl")
include("dpg.jl")
include("experience_replay.jl")

end # module


text:
# Test the scripts

module TestHeatEquationExample
    mktempdir() do dir
        cd(dir) do
            include(joinpath(@__DIR__, "../docs/src/literate/heat_equation.jl"))
        end
    end
end

module TestIncompressibleElasticity
    mktempdir() do dir
        cd(dir) do
            include(joinpath(@__DIR__, "../docs/src/literate/incompressible_elasticity.jl"))
        end
    end
end

module TestHyperElasticity
    mktempdir() do dir
        cd(dir) do
            include(joinpath(@__DIR__, "../docs/src/literate/hyperelasticity.jl"))
        end
    end
end


text:
# For exporting names in imported packages to scripts that use BeNGS
# maybe remove for final version
using Reexport
using Compat
using DPMeansClustering
@reexport using NextGenSeqUtils
using Distributions

include("consensus.jl")
include("clusterpipeline.jl")


text:
import ACE

import ACE: OneParticleBasis,
            symbols,
            indexrange,
            isadmissible,
            degree,
            get_index

import ACEbase: Discrete1pBasis

import JuLIP.Potentials: ZList, SZList, z2i, i2z

export Species1PBasis

"""
and auxiliary type to write a joint codebase for
   `Species1PBasis` and `Species1PBasisCtr`
"""
abstract type AbstractSpecies1PBasis{NZ} <: Discrete1pBasis{Bool} end

"""
`Species1PBasis` : todo write docs
"""
struct Species1PBasis{NZ} <: AbstractSpecies1PBasis{NZ}
   zlist::SZList{NZ}
end

zlist(basis::AbstractSpecies1PBasis) = basis.zlist

Species1PBasis(species) = Species1PBasis(ZList(species, static=true))

Base.length(basis::AbstractSpecies1PBasis) = length(basis.zlist)

function ACE.evaluate!(A, basis::Species1PBasis, X::AbstractState)
   fill!(A, false)
   A[z2i(basis.zlist, X.mu)] = true
   return A
end

ACE.valtype(::AbstractSpecies1PBasis, args...) = Bool

symbols(::Species1PBasis) = [:μ]

indexrange(basis::Species1PBasis) = Dict( :μ => Int.(basis.zlist.list) )

isadmissible(b, basis::Species1PBasis) = (AtomicNumber(b.μ) in basis.zlist.list)

degree(b, basis::AbstractSpecies1PBasis, args...) = 0

get_index(basis::Species1PBasis, b) = z2i(basis.zlist, AtomicNumber(b.μ))

ACE.get_spec(basis::Species1PBasis, i::Integer) = (mu = basis.zlist.list[i],)
ACE.get_spec(basis::Species1PBasis) = ACE.get_spec.(basis, 1:length(basis))

Base.rand(basis::AbstractSpecies1PBasis) = rand(basis.zlist.list)

write_dict(V::Species1PBasis) =
      Dict( "__id__" => "Species1PBasis",
            "zlist" => write_dict(V.zlist))

function read_dict(::Val{:Species1PBasis}, D::Dict)
   zlist = read_dict(D["zlist"])
   return Species1PBasis(zlist)
end


text:
module prml

using Distributions
using SpecialFunctions
using Statistics
using LinearAlgebra
using Plots

export BernoulliDist
export BetaDist
export DirichletDist
export CategoricalDist
export GaussianDist
export GaussianBayesMeanDist
export GaussianBayesVarDist
export GammaDist
export MultivariateGaussianDist

export pdf
export fitting
export draw

export PolynomialFeature
export GaussianFeature
export SigmoidalFeature
export transform

export LinearRegressor
export RidgeRegressor
export BayesianRegressor
export EmpiricalBayesianRegressor
export fitting
export log_evidence
export predict
export predictSampling

export AbstractLayer
export LinearLayer
export SigmoidLayer
export TanhLayer
export ReLULayer
export forward_propagation
export activate_derivative
export backward_propagation

export AbstractCostFunction
export SigmoidCrossEntropy
export SoftmaxCrossEntropy
export SumSquareError
export GaussianMixtureError
export delta
export cost

export NeuralNetwork
export fitting
export predict

export PolynomialKernel
export RBFKernel
export GaussianProcessRegressor
export fitting
export predict

include("distribution/distribution.jl")
include("feature/feature.jl")
include("linear/linear.jl")
include("neural_network/neural_network.jl")
include("kernel/kernel.jl")

end # module

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
rm(joinpath(@__DIR__, "Project.toml"); force=true)
rm(joinpath(@__DIR__, "Manifest.toml"); force=true)

using Pkg

Pkg.activate(@__DIR__)

if ARGS[1] == "full"
    pkgs = ["MadNLPHSL","MadNLPPardiso","MadNLPMumps","MadNLPKrylov"]
elseif ARGS[1] == "basic"
    pkgs = ["MadNLPMumps","MadNLPKrylov"]
elseif ARGS[1] == "cuda"
    pkgs = ["MadNLPGPU"]
else
    error("proper argument should be given - full or basic")
end

Pkg.develop(PackageSpec(path=joinpath(@__DIR__,"..")))
Pkg.develop(PackageSpec(path=joinpath(@__DIR__,"..","lib","MadNLPTests")))
Pkg.develop.([PackageSpec(path=joinpath(@__DIR__,"..","lib",pkg)) for pkg in pkgs])
Pkg.build()

if ARGS[1] == "full" || ARGS[1] == "basic"
    Pkg.test("MadNLP", coverage=true)
end
Pkg.test.(pkgs, coverage=true)


text:
# Tools for DataFrame Processing

defaultchkfn = (ismissing, isnothing, isnan);

"""
`filtern` works based on `filter` automatically drop rows where the corresponding specific column `col1` contains `NaN`, `missing` or `nothing`.

`filtern(col1, df)` is the shorthand of `filter(cols => fun, df::AbstractDataFrame)` where `fun` is for examing whether the column contains `NaN`, `missing` or `nothing`.
"""
function filtern!(col, df; chkfun=defaultchkfn)
    # see this thread: https://stackoverflow.com/questions/62789334/how-to-remove-drop-rows-of-nothing-and-nan-in-julia-dataframe
    return filter!(col => c1 -> !any(f -> f(c1), chkfun), df)
end

"""
`filtern(col1, col2, df)` automatically drop rows where either of `col1` and `col2`  contains `NaN`, `missing` or `nothing`.
"""
function filtern!(col1, col2, df; chkfun=defaultchkfn)
    return filter!([col1, col2] => (c1, c2) -> !any(f -> f(c1) || f(c2), chkfun), df)
end

"""
`(df_seg, col_indicator) = forsegplot(df, [:col1, :col2])` returns `df_seg` (a copy of dataframe `df`) and the column indicators (a named tuple), for `Gadfly.plot(df_seg, Geom.segment; col_indicator...)` plot in Gadfly.
Basically,
```julia
(df_seg, col_indicator) = forsegplot(df, [:col1, :col2]);
Gadfly.plot(df_seg, Geom.segment; col_indicator...);
```

is equivalently

```julia
x = df[:col1];
y = df[:col2];
df_seg = DataFrame(x0=x[1:end-1], x1=x[2:end], # segment's starts
                   y0=y[1:end-1], y1=y[2:end], # segment's ends

Gadfly.plot(df_seg, x=:x0,y=:y0,xend=:x1, yend=:y1, Geom.segment)
```

`forsegplot(df, [:col1, :col2, :col3])` does the a similar job but returns a dataframe with columns named `:x0`, `:y0`, `:z0`; higher dimension is not supported yet.
"""
function forsegplot(df::DataFrame, cols::Vector{<:Union{Symbol, String}}) #
    names0 = [:x0, :y0, :z0];
    names1 = [:x1, :y1, :z1];

    col_indicator = (x=:x0, y=:y0, xend=:x1, yend=:y1);
    if length(cols) > 2
        col_indicator = merge(col_indicator, (z=:z0, zend=:z1));
    end

    lenc = length(cols);
    names0 = names0[1:lenc];
    names1 = names1[1:lenc];

    values0 = [];
    values1 = [];

    for name in cols
        xi = df[!,name];
        push!(values0, xi[1:end-1]);
        push!(values1, xi[2:end]);
    end

    cb(v1, v2) = [permutedims(hcat(v1, v2))...]; # making a n by 2 matrix using v1 and v2, and transpose it and release every element that we can have something like :x0, :x1, :y0, :y1...

    allnames = (cb(names0, names1)...,); # must be a tuple for NamedTuple construction
    allvalues = cb(values0, values1);
    nt = NamedTuple{allnames}(allvalues);

    df_seg = DataFrame(nt);

    return df_seg, col_indicator

end


text:
@testset "Distance Fields" begin

@testset "SignedDistanceField" begin
    s1 = SignedDistanceField(HyperRectangle(Vec(1,2),Vec(3,4)), [1 2;3 4])
    @test typeof(s1) == GeometryTypes.SignedDistanceField{2,Int,Int}
    s1 = SignedDistanceField(HyperRectangle(Vec(1,2),Vec(3,4)), [1. 2.;3. 4.])
    @test typeof(s1) == GeometryTypes.SignedDistanceField{2,Int,Float64}

    # functional
    s2 = SignedDistanceField(HyperRectangle(Vec(0,0,0),Vec(1,1,1))) do v
        sqrt(sum(v.*v)) - 1 # sphere
    end
    @test size(s2) == (11, 11, 11)
    # functional
    s2 = SignedDistanceField(HyperRectangle(Vec(-1,-1),Vec(2,3))) do v
        sqrt(sum(v.*v)) - 1 # circle
    end
    @test size(s2) == (21, 31)
    @test size(s2, 1) == 21
    @test size(s2, 2) == 31

    @test HyperRectangle(s2) == HyperRectangle(Vec(-1,-1),Vec(2,3))
end

@testset "getindex" begin
    sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1),Vec(2,2))) do v
        sqrt(sum(v.*v)) - 1 # circle
    end

    @test sdf[1,1] == sdf[21*21] # by circle symmetry
    @test_throws BoundsError sdf[21*21+1]
end

end


text:
# Wrapper around the spglib library

import spglib_jll: libsymspg

mutable struct SpglibDataset
    spacegroup_number::Cint
    hall_number::Cint
    _international_symbol::NTuple{11,Cchar}
    _hall_symbol::NTuple{17,Cchar}
    _choice::NTuple{6,Cchar}
    _transformation_matrix::NTuple{9,Cdouble}
    _origin_shift::NTuple{3,Cdouble}
    n_operations::Cint
    _rotations::Ptr{Cint}
    _translations::Ptr{Cdouble}
    n_atoms::Cint
    _wyckoffs::Ptr{Cint}
    _site_symmetry_symbols::Ptr{NTuple{7,Cchar}}
    _equivalent_atoms::Ptr{Cint}
    _crystallographic_orbits::Ptr{Cint}
    _primitive_lattice::NTuple{9,Cdouble}
    _mapping_to_primitive::Ptr{Cint}
    n_std_atoms::Cint
    _std_lattice::NTuple{9,Cdouble}
    _std_types::Ptr{Cint}
    _std_positions::Ptr{Cdouble}
    _std_rotation_matrix::NTuple{9,Cdouble}
    _std_mapping_to_primitive::Ptr{Cint}
    _pointgroup_symbol::NTuple{6,Cchar}
end
function Base.getproperty(ds::SpglibDataset, name::Symbol)
    if name === :international_symbol
        letters = Char.(getfield(ds, :_international_symbol))
        fzero = findfirst(==('\0'), collect(letters))
        return join(letters[i] for i in 1:(isnothing(fzero) ? length(letters) : (fzero-1)))
    elseif name === :hall_symbol
        letters = Char.(getfield(ds, :_hall_symbol))
        fzero = findfirst(==('\0'), collect(letters))
        return join(letters[i] for i in 1:(isnothing(fzero) ? length(letters) : (fzero-1)))
    elseif name === :choice
        letters = Char.(getfield(ds, :_choice))
        fzero = findfirst(==('\0'), collect(letters))
        return join(letters[i] for i in 1:(isnothing(fzero) ? length(letters) : (fzero-1)))
    elseif name === :transformation_matrix
        return reshape(collect(getfield(ds, :_transformation_matrix)), (3,3))
    elseif name === :origin_shift
        return reshape(collect(getfield(ds, :_origin_shift)), 3)
    elseif name === :rotations
        return unsafe_wrap(Array, getfield(ds, :_rotations), (3,3,Int(getfield(ds, :n_operations))))
    elseif name === :translations
        return unsafe_wrap(Array, getfield(ds, :_translations), (3,Int(getfield(ds, :n_operations))))
    elseif name === :wyckoffs
        return unsafe_wrap(Array, getfield(ds, :_wyckoffs), Int(getfield(ds, :n_atoms)))
    elseif name === :site_symmetry_symbols
        a = unsafe_wrap(Array, getfield(ds, :_site_symmetry_symbols), getfield(ds, :n_atoms))
        ret = Vector{String}(undef, length(a))
        for i in 1:length(a)
            letters = Char.(a[i])
            fzero = findfirst(==('\0'), collect(letters))
            ret[i] = join(letters[j] for j in 1:(isnothing(fzero) ? length(letters) : (fzero-1)))
        end
        return ret
    elseif name === :equivalent_atoms
        return unsafe_wrap(Array, getfield(ds, :_equivalent_atoms), getfield(ds, :n_atoms))
    elseif name === :crystallographic_orbits
        return unsafe_wrap(Array, getfield(ds, :_crystallographic_orbits), getfield(ds, :n_atoms))
    elseif name === :primitive_lattice
        return reshape(collect(getfield(ds, :_primitive_lattice)), (3,3))
    elseif name === :mapping_to_primitive
        return unsafe_wrap(Array, getfield(ds, :_mapping_to_primitive), getfield(ds, :n_atoms))
    elseif name === :std_lattice
        return reshape(collect(getfield(ds, :_std_lattice)), (3,3))
    elseif name === :std_types
        return unsafe_wrap(Array, getfield(ds, :_std_types), getfield(ds, :n_std_atoms))
    elseif name === :std_positions
        return unsafe_wrap(Array, getfield(ds, :_std_positions), (3,Int(getfield(ds, :n_std_atoms))))
    elseif name === :std_rotation_matrix
        return reshape(collect(getfield(ds, :_std_rotation_matrix)), (3,3))
    elseif name === :std_mapping_to_primitive
        return unsafe_wrap(Array, getfield(ds, :_std_mapping_to_primitive), getfield(ds, :n_std_atoms))
    elseif name === :pointgroup_symbol
        letters = Char.(getfield(ds, :_pointgroup_symbol))
        fzero = findfirst(==('\0'), collect(letters))
        return join(letters[i] for i in 1:(isnothing(fzero) ? length(letters) : (fzero-1)))
    else
        getfield(ds, name)
    end
end

function get_spglib_dataset(net::CrystalNet)
    lattice = Matrix{Cdouble}(net.cell.mat') # transpose to account for row-major operations
    n = nv(net.graph)
    positions = Matrix{Cdouble}(undef, 3, n)
    types = Vector{Cint}(undef, n)
    symb_to_int = Dict{Symbol,Cint}()
    j = 1
    for i in 1:n
        positions[:,i] .= net.cell.mat * net.pos[i]
        j += ((types[i] = get!(symb_to_int, net.types[i], j)) == j)
    end
    ptr = ccall((:spg_get_dataset, libsymspg), Ptr{SpglibDataset},
                (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint, Cdouble),
                lattice, positions, types, n, 10*eps(Cdouble))
    dataset = unsafe_load(ptr)
    @assert dataset.n_atoms == n # otherwise the net is not minimal
    return dataset
end

function find_symmetries(net::CrystalNet{Rational{S}}) where S
    T = soft_widen(S)
    U = widen(T)
    lattice = Matrix{Cdouble}(LinearAlgebra.I, 3, 3) # positions are expressed in this basis

    I = sortperm(net.pos)
    uniquepos = SVector{3, Rational{S}}[net.pos[I[1]]]
    symb_to_int = Dict{Vector{Tuple{Symbol,Int}}, Cint}([(net.types[I[1]], 1)] => 1)
    last_types = [(net.types[I[1]], 1)]
    types = Cint[1]
    new_symbol = 2
    for i in 2:length(net.pos)
        j = I[i]
        if net.pos[j] == last(uniquepos)
            flag = true
            for k in 1:length(last_types)
                thistype, x = last_types[k]
                if thistype == net.types[j]
                    last_types[k] = (thistype, x+1)
                    flag = false
                    break
                end
            end
            if flag
                push!(last_types, (net.types[j], 1))
            end
            symb = get!(symb_to_int, sort!(last_types), new_symbol)
            types[end] = symb
            new_symbol += (symb == new_symbol)
        else
            push!(uniquepos, net.pos[j])
            last_types = [(net.types[j], 1)]
            symb = get!(symb_to_int, sort!(last_types), new_symbol)
            push!(types, symb)
            new_symbol += (symb == new_symbol)
        end
    end

    n = length(uniquepos)
    positions = Matrix{Cdouble}(undef, 3, n)
    den = 1
    for i in 1:n
        den = lcm(den, lcm((U∘denominator).(uniquepos[i])))
        positions[:,i] .= uniquepos[i]
    end

    len = ccall((:spg_get_multiplicity, libsymspg), Cint,
                (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint, Cdouble),
                lattice, positions, types, n, 100*eps(Cdouble))
    @assert len >= 1
    rotations = Array{Cint}(undef, 3, 3, len)
    translations = Array{Cdouble}(undef, 3, len)
    _len = ccall((:spg_get_symmetry, libsymspg), Cint,
                 (Ptr{Cint}, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint, Cdouble),
                 rotations, translations, len, lattice, positions, types, n, 100*eps(Cdouble))
    @assert _len == len
    # @show len
    den = lcm(len, den)
    # The first symmetry should always be the identity, which we can skip
    @assert isone(rotations[:,:,1])
    @assert iszero(translations[:,1])
    symmetries = SMatrix{3,3,Int,9}[rotations[:,:,1]]
    vmaps = Vector{Int}[collect(1:nv(net.graph))]
    floatpos = [float(x) for x in net.pos]
    hasmirror = false # whether a mirror plane exists or not. If so, the graph is achiral
    for i in 2:len
        rot = SMatrix{3,3,T,9}(transpose(rotations[:,:,i]))
        if det(rot) < 0
            @assert det(rot) == -1
            hasmirror = true
            # continue # No need to keep it as its mirror image will be kept
        end
        # @assert det(rot) == 1
        tr = SVector{3,Cdouble}(translations[:,i])
        trans = SVector{3,Rational{U}}(round.(U, den .* tr) .// den)
        vmap = check_valid_translation(net, trans, rot)
        if isnothing(vmap)
            trans = SVector{3,Rational{U}}(net.pos[last(findmin([norm(x .- tr) for x in floatpos]))])
            vmap = check_valid_translation(net, trans, rot)
            # if isnothing(vmap)
            #     @show translations[:,i]
            #     @show float(trans)
            #     @show trans
            #     @show rot
            #     println("net = CrystalNet(PeriodicGraph3D(\"", net.graph, "\"))")
            #     throw("Invalid symmetry") # This error is for debugging only, false positives without it being an error
            # end
            isnothing(vmap) && continue
        end
        @assert rot != LinearAlgebra.I # Otherwise the symmetry is a pure translation, so the net is not minimal
        push!(symmetries, rot)
        push!(vmaps, vmap::Vector{Int})
    end
    # @show vmaps
    # @show length(symmetries)
    return (symmetries, vmaps, hasmirror)
end