text:
struct Dihedral

   i :: Int64
   j :: Int64
   k :: Int64
   l :: Int64

   mult :: Int64
   kchi :: Vector{Float64}
   n :: Vector{Int64}
   delta :: Vector{Float64}

end


text:
using LossFunctions

lossFuncVec = [SquaredError(), AbsoluteError(), MinkowskiError(), SquaredPropError(), AbsolutePropError(), SquaredLogError(), AbsoluteLogError(), QLIKE()]
posInputIdent = ["squaredLogError", "absoluteLogError", "QLIKE"]

#Deterministic tests (don't bother for expressions involving logs)
loss(2, 5, SquaredError()) != 9 && error("Invalid det 0")
loss(2, 5, AbsoluteError()) != 3 && error("Invalid det 0")
loss(2, 5, MinkowskiError()) != 9 && error("Invalid det 0")
loss(8, 2, SquaredPropError()) != 9 && error("Invalid det 0")
loss(8, 2, AbsolutePropError()) != 3 && error("Invalid det 0")
#Loop over loss functions and simulate input for each type
N = 10
M = 5
for k = 1:length(lossFuncVec)
	#k = 6
	#Get current loss function
	lF = lossFuncVec[k]
	#Simulate data
	x = randn()
	y = randn()
	xVec = randn(N)
	yVec = randn(N)
	xMat = randn(N, M)
	yMat = randn(N, M)
	if any(string(lF) .== posInputIdent)
		x = abs(x)
		y = abs(y)
		xVec = abs(xVec)
		yVec = abs(yVec)
		xMat = abs(xMat)
		yMat = abs(yMat)
	end
	#Test loss on simulated numbers (note, reverse order tested)
	xyOut = loss(x, y, lF)
	typeof(xyOut) != typeof(x) && error("Invalid sim 0")
	xyOut = loss(xVec, yVec, lF)
	length(xyOut) != N && error("Invalid sim 1")
	xyOut = loss(xMat, yMat, lF)
	size(xyOut, 1) != N && error("Invalid sim 2")
	size(xyOut, 2) != M && error("Invalid sim 3")
	xyOut = loss(xVec, yMat, lF)
	size(xyOut, 1) != N && error("Invalid sim 4")
	size(xyOut, 2) != M && error("Invalid sim 5")
	xyOut = loss(xMat, yVec, lF)
	size(xyOut, 1) != N && error("Invalid sim 6")
	size(xyOut, 2) != M && error("Invalid sim 7")
end

ext: jl


text:
export divergence

function divergence{T}(X::DirectProductSpace{T})
    x = Space{T}[divergence(s) for s in X.factors]
    DirectProductSpace(x)
end

"""
    divergence(x)

Compute the divergence of a finite element space.
"""
function divergence(x::Space)
    ref = refspace(x)
    geo = geometry(x)
    els = elements(geo)
    fns = x.fns
    dvs = similar(fns)
    for (i,fn) in enumerate(fns)
        dvs[i] = similar(fns[i])
        for (j,sh) in enumerate(fn)
            el = els[sh.cellid]
            dvs[i][j] = divergence(ref, sh, el)
        end
    end
    divergence(x, geo, dvs)
end

"""
    curl(X)

Compute the curl of a finite element basis. The resulting set of functions might be linearly dependent because of the kernel of the curl operator.
"""
function curl(x::Space)
    ref = refspace(x)
    geo = geometry(x)
    els = elements(geo)
    crl = similar(x.fns)
    for (i,fn) in enumerate(x.fns)
        crl[i] = similar(x.fns[i], 0)
        for (j,sh) in enumerate(fn)
            el = els[sh.cellid]
            append!(crl[i], curl(ref, sh, el))
        end
    end
    curl(x, geo, crl)
end


text:
using MeshViz
using Meshes
using GeoTables
using ReferenceTests
using ImageIO
using Test

import CairoMakie as Mke

@testset "MeshViz.jl" begin
  # 1D point set
  p = PointSet(rand(Point1, 100))
  # TODO

  # 2D point set
  p = PointSet(rand(Point2, 100))
  viz(p)
  viz(p, color = :red)
  viz(p, color = 1:100)
  viz(p, color = 1:100, colormap = :inferno)

  # 3D point set
  p = PointSet(rand(Point3, 100))
  viz(p)
  viz(p, color = :red)
  viz(p, color = 1:100)
  viz(p, color = 1:100, colormap = :inferno)

  # 1D Cartesian grid
  d = CartesianGrid(10)
  # TODO

  # 2D Cartesian grid
  d = CartesianGrid(10,10)
  viz(d)
  viz(d, showfacets = true)
  viz(d, showfacets = true, facetcolor = :red)
  viz(d, color = 1:100)
  viz(d, color = 1:100, colormap = :inferno)

  # 3D Cartesian grid
  d = CartesianGrid(10,10,10)
  viz(d)
  viz(d, showfacets = true)
  viz(d, showfacets = true, facetcolor = :red)
  viz(d, color = 1:1000)
  viz(d, color = 1:1000, colormap = :inferno)

  # 2D N-gons
  t = Triangle((1.,0.), (2.,0.), (2.,1.))
  q = Quadrangle((0.,0.), (1.,0.), (1.,1.), (0.,1.))
  viz(t)
  viz(t, showboundary = false)
  viz(t, color = :orange)
  viz(t, color = :cyan, boundarycolor = :red)
  viz(q)
  viz(q, showboundary = false)
  viz(q, color = :orange)
  viz(q, color = :cyan, boundarycolor = :red)

  # 3D N-gons
  o = Octagon([(0.0,0.0,1.0), (0.5,-0.5,0.0), (1.0,0.0,0.0), (1.5,0.5,-0.5),
               (1.0,1.0,0.0), (0.5,1.5,0.0), (0.0,1.0,0.0), (-0.5,0.5,0.0)])
  viz(o)
  viz(o, showboundary = false)
  viz(o, color = :orange)
  viz(o, color = :cyan, boundarycolor = :red)

  # Polygonal areas
  p = PolyArea((0.,0.), (0.5,-1.5), (1.,0.), (1.5,0.5),
               (1.,1.), (0.5,1.5), (-0.5,0.5), (0.,0.))
  viz(p)
  viz(p, showboundary = false)
  viz(p, color = :orange)
  viz(p, color = :cyan, boundarycolor = :red)

  # Multi-geometries
  t = Triangle((1.,0.), (2.,0.), (2.,1.))
  q = Quadrangle((0.,0.), (1.,0.), (1.,1.), (0.,1.))
  m = Multi([t, q])
  viz(m)
  viz(m, showboundary = false)
  viz(m, color = :orange)
  viz(m, color = :cyan, boundarycolor = :red)

  # 2D boxes
  b = Box((0.,0.), (1.,1.))
  viz(b)
  viz(b, showboundary = false)
  viz(b, color = :orange)
  viz(b, color = :cyan, boundarycolor = :red)

  # 3D boxes
  b = Box((0.,0.,0.), (1.,1.,1.))
  viz(b)
  viz(b, color = :orange)

  # 3D spheres
  s = Sphere((0.,0.,0.), 1.)
  viz(s)
  viz(s, color = :orange)

  # Collections of geometries
  t = Triangle((1.,0.), (2.,0.), (2.,1.))
  q = Quadrangle((0.,0.), (1.,0.), (1.,1.), (0.,1.))
  d = Collection([t, q])
  viz(d)
  viz(d, showfacets = true)
  viz(d, showfacets = true, facetcolor = :red)
  viz(d, color = 1:2)
  viz(d, color = 1:2, colormap = :inferno)

  # Simple meshes
  s = Sphere((0.,0.,0.), 1.)
  m = discretize(s, RegularDiscretization(10))
  viz(m)
  viz(m, color = 1:nvertices(m))
  viz(m, color = 1:nelements(m))
  viz(m, color = 1:nelements(m), showfacets = true)
  viz(m, color = :orange, showfacets = true, facetcolor = :cyan)

  # 2D partitions
  g = CartesianGrid(10, 10)
  p = partition(g, PlanePartition((1.,1.)))
  viz(p)

  # 3D partitions
  g = CartesianGrid(10, 10, 10)
  p = partition(g, PlanePartition((1.,1.,1.)))
  viz(p)

  # Vector of points
  p = rand(Point2, 100)
  viz(p)
  viz(p, color = 1:100)
  viz(p, color = 1:100, colormap = :inferno)
  viz(p, size = 20)

  # Data over grid
  d = meshdata(CartesianGrid(10,10), etable = (z=rand(100),w=1:100))
  viz(d, variable = :z)
  viz(d, variable = :w)
  viz(d, colormap = :inferno)

  # Data over point set
  d = meshdata(PointSet(rand(Point2,100)), etable = (z=rand(100),))
  viz(d)
  viz(d, colormap = :inferno)
end


text:
module unit_cube_mode_examples
using FinEtools
using FinEtools.MeshExportModule

function unit_cube_modes()
    println("""
    Vibration modes of unit cube  of almost incompressible material.
    %
    Reference: Puso MA, Solberg J (2006) A stabilized nodally integrated
    tetrahedral. International Journal for Numerical Methods in
    Engineering 67: 841-867.
    """)
    t0 = time()

    E = 1*phun("PA");
    nu = 0.499;
    rho = 1*phun("KG/M^3");
    a = 1*phun("M"); b = a; h =  a;
    n1 = 10;# How many element edges per side?
    na =  n1; nb =  n1; nh  = n1;
    neigvs = 20                   # how many eigenvalues
    OmegaShift = (0.01*2*pi)^2;

    MR = DeforModelRed3D
    fens,fes  = H20block(a,b,h, na,nb,nh)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(3,2)), material)

    @time K  = stiffness(femm, geom, u)
    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(3,3)), material)
    @time M = mass(femm, geom, u)
    d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
    d = d - OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    mode = 17
    scattersysvec!(u, v[:,mode])
    File =  "unit_cube_modes.vtk"
    vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])

    true

end # unit_cube_modes

function unit_cube_modes_algo()
    println("""
    % Vibration modes of unit cube  of almost incompressible material.
    %
    % Reference: Puso MA, Solberg J (2006) A stabilized nodally integrated
    % tetrahedral. International Journal for Numerical Methods in
    % Engineering 67: 841-867.""")
    t0 = time()

    E = 1*phun("PA");
    nu = 0.499;
    rho= 1*phun("KG/M^3");
    a=1*phun("M"); b=a; h= a;
    n1=2;# How many element edges per side?
    na= n1; nb= n1; nh =n1;
    neigvs=20                   # how many eigenvalues
    omega_shift=(0.1*2*pi)^2;

    fens,fes =H20block(a,b,h, na,nb,nh)

    # Make the region
    MR = DeforModelRed3D
    material = MatDeforElastIso(MR, rho, E, nu, 0.0)
    region1 = FDataDict("femm"=>FEMMDeforLinear(MR, IntegData(fes, GaussRule(3,2)),
    material), "femm_mass"=>FEMMDeforLinear(MR, IntegData(fes, GaussRule(3,3)),
    material))

    # Make model data
    modeldata =  FDataDict(
    "fens"=> fens, "regions"=>  [region1],
    "omega_shift"=>omega_shift, "neigvs"=>neigvs)

    # Solve
    modeldata = FinEtools.AlgoDeforLinearModule.modal(modeldata)

    fs = modeldata["omega"]/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    modeldata["postprocessing"] = FDataDict("file"=>"unit_cube_mode",
    "mode"=>10)
    modeldata=FinEtools.AlgoDeforLinearModule.exportmode(modeldata)
    @async run(`"paraview.exe" $(modeldata["postprocessing"]["file"]*"1.vtk")`)

    true

end # unit_cube_modes_algo

function unit_cube_modes_export()
    println("""
    Vibration modes of unit cube  of almost incompressible material.

    This example EXPORTS the model to Abaqus.

    Reference: Puso MA, Solberg J (2006) A stabilized nodally integrated
    tetrahedral. International Journal for Numerical Methods in
    Engineering 67: 841-867.
    """)
    t0 = time()

    E = 1*phun("PA");
    nu = 0.499;
    rho = 1*phun("KG/M^3");
    a = 1*phun("M"); b = a; h =  a;
    n1 = 5;# How many element edges per side?
    na =  n1; nb =  n1; nh  = n1;
    neigvs = 20                   # how many eigenvalues
    OmegaShift = (0.01*2*pi)^2;

    MR = DeforModelRed3D
    fens,fes  = H20block(a,b,h, na,nb,nh)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    material=MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(3,2)), material)

    K =stiffness(femm, geom, u)
    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(3,3)), material)
    M =mass(femm, geom, u)
    d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
    d = d - OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    mode = 7
    scattersysvec!(u, v[:,mode])
    File =  "unit_cube_modes.vtk"
    vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])
    @async run(`"paraview.exe" $File`)

    AE = AbaqusExporter("unit_cube_modes_h20");
    # AE.ios = STDOUT;
    HEADING(AE, "Vibration modes of unit cube  of almost incompressible material.");
    COMMENT(AE, "The  first six frequencies are rigid body modes.");
    COMMENT(AE, "The  first nonzero frequency (7) should be around 0.26 Hz");
    PART(AE, "part1");
    END_PART(AE);
    ASSEMBLY(AE, "ASSEM1");
    INSTANCE(AE, "INSTNC1", "PART1");
    NODE(AE, fens.xyz);
    COMMENT(AE, "The hybrid form of the serendipity hexahedron is chosen because");
    COMMENT(AE, "the material is  nearly incompressible.");
    ELEMENT(AE, "c3d20rh", "AllElements", 1, connasarray(fes))
    ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
    SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements");
    END_INSTANCE(AE);
    END_ASSEMBLY(AE);
    MATERIAL(AE, "elasticity")
    ELASTIC(AE, E, nu)
    DENSITY(AE, rho)
    STEP_FREQUENCY(AE, neigvs)
    END_STEP(AE)
    close(AE)

    true

end # unit_cube_modes_export

function unit_cube_modes_msh8_algo()
    println("""
    % Vibration modes of unit cube  of almost incompressible material.
    % Mean-strain hexahedron.
    % Reference: Puso MA, Solberg J (2006) A stabilized nodally integrated
    % tetrahedral. International Journal for Numerical Methods in
    % Engineering 67: 841-867.""")
    t0 = time()

    E = 1*phun("PA");
    nu = 0.499;
    rho= 1*phun("KG/M^3");
    a=1*phun("M"); b=a; h= a;
    n1=8 # How many element edges per side?
    na= n1; nb= n1; nh =n1;
    neigvs=20                   # how many eigenvalues
    omega_shift=(0.1*2*pi)^2;

    fens,fes = H8block(a,b,h, na,nb,nh)

    # Make the region
    MR = DeforModelRed3D
    material = MatDeforElastIso(MR, rho, E, nu, 0.0)
    region1 = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegData(fes, GaussRule(3,2)),
    material), "femm_mass"=>FEMMDeforLinearMSH8(MR, IntegData(fes, GaussRule(3,3)),
    material))

    # Make model data
    modeldata =  FDataDict(
    "fens"=> fens, "regions"=>  [region1],
    "omega_shift"=>omega_shift, "neigvs"=>neigvs)

    # Solve
    modeldata = FinEtools.AlgoDeforLinearModule.modal(modeldata)

    fs = modeldata["omega"]/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    modeldata["postprocessing"] = FDataDict("file"=>"unit_cube_mode",
    "mode"=>10)
    modeldata=FinEtools.AlgoDeforLinearModule.exportmode(modeldata)
    @async run(`"paraview.exe" $(modeldata["postprocessing"]["file"]*"1.vtk")`)

    true

