text:
using Plots
include(joinpath(@__DIR__, "dpo.jl"))
include(joinpath(pwd(), "examples/direct_policy_optimization/simulate.jl"))

# Unpack trajectories
x̄, ū = unpack(z̄, prob)
x, u = unpack(z, prob_dpo.prob.prob.nom)
Θ = get_policy(z, prob_dpo)

# Simulation setup
model_sim = model
T_sim = 10 * T

W = Distributions.MvNormal(zeros(model_sim.n),
	Diagonal(5.0e-4 * ones(model_sim.n)))
w = rand(W, T_sim)

W0 = Distributions.MvNormal(zeros(model_sim.n),
	Diagonal(5.0e-4 * ones(model_sim.n)))
w0 = rand(W0, 1)

tf_nom = sum([ū[t][end] for t = 1:T-1])
t_nom = range(0, stop = tf_nom, length = T)
t_sim_nom = range(0, stop = tf_nom, length = T_sim)

tf_dpo = sum([u[t][end] for t = 1:T-1])
t_nom_dpo = range(0, stop = tf_dpo, length = T)
t_sim_nom_dpo = range(0, stop = tf_dpo, length = T_sim)

dt_sim_nom = tf_nom / (T_sim - 1)
dt_sim_dpo = tf_dpo / (T_sim - 1)

z_dpo, u_dpo, J_dpo, Jx_dpo, Ju_dpo = _simulate(
	model_sim,
	policy, Θ,
	x, u,
	Q, R,
	T_sim, u[1][end],
	vec(x[1] + w0), w,
	_norm = 2,
	ul = ul[1], uu = uu[1],
	u_idx = (1:model.m - 1))

jt = []
jtx = []
jtu = []

jd = []
jdx = []
jdu = []

n_trials = 100
for i = 1:n_trials
	W = Distributions.MvNormal(zeros(model_sim.n),
		Diagonal(δ * ones(model_sim.n)))
	w = rand(W, T_sim)

	W0 = Distributions.MvNormal(zeros(model_sim.n),
		Diagonal(δ0 * ones(model_sim.n)))
	w0 = rand(W0, 1)

	# Simulate
	z_tvlqr, u_tvlqr, J_tvlqr, Jx_tvlqr, Ju_tvlqr = _simulate(
		model_sim,
		policy, K,
	    x̄, ū,
		Q, R,
		T_sim, ū[1][end],
		vec(x̄[1] + w0), w,
		_norm = 2,
		ul = ul[1], uu = uu[1],
		u_idx = (1:model.m - 1))

	z_dpo, u_dpo, J_dpo, Jx_dpo, Ju_dpo = _simulate(
		model_sim,
		policy, Θ,
	    x, u,
		Q, R,
		T_sim, u[1][end],
		vec(x[1] + w0), w,
		_norm = 2,
		ul = ul[1], uu = uu[1],
		u_idx = (1:model.m - 1))

	push!(jt, J_tvlqr)
	push!(jtx, Jx_tvlqr)
	push!(jtu, Ju_tvlqr)

	push!(jd, J_dpo)
	push!(jdx, Jx_dpo)
	push!(jdu, Ju_dpo)
end

# average tracking performance

mean(jt)
mean(jd)

std(jt)
std(jd)

mean(jtx)
mean(jdx)

std(jtx)
std(jdx)

mean(jtu)
mean(jdu)

std(jtu)
std(jdu)

using RDataSets, StatsPlots
boxplot(["lqr" "dpo"], hcat(jt, jd), title="Tracking error (n = $n_trials)", label=["lqr" "dpo"], color = [:purple :orange])
boxplot(["lqr" "dpo"], hcat(jtx, jdx), title="State tracking error (n = $n_trials)", label=["lqr" "dpo"], color = [:purple :orange])
boxplot(["lqr" "dpo"], hcat(jtu, jdu), title="Control tracking error (n = $n_trials)", label=["lqr" "dpo"], color = [:purple :orange])

# plot(t_nom[1:end-1], hcat(ū...)[1:4,:]', label = "",
# 	color = :purple, width = 2.0)
# plot!(t_sim_nom[1:end-1], hcat(u_tvlqr...)', label = "", color = :black)
#
# plot(t_nom_dpo[1:end-1], hcat(u...)[1:4,:]', label = "",
# 	color = :orange, width = 2.0)
# plot!(t_sim_nom_dpo[1:end-1], hcat(u_dpo...)[1:4,:]', label = "", color = :black)

# Visualize
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

urdf = joinpath(pwd(), "models/biped/urdf/biped_left_pinned.urdf")
mechanism = parse_urdf(urdf, floating=false)
mvis = MechanismVisualizer(mechanism,
    URDFVisuals(urdf, package_path=[dirname(dirname(urdf))]), vis)

visualize!(mvis, model, x, Δt = u[1][end])
visualize!(mvis, model, z_dpo, Δt = dt_sim_dpo)

# visualize!(mvis, model, z_tvlqr, Δt = dt_sim_nom)
# visualize!(mvis, model, z_dpo, Δt = dt_sim_dpo)


text:
const se = 4.6051701859880916e+02
const se05 = 2.3025850929940458e+02
const ke = 1.0e-200
const ke05 = 1.0e-100
const keinv = 1.0e200
const ke05inv = 1.0e100

## TODO: Create PSP Struct

PSP_VDSx(VDS) = VDS^2/(sqrt(VDS^2+0.01)+0.1)
PSP_VSBx(VSBs,VDS) = VSBs+0.5*(VDS-PSP_VDSx(VDS))

MINA(x,y,a) = (x+y-sqrt((x-y)^2.0+a))/2.0
χ(y::Number) = y^2/(2.0+y^2)
χp(y::Number) = 4*y/(2+y^2)^2
χpp(y::Number) = (8-12*y^2)/(2+y^2)^3.0
ν(a,c) = a+c
μ1(a,c,τ) = (a+c)^2/τ + c^2/2 - a
σ1(a,c,τ,η) = a*ν(a,c)/(μ1(a,c,τ)+(c^2 /3-a)+c*(a+c)/μ1(a,c,τ)) + η
μ2(a,b,c,τ) = (a+c)^2/τ + c^2/2 - a*b
σ2(a,b,c,τ,η) = a*ν(a,c)/(μ2(a,b,c,τ)+(c^2/3-a*b)+c*(a+c)/μ2(a,b,c,τ)) + η

P3(u) =  (1.0 + (u) * (1.0 + 0.5 * ((u) * (1.0 + (u)/3.0))))
ϕx(φb) = 2*0.95*φb
aϕ(φb) = 2.5e-3*(2*φb)^2
bϕ(φb) = 2.5e-3*(2*φb)^2
ϕxs(φb) = MINA(ϕx(φb)-0.5*sqrt(bϕ(φb)),0,aϕ(φb))
ϕV(φb) = MINA(Vs,Vs+Vd,bϕ(φb))+ϕx(φb)
Vsbs(Vs,φb) = Vs-MINA(ϕV(φb),0,aϕ(φb))+ϕxs(φb)

function expl_high(x)
    if ((x) < se05)
            res       = exp(x)
    else
            res       =  ke05inv * P3((x) - se05)
    end
end

function expl_low(x)
    if ((x) > -se05)
            res       =  exp(x)
    else
            res       = ke05 / P3(-se05 - (x))
    end
    return res
end

function sigma2(a,b,c,tau,eta)
    nu           =  (a) + (c)
    mutau        =  (nu) * (nu) + (tau) * (0.5 * ((c) * (c)) - (a) * (b))
    y            =  (eta) + (a) * nu * (tau) / (mutau + (nu / mutau) * (tau) * (tau) * (c) * ((c) * (c) * 1/3.0 - (a) * (b)))
    return y
end

function sigma(a,c,tau,eta)
    nu           =  (a) + (c)
    mutau        =  nu * nu + (tau) * (0.5 * ((c) * (c)) - (a))
    y            =  (eta) + (a) * nu * (tau) / (mutau + (nu / mutau) * (tau) * (tau) * (c) * ((c) * (c) * 1/3.0 - (a)))
    return y
end

function γ(N,tox)
    Cox = ϵ₀*3.9/tox
    return sqrt(2*q*11.7*ϵ₀*N)/Cox
end

function γC(N,Cox)
    return sqrt(2*q*11.7*ϵ₀*N)/Cox
end

function SP_psi(SD::SemiconductorData,Vg,Vfb,tox,V=0)
    phit = kb*Temperature(SD)
    xg = (Vg-(Vfb))/phit;
    xn = (2.0*abs(Ei(SD))+V)/phit;
    G = γ(abs(C(SD)),tox)/sqrt(phit);
    xsub = xg + G^2.0/2 -G*(xg-1+G^2.0/4)^(1/2)
    eta = s(xsub,xn+3,5)
    a = (xg-eta)^2-G.^2*eta+G^2
    c = 2.0*(xg-eta)+G^2
    t = xn-eta+log(a/G^2)
    x0 = eta+sigmaSP(a,c,t)
    D0 = exp(x0-xn)
    p = 2*(xg-x0)+G.^2*(1+D0)
    q = (xg-x0)^2-G^2*(x0+D0-1)
    x = x0+2*q./(p+sqrt(p^2-2*(2-G^2*D0)*q))
    return phit*x
end

function PSP_psiS(Vg,Na,tox,Vfb,T,φb,Vsb=0.0)
    phit = kb*T
    φb *= 2.0
    φb = abs(φb)
    phib = φb
    xg = (Vg-Vfb)/phit
    Gf = γ(Na,tox)/sqrt(phit)
    inv_phit1 = 1/phit
    oneSixth = 1/6.0
    invSqrt2 = 1/sqrt(2.0)
    xi           =  1.0 + Gf * invSqrt2
    inv_xi       =  1.0 /  xi
    Ux           =  Vsb * inv_phit1
    xn         =  phib * inv_phit1 + Ux
    if (xn < se)
        delta     =  exp(-xn)
    else
        delta     = ke / P3(xn - se)
    end
    margin       =  1.0e-5 * xi

    inv_xi = 1/xi
    Gf2 = Gf^2
    inv_Gf2 = 1/Gf2
    if (abs(xg) <= margin)
        SP_S_temp1 =  inv_xi * inv_xi * oneSixth * invSqrt2
        sp         =  xg * inv_xi * (1.0 + xg * (1.0 - (delta)) * Gf * SP_S_temp1)
    else
        if (xg < -margin)
            SP_S_yg     = -xg
            SP_S_ysub   = 1.25 * (SP_S_yg * inv_xi)
            SP_S_eta    = 0.5 * (SP_S_ysub + 10 - sqrt((SP_S_ysub - 6.0) * (SP_S_ysub - 6.0) + 64.0))
            SP_S_temp   = SP_S_yg - SP_S_eta
            SP_S_a      = SP_S_temp * SP_S_temp + Gf2*(SP_S_eta + 1.0)
            SP_S_c      = 2.0 * SP_S_temp - Gf2
            SP_S_tau    = -SP_S_eta + log(SP_S_a * inv_Gf2)
            SP_S_y0 = sigma(SP_S_a, SP_S_c, SP_S_tau, SP_S_eta)
            SP_S_delta0 = expl_high(SP_S_y0)
            SP_S_delta1 = 1.0 / SP_S_delta0
            SP_S_temp   = 1.0 / (2.0 + SP_S_y0 * SP_S_y0)
            SP_S_xi0    = SP_S_y0 * SP_S_y0 * SP_S_temp
            SP_S_xi1    = 4.0 * (SP_S_y0 * SP_S_temp * SP_S_temp)
            SP_S_xi2    = (8.0 * SP_S_temp - 12.0 * SP_S_xi0) * SP_S_temp * SP_S_temp
            SP_S_temp   = SP_S_yg - SP_S_y0
            SP_S_temp1  = (delta) * SP_S_delta1
            SP_S_pC     = 2.0 * SP_S_temp + Gf2 * (SP_S_delta0 - 1.0 - SP_S_temp1 + (delta) * (1.0 - SP_S_xi1))
            SP_S_qC     = SP_S_temp * SP_S_temp - Gf2 * (SP_S_delta0 - SP_S_y0 - 1.0 + SP_S_temp1 + (delta) * (SP_S_y0 - 1.0 - SP_S_xi0))
            SP_S_temp   = 2.0 - Gf2 * (SP_S_delta0 + SP_S_temp1 - (delta) * SP_S_xi2)
            SP_S_temp   = SP_S_pC * SP_S_pC - 2.0 * (SP_S_qC * SP_S_temp)
            sp          = -SP_S_y0 - 2.0 * (SP_S_qC / (SP_S_pC + sqrt(SP_S_temp)))
        else
            SP_xg1    = 1.0 / (1.25 + Gf * 7.324648775608221e-001)
            SP_S_A_fac= (xi * 1.25 * SP_xg1 - 1.0) * SP_xg1
            SP_S_xbar = xg * inv_xi * (1.0 + SP_S_A_fac * xg)
            SP_S_temp = expl_low(-SP_S_xbar)
            SP_S_w    = 1.0 - SP_S_temp
            SP_S_x1   = xg + Gf2 * 0.5 - Gf * sqrt(xg + Gf2 * 0.25 - SP_S_w)
            SP_S_bx   = (xn) + 3.0
            SP_S_eta  = MINA(SP_S_x1, SP_S_bx, 5.0) - 0.5 * (SP_S_bx - sqrt(SP_S_bx * SP_S_bx + 5.0))
            SP_S_temp = xg - SP_S_eta
            SP_S_temp1= exp(-SP_S_eta)
            SP_S_temp2= 1.0 / (2.0 + SP_S_eta * SP_S_eta)
            SP_S_xi0  = SP_S_eta * SP_S_eta * SP_S_temp2
            SP_S_xi1  = 4.0 * (SP_S_eta * SP_S_temp2 * SP_S_temp2)
            SP_S_xi2  = (8.0 * SP_S_temp2 - 12.0 * SP_S_xi0) * SP_S_temp2 * SP_S_temp2
            SP_S_a    = max(1.0e-40, SP_S_temp * SP_S_temp - Gf2 * (SP_S_temp1 + SP_S_eta - 1.0 - (delta) * (SP_S_eta + 1.0 + SP_S_xi0)))
            SP_S_b    = 1.0 - 0.5 * (Gf2 * (SP_S_temp1 - (delta) * SP_S_xi2))
            SP_S_c    = 2.0 * SP_S_temp + Gf2 * (1.0 - SP_S_temp1 - (delta) * (1.0 + SP_S_xi1))
            SP_S_tau  = (xn) - SP_S_eta + log(SP_S_a / Gf2)
            SP_S_x0 = sigma2(SP_S_a, SP_S_b, SP_S_c, SP_S_tau, SP_S_eta)
            if (SP_S_x0 < se05)
                SP_S_delta0 = exp(SP_S_x0)
                SP_S_delta1 = 1.0 / SP_S_delta0
                SP_S_delta0 = (delta) * SP_S_delta0
            else
                if (SP_S_x0 > (xn) - se05)
                    SP_S_delta0 = exp(SP_S_x0 - (xn))
                    SP_S_delta1 = (delta) / SP_S_delta0
                else
                    SP_S_delta0 = ke05 / P3((xn) - SP_S_x0 - se05)
                    SP_S_delta1 = ke05 / P3(SP_S_x0 - se05)
                end
            end
            SP_S_temp   = 1.0 / (2.0 + SP_S_x0 * SP_S_x0)
            SP_S_xi0    = SP_S_x0 * SP_S_x0 * SP_S_temp
            SP_S_xi1    = 4.0 * (SP_S_x0 * SP_S_temp * SP_S_temp)
            SP_S_xi2    = (8.0 * SP_S_temp - 12.0 * SP_S_xi0) * SP_S_temp * SP_S_temp
            SP_S_temp   = xg - SP_S_x0
            SP_S_pC     = 2.0 * SP_S_temp + Gf2 * (1.0 - SP_S_delta1 + SP_S_delta0 - (delta) * (1.0 + SP_S_xi1))
            SP_S_qC     = SP_S_temp * SP_S_temp - Gf2 * (SP_S_delta1 + SP_S_x0 - 1.0 + SP_S_delta0 - (delta) * (SP_S_x0 + 1.0 + SP_S_xi0))
            SP_S_temp   = 2.0 - Gf2 * (SP_S_delta1 + SP_S_delta0 - (delta) * SP_S_xi2)
            SP_S_temp   = SP_S_pC * SP_S_pC - 2.0 * (SP_S_qC * SP_S_temp)
            sp          = SP_S_x0 + 2.0 * (SP_S_qC / (SP_S_pC + sqrt(SP_S_temp)))
        end
    end
    xs = sp
    Es = 0.0
    Ds = 0.0
    Ps = 0.0
    if xg > 0.0
        Es = exp(-xs)
        Ds = abs(1/Es-xs-1-χ(xs))*delta
        Ps = abs(xs-1+Es)
    end
    xgs = xg-xs <= 0 ? xg-xs : Gf*sqrt(Ds+Ps)
    psiss = phit*xs
    return psiss,Es,Ds,Ps,xgs
