text:
# ===============================
# Written by AAE
# TU Delft, Spring 2014
# simulkade.com
# ===============================

# ================================================================
# Changes:
#    2014-12-30 added 2D radial and 3D cylindrical grids
#    2014-12-31 many many bugs fixed
#    2015-01-10 extended to accept nonuniform grids
# ================================================================

function divergenceTerm(F::FaceValue)
d = F.domain.dimension
if d==1
  RHSdiv = divergenceTerm1D(F)
elseif d==1.5
  RHSdiv = divergenceTermCylindrical1D(F)
elseif d==2
  RHSdiv, RHSdivx, RHSdivy = divergenceTerm2D(F)
elseif d==2.5
  RHSdiv, RHSdivx, RHSdivy = divergenceTermCylindrical2D(F)
elseif d==2.8
  RHSdiv, RHSdivx, RHSdivy = divergenceTermRadial2D(F)
elseif d==3
  RHSdiv, RHSdivx, RHSdivy, RHSdivz = divergenceTerm3D(F)
elseif d==3.2
  RHSdiv, RHSdivx, RHSdivy, RHSdivz = divergenceTermCylindrical3D(F)
end
RHSdiv

end

# =============== Divergence 1D Term ============================
function divergenceTerm1D(F::FaceValue)
# This function calculates the divergence of a field
# using its face

# extract data from the mesh structure
Nx = F.domain.dims[1]
G = [1:Nx+2;]
DX = F.domain.cellsize.x[2:end-1]

# define the vector of cell index
row_index = reshape(G[2:Nx+1],Nx) # main diagonal

# compute the divergence
div_x = (F.xvalue[2:Nx+1]-F.xvalue[1:Nx])./DX

# define the RHS Vector
RHSdiv = zeros(Nx+2)

# assign the values of the RHS vector
RHSdiv[row_index] = reshape(div_x,Nx)

return RHSdiv

end

# =============== Divergence Cylindrical 1D Term ============================
function divergenceTermCylindrical1D(F::FaceValue)
# This function calculates the divergence of a field
# using its face

# extract data from the mesh structure
Nx = F.domain.dims[1]
G = [1:Nx+2;]
DX = F.domain.cellsize.x[2:end-1]
rp = F.domain.cellcenters.x
rf = F.domain.facecenters.x

# define the vector of cell index
row_index = reshape(G[2:Nx+1],Nx) # main diagonal

# reassign the east, west, north, and south flux vectors for the
# code readability
Fe = F.xvalue[2:Nx+1]
Fw = F.xvalue[1:Nx]
re = rf[2:Nx+1]
rw = rf[1:Nx]

# compute the divergence
div_x = (re.*Fe-rw.*Fw)./(DX.*rp)

# define the RHS Vector
RHSdiv = zeros(Nx+2)

# assign the values of the RHS vector
RHSdiv[row_index] = reshape(div_x,Nx)

return RHSdiv

end

# =============== Divergence 2D Term ============================
function divergenceTerm2D(F::FaceValue)
# This function calculates the divergence of a field
# using its face

# extract data from the mesh structure
Nx = F.domain.dims[1]
Ny = F.domain.dims[2]
G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)
DX = F.domain.cellsize.x[2:end-1]
DY = zeros( 1, Ny)
DY[:] = F.domain.cellsize.y[2:end-1]

# define the vector of cell index
row_index = reshape(G[2:Nx+1,2:Ny+1],Nx*Ny) # main diagonal

# reassign the east, west, north, and south flux vectors for the
# code readability
Fe = F.xvalue[2:Nx+1,:]
Fw = F.xvalue[1:Nx,:]
Fn = F.yvalue[:,2:Ny+1]
Fs = F.yvalue[:,1:Ny]

# compute the divergence
div_x = (Fe - Fw)./DX
div_y = (Fn - Fs)./DY

# define the RHS Vector
RHSdiv = zeros((Nx+2)*(Ny+2))
RHSdivx = zeros((Nx+2)*(Ny+2))
RHSdivy = zeros((Nx+2)*(Ny+2))

# assign the values of the RHS vector
RHSdiv[row_index] = reshape(div_x+div_y,Nx*Ny)
RHSdivx[row_index] = reshape(div_x,Nx*Ny)
RHSdivy[row_index] = reshape(div_y,Nx*Ny)

(RHSdiv, RHSdivx, RHSdivy)

end

# =============== Divergence 2D Cylindrical Term ============================
function divergenceTermCylindrical2D(F::FaceValue)
# This function calculates the divergence of a field
# using its face

# extract data from the mesh structure
Nr = F.domain.dims[1]
Nz = F.domain.dims[2]
G=reshape([1:(Nr+2)*(Nz+2);], Nr+2, Nz+2)
dr = F.domain.cellsize.x[2:end-1]
dz= zeros( 1, Nz)
dz[:] = F.domain.cellsize.y[2:end-1]
rp = F.domain.cellcenters.x
rf = F.domain.facecenters.x

# define the vector of cell index
row_index = reshape(G[2:Nr+1,2:Nz+1],Nr*Nz) # main diagonal

# reassign the east, west, north, and south flux vectors for the
# code readability
Fe = F.xvalue[2:Nr+1,:]
Fw = F.xvalue[1:Nr,:]
Fn = F.yvalue[:,2:Nz+1]
Fs = F.yvalue[:,1:Nz]
re = rf[2:Nr+1]
rw = rf[1:Nr]

# compute the divergence
div_x = (re.*Fe - rw.*Fw)./(dr.*rp)
div_y = (Fn - Fs)./dz

# define the RHS Vector
RHSdiv = zeros((Nr+2)*(Nz+2))
RHSdivx = zeros((Nr+2)*(Nz+2))
RHSdivy = zeros((Nr+2)*(Nz+2))

# assign the values of the RHS vector
RHSdiv[row_index] = reshape(div_x+div_y,Nr*Nz)
RHSdivx[row_index] = reshape(div_x,Nr*Nz)
RHSdivy[row_index] = reshape(div_y,Nr*Nz)

(RHSdiv, RHSdivx, RHSdivy)

end

# =============== Divergence 2D Radial Term ============================
function divergenceTermRadial2D(F::FaceValue)
# This function calculates the divergence of a field
# using its face

# extract data from the mesh structure
Nr = F.domain.dims[1]
Ntheta = F.domain.dims[2]
G=reshape([1:(Nr+2)*(Ntheta+2);], Nr+2, Ntheta+2)
dr = F.domain.cellsize.x[2:end-1]
dtheta= zeros( 1, Ntheta)
dtheta[:]= F.domain.cellsize.y[2:end-1]
rp = F.domain.cellcenters.x
rf = F.domain.facecenters.x

# define the vector of cell index
row_index = reshape(G[2:Nr+1,2:Ntheta+1],Nr*Ntheta) # main diagonal

# reassign the east, west, north, and south flux vectors for the
# code readability
Fe = F.xvalue[2:Nr+1,:]
Fw = F.xvalue[1:Nr,:]
Fn = F.yvalue[:,2:Ntheta+1]
Fs = F.yvalue[:,1:Ntheta]
re = rf[2:Nr+1]
rw = rf[1:Nr]

# compute the divergence
div_x = (re.*Fe-rw.*Fw)./(dr.*rp)
div_y = (Fn-Fs)./(dtheta.*rp)

# define the RHS Vector
RHSdiv = zeros((Nr+2)*(Ntheta+2))
RHSdivx = zeros((Nr+2)*(Ntheta+2))
RHSdivy = zeros((Nr+2)*(Ntheta+2))

# assign the values of the RHS vector
RHSdiv[row_index] = reshape(div_x+div_y,Nr*Ntheta)
RHSdivx[row_index] = reshape(div_x,Nr*Ntheta)
RHSdivy[row_index] = reshape(div_y,Nr*Ntheta)

(RHSdiv, RHSdivx, RHSdivy)

end

# =============== Divergence 3D Term ============================
function divergenceTerm3D(F::FaceValue)
# This function calculates the divergence of a field
# using its face

# extract data from the mesh structure
Nx = F.domain.dims[1]
Ny = F.domain.dims[2]
Nz = F.domain.dims[3]
G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)
dx = F.domain.cellsize.x[2:end-1]
dy = zeros( 1, Ny)
dy[:] = F.domain.cellsize.y[2:end-1]
dz = zeros( 1,1,Nz)
dz[:] = F.domain.cellsize.z[2:end-1]

# define the vector of cell index
row_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],Nx*Ny*Nz) # main diagonal

# reassign the east, west, north, and south flux vectors for the
# code readability
Fe = F.xvalue[2:Nx+1,:,:]
Fw = F.xvalue[1:Nx,:,:]
Fn = F.yvalue[:,2:Ny+1,:]
Fs = F.yvalue[:,1:Ny,:]
Ff = F.zvalue[:,:,2:Nz+1]
Fb = F.zvalue[:,:,1:Nz]

# compute the divergence
div_x = (Fe - Fw)./dx
div_y = (Fn - Fs)./dy
div_z = (Ff - Fb)./dz

# define the RHS Vector
RHSdiv = zeros((Nx+2)*(Ny+2)*(Nz+2))
RHSdivx = zeros((Nx+2)*(Ny+2)*(Nz+2))
RHSdivy = zeros((Nx+2)*(Ny+2)*(Nz+2))
RHSdivz = zeros((Nx+2)*(Ny+2)*(Nz+2))

# assign the values of the RHS vector
RHSdiv[row_index] = reshape(div_x+div_y+div_z,Nx*Ny*Nz)
RHSdivx[row_index] = reshape(div_x,Nx*Ny*Nz)
RHSdivy[row_index] = reshape(div_y,Nx*Ny*Nz)
RHSdivz[row_index] = reshape(div_z,Nx*Ny*Nz)

(RHSdiv, RHSdivx, RHSdivy, RHSdivz)

end

# =============== Divergence 3D Cylindrical Term ============================
function divergenceTermCylindrical3D(F::FaceValue)
# This function calculates the divergence of a field
# using its face

# extract data from the mesh structure
Nx = F.domain.dims[1]
Ny = F.domain.dims[2]
Nz = F.domain.dims[3]
G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)
dx = F.domain.cellsize.x[2:end-1]
dy = zeros( 1, Ny)
dy[:] = F.domain.cellsize.y[2:end-1]
dz = zeros( 1,1,Nz)
dz[:] = F.domain.cellsize.z[2:end-1]
rp = F.domain.cellcenters.x

# define the vector of cell index
row_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],Nx*Ny*Nz) # main diagonal

# reassign the east, west, north, and south flux vectors for the
# code readability
Fe = F.xvalue[2:Nx+1,:,:]
Fw = F.xvalue[1:Nx,:,:]
Fn = F.yvalue[:,2:Ny+1,:]
Fs = F.yvalue[:,1:Ny,:]
Ff = F.zvalue[:,:,2:Nz+1]
Fb = F.zvalue[:,:,1:Nz]

# compute the divergence
div_x = (Fe - Fw)./dx
div_y = (Fn - Fs)./(dy.*rp)
div_z = (Ff - Fb)./dz

# define the RHS Vector
RHSdiv = zeros((Nx+2)*(Ny+2)*(Nz+2))
RHSdivx = zeros((Nx+2)*(Ny+2)*(Nz+2))
RHSdivy = zeros((Nx+2)*(Ny+2)*(Nz+2))
RHSdivz = zeros((Nx+2)*(Ny+2)*(Nz+2))

# assign the values of the RHS vector
RHSdiv[row_index] = reshape(div_x+div_y+div_z,Nx*Ny*Nz)
RHSdivx[row_index] = reshape(div_x,Nx*Ny*Nz)
RHSdivy[row_index] = reshape(div_y,Nx*Ny*Nz)
RHSdivz[row_index] = reshape(div_z,Nx*Ny*Nz)

(RHSdiv, RHSdivx, RHSdivy, RHSdivz)

end

# ===================== Gradient Term ============================
function gradientTerm(phi)
  # calculates the gradient of a variable
  # the output is a face variable
  d=phi.domain.dimension
  if d==1 || d==1.5
    dx = 0.5*(phi.domain.cellsize.x[1:end-1]+phi.domain.cellsize.x[2:end])
    FaceValue(phi.domain,
      (phi.value[2:end]-phi.value[1:end-1])./dx,
      [1.0],
      [1.0])
  elseif d==2 || d==2.5
    dx = 0.5*(phi.domain.cellsize.x[1:end-1]+phi.domain.cellsize.x[2:end])
    Ny = phi.domain.dims[2]
    dy = zeros( 1, Ny+1)
    dy[:] = 0.5*(phi.domain.cellsize.y[1:end-1]+phi.domain.cellsize.y[2:end])
    FaceValue(phi.domain,
      (phi.value[2:end,2:end-1]-phi.value[1:end-1,2:end-1])./dx,
      (phi.value[2:end-1,2:end]-phi.value[2:end-1,1:end-1])./dy,
      [1.0])
  elseif d==2.8
    dx = 0.5*(phi.domain.cellsize.x[1:end-1]+phi.domain.cellsize.x[2:end])
    Ntheta = phi.domain.dims[2]
    dtheta = zeros( 1, Ntheta+1)
    dtheta[:] = 0.5*(phi.domain.cellsize.y[1:end-1]+phi.domain.cellsize.y[2:end])
    rp = phi.domain.cellcenters.x
    FaceValue(phi.domain,
      (phi.value[2:end,2:end-1]-phi.value[1:end-1,2:end-1])./dx,
      (phi.value[2:end-1,2:end]-phi.value[2:end-1,1:end-1])./(dtheta.*rp),
      [1.0])
  elseif d==3
    Ny = phi.domain.dims[2]
    Nz = phi.domain.dims[3]
    dx = 0.5*(phi.domain.cellsize.x[1:end-1]+phi.domain.cellsize.x[2:end])
    dy= zeros( 1, Ny+1)
    dy[:] = 0.5*(phi.domain.cellsize.y[1:end-1]+phi.domain.cellsize.y[2:end])
    dz= zeros( 1, 1, Nz+1)
    dz[:] = 0.5*(phi.domain.cellsize.z[1:end-1]+phi.domain.cellsize.z[2:end])
    FaceValue(phi.domain,
      (phi.value[2:end,2:end-1,2:end-1]-phi.value[1:end-1,2:end-1,2:end-1])./dx,
      (phi.value[2:end-1,2:end,2:end-1]-phi.value[2:end-1,1:end-1,2:end-1])./dy,
      (phi.value[2:end-1,2:end-1,2:end]-phi.value[2:end-1,2:end-1,1:end-1])./dz)
  elseif d==3.2
    Ntheta = phi.domain.dims[2]
    Nz = phi.domain.dims[3]
    dx = 0.5*(phi.domain.cellsize.x[1:end-1]+phi.domain.cellsize.x[2:end])
    dy= zeros( 1, Ntheta+1)
    dy[:] = 0.5*(phi.domain.cellsize.y[1:end-1]+phi.domain.cellsize.y[2:end])
    dz= zeros( 1, 1, Nz+1)
    dz[:] = 0.5*(phi.domain.cellsize.z[1:end-1]+phi.domain.cellsize.z[2:end])
    rp = phi.domain.cellcenters.x
    FaceValue(phi.domain,
      (phi.value[2:end,2:end-1,2:end-1]-phi.value[1:end-1,2:end-1,2:end-1])./dx,
      (phi.value[2:end-1,2:end,2:end-1]-phi.value[2:end-1,1:end-1,2:end-1])./(dy.*rp),
      (phi.value[2:end-1,2:end-1,2:end]-phi.value[2:end-1,2:end-1,1:end-1])./dz)

  end
end

