text:
"""
Spectral representation related functions
"""
module Spectral
export kernelT, kernelΩ, density, freq2Tau, freq2MatFreq
export kernelFermiT, kernelFermiΩ, kernelBoseT, kernelBoseΩ, fermiDirac, boseEinstein

"""
    kernelT(isFermi::Bool, symmetry::Symbol, τ::T, ω::T, β::T, regularized::Bool = false) where {T<:AbstractFloat}

Compute the imaginary-time kernel of different type.

# Arguments
- `isFermi`: fermionic or bosonic
- `symmetry`: :ph, :pha, or :none
- `τ`: the imaginary time, must be (-β, β]
- `ω`: energy
- `β`: the inverse temperature
- `regularized`: use regularized kernel or not
"""
@inline function kernelT(isFermi::Bool, symmetry::Symbol, τ::T, ω::T, β::T, regularized::Bool = false) where {T<:AbstractFloat}
    if symmetry == :none
        if regularized
            return isFermi ? kernelFermiT(τ, ω, β) : kernelBoseT_regular(τ, ω, β)
        else
            return isFermi ? kernelFermiT(τ, ω, β) : kernelBoseT(τ, ω, β)
        end
    elseif symmetry == :ph
        return isFermi ? kernelFermiT_PH(τ, ω, β) : kernelBoseT_PH(τ, ω, β)
    elseif symmetry == :pha
        return isFermi ? kernelFermiT_PHA(τ, ω, β) : kernelBoseT_PHA(τ, ω, β)
    else
        @error "Symmetry $symmetry is not implemented!"
    end
end
"""
    kernelT(isFermi::Bool, symmetry::Symbol, τGrid::AbstractVector{T}, ωGrid::AbstractVector{T}, β::T, regularized::Bool = false) where {T<:AbstractFloat}

Compute kernel with given τ and ω grids.
"""
function kernelT(isFermi::Bool, symmetry::Symbol, τGrid::AbstractVector{T}, ωGrid::AbstractVector{T}, β::T, regularized::Bool = false) where {T<:AbstractFloat}
    kernel = zeros(T, (length(τGrid), length(ωGrid)))
    for (τi, τ) in enumerate(τGrid)
        for (ωi, ω) in enumerate(ωGrid)
            kernel[τi, ωi] = kernelT(isFermi, symmetry, τ, ω, β, regularized)
        end
    end
    return kernel
end

"""
    kernelFermiT(τ, ω, β)

Compute the imaginary-time fermionic kernel.  Machine accuracy ~eps(g) is guaranteed``
```math
g(τ>0) = e^{-ωτ}/(1+e^{-ωβ}), g(τ≤0) = -e^{-ωτ}/(1+e^{ωβ})
```

# Arguments
- `τ`: the imaginary time, must be (-β, β]
- `ω`: frequency
- `β`: the inverse temperature
"""
@inline function kernelFermiT(τ::T, ω::T, β::T) where {T<:AbstractFloat}
    (-β < τ <= β) || error("τ=$τ must be (-β, β] where β=$β")
    # if τ == T(0.0)
    #     τ = -eps(T)
    # end
    if τ >= T(0.0)
        if ω > T(0.0)
            return exp(-ω * τ) / (1 + exp(-ω * β))
        else
            return exp(ω * (β - τ)) / (1 + exp(ω * β))
        end
    else
        if ω > T(0.0)
            return -exp(-ω * (τ + β)) / (1 + exp(-ω * β))
        else
            return -exp(-ω * τ) / (1 + exp(ω * β))
        end
    end
end

"""
    kernelBoseT(τ, ω, β)

Compute the imaginary-time bosonic kernel. Machine accuracy ~eps(g) is guaranteed``
```math
g(τ>0) = e^{-ωτ}/(1-e^{-ωβ}), g(τ≤0) = -e^{-ωτ}/(1-e^{ωβ})
```

# Arguments
- `τ`: the imaginary time, must be (-β, β]
- `ω`: frequency
- `β`: the inverse temperature
"""
@inline function kernelBoseT(τ::T, ω::T, β::T) where {T<:AbstractFloat}
    (-β < τ <= β) || error("τ must be (-β, β]")
    # if τ == T(0.0)
    #     τ = -eps(T)
    # end

    if τ >= T(0.0)
        if ω > T(0.0)
            # expm1(x)=exp(x)-1 fixes the accuracy for x-->0^+
            return exp(-ω * τ) / (-expm1(-ω * β))
        else
            return exp(ω * (β - τ)) / expm1(ω * β)
        end
    else
        if ω > T(0.0)
            return exp(-ω * (τ + β)) / (-expm1(-ω * β))
        else
            return exp(-ω * τ) / expm1(ω * β)
        end
    end
end

"""
    kernelBoseT_regular(τ, ω, β)

Compute the imaginary-time bosonic kernel with a regulator near ω=0. Machine accuracy ~eps(g) is guaranteed``
```math
g(τ>0) = e^{-ωτ}/(1+e^{-ωβ}), g(τ≤0) = e^{-ωτ}/(1+e^{ωβ})
```

# Arguments
- `τ`: the imaginary time, must be (-β, β]
- `ω`: frequency
- `β`: the inverse temperature
"""
@inline function kernelBoseT_regular(τ::T, ω::T, β::T) where {T<:AbstractFloat}
    (-β < τ <= β) || error("τ must be (-β, β]")
    # if τ == T(0.0)
    #     τ = -eps(T)
    # end

    if τ >= T(0.0)
        if ω > T(0.0)
            return exp(-ω * τ) / (1 + exp(-ω * β))
        else
            return exp(ω * (β - τ)) / (1 + exp(ω * β))
        end
    else
        if ω > T(0.0)
            return exp(-ω * (τ + β)) / (1 + exp(-ω * β))
        else
            return exp(-ω * τ) / (1 + exp(ω * β))
        end
    end
end

"""
    kernelFermiT_PH(τ, ω, β)

Compute the imaginary-time kernel for correlation function ``⟨O(τ)O(0)⟩``. Machine accuracy ~eps(C) is guaranteed``
```math
K(τ) = e^{-ω|τ|}+e^{-ω(β-|τ|)}
```

# Arguments
- `τ`: the imaginary time, must be (-β, β]
- `ω`: frequency, ω>=0
- `β`: the inverse temperature
"""
@inline function kernelFermiT_PH(τ::T, ω::T, β = T(1)) where {T<:AbstractFloat}
    (-β < τ <= β) || error("τ must be (0, β]")
    (ω >= 0) || error("ω must be >=0")
    τ = abs(τ)
    return exp(-ω * τ) + exp(-ω * (β - τ))
end

"""
    kernelBoseT_PH(τ, ω, β)

Compute the imaginary-time kernel for correlation function ``⟨O(τ)O(0)⟩``. Machine accuracy ~eps(C) is guaranteed``
```math
K(τ) = e^{-ω|τ|}+e^{-ω(β-|τ|)}
```

# Arguments
- `τ`: the imaginary time, must be (-β, β]
- `ω`: frequency, ω>=0
- `β`: the inverse temperature
"""
@inline function kernelBoseT_PH(τ::T, ω::T, β = T(1)) where {T<:AbstractFloat}
    (-β < τ <= β) || error("τ must be (0, β]")
    (ω >= 0) || error("ω must be >=0")
    τ = abs(τ)
    return exp(-ω * τ) + exp(-ω * (β - τ))
end

"""
    kernelFermiT_PHA(τ, ω, β)

Compute the imaginary-time kernel for correlation function ``⟨O(τ)O(0)⟩``. Machine accuracy ~eps(C) is guaranteed``
```math
K(τ) = e^{-ω|τ|}-e^{-ω(β-|τ|)}
```

# Arguments
- `τ`: the imaginary time, must be (0, β]
- `ω`: frequency, ω>=0
- `β`: the inverse temperature
"""
@inline function kernelFermiT_PHA(τ::T, ω::T, β::T) where {T<:AbstractFloat}
    (-β < τ <= β) || error("τ must be (-β, β]")
    (ω >= 0) || error("ω must be >=0")
    τ = abs(τ)
    return exp(-ω * τ) - exp(-ω * (β - τ))
end

"""
    kernelBoseT_PHA(τ, ω, β)

Compute the imaginary-time kernel for correlation function ``⟨O(τ)O(0)⟩``. Machine accuracy ~eps(C) is guaranteed``
```math
K(τ) = e^{-ω|τ|}-e^{-ω(β-|τ|)}
```

# Arguments
- `τ`: the imaginary time, must be (0, β]
- `ω`: frequency, ω>=0
- `β`: the inverse temperature
"""
@inline function kernelBoseT_PHA(τ::T, ω::T, β::T) where {T<:AbstractFloat}
    (-β < τ <= β) || error("τ must be (-β, β]")
    (ω >= 0) || error("ω must be >=0")
    τ = abs(τ)
    return exp(-ω * τ) - exp(-ω * (β - τ))
end

"""
    kernelΩ(isFermi::Bool, symmetry::Symbol, n::Int, ω::T, β::T, regularized::Bool = false) where {T<:AbstractFloat}

Compute the imaginary-time kernel of different type. Assume ``k_B T/\\hbar=1``

# Arguments
- `isFermi`: fermionic or bosonic
- `symmetry`: :ph, :pha, or :none
- `n`: index of the Matsubara frequency
- `ω`: energy
- `β`: the inverse temperature
- `regularized`: use regularized kernel or not
"""
@inline function kernelΩ(isFermi::Bool, symmetry::Symbol, n::Int, ω::T, β::T, regularized::Bool = false) where {T<:AbstractFloat}
    if symmetry == :none
        if regularized
            return isFermi ? kernelFermiΩ(n, ω, β) : kernelBoseΩ_regular(n, ω, β)
        else
            return isFermi ? kernelFermiΩ(n, ω, β) : kernelBoseΩ(n, ω, β)
        end
    elseif symmetry == :ph
        return isFermi ? kernelFermiΩ_PH(n, ω, β) : kernelBoseΩ_PH(n, ω, β)
    elseif symmetry == :pha
        return isFermi ? kernelFermiΩ_PHA(n, ω, β) : kernelBoseΩ_PHA(n, ω, β)
    else
        @error "Symmetry $symmetry  is not implemented!"
    end
end

"""
    kernelΩ(isFermi::Bool, symmetry::Symbol, nGrid::Vector{Int}, ωGrid::Vector{T}, β::T) where {T<:AbstractFloat}

Compute kernel matrix with given ωn (integer!) and ω grids.
"""
function kernelΩ(isFermi::Bool, symmetry::Symbol, nGrid::Vector{Int}, ωGrid::Vector{T}, β::T, regularized::Bool = false) where {T<:AbstractFloat}
    kernel = zeros(Complex{T}, (length(nGrid), length(ωGrid)))
    for (ni, n) in enumerate(nGrid)
        for (ωi, ω) in enumerate(ωGrid)
            kernel[ni, ωi] = kernelΩ(isFermi, symmetry, n, ω, β, regularized)
        end
    end
    return kernel
end

"""
    kernelFermiΩ(n::Int, ω::T, β::T) where {T <: AbstractFloat}

Compute the fermionic kernel with Matsubara frequency.
```math
g(iω_n) = -1/(iω_n-ω),
```
where ``ω_n=(2n+1)π/β``. The convention here is consist with the book "Quantum Many-particle Systems" by J. Negele and H. Orland, Page 95

# Arguments
- `n`: index of the Matsubara frequency
- `ω`: energy
- `β`: the inverse temperature
"""
@inline function kernelFermiΩ(n::Int, ω::T, β::T) where {T<:AbstractFloat}
    # fermionic Matsurbara frequency
    ω_n = (2 * n + 1) * π / β
    G = -1.0 / (ω_n * im - ω)
    return Complex{T}(G)
end

"""
    kernelBoseΩ(n::Int, ω::T, β::T) where {T <: AbstractFloat}

Compute the bosonic kernel with Matsubara frequency.
```math
g(iω_n) = -1/(iω_n-ω),
```
where ``ω_n=2nπ/β``. The convention here is consist with the book "Quantum Many-particle Systems" by J. Negele and H. Orland, Page 95

# Arguments
- `n`: index of the Matsubara frequency
- `ω`: energy
- `β`: the inverse temperature
"""
@inline function kernelBoseΩ(n::Int, ω::T, β::T) where {T<:AbstractFloat}
    # fermionic Matsurbara frequency
    ω_n = (2 * n) * π / β
    G = -1.0 / (ω_n * im - ω)
    if !isfinite(G)
        throw(DomainError(-1, "Got $G for the parameter $n, $ω and $β"))
    end
    return Complex{T}(G)
end

"""
    kernelBoseΩ_regular(n::Int, ω::T, β::T) where {T <: AbstractFloat}

Compute the bosonic kernel in Matsubara frequency with a regulartor near ω=0
```math
g(iω_n) = -(1-e^{-ωβ})/(1+e^{-ωβ})/(iω_n-ω),
```
where ``ω_n=2nπ/β``. The convention here is consist with the book "Quantum Many-particle Systems" by J. Negele and H. Orland, Page 95

# Arguments
- `n`: index of the Matsubara frequency
- `ω`: energy
- `β`: the inverse temperature
"""
@inline function kernelBoseΩ_regular(n::Int, ω::T, β::T) where {T<:AbstractFloat}
    # fermionic Matsurbara frequency
    ω_n = (2 * n) * π / β
    x = ω * β

    if n == 0 && abs(x) < 1.0e-5
        G = β * (1 - x / 2 + x^2 / 6) / (1 + exp(-x)) #β*(1-e^{-x})/x
    else
        # expm1(x)=exp(x)-1 fixes the accuracy for x-->0^+
        if ω > T(0.0)
            G = -1.0 / (ω_n * im - ω) * (-expm1(-x)) / (1 + exp(-x))
        else
            G = -1.0 / (ω_n * im - ω) * expm1(x) / (exp(x) + 1)
        end
    end
    if !isfinite(G)
        throw(DomainError(-1, "Got $G for the parameter $n, $ω and $β"))
    end
    return Complex{T}(G)
end

"""
    kernelFermiΩ_PH(n::Int, ω::T, β::T) where {T <: AbstractFloat}

Compute the Matsubara-frequency kernel for a correlator ``⟨O(τ)O(0)⟩_{iω_n}``.
```math
K(iω_n) = -\\frac{2iω_n}{ω^2+ω_n^2}(1+e^{-ωβ}),
```
where ``ω_n=(2n+1)π/β``. The convention here is consist with the book "Quantum Many-particle Systems" by J. Negele and H. Orland, Page 95

# Arguments
- `n`: index of the Matsubara frequency
- `ω`: energy
- `β`: the inverse temperature
"""
@inline function kernelFermiΩ_PH(n::Int, ω::T, β::T) where {T<:AbstractFloat}
    # Matsurbara-frequency correlator
    if ω < T(0.0)
        throw(DomainError("real frequency should be positive!"))
    end
    ω_n = (2n + 1) * π / β
    K = -2ω_n / (ω^2 + ω_n^2) * (1 + exp(-ω * β))
    if !isfinite(K)
        throw(DomainError(-1, "Got $K for the parameter $n, $ω and $β"))
    end
    return Complex{T}(T(0), K) #purely imaginary!
end