end

function PSP_psid(Vg,Na,tox,Vfb,T,φb;Vsb=0.0,Vdse=0.0)
    phit = kb*T
    x1 = 1.25
    φb *= 2.0
    xg = (Vg-Vfb)/phit
    G = γ(Na,tox)/sqrt(phit)
    ξ = 1+G/sqrt(2.0)
    xns = (φb+Vsb)/phit
    xnd = (φb+Vsb+Vdse)/phit
    kds = exp(-Vdse/phit)
    Δns = exp(-xns)
    Δnd = Δns*kds
    xmrg = 1e-5*ξ
    if xg > xmrg
        bx= xnd+3.0
        xg1 = x1+G*sqrt(exp(-x1)+x1-1)
        xbar = xg/ξ*(1+xg*(ξ*x1-xg1)/xg1^2)
        x0 = xg+G^2/2.0-G*sqrt(xg+G^2/4-1+exp(-xbar))
        η = MINA(x0,bx,5) - (bx-sqrt(bx^2+5.0))/2.0
        a = (xg-η)^2-G^2*(exp(-η)+η-1-Δnd*(η+1+χ(η)))
        b = 1-G^2/2.0*(exp(-η)-Δnd*χpp(η))
        c = 2*(xg-η)+G^2*(1-exp(-η)-Δnd*(1+χp(η)))
        τ = xnd-η+log(abs(a/G^2))
        y0 = σ2(a,b,c,τ,η)
        Δ0 = exp(y0)
        p = 2*(xg-y0)+G^2*(1-1/Δ0+Δnd*(Δ0-1-χp(y0)))
        q = (xg-y0)^2-G^2*(y0+1/Δ0-1+Δnd*(Δ0-y0-1-χ(y0)))
        xd = y0+2*q/(p+sqrt(p^2-2*q*(2-G^2*(1/Δ0+Δnd*(Δ0-χpp(y0))))))
    elseif xg <= xmrg
        xd = xg/ξ*(1+G*xg*(1-Δnd)/(ξ^2*6*sqrt(2)))
    else
        xg1 = x1+G*sqrt(exp(-x1)+x1-1)
        xbar = xg/ξ*(1+xg*(ξ*x1-xg1)/xg1^2)
        x0 = xg+G^2/2.0-G*sqrt(xg+G^2/4-1+exp(-xbar))
        bx = xns+3
        η=MINA(x0,bx,5.0)-(bx-sqrt(bx^2+5))/2.0
        a = (xg-η)^2-G^2*(exp(-η)+η-1-Δns*(η+1+χ(η)))
        b = 1-G^2/2.0*(exp(-η)-Δns*χpp(η))
        c = 2*(xg-η)+G^2*(1-exp(-η)-Δns*(1+χp(η)))
        τ = xns-η+log(a/G^2)
        y0 = σ2(a,b,c,τ,η)
        Δ0 = exp(y0)
        p = 2*(xg-y0)+G^2*(1-1/Δ0+Δns*(Δ0-1-χp(y0)))
        q = (xg-y0)^2-G^2*(y0+1/Δ0-1+Δns*(Δ0-y0-1-χ(y0)))
        xd = y0+2*q/(p+sqrt(p^2-2*q*(2-G^2*(1/Δ0+Δns*(Δ0-χpp(y0))))))
    end
    ψs,Es,Ds,Ps,xgs = PSP_psiS(Vg,Na,tox,Vfb,T,φb,Vsb)
    xs = ψs*phit
    xds = xd-xs
    if  xds < 1e-10
        p = 2*xgs+G^2*(1-Es+Δnd*(1/Es-1-χp(xs)))
        q = G^2*(1-kds)*Ds
        ξ = 1-G^2/2.0*(Es+Δnd*(1/Es-χpp(xs)))
        xds = 2*q/(p+sqrt(p^2-4*ξ*q))
        xd = xs+xds
    end
    Ed = exp(-xd)
    Dd = (1/Ed-xd-1-χ(xs))*Δnd
    psisd = phit*xd
    return psisd,Ed,Dd,xds
end

function PSP_VaricapPsi(Vg,Na,tox,Vfb,T,φb)
    #psisd,Ed,Dd,Δψ = PSP_psid(Vg,Na,tox,Vfb,T,φb)
    psiss,Es,Ds,Ps,xgs = PSP_psiS(Vg,Na,tox,Vfb,T,φb)
    return psiss
end

function PSP_ψm(Vg,Vd,Vs,Na,tox,Vfb,T,φb)
    psiss,Es,Ds,Ps,xgs = PSP_psiS(Vg,Na,tox,Vfb,T,φb,Vs)
    psisd,Ed,Dd,Pd,xds = PSP_psiS(Vg,Na,tox,Vfb,T,φb,Vd)
    phit = kb*T
    φb *= 2.0
    G = γ(Na,tox)/sqrt(phit)
    x_s = psiss/phit
    x_d = psisd/phit
    xg = (Vg-Vfb)/phit
    xds = x_d-x_s
    if xg > 0.0
        xm = (x_s+x_d)*0.5
        Em = sqrt(Es*Ed)
        Dbar = (Ds+Dd)*0.5
        Dm = Dbar + xds^2/8.0*(Em-2/G^2)
        Pm = xm-1.0+Em
        xgm = G*sqrt(Dm+Pm)

        return xm*phit,Em,Dm,Pm,xds,xgm
    else
        return 0,0,0,0,0,0
    end
end

PSP_qim(G,phit,Dm,Pm,xgm) = G^2*phit*Dm/(xgm+G/sqrt(Pm))
PSP_α(G,Em,Pm) = 1.0 + G*(1-Em)/2/sqrt(Pm)
PSP_qis(G,Ds,xgs,Ps,phit)= G^2*phit*Ds/(xgs+G*sqrt(Ps))
PSP_αs(G,Es,Ps) = 1 + G*(1-Es)/2/sqrt(Ps)
PSP_qbs(phit,G,Ps) = phit*G*sqrt(Ps)
PSP_μx(Vbsx,T,Xcor=0,stxcor=0) = 1+Xcor*Vdsx/(1+0.2*Xcor*Vbsx)
PSP_Eeff(qbs,qis,cox,esi,η=1/2)=1e-8*cox/esi*(qbs+η*qis)
#PSP_ρs(θr,qis) = θr
PSP_θμ(T,STTHEMU=1.5,THEMU=1.5)= THEMU*(T/300.0)^STTHEMU
PSP_μe(T,MUE=0.5,STMUE=0) = MUE*(T/300.0)^STMUE
PSP_Gmob(μe,Eeff,θμ,Cs,qbs,qis,θcs=2.0)=(1+(μe*Eeff)^θμ+Cs*(qbs/(qis+qbs))^θcs)
PSP_μ(μ0,T,Gmob,U0ST) = μ0*(T/300.0)^U0ST/Gmob
PSP_wsat(qis) = 100*qis/(100+qis)
PSP_θsat(T,Gmob,wsat,THESATG=0,THESAT=1,STTHETASAT=1.0) = THESAT*(T/300.0)^STTHETASAT/Gmob*(1+THESATG*wsat)
PSP_ϕ∞(qis,αs,phit)=qis/αs+phit
PSP_ysat(θsat,ϕ∞) = θsat*ϕ∞/sqrt(2)
PSP_za(ysat)=2/(1+sqrt(1+4*ysat))
PSP_ϕ0(ϕ∞,za,ysat) =ϕ∞*za*(1+0.86*za*ysat*(1-za^2*ysat)/(1+4*za^3*ysat^2))
PSP_asat(xgs,G) = xgs+G^2/2.0
PSP_ϕ2(phit,G,Ds,asat) = phit*0.98*G^2*Ds/(asat+sqrt(asat^2-0.98*G^2*Ds))
PSP_ϕsat(ϕ0,ϕ2) = 2*ϕ0*ϕ2/(ϕ0+ϕ2+sqrt((ϕ0+ϕ2)^2-3.95*ϕ0*ϕ2))

function PSP_θsat(Vg,Vs,Na,tox,Vfb,T,φb,Cox)
    psiss,Es,Ds,Ps,xgs = PSP_psiS(Vg,Na,tox,Vfb,T,φb,Vs)
    G = γ(Na,tox)/sqrt(phit)
    phit = kb*T
    qis =  PSP_qis(G,Ds,xgs,Ps,phit)
    αs = PSP_αs(G,Es,Ps)
    qbs = PSP_qbs(phit,G,Ps)
    Eeff = PSP_Eeff(qbs,qis,Cox,ϵ₀*11.7)
    μe = PSP_μe(T)
    θμ = PSP_θμ(T)
    CS = 0
    Gmob = PSP_Gmob(μe,Eff,θμ,CS,qbs,qis)
    wsat = PSP_wsat(qis)
    θsat = PSP_θsat(T,Gmob,wsat)
end

function PSP_Vdsat(Vg,Vs,Na,tox,Vfb,T,φb,Cox)
    psiss,Es,Ds,Ps,xgs = PSP_psiS(Vg,Na,tox,Vfb,T,φb,Vs)
    phit = kb*T
    G = γ(Na,tox)/sqrt(phit)

    qis =  PSP_qis(G,Ds,xgs,Ps,phit)
    αs = PSP_αs(G,Es,Ps)
    qbs = PSP_qbs(phit,G,Ps)
    Eeff = PSP_Eeff(qbs,qis,Cox,ϵ₀*11.7)
    μe = PSP_μe(T)
    θμ = PSP_θμ(T)
    CS = 0
    Gmob = PSP_Gmob(μe,Eeff,θμ,CS,qbs,qis)
    wsat = PSP_wsat(qis)
    θsat = PSP_θsat(T,Gmob,wsat)
    ϕ∞ = PSP_ϕ∞(qis,PSP_αs(G,Es,Ps),phit)
    ysat = PSP_ysat(θsat,ϕ∞)
    za = PSP_za(ysat)
    ϕ0 = PSP_ϕ0(ϕ∞,za,ysat)
    asat = PSP_asat(xgs,G)
    ϕ2 = PSP_ϕ2(phit,G,Ds,asat)
    ϕsat = PSP_ϕsat(ϕ0,ϕ2)
    Vdsat = ϕsat-phit*log(1+(ϕsat*(ϕsat-2*asat*phit))/G^2/Ds/phit/phit)
    return Vdsat
end

function PSP_Vdseff(Vg,Vd,Vs,Na,tox,Vfb,T,φb,Cox)
    Vdsat = PSP_Vdsat(Vg,Vs,Na,tox,Vfb,T,φb,Cox)
    AX = 3.0
    Vds = Vd-Vs
    return Vds/(1+(Vds/Vdsat)^AX)^(1/AX)
end

PSP_Gvsat(θsat,Δψ) = 0.5*(1+sqrt(1+2*(θsat*Δψ)^2))

PSP_T1(Vds,Vdse,Δψ,VP=0.05)=log(abs(1+(Vds-Δψ)/VP)/abs(1+(Vdse-Δψ)/VP))

PSP_ΔLL(T1,ALP=0.01) =ALP*T1

PSP_GΔL(ΔLL) = 1/(1+ΔLL+ΔLL^2)

function PSP_Gvsat(Vg,Vd,Vs,Na,tox,Vfb,T,φb,Cox,CFD,CF)
    Vdseff = PSP_Vdseff(Vg,Vd,Vs,Na,tox,Vfb,T,φb,Cox)
    Vds = Vd-Vs
    Vdsx = Vds^2/(sqrt(Vds^2+0.01)+0.1)
    ΔVgb = CF*2*(sqrt(1.0+CFD*Vdsx)-1.0)/CFD
    Vg = Vg+ΔVgb
    psiss,Es,Ds,Ps,xgs = PSP_psiS(Vg,Na,tox,Vfb,T,φb,Vs)
    psisd,Ed,Dd,Pd,xds = PSP_psiS(Vg,Na,tox,Vfb,T,φb,Vdseff)
    xds = (psisd-psiss)
    phit = kb*T
    G = γ(Na,tox)/sqrt(phit)
    qis =  PSP_qis(G,Ds,xgs,Ps,phit)
    αs = PSP_αs(G,Es,Ps)
    qbs = PSP_qbs(phit,G,Ps)
    Eeff = PSP_Eeff(qbs,qis,Cox,ϵ₀*11.7)
    μe = PSP_μe(T)
    θμ = PSP_θμ(T)
    CS = 0
    Gmob = PSP_Gmob(μe,Eeff,θμ,CS,qbs,qis)
    wsat = PSP_wsat(qis)
    θsat = PSP_θsat(T,Gmob,wsat)
    ϕ∞ = PSP_ϕ∞(qis,PSP_αs(G,Es,Ps),phit)
    ysat = PSP_ysat(θsat,ϕ∞)
    za = PSP_za(ysat)
    ϕ0 = PSP_ϕ0(ϕ∞,za,ysat)
    asat = PSP_asat(xgs,G)
    ϕ2 = PSP_ϕ2(phit,G,Ds,asat)
    ϕsat = PSP_ϕsat(ϕ0,ϕ2)
    return Gmob*PSP_GΔL(PSP_ΔLL(PSP_T1(Vds,Vdseff,xds))),Vdseff,psiss,Es,Ds,Ps,xgs,psisd,Ed,Dd,Pd,xds,G,phit
end

function PSP_ψms(Vg,Vd,Vs,Na,tox,Vfb,T,φb,Cox,CFD,CF)
    Gmob,Vdseff,psiss,Es,Ds,Ps,xgs,psisd,Ed,Dd,Pd,xds,G,phit = PSP_Gvsat(Vg,Vd,Vs,Na,tox,Vfb,T,φb,Cox,CFD,CF)
    φb *= 2.0
    x_s = psiss/phit
    x_d = psisd/phit
    xg = (Vg-Vfb)/phit
    xds = x_d-x_s
    if xg > 0.0
        xm = (x_s+x_d)*0.5
        Em = sqrt(Es*Ed)
        Dbar = (Ds+Dd)*0.5
        Dm = Dbar + xds^2/8.0*(Em-2/G^2)
        Pm = xm-1.0+Em
        xgm = G*sqrt(Dm+Pm)

        return xm*phit,Em,Dm,Pm,xds,xgm,Gmob
    else
        return 0,0,0,0,0,0,0
    end
end