end # unit_cube_modes_msh8_algo

function allrun()
    println("#####################################################")
    println("# unit_cube_modes ")
    unit_cube_modes()
    println("#####################################################")
    println("# unit_cube_modes_algo ")
    unit_cube_modes_algo()
    println("#####################################################")
    println("# unit_cube_modes_export ")
    unit_cube_modes_export()
    println("#####################################################")
    println("# unit_cube_modes_msh8_algo ")
    unit_cube_modes_msh8_algo()
    return true
end # function allrun

end # module unit_cube_mode_examples


text:
@testset "Geometries" begin
  @testset "Parametric dimension" begin
    @test paramdim(Segment) == 1
    @test paramdim(Triangle) == 2
    @test paramdim(Quadrangle) == 2
    @test paramdim(PolyArea) == 2
    @test paramdim(Hexahedron) == 3
    @test paramdim(Pyramid) == 3
    @test paramdim(Tetrahedron) == 3
  end
end


text:
module TestModularityUtilities
using Test
using EcologicalNetworks

# Perfectly modular bipartite network
A = [
  true true true false false false; true true true false false false;
  false false false true true true; false false false true true true
  ]
B = BipartiteNetwork(A)
U = convert(UnipartiteNetwork, B)

LB = Dict(zip(species(B), [1,1,2,2,1,1,1,2,2,2]))
LU = Dict(zip(species(U), [1,1,2,2,1,1,1,2,2,2]))

@test Q(B, LB) == 0.5
@test Q(U, LU) == 0.5
@test Qr(B, LB) == 1.0
@test Qr(U, LU) == 1.0

end


text:
# This file was generated, do not modify it. # hide
ℓ = symlogdensity(post)
f = codegen(post; ℓ=ℓ)
fast = @benchmark $f($(;x,λ), $(;y), $trace)

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
abstract type AbstractMissile <: AbstractEnv end

include("PointMass3DMissile.jl")
include("PursuerEvador3DMissile.jl")


text:
using PressureCoeffs
using Test

@testset "PressureCoeffs.jl" begin
    # Write your tests here.
end


text:
using NonuniformResampling1D
using Test
using Statistics
using Interpolations

@testset "NonuniformResampling1D Tests" begin
    include("range_utilities.test.jl")
    include("window_functions.test.jl")
    include("nuresample.test.jl")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ f13a03a0-24e6-11ec-2825-8123926abe5d
begin
	import Pkg
	Pkg.activate(joinpath(@__DIR__))
	using SimLib, Plots, PlutoUI, Statistics, LinearAlgebra, LaTeXStrings, XXZNumerics, SpinSymmetry, LsqFit, SparseArrays
	gr()
	LinearAlgebra.BLAS.set_num_threads(8)
	#cd(joinpath(@__DIR__, ".."))
end

# ╔═╡ 3eef0c5c-9ef5-4e0e-ab9b-0721e68f8413
using BenchmarkTools

# ╔═╡ 25cbd399-a199-4469-be1d-f402df522511
md"""
# Entanglement Entropy implementation
## Basics
Have density matrix $\rho = \sum_i p_i |i\rangle\langle i|$ on full system. Can divide system in subsystems $A$ and $B$.

### Partial trace
```math
\rho_B = \mathrm{Tr}_A \rho = \sum_{j} \sum_i \langle a_i |\langle b_j| \rho |a_i \rangle b_j \rangle |b_j\rangle\langle b_j| = \sum_{i,j} \langle a_i b_j| \rho |a_i b_j \rangle |b_j\rangle\langle b_j|
```
where $\{a_i\}$ ($\{b_i\}$) are a basis of $A$ ($B$).

### Von-Neumann entropy:
```math
S(\rho) = \mathrm{Tr} \rho \log_2 \rho = \sum_i p_i \log_2 p_i
```

### Entanglement entropy
Entanglement entropy between subsystem $A$ and $B$ is
```math
S_{A} = S(\rho_A) = S(\mathrm{Tr}_B \rho) = S(\rho_B) = S_B
```
→ Need to compute partial trace and compute eigenvalues.

→ Better way for pure states: SVD
```math
|\psi\rangle = \sum_{i,j} c_{i,j} |a_i,b_j\rangle \overset{SVD}{=} \sum_i \sigma_i |\tilde{a}_i\rangle|\tilde{b}_i\rangle
```
Here we basically *reshaped* the state vector into a matrix and used singular value decomposition on it. This allows easy computation of the eigenvalues of the reduced density matrix. Consider:
```math
\begin{align}
\rho_B = \mathrm{Tr}_A \rho &= \mathrm{Tr}_A \sum_{i,j} \sigma_i \sigma_j^\star |\tilde{a}_i\rangle|\tilde{b}_i\rangle \langle\tilde{a}_j |\langle\tilde{b}_j|\\
&= \sum_j |\sigma_j|^2 |\tilde{b}_j\rangle \langle \tilde{b}_j|
\end{align}
```
**Note:** This also shows that $S(\rho_A) = S(\rho_B)$.

**Note 2:** The SVD way is computationally much faster than partial trace + eigen decomposition and saves a lot of memory as well.
"""

# ╔═╡ 3f654c7b-ff3b-423c-91ad-8c942ca5ece9
md"""
## Additional symmetries
For this computation having $|\psi\rangle$ in a subspace of a symmetry does not make the computation more straight-forward as the *reshape* operation does not work too well in all symmetrized bases.

What follows is explicitly only for spin-1/2 chains of $N$ spins but can probably be generalized more or less easily.

### z-Block basis
Let's denote the eigenspaces of the magnetization $\hat{S}_z$ by $L_k^N$ (where $k$ counts the number of $\uparrow$ in that space) and states by $|k, i\rangle$. The second index just enumerates the different basis vectors in this sector. The whole Hilbert space $\mathcal{H}$ decomposes:
```math
\mathcal{H} = \bigoplus_{k=0}^{N} L^N_k
```

Assume $|\psi\rangle \in L_k$ and thus $|\psi\rangle = \sum_i c_i |k,i\rangle$ and we want to trace out the first $N_1$ spins of the chain.

In this case we are lucky because $L^N_k$ again decompose nicely into smaller magnetization eigenspaces:
```math
L_k^N = \bigoplus_{k_1=k_{min}}^{k_{max}} L^{N_1}_{k_1} \otimes L^{N-N_1}_{k-k_1}
```
with $k_{min} = \max\{0, N_1+k-N\}$ and $k_{max} = \min\{k, N_1\}$.

This means, that the reduced density matrix is actually block diagonal! And the computation of the Schmidt values can be simplified as well as for every $k$ we again can reshape the corresponding vectors into a rectangular matrix.

```math
|\psi\rangle = \sum_i c_i |k,i\rangle = \sum_{k_1=k_{min}}^{k_{max}} \sum_{i,j} c_{i,j} |k_1,i\rangle|k-k_1,j\rangle
```

**Implementation note:**
The only thing left is to determine the correct index $i \in L^N_k$ from $(i,j) \in L^{N_1}_{k_1} \otimes L^{N-N_1}_{k-k_1}$. I currently do this via a lookup table and some index math. I already have functions to compute the indices of $L^N_k \subset \mathcal{H}$.

So what I do is:
1. Compute the indices $I$ and $J$ of $L^{N_1}_{k_1} \subset L^N_k$ and $L^{N-N_1}_{k-k_1} \subset L^N_k$
2. This gives the indices of $L^{N_1}_{k_1} \otimes L^{N-N_1}_{k-k_1} \subset L^N_k$ by $I\otimes J = \{2^{N-N_1}i+j | i\in I, j\in J\}$.
3. Thus we have the subspace's indices $I\otimes J \subset L^N_k \subset \mathcal{H}$
4. This I map via lookup table to the actual indices of the state vector.
"""

# ╔═╡ fc4f5bbb-cc20-45ed-aa1a-928ae12b1566
md"""
### Parity (spin flip)

This is not as easy although we can find a decomposition. Let $P^N_\pm$ denote the symmetry sectors. We can see:
```math
P^N_+ = P^{N_1}_+ \otimes P^{N-N_1}_+ \oplus P^{N_1}_- \otimes P^{N-N_1}_-
```

Again we have 2 blocks not talking to each other that may be treated seperately. Finding the indices is a bit more involved this time around.

Consider a state $|s\rangle \in \mathcal{H} = |s_1, s_2\rangle$ and its flipped partner $|\bar{s}\rangle = |\bar{s}_1, \bar{s}_2\rangle$. Also consider the half-flipped state $|s^\star\rangle = |s_1, \bar{s}_2\rangle$ and its partner $|\bar{s}^\star\rangle = |\bar{s}_1, s_2\rangle$. Thus we have $|S\rangle = |s\rangle + |\bar{s}\rangle \in P^N_+$ and $|S^\star\rangle = |s^\star\rangle + |\bar{s}^\star\rangle \in P^N_+$. One can directly construct:
```math
\begin{align}
|t_+\rangle &= |S\rangle + |S^\star\rangle = (|s_1\rangle + |\bar{s}_1\rangle)\otimes (|s_2\rangle + |\bar{s}_2\rangle) \in P^{N_1}_+ \otimes P^{N-N_1}_+\\
|t_-\rangle &= |S\rangle - |S^\star\rangle = (|s_1\rangle - |\bar{s}_1\rangle)\otimes (|s_2\rangle - |\bar{s}_2\rangle) \in P^{N_1}_- \otimes P^{N-N_1}_-
\end{align}
```
"""

# ╔═╡ 119fb04d-219b-4a08-a50b-584e470f517b
N = 11

# ╔═╡ 6ad0d260-5284-41f2-b7d7-273023d9114c
H = xxzmodel(Symmetric((1 .+ (2*rand(N,N))) .^ 6),-0.7)

# ╔═╡ 942f4e4c-8f21-4a60-ae78-e183e003c748
eig = eigen!(Hermitian(Matrix(H)));

# ╔═╡ a29ef3cf-b4c1-4320-b580-f3cda39e7fac
eig.vectors[:, 5]

# ╔═╡ dd2a8d83-7251-402c-81ea-7b9e52138080
size(eig.vectors,1)

# ╔═╡ 0ea97d53-929e-4683-9742-3183a527abf4
k = 5

# ╔═╡ 4539ba8a-ec44-444d-8ad0-ccc07e7443f5
abs2.(svdvals(reshape(eig.vectors[:, k], 2^div(N-1,2), 2^(N-div(N-1,2)))))

# ╔═╡ f7145441-e14b-47dc-b8c4-c7d243b38c8e
entropy(abs2.(svdvals(reshape(eig.vectors[:, k], 2^div(N-1,2), 2^(N-div(N-1,2))))))

# ╔═╡ 6472d25b-94f7-42ee-8a9d-935b31748eb9
eig.values[:,1]

# ╔═╡ 9bf89fcc-81f9-40c2-9fa5-ddd2eea47941
function halfchain_entropy(state, N=nspins(state), N1=div(N,2))
	entropy(abs2.(svdvals(reshape(state, 2^N1, 2^(N-N1)))))
end

# ╔═╡ d9b06c60-56be-4128-9422-fbfcda7f53e8
plot(halfchain_entropy.(eachcol(eig.vectors)))

# ╔═╡ 2bf00bc0-1d51-4134-afe3-9d5cf897f495
halfchain_entropy.(eachcol(eig.vectors))

# ╔═╡ 0f9e0ab5-7c6e-4e7a-9e58-2185c0b55e4f
halfchain_entropy(eig.vectors[:, k])

# ╔═╡ 1274618f-d59e-4af7-ad6a-81e9c4b1e8bd
begin
	function partial_trace(vec)
		N = round(Int, log2(length(vec)))
		Nhalf = div(N,2)
		ρ = vec * vec'
		tmp = reshape(ρ, 2^Nhalf,2^(N-Nhalf), 2^Nhalf,2^(N-Nhalf))
		res = zeros(eltype(ρ), 2^Nhalf, 2^Nhalf)
		for i in 1:2^(N-Nhalf)
			res += tmp[:,i,:,i]
		end
		Hermitian(res)
	end
end

# ╔═╡ 6910adba-dea2-4095-b888-12bf5a259fd5
reverse(eigvals(partial_trace(eig.vectors[:,k])))

# ╔═╡ efc5f35f-bb74-40ab-b335-b91a1ac7b6b6
NB = N

# ╔═╡ bf58d3eb-46fe-41b8-819a-1827bb29ab51
kB = div(N-1,2)

# ╔═╡ 18820e48-e53b-4d06-ac56-a2a958737156
N1 = 3

# ╔═╡ ae78b8ca-2a42-401d-ac1e-1a0507d015a1
symm_basis = ZBlockBasis(N,kB)

# ╔═╡ 1a6fdac3-7f6f-4c7b-b379-977268b54707
ψsymm = normalize!(ones(basissize(symm_basis)));

# ╔═╡ ae734d90-65cf-4e37-9c6c-aa725d31131c
ψfull = Vector(sparsevec(SpinSymmetry._indices(symm_basis), ψsymm, 2^NB));

# ╔═╡ 44a313c0-3c1b-4769-9176-81acbeffb017
entropy(abs2.(svdvals(reshape(ψfull, 2^N1, 2^(NB-N1)))))

# ╔═╡ 1957e200-71fe-45e3-aee5-39b8eecba31a
SZ = sum(op_list(σz, NB))

# ╔═╡ fff4c13f-003c-49f8-968a-a11a984ab8b0
diag(SZ[SpinSymmetry._indices(symm_basis), SpinSymmetry._indices(symm_basis)])

# ╔═╡ 70d4e7b4-4149-4094-ae75-9c6e7eb097ea
SpinSymmetry._indices(symm_basis)

# ╔═╡ 90000a7c-e927-4de3-b7ab-f238bd5f7af5
indexLookup = let inds = SpinSymmetry._indices(symm_basis)
	sparsevec(inds, 1:length(inds))
end