# ===================== Gradient Term ============================
function gradientCellTerm(phi)
  # calculates the gradient of a variable
  # the output is a face variable
  d=phi.domain.dimension
  phi_face = linearMean(phi)
  if d==1 || d==1.5
    dx = phi.domain.cellsize.x[2:end-1]
    CellVector(phi.domain,
      (phi_face.xvalue[2:end] .- phi_face.xvalue[1:end-1])./dx,
      [1.0],
      [1.0])
  elseif d==2 || d==2.5
    dx = phi.domain.cellsize.x[2:end-1]
    Ny = phi.domain.dims[2]
    dy = zeros(1, Ny)
    dy[:] = phi.domain.cellsize.y[2:end-1]
    CellVector(phi.domain,
      (phi_face.xvalue[2:end, :] .- phi_face.xvalue[1:end-1, :])./dx,
      (phi_face.yvalue[:, 2:end] .- phi_face.yvalue[:, 1:end-1])./dy,
      [1.0])
  elseif d==2.8
    dx = phi.domain.cellsize.x[2:end-1]
    Ntheta = phi.domain.dims[2]
    dtheta = zeros(1, Ntheta)
    dtheta[:] = phi.domain.cellsize.y[2:end-1]
    rp = phi.domain.cellcenters.x
    CellVector(phi.domain,
      (phi_face.xvalue[2:end, :] .- phi_face.xvalue[1:end-1, :])./dx,
      (phi_face.yvalue[:, 2:end] .- phi_face.yvalue[:, 1:end-1])./(dtheta.*rp),
      [1.0])
  elseif d==3
    Ny = phi.domain.dims[2]
    Nz = phi.domain.dims[3]
    dx = phi.domain.cellsize.x[2:end-1]
    dy= zeros(1, Ny)
    dy[:] = phi.domain.cellsize.y[2:end-1]
    dz= zeros( 1, 1, Nz)
    dz[:] = phi.domain.cellsize.z[2:end-1]
    CellVector(phi.domain,
      (phi_face.xvalue[2:end, :, :] .- phi_face.xvalue[1:end-1, :, :])./dx,
      (phi_face.yvalue[:, 2:end, :] .- phi_face.yvalue[:, 1:end-1, :])./dy,
      (phi_face.zvalue[:, :, 2:end] .- phi_face.zvalue[:, :, 1:end-1])./dz)
  elseif d==3.2
    Ntheta = phi.domain.dims[2]
    Nz = phi.domain.dims[3]
    dx = phi.domain.cellsize.x[2:end-1]
    dy= zeros(1, Ntheta)
    dy[:] = phi.domain.cellsize.y[2:end-1]
    dz= zeros(1, 1, Nz)
    dz[:] = phi.domain.cellsize.z[2:end-1]
    rp = phi.domain.cellcenters.x
    CellVector(phi.domain,
      (phi_face.xvalue[2:end, :, :] .- phi_face.xvalue[1:end-1, :, :])./dx,
      (phi_face.yvalue[:, 2:end, :] .- phi_face.yvalue[:, 1:end-1, :])./(dy.*rp),
      (phi_face.zvalue[:, :, 2:end] .- phi_face.zvalue[:, :, 1:end-1])./dz)
  end
end


text:
using Bamberg, StaticArrays
using PkgBenchmark
using Compat

n = Int(10^4)