function PSP_Id(Vg,Vd,Vs,Na,tox,Vfb,T,φb,Cox,CFD,CF,BETN=7e-2,U0=5e-2)
    β = BETN*U0*Cox*(T/300.0)^(-1.0)
    phit = kb*T
    G = γ(Na,tox)/sqrt(phit)
    xm,Em,Dm,Pm,xds,xgm,Gmob = PSP_ψms(Vg,Vd,Vs,Na,tox,Vfb,T,φb,Cox,CFD,CF)
    φb *= 2.0
    qim = PSP_qim(G,phit,Dm,Pm,xgm)
    α = PSP_α(G,Em,Pm)
    return β*(qim + phit*α)*phit*xds/Gmob
end

function PSP_Id(Vg,Vd,Vs,MOSf::MOSStructure,CFD,CF,BETN=7e-2,U0=5e-2)
   return PSP_Id(Vg,Vd,Vs,abs(C(MOSf)),MOSf.tox,VFB(MOSf),Temperature(MOSf),ϕb(MOSf),Cox(MOSf),CFD,CF,BETN,U0)
end

function PSP_ID(Vg,Vd,Vs,Na,tox,Vfb,T,φb,Cox)
    phit = kb*T
    G = γ(Na,tox)/sqrt(phit)
    xm,Em,Dm,Pm,xds,xgm = PSP_ψm(Vg,Vd,Vs,Na,tox,Vfb,T,φb)
    φb *= 2.0
    qim = PSP_qim(G,phit,Dm,Pm,xgm)
    α = PSP_α(G,Em,Pm)
    return Cox*(qim + phit*α)*phit*xds
end

function PSP_ID(Vg,Vd,Vs,MOSf::MOSStructure)
   return PSP_ID(Vg,Vd,Vs,abs(C(MOSf)),MOSf.tox,VFB(MOSf),Temperature(MOSf),ϕb(MOSf),Cox(MOSf))
end

function qim_PSP(Vg,Vd,Vs,Na,tox,Vfb,T,φb)
    phit = kb*T
    G = γ(Na,tox)/sqrt(phit)
    xm,Em,Dm,Pm,xds,xgm = PSP_ψm(Vg,Vd,Vs,Na,tox,Vfb,T,φb)
    return qim(G,phit,Dm,Pm,xgm)
end

function α_PSP(Vg,Vd,Vs,Na,tox,Vfb,T,φb)
    phit = kb*T
    G = γ(Na,tox)/sqrt(phit)
    xm,Em,Dm,Pm,xds,xgm = PSP_ψm(Vg,Vd,Vs,Na,tox,Vfb,T,φb)
    return α(G,Em,Pm)

end

ext: jl


text:
@testset "791.custom-sort-string.jl" begin
    @test custom_sort_string("cba", "abcd") == "dcba"
end


text:
struct Label{Name}; end

@pure Label(name::Symbol) = Label{name}()

macro l_str(ex)
    return Label(Symbol(ex::String))
end

@pure (==)(::Label{Name}, ::Label{Name}) where {Name} = true
@pure (==)(::Label{Name1}, ::Label{Name2}) where {Name1, Name2} = false

show(io::IO, ::Label{Name}) where {Name} = print(io, "l\"$(Name)\"")
print(io::IO, ::Label{Name}) where {Name} = print(io, Name)


text:
"""
Computes the parameters of the inverse gamma distribution, given a mean and
standard deviation.
"""
function invgammamoment(m, s)
  v = s ^ 2
  a = (m / s) ^ 2 + 2
  b = m * (a - 1)
  return a, b
end

function generate_sample(; N::Int, gamma::Real, eta::AbstractVector{<:Real},
                         loc::AbstractVector{<:Real},
                         scale::AbstractVector{<:Real},
                         df::AbstractVector{<:Real},
                         skew::AbstractVector{<:Real})
  K = length(eta)
  Ninf= rand(Binomial(N, gamma))
  Nfinite = N - Ninf
  mm = MixtureModel([SkewT(loc[k], scale[k], df[k], skew[k]) for k in 1:K], eta)
  yfinite = rand(mm, Nfinite)
  return [yfinite; fill(-Inf, Ninf)], mm
end

function generate_samples(; NC, NT, gammaC, gammaT, etaC, etaT,
                          loc, scale, df, skew)
  K = length(loc)
  @assert K == length(etaC) == length(etaT)
  @assert K == length(scale) == length(df) == length(skew)
  yC, mmC = generate_sample(N=NC, gamma=gammaC, eta=etaC, loc=loc, scale=scale,
                            df=df, skew=skew)
  yT, mmT = generate_sample(N=NT, gamma=gammaT, eta=etaT, loc=loc, scale=scale,
                            df=df, skew=skew)
  return (yC=yC, yT=yT, gammaC=gammaC, gammaT=gammaT, etaC=etaC, etaT=etaT,
          loc=loc, scale=scale, df=df, skew=skew, mmC=mmC, mmT=mmT)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
# add all packages used here
packages = [
     "IJulia",
     "Gadfly",
     "DataFrames",
     "RCall",
     "TextAnalysis",
     "Query",
     "Vega",
     "Taro",
     "ProgressMeter",
     "Languages",
     "Requests",
     "Plotly",
     "Cairo",
     "Fontconfig"
 ]
# this function installs all required packages
function initApp()
    println("checking for packages...")
    # fetch all installed packages
    installedPkgs = Pkg.installed()
    # loop through each required packages and install of missing
    for pkg in packages
        if !haskey(installedPkgs, pkg)
                println("$pkg is missing in your system...")
                println("Installing $pkg...")
                Pkg.add(pkg)
        end
    end
    println(":) All Packages installed!")
end

# initilize app
initApp()

max_stars_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_issues_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_forks_repo_name: jmulloka/Breast-Cancer-Text-Mining


text:
#!/julia-1.1.1/bin/julia

using TSML
using TSML.ArgumentParsers

tsmlmain()

max_stars_repo_name: ppalmes/TSMLextra

max_issues_repo_name: ppalmes/TSMLextra

max_forks_repo_name: ppalmes/TSMLextra


text:
"""
    p=get_p(mesh,v,point,tet_idx=0)

Get pressure value `p` at point `point`in solution vector `v`.
The optional argument tet_idx accelarates finding the tetrahedron enclosing the point.
"""
function get_p(mesh::Mesh,v,point,tet_idx::Int=0)
    CT, order, smplx = initialize_CT(mesh,v,point,tet_idx)

    if order==:lin
        f=f1

    elseif order==:quad
        f=f2

    elseif order==:herm
        f=fh
    else
        return nothing
    end

    return sum(f(CT,point).*v[smplx])
end
"""
    get_n_grad_p(mesh::Mesh,v,point,n,tet_idx::Int=0)

Get directional derivative `dp` of teh pressure w.r.t. the direction `n` at
point `point`in solution vector `v`. The optional argument tet_idx accelarates
finding the tetrahedron enclosing the point.
"""
function get_n_grad_p(mesh::Mesh,v,point,n,tet_idx::Int=0)
    CT,order,smplx = initialize_CT(mesh,v,point,tet_idx)

    function gradsrc(J,n,x)
        if order==:lin
            return s43nv1rx(J,n,x)
        elseif order==:quad
            return s43nv2rx(J,n,x)
        elseif order==:herm
            return s43nvhrx(J,n,x)
        end
    end

    return sum(gradsrc(CT,n,point).*v[smplx])
end

function initialize_CT(mesh,v,point,tet_idx=0)
    if tet_idx==0
        tet_idx=find_tetrahedron_containing_point(mesh,point)
    end
    dim=length(v) #TODO:Blochify
    N_points=size(mesh.points,2)
    N_lines=length(mesh.lines)
    N_triangles=length(mesh.triangles)+length(mesh.int_triangles)
    if dim == N_points
        el_type=:lin
    elseif dim == N_points+N_lines
        el_type=:quad
    elseif dim == 4*N_points+N_triangles
        el_type=:herm
    else
        el_type=:unknown
    end
    smplx=aggregate_element(mesh,tet_idx, el_type)
    CT=CooTrafo(mesh.points[:,smplx[1:4]])

    return CT, el_type, smplx
end

##

function aggregate_element(mesh::Mesh, tet_idx, el_type=:lin)
    N_points=size(mesh.points)[2]
    if (el_type in (:quad,:herm) ) &&  length(mesh.lines)==0
        collect_lines!(mesh)
    end
    smplx=mesh.tetrahedra[tet_idx]
    if el_type==:lin
        tet=copy(smplx)

    elseif el_type==:quad
        tet=Array{UInt32}(undef,10)
        tet[1:4]=copy(smplx[:])
        tet[5]=get_line_idx(mesh,smplx[[1,2]])+N_points#find_smplx(mesh.lines,smplx[[1,2]])+N_points #TODO: type stability
        tet[6]=get_line_idx(mesh,smplx[[1,3]])+N_points
        tet[7]=get_line_idx(mesh,smplx[[1,4]])+N_points
        tet[8]=get_line_idx(mesh,smplx[[2,3]])+N_points
        tet[9]=get_line_idx(mesh,smplx[[2,4]])+N_points
        tet[10]=get_line_idx(mesh,smplx[[3,4]])+N_points

    elseif el_type==:herm
        tet=Array{UInt32}(undef,20)
        tet[1:4] = smplx[:].+0
        tet[5:8] = smplx[:].+N_points
        tet[9:12] = smplx[:].+2*N_points
        tet[13:16]= smplx[:].+3*N_points
        for (jdx,tria) in enumerate([smplx[[2,3,4]],smplx[[1,3,4]],smplx[[1,2,4]],smplx[[1,2,3]]])
            fcidx     =  find_smplx(mesh.triangles,tria)
            if fcidx !=0
                tet[16+jdx]   =  fcidx+4*N_points
            else
                fcidx=find_smplx(mesh.int_triangles,tria)
                if fcidx==0
                    println("Error, face not found!!!")
                    return nothing
                end
                tet[16+jdx]   =  fcidx+4*N_points+length(mesh.triangles)
            end
        end

    else
        println("Error: element order $(:el_type) not defined!")
        return nothing
    end

    return tet

end


text:
using CUDAdrv
using Base.Test

@test devcount() > 0

include("core.jl")


text:
using LabReports
using Test

@testset "LabReports.jl" begin
    folder = "fake_data"
    reference_folder = "reference_data"
    to_rename = joinpath("fake_data", "200", "200_C&D_3.4e-3")
    renamed = joinpath("fake_data", "200", "200_C&D_3.4e-3_D")

    data = @test_logs (:info, "Renamed $to_rename to $renamed") find_files(folder)
    reference_data = find_files(reference_folder)

    @test length(data) == length(reference_data) == 3
    @testset "Find $k" for k in keys(data)
        for (datafile, ref) in zip(data[k], reference_data[k])
            @test basename(datafile.filename) == basename(ref.filename)
            @test basename(datafile.savename) == basename(ref.savename)
            @test datafile.units == ref.units
            @test datafile.legend_units == ref.legend_units
            @test datafile.idx == ref.idx
        end
    end

    include("originlab.jl")

    ENV["UNITFUL_FANCY_EXPONENTS"] = false

    process_data("CV", data, select=(:Scan, ==, 2))

    process_data("EIS", data, select=(Symbol("-Phase (°)"), >, 0))

    process_data("C&D", data, insert_D=(0, 1.4, 0, 0, 0), continue_col=Symbol("Time (s)"))

    for ((root,dirs,files),(ref_root,ref_dirs,ref_files)) in zip(walkdir(folder), walkdir(reference_folder))
        @test dirs == ref_dirs
        @test length(files) == length(ref_files)
        @test files == ref_files
        @testset "File comparison for $file" for (file, ref_file) in zip(files, ref_files)
            f = read(joinpath(root, file), String)
            r = read(joinpath(ref_root, ref_file), String)
            @test f == r
        end
    end

    @test filevalue(data["CV"][1]) == "14000"
    @test files_with_val(data["CV"], "14000") == [data["CV"][1]]

    vals = Dict("15"=>Set(["2e-3","8.9e-5"]), "200"=>Set(["3.4e-3","8.9e-5","7.4e-2"]))
    @test filevalues(data["C&D"]) == vals

    include("analysis.jl")

    # Cleanup
    clear(folder, r".*\.dat")
    mv(renamed, to_rename)
end


text:
using DelaySSAToolkit
using Catalyst
# using Plots
################################################################################
# the following example illustrates how to deal with a oscillatory system
# d: X-> 0; k1/(1+Y^2): 0-> X; [trigger X-> Y after τ time;] k2*Y/(1+Y):  Y-> 0;

rn = @reaction_network begin
    1/(1+Y^2), 0 --> X
    1/(1+Y),   Y --> 0
end
jumpsys = convert(JumpSystem, rn, combinatoric_ratelaws = false)

u0 = [0,0]
de_chan0 = [[]]
tf = 400.
tspan = (0,tf)
τ = 20.
dprob = DiscreteProblem(jumpsys, u0, tspan)

delay_trigger_affect! = function (integrator, rng)
  append!(integrator.de_chan[1], τ)
end
delay_trigger = Dict(1=>delay_trigger_affect!)
delay_complete = Dict(1=>[2=>1, 1=>-1])
delay_interrupt = Dict()
delayjumpset = DelayJumpSet(delay_trigger, delay_complete, delay_interrupt)

djprob = DelayJumpProblem(jumpsys, dprob, DelayRejection(), delayjumpset, de_chan0, save_positions=(true,true))
sol = solve(djprob, SSAStepper(),seed=1234)

djprob2 = DelayJumpProblem(jumpsys, dprob, DelayDirect(), delayjumpset, de_chan0, save_positions = (true,true))
sol2 = solve(djprob2, SSAStepper(),seed=1234)

plot(sol.t,[sol.u[i][1] for i in 1:length(sol.u)])
plot!(sol.t,[sol.u[i][2] for i in 1:length(sol.u)])
scatter(sol2.t,[sol2.u[i][1] for i in 1:length(sol2.u)],legend = false)
scatter(sol2.t,[sol2.u[i][2] for i in 1:length(sol2.u)],legend = false)


text:
using BenchmarkTools
include("../../src/noise_src_jl/geometry.jl")

# Inputs
x = collect(range(0, 1000, length=100))
y = zeros(100)
z = collect(range(0, 100, length=100))
alpha = 10*ones(100)
gamma = 10*ones(100)
t_s = range(0, 50, length=100)
c_0 = 300*ones(100)
T_0 = 300*ones(100)
x_obs = [6500., 0., 1.2]
n_t = 100

# Settings struct
struct Settings
    dT::Float64
end
settings = Settings(0.)

# Compute geometry
@btime geometry(settings, x_obs, x, y, z, alpha, gamma, t_s, c_0, T_0)

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
"""
Higher-order clustering coefficient, which generalizes the traditional
clustering coefficients to closure of cliques larger than triangles.

Based on "Higher-order clustering in networks" by Hao Yin, Austin R. Benson, and
Jure Leskovec. The original paper may be found at https://arxiv.org/pdf/1704.03913.pdf.
Code from https://github.com/arbenson/HigherOrderClustering.jl.
"""

using LinearAlgebra
using SparseArrays
using StatsBase
using MatrixNetworks
include("cliques.jl")

"""
hoccf_data
----------

This is the data type returned by the function that computes higher-order clustering coefficients.

The field values are

order::Int64
    The order of the clustering coefficient (order = 2 is the classical
    clustering coefficient definition).

global_hoccf::Float64
    The global higher-order clustering coefficient.

avg_hoccf::Float64
    The average higher-order clustering coefficient (the mean is taken over
    nodes at the center of at least one wedge).

avg_hoccf2::Float64
     The average higher-order clustering coefficient, where the local clustering
    of a node not in any wedge is considered to be 0.

local_hoccfs::Vector{Float64}
    The vector of local higher-order clustering coefficients. If a node is not
    the center of at least one wedge, then the value is 0.

ho_wedge_counts::Vector{Int64}
    Higher-order wedge counts of nodes: ho_wedge_counts[v] is the number of higher-order
    wedges with node v at its center.

clique_counts::Vector{Int64}
    Clique counts of nodes: clique_counts[v] is the number of k-cliques containing
    node v where k = order + 1.
"""
struct hoccf_data
    order::Int64
    global_hoccf::Float64
    avg_hoccf::Float64
    avg_hoccf2::Float64
    local_hoccfs::Vector{Float64}
    ho_wedge_counts::Vector{Int64}
    clique_counts::Vector{Int64}