"""
    kernelBoseΩ_PH(n::Int, ω::T, β::T) where {T <: AbstractFloat}

Compute the Matsubara-frequency kernel for a correlator ``⟨O(τ)O(0)⟩_{iω_n}``.
```math
K(iω_n) = \\frac{2ω}{ω^2+ω_n^2}(1-e^{-ωβ}),
```
where ``ω_n=2nπ/β``. The convention here is consist with the book "Quantum Many-particle Systems" by J. Negele and H. Orland, Page 95

# Arguments
- `n`: index of the Matsubara frequency
- `ω`: energy
- `β`: the inverse temperature
"""
@inline function kernelBoseΩ_PH(n::Int, ω::T, β::T) where {T<:AbstractFloat}
    # Matsurbara-frequency correlator
    if ω < T(0.0)
        throw(DomainError("real frequency should be positive!"))
    end
    x = ω * β
    if n == 0 && x < 1.0e-5
        K = β * (2 - x + x^2 / 3) #2β*(1-e^{-x})/x
    else
        ω_n = 2n * π / β
        # expm1(x)=exp(x)-1 fixes the accuracy for x-->0^+
        K = 2ω / (ω^2 + ω_n^2) * (-expm1(-x))
    end
    if !isfinite(K)
        throw(DomainError(-1, "Got $K for the parameter $n, $ω and $β"))
    end
    return K
end

"""
    kernelFermiΩ_PHA(n::Int, ω::T, β::T) where {T <: AbstractFloat}

Compute the Matsubara-frequency kernel for a anormalous fermionic correlator with particle-hole symmetry.
```math
K(iω_n) = \\frac{2ω}{ω^2+ω_n^2}(1+e^{-ωβ}),
```
where ``ω_n=(2n+1)π/β``. The convention here is consist with the book "Quantum Many-particle Systems" by J. Negele and H. Orland, Page 95

# Arguments
- `n`: index of the fermionic Matsubara frequency
- `ω`: energy
- `β`: the inverse temperature
"""
@inline function kernelFermiΩ_PHA(n::Int, ω::T, β::T) where {T<:AbstractFloat}
    # Matsurbara-frequency correlator
    if ω < T(0.0)
        throw(DomainError("real frequency should be positive!"))
    end
    ω_n = (2n + 1) * π / β
    K = 2ω / (ω^2 + ω_n^2) * (1 + exp(-ω * β))
    if !isfinite(K)
        throw(DomainError(-1, "Got $K for the parameter $n, $ω and $β"))
    end
    return K
end

"""
    kernelBoseΩ_PHA(n::Int, ω::T, β::T) where {T <: AbstractFloat}

Compute the Matsubara-frequency kernel for a anormalous fermionic correlator with particle-hole symmetry.
```math
K(iω_n) = -\\frac{2iω_n}{ω^2+ω_n^2}(1-e^{-ωβ}),
```
where ``ω_n=2nπ/β``. The convention here is consist with the book "Quantum Many-particle Systems" by J. Negele and H. Orland, Page 95
# Arguments
- `n`: index of the fermionic Matsubara frequency
- `ω`: energy
- `β`: the inverse temperature
"""
@inline function kernelBoseΩ_PHA(n::Int, ω::T, β::T) where {T<:AbstractFloat}
    # Matsurbara-frequency correlator
    if ω < T(0.0)
        throw(DomainError("real frequency should be positive!"))
    end
    x = ω * β
    ω_n = 2n * π / β
    # expm1(x)=exp(x)-1 fixes the accuracy for x-->0^+
    K = -2ω_n / (ω^2 + ω_n^2) * (1 + exp(-ω * β))

    if n == 0
        K = T(0)
    else
        ω_n = 2n * π / β
        # expm1(x)=exp(x)-1 fixes the accuracy for x-->0^+
        K = -2ω_n / (ω^2 + ω_n^2) * (-expm1(-x))
    end
    if !isfinite(K)
        throw(DomainError(-1, "Got $K for the parameter $n, $ω and $β"))
    end
    return Complex{T}(T(0), K)
end

"""
    density(isFermi::Bool, ω, β)

Compute the imaginary-time kernel of different type. Assume ``k_B T/\\hbar=1``

# Arguments
- `isFermi`: fermionic or bosonic
- `ω`: energy
- `β`: the inverse temperature
"""
@inline function density(isFermi::Bool, ω::T, β::T) where {T<:AbstractFloat}
    return isFermi ? fermiDirac(ω, β) : boseEinstein(ω, β)
end

"""
fermiDirac(ω)

Compute the Fermi Dirac function. Assume ``k_B T/\\hbar=1``
```math
f(ω) = 1/(e^{ωβ}+1)
```

# Arguments
- `ω`: frequency
- `β`: the inverse temperature
"""
@inline function fermiDirac(ω::T, β::T) where {T<:AbstractFloat}
    x = ω * β
    if -T(50.0) < x < T(50.0)
        return 1.0 / (1.0 + exp(x))
    elseif x >= T(50.0)
        return exp(-x)
    else # x<=-50.0
        return 1.0 - exp(x)
    end
end

"""
boseEinstein(ω)

Compute the Fermi Dirac function. Assume ``k_B T/\\hbar=1``
```math
f(ω) = 1/(e^{ωβ}-1)
```

# Arguments
- `ω`: frequency
- `β`: the inverse temperature
"""
@inline function boseEinstein(ω::T, β::T) where {T<:AbstractFloat}
    # if -eps(T)<ω<eps(T)
    #     return 0.0
    # end
    n = 0.0
    x = ω * β
    if -T(50.0) < x < T(50.0)
        n = 1.0 / (exp(x) - 1.0)
    elseif x >= T(50.0)
        n = exp(-x)
    else # x<=-50.0
        n = -1.0 - exp(x)
    end
    if !isfinite(n)
        throw(DomainError(-1, "Got $n for the parameter $ω and $β"))
    end
    return n
end

end


text:
using BinaryBuilder
using SHA

# Define what we're downloading, where we're putting it
src_name = "Rmath-julia"
src_vers = "0.2.0"
src_url = "https://github.com/JuliaLang/Rmath-julia/archive/v$(src_vers).tar.gz"
src_hash = "087ada2913c5401c5772cde1606f9924dcb159f1c9d755630dcce350ef8036ac"

# First, download the source, store it in ./downloads/
src_path = joinpath(pwd(), "downloads", basename(src_url))
try mkpath(dirname(src_path)) end
download_verify(src_url, src_hash, src_path; verbose=true)

# Our build products will go into ./products
out_path = joinpath(pwd(), "products")
rm(out_path; force=true, recursive=true)
mkpath(out_path)

# Build Rmath-julia for all our platforms
products = Dict()
for platform in supported_platforms()
    target = platform_triplet(platform)

    # We build in a platform-specific directory
    build_path = joinpath(pwd(), "build", target)
    try mkpath(build_path) end

    cd(build_path) do
        # For each build, create a temporary prefix we'll install into
        temp_prefix() do prefix
            # Unpack the source into our build directory
            unpack(src_path, build_path; verbose=true)

            cd("./Rmath-julia-$(src_vers)") do
                # We expect this output from our build steps
                libRmath = LibraryProduct("./src", "libRmath")

                # The libRmath makefile is kind of dumb
                makevars = `fPIC=-fPIC`
                if startswith(String(platform), "win")
                    makevars = `$makevars OS=Windows_NT`
                end

                steps = [
                    `make clean`,
                    `make $(makevars) -j$(min(Sys.CPU_CORES + 1,8))`,
                ]

                dep = Dependency(src_name, [libRmath], steps, platform, prefix)
                build(dep; verbose=true)

                # Create the folders we're going to install into
                mkdir(libdir(prefix))

                # Once it's been built, manually copy `libRmath` to `prefix`
                # We have to do this because the Makefile doesn't install it.
                libRmath_path = locate(libRmath, platform=platform, verbose=true)
                cp(libRmath_path, joinpath(libdir(prefix), basename(libRmath_path)))
                # Also copy in the headers, because why not.
                cp("./include", includedir(prefix))
            end

            # Once we're built up, go ahead and package this prefix out
            tarball_path, tarball_hash = package(prefix, joinpath(out_path, src_name); platform=platform, verbose=true)
            products[target] = (basename(tarball_path), tarball_hash)
        end
    end

    # Finally, destroy the build_path
    rm(build_path; recursive=true)
end

# In the end, dump an informative message telling the user how to download/install these
info("Hash/filename pairings:")
for target in keys(products)
    filename, hash = products[target]
    println("    :$(platform_key(target)) => (\"\$bin_prefix/$(filename)\", \"$(hash)\"),")
end

ext: jl

max_stars_repo_name: staticfloat/Rmath-julia-builder

max_issues_repo_name: staticfloat/Rmath-julia-builder

max_forks_repo_name: staticfloat/Rmath-julia-builder


text:
#!/usr/bin/env julia
import vfld

vfld.main(ARGS)


text:
using ITensorInfiniteMPS
using ITensors

##############################################################################
# VUMPS/TDVP parameters
#

maxdim = 20 # Maximum bond dimension
cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension
max_vumps_iters = 10 # Maximum number of iterations of the VUMPS/TDVP algorithm at a fixed bond dimension
tol = 1e-5 # Precision error tolerance for outer loop of VUMPS or TDVP
outer_iters = 5 # Number of times to increase the bond dimension
time_step = -Inf # -Inf corresponds to VUMPS, finite time_step corresponds to TDVP
solver_tol = (x -> x / 100) # Tolerance for the local solver (eigsolve in VUMPS and exponentiate in TDVP)
multisite_update_alg = "parallel" # Choose between ["sequential", "parallel"]. Only parallel works with TDVP.
conserve_qns = true # Whether or not to conserve spin parity
nsite = 2 # Number of sites in the unit cell
localham_type = ITensor # Can choose `ITensor` or `MPO`

# Parameters of the transverse field Ising model
model_params = (J=1.0, h=0.9)

##############################################################################
# CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
#

model = Model"ising"()

function space_shifted(::Model"ising", q̃sz; conserve_qns=true)
  if conserve_qns
    return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
  else
    return [QN() => 2]
  end
end

# Shift the QNs by 1 so that the flux of the unit cell is zero
# and therefore the flux density of the uniform state is zero
# to avoid diverging flux in the thermodynamic limit.
space_ = fill(space_shifted(model, 1; conserve_qns=conserve_qns), nsite)
s = infsiteinds("S=1/2", nsite; space=space_)
initstate(n) = "↑"
ψ = InfMPS(s, initstate)

# Form the Hamiltonian
H = InfiniteSum{localham_type}(model, s; model_params...)

# Check translational invariance
@show norm(contract(ψ.AL[1:nsite]..., ψ.C[nsite]) - contract(ψ.C[0], ψ.AR[1:nsite]...))

vumps_kwargs = (
  tol=tol,
  maxiter=max_vumps_iters,
  solver_tol=solver_tol,
  multisite_update_alg=multisite_update_alg,
)
subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)

# Alternate steps of running VUMPS and increasing the bond dimension
@time for _ in 1:outer_iters
  println("\nIncrease bond dimension")
  global ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...)
  println("Run VUMPS with new bond dimension")
  global ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...)
end

# Check translational invariance
@show norm(contract(ψ.AL[1:nsite]..., ψ.C[nsite]) - contract(ψ.C[0], ψ.AR[1:nsite]...))

#
# Compare to DMRG
#

nsite_finite = 100
s_finite = siteinds("S=1/2", nsite_finite; conserve_szparity=true)
H_finite = MPO(model, s_finite; model_params...)
ψ_finite = randomMPS(s_finite, initstate)
@show flux(ψ_finite)
sweeps = Sweeps(10)
setmaxdim!(sweeps, maxdim)
setcutoff!(sweeps, cutoff)
energy_finite_total, ψ_finite = @time dmrg(H_finite, ψ_finite, sweeps)
@show energy_finite_total / nsite_finite

function energy_local(ψ1, ψ2, h::ITensor)
  ϕ = ψ1 * ψ2
  return (noprime(ϕ * h) * dag(ϕ))[]
end

energy_local(ψ1, ψ2, h::MPO) = energy_local(ψ1, ψ2, prod(h))

function ITensors.expect(ψ, o)
  return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[]
end

# Exact energy at criticality: 4/pi = 1.2732395447351628

n_finite = nsite_finite ÷ 2
orthogonalize!(ψ_finite, n_finite)
hn_finite = ITensor(model, s_finite[n_finite], s_finite[n_finite + 1]; model_params...)
energy_finite = energy_local(ψ_finite[n_finite], ψ_finite[n_finite + 1], hn_finite)
energy_infinite = energy_local(ψ.AL[1], ψ.AL[2] * ψ.C[2], H[(1, 2)])
@show energy_finite, energy_infinite
@show abs(energy_finite - energy_infinite)

energy_exact = reference(model, Observable("energy"); model_params...)
@show energy_exact

Sz1_finite = expect(ψ_finite[n_finite], "Sz")
orthogonalize!(ψ_finite, n_finite + 1)
Sz2_finite = expect(ψ_finite[n_finite + 1], "Sz")
Sz1_infinite = expect(ψ.AL[1] * ψ.C[1], "Sz")
Sz2_infinite = expect(ψ.AL[2] * ψ.C[2], "Sz")

@show Sz1_finite, Sz2_finite
@show Sz1_infinite, Sz2_infinite


text:
using OrdinaryDiffEq, Test, LinearAlgebra
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear
prob = prob_ode_2Dlinear
choice_function(integrator) = (Int(integrator.t<0.5) + 1)
alg_double = CompositeAlgorithm((Tsit5(),Tsit5()),choice_function)
alg_double2 = CompositeAlgorithm((Vern6(),Vern6()),choice_function)
alg_switch = CompositeAlgorithm((Tsit5(),Vern7()),choice_function)

@time sol1 = solve(prob_ode_linear,alg_double)
@time sol2 = solve(prob_ode_linear,Tsit5())
@test sol1.t == sol2.t
@test sol1(0.8) == sol2(0.8)

integrator1 = init(prob,alg_double2)
integrator2 = init(prob,Vern6())
solve!(integrator1)
solve!(integrator2)

@test integrator1.sol.t == integrator2.sol.t

sol = solve(prob,alg_switch)
@inferred DiffEqBase.__init(prob, alg_switch)
v = @inferred OrdinaryDiffEq.ode_interpolant(1.0, integrator1, integrator1.opts.save_idxs, Val{0})
@inferred OrdinaryDiffEq.ode_interpolant!(v, 1.0, integrator1, integrator1.opts.save_idxs, Val{0})
v = @inferred OrdinaryDiffEq.ode_extrapolant(1.0, integrator1, integrator1.opts.save_idxs, Val{0})
@inferred OrdinaryDiffEq.ode_extrapolant!(v, 1.0, integrator1, integrator1.opts.save_idxs, Val{0})

condition(u,t,integrator) = t==192.0
function affect!(integrator)
   integrator.u[1] += 14000
   integrator.u[2] += 14000
