text:
module PrimordialB

using  Dierckx, JLD, PyCall #, SpecialFunctions
export
	# Custom types
	FFTgrid,
	MatrixCls,
	LenseDecomp,
    Pix,

	# methods
	lense,
	sim_xk,
	radial_power,
	qu2eb, eb2qu,
	squash, squash!

# FFTW.set_num_threads(CPU_CORES)
# FFTW.set_num_threads(3)

#=  To lint this file run:
using Lint
lintfile("src/PrimordialB.jl")
=#

#=##########################################################

Custom type Definitions

=##############################################################

# ---- Holds grid, model and planned FFT parameters for the quadratic estimate.
immutable FFTgrid{dm, T}
	period::Float64
	nside::Int64
	deltx::Float64
	deltk::Float64
	nyq::Float64
	x::Array{Array{Float64,dm},1}
	k::Array{Array{Float64,dm},1}
	r::Array{Float64,dm}
	FFT::T  # saved plan for fast fft
end

#---- Holds the cls expanded out to the 2 d spectral matrices.
immutable MatrixCls{dm}
	cϕϕk::Array{Float64,dm}
	cϕψk::Array{Float64,dm}
	cψψk::Array{Float64,dm}
	cTTk::Array{Float64,dm}
	cTEk::Array{Float64,dm}
	cEEk::Array{Float64,dm}
	cBBk::Array{Float64,dm}
	cBB0k::Array{Float64,dm}
	cTTnoisek::Array{Float64,dm}
	cEEnoisek::Array{Float64,dm}
	cBBnoisek::Array{Float64,dm}
end

# ----- Holds the decomposition of the lensing displacements
immutable LenseDecomp
	indcol::Array{Int64,2}
	indrow::Array{Int64,2}
	rdisplx::Array{Float64,2}
	rdisply::Array{Float64,2}
	displx::Array{Float64,2}
	disply::Array{Float64,2}
	ϕk::Array{Complex{Float64},2}
	ψk::Array{Complex{Float64},2}
end

##########################################################

# Type constructors

##############################################################

# TODO: fix this up for general dimension
immutable Pix{period, nside} end