end

"""
`higher_order_ccfs`
======

Computes the global and average local l-th order clustering coefficients of a
network.

Arguments
---------
    - `A::MatrixNetwork`: The network to be analyzed
    - `l::Int64`: The order of the clustering coefficients

The authors define an l-wedge to be an l-clique with an extra node of degree 1 attached.
For example, this is a 3-wedge: https://i.imgur.com/Wc1pRzl.png
The clique node that is adjacent to the extra node is the "center" of the wedge.
The global l-th order clustering coefficient measures the ratio of (l+1)-cliques
to l-wedges: the probability that a randomly chosen l-wedge is "closed" into an
(l+1)-clique. The local l-th order clustering coefficient is defined similarly
as the ratio of the number of (l+1)-cliques a node is part of to the number of
l-wedges it is the center of.

The traditional clustering coefficients are a special case `hocc(A,2)` of the
above concept.
"""
function higher_order_ccfs(A::SparseMatrixCSC{T,Int64}, l::Int64) where T
    A = min.(A, 1)
    A -= Diagonal(A)
    n = size(A, 1)
    # Get clique counts
    clique_counts1 = kcliques(A, l)[1]
    clique_counts2 = kcliques(A, l + 1)[1]
    degs = vec(sum(A, dims=2))
    # normalize
    wedge_counts = (degs .- l .+ 1) .* clique_counts1
    nz_inds = findall(wedge_counts .> 0)
    local_hoccfs = zeros(Float64, n)
    local_hoccfs[nz_inds] = l .* clique_counts2[nz_inds] ./ wedge_counts[nz_inds]
    return hoccf_data(l, l .* sum(clique_counts2) / sum(wedge_counts),
                      mean(local_hoccfs[nz_inds]), mean(local_hoccfs),
                      local_hoccfs, wedge_counts, clique_counts2)
end


text:
import Base.==

"""

A multidimensional array together with additional axes attached to it.
Each axis is a one dimensional array, possibly a `DictArray`.
Use the function `larr` or a macro version `@larr` to create `LabeledArray`s, rather than call the constructor directly.

A `LabeledArray` consists of one main array, which we call the *base* array, and an axis array for each direction.

##### Constructors
* `LabeledArray(arr, axes)` creates a `LabeledArray` from a *base* array `arr` and a tuple `axes` of one dimensional arrays for axes.
* `LabeledArray(arr; kwargs...)` creates a `LabeledArray` from a *base* array `arr` and a keyword argument of the form `axisN=V` for some integer `N` for the axis direction and a one dimensional array `V`. If `N` in `axisN` is omitted, the smallest available number is automatically chosen.

"""
immutable LabeledArray{T,N,AXES<:Tuple,TN} <: AbstractLabeledArray{T,N,AXES,TN}
  data::TN
  axes::AXES

  LabeledArray(data::TN, axes::AXES) = begin
    if length(axes) != ndims(data)
      throw(ArgumentError("LabeledArray: the number of axes $(length(axes)) is not the same as the number of dimensions of data $(ndims(data))."))
    end
    if size(data) != map(length, axes)
      throw(ArgumentError("LabeledArray: the dimensions of data $(size(data)) do not match with the axes lengths $(map(length, axes))."))
    end
    for axis in axes
      if ndims(axis) != 1
        throw(ArgumentError("LabeledArray: axes should be 1 dimensional."))
      end
    end
    new(data, axes)
  end
end

(==)(arr1::LabeledArray, arr2::LabeledArray) = wrap_array(arr1.data)==wrap_array(arr2.data) && all(map(arr1.axes,arr2.axes) do x1,x2;wrap_array(x1)==wrap_array(x2) end)
(==)(::LabeledArray, ::DictArray) = false
(==)(::DictArray, ::LabeledArray) = false

Base.copy(arr::LabeledArray) = LabeledArray(copy(arr.data), map(copy, arr.axes))

"""

Default axis used when no axis is specified for a `LabeledArray`.
It behaves mostly as an array `[Nullable(1), Nullable(2), ...]`.
However, one notable exception is when using `@select`/`selct`/`extract`/`discard`/`getindex`/etc to choose, set or drop specific elements of a `LabeledArray`.
In this case, the result array of reduced size will again have a `DefaultAxis` of an appropriate size.

"""
immutable DefaultAxis <: AbstractArray{Nullable{Int}, 1}
  counts::Int
end
Base.size(axis::DefaultAxis) = (axis.counts,)
Base.linearindexing(::Type{DefaultAxis}) = Base.LinearFast()
Base.getindex(axis::DefaultAxis, i::Int) = Nullable(i)
Base.getindex(axis::DefaultAxis, i...) = begin
  labeledarray_mapnullable_inner((1:axis.counts)[i...])
end
#labeledarray_mapnullable_inner(i::Int) = Nullable(i)
labeledarray_mapnullable_inner(i::AbstractArray) = map(Nullable, i)

Base.start(::DefaultAxis) = 1
Base.next(arr::DefaultAxis, state) = (Nullable(state), state+1)
Base.done(arr::DefaultAxis, s) = s > arr.counts
Base.eltype(::Type{DefaultAxis}) = Nullable{Int}
Base.length(axis::DefaultAxis) = axis.counts
Base.ndims(::DefaultAxis) = 1
(==)(arr1::DefaultAxis, arr2::DefaultAxis) = arr1.counts == arr2.counts

LabeledArray{T,N,AXES<:Tuple}(data::AbstractArray{T,N}, axes::AXES) = begin
  LabeledArray{T,N,AXES,typeof(data)}(data, axes)
end
# don't know why this is necessary.
LabeledArray{K,N,VS,SV,AXES<:Tuple}(data::DictArray{K,N,VS,SV}, axes::AXES) = LabeledArray{LDict{K,SV},N,AXES,typeof(data)}(data, axes)
LabeledArray{T,N}(data::AbstractArray{T,N}; kwargs...) = begin
  axes = Any[DefaultAxis(x) for x in size(data)]
  axis_last_index = 0
  for (k, v) in kwargs
    kstr = string(k)
    if startswith(kstr, "axis")
      axisindex = if length(kstr) == 4
        axis_last_index += 1
        axis_last_index
      else
        axis_last_index = parse(Int, kstr[5:endof(kstr)])
        axis_last_index
      end
      axes[axisindex] = v
    end
  end
  axestuple = (axes...)
  LabeledArray(data, axestuple)
end

# used to calculate the appropriate number of columns/rows when showing an LabeledArray or a DictArray.
arr_width(x::DictArray) = length(x.data)
arr_width(x::AbstractArray) = 1

format_string_map = Dict{DataType,AbstractString}(Float64=>"%0.8g")

"""

`set_format_string!(T, format_string)`

Set the formatting string for type `T` to `format_string`.
By default, `Float64` is displayed with the type format string `"%0.8g"`, and this is the only one.

```julia

julia> @larr(a=rand(3,5),b=rand(1.0*1:10,3,5))
3 x 5 LabeledArray

  |1             |2              |3             |4              |5
--+--------------+---------------+--------------+---------------+--------------
  |a          b  |a          b   |a          b  |a          b   |a          b
--+--------------+---------------+--------------+---------------+--------------
1 |0.47740705 6  |0.29141232 10  |0.83865594 8  |0.99285171 10  |0.70529279 2
2 |0.39800798 4  |0.32770802 2   |0.6787478  3  |0.74177887 2   |0.19838391 9
3 |0.73599273 7  |0.53117299 4   |0.97372752 8  |0.58885829 4   |0.30529039 8

julia> DataCubes.set_format_string!(Float64, "%0.2f")
"%0.2f"

julia> @larr(a=rand(3,5),b=rand(1.0*1:10,3,5))
3 x 5 LabeledArray

  |1         |2         |3         |4          |5
--+----------+----------+----------+-----------+----------
  |a    b    |a    b    |a    b    |a    b     |a    b
--+----------+----------+----------+-----------+----------
1 |0.27 6.00 |0.29 7.00 |0.23 1.00 |0.68 2.00  |0.45 7.00
2 |0.95 6.00 |0.71 5.00 |0.61 1.00 |0.64 10.00 |0.10 7.00
3 |0.05 5.00 |0.52 7.00 |0.11 8.00 |0.52 6.00  |0.21 2.00
```

"""
set_format_string!{T}(::Type{T}, fmt::AbstractString) = (global format_string_map;format_string_map[T] = fmt)

# some functions to set the height and width of output LabeledArrays to console.
default_showsize = () -> ((height,width) = iosize_compat();(max(20,ceil(Int,height*0.5)),width-7))

"""

`set_default_showsize!()`

Set the show height and width limit function to be the default one.
The default version calculates the height and width based on the current console screen size.
This is used to set the print limits when showing a `DictArray` or a `LabeledArray`.

"""
set_default_showsize!() = (global show_size;show_size = default_showsize)
show_size = default_showsize

"""

`set_showsize!(height::Integer, width::Integer)`

Set the show height and width limits to be the constant `height` and `width`, respectively.
To get to the default behavior of adjusting to the current console screen size, use `set_default_showsize!()`.

"""
set_showsize!(height::Integer, width::Integer) = (global show_size;show_size = () -> (height, width))

"""

`set_showheight!(height::Integer)`

Set the show height limit to the constant `height`.
To get to the default behavior of adjusting to the current console screen size, use `set_default_showsize!()`.

"""
set_showheight!(height::Integer) = (origwidth=show_size()[2];global show_size;show_size = () -> (height, origwidth))

"""

`set_showwidth!(width::Integer)`

Set the show width limit to the constant `width`.
To get to the default behavior of adjusting to the current console screen size, use `set_default_showsize!()`.

"""
set_showwidth!(width::Integer) = (origheight=show_size()[1];global show_size;show_size = () -> (origheight, width))
toshow_alongrow = true

"""

`set_showalongrow!(alongrow::Bool)`

Determine how to show the fields of a `DictArray`.
By default (`alongrow` is `true`), the fields are shown from left to right.
If `alongrow=false`, the fields are shown from top to bottom.

##### Examples

```julia
julia> set_showalongrow!(true)
true

julia> darr(a=[1 2 3;4 5 6], b=[:x :y :z;:u :v :w])
2 x 3 DictArray

a b |a b |a b
----+----+----
1 x |2 y |3 z
4 u |5 v |6 w

julia> set_showalongrow!(false)
false

julia> darr(a=[1 2 3;4 5 6], b=[:x :y :z;:u :v :w])
2 x 3 DictArray

a |1 2 3
b |x y z
--+------
a |4 5 6
b |u v w
```

"""
set_showalongrow!(alongrow::Bool) = (global toshow_alongrow; toshow_alongrow = alongrow)

# some functions to set the height and width of output LabeledArrays to html.
default_dispsize = () -> (80,30)

"""

`set_default_dispsize!()`

Set the display height and width limit function to be the default one.
This is used to set the display limits when displaying a `DictArray` or a `LabeledArray`.

"""
set_default_dispsize!() = (global show_size;show_size = default_dispsize)
dispsize = default_dispsize

"""

`set_dispsize!(height::Integer, width::Integer)`

Set the display height and width limits to be the constant `height` and `width`, respectively.
To get to the default behavior, use `set_default_dispsize!()`.

"""
set_dispsize!(height::Integer, width::Integer) = (global dispsize;dispsize = () -> (height, width))

"""

`set_dispheight!(height::Integer)`

Set the display height limit to the constant `height`.
To get to the default behavior, use `set_default_dispsize!()`.

"""
set_dispheight!(height::Integer) = (origwidth=dispsize()[2];global dispsize;dispsize = () -> (height, origwidth))

"""

`set_dispwidth!(width::Integer)`

Set the display width limit to the constant `width`.
To get to the default behavior, use `set_default_dispsize!()`.

"""
set_dispwidth!(width::Integer) = (origheight=dispsize()[1];global dispsize;dispsize = () -> (origheight, width))
todisp_alongrow = true

"""

`set_dispalongrow!(alongrow::Bool)`

Determine how to display the fields of a `DictArray`.
By default (`alongrow` is `true`), the fields are displayed from left to right.
If `alongrow=false`, the fields are displayed from top to bottom.
See `set_showalongrow!` to see an analogous example when showing a `DictArray`.

"""
set_dispalongrow!(alongrow::Bool) = (global todisp_alongrow; todisp_alongrow = alongrow)

"""

`show(io::IO, table::LabeledArray [, indent=0; height=..., width=...])`

Show a LabeledArray.

##### Arguments

* `height` and `width`(optional, default set by show_size()): sets the maximum height and width to draw, beyond which the table will be cut.
* `alongrow`(optional, default set by `set_dispalongrow!`. `tru` by default): if `true`, the fields in the array will be displayed along the row in each cell. Otherwise, they will be stacked on top of each other.

"""
Base.show{N}(io::IO, table::LabeledArray{TypeVar(:T),N}, indent=0; height::Int=show_size()[1], width::Int=show_size()[2], alongrow::Bool=toshow_alongrow) = begin
  print(io, join(size(table), " x "))
  print(io, " LabeledArray")
  print(io, '\n')
  ndimstable = ndims(table)
  for index=1:last(size(table))
    axislabel = table.axes[end][index]
    coords = [fill(:,ndimstable-1);index]
    print(io, '\n')
    show_indent(io, indent)
    print(io, axislabel)
    print(io, ' ')
    print(io, '[')
    for i in 1:length(coords)-1
      print(io, ":,")
    end
    print(io, index)
    print(io, ']')
    print(io, '\n')
    show(io, getindex(table, coords...), indent+4; height=height, width=width, alongrow=alongrow)
  end
end

Base.show(io::IO, table::LabeledArray{TypeVar(:T),0}, indent=0; height::Int=show_size()[1], width::Int=show_size()[2], alongrow::Bool=toshow_alongrow) = begin
  println(io, "0 dimensional LabeledArray")
  show(io, table.data)
end

Base.show(io::IO, table::LabeledArray{TypeVar(:T),1}, indent=0; height::Int=show_size()[1], width::Int=show_size()[2], alongrow::Bool=toshow_alongrow) = begin
  print(io, join(size(table), " x "))
  print(io, " LabeledArray")
  print(io, '\n')
  (result, hlines, vlines) = if alongrow
    create_string_reprmat_alongrow(table; height=height, width=width)
  else
    create_string_reprmat_alongcol(table; height=height, width=width)
  end
  show_string_matrix(io, result, height, width, indent, hlines, vlines)
end

# create a string version of the matrix, properly truncated.
# returns a tuple of (string matrix, locations of horizontal lines, locations of vertical lines).
# location 10 of horizonal line means, e.g., there is a horizontal line to be inserted
# between the 9th and 10th rows.
create_string_reprmat_alongrow(table::LabeledArray{TypeVar(:T),1}; height::Int=show_size()[1], width::Int=show_size()[2]) = begin
  show_height,show_width = (height, width)
  tabledata = table.data
  tableaxes1 = table.axes[1]
  nrows = min(show_height, size(table,1) + 1)
  nkeys = arr_width(tableaxes1)
  ncols = max(nkeys, min(show_width, nkeys + arr_width(tabledata)))
  result = fill("", nrows, ncols)
  if isa(tabledata, DictArray)
    result[1,nkeys+1:ncols] = map(string, tabledata.data.keys[1:ncols-nkeys])
    for i in nkeys+1:ncols
      result[2:nrows,i] = map(cell_to_string, tabledata.data.values[i-nkeys][1:nrows-1])
    end
  else
    result[2:nrows,nkeys+1] = map(cell_to_string, tabledata[1:nrows-1])
  end
  if isa(tableaxes1, DictArray)
    result[1,1:nkeys] = map(string, tableaxes1.data.keys)
    for i in 1:nkeys
      result[2:nrows,i] = map(cell_to_string, tableaxes1.data.values[i][1:nrows-1])
    end
  else
    result[2:nrows,1] = map(cell_to_string, tableaxes1[1:nrows-1])
  end
  (result, [2], [nkeys+1])