end
A = [-0.027671669470584172 -0.0 -0.0 -0.0 -0.0 -0.0; -0.0 -0.05540281553537378 -0.0 -0.0 -0.0 -0.0; 0.011534597161021629 0.011933539591245327 -0.24891886153387743 0.023054812171672122 0.0 0.0; 0.0 0.0 0.17011732278405356 -0.023054812171672122 0.0 0.0; 0.01613707230956254 0.04346927594412846 0.03148193084515083 0.0 -1.5621055510998967e9 7.293040577236404; 0.0 0.0 0.0 0.0 1.559509231932001e9 -7.293040577236404]
prob = ODEProblem((du,u,p,t) -> mul!(du,A,u), zeros(6), (0.0, 1000), tstops=[192], callback=DiscreteCallback(condition, affect!));
sol = solve(prob, alg=AutoVern7(Rodas5()))
@test sol.t[end] == 1000.0


text:
import Base: convert, +, *, /, -, isless

###### Filesize algebra ######

export B, kB, MB, GB, TB

struct Bytes
    val::Float64
end

(*)(n::Number, b::Bytes) = Bytes(n*b.val)
(+)(a::Bytes, b::Bytes) = Bytes(a.val+b.val)
(-)(a::Bytes, b::Bytes) = Bytes(a.val-b.val)
isless(a::Bytes, b::Bytes) = a.val < b.val
(/)(a::Bytes, b::Number) = Bytes(ceil(UInt64, a.val/b))
(/)(a::Bytes, b::Bytes) = a.val/b.val

const B = Bytes(1)
const kB = 1024B
const MB = 1024kB
const GB = 1024MB
const TB = 1024GB

convert(::Type{Bytes}, n::Integer) = Bytes(n)

##### Base patches ######

_intersect(x,y) = intersect(x,y)

_intersect(::Colon, ::Colon) = Colon()
_intersect(::Colon, r) = r
_intersect(r, ::Colon) = r

##### Utility functions ######

"""
Utility function to divide the range `range` into `n` chunks
"""
function split_range(range::AbstractRange{T}, n) where T
    len = length(range)

    starts = len >= n ?
        round.(T, Base.range(first(range), stop=last(range)+1, length=n+1)) :
        vcat(collect(first(range):(last(range)+1)), zeros(T, n-len))

    map((x,y)->x:y, starts[1:end-1], starts[2:end] .- 1)
end

function split_range(r::AbstractRange{Char}, n)
    map((x) -> Char(first(x)):Char(last(x)), split_range(Int(first(r)):Int(last(r)), n))
end

"""
    split_range_interval(range, n)
split a range into pieces each of length `n` or lesser
"""
function split_range_interval(range, n)
    f = first(range)
    N = length(range)
    npieces = ceil(Int, N / n)

    ranges = Array{UnitRange}(npieces)
    for i=1:npieces
        ranges[i] = f:min(N, f+n-1)
        f += n
    end
    ranges
end

###### interface macro ######

getvarname(x::Expr) = x.args[1]
getvarname(x::Symbol) = x

"""
    @unimplemented fname(<args...>)

While it is nice to define generic function ad-hoc, it can sometimes
get confusing to figure out which method is missing.
`@interface` creates a function which errors out pointing
which method is missing.
"""
macro unimplemented(expr)
    @assert expr.head == :call

    fname = expr.args[1]
    args = expr.args[2:end]
    sig = string(expr)

    vars = map(x->getvarname(x), args)
    typs = Expr(:vect, map(x -> :(typeof($x)), vars)...,)

    :(function $(esc(fname))($(args...,))
        error(string("a method of ", $sig, " specialized to ", ($typs...,), " is missing"))
    end)
end

"""
Tree reduce
"""
function treereduce(f, xs)
    length(xs) == 1 && return xs[1]
    l = length(xs)
    m = div(l, 2)
    f(treereduce(f, xs[1:m]), treereduce(f, xs[m+1:end]))
end

const  ENABLE_DEBUG = true # Will remove code under @dbg

"""
Run a block of code only if DEBUG is true
"""
macro dbg(expr)
    if ENABLE_DEBUG
        esc(expr)
    end
end

if !isdefined(Base, :reduced_indices)
    function reduced_dims(x, dim)
        Base.reduced_dims(size(x), dim)
    end
else
    reduced_dims(x, dim) = Base.reduced_indices(axes(x), dim)
end

function treereducedim(op, xs::Array, dim::Int)
    l = size(xs, dim)
    colons = Any[Colon() for i=1:length(size(xs))]
    if dim > length(size(xs))
        return xs
    end
    ys = treereduce((x,y)->map(op, x,y), Any[(
        colons[dim] = [i];
        view(xs, colons...)
    ) for i=1:l])
    reshape(ys[:], reduced_dims(xs, dim))
end

function treereducedim(op, xs::Array, dim::Tuple)
    reduce((prev, d) -> treereducedim(op, prev, d), dim, init=xs)
end

function allslices(xs, n)
    idx = Any[Colon() for i in 1:ndims(xs)]
    [(idx[n] = j; view(xs, idx...)) for j in 1:size(xs, n)]
end

function treereduce_nd(fs, xs)
    n = ndims(xs)
    if n==1
        treereduce(fs[n], xs)
    else
        treereduce(fs[n], map(ys->treereduce_nd(fs, ys), allslices(xs, n)))
    end
end

function setindex(x::NTuple{N}, idx, v) where N
    map(ifelse, ntuple(x->idx === x, Val(N)), ntuple(x->v, Val(N)), x)
end

function showloc(f, argcount)
    args = ntuple(x->Any, argcount)
    ms = methods(f)
    if length(ms) == 1
        string(first(ms))
    else
        string(f)
    end
end


text:
using CSTParser, Test, CSTParser.Tokenize
using CSTParser: parse, remlineinfo!, span, headof, kindof, valof

@testset "CSTParser" begin
    include("spec.jl")
    include("parser.jl")
    include("interface.jl")
    include("display.jl")
    include("iterate.jl")
    CSTParser.check_base()
    include("errparse.jl")

end


text:
################################################################################
#  Copyright 2020, Marta Vanin, Tom Van Acker                                  #
################################################################################
# PowerModelsDistributionStateEstimation.jl                                    #
# An extention package of PowerModels(Distribution).jl for Static Power System #
# State Estimation.                                                            #
################################################################################

abstract type ConversionType end

struct SquareFraction<:ConversionType
    msr_id::Int64
    cmp_type::Symbol
    cmp_id::Int64
    bus_ind::Int64
    numerator::Array
    denominator::Array
end

struct Square<:ConversionType
    msr_id::Int64
    cmp_type::Symbol
    cmp_id::Int64
    bus_ind::Int64
    elements::Array
end

struct Multiplication<:ConversionType
    msr_sym::Symbol
    msr_id::Int64
    cmp_type::Symbol
    cmp_id::Int64
    bus_ind::Int64
    mult1::Array
    mult2::Array
end

struct Tangent<:ConversionType
    msr_id::Int64
    cmp_type::Symbol
    cmp_id::Int64
    numerator::Symbol
    denominator::Symbol
end

struct Fraction<:ConversionType
    msr_type::Symbol
    msr_id::Int64
    cmp_type::Symbol
    cmp_id::Int64
    bus_ind::Int64
    numerator::Array
    denominator::Symbol
end

struct MultiplicationFraction<:ConversionType
    msr_type::Symbol
    msr_id::Int64
    cmp_type::Symbol
    cmp_id::Int64
    bus_ind::Int64
    power::Array
    voltage::Array
end

function assign_conversion_type_to_msr(pm::_PMD.AbstractUnbalancedACPModel,i,msr::Symbol;nw=nw)
    cmp_id = _PMD.ref(pm, nw, :meas, i, "cmp_id")
    if msr == :cm
        msr_type = SquareFraction(i,:branch, cmp_id, _PMD.ref(pm,nw,:branch,cmp_id)["f_bus"], [:p, :q], [:vm])
    elseif msr == :cmg
        msr_type = SquareFraction(i,:gen, cmp_id, _PMD.ref(pm,nw,:gen,cmp_id)["gen_bus"], [:pg, :qg], [:vm])
    elseif msr == :cmd
        msr_type = SquareFraction(i,:load, cmp_id, _PMD.ref(pm,nw,:load,cmp_id)["load_bus"], [:pd, :qd], [:vm])
    elseif msr == :cr
        msr_type = Fraction(msr, i,:branch, cmp_id, _PMD.ref(pm,nw,:branch,cmp_id)["f_bus"], [:p, :q, :va], :vm)
    elseif msr == :crg
        msr_type = Fraction(msr, i,:gen, cmp_id, _PMD.ref(pm,nw,:gen,cmp_id)["gen_bus"], [:pg, :qg, :va], :vm)
    elseif msr == :crd
        msr_type = Fraction(msr, i,:load, cmp_id, _PMD.ref(pm,nw,:load,cmp_id)["load_bus"], [:pd, :qd, :va], :vm)
    elseif msr == :ci
        msr_type = Fraction(msr, i,:branch, cmp_id, _PMD.ref(pm,nw,:branch,cmp_id)["f_bus"], [:p, :q, :va], :vm)
    elseif msr == :cig
        msr_type = Fraction(msr, i,:gen, cmp_id, _PMD.ref(pm,nw,:gen,cmp_id)["gen_bus"], [:pg, :qg, :va], :vm)
    elseif msr == :cid
        msr_type = Fraction(msr, i,:load, cmp_id, _PMD.ref(pm,nw,:load,cmp_id)["load_bus"], [:pd, :qd, :va], :vm)
    else
       error("the chosen measurement $(msr) at $(_PMD.ref(pm, nw, :meas, i, "cmp")) $(_PMD.ref(pm, nw, :meas, i, "cmp_id")) is not supported and should be removed")
    end
    return msr_type
end

function assign_conversion_type_to_msr(pm::_PMD.AbstractUnbalancedACRModel,i,msr::Symbol;nw=nw)
    cmp_id = _PMD.ref(pm, nw, :meas, i, "cmp_id")
    if msr == :vm
        msr_type = Square(i,:bus, cmp_id, _PMD.ref(pm,nw,:bus,cmp_id)["index"], [:vi, :vr])
    elseif msr == :va
        msr_type = Tangent(i, :bus, cmp_id, :vi, :vr)
    elseif msr == :cm
        msr_type = SquareFraction(i,:branch, cmp_id, _PMD.ref(pm,nw,:branch,cmp_id)["f_bus"], [:p, :q], [:vi, :vr])
    elseif msr == :cmg
        msr_type = SquareFraction(i,:gen, cmp_id, _PMD.ref(pm,nw,:gen,cmp_id)["gen_bus"], [:pg, :qg], [:vi, :vr])
    elseif msr == :cmd
        msr_type = SquareFraction(i,:load, cmp_id, _PMD.ref(pm,nw,:load,cmp_id)["load_bus"], [:pd, :qd], [:vi, :vr])
    elseif msr == :cr
        msr_type = MultiplicationFraction(msr, i,:branch, cmp_id, _PMD.ref(pm,nw,:branch,cmp_id)["f_bus"], [:p, :q], [:vr, :vi])
    elseif msr == :crg
        msr_type = MultiplicationFraction(msr, i,:gen, cmp_id, _PMD.ref(pm,nw,:gen,cmp_id)["gen_bus"], [:pg, :qg], [:vr, :vi])
    elseif msr == :crd
        msr_type = MultiplicationFraction(msr, i,:load, cmp_id, _PMD.ref(pm,nw,:load,cmp_id)["load_bus"], [:pd, :qd], [:vr, :vi])
    elseif msr == :ci
        msr_type = MultiplicationFraction(msr, i,:branch, cmp_id, _PMD.ref(pm,nw,:branch,cmp_id)["f_bus"], [:p, :q], [:vr, :vi])
    elseif msr == :cig
        msr_type = MultiplicationFraction(msr, i,:gen, cmp_id, _PMD.ref(pm,nw,:gen,cmp_id)["gen_bus"], [:pg, :qg], [:vr, :vi])
    elseif msr == :cid
        msr_type = MultiplicationFraction(msr, i,:load, cmp_id, _PMD.ref(pm,nw,:load,cmp_id)["load_bus"], [:pd, :qd], [:vr, :vi])
    else
       error("the chosen measurement $(msr) at $(_PMD.ref(pm, nw, :meas, i, "cmp")) $(_PMD.ref(pm, nw, :meas, i, "cmp_id")) is not supported and will be ignored")
    end
    return msr_type
end

function assign_conversion_type_to_msr(pm::_PMD.AbstractUnbalancedIVRModel,i,msr::Symbol;nw=nw)
    cmp_id = _PMD.ref(pm, nw, :meas, i, "cmp_id")
    if msr == :vm
        msr_type = Square(i,:bus, cmp_id, _PMD.ref(pm,nw,:bus,cmp_id)["index"], [:vi, :vr])
    elseif msr == :va
        msr_type = Tangent(i, :bus, cmp_id, :vi, :vr)
    elseif msr == :cm
        msr_type = Square(i,:branch, cmp_id, _PMD.ref(pm,nw,:branch,cmp_id)["f_bus"], [:cr, :ci])
    elseif msr == :cmg
        msr_type = Square(i,:gen, cmp_id, _PMD.ref(pm,nw,:gen,cmp_id)["gen_bus"], [:crg, :cig])
    elseif msr == :cmd
        msr_type = Square(i,:load, cmp_id, _PMD.ref(pm,nw,:load,cmp_id)["load_bus"], [:crd, :cid])
    elseif msr == :ca
        msr_type = Tangent(i, :branch, cmp_id, :ci, :cr)
    elseif msr == :cag
        msr_type = Tangent(i, :gen, cmp_id, :cig, :crg)
    elseif msr == :cad
        msr_type = Tangent(i, :load, cmp_id, :cid, :crd)
    elseif msr == :p
        msr_type = Multiplication(msr, i,:branch, cmp_id, _PMD.ref(pm,nw,:branch,cmp_id)["f_bus"], [:cr, :ci], [:vr, :vi])
    elseif msr == :pg
        msr_type = Multiplication(msr, i,:gen, cmp_id, _PMD.ref(pm,nw,:gen,cmp_id)["gen_bus"], [:crg, :cig], [:vr, :vi])
    elseif msr == :pd
        msr_type = Multiplication(msr, i,:load, cmp_id, _PMD.ref(pm,nw,:load,cmp_id)["load_bus"], [:crd, :cid], [:vr, :vi])
    elseif msr == :q
        msr_type = Multiplication(msr, i,:branch, cmp_id, _PMD.ref(pm,nw,:branch,cmp_id)["f_bus"], [:cr, :ci], [:vr, :vi])
    elseif msr == :qg
        msr_type = Multiplication(msr, i,:gen, cmp_id, _PMD.ref(pm,nw,:gen,cmp_id)["gen_bus"], [:crg, :cig], [:vr, :vi])
    elseif msr == :qd
        msr_type = Multiplication(msr, i,:load, cmp_id, _PMD.ref(pm,nw,:load,cmp_id)["load_bus"], [:crd, :cid], [:vr, :vi])
    else
       error("the chosen measurement $(msr) at $(_PMD.ref(pm, nw, :meas, i, "cmp")) $(_PMD.ref(pm, nw, :meas, i, "cmp_id")) is not supported and should be removed")
    end
    return msr_type
end