@benchgroup "$n" begin

    gridpoints_fa = generateidealdata(n,"flat")
    @bench "findminnorm" findminnorm(gridpoints_fa, [n-1 1]')
    @bench "fnp" fnp(gridpoints_fa, [n-1 1]')

    gridpoints_na = generateidealdata(n,"nested")
    @bench "pointsvecnorm" pointsvecnorm(gridpoints_na, [n-1 1]')

    gridpoints_st = generateidealdata(n,"static-nested")
    @bench "pointsvecnorm-static" pointsvecnorm(gridpoints_st, SVector{2, Int}(n-1, 1))
    @bench "pointsvecnorm-static-opt" pointsvecnorm_opt(gridpoints_st, SVector{2, Int}(n-1, 1))
end

# res = benchmarkpkg("Bamberg")
# showall(res)

ext: jl

max_stars_repo_name: Affos/Bamberg

max_issues_repo_name: Affos/Bamberg

max_forks_repo_name: Affos/Bamberg


text:
using SFrontiers, DataFrames, Test

@testset verbose=true "Method of Moments" begin

	data= [-1.06966    1.53076     0.346522  1
			-1.55598   -1.54822     0.0296991  1
			-1.06309    0.0860199  -0.546911  1
			0.396344  -1.59922    -1.62234  1
			-0.367106  -1.31003     1.67005  1]

	 df = DataFrame(data, [:y, :x1, :x2, :cons])

    re1 = sfmodel_MoMTest(sftype(prod), sfdist(half),
                     @depvar(y), @frontier( cons, x1),
                     data = df,
                     ω = (0.5, 1, 2),
                     verbose=false )

	@test re1.coeff[1]  ≈ -0.31582 atol = 1e-5
	@test re1.jlms[1]   ≈  0.44173 atol = 1e-5
	@test re1.MoM_loglikelihood ≈ -4.74344 atol = 1e-5
end

@testset verbose=true "Part 1 (cross-section) of log-likelihood functions, marginal effect, jlms & bc index" begin

	data= [-1.06966    1.53076     0.346522  1
			-1.55598   -1.54822     0.0296991  1
			-1.06309    0.0860199  -0.546911  1
			0.396344  -1.59922    -1.62234  1
			-0.367106  -1.31003     1.67005  1]

	y = reshape(data[:,1], length(data[:,1]), 1)  # matrix
	x = z = w = q = data[:,2:3]
	v = data[:,4:4]
	PorC = 1
	nobs = 5
	dum2 = ()
	tvar = [1,2,3,1,2]
	ivar = [1,1,1,2,2]
	temp1 = [[1,2,3],[4,5]]
	temp2 = [3,2]
	idt = hcat(temp1, temp2)

	#--- N-HN ----

	posvec = (begx=1, endx=2, begz=0, endz=0, begq=0, endq=0, begw=3, endw=4, begv=5, endv=5)
	rho = ones(5,1); rho2=vec(rho);
	@test  LL_T(Half, y,x,(),(),w,v,1,5,posvec,rho, dum2, nothing) ≈ 10.32890 atol=1e-5

	jlms, bc = jlmsbc(Half, PorC, posvec, rho2, y, x, (), (), w, v, nothing)
    @test jlms[1,1] ≈ 2.27294 atol = 1e-5
    @test bc[1,1] ≈ 0.190717 atol = 1e-5

	nofvar = (nofobs=5, nofx=2, nofz=0, nofq=0,	nofw=2, nofv=1, nofpara=5, nofmarg = 0+0+2)
    sfmodel_spec(sftype(prod), sfdist(half), depvar(y), frontier(x),  σᵤ²(w), σᵥ²(v))
    aa, bb = get_marg(Half, posvec, nofvar, rho2, z, (), w)
    @test aa[1,1] ≈ 1.0199 atol=1e-5
	@test bb[1] ≈ 0.41615 atol = 1e-5

	# -- N-TN ----

	posvec = (begx=1, endx=2, begz=3, endz=4, begq=0, endq=0, begw=5, endw=6, begv=7, endv=7)
	rho = ones(7,1); rho2 = vec(rho);
	@test LL_T(Trun, y,x,z,(),w,v,1,5,posvec,rho, dum2, nothing)  ≈ 10.10842 atol=1e-5

	jlms, bc = jlmsbc(Trun, PorC, posvec, rho2, y, x, z, (), w, v, nothing)
    @test jlms[1,1] ≈ 2.72632 atol = 1e-5
    @test bc[1,1] ≈ 0.13462 atol = 1e-5

	nofvar = (nofobs=5, nofx=2, nofz=2, nofq=0,	nofw=2, nofv=1, nofpara=7, nofmarg = 2+0+2)
    sfmodel_spec(sftype(prod), sfdist(trun), depvar(y), frontier(x), μ(z), σᵤ²(z), σᵥ²(v))
    aa, bb = get_marg(Trun, posvec, nofvar, rho2, z, (), w)
    @test aa[1,1] ≈ 0.55182 atol=1e-5
	@test bb[1] ≈ 0.26175 atol = 1e-5

	# -- N-Expo ----

	posvec = (begx=1, endx=2, begz=0, endz=0, begq=0, endq=0, begw=3, endw=4, begv=5, endv=5)
	rho = ones(5,1); rho2=vec(rho)
	@test LL_T(Expo, y,x,(),(),w,v,1,5,posvec,rho, dum2, nothing) ≈ 10.66291 atol=1e-5

	jlms, bc = jlmsbc(Expo, PorC, posvec, rho2, y, x, (), (), w, v, nothing)
    @test jlms[1,1] ≈ 2.27578 atol = 1e-5
    @test bc[1,1] ≈ 0.207597 atol = 1e-5

	nofvar = (nofobs=5, nofx=2, nofz=0, nofq=0,	nofw=2, nofv=1, nofpara=5, nofmarg = 0+0+2)
    sfmodel_spec(sftype(prod), sfdist(e), depvar(y), frontier(x),  σᵤ²(w), σᵥ²(v))
    aa, bb = get_marg(Expo, posvec, nofvar, rho2, z, (), w)
    @test aa[1,1] ≈ 1.27825 atol=1e-5
	@test bb[1] ≈ 0.52157 atol = 1e-5

	# -- scaling ---

	posvec = (begx=1, endx=2, begz=3, endz=3, begq=4, endq=5, begw=6, endw=7, begv=8, endv=8)
	rho = ones(8,1); rho2=vec(rho);
	@test  LL_T(Trun_Scale, y,x,v,q,w,v,1,5,posvec,rho, dum2, nothing)  ≈ 11.93947 atol=1e-5

	jlms, bc = jlmsbc(Trun_Scale, PorC, posvec, rho2, y, x, v, q, w, v, nothing)
    @test jlms[1,1] ≈ 3.11162 atol = 1e-5
    @test bc[1,1] ≈ 0.11505 atol = 1e-5

	nofvar = (nofobs=5, nofx=2, nofz=1, nofq=2,	nofw=2, nofv=1, nofpara=8, nofmarg = 1+2+2)
	sfmodel_spec(sftype(prod), sfdist(trun_scale), depvar(y), frontier(x), hscale(q), μ(v), σᵤ²(w), σᵥ²(v))
    aa, bb = get_marg(Trun_Scale, posvec, nofvar, rho2, v, q, w)
    @test aa[1,1] ≈ 16.00399 atol=1e-5
	@test bb[1] ≈ 3.80746 atol = 1e-5

end

@testset verbose=true "Part 2 (panel) of log-likelihood functions, marginal effect, jlms & bc index" begin

	data= [-1.06966    1.53076     0.346522  1
			-1.55598   -1.54822     0.0296991  1
			-1.06309    0.0860199  -0.546911  1
			0.396344  -1.59922    -1.62234  1
			-0.367106  -1.31003     1.67005  1]

	y = reshape(data[:,1], length(data[:,1]), 1)  # matrix
	x = z = w = q = data[:,2:3]
	v = data[:,4:4]
	PorC = 1
	nobs = 5
	dum2 = ()
	tvar = [1,2,3,1,2]
	ivar = [1,1,1,2,2]
	temp1 = [[1,2,3],[4,5]]
	temp2 = [3,2]
	idt = hcat(temp1, temp2)

	# -- Wang and Ho 2010, half ----

	posvec = (begx=1, endx=2, begz=0, endz=0, begq=3, endq=4, begw=5, endw=5, begv=6, endv=6)
	rho = ones(6,1); rho2=vec(rho)
	@test LL_T(PFEWHH, data[:,1],x,(),q,(),(),1,5,posvec,rho, idt, nothing)  ≈ 5.97091 atol=1e-5

	jlms, bc = jlmsbc(PFEWHH, PorC, posvec, rho2, data[:,1], x, (), q, w, v, idt)
    @test jlms[1,1] ≈ 3.13039 atol = 1e-5
    @test bc[1,1] ≈ 0.13926 atol = 1e-5

	nofvar = (nofobs=5, nofx=2, nofz=0, nofq=2,	nofw=1, nofv=1, nofpara=6, nofmarg = 0+2+1)
	sfmodel_spec(sfpanel(TFE_WH2010), sftype(prod), sfdist(half), timevar(tvar), idvar(ivar),
	             depvar(y), frontier(x), hscale(q),  σᵤ²(v), σᵥ²(v))
    aa, bb = get_marg(PFEWHH, posvec, nofvar, rho2, v, q, v)

    @test aa[1,1] ≈ 8.59766 atol=1e-5
	@test bb[1] ≈ 2.33071 atol = 1e-5

	# --- Wang and Ho 2010, truncated ----

	posvec = (begx=1, endx=2, begz=3, endz=3, begq=4, endq=5, begw=6, endw=6, begv=7, endv=7)
	rho = ones(7,1); rho2=vec(rho)
	@test LL_T(PFEWHT, data[:,1],x,(),q,(),(),1,5,posvec,rho, idt, nothing)  ≈ 6.19169 atol=1e-5

	jlms, bc = jlmsbc(PFEWHT, PorC, posvec, rho2, data[:,1], x, v, q, w, v, idt)
    @test jlms[1,1] ≈ 3.31418 atol = 1e-5
    @test bc[1,1] ≈ 0.124392 atol = 1e-5

	nofvar = (nofobs=5, nofx=2, nofz=1, nofq=2,	nofw=1, nofv=1, nofpara=7, nofmarg = 1+2+1)
	sfmodel_spec(sfpanel(TFE_WH2010), sftype(prod), sfdist(trun), timevar(tvar), idvar(ivar),
	             depvar(y), frontier(x), hscale(q), μ(v), σᵤ²(v), σᵥ²(v))
    aa, bb = get_marg(PFEWHT, posvec, nofvar, rho2, v, q, v)

    @test aa[1,1] ≈ 11.44913 atol=1e-5
	@test bb[1] ≈ 3.10371 atol = 1e-5

	# --- time decay ----

	posvec = (begx=1, endx=2, begz=3, endz=4, begq=5, endq=6, begw=7, endw=7, begv=8, endv=8)
	rho = ones(8,1); rho2 = vec(rho)
	@test LL_T(PanDecay, y, x, z, q, w,v,1,5,posvec,rho, idt, nothing)  ≈ 11.79293 atol=1e-5

	jlms, bc = jlmsbc(PanDecay, PorC, posvec, rho2, y, x, z, q, w, v, idt)
    @test jlms[1,1] ≈ 3.282 atol = 1e-5
    @test bc[1,1] ≈ 0.10012 atol = 1e-5

	nofvar = (nofobs=5, nofx=2, nofz=2, nofq=2,	nofw=1, nofv=1, nofpara=8, nofmarg = 2+2+1)
	sfmodel_spec(sfpanel(TimeDecay), sftype(prod), sfdist(trun), timevar(tvar), idvar(ivar),
	             depvar(y), frontier(x), μ(z), gamma(q), σᵤ²(v), σᵥ²(v))
    aa, bb = get_marg(PanDecay, posvec, nofvar, rho2, z, q, v)
    @test aa[1,1] ≈ 4.38273 atol=1e-5
	@test bb[1] ≈ 1.04393 atol = 1e-5

	# --- Kumbhakar 1990 ----

	posvec = (begx=1, endx=2, begz=3, endz=4, begq=5, endq=6, begw=7, endw=7, begv=8, endv=8)
	rho = ones(8,1); rho2 = vec(rho)
	@test LL_T(PanKumb90, y, x, z, q, w,v,1,5,posvec,rho, idt, nothing)  ≈ 13.42116 atol=1e-5

	jlms, bc = jlmsbc(PanKumb90, PorC, posvec, rho2, y, x, z, q, w, v, idt)
    @test jlms[1,1] ≈ 0.23543 atol = 1e-5
    @test bc[1,1] ≈ 0.79868 atol = 1e-5

	nofvar = (nofobs=5, nofx=2, nofz=2, nofq=2,	nofw=1, nofv=1, nofpara=8, nofmarg = 2+2+1)
	sfmodel_spec(sfpanel(Kumbhakar1990), sftype(prod), sfdist(trun), timevar(tvar), idvar(ivar),
	             depvar(y), frontier(x), μ(z), gamma(q), σᵤ²(v), σᵥ²(v))
    aa, bb = get_marg(PanKumb90, posvec, nofvar, rho2, z, q, v)
    @test aa[1,1] ≈ 0.17797 atol=1e-5
	@test bb[1] ≈ 0.29235 atol = 1e-5

	# --- CSW 2014, half ---------

	posvec = (begx=1, endx=2, begz=0, endz=0, begq=0, endq=0, begw=3, endw=3, begv=4, endv=4)
	rho = ones(4,1); rho2=vec(rho)
	@test LL_T(PFECSWH, data[:,1], x, (), (), (), (),1,5,posvec,rho, idt, nothing) ≈ 8.61254 atol=1e-5

	jlms, bc = jlmsbc(PFECSWH, PorC, posvec, rho2, data[:,1], x, z, (), v, v, idt)
    @test jlms[1,1] ≈ 1.74963 atol = 1e-5
    @test bc[1,1] ≈ 0.26295 atol = 1e-5

    # This model does not allow exog determinants hence no marginal effect.

	# --- true random effect, half ----------

	posvec = (begx=1, endx=2, begz=0, endz=0, begq=3, endq=3, begw=4, endw=4, begv=5, endv=5)
	rho = ones(5,1); rho2=vec(rho)
	@test LL_T(PTREH, y, x, (), q, (), (),1,5,posvec,rho, idt, nothing) ≈ 11.70517 atol=1e-5

	jlms, bc = jlmsbc(PTREH, PorC, posvec, rho2, y, x, (), q, v, v, idt)
    @test jlms[1,1] ≈ 0 atol = 1e-5  # not available
    @test bc[1,1] ≈ 0.26956 atol = 1e-5

	# This model does not allow exog determinants hence no marginal effect.

	# --- true random effect, trun -----

	posvec = (begx=1, endx=2, begz=1, endz=1, begq=4, endq=4, begw=5, endw=5, begv=6, endv=6)
	rho = ones(6,1); rho2=vec(rho)
	@test LL_T(PTRET, y, x, v, q, (), (),1,5,posvec,rho, idt, nothing) ≈ 11.98406 atol=1e-5

	jlms, bc = jlmsbc(PTRET, PorC, posvec, rho2, y, x, v, q, v, v, idt)
    @test jlms[1,1] ≈ 0 atol = 1e-5  # not available
    @test bc[1,1] ≈ 0.18474 atol = 1e-5

	# This model does not allow exog determinants hence no marginal effect.

end

@testset verbose=true "full MLE: N-HN(z), boot (may take time)" begin

	data = [-0.930756     0.679107   -0.117138   1.0
			-0.611645     0.828413   -0.601254   1.0
			-1.31106     -0.353007    1.14228    1.0
			-0.00456725  -0.134854   -0.0886163  1.0
			-3.36723      0.586617    0.279466   1.0
			-1.49024      0.297336    0.111422   1.0
			1.33926      0.0649475  -0.357884   1.0
			0.0539396   -0.109017    0.473714   1.0
			-3.75177     -0.51421     0.300234   1.0
			-0.485357     1.57433    -0.762677   1.0
			-0.643492    -0.688907    1.42305    1.0
			0.675002    -0.762804    0.408387   1.0
			-2.48183      0.397482    0.588621   1.0
			-2.96432      0.81163    -0.296278   1.0
			-1.04406     -0.346355    0.691111   1.0
			-4.63685     -0.187573    0.506874   1.0
			-2.07605     -1.60726    -0.0569299  1.0
			-2.85586     -2.48079    -1.77102    1.0
			0.207866     2.27623     1.59062    1.0
			-1.38462      0.219693    1.39706    1.0]

	y = data[:,1];
	xvar = hcat(data[:,2], data[:,4]);
	_con = data[:,4];
	zvar = hcat(data[:,3], data[:,4]);

	myini = [0.5, 0.5, 0.5, log(3), log(1)]

	sfmodel_spec(sftype(prod), sfdist(half),
			depvar(y), frontier(xvar),
			sigma_u_2(zvar), sigma_v_2(_con))

	sfmodel_init(all_init(myini))

	sfmodel_opt(warmstart_solver(NelderMead()),
			warmstart_maxIT(200),
			main_solver(Newton()),
			main_maxIT(2000),
			tolerance(1e-8),
			verbose(false), banner(false)
			)

	res = sfmodel_fit();

	boot1, d1 = sfmodel_boot_marginal(result=res, R=10, seed=123, getBootData=true)

	@test res.coeff[1] ≈ 0.30200 atol=1e-5
	@test res.coeff[2] ≈ 0.13454 atol=1e-5
	@test res.marginal[1,1] ≈ -0.08340  atol=1e-5
	@test res.jlms[1] ≈ 1.21327  atol=1e-5
	@test res.bc[1] ≈ 0.37075 atol=1e-5
	@test boot1[1] ≈ 1.60067 atol=7e-5
	@test d1[1] ≈ -0.47043 atol=1e-5

    ci90 = sfmodel_CI(bootdata=d1, observed=res.marginal_mean, level=0.10, verbose=false);
    @test ci90[1][1] ≈ -0.38963 atol=1e-5

	pred1 = sfmodel_predict(@eq(frontier))
	@test pred1[1] ≈ 0.33963 atol=1e-5

end

@testset verbose=true "full MLE: N-Exp (may take time)" begin

	data = [-0.930756     0.679107   -0.117138   1.0
			-0.611645     0.828413   -0.601254   1.0
			-1.31106     -0.353007    1.14228    1.0
			-0.00456725  -0.134854   -0.0886163  1.0
			-3.36723      0.586617    0.279466   1.0
			-1.49024      0.297336    0.111422   1.0
			1.33926      0.0649475  -0.357884   1.0
			0.0539396   -0.109017    0.473714   1.0
			-3.75177     -0.51421     0.300234   1.0
			-0.485357     1.57433    -0.762677   1.0
			-0.643492    -0.688907    1.42305    1.0
			0.675002    -0.762804    0.408387   1.0
			-2.48183      0.397482    0.588621   1.0
			-2.96432      0.81163    -0.296278   1.0
			-1.04406     -0.346355    0.691111   1.0
			-4.63685     -0.187573    0.506874   1.0
			-2.07605     -1.60726    -0.0569299  1.0
			-2.85586     -2.48079    -1.77102    1.0
			0.207866     2.27623     1.59062    1.0
			-1.38462      0.219693    1.39706    1.0]

	y = data[:,1];
	xvar = hcat(data[:,2], data[:,4]);
	_con = data[:,4];
	zvar = hcat(data[:,3], data[:,4]);

	myini = [0.5, 0.5,  log(3), log(1)]

	sfmodel_spec(sftype(prod), sfdist(e),
			depvar(y), frontier(xvar),
			sigma_u_2(_con), sigma_v_2(_con))

	sfmodel_init(all_init(myini))

	sfmodel_opt(warmstart_solver(NelderMead()),
			warmstart_maxIT(200),
			main_solver(Newton()),
			main_maxIT(2000),
			tolerance(1e-8),
			verbose(false), banner(false)
			)

	res = sfmodel_fit();

	@test res.coeff[1] ≈ 0.34725 atol=1e-5
	@test res.coeff[2] ≈ -0.51886 atol=1e-5
	@test res.jlms[1] ≈ 0.67725 atol=1e-5
	@test res.bc[1] ≈ 0.57816 atol=1e-5

	pred1 = sfmodel_predict(@eq(frontier))
	@test pred1[1] ≈  -0.28303 atol=1e-5

end

@testset verbose=true "full MLE: N-HN (may take time)" begin

	data = [-0.930756     0.679107   -0.117138   1.0
			-0.611645     0.828413   -0.601254   1.0
			-1.31106     -0.353007    1.14228    1.0
			-0.00456725  -0.134854   -0.0886163  1.0
			-3.36723      0.586617    0.279466   1.0
			-1.49024      0.297336    0.111422   1.0
			1.33926      0.0649475  -0.357884   1.0
			0.0539396   -0.109017    0.473714   1.0
			-3.75177     -0.51421     0.300234   1.0
			-0.485357     1.57433    -0.762677   1.0
			-0.643492    -0.688907    1.42305    1.0
			0.675002    -0.762804    0.408387   1.0
			-2.48183      0.397482    0.588621   1.0
			-2.96432      0.81163    -0.296278   1.0
			-1.04406     -0.346355    0.691111   1.0
			-4.63685     -0.187573    0.506874   1.0
			-2.07605     -1.60726    -0.0569299  1.0
			-2.85586     -2.48079    -1.77102    1.0
			0.207866     2.27623     1.59062    1.0
			-1.38462      0.219693    1.39706    1.0]

	y = data[:,1];
	xvar = hcat(data[:,2], data[:,4]);
	_con = data[:,4];
	zvar = hcat(data[:,3], data[:,4]);

	myini = [0.5, 0.5,  log(3), log(1)]

	sfmodel_spec(sftype(prod), sfdist(half),
			depvar(y), frontier(xvar),
			sigma_u_2(_con), sigma_v_2(_con))

	sfmodel_init(all_init(myini))

	sfmodel_opt(warmstart_solver(NelderMead()),
			warmstart_maxIT(200),
			main_solver(Newton()),
			main_maxIT(2000),
			tolerance(1e-8),
			verbose(false), banner(false)
			)

	res = sfmodel_fit();

	@test res.coeff[1] ≈ 0.32107 atol=1e-5
	@test res.coeff[2] ≈ 0.11681 atol=1e-5
	@test res.jlms[1] ≈ 1.20146 atol=1e-5
	@test res.bc[1] ≈ 0.37468 atol=1e-5

	pred1 = sfmodel_predict(@eq(frontier))
	@test pred1[1] ≈  0.33485 atol=1e-5

end

@testset verbose=true "full MLE: N-TN(z) (may take time)" begin

	data = [-0.930756     0.679107   -0.117138   1.0
			-0.611645     0.828413   -0.601254   1.0
			-1.31106     -0.353007    1.14228    1.0
			-0.00456725  -0.134854   -0.0886163  1.0
			-3.36723      0.586617    0.279466   1.0
			-1.49024      0.297336    0.111422   1.0
			1.33926      0.0649475  -0.357884   1.0
			0.0539396   -0.109017    0.473714   1.0
			-3.75177     -0.51421     0.300234   1.0
			-0.485357     1.57433    -0.762677   1.0
			-0.643492    -0.688907    1.42305    1.0
			0.675002    -0.762804    0.408387   1.0
			-2.48183      0.397482    0.588621   1.0
			-2.96432      0.81163    -0.296278   1.0
			-1.04406     -0.346355    0.691111   1.0
			-4.63685     -0.187573    0.506874   1.0
			-2.07605     -1.60726    -0.0569299  1.0
			-2.85586     -2.48079    -1.77102    1.0
			0.207866     2.27623     1.59062    1.0
			-1.38462      0.219693    1.39706    1.0]

	y = data[:,1];
	xvar = hcat(data[:,2], data[:,4]);
	_con = data[:,4];
	zvar = hcat(data[:,3], data[:,4]);

	myini = [0.5, 0.5, 0.1,  0.5, log(3), log(1)]

	sfmodel_spec(sftype(prod), sfdist(trun), mu(_con),
			depvar(y), frontier(xvar),
			sigma_u_2(zvar), sigma_v_2(_con))

	sfmodel_init(all_init(myini))

	sfmodel_opt(warmstart_solver(NelderMead()),
			warmstart_maxIT(200),
			main_solver(Newton()), # may try warmstart_delta=0.2
			main_maxIT(2000),
			tolerance(1e-8),
			verbose(false),	banner(false)
			)

	res = sfmodel_fit();

	@test res.coeff[1] ≈ 0.29602 atol=1e-5
	@test res.coeff[2] ≈ 0.27351 atol=1e-5
	@test res.marginal[1,1] ≈ -0.07106 atol=1e-5
	@test res.jlms[1] ≈ 1.35874 atol=1e-5
	@test res.bc[1] ≈ 0.32483 atol=1e-5

	pred1 = sfmodel_predict(@eq(frontier))
	@test pred1[1] ≈  0.47454 atol=1e-5

end

@testset verbose=true "full MLE: N-TN (may take time) " begin

	data = [-0.930756     0.679107   -0.117138   1.0
			-0.611645     0.828413   -0.601254   1.0
			-1.31106     -0.353007    1.14228    1.0
			-0.00456725  -0.134854   -0.0886163  1.0
			-3.36723      0.586617    0.279466   1.0
			-1.49024      0.297336    0.111422   1.0
			1.33926      0.0649475  -0.357884   1.0
			0.0539396   -0.109017    0.473714   1.0
			-3.75177     -0.51421     0.300234   1.0
			-0.485357     1.57433    -0.762677   1.0
			-0.643492    -0.688907    1.42305    1.0
			0.675002    -0.762804    0.408387   1.0
			-2.48183      0.397482    0.588621   1.0
			-2.96432      0.81163    -0.296278   1.0
			-1.04406     -0.346355    0.691111   1.0
			-4.63685     -0.187573    0.506874   1.0
			-2.07605     -1.60726    -0.0569299  1.0
			-2.85586     -2.48079    -1.77102    1.0
			0.207866     2.27623     1.59062    1.0
			-1.38462      0.219693    1.39706    1.0]

		y = data[:,1];
		xvar = hcat(data[:,2], data[:,4]);
		_con = data[:,4];
		zvar = hcat(data[:,3], data[:,4]);

		myini = [0.5, 0.5, 0.1, log(3), log(1)]

		sfmodel_spec(sftype(prod), sfdist(trun), mu(_con),
				depvar(y), frontier(xvar),
				sigma_u_2(_con), sigma_v_2(_con))

		sfmodel_init(all_init(myini))

		sfmodel_opt(warmstart_solver(NelderMead()),
				warmstart_maxIT(200),
				main_solver(Newton()),
				main_maxIT(2000),
				tolerance(1e-8),
				verbose(false),	banner(false)
				)

		res = sfmodel_fit();

		@test res.coeff[1] ≈ 0.31051 atol=1e-5
		@test res.coeff[2] ≈ 0.31477 atol=1e-5
		@test res.jlms[1]  ≈ 1.41041 atol=1e-5
		@test res.bc[1]    ≈ 0.30956 atol=1e-5

		pred1 = sfmodel_predict(@eq(frontier))
		@test pred1[1] ≈  0.52564 atol=1e-5

end

@testset verbose=true "misc" begin

	a1 = sfmodel_MixTable(5)
	@test a1[1,3]  ≈ 10.371

end


text:
module hyporheicBiogeochemistry

using DifferentialEquations
using QuadGK
using InverseLaplace
using Turing
using CSV
using DataFrames
using MySQL

include("residenceTimeDistributions.jl")
include("reactionKinetics.jl")
include("mesocosmModel.jl")
include("prepData.jl")
include("statisticalModel.jl")

# travel time distributions
export E_powerLaw
export qCalc_powerLaw
# reaction kinetics
export foc
export fopac
# bayesian fitting
export prep_data_csv
export prep_data_db
export build_sam_model_lt
export build_sam_model_dde
export build_sam_model_dde2
export build_statistical_model

end


text:
# This file was generated, do not modify it. # hide
v = categorical(["high", "med", "low", "high", "med", "low"], ordered=true)

levels(v)


text:
module sequential

    mutable struct Sequential
        layers::Array{Any}
        add_layer::Any
        activate::Any
        initialize::Any
        update::Any
        num_layer::Int64

        loss::Float64

        default_input_shape::Tuple

        function Sequential()
            new(Any[Hidden_Input_Layer(Float32[])], Sequential_add_layer, activate_Sequential, init_Sequential, update_Sequential, 0, 0.0, ())
        end
    end

    mutable struct Hidden_Input_Layer
        output::Array{Float32}
    end
    mutable struct Hidden_Output_Layer
        δ::Array{Float32}
    end

    function Sequential_add_layer(model::Sequential, layer::Any; args...)
        push!(model.layers, layer(;input_shape=model.default_input_shape, args...))
        model.num_layer += 1
        model.default_input_shape = model.layers[end].output_shape
    end

    function init_Sequential(model::Sequential, mini_batch::Int64)
        if length(model.layers)==model.num_layer+1
            push!(model.layers, Hidden_Output_Layer(zeros(Float32, model.layers[end].output_shape..., mini_batch)))
        end
        for i in 2:length(model.layers)-1
            model.layers[i].initialize(model.layers[i], mini_batch)
        end
    end

    function activate_Sequential(model::Sequential, data::Array{Float32})
        model.layers[1].output = data
        for i in 2:length(model.layers)-1
            model.layers[i].activate(model.layers[i], model.layers[i-1].output)
        end
        return model.layers[end-1].output
    end

    function update_Sequential(model::Sequential, current_input_data::Array{Float32}, current_output_data::Array{Float32}, loss_function::Any, monitor::Any, optimizer::String, α::Float64, parameters...)
        model.activate(model, current_input_data)
        loss_function.prop!(model.layers[end].δ, model.layers[end-1].output, current_output_data)
        for i in length(model.layers)-1:-1:2
            model.layers[i].update(model.layers[i], optimizer, model.layers[i-1].output, model.layers[i+1].δ, α, parameters)
        end
        model.loss += monitor.func(model.layers[end-1].output, current_output_data)
    end
end


text:
@testset "match wave low frequency" begin

## Low frequency
    medium = Acoustic(2; ρ=1.0, c=1.0)
    ms = MultipleScattering

    volfracs = [0.1,0.5]
    species = [
        Specie(Particle(
            Acoustic(2; ρ=5., c=v), ms.Circle(0.002));
            volume_fraction=v
        )
    for v in volfracs]

    ω = 0.002
    θin = -0.2
    tol = 1e-7
    basis_order = 1

    k_eff_lows = map(species) do s
        med_eff = effective_medium(medium, [s])
        ω/med_eff.c
    end

    normal = [-1.0,0.0] # an outward normal to the surface
    materials = [Material(Halfspace(normal),s) for s in species]
    source = PlaneSource(medium, [cos(θin),sin(θin)])

    wave_eff_lows = [
        WaveMode(ω, k_eff_lows[i], source, materials[i]; tol = tol,
            basis_order=basis_order)
    for i in eachindex(species)]

    wave_effs_arr = [
        WaveModes(ω, source, materials[1];
            basis_order = basis_order,
            num_wavenumbers = 1,
            tol = tol,
            extinction_rescale = false),
        WaveModes(ω, source, materials[2];
            basis_order = basis_order,
            num_wavenumbers = 1,
            tol = tol,
            extinction_rescale = false)
    ]
    # for i in eachindex(species)]
    # Using the error causes a segmentation fault

    @test norm([ws[1].wavenumber for ws in wave_effs_arr] - k_eff_lows) < tol

    match_ws = [
        MatchPlaneWaveMode(ω, source, materials[i]; basis_order=basis_order,
            max_size = 60,
            tol = tol,
            wave_effs = wave_effs_arr[i])
    for i in eachindex(species)]

   @test maximum(match_error(match_ws[i],materials[i].shape) for i in eachindex(species)) < tol

    map(eachindex(species)) do i
        sx = abs(real(match_ws[i].PlaneWaveModes[1].wavenumber))
        x = LinRange(match_ws[i].x_match[end],2pi / sx,200)
        avg_low = DiscretePlaneWaveMode(x, wave_eff_lows[i], materials[i].shape)
        avg = DiscretePlaneWaveMode(x, match_ws[i].PlaneWaveModes, materials[i].shape)
        @test norm(avg.amplitudes[:] - avg_low.amplitudes[:]) < tol
        @test maximum(abs.(avg.amplitudes[:] - avg_low.amplitudes[:])) < tol
    end
end

max_stars_repo_name: arturgower/EffectiveWaves


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

using Markdown
using InteractiveUtils

# ╔═╡ d01a9b4f-8b55-4607-abb6-717d227fcd48
begin
	using PlutoUI, LinearAlgebra
	PlutoUI.TableOfContents(aside=true)
end

# ╔═╡ 479a40d9-b81e-442b-9962-f972b110a4dd
# Pkg.checkout("SpecialMatrices")
using SpecialMatrices

# ╔═╡ 574e86bb-159c-4141-8d8a-21bdcc9b5304
# Plot the eigenvalues (singular values) and left singular vectors
using Plots

# ╔═╡ 5eb73af5-f78a-4811-83a5-ac39063a4516
# pkg> add Arrowhead#master
using Arrowhead

# ╔═╡ 5d95dc2c-bf94-4b13-b9d5-b7b261e86cf6
md"""
# Algorithms for Structured Matrices

For matrices with some special structure, it is possible to derive versions of algorithms which are faster and/or more accurate than the standard algorithms.

__Prerequisites__

The reader should be familiar with concepts of eigenvalues and eigen vectors, singular values and singular vectors, related perturbation theory, and algorithms.

__Competences__

The reader should be able to recognise matrices which have rank-revealing decomposition and apply adequate algorithms, and to apply forward stable algorithms to arrowhead and diagonal-plus-rank-one matrices.
"""

# ╔═╡ e3e43840-d0a3-4cde-9b1e-5785759912b2
md"""
# Rank revealing decompositions

For more details, see [Z. Drmač, Computing Eigenvalues and Singular Values to High Relative Accuracy](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and [J. Demmel et al, Computing the singular value decomposition with high relative accuracy](http://www.sciencedirect.com/science/article/pii/S0024379599001342), and the references therein.

Let $A\in\mathbb{R}^{m\times n}$ with $\mathop{\mathrm{rank}}(A)=n$ (therefore, $m\geq n$) and $A=U\Sigma V^T$ its thin SVD.

## Definitions

Let $A\in\mathbb{R}^{m\times n}$.

The singular values of $A$ are (__perfectly__) __well determined to high relative accuracy__ if changing any entry $A_{kl}$ to $\theta A_{kl}$, $\theta \neq 0$, causes perturbations in singular values bounded by

$$
\min\{|\theta|,1/|\theta|\}\sigma_j \leq\tilde \sigma_j \leq
\max\{|\theta|,1/|\theta|\}\sigma_j,\quad  \forall j.$$

The __sparsity pattern__ of $A$, $Struct(A)$, is the set of indices for which $A_{kl}$ is permitted to be non-zero.

The __bipartite graph__ of the sparsity pattern $S$, $\mathcal{G}(S)$, is the graph with vertices partitioned into row vertices $r_1,\ldots,r_m$ and column vertices $c_1,\ldots,c_n$, where $r_k$ and $c_l$ are connected if and only if $(k,l)\in S$.

If $\mathcal{G}(S)$ is acyclic, matrices with sparsity pattern $S$ are __biacyclic__.

A decomposition $A=XDY^T$ with diagonal matrix $D$ is called a __rank revealing decomposition__ (RRD) if $X$ and $Y$ are full-column rank well-conditioned matrices.

__Hilbert matrix__ is a square matrix $H$ with elements $H_{ij}=\displaystyle\frac{1}{i+j-1}$.

__Hankel matrix__ is a square matrix with constant elements along skew-diagonals.

__Cauchy matrix__ is an $m\times n$ matrix $C$ with elements $C_{ij}=\displaystyle\frac{1}{x_i+y_j}$ with $x_i+y_j\neq 0$ for all $i,j$.
"""

# ╔═╡ 03f797ac-8764-4688-9e7e-e144cafb3b4c
md"""
## Facts

1. The singular values of $A$ are perfectly well determined to high relative accuracy if and only if the bipartite graph $\mathcal{G}(S)$ is acyclic (forest of trees). Examples are bidiagonal and arrowhead matrices. Sparsity pattern $S$ of acyclic bipartite graph allows at most $m+n-1$ nonzero entries. A bisection algorithm computes all singular values of biacyclic matrices to high relative accuracy.

2. An RRD of $A$ can be given or computed to high accuracy by some method. Typical methods are Gaussian elimination with complete pivoting or QR factorization with complete pivoting.

3. Let $\hat X \hat D \hat Y^T$ be the computed RRD of $A$ satisfying $|D_{jj}-\hat D_{jj}| \leq O(\varepsilon)|D_{jj}|$, $\| X-\hat X\|\leq O(\varepsilon) \|X\|$, and $\| Y-\hat Y\|\leq O(\varepsilon) \|Y\|$. The following algorithm computes the EVD of $A$ with high relative accuracy:

    1. Perform QR factorization with pivoting to get $\hat X\hat D=QRP$, where $P$ is a permutation matrix. Thus $A=QRP\hat Y^T$.
    2. Multiply $W=RP\hat Y^T$ (_NOT_ Strassen's multiplication). Thus $A=QW$ and $W$ is well-scaled from the left.
    3. Compute the SVD of $W^T=V\Sigma^T \bar U^T$ using one-sided Jacobi method. Thus $A=Q\bar U \Sigma V^T$.
    4. Multiply $U=Q\bar U$. Thus $A=U\Sigma V^T$ is the computed SVD of $A$.

4. Let $R=D'R'$, where $D'$ is such that the _rows_ of $R'$ have unit norms.  Then the following error bounds hold:

$$
\frac{|\sigma_j-\tilde\sigma_j|}{\sigma_j}\leq O(\varepsilon \kappa(R')\cdot \max\{\kappa(X),\kappa(Y)\})\leq
O(\varepsilon n^{3/2}\kappa(X)\cdot \max\{\kappa(X),\kappa(Y)\}).$$

5. Hilbert matrix is Hankel matrix and Cauchy matrix, it is symmetric positive definite and _very_ ill-conditioned.

6. Every sumbatrix of a Cauchy matrix is itself a Cauchy matrix.

7. Determinant of a square Cauchy matrix is

$$
\det(C)=\frac{\prod_{1\leq i<j\leq n}(x_j-x_i)(y_j-y_i)}
{\prod_{1\leq i,j\leq n} (x_i+y_j)}.$$

It is computed with elementwise high relative accuracy.

8. Let $A$ be square and nonsingular and let $A=LDR$ be its decomposition with diagonal $D$, lower unit-triangular $L$, and upper unit-triangular $R$. The closed formulas using quotients of minors are (see [A. S. Householder, The Theory of Matrices in Numerical Analysis](https://books.google.hr/books?id=hCre109IpRcC&printsec=frontcover&hl=hr&source=gbs_ge_summary_r&cad=0#v=onepage&q&f=false)):

$$
\begin{aligned}
D_{11}&=A_{11}, \\
D_{jj}&=\frac{\det(A_{1:j,1:j})}{\det(A_{1:j-1,1:j-1})}, \quad j=2,\ldots,n, \\
L_{jj}&=1, \\
L_{ij}&=\frac{\det(A_{[1,2,\ldots,j-1,i],[1:j]}\, )}
{\det(A_{1:j,1:j})}, \quad j < i, \\
R_{jj}&=1, \\
R_{ji}&=\frac{\det(A_{[1,2,\ldots,j],[1,2, \ldots,j-1,i]}\, )}
{\det(A_{1:j,1:j})}, \quad i > j,
\end{aligned}$$

"""

# ╔═╡ f280b119-76a7-4ee8-b6fd-608d977af0c6
md"""
## Examples

### Positive definite matrix

Let $A=DA_S D$ be strongly scaled symmetric positive definite matrix. Then Cholesky factorization with complete (diagonal) pivoting is an RRD. Consider the following three step algorithm:

1. Compute $P^T A P=LL^T$ (_Cholesky factorization with complete pivoting_).
2. Compute the $L=\bar U\Sigma V^T$ (_one-sided Jacobi, V is not needed_).
3. Set $\Lambda=\Sigma^2$ and $U=P\bar U$. Thus $A=U\Lambda U^T$ is an EVD of $A$.

The Cholesky factorization with pivoting can be implemented very fast with block algorithm (see [C. Lucas, LAPack-Style Codes for Level 2 and 3 Pivoted Cholesky Factorizations](http://www.netlib.org/lapack/lawnspdf/lawn161.pdf)).

The eigenvalues $\tilde \lambda_j$ computed using the above algorithm satisfy relative error bounds:

$$
\frac{|\lambda_j-\tilde\lambda_j|}{\lambda_j} \leq O(n\varepsilon \|A_S\|_2^{-1}).$$

"""

# ╔═╡ 28c1a9e7-6c65-4184-b41b-b5cfd17645b5
function JacobiR(A₁::AbstractMatrix)
    A=deepcopy(A₁)
    m,n=size(A)
    T=typeof(A[1,1])
    V=Matrix{T}(I,n,n)
    # Tolerance for rotation
    tol=√n*eps(T)
    # Counters
    p=n*(n-1)/2
    sweep=0
    pcurrent=0
    # First criterion is for standard accuracy, second one is for relative accuracy
    # while sweep<30 && vecnorm(A-diagm(diag(A)))>tol
    while sweep<30 && pcurrent<p
        sweep+=1
        # Row-cyclic strategy
        for i = 1 : n-1
            for j = i+1 : n
                # Compute the 2 x 2 sumbatrix of A'*A
                # F=A[:,[i,j]]'*A[:,[i,j]]
				F=view(A,:,[i,j])'*view(A,:,[i,j])
                # Check the tolerance - the first criterion is standard,
                # the second one is for relative accuracy
                # if A[i,j]!=zero(T)
                #
                if abs(F[1,2])>tol*√(F[1,1]*F[2,2])
                    # Compute c and s
                    τ=(F[2,2]-F[1,1])/(2*F[1,2])
                    t=sign(τ)/(abs(τ)+√(1+τ^2))
                    c=1/√(1+t^2)
                    s=c*t
                    G=LinearAlgebra.Givens(i,j,c,s)
                    # A*=G'
                    rmul!(A,G)
                    # V*=G'
                    rmul!(V,G)
                    pcurrent=0
                else
                    pcurrent+=1
                end
            end
        end
    end
    σ=[norm(A[:,k]) for k=1:n]
    for k=1:n
        A[:,k]./=σ[k]
    end
    # A, σ, V
    SVD(A,σ,adjoint(V))
end

# ╔═╡ bc8b94b7-7e20-4cd3-be68-7e9152fe6d7b
begin
	n=20
	import Random
	Random.seed!(421)
	B=randn(n,n)
	# Scaled matrix
	As=Matrix(Symmetric(B'*B))
	# Scaling
	D₀=exp.(50*(rand(n).-0.5))
	# Parentheses are necessary!
	A=[As[i,j]*(D₀[i]*D₀[j]) for i=1:n, j=1:n]
	issymmetric(A), cond(As), cond(A)
end

# ╔═╡ 66dedda0-24a0-48a0-9286-4ce688e5da72
# ?cholesky;

# ╔═╡ cc924252-137a-4241-b178-2eabf653ff71
md"""
We will not use the Cholesky factorization with complete pivoting. Instead, we will just sort the diagonal of $A$ in advance, which is sufficient for this example.

_Write the function for Cholesky factorization with complete pivoting as an excercise._
"""

# ╔═╡ c41e4534-6b29-4cbe-9b17-2c69c89e5570
# ?sortperm;

# ╔═╡ bcee7b90-2dd4-47b0-a781-67d35962d5f2
begin
	p=sortperm(diag(A), rev=true)
	L=cholesky(A[p,p])
end

# ╔═╡ f471d60f-1c70-4bd5-bb67-2bb17c18f3f8
U,σ,V=JacobiR(Matrix(L.L));

# ╔═╡ c474702a-8af8-4640-9f45-bca748ab3952
begin
	U₁=U[invperm(p),:]
	λ=σ.^2
end

# ╔═╡ 25143fe5-3965-468a-8cb1-7c3e8e8027ea
U'*A[p,p]*U

# ╔═╡ 86dc734d-f45d-491f-9f80-7d958e642fbd
# Due to large condition number, this is not
# as accurate as expected
Ξ=U₁'*A*U₁

# ╔═╡ 001f801a-06d2-48be-a5ad-4c07f7a92890
# Orthogonality
norm(U₁'*U₁-I)

# ╔═╡ 5196b508-5aad-4811-bfc0-1c1fa00eec37
begin
	DΞ=sqrt.(diag(Ξ))
	Ξs=[Ξ[i,j]/(DΞ[i]*DΞ[j]) for i=1:n, j=1:n]
end

# ╔═╡ bcfdae20-040f-492e-bc67-60cf5f420cc0
begin
	K=U₁*Diagonal(σ)
	K'*K
end

# ╔═╡ 3ad5b60e-3fb9-443e-90eb-40b2fd78d9a3
# Explain why is the residual so large.
norm(A*U₁-U₁*Diagonal(λ))

# ╔═╡ 3e863b09-dc68-4560-adeb-e56ab22d9afd
# Relative residual is percfect
norm(A*U₁-U₁*Diagonal(λ))/norm(A)

# ╔═╡ a5555731-c229-43de-9d88-cc3b91d7db67
[λ sort(eigvals(A),rev=true)]

# ╔═╡ c3fcbaf1-b98f-401f-9244-aad632805ecb
md"""
### Hilbert matrix

We need the newest version of the package
[SpecialMatrices.jl](https://github.com/JuliaMatrices/SpecialMatrices.jl).
"""

# ╔═╡ 8d7e4960-9b04-4f5e-b2e6-acce80a6812f
varinfo(SpecialMatrices)

# ╔═╡ 80d64fc6-33e7-4069-a6f4-8e85369dad9f
C=Cauchy([1,2,3,4,5],[0,1,2,3,4])

# ╔═╡ 804c4972-bf3d-4284-8426-b50b3ea5bb1b
H=Hilbert(5)

# ╔═╡ e9bf2b02-a60d-495b-bbdb-d87e31d8d9e0
Hf=Matrix(H)

# ╔═╡ c113b7fe-87b9-454c-aeb9-166af79bbe61
begin
	# Exact formula for the determinant of a Cauchy matrix from Fact 7.
	import LinearAlgebra.det
	function det(C::Cauchy{T}) where T
	    n=length(C.x)
	    F=triu([(C.x[j]-C.x[i])*(C.y[j]-C.y[i]) for i=1:n, j=1:n],1)
	    num=prod(F[findall(!iszero,F)])
	    den=prod([(C.x[i]+C.y[j]) for i=1:n, j=1:n])
	    if all(isinteger,C.x)&all(isinteger,C.y)
	        return num//den
	    else
	        return num/den
	    end
	end
end

# ╔═╡ c6c0e7b2-400c-478e-90af-496659fed9c4
# This is exact
det(Hf)

# ╔═╡ d057146d-4a5c-415f-9332-9b3bb7e82fc4
det(C)

# ╔═╡ 66115e0a-a5c8-4a66-958b-27d9a7865855
md"""
Compute componentwise highly accurate $A=LDL ^T$ factorization of a Hilbert (Cauchy) matrix. Using `Rational` numbers gives high accuracy.
"""

# ╔═╡ 43656da7-93db-4de3-b420-098f64dfba85
# Exact LDLT factorization from Fact 8, no pivoting.
function LDLT(C::Cauchy)
    n=length(C.x)
    T=typeof(C.x[1])
    D=Array{Rational{T}}(undef,n)
    L=Matrix{Rational{T}}(I,n,n)
    δ=[det(Cauchy(C.x[1:j],C.y[1:j])) for j=1:n]
    D[1]=map(Rational{T},C[1,1])
    D[2:n]=δ[2:n]./δ[1:n-1]
    for i=2:n
        for j=1:i-1
            L[i,j]=det(Cauchy( C.x[[1:j-1;i]], C.y[1:j])) / δ[j]
        end
    end
    L,D
end

# ╔═╡ d0b35d9b-8b87-4831-846c-b63fe8912ec7
L₁,D₁=LDLT(C)

# ╔═╡ d25dcaa1-2448-49ec-b863-b224edc4ec2c
L₁*Diagonal(D₁)*L₁' # -Matrix(H)

# ╔═╡ d668a5a2-8060-48ef-adf0-6be0b1c3859e
# L*D*L' is an RRD
cond(L₁)

# ╔═╡ a4a40d52-94fc-4ccc-8e35-09c75c5de755
cond(C)

# ╔═╡ dbfbbb00-ab00-4c85-a4ce-32468f18088a
md"""
We now compute the accurate EVD of the Hilbert matrix of order $n=100$. We cannot use the function `LDLT()` since the _computation of determinant causes overflow_ and _there is no pivoting_. Instead, we use Algorithm 3 from [J. Demmel, Computing the singular value decomposition with high relative accuracy, SIAM J. Matrix Anal. Appl, 21 (1999) 562-580](http://www.netlib.org/lapack/lawnspdf/lawn130.pdf).
"""

# ╔═╡ 0176739c-2846-4b07-a8d3-68454c1251a6
function GECP(C::Cauchy)
    n=length(C.x)
    G=Matrix(C)
    x=copy(C.x)
    y=copy(C.y)
    pr=collect(1:n)
    pc=collect(1:n)
    # Find the maximal element
    for k=1:n-1
        i,j=Tuple(argmax(abs.(G[k:n,k:n])))
        i+=k-1
        j+=k-1
        if i!=k || j!=k
            G[[i,k],:]=G[[k,i],:]
            G[:, [j,k]]=G[:, [k,j]]
            x[[k,i]]=x[[i,k]]
            y[[k,j]]=y[[j,k]]
            pr[[i,k]]=pr[[k,i]]
            pc[[j,k]]=pc[[k,j]]
        end
        for r=k+1:n
            for s=k+1:n
                G[r,s]=G[r,s]*(x[r]-x[k])*(y[s]-y[k])/
                ((x[k]+y[s])*(x[r]+y[k]))
            end
        end
        G=Matrix(Symmetric(G))
    end
    D=diag(G)
    X=tril(G,-1)*Diagonal(1.0./D)+I
    Y=Diagonal(1.0./D)*triu(G,1)+I
    X,D,Y', pr,pc
end

# ╔═╡ 8b844cf7-c574-45e9-b682-fa9cc6e9cb73
X,D,Y,pr,pc=GECP(C)

# ╔═╡ b3877ea1-d479-4d08-af8c-26e17de77106
# Check
norm(X*Diagonal(D)*Y'-Matrix(C)[pr,pc]),
norm(X[invperm(pr),:]*Diagonal(D)*Y[invperm(pc),:]'-C)

# ╔═╡ 9bff72d8-68b6-41b6-a954-ee39f90ec7b0
begin
	# Now the big test.
	n₂=100
	H₂=Hilbert(n₂)
	C₂=Cauchy(collect(1:n₂), collect(0:n₂-1))
end

# ╔═╡ 72390bf0-c5d6-46de-b8d0-0bee2cbb0af7
md"""
We need a function to compute RRD from `GECP()`
"""

# ╔═╡ 46656cee-2202-4eb9-9725-1e3e3af4df42
function RRD(C::Cauchy)
    X,D,Y,pr,pc=GECP(C)
    X[invperm(pr),:], D, Y[invperm(pc),:]
end

# ╔═╡ 6b544eaf-858a-40e5-b59b-75cfa2237a6f
X₂,D₂,Y₂=RRD(C₂);

# ╔═╡ 8eb2d15e-53c0-40ab-98ce-7b0bbf4d6cd1
# Check
norm((X₂*Diagonal(D₂)*Y₂')-C₂)

# ╔═╡ c7fd604c-9e63-4ac0-8839-82371f968ba7
cond(C₂)

# ╔═╡ 23df1242-1b43-4d52-a8ed-1b12b0d5d4b9
# Is this RRD? here X=Y
cond(X₂), cond(Y₂)

# ╔═╡ 14374962-2eea-4c8f-9390-73ef886c25d9
# Algorithm from Fact 3
function RRDSVD(X,D,Y)
    Q,R,p=qr(X*Diagonal(D),Val(true))
    W=R[:,p]*Y'
    V,σ,U₁=JacobiR(W')
    U=Q*U₁
    U,σ,V
end

# ╔═╡ d4a86d03-0afe-4b1d-8828-52eae055aa9f
U₂,σ₂,V₂=RRDSVD(X₂,D₂,Y₂);

# ╔═╡ 1dfe86bf-7553-444b-a94e-340acccf5375
# Residual and orthogonality
norm(Matrix(C₂)*V₂-U₂*Diagonal(σ₂)), norm(U₂'*U₂-I), norm(V₂'*V₂-I)

# ╔═╡ cc123376-b657-4d8f-88e0-0f7577058c2b
# Observe the differences!!
[sort(σ₂) sort(svdvals(C₂)) sort(eigvals(Matrix(C₂)))]

# ╔═╡ 337ade70-accb-417b-b655-9640ed61b375
plot(σ₂,yscale = :log10,legend=false, title="Singular values of Hilbert matrix")

# ╔═╡ f7556414-693d-4d46-889b-ed6b091a235e
begin
	# Better spy
	# some options :bluesreds,clim=(-1.0,1.0)
	import Plots.spy
	spy(A)=heatmap(A, yflip=true, color=:bluesreds, aspectratio=1)
end

# ╔═╡ 9ee70a5d-59f7-460c-92af-8872966cec40
spy(U₂)

# ╔═╡ abb64f41-6817-49e9-9d44-81ef0ce32934
md"""
# Symmetric arrowhead and DPR1 matrices

For more details, see
[N. Jakovčević Stor, I. Slapničar and J. Barlow, Accurate eigenvalue decomposition of real symmetric arrowhead matrices and applications](https://arxiv.org/abs/1302.7203) and [N. Jakovčević Stor, I. Slapničar and J. Barlow, Forward stable eigenvalue decomposition of rank-one modifications of diagonal matrices](https://arxiv.org/abs/1405.7537).
"""

# ╔═╡ b24b1c96-c73e-4c8e-83af-f35e2b9df304
md"""
## Definitions

An __arrowhead matrix__ is a real symmetric matrix of order $n$ of the form $A=\begin{bmatrix} D & z \\ z^{T} & \alpha \end{bmatrix}$, where $D=\mathop{\mathrm{diag}}(d_{1},d_{2},\ldots ,d_{n-1})$, $z=\begin{bmatrix} \zeta _{1} & \zeta _{2} & \cdots & \zeta _{n-1} \end{bmatrix}^T$ is a vector, and $\alpha$ is a scalar.

An arrowhead matrix is __irreducible__ if $\zeta _{i}\neq 0$ for all $i$ and $d_{i}\neq d_{j}$ for all $i\neq j$.

A __diagonal-plus-rank-one matrix__ (DPR1 matrix) is a real symmetric matrix of order $n$ of the form $A= D +\rho z z^T$, where $D=\mathop{\mathrm{diag}}(d_{1},d_{2},\ldots ,d_{n})$, $z=\begin{bmatrix} \zeta _{1} & \zeta _{2} & \cdots & \zeta _{n} \end{bmatrix}^T$ is a vector, and $\rho \neq 0$ is a scalar.

A DPR1 matrix is __irreducible__ if $\zeta _{i}\neq 0$ for all $i$ and $d_{i}\neq d_{j}$ for all $i\neq j$.
"""

# ╔═╡ 5e068291-cb05-4fd5-acaa-8e5d766e375a
md"""
## Facts on arrowhead matrices

Let $A$ be an arrowhead matrix of order $n$ and let $A=U\Lambda U^T$ be its EVD.

1. If $d_i$ and $\lambda_i$ are nonincreasingy ordered, the Cauchy Interlace Theorem implies

$$\lambda _{1}\geq d_{1}\geq \lambda _{2}\geq d_{2}\geq \cdots \geq d_{n-2}\geq\lambda
_{n-1}\geq d_{n-1}\geq \lambda _{n}.$$

2. If $\zeta _{i}=0$ for some $i$, then $d_{i}$ is an eigenvalue whose corresponding eigenvector is the $i$-th unit vector, and we can reduce the size of the problem by deleting the $i$-th row and column of the matrix. If $d_{i}=d_{j}$, then $d_{i}$ is an eigenvalue of $A$ (this follows from the interlacing property) and we can reduce the size of the problem by annihilating $\zeta_j$ with a Givens rotation in the $(i,j)$-plane.

3. If $A$ is irreducible, the interlacing property holds with strict inequalities.

4. The eigenvalues of $A$ are the zeros of the __Pick function__ (also, _secular equation_)

$$
f(\lambda )=\alpha -\lambda -\sum_{i=1}^{n-1}\frac{\zeta _{i}^{2}}{%
d_{i}-\lambda }=\alpha -\lambda -z^{T}(D-\lambda I)^{-1}z,$$

and the corresponding eigenvectors are

$$
U_{:,i}=\frac{x_{i}}{\left\Vert x_{i}\right\Vert _{2}},\quad
x_{i}=\begin{bmatrix}
\left( D-\lambda _{i}I\right) ^{-1}z \\
-1%
\end{bmatrix},
\quad i=1,\ldots ,n.$$

5. Let $A$ be irreducible and nonsingular. If $d_i\neq 0$ for all $i$, then $A^{-1}$ is a DPR1 matrix

$$
A^{-1}=\begin{bmatrix} D^{-1} &  \\ & 0 \end{bmatrix} + \rho uu^{T},$$

where $u=\begin{bmatrix} D^{-1}z  \\ -1 \end{bmatrix}$, and $\rho =\displaystyle\frac{1}{\alpha-z^{T}D^{-1}z}$. If $d_i=0$, then $A^{-1}$ is a permuted arrowhead matrix,

$$
A^{-1}\equiv
\begin{bmatrix}
D_{1} & 0 & 0 & z_{1} \\
0 & 0 & 0 & \zeta _{i} \\
0 & 0 & D_{2} & z_{2} \\
z_{1}^{T} & \zeta _{i} & z_{2}^{T} & \alpha
\end{bmatrix}^{-1}
= \begin{bmatrix}
D_{1}^{-1} & w_{1} & 0 & 0 \\
w_{1}^{T} & b & w_{2}^{T} & 1/\zeta _{i} \\
0 & w_{2} & D_{2}^{-1} & 0 \\
0 & 1/\zeta _{i} & 0 & 0
\end{bmatrix},$$

where

$$
\begin{aligned}
w_{1}&=-D_{1}^{-1}z_{1}\displaystyle\frac{1}{\zeta _{i}},\\
w_{2}&=-D_{2}^{-1}z_{2}\displaystyle\frac{1}{\zeta _{i}},\\
b&= \displaystyle\frac{1}{\zeta _{i}^{2}}\left(-\alpha +z_{1}^{T}D_{1}^{-1}z_{1}+z_{2}^{T}D_{2}^{-1}z_{2}\right).
\end{aligned}$$

6. The algorithm based on the following approach computes all eigenvalues and _all components_ of the corresponding eigenvectors in a forward stable manner to almost full accuracy in $O(n)$ operations per eigenpair:

    1. Shift the irreducible $A$ to $d_i$ which is closer to $\lambda_i$ (one step of bisection on $f(\lambda)$).
    2. Invert the shifted matrix.
    3. Compute the absolutely largest eigenvalue of the inverted shifted matrix and the corresponding eigenvector.

7. The algorithm is implemented in the package [Arrowhead.jl](https://github.com/ivanslapnicar/Arrowhead.jl). In certain cases, $b$ or $\rho$ need to be computed with extended precision. For this, we use the functions from file [DoubleDouble.jl](https://github.com/ivanslapnicar/Arrowhead.jl/blob/master/src/DoubleDouble.jl), originally from the the package [DoubleDouble.jl](https://github.com/simonbyrne/DoubleDouble.jl).
"""

# ╔═╡ 6d7f330f-3c88-44bd-aff9-65daa7f1ea1c
md"
## Examples

### Extended precision arithmetic
"

# ╔═╡ eea30449-2410-448d-885f-e27b7aa657c0
begin
	# Extended precision arithmetic
	a=2.0
	b=3.0
	√a
end

# ╔═╡ 21f8e11a-09fe-49a9-91cd-b118876910a8
√BigFloat(a)

# ╔═╡ 6d484719-e545-4742-9e46-515525af8244
# Double numbers according to Dekker, 1971
ad=Arrowhead.Double(a)

# ╔═╡ 310ace86-6193-43ee-bdac-fe1a79a2cb26
bd=Arrowhead.Double(b)

# ╔═╡ be91a6f1-2594-4084-8b23-312c51f35553
roota=√ad

# ╔═╡ 5762ccae-2a28-4041-8a02-4bec918e192a
rootb=√bd

# ╔═╡ 4c1a73a1-bdfe-4939-8282-54b7fa193e5a
# 30 digits should match
BigFloat(roota.hi)+BigFloat(roota.lo)

# ╔═╡ f5b83025-64e7-44d1-a841-f9a7fc8614fb
√BigFloat(a)*√BigFloat(b)

# ╔═╡ 89e578d9-7608-42f2-a2f0-0ec2c89bd510
rootab=roota*rootb

# ╔═╡ 24f1cb57-ed20-4399-996f-80294c7b1bb2
BigFloat(rootab.hi)+BigFloat(rootab.lo)

# ╔═╡ 8339adb3-b0fa-45e4-b107-62954547345b
md"""
### Random arrowhead matrix
"""

# ╔═╡ 321070d6-c0c9-482d-9a5f-0e5a4f7b522f
varinfo(Arrowhead)

# ╔═╡ 9c207955-174d-474f-be95-061c71761023
methods(GenSymArrow)

# ╔═╡ bd8ea662-f4c2-4f2d-943b-b1a5f59cbe74
begin
	n₃=8
	A₃=GenSymArrow(n₃,n₃)
end

# ╔═╡ 9fe07467-4033-4374-b4f2-b0ceecf91a03
# Elements of the type SymArrow
A₃.D, A₃.z, A₃.a, A₃.i

# ╔═╡ a3ba9c15-1d55-4e2b-b32e-4202fbeca671
E₃,info₃=eigen(A₃)

# ╔═╡ 34d1bc3f-b703-4b31-887e-36d036e5c3f9
@which eigen(A₃)

# ╔═╡ 8677d8b8-f68b-4013-864c-89902ffff8fd
# Residual and orthogonality
norm(A₃*E₃.vectors-E₃.vectors*Diagonal(E₃.values)),
norm(E₃.vectors'*E₃.vectors-I)

# ╔═╡ 360829ce-2790-4d70-8699-687650fc51b4
begin
	# Timings - notice the O(n^2)
	@time eigen(GenSymArrow(1000,1000))
	@time eigen(GenSymArrow(2000,2000))
	1
end

# ╔═╡ 9f2a5e5d-a705-4d65-8ebf-6b6d9d548999
md"""
### Numerically demanding matrix
"""

# ╔═╡ 19d95a91-8c03-4b12-a8d3-46a6861991a1
A₄=SymArrow( [ 1e10+1.0/3.0, 4.0, 3.0, 2.0, 1.0 ],
    [ 1e10 - 1.0/3.0, 1.0, 1.0, 1.0, 1.0 ], 1e10, 6 )

# ╔═╡ 4a86007c-1137-4964-be3c-da51ab7fca3c
begin
	E₄,info₄=eigen(A₄)
	[sort(E₄.values) sort(eigvals(Matrix(A₄))) sort(E₄.values)-sort(eigvals(Matrix(A₄)))]
end

# ╔═╡ 33e70165-1fab-4996-99e6-a5d27bf976c5
# Residual and orthogonality
norm(A₄*E₄.vectors-E₄.vectors*Diagonal(E₄.values)),
norm(E₄.vectors'*E₄.vectors-I)

# ╔═╡ 65e00890-d4d7-4d36-a359-4606378401b7
md"""
## Facts on DPR1 matrices

The properties of DPR1 matrices are very similar to those of arrowhead matrices. Let $A$ be a DPR1 matrix of order $n$ and let $A=U\Lambda U^T$ be its EVD.

1. If $d_i$ and $\lambda_i$ are nonincreasingy ordered and $\rho>0$, then

$$\lambda _{1}\geq d_{1}\geq \lambda _{2}\geq d_{2}\geq \cdots \geq d_{n-2}\geq\lambda
_{n-1}\geq d_{n-1}\geq \lambda _{n}\geq d_n.$$

If $A$ is irreducible, the inequalities are strict.

2. Fact 2 on arrowhead matrices holds.

3. The eigenvalues of $A$ are the zeros of the __secular equation__

$$
f(\lambda )=1+\rho\sum_{i=1}^{n}\frac{\zeta _{i}^{2}}{d_{i}-\lambda }
=1 +\rho z^{T}(D-\lambda I)^{-1}z=0,$$

and the corresponding eigenvectors are

$$
U_{:,i}=\frac{x_{i}}{\left\Vert x_{i}\right\Vert _{2}},\quad
x_{i}=( D-\lambda _{i}I) ^{-1}z.$$

4. Let $A$ be irreducible and nonsingular. If $d_i\neq 0$ for all $i$, then

$$
A^{-1}=D^{-1} +\gamma uu^{T},\quad  u=D^{-1}z, \quad \gamma =-\frac{\rho}{1+\rho z^{T}D^{-1}z},$$

is also a DPR1 matrix. If $d_i=0$, then $A^{-1}$ is a permuted arrowhead matrix,

$$
A^{-1}\equiv \left(\begin{bmatrix} D_{1} & 0 & 0 \\  0 & 0 & 0  \\  0 & 0 & D_{2} \end{bmatrix}
+\rho \begin{bmatrix} z_{1} \\ \zeta _{i} \\ z_{2}
\end{bmatrix}
\begin{bmatrix}
z_{1}^{T} & \zeta _{i} & z_{2}^{T}
\end{bmatrix}\right)^{-1}=
\begin{bmatrix}
D_{1}^{-1} & w_{1} & 0 \\
w_{1}^{T} & b & w_{2}^{T} \\
0 & w_{2} & D_{2}^{-1}
\end{bmatrix},$$

where

$$
\begin{aligned}
w_{1}&=-D_{1}^{-1}z_{1}\displaystyle\frac{1}{\zeta _{i}},\\
w_{2}&=-D_{2}^{-1}z_{2}\displaystyle\frac{1}{\zeta _{i}},\\
b &=\displaystyle\frac{1}{\zeta _{i}^{2}}\left(
\frac{1}{\rho}+z_{1}^{T}D_{1}^{-1}z_{1}+z_{2}^{T}D_{2}^{-1}z_{2}\right).
\end{aligned}$$

5. The algorithm based on the same approach as above, computes all eigenvalues and all components of the corresponding eigenvectors in a forward stable manner to almost full accuracy in $O(n)$ operations per eigenpair. The algorithm is implemented in the package `Arrowhead.jl`. In certain cases, $b$ or $\gamma$ need to be computed with extended precision.
"""

# ╔═╡ 6b3c7336-f1fc-428b-8b1e-94546b6622e8
md"""

## Examples

### Random DPR1 matrix
"""

# ╔═╡ 71d0022d-3cd8-488c-80b8-4c3c958ed8fa
begin
	n₅=8
	A₅=GenSymDPR1(n₅)
end

# ╔═╡ 4655eb27-0ae8-46db-9deb-575c89d08e7e
# Elements of the type SymDPR1
A₅.D, A₅.u, A₅.r

# ╔═╡ 6642fa72-8b1b-40b2-87f2-91199a96d7f9
begin
	E₅,info₅=eigen(A₅)
	norm(A₅*E₅.vectors-E₅.vectors*Diagonal(E₅.values)),
	norm(E₅.vectors'*E₅.vectors-I)
end

# ╔═╡ 827bb865-c361-4a98-a6d6-90dca425eddc
md"""
### Numerically demanding DPR1 matrix
"""

# ╔═╡ 56c0c2ec-220f-4e04-a0cf-a144d89dd225
# Choose one
A₆=SymDPR1( [ 1e10, 5.0, 4e-3, 0.0, -4e-3,-5.0 ], [ 1e10, 1.0, 1.0, 1e-7, 1.0,1.0 ], 1.0 )

# ╔═╡ 31114bf7-450b-458d-9a1d-c22123f5f291
begin
	E₆,info₆=eigen(A₆)
	[sort(E₆.values) sort(eigvals(Matrix(A₆)))]
end

# ╔═╡ ad305a98-6cdc-424e-8f06-386b801b53e5
# Residual and orthogonality
norm(A₆*E₆.vectors-E₆.vectors*Diagonal(E₆.values)),
norm(E₆.vectors'*E₆.vectors-I)

# ╔═╡ Cell order:
# ╟─d01a9b4f-8b55-4607-abb6-717d227fcd48
# ╟─5d95dc2c-bf94-4b13-b9d5-b7b261e86cf6
# ╟─e3e43840-d0a3-4cde-9b1e-5785759912b2
# ╟─03f797ac-8764-4688-9e7e-e144cafb3b4c
# ╟─f280b119-76a7-4ee8-b6fd-608d977af0c6
# ╠═28c1a9e7-6c65-4184-b41b-b5cfd17645b5
# ╠═bc8b94b7-7e20-4cd3-be68-7e9152fe6d7b
# ╠═66dedda0-24a0-48a0-9286-4ce688e5da72
# ╟─cc924252-137a-4241-b178-2eabf653ff71
# ╠═c41e4534-6b29-4cbe-9b17-2c69c89e5570
# ╠═bcee7b90-2dd4-47b0-a781-67d35962d5f2
# ╠═f471d60f-1c70-4bd5-bb67-2bb17c18f3f8
# ╠═c474702a-8af8-4640-9f45-bca748ab3952
# ╠═25143fe5-3965-468a-8cb1-7c3e8e8027ea
# ╠═86dc734d-f45d-491f-9f80-7d958e642fbd
# ╠═001f801a-06d2-48be-a5ad-4c07f7a92890
# ╠═5196b508-5aad-4811-bfc0-1c1fa00eec37
# ╠═bcfdae20-040f-492e-bc67-60cf5f420cc0
# ╠═3ad5b60e-3fb9-443e-90eb-40b2fd78d9a3
# ╠═3e863b09-dc68-4560-adeb-e56ab22d9afd
# ╠═a5555731-c229-43de-9d88-cc3b91d7db67
# ╟─c3fcbaf1-b98f-401f-9244-aad632805ecb
# ╠═479a40d9-b81e-442b-9962-f972b110a4dd
# ╠═8d7e4960-9b04-4f5e-b2e6-acce80a6812f
# ╠═80d64fc6-33e7-4069-a6f4-8e85369dad9f
# ╠═804c4972-bf3d-4284-8426-b50b3ea5bb1b
# ╠═e9bf2b02-a60d-495b-bbdb-d87e31d8d9e0
# ╠═c6c0e7b2-400c-478e-90af-496659fed9c4
# ╠═c113b7fe-87b9-454c-aeb9-166af79bbe61
# ╠═d057146d-4a5c-415f-9332-9b3bb7e82fc4
# ╟─66115e0a-a5c8-4a66-958b-27d9a7865855
# ╠═43656da7-93db-4de3-b420-098f64dfba85
# ╠═d0b35d9b-8b87-4831-846c-b63fe8912ec7
# ╠═d25dcaa1-2448-49ec-b863-b224edc4ec2c
# ╠═d668a5a2-8060-48ef-adf0-6be0b1c3859e
# ╠═a4a40d52-94fc-4ccc-8e35-09c75c5de755
# ╟─dbfbbb00-ab00-4c85-a4ce-32468f18088a
# ╠═0176739c-2846-4b07-a8d3-68454c1251a6
# ╠═8b844cf7-c574-45e9-b682-fa9cc6e9cb73
# ╠═b3877ea1-d479-4d08-af8c-26e17de77106
# ╠═9bff72d8-68b6-41b6-a954-ee39f90ec7b0
# ╟─72390bf0-c5d6-46de-b8d0-0bee2cbb0af7
# ╠═46656cee-2202-4eb9-9725-1e3e3af4df42
# ╠═6b544eaf-858a-40e5-b59b-75cfa2237a6f
# ╠═8eb2d15e-53c0-40ab-98ce-7b0bbf4d6cd1
# ╠═c7fd604c-9e63-4ac0-8839-82371f968ba7
# ╠═23df1242-1b43-4d52-a8ed-1b12b0d5d4b9
# ╠═14374962-2eea-4c8f-9390-73ef886c25d9
# ╠═d4a86d03-0afe-4b1d-8828-52eae055aa9f
# ╠═1dfe86bf-7553-444b-a94e-340acccf5375
# ╠═cc123376-b657-4d8f-88e0-0f7577058c2b
# ╠═574e86bb-159c-4141-8d8a-21bdcc9b5304
# ╠═337ade70-accb-417b-b655-9640ed61b375
# ╠═f7556414-693d-4d46-889b-ed6b091a235e
# ╠═9ee70a5d-59f7-460c-92af-8872966cec40
# ╟─abb64f41-6817-49e9-9d44-81ef0ce32934
# ╟─b24b1c96-c73e-4c8e-83af-f35e2b9df304
# ╟─5e068291-cb05-4fd5-acaa-8e5d766e375a
# ╟─6d7f330f-3c88-44bd-aff9-65daa7f1ea1c
# ╠═5eb73af5-f78a-4811-83a5-ac39063a4516
# ╠═eea30449-2410-448d-885f-e27b7aa657c0
# ╠═21f8e11a-09fe-49a9-91cd-b118876910a8
# ╠═6d484719-e545-4742-9e46-515525af8244
# ╠═310ace86-6193-43ee-bdac-fe1a79a2cb26
# ╠═be91a6f1-2594-4084-8b23-312c51f35553
# ╠═5762ccae-2a28-4041-8a02-4bec918e192a
# ╠═4c1a73a1-bdfe-4939-8282-54b7fa193e5a
# ╠═f5b83025-64e7-44d1-a841-f9a7fc8614fb
# ╠═89e578d9-7608-42f2-a2f0-0ec2c89bd510
# ╠═24f1cb57-ed20-4399-996f-80294c7b1bb2
# ╟─8339adb3-b0fa-45e4-b107-62954547345b
# ╠═321070d6-c0c9-482d-9a5f-0e5a4f7b522f
# ╠═9c207955-174d-474f-be95-061c71761023
# ╠═bd8ea662-f4c2-4f2d-943b-b1a5f59cbe74
# ╠═9fe07467-4033-4374-b4f2-b0ceecf91a03
# ╠═a3ba9c15-1d55-4e2b-b32e-4202fbeca671
# ╠═34d1bc3f-b703-4b31-887e-36d036e5c3f9
# ╠═8677d8b8-f68b-4013-864c-89902ffff8fd
# ╠═360829ce-2790-4d70-8699-687650fc51b4
# ╟─9f2a5e5d-a705-4d65-8ebf-6b6d9d548999
# ╠═19d95a91-8c03-4b12-a8d3-46a6861991a1
# ╠═4a86007c-1137-4964-be3c-da51ab7fca3c
# ╠═33e70165-1fab-4996-99e6-a5d27bf976c5
# ╟─65e00890-d4d7-4d36-a359-4606378401b7
# ╟─6b3c7336-f1fc-428b-8b1e-94546b6622e8
# ╠═71d0022d-3cd8-488c-80b8-4c3c958ed8fa
# ╠═4655eb27-0ae8-46db-9deb-575c89d08e7e
# ╠═6642fa72-8b1b-40b2-87f2-91199a96d7f9
# ╟─827bb865-c361-4a98-a6d6-90dca425eddc
# ╠═56c0c2ec-220f-4e04-a0cf-a144d89dd225
# ╠═31114bf7-450b-458d-9a1d-c22123f5f291
# ╠═ad305a98-6cdc-424e-8f06-386b801b53e5

max_stars_repo_path: src/Module B - Eigenvalue and Singular Value Decompositions/L7 Algorithms for Structured Matrices.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
module PPInterpolation

include("CubicPP.jl")
include("QuinticHermiteSpline.jl")
end # module


text:
using .PhysicalConstants: c, ca40_qubit_transition_frequency

export Laser

"""
    Laser(;E=0, Δ=0, ϵ=(x̂+ŷ)/√2, k=ẑ, ϕ=0, λ=729.147e-9, pointing::Array{Tuple{Int,Real}})

The physical parameters defining laser light.
**args**
* `E::Union{Function,Real}`: magnitude of the E-field in V/m
* `Δ`: static detuning from f = c/λ in [Hz]
* `ϵ::NamedTuple`: (ϵ.x, ϵ.y, ϵ.z), polarization direction, requires norm of 1
* `k::NamedTuple`: (k.x, k.y, k.z), propagation direction, requires norm of 1
* `ϕ::Union{Function,Real}`: time-dependent phase. of course, this can also be used to model a
    time-dependent detuning. Units are in radians. Note: if this is set to a function of time,
    then when constructing a Hamiltonian with the `hamiltonian` function, the units of time
    will be as specified by the `timescale` keyword argument.
* `λ::Real`: the wavelength of the laser in meters
* `pointing`: an array of `Tuple{Int,Real}` for describing ion-laser pointing configuration.
    (first element of the tuple is the index for an ion and the second element is the scaling
    factor for the laser's Efield which must be between 0 and 1).
"""
mutable struct Laser
    E::Function
    Δ::Real
    ϵ::NamedTuple{(:x,:y,:z)}
    k::NamedTuple{(:x,:y,:z)}
    ϕ::Function
    λ::Real
    pointing::Vector
    function Laser(;
            E::TE=0, Δ=0, ϵ=(x̂+ŷ)/√2, k=ẑ, ϕ::Tϕ=0, λ=c/ca40_qubit_transition_frequency,
            pointing=Array{Tuple{Int,<:Real}}(undef, 0)
        ) where {TE, Tϕ}
        rtol = 1e-6
        @assert isapprox(norm(ϵ), 1, rtol=rtol) "!(|ϵ| = 1)"
        @assert isapprox(norm(k), 1, rtol=rtol) "!(|k| = 1)"
        @assert isapprox(ndot(ϵ, k), 0, rtol=rtol) "!(ϵ ⟂ k)"
        a = pointing
        (ion_num, scaling) = map(x->getfield.(a, x), fieldnames(eltype(a)))
        @assert length(ion_num) == length(unique(ion_num)) (
                "a laser is pointing at the same ion twice"
        )
        for s in scaling
            @assert 0 <= s <= 1 "must have s ∈ [0,1]"
        end
        TE <: Number ?  Et(t) = E : Et = E
        Tϕ <: Number ? ϕt(t) = ϕ : ϕt = ϕ
        new(Et, Δ, ϵ, k, ϕt, λ, pointing)
    end
    # for copying
    Laser(E, Δ, ϵ, k, ϕ, λ, pointing) = new(E, Δ, ϵ, k, ϕ, λ, pointing)
end

function Base.:(==)(L1::Laser, L2::Laser)
    for field in fieldnames(Laser)
        if getfield(L1, field) != getfield(L2, field)
            return false
        end
    end
    true
end

function Base.print(L::Laser)
    println("λ: ", L.λ, " m")
    println("Δ: ", L.Δ, " Hz")
    println("ϵ̂: ", "(x=$(L.ϵ.x), y=$(L.ϵ.y), z=$(L.ϵ.z))")
    println("k̂: ", "(z=$(L.k.x), y=$(L.k.y), z=$(L.k.z))")
    println("E(t=0): ", "$(L.E(0.0)) V/m")
    println("ϕ(t=0): ", "$(L.ϕ(0.0)) ⋅ 2π")
end

function Base.setproperty!(L::Laser, s::Symbol, v::Tv) where{Tv}
    rtol = 1e-6
    if s == :ϵ
        @assert isapprox(norm(v), 1, rtol=rtol) "!(|ϵ| = 1)"
        if ! isapprox(ndot(L.k, v), 0, rtol=rtol)
            @warn "!(ϵ ⟂ k)"
        end
    elseif s == :k
        @assert isapprox(norm(v), 1, rtol=rtol) "!(|k| = 1)"
        if ! isapprox(ndot(v, L.ϵ), 0, rtol=rtol)
            @warn "!(ϵ ⟂ k)"
        end
    elseif s == :pointing
        b = Tv <: Vector{Tuple{Int64,Float64}} || Tv <: Vector{Tuple{Int64,Int64}}
        @assert b  "type != Vector{Tuple{Int,Real}}"
        (ion_num, scaling) = map(x->getfield.(v, x), fieldnames(eltype(v)))
        @assert length(ion_num) == length(unique(ion_num)) (
                "a laser is pointing at the same ion twice"
        )
        for s in scaling
            @assert 0 <= s <= 1 "must have s ∈ [0,1]"
        end
    elseif s == :E || s == :ϕ
        Tv <: Number ? vt(t) = v : vt = v
        Core.setproperty!(L, s, vt)
        return
    end
    Core.setproperty!(L, s, v)
end


text:
include((@__DIR__)*"/../preamble.jl")

include("reflecting_model/model_def.jl")
include("default_params.jl")

pth = mkpath((@__DIR__)*"/hitting_times/data/errors")

for ic_string in (ic,)
    ks_error_lwr = DataFrame(DG_limit=[])
    l1_error_lwr = DataFrame(DG_limit=[])
    ks_error_upr = DataFrame(DG_limit=[])
    l1_error_upr = DataFrame(DG_limit=[])
    for order in orders
        ks_error_row_lwr = zeros(1)
        l1_error_row_lwr = zeros(1)
        ks_error_row_upr = zeros(1)
        l1_error_row_upr = zeros(1)
        # for (c,model_string) in [:dg;]
        c=1
            ks_error_samples = zeros(n_boot)
            l1_error_samples = zeros(n_boot)
            for n in 1:n_boot
                cdf_boot = CSV.read("hitting_times_model/hitting_times/data/"*
                    ic_string*"/bootstrap/bootsample_"*string(n)*".csv",DataFrame)

                approx_data = CSV.read((@__DIR__)*"/hitting_times/data/"*
                    ic_string*"/order_"*string(order)*"_model_dg1.csv",DataFrame)

                ks_error_samples[n] = max(
                    maximum(abs.(approx_data.phase_1-cdf_boot.phase_1)),
                    maximum(abs.(approx_data.phase_2-cdf_boot.phase_2)),
                )
                l1_error_samples[n] = sum(abs.(approx_data.phase_1[1:end-1]+approx_data.phase_1[2:end]-
                        cdf_boot.phase_1[1:end-1]-cdf_boot.phase_1[2:end]).*diff(approx_data.t))/2 +
                    sum(abs.(approx_data.phase_2[1:end-1]+approx_data.phase_2[2:end]-
                        cdf_boot.phase_2[1:end-1]-cdf_boot.phase_2[2:end]).*diff(approx_data.t))/2
            end
            ks_q5 = quantile(ks_error_samples,qtiles[1])
            ks_q95 = quantile(ks_error_samples,qtiles[2])
            l1_q5 = quantile(l1_error_samples,qtiles[1])
            l1_q95 = quantile(l1_error_samples,qtiles[2])
            ks_error_row_lwr[c] = ks_q5
            ks_error_row_upr[c] = ks_q95
            l1_error_row_lwr[c] = l1_q5
            l1_error_row_upr[c] = l1_q95
        # end
        push!(ks_error_lwr,ks_error_row_lwr)
        push!(ks_error_upr,ks_error_row_upr)
        push!(l1_error_lwr,l1_error_row_lwr)
        push!(l1_error_upr,l1_error_row_upr)
    end
    CSV.write(pth*"/ks_lwr_"*ic_string*".csv",ks_error_lwr)
    CSV.write(pth*"/ks_upr_"*ic_string*".csv",ks_error_upr)
    CSV.write(pth*"/l1_lwr_"*ic_string*".csv",l1_error_lwr)
    CSV.write(pth*"/l1_upr_"*ic_string*".csv",l1_error_upr)
end


text:
##### Weighted var & std

## var

function Base.varzm(v::RealArray, wv::WeightVec)
    n = length(v)
    length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths."))
    w = values(wv)
    s = 0.0
    for i = 1:n
        @inbounds s += w[i] * abs2(v[i])
    end
    return s / sum(wv)
end

Base.varm(v::RealArray, m::Real, wv::WeightVec) = _moment2(v, m, wv)

function Base.var(v::RealArray, wv::WeightVec; mean=nothing)
    mean == 0 ? Base.varzm(v, wv) :
    mean == nothing ? varm(v, Base.mean(v, wv), wv) :
    varm(v, mean, wv)
end

## var along dim

Base.varzm!(R::AbstractArray, A::RealArray, wv::WeightVec, dim::Int) =
    scale!(_wsum_general!(R, Abs2Fun(), A, values(wv), dim, true), inv(sum(wv)))

Base.varm!(R::AbstractArray, A::RealArray, M::RealArray, wv::WeightVec, dim::Int) =
    scale!(_wsum_centralize!(R, Abs2Fun(), A, M, values(wv), dim, true), inv(sum(wv)))

function var!(R::AbstractArray, A::RealArray, wv::WeightVec, dim::Int; mean=nothing)
    if mean == 0
        Base.varzm!(R, A, wv, dim)
    elseif mean == nothing
        Base.varm!(R, A, Base.mean(A, wv, dim), wv, dim)
    else
        # check size of mean
        for i = 1:ndims(A)
            dA = size(A,i)
            dM = size(mean,i)
            if i == dim
                dM == 1 || throw(DimensionMismatch("Incorrect size of mean."))
            else
                dM == dA || throw(DimensionMismatch("Incorrect size of mean."))
            end
        end
        Base.varm!(R, A, mean, wv, dim)
    end
end

Base.varm(A::RealArray, M::RealArray, wv::WeightVec, dim::Int) =
    Base.varm!(Array(Float64, Base.reduced_dims(A, dim)), A, M, wv, dim)

Base.var(A::RealArray, wv::WeightVec, dim::Int; mean=nothing) =
    var!(Array(Float64, Base.reduced_dims(A, dim)), A, wv, dim; mean=mean)

## std

Base.stdm(v::RealArray, m::Real, wv::WeightVec) = sqrt(varm(v, m, wv))
Base.std(v::RealArray, wv::WeightVec; mean=nothing) = sqrt(var(v, wv; mean=mean))

Base.stdm(v::RealArray, m::RealArray, dim::Int) = Base.sqrt!(varm(v, m, dim))
Base.stdm(v::RealArray, m::RealArray, wv::WeightVec, dim::Int) = sqrt(varm(v, m, wv, dim))
Base.std(v::RealArray, wv::WeightVec, dim::Int; mean=nothing) = sqrt(var(v, wv, dim; mean=mean))

##### Fused statistics

mean_and_var(A::RealArray) = (m = mean(A); (m, varm(A, m)))
mean_and_std(A::RealArray) = (m = mean(A); (m, stdm(A, m)))

mean_and_var(A::RealArray, wv::WeightVec) = (m = mean(A, wv); (m, varm(A, m, wv)))
mean_and_std(A::RealArray, wv::WeightVec) = (m = mean(A, wv); (m, stdm(A, m, wv)))

mean_and_var(A::RealArray, dim::Int) = (m = mean(A, dim); (m, varm(A, m, dim)))
mean_and_std(A::RealArray, dim::Int) = (m = mean(A, dim); (m, stdm(A, m, dim)))

mean_and_var(A::RealArray, wv::WeightVec, dim::Int) = (m = mean(A, wv, dim); (m, varm(A, m, wv, dim)))
mean_and_std(A::RealArray, wv::WeightVec, dim::Int) = (m = mean(A, wv, dim); (m, stdm(A, m, wv, dim)))

##### General central moment

function _moment2(v::RealArray, m::Real)
    n = length(v)
    s = 0.0
    for i = 1:n
        @inbounds z = v[i] - m
        s += z * z
    end
    s / n
end

function _moment2(v::RealArray, m::Real, wv::WeightVec)
    n = length(v)
    s = 0.0
    w = values(wv)
    for i = 1:n
        @inbounds z = v[i] - m
        @inbounds s += (z * z) * w[i]
    end
    s / sum(wv)
end

function _moment3(v::RealArray, m::Real)
    n = length(v)
    s = 0.0
    for i = 1:n
        @inbounds z = v[i] - m
        s += z * z * z
    end
    s / n
end

function _moment3(v::RealArray, m::Real, wv::WeightVec)
    n = length(v)
    s = 0.0
    w = values(wv)
    for i = 1:n
        @inbounds z = v[i] - m
        @inbounds s += (z * z * z) * w[i]
    end
    s / sum(wv)
end

function _moment4(v::RealArray, m::Real)
    n = length(v)
    s = 0.0
    for i = 1:n
        @inbounds z = v[i] - m
        s += abs2(z * z)
    end
    s / n
end

function _moment4(v::RealArray, m::Real, wv::WeightVec)
    n = length(v)
    s = 0.0
    w = values(wv)
    for i = 1:n
        @inbounds z = v[i] - m
        @inbounds s += abs2(z * z) * w[i]
    end
    s / sum(wv)
end

function _momentk(v::RealArray, k::Int, m::Real)
    n = length(v)
    s = 0.0
    for i = 1:n
        @inbounds z = v[i] - m
        s += (z ^ k)
    end
    s / n
end

function _momentk(v::RealArray, k::Int, m::Real, wv::WeightVec)
    n = length(v)
    s = 0.0
    w = values(wv)
    for i = 1:n
        @inbounds z = v[i] - m
        @inbounds s += (z ^ k) * w[i]
    end
    s / sum(wv)
end

function moment(v::RealArray, k::Int, m::Real)
    k == 2 ? _moment2(v, m) :
    k == 3 ? _moment3(v, m) :
    k == 4 ? _moment4(v, m) :
    _momentk(v, k, m)
end

function moment(v::RealArray, k::Int, m::Real, wv::WeightVec)
    k == 2 ? _moment2(v, m, wv) :
    k == 3 ? _moment3(v, m, wv) :
    k == 4 ? _moment4(v, m, wv) :
    _momentk(v, k, m, wv)
end

moment(v::RealArray, k::Int) = moment(v, k, mean(v))
moment(v::RealArray, k::Int, wv::WeightVec) = moment(v, k, mean(v, wv), wv)

##### Skewness and Kurtosis

# Skewness
# This is Type 1 definition according to Joanes and Gill (1998)
function skewness(v::RealArray, m::Real)
    n = length(v)
    cm2 = 0.0   # empirical 2nd centered moment (variance)
    cm3 = 0.0   # empirical 3rd centered moment
    for i = 1:n
        @inbounds z = v[i] - m
        z2 = z * z

        cm2 += z2
        cm3 += z2 * z
    end
    cm3 /= n
    cm2 /= n
    return cm3 / sqrt(cm2 * cm2 * cm2)  # this is much faster than cm2^1.5
end

function skewness(v::RealArray, wv::WeightVec, m::Real)
    n = length(v)
    length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths."))
    cm2 = 0.0   # empirical 2nd centered moment (variance)
    cm3 = 0.0   # empirical 3rd centered moment
    w = values(wv)

    @inbounds for i = 1:n
        x_i = v[i]
        w_i = w[i]
        z = x_i - m
        z2w = z * z * w_i
        cm2 += z2w
        cm3 += z2w * z
    end
    sw = sum(wv)
    cm3 /= sw
    cm2 /= sw
    return cm3 / sqrt(cm2 * cm2 * cm2)  # this is much faster than cm2^1.5
end

skewness(v::RealArray) = skewness(v, mean(v))
skewness(v::RealArray, wv::WeightVec) = skewness(v, wv, mean(v, wv))

# (excessive) Kurtosis
# This is Type 1 definition according to Joanes and Gill (1998)
function kurtosis(v::RealArray, m::Real)
    n = length(v)
    cm2 = 0.0  # empirical 2nd centered moment (variance)
    cm4 = 0.0  # empirical 4th centered moment
    for i = 1:n
        @inbounds z = v[i] - m
        z2 = z * z
        cm2 += z2
        cm4 += z2 * z2
    end
    cm4 /= n
    cm2 /= n
    return (cm4 / (cm2 * cm2)) - 3.0
end

function kurtosis(v::RealArray, wv::WeightVec, m::Real)
    n = length(v)
    length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths."))
    cm2 = 0.0  # empirical 2nd centered moment (variance)
    cm4 = 0.0  # empirical 4th centered moment
    w = values(wv)

    @inbounds for i = 1 : n
        x_i = v[i]
        w_i = w[i]
        z = x_i - m
        z2 = z * z
        z2w = z2 * w_i
        cm2 += z2w
        cm4 += z2w * z2
    end
    sw = sum(wv)
    cm4 /= sw
    cm2 /= sw
    return (cm4 / (cm2 * cm2)) - 3.0
end

kurtosis(v::RealArray) = kurtosis(v, mean(v))
kurtosis(v::RealArray, wv::WeightVec) = kurtosis(v, wv, mean(v, wv))


text:
# wrapper `TransformedTargetModel`
# "TT" is shorthand for "TransformedTargetModel"

# # TYPES

const TT_SUPPORTED_ATOMS = (
    :Probabilistic,
    :ProbabilisticUnsupervisedDetector,
    :ProbabilisticSupervisedDetector,
    :Deterministic,
    :DeterministicUnsupervisedDetector,
    :DeterministicSupervisedDetector,
    :Interval)

# Each supported atomic type gets its own wrapper:

const TT_TYPE_GIVEN_ATOM =
    Dict(atom =>
         Symbol("TransformedTargetModel$atom") for atom in TT_SUPPORTED_ATOMS)

# ...which must have appropriate supertype:

const TT_SUPER_GIVEN_ATOM =
    Dict(atom =>
         Symbol("$(atom)Composite") for atom in TT_SUPPORTED_ATOMS)

# The type definitions:

for From in TT_SUPPORTED_ATOMS
    New = TT_TYPE_GIVEN_ATOM[From]
    To  = TT_SUPER_GIVEN_ATOM[From]
    ex = quote
        mutable struct $New{M <: $From} <: $To
            model::M
            target   # Unsupervised or callable
            inverse  # callable or `nothing`
            cache
        end
    end
    eval(ex)
end

# dict whose keys and values are now types instead of symbols:
const tt_type_given_atom = Dict()
for atom in TT_SUPPORTED_ATOMS
    atom_str = string(atom)
    type = TT_TYPE_GIVEN_ATOM[atom]
    @eval(tt_type_given_atom[$atom] = $type)
end

# not exported:
const TT_TYPES = values(tt_type_given_atom)
const TT_TYPE_EXS = values(TT_TYPE_GIVEN_ATOM)
const SomeTT = Union{TT_TYPES...}
const TTSupported = Union{keys(tt_type_given_atom)...}

# # CONSTRUCTOR

const ERR_MODEL_UNSPECIFIED = ArgumentError(
    "Expecting atomic model as argument. None specified. ")
const ERR_TOO_MANY_ARGUMENTS = ArgumentError(
    "At most one non-keyword argument, a model, allowed. ")
const PRETTY_TT_SUPPORT_OPTIONS =
    join([string("`", opt, "`") for opt in TT_SUPPORTED_ATOMS],
         ", ",
         ", and ")
const err_tt_unsupported(model) = ArgumentError(
    "Only these model supertypes support wrapping as in "*
    "`TransformedTarget`: $PRETTY_TT_SUPPORT_OPTIONS.\n"*
    "Model provided has type `$(typeof(model))`. ")
const WARN_IDENTITY_INVERSE =
    "Model being wrapped is not a deterministic predictor. "*
    "Setting `inverse=identity` to suppress inverse transformations "*
    "of predictions. "
const WARN_MISSING_INVERSE =
    "Specified `target` is not a model instance or type "*
    "and so is assumed callable (eg, is a function). "*
    "I am setting `inverse=identity` as no `inverse` specified. This means "*
    "predictions of the (semi)supervised model will be "*
    "returned on a scale different from the training target. "

"""
    TransformedTargetModel(model; target=nothing, inverse=nothing, cache=true)

Wrap the supervised or semi-supervised `model` in a transformation of
the target variable.

Here `target` one of the following:

- The `Unsupervised` model that is to transform the training target.
  By default (`inverse=nothing`) the parameters learned by this
  transformer are also used to inverse-transform the predictions of
  `model`, which means `target` must implement the `inverse_transform`
  method. If this is not the case, specify `inverse=identity` to
  suppress inversion.

- A callable object for transforming the target, such as `y ->
  log.(y)`. In this case a callable `inverse`, such as `z -> exp.(z)`,
  should be specified.

Specify `cache=false` to prioritize memory over speed, or to guarantee data
anonymity.

Specify `inverse=identity` if `model` is a probabilistic predictor, as
inverse-transforming sample spaces is not supported. Alternatively,
replace `model` with a deterministic model, such as `Pipeline(model,
y -> mode.(y))`.

### Examples

A model that normalizes the target before applying ridge regression,
with predictions returned on the original scale:

```
@load RidgeRegressor pkg=MLJLinearModels
model = RidgeRegressor()
tmodel = TransformedTargetModel(model, target=Standardizer())
```

A model that applies a static `log` transformation to the data, again
returning predictions to the original scale:

```
tmodel2 = TransformedTargetModel(model, target=y->log.(y), inverse=z->exp.(y))
```

"""
function TransformedTargetModel(args...;
                           model=nothing,
                           target=nothing,
                           inverse=nothing,
                           cache=true)
    length(args) < 2 || throw(ERR_TOO_MANY_ARGUMENTS)
    if length(args) === 1
        atom = first(args)
        model === nothing ||
            @warn "Using `model=$atom`. Ignoring specification `model=$model`. "
    else
        model === nothing && throw(ERR_MODEL_UNSPECIFIED)
        atom = model
    end
    atom isa TTSupported || throw(err_tt_unsupported(atom))
    target === nothing && throw(ERR_TARGET_NOT_SPECIFIED)

    metamodel =
        tt_type_given_atom[MMI.abstract_type(atom)](atom,
                                                    target,
                                                    inverse,
                                                    cache)
    message = clean!(metamodel)
    isempty(message) || @warn message
    return metamodel
end

_is_model_type(m) = m isa Type && m <: Model

function clean!(model::SomeTT)
    message = ""
    if prediction_type(model.model) !== :deterministic &&
        model.inverse != identity
        model.inverse = identity
        message *= WARN_IDENTITY_INVERSE
    end
    if !(model.target isa Model) &&
        !_is_model_type(model.target) && model.inverse === nothing
        model.inverse = identity
        message *= WARN_MISSING_INVERSE
    end
    return message
end

# # FIT METHOD

function MMI.fit(model::SomeTT, verbosity, X, y, other...)

    _target = model.target
    inverse = model.inverse
    atom = model.model
    cache = model.cache

    Xs = source(X)
    ys = source(y)
    others = source.(other)

    target = _is_model_type(_target) ? _target() : _target

    if target isa Model
        if target isa Static
            unsupervised_mach = machine(target, cache=cache)
        else
            unsupervised_mach = machine(target, ys, cache=cache)
        end
        z = transform(unsupervised_mach, ys)
    else
        z = node(target, ys)
    end

    supervised_mach = machine(atom, Xs, z, cache=cache)
    zhat = predict(supervised_mach, Xs)

    yhat = if target isa Model && inverse != identity
        inverse_transform(unsupervised_mach, zhat)
    else
        node(inverse, zhat)
    end

    # in case the atomic mode implements `transform`:
    W = transform(supervised_mach, Xs)

    network_mach =  machine(MMI.abstract_type(atom)(),
                            Xs,
                            ys,
                            others...;
                            predict=yhat,
                            transform=W)

    return!(network_mach, model, verbosity)
end

# # TRAINING LOSSES

function training_losses(model::SomeTT, tt_report)
    mach = first(tt_report.machines)
    return training_losses(mach)
end

## MODEL TRAITS

MMI.package_name(::Type{<:SomeTT}) = "MLJBase"
MMI.package_license(::Type{<:SomeTT}) = "MIT"
MMI.package_uuid(::Type{<:SomeTT}) = "a7f614a8-145f-11e9-1d2a-a57a1082229d"
MMI.is_wrapper(::Type{<:SomeTT}) = true
MMI.package_url(::Type{<:SomeTT}) =
    "https://github.com/JuliaAI/MLJBase.jl"

for New in TT_TYPE_EXS
    quote
        MMI.iteration_parameter(::Type{<:$New{M}}) where M =
            MLJBase.prepend(:model, iteration_parameter(M))
    end |> eval
    for trait in [:input_scitype,
        :output_scitype,
        :target_scitype,
        :fit_data_scitype,
        :predict_scitype,
        :transform_scitype,
        :inverse_transform_scitype,
        :is_pure_julia,
        :supports_weights,
        :supports_class_weights,
        :supports_online,
        :supports_training_losses,
        :is_supervised,
        :prediction_type]
        quote
            MMI.$trait(::Type{<:$New{M}}) where M = MMI.$trait(M)
        end |> eval
    end
end


text:
using Test, Publish

module CellImportedModule
using DataFrames
export imported_dataframe
const imported_dataframe = DataFrame(a=1:4, b=2:5)
end

@testset "Publish" begin
    @testset "Cells" begin
        path = joinpath(@__DIR__, "_projects/cell_tests/Project.toml")
        @test html(path; globals=Dict("publish" => Dict("cell-imports" => [CellImportedModule]))) == path
        @test pdf(path; globals=Dict("publish" => Dict("cell-imports" => [CellImportedModule]))) == path
        mktempdir() do dir
            cd(dir) do
                @test deploy(path, "deploy"; versioned=false, globals=Dict("publish" => Dict("cell-imports" => [CellImportedModule]))) == path
                cd("deploy") do
                    x1 = repr("text/plain", [1, 2])
                    x2 = repr("text/plain", [2, 2])
                    @test occursin(x1, read("README.html", String))
                    @test occursin(x2, read("README.html", String))
                end
            end
        end
    end
    @testset "MIME type display" begin
        path = joinpath(@__DIR__, "_projects/mime_tests/Project.toml")
        @test html(path) == path
        @test pdf(path) == path
    end
    @testset "User Templates" begin
        path = joinpath(@__DIR__, "_projects/user_templates/Project.toml")
        @test html(path) == path
        mktempdir() do dir
            cd(dir) do
                @test deploy(path, "deploy"; versioned=false) == path
                cd("deploy") do
                    @test isfile("README.html")
                    @test isfile("custom.template")
                    @test isfile("search.html")
                    @test isfile("search.json")
                    @test isfile("index.html")
                    @test occursin("<!--user_templates-->", read("README.html", String))
                    @test occursin("<!--user_templates-->", read("search.html", String))
                end
            end
        end
    end
    @testset "Custom Themes" begin
        path = joinpath(@__DIR__, "_projects/theme_test/Project.toml")
        @test html(path) == path
        @test pdf(path) == path
        mktempdir() do dir
            cd(dir) do
                @test deploy(path, "deploy"; versioned=false) == path
                file = joinpath("deploy", "README.html")
                cd("deploy") do
                    @test isfile("README.html")
                    @test isfile("test.css")
                    @test isfile("test.js")
                    @test isfile("search.html")
                    @test isfile("search.json")
                    @test isfile("index.html")
                    @test occursin("<!--test-template-->", read("README.html", String))
                    @test occursin("<!--test-template-->", read("search.html", String))
                end
            end
        end
    end
    @testset "Integration" begin
        version = Publish.Project(Publish).env["version"]
        @test html(Publish) == Publish
        # @test pdf(Publish) == Publish
        mktempdir() do dir
            cd(dir) do
                mkdir("project")
                setup("project"; pkg=false)
                @test isfile(joinpath("project", "Project.toml"))
                @test isfile(joinpath("project", "README.md"))

                @test deploy(Publish, "deploy") == Publish
                @test isfile(joinpath("deploy", version, "index.html"))
                rm("deploy"; recursive=true)

                @test deploy(Publish, "deploy"; named=true) == Publish
                @test isfile(joinpath("deploy", "Publish", version, "index.html"))
                rm("deploy"; recursive=true)

                @test deploy(Publish, "deploy"; versioned=false) == Publish
                @test isfile(joinpath("deploy", "index.html"))
                rm("deploy"; recursive=true)

                @test deploy(Publish, "deploy"; named=true, versioned=false) == Publish
                @test isfile(joinpath("deploy", "Publish", "index.html"))
                rm("deploy"; recursive=true)

                # @test deploy(Publish, "deploy", pdf) == Publish
                # @test isfile(joinpath("deploy", version, "Publish.pdf"))
                # rm("deploy"; recursive=true)

                # @test deploy(Publish, "deploy", pdf, html) == Publish
                # @test isfile(joinpath("deploy", version, "Publish.pdf"))
                # @test isfile(joinpath("deploy", version, "index.html"))
                # rm("deploy"; recursive=true)
            end
        end
    end
    @testset "Virtual Documents" begin
        ex = Publish.Experimental
        p = ex.Project(
            name = "virtual_doc_test",
            publish = (;
                tectonic = (; args = ["-w", "https://ttassets.z13.web.core.windows.net/tlextras-2020.0r0.tar"])
            ),
            ex.Page(
                "# Virtual Document Tests",
                ex.Table(
                    CellImportedModule.imported_dataframe,
                    caption = "imported dataframe",
                    pretty_table = (nosubheader = true,),
                ),
            ),
        )

        deploy(p, "vdoctest", pdf)
        @test isdir("vdoctest")
        isdir("vdoctest") && rm("vdoctest"; force = true, recursive = true)

        deploy(p, "vdoctest", pdf; clean = true)
        @test isdir("vdoctest")
        @test isfile("vdoctest/virtual_doc_test.pdf")
        isdir("vdoctest") && rm("vdoctest"; force = true, recursive = true)
    end
end


text:
@enum Polarisation TE TM

struct Mu
    muxx::Matrix{ComplexF64}
    muyy::Matrix{ComplexF64}
    muzz::Matrix{ComplexF64}
    muxy::Matrix{ComplexF64}
    muyx::Matrix{ComplexF64}
end
struct Eps
    epsxx::Matrix{ComplexF64}
    epsyy::Matrix{ComplexF64}
    epszz::Matrix{ComplexF64}
    epsxy::Matrix{ComplexF64}
    epsyx::Matrix{ComplexF64}
end

struct Geometry
    a1::Array{<:Real,1}
    a2::Array{<:Real,1}
    ep::Eps
    mu::Mu
end

function MaterialTensor(mat::Union{AbstractVecOrMat, AbstractFloat, Int})
    if typeof(mat) <: Int || typeof(mat) <: AbstractFloat
        return Matrix(mat*I, 3, 3)
    elseif typeof(mat) <: AbstractVector && size(mat) == (3, )
        return diagm(mat)
    elseif typeof(mat) <: AbstractMatrix && size(mat) == (3, 3)
        return mat
    else
        throw("only support 3X3 Matrix, 3X1 Vector or a number")
    end
end

function Geometry(epf::Function, muf::Function, a1::Array{<:Real,1}, a2::Array{<:Real,1}, d1::Real, d2::Real)
    epfT(x, y) = MaterialTensor(epf(x, y))
    mufT(x, y) = MaterialTensor(muf(x, y))
    P = ceil(Int, norm(a1)/(d1 - 1e-6))
    Q = ceil(Int, norm(a2)/(d2 - 1e-6))
    ps = range(-0.5, stop=0.5, length=P)[1:end-1] #why end-1?
    qs = range(-0.5, stop=0.5, length=Q)[1:end-1]
    ps = ps .+ step(ps)/2 # ensure that p and q are centered on zero
    qs = qs .+ step(qs)/2
    xs = [p*a1[1]+q*a2[1] for p in ps, q in qs]
    ys = [p*a1[2]+q*a2[2] for p in ps, q in qs]

    ep = epfT.(xs, ys)
    mu = mufT.(xs, ys)
    epsij(i, j) = [ep[x, y][i, j] for (x, y) in Iterators.product(1:size(ep)[1], 1:size(ep)[2])]
    muij(i, j) = [mu[x, y][i, j] for (x, y) in Iterators.product(1:size(mu)[1], 1:size(mu)[2])]
    epsxx, epsyy, epszz, epsxy, epsyx = epsij(1, 1), epsij(2, 2), epsij(3, 3), epsij(1, 2), epsij(2, 1)
    muxx, muyy, muzz, muxy, muyx = muij(1, 1), muij(2, 2), muij(3, 3), muij(1, 2), muij(2, 1)

    ep = Eps(epsxx, epsyy, epszz, epsxy, epsyx)
    mu = Mu(muxx, muyy, muzz, muxy, muyx)
    return Geometry(a1, a2, ep, mu)
end

function Geometry(epf::Function, muf::Function, a1_deg::Real, a2_deg::Real, d1::Real, d2::Real)
    a1 = [cosd(a1_deg), sind(a1_deg)]
    a2 = [cosd(a2_deg), sind(a2_deg)]
    return Geometry(epf, muf, a1, a2, d1, d2)
end

function Geometry(epf::Function, muf::Function, a1::Array{<:Real,1}, a2::Array{<:Real,1}, d1::Real, d2::Real, polarisation::Polarisation)
    u1, u2, u3 = transpose(vcat(a1, 0)), transpose((vcat(a2, 0))), [0 0 1]
    g=[u1*u1' u1*u2' u1*u3';u2*u1' u2*u2' u2*u3';u3*u1' u3*u2' u3*u3']
    Area = abs(u1*cross(vec(u2), vec(u3)))
    if polarisation == TE
        epfTE(x, y) = inv(MaterialTensor(epf(x, y))*inv(g)*Area)
        mufTE(x, y) = MaterialTensor(muf(x, y))*inv(g)*Area
        return Geometry(epfTE, mufTE, a1, a2, d1, d2)
    else
        epfTM(x, y) = MaterialTensor(epf(x, y))*inv(g)*Area
        mufTM(x, y) = inv(MaterialTensor(muf(x, y))*inv(g)*Area)
        return Geometry(epfTM, mufTM, a1, a2, d1, d2)
    end
end

function Geometry(epf::Function, muf::Function, a1_deg::Real, a2_deg::Real, d1::Real, d2::Real, polarisation::Polarisation)
    a1 = [cosd(a1_deg), sind(a1_deg)]
    a2 = [cosd(a2_deg), sind(a2_deg)]
    return Geometry(epf, muf, a1, a2, d1, d2, polarisation)
end


text:
using VectorSphericalWaves
using ProfileView
using Zygote

m, n, kr_r, kr_i, θ, ϕ = 6, 10, 1e5, 0, 0.2, 0.3
kr = complex(kr_r, kr_i)

θ = 0.4
@profiler for i=1:1e6; VectorSphericalWaves.πₘₙ(m, n, θ); end

@time for i=1:1e5; VectorSphericalWaves.πₘₙ(m, n, θ); end

Zygote.gradient(VectorSphericalWaves.πₘₙ, m, n, θ)

@time for i=1:1e2; N_mn_wave(m, n, complex(kr_r, kr_i), θ, ϕ; kind="regular"); end

using VectorSphericalWaves
m, n, kr_r, kr_i, θ, ϕ = 6, 10, 1e5, 0, 0.2, 0.3 ; kr = complex(kr_r, kr_i); θ = rand(1000);
N_mn_wave(m, n, kr, θ[1], ϕ; kind="regular")
M_mn_wave(m, n, kr, θ[1], ϕ; kind="regular")
@time for n = 1 : 6
    for m = -n:n
        N_mn_wave.(m, n, kr, θ, ϕ; kind="regular")
        M_mn_wave.(m, n, kr, θ, ϕ; kind="regular")
        N_mn_wave.(m, n, kr, θ, ϕ; kind="irregular")
        M_mn_wave.(m, n, kr, θ, ϕ; kind="irregular")
    end
end


text:
#! /usr/bin/env julia

gear_ratio = 37/13

steps_per_rev = 400*16

bolt_diameter = 7.5

steps_per_bolt_rev = gear_ratio * steps_per_rev # steps/rev

bolt_circumference = bolt_diameter * pi # mm/rev

steps_per_mm = steps_per_bolt_rev/bolt_circumference

println("Steps per mm of filament: $steps_per_mm")

max_stars_repo_name: sjkelly/JuliaPlayground

max_issues_repo_name: sjkelly/JuliaPlayground

max_forks_repo_name: sjkelly/JuliaPlayground


text:
using .URIParser

function URIParser.URI(p::AbstractPath; query="", fragment="")
    if isempty(p.root)
        throw(ArgumentError("$p is not an absolute path"))
    end

    b = IOBuffer()
    print(b, "file://")

    if !isempty(p.drive)
        print(b, "/")
        print(b, p.drive)
    end

    for s in p.segments
        print(b, "/")
        print(b, URIParser.escape(s))
    end

    return URIParser.URI(URIParser.URI(String(take!(b))); query=query, fragment=fragment)
end