# ╔═╡ c1aac45a-ca0a-4188-af09-b082ae14cd7c
let res = []
	for k in reverse(0:min(N1,kB))
		let b1 = SpinSymmetry._indices(ZBlockBasis(N1,k)),
			b2 = SpinSymmetry._indices(ZBlockBasis(NB-N1, kB-k))
			push!(res, vec(b1 .+  2^N1 * (b2 .- 1)'))
		end
	end
	vcat(res...)
end

# ╔═╡ b121dd76-ecde-415b-8d29-e92c452b305b
let inds = SpinSymmetry._indices(symm_basis),
	v = Dict((i, ind) for (ind, i) in enumerate(inds))
	@btime $v[118]
end

# ╔═╡ db4667e3-0ed5-430b-b9e2-b245608d8d7b
v = Dict((i, ind) for (ind, i) in enumerate(SpinSymmetry._indices(symm_basis)))

# ╔═╡ 5db25eaa-8fe0-45ed-9029-a42ebc968c63
let v = Vector(indexLookup),
	k = [56,111], res = zeros(Int, length(k))
	@btime $res .= getindex.(Ref($v), $k)
end

# ╔═╡ b9928ff9-8c98-4498-b1d3-3597357a1ffd
let v = Dict((i, ind) for (ind, i) in enumerate(SpinSymmetry._indices(symm_basis))),
	k = [56,111], res = zeros(Int, length(k))
	@btime $res .= getindex.(Ref($v), $k)
end

# ╔═╡ 1e976790-2e7c-4e43-a222-a1e0d13c4f7f
let v = indexLookup,
	k = [56,111], res = zeros(Int, length(k))
	@btime $res .= getindex.(Ref($v), $k)
end

# ╔═╡ ea33e44e-65fd-46a7-b077-cfac2a0ef768
allinds =
	[SpinSymmetry._indices(ZBlockBasis(N1,k)) .+
		2^N1 .* (SpinSymmetry._indices(ZBlockBasis(NB-N1, kB-k)) .- 1)'
		 for k in max(0,N1+kB-NB):min(N1,kB)]

# ╔═╡ 3c1b94c5-2279-40bf-8b69-3b4dac8ea490
reducedρ2 = let res = [], indexLookup = Vector(indexLookup)
	for indmat in allinds
		l = size(indmat, 1)
		temp = zeros(eltype(ψsymm), l, l)
		vals = zeros(eltype(ψsymm), l)
		for inds in eachcol(indmat)
			vals .= getindex.(Ref(ψsymm), getindex.(Ref(indexLookup), inds))
			temp = mul!(temp, vals, vals', 1, 1)
		end
		push!(res, temp)
	end
	res
end

# ╔═╡ 5c262eea-0622-4046-a60d-432f24daf651
eigvals.(reducedρ2)

# ╔═╡ 7cc24ef6-a45b-4f87-a17f-91f7b196f5ad
svdvals(ψsymm[getindex.(Ref(indexLookup), allinds[2])]).^2

# ╔═╡ a84a4a13-ba18-4697-a643-817147e597bd
typeof(allinds)

# ╔═╡ 49481125-553b-48df-907b-ef98051cac47
sum(length, allinds)

# ╔═╡ 7f0db04e-5701-4978-bc8c-e98753bbb749
reducedρ = let res = []
	for k in reverse(0:min(N1,kB))
		let b1 = SpinSymmetry._indices(ZBlockBasis(N1,k)),
			b2 = 2^N1 * (SpinSymmetry._indices(ZBlockBasis(NB-N1, kB-k)) .- 1),
			temp = zeros(eltype(ψsymm), length(b2), length(b2))
			for indA in b1
				fullBasisInds = indA .+ b2
				symmBasisInds = indexLookup[fullBasisInds]
				vals = ψsymm[symmBasisInds]
				temp += vals * vals'
			end
			push!(res, temp)
		end
	end
	res
end

# ╔═╡ 2ba4c3c4-ec8a-4791-b940-4251fa9ae145
entropy(vcat(svdvals.(reducedρ)...))

# ╔═╡ 8cfdd892-56c1-49f8-a993-2f426f9ec690
function entanglement_entropy_zblock(state, N, kB, N1)
	(N == N1 || N1 == 0) && return 0
	(N == k  ||  k == 0) && return 0
	reducedρ = []
	for k in reverse(0:min(N1,kB))
		let b1 = SpinSymmetry._indices(ZBlockBasis(N1,k)),
			b2 = 2^N1 * (SpinSymmetry._indices(ZBlockBasis(NB-N1, kB-k)) .- 1),
			temp = zeros(eltype(ψsymm), length(b2), length(b2))
			for indA in b1
				fullBasisInds = indA .+ b2
				symmBasisInds = indexLookup[fullBasisInds]
				vals = state[symmBasisInds]
				temp += vals * vals'
			end
			push!(reducedρ, temp)
		end
	end
	entropy(vcat(eigvals!.(reducedρ)...))
end

# ╔═╡ 4b3461ca-73e6-4c1a-819a-53ab3774a6f4
entanglement_entropy_zblock(ψsymm, N, kB, N1)

# ╔═╡ c1343cbc-37b9-446e-9035-57983934c99b
Hsymm = symmetrize_operator(H, ZBlockBasis(N, kB))

# ╔═╡ 2cf277f4-38f9-43c6-90b4-20c1a0c5a5df
eigsymm = eigen!(Hermitian(Matrix(Hsymm)));

# ╔═╡ 98f5eff9-d71b-4a2d-ab88-a88b9fffa944
plot(entanglement_entropy_zblock.(eachcol(eigsymm.vectors), N, div(N-1,2), kB))

# ╔═╡ 79893553-b8b2-4d90-80f8-3b0d01d5096d

# ╔═╡ Cell order:
# ╠═25cbd399-a199-4469-be1d-f402df522511
# ╠═3f654c7b-ff3b-423c-91ad-8c942ca5ece9
# ╠═fc4f5bbb-cc20-45ed-aa1a-928ae12b1566
# ╠═f13a03a0-24e6-11ec-2825-8123926abe5d
# ╠═119fb04d-219b-4a08-a50b-584e470f517b
# ╠═6ad0d260-5284-41f2-b7d7-273023d9114c
# ╠═942f4e4c-8f21-4a60-ae78-e183e003c748
# ╠═a29ef3cf-b4c1-4320-b580-f3cda39e7fac
# ╠═dd2a8d83-7251-402c-81ea-7b9e52138080
# ╠═0ea97d53-929e-4683-9742-3183a527abf4
# ╠═4539ba8a-ec44-444d-8ad0-ccc07e7443f5
# ╠═f7145441-e14b-47dc-b8c4-c7d243b38c8e
# ╠═d9b06c60-56be-4128-9422-fbfcda7f53e8
# ╠═6472d25b-94f7-42ee-8a9d-935b31748eb9
# ╠═2bf00bc0-1d51-4134-afe3-9d5cf897f495
# ╠═0f9e0ab5-7c6e-4e7a-9e58-2185c0b55e4f
# ╠═9bf89fcc-81f9-40c2-9fa5-ddd2eea47941
# ╠═6910adba-dea2-4095-b888-12bf5a259fd5
# ╠═1274618f-d59e-4af7-ad6a-81e9c4b1e8bd
# ╠═efc5f35f-bb74-40ab-b335-b91a1ac7b6b6
# ╠═bf58d3eb-46fe-41b8-819a-1827bb29ab51
# ╠═18820e48-e53b-4d06-ac56-a2a958737156
# ╠═1a6fdac3-7f6f-4c7b-b379-977268b54707
# ╠═ae734d90-65cf-4e37-9c6c-aa725d31131c
# ╠═44a313c0-3c1b-4769-9176-81acbeffb017
# ╠═2ba4c3c4-ec8a-4791-b940-4251fa9ae145
# ╠═ae78b8ca-2a42-401d-ac1e-1a0507d015a1
# ╠═1957e200-71fe-45e3-aee5-39b8eecba31a
# ╠═fff4c13f-003c-49f8-968a-a11a984ab8b0
# ╠═70d4e7b4-4149-4094-ae75-9c6e7eb097ea
# ╠═90000a7c-e927-4de3-b7ab-f238bd5f7af5
# ╠═c1aac45a-ca0a-4188-af09-b082ae14cd7c
# ╠═3eef0c5c-9ef5-4e0e-ab9b-0721e68f8413
# ╠═b121dd76-ecde-415b-8d29-e92c452b305b
# ╠═db4667e3-0ed5-430b-b9e2-b245608d8d7b
# ╠═5db25eaa-8fe0-45ed-9029-a42ebc968c63
# ╠═b9928ff9-8c98-4498-b1d3-3597357a1ffd
# ╠═1e976790-2e7c-4e43-a222-a1e0d13c4f7f
# ╠═ea33e44e-65fd-46a7-b077-cfac2a0ef768
# ╠═3c1b94c5-2279-40bf-8b69-3b4dac8ea490
# ╠═5c262eea-0622-4046-a60d-432f24daf651
# ╠═7cc24ef6-a45b-4f87-a17f-91f7b196f5ad
# ╠═a84a4a13-ba18-4697-a643-817147e597bd
# ╠═49481125-553b-48df-907b-ef98051cac47
# ╠═7f0db04e-5701-4978-bc8c-e98753bbb749
# ╠═8cfdd892-56c1-49f8-a993-2f426f9ec690
# ╠═4b3461ca-73e6-4c1a-819a-53ab3774a6f4
# ╠═c1343cbc-37b9-446e-9035-57983934c99b
# ╠═2cf277f4-38f9-43c6-90b4-20c1a0c5a5df
# ╠═98f5eff9-d71b-4a2d-ab88-a88b9fffa944
# ╠═79893553-b8b2-4d90-80f8-3b0d01d5096d


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("World of machine learning")

# # IMDB Movie Review Sentiment Analysis Demo

#-

# ## Setup
# * Adds required packages to Julia.
# * Loads the data and a pretrained model.
# * See [07.imdb.ipynb](https://github.com/denizyuret/Knet.jl/blob/v1.1.1/tutorial/07.imdb.ipynb) to see how the model was trained.

cd(datapath("imdb"))
include("imdb.jl")

# ## Sample Data
# The model was trained using 25000 movie reviews such as the following (shift-ENTER to see random example)
# Each review was tokenized, lowercased, truncated to max 150 words and a 30,000 word vocabulary.

r = rand(1:length(xtrn))
println(reviewstring(xtrn[r],ytrn[r]))
flush(stdout)

# ## Test
# We test the model on 25000 never before seen reviews on the test set (shift-ENTER to see random example)
# The test accuracy is around 86%

r = rand(1:length(xtst))
println(reviewstring(xtst[r],ytst[r]))
println(predictstring(xtst[r]))
flush(stdout)

# ## User Input
# In this cell you can enter your own review and let the model guess the sentiment

userinput = readline(stdin)
words = split(lowercase(userinput))
ex = [get!(imdbdict,wr,UNK) for wr in words]
ex[ex.>VOCABSIZE].=UNK
println(predictstring(ex))

#-

max_stars_repo_path: Courses/World of machine learning/imdbdemo.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
###########################################################################
#                             Data Structures                             #
###########################################################################

# A new type specifically for spatial branch and bound, essentially SDPWRMPowerModel
mutable struct NodeWRMPowerModel <: PM.AbstractWRMModel PM.@pm_fields end

# A new branch type, mapping (Wij, Tij) to their lower and upper bounds
mutable struct ComplexVariableBranch <: BB.AbstractBranch
    lb::Dict{Tuple{JuMP.VariableRef, JuMP.VariableRef},Real} # Lij
    ub::Dict{Tuple{JuMP.VariableRef, JuMP.VariableRef},Real} # Uij
end

# This computes coefficients for cuts (3a), (3b) in the paper
function compute_π(bounds::NTuple{6, <:Real})::NTuple{5, <:Real}
    function sigmoid(x::Real)::Real
        x == 0 ? res = 0 : res = (sqrt(1 + x^2) - 1) / x
        return res
    end
    (Lii, Uii, Ljj, Ujj, Lij, Uij) = bounds
    π0 = -sqrt(Lii * Ljj * Uii * Ujj)
    π1 = -sqrt(Ljj * Ujj)
    π2 = -sqrt(Lii * Uii)
    π3 = (sqrt(Lii) + sqrt(Uii)) * (sqrt(Ljj) + sqrt(Ujj)) * (1 - sigmoid(Lij) * sigmoid(Uij)) / (1 + sigmoid(Lij) * sigmoid(Uij))
    π4 = (sqrt(Lii) + sqrt(Uii)) * (sqrt(Ljj) + sqrt(Ujj)) * (sigmoid(Lij) + sigmoid(Uij)) / (1 + sigmoid(Lij) * sigmoid(Uij))
    return (π0, π1, π2, π3, π4)
end

# This is the main branch type used in our implementation
# (i,j) should always be one of the lines in the power system
mutable struct SpatialBCBranch <: BB.AbstractBranch
    i::Int
    j::Int
    wii::JuMP.VariableRef
    wjj::JuMP.VariableRef
    wr::JuMP.VariableRef
    wi::JuMP.VariableRef
    mod_branch::Union{Nothing, ComplexVariableBranch, BB.VariableBranch} # this captures either a changed bound for Wij and Tij, or a changed bound for Wii/Wjj
    bounds::NTuple{6, <:Real} # Lii, Uii, Ljj, Ujj, Lij, Uij
    valid_ineq_coeffs::NTuple{5, <:Real} # π coefficients
end

# create a SpatialBCBranch, in which mod_branch is nothing (because the exact matrix entry to branch is not decided yet)
function create_sbc_branch(i::Int64, j::Int64, prev_node::BB.AbstractNode)
    root = find_root(prev_node)
    pm = root.auxiliary_data["PM"]
    bus_ids = PM.ids(pm, :bus)
    lookup_w_index = Dict((bi,i) for (i,bi) in enumerate(bus_ids))
    widx_i = lookup_w_index[i]
    widx_j = lookup_w_index[j]
    wii = PM.var(pm, :WR)[widx_i,widx_i]
    wjj = PM.var(pm, :WR)[widx_j,widx_j]
    wr = PM.var(pm, :WR)[widx_i,widx_j]
    wi = PM.var(pm, :WI)[widx_i,widx_j]
    bounds = get_LU_from_branches(prev_node, i, j)
    πs = compute_π(bounds)
    return SpatialBCBranch(i, j, wii, wjj, wr, wi, nothing, bounds, πs)
end

# This function tracks the branches and update the bounds for selected entries (i,j)
function get_LU_from_branches(node::BB.AbstractNode, i::Int64, j::Int64)::NTuple{6, <:Real}
    pnode = node
    LU = [NaN for i in 1:6]
    while !isnothing(pnode.parent)
        sbc_branch = pnode.branch
        mod_branch = sbc_branch.mod_branch
        if sbc_branch.i == i && sbc_branch.j == j
            for i in eachindex(LU)
                if isnan(LU[i]) LU[i] = sbc_branch.bounds[i] end
            end
        elseif sbc_branch.i == i
            if isnan(LU[1]) LU[1] = sbc_branch.bounds[1] end
            if isnan(LU[2]) LU[2] = sbc_branch.bounds[2] end
        elseif sbc_branch.i == j
            if isnan(LU[3]) LU[3] = sbc_branch.bounds[1] end
            if isnan(LU[4]) LU[4] = sbc_branch.bounds[2] end
        elseif sbc_branch.j == i
            if isnan(LU[1]) LU[1] = sbc_branch.bounds[3] end
            if isnan(LU[2]) LU[2] = sbc_branch.bounds[4] end
        elseif sbc_branch.j == j
            if isnan(LU[3]) LU[3] = sbc_branch.bounds[3] end
            if isnan(LU[4]) LU[4] = sbc_branch.bounds[4] end
        end
        pnode = pnode.parent
    end
    if isnan(LU[1]) LU[1] = pnode.auxiliary_data["Lii"][i] end
    if isnan(LU[2]) LU[2] = pnode.auxiliary_data["Uii"][i] end
    if isnan(LU[3]) LU[3] = pnode.auxiliary_data["Lii"][j] end
    if isnan(LU[4]) LU[4] = pnode.auxiliary_data["Uii"][j] end
    if isnan(LU[5]) LU[5] = pnode.auxiliary_data["Lij"][(i,j)] end
    if isnan(LU[6]) LU[6] = pnode.auxiliary_data["Uij"][(i,j)] end

    return Tuple(LU)
end

max_stars_repo_name: kibaekkim/ComplexBC

max_issues_repo_name: kibaekkim/ComplexBC

max_forks_repo_name: kibaekkim/ComplexBC


text:
using Oceananigans.Grids: underlying_left_halo_indices, underlying_right_halo_indices, topology
using Oceananigans: location

"""
    zero_halo_regions!(fields::Union{Tuple, NamedTuple})

Zero out the halo regions of each field in `fields`.
"""
function zero_halo_regions!(fields::Union{Tuple, NamedTuple}; kwargs...)

    for field in fields
        zero_halo_regions!(field; kwargs...)
    end

    return nothing
end

"""
    zero_halo_regions!(field)

Zero out the halo regions of `field`.
"""
zero_halo_regions!(field; dims=(1, 2, 3)) = zero_halo_regions!(parent(field), location(field), field.grid, dims=dims)

"""
    zero_halo_regions!(underlying_data, loc, grid)

Zero out the halo regions of the `underlying_data` of a field
at `loc`ation on `grid`.
"""
function zero_halo_regions!(underlying_data, loc, grid; dims=(1, 2, 3))

    dims = tuple(dims...)
    Lx, Ly, Lz = loc
    Tx, Ty, Tz = topology(grid)

    if 1 ∈ dims
        zero_west_halo!(underlying_data,  underlying_left_halo_indices(Lx, Tx, grid.Nx, grid.Hx))
        zero_east_halo!(underlying_data, underlying_right_halo_indices(Lx, Tx, grid.Nx, grid.Hx))
    end

    if 2 ∈ dims
        zero_south_halo!(underlying_data,  underlying_left_halo_indices(Ly, Ty, grid.Ny, grid.Hy))
        zero_north_halo!(underlying_data, underlying_right_halo_indices(Ly, Ty, grid.Ny, grid.Hy))
    end

    if 3 ∈ dims
           zero_top_halo!(underlying_data,  underlying_left_halo_indices(Lz, Tz, grid.Nz, grid.Hz))
        zero_bottom_halo!(underlying_data, underlying_right_halo_indices(Lz, Tz, grid.Nz, grid.Hz))
    end

    return nothing
end

  @inline zero_west_halo!(d, halo_indices) = @views @. d[halo_indices, :, :] = 0
 @inline zero_south_halo!(d, halo_indices) = @views @. d[:, halo_indices, :] = 0
@inline zero_bottom_halo!(d, halo_indices) = @views @. d[:, :, halo_indices] = 0

 @inline zero_east_halo!(d, halo_indices) = @views @. d[halo_indices, :, :] = 0
@inline zero_north_halo!(d, halo_indices) = @views @. d[:, halo_indices, :] = 0
  @inline zero_top_halo!(d, halo_indices) = @views @. d[:, :, halo_indices] = 0


text:
using StaticArrays
using Parameters
using OrdinaryDiffEq
using Test

@inbounds @inline function ż(z, p, t)
    @unpack A, B, D = p
    p₀, p₂ = z[1:2]
    q₀, q₂ = z[3:4]

    return SVector{4}(
        -A * q₀ - 3 * B / √2 * (q₂^2 - q₀^2) - D * q₀ * (q₀^2 + q₂^2),
        -q₂ * (A + 3 * √2 * B * q₀ + D * (q₀^2 + q₂^2)),
        A * p₀,
        A * p₂
  )
end

condition(u, t, integrator) = u
affect!(integrator) = nothing
cbf(idx) = ContinuousCallback(condition,
    affect!, nothing, save_positions=(false, true), idxs=idx)
z0 = SVector{4}(7.1989885061904335, -0.165912283356219, 0., -3.63534900748947)

tspan=(0.,300.)
prob=ODEProblem(ż, z0, tspan, (A=1, B=0.55, D=0.4), callback=cbf(3))
sol=solve(prob, Vern9(), abstol=1e-14, reltol=1e-14,
    save_everystep=false, save_start=false, save_end=false, maxiters=1e6)

@test length(sol) == 126

prob=ODEProblem(ż, z0, (0,400.), (A=1, B=0.55, D=0.4), callback=cbf(3))
sol=solve(prob, Vern9(), abstol=1e-14, reltol=1e-14, save_everystep=false, save_start=false, save_end=false, maxiters=2e4)

@test length(sol) == 148

prob=ODEProblem(ż, z0, (0,5000.), (A=1, B=0.55, D=0.4), callback=cbf(3))
sol=solve(prob, Vern9(), abstol=1e-14, reltol=1e-14, save_everystep=false, save_start=false, save_end=false, maxiters=1e6)

@test length(sol) == 1883

using ParameterizedFunctions
f = @ode_def BallBounce begin
  dy =  v
  dv = -g
end g
function condition(u,t,integrator) # Event when event_f(u,t) == 0
  u[1]
end
function affect!(integrator)
  integrator.u[2] = -integrator.u[2]
end
cb2 = ContinuousCallback(condition,affect!)
tspan = (0.0,100000.0)
u0 = [50.0,0.0]
p = 9.8
prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Tsit5(),callback=cb2)
@test minimum(sol') > -40


text:
using TropicalTensors, TropicalTensors.Reversible
using Yao, TropicalYao
using Test
using Viznet
using NiLang, NiLang.AD
using Random
using TropicalTensors.Reversible: cachesize_A, cachesize_B, cachesize_largemem

@testset "spinglass" begin
    Random.seed!(2)
    lt = rand_maskedsquare(8, 10, 0.8)
    sg = Spinglass(lt, rand(Int32[-1,1], length(sgbonds(lt))), zeros(Int32, length(lt)))
    Lx, Ly = size(sg.lattice)
    reg = ArrayReg(ones(Tropical{Int32}, 1<<(Ly+2)))
    STK = stack4reg(reg, cachesize_largemem(lt))
    @test isolve_largemem(Int32(0), sg, reg, STK)[1] == solve(sg).n
    A = stack4reg(reg, cachesize_A(lt))
    B = stack4reg(reg, cachesize_B(lt))
    @test check_inv(isolve, (Int32(0), sg, reg, A, B))
    reg = ArrayReg(ones(Tropical{Int32}, 1<<(Ly+2)))
    A = stack4reg(reg, cachesize_A(lt))
    B = stack4reg(reg, cachesize_B(lt))
    @test isolve(Int32(0), sg, reg, A, B)[1] == solve(sg).n
end

@testset "optconfig" begin
    Random.seed!(5)
    sg = rand_spinglass(Int64, rand_maskedsquare(9, 7, 0.8); jt=Randpm(), ht=Zero(), seed=2)
    optc = opt_config(sg)
    optc2 = opt_config_largemem(sg)
    @test optc.grad_Js == optc2.grad_Js
    @test optc.grad_hs == optc2.grad_hs
    @test sum(sg.Js .* optc.grad_Js) == optc.eng
    @test sum(sg.Js .* optc2.grad_Js) == optc2.eng
    @test optc.eng == optc2.eng
end

@testset "fallback" begin
    sg = Spinglass(MaskedSquareLattice(ones(Bool,2,2)), [-1,1,1,1,1,-1], [0,0,0,0])
    res = opt_config_largemem(sg)
    @test res.eng == 2
end


text:
using mopkg
using Test

@testset "mopkg.jl" begin
	x, f = 1000.0, x -> x^4 - 24x^3 + 7x^2 - 11x + 420
	@test minimum_nr(f, x) == (-32524.475393806228, 17.812172762347867)
	x, f = 500.0, x -> x^4 - 7x^3 + 7x^2 - 77x + 77
	@test minimum_nr(f, x) == (-387.5919553775244, 5.277905093806142)
end


text:
using Documenter, PowerModelsDistributionStateEstimation

makedocs(
    modules     = [PowerModelsDistributionStateEstimation],
    format      = Documenter.HTML(mathengine = Documenter.MathJax()),
    sitename    = "PowerModelsDistributionStateEstimation.jl",
    authors     = "Marta Vanin, and Tom Van Acker",
    pages       = [
              "Home"    => "index.md",
              "Manual"  => [
                            "Getting Started"                         => "quickguide.md",
                            "Input Data Format"                       => "input_data_format.md",
                            "Measurements and Conversions"            => "measurements.md",
                            "State Estimation Criteria"               => "se_criteria.md",
                            "Bad Data Detection and Identification"   => "bad_data.md",
                            ],
             "Library " => [
                            "Power Flow Formulations"        => "formulations.md",
                            "Problem Specifications"         => "problems.md",
                           ],
                 ]
)

deploydocs(
     repo = "github.com/Electa-Git/PowerModelsDistributionStateEstimation.jl.git"
)


text:
mutable struct ProjMPOSum
  pm::Vector{ProjMPO}
end

ProjMPOSum(mpos::Vector{MPO}) = ProjMPOSum([ProjMPO(M) for M in mpos])

ProjMPOSum(Ms::MPO...) = ProjMPOSum([Ms...])

nsite(P::ProjMPOSum) = nsite(P.pm[1])

Base.length(P::ProjMPOSum) = length(P.pm[1])

function product(P::ProjMPOSum,
                 v::ITensor)::ITensor
  Pv = product(P.pm[1],v)
  for n=2:length(P.pm)
    Pv += product(P.pm[n],v)
  end
  return Pv
end

function Base.eltype(P::ProjMPOSum)
  elT = eltype(P.pm[1])
  for n=2:length(P.pm)
    elT = promote_type(elT,eltype(P.pm[n]))
  end
  return elT
end

(P::ProjMPOSum)(v::ITensor) = product(P,v)

Base.size(P::ProjMPOSum) = size(P.pm[1])

function position!(P::ProjMPOSum,psi::MPS,pos::Int)
  for M in P.pm
    position!(M,psi,pos)
  end
end

function noiseterm(P::ProjMPOSum,
                   phi::ITensor,
                   b::Int,
                   dir::String)
  nt = noiseterm(P.pm[1],phi,b,dir)
  for n=2:length(P.pm)
    nt += noiseterm(P.pm[n],phi,b,dir)
  end
  return nt
end


text:
# ==============================================================================
# This file fixes a some issues found in Julia.
# ==============================================================================

# ==============================================================================
# https://github.com/JuliaLang/julia/issues/33332
# Fix bug in Sparse Matrix in Julia. They used vals1[j1] != 0 instead of using
# iszero(vals1[j1]) idem for vals2[j2]

function Base.:(==)(A1::SparseMatrixCSC, A2::SparseMatrixCSC)
    size(A1) != size(A2) && return false
    vals1, vals2 = nonzeros(A1), nonzeros(A2)
    rows1, rows2 = rowvals(A1), rowvals(A2)
    m, n = size(A1)
    @inbounds for i = 1:n
        nz1,nz2 = nzrange(A1,i), nzrange(A2,i)
        j1,j2 = first(nz1), first(nz2)
        # step through the rows of both matrices at once:
        while j1 <= last(nz1) && j2 <= last(nz2)
            r1,r2 = rows1[j1], rows2[j2]
            if r1==r2
                vals1[j1]!=vals2[j2] && return false
                j1+=1
                j2+=1
            else
                if r1<r2
                    !iszero(vals1[j1]) && return false
                    j1+=1
                else
                    !iszero(vals2[j2]) && return false
                    j2+=1
                end
            end
        end
        # finish off any left-overs:
        for j = j1:last(nz1)
            !iszero(vals1[j]) && return false
        end
        for j = j2:last(nz2)
            !iszero(vals2[j]) && return false
        end
    end
    return true
end

# ==============================================================================
# https://github.com/JuliaLang/julia/issues/33332 fixed in Julia 1.3
# Fix bug of sparse[i,j] = mp0 confusing with sparse[i,j] = 0

getcolptr(S::SparseMatrixCSC) = getfield(S, :colptr)

function _insert!(v::Vector, pos::Integer, item, nz::Integer)
    if nz > length(v)
        insert!(v, pos, item)
    else # nz < length(v)
        Base.unsafe_copyto!(v, pos+1, v, pos, nz - pos)
        v[pos] = item
        v
    end
end

function _setindex_scalar!(A::SparseMatrixCSC{MP,Ti}, _v, _i::Integer, _j::Integer) where {Tv,Ti<:Integer}
    v = convert(MP, _v)
    i = convert(Ti, _i)
    j = convert(Ti, _j)
    if !((1 <= i <= size(A, 1)) & (1 <= j <= size(A, 2)))
        throw(BoundsError(A, (i,j)))
    end
    coljfirstk = Int(getcolptr(A)[j])
    coljlastk = Int(getcolptr(A)[j+1] - 1)
    searchk = searchsortedfirst(rowvals(A), i, coljfirstk, coljlastk, Base.Order.Forward)
    if searchk <= coljlastk && rowvals(A)[searchk] == i
        nonzeros(A)[searchk] = v
        return A
    end
    if !iszero(v)
        nz = getcolptr(A)[size(A, 2)+1]
        !isbitstype(Ti) || nz < typemax(Ti) ||
            throw(ArgumentError("nnz(A) going to exceed typemax(Ti) = $(typemax(Ti))"))
        _insert!(rowvals(A), searchk, i, nz)
        _insert!(nonzeros(A), searchk, v, nz)
        @simd for m in (j + 1):(size(A, 2) + 1)
            @inbounds getcolptr(A)[m] += Ti(1)
        end
    end
    return A
end

@inline Base.:setindex!(A::SparseMatrixCSC{MP,Ti}, _v, _i::Integer, _j::Integer) where {Tv,Ti<:Integer} =
    _setindex_scalar!(A, _v, _i, _j)

# Used by SimpleWeightedDiGraph
function SparseMatrixCSC{MP,Ti}(M::StridedMatrix) where {Tv,Ti}
    nz = count(!iszero, M)
    colptr = zeros(Ti, size(M, 2) + 1)
    nzval = Vector{MP}(undef, nz)
    rowval = Vector{Ti}(undef, nz)
    colptr[1] = 1
    cnt = 1
    @inbounds for j in 1:size(M, 2)
        for i in 1:size(M, 1)
            v = M[i, j]
            if !iszero(v)
                rowval[cnt] = i
                nzval[cnt] = v
                cnt += 1
            end
        end
        colptr[j+1] = cnt
    end
    return SparseMatrixCSC(size(M, 1), size(M, 2), colptr, rowval, nzval)
end

# ==============================================================================
# Since Julia 1.4.x the matrix product Max-Plus sparse * full or full * sparse
# is no longer working (while sparse * sparse keep working)
Base.:(*)(A::ArrMP, S::SpaMP) = A * full(S)
Base.:(*)(S::SpaMP, A::ArrMP) = full(S) * A
Base.:(*)(A::ArrMI, S::SpaMI) = A * full(S)
Base.:(*)(S::SpaMI, A::ArrMI) = full(S) * A

# ==============================================================================
# Since Julia 1.6.x sparse matrices are displayed like dense matrices but with
# dots instead of zero elements. We prefer the older Julia 1.5 display because
# more compact, so we force older display but only for sparse Max-Plus only and
# for sparse matrices of other type because our purpose is not to change the
# Julia behavior but we cherry pick good behavior for our toolbox.

function Base.show(io::IO, ::MIME"text/plain", S::SpaMP)
    xnnz = nnz(S)
    m, n = size(S)
    print(io, m, "×", n, " Max-Plus sparse matrix with ", xnnz, " stored ",
              xnnz == 1 ? "entry" : "entries")
    if xnnz != 0
        print(io, ":")
        show(IOContext(io, :typeinfo => eltype(S)), S)
    end
end

function Base.show(io::IO, S::SpaMP)
    Base.show(convert(IOContext, io), S::SpaMP)
end

function Base.show(io::IOContext, S::SpaMP)
    nnz(S) == 0 && return show(io, MIME("text/plain"), S)

    ioc = IOContext(io, :compact => true)
    function _format_line(r, col, padr, padc)
        print(ioc, "\n  [", rpad(rowvals(S)[r], padr), ", ", lpad(col, padc), "]  =  ")
        if isassigned(nonzeros(S), Int(r))
            show(ioc, nonzeros(S)[r])
        else
            print(ioc, Base.undef_ref_str)
        end
    end

    function _get_cols(from, to)
        idx = eltype(getcolptr(S))[]
        c = searchsortedlast(getcolptr(S), from)
        for i = from:to
            while i == getcolptr(S)[c+1]
                c +=1
            end
            push!(idx, c)
        end
        idx
    end

    rows = displaysize(io)[1] - 4 # -4 from [Prompt, header, newline after elements, new prompt]
    if !get(io, :limit, false) || rows >= nnz(S) # Will the whole matrix fit when printed?
        cols = _get_cols(1, nnz(S))
        padr, padc = ndigits.((maximum(rowvals(S)[1:nnz(S)]), cols[end]))
        _format_line.(1:nnz(S), cols, padr, padc)
    else
        if rows <= 2
            print(io, "\n  \u22ee")
            return
        end
        s1, e1 = 1, div(rows - 1, 2) # -1 accounts for \vdots
        s2, e2 = nnz(S) - (rows - 1 - e1) + 1, nnz(S)
        cols1, cols2 = _get_cols(s1, e1), _get_cols(s2, e2)
        padr = ndigits(max(maximum(rowvals(S)[s1:e1]), maximum(rowvals(S)[s2:e2])))
        padc = ndigits(cols2[end])
        _format_line.(s1:e1, cols1, padr, padc)
        print(io, "\n  \u22ee")
        _format_line.(s2:e2, cols2, padr, padc)
    end
    return
end


text:
using Base.Test

@testset "generated_weights" begin
    include("generated_weights/fc+softmax.jl")
    include("generated_weights/mfc+mfc+softmax.jl")
    include("generated_weights/conv+fc+softmax.jl")

    if Pkg.installed("Gurobi") != nothing
        # Skip these tests if Gurobi is not installed.
        include("generated_weights/conv+conv+fc+softmax.jl")
    end
end


text:
#=

  From the WebPPL model (written months ago):
  """
  Today I changed the last of two bulbs which were supposed to hold for
  2000h (they are Halogen bulbs). The first went out some day ago,
  say 25h ago.

  What is the probability of this to happen, i.e. that both bulbs went out
  in a 25h period, if we assume that the time for light bulbs are an
  Exponential distribution?
  """

  Mathematica code:
    Probability[Abs[a - b] <= c, {a, b} \[Distributed]
      ProductDistribution[ExponentialDistribution[p],
      ExponentialDistribution[p]]]

  1 - E^(-c p) (if c > 0)

  -> 0.0124222

  Here's the summary of the three models: exponential, gaussian, and poisson,
  where 1 is the probability that both bulbs went out in the range of 25h.
  I guess that the exponential model is the most realistic of these.

  Model 1: Exponential
  Distributions of variable p (num:0)
  0.00000 =>    9881  (0.988100)
  1.00000 =>     119  (0.011900)

  Model 2: Gaussian
  Distributions of variable p (num:0)
  0.00000 =>    8639  (0.863900)
  1.00000 =>    1361  (0.136100)

  Model 3: Poisson
  Distributions of variable p (num:0)
  0.00000 =>    6851  (0.685100)
  1.00000 =>    3149  (0.314900)

  Model 4: Weibull (experimental)
  Distributions of variable p (num:0)
  0.00000 =>    9880  (0.988000)
  1.00000 =>     120  (0.012000)

  Note: 2000h is 83.33 days.
  1/0.011900 (from the exponential model) ~ 84.03: This happens once each 84'th day...

  Cf ~/webppl/light_bulbs.jl

=#

using Turing
include("jl_utils.jl")

# Exponential
@model function light_bulbs_exp()
    lambda = 2000
    b1 ~ Exponential(lambda)
    b2 ~ Exponential(lambda)
    # The differences in failures
    diff ~ Dirac(abs(b1-b2))
    # Probability that the difference between the events is <= 25
    p ~ Dirac(diff <= 25)
end

println("Model 1: Exponential")
model = light_bulbs_exp()

# chns = sample(model, Prior(), 100_000)
# chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(15), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 100_000)

display(chns)
show_var_dist_pct(chns, :p)

# Gaussian distribution: 2000h, standard deviation 100h
@model function light_bulbs_gauss()
    mu = 2000;
    stdev = 100;
    b1 ~ Normal(mu,stdev)
    b2 ~ Normal(mu,stdev)
    diff ~ Dirac(abs(b1-b2))
    p ~ Dirac(diff <= 25)
end

println("Model 2: Gaussian")
model = light_bulbs_gauss()

# chns = sample(model, Prior(), 100_000)
# chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(15), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 100_000)

display(chns)
show_var_dist_pct(chns, :p)

# Poisson distribution: 1/2000h
@model function light_bulbs_poisson()
    lambda = 2000
    b1 ~ Poisson(lambda)
    b2 ~ Poisson(lambda)
    diff ~ Dirac(abs(b1-b2))
    p ~ Dirac(diff <= 25)
end

println("Model 3: Poisson")
model = light_bulbs_poisson()

# chns = sample(model, Prior(), 100_000)
# chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(15), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 100_000)

display(chns)
show_var_dist_pct(chns, :p)

# Weibul distribution
@model function light_bulbs_weibul()
    shape = 1
    scale = 2000
    b1 ~ Weibull(shape,scale)
    b2 ~ Weibull(shape,scale)
    diff ~ Dirac(abs(b1-b2))
    p ~ Dirac(diff <= 25)
end

println("Model 4: Weibull (experimental)")
model = light_bulbs_weibul()

# chns = sample(model, Prior(), 100_000)
# chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(15), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 100_000)

display(chns)
show_var_dist_pct(chns, :p)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#   This file is part of Fatou.jl. It is licensed under the MIT license
#   Copyright (C) 2019 Michael Reed

PyPlot.plot(K::FilledSet;c::String="",bare::Bool=false) = PyPlot.imshow(K;cmap=c,bare=bare)

function PyPlot.imshow(K::FilledSet;cmap::String="",bare::Bool=false)
    PyPlot.figure() # plot figure using imshow based in input preferences
    isempty(cmap) && (cmap = K.meta.cmap)
    isempty(cmap) ? PyPlot.imshow(K.meta.iter ? K.iter : K.mix, extent=bounds(K)) :
        PyPlot.imshow(K.meta.iter ? K.iter : K.mix, cmap=cmap, extent=bounds(K))
    PyPlot.tight_layout()
    !bare && PyPlot.title(K)
end

function PyPlot.title(K::FilledSet) # annotate title using LaTeX
    text,t = "f:z\\mapsto $(rdpm(Algebra.latex(K.meta.E))),\\,",LaTeXString(typeplot(K))
    if K.meta.newt # annotate y-axis with Newton's method
        PyPlot.title(latexstring("$text m = $(K.meta.m), ")*t)
        PyPlot.ylabel(L"Fatou\,set:\,"*L"z\,↦\,z-m\,×\,f(z)\,/\,f\,'(z)")
    else
        PyPlot.title(latexstring(text)*t)
    end
    PyPlot.tight_layout()
    PyPlot.colorbar()
end

function orbit(E,f::Function,bi::Matrix{Float64},orb::Int=0,depth::Int=1,incr::Int=384)
    x,N,N2,orbit,bis = real_orb(E,f,bi,orb,depth,incr)
    # prepare for next figure
    PyPlot.figure()
    # plot background lines
    PyPlot.plot(x[:],N[:,1],"k--",x[:],N[:,2])
    # plot orbit cobweb path
    PyPlot.plot(orbit[:,1],orbit[:,2],"r")
    # plot f^2,f^3,f^4,...
    for h ∈ 3:depth+1
        PyPlot.plot(x,N[:,h],lw=1)
    end
    if ~(orb == 0)
        PyPlot.plot(range(bi[1],stop=bi[2],length=length(N2)),N2[:],"gray",marker="x",linestyle=":",lw=1)
        funs = [latexstring("\\phi(x_{0:$orb})")]
        funt = ", IC: \$ x_0 = $(bis[3])\$, \$ n\\in0:$orb\$"
    else
        funs = []
        funt = ""
    end
    # trim graph
    d=1.07
    PyPlot.xlim(bi[1],bi[2])
    PyPlot.ylim(minimum([d*minimum(N[:,2]),0]),maximum([d*maximum(N[:,2]),0]))
    # set title
    fune = rdpm(Algebra.latex(E))
    PyPlot.title(latexstring("\$ x \\mapsto $fune\$$funt"))
    # set legend
    PyPlot.legend(vcat([L"$y=x$",L"$\phi(x)$",L"(x_n,\phi(x_n))"],[latexstring("\\phi^{$x}(x)") for x ∈ 2:depth],funs))
    PyPlot.tight_layout()
    return
end


text:
export
	benchmark_solve!,
	shift_fill!

function interp_rows(N::Int,tf::Float64,X::AbstractMatrix)::Matrix
    n,N1 = size(X)
    t1 = range(0,stop=tf,length=N1)
    t2 = collect(range(0,stop=tf,length=N))
    X2 = zeros(n,N)
    for i = 1:n
        interp_cubic = CubicSplineInterpolation(t1, X[i,:])
        X2[i,:] = interp_cubic(t2)
    end
    return X2
end

function ispossemidef(A)
	eigs = eigvals(A)
	if any(real(eigs) .< 0)
		return false
	else
		return true
	end
end

function convertInf!(A::VecOrMat{Float64},infbnd=1.1e20)
    infs = isinf.(A)
    A[infs] = sign.(A[infs])*infbnd
    return nothing
end

function set_logger()
    if !(global_logger() isa SolverLogger)
        global_logger(default_logger(true))
    end
end

function benchmark_solve!(solver; samples=10, evals=10)
    Z0 = deepcopy(get_trajectory(solver))
    solver.opts.verbose = false
    b = @benchmark begin
        initial_trajectory!($solver,$Z0)
        solve!($solver)
    end samples=samples evals=evals
    return b
end

function benchmark_solve!(solver, data::Dict; samples=10, evals=10)
	b = benchmark_solve!(solver, samples=samples, evals=evals)

   # Run stats
   push!(data[:time], time(median(b))*1e-6)  # ms
   push!(data[:iterations], iterations(solver))
   push!(data[:cost], cost(solver))
   return b
end

function shift_fill!(A::Vector{<:SVector}, n=1)
	N = length(A)
	@inbounds for k = n+1:N
		A[k-n] = A[k]
	end
	a_last = A[N-n]
	@inbounds for k = N-n:N
		A[k] = a_last
	end
	return nothing
end


text:
module test_stdlib_distributed

using Test
using Distributed

@test nprocs() >= 1
@test nworkers() >= 1

if myid() == 1
    @test workers() ⊆ procs()
    ws = addprocs(2)
    @test ws ⊆ workers()
    @test nprocs() == nworkers() + 1

    id = workers()[1]
    f = Future(id)
    put!(f, :OK)
    @test id > 1
    @test fetch(f) == :OK

    t = rmprocs(ws)
    wait(t)

    @test nprocs() >= 1
    @test nworkers() >= 1
end

end # module test_stdlib_distributed


text:
using Test, Factorions

@test is_factorion(145)
@test find_factorions(6) == Set([1, 2, 25, 26])


text:
read_csv(path; kws...) = read_csv(path, DataFrame; kws...)

# function read_csv(path, sink::Type{<:AbstractMatrix{T}}; delim=nothing, kws...) where T
#     x = delim === nothing ? readdlm(path, T; kws...) : readdlm(path, delim, T; kws...)
#     return x
# end

function read_csv(path, sink::Type{A}; kws...) where A <: AbstractMatrix
    A(read_csv(path; kws...))
end

function read_csv(path, sink::Type{<:DataFrame}; kws...)
    return CSV.read(path, sink; kws...)
end

function parse_pystring(s::AbstractString)
    s == "False" && return false
    s == "True" && return true
    s == "None" && return nothing
    try return parse(Int, s); catch; end
    try return parse(Float64, s); catch; end
    return s
end

function restrict_array_type(res::AbstractArray)
    # attempt conversion
    if all(x -> x isa Integer, res)
        return Int.(res)
    elseif all(x -> x isa AbstractFloat, res)
        return Float32.(res)
    elseif all(x -> x isa String, res)
        return String.(res)
    else
        return res
    end
end

function df_to_matrix(df::AbstractDataFrame)
    x = Matrix(df)
    if size(x, 2) == 1
        return reshape(x, 1, size(x, 1))
    else
        return permutedims(x, (2, 1))
    end
end

bytes_to_type(::Type{UInt8}, A::Array{UInt8}) = A
bytes_to_type(::Type{N0f8}, A::Array{UInt8}) = reinterpret(N0f8, A)
bytes_to_type(::Type{T}, A::Array{UInt8}) where T<:Integer = convert(Array{T}, A)
bytes_to_type(::Type{T}, A::Array{UInt8}) where T<:AbstractFloat = A ./ T(255)
bytes_to_type(::Type{T}, A::Array{UInt8}) where T<:Number  = convert(Array{T}, reinterpret(N0f8, A))

# PIRACY # see https://github.com/JuliaML/MLUtils.jl/issues/67
MLUtils.numobs(x::AbstractDataFrame) = size(x, 1)
MLUtils.getobs(x::AbstractDataFrame, i) = x[i, :]

function clean_nt(nt::NamedTuple)
    res = (; (p for p in  pairs(nt) if p[2] !== nothing)...)
    if isempty(res)
        return nothing
    else
        return res
    end
end

"""
    convert2image(d, i)
    convert2image(d, x)
    convert2image(DType, x)

Convert the observation(s) `i` from dataset `d` to image(s).
It can also convert a numerical array `x`.

# Examples

```julia-repl
julia> using MLDatasets: MNIST

julia> using ImageInTerminal

julia> d = MNIST()

julia> i = 1:2;

julia> convert2image(d, i) # You should see 2 images in the terminal

julia> x = d[1].features;

julia> convert2image(MNIST, x) # or convert2image(d, x)
```
"""
function convert2image end

convert2image(d::SupervisedDataset, i::Integer) =
    convert2image(typeof(d), d[i].features)
convert2image(d::SupervisedDataset, i::AbstractVector) =
    convert2image(typeof(d), d[i].features)
convert2image(d::SupervisedDataset, x::AbstractArray) =
    convert2image(typeof(d), x)


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

using DynamicHMCModels

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

ProjDir = @__DIR__
cd(ProjDir)

# Read the milk data

df = DataFrame(CSV.read(joinpath("..", "..", "data", "milk.csv"), delim=';'))
df = filter(row -> !(row[:neocortex_perc] == "NA"), df)
#df[:, :kcal_per_g] = convert(Vector{Float64}, df[:, :kcal_per_g])
df[!, :log_mass] = log.(convert(Vector{Float64}, df[:, :mass]))

# Define the model struct

Base.@kwdef mutable struct MilkModel{Ty <: AbstractVector, Tx <: AbstractMatrix}
    "Observations."
    y::Ty
    "Covariates"
    x::Tx
end

# Write a function to return properly dimensioned transformation.

function make_transformation(model::MilkModel)
  as((β = as(Array, size(model.x, 2)), σ = asℝ₊))
end

# Instantiate the model with data and inits.

x = hcat(ones(size(df, 1)), df[:, :log_mass]);
model = MilkModel(;y=df[:, :kcal_per_g], x=x)

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

function (model::MilkModel)(θ)
    @unpack y, x, = model   # extract the data
    @unpack β, σ = θ            # works on the named tuple too
    ll = 0.0
    ll += logpdf(Normal(0, 100), x[1])
    ll += logpdf(Normal(0, 1), x[2])
    ll += logpdf(TDist(1.0), σ)
    ll += loglikelihood(Normal(0, σ), y .- x*β)
    ll
end

println()
model((β = [1.0, 2.0], σ = 1.0))
println()

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

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P);

# Tune and sample.

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)

println()
DynamicHMC.Diagnostics.EBFMI(results.tree_statistics) |> display
println()

DynamicHMC.Diagnostics.summarize_tree_statistics(results.tree_statistics) |> display
println()

a3d = Array{Float64, 3}(undef, 1000, 3, 1);
for j in 1:1
  for i in 1:1000
    a3d[i, 1:2, j] = values(posterior[i].β)
    a3d[i, 3, j] = values(posterior[i].σ)
  end
end

pnames = ["β[1]", "β[2]", "σ"]
sections =   Dict(
  :parameters => pnames,
)
chns = create_mcmcchains(a3d, pnames, sections, start=1);
chns = set_names(chns, Dict("β[1]" => "α", "β[2]" => "β"))

stan_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
          Mean         SD        Naive SE       MCSE      ESS
    a  0.70472876 0.057040655 0.00090189195 0.0011398893 1000
   bm -0.03150330 0.023642759 0.00037382484 0.0004712342 1000
sigma  0.18378372 0.039212805 0.00062000888 0.0011395979 1000

Quantiles:
          2.5%       25.0%       50.0%        75.0%       97.5%
    a  0.59112968  0.66848775  0.70444950  0.741410500 0.81915225
   bm -0.07729257 -0.04708425 -0.03104865 -0.015942925 0.01424901
sigma  0.12638780  0.15605950  0.17800600  0.204319250 0.27993590
";

describe(chns)

# End of `05/5.6d.jl`


text:
struct CeedError <: Exception
    fname::String
    lineno::Int
    func::String
    ecode::Int
    message::String
end

# COV_EXCL_START
function Base.showerror(io::IO, e::CeedError)
    println(io, "libCEED error code ", e.ecode, " in ", e.func)
    println(io, e.fname, ':', e.lineno, '\n')
    println(io, e.message)
end
# COV_EXCL_STOP

function handle_ceed_error(
    ceed::C.Ceed,
    c_fname::Cstring,
    lineno::Cint,
    c_func::Cstring,
    ecode::Cint,
    c_format::Cstring,
    args::Ptr{Cvoid},
)
    c_message = ccall(
        (:CeedErrorFormat, C.libceed),
        Cstring,
        (C.Ceed, Cstring, Ptr{Cvoid}),
        ceed,
        c_format,
        args,
    )
    fname = unsafe_string(c_fname)
    func = unsafe_string(c_func)
    message = unsafe_string(c_message)
    throw(CeedError(fname, lineno, func, ecode, message))
end

mutable struct Ceed
    ref::RefValue{C.Ceed}
end

"""
    Ceed(spec="/cpu/self")

Wraps a libCEED `Ceed` object, created with the given resource specification string.
"""
function Ceed(spec::AbstractString="/cpu/self")
    obj = Ceed(Ref{C.Ceed}())
    C.CeedInit(spec, obj.ref)
    ehandler = @cfunction(
        handle_ceed_error,
        Cint,
        (C.Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid})
    )
    C.CeedSetErrorHandler(obj.ref[], ehandler)
    finalizer(obj) do x
        # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x))
        destroy(x)
    end
    return obj