function assign_conversion_type_to_msr(pm::_PMD.LinDist3FlowPowerModel,i,msr::Symbol;nw=nw)
    cmp_id = _PMD.ref(pm, nw, :meas, i, "cmp_id")
    if msr == :vm
        msr_type = Square(i,:bus, cmp_id, _PMD.ref(pm,nw,:bus,cmp_id)["index"], [:w])
    elseif msr == :cm
        msr_type = SquareFraction(i,:branch, cmp_id, _PMD.ref(pm,nw,:branch,cmp_id)["f_bus"], [:p, :q], [:w])
    elseif msr == :cmg
        msr_type = SquareFraction(i,:gen, cmp_id, _PMD.ref(pm,nw,:gen,cmp_id)["gen_bus"], [:pg, :qg], [:w])
    elseif msr == :cmd
        msr_type = SquareFraction(i,:load, cmp_id, _PMD.ref(pm,nw,:load,cmp_id)["load_bus"], [:pd, :qd], [:w])
    else
       error("the chosen measurement $(msr) at $(_PMD.ref(pm, nw, :meas, i, "cmp")) $(_PMD.ref(pm, nw, :meas, i, "cmp_id")) is not supported and should be removed")
    end
    return msr_type
end

function assign_conversion_type_to_msr(pm::_PMD.SDPUBFPowerModel,i,msr::Symbol;nw=nw)
   error("Currently only a limited amount of measurement types is supported for the SDP model, $(msr) is not available")
end

function no_conversion_needed(pm::_PMD.AbstractUnbalancedACPModel, msr_var::Symbol)
  return msr_var ∈ [:vm, :va, :pd, :qd, :pg, :qg, :p, :q]
end

function no_conversion_needed(pm::_PMD.AbstractUnbalancedACRModel, msr_var::Symbol)
  return msr_var ∈ [:vr, :vi, :pd, :qd, :pg, :qg, :p, :q]
end

function no_conversion_needed(pm::_PMD.AbstractUnbalancedIVRModel, msr_var::Symbol)
  return msr_var ∈ [:vr, :vi, :cr, :ci, :crg, :cig, :crd, :cid]
end

function no_conversion_needed(pm::_PMD.AbstractLPUBFModel, msr_var::Symbol)
    return msr_var ∈ [:w, :pd, :qd, :pg, :qg, :p, :q]
end

function no_conversion_needed(pm::_PMD.SDPUBFPowerModel, msr_var::Symbol)
    return msr_var ∈ [:w, :pd, :qd, :pg, :qg]
end

function create_conversion_constraint(pm::_PMD.AbstractUnbalancedPowerModel, original_var, msr::SquareFraction; nw=nw)

    new_var_num = []
    for nvn in msr.numerator
        if occursin("v", String(nvn)) && msr.cmp_type != :bus
            push!(new_var_num, _PMD.var(pm, nw, nvn, msr.bus_ind))
        elseif msr.cmp_type == :branch
            push!(new_var_num, _PMD.var(pm, nw, nvn, (msr.cmp_id, msr.bus_ind, _PMD.ref(pm,nw,:branch,msr.cmp_id)["t_bus"])))
        else
            push!(new_var_num, _PMD.var(pm, nw, nvn, msr.cmp_id))
        end
    end

    msr.cmp_type == :branch ? id = (msr.cmp_id,  msr.bus_ind, _PMD.ref(pm,nw,:branch,msr.cmp_id)["t_bus"]) : id = msr.cmp_id
    conn = get_active_connections(pm, nw, msr.cmp_type, msr.cmp_id)

    if typeof(pm) <: _PMD.LinDist3FlowPowerModel

        new_var_den = [_PMD.var(pm, nw, msr.denominator, msr.cmp_id)]

        JuMP.@NLconstraint(pm.model, [c in conn],
            original_var[id][c]^2 == (sum( n[c]^2 for n in new_var_num ))/
                       (sum( d[c] for d in new_var_den))
            )
    else
        new_var_den = []
        for nvd in msr.denominator
            if !isa(nvd, Symbol) #only case is when I have an array of ones
                push!(new_var_den, nvd)
            elseif occursin("v", String(nvd)) && msr.cmp_type != :bus
                push!(new_var_den, _PMD.var(pm, nw, nvd, msr.bus_ind))
            else
                push!(new_var_den, _PMD.var(pm, nw, nvd, msr.cmp_id))
            end
        end

        JuMP.@NLconstraint(pm.model, [c in conn],
            original_var[id][c]^2 == (sum( n[c]^2 for n in new_var_num ))/
                       (sum( d[c]^2 for d in new_var_den))
            )
    end
end

function create_conversion_constraint(pm::_PMD.AbstractUnbalancedPowerModel, original_var, msr::Square; nw=nw)

    new_var = []
    for nvn in msr.elements
        if msr.cmp_type == :branch
            push!(new_var, _PMD.var(pm, nw, nvn, (msr.cmp_id, _PMD.ref(pm,nw,:branch,msr.cmp_id)["f_bus"], _PMD.ref(pm,nw,:branch,msr.cmp_id)["t_bus"])))
        else
            push!(new_var, _PMD.var(pm, nw, nvn, msr.cmp_id))
        end
    end

    msr.cmp_type == :branch ? id = (msr.cmp_id,  _PMD.ref(pm,nw,:branch,msr.cmp_id)["f_bus"], _PMD.ref(pm,nw,:branch,msr.cmp_id)["t_bus"]) : id = msr.cmp_id
    conn = get_active_connections(pm, nw, msr.cmp_type, msr.cmp_id)

    if typeof(pm) <: _PMD.LinDist3FlowPowerModel
        JuMP.@constraint(pm.model, [c in conn],
            original_var[id][c]^2 == (sum( n[c] for n in new_var ))
            )
    else
        JuMP.@constraint(pm.model, [c in conn],
            original_var[id][c]^2 == (sum( n[c]^2 for n in new_var ))
            )
    end
end

function create_conversion_constraint(pm::_PMD.AbstractUnbalancedPowerModel, original_var, msr::Multiplication; nw=nw)

    m1 = []
    m2 = []

    for m in msr.mult1
        if occursin("v", String(m)) && msr.cmp_type != :bus
            push!(m1, _PMD.var(pm, nw, m, msr.bus_ind))
        elseif msr.cmp_type == :branch
            push!(m1, _PMD.var(pm, nw, m, (msr.cmp_id, msr.bus_ind, _PMD.ref(pm, nw, :branch,msr.cmp_id)["t_bus"])))
        else
            push!(m1, _PMD.var(pm, nw, m, msr.cmp_id))
        end
    end

    for mm in msr.mult2
        if occursin("v", String(mm)) && msr.cmp_type != :bus
            push!(m2, _PMD.var(pm, nw, mm, msr.bus_ind))
        elseif msr.cmp_type == :branch
            push!(m2, _PMD.var(pm, nw, mm, (msr.cmp_id, msr.bus_ind, _PMD.ref(pm, nw, :branch,msr.cmp_id)["t_bus"])))
        else
            push!(m2, _PMD.var(pm, nw, mm, msr.cmp_id))
        end
    end
    msr.cmp_type == :branch ? id = (msr.cmp_id,  msr.bus_ind, _PMD.ref(pm,nw,:branch,msr.cmp_id)["t_bus"]) : id = msr.cmp_id
    conn = get_active_connections(pm, nw, msr.cmp_type, msr.cmp_id)

    if occursin("p", String(msr.msr_sym))
        JuMP.@constraint(pm.model, [c in conn],
            original_var[id][c] == m1[1][c]*m2[1][c]+m1[2][c]*m2[2][c]
            )
    elseif occursin("q", String(msr.msr_sym))
        JuMP.@constraint(pm.model, [c in conn],
            original_var[id][c] == -m1[2][c]*m2[1][c]+m1[1][c]*m2[2][c]
            )
    end
end

function create_conversion_constraint(pm::_PMD.AbstractUnbalancedPowerModel, original_var, msr::Tangent; nw=nw)
    #TODO for v0.2.0 this needs to be general to every distribution or we need to provide an exception
    @warn "Performing a Tangent conversion only makes sense for Normal distributions and is in general not advised"
    conn = get_active_connections(pm, nw, msr.cmp_type, msr.cmp_id)
    for c in conn
        if _PMD.ref(pm, nw, :meas, msr.msr_id, "dst")[c] != 0.0

            μ_tan = tan(_DST.mean(_PMD.ref(pm, nw, :meas, msr.msr_id, "dst")[c]))
            σ_tan = abs(sec(μ_tan)*_DST.std(_PMD.ref(pm, nw, :meas, msr.msr_id, "dst")[c]))
            _PMD.ref(pm, nw, :meas, msr.msr_id, "dst")[c] = _DST.Normal( μ_tan, σ_tan )

            if msr.cmp_type == :branch
                num = _PMD.var(pm, nw, msr.numerator, (msr.cmp_id, _PMD.ref(pm, nw, :branch,msr.cmp_id)["f_bus"], _PMD.ref(pm, nw, :branch,msr.cmp_id)["t_bus"]))
                den = _PMD.var(pm, nw, msr.denominator, (msr.cmp_id, _PMD.ref(pm, nw, :branch,msr.cmp_id)["f_bus"], _PMD.ref(pm, nw, :branch,msr.cmp_id)["t_bus"]))
            else
                num = _PMD.var(pm, nw, msr.numerator, msr.cmp_id)
                den = _PMD.var(pm, nw, msr.numerator, msr.cmp_id)
            end
            msr.cmp_type == :branch ? id = (msr.cmp_id, _PMD.ref(pm,nw,:branch,msr.cmp_id)["f_bus"], _PMD.ref(pm,nw,:branch,msr.cmp_id)["t_bus"]) : id = msr.cmp_id
            JuMP.@NLconstraint(pm.model,
                original_var[id][c]*den[c] == num[c]
                )
        end
    end
end

function create_conversion_constraint(pm::_PMD.AbstractUnbalancedPowerModel, original_var, msr::Fraction; nw=nw)
    num = []
    for n in msr.numerator
        if occursin("v", String(n))
            push!(num, _PMD.var(pm, nw, n, msr.bus_ind))
        elseif !occursin("v", String(n)) && msr.cmp_type != :branch
            push!(num, _PMD.var(pm, nw, n, msr.cmp_id))
        else
            push!(num, _PMD.var(pm, nw, n, (msr.cmp_id, _PMD.ref(pm, nw, :branch,msr.cmp_id)["f_bus"], _PMD.ref(pm, nw, :branch,msr.cmp_id)["t_bus"])))
        end
    end

    den = _PMD.var(pm, nw, msr.denominator, msr.bus_ind)
    msr.cmp_type == :branch ? id = (msr.cmp_id,  msr.bus_ind, _PMD.ref(pm,nw,:branch,msr.cmp_id)["t_bus"]) : id = msr.cmp_id
    conn = get_active_connections(pm, nw, msr.cmp_type, msr.cmp_id)

    if occursin("r", String(msr.msr_type))
        JuMP.@NLconstraint(pm.model, [c in conn],
            original_var[id][c]*den[c] == num[1][c]*cos(num[3][c])+num[2][c]*sin(num[3][c])
            )
    elseif occursin("i", String(msr.msr_type))
        JuMP.@NLconstraint(pm.model, [c in conn],
            original_var[id][c]*den[c] == -num[2][c]*cos(num[3][c])+num[1][c]*sin(num[3][c])
            )
    else
        error("wrong measurement association")
    end
end

function create_conversion_constraint(pm::_PMD.AbstractUnbalancedPowerModel, original_var, msr::MultiplicationFraction; nw=nw)

    p = []
    v = []
    for pw in msr.power
        if msr.cmp_type == :branch
            push!(p, _PMD.var(pm, nw, pw, (msr.cmp_id, _PMD.ref(pm, nw, :branch,msr.cmp_id)["f_bus"], _PMD.ref(pm, nw, :branch,msr.cmp_id)["t_bus"])))
        else
            push!(p, _PMD.var(pm, nw, pw, msr.cmp_id))
        end
    end
    for vl in msr.voltage
        push!(v, _PMD.var(pm, nw, vl, msr.bus_ind))
    end

    msr.cmp_type == :branch ? id = (msr.cmp_id,  _PMD.ref(pm, nw, :branch,msr.cmp_id)["f_bus"], _PMD.ref(pm,nw,:branch,msr.cmp_id)["t_bus"]) : id = msr.cmp_id
    conn = get_active_connections(pm, nw, msr.cmp_type, msr.cmp_id)

    if occursin("cr", string(msr.msr_type))
        JuMP.@NLconstraint(pm.model, [c in conn],
            original_var[id][c] == (p[1][c]*v[1][c]+p[2][c]*v[2][c])/(v[1][c]^2+v[2][c]^2)
            )
    elseif occursin("ci", string(msr.msr_type))
        JuMP.@NLconstraint(pm.model, [c in conn],
            original_var[id][c] == (-p[2][c]*v[1][c]+p[1][c]*v[2][c])/(v[1][c]^2+v[2][c]^2)
            )
    end
end


text:
function solver(
    Nt,
    dt,
    L,
    N,
    phi,
    x,
    x_dot,
    parameter_phi,
    parameter_impurity,
    parameter_V,
    phi_fct,
    impurity_fct,
    V_fct,
    times_save,
)
    V = [zero(phi) for _ = 1:length(size(phi))+1]
    xarray = [i for i in range(-L / 2.0, L / 2.0, length = N)]
    update_V = ((V, x, t) -> V_fct(V, x, t, xarray, parameter_V...))
    phi_step = ((phi, V, t) -> phi_fct(phi, V, t, dt, parameter_phi...))
    x_step = (
        (phi, V, x, xdot, t) ->
            impurity_fct(phi, V, x, xdot, t, dt, parameter_impurity...)
    )
    phi_r = [zero(phi) for _ = 1:length(times_save)]
    x_r = [zero(x) for _ = 1:length(times_save)]
    x_dot_r = [zero(x_dot) for _ = 1:length(times_save)]
    j = 1
    for i = 1:Nt
        update_V(V, x, i * dt)
        x_step(phi, V, x, x_dot, i * dt)
        phi_step(phi, V, i * dt)
        #push!(o,observables(phi, x, xdot, i*dt)) # this is slow think about it!
        #push!(o,observables(phi, x , xdot,t))
        if i in times_save
            phi_r[j] = deepcopy(phi)
            x_r[j] = deepcopy(x)
            x_dot_r[j] = deepcopy(x_dot)
            j += 1
        end
    end
    return phi_r, x_r, x_dot_r
end

if abspath(PROGRAM_FILE) == @__FILE__
    include("helper.jl")
    include("SGPE.jl")
    include("NewtonImp.jl")
    using BenchmarkTools
    phi = ones(100, 100) + 0im * ones(100, 100)
    V = [
        ones(100, 100) + 0im * ones(100, 100),
        ones(100, 100) + 0im * ones(100, 100),
        ones(100, 100) + 0im * ones(100, 100),
    ]
    xarray = [i for i in range(-5 / 2.0, 5 / 2.0, length = 100)]
    stencil_bc = tuple_stencil_bc(phi)
    dt = 0.1
    lambda = 10.0
    m_b = 1.0
    mu = 0.1
    g = 1.0
    xstep = xarray[2] - xarray[1]
    T = 0.0
    M = 1.0
    re = 1.0
    parameter_phi = [lambda, m_b, mu, g, stencil_bc, xstep, T]
    parameter_impurity = [M, lambda, xstep]
    parameter_V = [re]
    N = 100
    L = 5.0
    x = [0.0, 0.0]
    x_dot = [0.0, 0.0]