end

create_string_reprmat_alongcol(table::LabeledArray{TypeVar(:T),1}; height::Int=show_size()[1], width=show_size()[2]) = begin
  show_height,show_width = (height, width)
  tabledata = table.data
  tableaxes1 = table.axes[1]
  widthtabledata = arr_width(tabledata)
  nrows = min(show_height, 1+widthtabledata*length(table))
  nkeys = arr_width(tableaxes1)
  ncols = max(nkeys, min(show_width, nkeys + 2))
  result = fill("", nrows, ncols)
  if isa(tabledata, DictArray)
    for i in 2:nrows
      result[i,nkeys+1] = string(tabledata.data.keys[1+(i-2) % widthtabledata]) #string(data.keys[i loop_key_index])
      result[i,nkeys+2] = cell_to_string(tabledata.data.values[1+(i-2) % widthtabledata][1+div(i-2,widthtabledata)])
    end
  else
    result[2:nrows,nkeys+1] = map(cell_to_string, tabledata[1:nrows-1])
  end
  # display axis 1.
  if isa(tableaxes1, DictArray)
    result[1,1:nkeys] = map(string, tableaxes1.data.keys)
    for i in 2:nrows
      for j in 1:nkeys
        if (i-2) % widthtabledata == 0
          result[i,j] = cell_to_string(tableaxes1.data.values[j][1+div(i-2,widthtabledata)])
        end
      end
    end
  else
    for i in 2:nrows
      result[i,1] = cell_to_string(tableaxes1[1+div(i-2,widthtabledata)])
    end
  end
  (result, collect(2+widthtabledata*(0:1+fld(nrows,widthtabledata))), [nkeys+1,nkeys+2])
end

Base.show(io::IO, table::LabeledArray{TypeVar(:T),2}, indent; height::Int=show_size()[1], width::Int=show_size()[2], alongrow=toshow_alongrow) = begin
  print(io, join(size(table), " x "))
  print(io, " LabeledArray")
  print(io, '\n')
  (result, hlines, vlines) = if alongrow
    create_string_reprmat_alongrow(table; height=height, width=width)
  else
    create_string_reprmat_alongcol(table; height=height, width=width)
  end
  show_string_matrix(io, result, height, width, indent, hlines, vlines)
end

# create a string version of the matrix, properly truncated.
# returns a tuple of (string matrix, locations of horizontal lines, locations of vertical lines).
# location 10 of horizonal line means, e.g., there is a horizontal line to be inserted
# between the 9th and 10th rows.
create_string_reprmat_alongrow(table::LabeledArray{TypeVar(:T),2}; height::Int=show_size()[1], width::Int=show_size()[2]) = begin
  show_height,show_width = (height, width)
  tabledata = table.data
  (tableheight,tablewidth) = size(tabledata)
  widthtabledata = arr_width(tabledata)
  tableaxes1 = table.axes[1]
  tableaxes2 = table.axes[2]
  widthtableaxes2 = arr_width(tableaxes2)
  nrows = min(show_height, tableheight + widthtableaxes2 + 1)
  nkeys = arr_width(tableaxes1)
  ncols = max(nkeys, min(show_width, nkeys + length(tableaxes2)*widthtabledata))
  result = fill("", nrows, ncols)
  # fill in axis2 above data columns.
  if isa(tableaxes2, DictArray)
    result[1:widthtableaxes2,nkeys] = map(string, tableaxes2.data.keys)
    for i in 1:min(show_height, widthtableaxes2)
      for j in 1:length(tableaxes2)
        ind = 1+nkeys+widthtabledata*(j-1)
        if ind <= ncols
          result[i,ind] = cell_to_string(tableaxes2.data.values[i][j])
        else
          break
        end
      end
    end
  else
    for i in 1:length(tableaxes2)
      ind = 1+nkeys+widthtabledata*(i-1)
      if ind <= ncols
        result[1,ind] = cell_to_string(tableaxes2[i])
      end
    end
  end
  # fill in the key columns.
  if isa(tabledata, DictArray)
    if widthtableaxes2+1 <= nrows
      for i in 1:length(tableaxes2)
        for j in 1:widthtabledata
          ind = nkeys+widthtabledata*(i-1)+j
          if ind>ncols break end
          result[widthtableaxes2+1,ind] = string(tabledata.data.keys[j])
        end
      end
      for i in widthtableaxes2+2:nrows
        for j in 1:length(tableaxes2)
          for k in 1:widthtabledata
            ind = nkeys+widthtabledata*(j-1)+k
            if ind>ncols break end
            result[i,ind] = cell_to_string(tabledata.data.values[k][i-widthtableaxes2-1,j])
          end
        end
      end
    end
  else
    result[widthtableaxes2+2:nrows,nkeys+1:ncols] = map(cell_to_string, tabledata[1:nrows-1-widthtableaxes2,1:ncols-nkeys])
  end
  if isa(tableaxes1, DictArray)
    result[widthtableaxes2+1,1:nkeys] = map(string, tableaxes1.data.keys)
    for i in 1:nkeys
      result[widthtableaxes2+2:nrows,i] = map(cell_to_string, tableaxes1.data.values[i][1:nrows-widthtableaxes2-1])
    end
  else
    result[widthtableaxes2+2:nrows,1] = map(cell_to_string, tableaxes1[1:nrows-1-widthtableaxes2])
  end
  (result, [widthtableaxes2+1,widthtableaxes2+2], map(x->nkeys+1+widthtabledata*(x-1), 1:ncols))
end

create_string_reprmat_alongcol(table::LabeledArray{TypeVar(:T),2}; height::Int=show_size()[1], width::Int=show_size()[2]) = begin
  show_height,show_width = (height, width)
  tabledata = table.data
  (tableheight,tablewidth) = size(tabledata)
  widthtabledata = arr_width(tabledata)
  tableaxes1 = table.axes[1]
  tableaxes2 = table.axes[2]
  widthtableaxes2 = arr_width(tableaxes2)
  nrows = min(show_height, widthtabledata*tableheight + widthtableaxes2)
  nkeys = arr_width(tableaxes1)
  ncols = max(nkeys, min(show_width, nkeys + length(tableaxes2)+1))
  result = fill("", nrows, ncols)
  # fill in axis2 above data columns.
  if isa(tableaxes2, DictArray)
    result[1:widthtableaxes2,nkeys+1] = map(string, tableaxes2.data.keys)
    for i in 1:min(show_height, widthtableaxes2)
      for j in 1:length(tableaxes2)
        ind = 1 + nkeys + j
        if ind <= ncols
          result[i,ind] = cell_to_string(tableaxes2.data.values[i][j])
        else
          break
        end
      end
    end
  else
    for i in 1:length(tableaxes2)
      ind = 1 + nkeys + i
      if ind <= ncols
        result[1,ind] = cell_to_string(tableaxes2[i])
      end
    end
  end
  # fill in the key columns.
  if isa(tabledata, DictArray)
    for j in widthtableaxes2+1:nrows
      result[j,nkeys+1] = string(tabledata.data.keys[1+(j-widthtableaxes2-1) % widthtabledata])
    end
    for i in widthtableaxes2+1:nrows
      for j in 1:length(tableaxes2)
        ind = nkeys + 1 + j
        if ind>ncols break end
        result[i,ind] = cell_to_string(tabledata.data.values[1+(i-widthtableaxes2-1) % widthtabledata][div(i-widthtableaxes2-1, widthtabledata)+1,j])
      end
    end
  else
    result[widthtableaxes2+1:nrows,nkeys+2:ncols] = map(cell_to_string, tabledata[1:nrows-widthtableaxes2,1:ncols-nkeys-1])
  end
  # fill in the axis 1 field names and values.
  if isa(tableaxes1, DictArray)
    result[widthtableaxes2,1:nkeys] = map(string, tableaxes1.data.keys)
    for i in 1:nkeys
      for j in widthtableaxes2+1:nrows
        if (j-widthtableaxes2-1) % widthtabledata == 0
          result[j,i] = cell_to_string(tableaxes1.data.values[i][1+div(j-1-widthtableaxes2, widthtabledata)])
        end
      end
    end
  else
    for j in widthtableaxes2+1:nrows
      if (j-widthtableaxes2-1) % widthtabledata == 0
        result[j,1] = cell_to_string(tableaxes1[1+div(j-1-widthtableaxes2, widthtabledata)])
      end
    end
  end
  (result, collect(widthtableaxes2+1+widthtabledata*(0:1+fld(nrows-1,widthtabledata))), [nkeys+1,nkeys+2])
end

show_indent(io::IO, indent) = print(io, repeat(" ", indent))

show_string_matrix(io::IO, m::Matrix, show_height, show_width, indent, hline_pos=[], vline_pos=[]) = begin
  # calculate maximum width for the displayed cells for each column.
  msize = size(m)
  max_widths = zeros(Int, msize[2])
  for j in 1:msize[1]
    for i in 1:msize[2]
      max_widths[i] = max(max_widths[i], length(m[j,i]))
    end
  end
  # calculate total width in character.
  max_i = 0
  nchar_sofar = indent
  istrunc = false
  for i in 1:length(max_widths)
    if nchar_sofar <= show_width
      nchar_sofar += max_widths[i] + (i in vline_pos ? 2 : 1)
      max_i +=1
    else
      istrunc = true
      break
    end
  end
  j_count = 0
  print(io, '\n')
  for j in 1:msize[1]
    if j in hline_pos
      show_indent(io, indent)
      running_nchars = indent
      for i in 1:max_i
        if i in vline_pos
          print(io, '+')
          running_nchars += 1
        end
        linelength = max(0, min(show_width-running_nchars, max_widths[i]+1))
        print(io, repeat("-", linelength))
        running_nchars += linelength
      end
      print(io, '\n')
      j_count += 1
    end
    show_indent(io, indent)
    running_nchars = indent
    for i in 1:max_i
      if i in vline_pos
        print(io, '|')
        running_nchars += 1
      end
      if max_widths[i] + running_nchars > show_width
        # it could be a unicode character.
        if running_nchars >= show_width-2
          print(io, "..")
          running_nchars += 2
        else
          for c in m[j,i]
            if running_nchars >= show_width-2
              print(io, "..")
              running_nchars += 2
              break
            end
            print(io, c)
            running_nchars += 1
          end
        end
        print(io, repeat(" ", max(0, show_width - running_nchars)))
        running_nchars += max(0, show_width - running_nchars)
      else
        print(io, rpad(m[j,i], max_widths[i]))
        running_nchars += max_widths[i]
      end
      print(io, ' ')
      running_nchars += 1
    end
    if istrunc print(io, "...") end
    print(io, '\n')
    j_count += 1
    if j_count > show_height
      show_indent(io, indent)
      running_nchars = indent
      for i in 1:max_i
        if i in vline_pos
          print(io, '|')
          running_nchars = 1
        end
        expected_running_nchars = running_nchars + max_widths[i] + 1
        if expected_running_nchars >= show_width
          print(io, ":")
          break
        else
          print(io, rpad(":", max_widths[i]))
          running_nchars = expected_running_nchars
          print(io, ' ')
        end
      end
      print(io, '\n')
      break
    end
  end
end

Base.eachindex(arr::LabeledArray) = eachindex(arr.data)
Base.start(arr::LabeledArray) = begin
  iter = eachindex(arr.data)
  i = 1
  axes = map(arr.axes) do axis
    a = BroadcastAxis(axis, arr.data, i)
    i += 1
    a
  end
  (axes, iter, start(iter))
end
Base.next(arr::LabeledArray, state) = begin
  nextelem, nextstate = next(state[2], state[3])
  axesvalues = map(state[1]) do axis
    axis[nextelem]
  end
  nextvalue = arr[nextelem]
  (nextvalue, (state[1], state[2], nextstate))
end
Base.done(arr::LabeledArray, state) = done(state[2], state[3])
Base.linearindexing{T,N,AXES<:Tuple,TN}(::Type{LabeledArray{T,N,AXES,TN}}) = Base.LinearSlow()
Base.eltype{T,N,AXES<:Tuple,TN}(::Type{LabeledArray{T,N,AXES,TN}}) = T
Base.length(table::LabeledArray) = length(table.data)
Base.getindex(arr::LabeledArray, arg::Symbol) = selectfield(arr, arg)
Base.getindex(arr::LabeledArray, arg::Symbol, args::Symbol...) = [selectfield(arr, a) for a in [arg;args...]]
Base.getindex{N}(arr::LabeledArray, args::Tuple{N,Symbol}) = map(a->selectfield(arr, a), args)
Base.getindex(arr::LabeledArray, args::AbstractVector{Symbol}) = selectfields(arr, args...)
Base.getindex(arr::LabeledArray, indices::CartesianIndex) = getindex(arr, indices.I...)
Base.getindex(table::LabeledArray, indices...) = begin
  getindex_labeledarray_nonscalar_indexing(table, indices)
end

getindex_labeledarray_nonscalar_indexing(table::LabeledArray, indices) = begin
  dataelem = getindex(table.data, indices...)
  zippedindices = zip(table.axes, indices)
  filtered_zippedindices = [filter(elem->!isa(elem[2], Real), zip(table.axes, indices))...]
  axeselem = [wrap_scalar_inner(index, axis_index(axis,index)) for (axis, index) in zippedindices]
  ndimstable = ndims(table)
  v04ndims = ndimstable
  for i in ndimstable:-1:1
    if !isa(axeselem[i], Real)
      v04ndims = i
      break
    end
  end
  if ndims(dataelem) == v04ndims
    # this is the expected behavior for v0.4:
    LabeledArray(dataelem, (axeselem[1:v04ndims]...))
  else
    # this is the expected behavior for v0.5, if one direction before the last non integer index is an integer index.
    LabeledArray(dataelem, ([axis_index(axis, index) for (axis, index) in filtered_zippedindices]...))
  end
end

wrap_scalar_inner(index, elem) = elem
wrap_scalar_inner(index::Real, elem) = [elem]

axis_index(axis, index) = axis[index]
axis_sub(axis, index) = sub(axis, index)
axis_slice(axis, index) = slice(axis, index)

axis_index(axis::DefaultAxis, i::Int) = i
axis_index(axis::DefaultAxis, ::Colon) = axis
axis_index(axis::DefaultAxis, index::AbstractArray{Bool}) = DefaultAxis(length(find(index)))
axis_index(axis::DefaultAxis, index) = DefaultAxis(length(index))
axis_sub(axis::DefaultAxis, index) = axis_index(axis, index)
axis_slice(axis::DefaultAxis, index) = axis_index(axis, index)

Base.getindex(table::LabeledArray, indices::Int) =
  getindex_inner(table, ind2sub(table, indices))
Base.getindex(table::LabeledArray, indices::Int...) =
  getindex_inner(table, indices)
getindex_inner{N}(table::LabeledArray, indices::NTuple{N,Int}) = getindex(table.data, indices...)

# returns nested pairs axis1 => (axis2 => ... => (axisN => value)).
getindexpair(table::LabeledArray, indices::Int) =
  getindexpair_inner(table, ind2sub(table, indices))
getindexpair(table::LabeledArray, indices::Int...) =
  getindexpair_inner(table, indices)
getindexpair_inner{N}(table::LabeledArray{TypeVar(:T),N}, indices::NTuple{N,Int}) = begin
  tableaxes = table.axes
  foldr((x,acc)->Pair(x,acc),
        getindex(table.data, indices...),
        [x[y] for (x,y) in zip(tableaxes,indices)])