@generated function FFTgrid{period, nside}(p::Pix{period, nside})
    deltx     = period / nside
	deltk     = 2π / period
    nyq       = 2π / (2deltx)
    k_side = ifftshift(-nside÷2:(nside-1)÷2) * deltk
    x_side = ifftshift(-nside÷2:(nside-1)÷2) * deltx
    zrs    = zeros(nside, nside)
    k      = [k_side.' .+ zrs, k_side .+ zrs]#[reshape(k_side, 1, nside), reshape(k_side, nside, 1)]
    x      = [x_side.' .+ zrs, x_side .+ zrs]#[reshape(x_side, 1, nside), reshape(x_side, nside, 1)]
    r      = sqrt.(k[1].^2 + k[2].^2)
    FFT    =  ((deltx^2) / (2π)) * plan_fft(rand(nside,nside); flags=FFTW.ESTIMATE, timelimit=5)
    return :(FFTgrid{2, typeof($(FFT))}($period, $nside, $deltx, $deltk, $nyq, $x, $k, $r, $(FFT)))
end

function MatrixCls{dm,T}(g::FFTgrid{dm,T}, cls; σTTrad=0.0, σEErad=0.0,  σBBrad=0.0, beamFWHM=0.0)
	cϕϕk = cls_to_cXXk(cls[:ell], cls[:ϕϕ], g.r)
	cϕψk = cls_to_cXXk(cls[:ell], cls[:ϕψ], g.r)
	cψψk = cls_to_cXXk(cls[:ell], cls[:ψψ], g.r)
	cTTk = cls_to_cXXk(cls[:ell], cls[:tt], g.r)
	cTEk = cls_to_cXXk(cls[:ell], cls[:te], g.r)
	cEEk = cls_to_cXXk(cls[:ell], cls[:ee], g.r)
	cBBk = cls_to_cXXk(cls[:ell], cls[:bb], g.r)
	cBB0k= cls_to_cXXk(cls[:ell], cls[:bb0], g.r)
	cTTnoisek = cNNkgen(g.r; σrad=σTTrad, beamFWHM=beamFWHM)
	cEEnoisek = cNNkgen(g.r; σrad=σEErad, beamFWHM=beamFWHM)
	cBBnoisek = cNNkgen(g.r; σrad=σBBrad, beamFWHM=beamFWHM)
	MatrixCls{dm}(cϕϕk, cϕψk, cψψk, cTTk, cTEk, cEEk, cBBk, cBB0k, cTTnoisek, cEEnoisek, cBBnoisek)
end

function LenseDecomp(ϕk, ψk, g)
	#displx = real(g.FFT \ (im .* g.k[1] .* ϕk) +  g.FFT \ (im .* g.k[2] .* ψk))
	#disply = real(g.FFT \ (im .* g.k[2] .* ϕk) -  g.FFT \ (im .* g.k[1] .* ψk))
	displx, disply = LenseDecomp_helper1(ϕk, ψk, g)

	row, col  = size(g.x[1])
	indcol    = Array(Int64, row, col)
	indrow	  = Array(Int64, row, col)
	rdisplx   = Array(Float64, row, col)
	rdisply   = Array(Float64, row, col)
	@inbounds for j = 1:col, i = 1:row
		round_displx_deltx = round(Int64, displx[i,j]/g.deltx)
		round_disply_deltx = round(Int64, disply[i,j]/g.deltx)
	    indcol[i,j]  = indexwrap(j + round_displx_deltx, col)
	    indrow[i,j]  = indexwrap(i + round_disply_deltx, row)
		rdisplx[i,j] = displx[i,j] - g.deltx * round_displx_deltx
	    rdisply[i,j] = disply[i,j] - g.deltx * round_disply_deltx
	end
	return LenseDecomp(indcol, indrow, rdisplx, rdisply, displx, disply, copy(ϕk), copy(ψk))
end
function LenseDecomp_helper1(ϕk, ψk, g)
	tmpdxk = Array(Complex{Float64}, size(g.r))
	tmpdyk = Array(Complex{Float64}, size(g.r))
	@inbounds @simd for i in eachindex(tmpdxk, tmpdyk)
		tmpdxk[i] = complex(im * g.k[1][i] * ϕk[i] + im * g.k[2][i] * ψk[i])
		tmpdyk[i] = complex(im * g.k[2][i] * ϕk[i] - im * g.k[1][i] * ψk[i])
	end
	displx = real(g.FFT \ tmpdxk)
	disply = real(g.FFT \ tmpdyk)
	return displx, disply
end

LenseDecomp(len::LenseDecomp, g) = LenseDecomp(copy(len.ϕk), copy(len.ψk), g)

##########################################################

# Helper functions for the type constructors

##############################################################

indexwrap(ind::Int64, uplim)  = mod(ind - 1, uplim) + 1

function cls_to_cXXk{dm}(ell, cxxls, r::Array{Float64, dm})
		spl = Dierckx.Spline1D(ell, cxxls; k=1, bc="zero", s=0.0)
		rtn = squash(map(spl, r))::Array{Float64, dm}
		rtn[r.==0.0] = 0.0
		return squash(map(spl, r))::Array{Float64, dm}
end

function cNNkgen{dm}(r::Array{Float64,dm}; σrad=0.0, beamFWHM=0.0)
	beamSQ = exp(- (beamFWHM ^ 2) * abs2(r) ./ (8 * log(2)) )    # my misunderstanding ? anyway we alway use zero beam here
    rtn = ones(size(r)) .* abs2(σrad) ./ beamSQ
    rtn[r.==0.0] = 0.0
    return rtn
end

function getgrid{T}(g::FFTgrid{2,T})
	xco_side, kco_side = getxkside(g)
	kco1, kco2 = meshgrid(kco_side, kco_side)
	xco1, xco2 = meshgrid(xco_side, xco_side)
	kco    = Array{Float64,2}[kco1, kco2]
	xco    = Array{Float64,2}[xco1, xco2]
	return xco, kco
end

function getxkside{dm,T}(g::FFTgrid{dm,T})
	deltx    = g.period / g.nside
	deltk    = 2π / g.period
	xco_side = zeros(g.nside)
	kco_side = zeros(g.nside)
	for j in 0:(g.nside-1)
		xco_side[j+1] = (j < g.nside/2) ? (j*deltx) : (j*deltx - g.period)
		kco_side[j+1] = (j < g.nside/2) ? (j*deltk) : (j*deltk - 2*π*g.nside/g.period)
	end
	xco_side, kco_side
end

function meshgrid(side_x,side_y)
    	nx = length(side_x)
    	ny = length(side_y)
    	xt = repmat(vec(side_x).', ny, 1)
    	yt = repmat(vec(side_y)  , 1 , nx)
    	return xt, yt
end

###############################

# local likelihood functions

###############################
include("likewhite.jl")

###############################

# functions for simulation code

###############################
include("simulation_src.jl")

###############################

# the quadratic estimate

###############################

""" EB quadratic estimate """
function eb_quad_est(dek, dbk, cEEfid, cBBfid, cEEobsfid, cBBobsfid, g, Aℓ)
    CX = 1./cEEobsfid
    CY = 1./cBBobsfid
    estk  = -im .* g.k[1] .* sin_conv(im .* g.k[1] .* cEEfid .* dek, dbk, CX, CY, g)
    estk += -im .* g.k[2] .* sin_conv(im .* g.k[2] .* cEEfid .* dek, dbk, CX, CY, g)
    estk += -im .* g.k[1] .* sin_conv(dek,  im .* g.k[1] .* cBBfid .* dbk, CX, CY, g)
    estk += -im .* g.k[2] .* sin_conv(dek,  im .* g.k[2] .* cBBfid .* dbk, CX, CY, g)
    rtnk = Aℓ .* estk
    rtnx = real(g.FFT \ rtnk)
    return rtnk, rtnx
end

""" Aℓ for EB quadratic estimate"""
function getAℓ(cEEtru, cBBtru, cEEfid, cBBfid, cEEobsfid, cBBobsfid, g)
    tmp = zeros(Complex{Float64}, size(cEEtru))
    ons = ones(size(tmp))
    CX = 1./cEEobsfid
    CY = 1./cBBobsfid
    for p = 1:2, q = 1:2
        tmp += -g.k[p] .* g.k[q] .* sin2_conv( g.k[p] .* cEEfid, g.k[q] .* cBBtru, CX, CY, g)
        tmp += -g.k[p] .* g.k[q] .* sin2_conv( g.k[p] .* cEEtru, g.k[q] .* cBBfid, CX, CY, g)
        tmp +=  g.k[p] .* g.k[q] .* sin2_conv( g.k[p] .* g.k[q] .* cEEfid .* cEEtru, ons, CX, CY, g)
        tmp +=  g.k[p] .* g.k[q] .* sin2_conv( ons, g.k[p] .* g.k[q] .* cBBfid .* cBBtru, CX, CY, g)
    end
    tmp .*= 1 / 2 / π
    rtnk  = 1 ./ real(tmp)
    return squash(rtnk)
end

""" Aℓ for EB quadratic estimate"""
function getN0ℓ(cEEfid, cBBfid, cEEobsfid, cBBobsfid, cEEobstru, cBBobstru, cEBobstru, g, Aℓ)
    tmp = zeros(Complex{Float64}, size(cEEfid))
    ons = ones(size(tmp))

    CX = cEEobstru./abs2(cEEobsfid)
    CY = cBBobstru./abs2(cBBobsfid)
    for p = 1:2, q = 1:2
        tmp += -g.k[p] .* g.k[q] .* sin2_conv( g.k[p] .* cEEfid, g.k[q] .* cBBfid, CX, CY, g)
        tmp += -g.k[p] .* g.k[q] .* sin2_conv( g.k[p] .* cEEfid, g.k[q] .* cBBfid, CX, CY, g)
        tmp +=  g.k[p] .* g.k[q] .* sin2_conv( g.k[p] .* g.k[q] .* cEEfid .* cEEfid, ons, CX, CY, g)
        tmp +=  g.k[p] .* g.k[q] .* sin2_conv( ons, g.k[p] .* g.k[q] .* cBBfid .* cBBfid, CX, CY, g)
    end
    CX = cEBobstru ./ cEEobsfid./ cBBobsfid
    CY = cEBobstru ./ cEEobsfid./ cBBobsfid
    for p = 1:2, q = 1:2
        tmp +=  g.k[p] .* g.k[q] .* sin2_conv( g.k[p] .* cEEfid, g.k[q] .* cEEfid, CX, CY, g)
        tmp +=  g.k[p] .* g.k[q] .* sin2_conv( g.k[p] .* cBBfid, g.k[q] .* cBBfid, CX, CY, g)
        tmp += -g.k[p] .* g.k[q] .* sin2_conv( ons, g.k[p] .* g.k[q] .* cEEfid .* cBBfid, CX, CY, g)
        tmp += -g.k[p] .* g.k[q] .* sin2_conv( g.k[p] .* g.k[q] .* cEEfid .* cBBfid,  ons, CX, CY, g)
    end
    tmp .*= 1 / 2 / π
    rtnk  = abs2(Aℓ) .* real(tmp)
    return squash(rtnk)
end

function clenEBk(mCls, g)
    φ_k  = angle(g.k[1] + im * g.k[2])
    cQQk = cos(2φ_k).^2 .* mCls.cEEk + sin(2φ_k).^2 .* mCls.cBBk
    cUUk = sin(2φ_k).^2 .* mCls.cEEk + cos(2φ_k).^2 .* mCls.cBBk
    cQUk = cos(2φ_k) .* sin(2φ_k) .* (mCls.cEEk - mCls.cBBk)
    δ1Qδ1Qk = clenEBk_helper(cQQk, mCls.cϕϕk, g)
    δ1Uδ1Uk = clenEBk_helper(cUUk, mCls.cϕϕk, g)
    δ1Qδ1Uk = clenEBk_helper(cQUk, mCls.cϕϕk, g)
    rtnk   =  cos(2φ_k) .* sin(2φ_k) .* (δ1Uδ1Uk - δ1Qδ1Qk)
    rtnk  +=  (cos(2φ_k).^2 - sin(2φ_k).^2) .* δ1Qδ1Uk
    return real(rtnk)
end
# (cXYk, cϕϕk, g) → ∫dk̲̲ kp * kq * (k+l)p * (k+l)q * cXYk * cϕϕk
function clenEBk_helper(cXYk, cϕϕk, g)
    rtnk = zeros(Complex{Float64}, size(cϕϕk))
    for p = 1:2, q = 1:2
        rtnk += conv( g.k[p].*g.k[q].*cXYk,  g.k[p].*g.k[q].*cϕϕk, g)
    end
    rtnk .*= 1 / 2 / π
    return rtnk
end

function δ01E_δ01B(cEEk, cBBk, ϕk, g)
    φ_k     = angle(g.k[1] + im * g.k[2])
    ex, ek  = sim_xk(cEEk, g)
    bx, bk  = sim_xk(cBBk, g)
    qk, uk, qx, ux  = PrimordialB.eb2qu(ek, bk, g)
    δ0Ek = copy(ek)
    δ0Bk = copy(bk)
    δ1Ek = zeros(Complex{Float64}, size(g.x[1]))
    δ1Bk = zeros(Complex{Float64}, size(g.x[1]))
    for p = 1:2
        δ1Ek +=  - cos(2φ_k).*conv(im.*g.k[p].*qk, im.*g.k[p].*ϕk, g)  - sin(2φ_k).*conv(im.*g.k[p].*uk, im.*g.k[p].*ϕk, g)
        δ1Bk +=    sin(2φ_k).*conv(im.*g.k[p].*qk, im.*g.k[p].*ϕk, g)  - cos(2φ_k).*conv(im.*g.k[p].*uk, im.*g.k[p].*ϕk, g)
    end
    return δ0Ek, δ1Ek, δ0Bk, δ1Bk
end

function simulateEB_N1(itrs, cEEtru, cBBtru, cEEfid, cBBfid, cEEobsfid, cBBobsfid, mCls, g, Aℓ)
    N1k = zeros(Float64, size(g.x[1]))
    for cntr = 1:itrs
        ϕx, ϕk  = sim_xk(mCls.cϕϕk, g)
        δ0Ek_3, δ1Ek_3, δ0Bk_3, δ1Bk_3 = δ01E_δ01B(cEEtru, cBBtru, ϕk, g)
        δ0Ek_4, δ1Ek_4, δ0Bk_4, δ1Bk_4 = δ01E_δ01B(cEEtru, cBBtru, ϕk, g)
        ϕestk_3, = eb_quad_est(δ0Ek_4, δ1Bk_3, cEEfid, cBBfid, cEEobsfid, cBBobsfid, g, Aℓ)
        ϕestk_4, = eb_quad_est(δ0Ek_3, δ1Bk_4, cEEfid, cBBfid, cEEobsfid, cBBobsfid, g, Aℓ)
        N1k += ϕestk_3 .* conj(ϕestk_4) .* (g.deltk^2)
    end
    return real(N1k./ itrs)
end

""" (Xk, Yk, g) → ∫dk̲ X_{k+l} * Y_{-k}"""
function conv(Xk, Yk, g)
    Xx   = (g.FFT \ Xk)
    Yx   = (g.FFT \ Yk)
    rtnk = g.FFT * (Xx .* Yx)
    return  rtnk
end

""" (Xk, Yk, filtk, g) → ∫dk̲ im(k+l)X_{k+l} * im(-k)Y_{-k}  """
function ∂conv(Xk, Yk, g)
    rtnk = zeros(Complex{Float64},size(Xk))
    for p in 1:2
        Xx     = g.FFT \ (im .* g.k[p] .* Xk)
        Yx     = g.FFT \ (im .* g.k[p] .* Yk)
        rtnk .+= g.FFT * (Xx .* Yx)
    end
    return  rtnk
end

""" (Xk, Yk, g) → ∫dk̲ sin(2φ_{k+l} - 2φ_{k}) * (X_{k+l} * CX_{k+l}) * conj(Y_k * CY_{k})
where CX and CY are real """
function sin_conv(Xk, Yk, CXk, CYk, g)
    φ_k     = angle(g.k[1] + im * g.k[2])
    expi2φk = exp(im * 2 * φ_k)
    hatXx   = g.FFT \ (squash(Xk .* CXk) .* expi2φk)
    hatYx   = g.FFT \ (squash(Yk .* CYk) .* expi2φk)
    return  g.FFT * imag(hatXx .* conj(hatYx))
end

""" (Xk, Yk, g) → ∫dk̲ sin²(2φ_{k+l} - 2φ_{k}) * (X_{k+l} * CX_{k+l}) * conj(Y_k * CY_{k})
where CX and CY are real """
function sin2_conv(Xk, Yk, CXk, CYk,  g)
    φ_k     = angle(g.k[1] + im * g.k[2])
    expi4φk = exp(im * 4 * φ_k)
    Xx      = g.FFT \ squash(Xk .* CXk)
    Yx      = g.FFT \ squash(Yk .* CYk)
    checkXx = g.FFT \ (squash(Xk .* CXk) .* expi4φk)
    checkYx = g.FFT \ (squash(Yk .* CYk) .* expi4φk)
    rtn     = g.FFT * (Xx .* Yx - real(checkXx .* conj(checkYx)))
    rtn   .*= 0.5
    return  rtn
end

###############################

# De-lensing by treating observed E as unlensed template (i.e. Blakes method)

###############################

function full_blakes_method(lmin, lmax, dek, dbk, mCls, rng, ρl, g, ϕestk)
    mat_clnbk  = ResB_matpwr(mCls, ρl, g)
    δ1residk   = dbk - estδ1BfromEk(dek, ϕestk, ρl, mCls, g)
    rll        = lminlmax_loglike_marg_profile(rng, δ1residk, lmin, lmax, g, mCls, mat_clnbk)
    return rll
end

# new
function δ1BfromEk(_ek, _ϕk, g ) # was LenB1
	zr = zeros(size(g.r))
    qk, uk,    = PrimordialB.eb2qu(_ek, zr, g)
    δ1qk       = PrimordialB.∂conv(qk, _ϕk, g)
    δ1uk       = PrimordialB.∂conv(uk, _ϕk, g)
    δ1E_fromEk, δ1B_fromEk,   = PrimordialB.qu2eb(δ1qk, δ1uk, g)
    return δ1B_fromEk
end

# new
function estδ1BfromEk(dek, ϕestk, ρl, mCls, g)  # was LenB_cib
    μek  = dek .* squash(mCls.cEEk./(mCls.cEEk + mCls.cEEnoisek))
    μϕk  = ϕestk .* squash(ρl.^2)
    estδ1B_fromEk =  δ1BfromEk(μek, μϕk, g)
    return estδ1B_fromEk
end

function ResB_matpwr(mCls, ρl, g)
    mCls_ee_eff   = mCls.cEEk .* squash(mCls.cEEk ./ (mCls.cEEk + mCls.cEEnoisek))
    mCls_ϕϕ_eff   = mCls.cϕϕk .* squash(ρl.^2)
    mCl_lnbb      = lnB_matpwr(mCls.cEEk, mCls.cϕϕk, g)
    mCl_lnbb_cib  = lnB_matpwr(mCls_ee_eff, mCls_ϕϕ_eff, g)
    mCl_bb_res    = mCl_lnbb - mCl_lnbb_cib
    return mCl_bb_res
end
function lnB_matpwr(mCls_ee, mCls_ϕϕ, g)
    φ_l    = angle(g.k[1] + im * g.k[2])
    cosφ2  = cos(2φ_l)
    sinφ2  = sin(2φ_l)
    cosφ4  = cos(4φ_l)
    sinφ4  = sin(4φ_l)

    unt_t1 = lnB_pwr_helper(g.r.^2 .* mCls_ee,          g.r.^2 .*          mCls_ϕϕ, g, opt=:unit)
    unt_t2 = lnB_pwr_helper(g.r.^2 .* cosφ2 .* mCls_ee, g.r.^2 .* cosφ2 .* mCls_ϕϕ, g, opt=:unit)
    unt_t3 = lnB_pwr_helper(g.r.^2 .* sinφ2 .* mCls_ee, g.r.^2 .* sinφ2 .* mCls_ϕϕ, g, opt=:unit)

    cos_t1 = lnB_pwr_helper(g.r.^2 .* mCls_ee,          g.r.^2 .*          mCls_ϕϕ.* cosφ4, g, opt=:cos).* (-cosφ4)
    cos_t2 = lnB_pwr_helper(g.r.^2 .* cosφ2 .* mCls_ee, g.r.^2 .* cosφ2 .* mCls_ϕϕ.* cosφ4, g, opt=:cos).* (-cosφ4)
    cos_t3 = lnB_pwr_helper(g.r.^2 .* sinφ2 .* mCls_ee, g.r.^2 .* sinφ2 .* mCls_ϕϕ.* cosφ4, g, opt=:cos).* (-cosφ4)
    cos_t4 = lnB_pwr_helper(g.r.^2 .* mCls_ee,          g.r.^2 .*          mCls_ϕϕ.* sinφ4, g, opt=:cos).* (-sinφ4)
    cos_t5 = lnB_pwr_helper(g.r.^2 .* cosφ2 .* mCls_ee, g.r.^2 .* cosφ2 .* mCls_ϕϕ.* sinφ4, g, opt=:cos).* (-sinφ4)
    cos_t6 = lnB_pwr_helper(g.r.^2 .* sinφ2 .* mCls_ee, g.r.^2 .* sinφ2 .* mCls_ϕϕ.* sinφ4, g, opt=:cos).* (-sinφ4)

    sin_t1 = lnB_pwr_helper(g.r.^2 .* mCls_ee,          g.r.^2 .*          mCls_ϕϕ.* sinφ4, g, opt=:sin).* cosφ4
    sin_t2 = lnB_pwr_helper(g.r.^2 .* cosφ2 .* mCls_ee, g.r.^2 .* cosφ2 .* mCls_ϕϕ.* sinφ4, g, opt=:sin).* cosφ4
    sin_t3 = lnB_pwr_helper(g.r.^2 .* sinφ2 .* mCls_ee, g.r.^2 .* sinφ2 .* mCls_ϕϕ.* sinφ4, g, opt=:sin).* cosφ4
    sin_t4 = lnB_pwr_helper(g.r.^2 .* mCls_ee,          g.r.^2 .*          mCls_ϕϕ.* cosφ4, g, opt=:sin).* (-sinφ4)
    sin_t5 = lnB_pwr_helper(g.r.^2 .* cosφ2 .* mCls_ee, g.r.^2 .* cosφ2 .* mCls_ϕϕ.* cosφ4, g, opt=:sin).* (-sinφ4)
    sin_t6 = lnB_pwr_helper(g.r.^2 .* sinφ2 .* mCls_ee, g.r.^2 .* sinφ2 .* mCls_ϕϕ.* cosφ4, g, opt=:sin).* (-sinφ4)

    mCl = 0.25*(unt_t1 + unt_t2 + unt_t3
                +sin_t1 + sin_t2 + sin_t3 + sin_t4 + sin_t5 + sin_t6
                +cos_t1 + cos_t2 + cos_t3 + cos_t4 + cos_t5 + cos_t6)/(2pi)
    return mCl
end
function lnB_pwr_helper(Ak::Array{Float64,2}, Bk::Array{Float64,2}, g; opt = :sin)
    φ4_l    = 4angle(g.k[1] + im * g.k[2])
    Ak_hat  = Ak .* exp(im * φ4_l)
    Bk_hat  = Bk .* exp(im * φ4_l)

    if opt == :unit
        Ax = g.FFT \ Ak
        Bx = g.FFT \ Bk
        return real( g.FFT * (Ax .* conj(Bx)) )
    elseif opt == :sin
        Ax_hat  = g.FFT \ Ak_hat
        Bx_hat  = g.FFT \ Bk_hat
        return real( g.FFT * imag(Ax_hat .* conj(Bx_hat)) )
    elseif opt == :cos
        Ax_hat  = g.FFT \ Ak_hat
        Bx_hat  = g.FFT \ Bk_hat
        return real( g.FFT * real(Ax_hat .* conj(Bx_hat)) )
    else
        println("Wrong Option")
    end
end

##########################################################

# Lensing functions

##############################################################

""" Lense qx, ux:  `rqx, rux = lense(qx, ux, len, g, order=2, qk=g.FFT*qx, uk=g.FFT*ux)` """
function lense{T}(
			qx::Matrix{Float64},
			ux::Matrix{Float64},
			len,
			g::FFTgrid{2,T},
			order::Int64 = 2,
			qk::Matrix{Complex{Float64}} = g.FFT * qx,
			uk::Matrix{Complex{Float64}} = g.FFT * ux,
	)
	rqx, rux  = intlense(qx, ux, len)  # <--- return values
	@inbounds for n in 1:order, α₁ in 0:n
		# kα   = im ^ n .* g.k[1] .^ α₁ .* g.k[2] .^ (n - α₁)
		kα  = intlense_helper1(n, α₁, g.k[1], g.k[2])

		∂α_qx = real(g.FFT \ (kα .* qk))
		∂α_ux = real(g.FFT \ (kα .* uk))
		∂α_qx, ∂α_ux  = intlense(∂α_qx, ∂α_ux, len)

		# xα   = len.rdisplx .^ α₁ .* len.rdisply .^ (n - α₁) ./ factorial(α₁) ./ factorial(n - α₁)
		# rqx += xα .* ∂α_qx
		# rux += xα .* ∂α_ux
		intlense_helper2!(rqx, rux, n, α₁, len.rdisplx, len.rdisply, ∂α_qx, ∂α_ux)
    end
    return rqx, rux
end
function intlense_helper1(n, α₁, k1, k2)
	rtn  = Array(Complex{Float64}, size(k1))
	imn, nmα₁  = im ^ n, n - α₁
	@inbounds @simd for i in eachindex(rtn)
		rtn[i] = complex(imn * k1[i] ^ α₁ * k2[i] ^ nmα₁)
	end
	return rtn
end
function intlense_helper2!(rqx, rux, n, α₁, rx, ry, ∂qx, ∂ux)
	fα₁, fnmα₁, nmα₁ = factorial(α₁), factorial(n - α₁), n - α₁
	@inbounds @simd for i in eachindex(rqx, rux)
		xα      = rx[i] ^ α₁ * ry[i] ^ nmα₁ / fα₁ / fnmα₁
		rqx[i] += xα * ∂qx[i]
		rux[i] += xα * ∂ux[i]
	end
	return nothing
end

function intlense(qx, ux, len)
	rqx  = similar(qx)
	rux  = similar(ux)
    @inbounds for i in eachindex(rqx, rux)
            rqx[i] = qx[len.indrow[i], len.indcol[i]]
            rux[i] = ux[len.indrow[i], len.indcol[i]]
    end
    return rqx, rux
end

##########################################################

# Miscellaneous functions

##############################################################
classy = try
	pyimport(:classy)
end
function class(;ϕscale  = 1.0, ψscale = 0.0, lmax = 6_000, r = 0.005, r0 = 1.0, omega_b = 0.0224567, omega_cdm =0.118489, tau_reio = 0.128312, theta_s = 0.0104098, logA_s_1010 = 3.29056, n_s = 0.968602)
	cosmo = classy[:Class]()
	cosmo[:struct_cleanup]()
	cosmo[:empty]()
	params              = Dict(
	   		"output"        => "tCl, pCl, lCl",
	   		"modes"         => "s,t",
	   		"lensing"       => "yes",
			"l_max_scalars" => lmax + 500,
			"l_max_tensors" => lmax + 500, #lmax + 500,
	      	"omega_b"       => omega_b,
	    	"omega_cdm"     => omega_cdm,
	      	"tau_reio"      => tau_reio,
	      	"100*theta_s"   => 100*theta_s,
	      	"ln10^{10}A_s"  => logA_s_1010,
	      	"n_s"           => n_s,
			"r"             => r,
	        #"k_pivot"      => 0.05,
			#"k_step_trans" => 0.1, # 0.01 for super high resolution
	   		#"l_linstep"    => 10, # 1 for super high resolution
	   		)
	cosmo[:set](params)
	cosmo[:compute]()
	cls_ln              = cosmo[:lensed_cl](lmax)
	cls                 = cosmo[:raw_cl](lmax)
	rtn                 = Dict{Symbol, Array{Float64,1}}(
			:ell        => cls["ell"],
			:ln_tt      => cls_ln["tt"] * (10^6 * cosmo[:T_cmb]()) ^ 2,
			:ln_ee      => cls_ln["ee"] * (10^6 * cosmo[:T_cmb]()) ^ 2,
			:ln_bb      => cls_ln["bb"] * (10^6 * cosmo[:T_cmb]()) ^ 2,
			:ln_te      => cls_ln["te"] * (10^6 * cosmo[:T_cmb]()) ^ 2,
			:ln_tϕ      => cls_ln["tp"] * (10^6 * cosmo[:T_cmb]()),
			:tt 	    => cls["tt"] * (10^6 * cosmo[:T_cmb]()) ^ 2,
			:ee 	    => cls["ee"] * (10^6 * cosmo[:T_cmb]()) ^ 2,
			:bb 	    => cls["bb"] * (10^6 * cosmo[:T_cmb]()) ^ 2,
			:te 	    => cls["te"] * (10^6 * cosmo[:T_cmb]()) ^ 2,
			:tϕ 	    => cls["tp"] * (10^6 * cosmo[:T_cmb]()),
			:ϕϕ         => ϕscale.*cls["pp"],
			:ϕψ         => 0.0.*cls["pp"],
			:ψψ         => ψscale.*cls["pp"],
			:bb0        => cls["bb"] * (10^6 * cosmo[:T_cmb]()) ^ 2 * (r0 / r),
		)
	return rtn
end

""" Convert qu to eb:  `ek, bk, ex, bx = qu2eb(qk, uk, g)` """
function qu2eb(qk, uk, g)
	φ2_l = 2angle(g.k[1] + im * g.k[2])
	ek   = - qk .* cos(φ2_l) - uk .* sin(φ2_l)
	bk   =   qk .* sin(φ2_l) - uk .* cos(φ2_l)
	ex   = real(g.FFT \ ek)
	bx   = real(g.FFT \ bk)
	return ek, bk, ex, bx
end

""" Convert eb to qu: `qk, uk, qx, ux = eb2qu(ek, bk, g)` """
function eb2qu(ek, bk, g)
	φ2_l = 2angle(g.k[1] + im * g.k[2])
	qk   = - ek .* cos(φ2_l) + bk .* sin(φ2_l)
	uk   = - ek .* sin(φ2_l) - bk .* cos(φ2_l)
    qx   = real(g.FFT \ qk)
	ux   = real(g.FFT \ uk)
	return qk, uk, qx, ux
end

function ekbk2_lense_cex_sex_cbx_sbx(ek, bk, len, g, order::Int64 = 7)  # scaled b mode to 1
	cex, sex, cbx, sbx = ekbk2_cex_sex_cbx_sbx(ek, bk, g)
	ln_cex, ln_sex         = PrimordialB.lense(cex, sex, len, g, order)
	ln_sbx, ln_cbx         = PrimordialB.lense(sbx, cbx, len, g, order)
	return ln_cex, ln_sex, ln_cbx, ln_sbx
end

function ekbk2_cex_sex_cbx_sbx(ek, bk, g)
	zr = zeros(size(g.r))
	cek, sek, cex, sex = PrimordialB.eb2qu(ek, zr, g)
	sbk, cbk, sbx, cbx = PrimordialB.eb2qu(zr, bk, g)
	return cex, sex, cbx, sbx
end

function pixelsize_to_muK(pixel_side_length, nside)
    fsky_percent = 100 .* (pixel_side_length .* nside .* π ./ (180*60)) .^2 ./ (4π)
    μKrange      = √(fsky_percent ./ 4.0)
    return μKrange
end

""" kbins, rtnk  = radial_power(fk, smooth, g) """
function radial_power(fk, smooth, g)
	rtnk = Float64[]
	dk = g.deltk
	dm =  length(g.x)
	kbins = collect((smooth*dk):(smooth*dk):(g.nyq))
	for wavenumber in kbins
		indx = (wavenumber-smooth*dk) .< g.r .<= (wavenumber+smooth*dk)
		push!(rtnk, sum(abs2(fk[indx]).* (dk.^dm)) / sum(indx))
	end
	return kbins, rtnk
end

""" kbins, rtnk  = radial_sum(fk, smooth, g) """
function radial_ave(fk, smooth, g)
	rtnk  = Float64[]
	dk    = g.deltk
	kbins = collect((smooth*dk):(smooth*dk):(g.nyq))
	for wavenumber in kbins
		indx = (wavenumber-smooth*dk) .< g.r .<= (wavenumber+smooth*dk)
		push!(rtnk, sum(fk[indx]) / sum(indx))
	end
	return kbins, rtnk
end

# -------- Simulate a mean zero Gaussian random field in the pixel domain given a spectral density.
function sim_xk{dm, T}(cXXk::Array{Float64,dm}, g::FFTgrid{dm, T})
	wx, wk = white_wx_wk(g)
	zk = √(cXXk) .* wk
	zx = real(g.FFT \ zk)
	return zx, zk
end

# ----- white noise
function white_wx_wk{dm, T}(g::FFTgrid{dm, T})
	dx  = g.deltx ^ dm
	wx = randn(size(g.r)) ./ √(dx)
	wk = g.FFT * wx
	return wx, wk
end

squash{T<:Number}(x::T)  = isnan(abs(x)) ? zero(T) : isfinite(abs(x)) ? x : zero(T)
function squash!{dm,T}(x::Array{T,dm}, mask::BitArray{dm}=trues(size(x)))
	@inbounds @simd for i in eachindex(x)
		if !mask[i] || !isfinite(x[i]) || isnan(x[i])
			x[i] = zero(T)
		end
	end
	return x
end
function squash{dm,T}(x::Array{T,dm}, mask::BitArray{dm}=trues(size(x)))
	y = copy(x)
	return squash!(y, mask)
end

end # module


text:
f = PetriCospan([1, 2], Petri.Model(1:4, [(Dict(1=>1), Dict(3=>1)), (Dict(2=>1), Dict(4=>1))]), [3, 4])

g = PetriCospan([1,2], Petri.Model(1:3, [(Dict(1=>1, 2=>1), Dict(3=>1))]), [3])

h = f ⋅ g

h′ = PetriCospan([1,2], Petri.Model(1:5, [(Dict(1=>1), Dict(3=>1)), (Dict(2=>1), Dict(4=>1)), (Dict(3=>1, 4=>1), Dict(5=>1))]), [5])

h_id = h ⋅ id(PetriCospanOb(1))

@test dom(f) == PetriCospanOb(2)
@test codom(f) == PetriCospanOb(2)
@test dom(g) == codom(f)
@test codom(g) == PetriCospanOb(1)
@test dom(h) == dom(f)
@test codom(h) == codom(g)

compare_petricospan(h, h′)
compare_petricospan(h, h_id)

ext: jl


text:
const PCG_LIST = (
    (:Oneseq, UInt16, :XSH_RS, UInt8),
    (:Oneseq, UInt32, :XSH_RS, UInt16),
    (:Oneseq, UInt64, :XSH_RS, UInt32),
    (:Oneseq, UInt128, :XSH_RS, UInt64),
    (:Unique, UInt16, :XSH_RS, UInt8),
    (:Unique, UInt32, :XSH_RS, UInt16),
    (:Unique, UInt64, :XSH_RS, UInt32),
    (:Unique, UInt128, :XSH_RS, UInt64),
    (:Setseq, UInt16, :XSH_RS, UInt8),
    (:Setseq, UInt32, :XSH_RS, UInt16),
    (:Setseq, UInt64, :XSH_RS, UInt32),
    (:Setseq, UInt128, :XSH_RS, UInt64),
    (:MCG, UInt16, :XSH_RS, UInt8),
    (:MCG, UInt32, :XSH_RS, UInt16),
    (:MCG, UInt64, :XSH_RS, UInt32),
    (:MCG, UInt128, :XSH_RS, UInt64),

    (:Oneseq, UInt16, :XSH_RR, UInt8),
    (:Oneseq, UInt32, :XSH_RR, UInt16),
    (:Oneseq, UInt64, :XSH_RR, UInt32),
    (:Oneseq, UInt128, :XSH_RR, UInt64),
    (:Unique, UInt16, :XSH_RR, UInt8),
    (:Unique, UInt32, :XSH_RR, UInt16),
    (:Unique, UInt64, :XSH_RR, UInt32),
    (:Unique, UInt128, :XSH_RR, UInt64),
    (:Setseq, UInt16, :XSH_RR, UInt8),
    (:Setseq, UInt32, :XSH_RR, UInt16),
    (:Setseq, UInt64, :XSH_RR, UInt32),
    (:Setseq, UInt128, :XSH_RR, UInt64),
    (:MCG, UInt16, :XSH_RR, UInt8),
    (:MCG, UInt32, :XSH_RR, UInt16),
    (:MCG, UInt64, :XSH_RR, UInt32),
    (:MCG, UInt128, :XSH_RR, UInt64),

    (:Oneseq, UInt8, :RXS_M_XS, UInt8),
    (:Oneseq, UInt16, :RXS_M_XS, UInt16),
    (:Oneseq, UInt32, :RXS_M_XS, UInt32),
    (:Oneseq, UInt64, :RXS_M_XS, UInt64),
    (:Oneseq, UInt128, :RXS_M_XS, UInt128),
    (:Unique, UInt16, :RXS_M_XS, UInt16),
    (:Unique, UInt32, :RXS_M_XS, UInt32),
    (:Unique, UInt64, :RXS_M_XS, UInt64),
    (:Unique, UInt128, :RXS_M_XS, UInt128),
    (:Setseq, UInt8, :RXS_M_XS, UInt8),
    (:Setseq, UInt16, :RXS_M_XS, UInt16),
    (:Setseq, UInt32, :RXS_M_XS, UInt32),
    (:Setseq, UInt64, :RXS_M_XS, UInt64),
    (:Setseq, UInt128, :RXS_M_XS, UInt128),

    (:Oneseq, UInt64, :XSL_RR, UInt32),
    (:Oneseq, UInt128, :XSL_RR, UInt64),
    (:Unique, UInt64, :XSL_RR, UInt32),
    (:Unique, UInt128, :XSL_RR, UInt64),
    (:Setseq, UInt64, :XSL_RR, UInt32),
    (:Setseq, UInt128, :XSL_RR, UInt64),
    (:MCG, UInt64, :XSL_RR, UInt32),
    (:MCG, UInt128, :XSL_RR, UInt64),

    (:Oneseq, UInt64, :XSL_RR_RR, UInt64),
    (:Oneseq, UInt128, :XSL_RR_RR, UInt128),
    (:Unique, UInt64, :XSL_RR_RR, UInt64),
    (:Unique, UInt128, :XSL_RR_RR, UInt128),
    (:Setseq, UInt64, :XSL_RR_RR, UInt64),
    (:Setseq, UInt128, :XSL_RR_RR, UInt128)
)

let p() = begin
        s = "The list of all the parameter combinations that can be used for PCG.\n\n"
        s *= "|Stream variation|State Type|Method Type|Output Type|\n"
        s *= "|---|---|---|---|\n"
        for (pcg_type, uint_type, method, output_type) in PCG_LIST
            s *= "|`PCGState$pcg_type`|`$uint_type`|`PCG_$method`|`$output_type`|\n"
        end
        s
    end
    @doc p() PCG_LIST
end


text:
#!/usr/bin/env julia

using Luxor

using Random
Random.seed!(42)

using Test

function offset_poly_test(fname)
    Drawing("A2", fname)
    origin()
    background("white")
    setline(1.5)

    tiles = Tiler(Luxor.current_width(), Luxor.current_height(), 6, 6, margin=20)
    randomoffset = 18
    for (pos, n) in tiles
        gsave()
        translate(pos)
        radius =  tiles.tilewidth/4
        p = star(O, radius, 5, 0.25, 0, vertices=true)
        sethue("red")
        setdash("dot")
        poly(p, :stroke, close=true)
        setdash("solid")
        sethue("black")
        randomoffset -= 1
        poly(offsetpoly(p, randomoffset), :stroke, close=true)
        text(string(randomoffset), O.x, O.y + tiles.tilewidth/2, halign=:center)
        grestore()
    end
    @test finish() == true
end

fname="offsetpolys.pdf"
offset_poly_test(fname)
println("...finished test: output in $(fname)")


text:
#### InitParams

export Params

"""
    Params(::Case)

Initialize stand-alone input parameters to
solve the EDMF equations for a given case.
"""
function Params end

function Params(param_set, ::BOMEX)
    FT = Float64
    IT = Int
    params = Dict{Symbol, Any}()

    #####
    ##### TODO: Parameters that need to be added to CLIMAParameters
    #####
    params[:k_Karman] = 0.4 # "Von Karman constant (unit-less)"

    #####
    ##### Filter parameters
    #####
    params[:param_set] = param_set
    params[:N_subdomains] = 3

    #####
    ##### IO parameters
    #####

    params[:plot_single_fields] = true
    params[:export_frequency] = 200

    #####
    ##### Filter parameters
    #####

    params[:a_bounds] = [1e-3, 1 - 1e-3]                                  # filter for a
    params[:w_bounds] = [0.0, 10000.0]                                  # filter for w
    params[:q_bounds] = [0.0, 1.0]                                      # filter for q

    #####
    ##### Space and time discretization
    #####

    params[:UniformGridParams] =
        UniformGridParams{FT, IT}(; z_min = 0.0, z_max = 3000.0, n_elems = 75)
    grid_params = params[:UniformGridParams]
    params[:Δz] = (grid_params.z_max - grid_params.z_min) / grid_params.n_elems

    params[:Δt] = FT(20.0)
    params[:Δt_min] = FT(20.0)
    params[:t_end] = FT(21600.0)
    params[:CFL] = FT(0.8)

    #####
    ##### External conditions
    #####

    # Looks good
    params[:ForcingType] = NoForcing()

    # Looks okay
    # params[:ForcingType] = StandardForcing(apply_subsidence=true,
    #                                        apply_coriolis=true,
    #                                        coriolis_param=FT(0.376e-4))
    #####
    ##### Physical models
    #####

    params[:EntrDetrModel] = BOverW2{FT}(1, 1)
    # Looks okay
    params[:MixingLengthModel] = ConstantMixingLength{FT}(100)
    # Looks okay
    # params[:MixingLengthModel]      = SCAMPyMixingLength{FT}(StabilityDependentParam{FT}(2.7,-100.0),
    #                                                          StabilityDependentParam{FT}(-1.0,-0.2))

    # Getting NaNs for TKE and other fields. Something needs to be fixed
    # params[:MixingLengthModel]      = IgnaciosMixingLength(StabilityDependentParam{FT}(2.7,-100.0),
    #                                                        StabilityDependentParam{FT}(-1.0,-0.2),
    #                                                        0.1, 0.12, 0.4, 40/13, 0.74)

    params[:EddyDiffusivityModel] = SCAMPyEddyDiffusivity{FT}(0.1)

    params[:PressureModel] = SCAMPyPressure{FT}(;
        buoy_coeff = FT(1.0 / 3.0),
        drag_coeff = FT(0.375),
        plume_spacing = FT(500.0),
    )
    params[:prandtl_number] = FT(1.0)
    params[:tke_diss_coeff] = FT(2.0)

    params[:f_coriolis] = 0.376e-4                                      # coriolis force coefficient
    params[:f_c] = 0                                                    # buoyancy gradient factor

    params[:zrough] = 1.0e-4                                            # surface roughness

    params[:SurfaceModel] = SurfaceFixedFlux{FT}(
        param_set;
        T = FT(300.4),
        P = FT(1.015e5),
        q_tot = FT(0.02245),
        ustar = FT(0.28),
        windspeed_min = FT(0.0),
        tke_tol = FT(0.01),
        area = FT(0.1),
    )

    params[:inversion_height] = [1.0 for i in 1:params[:N_subdomains]]  # inversion height
    params[:Ri_bulk_crit] = 0.0                                         # inversion height parameters
    _molmass_ratio::FT = FT(molmass_ratio(param_set))
    _grav::FT = FT(grav(param_set))
    params[:bflux] = (
        _grav * (
            (
                8.0e-3 +
                (_molmass_ratio - 1) * (299.1 * 5.2e-5 + 22.45e-3 * 8.0e-3)
            ) / (299.1 * (1.0 + (_molmass_ratio - 1) * 22.45e-3))
        )
    )
    params[:cq] = 0.001133                                              # Some surface parameter in SCAMPy
    params[:ch] = 0.001094                                              # Some surface parameter in SCAMPy
    params[:cm] = 0.001229                                              # Some surface parameter in SCAMPy
    params[:grid_adjust] =
        (log(20.0 / params[:zrough]) / log(params[:Δz] / 2 / params[:zrough]))^2 # Some surface parameter in SCAMPy
    params[:cq] *= params[:grid_adjust]                                 # Some surface parameter in SCAMPy
    params[:ch] *= params[:grid_adjust]                                 # Some surface parameter in SCAMPy
    params[:cm] *= params[:grid_adjust]                                 # Some surface parameter in SCAMPy

    return params
end


text:
"""
    link!(p, q [, length])
    link!(p, q, branch)

Connect node `p` to node `q` with a new branch. The new branch will have no length if `length` is not given. Alternatively, an existing `Branch` object can be passed instead of creating a new one.

Returns the branch.
"""
function link!(p::AbstractNode, q::AbstractNode, br::AbstractBranch)
    p ≡ q && throw(InvalidTopology("`p` and `q` are the same node"))
    link_p2q = NodeLink(q, br)
    link_q2p = NodeLink(p, br)

    push!(p.links, link_p2q)
    push!(q.links, link_q2p)

    return br
end # function link!

link!(p::AbstractNode, q::AbstractNode) = link!(p, q, Branch())

link!(p::AbstractNode, q::AbstractNode, length::Float64) =
    link!(p, q, Branch(length))

"""
    unlink!(p, q)

Dissociate the nodes `p` and `q` and return the branch that was connecting them.
"""
function unlink!(p::AbstractNode, q::AbstractNode)
    for (i, link) in enumerate(p.links)
        if link.to ≡ q
            for (j, link) in enumerate(q.links)
                if link.to ≡ p
                    @inbounds br = p.links[i].branch
                    deleteat!(p.links, i)
                    deleteat!(q.links, j)

                    return br
                end
            end
        end
    end

    throw(InvalidTopology("the two nodes are not linked."))
end # function unlink!

"""
    graft!(p, q, r; prop_brlength=0.5)
    graft!(p, q, r, branch_pq, branch_qr)

Insert node `p` between the nodes `q` and `r`, where `q` and `r` are neighbours.

The length of the branch between `q` and `r` is redistributed between the new branches (between `p` and `q` and between `q` and `r`) according to the `prop_brlength` parameter, which gives the proportion of the length of the branch between `q` and `r` that is to be assigned to the branch between `p` and `q`. The branch between `q` and `r` is assigned the remainder of the original branch length.

Alternatively, branch objects can be passed to be used as the new branches.
"""
function graft!(
    p::AbstractNode,
    q::AbstractNode,
    r::AbstractNode,
    br_qp::AbstractBranch,
    br_pr::AbstractBranch
    )

    unlink!(q, r)
    link!(q, p, br_qp)
    link!(p, r, br_pr)

    return nothing
end # function graft!

function graft!(
    p::AbstractNode,
    q::AbstractNode,
    r::AbstractNode,
    prop_qp::Float64 = 0.5;
    f::Union{Function, DataType} = Branch
    )

    @assert 0.0 ≤ prop_qp ≤ 1.0
    br0 = unlink!(q, r)
    if br0.length ≠ nothing
        link!(q, p, f(br0.length * prop_qp))
        link!(p, r, f(br0.length * (1.0 - prop_qp)))
    else
        link!(q, p, f())
        link!(p, r, f())
    end

    return nothing
end # function graft!

"""
    pluck!(p, q, r [, constructor])

Disconnect node `p` from its neighbours `q` and `r`, and connect `q` and `r` to each other. An
"""
function pluck!(
    p::AbstractNode,
    q::AbstractNode,
    r::AbstractNode,
    f::Union{Function, DataType} = Branch
    )
    br_pq = unlink!(p, q)
    br_pr = unlink!(p, r)

    if isnothing(br_pq.length) || isnothing(br_pr.length)
        br_qr = link!(q, r)
    else
        br_qr = link!(q, r, f(br_pq.length + br_pr.length))
    end

    sizehint!(br_qr.datablocks, length(br_pq.datablocks))

    return nothing
end # function pluck!


text:
#######################
###### Variables ######
#######################

# Variable names must begin with a letter (A-Z or a-z),
# underscore, or a subset of Unicode code points greater than 00A0

# Assign the value 10 to the variable x
x = 10
# println(x)

# Doing math with x's value
x = x + 1
# println(x)

# Reassign x's value
x = 1 + 1
# println(x)

# You can assign values of other types, like strings of text
x = "Hello World!"
# println(x)

# Variable names are case-sensitive, and have no semantic meaning
δ = 0.00001
# println(δ)

안녕하세요 = "Hello"
# println(안녕하세요)

# You can type many Unicode math symbols by typing the backslashed LaTeX symbol name followed by tab.
# For example, the variable name δ can be entered by typing \delta-tab,
# or even α̂⁽²⁾ by \alpha-tab-\hat- tab-\^(2)-tab.

###################################
###### Stylistic Conventions ######
###################################

# 1) Names of variables are in lower case.
# 2) Word separation can be indicated by underscores ('_'),
#    but use of underscores is discouraged unless the name would be hard to read otherwise.
# 3) Names of Types and Modules begin with a capital letter and word separation is shown
#    with upper camel case instead of underscores.
# 4) Names of functions and macros are in lower case, without underscores.
# 5) Functions that write to their arguments have names that end in !.
#    These are sometimes called "mutating" or "in-place" functions because they are intended
#    to produce changes in their arguments after the function is called, not just return a value.

max_stars_repo_name: mauricioriva/learn-julia

max_issues_repo_name: mauricioriva/learn-julia

max_forks_repo_name: mauricioriva/learn-julia


text:
type EnergiesT
  Total::Float64
  Kinetic::Float64
  Ionic::Float64
  Hartree::Float64
  XC::Float64
end

function print_Energies( Energies::EnergiesT )
  @printf("Total   energy: %18.10f\n", Energies.Total )
  @printf("Kinetic energy: %18.10f\n", Energies.Kinetic )
  @printf("Ionic   energy: %18.10f\n", Energies.Ionic )
  @printf("Hartree energy: %18.10f\n", Energies.Hartree )
  @printf("XC      energy: %18.10f\n", Energies.XC )
end


text:
using QuBase.bases

m = [1+1im 2+2im 3+3im 4+4im;
     5+5im 6+6im 7+7im 8+8im;
     9+9im 10+10im 11+11im 12+12im;
     13+13im 14+14im 15+15im 16+16im]

v = [1+1im, 2+2im, 3+3im, 4+4im]

qm = QuArray(m, FiniteBasis(4,1), FiniteBasis(2,2))
qv = QuArray(v)

# operator ⊗ operator
qm_qm = tensor(qm, qm)
@assert rawcoeffs(qm_qm) == kron(m, m)
@assert bases(qm_qm) == (tensor(bases(qm, 1), bases(qm, 1)), tensor(bases(qm, 2), bases(qm, 2)))

qmc_qm = tensor(qm', qm)
@assert rawcoeffs(qmc_qm) == kron(m', m)
@assert bases(qmc_qm) == (tensor(bases(qm, 2), bases(qm, 1)), tensor(bases(qm, 1), bases(qm, 2)))

qm_qmc = tensor(qm, qm')
@assert rawcoeffs(qm_qmc) == kron(m, m')
@assert bases(qm_qmc) == (tensor(bases(qm, 1), bases(qm, 2)), tensor(bases(qm, 2), bases(qm, 1)))

qmc_qmc = tensor(qm', qm')
@assert rawcoeffs(qmc_qmc) == kron(m, m)
@assert bases(qmc_qmc) == (tensor(bases(qm, 2), bases(qm, 2)), tensor(bases(qm, 1), bases(qm, 1)))

# state ⊗ state
qv_qv = tensor(qv, qv)
@assert rawcoeffs(qv_qv) == kron(v, v)
@assert bases(qv_qv) == (tensor(bases(qv, 1), bases(qv, 1)),)

qv_qvc = tensor(qv, qv')
@assert rawcoeffs(qv_qvc) == kron(v, v')
@assert bases(qv_qvc) == (bases(qv, 1), bases(qv, 1))

qvc_qv = tensor(qv', qv)
@assert rawcoeffs(qvc_qv) == kron(v', v)
@assert bases(qvc_qv) == (bases(qv, 1), bases(qv, 1))

qvc_qvc = tensor(qv', qv')
@assert rawcoeffs(qvc_qvc) == kron(v, v)
@assert bases(qvc_qvc) == (tensor(bases(qv, 1), bases(qv, 1)),)

# n-ary version of tensor
a = QuArray(rand(Complex{Float64},4), FiniteBasis(2,2))
b = QuArray(rand(Complex{Float64},1), FiniteBasis(1,1))
c = QuArray(rand(Complex{Float64},3))

abc = tensor(a,b,c)
ab_c = tensor(tensor(a, b), c)
a_bc = tensor(a, tensor(b, c))

@assert bases(abc) == bases(ab_c)
@assert bases(abc) == bases(a_bc)
@test_approx_eq coeffs(abc) coeffs(ab_c)
@test_approx_eq coeffs(abc) coeffs(a_bc)


text:
#int64ops.jl
#collates all of the relevant int64 operations

#bring in some important UInt64 bitwise methods
include("./i64o-helpers.jl")
include("./i64o-constants.jl")      #functions generating int64 constants
include("./i64o-bitpatterns.jl")    #boolean functions reporting on bit patterns
include("./i64o-clzctz.jl")         #leading_zeros and trailing_zeros methods.
include("./i64o-masks.jl")          #mask generating operators
include("./i64o-shifts.jl")         #left and right shift operators
include("./i64o-comparison.jl")     #chaining comparison operators
include("./i64o-utilities.jl")      #other int64 utilities that help unums.
include("./i64o-math.jl")

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
"""
    static_analysis(assembly; kwargs...)

Perform a static analysis of the system of nonlinear beams contained in
`assembly`. Return the resulting system and a flag indicating whether the
iteration procedure converged.

# Keyword Arguments
 - `prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}()`: Dictionary
     holding PrescribedConditions composite types for the points in `keys(prescribed_conditions)`
 - `distributed_loads = Dict{Int,DistributedLoads{Float64}}()`: Dictionary holding
     DistributedLoads composite types for the beam elements in `keys(distributed_loads)`
 - `linear = false`: Set to `true` for a linear analysis
 - `method = :newton`: Method (as defined in NLsolve) to solve nonlinear system of equations
 - `linesearch = LineSearches.BackTracking()`: Line search used to solve nonlinear system of equations
 - `ftol = 1e-9`: tolerance for solving nonlinear system of equations
 - `iterations = 1000`: maximum iterations for solving the nonlinear system of equations
 - `nstep = 1`: Number of time steps. May be used in conjunction with time varying
     prescribed conditions and distributed loads to gradually increase
    displacements/loads.
"""
function static_analysis(assembly;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads = Dict{Int,DistributedLoads{Float64}}(),
    linear = false,
    method = :newton,
    linesearch = BackTracking(),
    ftol = 1e-9,
    iterations = 1000,
    nstep = 1,
    )

    static = true

    system = System(assembly, keys(prescribed_conditions), static)

    return static_analysis!(system, assembly;
        prescribed_conditions = prescribed_conditions,
        distributed_loads = distributed_loads,
        linear = linear,
        method = method,
        linesearch = linesearch,
        ftol = ftol,
        iterations = iterations,
        nstep = nstep)
end

"""
    static_analysis!(system, assembly; kwargs...)

Pre-allocated version of `static_analysis`.  Uses the state variables stored in
`system` as an initial guess for iterating.
"""
function static_analysis!(system, assembly;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads = Dict{Int,DistributedLoads{Float64}}(),
    linear = false,
    method = :newton,
    linesearch = BackTracking(),
    ftol = 1e-9,
    iterations = 1000,
    nstep = 1)

    # check to make sure system is static
    @assert system.static == true

    # unpack pre-allocated storage and pointers
    x = system.x
    F = system.r
    J = system.K
    irow_pt = system.irow_pt
    irow_beam = system.irow_beam
    irow_beam1 = system.irow_beam1
    irow_beam2 = system.irow_beam2
    icol_pt = system.icol_pt
    icol_beam = system.icol_beam
    current_step_ref = system.current_step

    converged = true
    for istep = 1:nstep

        # update current time step
        current_step_ref[] = istep

        # solve the system of equations
        f! = (F, x) -> system_residual!(F, x, assembly, prescribed_conditions,
            distributed_loads, istep, irow_pt, irow_beam, irow_beam1,
            irow_beam2, icol_pt, icol_beam)

        j! = (J, x) -> system_jacobian!(J, x, assembly, prescribed_conditions,
            distributed_loads, istep, irow_pt, irow_beam, irow_beam1, irow_beam2,
            icol_pt, icol_beam)

        if linear
            # linear analysis
            x .= 0.0
            f!(F, x)
            j!(J, x)
            x .-= J\F
        else
            # nonlinear analysis
            df = NLsolve.OnceDifferentiable(f!, j!, x, F, J)

            result = NLsolve.nlsolve(df, x,
                linsolve=(x, A, b) -> ldiv!(x, lu(A), b),
                method=method,
                linesearch=linesearch,
                ftol=ftol,
                iterations=iterations)

            # update the solution
            x .= result.zero

            # update convergence flag
            converged = result.f_converged
        end
    end

    return system, converged
end

"""
    steady_state_analysis(assembly; kwargs...)

Perform a steady-state analysis for the system of nonlinear beams contained in
`assembly`.  Return the resulting system and a flag indicating whether the
iteration procedure converged.

# Keyword Arguments
 - `prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}()`: Dictionary
     holding PrescribedConditions composite types for the points in `keys(prescribed_conditions)`
 - `distributed_loads = Dict{Int,DistributedLoads{Float64}}()`: Dictionary holding
     DistributedLoads composite types for the beam elements in `keys(distributed_loads)`
 - `linear = false`: Set to `true` for a linear analysis
 - `method = :newton`: Method (as defined in NLsolve) to solve nonlinear system of equations
 - `linesearch = LineSearches.BackTracking()`: Line search used to solve nonlinear system of equations
 - `ftol = 1e-9`: tolerance for solving nonlinear system of equations
 - `iterations = 1000`: maximum iterations for solving the nonlinear system of equations
 - `nstep = 1`: Number of time steps. May be used in conjunction with time varying
     prescribed conditions, distributed loads, and global motion to gradually
     increase displacements/loads.
 - `origin = zeros(3)`: Global frame origin
 - `linear_velocity = fill(zeros(3), nstep)`: Global frame linear velocity for each time step.
 - `angular_velocity = fill(zeros(3), nstep)`: Global frame angular velocity for each time step
"""
function steady_state_analysis(assembly;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads = Dict{Int,DistributedLoads{Float64}}(),
    linear = false,
    method = :newton,
    linesearch = BackTracking(),
    ftol = 1e-9,
    iterations = 1000,
    nstep = 1,
    origin = (@SVector zeros(3)),
    linear_velocity = fill((@SVector zeros(3)), nstep),
    angular_velocity = fill((@SVector zeros(3)), nstep),
    )

    static = false

    system = System(assembly, keys(prescribed_conditions), static)

    return steady_state_analysis!(system, assembly;
        prescribed_conditions = prescribed_conditions,
        distributed_loads = distributed_loads,
        linear = linear,
        method = method,
        linesearch = linesearch,
        ftol = ftol,
        iterations = iterations,
        nstep = nstep,
        origin = origin,
        linear_velocity = linear_velocity,
        angular_velocity = angular_velocity
        )
end

"""
    steady_state_analysis!(system, assembly; kwargs...)

Pre-allocated version of `steady_state_analysis`.  Uses the state variables stored in
`system` as an initial guess for iterating.
"""
function steady_state_analysis!(system, assembly;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads = Dict{Int,DistributedLoads{Float64}}(),
    linear = false,
    method = :newton,
    linesearch = BackTracking(),
    ftol = 1e-9,
    iterations = 1000,
    nstep = 1,
    origin = (@SVector zeros(3)),
    linear_velocity = fill((@SVector zeros(3)), nstep),
    angular_velocity = fill((@SVector zeros(3)), nstep),
    )

    # check to make sure the simulation is dynamic
    @assert system.static == false

    # unpack pointers to pre-allocated storage
    x = system.x
    F = system.r
    J = system.K
    irow_pt = system.irow_pt
    irow_beam = system.irow_beam
    irow_beam1 = system.irow_beam1
    irow_beam2 = system.irow_beam2
    icol_pt = system.icol_pt
    icol_beam = system.icol_beam
    current_step_ref = system.current_step

    # process global frame motion inputs (in case only one value was provided)
    x0 = SVector{3}(origin)
    if eltype(linear_velocity) <: Number
        v0 = fill(SVector{3}(linear_velocity), nstep)
    else
        v0 = SVector{3}.(linear_velocity)
    end
    if eltype(angular_velocity) <: Number
        ω0 = fill(SVector{3}(angular_velocity), nstep)
    else
        ω0 = SVector{3}.(angular_velocity)
    end

    converged = true
    for istep = 1:nstep

        # update current time step
        current_step_ref[] = istep

        f! = (F, x) -> system_residual!(F, x, assembly, prescribed_conditions,
            distributed_loads, istep, irow_pt, irow_beam, irow_beam1, irow_beam2,
            icol_pt, icol_beam, x0, v0[istep], ω0[istep])

        j! = (J, x) -> system_jacobian!(J, x, assembly, prescribed_conditions,
            distributed_loads, istep, irow_pt, irow_beam, irow_beam1, irow_beam2,
            icol_pt, icol_beam, x0, v0[istep], ω0[istep])

        # solve the system of equations
        if linear
            # linear analysis
            x .= 0.0
            f!(F, x)
            j!(J, x)
            x .-= J\F
        else
            # nonlinear analysis
            df = NLsolve.OnceDifferentiable(f!, j!, x, F, J)

            result = NLsolve.nlsolve(df, x,
                linsolve=(x, A, b) -> ldiv!(x, lu(A), b),
                method=method,
                linesearch=linesearch,
                ftol=ftol,
                iterations=iterations)

            # update the solution
            x .= result.zero

            # update the convergence flag
            convergence = result.f_converged
        end
    end

    return system, converged
end

"""
    eigenvalue_analysis(assembly; kwargs...)

Compute the eigenvalues and eigenvectors of the system of nonlinear beams
contained in `assembly`.  Return the modified system, eigenvalues, eigenvectors,
and a convergence flag indicating whether the corresponding steady-state analysis
converged.

# Keyword Arguments
 - `prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}()`: Dictionary
     holding PrescribedConditions composite types for the points in `keys(prescribed_conditions)`
 - `distributed_loads = Dict{Int,DistributedLoads{Float64}}()`: Dictionary holding
     DistributedLoads composite types for the beam elements in `keys(distributed_loads)`
 - `linear = false`: Set to `true` for a linear analysis
 - `method = :newton`: Method (as defined in NLsolve) to solve nonlinear system of equations
 - `linesearch = LineSearches.BackTracking()`: Line search used to solve nonlinear system of equations
 - `ftol = 1e-9`: tolerance for solving nonlinear system of equations
 - `iterations = 1000`: maximum iterations for solving the nonlinear system of equations
 - `nstep = 1`: Number of time steps. May be used in conjunction with time varying
     prescribed conditions, distributed loads, and global motion to gradually
     increase displacements/loads.
 - `origin = zeros(3)`: Global frame origin
 - `linear_velocity = fill(zeros(3), nstep)`: Global frame linear velocity for each time step.
 - `angular_velocity = fill(zeros(3), nstep)`: Global frame angular velocity for each time step
 - `nev = 6`: Number of eigenvalues to compute
"""
function eigenvalue_analysis(assembly;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads = Dict{Int,DistributedLoads{Float64}}(),
    method = :newton,
    linear = false,
    linesearch = BackTracking(),
    ftol = 1e-9,
    iterations = 1000,
    nstep = 1,
    origin = (@SVector zeros(3)),
    linear_velocity = fill((@SVector zeros(3)), nstep),
    angular_velocity = fill((@SVector zeros(3)), nstep),
    nev = 6
    )

    static = false

    system = System(assembly, keys(prescribed_conditions), static)

    return eigenvalue_analysis!(system, assembly;
        prescribed_conditions = prescribed_conditions,
        distributed_loads = distributed_loads,
        linear = linear,
        method = method,
        linesearch = linesearch,
        ftol = ftol,
        iterations = iterations,
        nstep = nstep,
        origin = origin,
        linear_velocity = linear_velocity,
        angular_velocity = angular_velocity,
        nev = nev
        )
end

"""
    eigenvalue_analysis!(system, assembly; kwargs...)

Pre-allocated version of `eigenvalue_analysis`.  Uses the state variables stored in
`system` as an initial guess for iterating to find the steady state solution.
"""
function eigenvalue_analysis!(system, assembly;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads = Dict{Int,DistributedLoads{Float64}}(),
    linear = false,
    method = :newton,
    linesearch = BackTracking(),
    ftol = 1e-9,
    iterations = 1000,
    nstep = 1,
    origin = (@SVector zeros(3)),
    linear_velocity = fill((@SVector zeros(3)), nstep),
    angular_velocity = fill((@SVector zeros(3)), nstep),
    nev = 6,
    )

    # perform steady state analysis (if nonlinear)
    if linear
        system.x .= 0.0
        converged = true
    else
        system, converged = steady_state_analysis!(system, assembly;
            prescribed_conditions = prescribed_conditions,
            distributed_loads = distributed_loads,
            linear = linear,
            method = method,
            linesearch = linesearch,
            ftol = ftol,
            iterations = iterations,
            nstep = nstep,
            origin = origin,
            linear_velocity = linear_velocity,
            angular_velocity = angular_velocity
            )
    end

    # unpack state vector, stiffness, and mass matrices
    x = system.x # populated during steady state solution
    K = system.K # needs to be updated
    M = system.M # still needs to be populated

    # also unpack system indices
    irow_pt = system.irow_pt
    irow_beam = system.irow_beam
    irow_beam1 = system.irow_beam1
    irow_beam2 = system.irow_beam2
    icol_pt = system.icol_pt
    icol_beam = system.icol_beam

    # get current time step
    istep = system.current_step[]

    # process global frame motion inputs (in case only one value was provided)
    x0 = SVector{3}(origin)
    if eltype(linear_velocity) <: Number
        v0 = SVector{3}(linear_velocity)
    else
        v0 = SVector{3}(linear_velocity[istep])
    end
    if eltype(angular_velocity) <: Number
        ω0 = SVector{3}(angular_velocity)
    else
        ω0 = SVector{3}(angular_velocity[istep])
    end

    # solve for the system stiffness matrix
    K = system_jacobian!(K, x, assembly, prescribed_conditions,
        distributed_loads, istep, irow_pt, irow_beam,
        irow_beam1, irow_beam2, icol_pt, icol_beam, x0, v0, ω0)

    # solve for the system mass matrix
    M = system_mass_matrix!(M, x, assembly, irow_pt, irow_beam, irow_beam1,
        irow_beam2, icol_pt, icol_beam)

    # construct linear map
    T = eltype(system)
    nx = length(x)
    Kfact = lu(K)
    f! = (b, x) -> ldiv!(b, Kfact, M*x)
    fc! = (b, x) -> mul!(b, M', Kfact'\x)
    A = LinearMap{T}(f!, fc!, nx, nx; ismutating=true)

    # compute eigenvalues and eigenvectors
    λ, V = partialeigen(partialschur(A; nev=min(nx,nev), which=LM())[1])

    # sort eigenvalues by magnitude
    perm = sortperm(λ, by=(λ)->(abs(λ),imag(λ)), rev=true)
    λ .= λ[perm]
    V .= V[:,perm]

    # eigenvalues are actually 1/λ, no modification necessary for eigenvectors
    λ .= 1 ./ λ

    return system, λ, V, converged
end

"""
    time_domain_analysis(assembly, dt; kwargs...)

Perform a time-domain analysis for the system of nonlinear beams contained in
`assembly`.  Return the final system, a post-processed solution history, and a
convergence flag indicating whether the iterations converged for each time step.

# Keyword Arguments
 - `prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}()`: Dictionary
     holding PrescribedConditions composite types for the points in `keys(prescribed_conditions)`
 - `distributed_loads = Dict{Int,DistributedLoads{Float64}}()`: Dictionary holding
     DistributedLoads composite types for the beam elements in `keys(distributed_loads)`
 - `linear = false`: Set to `true` for a linear analysis
 - `method = :newton`: Method (as defined in NLsolve) to solve nonlinear system of equations
 - `linesearch = LineSearches.BackTracking()`: Line search used to solve nonlinear system of equations
 - `ftol = 1e-9`: tolerance for solving nonlinear system of equations
 - `iterations = 1000`: maximum iterations for solving the nonlinear system of equations
 - `nstep = 1`: The total length of the time vector
 - `origin = zeros(3)`: Global frame origin
 - `linear_velocity = fill(zeros(3), nstep)`: Global frame linear velocity for each time step.
 - `angular_velocity = fill(zeros(3), nstep)`: Global frame angular velocity for each time step.
 - `u0=fill(zeros(3), length(assembly.elements))`: initial displacment of each beam element,
 - `theta0=fill(zeros(3), length(assembly.elements))`: initial angular displacement of each beam element,
 - `udot0=fill(zeros(3), length(assembly.elements))`: initial time derivative with respect to `u`
 - `thetadot0=fill(zeros(3), length(assembly.elements))`: initial time derivative with respect to `theta`
 - `save=1:nstep`: Steps at which to save the time history
"""
function time_domain_analysis(assembly, dt;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads = Dict{Int,DistributedLoads{Float64}}(),
    linear = false,
    method = :newton,
    linesearch = BackTracking(),
    ftol = 1e-9,
    iterations = 1000,
    nstep = 1,
    origin = (@SVector zeros(3)),
    linear_velocity = fill((@SVector zeros(3)), nstep),
    angular_velocity = fill((@SVector zeros(3)), nstep),
    u0 = fill((@SVector zeros(3)), length(assembly.elements)),
    theta0 = fill((@SVector zeros(3)), length(assembly.elements)),
    udot0 = fill((@SVector zeros(3)), length(assembly.elements)),
    thetadot0 = fill((@SVector zeros(3)), length(assembly.elements)),
    save = 1:nstep
    )

    static = false

    system = System(assembly, keys(prescribed_conditions), static)

    return time_domain_analysis!(system, assembly, dt;
        prescribed_conditions = prescribed_conditions,
        distributed_loads = distributed_loads,
        linear = linear,
        method = method,
        linesearch = linesearch,
        ftol = ftol,
        iterations = iterations,
        nstep = nstep,
        origin = origin,
        linear_velocity = linear_velocity,
        angular_velocity = angular_velocity,
        u0 = u0,
        theta0 = theta0,
        udot0 = udot0,
        thetadot0 = thetadot0,
        save = 1:nstep
        )
end

"""
    time_domain_analysis!(system, assembly, dt; kwargs...)

Pre-allocated version of `time_domain_analysis`.  Uses the state variables stored
in `system` as an initial guess for iterating.
"""
function time_domain_analysis!(system, assembly, dt;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads = Dict{Int,DistributedLoads{Float64}}(),
    linear = false,
    method = :newton,
    linesearch = BackTracking(),
    ftol = 1e-9,
    iterations = 1000,
    nstep = 1,
    origin = (@SVector zeros(3)),
    linear_velocity = fill((@SVector zeros(3)), nstep),
    angular_velocity = fill((@SVector zeros(3)), nstep),
    u0 = fill((@SVector zeros(3)), length(assembly.elements)),
    theta0 = fill((@SVector zeros(3)), length(assembly.elements)),
    udot0 = fill((@SVector zeros(3)), length(assembly.elements)),
    thetadot0 = fill((@SVector zeros(3)), length(assembly.elements)),
    save = 1:nstep
    )

    # check to make sure the simulation is dynamic
    @assert system.static == false

    # unpack pre-allocated storage and pointers for system
    x = system.x
    F = system.r
    J = system.K
    irow_pt = system.irow_pt
    irow_beam = system.irow_beam
    irow_beam1 = system.irow_beam1
    irow_beam2 = system.irow_beam2
    icol_pt = system.icol_pt
    icol_beam = system.icol_beam
    udot_init = system.udot_init
    θdot_init = system.θdot_init
    CtCabPdot_init = system.CtCabPdot_init
    CtCabHdot_init = system.CtCabHdot_init
    current_step_ref = system.current_step

    nbeam = length(assembly.elements)

    # process global frame motion inputs (in case only one value was provided)
    x0 = SVector{3}(origin)
    if eltype(linear_velocity) <: Number
        v0 = fill(SVector{3}(linear_velocity), nstep)
    else
        v0 = SVector{3}.(linear_velocity)
    end
    if eltype(angular_velocity) <: Number
        ω0 = fill(SVector{3}(angular_velocity), nstep)
    else
        ω0 = SVector{3}.(angular_velocity)
    end

    # --- Initial Condition Run --- #

    # update current time step
    current_step_ref[] = 1

    # construct residual and jacobian functions
    f! = (F, x) -> system_residual!(F, x, assembly, prescribed_conditions,
        distributed_loads, 1, irow_pt, irow_beam, irow_beam1, irow_beam2,
        icol_pt, icol_beam, x0, v0[1], ω0[1], u0, theta0, udot0, thetadot0)

    j! = (J, x) -> system_jacobian!(J, x, assembly, prescribed_conditions,
        distributed_loads, 1, irow_pt, irow_beam, irow_beam1, irow_beam2,
        icol_pt, icol_beam, x0, v0[1], ω0[1], u0, theta0, udot0, thetadot0)

    # solve system of equations
    if linear
        # linear analysis
        x .= 0.0
        f!(F, x)
        j!(J, x)
        x .-= J\F
    else
        # nonlinear analysis
        df = OnceDifferentiable(f!, j!, x, F, J)

        result = NLsolve.nlsolve(df, x,
            linsolve=(x, A, b) -> ldiv!(x, lu(A), b),
            method=method,
            linesearch=linesearch,
            ftol=ftol,
            iterations=iterations)

        x .= result.zero
    end

    converged = result.f_converged

    # --- End Initial Condition Run --- #

    # now set up for the time-domain run
    for ibeam = 1:nbeam
        icol = icol_beam[ibeam]
        # calculate udot_init
        udot_init[ibeam] = 2/dt*u0[ibeam] + udot0[ibeam]
        # calculate θdot_init
        θdot_init[ibeam] = 2/dt*theta0[ibeam] + thetadot0[ibeam]
        # extract rotation parameters
        C = get_C(theta0[ibeam])
        Cab = assembly.elements[ibeam].Cab
        CtCab = C'*Cab
        # calculate CtCabPdot_init
        CtCabP = CtCab*SVector{3}(x[icol+12], x[icol+13], x[icol+14])
        CtCabPdot = SVector{3}(x[icol], x[icol+1], x[icol+2])
        CtCabPdot_init[ibeam] = 2/dt*CtCabP + CtCabPdot
        # calculate CtCabHdot_init
        CtCabH = CtCab*SVector{3}(x[icol+15], x[icol+16], x[icol+17])
        CtCabHdot = SVector{3}(x[icol+3], x[icol+4], x[icol+5])
        CtCabHdot_init[ibeam] = 2/dt*CtCabH + CtCabHdot
        # insert initial conditions for time-domain analysis
        x[icol:icol+2] .= u0[ibeam]
        x[icol+3:icol+5] .= theta0[ibeam]
    end

    # initialize storage for each time step
    isave = 1
    history = Vector{AssemblyState{eltype(system)}}(undef, length(save))

    # add initial state to the solution history (if it should be saved)
    if 1 in save
        history[isave] = AssemblyState(system, assembly, prescribed_conditions = prescribed_conditions)
        isave += 1
    end

    # --- Begin Time Domain Simulation --- #

    for istep = 2:nstep

        # update current time step
        current_step_ref[] = istep

        # solve for the state variables at the next time step
        f! = (F, x) -> system_residual!(F, x, assembly, prescribed_conditions,
            distributed_loads, istep, irow_pt, irow_beam, irow_beam1, irow_beam2,
            icol_pt, icol_beam, x0, v0[istep], ω0[istep], udot_init, θdot_init, CtCabPdot_init,
            CtCabHdot_init, dt)

        j! = (J, x) -> system_jacobian!(J, x, assembly, prescribed_conditions,
            distributed_loads, istep, irow_pt, irow_beam, irow_beam1, irow_beam2,
            icol_pt, icol_beam, x0, v0[istep], ω0[istep], udot_init, θdot_init, CtCabPdot_init,
            CtCabHdot_init, dt)

        # solve system of equations
        if linear
            # linear analysis
            x .= 0.0
            f!(F, x)
            j!(J, x)
            x .-= J\F
        else
            df = OnceDifferentiable(f!, j!, x, F, J)

            result = NLsolve.nlsolve(df, x,
                linsolve=(x, A, b) -> ldiv!(x, lu(A), b),
                method=method,
                linesearch=linesearch,
                ftol=ftol,
                iterations=iterations)

            x .= result.zero
        end

        # add state to history
        if istep in save
            history[isave] = AssemblyState(system, assembly, prescribed_conditions = prescribed_conditions)
            isave += 1
        end

        # stop early if unconverged
        if !result.f_converged
            converged = false
            break
        end

        # exit loop if done iterating
        if istep == nstep
            break
        end

        # update time derivative terms for the next time step
        for ibeam = 1:nbeam
            icol = icol_beam[ibeam]
            # calculate udot for next time step
            u = SVector(x[icol], x[icol+1], x[icol+2])
            udot = 2/dt*u - udot_init[ibeam]
            udot_init[ibeam] = 2/dt*u + udot
            # calculate θdot for next time step
            θ = SVector(x[icol+3], x[icol+4], x[icol+5])
            θdot = 2/dt*θ - θdot_init[ibeam]
            θdot_init[ibeam] = 2/dt*θ + θdot
            # extract rotation parameters
            C = get_C(θ)
            Cab = assembly.elements[ibeam].Cab
            CtCab = C'*Cab
            # calculate CtCabPdot for next time step
            CtCabP = CtCab*SVector(x[icol+12], x[icol+13], x[icol+14])
            CtCabPdot = 2/dt*CtCabP - CtCabPdot_init[ibeam]
            CtCabPdot_init[ibeam] = 2/dt*CtCabP + CtCabPdot
            # calculate CtCabHdot for next time step
            CtCabH = CtCab*SVector(x[icol+15], x[icol+16], x[icol+17])
            CtCabHdot = 2/dt*CtCabH - CtCabHdot_init[ibeam]
            CtCabHdot_init[ibeam] = 2/dt*CtCabH + CtCabHdot
        end
    end

    # --- End Time Domain Simulation --- #

    return system, history, converged
end

# this function is currently unused
function linsolve!(x, A, b)

    # get scaling vectors
    r, s = linf_norm_scaling(A)

    # scale matrix and vector
    A .= A .* r .* s'
    b .= b .* r

    # do linear solve
    ldiv!(x, lu(A), b)

    # remove scaling from result
    x .= x .* s

    # remove scaling from matrix and vector
    A .= A ./ r ./ s'
    b .= b ./ r

    return x
end


text:
include("../common.jl")

compiler_target = nothing
try
    global compiler_target = triplet(platform_key(ARGS[end]))
    if typeof(compiler_target) <:UnknownPlatform
        error()
    end
catch
    error("This is not a typical build_tarballs.jl!  Must provide exactly one platform as the last argument!")
end

name = "BaseCompilerShard-$(compiler_target)"
version = v"2019.04.15"

# Refresh cache of cmake toolchains
cd(joinpath(@__DIR__, "cmake_toolchains")) do
    run(`bash build_toolchains.sh`)
end

# We're going to install things into this prefix, then package it up as a .tar.gz
temp_prefix() do prefix
    # Install KernelHeaders
    headers_path = find_tarball("KernelHeaders", "KernelHeaders.*$(compiler_target)").url
    @info("Unpacking $(headers_path)")
    unpack(headers_path, prefix.path)

    # Install Binutils
    #binutils_path = find_tarball("Binutils", "Binutils-$(compiler_target)").url
    #@info("Unpacking $(binutils_path)")
    #unpack(binutils_path, prefix.path)

    # Install Libc
    libc_path = ""
    if occursin("-gnu", compiler_target)
        # If we're a `-gnu` target, install Glibc.  We install version 2.12.2 on
        # x86_64 and i686, 2.19 on armv7l and aarch64, and 2.25 on ppc64le. These
        # versions are chosen through the time-tested method of "pick the earliest
        # one that works", so as to maximize compatibility.  We also do our best
        # to match versions of official Julia binaries.
        glibc_version_dict = Dict(
            :x86_64 => "2.12.2",
            :i686 => "2.12.2",
            :aarch64 => "2.19",
            :armv7l => "2.19",
            :powerpc64le => "2.25",
        )
        glibc_version = glibc_version_dict[arch(platform_key(compiler_target))]

        libc_path = find_tarball("Glibc", "Glibc.*$(glibc_version).*$(compiler_target)").url
    elseif occursin("-musl", compiler_target)
        libc_path = find_tarball("Musl", "Musl.*$(compiler_target)").url
    elseif occursin("-mingw", compiler_target)
        libc_path = find_tarball("Mingw", "Mingw.*$(compiler_target)").url
    elseif occursin("-freebsd", compiler_target)
        libc_path = find_tarball("FreeBSDLibc", "FreeBSDLibc.*").url
    elseif occursin("-darwin", compiler_target)
        libc_path = find_tarball("MacOSLibc", "MacOSLibc.*").url
    else
        error("Don't know how to install a Libc for $(compiler_target)!")
    end
    @info("Unpacking $(libc_path)")
    unpack(libc_path, prefix.path)

    # Deploy our cmake toolchain files
    toolchain_dir = joinpath(@__DIR__, "cmake_toolchains", "$(compiler_target)")
    for f in readdir(toolchain_dir)
        cp(joinpath(toolchain_dir, f), joinpath(prefix.path, f))
    end

    # We create a link from /opt/${target}/${target}/sys-root/usr/local/lib to /workspace/destdir/lib
    # This is the most reliable way for our sysroot'ed compilers to find destination libraries so far.
    usr_local_path = joinpath(prefix.path, "$(compiler_target)", "sys-root", "usr", "local")
    mkpath(usr_local_path)
    symlink("/workspace/destdir/lib", joinpath(usr_local_path, "lib"))
    symlink("/workspace/destdir/lib64", joinpath(usr_local_path, "lib64"))

    # Do a little bit of cleanup; at this point we don't care about logs, manifests, etc...
    rm(joinpath(prefix, "logs"); recursive=true, force=true)
    rm(joinpath(prefix, "manifests"); recursive=true, force=true)

    # Package this prefix up nice and tight
    mkpath(joinpath(pwd(), "products"))
    out_path = joinpath(pwd(), "products", name)
    package(prefix, out_path, version; platform=Linux(:x86_64, :glibc), verbose=true, force=true)
end

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
N: universe dimension
D: manifold dimension
T: coordinate type
"""
struct RectangleChart{N,D,T}
  origin::SVector{N,T}
  tangents::SVector{D,SVector{N,T}}
  normal::SVector{N,T}
end

function RectangleChart(o, tgs)
  n = tgs[1] × tgs[2]
  n = normalize(n)
  RectangleChart(o, tgs, n)
end

struct RectangleNBD{N,D,T}
  chart::RectangleChart{N,D,T}
  parm::SVector{D,T}
  cart::SVector{N,T}
end

cartesian(nbd::RectangleNBD) = nbd.cart

function neighborhood(rect::RectangleChart{3,2}, uv)
  o = rect.origin
  t1, t2 = rect.tangents
  u1, u2 = uv
  cart = o + u1*t1 + u2*t2
  RectangleNBD(rect, SVector(uv), cart)
end

getindex(nbd::RectangleNBD, i) = nbd.cart[i]
normal(nbd::RectangleNBD) = nbd.chart.normal

function quadpoints(chart::RectangleChart{3,2}, rule)

  T = eltype(chart.origin)

  o = chart.origin
  t1 = chart.tangents[1]
  t2 = chart.tangents[2]
  j = norm(t1 × t2)
  Q, P = rule
  U1, W1 = legendre(Q, zero(T), one(T))
  U2, W2 = legendre(P, zero(T), one(T))

  #pts = [o + u1 * t1 + u2 * t2 for u1 in U1 for u2 in U2]
  pts = [neighborhood(chart, (u1,u2)) for u1 in U1 for u2 in U2]
  wts = [j*w1*w2 for w1 in W1 for w2 in W2]

  collect(zip(pts,wts))
end


text:
using Clustering

function clustering(train_x::XData, mix_num::Int64; maxiter=200)::XBatches{<:Bool}
    if mix_num == 1
        return convert(XBatches, train_x)
    end

    n = num_examples(train_x)
    X = feature_matrix(train_x)'

    println("Running K-means clustering algorithm with num of components $mix_num, maximum iterations $maxiter")
    R = kmeans(X, mix_num; maxiter=maxiter)
    @assert nclusters(R) == mix_num
    a = assignments(R)

    clustered_train_x = Vector{PlainXData{Bool,BitMatrix}}()
    for k in 1 : mix_num
        push!(clustered_train_x, XData(convert(BitMatrix, X[:, findall(x -> x == k, a)]')))
    end
    return clustered_train_x
end


text:
module TestMixtape

using Test
using Mixtape
using CodeInfoTools
import Mixtape: CompilationContext, transform, allow, debug
using MacroTools

include("dynamic_overlay.jl")
include("rand_swap.jl")
include("insert_state.jl")
include("recursion.jl")
include("invalidation.jl")

EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")

for ex in readdir(EXAMPLES_DIR)
    @testset "$ex" begin
        # Just check the examples run, for now.
        include(joinpath(EXAMPLES_DIR, ex))
    end
end

end # module


text:
"""
    _get_cubeindex(iso_vals, iso)

given `iso_vals` and iso, return an 8 bit value corresponding
to each corner of a cube. In each bit position,
0 indicates in the isosurface and 1 indicates outside the surface,
where the sign convention indicates negative inside the surface
"""
@inline function _get_cubeindex(iso_vals, iso)
    cubeindex = iso_vals[1] < iso ? 0x01 : 0x00
    iso_vals[2] < iso && (cubeindex |= 0x02)
    iso_vals[3] < iso && (cubeindex |= 0x04)
    iso_vals[4] < iso && (cubeindex |= 0x08)
    iso_vals[5] < iso && (cubeindex |= 0x10)
    iso_vals[6] < iso && (cubeindex |= 0x20)
    iso_vals[7] < iso && (cubeindex |= 0x40)
    iso_vals[8] < iso && (cubeindex |= 0x80)
    cubeindex
end

"""
    _get_cubeindex_pos(iso_vals, iso)

given `iso_vals` and iso, return an 8 bit value corresponding
to each corner of a cube. In each bit position,
0 indicates in the isosurface and 1 indicates outside the surface,
where the sign convention indicates positive inside the surface
"""
@inline function _get_cubeindex_pos(iso_vals, iso)
    cubeindex = iso_vals[1] > iso ? 0x01 : 0x00
    iso_vals[2] > iso && (cubeindex |= 0x02)
    iso_vals[3] > iso && (cubeindex |= 0x04)
    iso_vals[4] > iso && (cubeindex |= 0x08)
    iso_vals[5] > iso && (cubeindex |= 0x10)
    iso_vals[6] > iso && (cubeindex |= 0x20)
    iso_vals[7] > iso && (cubeindex |= 0x40)
    iso_vals[8] > iso && (cubeindex |= 0x80)
    cubeindex
end

"""
    no_triangles(cubeindex)

Called after `_get_cubeindex`. Determines if a voxel index has triangles.
"""
@inline function _no_triangles(cubeindex::UInt8)
    cubeindex == 0x00 || cubeindex == 0xff
end

#
# General isosurface docstring
#

@doc """

    function isosurface(sdf::AbstractArray{T, 3}, [ method::AbstractMeshingAlgorithm ],
                         [ VertType = SVector{3,Float64} ], [ FaceType} = SVector{3, Int} ] ;
                         origin = SVector(-1.0,-1.0,-1.0), widths = SVector(2.0,2.0,2.0))

    function isosurface(f::Function, [ method::AbstractMeshingAlgorithm ],
                         [ VertType = SVector{3,Float64} ], [FaceType = SVector{3, Int} ] ;
                         origin = SVector(-1.0,-1.0,-1.0), widths = SVector(2.0,2.0,2.0)
                         samples=(24,24,24))`

`isosurface` is the general interface to all isosurface extraction algorithms.

Returns: (Vector{VertType}, Vector{FaceType})

Defaults:
- VertType = SVector{3,Float64} (positional)
- FaceType = SVector{3, Int} ] ; (positional)
- origin = SVector(-1.0,-1.0,-1.0) (keyword)
- widths = SVector(2.0,2.0,2.0) (keyword)
- samples=(24,24,24) (keyword, function sampling only)

`method` must be an instance of an `AbstractMeshingAlgorithm`, e.g.:
- MarchingCubes()
- MarchingTetrahedra()
- NaiveSurfaceNets()

If `isosurface` is called without a specified algorithm, it will default to MarchingCubes.

If a subtype of `AbstractArray` is specified, the mesh will by default be centered at the origin between
(-1,1) in each axis. This may be overridden by specifying a new origin and widths for the axis-aligned bounding box
using keywords of the same names. For example if we want our vertices in the range of (0,1), we can specify `origin=SVector(0,0,0)`
and `widths = SVector(1,1,1)`.

If a function is specified, it will be uniformly sampled in each axis by the amount specified in `samples`.
The function is called using the specifed `VertType`.

Performance Tips:
- ensure `VertType`, `origin`, and `widths` are all of the same type
- ensure the element type of `VertType` is the same as the specified isolevel

See also:
- MarchingCubes
- MarchingTetrahedra
- NaiveSurfaceNets

"""
function isosurface(A; kwargs...)
    isosurface(A, MarchingCubes(); kwargs...)
end


text:
module AST

export emit

"Correspondence between Julia types and Rust types from libc"
const RUST_TYPES = Dict{DataType, String}(
    Int8 => "int8_t",
    Int16 => "int16_t",
    Int32 => "int32_t",
    Int64 => "int64_t",
    UInt8 => "uint8_t",
    UInt16 => "uint16_t",
    UInt32 => "uint32_t",
    UInt64 => "uint16_t",
    Float32 => "c_float",
    Float64 => "c_double",
)

"Correspondence between Julia functions and Rust functions"
const RUST_FUNCTIONS = Dict{Symbol, String}(
    :+ => "add",
    :div => "div",
    :* => "mul",
    :% => "rem",
    :- => "sub",
    :~ => "not",
    :& => "bitand",
    :| => "bitor",
    :$ => "bitxor",
    :<< => "shl",
    :>> => "shr", # I believe this is the correct shift
)

function emit(li::LambdaInfo, sn::SlotNumber)::String
    emit(li, li.slotnames[sn.id])
end

"Generates a valid Rust identifier from a Julia Symbol"
function emit(li::LambdaInfo, s::Symbol)::String
    replace(string(s), r"[#0-9]", "")
end

function emit(li::LambdaInfo, t::DataType)::String
    RUST_TYPES[t]
end

"This will only handle global references to functions"
function emit(li::LambdaInfo, f::GlobalRef)::String
    RUST_FUNCTIONS[f.name]
end

function emit(li::LambdaInfo, ex::Expr)::String
    result = ""

    if ex.head === :call
        if length(ex.args) == 1
            result = "$(emit(li, ex.args[1]))()"
        else
            io = IOBuffer()
            print(io, "(", emit(li, ex.args[2]), ").", emit(li, ex.args[1]), "(")
            join(io, map(ex.args[3:end]) do arg
                emit(li, arg)
            end, ", ")
            print(io, ")")
            result = String(io)
        end

        if isleaftype(ex.typ)
            result_type = emit(li, ex.typ)
            result = "($result) as $result_type"
        end
    elseif ex.head === :return
        if isleaftype(ex.typ)
            result_type = emit(li, ex.typ)
            result = "return $(emit(li, ex.args[1])) as $result_type"
        else
            result = "return $(emit(li, ex.args[1]))"
        end
    elseif ex.head === :line
        result = ""
    else
        error("Unhandled expression type $(ex.head)")
    end

    if isleaftype(ex.typ)
        result_type = emit(li, ex.typ)
        result = "($result) as $result_type"
    end

    return result
end

function emit(li::LambdaInfo)::String
    ast = Base.uncompressed_ast(li)

    io = IOBuffer()
    print(io, "fn ", emit(li, li.def.name), "(")
    join(io, map(li.slotnames[2:end], li.slottypes[2:end]) do name, typ
        "$(emit(li, name)): $(emit(li, typ))"
    end, ", ")
    print(io, ") -> ", emit(li, li.rettype), " { ")
    for ex in ast
        println(io, "    ", emit(li, ex))
    end
    println(io, "}")

    return String(io)
end

function generate_rust(func::Function, args::ANY)::Tuple{String, String}
    li = code_typed(func, args; optimize=false)[1]

    if !li.inferred
        error("Unable to infer types. Try adding lots of type assertions!")
    end

    emit(li, li.def.name), emit(li)
end

end


text:
module sl_model

using DifferentialEquations, Flux, Optim, DiffEqFlux, Plots

model_params=[0.1, 6.28, 0.1, 1.0];
all_ic=[[0.0,0.01], [π,0.2], [π/2.0, 0.1]];
tend=15.0
tstart=0.0
sampling=0.1

function sl_oscillator(du,u,p,t)
    γ, ω0, β, K=p
    ψ,R=u
    du[1]=ω0+K/2.0*sin(β)*(1+R^2)
    du[2]=-γ*R+K/2.0*cos(β)*R*(1-R^2)
end

function InitPlot()

    mycolor=[:blue, :red, :green, :purple, :black, :cyan, :orange]
    odedata=[solve(ODEProblem(sl_oscillator,η,(tstart,tend), model_params), Tsit5()) for η in all_ic];
    pl = plot(odedata[1],proj=:polar,vars=(1,2), color=:blue, legend=false)

    for k in 2:length(all_ic)
        plot!(pl, odedata[k], color=mycolor[k], vars=(1,2), proj=:polar)
    end
    title!(pl, "Oscillator Model")
    xlabel!(pl, "X")
    ylabel!(pl, "Y")
    savefig("SL_init.png")
    display(pl)

end

function predict_adjoint(param) # Our 1-layer neural network
    prob=ODEProblem(sl_oscillator,[0.1,0.1],(tstart,tend), model_params)
    Array(concrete_solve(prob,Tsit5(),param[1:2],param[3:end],saveat=tstart:sampling:tend))
end

# Generate some data to fit, and add some noise to it
data=predict_adjoint([0.0,0.08,0.1, 6.28, 0.1, 1.0])
σN=0.05
data+=σN*randn(size(data))
data=abs.(data)

function loss_adjoint(param)
  prediction = predict_adjoint(param)
  loss = sum(abs2,prediction - data)
  loss
end

function train_model(;pguess=[0.1,0.1,0.15, 6.0, 0.2, 0.8])
    println("The initial loss function is $(loss_adjoint(pguess)[1])")
    resinit=DiffEqFlux.sciml_train(loss_adjoint,pguess,ADAM(),maxiters=3000)
    res = DiffEqFlux.sciml_train(loss_adjoint,resinit.minimizer,BFGS(initial_stepnorm = 0.0001))
    println("The parameters are $(res.minimizer) with final loss value $(res.minimum)")
    return(res)
end

function plotFit(param)

    validationPlot(param, param[1:2])

end

function validationPlot(param, ic)

    lay = @layout [ a{0.6w} [b; c]]

    sol_fit=solve(ODEProblem(sl_oscillator,ic,(0.0,15.0),param[3:end]), Tsit5(), saveat=tstart:0.01:tend)
    sol_data=solve(ODEProblem(sl_oscillator,ic,(0.0,15.0),model_params), Tsit5(), saveat=tstart:0.01:tend)
    pl1=plot(sol_fit,lw=2, legend=false, vars=(1,2), color=:green, proj = :polar)
    plot!(pl1, sol_data, color=:blue, proj=:polar, vars=(1,2))

    #radius plot
    pl2=plot(sol_fit,lw=2, legend=false, vars=(0,2), color=:green)
    plot!(pl2, sol_data, color=:blue, vars=(0,2), lw=2)

    #angle plot
    pl3=plot(sol_fit.t,sin.(sol_fit[1,:]), lw=2, legend=false, color=:green)
    plot!(pl3, sol_data.t, sin.(sol_data[1,:]), color=:blue)

    pl=plot(pl1, pl2, pl3, layout=lay)
    savefig("Sl_valid_2.png")
    display(pl)

end

end #module

sl_model.InitPlot()
res=sl_model.train_model()
sl_model.plotFit(res.minimizer)
sl_model.validationPlot(res.minimizer, [π, 0.6])

max_stars_repo_name: khannay/FittingParamsDiffEqFlux

max_issues_repo_name: khannay/FittingParamsDiffEqFlux

max_forks_repo_name: khannay/FittingParamsDiffEqFlux


text:
# This file was generated, do not modify it. # hide
ypred = predict_mode(mach, X)
misclassification_rate(ypred, y)


text:
module PDiscRefFEs

using Gridap

export PDiscRefFE

import Gridap: shfbasis
import Gridap: polytope
import Gridap: nfacedofs
import Gridap: dofbasis

struct PDiscRefFE{D,T} <: RefFE{D,T}
  reffe_tet::LagrangianRefFE{D,T}
  polytope::Polytope{D}
  nfacedofs::Vector{Vector{Int}}
end

function PDiscRefFE(::Type{T},D::Integer,order::Integer) where T

  extrusion = fill(HEX_AXIS,D)
  polytope = Polytope(extrusion...)

  PDiscRefFE(T,polytope,order)

end

function PDiscRefFE(
  ::Type{T}, polytope::Polytope{D}, order::Integer) where {T,D}

  extrusion_tet = fill(TET_AXIS,D)
  polytope_tet = Polytope(extrusion_tet...)

  reffe_tet = LagrangianRefFE(T,polytope_tet,order)

  # We assign all DOFs to the interior of the cell.
  # It is not possible to construct continuous spaces
  # with this RefFE
  nnfaces = length(polytope.nfaces)
  ndofs = length(reffe_tet.shfbasis)
  nfacedofs = [Int[] for i in 1:nnfaces]
  nfacedofs[end] = collect(1:ndofs)

  PDiscRefFE{D,T}(reffe_tet,polytope,nfacedofs)

end

dofbasis(this::PDiscRefFE) = this.reffe_tet.dofbasis

polytope(this::PDiscRefFE) = this.polytope

shfbasis(this::PDiscRefFE) = this.reffe_tet.shfbasis

nfacedofs(this::PDiscRefFE) = this.nfacedofs

end # module


text:
@testset "50.powx-n.jl" begin
    @test my_pow(2.0, 0) ≈ 1.0
    @test all(my_pow(2.0, i) ≈ 2.0^i for i in 1:20)
    @test all(my_pow(-2.0, i) ≈ (-2.0)^i for i in 1:20)
end


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

include("../common.jl")
name = "oneAPI_Level_Zero_Loader"

# Bash recipe for building across all platforms
script = raw"""
cd level-zero
install_license LICENSE

mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DCMAKE_BUILD_TYPE=Release ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc")
]
platforms = expand_cxxstring_abis(platforms)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("OpenCL_Headers_jll"),

    # The oneAPI Level Zero loader expects an implementation of its spec to be available.
    # That means we need to dlopen the implementation, e.g. Intel's NEO, before the loader!
    # To make sure that implementation sits at the same API level of the loader, we have
    # both packages depend on the oneAPI Level Zero Headers package.
    #
    # Users of these packages, e.g. oneAPI.jl, should only depend on the loader and an
    # implementation, but will need to make sure to load one before the other.
    Dependency("oneAPI_Level_Zero_Headers_jll"; compat="=$api_version"),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"5")

# bump counter: 3


text:
module MaxHelpingHandAmbienceVolumeTrigger

using ..Ahorn, Maple

@mapdef Trigger "MaxHelpingHand/AmbienceVolumeTrigger" AmbienceVolumeTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    from::Number=0.0, to::Number=0.0, direction::String="NoEffect")

const placements = Ahorn.PlacementDict(
    "Ambience Volume (max480's Helping Hand)" => Ahorn.EntityPlacement(
        AmbienceVolumeTrigger,
        "rectangle"
    )
)

function Ahorn.editingOptions(trigger::AmbienceVolumeTrigger)
    return Dict{String, Any}(
        "direction" => Maple.trigger_position_modes
    )
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
# This is the R2-IBEA weight vector generation algorithm as described in:
# Dung H. Phan and Junichi Suzuki, "R2-IBEA: R2 Indicator Based Evolutionary
# Algorithm for Multiobjective Optimization", 2013.

function generate_random_weight_vector(m)

  x = sort(rand(m-1))

  w = zeros(m)
  w[1] = x[1]
  for i in 2:(m-1)
    w[i] = x[i] - x[i-1]
  end
  w[m] = 1 - x[m-1]

  w

end

function generate_weight_vectors(numobjectives, numvectors, maxnum_iterations = 10000)

  # The final returned set of vectors will be returned as a
  # numobjectives * numvectors matrix but we need room for one more column
  # to keep the last candidate vector to be added.
  generated_vectors = zeros(numobjectives, numvectors+1) # W
  numgenerated = 0

  for t in 1:maxnum_iterations

    numgenerated += 1
    generated_vectors[:, numgenerated] = generate_random_weight_vector(numobjectives)

    if numgenerated > numvectors

      i = index_to_column_with_min_hypervolume_contribution(generated_vectors)
      generated_vectors[:, i] = generated_vectors[:, numgenerated]
      numgenerated -= 1

    end
  end

  # Delete the additional column and return the correct number of vectors
  generated_vectors[:, 1:numobjectives]

end

function index_to_column_with_min_hypervolume_contribution(vectors)

end

#
# Calculate the approximate hypervolume contribution of each of the L
# M-dimensional fitness vectors given as a M*L matrix, v. The calculation is
# via monte-carlo simulation.
#
function approximate_hypervolume_monte_carlo(v::Array{Float64, 2};
  upperbound = maximum(v, 2), numsamples = 100000)

  M, L = size(v)

  # Get the lower bound point and the delta, i.e. size of box between lower and upper bounds.
  lowerbound = minimum(v, 2)
  delta = upperbound - lowerbound

  # Alloc mem once so we can write into it later for speed.
  sampledpoint = zeros(M)

  # Count the number of sampled points that are dominated by the given points.
  num_dominated = 0
  for i in 1:numsamples

    # Sample in the box from lowerbound to upperbound.
    for m in 1:M
      sampledpoint[m] = lowerbound[m] + rand() * delta[m]
    end

    # Compare sampled point to points on the front and count if it is dominated.
    for frontpoint in 1:L
      if dominates(v[:, frontpoint], sampledpoint)
        num_dominated += 1
        break
      end
    end

  end

  # Approximate hypervolume based on how many of the sampled points were dominated.
  return prod(delta) * num_dominated / numsamples

end

#
# Return true iff points a dominates point b in the Pareto sense.
#
function dominates(a, b)

  dominates_in_at_least_one = false

  for i in 1:length(a)
    if a[i] > b[i]
      return false
    elseif a[i] < b[i]
      dominates_in_at_least_one = true
    end
  end

  return dominates_in_at_least_one

end


text:
# Miscellaneous utilities

"""
```
crossover(x::Array{T}, y::Array{T}) where {T<:Real}
```

Find where `x` crosses over `y` (returns boolean vector where crossover occurs)
"""
function crossover(x::Array{T}, y::Array{T}) where {T<:Real}
    @assert size(x,1) == size(y,1)
    out = falses(size(x))
    @inbounds for i in 2:size(x,1)
        out[i] = ((x[i] > y[i]) && (x[i-1] < y[i-1]))
    end
    return out
end

"""
```
crossunder(x::Array{T}, y::Array{T}) where {T<:Real}
```

Find where `x` crosses under `y` (returns boolean vector where crossunder occurs)
"""
function crossunder(x::Array{T}, y::Array{T}) where {T<:Real}
    @assert size(x,1) == size(y,1)
    out = falses(size(x))
    @inbounds for i in 2:size(x,1)
        out[i] = ((x[i] < y[i]) && (x[i-1] > y[i-1]))
    end
    return out
end

"""
```
wilder_sum(x::Vector{T}; n::Int=10)::Vector{T}
```

Welles Wilder summation of an array
"""
function wilder_sum(x::Vector{T}; n::Int=10)::Vector{T} where {T<:Real}
    @assert n<size(x,1) && n>0 "Argument n is out of bounds."
    nf = float(n)  # type stability -- all arithmetic done on floats
    out = zeros(size(x))
    out[1] = x[1]
    @inbounds for i = 2:size(x,1)
        out[i] = x[i] + out[i-1]*(nf-1.0)/nf
    end
    return out
end
wilder_sum(X::Matrix; n::Int=10)::Matrix = hcat((wilder_sum(X[:,j], n=n) for j in 1:size(X,2))...)

"""
(Adapted from StatsBase: https://raw.githubusercontent.com/JuliaStats/StatsBase.jl/master/src/scalarstats.jl)

Compute the mode of an arbitrary array::Array{T}
"""
function mode(a::AbstractArray{T}) where {T<:Real}
    isempty(a) && error("mode: input array cannot be empty.")
    cnts = Dict{T,Int}()
    # first element
    mc = 1
    mv = a[1]
    cnts[mv] = 1
    # find the mode along with table construction
    @inbounds for i = 2 : length(a)
        x = a[i]
        if haskey(cnts, x)
            c = (cnts[x] += 1)
            if c > mc
                mc = c
                mv = x
            end
        else
            cnts[x] = 1
            # in this case: c = 1, and thus c > mc won't happen
        end
    end
    return mv
end

"""
```
diffn(x::Vector{T}; n::Int=1)::Vector{T} where {T<:Real}
diffn(X::Matrix; n::Int=1)::Matrix = hcat([diffn(X[:,j], n=n) for j in 1:size(X,2)]...)
```

Lagged differencing
"""
function diffn(x::Vector{T}; n::Int=1)::Vector{T} where {T<:Real}
    @assert n<size(x,1) && n>0 "Argument n out of bounds."
    dx = zeros(size(x))
    dx[1:n] .= NaN
    @inbounds for i=n+1:size(x,1)
        dx[i] = x[i] - x[i-n]
    end
    return dx
end
diffn(X::Matrix; n::Int=1)::Matrix = hcat([diffn(X[:,j], n=n) for j in 1:size(X,2)]...)


text:
using SurrealNumbers
using Test

lr_eq(x, l, r) = issetequal(left(x), l) && issetequal(right(x), r)
lr_eq(x1, x2) = lr_eq(x1, left(x2), right(x2))

@testset "util" begin include("util.jl") end
@testset "sets" begin include("sets.jl") end
@testset "SurrealForm" begin include("SurrealForm.jl") end
@testset "builtins" begin include("builtins.jl") end


text:
module RadioButtons

Base.Experimental.@optlevel 1

using Redux
using CImGui

# actions
abstract type AbstractRadioButtonAction <: AbstractSyncAction end

get_label(a::AbstractRadioButtonAction) = a.label

"""
    Rename(label, new_label)
Change radio button's label to `new_label`.
Note that renaming may also change the identifier, please refer to `help?> CImGui.PushID`.
"""
struct Rename <: AbstractRadioButtonAction
    label::String
    new_label::String
end

"""
    SetTriggeredTo(label, is_triggered)
Update widget's `is_triggered` state.
"""
struct SetTriggeredTo <: AbstractRadioButtonAction
    label::String
    is_triggered::Bool
end

"""
    SetActiveTo(label, is_active)
Set radio button's state to `is_active`.
"""
struct SetActiveTo <: AbstractRadioButtonAction
    label::String
    is_active::Bool
end

"""
    AddButton(label::AbstractString, is_active = false)
Add a [`RadioButton`](@ref).
"""
struct AddButton <: AbstractRadioButtonAction
    label::String
    is_active::Bool
end
AddButton(label::AbstractString) = AddButton(label, false)

"""
    DeleteButton(label::AbstractString)
Delete the [`RadioButton`](@ref).
"""
struct DeleteButton <: AbstractRadioButtonAction
    label::String
end

"""
    OnlySetThisOneToActive(label::AbstractString)
Clear all the `is_active` states in the button list and only set the `label`'s `is_active` state to active.
"""
struct OnlySetThisOneToActive <: AbstractRadioButtonAction
    label::String
end

# state
abstract type AbstractRadioButtonState <: AbstractImmutableState end

get_label(s::AbstractRadioButtonState) = s.label

"""
    RadioButtons.State(label::AbstractString, is_active = false, is_triggered = false)
A radio button state which contains a label a.k.a the identifier, a flag value `is_active` and
a flag value `is_triggered` that records the state of the latest poll events.
"""
struct State <: AbstractRadioButtonState
    label::String
    is_active::Bool
    is_triggered::Bool
end
State(label::AbstractString, is_active = false) = State(label, is_active, false)

is_active(s::State) = s.is_active
is_triggered(s::State) = s.is_triggered

# reducers
reducer(state::AbstractState, action::AbstractAction) = state
reducer(state::Vector{<:AbstractState}, action::AbstractAction) = state
reducer(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

reducer(s::State, a::Rename) = State(a.new_label, s.is_active, s.is_triggered)
reducer(s::State, a::SetActiveTo) = State(s.label, a.is_active, s.is_triggered)
reducer(s::State, a::SetTriggeredTo) = State(s.label, s.is_active, a.is_triggered)

reducer(s::Dict{String,<:AbstractRadioButtonState}, a::AbstractRadioButtonAction) =
    Dict(k => (get_label(v) == get_label(a) ? reducer(v, a) : v) for (k, v) in s)

reducer(s::Vector{<:AbstractRadioButtonState}, a::AbstractRadioButtonAction) = map(s) do s
    get_label(s) == get_label(a) ? reducer(s, a) : s
end

reducer(s::Vector{<:AbstractRadioButtonState}, a::OnlySetThisOneToActive) = map(s) do s
    get_label(s) == get_label(a) ? State(s.label, true, s.is_triggered) : State(s.label, false, s.is_triggered)
end

reducer(s::Vector{<:AbstractRadioButtonState}, a::AddButton) =
    [s..., State(a.label, a.is_active, false)]
reducer(s::Vector{<:AbstractRadioButtonState}, a::DeleteButton) =
    filter(s -> get_label(s) !== get_label(a), s)

# helper
"""
    RadioButton(store::AbstractStore, get_state=Redux.get_state) -> Bool
Return `true` when triggered.
`get_state` is a router function that tells how to find the target state from `store`. `chain_action` is for chaining upstream actions.
"""
function RadioButton(store::AbstractStore, get_state=Redux.get_state)
    s = get_state(store)
    is_triggered = CImGui.RadioButton(get_label(s), is_active(s))
    dispatch!(store, SetTriggeredTo(get_label(s), is_triggered))
    is_triggered && dispatch!(store, SetActiveTo(get_label(s), true))
    return is_triggered
end

end # module


text:
using Electron
using URIParser
using FilePaths
using Test

Electron.prep_test_env()

@testset "Electron" begin

@testset "Core" begin

w = Window(URI("file://test.html"))

a = applications()[1]

@test isa(w, Window)

@test length(applications()) == 1
@test length(windows(a)) == 1

res = run(w, "Math.log(Math.exp(1))")

@test res == 1

@test_throws ErrorException run(w, "syntaxerror")

res = run(a, "Math.log(Math.exp(1))")

@test res ==1

close(w)
@test length(applications()) == 1
@test isempty(windows(a))

w2 = Window(join(@__PATH__, p"test.html"))

toggle_devtools(w2)

close(a)
@test length(applications()) == 1
@test length(windows(a)) == 0

sleep(1)
@test isempty(applications())
@test isempty(windows(a))

w3 = Window(Dict("url" => string(URI("file://test.html"))))

w4 = Window(URI("file://test.html"), options=Dict("title" => "Window title"))

w5 = Window("", options=Dict("title" => "Window title"))

a2 = applications()[1]

w6 = Window(a2, "", options=Dict("title" => "Window title"))

w7 = Window(a2)

run(w7, "sendMessageToJulia('foo')")

@test take!(msgchannel(w7)) == "foo"

load(w7, "bar")

run(w7, "sendMessageToJulia(window.document.documentElement.innerHTML)")

@test occursin("bar", take!(msgchannel(w7)))

load(w7, join(@__PATH__, p"test.html"))
load(w7, URI(join(@__PATH__, p"test.html")))

close(w7)

close(w3)
close(w4)
close(w5)
close(w6)
close(a2)

end # testset "Core"

@testset "ElectronAPI" begin
    win = Window()

    @test (ElectronAPI.setBackgroundColor(win, "#000"); true)
    @test ElectronAPI.isFocused(win) isa Bool

    bounds = ElectronAPI.getBounds(win)
    boundskeys = ["width", "height", "x", "y"]
    @test Set(boundskeys) <= Set(keys(bounds))
    @test all(isa.(get.(Ref(bounds), boundskeys, nothing), Real))

    close(win)
end

for app in applications()
    close(app)
end

end