end

max_stars_repo_name: phyjonas/ImpuritySGPE

max_issues_repo_name: phyjonas/ImpuritySGPE

max_forks_repo_name: phyjonas/ImpuritySGPE


text:
@testset "Metaevent" begin
    validtestvalues = [
        # Structure: dT (variable length), begin byte (0xFF), type byte, length (variable length type), data
        # dT is a variable length value which is handled by readvariablelength and not readsysexevent. We don't need to test it here, which
        # is why it's set to 0 in each case.
        ([0x00, 0xFF, 0x0F, 0x04, 0x11, 0x21, 0x53, 0x1F], MIDI.MetaEvent(0, 0x0F, [0x11, 0x21, 0x53, 0x1F])),
        ([0x00, 0xFF, 0x0F, 0x06, 0x25, 0x61, 0x23, 0x5B, 0x00, 0x02], MIDI.MetaEvent(0, 0x0F, [0x25, 0x61, 0x23, 0x5B, 0x00, 0x02]))
    ]

    invalidtestvalues = [
        ([0x00, 0xF0, 0x05, 0x11, 0xF7, 0x53, 0x1F], ErrorException)
    ]

    @testset "it should correctly generate a metaevent from raw data" begin
        for (input, output) in validtestvalues
            result = MIDI.readmetaevent(Int(input[1]), IOBuffer(input[2:length(input)]))
            @test result.dT == output.dT && result.data == output.data
        end
    end

    @testset "it should corrrectly write a MetaEvent to raw data" begin
        for (output, input) in validtestvalues
            buf = IOBuffer()
            MIDI.writeevent(buf, input)
            @test take!(buf) == output
        end
    end

    @testset "it should fail when invalid raw data is provided" begin
        for (input, errtype) in invalidtestvalues
            @test_throws errtype MIDI.readsysexevent(Int(input[1]), IOBuffer(input[2:length(input)]))
        end
    end
end


text:
import DataFrames

function disallowmissing(x::AbstractArray{Union{Missing, T}, N})::Array{T, N} where T where N
    result::Array{T, N} = convert(Array{T, N}, x)
    return result
end

"""
"""
function convert_value_to_missing! end

function convert_value_to_missing!(
        df::DataFrames.AbstractDataFrame,
        value,
        column_names::AbstractArray{Symbol},
        )::Nothing
    function f(old_vector::AbstractVector)::Vector{Any}
        new_vector = Vector{Any}(length(old_vector))
        for i = 1:length(old_vector)
            old_vector_ith_element = old_vector[i]
            if DataFrames.ismissing(old_vector_ith_element)
                new_vector[i] = DataFrames.missing
            elseif old_vector_ith_element == value
                new_vector[i] = DataFrames.missing
            else
                new_vector[i] = old_vector_ith_element
            end
        end
        return new_vector
    end
    transform_columns!(df, f, column_names)
    return nothing
end

function convert_value_to_missing!(
        df::DataFrames.AbstractDataFrame,
        value,
        column_name::Symbol,
        )::Nothing
    convert_value_to_missing!(df, value, Symbol[column_name])
    return nothing
end


text:
using LinearAlgebra, IntervalSets, UnPack
import ClimaCore:
    Domains, Topologies, Meshes, Spaces, Geometry, Operators, Fields

using Test

using StaticArrays, LinearAlgebra

function rotational_field(space, α0 = 45.0)
    coords = Fields.coordinate_field(space)
    map(coords) do coord
        ϕ = coord.lat
        λ = coord.long
        uu = (cosd(α0) * cosd(ϕ) + sind(α0) * cosd(λ) * sind(ϕ))
        uv = -sind(α0) * sind(λ)
        Geometry.UVVector(uu, uv)
    end
end

@testset "Spherical geometry properties" begin
    # test different combinations of odd/even to ensure pole is correctly
    # handled
    for Ne in (4, 5), Nq in (4, 5)
        FT = Float64
        radius = FT(3)

        domain = Domains.SphereDomain(radius)
        mesh = Meshes.EquiangularCubedSphere(domain, Ne)
        grid_topology = Topologies.Topology2D(mesh)
        quad = Spaces.Quadratures.GLL{Nq}()
        space = Spaces.SpectralElementSpace2D(grid_topology, quad)

        @test sum(ones(space)) ≈ 4pi * radius^2 rtol = 1e-3

        for α in [0.0, 45.0, 90.0]
            div = Operators.Divergence()
            u = rotational_field(space, α)
            divu = Spaces.weighted_dss!(div.(u))
            @test norm(divu) < 1e-2

            # test dss on UVcoordinates
            uu = Spaces.weighted_dss!(copy(u))
            @test norm(uu .- u) < 1e-14

            uᵢ = Geometry.transform.(Ref(Geometry.Covariant12Axis()), u)
            uuᵢ = Spaces.weighted_dss!(copy(uᵢ))
            @test norm(uuᵢ .- uᵢ) < 1e-14

            uⁱ = Geometry.transform.(Ref(Geometry.Contravariant12Axis()), u)
            uuⁱ = Spaces.weighted_dss!(copy(uⁱ))
            @test norm(uuⁱ .- uⁱ) < 1e-14
        end
    end
end


text:
function algebry(rel::Parallelism, coords)
    A = coords[rel.l.l]
    B = coords[rel.l.r]
    C = coords[rel.r.l]
    D = coords[rel.r.r]
    return ((B.x - A.x)*(D.y - C.y) - (B.y - A.y)*(D.x - C.x),)
end

function algebry(rel::Perpendicularity, coords)
    A = coords[rel.l.l]
    B = coords[rel.l.r]
    C = coords[rel.r.l]
    D = coords[rel.r.r]
    return ((B.x - A.x)*(D.x - C.x) + (B.y - A.y)*(D.y - C.y),)
end

function algebry(rel::Congruence{Segment}, coords)
    A = coords[rel.l.l]
    B = coords[rel.l.r]
    C = coords[rel.r.l]
    D = coords[rel.r.r]
    return ((B.x-A.x)^2 + (B.y-A.y)^2 - (D.x-C.x)^2 - (D.y-C.y)^2,)
end

function algebry(rel::Parallelogram, coords)
    A = coords[rel.A]
    B = coords[rel.B]
    C = coords[rel.C]
    D = coords[rel.D]
    return (B.x + D.x - A.x - C.x, B.y + D.y - A.y - C.y)
end

function algebry(rel::Midpoint, coords)
    M = coords[rel.M]
    A = coords[rel.seg.l]
    B = coords[rel.seg.r]
    return (2*M.x-A.x-B.x, 2*M.y-A.y-B.y)
end

function algebry(rel::Appartenence{Segment}, coords)
    P = coords[rel.l]
    A = coords[rel.r.l]
    B = coords[rel.r.r]
    return ((B.x - A.x)*(P.y - A.y) - (B.y - A.y)*(P.x - A.x),)
end

function algebry(rel::Intersection{Segment, Segment}, coords)
    return (algebry(Appartenence(rel.P, rel.l), coords)[1], algebry(Appartenence(rel.P, rel.r), coords)[1])
end

function algebry(rel::Projection, coords)
    return (algebry(Appartenence(rel.H, rel.seg), coords)[1], algebry(Perpendicularity(Segment(rel.P, rel.H), rel.seg), coords)[1])
end

function algebry(rel::Appartenence{Circle}, coords)
    OP = Segment(rel.r.O, rel.l)
    OA = Segment(rel.r.O, rel.r.A)
    return algebry(Congruence(OA, OP), coords)
end

function algebry(rel::CircleThreePoints, coords)
    OA = Segment(rel.O, rel.A)
    OB = Segment(rel.O, rel.B)
    OC = Segment(rel.O, rel.C)
    return (algebry(Congruence(OA, OB), coords)[1], algebry(Congruence(OA, OC), coords)[1])
end


text:
# %% Let everyone in the world play the St. Petersburg game
using Distributions
using DelimitedFiles
include("../util.jl")
# %%
# TODO Cache

g = Geometric(0.5)

population = Int(7.9e9)

function computeWinnings(N)
    winnings = zeros(Int, round(Int, log2(N)+5))
    for _ in 1:N
        w = rand(g)+1
        winnings[w] += 1
    end
    return winnings
end

@time winnings = computeWinnings(population)

# %%
writedlm("../../data/world.csv", winnings)
write(ARGS[1], table)

# %% Write to table
function commas(num::Integer)
    str = string(num)
    return replace(str, r"(?<=[0-9])(?=(?:[0-9]{3})+(?![0-9]))" => ",")
end

data = readdlm("../../data/world.csv") |> vec .|> x -> round(x, sigdigits=3) .|> Int64 .|> commas
payout = 1:length(data) .|> x -> 2^x .|> x -> round(x, sigdigits=3) .|> Int64 .|> commas
# %%
table = mat2latex([payout[1:34] data[1:34]])
write("../../tex/world.tex", table)
# %%

max_stars_repo_name: akoen/petersburg

max_issues_repo_name: akoen/petersburg

max_forks_repo_name: akoen/petersburg


text:
# IO functions

export AIO
"""
    struct AIO

Holds a C pointer `adios2_io *`.
"""
struct AIO
    ptr::Ptr{Cvoid}
    adios::Adios
    AIO(ptr::Ptr{Cvoid}, adios::Adios) = new(ptr, adios)
end

function Base.show(io::IO, ::MIME"text/plain", aio::AIO)
    return print(io, "AIO(0x$(string(UInt(aio.ptr); base=16)))")
end

export define_variable
"""
    variable = define_variable(io::AIO, name::AbstractString, type::Type,
                         shape::Union{Nothing,NTuple{N,Int} where N,CartesianIndex}=nothing,
                         start::Union{Nothing,NTuple{N,Int} where N,CartesianIndex}=nothing,
                         count::Union{Nothing,NTuple{N,Int} where N,CartesianIndex}=nothing,
                         constant_dims::Bool=false)
    variable::Union{Nothing,Variable}

Define a variable within `io`.

# Arguments
- `io`: handler that owns the variable
- `name`: unique variable identifier
- `type`: primitive type
- `ndims`: number of dimensions
- `shape`: global dimension
- `start`: local offset
- `count`: local dimension
- `constant_dims`: `true`: shape, start, count won't change; `false`:
  shape, start, count will change after definition
"""
function define_variable(io::AIO, name::AbstractString, type::Type,
                         shape::LocalValue)
    ptr = ccall((:adios2_define_variable, libadios2_c), Ptr{Cvoid},
                (Ptr{Cvoid}, Cstring, Cint, Csize_t, Ptr{Csize_t}, Ptr{Csize_t},
                 Ptr{Csize_t}, Cint), io.ptr, name, adios_type(type), 1,
                Ref(local_value_dim), C_NULL, C_NULL, true)
    return ptr == C_NULL ? nothing : Variable(ptr, io.adios)
end
function define_variable(io::AIO, name::AbstractString, type::Type,
                         shape::Union{Nothing,NTuple{N,Int} where N,
                                      CartesianIndex}=nothing,
                         start::Union{Nothing,NTuple{N,Int} where N,
                                      CartesianIndex}=nothing,
                         count::Union{Nothing,NTuple{N,Int} where N,
                                      CartesianIndex}=nothing;
                         constant_dims::Bool=false)
    ndims = max(length(maybe(shape, ())), length(maybe(start, ())),
                length(maybe(count, ())))
    @assert all(x -> x ≡ nothing || length(x) == ndims, [shape, start, count])
    ptr = ccall((:adios2_define_variable, libadios2_c), Ptr{Cvoid},
                (Ptr{Cvoid}, Cstring, Cint, Csize_t, Ptr{Csize_t}, Ptr{Csize_t},
                 Ptr{Csize_t}, Cint), io.ptr, name, adios_type(type), ndims,
                shape ≡ nothing ? C_NULL : Csize_t[reverse(Tuple(shape))...],
                start ≡ nothing ? C_NULL : Csize_t[reverse(Tuple(start))...],
                count ≡ nothing ? C_NULL : Csize_t[reverse(Tuple(count))...],
                constant_dims)
    return ptr == C_NULL ? nothing : Variable(ptr, io.adios)
end

function define_variable(io::AIO, name::AbstractString, var::AdiosType)
    return define_variable(io, name, typeof(var), LocalValue())
end
function define_variable(io::AIO, name::AbstractString,
                         arr::AbstractArray{<:AdiosType})
    return define_variable(io, name, eltype(arr), nothing, nothing, size(arr);
                           constant_dims=true)
end

export inquire_variable
"""
    variable = inquire_variable(io::AIO, name::AbstractString)
    variable::Union{Nothing,Variable}

Retrieve a variable handler within current `io` handler.

# Arguments
- `io`: handler to variable `io` owner
- `name`: unique variable identifier within `io` handler
"""
function inquire_variable(io::AIO, name::AbstractString)
    ptr = ccall((:adios2_inquire_variable, libadios2_c), Ptr{Cvoid},
                (Ptr{Cvoid}, Cstring), io.ptr, name)
    return ptr == C_NULL ? nothing : Variable(ptr, io.adios)
end

export inquire_all_variables
"""
    variables = inquire_all_variables(io::AIO)
    variables::Union{Nothing,Vector{Variable}}

Returns an array of variable handlers for all variable present in the
`io` group.

# Arguments
- `io`: handler to variables io owner
"""
function inquire_all_variables(io::AIO)
    c_variables = Ref{Ptr{Ptr{Cvoid}}}()
    size = Ref{Csize_t}()
    err = ccall((:adios2_inquire_all_variables, libadios2_c), Cint,
                (Ref{Ptr{Ptr{Cvoid}}}, Ref{Csize_t}, Ptr{Cvoid}), c_variables,
                size, io.ptr)
    Error(err) ≠ error_none && return nothing
    variables = Array{Variable}(undef, size[])
    for n in 1:length(variables)
        ptr = unsafe_load(c_variables[], n)
        @assert ptr ≠ C_NULL
        variables[n] = Variable(ptr, io.adios)
    end
    free(c_variables[])
    return variables
end

export inquire_group_variables
"""
    vars = inquire_group_variables(io::AIO, full_prefix::AbstractString)
    vars::Vector{String}

List all variables in the group `full_prefix`.
"""
function inquire_group_variables(io::AIO, full_prefix::AbstractString)
    c_variables = Ref{Ptr{Ptr{Cvoid}}}(C_NULL)
    size = Ref{Csize_t}()
    err = ccall((:adios2_inquire_group_variables, libadios2_c), Cint,
                (Ref{Ptr{Ptr{Cvoid}}}, Cstring, Ref{Csize_t}, Ptr{Cvoid}),
                c_variables, full_prefix, size, io.ptr)
    Error(err) ≠ error_none && return nothing
    variables = Array{Variable}(undef, size[])
    for n in 1:length(variables)
        ptr = unsafe_load(c_variables[], n)
        @assert ptr ≠ C_NULL
        variables[n] = Variable(ptr, io.adios)
    end
    free(c_variables[])
    return variables
end