end

Base.size{N}(table::LabeledArray{TypeVar(:T),N}) = size(table.data)::NTuple{N,Int}
Base.transpose(table::LabeledArray{TypeVar(:T),2}) = LabeledArray(transpose(table.data), (table.axes[2],table.axes[1]))
Base.permutedims(table::LabeledArray, perms) = begin
  lenaxes = length(table.axes)
  if length(perms) != lenaxes
    throw(ArgumentError("LabeledArray permutedims: the number of axes is not the same as the length of the permutations."))
  end
  newdata = permutedims(table.data, perms)
  newaxes = cell(lenaxes)
  for i=1:lenaxes
    newaxes[i] = table.axes[perms[i]]
  end
  LabeledArray(newdata,(newaxes...))
end
Base.setindex!(table::LabeledArray, v, args...) = setindex!(table.data, v, args...)
Base.Multimedia.writemime(io::IO, ::MIME"text/plain", table::LabeledArray) = show(io, table)
Base.sub(table::LabeledArray, indices::Union{Colon, Int64, AbstractArray{TypeVar(:T),1}}...) = begin
  dataelem = sub(table.data, indices...)
  axeselem = [axis_sub(axis, isa(index,Int) ? (index:index) : index) for (axis, index) in zip(table.axes,indices)]
  newN = ndims(dataelem)
  LabeledArray(dataelem, (axeselem[1:newN]...))
end
Base.slice(table::LabeledArray, indices::Union{Colon, Int64, AbstractArray{TypeVar(:T),1}}...) = begin
  dataelem = slice(table.data, indices...)
  axeselem = [axis_slice(axis, index) for (axis, index) in [filter(elem->!isa(elem[2], Real), zip(table.axes, indices))...]]
  newN = ndims(dataelem)
  LabeledArray(dataelem, (axeselem[1:newN]...))
end

# select the field with name from table. If the name belongs to an axis,
# the axis array is broadcast to the full matrix shape.
selectfield(table::LabeledArray, name) = begin
  if isa(table.data, DictArray) && name in table.data.data.keys
    return selectfield(table.data, name)
  end
  counter = 1
  for axis in table.axes
    if isa(axis, DictArray) && name in axis.data.keys
      return AbstractArrayWrapper(BroadcastAxis(selectfield(axis, name), table.data, counter))
    end
    counter += 1
  end
  throw(ArgumentError(string("cannot find field ", name, " in the data or axes.")))
end
selectfields(arr::LabeledArray, fields...) = DictArray(LDict(map(n->n=>selectfield(arr,n), fields)...))

# broadcast an axis so that they can be treated as an N dimensional array for LabeledArray{T,N}.
immutable BroadcastAxis{T,N,V,W} <: AbstractArray{T,N}
  axis::V #AbstractArray{T,1}
  base::W #AbstractArray{TypeVar(:U),N}
  index::Int
end
(==)(arr1::BroadcastAxis, arr2::BroadcastAxis) =
  arr1.index==arr2.index && size(arr1.base)==size(arr2.base) && arr1.axis==arr2.axis
BroadcastAxis{T,N}(axis::AbstractArray{T,1}, base::AbstractArray{TypeVar(:U),N}, index::Int) = BroadcastAxis{T,N,typeof(axis),typeof(base)}(axis, base, index)
BroadcastAxis(axis::DictArray, base::AbstractArray, index::Int) = create_dictarray_nocheck(mapvalues(v->BroadcastAxis(v, base, index), peel(axis)))
Base.size(axis::BroadcastAxis) = size(axis.base)
Base.linearindexing(::Type{BroadcastAxis}) = Base.LinearSlow()
Base.getindex(axis::BroadcastAxis, arg::CartesianIndex) = getindex(axis.axis, arg[axis.index])
Base.getindex(axis::BroadcastAxis, args...) = begin
  newaxis = getindex(axis.axis, args[axis.index])
  if isa(newaxis, AbstractArray)
    newbase = sub(axis.base, args...)
    BroadcastAxis(newaxis, newbase, axis.index)
  else
    # this is definitely not one point. it is one coord along the axis index, but ranges along some other directions.
    # I cannot think of an alternative. Let's do a copy.
    fill(newaxis, size(sub(axis.base, args...)))
  end
end
Base.getindex(axis::BroadcastAxis, index::Int...) = getindex(axis.axis, index[axis.index])
Base.getindex(axis::BroadcastAxis, index::Int) = begin
  getindex(axis.axis, ind2sub(axis.base, index)[axis.index])
end
#getindexvalue{T}(axis::BroadcastAxis, ::Type{T}, args...) = getindexvalue(axis, args...)::T
getindexvalue(axis::BroadcastAxis, arg::CartesianIndex) = getindexvalue(axis.axis, arg[axis.index])
getindexvalue(axis::BroadcastAxis, args...) = begin
  newaxis = getindexvalue(axis.axis, args[axis.index])
  if isa(newaxis, AbstractArray)
    newbase = sub(axis.base, args...)
    BroadcastAxis(newaxis, newbase, axis.index)
  else
    # this is definitely not one point. it is one coord along the axis index, but ranges along some other directions.
    # I cannot think of an alternative. Let's do a copy.
    fill(newaxis, size(sub(axis.base, args...)))
  end
end
getindexvalue(axis::BroadcastAxis, index::Int...) = getindexvalue(axis.axis, index[axis.index])
getindexvalue(axis::BroadcastAxis, index::Int) = begin
  getindexvalue(axis.axis, ind2sub(axis.base, index)[axis.index])
end
# TODO let's make sub/slice for one dimensional indexing.
Base.sub(arr::BroadcastAxis, args::Union{Colon,Int,AbstractVector}...) = BroadcastAxis(arr.axis[args[arr.index]], sub(arr.base, args...), arr.index)
Base.slice(arr::BroadcastAxis, args::Union{Colon,Int,AbstractVector}...) = BroadcastAxis(arr.axis[args[arr.index]], slice(arr.base, args...), arr.index)
Base.sub(arr::BroadcastAxis, args::Tuple{Vararg{Union{Colon,Int,AbstractVector}}})= sub(arr, args...)
Base.slice(arr::BroadcastAxis, args::Tuple{Vararg{Union{Colon,Int,AbstractVector}}}) = slice(arr, args...)

"""

returns all field names for LabeledArray or DictArray. Returns an empty array for other types of arrays.

"""
allfieldnames(table::LabeledArray) = simplify_array(unique(vcat(vcat(map(allfieldnames, table.axes)...),allfieldnames(table.data))))
allfieldnames(table::AbstractArray) = Array(Any, 0)
cell_to_string(cell::Nullable) = cell.isnull ? "" : string(cell.value)
cell_to_string{T}(cell::Nullable{T}) = cell.isnull ? "" : haskey(format_string_map,T) ? sprintf1(format_string_map[T], cell.value) : string(cell.value)
cell_to_string{K,T}(cell::LDict{K,Nullable{T}}) = string(cell)
cell_to_string(cell::AbstractArray) = string(map(x->string(cell_to_string(x)," "), cell)...)

isna(arr::LabeledArray) = isna(peel(arr))
isna(arr::LabeledArray, indices...) = isna(peel(arr), indices...)

"""

`cat(catdim::Integer, arrs::LabeledArray...)`

Concatenate `LabeledArray`s `arrs` along the `catdim` direction.
The base of each element of `arrs` are concatenated and become the new base of the return `LabeledArray`.
The axes of each of `arrs` not along the `catdim` direction should be identical. Otherwise, there will be an error.
The axis along the `catdim` direction will be the concatenation of the axis of each of `arrs` along that direction.

```julia
julia> t1 = larr(a=[1 2 3;4 5 6], axis1=[:x,:y], axis2=["A","B","C"])
t2 x 3 LabeledArray

  |  |A B C
--+--+------
x |a |1 2 3
--+--+------
y |a |4 5 6

julia> t2 = larr(a=[11 12 13;14 15 16], axis1=[:x,:y], axis2=["D","E","F"])
2 x 3 LabeledArray

  |  |D  E  F
--+--+---------
x |a |11 12 13
--+--+---------
y |a |14 15 16

julia> cat(2, t1, t2)
2 x 6 LabeledArray

  |  |A B C D  E  F
--+--+---------------
x |a |1 2 3 11 12 13
--+--+---------------
y |a |4 5 6 14 15 16
```

"""
Base.cat(catdim::Integer, arr1::LabeledArray, arrs::LabeledArray...) = begin
  if isempty(arrs)
    return arr1
  end
  newdata = cat(catdim, arr1.data, map(x->x.data, arrs)...)
  arrlist = Any[arr1,arrs...]
  newaxes = cell(max(map(ndims, arrlist)...))
  for arr in arrlist
    for i in 1:length(arr.axes)
      if i != catdim && isdefined(newaxes, i) && arr.axes[i] != newaxes[i]
        throw(ArgumentError("array dimensions in the arguments for cat do not match."))
      else
        newaxes[i] = arr.axes[i]
      end
    end
  end
  newaxes[catdim] = if all(arr->ndims(arr)<catdim || isa(arr.axes[catdim], DefaultAxis), arrlist) #isa(newaxes[catdim], DefaultAxis)
    DefaultAxis(sum(map(arr->ndims(arr)<catdim ? 1 : length(arr.axes[catdim]), arrlist)))
  else
    maxsizearr = arrlist[indmax(map(ndims, arrlist))]
    null_oneelem_array = create_naelem_array(maxsizearr.axes[catdim])
    cat(1, map(arr->ndims(arr)<catdim ? null_oneelem_array : arr.axes[catdim], arrlist)...)
  end
  LabeledArray(newdata, (newaxes...))
end

create_naelem_array(arr::DictArray) = darr(mapvalues(arr.data) do v
  eltype(v)[Nullable{eltype(eltype(v))}()]
end...)
create_naelem_array(arr::AbstractArray) = eltype(arr)[Nullable{eltype(eltype(arr))}()]

Base.vcat(arr1::LabeledArray, arrs::LabeledArray...) = cat(1, arr1, arrs...)
Base.hcat(arr1::LabeledArray, arrs::LabeledArray...) = cat(2, arr1, arrs...)

Base.repeat(arr::LabeledArray; inner::Array{Int}=ones(Int,ndims(arr)), outer::Array{Int}=ones(Int,ndims(arr))) = begin
  newdata = repeat(arr.data, inner=inner, outer=outer)
  newaxes = ntuple(ndims(arr)) do d
    repeat(arr.axes[d], inner=collect(inner[d]), outer=collect(outer[d]))
  end
  additional_axes = map(n->DefaultAxis(n), size(newdata)[length(newaxes)+1:end])
  LabeledArray(newdata, (newaxes..., additional_axes...))
end

Base.Multimedia.writemime{N}(io::IO,
                             ::MIME"text/html",
                             table::LabeledArray{TypeVar(:T),N};
                             height::Int=dispsize()[1],
                             width::Int=dispsize()[2],
                             alongrow::Bool=todisp_alongrow) = begin
  ndimstable = ndims(table)
  print(io, join(size(table), " x "))
  print(io, " LabeledArray")
  print(io, '\n')
  print(io, "<ul>")
  for index=1:last(size(table))
    axislabel = table.axes[end][index]
    coords = [fill(:,ndimstable-1);index]
    print(io, "<li>")
    print(io, axislabel)
    print(io, ' ')
    print(io, '[')
    for i in 1:length(coords)-1
      print(io, ":,")
    end
    print(io, index)
    print(io, ']')
    print(io, '\n')
    Base.Multimedia.writemime(io, MIME("text/html"), getindex(table, coords...); height=height, width=width, alongrow=alongrow)
    print(io, "</li>")
  end
  print(io, "</ul>")
end

Base.Multimedia.writemime(io::IO, ::MIME"text/html", table::LabeledArray{TypeVar(:T),0}; height::Int=dispsize()[1], width::Int=dispsize()[2], alongrow::Bool=todisp_alongrow) = begin
  print(io, "0 dimensional LabeledArray")
end

Base.Multimedia.writemime(io::IO, ::MIME"text/html", table::Union{LabeledArray{TypeVar(:T),1},LabeledArray{TypeVar(:T),2}}; height::Int=dispsize()[1], width::Int=dispsize()[2], alongrow::Bool=todisp_alongrow) = begin
  print(io, join(size(table), " x "))
  print(io, " LabeledArray")
  print(io, '\n')
  (result, hlines, vlines) = if alongrow
    create_string_reprmat_alongrow(table; height=height, width=width)
  else
    create_string_reprmat_alongcol(table; height=height, width=width)
  end
  print_string_reprmat_to_html_table(io, result, height, width, hlines, vlines)
end

# an internal function to convert a string represented matrix into an html table string printed into an output io).
print_string_reprmat_to_html_table{T<:AbstractString}(io::IO, strrep::AbstractArray{T,2}, height, width, hlines=[], vlines=[]) = begin
  buffer = 1
  print(io, """<table border="0" style="border:none; width:100%; border-collapse:collapse">""")
  for row in 1:size(strrep,1)
    if row > height-buffer
      print(io, """<tr style="border:0">""")
      for col in 1:size(strrep,2)
        if col > width-buffer
          break
        end
        print(io, """<td style="border:1px solid #000; border-width:0 0 0 0">:</td>""")
      end
      print(io, "</tr>")
      break
    end
    print(io, """<tr style="border:0">""")
    for col in 1:size(strrep,2)
      if col > width-buffer
        print(io, """<td style="border:1px solid #000; border-width:0 0 0 0">...</td>""")
        break
      end
      ishline = row in hlines
      isvline = col in vlines
      if ishline && isvline
        print(io, """<td style="border:1px solid #ccc;border-width:1px 0px 0px 1px">""")
      end
      if ishline && !isvline
        print(io, """<td style="border:1px solid #ccc;border-width:1px 0px 0px 0px">""")
      end
      if !ishline && isvline
        print(io, """<td style="border:1px solid #ccc;border-width:0px 0px 0px 1px">""")
      end
      if !ishline && !isvline
        print(io, """<td style="border:1px solid #ccc;border-width:0px 0px 0px 0px">""")
      end
      print(io, strrep[row, col])
      print(io, "</td>")
    end
    print(io, "</tr>")
   end
  print(io, """</table>""")
end

Base.convert(::Type{DictArray}, larr::LabeledArray) = selectfields(larr, allfieldnames(larr)...)
Base.convert(::Type{LabeledArray}, larr::LabeledArray) = larr

Base.map(f::Function, arr::LabeledArray) = mapslices(f, arr, [])
Base.reducedim(f::Function, arr::LabeledArray, dims, initial) = begin
  mapslices(arr, dims) do slice
    reduce(f, initial, slice)
  end
end

