text:
using  StaticArrays

"""
	RealVector(x1::Float64, x2::Float64, x3::Float64)

Static Float64 vector with 3 elements.
"""
const RealVector = SArray{Tuple{3},Float64,1,3}
const RealMatrix = SArray{Tuple{3,3},Float64,2,9}

"""
	VECX

Static cartesian basis vector in the X direction. Equivalent to `RealVector(1.,0.,0.)`
"""
const VECX = RealVector(1., 0., 0.)

"""
	VECY

Static cartesian basis vector in the Y direction. Equivalent to `RealVector(0.,1.,0.)`.
"""
const VECY = RealVector(0., 1., 0.)

"""
	VECZ

Static cartesian basis vector in the Z direction. Equivalent to `RealVector(0.,0.,1.)`.
"""
const VECZ = RealVector(0., 0., 1.)

"""
	VEC0

Static zero vector. Equivalent to `zero(RealVector)`.
"""
const VEC0 = zero(RealVector)

"""
	MAT0

Static 3x3 zero matrix.
"""
const MAT0 = zero(RealMatrix)

"""
	MAT1

Static 3x3 identity matrix.
"""
const MAT1 = RealMatrix(1., 0., 0.,
                        0., 1., 0.,
						0., 0., 1.)

#backward compatibility
function Vec2(x1::Float64, x2::Float64)::RealVector
	return RealVector(x1, x2, 0.)
end

#backward compatibility
function Mat2(x11::Float64, x12::Float64, x21::Float64, x22::Float64)
	return RealMatrix(x11,x12, 0.,
	                  x21,x22, 0.,
					   0., 0., 1.)
end

#backward compatibility
const Vec3 = RealVector
const Mat3 = RealMatrix

"""
	AbstractParticle

Abstract supertype for smoothed particles. Any structure with `AbstractParticle`
supertype is expected to:

* be mutable
* have field `x::RealVector` (particle position)
"""
abstract type AbstractParticle end

"""
    Shape

Supertype for geometrical shapes.
"""
abstract type Shape end

#cell for storing particle indeces
mutable struct Cell
	entries::Vector{Int64}
	lock::Threads.ReentrantLock
	Cell() = new(Int64[], Threads.ReentrantLock())
end

Base.size(c::Cell) = (length(c.entries),)
Base.IndexStyle(::Type{<:Cell}) = IndexLinear()
Base.getindex(c::Cell, i::Int64) = getindex(c.entries, i)
Base.setindex!(c::Cell, val::Any, i::Int64) = setindex!(c.entries, val, i)

"""
	ParticleSystem(T::Type, domain::Shape, h::Float64)

Struct that contains all vital information about the simulation.
The constructor specifies that:
* the simulation will use particles of type `T <: AbstractParticle`,
* Particles outside of the 'domain' can be disregarded (and will be automatically removed).
* Particles are considered neighbours if their distance is less than `h`.
Please, do not make 'domain' unnecessarily large (has negative impact on performance).
"""
struct ParticleSystem{T <: AbstractParticle}
	h::Float64
	domain::Shape

	#values for computing location keys
	key_phase::NTuple{3, Int64}
	key_lim::NTuple{3, Int64}
	key_max::Int64
	key_diff::Vector{Int64}

	particles::Vector{T}
	cell_list::Vector{Cell}
	removal_cell::Cell

	ParticleSystem(T::DataType, domain::Shape, h::Float64) =
	begin
			@assert(h  > 0.0, "invalid ParticleSystem declaration! (h must be a positive float)")
			@assert(T <: AbstractParticle, "invalid ParticleSystem declaration! ("*string(T)*" is not an AbstractParticle subtype)")
			@assert(hasfield(T, :x) && (attribute_type(T, :x) == RealVector), "invalid ParticleSystem declaration! (particles must have a field `x::RealVector`)")

			box = boundarybox(domain)
			x_min = (box.x1_min, box.x2_min, box.x3_min)
			x_max = (box.x1_max, box.x2_max, box.x3_max)
			key_phase = Int64.(floor.(x_min./h))
			key_lim = Int64.(floor.(x_max./h)) .- key_phase .+ 1
			key_max = prod(key_lim)
			key_diff = Int64[]
			if key_lim[3] == 1
				@info("2D sim")
				#simulation in 2d
				for di in -1:1, dj in -1:1
					push!(key_diff, di + key_lim[1]*dj)
				end
			else
				@info("3D sim")
				#simulation in 3d
				for di in -1:1, dj in -1:1, dk in -1:1
					push!(key_diff, di + key_lim[1]*(dj + key_lim[2]*dk))
				end
			end
			particles = Vector{T}()
			cell_list = [Cell() for _ in 1:key_max]
			removal_cell = Cell()
			return new{T}(
				h, box,
				key_phase, key_lim, key_max, key_diff,
				particles, cell_list, removal_cell
			)
	end
end

get_particle_type(::ParticleSystem{T}) where T = T

#Find the cell index for given coordinate x::RealVector
function find_key(sys::ParticleSystem, x::RealVector)::Int64
	try
		i = 1 + Int64(floor(x[1]/sys.h)) - sys.key_phase[1]
		j = 1 + Int64(floor(x[2]/sys.h)) - sys.key_phase[2]
		k = 1 + Int64(floor(x[3]/sys.h)) - sys.key_phase[3]
		return i + sys.key_lim[1]*(j-1) + sys.key_lim[1]*sys.key_lim[2]*(k-1)
	catch #Int64(x) fails when x is nan or infinity
		return -1
	end
end

"""
	ParticleField(sys::ParticleSystem, varS::Symbol)

Creates an abstract array whose ``n``-th element is the value of scalar `varS`
of ``n``-th particle in `sys`.

!!! warning "Warning"
    The indentity of ``n``-th particle in `ParticleSystem` may change
    when particles are added or removed.
"""
struct ParticleField <: AbstractArray{Float64, 1}
	sys::ParticleSystem
	varS::Symbol
end

Base.size(f::ParticleField) = (length(f.sys.particles),)
Base.getindex(f::ParticleField, i::Int64) = getproperty(f.sys.particles[i], f.varS)
Base.setindex!(f::ParticleField, val::Any, k::Int64) = setproperty!(f.sys.particles[k], f.varS, val)

function attribute_type(type::DataType, var::Symbol)
    ind = findfirst(s -> s == var, fieldnames(type))
    if typeof(ind) == Nothing
        throw("Variable "*string(var)*
            " does not exist!")
    end
    return fieldtypes(type)[ind]
end


text:
include("imports.jl")

@testset ExtendedTestSet "changehrefextension" begin
    ext = "html"
    @test changehrefextension("./hello.md", ext) == "./hello.md.html"
    @test changehrefextension("#", ext) == "#"
    @test changehrefextension("#id", ext) == "#id"
    @test changehrefextension("/root/bla#id", ext) == "/root/bla.html#id"
    @test_broken changehrefextension("/root/bla.blub.bab#id", ext) == "/root/bla.blub.bab.html#id"
end

ext: jl


text:
using BasisFunctions, BasisFunctions.Test, DomainSets, DoubleFloats

using Test

types = (Float64, LargeFloat)

#####
# Orthogonal polynomials
#####