export define_attribute
"""
    attribute = define_attribute(io::AIO, name::AbstractString, value)
    attribute::Union{Nothing,Attribute}

Define an attribute value inside `io`.
"""
function define_attribute(io::AIO, name::AbstractString, value::AdiosType)
    T = typeof(value)
    if T <: AbstractString
        ptr = ccall((:adios2_define_attribute, libadios2_c), Ptr{Cvoid},
                    (Ptr{Cvoid}, Cstring, Cint, Cstring), io.ptr, name,
                    adios_type(T), value)
    else
        ptr = ccall((:adios2_define_attribute, libadios2_c), Ptr{Cvoid},
                    (Ptr{Cvoid}, Cstring, Cint, Ptr{Cvoid}), io.ptr, name,
                    adios_type(T), Ref(value))
    end
    return ptr == C_NULL ? nothing : Attribute(ptr, io.adios)
end

export define_attribute_array
"""
    attribute = define_attribute_array(io::AIO, name::AbstractString,
                                       values::AbstractVector)
    attribute::Union{Nothing,Attribute}

Define an attribute array inside `io`.
"""
function define_attribute_array(io::AIO, name::AbstractString,
                                values::AbstractVector{<:AdiosType})
    T = eltype(values)
    if T <: AbstractString
        cvalues = pointer.(values)
        ptr = ccall((:adios2_define_attribute_array, libadios2_c), Ptr{Cvoid},
                    (Ptr{Cvoid}, Cstring, Cint, Ptr{Ptr{Cchar}}, Csize_t),
                    io.ptr, name, adios_type(T), cvalues, length(values))
        # Use `values` again to ensure it is not GCed too early
        cvalues = pointer.(values)
    else
        ptr = ccall((:adios2_define_attribute_array, libadios2_c), Ptr{Cvoid},
                    (Ptr{Cvoid}, Cstring, Cint, Ptr{Cvoid}, Csize_t), io.ptr,
                    name, adios_type(T), values, length(values))
    end
    return ptr == C_NULL ? nothing : Attribute(ptr, io.adios)
end

export define_variable_attribute
"""
    attribute = define_variable_attribute(io::AIO, name::AbstractString, value,
                                          variable_name::AbstractString,
                                          separator::AbstractString="/")
    attribute::Union{Nothing,Attribute}

Define an attribute single value associated to an existing variable by
its name.
"""
function define_variable_attribute(io::AIO, name::AbstractString,
                                   value::AdiosType,
                                   variable_name::AbstractString,
                                   separator::AbstractString="/")
    T = typeof(value)
    if T <: AbstractString
        ptr = ccall((:adios2_define_variable_attribute, libadios2_c),
                    Ptr{Cvoid},
                    (Ptr{Cvoid}, Cstring, Cint, Cstring, Cstring, Cstring),
                    io.ptr, name, adios_type(T), value, variable_name,
                    separator)
    else
        ptr = ccall((:adios2_define_variable_attribute, libadios2_c),
                    Ptr{Cvoid},
                    (Ptr{Cvoid}, Cstring, Cint, Ptr{Cvoid}, Cstring, Cstring),
                    io.ptr, name, adios_type(T), Ref(value), variable_name,
                    separator)
    end
    return ptr == C_NULL ? nothing : Attribute(ptr, io.adios)
end

export define_variable_attribute_array
"""
    attribute = define_variable_attribute_array(io::AIO, name::AbstractString,
                                                values::AbstractVector,
                                                variable_name::AbstractString,
                                                separator::AbstractString="/")
    attribute::Union{Nothing,Attribute}

Define an attribute array associated to an existing variable by its
name.
"""
function define_variable_attribute_array(io::AIO, name::AbstractString,
                                         values::AbstractVector{<:AdiosType},
                                         variable_name::AbstractString,
                                         separator::AbstractString="/")
    T = eltype(values)
    if T <: AbstractString
        cvalues = pointer.(values)
        ptr = ccall((:adios2_define_variable_attribute_array, libadios2_c),
                    Ptr{Cvoid},
                    (Ptr{Cvoid}, Cstring, Cint, Ptr{Ptr{Cchar}}, Csize_t,
                     Cstring, Cstring), io.ptr, name, adios_type(T), cvalues,
                    length(values), variable_name, separator)
        # Use `values` again to ensure it is not GCed too early
        cvalues = pointer.(values)
    else
        ptr = ccall((:adios2_define_variable_attribute_array, libadios2_c),
                    Ptr{Cvoid},
                    (Ptr{Cvoid}, Cstring, Cint, Ptr{Cvoid}, Csize_t, Cstring,
                     Cstring), io.ptr, name, adios_type(T), values,
                    length(values), variable_name, separator)
    end
    return ptr == C_NULL ? nothing : Attribute(ptr, io.adios)
end

export inquire_attribute
"""
    attribute = inquire_attribute(io::AIO, name::AbstractString)
    attribute::Union{Nothing,Attribute}

Return a handler to a previously defined attribute by name.
"""
function inquire_attribute(io::AIO, name::AbstractString)
    ptr = ccall((:adios2_inquire_attribute, libadios2_c), Ptr{Cvoid},
                (Ptr{Cvoid}, Cstring), io.ptr, name)
    return ptr == C_NULL ? nothing : Attribute(ptr, io.adios)
end

export inquire_variable_attribute
"""
    attribute = inquire_variable_attribute(io::AIO, name::AbstractString,
                                           variable_name::AbstractString,
                                           separator::AbstractString="/")
    attribute::Union{Nothing,Attribute}

Return a handler to a previously defined attribute by name.
"""
function inquire_variable_attribute(io::AIO, name::AbstractString,
                                    variable_name::AbstractString,
                                    separator::AbstractString="/")
    ptr = ccall((:adios2_inquire_variable_attribute, libadios2_c), Ptr{Cvoid},
                (Ptr{Cvoid}, Cstring, Cstring, Cstring), io.ptr, name,
                variable_name, separator)
    return ptr == C_NULL ? nothing : Attribute(ptr, io.adios)
end

export inquire_all_attributes
"""
    attributes = inquire_all_attributes(io::AIO)
    attributes::Union{Nothing,Vector{Attribute}}

Return an array of attribute handlers for all attribute present in the
io group.
"""
function inquire_all_attributes(io::AIO)
    c_attributes = Ref{Ptr{Ptr{Cvoid}}}()
    size = Ref{Csize_t}()
    err = ccall((:adios2_inquire_all_attributes, libadios2_c), Cint,
                (Ref{Ptr{Ptr{Cvoid}}}, Ref{Csize_t}, Ptr{Cvoid}), c_attributes,
                size, io.ptr)
    Error(err) ≠ error_none && return nothing
    attributes = Array{Attribute}(undef, size[])
    for n in 1:length(attributes)
        ptr = unsafe_load(c_attributes[], n)
        @assert ptr ≠ C_NULL
        attributes[n] = Attribute(ptr, io.adios)
    end
    free(c_attributes[])
    return attributes
end

export inquire_group_attributes
"""
    vars = inquire_group_attributes(io::AIO, full_prefix::AbstractString)
    vars::Vector{String}

List all attributes in the group `full_prefix`.
"""
function inquire_group_attributes(io::AIO, full_prefix::AbstractString)
    c_attributes = Ref{Ptr{Ptr{Cvoid}}}(C_NULL)
    size = Ref{Csize_t}()
    err = ccall((:adios2_inquire_group_attributes, libadios2_c), Cint,
                (Ref{Ptr{Ptr{Cvoid}}}, Cstring, Ref{Csize_t}, Ptr{Cvoid}),
                c_attributes, full_prefix, size, io.ptr)
    Error(err) ≠ error_none && return nothing
    attributes = Array{Attribute}(undef, size[])
    for n in 1:length(attributes)
        ptr = unsafe_load(c_attributes[], n)
        @assert ptr ≠ C_NULL
        attributes[n] = Attribute(ptr, io.adios)
    end
    free(c_attributes[])
    return attributes
end

export inquire_subgroups
"""
    groups = inquire_subgroups(io::AIO, full_prefix::AbstractString)
    groups::Vector{String}

List all subgroups in the group `full_prefix`.
"""
function inquire_subgroups(io::AIO, full_prefix::AbstractString)
    c_subgroups = Ref{Ptr{Ptr{Cchar}}}()
    size = Ref{Csize_t}()
    err = ccall((:adios2_inquire_subgroups, libadios2_c), Cint,
                (Ref{Ptr{Ptr{Cchar}}}, Cstring, Ref{Csize_t}, Ptr{Cvoid}),
                c_subgroups, full_prefix, size, io.ptr)
    Error(err) ≠ error_none && return nothing
    subgroups = Array{String}(undef, size[])
    for n in 1:length(subgroups)
        ptr = unsafe_load(c_subgroups[], n)
        @assert ptr ≠ C_NULL
        subgroups[n] = unsafe_string(ptr)
        free(ptr)
    end
    free(c_subgroups[])
    return subgroups
end

"""
    engine = open(io::AIO, name::AbstractString, mode::Mode)
    engine::Union{Nothing,Engine}

 Open an Engine to start heavy-weight input/output operations.

In MPI version reuses the communicator from [`adios_init_mpi`](@ref).
MPI Collective function as it calls `MPI_Comm_dup`.

# Arguments
- `io`: engine owner
- `name`: unique engine identifier
- `mode`: `mode_write`, `mode_read`, `mode_append` (not yet supported)
"""
function Base.open(io::AIO, name::AbstractString, mode::Mode)
    ptr = ccall((:adios2_open, libadios2_c), Ptr{Cvoid},
                (Ptr{Cvoid}, Cstring, Cint), io.ptr, name, mode)
    return ptr == C_NULL ? nothing : Engine(ptr, io.adios)
end

export engine_type
"""
    type = engine_type(io::AIO)
    type::Union{Nothing,String}

Return engine type string.
"""
function engine_type(io::AIO)
    size = Ref{Csize_t}()
    err = ccall((:adios2_engine_type, libadios2_c), Cint,
                (Ptr{Cchar}, Ref{Csize_t}, Ptr{Cvoid}), C_NULL, size, io.ptr)
    Error(err) ≠ error_none && return nothing
    type = '\0'^size[]
    err = ccall((:adios2_engine_type, libadios2_c), Cint,
                (Ptr{Cchar}, Ref{Csize_t}, Ptr{Cvoid}), type, size, io.ptr)
    Error(err) ≠ error_none && return nothing
    return type
end

export get_engine
"""
    engine = get_engine(io::AIO, name::AbstractString)
    engine::Union{Nothing,Engine}
"""
function get_engine(io::AIO, name::AbstractString)
    ptr = ccall((:adios2_get_engin, libadios2_c), Ptr{Cvoid},
                (Ptr{Cvoid}, Cstring), io.ptr, name)
    return ptr == C_NULL ? nothing : Engine(ptr, io.adios)
end


text:
"""
    MIPVerify(optimizer)

MIPVerify computes maximum allowable disturbance using mixed integer linear programming.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: hyperrectangle
3. Output: halfspace

# Return
`AdversarialResult`

# Method
MILP encoding. Use presolve to compute a tight node-wise bounds first.
Default `optimizer` is `GLPKSolverMIP()`.

# Property
Sound and complete.

# Reference

V. Tjeng, K. Xiao, and R. Tedrake,
["Evaluating Robustness of Neural Networks with Mixed Integer Programming,"
*ArXiv Preprint ArXiv:1711.07356*, 2017.](https://arxiv.org/abs/1711.07356)

[https://github.com/vtjeng/MIPVerify.jl](https://github.com/vtjeng/MIPVerify.jl)
"""
@with_kw struct MIPVerify{O<:AbstractMathProgSolver}
    optimizer::O
end

function solve(solver::MIPVerify, problem::Problem)
    model = Model(solver = solver.optimizer)
    neurons = init_neurons(model, problem.network)
    deltas = init_deltas(model, problem.network)
    add_complementary_set_constraint!(model, problem.output, last(neurons))
    bounds = get_bounds(problem)
    encode_mip_constraint!(model, problem.network, bounds, neurons, deltas)
    o = max_disturbance!(model, first(neurons) - problem.input.center)
    status = solve(model, suppress_warnings = true)
    if status == :Infeasible
        return AdversarialResult(:SAT)
    end
    if getvalue(o) >= minimum(problem.input.radius)
        return AdversarialResult(:SAT)
    else
        return AdversarialResult(:UNSAT, getvalue(o))
    end
end


text:
abstract type Shape{T} end

mutable struct Mesh{T} <: Shape{T}
    bodyids::Vector{Int64}
    m::T
    J::SMatrix{3,3,T,9}

    xoffset::SVector{3,T}
    qoffset::UnitQuaternion{T}

    path::String
    color::RGBA

    function Mesh(path::String, m::Real, J::AbstractMatrix;
            color = RGBA(0.75, 0.75, 0.75), xoffset::AbstractVector = zeros(3), qoffset::UnitQuaternion = one(UnitQuaternion)
        )
        T = promote_type(eltype.((m, J, xoffset, qoffset))...)

        bodyids = Int64[]
        new{T}(bodyids, m, J, xoffset, qoffset, path, color)
    end
end

mutable struct Box{T} <: Shape{T}
    bodyids::Vector{Int64}
    m::T
    J::SMatrix{3,3,T,9}

    xoffset::SVector{3,T}
    qoffset::UnitQuaternion{T}

    xyz::SVector{3,T}
    color::RGBA

    function Box(x::Real, y::Real, z::Real, m::Real;
            color = RGBA(0.75, 0.75, 0.75), xoffset::AbstractVector = zeros(3), qoffset::UnitQuaternion = one(UnitQuaternion)
        )
        T = promote_type(eltype.((x, y, z, m, xoffset, qoffset))...)

        J = 1 / 12 * m * diagm([y^2 + z^2;x^2 + z^2;x^2 + y^2])
        bodyids = Int64[]

        new{T}(bodyids, m, J, xoffset, qoffset, [x;y;z], color)
    end
end

mutable struct Cylinder{T} <: Shape{T}
    bodyids::Vector{Int64}
    m::T
    J::SMatrix{3,3,T,9}

    xoffset::SVector{3,T}
    qoffset::UnitQuaternion{T}

    rh::SVector{2,T}
    color::RGBA

    # Cylinder points in the z direction
    function Cylinder(r::Real, h::Real, m::Real;
            color = RGBA(0.75, 0.75, 0.75), xoffset::AbstractVector = zeros(3), qoffset::UnitQuaternion = one(UnitQuaternion)
        )
        T = promote_type(eltype.((r, h, m, xoffset, qoffset))...)

        J = 1 / 2 * m * diagm([r^2 + 1 / 6 * h^2;r^2 + 1 / 6 * h^2;r^2])
        bodyids = Int64[]

        new{T}(bodyids, m, J, xoffset, qoffset, [r;h], color)
    end
end

mutable struct Sphere{T} <: Shape{T}
    bodyids::Vector{Int64}
    m::T
    J::SMatrix{3,3,T,9}

    xoffset::SVector{3,T}
    qoffset::UnitQuaternion{T}

    r::T
    color::RGBA

    function Sphere(r::Real, m::Real;
            color = RGBA(0.75, 0.75, 0.75), xoffset::AbstractVector = zeros(3), qoffset::UnitQuaternion = one(UnitQuaternion)
        )
        T = promote_type(eltype.((r, m, xoffset, qoffset))...)

        J = 2 / 5 * m * diagm([r^2 for i = 1:3])
        bodyids = Int64[]

        new{T}(bodyids, m, J, xoffset, qoffset, r, color)
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, shape::Shape{T}) where {T}
    summary(io, shape)
    println(io,"")
    println(io," bodyids: "*string(shape.bodyids))
    println(io," m:       "*string(shape.m))
    println(io," J:       "*string(shape.J))
    println(io," xoffset: "*string(shape.xoffset))
    println(io," qoffset: "*string(shape.qoffset))
    println(io," color:   "*string(shape.color))