"""

`mapslices(f::Function, arr::LabeledArray, dims)`

Apply the function `f` to each slice of `arr` specified by `dims`. `dims` is a vector of integers along which direction to reduce.

* If `dims` includes all dimensions, `f` will be applied to the whole `arr`.
* If `dims` is empty, `mapslices` is the same as `map`.
* Otherwise, `f` is applied to each slice spanned by the directions.

##### Return

Return a dimensionally reduced array along the directions in `dims`.
If the return value of `f` is an `LDict`, the return value of the corresponding `mapslices` is a `DictArray`.
If the return valud of `f` is an AbstractArray, the return value of the corresponding `mapslices` will be an array of the same type, but the dimensions have been extended. If the original dimensions are `N1 x N2 x ... x Nn` and `f` returns an array of dimensions `M1 x ... x Mm`, the return array will be of dimensions `M1 x ... x Mm x N1 x ... x Nn` with those `Ni`s omitted that belong to `dims`. You can think of this as the usual mapslices, followed by `cat`. See the last two examples below.
Otherwise, the return value is an `Array`.

```julia
julia> mapslices(d->d[:a] .* 2,larr(a=[1 2 3;4 5 6],b=[10 11 12;13 14 15],axis1=darr(k=[:X,:Y]),axis2=['A','B','C']),[1])
3 LabeledArray

  |
  --+---------------------------
  A |[Nullable(2),Nullable(8)]
  B |[Nullable(4),Nullable(10)]
  C |[Nullable(6),Nullable(12)]

julia> mapslices(d->d[:a] .* 2,larr(a=[1 2 3;4 5 6],b=[10 11 12;13 14 15],axis1=darr(k=[:X,:Y]),axis2=['A','B','C']),[2])
  2 LabeledArray

  k |
  --+----------------------------------------
  X |[Nullable(2),Nullable(4),Nullable(6)]
  Y |[Nullable(8),Nullable(10),Nullable(12)]

# note how the order of axes changed. the newly generated directions always preced the existing ones.
julia> mapslices(x->msum(x), larr(axis1=[:A,:B,:C,:D,:E],axis2=['X','Y','Z'],a=reshape(1:15,5,3), b=1.0*reshape(1:15,5,3)), [2])
3 x 5 LabeledArray

  |A       |B       |C       |D       |E
--+--------+--------+--------+--------+--------
  |a  b    |a  b    |a  b    |a  b    |a  b
--+--------+--------+--------+--------+--------
X |1  1.0  |2  2.0  |3  3.0  |4  4.0  |5  5.0
Y |7  7.0  |9  9.0  |11 11.0 |13 13.0 |15 15.0
Z |18 18.0 |21 21.0 |24 24.0 |27 27.0 |30 30.0

```

"""
Base.mapslices(f::Function, arr::LabeledArray, dims::AbstractVector) = mapslices_darr_larr(f, arr, dims)
Base.mapslices(f::Function, arr::LabeledArray, dims::Int...) = mapslices_darr_larr(f, arr, [dims...])

recursive_pair_types{K1,K2,V}(::Type{Pair{K1,Pair{K2,V}}}) = DataType[K1, recursive_pair_types(Pair{K2,V})...]
recursive_pair_types{K,V}(::Type{Pair{K,V}}) = DataType[K, V]

"""

`create_dict(::LabeledArray)`

Create a nested `Dict` from a `LabeledArray`.

##### Examples

```julia
julia> t = larr(a=[1 2;3 4], axis1=[:x,:y], axis2=["A","B"])
2 x 2 LabeledArray

  |  |A B
--+--+----
x |a |1 2
--+--+----
y |a |3 4

julia> create_dict(t)
Dict{Nullable{Symbol},Dict{Nullable{ASCIIString},DataCubes.LDict{Symbol,Nullable{Int64}}}} with 2 entries:
  Nullable(:y) => Dict(Nullable("B")=>DataCubes.LDict(:a=>Nullable(4)),Nullable("A")=>DataCubes.LDict(:a=>Nullable(3)))
  Nullable(:x) => Dict(Nullable("B")=>DataCubes.LDict(:a=>Nullable(2)),Nullable("A")=>DataCubes.LDict(:a=>Nullable(1)))
```

"""
function create_dict end

create_dict{T,N,AXES<:Tuple,TN}(arr::LabeledArray{T,N,AXES,TN}) = begin
  kvtype = foldr((x,acc)->Pair{x,acc}, T, map(eltype, arr.axes))
  typevec = recursive_pair_types(kvtype)
  result = foldr((x,acc)->Dict{x,acc}, typevec)()
  for i in eachindex(arr)
    if !isna(arr, i)
      insert_dict!(typevec, result, getindexpair(arr, i))
    end
  end
  result
end

insert_dict!{K,KK,DV,PV}(typevec, d::Dict{K,Dict{KK,DV}}, elem::Pair{K,PV}) = begin
  k = elem.first
  v = elem.second
  if !haskey(d, k)
    d[k] = Dict{KK,DV}()
  end
  insert_dict!(typevec[2:end], d[k], v)
end
insert_dict!{K,V}(typevec, d::Dict{K,V}, elem::Pair{K,V}) = begin
  k = elem.first
  v = elem.second
  d[k] = v
end

# to avoid ambiguity with conventional use of reverse, specifically prohobit the integer argument.
Base.reverse(arr::LabeledArray, dummy::Int) = error("not yet implemented.")
Base.reverse(arr::LabeledArray) = Base.reverse(arr, [1])
# any iterable dims can be an argument.
Base.reverse(arr::LabeledArray, dims) = begin
  coords = ntuple(ndims(arr)) do d
    if d in dims
      size(arr, d):-1:1
    else
      Colon()
    end
  end
  getindex(arr, coords...)
end
Base.flipdim(arr::LabeledArray, dims::Int...) = reverse(arr, dims)

Base.reshape(arr::LabeledArray, dims::Tuple{Vararg{Int}}) = reshape(arr, dims...)

"""

`reshape(arr::LabeledArray, dims...)`

Reshape `arr` into different sizes, if there is no ambiguity.
This means you can collapse several contiguous directions into one direction,
in which case all axes belonging to collapsing direction will be concatenated.
For other case, sometimes it is possible to disambiguate the axis position.
But in general, the result is either an error or an undetermined result
as long as the axis positions are concerned.

##### Examples

```julia
julia> t = larr(a=[1 2 3;4 5 6], axis1=[:x,:y], axis2=["A","B","C"])
2 x 3 LabeledArray

  |A |B |C
--+--+--+--
  |a |a |a
--+--+--+--
x |1 |2 |3
y |4 |5 |6

1 x 6 LabeledArray

x1 |x |y |x |y |x |y
x2 |A |A |B |B |C |C
---+--+--+--+--+--+--
   |a |a |a |a |a |a
---+--+--+--+--+--+--
1  |1 |4 |2 |5 |3 |6

julia> reshape(t, 6, 1)
6 x 1 LabeledArray

      |1
------+--
x1 x2 |a
------+--
x  A  |1
y  A  |4
x  B  |2
y  B  |5
x  C  |3
y  C  |6

julia> reshape(t, 6)
6 LabeledArray

x1 x2 |a
------+--
x  A  |1
y  A  |4
x  B  |2
y  B  |5
x  C  |3
y  C  |6

julia> reshape(t, 3,2)
ERROR: ArgumentError: dims (3,2) are inconsistent.
```

"""
Base.reshape(arr::LabeledArray, dims::Int...) = begin
  newdata = reshape(arr.data, dims...)
  arraxes = arr.axes
  srcdims = [1;cumprod(collect(size(arr)))]
  tgtdims = [1;cumprod(collect(dims))]
  newaxes = Any[DefaultAxis(d) for d in dims]

  tracker = []
  for d in 1:length(arr.axes)
    if !isa(arraxes[d], DefaultAxis)
      startdims = srcdims[d]
      enddims = srcdims[d+1]
      # find where startdims, enddims lie in the tgtdims.
      for p in 1:length(dims)
        tgtdim = tgtdims[p]
        ntgtdim = tgtdims[p+1]
        if tgtdim<=startdims && ntgtdim>=enddims
          inner_dims = divrem(startdims, tgtdim)
          outer_dims = divrem(ntgtdim, enddims)
          if inner_dims[2] != 0 || outer_dims[2] != 0
            throw(ArgumentError("dims $(dims) are inconsistent."))
          end
          axistoadd = repeat(arraxes[d], inner=collect(inner_dims[1]), outer=collect(outer_dims[1]))
          if isa(newaxes[p], DefaultAxis)
            newaxes[p] = axistoadd
          else
            newaxes[p] = DictArray(merge(map(x->isa(x,DictArray) ? x.data : LDict(create_additional_fieldname(arr,tracker)=>x), (newaxes[p], axistoadd))...))
          end
        end
      end
    end
  end
  LabeledArray(newdata, (newaxes...))
end

"""

##### Description

reorders the fields such that the first field names are `names`.
The rest field names are placed sequentially after that.

"""
function reorder end

"""

##### Description

renames the fields such that the first field names are `names`.
The rest field names remain the same.

"""
function rename end

"""

`reorder(arr::LabeledArray, ks...)`

Reorder the field names of the base of `arr` so that the first few field names are `ks`.
The base of `arr` is expected to be a `DictArray`.

##### Return
A new `LabeledArray` whose base fields are shuffled from `arr` so that the first few field names are `ks`.

"""
reorder(arr::LabeledArray, ks...) = LabeledArray(reorder(arr.data, ks...), arr.axes)

"""

`rename(arr::LabeledArray, ks...)`

Rename the field names of the base of `arr` so that the first few field names are `ks`.

##### Return
A new `LabeledArray` whose first few field names of the base of `arr` are `ks`.

"""
rename(arr::LabeledArray, ks...) = LabeledArray(rename(arr.data, ks...), arr.axes)

"""

`@larr(...)`

Create a `LabeledArray`. The arguments `...` can be one of the following:

##### Arguments

* `k=>v` creates a field using array `v` with field name `k` for the base of the return `LabeledArray`. `k` can be an arbitrary type. If the element type of `v` is not `Nullable`, each element will be wrapped by `Nullable`. If `NA` is provided as an element, it is translated as `Nullable{T}()` for an appropriate type `T`.
* `k=v` creates a field using array `v` of the base with field name `:k`.
* There can be at most one non pair type argument, which will be converted into a `LabeledArray` and other pair arguments will update it.
* `axisN[...]` for some integer `N`: this creates an axis along the `N`th direction. If `...` are either keywords or pairs, those are used to create a `DictArray`. Otherwise, an array will be created using `...`. If `N` in `axisN` is omitted, the smallest available number is automatically chosen.

##### Examples

```julia
julia> t = @larr(a=[1 NA;3 4;NA NA],:b=>[1.0 1.5;:sym 'a';"X" "Y"],c=1,axis1[:U,NA,:W],axis[r=['m','n']])
3 x 2 LabeledArray

r |m       |n
--+--------+--------
  |a b   c |a b   c
--+--------+--------
U |1 1.0 1 |  1.5 1
  |3 sym 1 |4 a   1
W |  X   1 |  Y   1

julia> @larr(t, c=[NA NA;3 4;5 6], :d=>:X, axis1[k=["g","h","i"]])
3 x 2 LabeledArray

r |m         |n
--+----------+----------
k |a b   c d |a b   c d
--+----------+----------
g |1 1.0   X |  1.5   X
h |3 sym 3 X |4 a   4 X
i |  X   5 X |  Y   6 X
```

"""
macro larr(args...)
  data_pairs = Any[]
  axes_pairs = Any[]
  template = Nullable()
  for arg in args
    if :head in fieldnames(arg) && (arg.head == :kw || arg.head == :(=>))
      key = arg.head==:kw ? quote_symbol(arg.args[1]) : arg.args[1]
      value = arg.args[2]
      push!(data_pairs, quote
        $key => @nalift($(esc(value)))
      end)
    elseif :head in fieldnames(arg) && arg.head ==:ref && startswith(string(arg.args[1]), "axis") &&
           all(Bool[is_kw(a) for a in arg.args[2:end]])
      push!(axes_pairs, arg.args[1] => read_kws(arg.args[2:end]))
    elseif :head in fieldnames(arg) && arg.head == :ref && startswith(string(arg.args[1]), "axis") && length(arg.args) > 1
      if length(arg.args) != 2
        valueexpr = Expr(:vect, map(esc, arg.args[2:end])...)
        push!(axes_pairs, arg.args[1] => (verbatim_magic_symbol, quote
          @nalift($valueexpr)
        end))
        #throw(ArgumentError("if not a keyword/pair argument, the length has to be 1, consisting of an array for that axis."))
      else
        push!(axes_pairs, arg.args[1] => (verbatim_magic_symbol, quote
          @nalift($(esc(arg.args[2])))
        end))
      end
    elseif template.isnull #length(args) == 1
      template = Nullable(arg)
    else
      throw(ArgumentError("not a key=value type argument or an array argument $arg, or two or more base LabeledArrays are provided."))
    end
  end
  dataexp = if template.isnull && isempty(data_pairs)
    throw(ArgumentError("neither a base LabeledArray or a fieldname=>array pair provided."))
  elseif isempty(data_pairs)
    nothing
  else
    Expr(:call, :darr, Expr(:..., Expr(:vect, data_pairs...)))
  end
  axesexp = map(axes_pairs) do axispair
    axisname = axispair[1]
    axisvalue = axispair[2]
    axisexpr = if isa(axisvalue, Tuple) && length(axisvalue) == 2 && axisvalue[1] == verbatim_magic_symbol
      Expr(:kw, axisname, axisvalue[2])
    else
      Expr(:kw, axisname, Expr(:call, :darr, Expr(:..., Expr(:vect, axisvalue...))))
    end
  end
  if template.isnull
    Expr(:call, :larr, dataexp, axesexp...)
  else
    Expr(:call, :update_larr, quote @nalift($(esc(template.value))) end, dataexp, axesexp...)
  end
end

larr_inner(kwargs...) = begin
  data = []
  axes = []
  for (k,v) in kwargs
    if startswith(string(k), "axis")
      push!(axes, k=>nalift(v))
    else
      push!(data, k=>v)
    end
  end
  LabeledArray(darr(data...); axes...)
end

"""

`larr(...)`

Create a `LabeledArray`. The arguments `...` can be one of the following:

##### Arguments

* `k=>v` creates a field using array `v` with field name `k` for the underlying base `DictArray`. `k` can be an arbitrary type. If the element type of `v` is not `Nullable`, each element will be wrapped by `Nullable`. If you want to manually provide a `Nullable` array with `Nullable{T}()` elements in it, the macro version `@larr` may be more convenient to use. Note that this type of argument precedes the keyword type argument in the return `LabeledArray`, as shown in Examples below.
* `k=v` creates a field using an array `v` with field name `:k` for the underlying base `DictArray`.
* The keyword `axisN=v` for some integer `N` and an array `v` is treated specially. This will create the `N`th axis using the array `v`. If `N` in `axisN` is omitted, the smallest available number is automatically chosen.
* There can be at most one non pair type argument, which will be converted into a `LabeledArray` and other pair arguments will update it.
Especially, if the non pair type argument is an array of `LDict`, it will be converted into a `DictArray`.

##### Examples

```julia
julia> t = larr(a=[1 2;3 4;5 6],:b=>[1.0 1.5;:sym 'a';"X" "Y"],c=1,axis=[:U,:V,:W],axis2=darr(r=['m','n']))
3 x 2 LabeledArray

r |m       |n
--+--------+--------
  |b   a c |b   a c
--+--------+--------
U |1.0 1 1 |1.5 2 1
V |sym 3 1 |a   4 1
W |X   5 1 |Y   6 1

julia> larr(t, c=[1 2;3 4;5 6], :d=>:X, axis1=darr(k=["g","h","i"]))
3 x 2 LabeledArray

r |m         |n
--+----------+----------
k |b   a c d |b   a c d
--+----------+----------
g |1.0 1 1 X |1.5 2 2 X
h |sym 3 3 X |a   4 4 X
i |X   5 5 X |Y   6 6 X
```

"""
function larr end

larr0(arr::LabeledArray, kwargs...) = begin
  data = []
  newaxes = Any[arr.axes...]
  axis_count = 0
  for (k,v) in kwargs
    kstr = string(k)
    if startswith(kstr, "axis")
      axisindex = if length(kstr) == 4
        axis_count += 1
        axis_count
      else
        axis_count = parse(Int, kstr[5:endof(kstr)])
        axis_count
      end
      newaxes[axisindex] = nalift(v)
    else
      push!(data, k=>nalift(v))
    end
  end
  LabeledArray(combine_data_inner(arr.data, LDict(data...)), (newaxes...))
end

combine_data_inner(arr1::DictArray, ldict::LDict) = darr(merge(arr1.data, ldict)...)
combine_data_inner(arr1::AbstractArray, ldict::LDict) = if isempty(ldict)
  arr1