function test_chebyshevT(T)
    println("- ChebyshevT polynomials ($T)")
    bc = ChebyshevT{T}(20)
    test_ops_generic(bc)
    x1 = T(4//10)
    @test bc[4](x1) ≈ cos(3*acos(x1))

    @test support(bc) == ChebyshevInterval{T}()

    n1 = 160
    b1 = ChebyshevT{T}(n1)

    n1 = 160
    b1 = ChebyshevT{T}(n1)
    f = exp
    e = approximate(b1, exp)
    x0 = T(1//2)
    @test abs(e(T(x0))-f(x0)) < sqrt(eps(T))

    @test hastransform(b1)
    @test hastransform(b1, interpolation_grid(b1))
    @test hastransform(b1, ChebyshevExtremae(n1))
    @test hastransform(b1, ChebyshevNodes(n1))
    @test ~hastransform(b1, PeriodicEquispacedGrid(n1, -1, 1))
    @test ~hastransform(b1, ChebyshevExtremae(n1+5))
    @test ~hastransform(b1, ChebyshevNodes(n1+5))
    @test BasisFunctions.jacobi_α(b1) ≈ -.5 ≈ BasisFunctions.jacobi_β(b1)
end

function test_chebyshevU(T)
    println("- ChebyshevU polynomials ($T)")
    n1 = 20
    b1 = ChebyshevU{T}(n1)
    test_ops_generic(b1)
    x1 = T(4//10)
    @test b1[4](x1) ≈ sin(4*acos(x1))/sin(acos(x1))

    @test support(b1) == ChebyshevInterval{T}()
    @test BasisFunctions.jacobi_α(b1) ≈ .5 ≈ BasisFunctions.jacobi_β(b1)

    test_orthogonality_orthonormality(b1, true, false, BasisFunctions.ChebyshevUWeight{T}())
    test_orthogonality_orthonormality(b1, true, false, gauss_rule(b1))
    test_orthogonality_orthonormality(b1, true, false, gauss_rule(resize(b1,length(b1)+1)))
end

function test_legendre(T)
    println("- Legendre polynomials ($T)")
    bl = Legendre{T}(10)
    test_ops_generic(bl)
    x1 = T(4//10)
    @test abs(bl[6](x1) - 0.27064) < 1e-5

    test_orthogonality_orthonormality(bl, true, false, BasisFunctions.LegendreWeight{T}())
    test_orthogonality_orthonormality(bl, true, false, gauss_rule(bl))
    test_orthogonality_orthonormality(bl, true, false, gauss_rule(resize(bl,length(bl)+1)))
end

function test_laguerre(T)
    println("- Laguerre polynomials ($T)")
    bl = Laguerre(10, T(1//3))
    test_ops_generic(bl)
    x1 = T(4//10)
    @test abs(bl[6](x1) + 0.08912346) < 1e-5

    test_orthogonality_orthonormality(bl, true, false, BasisFunctions.LaguerreWeight(T(1//3)))
    test_orthogonality_orthonormality(Laguerre(10,T(0)), true, true, BasisFunctions.LaguerreWeight(T(0)))
    test_orthogonality_orthonormality(bl, true, false, gauss_rule(bl))
    test_orthogonality_orthonormality(bl, true, false, gauss_rule(resize(bl,length(bl)+1)))
end

function test_hermite(T)
    println("- Hermite polynomials ($T)")
    bh = Hermite{T}(10)
    test_ops_generic(bh)
    x1 = T(4//10)
    @test abs(bh[6](x1) - 38.08768) < 1e-5

    test_orthogonality_orthonormality(bh, true, false, BasisFunctions.HermiteWeight{T}())
    test_orthogonality_orthonormality(bh, true, false, gauss_rule(bh))
    test_orthogonality_orthonormality(bh, true, false, gauss_rule(resize(bh,length(bh)+1)))
end

function test_jacobi(T)
    println("- Jacobi polynomials ($T)")
    bj = Jacobi(10, T(2//3), T(3//4))
    test_ops_generic(bj)
    x1 = T(4//10)
    @test abs(bj[6](x1) - 0.335157) < 1e-5
    @test BasisFunctions.jacobi_α(bj) ≈ T(2//3)
    @test BasisFunctions.jacobi_β(bj) ≈ T(3//4)

    test_orthogonality_orthonormality(bj, true, false, BasisFunctions.JacobiWeight(T(2//3), T(3//4)))
    test_orthogonality_orthonormality(bj, true, false, gauss_rule(bj))
    test_orthogonality_orthonormality(bj, true, false, gauss_rule(resize(bj,length(bj)+1)))
end

function test_ops_generic(ops)
    T = codomaintype(ops)
    tol = test_tolerance(T)

    x = fixed_point_in_domain(ops)
    z1 = BasisFunctions.unsafe_eval_element(ops, length(ops), x)
    z2 = recurrence_eval(ops, length(ops), x)
    @test abs(z1-z2) < tol
    a,b = monic_recurrence_coefficients(ops)
    z3 = monic_recurrence_eval(a, b, length(ops), x)
    γ = leading_order_coefficient(ops, length(ops))
    @test abs(z1 - γ*z3) < tol

    d1 = BasisFunctions.unsafe_eval_element_derivative(ops, length(ops), x, 1)
    d2 = recurrence_eval_derivative(ops, length(ops), x)
    @test abs(d1-d2) < tol

    # if prectype(ops) == Float64
        # We only do these tests for Float64 because eig currently does not support BigFloat
        r = roots(ops)
        @test maximum(abs.(BasisFunctions.unsafe_eval_element.(Ref(ops), length(ops)+1, r))) < 100tol

        m = gauss_rule(ops)
        x = points(m)
        w = BasisFunctions.weights(m)
        @test abs(sum(w) - first_moment(ops)) < tol
    # end
end

# for T in types
for T in (Float64, LargeFloat)
    @testset "Orthogonal polynomials ($T)" begin
        println("Classical orthogonal polynomials ($T):")
        test_chebyshevT(T)
        test_chebyshevU(T)
        test_legendre(T)
        test_laguerre(T)
        test_hermite(T)
        test_jacobi(T)
    end
    println()
end

@testset "Orthogonality of orthogonal polynomials" begin
    OPSs = [ChebyshevT, ChebyshevU, Legendre, Hermite, Jacobi, Laguerre]
    for ops in OPSs, n in (5,6), T in (Float64,LargeFloat)
        B = ops{T}(n)
        test_orthogonality_orthonormality(B, gauss_rule(B))
        test_orthogonality_orthonormality(B, gauss_rule(resize(B,2n)))
        test_orthogonality_orthonormality(B, gauss_rule(resize(B,n-1)))
    end
end


text:
using VectorAutoregressions
using Base.Test

#-----------Set base-path------------------------------------------
path = Pkg.dir("VectorAutoregressions")

#--------------------Test VAR model---------------------------------------------------------
# comparison against http://www-personal.umich.edu/~lkilian/figure9_1_chol.zip
y = readdlm(joinpath(path,"test","cholvar_test_data.csv"),',')

V = VAR(y, 4, true)

@test V.p == 4
@test V.Σ ≈ [
        312.524637669663 0.773643992117187 0.919341486304653
        0.773643992117187 0.0515304724687740 0.0148802276681216
        0.919341486304653 0.0148802276681216 0.557021891337960
        ]

@test V.β[:,1:4] ≈ [
        -0.680676727884458 -0.00637960178073601 0.936519611007488 1.75104737240770
        -0.0152242907833718 0.00177981748889471 0.599086054531876 0.0281955980964510
        0.436329751995733 -0.00236721376013309 -0.218737236527547 0.332361646200618
        ]

@test V.ϵ[:,1:4] ≈ [
        70.3934221013155 -0.443580106581005 4.21164200500382 1.22299717886091
        0.189207494365181 0.287694213161629 0.692414550347607 0.405418261036420
        -1.39379479232130 0.277340573620367 -0.573151724739961 -0.395501452903569
        ]

mIRFa = IRFs_a(V,4,true)

@test [mIRFa.IRF[1,1:4]'; mIRFa.IRF[4,1:4]'; mIRFa.IRF[7,1:4]'] ≈ [
        17.6783663744607 0.0192642533976305 -2.85156060358692 0.594909295963393
        0.0437621879606954 0.0591478592307352 0.113331057827868 0.109932348126692
        0.0520037579735191 -0.0341368375885218 -0.111284843056732 -0.0491881671489620
        ]

mIRFb = IRFs_b(V,4,10,true)

@test [mIRFb.IRF[1,1:4]'; mIRFb.IRF[4,1:4]'; mIRFb.IRF[7,1:4]'] ≈ [
        17.6783663744607 0.0192642533976305 -2.85156060358692 0.594909295963393
        0.0437621879606954 0.0591478592307352 0.113331057827868 0.109932348126692
        0.0520037579735191 -0.0341368375885218 -0.111284843056732 -0.0491881671489620
        ]

#--------------------Test BVAR model---------------------------------------------------------
# comparison against http://cremfi.econ.qmul.ac.uk/outgoing/bvar.zip
include(joinpath(path,"src","bvar.jl"))
y = readdlm(joinpath(path,"test","bvar_data.csv"), ',')
y = y[:,1:3]
prior = Hyperparameter()
mForecast = fit_bvar(y,prior)

@test isapprox([0.8510 1.4081 2.2570 2.3415 2.4622 2.5835 2.6867 2.5790 2.5897 2.5767],median(mForecast[:,:,1],1); atol = 1)
@test isapprox([1.9614 2.2587 1.8328 1.8745 2.0870 2.2014 2.3303 2.5225 2.5453 2.59387],median(mForecast[:,:,2],1); atol = 1)
@test isapprox([-0.3827 -0.2272 -0.1532 -0.0735 0.0784 0.1620 0.3764 0.5247 0.7264 0.8861],median(mForecast[:,:,3],1); atol = 1)

#--------------------Test Local Projection IRFs---------------------------------------------------------
# comparison against Kilian and Kim (2011) generated data from VAR(12)

#-----------Load data----------------------------------------------
y      = readcsv(joinpath(path,"test","lp_data.csv"))
irfl   = readcsv(joinpath(path,"test","lp_test_lp_chol_irf.csv"))
stdl   = readcsv(joinpath(path,"test","lp_test_lp_std.csv"))
stdv   = readcsv(joinpath(path,"test","lp_test_var_stdv.csv"))
COVsig = readcsv(joinpath(path,"test","lp_test_var_covsig.csv"))
rfirfl = readcsv(joinpath(path,"test","lp_test_lp_rf_irf.csv"))
mlag   = readcsv(joinpath(path,"test","lp_test_lp_laglength.csv"))

#-----------Hyperparameter-----------------------------------------
const p = 12    # lag length
const H = 24    # horizon
const intercept = true

#-----------Reduced form local projection IRFs-------------------
RF_IRFs = IRFs_localprojection(y, p, H)
mRFIRFs,CI = RF_IRFs.IRF, RF_IRFs.CI

#-----------Structural local projection IRFs-------------------
# Using a VAR(pbar) as auxiliary model for cholesky identification
V = VAR(y,p,intercept)
A0inv = V.Σ |> λ -> cholfact(λ)[:L] |> full
mStd,mCov_Σ = irf_ci_asymptotic(V, H, V.inter)

IRF = IRFs_localprojection(y, p, H, A0inv, mCov_Σ)
mIRF,CI = IRF.IRF, IRF.CI
CIl,CIh = CI.CIl, CI.CIh

@test isapprox(mRFIRFs,rfirfl)
@test isapprox(A0inv,[0.721465    0.0          0.0        0.0
               -0.0842337   2.68627      0.0        0.0
                2.61134     0.724746    26.8394     0.0
                0.0886017  -0.00940385   0.0693017  0.486568],
                atol = 0.0001)
@test isapprox(V.Σ, [0.520512   -0.0607717    1.88399   0.063923
              -0.0607717   7.22317      1.7269   -0.0327246
               1.88399     1.7269     727.695     2.08457
               0.063923   -0.0327246    2.08457   0.249489],
               atol = 0.0001)
@test isapprox(mStd, stdv, atol = 0.1)
@test isapprox(mCov_Σ, COVsig, atol = 0.1)
@test isapprox(mIRF, irfl)
@test isapprox(CIl,(irfl - 1.96*stdl),atol = 0.1)
@test isapprox(CIh,(irfl + 1.96*stdl),atol = 0.1)

#-----------Test LP with different lag-length at different horizon-------
const vP = p*ones(Int64,H) # vector of lag-length

#-----------Reduced form local projection IRFs-------------------
RF_IRFs = IRFs_localprojection(y, vP, H)
mRFIRFs,CI = RF_IRFs.IRF, RF_IRFs.CI

#-----------Structural local projection IRFs-------------------
# Using a VAR(pbar) as auxiliary model for cholesky identification
V = VAR(y,p,intercept)
A0inv = V.Σ |> λ -> cholfact(λ)[:L] |> full
mStd,mCov_Σ = irf_ci_asymptotic(V, H, V.inter)

IRF = IRFs_localprojection(y, vP, H, A0inv, mCov_Σ)
mIRF,CI = IRF.IRF, IRF.CI
CIl,CIh = CI.CIl, CI.CIh

@test isapprox(mRFIRFs,rfirfl)
@test isapprox(mIRF, irfl)
@test isapprox(CIl,(irfl - 1.96*stdl),atol = 0.1)
@test isapprox(CIh,(irfl + 1.96*stdl),atol = 0.1)

#-----------Test LP lag-length selecion procedure----------------
const pbar = 12      # max order of lag to test

#-----------Select lag-length with AIC, BIC, AICC, HQC-----------
mVARlag = (var_lagorder(y,pbar,ic) for ic in ["aic","bic","aicc","hqc"]) |> λ -> hcat(collect(λ)...)
mLPlag  = (lp_lagorder(y,pbar,H,ic) for ic in ["aic","bic","aicc","hqc"]) |> λ -> hcat(collect(λ)...)

@test mLPlag == mlag
@test mVARlag == [12 2 12 3]


text:
#
# Entropy-based cross frequency coupling
#
# See Hurtado, J. M., Rubchinsky, L. L., & Sigvardt, K. A. (2004).
# Statistical Method for Detection of Phase-Locking Episodes in Neural
# Oscillations. Journal of Neurophysiology, 91(4), 1883–1898.
# doi:10.1152/jn.00853.2003
#
# Tort, A. B. L., Kramer, M. A., Thorn, C., Gibson, D. J., Kubota, Y.,
# Graybiel, A. M., & Kopell, N. J. (2008). Dynamic cross-frequency
# couplings of local field potential oscillations in rat striatum and
# hippocampus during performance of a T-maze task. Proceedings of the
# National Academy of Sciences, 105(51), 20517–20522.
# doi:10.1073/pnas.0810524105
struct HurtadoModulationIndex <: PairwiseStatistic
    nbins::UInt8
    hmax::Float64

    HurtadoModulationIndex(nbins) = new(nbins, log(nbins))
end
HurtadoModulationIndex() = HurtadoModulationIndex(18)
Base.eltype(::HurtadoModulationIndex, X::AbstractArray{Complex{T}}) where {T<:Real} = T

allocwork(t::HurtadoModulationIndex, X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}}) where {T<:Real} =
    (Array(UInt8, size(X, 1), size(X, 2)), Array(Int32, t.nbins, nchannels(X)),
     Array(T, size(X, 1)), Array(Float64, t.nbins))

function binphases!(phasebin::Matrix{UInt8}, ninbin::Matrix{Int32}, X::AbstractVecOrMat{Complex{T}}) where T<:Real
    fill!(ninbin, UInt8(0))
    m = size(ninbin, 1)/2pi
    for j = 1:size(X, 2), i = 1:size(X, 1)
        @inbounds bin = phasebin[i, j] = ceil(UInt8, (angle(X[i, j])+pi)*m)
        @inbounds ninbin[bin, j] += 1
    end
end

function sumphases!(meanamp::Vector{Float64}, phasebin::Matrix{UInt8}, amp::Vector{T}, iphase::Int) where T<:Real
    # Empty phase bins
    for i = 1:length(meanamp)
        @inbounds meanamp[i] = 0
    end

    # Compute sum of amplitudes for each phase bin
    for i = 1:size(phasebin, 1)
        @inbounds meanamp[phasebin[i, iphase]] += amp[i]
    end
end

function finish_mi!(meanamp::Vector{Float64}, ninbin::AbstractVecOrMat{Int32}, iphase::Int, hmax)
    # Normalize to mean amplitudes
    d = 0.
    for ibin = 1:length(meanamp)
        @inbounds d += (meanamp[ibin] /= ninbin[ibin, iphase])
    end
    d = inv(d)

    # Compute h
    h = 0.
    for ibin = 1:length(meanamp)
        @inbounds p = meanamp[ibin]*d
        h -= p*log(p)
    end

    # Compute modulation index
    (hmax - h)/hmax
end

# X is phase, Y is amplitude
function computestat!(t::HurtadoModulationIndex, out::AbstractMatrix{T},
                      work::Tuple{Matrix{UInt8}, Matrix{Int32}, Vector{T}, Vector{Float64}},
                      X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}}) where T<:Real
    nbins = t.nbins
    phasebin, ninbin, amp, meanamp = work
    chkinput(out, X, Y)
    size(ninbin, 1) == length(meanamp) == nbins || throw(ArgumentError("invalid work object"))

    binphases!(phasebin, ninbin, X)

    # Compute phase-amplitude coupling
    for k = 1:size(Y, 2)
        # Compute amplitudes
        for i = 1:size(X, 1)
            @inbounds amp[i] = abs(Y[i, k])
        end

        for j = 1:size(X, 2)
            sumphases!(meanamp, phasebin, amp, j)
            out[j, k] = finish_mi!(meanamp, ninbin, j, t.hmax)
        end
    end

    out
end

allocwork(t::AbstractJackknifeSurrogates{HurtadoModulationIndex}, X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}}=X) where {T<:Real} =
    (Array(UInt8, size(X, 1), size(X, 2)), Array(Int32, t.transform.nbins, nchannels(X)),
     Array(T, size(X, 1)), Array(Float64, t.transform.nbins), Array(Int32, t.transform.nbins), Array(Float64, t.transform.nbins))
function computestat!(t::AbstractJackknifeSurrogates{HurtadoModulationIndex}, out::JackknifeSurrogatesOutput,
                      work::Tuple{Matrix{UInt8}, Matrix{Int32}, Vector{T}, Vector{Float64}, Vector{Int32}, Vector{Float64}},
                      X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}}) where T<:Real
    trueval = out.trueval
    surrogates = out.surrogates
    nbins = t.transform.nbins
    phasebin, ninbin, amp, meanamp, ninbin2, meanamp2 = work
    chkinput(trueval, X, Y)
    size(ninbin, 1) == length(meanamp) == nbins || throw(ArgumentError("invalid work object"))
    ntrials(X) % jnn(t) == 0 || throw(DimensionMismatch("ntrials not evenly divisible by $(jnn(t))"))
    size(out.surrogates, 1) == div(ntrials(X), jnn(t)) || throw(DimensionMismatch("invalid output size"))

    binphases!(phasebin, ninbin, X)

    for k = 1:size(Y, 2)
        # Compute amplitudes
        for i = 1:size(X, 1)
            @inbounds amp[i] = abs(Y[i, k])
        end

        for j = 1:size(X, 2)
            sumphases!(meanamp, phasebin, amp, j)

            # Compute jackknifes
            for i = 1:size(surrogates, 1)
                for m = 1:nbins
                    @inbounds meanamp2[m] = meanamp[m]
                    @inbounds ninbin2[m] = ninbin[m, j]
                end
                for idel = (i-1)*jnn(t)+1:i*jnn(t)
                    @inbounds meanamp2[phasebin[idel, j]] -= amp[idel]
                    @inbounds ninbin2[phasebin[idel, j]] -= 1
                end
                surrogates[i, j, k] = finish_mi!(meanamp2, ninbin2, 1, t.transform.hmax)
            end

            # Compute true modulation index
            trueval[j, k] = finish_mi!(meanamp, ninbin, j, t.transform.hmax)
        end
    end

    out
end


text:
module PrepareWindow

export displayAll,createAndInitShaderProgram

using ModernGL, GeometryTypes, GLFW
using  ..PrepareWindowHelpers
using   ..OpenGLDisplayUtils,  ..DataStructs,Logging
using  ..ShadersAndVerticies,  ..ForDisplayStructs, ..ShadersAndVerticiesForText,  ..ModernGlUtil

"""
preparing all for displaying the images and responding to mouse and keyboard input
	listOfTexturesToCreate- list of texture specifications needed to for example create optimal shader
	calcDimsStruct - holds important data about verticies, textures dimensions etc.
"""
function displayAll(listOfTexturesToCreate::Vector{TextureSpec}
					,calcDimsStruct::CalcDimsStruct)
	# atomic variable that is enabling stopping async loop of event listening in order to enable othe actions with GLFW context
	stopListening = Threads.Atomic{Bool}(0)
	stopListening[]=false

	if(Threads.nthreads()==1)
		println("increase number of available threads look into https://docs.julialang.org/en/v1/manual/multi-threading/  or modify for example in vs code extension")
    end
    # Create the window. This sets all the hints and makes the context current.
	window = initializeWindow(calcDimsStruct.windowWidth,calcDimsStruct.windowHeight)

   	# The shaders
	println(createcontextinfo())
	gslsStr = get_glsl_version_string()

	vertex_shader = createVertexShader(gslsStr)

	masks = filter(textSpec-> !textSpec.isMainImage ,listOfTexturesToCreate)
	someExampleMask = masks[1]
	someExampleMaskB = masks[2]
	@info "masks set for subtraction $(someExampleMask.name)" someExampleMaskB.name
	fragment_shader_main,shader_program= createAndInitShaderProgram(vertex_shader,listOfTexturesToCreate,someExampleMask,someExampleMaskB,gslsStr  )

	glUseProgram(shader_program)

##for control of text display
	fragment_shader_words = ShadersAndVerticiesForText.createFragmentShader(gslsStr)
	shader_program_words = glCreateProgram()
	glAttachShader(shader_program_words, fragment_shader_words)
	glAttachShader(shader_program_words, vertex_shader)

	vbo_words = Ref(GLuint(1))   # initial value is irrelevant, just allocate space
    glGenBuffers(1, vbo_words)
##for control of text display

###########buffers
	#create vertex buffer
	createVertexBuffer()
	# Create the Vertex Buffer Objects (VBO)
	vbo = createDAtaBuffer(calcDimsStruct.mainImageQuadVert)

	# Create the Element Buffer Object (EBO)
	ebo = createElementBuffer( ShadersAndVerticies.elements)
	############ how data should be read from data buffer
	encodeDataFromDataBuffer()
	#capturing The data from GLFW
	controllWindowInput(window)

#loop that enables reacting to mouse and keyboards inputs  so every 0.1 seconds it will check GLFW weather any new events happened
	t = @task begin;
		while(!GLFW.WindowShouldClose(window))
		sleep(0.005);
		if(!stopListening[])
		 # Poll for and process events
		  GLFW.PollEvents()
			end
		end
		end
	schedule(t)

return (window,vertex_shader,fragment_shader_main ,shader_program,stopListening,vbo,ebo,fragment_shader_words,vbo_words,shader_program_words,gslsStr)

end# displayAll

"""
On the basis of information from listOfTexturesToCreate it creates specialized shader program
"""
function createAndInitShaderProgram(vertex_shader::UInt32
	,listOfTexturesToCreate::Vector{TextureSpec}
	,maskToSubtractFrom::TextureSpec
	,maskWeAreSubtracting ::TextureSpec
	,gslsStr::String)::Tuple{UInt32, UInt32}
			fragment_shader = createFragmentShader(gslsStr,listOfTexturesToCreate,maskToSubtractFrom,maskWeAreSubtracting)
			shader_program = glCreateProgram()
			glAttachShader(shader_program, fragment_shader)
			glAttachShader(shader_program, vertex_shader)
			glLinkProgram(shader_program)
return fragment_shader,shader_program

end#createShaderProgram

end #PreperWindow


text:
function get_arg_1(args...)
  if length(args) > 1
    return args[1]
  else
    return false
  end
end

println(get_arg_1(1,2,3))
println(get_arg_1([21,23]...))
println(get_arg_1([3:10;]...))


text:
using Surrogates

#1D
obj = x -> sin(x) + sin(x)^2 + sin(x)^3
lb = 0.0
ub = 10.0
x = sample(5,lb,ub,LowDiscrepancySample(2))
y = obj.(x)
p = 3.5
InverseDistance = InverseDistanceSurrogate(x,y,p,lb,ub)
prediction = InverseDistance(5.0)
add_point!(InverseDistance,5.0,-0.91)
add_point!(InverseDistance,[5.1,5.2],[1.0,2.0])

#ND
lb = [0.0,0.0]
ub = [10.0,10.0]
n = 100
x = sample(n,lb,ub,SobolSample())
f = x -> x[1]*x[2]^2
y = f.(x)
p = 3.0
InverseDistance = InverseDistanceSurrogate(x,y,p,lb,ub)
prediction = InverseDistance((1.0,2.0))
add_point!(InverseDistance,(5.0,3.4),-0.91)
add_point!(InverseDistance,[(5.1,5.2),(5.3,6.7)],[1.0,2.0])


text:
using Flux, MLDataPattern
using Flux: onehotbatch
include("spectrograms.jl")

# Image size constants
img_n = 150
img_m = 100

# Definition of a dataset with training and testing partitions
struct DataSet
  X_train::Array{Float32, 2}
  Y_train::Flux.OneHotMatrix
  X_test::Array{Float32, 2}
  Y_test::Flux.OneHotMatrix
end

# Split the data
function create_data_splits(n=img_n, m=img_m)
  @info("Generating spectrograms of size $(n)x$(m)...")

  # Read labels from source CSV file
  file_to_class = gen_label_dict("train/train.csv")
  labels = simple_labels(file_to_class)

  # Create spectrograms with values between 0 and 1 (via sigmoid)
  spec_base = create_spectrogram_base(file_to_class, n, m)

  # Format dataset
  X = hcat(float.(reshape.(spec_base, :))...) |> gpu
  XY = vcat(X, labels')

  get_label(x) = convert(Int, x[end])

  @info("Generating train and test splits...")

  # Stratified split: 80% for training and 20% for testing
  # Encode Ys with one-hot encoding
  train_data, test_data = stratifiedobs(get_label, XY, 0.8)
  X_train = train_data[1:end-1,:]
  Y_train = onehotbatch(convert(Array{Int}, train_data[end,:]), 0:9)
  X_test = test_data[1:end-1,:]
  Y_test = onehotbatch(convert(Array{Int}, test_data[end,:]), 0:9)

  @info("Done")

  return DataSet(X_train, Y_train, X_test, Y_test)
end


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

@doc raw"""DEPRECATED 1.9 - This group version of IPBlock is deprecated by networking/v1/IPBlock. IPBlock describes a particular CIDR (Ex. \"192.168.1.1/24\") that is allowed to the pods matched by a NetworkPolicySpec's podSelector. The except entry describes CIDRs that should not be included within this rule.

    IoK8sApiExtensionsV1beta1IPBlock(;
        cidr=nothing,
        except=nothing,
    )

    - cidr::String : CIDR is a string representing the IP Block Valid examples are \"192.168.1.1/24\"
    - except::Vector{String} : Except is a slice of CIDRs that should not be included within an IP Block Valid examples are \"192.168.1.1/24\" Except values will be rejected if they are outside the CIDR range
"""
mutable struct IoK8sApiExtensionsV1beta1IPBlock <: SwaggerModel
    cidr::Any # spec type: Union{ Nothing, String } # spec name: cidr
    except::Any # spec type: Union{ Nothing, Vector{String} } # spec name: except

    function IoK8sApiExtensionsV1beta1IPBlock(;cidr=nothing, except=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1IPBlock, Symbol("cidr"), cidr)
        setfield!(o, Symbol("cidr"), cidr)
        validate_property(IoK8sApiExtensionsV1beta1IPBlock, Symbol("except"), except)
        setfield!(o, Symbol("except"), except)
        o
    end
end # type IoK8sApiExtensionsV1beta1IPBlock

const _property_map_IoK8sApiExtensionsV1beta1IPBlock = Dict{Symbol,Symbol}(Symbol("cidr")=>Symbol("cidr"), Symbol("except")=>Symbol("except"))
const _property_types_IoK8sApiExtensionsV1beta1IPBlock = Dict{Symbol,String}(Symbol("cidr")=>"String", Symbol("except")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1IPBlock }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1IPBlock))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1IPBlock }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1IPBlock[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1IPBlock }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1IPBlock[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1IPBlock)
    (getproperty(o, Symbol("cidr")) === nothing) && (return false)
    true
end

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


text:
module AllenCahn

using SparseArrays
using LinearAlgebra

include("structures.jl")
export AllenCahnProblem1D, NeumannBC, PeriodicBC, BackwardEulerMethod, CrankNicolsonMethod

include("constructors.jl")
export AllenCahnProblem1D

include("assembly.jl")
export assemble_system!

include("solve.jl")
export solve

end #module


text:
using StructuredMatrices: binomial2

### Data layout
"""
Intention is for default data layout of CorrCholesky to be
Diagonals [aligned]
Lower sub-diagonal triangle (packed column major)
log(Diagonals) [aligned]
inv(Diagonals) [aligned] -- OPTIONAL, default absent
"""
mutable struct CorrCholesky{M,T,L} <: StructuredMatrices.AbstractMutableLowerTriangularMatrix{M,T,L}
    data::NTuple{L,T}
    CorrCholesky{M,T,L}(::UndefInitializer) where {M,T,L} = new{M,T,L}()
end
struct PtrCorrCholesky{M,T,L} <: StructuredMatrices.AbstractMutableLowerTriangularMatrix{M,T,L}
    ptr::Ptr{T}
end
"""
Intention is for default data layout of CorrCholesky to be
Diagonals [aligned]
Lower sub-diagonal triangle (packed column major)
log(Diagonals) [aligned]
inv(Diagonals) [aligned]
"""
mutable struct CovarCholesky{M,T,L} <: StructuredMatrices.AbstractMutableLowerTriangularMatrix{M,T,L}
    data::NTuple{L,T}
    CovarCholesky{M,T,L}(::UndefInitializer) where {M,T,L} = new{M,T,L}()
end
struct PtrCovarCholesky{M,T,L} <: StructuredMatrices.AbstractMutableLowerTriangularMatrix{M,T,L}
    ptr::Ptr{T}
end

const AbstractCorrCholesky{M,T,L} = Union{CorrCholesky{M,T,L},PtrCorrCholesky{M,T,L}}
const AbstractCovarCholesky{M,T,L} = Union{CovarCholesky{M,T,L},PtrCovarCholesky{M,T,L}}

@generated function PtrCovarCholesky{M,T}(sp::StackPointer) where {M,T}
    L = VectorizationBase.align(VectorizationBase.align(VectorizationBase.align(binomial2(M+1), T) + M, T) + M, T)
    quote
        $(Expr(:meta,:inline))
        sp + $(L*sizeof(T)), PtrCovarCholesky{$M,$T,$L}(pointer(sp, T))
    end
end

@generated function ReverseDiffExpressionsBase.alloc_adjoint(C::Union{<:AbstractCorrCholesky{M,T},<:AbstractCovarCholesky{M,T}}) where {M,T}
    L = PaddedMatrices.calc_padding(StructuredMatrices.binomial2(M+1), T)
    quote
        $(Expr(:meta,:inline))
        StructuredMatrices.PtrLowerTriangularMatrix{$M,$T,$L}(SIMDPirates.alloca(Val{$L}(),$T))
    end
end
@inline ReverseDiffExpressionsBase.alloc_adjoint(sp::StackPointer, C::AbstractCorrCholesky{M,T}) where {M,T} = StructuredMatrices.PtrLowerTriangularMatrix{M,T}(sp)
@inline ReverseDiffExpressionsBase.alloc_adjoint(sp::StackPointer, C::AbstractCovarCholesky{M,T}) where {M,T} = StructuredMatrices.PtrLowerTriangularMatrix{M,T}(sp)

function caches_logdiag(M, L, ::Type{T} = Float64) where {T}
    offset_length = VectorizationBase.align(binomial2(M+1),T)
    L >= offset_length + M, offset_length
end
function caches_invdiag(M, L, ::Type{T} = Float64) where {T}
    offset_length = VectorizationBase.align(VectorizationBase.align(binomial2(M+1),T) + M, T)
    L >= offset_length + M, offset_length
end
@inline logdiag(C::StructuredMatrices.AbstractDiagTriangularMatrix{M,T,L}) where {M,T,L} = PaddedMatrices.LazyMap(SLEEFPirates.log, PtrVector{M,T,M,true}(pointer(C)))
@inline invdiag(C::StructuredMatrices.AbstractDiagTriangularMatrix{M,T,L}) where {M,T,L} = PaddedMatrices.LazyMap(SIMDPirates.vinv, PtrVector{M,T,M,true}(pointer(C)))
@generated function logdiag(C::Union{<:AbstractCorrCholesky{M,T,L},<:AbstractCovarCholesky{M,T,L}}) where {M,T,L}
    c, offset = caches_logdiag(M,L)
    if c
        quote
            $(Expr(:meta,:inline))
            PtrVector{$M,$T,$(PaddedMatrices.calc_padding(M,T)),false}(pointer(C) + $(offset*sizeof(T)))
        end
    else
        quote
            $(Expr(:meta,:inline))
            d = PtrVector{$M,$T,$M,true}(pointer(C))
            PaddedMatrices.LazyMap(SLEEFPirates.log, d)
        end
    end
end
@generated function invdiag(C::Union{<:AbstractCorrCholesky{M,T,L},<:AbstractCovarCholesky{M,T,L}}) where {M,T,L}
    c, offset = caches_invdiag(M, L)
    if c
        quote
            $(Expr(:meta,:inline))
            PtrVector{$M,$T,$(PaddedMatrices.calc_padding(M,T)),false}(pointer(C) + $(offset*sizeof(T)))
        end
    else
        quote
            $(Expr(:meta,:inline))
            d = PtrVector{$M,$T,$M,true}(pointer(C))
            PaddedMatrices.LazyMap(SIMDPirates.vinv, d)
        end
    end
end

@generated PaddedMatrices.type_length(::Type{<: AbstractCorrCholesky{M}}) where {M} = StructuredMatrices.binomial2(M+1)
@generated PaddedMatrices.type_length(::AbstractCorrCholesky{M}) where {M} = StructuredMatrices.binomial2(M+1)
@generated PaddedMatrices.param_type_length(::Type{<: AbstractCorrCholesky{M}}) where {M} = StructuredMatrices.binomial2(M)
@generated PaddedMatrices.param_type_length(::AbstractCorrCholesky{M}) where {M} = StructuredMatrices.binomial2(M)

@inline Base.pointer(L::Union{CorrCholesky{M,T},CovarCholesky{M,T}}) where {M,T} = Base.unsafe_convert(Ptr{T}, pointer_from_objref(L))
@inline Base.pointer(L::Union{PtrCorrCholesky,PtrCovarCholesky}) = L.ptr
@inline Base.unsafe_convert(::Type{Ptr{T}}, L::PtrCorrCholesky{M,T}) where {M,T} = L.ptr

abstract type AbstractCholeskyConstraintAdjoint{P,T,L} <: AbstractArray{T,4} end
mutable struct CholeskyConstraintAdjoint{P,T,L} <: AbstractCholeskyConstraintAdjoint{P,T,L}
    data::NTuple{L,T}
    CholeskyConstraintAdjoint{P,T,L}(::UndefInitializer) where {P,T,L} = new()
end
struct PtrCholeskyConstraintAdjoint{P,T,L} <: AbstractCholeskyConstraintAdjoint{P,T,L}
    ptr::Ptr{T}
end
@inline VectorizationBase.vectorizable(A::Union{PtrCorrCholesky,PtrCovarCholesky}) = VectorizationBase.Pointer(A.ptr)
@inline Base.pointer(A::PtrCholeskyConstraintAdjoint) = A.ptr
@inline Base.unsafe_convert(::Type{Ptr{T}}, A::PtrCholeskyConstraintAdjoint{P,T}) where {P,T} = A.ptr
Base.size(::AbstractCholeskyConstraintAdjoint{P}) where {P} = (P,P,P-1,P-1)

@inline function Base.getindex(adj::CholeskyConstraintAdjoint{P,T,L}, i::Integer) where {P,T,L}
    @boundscheck i > L && PaddedMatrices.ThrowBoundsError("i = $i > $L.")
    @inbounds adj.data[i]
end
@inline function Base.getindex(A::Union{PtrCorrCholesky{P,T,L},PtrCholeskyConstraintAdjoint{P,T,L}}, i::Integer) where {P,T,L}
    @boundscheck i > L && PaddedMatrices.ThrowBoundsError("i = $i > $L.")
    VectorizationBase.load(A.ptr + (i-1)*sizeof(T))
end
@inline function Base.setindex!(A::Union{PtrCorrCholesky{P,T,L},PtrCholeskyConstraintAdjoint{P,T,L}}, v::T, i::Integer) where {P,T,L}
    @boundscheck i > L && PaddedMatrices.ThrowBoundsError("i = $i > $L.")
    VectorizationBase.store!(A.ptr + (i-1)*sizeof(T), v)
end
@inline function Base.getindex(adj::AbstractCholeskyConstraintAdjoint{P,T,L}, i::CartesianIndex{4}) where {P,T,L}
    adj[i[1],i[2],i[3],i[4]]
end
function Base.getindex(adj::AbstractCholeskyConstraintAdjoint{Mp1,T,L}, i, j, k, l) where {Mp1,T,L}
    # Welcome to Branchapalooza!!!
    kp1 = k + 1
    M = Mp1 - 1
    @boundscheck max(i,j,kp1,l+1) > Mp1 && PaddedMatrices.ThrowBoundsError("At least one of (i,j,k,l) = ($i,$j,$k,$l) > ($Mp1,$Mp1,$M,$M)")
    # Lets try and skip trudging through the loop
    j > i && return zero(T)
    l > k && return zero(T)
    i == kp1 || return zero(T)
    i == j == 1 && return zero(T)
    # note order matters.
    # if i == 1, return zero(T), or if
    # if i != k + 1, return zero(T)
    l == j == 1 && return one(T)# + 10l + 100j
    if i == j == kp1
        ind = 0
        # This for loop can be removed
        # I'll get around to it
        # but this function is not performance critical.
        for mp ∈ 1:M
            if mp == l
                return @inbounds adj[ind + i - mp]
            else
                ind += Mp1 - mp
            end
        end
    end
    ind = StructuredMatrices.binomial2(Mp1)
    for mc ∈ 2:M # columns of output ()
        mc_equal_j = mc == j
        for mp ∈ 1:mc #
            mp_equal_l = mp == l
            for mr ∈ mc+1:Mp1 # iter through columns of z
                ind += 1
                if mc_equal_j && mp_equal_l && (mr == i)
                    return @inbounds adj[ind]
                end
            end
        end
    end
    zero(T)
end

function corr_cholesky_adjoint_mul_quote(Mp1, T, add::Bool)#,sp::Bool = false)
    M = Mp1 - 1
    q = quote end
    StructuredMatrices.load_packed_L_quote!(q.args, Mp1, :∂target_∂L, :t)
    # first, we initialize all outputs with the diagonal pass
    # for the partials with respect to column 1 of the inputs,
    # we also add the identity matrix part of the jacobian
    # (ie, partials with respect to the first column are "1")
    for mc ∈ 1:M
        push!(q.args, :(
            $(Symbol(:∂lkj_∂z_, mc, :_, 1)) = $(PaddedMatrices.sym(:∂target_∂L, mc+1, 1)) +
                                        $(PaddedMatrices.sym(:∂target_∂L, mc+1, mc+1)) * adj[$mc]
                        )
        )
    end
    # partials with respec to remaining columns of the inputs
    ind = M
    for mp ∈ 2:M
        for mc ∈ mp:M
            ind += 1
            push!(q.args, :($(Symbol(:∂lkj_∂z_, mc, :_, mp)) = $(PaddedMatrices.sym(:∂target_∂L, mc+1, mc+1)) * adj[$ind] ))
        end
    end
    # now that all partials are initialized, we update them.
    for mc ∈ 2:M
        for mp ∈ 1:mc
            for mr ∈ mc+1:Mp1 # iter through columns of z
                ind += 1
                push!(q.args, :($(Symbol(:∂lkj_∂z_, mr-1, :_, mp)) += $(PaddedMatrices.sym(:∂target_∂L, mr, mc)) * adj[$ind] ))
            end
        end
    end
    lkj_l = StructuredMatrices.binomial2(Mp1)
    lkj_l_full = PaddedMatrices.calc_padding(lkj_l, T)
    push!(q.args, :(sptr = pointer(out)))
    ind = 0
    size_T = sizeof(T)
    for pc ∈ 1:M
        for pr ∈ pc:M
            sptrind = :(sptr + $size_T*$ind)
            ∂lkjsym = PaddedMatrices.sym(:∂lkj_∂z, pr, pc)
            if add
                push!(q.args, :(VectorizationBase.store!($sptrind, Base.FastMath.add_fast(VectorizationBase.load($sptrind), $∂lkjsym))))
            else
                push!(q.args, :(VectorizationBase.store!($sptrind, $∂lkjsym)))
            end
            ind += 1
        end
    end
    push!(q.args, :out)
    quote
        @fastmath @inbounds begin
            $q
        end
    end
end

# @generated function LinearAlgebra.mul!(
#     c::AbstractLowerTriangularMatrix{Mp1,T},
#     t::LinearAlgebra.Adjoint{T,<:PaddedMatrices.AbstractFixedSizeVector{BP,T,BPL}},
#     adj::AbstractCholeskyConstraintAdjoint{Mp1,T,L}
# ) where {Mp1,T,L,BP,BPL}
#     corr_cholesky_adjoint_mul_quote(Mp1,T)
# end

@generated function LinearAlgebra.mul!(
    out::StructuredMatrices.AbstractLowerTriangularMatrix{M,T},
    t::StructuredMatrices.AbstractLowerTriangularMatrix{Mp1,T},
    adj::AbstractCholeskyConstraintAdjoint{Mp1,T,L}
) where {Mp1,T,L,M}
    corr_cholesky_adjoint_mul_quote(Mp1,T,false)
end
@generated function PaddedMatrices.muladd!(
    out::StructuredMatrices.AbstractLowerTriangularMatrix{M,T},
    t::StructuredMatrices.AbstractLowerTriangularMatrix{Mp1,T},
    adj::AbstractCholeskyConstraintAdjoint{Mp1,T,L}
) where {Mp1,T,L,M}
    corr_cholesky_adjoint_mul_quote(Mp1,T,true)
end

function add_calc_logdet_expr!(q, M, ::Type{T}, logdetsym) where {T}
    Mp1 = M + 1; Mm1 = M - 1
    if M < 2
        push!(q.args, :($logdetsym = zero(T)))
    elseif Mp1 < 10
        push!(q.args, :($logdetsym = log( $(Expr(:call, :(Base.FastMath.mul_fast), [Symbol(:ljp_,m) for m ∈ 3:Mp1]...)) ) ))
    else
        W, Wshift = VectorizationBase.pick_vector_width_shift(Mm1, T)
        ljp_reps = Mm1 >>> Wshift
        ljp_rem = Mm1 & (W - 1)
        mincr = 2
        push!(q.args, :($logdetsym = SLEEFPirates.log($(Expr(:tuple, [:(Core.VecElement{$T}($(Symbol(:ljp_,w+mincr)))) for w ∈ 1:W]...)))))
        for mr ∈ 1:ljp_reps
            mincr += W
            push!(q.args, :($logdetsym = vadd($logdetsym, SLEEFPirates.log($(Expr(:tuple, [:(Core.VecElement{$T}($(Symbol(:ljp_,w+mincr)))) for w ∈ 1:W]...))))))
        end
        if ljp_rem > 0
            mincr += W
            push!(q.args, :($logdetsym = vadd($logdetsym, SLEEFPirates.log($(Expr(:tuple, [:(Core.VecElement{$T}($(w > Mp1 ? one(T) : Symbol(:ljp_,w)))) for w ∈ 1+mincr:W+mincr]...))))))
        end
    end
    nothing
end

function store_symbol_vector!(q::Expr, sv::Vector, out::Symbol, ::Type{T}) where {T}
    num_syms = length(sv)
    W, Wshift = VectorizationBase.pick_vector_width_shift(num_syms, T)
    i = 0
    ptr_out = gensym(Symbol(:ptr_, out))
    push!(q.args, Expr(:(=), ptr_out, Expr(:call, :pointer, out)))
    for _ ∈ 1:(num_syms >>> Wshift)
        push!(q.args, :(vstore!($ptr_out, $(Expr(:tuple, [:(Core.VecElement($(sv[i+w]))) for w ∈ 1:W]...)))))
        push!(q.args, :($ptr_out += $(sizeof(T)*W)))
        i += W
    end
    rem = num_syms & (W - 1)
    if rem > 0
        push!(q.args, :(vstore!($ptr_out, $(Expr(:tuple, [w > rem ? Core.VecElement(zero(T)) : :(Core.VecElement($(sv[i+w]))) for w ∈ 1:W]...)), $(VectorizationBase.mask(T, rem)))))
    end
    nothing
end

"""
Generates the quote for the constraining transformation from z ∈ (-1,1) to _Correlation_Cholesky
without taking the derivative of the expression.
"""
function constrain_lkj_factor_quote(
    Mp1::Int, L::Int, LL::Int, ::Type{T}, lkjsym::Symbol, zsym::Symbol
) where {T}
    # @show L
    M = Mp1 - 1
    # M = (Int(sqrt(1 + 8L))-1)>>>1
    # Mp1 = M+1
    q = quote x_1_1 = one($T) end
    for m ∈ 1:M
        xm = Symbol(:x_, m+1, :_, 1)
        push!(q.args, :( $xm = $zsym[$m]) )
        push!(q.args, :($(Symbol(:Omx²_, m+1, :_, 2)) = 1 - $xm * $xm))
    end
    ljnum = 0
    zind = M
    for mc ∈ 2:Mp1
        push!(q.args, :($(Symbol(:x_, mc, :_, mc)) = sqrt($(Symbol(:Omx²_, mc, :_, mc)))))
        for mr ∈ mc+1:Mp1
            xm = Symbol(:x_, mr, :_, mc)
            Omx²m = Symbol(:Omx²_, mr, :_, mc)
            ljnum += 1
            zind += 1
            lj = Symbol(:lj_, mr, :_, mc)
            # We're taking the derivative with respect to 0.5log(Omx²m) for each of the Omx²m encountered in this loop.
            push!(q.args, :($lj = sqrt($Omx²m) ))
            push!(q.args, :($xm = $zsym[$zind] * $lj ))
            push!(q.args, :($(Symbol(:Omx²_, mr, :_, mc+1)) = $Omx²m - $xm * $xm ))
            if mc == 2
                push!(q.args, :( $(Symbol(:ljp_, mr)) = $lj))
            else
                push!(q.args, :( $(Symbol(:ljp_, mr)) *= $lj))
            end
        end
    end
    logdetsym = gensym(:logdet)
    add_calc_logdet_expr!(q, M, T, logdetsym)
    lkjsyms = Symbol[]
    # i = 0
    for mc ∈ 1:M+1
        # i += 1
        # push!(q.args, :($lkjsym[$i] = $(Symbol(:x_, mc, :_, mc))))
        push!(lkjsyms, Symbol(:x_, mc, :_, mc))
    end
    for mc ∈ 1:M+1
        for mr ∈ mc+1:M+1
            # i += 1
            # push!(q.args, :($lkjsym[$i] = $(Symbol(:x_, mr, :_, mc))))
            push!(lkjsyms, Symbol(:x_, mr, :_, mc))
        end
    end
    store_symbol_vector!(q, lkjsyms, lkjsym, T)
    if first(caches_logdiag(Mp1, LL, T)) # store log of diagonals
        lkj_length_triangle = VectorizationBase.align(StructuredMatrices.binomial2(Mp1+1), T)
        msym = gensym(:m)
        W = VectorizationBase.pick_vector_width(M, T)
        logdiag_quote = if Mp1 % W == 1
            quote
                $lkjsym[$(lkj_length_triangle + 1)] = zero($T)
                LoopVectorization.@vvectorize_unsafe $T for $msym ∈ 1:$M
                    $lkjsym[$msym+$(lkj_length_triangle+1)] = SLEEFPirates.log($lkjsym[$msym+1])
                end
            end
        else
            quote
                LoopVectorization.@vvectorize_unsafe $T for $msym ∈ 1:$(PaddedMatrices.calc_padding(Mp1,T))
                    $lkjsym[$msym+$lkj_length_triangle] = SLEEFPirates.log($lkjsym[$msym])
                end
            end
        end
        push!(q.args, macroexpand(LoopVectorization, logdiag_quote))
    end
    quote
        # $(Expr(:meta,:inline))
        @fastmath @inbounds begin
            $q
        end
        $logdetsym
    end
end

function lkj_adjoint_length(M)
    Ladj = 0
    for mp ∈ 1:M, mc ∈ mp+1:M+1
        Ladj += 1
    end
    for mc ∈ 2:M, mp ∈ 1:mc, mr ∈ mc+1:M+1
        Ladj += 1
    end
    Ladj
end

∂lkj_sym(i, j, k) = Symbol(:∂x_, i, :_, j, :_∂z_, k)

"""
This function also generates code to calculate the gradient with respect to the log determinant of the Jacobian,
as well as the full Jacobian of the constraining transformation.
"""
function constrain_lkj_factor_jac_quote(
    Mp1::Int, L::Int, LL::Int, ::Type{T}, lkjsym::Symbol, ∂logdetsym::Symbol, jacobiansym::Symbol, zsym::Symbol
) where {T}
    # @show L
    M = Mp1 - 1
    # partial = true
    # Storage is the subdiagonal lower triangle
    # (the diagonal itself is determined by the correlation matrix constraint)
    #
    # The correlation matrix is (M+1) x (M+1)
    # our subdiagonal triangle is M x M
    #
    # We denote the (-1, 1) constrained parameters "z"
    # while "x" is the (M+1) x (M+1) cholesky factor of the correlation matrix.
    q = quote x_1_1 = one($T) end
    for m ∈ 1:M
        xm = Symbol(:x_, m+1, :_, 1)
        Omx²m = Symbol(:Omx²_, m+1, :_2)
        push!(q.args, :( $xm = $zsym[$m]) )
        push!(q.args, :( $Omx²m = one($T) - $xm * $xm ) )
        push!(q.args, :($(Symbol(:∂x_, m+1, :_1_∂_1)) = one($T) ) )
    end
    ljnum = 0
    zind = M
    for mc ∈ 2:Mp1
        xm = Symbol(:x_, mc, :_, mc)
        # Note that the partials terminate at 1 less, because "z" only fills the
        # lower subdiagonal triangle. Iteration stops 2 early, as these are ones that reference
        # previous calculations. Final is the first partial for that z-term
        # (that z term would be the furthest to the right in its given row)
        xm_prev = Symbol(:x_, mc, :_, mc-1) # update based on that z term
        Omx²m = Symbol(:Omx²_, mc, :_, mc); Omx²_previous = Symbol(:Omx²_, mc, :_, mc-1)
        mc > 2 && push!(q.args, :($Omx²m = $Omx²_previous - $xm_prev * $xm_prev))
        push!(q.args, :($xm = sqrt($Omx²m)))
        for mc_nested ∈ 1:mc-2
            ∂Omx²_previous = Symbol(:∂, Omx²_previous, :_, mc_nested)
            push!(q.args, :($(Symbol(:∂, Omx²m, :_, mc_nested)) = $∂Omx²_previous -2*$xm_prev * $(Symbol(:∂, xm_prev, :_∂_, mc_nested))   ))
        end
        let mc_nested = mc - 1 # initialize
            push!(q.args, :($(Symbol(:∂, Omx²m, :_, mc_nested)) = -2*$xm_prev * $(Symbol(:∂, xm_prev, :_∂_, mc_nested)) ))
        end
        # Now determin ∂x/∂z
        for mc_nested ∈ 1:mc-1 # partial is a function of previous values
            ∂xm = Symbol(:∂, xm, :_∂_, mc_nested)
            ∂Omx²m = Symbol(:∂, Omx²m, :_, mc_nested)
            push!(q.args, :($∂xm = $(T(0.5)) * $xm * $∂Omx²m / $Omx²m ) )
        end
        for mr ∈ mc+1:Mp1
            xm = Symbol(:x_, mr, :_, mc)
            xm_prev = Symbol(:x_, mr, :_, mc-1)
            Omx²_previous = Symbol(:Omx²_, mr, :_, mc-1)
            Omx²m = Symbol(:Omx²_, mr, :_, mc)
            ljnum += 1
            zind += 1
            lj = Symbol(:lj_, mr, :_, mc)
            mc > 2 && push!(q.args, :($Omx²m = $Omx²_previous - $xm_prev * $xm_prev ))
            # We're taking the derivative with respect to 0.5log(Omx²m) for each of the Omx²m encountered in this loop.
            push!(q.args, :($lj = sqrt($Omx²m) ))
            push!(q.args, :($xm = $zsym[$zind] * $lj ))
            if mc == 2
                push!(q.args, :( $(Symbol(:ljp_, mr)) = $lj))
            else
                push!(q.args, :( $(Symbol(:ljp_, mr)) *= $lj))
            end
            # Now we calculate partials with respect to Omx²m from this collomn
            for mc_nested ∈ 1:mc-2
                ∂Omx²_previous = Symbol(:∂, Omx²_previous, :_, mc_nested)
                push!(q.args, :($(Symbol(:∂, Omx²m, :_, mc_nested)) = $∂Omx²_previous -2*$xm_prev * $(Symbol(:∂, xm_prev, :_∂_, mc_nested))   ))
                push!(q.args, :($(Symbol(:∂ljp_, mr, :_, mc_nested)) += $(T(0.5)) * $(Symbol(:∂, Omx²m, :_, mc_nested)) / $Omx²m ))
            end
            let mc_nested = mc - 1 # initialize
                push!(q.args, :($(Symbol(:∂, Omx²m, :_, mc_nested)) = -2*$xm_prev * $(Symbol(:∂, xm_prev, :_∂_, mc_nested)) ))
                push!(q.args, :($(Symbol(:∂ljp_, mr, :_, mc_nested)) = $(T(0.5)) * $(Symbol(:∂, Omx²m, :_, mc_nested)) / $Omx²m ))
            end
            for mc_nested ∈ 1:mc-1 # partial is a function of previous values
                ∂xm = Symbol(:∂, xm, :_∂_, mc_nested)
                ∂Omx²m = Symbol(:∂, Omx²m, :_, mc_nested)
                push!(q.args, :($∂xm = $(T(0.5)) * $xm * $∂Omx²m / $Omx²m ) )
            end
            let mc_nested = mc # new partial
                ∂xm = Symbol(:∂, xm, :_∂_, mc_nested)
                push!(q.args, :($∂xm = $lj ) )
            end
        end
    end
    logdetsym = gensym(:logdet)
    add_calc_logdet_expr!(q, M, T, logdetsym)
    lkj_length_triangle = VectorizationBase.align(binomial2(Mp1+1),T)
    #lkj_length = lkj_length_triangle + VectorizationBase.align(Mp1, T)
    # lkjsym = gensym(:CorrCholesky)
    # i = 0
    lkjsyms = Symbol[]
    for mc ∈ 1:Mp1
        # i += 1
        # push!(q.args, :($lkjsym[$i] = $(Symbol(:x_, mc, :_, mc))))
        push!(lkjsyms, Symbol(:x_, mc, :_, mc))
    end
    for mc ∈ 1:Mp1
        for mr ∈ mc+1:Mp1
            # i += 1
            # push!(q.args, :($lkjsym[$i] = $(Symbol(:x_, mr, :_, mc))))
            push!(lkjsyms, Symbol(:x_, mr, :_, mc))
        end
    end
    store_symbol_vector!(q, lkjsyms, lkjsym, T)
    W = VectorizationBase.pick_vector_width(M, T)
    msym = gensym(:m)
    logdiag_quote = if Mp1 % W == 1
        quote  # If the remainder is 1, we can save a loop iteration by simply assigning zero to the first element
            $lkjsym[$(1 + lkj_length_triangle)] = zero($T)
            LoopVectorization.@vvectorize_unsafe $T for $msym ∈ 1:$M
                $lkjsym[$msym + $(1+lkj_length_triangle)] = SLEEFPirates.log($lkjsym[$msym+1])
            end
        end
    else
        quote
            LoopVectorization.@vvectorize_unsafe $T for $msym ∈ 1:$(PaddedMatrices.calc_padding(Mp1, T))
                $lkjsym[$msym + $lkj_length_triangle] = SLEEFPirates.log($lkjsym[$msym])
            end
        end
    end
    push!(q.args, macroexpand(LoopVectorization, logdiag_quote))
    bin2M = binomial2(Mp1)
    ∂logdetsyms = Union{T,Symbol}[]
    # i = 1
    for mc ∈ 1:M
        push!(∂logdetsyms, zero(T))
              # push!(q.args, :($∂logdetsym[$i] = zero($T)))
        # i += 1
        for mr ∈ mc+2:M+1
            push!(∂logdetsyms, Symbol(:∂ljp_, mr, :_, mc))
            # push!(q.args, :($∂logdetsym[$i] = $(Symbol(:∂ljp_, mr, :_, mc))))
            # i += 1
        end
    end
    store_symbol_vector!(q, ∂logdetsyms, ∂logdetsym, T)
    # Ladj = DistributionParameters.lkj_adjoint_length(M)
    # diagonal block of lkj
    # i = 1
    jacsyms = Symbol[]
    for mp ∈ 1:M
        for mc ∈ mp+1:Mp1
            push!(jacsyms, Symbol(:∂x_, mc, :_, mc, :_∂_, mp))
            # push!(q.args, :($jacobiansym[$i] = $(Symbol(:∂x_, mc, :_, mc, :_∂_, mp))))
            # i +=1
        end
    end
    for mc ∈ 2:M
        for mp ∈ 1:mc
            for mr ∈ mc+1:Mp1 # iter through columns of z
                push!(jacsyms, Symbol(:∂x_, mr, :_, mc, :_∂_, mp))
                # push!(q.args, :($jacobiansym[$i] = $(Symbol(:∂x_, mr, :_, mc, :_∂_, mp))))
                # i += 1
            end
        end
    end
    store_symbol_vector!(q, jacsyms, jacobiansym, T)
# Return quote, followed by 4 outputs of the function
# four outputs are:
# 1. constrained
# 2. logdeterminant of constraining transformation
# 3. gradient with respect to log determinant
# 4. Jacobian of the constraining transformation
# constrain_q, constrained_expr, logdetsym, logdetgrad, jacobian
    quote
        # $(Expr(:meta,:inline))
        @fastmath @inbounds begin
            $q
        end
        $logdetsym
    end#, lkjsym, logdetsym, ∂logdetsym, jacobiansym
end

@generated function lkj_constrain!(
    lkjcorrchol::StructuredMatrices.AbstractMutableLowerTriangularMatrix{Mp1,T,LL},
    zlkj::PaddedMatrices.AbstractFixedSizeVector{L,T}
) where {T,L,Mp1,LL}
    constrain_lkj_factor_quote(Mp1, L, LL, T, :lkjcorrchol, :zlkj)
end
@generated function ∂lkj_constrain!(
    lkjcorrchol::StructuredMatrices.AbstractMutableLowerTriangularMatrix{Mp1,T,LL},
    lkjgrad::StructuredMatrices.AbstractMutableLowerTriangularMatrix{M,T},
    lkjjacobian::AbstractCholeskyConstraintAdjoint{Mp1,T,Ladj},
    zlkj::PaddedMatrices.AbstractFixedSizeVector{L,T}
) where {L,T,Mp1,M,LL,Ladj}
# ) where {L,T,Mp1,LL,Ladj,M}
    constrain_lkj_factor_jac_quote(Mp1, L, LL, T, :lkjcorrchol, :lkjgrad, :lkjjacobian, :zlkj)
end

# @generated function lkj_constrain(
#     sp::StackPointer,
#     zlkj::PaddedMatrices.AbstractFixedSizeVector{L,T},
#     ::Val{Align} = Val{true}()
# ) where {L,T,Align}
#     lkjconstrain_q, lkjconstrained_expr, lkjlogdetsym = constrain_lkj_factor_quote(L, T, :zlkj, true, Align)
#     quote
#         # Inlined because of:
#         # https://github.com/JuliaLang/julia/issues/32414
#         # Stop forcing inlining when the issue is fixed.
#         $(Expr(:meta,:inline))
#         $lkjconstrain_q
#         sp, ($lkjconstrained_expr, $lkjlogdetsym)
#     end
# end
# @generated function ∂lkj_constrain(
#     sp::StackPointer,
#     zlkj::PaddedMatrices.AbstractFixedSizeVector{L,T}
# ) where {L,T}
#     lkjconstrain_q, lkjconstrained_expr, lkjlogdetsym, lkjlogdetgrad, lkjjacobian = constrain_lkj_factor_jac_quote(L, T, :zlkj, true)
#     quote
#         # Inlined because of:
#         # https://github.com/JuliaLang/julia/issues/32414
#         # Stop forcing inlining when the issue is fixed.
#         $(Expr(:meta,:inline))
#         $lkjconstrain_q
#         sp, ($lkjconstrained_expr, $lkjlogdetsym, $lkjlogdetgrad, $lkjjacobian)
#     end
# end

function load_parameter!(
    first_pass::Vector{Any}, second_pass::Vector{Any}, out::Symbol,
    ::Type{<: AbstractCorrCholesky{M,T}}, partial::Bool,
    m::Module, sp::Union{Symbol,Nothing}, logjac::Bool = true, exportparam::Bool = false
) where {M,T}
    θ = Symbol("##θparameter##")
    ∂θ = Symbol("##∂θparameter##")
    ninvlogitout = gensym(out)
    invlogitout = gensym(out)
    ∂invlogitout = gensym(out)
    use_sptr = sp isa Symbol
    W, Wshift = VectorizationBase.pick_vector_width_shift(M, T)
    sumθᵢ = gensym(:sumθᵢ)
    N = (M * (M-1)) >>> 1
    Wm1 = W - 1
    rem = N & Wm1
    L = (N + Wm1) & ~Wm1
    zsym = gensym(:z) # z ∈ (-1, 1)
    zsymoffset = VectorizationBase.align(binomial2(M+1)*sizeof(T))
    if partial
        zsymoffset += VectorizationBase.align(max(2L,binomial2(M-1) + lkj_adjoint_length(M-1))*sizeof(T))
    end
    limitlife = use_sptr && !exportparam
    # zsym will be discarded, so we allocate it behind all the data we actually return.
    if use_sptr && !exportparam
        if partial
            push!(first_pass, :(($sp,$invlogitout) = $m.PtrVector{$L,$T}($sp)))
            push!(first_pass, :(($sp,$∂invlogitout) = $m.PtrVector{$L,$T}($sp)))
            limitlife && push!(first_pass, :($m.lifetime_start!($invlogitout); $m.lifetime_start!($∂invlogitout)))
        end
        # push!(first_pass, :( $zsym = $m.PtrVector{$N,$T}(pointer($sp + $zsymoffset,$T)) ))
        # push!(first_pass, :( ($sp,$zsym) = $m.PtrVector{$N,$T}($sp)))
        # limitlife && push!(first_pass, :($m.lifetime_start!($zsym)))
        push!(first_pass, :( $zsym = $m.PtrVector{$N,$T}($m.SIMDPirates.alloca(Val{$N}(),$T)) ))
    else
        push!(first_pass, :( $zsym = $m.PtrVector{$N,$T}($m.SIMDPirates.alloca(Val{$N}(),$T)) ))
        if partial
            push!(first_pass, :($invlogitout = $m.PtrVector{$L,$T}($m.alloca(Val{$L}(),$T))))
            push!(first_pass, :($∂invlogitout = $m.PtrVector{$L,$T}($m.alloca(Val{$L}(),$T))))
        end
    end
    i = gensym(:i)
    loop_body = quote $ninvlogitout = one($T) / (one($T) + $m.SLEEFPirates.exp($θ[$i])) end
    if partial
        push!(loop_body.args, :($invlogitout[$i] = one($T) - $ninvlogitout))
        push!(loop_body.args, :($∂invlogitout[$i] = $ninvlogitout * $invlogitout[$i]))
        push!(loop_body.args, :($zsym[$i] = $m.SIMDPirates.vmuladd($(T(-2)), $ninvlogitout, one($T))))
        logjac && push!(loop_body.args, :(target = $m.vadd(target, $m.SLEEFPirates.log($∂invlogitout[$i]))))
    else
        push!(loop_body.args, :($invlogitout = one($T) - $ninvlogitout))
        push!(loop_body.args, :($∂invlogitout = $ninvlogitout * $invlogitout))
        push!(loop_body.args, :($zsym[$i] = $m.SIMDPirates.vmuladd($(T(-2)), $ninvlogitout, one($T))))
        logjac && push!(loop_body.args, :(target = $m.vadd(target, $m.SLEEFPirates.log($∂invlogitout))))
    end
    vloop_quote = quote
        LoopVectorization.@vvectorize_unsafe $T $((m)) for $i ∈ 1:$N
            $loop_body
        end
    end
    push!(first_pass, macroexpand(m, vloop_quote))
    lkjlogdetsym = gensym(:lkjlogdetsym)
    if exportparam
        @assert partial == false
        lkj_length = binomial2(M+1)
    else
        lkj_length_triangle = VectorizationBase.align(binomial2(M+1),T)
        lkj_length = lkj_length_triangle + VectorizationBase.align(M, T)
    end
    if partial
        seedlkj = gensym(:seedlkj)
        lkjlogdetgradsym = gensym(:lkjlogdetgrad)
        lkjjacobiansym = gensym(:lkjjacobian)
        adjout = adj(out)
        push!(second_pass, :($m.muladd!($seedlkj, $adjout, $lkjjacobiansym)))
        loopq = quote
            LoopVectorization.@vvectorize_unsafe $T $((m)) for $i ∈ 1:$N
                $seedlkj[$i] = $(one(T)) - $(T(2)) * ( ($invlogitout)[$i] - (($seedlkj)[$i]) * ($∂invlogitout)[$i] )
            end
        end
        push!(second_pass, macroexpand(m, loopq))
        if limitlife
            for d ∈ (invlogitout, ∂invlogitout, adjout, out, lkjjacobiansym)
                push!(second_pass, :($m.lifetime_end!($d)))
            end
        end
        Ladj = DistributionParameters.lkj_adjoint_length(M - 1)
        if use_sptr #out
            push!(first_pass, :($lkjjacobiansym = $m.DistributionParameters.PtrCholeskyConstraintAdjoint{$M,$T,$Ladj}(pointer($sp,$T))))
            push!(first_pass, :($sp += $(VectorizationBase.align(Ladj*sizeof(T)))))
            push!(first_pass, :($out = $m.DistributionParameters.PtrCorrCholesky{$M,$T,$lkj_length}(pointer($sp,$T))))
            # push!(first_pass, Expr(:(=), out, :(PtrCorrCholesky{$M,$T,$lkj_length}(pointer($sp,$T)))))
            push!(first_pass, :($sp += $(VectorizationBase.align(sizeof(T)*lkj_length))))
        else
            push!(first_pass, :($lkjjacobiansym = $m.DistributionParameters.CholeskyConstraintAdjoint{$M,$T,$Ladj}(undef)))
            push!(first_pass, :($out = $m.DistributionParameters.CorrCholesky{$M,$T,$lkj_length}(undef)))
            # push!(first_pass, Expr(:(=), out, :(CorrCholesky{$M,$T,$lkj_length}(undef))))
        end
        push!(first_pass, :($seedlkj = $m.DistributionParameters.StructuredMatrices.PtrLowerTriangularMatrix{$(M-1),$T,$N}(pointer($∂θ))))
        push!(first_pass, :($lkjlogdetsym = $m.DistributionParameters.∂lkj_constrain!($out, $seedlkj, $lkjjacobiansym, $zsym)))
        push!(first_pass, :($adjout = $m.alloc_adjoint($out)))
        push!(first_pass, :($m.lifetime_start!($adjout)))
        push!(first_pass, :($θ += $N))
        push!(first_pass, :($∂θ += $N))
    else
        allocq = if use_sptr
            push!( first_pass, :( $out = $m.DistributionParameters.PtrCorrCholesky{$M,$T,$lkj_length}(pointer($sp,$T)) ) )
            # push!( first_pass, Expr(:(=), out, :(PtrCorrCholesky{$M,$T,$lkj_length}(pointer($sp,$T)) ) ))
            spoffset = sizeof(T)*lkj_length
            if !exportparam
                spoffset = VectorizationBase.align(spoffset)
            end
            push!( first_pass, :( $sp += $spoffset ) )
        else
            push!(first_pass, :( $out = $m.DistributionParameters.CorrCholesky{$M,$T,$lkj_length}(undef) ) )
            # push!(first_pass, Expr(:(=), out, :(CorrCholesky{$M,$T,$lkj_length}(undef) ) ))
        end
        push!(first_pass, :($lkjlogdetsym = $m.DistributionParameters.lkj_constrain!($out, $zsym)))
        push!(first_pass, :($θ += $N))
        limitlife && push!(second_pass, :($m.lifetime_end!($out)))
    end
    limitlife && push!(first_pass, :($m.lifetime_end!($zsym)))
    logjac && push!(first_pass, :(target = $m.vadd(target, $lkjlogdetsym)))
    nothing
end

function load_parameter!(
    first_pass::Vector{Any}, second_pass::Vector{Any}, out::Symbol,
    ::Type{<: AbstractCorrCholesky{M}}, partial::Bool = false,
    m::Module = DistributionParameters, sp::Union{Symbol,Nothing} = nothing,
    logjac::Bool = true, exportparam::Bool = false
) where {M}
    load_parameter!(first_pass, second_pass, out, CorrCholesky{M,Float64}, partial, m, sp, logjac, exportparam)
end

function parameter_names(::Type{<:StructuredMatrices.AbstractLowerTriangularMatrix{M}}, s::Symbol) where {M}
    ss = strip_hashtags(s)
    names = Vector{String}(undef, binomial2(M+1))
    for m ∈ 1:M
        names[m] = "$ss[$m,$m]"
    end
    ind = M
    for c ∈ 1:M-1, r ∈ c+1:M
        ind += 1
        names[ind] = "$ss[$r,$c]"
    end
    names::Vector{String}
end
parameter_names(x::LT) where {LT <: StructuredMatrices.AbstractLowerTriangularMatrix} = param_names(LT)

@generated function LinearAlgebra.mul!(
    σL::AbstractCovarCholesky{M,T,MTV},
    σ::Diagonal{T,<:AbstractFixedSizeVector{M,T,MV}},
    L::AbstractCorrCholesky{M,T,MTR}
) where {M,MTV,MTR,MV,T}
    q = quote
        σLpt = StructuredMatrices.PtrLowerTriangularMatrix{$M,Float64,$MTR}(pointer(σL))
        Lpt = StructuredMatrices.PtrLowerTriangularMatrix{$M,Float64,$MTR}(pointer(L))
        mul!(σLpt, σ, L)
    end
    σLcaches_ld = first(caches_logdiag(M, MTV, T))
    # tri_length = VectorizationBase.align(binomial2(M+1), T)
    # If σL does not have space, we don't calculate it.
    if !σLcaches_ld
        push!(q.args, :σL)
        return q
    end
    # We need at least one loop
    loopbody = quote
        logdiag_σL[m] = logσ[m] + logdiag_L[m]
    end
    # Do we also calculate invdiag?
    σLcaches_inv = first(caches_invdiag(M, MTV, T))
    if σLcaches_inv
        push!(q.args, :(invdiag_σL = invdiag(σL)))
        # push!(q.args, :(invdiag_L = invdiag(L)))
        push!(loopbody.args, :( invdiag_σL[m] = inv(σL[m]) ))
    end
    calc_logdiag_q = quote
        logdiag_σL = logdiag(σL)
        logdiag_L = logdiag(L)
        logσ = LazyMap(SLEEFPirates.log, σ.diag)
        @vvectorize $T for m ∈ 1:$MV
            $loopbody
        end
        σL
    end
    push!(q.args, calc_logdiag_q)
    q
end

@inline function Base.:*(
    sp::StackPointer,
    σ::Diagonal{T,<:AbstractFixedSizeVector{M,T,MV}},
    L::AbstractCorrCholesky{M,T,MTR}
) where {M,T,MV,MTR}
    (sp, σL) = PtrCovarCholesky{M,T}(sp)
    (sp, mul!(σL, σ, L))
end


text:
using LinearAlgebraicRepresentation
LAR = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

function randomshapes(n::Int=100)
	V, EV = LAR.randomcuboids(n, .4)
	V = GL.normalize2(V, flag=true)
	GL.VIEW([ GL.GLLines(V, EV, GL.COLORS[1]) ])

	W, EW = LAR.fraglines(1.5,1.5,1.5)((V,EV))
	GL.VIEW([ GL.GLLines(W, EW, GL.COLORS[1]) ])

	model = Lar.Model(V)
	Lar.addModelCells!(model, 1, EV)
	model = Lar.Arrangement.planar_arrangement(model)

	FVs = Lar.triangulate2D(model)
	GL.VIEW(GL.GLExplode(model.G, FVs, 1.2,1.2,1.2,99,1));
	GL.VIEW(GL.GLExplode(model.G, FVs, 1,1,1,99,1));
end


text:
module ArrayInterfaceBlockBandedMatrices

using ArrayInterfaceCore
using ArrayInterfaceBandedMatrices
using BlockArrays
using BlockBandedMatrices

struct BlockBandedMatrixIndex <: ArrayInterfaceCore.MatrixIndex
    count::Int
    refinds::Array{Int,1}
    refcoords::Array{Int,1}# storing col or row inds at ref points
    isrow::Bool
end
Base.firstindex(i::BlockBandedMatrixIndex) = 1
Base.lastindex(i::BlockBandedMatrixIndex) = i.count
Base.length(i::BlockBandedMatrixIndex) = lastindex(i)
function BlockBandedMatrixIndex(nrowblock, ncolblock, rowsizes, colsizes, l, u)
    blockrowind = ArrayInterfaceBandedMatrices.BandedMatrixIndex(nrowblock, ncolblock, l, u, true)
    blockcolind = ArrayInterfaceBandedMatrices.BandedMatrixIndex(nrowblock, ncolblock, l, u, false)
    sortedinds = sort(
        [(blockrowind[i], blockcolind[i]) for i = 1:length(blockrowind)],
        by=x -> x[1],
    )
    sort!(sortedinds, by=x -> x[2], alg=InsertionSort)# stable sort keeps the second index in order
    refinds = Array{Int,1}()
    refrowcoords = Array{Int,1}()
    refcolcoords = Array{Int,1}()
    rowheights = pushfirst!(copy(rowsizes), 1)
    cumsum!(rowheights, rowheights)
    blockheight = 0
    blockrow = 1
    blockcol = 1
    currenti = 1
    lastrowind = sortedinds[1][1] - 1
    lastcolind = sortedinds[1][2]
    for ind in sortedinds
        rowind, colind = ind
        if colind == lastcolind
            if rowind > lastrowind
                blockheight += rowsizes[rowind]
            end
        else
            for j = blockcol:blockcol+colsizes[lastcolind]-1
                push!(refinds, currenti)
                push!(refrowcoords, blockrow)
                push!(refcolcoords, j)
                currenti += blockheight
            end
            blockcol += colsizes[lastcolind]
            blockrow = rowheights[rowind]
            blockheight = rowsizes[rowind]
        end
        lastcolind = colind
        lastrowind = rowind
    end
    for j = blockcol:blockcol+colsizes[lastcolind]-1
        push!(refinds, currenti)
        push!(refrowcoords, blockrow)
        push!(refcolcoords, j)
        currenti += blockheight
    end
    push!(refinds, currenti)# guard
    push!(refrowcoords, -1)
    push!(refcolcoords, -1)
    rowindobj = BlockBandedMatrixIndex(currenti - 1, refinds, refrowcoords, true)
    colindobj = BlockBandedMatrixIndex(currenti - 1, refinds, refcolcoords, false)
    rowindobj, colindobj
end
Base.@propagate_inbounds function Base.getindex(ind::BlockBandedMatrixIndex, i::Int)
    @boundscheck 1 <= i <= ind.count || throw(BoundsError(ind, i))
    p = 1
    while i - ind.refinds[p] >= 0
        p += 1
    end
    p -= 1
    _i = i - ind.refinds[p]
    if ind.isrow
        return ind.refcoords[p] + _i
    else
        return ind.refcoords[p]
    end
end

function ArrayInterfaceCore.findstructralnz(x::BlockBandedMatrices.BlockBandedMatrix)
    l, u = BlockBandedMatrices.blockbandwidths(x)
    nrowblock = BlockBandedMatrices.blocksize(x, 1)
    ncolblock = BlockBandedMatrices.blocksize(x, 2)
    rowsizes = BlockArrays.blocklengths(axes(x, 1))
    colsizes = BlockArrays.blocklengths(axes(x, 2))
    return BlockBandedMatrixIndex(
        nrowblock,
        ncolblock,
        rowsizes,
        colsizes,
        l,
        u,
    )
end
struct BandedBlockBandedMatrixIndex <: ArrayInterfaceCore.MatrixIndex
    count::Int
    refinds::Array{Int,1}
    refcoords::Array{Int,1}# storing col or row inds at ref points
    reflocalinds::Array{ArrayInterfaceBandedMatrices.BandedMatrixIndex,1}
    isrow::Bool
end
Base.firstindex(i::BandedBlockBandedMatrixIndex) = 1
Base.lastindex(i::BandedBlockBandedMatrixIndex) = i.count
Base.length(i::BandedBlockBandedMatrixIndex) = lastindex(i)
Base.@propagate_inbounds function Base.getindex(ind::BandedBlockBandedMatrixIndex, i::Int)
    @boundscheck 1 <= i <= ind.count || throw(BoundsError(ind, i))
    p = 1
    while i - ind.refinds[p] >= 0
        p += 1
    end
    p -= 1
    _i = i - ind.refinds[p] + 1
    ind.reflocalinds[p][_i] + ind.refcoords[p] - 1
end

function BandedBlockBandedMatrixIndex(
    nrowblock,
    ncolblock,
    rowsizes,
    colsizes,
    l,
    u,
    lambda,
    mu,
)
    blockrowind = ArrayInterfaceBandedMatrices.BandedMatrixIndex(nrowblock, ncolblock, l, u, true)
    blockcolind = ArrayInterfaceBandedMatrices.BandedMatrixIndex(nrowblock, ncolblock, l, u, false)
    sortedinds = sort(
        [(blockrowind[i], blockcolind[i]) for i = 1:length(blockrowind)],
        by=x -> x[1],
    )
    sort!(sortedinds, by=x -> x[2], alg=InsertionSort)# stable sort keeps the second index in order
    rowheights = pushfirst!(copy(rowsizes), 1)
    cumsum!(rowheights, rowheights)
    colwidths = pushfirst!(copy(colsizes), 1)
    cumsum!(colwidths, colwidths)
    currenti = 1
    refinds = Array{Int,1}()
    refrowcoords = Array{Int,1}()
    refcolcoords = Array{Int,1}()
    reflocalrowinds = Array{ArrayInterfaceBandedMatrices.BandedMatrixIndex,1}()
    reflocalcolinds = Array{ArrayInterfaceBandedMatrices.BandedMatrixIndex,1}()
    for ind in sortedinds
        rowind, colind = ind
        localrowind =
            ArrayInterfaceBandedMatrices.BandedMatrixIndex(rowsizes[rowind], colsizes[colind], lambda, mu, true)
        localcolind =
            ArrayInterfaceBandedMatrices.BandedMatrixIndex(rowsizes[rowind], colsizes[colind], lambda, mu, false)
        push!(refinds, currenti)
        push!(refrowcoords, rowheights[rowind])
        push!(refcolcoords, colwidths[colind])
        push!(reflocalrowinds, localrowind)
        push!(reflocalcolinds, localcolind)
        currenti += localrowind.count
    end
    push!(refinds, currenti)
    push!(refrowcoords, -1)
    push!(refcolcoords, -1)
    rowindobj = BandedBlockBandedMatrixIndex(
        currenti - 1,
        refinds,
        refrowcoords,
        reflocalrowinds,
        true,
    )
    colindobj = BandedBlockBandedMatrixIndex(
        currenti - 1,
        refinds,
        refcolcoords,
        reflocalcolinds,
        false,
    )
    rowindobj, colindobj
end

function ArrayInterfaceCore.findstructralnz(x::BlockBandedMatrices.BandedBlockBandedMatrix)
    l, u = BlockBandedMatrices.blockbandwidths(x)
    lambda, mu = BlockBandedMatrices.subblockbandwidths(x)
    nrowblock = BlockBandedMatrices.blocksize(x, 1)
    ncolblock = BlockBandedMatrices.blocksize(x, 2)
    rowsizes = BlockArrays.blocklengths(axes(x, 1))
    colsizes = BlockArrays.blocklengths(axes(x, 2))
    return BandedBlockBandedMatrixIndex(
        nrowblock,
        ncolblock,
        rowsizes,
        colsizes,
        l,
        u,
        lambda,
        mu,
    )
end

ArrayInterfaceCore.has_sparsestruct(::Type{<:BlockBandedMatrices.BlockBandedMatrix}) = true
ArrayInterfaceCore.has_sparsestruct(::Type{<:BlockBandedMatrices.BandedBlockBandedMatrix}) = true
ArrayInterfaceCore.isstructured(::Type{<:BlockBandedMatrices.BlockBandedMatrix}) = true
ArrayInterfaceCore.isstructured(::Type{<:BlockBandedMatrices.BandedBlockBandedMatrix}) = true
ArrayInterfaceCore.fast_matrix_colors(::Type{<:BlockBandedMatrices.BlockBandedMatrix}) = true
ArrayInterfaceCore.fast_matrix_colors(::Type{<:BlockBandedMatrices.BandedBlockBandedMatrix}) = true

function ArrayInterfaceCore.matrix_colors(A::BlockBandedMatrices.BlockBandedMatrix)
    l, u = BlockBandedMatrices.blockbandwidths(A)
    blockwidth = l + u + 1
    nblock = BlockBandedMatrices.blocksize(A, 2)
    cols = BlockArrays.blocklengths(axes(A, 2))
    blockcolors = ArrayInterfaceCore._cycle(1:blockwidth, nblock)
    # the reserved number of colors of a block is the maximum length of columns of blocks with the same block color
    ncolors = [maximum(cols[i:blockwidth:nblock]) for i = 1:blockwidth]
    endinds = cumsum(ncolors)
    startinds = [endinds[i] - ncolors[i] + 1 for i = 1:blockwidth]
    colors = [
        (startinds[blockcolors[i]]:endinds[blockcolors[i]])[1:cols[i]]
        for i = 1:nblock
    ]
    return reduce(vcat, colors)
end

function ArrayInterfaceCore.matrix_colors(A::BlockBandedMatrices.BandedBlockBandedMatrix)
    l, u = BlockBandedMatrices.blockbandwidths(A)
    lambda, mu = BlockBandedMatrices.subblockbandwidths(A)
    blockwidth = l + u + 1
    subblockwidth = lambda + mu + 1
    nblock = BlockBandedMatrices.blocksize(A, 2)
    cols = BlockArrays.blocklengths(axes(A, 2))
    blockcolors = ArrayInterfaceCore._cycle(1:blockwidth, nblock)
    # the reserved number of colors of a block is the min of subblockwidth and the largest length of columns of blocks with the same block color
    ncolors = [
        min(subblockwidth, maximum(cols[i:blockwidth:nblock]))
        for i = 1:min(blockwidth, nblock)
    ]
    endinds = cumsum(ncolors)
    startinds = [endinds[i] - ncolors[i] + 1 for i = 1:min(blockwidth, nblock)]
    colors = [
        ArrayInterfaceCore._cycle(startinds[blockcolors[i]]:endinds[blockcolors[i]], cols[i])
        for i = 1:nblock
    ]
    return reduce(vcat, colors)
end

end # module


text:
using Test
using OpenQASM
using YaoTargetQASM

module QeLib

using YaoTargetQASM
qelib = joinpath(pkgdir(YaoTargetQASM), "test", "qelib1.inc")

qasm"""OPENQASM 2.0;
include "$qelib";

gate post q {x q;}
"""

end

ast = code_qasm(typeof(QeLib.s()))
@test ast ≈ OpenQASM.parse("""
OPENQASM 2.0;
qreg qreg_1[1];
rz(0) qreg_1;
ry(0) qreg_1;
rz(1.5707963267948966) qreg_1;
""")

ast = code_qasm(typeof(QeLib.x()))
@test ast ≈ OpenQASM.parse("""
OPENQASM 2.0;
qreg qreg_1[1];
rz(pi) qreg_1;
ry(0) qreg_1;
rz(pi) qreg_1;
""")

ast = code_qasm(typeof(QeLib.y()))
@test ast ≈ OpenQASM.parse("""
OPENQASM 2.0;
qreg qreg_1[1];
rz(pi) qreg_1;
ry(1.5707963267948966) qreg_1;
rz(1.5707963267948966) qreg_1;
""")

ast = code_qasm(typeof(QeLib.z()))
@test ast ≈ OpenQASM.parse("""
OPENQASM 2.0;
qreg qreg_1[1];
rz(0) qreg_1;
ry(0) qreg_1;
rz(pi) qreg_1;
""")

ast = code_qasm(typeof(QeLib.h()))
@test ast ≈ OpenQASM.parse("""
OPENQASM 2.0;
qreg qreg_1[1];
rz(1.5707963267948966) qreg_1;
ry(0) qreg_1;
rz(pi) qreg_1;
""")

ast = code_qasm(typeof(QeLib.rx(2.0)); toplevel=false)
@test ast ≈ OpenQASM.parse_gate("""
gate rx(theta) qreg_1 {
  rz(theta) qreg_1;
  ry(-1.5707963267948966) qreg_1;
  rz(1.5707963267948966) qreg_1;
}
""")

ast = code_qasm(typeof(QeLib.ccx()))
@test ast ≈ OpenQASM.parse("""
OPENQASM 2.0;
qreg qreg_1[3];
h qreg_1[2];
CX qreg_1[1], qreg_1[2];
rz(0) qreg_1[2];
ry(0) qreg_1[2];
rz(-0.7853981633974483) qreg_1[2];
CX qreg_1[0], qreg_1[2];
rz(0) qreg_1[2];
ry(0) qreg_1[2];
rz(0.7853981633974483) qreg_1[2];
CX qreg_1[1], qreg_1[2];
rz(0) qreg_1[2];
ry(0) qreg_1[2];
rz(-0.7853981633974483) qreg_1[2];
CX qreg_1[0], qreg_1[2];
rz(0) qreg_1[1];
ry(0) qreg_1[1];
rz(0.7853981633974483) qreg_1[1];
rz(0) qreg_1[2];
ry(0) qreg_1[2];
rz(0.7853981633974483) qreg_1[2];
h qreg_1[2];
CX qreg_1[0], qreg_1[1];
rz(0) qreg_1[0];
ry(0) qreg_1[0];
rz(0.7853981633974483) qreg_1[0];
rz(0) qreg_1[1];
ry(0) qreg_1[1];
rz(-0.7853981633974483) qreg_1[1];
CX qreg_1[0], qreg_1[1];
""")

ast = code_qasm(typeof(QeLib.crz(2.0)); toplevel=false)
@test ast ≈ OpenQASM.parse_gate("""
gate crz(lambda) qreg_2, qreg_1 {
  rz(0) qreg_1;
  ry(0) qreg_1;
  rz(lambda/2.0) qreg_1;
  CX qreg_2, qreg_1;
  rz(0) qreg_1;
  ry(0) qreg_1;
  rz(-(lambda)/2.0) qreg_1;
  CX qreg_2, qreg_1;
}
""")

module Adder3

using YaoTargetQASM
qelib = joinpath(pkgdir(YaoTargetQASM), "test", "qelib1.inc")

const adder3 = qasm"""
OPENQASM 2.0;
include "$qelib";
gate majority a,b,c
{
  cx c,b;
  cx c,a;
  ccx a,b,c;
}
gate unmaj a,b,c
{
  ccx a,b,c;
  cx c,a;
  cx a,b;
}
qreg cin[1];
qreg a[4];
qreg b[4];
qreg cout[1];
creg ans[5];
// set input states
x a[0]; // a = 0001
x b;    // b = 1111
// add a to b, storing result in b
majority cin[0],b[0],a[0];
majority a[0],b[1],a[1];
majority a[1],b[2],a[2];
majority a[2],b[3],a[3];
cx a[3],cout[0];
unmaj a[2],b[3],a[3];
unmaj a[1],b[2],a[2];
unmaj a[0],b[1],a[1];
unmaj cin[0],b[0],a[0];
measure b[0] -> ans[0];
measure b[1] -> ans[1];
measure b[2] -> ans[2];
measure b[3] -> ans[3];
measure cout[0] -> ans[4];
"""
end


text:
include("../src/PlutoClient.jl")

import DandelionWebSockets: getrequestheaders, HTTPHandshakeLogic

notebook_id, cookie = ARGS

getrequestheaders(h::HTTPHandshakeLogic) = [
        "Sec-WebSocket-Version" => "13",
        "Upgrade" => "websocket",
        "Connection" => "Upgrade",
        "Sec-WebSocket-Key" => h.key,
        "Cookie" => "secret=$(cookie)"]

client = PlutoClient(notebook_id=notebook_id)
ok = connect(client, "127.0.0.1:1235")
if !ok
        @error "Failed to connect to the notebook"
        exit(1)
end

while true
        stats = Dict(st.first => compute_percentage(st.second, st.first) for st in get_cpu_stats())
        set_bond_value(handler, "bonds/a", stats)
        sleep(4)
end

take!(client.stop_channel)


text:
export richardson_lucy_iterative

"""
    richardson_lucy_iterative(measured, psf; <keyword arguments>)

Classical iterative Richardson-Lucy iteration scheme for deconvolution.
`measured` is the measured array and `psf` the point spread function.
Converges slower than the optimization approach of `deconvolution`

# Keyword Arguments
- `regularizer=GR()`: A regularizer function. Can be exchanged
- `λ=0.05`: A float indicating the total weighting of the regularizer with
    respect to the global loss function
- `iterations=100`: Specifies number of iterations.
- `progress`: if not `nothing`, the progress will be monitored in a summary dictionary as obtained by
              DeconvOptim.options_trace_deconv()

# Example
```julia-repl
julia> using DeconvOptim, TestImages, Colors, Noise;

julia> img = Float32.(testimage("resolution_test_512"));

julia> psf = Float32.(generate_psf(size(img), 30));

julia> img_b = conv(img, psf);

julia> img_n = poisson(img_b, 300);

julia> @time res = richardson_lucy_iterative(img_n, psf);
```
"""
function richardson_lucy_iterative(measured, psf;
                                   regularizer=GR(),
                                   λ=0.05,
                                   iterations=100,
                                   conv_dims=1:ndims(psf),
                                   progress = nothing)

    otf, conv_temp = plan_conv(measured, psf, conv_dims)
    otf_conj = conj.(otf)
    # initializer
    rec = abs.(conv_temp(measured, otf))#ones(eltype(measured), size(measured))

    # buffer for gradient
    # we need Base.invokelatest because of world age issues with generated
    # regularizers
    buffer_grad =  let
        if !isnothing(regularizer)
            Base.invokelatest(gradient, regularizer, rec)[1]
        else
            nothing
        end
    end

    ∇reg(x) = buffer_grad .= Base.invokelatest(gradient, regularizer, x)[1]

    buffer = copy(measured)

    iter_without_reg(rec) = begin
        buffer .= measured ./ (conv_temp(rec, otf))
        conv_temp(buffer, otf_conj)
    end
    iter_with_reg(rec) = buffer .= (iter_without_reg(rec) .- λ .* ∇reg(rec))

    iter = isnothing(regularizer) ? iter_without_reg : iter_with_reg

    loss(myrec) = begin
        fwd = conv_temp(myrec, otf)
        return sum(fwd .- measured .* log.(fwd))
    end

    tmp_time = 0.0
    if progress !== nothing
        record_progress!(progress, rec, 0, loss(rec), 0.0, 1.0)
        tmp_time=time()
    end
    code_time = 0.0
    for i in 1:iterations
        rec .*= iter(rec)
        if progress !== nothing
            # do not count the time for evaluating the loss here.
            code_time += time() .- tmp_time
            record_progress!(progress, copy(rec), i, loss(rec), code_time, 1.0)
            tmp_time=time()
        end
    end

    return rec
end


text:
# An example of how the solution for the wave equation
# using a gaussian initial condition converges

using ScalarWaveEquation
using PyPlot

function get_gaussian(N::Int)
	# returns a normalized Gaussian distribution with (mean, std) = (0.5, 0.1)
	mean = 0.5
	std  = 0.1

	xdata      = range(0.0, stop = 1.0, length = N)
	exponent   = -0.5 * (xdata.-mean) .* (xdata.-mean) / (std*std)
	norm_coeff = std*sqrt(2*pi)

	return exp.(exponent)/norm_coeff
end

function get_gaussian_solution(Nx::Int, Nt::Int, Ttot::T) where {T}
	# returns a FunArr2D instance for the solution to the gaussian
	# initial condition, performed at difference resolutions

	# initial conditions are accepted in the form of FunArr instances
	# the x-range is assumed to be [0, 1]
	init_positions  = get_gaussian(Nx)
	init_velocities = zeros(Nx)
	ufun  = FunArr{Float64, Float64}(init_positions)
	utfun = FunArr{Float64, Float64}(init_velocities)

	# use the method of lines to calculate u(x, t) over [0, 1]x[0, Ttot]
	u2fun, energies = scalarwave_and_energy(ufun, utfun, Nt, Ttot)

	return u2fun
end

# get the solution at different resolutions for the same duration
u2fun1 = get_gaussian_solution(100, 150, 0.2)
u2fun2 = get_gaussian_solution(200, 300, 0.2)
u2fun3 = get_gaussian_solution(300, 450, 0.2)
u2fun4 = get_gaussian_solution(400, 600, 0.2)
u2fun5 = get_gaussian_solution(500, 750, 0.2)
u2fun6 = get_gaussian_solution(600, 900, 0.2)

# get the integrated squared difference between solutions
# for successive resolutions, for a 2000x2000 spacing
area12 = integrate_difference(u2fun1, u2fun2, 2000, 2000)
area23 = integrate_difference(u2fun2, u2fun3, 2000, 2000)
area34 = integrate_difference(u2fun3, u2fun4, 2000, 2000)
area45 = integrate_difference(u2fun4, u2fun5, 2000, 2000)
area56 = integrate_difference(u2fun5, u2fun6, 2000, 2000)

# the differences between increasing resolutions decreases
println("The integrated difference between 100x150 and 200x300 is $area12")
println("The integrated difference between 200x300 and 300x450 is $area23")
println("The integrated difference between 300x450 and 400x600 is $area34")
println("The integrated difference between 400x600 and 500x750 is $area45")
println("The integrated difference between 500x750 and 600x900 is $area56")


text:
function Base.getproperty(monitor::Monitor, name::Symbol)
	# TODO for GLFW 3.3: contentscale is GetMonitorContentScale
	if name == :isprimary
		GetPrimaryMonitor() == monitor
	elseif name == :name
		GetMonitorName(monitor)
	elseif name == :physicalsize
		GetMonitorPhysicalSize(monitor)
	elseif name == :position
		GetMonitorPos(monitor)
	elseif name == :supported_videomodes
		GetVideoModes(monitor)
	elseif name == :videomode
		GetVideoMode(monitor)
	else
		getfield(monitor, name)
	end
end


text:
################################################################################
#
#    This file contains types for a variety of empirical atomic potentials
#
################################################################################

abstract type ArbitraryPotential end
abstract type EmpiricalPotential <: ArbitraryPotential end
abstract type MixedPotential <: ArbitraryPotential end

################################################################################
# InteratomicPotentials API default generic implmentations
################################################################################
include("ArbitraryPotential/arbitrary_potential.jl")
export ArbitraryPotential

################################################################################
# Empirical Potentials
################################################################################
include("EmpiricalPotentials/empirical_potentials.jl")
export EmpiricalPotential

################################################################################
# Mixed Potentials
################################################################################
include("MixedPotentials/mixed_potentials.jl")
export MixedPotential


text:
"""
get_edges(g,node_x,node_y)
g is a SimpleWeightedGraph, while node_x and node_y are
arrays containing the position of the nodes.

Return a DataFrame containing the edge information,
such as position, pairs and weights.
"""
function get_edges(g,node_x,node_y)
    edges_p1 = []
    edges_p2 = []
    edges_w  = Float64[]
    for i in edges(g)
        push!(edges_p1, i.src)
        push!(edges_p1, i.src)
        push!(edges_p2, i.dst)
        push!(edges_p2, i.dst)
        try
            push!(edges_w, i.weight)
            push!(edges_w, i.weight)
        catch e
            push!(edges_w,1)
            push!(edges_w,1)
        end
    end

    edges_node  = Int[]
    edges_pairs = []
    for i in 1:size(edges_p1)[1]
        if i % 2 == 0
            push!(edges_node,edges_p2[i])
        else
            push!(edges_node,edges_p1[i])
        end
        push!(edges_pairs,[edges_p1[i],edges_p2[i]])
    end
    graph_edges = DataFrames.DataFrame(edges_x = node_x[edges_node], edges_y = node_y[edges_node],
                       pairs = edges_pairs, ew = edges_w, node = edges_node)
    return graph_edges
  end


text:
# ---
# title: 780. Reaching Points
# id: problem780
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math
# link: <https://leetcode.com/problems/reaching-points/description/>
# hidden: true
# ---
#
# A move consists of taking a point `(x, y)` and transforming it to either `(x,
# x+y)` or `(x+y, y)`.
#
# Given a starting point `(sx, sy)` and a target point `(tx, ty)`, return `True`
# if and only if a sequence of moves exists to transform the point `(sx, sy)` to
# `(tx, ty)`. Otherwise, return `False`.
#
#
#
#     **Examples:**
#     Input: sx = 1, sy = 1, tx = 3, ty = 5
#     Output: True
#     Explanation:
#     One series of moves that transforms the starting point to the target is:
#     (1, 1) -> (1, 2)
#     (1, 2) -> (3, 2)
#     (3, 2) -> (3, 5)
#
#     Input: sx = 1, sy = 1, tx = 2, ty = 2
#     Output: False
#
#     Input: sx = 1, sy = 1, tx = 1, ty = 1
#     Output: True
#
#
#
# **Note:**
#
#   * `sx, sy, tx, ty` will all be integers in the range `[1, 10^9]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
####################################################################
#  JSON schema definition and parsing
####################################################################

import Base: setindex!, getindex, haskey, show

## transforms escaped characters in JPaths back to their original value
function unescapeJPath(raw::String)
  ret = replace(raw, "~0" => "~")
  ret = replace(ret, "~1" => "/")

  m = match(r"%([0-9A-F]{2})", ret)
  if m != nothing
    repls = Dict()
    for c in m.captures
      # c = first(m.captures)
      haskey(repls, String("%$c")) && continue
      repls[String("%$c")] = "$(Char(Meta.parse("0x" * c)))"
    end

    for (k,v) in repls
      ret = replace(ret, k => v)
    end
  end

  ret
end

## decomposes an URI into its elements
function todict(uri::HTTP.URI)
  pdict = Dict( fn => getfield(uri, fn) for fn in fieldnames(HTTP.URI))
  delete!(pdict, :uri)
end

## removes the fragment part of an URI
function rmfragment(uri::HTTP.URI)
  els = todict(uri)
  delete!(els, :fragment)
  HTTP.URI(;els...)
end

## calculates the full id of the explored item using the parent URI in
#   'id0', and the 'id' property in 's'
function updateid(id0::HTTP.URI, s::String)
  id2 = HTTP.URI(s)
  id2.scheme != "" && return id2

  els = todict(id0)
  if id2.path != ""  # replace path of id0
    oldpath = match(r"^(.*/).*$", id0.path)
    if oldpath == nothing
      els[:path] = id2.path
    else
      els[:path] = oldpath.captures[1] * id2.path
    end
  end

  els[:fragment] = id2.fragment

  HTTP.URI(;els...)
end

## constructs the map of ids to schema elements
mkidmap!(map, x, id0::HTTP.URI) = nothing
mkidmap!(map, v::Vector, id0::HTTP.URI) = foreach(e -> mkidmap!(map,e,id0), v)
function mkidmap!(map::Dict, el::Dict, id0::HTTP.URI)
  if haskey(el, "id") && isa(el["id"], String) # draft 04
    id0 = updateid(id0, el["id"])
    map[string(id0)] = el
  elseif haskey(el, "\$id") && isa(el["\$id"], String) # draft 06+
    id0 = updateid(id0, el["\$id"])
    map[string(id0)] = el
  end

  for (k,v) in el
    mkidmap!(map, v, id0)
  end
end

## searches the element refered to by JSPointer/fragment in path in the schema s
function findelement(s, path)
  for el in split(path, "/")
    realel = unescapeJPath(String(el))
    (realel == "") && continue
    if s isa Dict
      haskey(s, realel) || error("missing property '$realel' in $s")
      s = s[realel]
    elseif s isa Vector
      idx = Meta.parse(realel)
      (idx isa Int) || error("expected numeric array index instead of '$realel'")
      idx += 1
      (length(s) < idx) && error("item index $(idx-1) larger than array $s")
      s = s[idx]
    else
      error("unmanaged type in ref resolution $(typeof(s)) - $s")
    end
  end
  s
end

## fetch remote schema
function getremoteschema(uri::HTTP.URI)
  r = HTTP.request("GET", uri)
  (r.status != 200) && error("remote ref $uri not found")
  Schema(JSON.parse(String(r.body))) # process remote ref
end

function findref(id0, idmap, path::String, parentFileDirectory::String)
  # path refers to root
  (path in ["", "#"]) && return idmap[string(id0)]

  # path is a JPointer
  (length(path) > 1 ) && (path[1:2] == "#/") &&
    return findelement(idmap[string(id0)], path[3:end])

  # path is a URI
  uri = updateid(id0, path) # fullRefURI(HTTP.URI(path), id0)
  uri2 = rmfragment(uri) # without JPointer

  isFileUri = startswith(uri2.scheme, "file") || isempty(uri2.scheme)
  # normalize a file path to an absolute path so creating a key is consistent
  if isFileUri
    if !isabspath(uri2.path)
      uri2 = HTTP.URIs.merge(uri2; path = abspath(joinpath(parentFileDirectory, uri2.path)))
    end
  end

  if !haskey(idmap, string(uri2))  # if not referenced already, fetch remote ref, add to idmap
    if startswith(uri2.scheme, "http")
      @info("fetching remote ref $(uri2)")
      idmap[string(uri2)] = getremoteschema(uri2).data
    elseif isFileUri
      @info("loading local ref $(uri2)")
      idmap[string(uri2)] = Schema(JSON.parsefile(uri2.path); parentFileDirectory = dirname(uri2.path)).data
    end
  end

  findelement(idmap[string(uri2)], uri.fragment)
end

# finds recursively all "$ref" and resolve their path
resolverefs!(s, id0, idmap, parentFileDirectory) = nothing
resolverefs!(s::Vector, id0, idmap, parentFileDirectory) = foreach(e -> resolverefs!(e, id0, idmap, parentFileDirectory), s)
function resolverefs!(s::Dict, id0, idmap, parentFileDirectory::String)
  if haskey(s, "id") && isa(s["id"], String) # draft 04
    id0 = updateid(id0, s["id"])
  elseif haskey(s, "\$id") && isa(s["\$id"], String) # draft 06+
    id0 = updateid(id0, s["\$id"])
  end

  for (k,v) in s
    if (k == "\$ref") && (v isa String)
      # This ref has not been resolved yet (otherwise it would not be a String)
      # We will replace the path string with the schema element pointed at, thus marking it as
      # resolved. This should prevent infinite recursions caused by self referencing
      s["\$ref"] = findref(id0, idmap, v, parentFileDirectory)
    else
      resolverefs!(v, id0, idmap, parentFileDirectory)
    end
  end
end

################################################################################
#  Schema struct definition
################################################################################

"""
`Schema(sch::String)`

Create a schema for document validation. `sch` should be a String containing a
valid JSON.

`Schema(sch::Dict)`

Create a schema but with `sch` being a parsed JSON created with `JSON.parse()` or
`JSON.parsefile()`.

## Example

```julia
julia> myschema = Schema("
  {
    \"properties\": {
      \"foo\": {},
      \"bar\": {}
    },
    \"required\": [\"foo\"]
  }
  ")

julia> sch = JSON.parsefile(filepath)
julia> myschema = Schema(sch)
```
"""
struct Schema
  data::Union{Dict{String, Any}, Bool}

  function Schema(spec0::Bool; idmap0=Dict{String, Any}())
    new(spec0)
  end

  function Schema(sp::String; idmap0=Dict{String, Any}(), parentFileDirectory::String)
    Schema(JSON.parse(sp), idmap0=idmap0, parentFileDirectory=parentFileDirectory)
  end

  function Schema(spec0::Dict; idmap0=Dict{String, Any}(), parentFileDirectory::String = abspath("."))
    spec  = deepcopy(spec0)
    idmap = deepcopy(idmap0)

    # construct dictionary of 'id' properties to resolve references later
    id0 = HTTP.URI()
    idmap[string(id0)] = spec
    mkidmap!(idmap, spec, id0)

    # resolve all refs to the corresponding schema elements
    resolverefs!(spec, id0, idmap, parentFileDirectory)

    new(spec)
  end
end

setindex!(x::Schema, val, key) = setindex!(x.data, val, key)
getindex(x::Schema, key) = getindex(x.data, key)
haskey(x::Schema, key) = haskey(x.data, key)

function show(io::IO, sch::Schema)
  show(io, Schema)
end