end
destroy(c::Ceed) = C.CeedDestroy(c.ref) # COV_EXCL_LINE
Base.getindex(c::Ceed) = c.ref[]

Base.show(io::IO, ::MIME"text/plain", c::Ceed) = ceed_show(io, c, C.CeedView)

"""
    getresource(c::Ceed)

Returns the resource string associated with the given [`Ceed`](@ref) object.
"""
function getresource(c::Ceed)
    res = Ref{Ptr{Cchar}}()
    C.CeedGetResource(c[], res)
    unsafe_string(Cstring(res[]))
end

"""
    isdeterministic(c::Ceed)

Returns true if backend of the given [`Ceed`](@ref) object is deterministic, and false
otherwise.
"""
function isdeterministic(c::Ceed)
    isdet = Ref{Bool}()
    C.CeedIsDeterministic(c[], isdet)
    isdet[]
end

"""
    get_preferred_memtype(c::Ceed)

Returns the preferred [`MemType`](@ref) (either `MEM_HOST` or `MEM_DEVICE`) of the given
[`Ceed`](@ref) object.
"""
function get_preferred_memtype(c::Ceed)
    mtype = Ref{MemType}()
    C.CeedGetPreferredMemType(c[], mtype)
    mtype[]
end

"""
    iscuda(c::Ceed)

Returns true if the given [`Ceed`](@ref) object has resource `"/gpu/cuda/*"` and false
otherwise.
"""
function iscuda(c::Ceed)
    res_split = split(getresource(c), "/")
    length(res_split) >= 3 && res_split[3] == "cuda"