else
  throw(ArgumentError("cannot combine an array of type $(typeof(arr1)) with a dict of type $(typeof(ldict))"))
end

larr(arr::AbstractArray, pairs...; kwargs...) = if isempty(kwargs) && isempty(pairs)
  convert(LabeledArray, nalift(arr))
else
  larr(convert(LabeledArray, nalift(arr)), pairs..., kwargs...)
end
larr(arr::LabeledArray, pairs...; kwargs...) = larr0(arr, pairs..., kwargs...)
larr(arr::LabeledArray) = arr
larr(args...; kwargs...) = larr_inner(args..., kwargs...)

Base.convert{K,V,N}(::Type{LabeledArray}, arr::LabeledArray{LDict{K,V},N}) = LabeledArray(convert_to_dictarray_if_possible(peel(arr)), arr.axes)
Base.convert{K,V,N}(::Type{LabeledArray}, arr::LabeledArray{Nullable{LDict{K,V}},N}) = LabeledArray(convert_to_dictarray_if_possible(peel(arr)), arr.axes)

# for some reason, this generic verion is gtakeed up instead of the LDict specialized versions.
Base.convert(::Type{LabeledArray}, arr::DictArray) = LabeledArray(arr)
Base.convert(::Type{LabeledArray}, arr::AbstractArray) = if eltype(arr) <: LDict
  LabeledArray(convert_to_dictarray_if_possible(arr))
elseif eltype(arr) <: Nullable && eltype(eltype(arr)) <: LDict
  LabeledArray(convert_to_dictarray_if_possible(arr))
else
  LabeledArray(arr)
end

update_larr(base, data; kwargs...) = update_larr(convert(LabeledArray, base), data;kwargs...)
update_larr(base::LabeledArray, data; kwargs...) = begin
  newdata = if data == nothing
    base.data
  else
    darr(merge(base.data.data, data.data)...)
  end
  newaxes = if isempty(kwargs)
    base.axes
  else
    axes = cell(length(base.axes))
    copy!(axes, base.axes)
    for (k,v) in kwargs
      kstr = string(k)
      if startswith(kstr, "axis")
        axisindex = parse(Int, kstr[5:endof(kstr)])
        axes[axisindex] = v
      end
    end
    axes
  end
  LabeledArray(newdata, (newaxes...))
end

"""

`merge(::LabeledArray, ::LabeledArray)`

Merge two `LabeledArrays`. The axes set of the two should be identical.
The bases are merged together and the common axes set is used.

##### Examples

```julia
julia> merge(larr(a=[1,2,3],b=[:x,:y,:z],axis1=[:a,:b,:c]),larr(c=[4,5,6],b=[:m,:n,:p],axis1=[:a,:b,:c]))
3 LabeledArray

  |a b c
--+------
a |1 m 4
b |2 n 5
c |3 p 6
```

"""
Base.merge(arr1::LabeledArray, arr2::LabeledArray) = begin
  @assert(pickaxis(arr1) == pickaxis(arr2))
  LabeledArray(merge(peel(arr1), peel(arr2)), pickaxis(arr1))
end
Base.merge(arr1::LabeledArray, arr2::LabeledArray, args...;kwargs...) = begin
  merge(merge(arr1, arr2), args...;kwargs...)
end
Base.merge(arr::LabeledArray) = arr

"""

`merge(::LabeledArray, ::DictArray...)`
`merge(::LabeledArray, args...)`

Merge the base of the `LabeledArray` and either the rest `DictArray`s, or a `DictArray` created by `args...`.
Together with the axes set of the input `LabeledArray`, return a new `LabeledArray`.

##### Examples

```julia
julia> merge(larr(a=[1,2,3],b=[:x,:y,:z],axis1=[:a,:b,:c]),darr(c=[4,5,6],b=[:m,:n,:p]),darr(a=["X","Y","Z"]))
3 LabeledArray

  |a b c
--+------
a |X m 4
b |Y n 5
c |Z p 6
```

"""
Base.merge(arr1::LabeledArray, arr2::DictArray) = LabeledArray(merge(peel(arr1), arr2), pickaxis(arr1))
Base.merge(arr1::DictArray, arr2::LabeledArray) = LabeledArray(merge(arr1, peel(arr2)), pickaxis(arr2))
Base.merge(arr1::LabeledArray, args...;kwargs...) = merge(arr1, darr(TableSize(arr1), args...;kwargs...))
Base.merge(arr1::LabeledArray, arr2::DictArray, args...;kwargs...) = merge(merge(arr1, arr2), args...;kwargs...)
Base.merge(arr1::DictArray, arr2::LabeledArray, args...;kwargs...) = merge(merge(arr1, arr2), args...;kwargs...)

Base.similar(arr::LabeledArray) = LabeledArray(similar(peel(arr)), arr.axes)
Base.similar{T}(arr::LabeledArray, ::Type{T}) = LabeledArray(similar(peel(arr), T), arr.axes)
Base.similar(arr::LabeledArray, args...) = similar(peel(arr), args...)


text:
module Geomechanics

using RecipesBase
using Plots
using Colors
using DataFrames

export C, Shmin, Shmax, FaultStress

include("Zobackogram.jl")
include("TriplotRecipe.jl")
include("FaultStress.jl")
include("Mohr3.jl")

end # module


text:
# utilility functions for formatting API responses

function to_table() end
function pivot_table() end
function order_table() end


text:
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))
using3D()

function crazy_towel(x, p, n)
    x1, x2, x3, x4, x5 = x
    SVector(
        3.8*x1*(1-x1) - 0.05*(x2+0.35)*(1-2*x3),
        0.1*( (x2+0.35)*(1-2*x3) - 1 )*(1 - 1.9*x1),
        3.78*x3*(1-x3) + 0.5*x2*x4,
        2.0*x4*(1-x4*x3) - 0.2x3,
        0.2*( x5^2*(1-2*x1) - 1 )*(1 - 0.9*x4),
    )
end

u0 = -0.01rand(5)
p0 = rand(3)
ds = DiscreteDynamicalSystem(crazy_towel, u0, p0)

tr = trajectory(ds, 2000; Ttr = 10)
x, y, z, w, e = columns(tr)
scatter3D(z, w, e)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
module LinearRationalExpectations

include("linear_rational_expectations.jl")
export LinearRationalExpectationsWs,
    LinearRationalExpectationsResults, first_order_solver!,
    LinearRationalExpectationsOptions, CyclicReductionOptions,
    GeneralizedSchurOptions

include("extended_lyapunov.jl")
export LyapdWs, extended_lyapd!, extended_lyapd_core!, is_stationary

include("variance.jl")
export VarianceWs, compute_variance!, correlation, variance_decomposition,
    variance_decomposition!, autocovariance!, autocorrelation!

end


text:
using BenchmarkTools: @btime

input = open("a1.txt") do file
  [parse(Int32, line) for line in eachline(file)]
end

ascending = @btime sort(input)
numbers = @btime Set(ascending)

function find_2020_complement()
  for (i,n) in enumerate(ascending)
    c = 2020 - n
    if c <= n && ascending[i+1] != c
      break
    end
    if in(c, numbers)
      return (n,c)
    end
  end
end

part1 = @btime find_2020_complement()
show(part1)
show(part1[1]*part1[2])
print("\n")

function find_2020_triple()
  for (i,n) in enumerate(ascending)
    for j in i+1:length(ascending)
      c = 2020 - n - ascending[j]
      if c <= ascending[j] && ascending[j+1] != c
        break
      end
      if in(c, numbers)
        return (n, ascending[j], c)
      end
    end
  end
end

part2 = @btime find_2020_triple()
show(part2)
show(part2[1]*part2[2]*part2[3])

function find_2020_triple_brute()
  for i in 1:length(ascending)
    for j in i+1:length(ascending)
      for k in j+1:length(ascending)
        ascending[i] + ascending[j] + ascending[k] == 2020 && return (ascending[i], ascending[j], ascending[k])
      end
    end
  end
end

@btime find_2020_triple_brute()


text:
################################################################################
# assumes A is a matrix of matrices or factorizations each of which has size d by d
# IDEA: add block matrix inverse, so far only interested in enabling efficient multiply for iterative methods
"""
```
BlockFactorization{T} <: Factorization{T}
```
Wrapper type for a matrix of matrices or factorizations `A` with `eltype(A) = T`,
which acts as if component matrices were concatenated to form a single matrix `B`
with `eltype(B) = T`.
"""
struct BlockFactorization{T, AT, U, V} <: Factorization{T}
    A::AT
    nindices::U
    mindices::V
end
function BlockFactorization(A::AbstractMatrix, nindices, mindices)
    T = block_eltype(A)
    BlockFactorization{T, typeof(A), typeof(nindices), typeof(mindices)}(A, nindices, mindices)
end
const StridedBlockFactorization = BlockFactorization{<:Any, <:Any, <:StepRange, <:StepRange}
const BlockDiagonalFactorization = BlockFactorization{<:Any, <:Diagonal, <:Any, <:Any}

# calculates the element type of the BlockFactorization with blocks given by A
block_eltype(A::AbstractMatrix{<:Number}) = eltype(A)
block_eltype(A::AbstractMatrix{<:AbstractMatrix}) = eltype(eltype(A))
block_eltype(A::AbstractMatrix{<:Factorization}) = eltype(eltype(A))
function block_eltype(A::AbstractMatrix)
    T = Union{} # bottom type
    for Ai in A
        T = promote_type(T, eltype(Ai))
    end
    return T
end

function BlockFactorization(A::AbstractMatrix, di::Int, dj::Int = di)
    BlockFactorization(A, 1:di:di*size(A, 1)+1, 1:dj:dj*size(A, 2)+1)
end

# strided = true assumes that A has strided block indices, i.e. every element has the same size
function BlockFactorization(A::AbstractMatrix; isstrided::Bool = false)
    if isstrided
        di, dj = size(A[1, 1])
        nind, mind = 1:di:di*size(A, 1)+1, 1:dj:dj*size(A, 2)+1
    else
        n, m = size(A)
        nind, mind = zeros(Int, n+1), zeros(Int, m+1)
        nind[1], mind[1] = 1, 1
        for i in 1:n
            ni = size(A[i, 1], 1)
            all(==(ni), (size(A[i, j], 1) for j in 1:m)) || throw(DimensionMismatch("elements in column $i do not all have the same size"))
            nind[i+1] = nind[i] + ni
        end
        for j in 1:m
            mj = size(A[1, j], 2)
            all(==(mj), (size(A[i, j], 2) for i in 1:n)) || throw(DimensionMismatch("elements in row $j do not all have the same size"))
            mind[j+1] = mind[j] + mj
        end
        nind, mind = strided(nind), strided(mind)
    end
    return BlockFactorization(A, nind, mind)
end

# if possible, convert x to a range, otherwise return x unchanged
function strided(x::AbstractVector)
    length(x) < 2 && return x
    d = x[2] - x[1]
    isstrided = true
    for i in 2:length(x)-1
        if !(x[i+1] - x[i] ≈ d)
            isstrided = false
            break
        end
    end
    isstrided ? range(1, step = d, length = length(x)) : x
end

# allowing AT to be Matrix of Vectors would not work well with blockmul!
# better to cast to Matrix of Matrices
function BlockFactorization(A::AbstractMatrix{<:AbstractVector})
    B = [reshape(aij, :, 1) for aij in A]
    BlockFactorization(B)
end

Base.size(B::BlockFactorization, i::Int) = (1 ≤ i ≤ 2) ? size(B)[i] : 1
Base.size(B::BlockFactorization) = B.nindices[end]-1, B.mindices[end]-1
function Base.Matrix(B::BlockFactorization)
    C = zeros(eltype(B), size(B))
    @threads for j in 1:size(B.A, 2)
        jd = B.mindices[j] : B.mindices[j+1]-1
        for i in 1:size(B.A, 1)
            id = B.nindices[i] : B.nindices[i+1]-1
            M = Matrix(B.A[i, j])
            C[id, jd] .= M
        end
    end
    return C
end

# fallback for now
# Base.getindex(B::BlockFactorization, i::Int, j::Int) = Matrix(B)[i, j]
function Base.getindex(B::BlockFactorization, i::Int, j::Int)
    ni = findlast(≤(i), B.nindices)
    ni = isnothing(ni) ? 1 : ni
    nj = findlast(≤(j), B.mindices)
    nj = isnothing(nj) ? 1 : nj
    ri = i - B.nindices[ni] + 1
    rj = j - B.mindices[nj] + 1
    return B.A[ni, nj][ri, rj]
end

# IDEA: more efficient if neccesary
# function Base.getindex(B::StridedBlockFactorization, i::Int, j::Int)
#     println("strided")
#     ni, nj = mod1(i, B.nindices.step), mod1(j, B.mindices.step)
#     ri, rj = rem(i, B.nindices.step)+1, rem(j, B.mindices.step)+1
#     return B.A[ni, nj][ri, rj]
# end

function LinearAlgebra.diag(B::BlockFactorization)
    n = minimum(size(B))
    d = zeros(eltype(B), n)
    @threads for i in 1:n
        d[i] = B[i, i]
    end
    return d
end

############################## matrix multiplication ###########################
function Base.:*(B::BlockFactorization, x::AbstractVector)
    T = promote_type(eltype(B), eltype(x))
    y = zeros(T, size(B, 1))
    mul!(y, B, x)
end
function Base.:*(B::BlockFactorization, X::AbstractMatrix)
    T = promote_type(eltype(B), eltype(X))
    Y = zeros(T, size(B, 1), size(X, 2))
    mul!(Y, B, X)
end

function LinearAlgebra.mul!(y::AbstractVector, B::BlockFactorization, x::AbstractVector, α::Real = 1, β::Real = 0)
    xx = [@view x[B.mindices[i] : B.mindices[i+1]-1] for i in 1:length(B.mindices)-1]
    yy = [@view y[B.nindices[i] : B.nindices[i+1]-1] for i in 1:length(B.nindices)-1]
    blockmul!(yy, B.A, xx, α, β)
    return y
end

function LinearAlgebra.mul!(Y::AbstractMatrix, B::BlockFactorization, X::AbstractMatrix, α::Real = 1, β::Real = 0)
    XX = [@view X[B.mindices[i] : B.mindices[i+1]-1, :] for i in 1:length(B.mindices)-1]
    YY = [@view Y[B.nindices[i] : B.nindices[i+1]-1, :] for i in 1:length(B.nindices)-1]
    blockmul!(YY, B.A, XX, α, β)
    return Y
end

# carries out multiplication for general BlockFactorization
function blockmul!(y::AbstractVecOfVecOrMat, G::AbstractMatrix, x::AbstractVecOfVecOrMat, α::Real = 1, β::Real = 0)
    @threads for i in eachindex(y)
        @. y[i] = β * y[i]
        for j in eachindex(x)
            mul!(y[i], G[i, j], x[j], α, 1)
        end
    end
    return y
end

################## specialization for diagonal block factorizations ############
# const DiagonalBlockFactorization{T} = BlockFactorization{T, <:Diagonal}
# carries out multiplication for Diagonal BlockFactorization
function blockmul!(y::AbstractVecOfVecOrMat, G::Diagonal,
                   x::AbstractVecOfVecOrMat, α::Real = 1, β::Real = 0)
    @threads for i in eachindex(y)
        mul!(y[i], G[i, i], x[i], α, β)
    end
    return y
end


text:
using Autocorrelations
using Autocorrelations: Sokal, fake_history
using Test

@testset "Autocorrelations.jl" begin
    for _ in 1:10
        history = fake_history(length=1e5, autocortime=10)
        @test effsize(history) ≈ 1e4 rtol=1e-1
        @test autocortime(history) ≈ 10 rtol=1e-1
    end
end