end


text:
function check_error(err)
    if err[] != C_NULL
        funcname = unsafe_string(rs2_get_failed_function(err[]))
        funcargs = unsafe_string(rs2_get_failed_args(err[]))
        message = unsafe_string(rs2_get_error_message(err[]))
        println("rs_error was raised when calling $funcname($funcargs):")
        println("    $message")
    end
end

function print_device_info(dev)
    err = Ref{Ptr{rs2_error}}(0)
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_NAME, err)
    check_error(err)
    println("\nUsing device 0, an $(unsafe_string(device_info))")
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_SERIAL_NUMBER, err)
    check_error(err)
    println("    Serial number: $(unsafe_string(device_info))")
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_FIRMWARE_VERSION, err)
    check_error(err)
    println("    Firmware version: $(unsafe_string(device_info))")
end


text:
using Random

try
    using Primes
catch
    import Pkg
    Pkg.add("Primes")
end

function rndpw2(pw2)
    @assert pw2>=5
    r = rand()
    s = 0.5
        if pw2==5
            r = r/2
            s = 0
        end
    pw2 = pw2 + r - s
    b = Int(round(2^pw2))
end

function rndf3(n)
    while true
        n1 = rndpw2(n)
        fc = factor(Vector, n1)
        length(fc)>=3 && return fc
    end
end

function sizes2(vi)
    vi = shuffle(vi)
    l = length(vi)
    dv2 = rand(1:l-1)
    vi1 = vi[begin:dv2]
    vi2 = vi[dv2+1:end]
    return (prod(vi1), prod(vi2))
end

function sizes3(vi)
    vi = shuffle(vi)
    l = length(vi)
    while true
        dv2 = rand(1:l-1)
        dv3 = rand(1:l-1)
        dv2 != dv3 && break
    end
    if dv2>dv3
        dv2, dv3 = dv3, dv2
    end
    vi1 = vi[begin:dv2]
    vi2 = vi[dv2+1:dv3]
    vi3 = vi[dv3+1:end]
    return (prod(vi1), prod(vi2), prod(vi3))
end

# function testfunc(p; disp=true)
#     n = rndf3(p)
#     k = prod(n)
#     if disp
#         println(k, " ", n)
#     end
#     s2 = sizes2(n)
#     s3 = sizes3(n)
#     @assert prod(s2)==k
#     @assert prod(s3)==k
# end
#
# function testloop(it)
#     for i in 1:it
#         p = rand(5:25)
#         testfunc(p, disp=false)
#     end
# end
"""
    sizes123(; pw2)

Fancy way to make up 1D, 2D and 3D array sizes such that their
number of elements is equal. N of elements is random and approx 2^(pw2+-0.5).
Factorization of N into 2 and 3 terms is random as well. s1=N, s2&s3 are
corresponding factorizations as tuples.
"""
function sizes123(; pw2)
    n = rndf3(pw2)
    s1 = prod(n)
    s2 = sizes2(n)
    s3 = sizes3(n)
    @assert s1==prod(s2)==prod(s3)
    return (;sizes=(;s1, s2, s3))
end

#-----------------

function loopbck2(; returnarrs=true, testarr)

    bigarrs = (; testarr)

    el_type = string(eltype(testarr))
    arsize = size(testarr)

    ret = (; el_type, arsize)
    if returnarrs
        ret = merge((; bigarrs), ret)
    end
    return ret
end

fns = (; sizes123, loopbck2)


text:
## Exercise 5-2
## The function time returns the current Greenwich Mean Time in seconds since “the epoch”, which is an arbitrary time used as a reference point. On UNIX systems, the epoch is 1 January 1970.

## julia> time()
## 1.554236060271472e9

## Write a script that reads the current time and converts it to a time of day in hours, minutes, and seconds, plus the number of days since the epoch.
println("Ans: ")
using Printf

curr_time = time()

seconds = curr_time % 60
h = 60 * 60
minutes = ((curr_time - seconds) % (h))/h
d = 60 * 60 * 24
hours = ((curr_time - (seconds + minutes)) % (d)) / d
days = (curr_time - (seconds + minutes + hours)) / (60 * 60 * 24)

@printf("%d Days, %d Hours, %d Minutes, %d Seconds\n", days, hours, minutes, seconds)

println("End.")


text:
#jl #! format: off
# Using PowerSystems to calculate network matrices

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# PowerSystems.jl supports the calculation of several different matrix representations of
# power system networks. This example demonstrates how to use PowerSystems.jl to calculate:
#  - Y bus
#  - Power transfer distribution factor (PTDF)
#  - Line outage distribution  factor (LODF)

# ### Dependencies
# Let's use a dataset from the [tabular data parsing example](https://nbviewer.jupyter.org/github/NREL-SIIP/SIIPExamples.jl/blob/master/notebook/2_PowerSystems_examples/parse_matpower.ipynb)
using SIIPExamples

pkgpath = pkgdir(SIIPExamples)
include(joinpath(pkgpath, "test", "2_PowerSystems_examples", "02_parse_matpower.jl"))

# ### Ybus
ybus = Ybus(sys)

# ### PTDF
ptdf = PTDF(sys)

# ### LODF
lodf = LODF(sys)

# ### Indexing
# Note that the axes of these matrices that correspond to buses are indexed by bus number
# (::Int64) while the branch axes are indexed by branch name (::String). You can access
# specific elements of the matrices as follows:

ptdf["bus3-bus4-i_6", 3]

# Additionally, PowerSystems provides accessors to the network matrices that take `Componets`
# as arguments so that you can pass references to the components themselves rather than the
# name or number. For example:
buses = collect(get_components(Bus, sys))
ybus[buses[1], buses[2]]

# If you would instead like to index by bus name, something like the following should work:
busname2num = get_components(Bus, sys) |> (c -> Dict(zip(get_name.(c), get_number.(c))))
ptdf["bus3-bus4-i_6", busname2num["bus3"]]

max_stars_repo_name: NREL-SIIP/Examples

max_issues_repo_name: NREL-SIIP/Examples

max_forks_repo_name: NREL-SIIP/Examples


text:
"""
    WSVarScoreTestBaseObs

A base per-observation object for the score test of within-subject variance
linear mixed model data instance without information on X1 or W1.
Contains base variables for testing
H0: β1 = 0 and τ1 = 0, H1: β1 ≠ 0 or τ1 ≠ 0,
for the full model of WiSER (with parameters β = [β1, β2],
     τ = [τ1, τ2], and Lγ),
We make use of the fitted null model.
"""
struct WSVarScoreTestBaseObs{T <: BlasReal}
    nullObs       :: WSVarLmmObs{T} # is it redundant?
    n :: Int
    p :: Int
    q :: Int
    l :: Int

    # working arrays
    A_21_β2β1_pre   :: Matrix{T} # becomes A_21_β2β1 when right-multiplied by X1
    A_21_τ2τ1_pre   :: Matrix{T} # becomes A_21_τ2τ1 when right-multiplied by W1
    A_21_Lγτ1_pre   :: Matrix{T} # becomes A_21_Lγτ1 when right-multiplied by W1

    A_21_β2β1_rowsum :: Vector{T} # for time-invariant X1, length p
    A_21_τ2τ1_rowsum :: Vector{T} # for time-invariant W1, length l
    A_21_Lγτ1_rowsum :: Vector{T} # for time-invariant W1, length q◺
end

function WSVarScoreTestBaseObs(nullObs::WSVarLmmObs{T}) where T <: BlasReal
    n, p, q, l = size(nullObs.Xt, 2), size(nullObs.Xt, 1),
        size(nullObs.Zt, 1), size(nullObs.Wt, 1)
    q◺ = ◺(q)
    A_21_β2β1_pre = Matrix{T}(undef, p, n)
    A_21_τ2τ1_pre = Matrix{T}(undef, l, n)
    A_21_Lγτ1_pre = Matrix{T}(undef, q◺, n)

    # construct A_21_β2β1_pre
    # X2t * Vinv, Vinv = Dinv - UUt
    X2t_U = nullObs.Xt * transpose(nullObs.Ut)
    mul!(A_21_β2β1_pre, -X2t_U, nullObs.Ut)
    @inbounds @simd for j in 1:n
        for i in 1:p
            A_21_β2β1_pre[i, j] += nullObs.Xt[i, j] *
                nullObs.Dinv[j] # first term
        end
    end

    # construct A_21_τ2τ1_pre
    # W2t * D * Vinv .* Vinv * D, this is no longer symmetric.
    # compute W2t * D * Vinv .* Vinv
    mul!(A_21_τ2τ1_pre, nullObs.Wt_D_Ut_kr_Utt, nullObs.Ut_kr_Ut) # third term
    @inbounds @simd for j in 1:n
        for i in 1:l
            # first term
            A_21_τ2τ1_pre[i, j] += nullObs.Wt_D_Dinv[i, j] * nullObs.Dinv[j]
            # second term
            A_21_τ2τ1_pre[i, j] += -2 * nullObs.Wt_D_sqrtdiagDinv_UUt[i, j] *
                nullObs.sqrtDinv_UUt[j]
        end
    end
    # right-multiply by D.
    @inbounds @simd for j in 1:n
        for i in 1:l
            A_21_τ2τ1_pre[i, j] = A_21_τ2τ1_pre[i, j] * nullObs.expwτ[j]
        end
    end

    # construct A_21_Lγτ1_pre
    # 2 * Cq' * (L'Z'(V^-1) ⊙ Z'(V^-1)) * Diagonal(expwτ)
    # nullObs.storage_q◺n is always Cq' * (L'Z'(V^-1) ⊙ Z'(V^-1)).
    A_21_Lγτ1_pre .=  2 * nullObs.storage_q◺n
    @inbounds @simd for j in 1:n
        for i in 1:q◺
            A_21_Lγτ1_pre[i, j] = A_21_Lγτ1_pre[i, j] * nullObs.expwτ[j]
        end
    end

    A_21_β2β1_rowsum = reshape(sum(A_21_β2β1_pre; dims=2), :)
    A_21_τ2τ1_rowsum = reshape(sum(A_21_τ2τ1_pre; dims=2), :)
    A_21_Lγτ1_rowsum = reshape(sum(A_21_Lγτ1_pre; dims=2), :)

    WSVarScoreTestBaseObs{T}(nullObs, n, p, q, l, A_21_β2β1_pre, A_21_τ2τ1_pre,
        A_21_Lγτ1_pre, A_21_β2β1_rowsum, A_21_τ2τ1_rowsum, A_21_Lγτ1_rowsum)
end

"""
    WSVarScoreTestObs

A base per-observation object for the score test of
within-subject variance linear mixed model data instance.
H0: β1 = 0 and τ1 = 0, H1: β1 ≠ 0 or τ1 ≠ 0,
for the full model of WiSER (with parameters β = [β1, β2], τ = [τ1, τ2], and Lγ).
We make use of the fitted null model.
"""
struct WSVarScoreTestObs{T <: BlasReal}
    # data
    testbaseobs         :: WSVarScoreTestBaseObs{T}
    r_X1                :: Int  # number of test variables in X1
    r_W1                :: Int  # number of test variables in W1
    r                   :: Int  # number of total test variables, = r_X1 + r_W1.
    #X1t                 :: Matrix{T}        # test variables for X
    #W1t                 :: Matrix{T}        # test variables for W

    # working arrays
    ψ_1                 :: Vector{T}
    ψ_β1                :: AbstractVector{T}
    ψ_τ1                :: AbstractVector{T}
    # ψ_β2 = testbaseobs.nullObs.∇β
    # ψ_τ2 = testbaseobs.nullObs.∇τ
    # ψ_Lγ = vech(testbaseobs.nullObs.∇Lγ),
    # be consistent with order of variables in WiSER.sandwich!()

    A_21_β2β1   :: Matrix{T} # p x r_X1, testbaseobs.A_21_β2β1_pre * X1
    A_21_τ2τ1   :: Matrix{T} # l x r_W1, testbaseobs.A_21_τ2τ1_pre * W1
    A_21_Lγτ1   :: Matrix{T} # q◺ x r_W1, testbaseobs.A_21_Lγτ1_pre * W1
end

function WSVarScoreTestObs(testbaseobs::WSVarScoreTestBaseObs{T},
    r_X1, r_W1
    #X1obs::AbstractMatrix{T}, W1obs::AbstractMatrix{T}
) where {T <: BlasReal}
    n, p, q, l = testbaseobs.n, testbaseobs.p, testbaseobs.q, testbaseobs.l
    q◺ = ◺(q)
    #r_X1 = size(X1obs, 2)
    #r_W1 = size(W1obs, 2)
    r = r_X1 + r_W1
    #X1t = Matrix{T}(undef, r_X1, n)
    #W1t = Matrix{T}(undef, r_W1, n)

    ψ_1 = Vector{T}(undef, r)
    ψ_β1 = @view(ψ_1[1:r_X1])
    ψ_τ1 = @view(ψ_1[r_X1+1:end])

    A_21_β2β1 = Matrix{T}(undef, p, r_X1)
    A_21_τ2τ1 = Matrix{T}(undef, l, r_W1)
    A_21_Lγτ1 = Matrix{T}(undef, q◺, r_W1)

    WSVarScoreTestObs{T}(testbaseobs, r_X1, r_W1, r, #X1t, W1t,
        ψ_1, ψ_β1, ψ_τ1, A_21_β2β1, A_21_τ2τ1, A_21_Lγτ1)
end

"""
    update!(testobs::WSVarScoreTestObs, X1obs, W1obs)

Updates WSVarScoreTestObs based on the test data.
"""
function update!(testobs::WSVarScoreTestObs,
    X1obs::Union{Nothing, AbstractVecOrMat{T}},
    W1obs::Union{Nothing, AbstractVecOrMat{T}}
) where {T <: BlasReal}
    testbaseobs = testobs.testbaseobs
    if testobs.r_X1 > 0
        mul!(testobs.ψ_β1,
            transpose(X1obs),
            reshape(testbaseobs.nullObs.Dinv_r -
                transpose(testbaseobs.nullObs.rt_UUt), :
            )
        )
        mul!(testobs.A_21_β2β1, testbaseobs.A_21_β2β1_pre, X1obs)
    end
    if testobs.r_W1 > 0
        mul!(testobs.ψ_τ1, -transpose(W1obs), testbaseobs.nullObs.diagDVRV)
        mul!(testobs.A_21_τ2τ1, testbaseobs.A_21_τ2τ1_pre, W1obs)
        mul!(testobs.A_21_Lγτ1, testbaseobs.A_21_Lγτ1_pre, W1obs)
    end
end