end

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
"""
vtkUnsGridPlot: function for writing vtk unstructured grid files and magnitudes

"""
function vtkUnsGridPlot( connec, nodes, nodalMagnitudes, filename )

    cells = MeshCell{VTKCellType, Vector{Int64}}[]
    for i in (1:size(connec,1) )
      push!(cells, MeshCell( VTKCellTypes.VTK_TETRA, connec[i] ) )
    end

    vtk_grid( filename, Array( nodes' ), cells ) do vtk
        vtk["magnitude", VTKPointData()] = nodalMagnitudes
    end

end

"""
vtkStrGridPlot: function for writing vtk structured grid files and magnitudes
 - nodalMagnitudes: vector of Any with [ magnitudes1, magnitudes2, ...]
"""
function vtkStrGridPlot( grid, nodalMagnitudes, filename )

  x, y, z = LinRange.( grid.startVox, grid.endVox, grid.voxelNums )

  # vtk_grid( filename, x, y, z,  append = false, ascii=true ) do vtk
  vtk_grid( filename, x, y, z ) do vtk
    for i in (1:length(nodalMagnitudes))
      vtk[ string("magnitude_",i), VTKPointData() ] = nodalMagnitudes[i]
    end
  end

end


text:
#=
	we make heavy use of custom type aliases. Have a look at the first lines of quickPIV.jl,
	to see the type abreviations used throughout the file "III", "I", "F"
=#

""" PIV PARAMETER OBJECT """

mutable struct PIVParameters

         corr::CORRTYPES # cross-correlation algorithm
    interSize::III       # interrogation area/volume size in pixels/voxels
 searchMargin::III       # search margin in pixels/voxels
      overlap::III       # overlap in pixels/voxels between adjacent IA/IV
        mpass::I         # multi-pass depth
        width::I         # exclusion radius around 1st peak for finding 2nd peak
         peak::String    # one of "", "gaussian", "centroid"
     sigNoise::String    # one of "", "peak2peak"
      filtFun::Function  # Filtering function to detect IA/IV corresponding to background
    threshold::F         # skip cross correlation if filtFun( interArea ) < threshold
        units::Tuple{DU,DU,DU,TU} # spatial and temporal units of the analyzed data
         args::Tuple{III,III,III,I,I,String,String,Function,F}
end

function parseCorr( corr::String )
    if     lowercase(corr) == "zncc"
		return   ZNCC();
	elseif lowercase(corr) == "fft"
		return    FFT();
	elseif lowercase(corr) == "nsqecc" || lowercase(corr) == "diff"
		return NSQECC();
	else
		return NSQECC();
	end
end
toIII( input::I   ) = ( input, input, input )
toIII( input::II  ) = ( input[1], input[2], 0 )
toIII( input::III ) = input

# convenient construction that takes care of filling up "args"
function PIVParameters( corr::CORRTYPES, interSize::UnI, searchMargin::UnI, ovp::UnI, mpass::I,
                        width::I, peak::String, SN::String, ffn::Function, th, uns::UNTS )

	th_f = th + 0.0; # converting to float
	args = (toIII(interSize), toIII(searchMargin), toIII(ovp), mpass, width, peak, SN, ffn, th_f)

	return PIVParameters( corr, args..., uns, args )
end

# constructor with default values
function setPIVParameters(;corr = "nsqecc", interSize = 32, searchMargin = 0, overlap = 0,
                           mpass = 1, width = 2, peak = "gaussian", sigNoise = "",
                           filtFun = (x)->maxval(x), threshold = -1.0, units = (1m,1m,1m,1second))

    return PIVParameters( parseCorr(corr), interSize, searchMargin, overlap, mpass,
                          width, peak, sigNoise, filtFun, threshold, units )
end

function checkPIVParameters( dsize::NTuple{N,I}, params ) where {N}

    !( N !== 2  ||  N !== 3 ) && error("Your input data isn't 2D nor 3D.")

	any( params.interSize[1:N] .> dsize ) && error("interSize is bigger than the data")
    any( params.interSize .< 2 ) && error("interSize = 1 makes little sense.")
    any( params.overlap .>= params.interSize ) && error("Overlap must be smaller than interSize")
end

""" SYNTHETIC PARTICLES PARAMETER OBJECT """

mutable struct synthParameters
      n::I    # total number of particles, ignored when dens > 0
   dens::I    # particle density ( nºparticles/IA )
      d::I    # depth. If depth > 1, generates a volume
      w::I    # width of generated image/volume
      h::I    # height of generated image/volume
     i0::F    # max intensity of particles
     dt::F    # standard deviation of gaussian around each particles
     th::F    # standard deviation of gaussian around lightsheet
    err::F    # determines rendering radius around each pixel
      z::I    # light-sheet position. Default z = 0
  noise::I    # varible for controlling noise. Not really implemented.
    rad::I    # overwrites the rendering radius computed from err.
   mode::Bool # switches between ( x, y, z ) or ( y, x, z ) coordinates
   args::Tuple{I,I,I,I,I,F,F,F,F,I,I,I,Bool}
end

# convenience constructor
function synthParameters( n::I, dens::I, d::I, w::I, h::I, i0::F, dt::F, th::F, err::F, z::I,
                          noise::I, rad::I, mode::Bool )

	args = (n, dens, d, w, h, i0, dt, th, err, z, noise, rad, mode)
	return synthParameters( args..., args )
end

# constructor with default values
function setSynthParameters(; n=0::I, dens=0::I, d=1::I, w=300::I, h=300::I,
                              z=0::I, i0=255.0::F, dt=3.0::F, err=0.1::F,
                              th=10.0::F, noise=0::I, rad=0::I, mode="whd"::S )

    return synthParameters( n, dens, d, w, h, i0, dt, th, err, z, noise, rad, mode=="whd" )
end

""" TRANSFORMATION PARAMETERS """

mutable struct transformParameters{N,M}
      kind::Symbol      # currently ony :translation is supported
     means::NTuple{N,F} # mean value of transform in x, y [ and z ] if N = 3
      vars::NTuple{N,F} # variance component in covariance matrix
 cov_ratio::NTuple{M,F} # covariance components in covariance matrix. -1 < 1
      args::Tuple{Symbol,NTuple{N,F},NTuple{N,F},NTuple{M,F}}
end

TP2D = transformParameters{2,1}
TP3D = transformParameters{3,3}

# convenience constructor
function transformParameters( kind::Symbol, means::NTuple{N,F}, vars::NTuple{N,F},
                              covs::NTuple{M,F} ) where {N,M}

	args = ( kind, means, vars, covs )
	return transformParameters{N,M}( args..., args )
end

# constructor with default values
function setTransformParameters(; kind=:translation, means=(1.0, 1.0), vars=(eps(Float64),
								  eps(Float64)), cov_ratio=(0.0,) )

    return transformParameters( kind, means, vars, cov_ratio );
end

to2DTransform( t2D::TP2D ) = t2D

function to2DTransform( t3D::TP3D )
    means     = (  t3D.means[1], t3D.means[2] );
    vars      = (  t3D.vars[1] ,  t3D.vars[2] );
    cov_ratio = (  t3D.cov_ratio[1], );

    return transformParameters{2,1}( t3D.kind, means, vars, cov_ratio );
end

function to3DTransform( t2D::TP2D )
    means     = (   t2D.means[1]  ,   t2D.means[2]  , (t2D.means[1] + t2D.means[2])/2 );
    vars      = (    t2D.vars[1]  ,    t2D.vars[2]  , ( t2D.vars[1] +  t2D.vars[2])/2 );
    cov_ratio = ( t2D.cov_ratio[1], t2D.cov_ratio[1], t2D.cov_ratio[1] );

    return transformParameters{3,3}( t2D.kind, means, vars, cov_ratio );
end

to3DTransform( t3D::TP3D ) = t3D

""" EVALUATION PARAMETER OBJECT """

mutable struct metaParameters
 metaloop::I      # determines the number of iterations between "min" and "max" values
 variable::Symbol # specifies what variable will be changed each iteration.
      min::F      # minimum value of changing variable
      max::F      # maximum value of changing variable
  repeats::I      # number of measurements for computing bias and error
     args::Tuple{I,Symbol,F,F,I}
end

function metaParameters( metaloop::I, variable::Symbol, min::F, max::F, repeats::I )

	args = (metaloop, variable, min, max, repeats )
	return metaParameters( args..., args )
end

# short type aliases, to make function signatures shorter
PP   = PIVParameters;
SP   = synthParameters;
TP   = Union{TP2D,TP3D}
MP   = metaParameters;

"""
    Creates a set of parameters for directing PIV evaluation. Accepted parameters are:
        variable (  Symbol ): a symbol indicating what variable should be evaluated. Can be any parameters
                              in synthParameters, PIVParameters or transformParameters.
        metaloop (  Int64  ): number of iterations
            min  ( Float64 ): min value of changing variable
            max  ( Float64 ): max value of changing variable
        repeats  (  Int64  ): number of repetitions
"""
function setMetaParameters( ; metaloop=0::I, var=Symbol("")::Symbol, min=0.0::F,  max=0.0::F, repeats=1::I )

    return metaParameters( metaloop, var, min,  max, repeats )
end

# function to update parameters, and also update the .args tuple, which was created to pass
# parameters by tuple deconstruction "parameter.args..."

function updateParameters!( params, field, newVal )
    hasfield = false;
    fieldidx =     0;
    for e in fieldnames( typeof(params) )
        fieldidx += 1;
        if ( e == field )
            hasfield = true;
            setproperty!( params, e, newVal )
            break
        end
    end

    ( !hasfield ) && return nothing;

    if typeof(params) == PIVParameters
        fieldidx -= 1 # ignore corr parameter
    end
    newargs = [ arg for arg in params.args ]
    newargs[ fieldidx ] = newVal;
    setproperty!( params, :args, Tuple(newargs) )
end

max_issues_repo_name: ErichZimmer/quickPIV

max_forks_repo_name: ErichZimmer/quickPIV


text:
function try_decode_jwt(contents::AbstractString)
    try
        jwt_decoded = JSONWebTokens.decode(JSONWebTokens.None(), contents)
        return true, jwt_decoded
    catch
    end
    return false, Dict{String, Any}()
end


text:
include("SortSparse.jl")
include("KoLesky.jl")
include("CovFuncs.jl")
using LinearAlgebra
using Distributions
using Distances
using Random
using Profile
using JLD
using MAT
using Match
Random.seed!(123)

NTrain = 2000
NTest = 100
N = NTest + NTrain
d =  2
y = randn(N)
x = randn(d, N)

setSplit = randperm(N)
trainInds = setSplit[1:NTrain]
testInds = setSplit[(NTrain + 1) : end]

xTrain = x[ :, trainInds ]
yTrain = y[ trainInds ]
xTest = x[:, testInds ]
yTest = y[ testInds ]

ρ = 2.5
λ = 1.5

l = 0.5

maternOrder = "12"
covfunc12!(r) = matern12!(r,l)
covfunc32!(r) = matern32!(r,l)
covfunc52!(r) = matern52!(r,l)

# Make sure to avoid NANs
@match maternOrder begin
  "12" => (covfunc! = covfunc12!)
  "32" => (covfunc! = covfunc32!)
  "52" => (covfunc! = covfunc52!)
end

@time colptr, rowval, P, revP, distances = sortSparseRev(xTrain,ρ, 1)
xTrain = xTrain[:, P]
yTrain = yTrain[P]
x = hcat(xTrain, xTest)

skeletons = construct_skeletons(colptr, rowval, distances, λ)
# TODO: Remove again
# skeletons = singleParents(skeletons)

colptr = 0
length(skeletons)
rowval = 0
GC.gc()

#Comparison to full model
LinearAlgebra.BLAS.set_num_threads(8)
@time begin
  Θ = covfunc!(pairwise(Euclidean(), x, dims=2))
  μFull = Θ[(NTrain + 1) : end, 1 : NTrain] * ( Θ[1 : NTrain, 1 : NTrain] \ yTrain)
end

Θ = 0
GC.gc()
# @show mean(abs.(yTest - μFull))
# @show mean(abs.(yTest - μFull)) / std(yTest)

@time μ, σ = predict(skeletons, (NTrain+1): N, yTrain, x, covfunc!)

# @show mean(abs.(yTest - μ)) / std(yTest)
# @show mean(abs.(yTest - μ))

@show norm(μFull - μ) / norm(μFull)
@show norm(μFull - μ)


text:
import Nemo

function read_SEAIJRC()
    @info "Reading SEAIJRC"
    polys = []
    apath = "C:\\data\\projects\\gbdata\\gbSEAIJRC"
    open(apath, "r") do f
        head = "b, alpha, g2, Ninv, k, g1, q, r, sat_aux1, sat_aux2"
        body = readline(f)
        @info "First line" head
        headv = map(x -> Symbol(x), split(head, ", "))
        e = Meta.parse(join(headv, ","))
        R, xs = eval(:((R, $e) = Nemo.PolynomialRing(Nemo.QQ, $headv)))
        @info "Created ring" R xs
        rrr = split(body, ",")
        for (i, s) in enumerate(rrr)
            @info "Reading poly $i/$(length(rrr))"
            poly = eval(Meta.parse(s))
            push!(polys, poly)
        end
    end
    @info "Read" length(polys)
    R.(polys)
end

function read_SIWR()
    @info "Reading SIWR"
    polys = []
    apath = "C:\\data\\projects\\gbdata\\gbSIWR"
    open(apath, "r") do f
        head = readline(f)
        body = readline(f)
        headv = map(x -> Symbol(x), split(head, ", "))
        e = Meta.parse(join(headv, ","))
        R, xs = eval(:((R, $e) = Nemo.PolynomialRing(Nemo.QQ, $headv)))
        @info "Created ring" R xs
        rrr = split(body, ",")
        for (i, s) in enumerate(rrr)
            @info "Reading poly $i/$(length(rrr))"
            poly = eval(Meta.parse(s))
            push!(polys, poly)
        end
    end
    @info "Read" length(polys)
    R.(polys)
end

function read_MAPK()
    @info "Reading MAPK"
    polys = []
    apath = "../../MAPK"
    # apath = "MAPK"
    open(apath, "r") do f
        head = "c0001, a10, gamma1000, alpha10, b00, beta11, c0111, beta10, alpha11, beta01, alpha01, gamma1100, c0011, c0010, b10, gamma1101, a00, b01, c1011, a01, gamma1110, gamma0100, sat_aux1, sat_aux2, sat_aux3, sat_aux4, sat_aux5"
        body = readline(f)
        @info "First line" head
        headv = map(x -> Symbol(x), split(head, ", "))
        e = Meta.parse(join(headv, ","))
        R, xs = eval(:((R, $e) = Nemo.PolynomialRing(Nemo.QQ, $headv)))
        @info "Created ring" R xs
        rrr = split(body, ",")
        l = length(rrr)
        for i in 1:length(rrr)
            @info "Reading poly $i/$(length(rrr))"
            str = map(strip, split(rrr[1], "+"))
            poly = R(0)
            for s1 in str
                # println(s1)
                # println("_________________")
                if occursin("-", s1)
                    sp1 = split(s1, "-")
                    if startswith(s1, "-")
                        sp1 = sp1[2:end]
                        poly -= eval(Meta.parse(sp1[1]))
                    else
                        poly += eval(Meta.parse(sp1[1]))
                    end
                    for s2 in sp1[2:end]
                        poly -= eval(Meta.parse(s2))
                    end
                else
                    poly += eval(Meta.parse(s1))
                end
            end
            push!(polys, poly)
            popfirst!(rrr)
        end
    end
    @info "Read" length(polys)
    R.(polys)
end

function read_MAPK_include()
    @info "Reading MAPK"
    polys = []
    apath = "../../MAPK"

    head = "c0001, a10, gamma1000, alpha10, b00, beta11, c0111, beta10, alpha11, beta01, alpha01, gamma1100, c0011, c0010, b10, gamma1101, a00, b01, c1011, a01, gamma1110, gamma0100, sat_aux1, sat_aux2, sat_aux3, sat_aux4, sat_aux5"
    @info "First line" head
    headv = map(x -> Symbol(x), split(head, ", "))
    e = Meta.parse(join(headv, ","))
    R, xs = eval(:((R, $e) = Nemo.PolynomialRing(Nemo.QQ, $headv)))

    collect(include(apath))
end

function parse_include(filename)
    head = "x, y, z"
    headv = map(x -> Symbol(x), split(head, ", "))
    e = Meta.parse("x, y, z")
    R, xs = eval(:((R, $e) = Nemo.PolynomialRing(Nemo.QQ, $headv)))

    include(abspath("benchmark\\data\\$filename"))
end


text:
macro yao_code_lowered(ex)
    esc(yao_code_lowered_m(ex))
end

macro yao_code_lowered(ctrl, ex)
    esc(yao_code_lowered_m(ctrl, ex))
end

macro yao_code_typed(ex)
    esc(yao_code_typed_m(ex))
end

macro yao_code_typed(ctrl, ex)
    esc(yao_code_typed_m(ctrl, ex))
end

function yao_code_lowered_m(ex)
    @match ex begin
        # apply(reg, gate, locs)
        :($locs => $gate) => quote
            Base.code_lowered(
                $Intrinsics.apply,
                ($YaoCompiler.AnyReg, typeof($gate), typeof($Locations($(locs)))),
            )
        end
        # apply(reg, gate)
        _ => quote
            Base.code_lowered($Intrinsics.apply, ($YaoCompiler.AnyReg, typeof($ex)))
        end
    end
end

function yao_code_lowered_m(ctrl, ex)
    @match ex begin
        # apply(reg, gate, locs, ctrl)
        :($locs => $gate) => quote
            Base.code_lowered(
                $Intrinsics.apply,
                (
                    $YaoCompiler.AnyReg,
                    typeof($gate),
                    typeof($Locations($(locs))),
                    typeof($CtrlLocations($ctrl)),
                ),
            )
        end
        _ => error("expect a locs => gate expression")
    end
end

function yao_code_typed_m(ex)
    @match ex begin
        # apply(reg, gate, locs)
        :($locs => $gate) => quote
            $Base.code_typed(
                $Intrinsics.apply,
                ($YaoCompiler.AnyReg, typeof($gate), typeof($Locations($(locs))));
                interp = $(YaoInterpreter()),
            )
        end
        # apply(reg, gate)
        _ => quote
            $Base.code_typed(
                $Intrinsics.apply,
                ($YaoCompiler.AnyReg, typeof($ex));
                interp = $(YaoInterpreter()),
            )
        end
    end
end

function yao_code_typed_m(ctrl, ex)
    @match ex begin
        # apply(reg, gate, locs, ctrl)
        :($locs => $gate) => quote
            Base.code_typed(
                $Intrinsics.apply,
                (
                    $YaoCompiler.AnyReg,
                    typeof($gate),
                    typeof($Locations($(locs))),
                    typeof($CtrlLocations($ctrl)),
                );
                interp = $(YaoInterpreter()),
            )
        end
        _ => error("expect a locs => gate expression")
    end
end


text:
# This is from https://github.com/JuliaGPU/CUDAapi.jl/pull/84/files

if has_cuda()
    try
        using CuArrays: CuArrays, CuPtr, unsafe_free!, usage_limit
        import CuArrays: CuArray
    catch ex
        @warn "CUDA is installed, but CuArrays.jl fails to load" exception=(ex,catch_backtrace())
    end
end

### Use CuArrays kernels as fallback for undefined KnetArray operations.

import Base: getindex, setindex!, permutedims, permutedims!, cat, hcat, vcat

# Extend function CuArray to create a memory shared CuArray from KnetArray:
# Avoid the cu function as it changes eltype to Float32
function CuArray(x::KnetArray{T}) where {T}
    p = CuPtr{T}(UInt(x.ptr.ptr))
    Base.unsafe_wrap(CuArray{T}, p, size(x); own=false)
end

# Based on _unsafe_getindex, multidimensional.jl:679, julia 1.2.0
function getindex(A::KnetArray, I...)
    I = Base.to_indices(A, I)
    shape = Base.index_shape(I...)
    B = similar(A, length.(shape))
    Base._unsafe_getindex!(CuArray(B), CuArray(A), I...)
    return B
end

function setindex!(A::KnetArray, B, I...)
    if B isa KnetArray || B isa AbstractArray
        B = CuArray(B)
    end
    setindex!(CuArray(A), B, I...)
    return A
end

permutedims!(y::KnetArray, x::KnetArray, perm) = (permutedims!(CuArray(y), CuArray(x), perm); y)

# Based on permutedims, multidimensional.jl:1334, julia 1.2.0
function permutedims(B::KnetArray,perm)
    dimsB = size(B)
    ndimsB = length(dimsB)
    (ndimsB == length(perm) && isperm(perm)) || throw(ArgumentError("no valid permutation of dimensions"))
    dimsP = ntuple(i->dimsB[perm[i]], ndimsB)::typeof(dimsB)
    P = similar(B, dimsP)
    permutedims!(P,B,perm)
end

#permutedims(x::KnetMatrix)=permutedims(x,(2,1))  # CuArrays is %10 faster but has startup cost
permutedims(x::KnetMatrix)=_transpose(x)          # cuDNN is %10 slower but no startup cost
permutedims(x::KnetVector)=copy(reshape(x,1,:))

using Base: dims2cat, cat_shape, __cat

# vcat(X::KnetArray...)=cat(X...; dims=Val(1)) # karray.jl version is 30%-80% faster
hcat(X::KnetArray...)=cat(X...; dims=Val(2))   # This should only kick in for dims > 2, karray.jl 1/2D versions are 100% faster

# Based on _cat_t, abstractarray.jl:1439, julia 1.2.0
function cat(X::KnetArray{T}...; dims) where {T}
    catdims = dims2cat(dims)
    # catdims == (true,) && return vcat_old(X...) # 10-30% faster
    shape = cat_shape(catdims, (), map(size, X)...)
    # length(shape) <= 2 && catdims == (false,true) && return hcat_old(X...) # 50% faster
    A = similar(X[1], T, shape) # cat_similar(X[1], T, shape)
    if T <: Number && count(!iszero, catdims) > 1
        fill!(A, zero(T))
    end
    __cat(CuArray(A), shape, catdims, map(CuArray,X)...)
    return A
end

# Must be careful with memory management, for now we will let Knet manage memory.
# use CuArray(x) with overwriting kernels only.
# use the following with caution.

if has_cuda()
    function KnetArray(x::CuArray{T,N}) where {T,N}
        p = Base.bitcast(Knet.Cptr, x.ptr)
        k = Knet.KnetPtr(p, sizeof(x), gpu(), x)
        KnetArray{T,N}(k, size(x))
    end
end

# Testing the CuArrays allocator: set Knet.cuallocator()=true to use this
function KnetPtrCu(len::Int)
    c = CuArray{UInt8}(undef, len)
    p = convert(Cptr, convert(UInt, Base.unsafe_convert(CuPtr{UInt8}, Base.cconvert(CuPtr{UInt8}, c))))
    kp = KnetPtr(p, len, gpu(), c)
    finalizer(freeKnetPtrCu, kp)
end

function freeKnetPtrCu(p::KnetPtr)
    # GC.gc comes here directly, manual calls come through freeKnetPtr()
    # After a manual call, GC.gc may call the finalizer again, avoid double free
    if p.parent isa CuArray
        unsafe_free!(p.parent)
        p.ptr, p.parent = C_NULL, nothing
    elseif p.parent isa Nothing
        @assert p.ptr == C_NULL
        # already freed, do nothing
    elseif p.parent isa KnetPtr
        # subarray, do nothing
    else
        error("Bad parent pointer $(typeof(p.parent))")
    end
end

# argmax, argmin etc. Fixes https://github.com/denizyuret/Knet.jl/issues/368.
# Two options: argmax(Array(KnetArray)) vs argmax(CuArray)
# Experiments: 10x10, 100x100, 1000x1000 with no dims, dims=1, dims=2
# With no dims, CuArrays is better for 100x100, 1000x1000.
# With all others, KnetArray is better.

import Base: argmax, argmin, findmax, findmin
# TODO: try this again after https://github.com/JuliaGPU/CuArrays.jl/issues/304 is resolved
# argmaxarray(x,d)=((d===:) && length(x) > 4096 ? CuArray(x) : Array(x))
argmaxarray(x,d)=Array(x)
argmax(x::KnetArray; dims=:)=argmax(argmaxarray(x,dims); dims=dims)
argmin(x::KnetArray; dims=:)=argmin(argmaxarray(x,dims); dims=dims)
findmax(x::KnetArray; dims=:)=findmax(argmaxarray(x,dims); dims=dims)
findmin(x::KnetArray; dims=:)=findmin(argmaxarray(x,dims); dims=dims)

# Issue #108:Element-wise power of KnetArray give NaN results #108
# This is a bug with CUDA giving NaN for integer powers of negative numbers (powf is broken)

import Base.Broadcast: broadcasted

function broadcasted(::typeof(^),a::KnetArray{T},s::Number) where T
    b = similar(a)
    ca = CuArray(a)
    cb = CuArray(b)
    cb .= ca .^ T(s)
    return b
end

function broadcasted(::typeof(^),s::Number,a::KnetArray{T}) where T
    b = similar(a)
    ca = CuArray(a)
    cb = CuArray(b)
    cb .= T(s) .^ ca
    return b
end


text:
module actor_world_test
using Test
include("../src/ActorsWorld.jl"); using .ActorsWorld

@testset "Actors-World" begin
    @testset "Basic" begin
        or = Orientation(ActorsWorld.DIRECTIONS[1])
        @test orientation_rotate(or,Val{true}).value == ActorsWorld.DIRECTIONS[2]
        @test orientation_rotate(or,Val{false}).value == ActorsWorld.DIRECTIONS[4]
        @test orientation_rotate(Orientation(ActorsWorld.DIRECTIONS[2]),Val{true}).value == ActorsWorld.DIRECTIONS[3]
        @test orientation_rotate(Orientation(ActorsWorld.DIRECTIONS[3]),Val{true}).value == ActorsWorld.DIRECTIONS[4]
        @test orientation_rotate(Orientation(ActorsWorld.DIRECTIONS[4]),Val{true}).value == ActorsWorld.DIRECTIONS[1]
        @test_throws InvalidDirectionError Orientation(:foo)
        lo = Location(1,2)
        wo = World(10,13)
        @test location_move(lo,Orientation(ActorsWorld.DIRECTIONS[2])).x == 2
        @test location_fix_ooBound(wo,Location(0,15)) == Location(10,2)
    end
    @testset "World generation" begin
        wo = World(10,12)
        @test wo.size.height == 12
        @test wo.size.width == 10
    end
    @testset "Actors" begin
        kara = Actor_Definition(
            moveable=true,
            turnable=true
        )
        leaf = Actor_Definition(
            layer=0
        )
        wo = World(10,10)
        ac1 = actor_create!(wo,kara,Location(1,3),Orientation(ActorsWorld.DIRECTIONS[1]))
        ac1_c = copy(ac1)
        @test ac1_c.actor_definition == ac1.actor_definition
        @test ac1_c.location == ac1.location
        @test ac1_c.orientation == ac1.orientation
        @test ac1 == wo.actors[1]
        actor_delete!(wo,ac1)
        @test length(wo.actors) == 0
        ac2 = actor_create!(wo,kara,Location(1,3),Orientation(ActorsWorld.DIRECTIONS[1]))
        ac3 = actor_create!(wo,leaf,Location(1,3),Orientation(ActorsWorld.DIRECTIONS[1]))
        @test length(get_actors_at_location(wo,Location(1,3))) == 2
        @test length(get_actors_at_location(wo,Location(2,3))) == 0
        @test_throws LocationFullError actor_create!(wo,kara,Location(1,3),
                                                     Orientation(ActorsWorld.DIRECTIONS[1]))
        k2 = actor_create!(wo,kara,Location(5,5),Orientation(ActorsWorld.DIRECTIONS[1]))
        @test_throws LocationFullError actor_create!(wo,kara,Location(5,5),Orientation(ActorsWorld.DIRECTIONS[1]))
        actor_delete!(wo,k2)
        @test_throws ActorNotFound actor_delete!(wo,k2)

        actor_moveto!(wo,ac2,Location(10,10))
        @test ac2.location == Location(10,10)
        actor_create!(wo,kara,Location(5,5),Orientation(ActorsWorld.DIRECTIONS[1]))
        @test_throws LocationFullError actor_moveto!(wo,ac2,Location(5,5))
        wo_c = copy(wo)
        @test length(wo_c.actors) == length(wo.actors)
        @test wo_c.size == wo.size
    end
    @testset "World Boundaries" begin
        ka = Actor_Definition(
            moveable=true,
            turnable=true
        )
        wo = World(2,5)
        ac = actor_create!(
            wo,ka,Location(1,1),Orientation(ActorsWorld.DIRECTIONS[1])
        )
        @test_throws LocationOutsideError actor_create!(
            wo,ka,Location(1,10),Orientation(ActorsWorld.DIRECTIONS[1])
        )
        actor_move!(wo,ac,ActorsWorld.DIRECTIONS[1]) # 1,2
        actor_move!(wo,ac,ActorsWorld.DIRECTIONS[1]) # 1,3
        actor_move!(wo,ac,ActorsWorld.DIRECTIONS[1]) # 1,4
        actor_move!(wo,ac,ActorsWorld.DIRECTIONS[1]) # 1,5
        @test ac.location == Location(1,5)
        actor_move!(wo,ac,ActorsWorld.DIRECTIONS[1]) # 1,1
        @test ac.location == Location(1,1)
        actor_move!(wo,ac,ActorsWorld.DIRECTIONS[1]) # 1,2
        @test ac.location == Location(1,2)
        actor_rotate!(ac,true)
        @test ac.orientation == Orientation(ActorsWorld.DIRECTIONS[2])
        actor_move!(wo,ac,ActorsWorld.DIRECTIONS[2]) # 2,2
        actor_move!(wo,ac,ActorsWorld.DIRECTIONS[2]) # 1,2
        @test ac.location == Location(1,2)
    end
    @testset "Moving other Actors" begin
        kara_2 = Actor_Definition(
            moveable=true,
            turnable=true
        )
        mushroom = Actor_Definition(
            moveable=true,
        )
        tree = Actor_Definition(
            moveable=false
        )
        wo = World(2,5)
        ac_kara_2 = actor_create!(
            wo,kara_2,Location(1,1),Orientation(ActorsWorld.DIRECTIONS[1])
        )
        ac_mushroom = actor_create!(
            wo,mushroom,Location(1,2),Orientation(ActorsWorld.DIRECTIONS[1])
        )
        actor_move!(wo,ac_kara_2,ActorsWorld.DIRECTIONS[1])
        @test ac_mushroom.location.y == 3
        actor_move!(wo,ac_kara_2,ActorsWorld.DIRECTIONS[3])
        ac_mushroom2 = actor_create!(
            wo,mushroom,Location(1,2),Orientation(ActorsWorld.DIRECTIONS[1])
        )
        actor_create!(
            wo,tree,Location(2,1),Orientation(ActorsWorld.DIRECTIONS[1])
        )
        @test_throws ActorInvalidMultipleMovementError actor_move!(wo,ac_kara_2,ActorsWorld.DIRECTIONS[1])
        @test_throws ActorInvalidMovementError actor_move!(wo,ac_kara_2,ActorsWorld.DIRECTIONS[2])
    end
    @testset "Putting and Picking" begin
        wo = World(1,10)
        leaf = Actor_Definition(
            layer=0,
            grabable=true
        )
        leaf_2 = Actor_Definition(
            layer=-1,
            grabable=true
        )
        kara_p = Actor_Definition()
        kara_ac = actor_create!(wo,kara_p,Location(1,1),Orientation(ActorsWorld.DIRECTIONS[1]))
        actor_putdown!(wo,kara_ac,leaf)
        @test wo.actors[2].actor_definition == leaf
        actor_pickup!(wo,kara_ac)
        @test length(wo.actors) == 1
        @test_throws ActorGrabNotFoundError actor_pickup!(wo,kara_ac)
        @test_throws ActorInvalidGrabError actor_putdown!(wo,kara_ac,kara_p)
        @test_throws ActorNotPlaceableError actor_putdown!(wo,kara_ac,leaf_2)
    end
    @testset "Sensors" begin
        wo = World(10,10)
        acd = Actor_Definition()
        acdp = Actor_Definition(layer=0)
        ac = actor_create!(wo,acd,Location(5,5),Orientation(:NORTH))
        actor_create!(wo,acd,Location(6,5),Orientation(:NORTH))
        actor_create!(wo,acd,Location(5,6),Orientation(:NORTH))
        actor_create!(wo,acdp,Location(5,5),Orientation(:NORTH))
        @test is_actor_definition_left(wo,ac,acd) == false
        @test is_actor_definition_right(wo,ac,acd) == true
        @test is_actor_definition_front(wo,ac,acd) == true
        @test is_actor_definition_right(wo,ac,Actor_Definition(moveable=true)) == false
        @test is_actor_definition_here(wo,ac,acdp) == true
    end
end
end


text:
using StaticArrays
using IntervalArithmetic
import IntervalArithmetic: interval

using LinearAlgebra

import Base: +, *, ^, -, sqrt, inv

"""
Affine form with center `c`, affine components `γ` and error `Δ`.

Variant where Δ is an interval
"""
struct AFF{N,T<:AbstractFloat}
    affine::AF{N,T}
    range::Interval{T}
end

interval(x::AFF) = x.range

function Base.show(io::IO, C::AFF{N,T}) where {N,T}
    print(io, "affine=", C.affine, "; range=", C.range)
end

# ==(C::Affine, D::Affine) = C.c == D.c && C.γ == D.γ

"""
Make an `AF` based on an interval, which is number `i` of `n` total variables.
"""
function AFF(X::Interval{T}, n, i) where {T}
    return AFF(AF(X, n, i), X)
end

for op in (:+, :*, :-)
    @eval function $op(x::AFF, y::AFF)
        affine = $op(x.affine, y.affine)

        range = $op(x.range, y.range)

        range = range ∩ interval(affine)

        return AFF(affine, range)
    end
end

for op in (:sqrt, :inv)
    @eval function $op(x::AFF)
        affine = $op(x.affine, x.range)

        range = $op(x.range)

        range = range ∩ interval(affine)

        return AFF(affine, range)
    end
end

*(x::AFF, α::Real) = AFF(α*x.affine, α*x.range)
*(α::Real, x::AFF) = x * α

+(x::AFF, α::Real) = AFF(α+x.affine, α + x.range)
+(α::Real, x::AFF) = x + α

-(x::AFF) = AFF(-x.affine, -x.range)
-(x::AFF, α::Real) = AFF(x.affine - α, x.range - α)
-(α::Real, x::AFF) = α + (-x)

/(x::AFF, α::Real) = AFF(x.affine/α, x.range/α)

function ^(x::AFF, n::Integer)

    invert = false

    if n < 0
        invert = true
        n = -n
        # @show n
    end

    result = Base.power_by_squaring(x, n)

    if invert
        result = inv(result)
    end

    return result
end
Base.literal_pow(::typeof(^), x::AFF, ::Val{p}) where {T,p} = x^p

# x = AF{2,Float64}(0.0, SVector(1.0, 0.0), 0..0)
# y = AF{2,Float64}(0.0, SVector(0.0, 1.0), 0..0)

#
# x = AF(3..5, 2, 1)
# y = AF(2..4, 2, 2)
#
# 3-x
# interval(3-x)
#
# x + y
#
#
# interval(x+y)
#
# x * y
# interval(x * y)
#
# interval(x * y)
# interval(x) * interval(y)
#
# z = AF(-1..1, 1, 1)
# z^2
# interval(z^2)
#
# using Polynomials
#
# p = Poly([-3, 1])
# p2 = p^8
#
# x = 4 ± 1e-4
# y = AF(x, 1, 1)
#
# interval(y)
# interval(p2(x))
# interval(p2(y))
#
# @time interval(p2(y))
#
#
# f( (x, y) ) = [x^3 + y, (x - y)^2]
#
# X = IntervalBox(-1..1, -1..1)
#
# f(X)
#
# xx = AF(X[1], 2, 1)
# yy = AF(X[2], 2, 2)
#
# interval.(f((xx, yy)))
#
# f(X)
#
#
#
#
# x = AF(4..6, 1, 1)    # example from Messine
# f(x) = x * (10 - x)
#
# f(x)
# interval(f(x))
#
# interval(10*x - x^2)