"""
    WSVarScoreTest

A base object for the score test of within-subject vLMM data instance.
H0: β1 = 0 and τ1 = 0, H1: β1 ≠ 0 or τ1 ≠ 0,
for the full model of WiSER (with parameters β = [β1, β2], τ = [τ1, τ2], and Lγ).
We make use of the fitted null model.
"""
struct WSVarScoreTest{T <: BlasReal}
    nullmodel   :: WSVarLmmModel{T}
    testobsvec :: Vector{WSVarScoreTestObs{T}}

    # dimensions
    p          :: Int       # #mean parameters in linear regression
    q          :: Int       # #random effects
    l          :: Int       # #parameters for modeling WS variability
    m          :: Int       # #individuals/clusters
    nsum       :: Int       # #observations (summed across individuals)
    r_X1                :: Int      # #test variables in X1
    r_W1                :: Int      # #test variables in W1
    r                   :: Int      # #total test variables, = r_X1 + r_W1.

    # working arrays
    ψ_1             :: Vector{T}        # length-r vector, sum_i testobs.ψ_1
    ψ_1obs          :: Matrix{T}        # r x m
    ψ_2obs          :: Matrix{T}        # (p + l + q◺) x m
    B_11            :: Matrix{T}        # r x r matrix.
    B_21            :: Matrix{T}        # (p + l + q◺) x r matrix.
    A_21            :: Matrix{T}        # (p + l + q◺) x r matrix.
    # B_22 = nullmodel.B
    # inv(A_22) = nullmodel.Ainv
    AinvBAinv   :: AbstractMatrix{T} # r x r.
    tmp_sr      :: AbstractMatrix{T} # p + l + q◺ x r.
    tmp_srx1    :: AbstractMatrix{T} # p + l + q◺ x r_X1.
    tmp_srw1    :: AbstractMatrix{T} # p + l + q◺ x r_W1.
    tmp_rr      :: AbstractMatrix{T} # r x r.
    tmp_rx1rx1  :: AbstractMatrix{T} # r_X1 x r_X1.
    tmp_rw1rw1  :: AbstractMatrix{T} # r_W1 x r_W1.
    tmp_r       :: AbstractVector{T}
    tmp_rx1     :: AbstractVector{T}
    tmp_rw1     :: AbstractVector{T}
end

function WSVarScoreTest(nullmodel::WSVarLmmModel{T},
    r_X1::Int, r_W1::Int;
    testbaseobsvec::Union{Vector{WSVarScoreTestBaseObs{T}}, Nothing} = nothing
) where {T <: BlasReal}

    @assert r_X1 >= 0 && r_W1 >= 0
    @assert r_X1 > 0 || r_W1 > 0
    if testbaseobsvec === nothing
        testbaseobsvec = [WSVarScoreTestBaseObs(obs) for obs in nullmodel.data]
    end
    testobsvec = [WSVarScoreTestObs(testbaseobs, r_X1, r_W1) for
        testbaseobs in testbaseobsvec]

    @assert nullmodel.isfitted[1] "Please fit the model first."

    p, q, l, m, nsum = nullmodel.p, nullmodel.q,
                        nullmodel.l, nullmodel.m, nullmodel.nsum
    q◺ = ◺(q)
    r = r_X1 + r_W1

    ψ_1 = Vector{T}(undef, r)
    ψ_1obs = Matrix{T}(undef, r, m)
    ψ_2obs = Matrix{T}(undef, p + l + q◺, m)
    fill!(ψ_2obs, zero(T))
    for (i, obs) in enumerate(nullmodel.data)
        ψ_β2 = @view ψ_2obs[1:p, i]
        ψ_τ2 = @view ψ_2obs[(p + 1):(p + l), i]
        ψ_Lγ = @view ψ_2obs[(p + l + 1):end, i]
        ψ_β2 .= obs.∇β
        ψ_τ2 .= obs.∇τ
        offset = 1
        @inbounds for j in 1:q, i in j:q
            ψ_Lγ[offset] = obs.∇Lγ[i, j]
            offset += 1
        end
    end

    B_11 = Matrix{T}(undef, r, r)
    B_21 = Matrix{T}(undef, p + l + q◺, r)
    A_21 = Matrix{T}(undef, p + l + q◺, r)

    AinvBAinv = nullmodel.Ainv * nullmodel.B * nullmodel.Ainv
    tmp_sr      = Matrix{T}(undef, p + l + q◺, r)
    tmp_srx1    = Matrix{T}(undef, p + l + q◺, r_X1)
    tmp_srw1    = Matrix{T}(undef, p + l + q◺, r_W1)
    tmp_rr      = Matrix{T}(undef, r, r)
    tmp_rx1rx1  = Matrix{T}(undef, r_X1, r_X1)
    tmp_rw1rw1  = Matrix{T}(undef, r_W1, r_W1)
    tmp_r       = Vector{T}(undef, r)
    tmp_rx1     = Vector{T}(undef, r_X1)
    tmp_rw1     = Vector{T}(undef, r_W1)

    WSVarScoreTest{T}(nullmodel, testobsvec, p, q, l, m, nsum, r_X1, r_W1, r,
        ψ_1, ψ_1obs, ψ_2obs, B_11, B_21, A_21,
        AinvBAinv, tmp_sr, tmp_srx1, tmp_srw1, tmp_rr, tmp_rx1rx1, tmp_rw1rw1,
        tmp_r, tmp_rx1, tmp_rw1
    )
end

"""
    test!(st::WSVarScoreTest, X1Vec, W1Vec)

Performs the score test, and returns the three p-values on the
time-variant test data
"""
function test!(st::WSVarScoreTest,
    X1vec::Union{Nothing, Vector{<:AbstractVecOrMat{T}}},
    W1vec::Union{Nothing, Vector{<:AbstractVecOrMat{T}}}) where {T <: BlasReal}

    if X1vec === nothing
        @assert st.r_X1 == 0
        X1vec = repeat([nothing], length(st.testobsvec))
    elseif W1vec === nothing
        @assert st.r_W1 == 0
        W1vec = repeat([nothing], length(st.testobsvec))
    end

    nm = st.nullmodel
    p, q, l, m, nsum = nm.p, nm.q, nm.l, nm.m, nm.nsum
    r_X1, r_W1, r = st.r_X1, st.r_W1, st.r

    # update testobsvec first
    for (testobs, X1obs, W1obs) in zip(st.testobsvec, X1vec, W1vec)
        if r_X1 > 0
            @assert size(X1obs, 1) == testobs.testbaseobs.n
            @assert size(X1obs, 2) == st.r_X1
        end
        if r_W1 > 0
            @assert size(W1obs, 1) == testobs.testbaseobs.n
            @assert size(W1obs, 2) == st.r_W1
        end
        update!(testobs, X1obs, W1obs)
    end

    # build ψ_1: sum_i testobs.ψ_1
    fill!(st.ψ_1, zero(T))
    for (i, testobs) in enumerate(st.testobsvec)
        st.ψ_1obs[:, i] .= testobs.ψ_1
        st.ψ_1 .+= testobs.ψ_1
    end

    # build B_11: using BLAS.syr!()
    fill!(st.B_11, zero(T))
    BLAS.syrk!('U', 'N', one(T), st.ψ_1obs, zero(T), st.B_11)
    copytri!(st.B_11, 'U')
    lmul!(one(T) / m, st.B_11)

    # build B_21
    mul!(st.B_21, st.ψ_2obs, transpose(st.ψ_1obs), one(T) / m, zero(T))

    # build A_21: 1/m sum_i Ai_21.
    fill!(st.A_21, zero(T))
    A_21_β2β1 = @view st.A_21[1         : p    , 1        : r_X1]
    A_21_τ2τ1 = @view st.A_21[p + 1     : p + l, r_X1 + 1 : r]
    A_21_Lγτ1 = @view st.A_21[p + l + 1 : end  , r_X1 + 1 : r]
    if r_X1 > 0
        for testobs in st.testobsvec
            A_21_β2β1 .+= testobs.A_21_β2β1
        end
    end
    if r_W1 > 0
        for testobs in st.testobsvec
            A_21_τ2τ1 .+= testobs.A_21_τ2τ1
            A_21_Lγτ1 .+= testobs.A_21_Lγτ1
        end
    end
    lmul!(1 / m, st.A_21)

    pvalues!(st)
end


text:
# This is an GCV extension for automatic rules transition with Gaussian Nodes

@rule NormalMeanVariance(:μ, Marginalisation) (m_out::ExponentialLinearQuadratic, m_v::PointMass, meta::Any) = begin
    m_out_mean, m_out_var = mean_var(m_out)
    return @call_rule NormalMeanVariance(:μ, Marginalisation) (m_out = NormalMeanVariance(m_out_mean, m_out_var), m_v = m_v, meta = meta)
end

@rule NormalMeanVariance(:μ, Marginalisation) (m_out::ExponentialLinearQuadratic, q_v::Any, meta::Any) = begin
    m_out_mean, m_out_var = mean_var(m_out)
    return @call_rule NormalMeanVariance(:μ, Marginalisation) (m_out = NormalMeanVariance(m_out_mean, m_out_var), q_v = q_v, meta = meta)
end

@rule NormalMeanPrecision(:μ, Marginalisation) (m_out::ExponentialLinearQuadratic, m_τ::PointMass, meta::Any) = begin
    m_out_mean, m_out_var = mean_var(m_out)
    return @call_rule NormalMeanPrecision(:μ, Marginalisation) (m_out = NormalMeanVariance(m_out_mean, m_out_var), m_τ = m_τ, meta = meta)
end

@rule NormalMeanPrecision(:μ, Marginalisation) (m_out::ExponentialLinearQuadratic, q_τ::Any, meta::Any) = begin
    m_out_mean, m_out_var = mean_var(m_out)
    return @call_rule NormalMeanPrecision(:μ, Marginalisation) (m_out = NormalMeanVariance(m_out_mean, m_out_var), q_τ = q_τ, meta = meta)
end

@marginalrule NormalMeanVariance(:out_μ_v) (m_out::ExponentialLinearQuadratic, m_μ::UnivariateNormalDistributionsFamily, m_v::PointMass, meta::Any) = begin
    m_out_mean, m_out_var = mean_var(m_out)
    return @call_marginalrule NormalMeanVariance(:out_μ_v) (m_out = NormalMeanVariance(m_out_mean, m_out_var), m_μ = m_μ, m_v = m_v, meta = meta)
end

@marginalrule NormalMeanVariance(:out_μ) (m_out::ExponentialLinearQuadratic, m_μ::UnivariateNormalDistributionsFamily, q_v::Any, meta::Any) = begin
    m_out_mean, m_out_var = mean_var(m_out)
    return @call_marginalrule NormalMeanVariance(:out_μ) (m_out = NormalMeanVariance(m_out_mean, m_out_var), m_μ = m_μ, q_v = q_v, meta = meta)
end

@marginalrule NormalMeanPrecision(:out_μ_τ) (m_out::ExponentialLinearQuadratic, m_μ::UnivariateNormalDistributionsFamily, m_τ::Any, meta::Any) = begin
    m_out_mean, m_out_var = mean_var(m_out)
    return @call_marginalrule NormalMeanPrecision(:out_μ_τ) (m_out = NormalMeanVariance(m_out_mean, m_out_var), m_μ = m_μ, m_τ = m_τ, meta = meta)
end

@marginalrule NormalMeanPrecision(:out_μ) (m_out::ExponentialLinearQuadratic, m_μ::UnivariateNormalDistributionsFamily, q_τ::Any, meta::Any) = begin
    m_out_mean, m_out_var = mean_var(m_out)
    return @call_marginalrule NormalMeanPrecision(:out_μ) (m_out = NormalMeanVariance(m_out_mean, m_out_var), m_μ = m_μ, q_τ = q_τ, meta = meta)
end


text:
using Markdown
using JSServe, Observables
using JSServe: Session, evaljs, linkjs
using JSServe: @js_str, onjs, Button, Slider, Asset
using WGLMakie, AbstractPlotting
using JSServe.DOM

app = App() do
    cmap_button = Button("change colormap")
    algorithm_button = Button("change algorithm")
    algorithms = ["mip", "iso", "absorption"]
    algorithm = Observable(first(algorithms))
    dropdown_onchange = js"JSServe.update_obs($algorithm, this.options[this.selectedIndex].text);"
    algorithm_drop = DOM.select(DOM.option.(algorithms); class="bandpass-dropdown", onclick=dropdown_onchange)

    data_slider = Slider(LinRange(1f0, 10f0, 100))
    iso_value = Slider(LinRange(0f0, 1f0, 100))
    N = 100
    slice_idx = Slider(1:N)

    signal = map(Observables.async_latest(data_slider.value)) do α
        a = -1; b = 2
        r = LinRange(-2, 2, N)
        z = ((x,y) -> x + y).(r, r') ./ 5
        me = [z .* sin.(α .* (atan.(y ./ x) .+ z.^2 .+ pi .* (x .> 0))) for x=r, y=r, z=r]
        return me .* (me .> z .* 0.25)
    end

    slice = map(signal, slice_idx) do x, idx
        view(x, :, idx, :)
    end
    fig = Figure()

    vol = volume(fig[1,1], signal; algorithm=map(Symbol, algorithm), ambient=Vec3f0(0.8), isovalue=iso_value)

    colormaps = collect(AbstractPlotting.all_gradient_names)
    cmap = map(cmap_button) do click
        return colormaps[rand(1:length(colormaps))]
    end

    heat = heatmap(fig[1, 2], slice, colormap=cmap)

    dom = md"""
    # More MD

    [Github-flavored Markdown info page](http://github.github.com/github-flavored-markdown/)

    [![Build Status](https://travis-ci.com/SimonDanisch/JSServe.jl.svg?branch=master)](https://travis-ci.com/SimonDanisch/JSServe.jl)

    Thoughtful example
    ======

    Alt-H2
    ------

    *italic* or **bold**

    Combined emphasis with **asterisks and _underscores_**.

    1. First ordered list item
    2. Another item
        * Unordered sub-list.
    1. Actual numbers don't matter, just that it's a number
        1. Ordered sub-list

    * Unordered list can use asterisks

    Inline `code` has `back-ticks around` it.
    ```julia
    test("haha")
    ```

    ---
    # JSServe

    [![Build Status](https://travis-ci.com/SimonDanisch/JSServe.jl.svg?branch=master)](https://travis-ci.com/SimonDanisch/JSServe.jl)
    [![Build Status](https://ci.appveyor.com/api/projects/status/github/SimonDanisch/JSServe.jl?svg=true)](https://ci.appveyor.com/project/SimonDanisch/JSServe-jl)
    [![Codecov](https://codecov.io/gh/SimonDanisch/JSServe.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/SimonDanisch/JSServe.jl)
    [![Build Status](https://travis-ci.com/SimonDanisch/JSServe.jl.svg?branch=master)](https://travis-ci.com/SimonDanisch/JSServe.jl)

    | Tables        | Are           | Cool  |
    | ------------- |:-------------:| -----:|
    | col 3 is      | right-aligned | $1600 |
    | col 2 is      | centered      |   $12 |
    | zebra stripes | are neat      |    $1 |

    > Blockquotes are very handy in email to emulate reply text.
    > This line is part of the same quote.

    # Plots:

    $(DOM.div("data param", data_slider))

    $(DOM.div("iso value", iso_value))

    $(DOM.div("y slice", slice_idx))

    $(algorithm_drop)

    $(cmap_button)

    ---

    $(fig.scene)

    ---
    """
    return JSServe.DOM.div(JSServe.MarkdownCSS, JSServe.Styling, dom)
end

display(app)