text:
#update the stomatal conductance for empirical stomatal models
function update_gsw!(clayer::CanopyLayer{FT},
            sm::EmpiricalStomatalModel{FT},
            photo_set::C3ParaSet{FT},
            envir::AirLayer{FT},
            Δt::FT,
			betafac::FT;
            τ::FT = FT(600)
) where {FT<:AbstractFloat}
(
    # unpack values
    @unpack g_sw, n_leaf, ps = clayer;

    # calculate steady state values
    # assume τ = 10 minutes
    for _iLF in 1:n_leaf
        _gsw_ss = max(sm.g0,
                      stomatal_conductance(sm, clayer, envir, betafac, _iLF));
        g_sw[_iLF] += (_gsw_ss - g_sw[_iLF]) / τ * Δt;
    end;

    return nothing
)
end

#update the stomatal conductance for Wang optimal stomatal model
function update_gsw!(clayer::CanopyLayer{FT},
            sm::OSMWang{FT},
            photo_set::C3ParaSet{FT},
            envir::AirLayer{FT},
            Δt::FT;
            τ::FT = FT(1e-6)
) where {FT<:AbstractFloat}
(
    # unpack values
    @unpack An, ec, g_bc, g_bw, g_lw, g_m, g_max, g_min, g_sw, n_leaf, p_sat,
            ps = clayer;
    @unpack p_atm, p_H₂O = envir;

    # update g_sw
    for iLF in 1:n_leaf
        _gsw = g_sw[iLF] .+ FT(0.001);
        _glw = 1 / ( 1/_gsw + 1/g_bw[iLF] );
        _gsc = _gsw / FT(1.6);
        _glc = 1 / ( 1/_gsc + 1/g_m[iLF] + 1/g_bc[iLF] );
        leaf_photosynthesis!(photo_set, ps, envir, GCO₂Mode(), _glc);
        _∂A   = ps.An - An[iLF];
        _e0   = g_lw[iLF] * (p_sat - p_H₂O) / p_atm;
        _e1   = _glw * (p_sat - p_H₂O) / p_atm;
        _∂E   = _e1 - _e0;
        _∂A∂E = _∂A / _∂E;
        _∂Θ∂E = FT(max(0.1, An[iLF]) / max(ec - _e0, 1e-7));

        # ensure that dgsw does not change too rapidly
        _Δgsw = (_∂A∂E - _∂Θ∂E) * τ * Δt;
        if _Δgsw > 0
            _Δgsw = min(_Δgsw, (g_max-g_sw[iLF]) / 4);
        else
            _Δgsw = max(_Δgsw, (g_min-g_sw[iLF]) / 4);
        end;
        g_sw[iLF] += _Δgsw;
    end;

    return nothing
)
end

#version of gas exchange functions from Land
#modified by skipping fluorescence calculations

function gas_exchange_nofluor!(
            photo_set::C3ParaSet{FT},
            canopyi::CanopyLayer{FT},
            envir::AirLayer{FT},
            drive::GswDrive,
            ind::Int
) where {FT<:AbstractFloat}
    @unpack g_max, g_min = canopyi;

    # update the conductances
    if canopyi.g_sw[ind] > g_max
        canopyi.g_sw[ind] = g_max;
    end

    if canopyi.g_sw[ind] < g_min
        canopyi.g_sw[ind] = g_min;
    end

    canopyi.g_lw[ind] = 1 / ( 1 / canopyi.g_sw[ind] +
                              1 / canopyi.g_bw[ind] );
    canopyi.g_sc[ind] = canopyi.g_sw[ind] / FT(1.6);
    canopyi.g_lc[ind] = 1 / ( FT(1.6) / canopyi.g_sw[ind] +
                              1 / canopyi.g_m[ind] +
                              1 / canopyi.g_bc[ind] );

    # update the photosynthetic rates
    leaf_photosynthesis!(photo_set, canopyi.ps, envir, GCO₂Mode(),
                         canopyi.g_lc[ind]);

    canopyi.Ac[ind] = canopyi.ps.Ac;
    canopyi.Aj[ind] = canopyi.ps.Aj;
    canopyi.Ap[ind] = canopyi.ps.Ap;
    canopyi.Ag[ind] = canopyi.ps.Ag;
    canopyi.An[ind] = canopyi.ps.An;
    try
    	canopyi.ϕs[ind] = canopyi.ps.ϕs;
	catch found_error
        canopyi.φs[ind] = canopyi.ps.φs;
    end
    # update the pressures
    canopyi.p_i[ind] = canopyi.ps.p_i;
    canopyi.p_s[ind] = canopyi.ps.p_s;

    return nothing
end

function gas_exchange_new!(
            photo_set::C3ParaSet{FT},
            canopyi::CanopyLayer{FT},
            envir::AirLayer{FT},
            drive::GswDrive
) where {FT<:AbstractFloat}
    # update the conductances for each "leaf"
    for i in eachindex(canopyi.g_sw)
        canopyi.ps.APAR = canopyi.APAR[i];
        Photosynthesis.leaf_ETR!(photo_set, canopyi.ps);
        gas_exchange_nofluor!(photo_set, canopyi, envir, drive, i);
    end

    return nothing
end

function update_LAI!(node::SPACMono{FT}, lai::FT) where {FT<:AbstractFloat}
    # 1. Update the general information of the SPAC
    node.la = node.ga * lai;

    # 2. Update the LAI for canopy radiation model
    node.canopy_rt.LAI = lai;

    # 3. Update the LA and LAI for stomatal and photosynthesis models
    for _iPS in node.plant_ps
        _iPS.LA   = node.la / node.n_canopy;
        _iPS.LAI  = lai / node.n_canopy;
        _iPS.tLAI = lai;
    end

    # 4. Update the LA for hydraulic model
    #    TODO abstract this to avoid memory allocation
    for _iHS in node.plant_hs.leaves
        _iHS.area = node.la / node.n_canopy;
    end

    return nothing
end

function update_VJR!(node::SPACMono{FT}, ratio::FT) where {FT<:AbstractFloat}

    # TODO change the ratio accordingly to photo_set
    # TODO add another ratio V2J in photo_set
    # Update Vcmax25, Jmax25 (1.67 Vcmax), and Rd25 (0.015 Vcmax)
    for _iPS in node.plant_ps
        _iPS.ps.Vcmax   = _iPS.ps.Vcmax25WW * ratio;
        _iPS.ps.Vcmax25 = _iPS.ps.Vcmax25WW * ratio;
        _iPS.ps.Jmax    = _iPS.ps.Jmax25WW * ratio;
        _iPS.ps.Jmax25  = _iPS.ps.Jmax25WW * ratio;
        _iPS.ps.Rd      = _iPS.ps.Rd25WW * ratio;
        _iPS.ps.Rd25    = _iPS.ps.Rd25WW * ratio;
    end

    return nothing

end

#the following functions update the soil moisture given the flux into the roots
#while accounting for rain entering the top layer, vertical drainage between the layers,
#and free drainage runoff out of the bottom layer if slope_index is greater than zero
#set slope_index to zero to assume no runoff
#ksoil max is in mm/s
function do_soil_drain!(node::SPACMono{FT}, rain_in::FT, deltaT::FT, k_soil_max::FT, slope_index::FT)
	layer_thick = -1*diff(node.soil_bounds)*1000

	runoff = FT(0);

	for i_root in eachindex(node.plant_hs.roots)
		soil_layerI = node.plant_hs.root_index_in_soil[i_root]

		out_layer = node.plant_hs.roots[i_root].flow /KG_H_2_MOL_S * FT(1/node.ga*deltaT/(60*60));

		node.swc[soil_layerI] -= out_layer / layer_thick[soil_layerI]
		node.swc[soil_layerI] = max(node.plant_hs.roots[soil_layerI].sh.Θr,min(node.swc[soil_layerI],node.plant_hs.roots[soil_layerI].sh.Θs));
	end

	smc_rain = node.swc[1] + rain_in /layer_thick[1]
	smc_rain = max(min(smc_rain, node.plant_hs.roots[1].sh.Θs),node.plant_hs.roots[1].sh.Θr);
	node.swc[1] = smc_rain

	k_soil = zeros(length(node.swc))
	p_soil_tot = zeros(length(node.swc))
	midpoints = zeros(length(node.swc))

	for i_soil in eachindex(node.swc)
		k_soil[i_soil] = soil_k_ratio_swc(node.plant_hs.roots[1].sh, node.swc[i_soil])
		midpoint = (node.soil_bounds[i_soil] + node.soil_bounds[i_soil+1]) /2 * 1000
		midpoints[i_soil] = midpoint;
		p_soil_tot[i_soil] = soil_p_25_swc(node.plant_hs.roots[1].sh, node.swc[i_soil])*mpa2mm + midpoint
	end

	k_soil *= k_soil_max

	p_diff = -1*diff(p_soil_tot) #in mm
	layer_len = -1*diff(midpoints)
	q_down = k_soil[1:length(k_soil)-1] .* p_diff ./ layer_len * deltaT #mm/s * mm / mm * s

	#drainage = k_soil*FT(0.15*6)*deltaT
	drainage = k_soil*slope_index*deltaT
	#drainage = 0

	runoff += drainage[length(k_soil)]

	node.swc[1:(length(k_soil)-1)] -= q_down ./ layer_thick[1:(length(k_soil)-1)]
	node.swc[2:length(k_soil)] += q_down ./ layer_thick[2:length(k_soil)]

	node.swc[length(k_soil)] -= drainage[length(k_soil)] ./ layer_thick[length(k_soil)]
	#node.swc -= drainage ./ layer_thick

	for i_soil in eachindex(node.swc)
		node.swc[i_soil] = max(min(node.swc[i_soil], node.plant_hs.roots[1].sh.Θs),node.plant_hs.roots[1].sh.Θr);
	end

	#return runoff
end

#this function is the same as the previous one, but for use with non steady state plant hydraulics
#because in that case the root water uptake is stored as root.q_in instead of root.flow
function do_soil_nss_drain!(node::SPACMono{FT}, rain_in::FT, deltaT::FT, k_soil_max::FT, slope_index::FT)
	layer_thick = -1*diff(node.soil_bounds)*1000

	runoff = FT(0);

	for i_root in eachindex(node.plant_hs.roots)
		soil_layerI = node.plant_hs.root_index_in_soil[i_root]

		out_layer = node.plant_hs.roots[i_root].q_in /KG_H_2_MOL_S * FT(1/node.ga*deltaT/(60*60));

		node.swc[soil_layerI] -= out_layer / layer_thick[soil_layerI]
		node.swc[soil_layerI] = max(node.plant_hs.roots[i_root].sh.Θr,min(node.swc[soil_layerI],node.plant_hs.roots[i_root].sh.Θs));
	end

	smc_rain = node.swc[1] + rain_in /layer_thick[1]
	smc_rain = max(min(smc_rain, node.plant_hs.roots[1].sh.Θs),node.plant_hs.roots[1].sh.Θr);
	node.swc[1] = smc_rain

	k_soil = zeros(length(node.swc))
	p_soil_tot = zeros(length(node.swc))
	midpoints = zeros(length(node.swc))

	for i_soil in eachindex(node.swc)
		k_soil[i_soil] = soil_k_ratio_swc(node.plant_hs.roots[1].sh, node.swc[i_soil])
		midpoint = (node.soil_bounds[i_soil] + node.soil_bounds[i_soil+1]) /2 * 1000
		midpoints[i_soil] = midpoint;
		p_soil_tot[i_soil] = soil_p_25_swc(node.plant_hs.roots[1].sh, node.swc[i_soil])*mpa2mm + midpoint
	end

	k_soil *= k_soil_max

	p_diff = -1*diff(p_soil_tot) #in mm
	layer_len = -1*diff(midpoints)
	q_down = k_soil[1:length(k_soil)-1] .* p_diff ./ layer_len * deltaT #mm/s * mm / mm * s

	#drainage = k_soil*FT(0.15*6)*deltaT
	drainage = k_soil*slope_index*deltaT
	#drainage = 0

	runoff += drainage[length(k_soil)]

	node.swc[1:(length(k_soil)-1)] -= q_down ./ layer_thick[1:(length(k_soil)-1)]
	node.swc[2:length(k_soil)] += q_down ./ layer_thick[2:length(k_soil)]

	node.swc[length(k_soil)] -= drainage[length(k_soil)] ./ layer_thick[length(k_soil)]
	#node.swc -= drainage ./ layer_thick

	for i_soil in eachindex(node.swc)
		node.swc[i_soil] = max(min(node.swc[i_soil], node.plant_hs.roots[1].sh.Θs),node.plant_hs.roots[1].sh.Θr);
	end

	return runoff #in units of mm
end

function do_soil_nss_drain2!(node::SPACMono{FT}, rain_in::FT, deltaT::FT, k_soil_max::FT, slope_index::FT)
	layer_thick = -1*diff(node.soil_bounds)*1000

	runoff = FT(0);

	for i_root in eachindex(node.plant_hs.roots)
		soil_layerI = node.plant_hs.root_index_in_soil[i_root]

		out_layer = node.plant_hs.roots[i_root].q_in * FT(1/node.ga) * FT(18.02/1000) * deltaT;
        #mol/s over basal area to mol/s over ground area to mm/s to mm

		node.swc[soil_layerI] -= 1*out_layer / layer_thick[soil_layerI]
		node.swc[soil_layerI] = max(node.plant_hs.roots[i_root].sh.Θr,min(node.swc[soil_layerI],node.plant_hs.roots[i_root].sh.Θs));
	end

	smc_rain = node.swc[1] + 1*rain_in/layer_thick[1] #mm/timestep
	smc_rain = max(min(smc_rain, node.plant_hs.roots[1].sh.Θs),node.plant_hs.roots[1].sh.Θr);
	node.swc[1] = smc_rain

	k_soil = zeros(length(node.swc))
	p_soil_tot = zeros(length(node.swc))
	midpoints = zeros(length(node.swc))

	for i_soil in eachindex(node.swc)
		k_soil[i_soil] = soil_k_ratio_swc(node.plant_hs.roots[1].sh, node.swc[i_soil])
		midpoint = (node.soil_bounds[i_soil] + node.soil_bounds[i_soil+1]) /2 * 1000
		midpoints[i_soil] = midpoint;
		p_soil_tot[i_soil] = soil_p_25_swc(node.plant_hs.roots[1].sh, node.swc[i_soil])*mpa2mm + midpoint
        #pressure head in mm
    end

	k_soil *= k_soil_max*1000 #m/s to mm/s

	p_diff = -1*diff(p_soil_tot) #in mm
	layer_len = -1*diff(midpoints)
	q_down = k_soil[1:length(k_soil)-1] .* p_diff ./ layer_len * deltaT #mm/s * mm / mm * s = mm

	#drainage = k_soil*FT(0.15*6)*deltaT
    drainage = k_soil[end]*mpa2mm*(soil_p_25_swc(node.plant_hs.roots[1].sh, node.swc[end])- FT(-0.75))/FT(1000);

	#drainage = 0

	runoff += drainage

	node.swc[1:(length(k_soil)-1)] -= q_down ./ layer_thick[1:(length(k_soil)-1)]
	node.swc[2:length(k_soil)] += q_down ./ layer_thick[2:length(k_soil)]

	node.swc[length(k_soil)] -= 1*drainage*deltaT ./ layer_thick[length(k_soil)]
	#node.swc -= drainage ./ layer_thick

	for i_soil in eachindex(node.swc)
		node.swc[i_soil] = max(min(node.swc[i_soil], node.plant_hs.roots[1].sh.Θs),node.plant_hs.roots[1].sh.Θr);
	end

	return runoff #in units of mm
end

function do_soil_nss_drain3!(node::SPACMono{FT}, rain_in::FT, deltaT::FT, k_soil_max::FT, drainage_pot::FT)
	layer_thick = -1*diff(node.soil_bounds)*1000

	runoff = FT(0);

    layer_flux = zeros(FT,length(layer_thick));
	for i_root in eachindex(node.plant_hs.roots)
		soil_layerI = node.plant_hs.root_index_in_soil[i_root]

		layer_flux[soil_layerI] -= node.plant_hs.roots[i_root].q_in * FT(1/node.ga) * FT(18.02/1000); #* deltaT;
  	end

    layer_flux[1] += rain_in/deltaT;

	k_soil = zeros(length(node.swc))
	p_soil_tot = zeros(length(node.swc))
	midpoints = zeros(length(node.swc))

	for i_soil in eachindex(node.swc)
		k_soil[i_soil] = soil_k_ratio_swc(node.plant_hs.roots[1].sh, node.swc[i_soil])
		midpoint = (node.soil_bounds[i_soil] + node.soil_bounds[i_soil+1]) /2 * 1000
		midpoints[i_soil] = midpoint;
		p_soil_tot[i_soil] = soil_p_25_swc(node.plant_hs.roots[1].sh, node.swc[i_soil])*mpa2mm + midpoint
        #pressure head in mm
    end

	k_soil *= k_soil_max*1000 #m/s to mm/s

	p_diff = -1*diff(p_soil_tot) #in mm
	layer_len = -1*diff(midpoints)
	q_down = k_soil[1:length(k_soil)-1] .* p_diff ./ layer_len;# * deltaT #mm/s * mm / mm * s = mm

	#drainage = k_soil*FT(0.15*6)*deltaT
    drainage = k_soil[end]*mpa2mm*(soil_p_25_swc(node.plant_hs.roots[1].sh, node.swc[end])- drainage_pot)/FT(1000);

	#drainage = 0

	runoff += drainage

	layer_flux[1:(length(k_soil)-1)] -= q_down #./ layer_thick[1:(length(k_soil)-1)]
	layer_flux[2:length(k_soil)] += q_down #./ layer_thick[2:length(k_soil)]

	layer_flux[length(k_soil)] -= 1*drainage #*deltaT ./ layer_thick[length(k_soil)]
	#node.swc -= drainage ./ layer_thick
    node.swc += layer_flux*deltaT ./ layer_thick;

	for i_soil in eachindex(node.swc)
		node.swc[i_soil] = max(min(node.swc[i_soil], node.plant_hs.roots[1].sh.Θs),node.plant_hs.roots[1].sh.Θr);
	end

	return runoff #in units of mm
end


text:
export pmcmc, smc

immutable Step
    thunk::Function
    score::Float64
end

type Particle
    path::Vector{Step}
    value
end

Particle(thunk) = Particle([Step(thunk, 0)], nothing)

type PMCMC <: Ctx
    numparticles::Int64
    thunk::Function
    currentindex::Int64
    retainedparticle::Union(Nothing,Particle)
    particles::Vector{Particle}
    function PMCMC(numparticles, thunk)
        ctx = new(numparticles, thunk, 1, nothing)
        resetparticles!(ctx)
        ctx
    end
end

function resetparticles!(ctx::PMCMC)
    ctx.particles = [Particle(ctx.thunk) for _ in 1:ctx.numparticles]
end

function sample(s::Store, k::Function, address, e::ERP, ::PMCMC)
    k(s,rand(e))
end

function factor(s::Store, k::Function, address, score, ctx::PMCMC)
    push!(ctx.particles[ctx.currentindex].path, Step(()->k(s,nothing), score))
    if ctx.currentindex < length(ctx.particles)
        ctx.currentindex += 1
    else
        @assert length(unique([length(p.path) for p in ctx.particles])) == 1
        ctx.particles = resample(ctx.particles, ctx.retainedparticle)
        ctx.currentindex = 1
    end
    ctx.particles[ctx.currentindex].path[end].thunk()
end

function resample(particles, weights, n)
    @assert length(particles) == length(weights)
    ps = weights/sum(weights)
    Particle[deepcopy(particles[rand(ps)]) for _ in 1:n]
end

function resample(particles, ::Nothing)
    weights = Float64[exp(p.path[end].score) for p in particles]
    resample(particles, weights, length(particles))
end

function resample(particles, retainedparticle::Particle)
    currentstep = length(particles[1].path)
    truncparticle = Particle(retainedparticle.path[1:currentstep], retainedparticle.value)
    allparticles = [particles, truncparticle]
    @assert length(unique([length(p.path) for p in allparticles])) == 1
    weights = Float64[exp(p.path[end].score) for p in allparticles]
    resample(allparticles, weights, length(particles))
end

function pmcmcexit(store::Store, value)
    ctx.particles[ctx.currentindex].value = value
    if ctx.currentindex < length(ctx.particles)
        ctx.currentindex += 1
        ctx.particles[ctx.currentindex].path[end].thunk()
    else
        ctx.currentindex = 1
    end
end

function pmcmc(store::Store, k::Function, address, comp::Function, numiterations, numparticles)
    global ctx
    counts = Dict{Any,Int64}()
    ctxold, ctx = ctx, PMCMC(numparticles, ()->comp(store,pmcmcexit,address))
    try
        for i in 1:numiterations
            # Check some invariants.
            @assert ctx.currentindex == 1
            @assert (i==1?==:!=)(ctx.retainedparticle, nothing)
            @assert all(p->length(p.path)==1 && p.value == nothing, ctx.particles)
            trampoline(ctx.particles[ctx.currentindex].path[end].thunk)
            ctx.retainedparticle = ctx.particles[1]
            for p in ctx.particles
                counts[p.value] = get(counts, p.value, 0) + 1
            end
            resetparticles!(ctx)
        end
    finally
        ctx = ctxold
    end
    k(store, Discrete(counts))
end

# PMCMC uses plain SMC (i.e. no retained particle) for the first
# iteration.
smc(store::Store, k::Function, a, comp::Function, n) = pmcmc(store,k,a,comp,1,n)

ext: jl


text:
@testset "Simplicial Complex" begin
    cplx = SimplicialComplex(Simplex(1,2,3), Simplex(2,4), Simplex(4))

    @test celltype(cplx) == Simplex{Int}
    @test length(cells(cplx)) == 3
    cs = cells(cplx,2)
    @test cs !== nothing
    @test length(cs) == 1
    cs = cells(cplx,3)
    @test cs === nothing
    @test size(cplx) == (4,4,1)
    @test size(cplx,-1) == 0
    @test size(cplx,0) == 4
    @test size(cplx,1) == 4
    @test size(cplx,2) == 1
    @test size(cplx,3) == 0
    @test cplx[Simplex(4)] == hash(Simplex(4))
    @test position(cplx, Simplex(4)) == 4
    @test cplx[Simplex(10)] === nothing
    @test position(cplx, Simplex(10)) === nothing

    @test length(simplices(cplx)) == 9
    scitr = simplices(cplx, 1)
    @test length(scitr) == 4
    for s in simplices(cplx, 2)
        @test s == Simplex(1,2,3)
    end

    cplx = SimplicialComplex(Simplex(1,2,3),
                       Simplex(2,4),
                       Simplex(3,4),
                       Simplex(5,4),
                       Simplex(6))

    # boundary
    qch = Chain(2, [1], [hash(Simplex(1,2,3))])
    resch = Chain(1, [-1, 1, 1], map(hash, [Simplex(1,2), Simplex(2,3), Simplex(1,3)]))
    @testset "Complex boundary chain" for (a,b) in zip(boundary(cplx, qch), resch)
        @test a == b
    end

    # ∂∘∂(c) == 0
    @test iszero(boundary(cplx, boundary(cplx, hash(Simplex(1,2,3)), 2)))

    @test convert(Matrix, boundary(cplx, 2, Int)) == [1 -1 1 0 0 0]'

    # coboundary
    push!(cplx, Simplex(1,3,5), recursive=true)
    push!(cplx, Simplex(1,2,3,5), recursive=true)

    qch = Chain(1, [1], [hash(Simplex(1,3))])
    resch = Chain(2, [-1, 1], map(hash,[Simplex(1,3,5), Simplex(1,2,3)]))
    @testset "Complex coboundary chain" for (a,b) in zip(coboundary(cplx, qch), resch)
        @test a == b
    end

    # δ∘δ(c) == 0
    @test iszero(coboundary(cplx, coboundary(cplx, hash(Simplex(1,3)), 1)))

    cplx = SimplicialComplex(Simplex{Int})

    splx = push!(cplx, Simplex(1))
    @test position(cplx, splx[1]) == 1
    @test dim(splx[1]) == 0
    @test sum(size(cplx)) == 1
    @test values(cells(cplx, 0)[1]) == [1]

    splx = push!(cplx, Simplex(1,2))
    @test splx[1] == Simplex(1,2)
    @test dim(splx[1]) == 1
    @test sum(size(cplx)) == 2
    @test sort!(values(cells(cplx, 1)[1])) == [1,2]

    cplx = SimplicialComplex(Simplex{Int})
    splxs = push!(cplx, Simplex(1, 2), recursive=true)
    @test splxs[1] == Simplex(1,2)
    @test sort!(values(splxs[1])) == [1, 2]
    @test splxs[2] == Simplex(2)
    @test values(splxs[2]) == [2]
    @test splxs[3] == Simplex(1)
    @test values(splxs[3]) == [1]
    @test sum(size(cplx)) == 3
    @test size(cplx, 0) == 2
    @test size(cplx, 1) == 1

    splxs = push!(cplx, Simplex(1,3), recursive=true)
    @test splxs[1] == Simplex(1,3)
    @test sort!(values(splxs[1])) == [1,3]
    @test splxs[2] == Simplex(3)
    @test values(splxs[2]) == [3]
    @test sum(size(cplx)) == 5
    @test size(cplx, 0) == 3
    @test size(cplx, 1) == 2

    @test Simplex(1,3) in cplx
    @test first(cells(cplx,1)) in cplx
    @test !(Simplex(2,3) in cplx)
    @test !(Simplex(2,3,3) in cplx)

    io = IOBuffer()
    write(io, cplx)
    seekstart(io)
    cplx2 = read(io, SimplicialComplex{Simplex{Int}})
    @test size(cplx2, 0) == 3
    @test size(cplx2, 1) == 2

    ch = boundary(cplx2, 1, 0, Int)
    @test dim(ch) == -1
end


text:
using CoolTensors
using Test

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


text:
#=
   -x, abs(x), inv(x),
   sqrt(x), invsqrt(x)
   x+y, x-y, x*y, x/y, hypot(x,y)
   muladd(a,b,c), fma(a,b,c)
=#

function signbit{P}(x::ArbFloat{P})
    0 != ccall(@libarb(arb_is_negative), Int, (Ptr{ArbFloat},), &x)
end

for (op,cfunc) in ((:-,:arb_neg), (:abs, :arb_abs), (:sign, :arb_sgn))
  @eval begin
    function ($op){P}(x::ArbFloat{P})
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}), &z, &x)
      z
    end
  end
end

for (op,cfunc) in ((:inv, :arb_inv), (:sqrt, :arb_sqrt), (:invsqrt, :arb_rsqrt))
  @eval begin
    function ($op){P}(x::ArbFloat{P})
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, P)
      z
    end
  end
end

for (op,cfunc) in ((:+,:arb_add), (:-, :arb_sub), (:*, :arb_mul), (:/, :arb_div), (:hypot, :arb_hypot))
  @eval begin
    function ($op){P}(x::ArbFloat{P}, y::ArbFloat{P})
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, &y, P)
      z
    end
    ($op){T<:AbstractFloat,P}(x::ArbFloat{P}, y::T) = ($op)(x, convert(ArbFloat{P}, y))
    ($op){T<:AbstractFloat,P}(x::T, y::ArbFloat{P}) = ($op)(convert(ArbFloat{P}, x), y)
  end
end

for (op,cfunc) in ((:+,:arb_add_si), (:-, :arb_sub_si), (:*, :arb_mul_si), (:/, :arb_div_si))
  @eval begin
    function ($op){P}(x::ArbFloat{P}, y::Int)
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Int, Int), &z, &x, y, P)
      z
    end
  end
end

(+){P}(x::Int, y::ArbFloat{P}) = (+)(y,x)
(-){P}(x::Int, y::ArbFloat{P}) = -((-)(y,x))
(*){P}(x::Int, y::ArbFloat{P}) = (*)(y,x)
(/){P}(x::Int, y::ArbFloat{P}) = (/)(ArbFloat{P}(x),y)

(+){P}(x::ArbFloat{P}, y::Integer) = (+)(x, convert(ArbFloat{P}, y))
(-){P}(x::ArbFloat{P}, y::Integer) = (-)(x, convert(ArbFloat{P}, y))
(*){P}(x::ArbFloat{P}, y::Integer) = (*)(x, convert(ArbFloat{P}, y))
(/){P}(x::ArbFloat{P}, y::Integer) = (/)(x, convert(ArbFloat{P}, y))

(+){P}(x::Integer, y::ArbFloat{P}) = (+)(convert(ArbFloat{P}, x), y)
(-){P}(x::Integer, y::ArbFloat{P}) = -((-)(y, convert(ArbFloat{P}, x)))
(*){P}(x::Integer, y::ArbFloat{P}) = (*)(convert(ArbFloat{P},x), y)
(/){P}(x::Integer, y::ArbFloat{P}) = (/)(convert(ArbFloat{P},x), y)

for (op,cfunc) in ((:addmul,:arb_addmul), (:submul, :arb_submul))
  @eval begin
    function ($op){P}(w::ArbFloat{P}, x::ArbFloat{P}, y::ArbFloat{P})
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &w, &x, &y, P)
      z
    end
  end
end

muladd{P}(a::ArbFloat{P}, b::ArbFloat{P}, c::ArbFloat{P}) = addmul(c,a,b)
mulsub{P}(a::ArbFloat{P}, b::ArbFloat{P}, c::ArbFloat{P}) = addmul(-c,a,b)


text:
# set LD_LIBRARY_PATH for matplotlib:
#
# ENV["LD_LIBRARY_PATH"] = "\$HOME/.julia/v0.6/Conda/deps/usr/lib:\$LD_LIBRARY_PATH"
# ENV["LD_LIBRARY_PATH"] = Pkg.dir() * "/Conda/deps/usr/lib"
# run(`bash -c 'export LD_LIBRARY_PATH=\$HOME/.julia/v0.6/Conda/deps/usr/lib:\$LD_LIBRARY_PATH'`)

# ld_load_path = Pkg.dir() * "/Conda/deps/usr/lib"
# import Base._setenv
# _setenv("LD_LIBRARY_PATH", ld_load_path)

# ENV["MPLBACKEND"]="agg"
# using PyCall
# py"""
# import os
# os.environ["LD_LIBRARY_PATH"]=$(Pkg.dir()) + "/Conda/deps/usr/lib"
# os.putenv("LD_LIBRARY_PATH", $(Pkg.dir()) + "/Conda/deps/usr/lib")
# """
#
# println("Begin testing: LD_LIBRARY_PATH: $(ENV["LD_LIBRARY_PATH"])")
# println("OLA")
#
# do the tests only, if matplotlib is setup correctly!
#
#if haskey(ENV, "LD_LIBRARY_PATH") && contains(ENV["LD_LIBRARY_PATH"], "Conda")

# ENV["PYTHON"]=""

# using Pkg
# Pkg.add("Conda")
# Pkg.update()
# using Conda
# Conda.update()
# Conda.add("matplotlib")
# Pkg.add("PyCall")
# Pkg.build("PyCall")
# Pkg.add("PyPlot");

using SOM
using Test

# Tests are using Iris dataset:
#
using RDatasets

iris = dataset("datasets", "iris")
train = iris[:,1:4]

include("testFuns.jl")

# test hexagonal, rectangular and spherical training:
#
@test testTrain(train, :hexagonal, toroidal = false)
@test testTrain(train, :rectangular, toroidal = false)
@test testTrain(train, :spherical, toroidal = false)

# test toroidal:
#
@test testTrain(train, :rectangular, toroidal = true)

# test normalisation:
#
@test testTrain(train, :rectangular, normType = :minmax)
@test testTrain(train, :rectangular, normType = :zscore)
@test testTrain(train, :rectangular, normType = :none)

# test convertTrainingData:
#
@test testTrain(convert(Array{Float64}, train), :rectangular)

# test kernels:
#
@test testTrain(train, :rectangular, kernel = gaussianKernel)
@test testTrain(train, :rectangular, kernel = bubbleKernel)

# test visual:
#
@test testVisual(train, :rectangular)

# test frequencies:
#
@test testFreqs(train, iris, :Species)

# test plots:
#
if MPL_INSTALLED
    @test testDensityPlot(train, :hexagonal)
    @test testDensityPlot(train, :spherical)

    @test testOtherDensityPlot(train[1:100,:], train[101:150,:])

    @test testClassesPlot(train, iris, :rectangular)
    @test testClassesPlot(train, iris, :hexagonal)
    @test testClassesPlot(train, iris, :spherical)

    # plot to file:
    #
    @test testClassesToFile(train, iris)
end
#
#eof.


text:
fileid = parse(Int64, ARGS[1])
simregexp = Regex("postsim.*_.*$(fileid).*.bson")

using GPMTD
using SparseProbVec

using BSON
using BSON: @load

using RCall
using StatsBase

## for compatibility with the imported types
using Random
using Random123
using PDMats
using LinearAlgebra
using Statistics
using StatsBase
using Dates

files = filter(x -> occursin(simregexp, x), readdir("./postsim"))
nfiles = length(files)

i = 1
for i in 1:nfiles
    file_now = files[i]
    println("starting file: $(file_now)")

    BSON.@load "./postsim/$(file_now)" model sims init mesg
    nsim = length(sims)

    ### Send to R
    R"rm(list=ls())"
    iter = model.state.iter
    accptr =  model.state.accpt ./ model.state.iter
    L = model.L
    n = model.n
    priorinfo = "$(model.prior)"
    initinfo = "$(init)"

    sims_llik = [ sims[i][:llik] for i=1:nsim ]
    # sims_Scounts = permutedims( hcat( [ StatsBase.counts( sims[i][:ζ], 0:L ) for i = 1:nsim ]... ) )

    sims_lam = permutedims(hcat([ exp.(sims[i][:lλ]) for i=1:nsim ]...))
    sims_int_mu = [ sims[i][:intercept].μ for i=1:nsim ]
    sims_int_sig2 = [ sims[i][:intercept].σ2 for i=1:nsim ]

    @rput mesg nsim iter accptr L n priorinfo initinfo sims_llik sims_lam sims_int_mu sims_int_sig2

    R"ls()"
    R"save.image(file=paste0('postsim/mcmc_', $(mesg), '.rda'))"

    cp("./postsim/$(file_now)", "./postsimB/$(file_now)")
    rm("./postsim/$(file_now)")

    println("completed file number: $(i) of $(nfiles) \n\n")
end


text:
import Base: GenericIOBuffer

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tokenize.Tokens.iskeyword, (Tokenize.Tokens.Kind,))
    precompile(Tokenize.Tokens.isliteral, (Tokenize.Tokens.Kind,))
    precompile(Tokenize.Tokens.isoperator, (Tokenize.Tokens.Kind,))
    precompile(Tokenize.Tokens.Token, (Tokenize.Tokens.Kind,Tuple{Int,Int},Tuple{Int,Int},Int,Int,String))
    precompile(Tokenize.Tokens.Token, ())
    precompile(Tokenize.Tokens.kind, (Tokenize.Tokens.Token,))
    precompile(Tokenize.Tokens.startpos, (Tokenize.Tokens.Token,))
    precompile(Tokenize.Tokens.endpos, (Tokenize.Tokens.Token,))
    precompile(Tokenize.Tokens.untokenize, (Tokenize.Tokens.Token,))
    precompile(Tokenize.Tokens.untokenize, (Tokenize.Tokens.RawToken,String))
    precompile(Tokenize.Tokens.untokenize, (Array{Tokenize.Tokens.Token, 1},))
    precompile(Tokenize.Tokens.untokenize, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))

    precompile(Tokenize.Lexers.is_cat_id_start, (Char, Int32,))
    precompile(Tokenize.Lexers.is_identifier_char, (Char,))
    precompile(Tokenize.Lexers.is_identifier_start_char, (Char,))
    precompile(Tokenize.Lexers.peekchar, (GenericIOBuffer{Array{UInt8, 1}},))
    precompile(Tokenize.Lexers.dpeekchar, (GenericIOBuffer{Array{UInt8, 1}},))
    precompile(Tokenize.Lexers.readchar, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.readchar, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.RawToken},))
    precompile(Tokenize.Lexers.next_token, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))

    precompile(Tokenize.Lexers.ishex, (Char,))
    precompile(Tokenize.Lexers.isbinary, (Char,))
    precompile(Tokenize.Lexers.isoctal, (Char,))
    precompile(Tokenize.Lexers.iswhitespace, (Char,))
    precompile(Tokenize.Lexers.Lexer, (String,))
    precompile(Tokenize.Lexers.Lexer, (String,Type{Tokenize.Tokens.Token}))
    precompile(Tokenize.Lexers.Lexer, (String,Type{Tokenize.Tokens.RawToken}))
    precompile(Tokenize.Lexers.Lexer, (GenericIOBuffer{Array{UInt8, 1}},Type{Tokenize.Tokens.Token}))
    precompile(Tokenize.Lexers.Lexer, (GenericIOBuffer{Array{UInt8, 1}},Type{Tokenize.Tokens.RawToken}))
    precompile(Tokenize.Lexers.tokenize, (String,))

    precompile(Tokenize.Lexers.iterate, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.iterate, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.RawToken},))
    precompile(Tokenize.Lexers.iterate, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token}, Bool,))
    precompile(Tokenize.Lexers.iterate, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.RawToken}, Bool,))
    precompile(Tokenize.Lexers.iterate, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token}, Bool,))
    precompile(Tokenize.Lexers.iterate, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.RawToken}, Bool,))
    precompile(Tokenize.Lexers.startpos, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.startpos, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.RawToken},))
    precompile(Tokenize.Lexers.startpos!, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},Int))
    precompile(Tokenize.Lexers.startpos!, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.RawToken},Int))
    precompile(Tokenize.Lexers.start_token!, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.start_token!, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.RawToken},))

    precompile(Tokenize.Lexers.lex_greater, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_prime, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_digit, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},Tokenize.Tokens.Kind))
    precompile(Tokenize.Lexers.lex_identifier, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token}, Char,))
    precompile(Tokenize.Lexers.lex_less, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_forwardslash, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_minus, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_xor, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_equal, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_bar, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_quote, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_plus, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_dot, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_exclaim, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_colon, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_percent, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_comment, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_comment, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},Bool))
    precompile(Tokenize.Lexers.lex_division, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_circumflex, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_backslash, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_star, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))
    precompile(Tokenize.Lexers.lex_amper, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))

    precompile(Tokenize.Lexers.lex_whitespace, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},))

    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token}, Char,))
    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token}, String,))
    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},typeof( Base.isdigit),))
    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},typeof( Tokenize.Lexers.iswhitespace),))
    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},typeof( Tokenize.Lexers.is_identifier_char),))
    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token},typeof(Tokenize.Lexers.ishex),))
    precompile(Tokenize.Lexers.accept_batch, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token}, typeof(Tokenize.Lexers.iswhitespace),))
    precompile(Tokenize.Lexers.accept_batch, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.Token}, typeof(Tokenize.Lexers.isdigit),))

    precompile(Tokenize.Lexers.accept_batch, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.RawToken}, typeof(Tokenize.Lexers.iswhitespace),))
    precompile(Tokenize.Lexers.accept_batch, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.RawToken}, typeof(Tokenize.Lexers.isdigit),))
    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}},Tokenize.Tokens.RawToken}, Char,))

    precompile(Tokenize.Lexers.readchar, (GenericIOBuffer{Array{UInt8, 1}},))
end


text:
module Dash
import HTTP, JSON2
using MacroTools
include("ComponentPackages.jl")
include("ComponentMetas.jl")
include("Components.jl")
include("Front.jl")

import .ComponentPackages
import .Front
using .ComponentMetas
using .Components

export dash, Component, Front, @use, <|, @callid_str, CallbackId, callback!,
run_server, PreventUpdate, no_update, @wildprop

ComponentPackages.@reg_components()
include("app.jl")
include("handlers.jl")
include("utils.jl")

@doc """
    module Dash

Julia backend for [Plotly Dash](https://github.com/plotly/dash)

# Examples
```julia

using Dash
app = dash("Test", external_stylesheets=["https://codepen.io/chriddyp/pen/bWLwgP.css"]) do
    html_div() do
        dcc_input(id="graphTitle", value="Let's Dance!", type = "text"),
        html_div(id="outputID"),
        dcc_graph(id="graph",
            figure = (
                data = [(x = [1,2,3], y = [3,2,8], type="bar")],
                layout = Dict(:title => "Graph")
            )
        )

    end
end
callback!(app, callid"{graphTitle.type} graphTitle.value => outputID.children") do type, value
    "You've entered: '\$(value)' into a '\$(type)' input control"
end
callback!(app, callid"graphTitle.value => graph.figure") do value
    (
        data = [
            (x = [1,2,3], y = abs.(randn(3)), type="bar"),
            (x = [1,2,3], y = abs.(randn(3)), type="scatter", mode = "lines+markers", line = (width = 4,))
        ],
        layout = (title = value,)
    )
end
handle = make_handler(app, debug = true)
run_server(handle, HTTP.Sockets.localhost, 8080)
```

# Available components
$(ComponentPackages.component_doc_list())
""" Dashboards

"""
    run_server(app::DashApp, host = HTTP.Sockets.localhost, port = 8080; debug::Bool = false)

Run Dash server

#Arguments
- `app` - Dash application
- `host` - host
- `port` - port
- `debug::Bool = false` - Enable/disable all the dev tools

#Examples
```jldoctest
julia> app = dash("Test") do
    html_div() do
        html_h1("Test Dashboard")
    end
end
julia>
julia> run_server(handler,  HTTP.Sockets.localhost, 8080)
```

"""
function run_server(app::DashApp, host = HTTP.Sockets.localhost, port = 8080; debug = false)
    handler = make_handler(app, debug = debug);
    HTTP.serve(handler, host, port)
end

end # module


text:
# Quick example

using Revise, BinaryTraits

# Let's define two ability traits
abstract type Ability end
@trait Fly as Ability
@trait Swim as Ability

# Define interface contracts for the type.
@implement CanFly by liftoff()
@implement CanFly by fly(direction::Float64, altitude::Float64)
@implement CanFly by speed()::Float64

# Define a data type and assign it traits.
struct Crane end
@assign Crane with CanFly,CanSwim

# Check compliance.
@check(Crane)
#=
julia> @check(Crane)
┌ Warning: Missing implementation: FlyTrait: CanFly ⇢ speed(::Crane)::Float64
└ @ BinaryTraits ~/.julia/dev/BinaryTraits/src/interface.jl:200
┌ Warning: Missing implementation: FlyTrait: CanFly ⇢ liftoff(::Crane)::Any
└ @ BinaryTraits ~/.julia/dev/BinaryTraits/src/interface.jl:200
┌ Warning: Missing implementation: FlyTrait: CanFly ⇢ fly(::Crane, ::Float64, ::Float64)::Any
└ @ BinaryTraits ~/.julia/dev/BinaryTraits/src/interface.jl:200
❌ Crane is missing these implementations:
1. FlyTrait: CanFly ⇢ speed(::<Type>)::Float64
2. FlyTrait: CanFly ⇢ liftoff(::<Type>)::Any
3. FlyTrait: CanFly ⇢ fly(::<Type>, ::Float64, ::Float64)::Any
=#

# What about composite traits?
@trait FlySwim with CanFly,CanSwim

# Define a new type and assign with composite trait
struct Swan end
@assign Swan with CanFlySwim

# Check compliance. It should drill down to figure out required interface contracts.
@check(Swan)
#=
julia> @check(Swan)
┌ Warning: Missing implementation: FlyTrait: CanFly ⇢ fly(::Swan, ::Float64, ::Float64)::Any
└ @ BinaryTraits ~/.julia/dev/BinaryTraits/src/interface.jl:77
┌ Warning: Missing implementation: FlyTrait: CanFly ⇢ speed(::Swan)::Float64
└ @ BinaryTraits ~/.julia/dev/BinaryTraits/src/interface.jl:77
┌ Warning: Missing implementation: FlyTrait: CanFly ⇢ liftoff(::Swan)::Any
└ @ BinaryTraits ~/.julia/dev/BinaryTraits/src/interface.jl:77
❌ Swan is missing these implementations:
1. FlyTrait: CanFly ⇢ fly(::<Type>, ::Float64, ::Float64)::Any
2. FlyTrait: CanFly ⇢ speed(::<Type>)::Float64
3. FlyTrait: CanFly ⇢ liftoff(::<Type>)::Any
=#


text:
using LatexTables

# Example 1
body        = rand(4,5)
body[:,3]   = 0.3*rand(4)
body[:,end] = 0.7.+0.3*rand(4)

s1 = table_to_tex(body;
    highlight_max_row=true,
    highlight_min_row=true,
    caption="And the winner is ... swan. Sorry penguin",
    leading_col=["We", "all", "like", "swans"],
    header=["Header", "Header2", "Penguin", "Header4", "Swan"],
)

print(s1)
make_tex_file("test1.tex", s1; dir_name="output", compile=true)

# Example 2
n    = 5
x1   = rand(n)
x2   = 1 .+ rand(n)
body = hcat(collect(1:n), x1, x2, x1./x2, [x./y >= 0.5 ? "Yes" : "No" for (x,y) in zip(x1,x2)])

s2 = table_to_tex(body;
    col_format=["d", "", "", "1p", ""],
    header=["Employee ID", "Result", "Target", "Performance", "Satisfaction"],
    alignment="lrrrr",
)

print(s2)
make_tex_file("test2.tex", s2; dir_name="output", compile=true)


text:
function soil(ebm::EBM)

    a = ebm.soil_a
    b = ebm.soil_b
    c = ebm.soil_c
    Gs = ebm.soil_Gs
    rhs = ebm.soil_rhs
    dTs = ebm.soil_dTs

    for k = 1:(ebm.Nsoil - 1)
        Gs[k] = 2.0 / (ebm.Dzsoil[k] / ebm.ksoil[k] + ebm.Dzsoil[k + 1] / ebm.ksoil[k + 1])
    end
    a[1] = 0.0
    b[1] = ebm.csoil[1] + Gs[1] * ebm.dt
    c[1] = - Gs[1] * ebm.dt
    rhs[1] = (ebm.Gsoil - Gs[1] * (ebm.Tsoil[1] - ebm.Tsoil[2])) * ebm.dt
    for k = 2:(ebm.Nsoil - 1)
        a[k] = c[k - 1]
        b[k] = ebm.csoil[k] + (Gs[k - 1] + Gs[k]) * ebm.dt
        c[k] = - Gs[k] * ebm.dt
        rhs[k] = Gs[k - 1] * (ebm.Tsoil[k - 1] - ebm.Tsoil[k]) * ebm.dt + Gs[k] * (ebm.Tsoil[k + 1] - ebm.Tsoil[k]) * ebm.dt
    end
    k = ebm.Nsoil
    Gs[k] = ebm.ksoil[k] / ebm.Dzsoil[k]
    a[k] = c[k - 1]
    b[k] = ebm.csoil[k] + (Gs[k - 1] + Gs[k]) * ebm.dt
    c[k] = 0.0
    rhs[k] = Gs[k - 1] * (ebm.Tsoil[k - 1] - ebm.Tsoil[k]) * ebm.dt
    tridiag(ebm.Nsoil, ebm.soil_gamma, a, b, c, rhs, dTs)
    for k = 1:ebm.Nsoil
        ebm.Tsoil[k] = ebm.Tsoil[k] + dTs[k]
    end

    return  nothing

end


text:
# test for the Utils module

"""
  Test calcBCNormal reverse mode.  This function works for any mesh dimension
"""
function test_calcBCNormal(mesh, eqn)
  dim = mesh.dim
  dxidx = rand(dim, dim)
  nrm = rand(dim)

  dxidx_c = zeros(Complex128, size(dxidx))
  copy!(dxidx_c, dxidx)
  nrm_out = zeros(Complex128, size(dxidx, 1))
  h = 1e-20
  pert = Complex128(0, h)

  jac = zeros(dim, dim*dim)
  for i=1:length(dxidx)
    dxidx_c[i] += pert
    calcBCNormal(eqn.params, dxidx_c, nrm, nrm_out)
    for j=1:length(nrm_out)
      jac[j, i] = imag(nrm_out[j])/h
    end
    dxidx_c[i] -= pert
  end

  nrm2_bar = zeros(dim)
  dxidx_bar = zeros(dim, dim)
  jac2 = zeros(jac)
  # compute reverse mode
  for i=1:length(nrm2_bar)
    nrm2_bar[i] = 1
    calcBCNormal_revm(eqn.params, dxidx, nrm, nrm2_bar, dxidx_bar)
    for j=1:length(dxidx_bar)
      jac2[i, j] = dxidx_bar[j]
    end

    fill!(dxidx_bar, 0.0)
    nrm2_bar[i] = 0
  end

  @test isapprox( norm(jac - jac2), 0.0) atol=1e-13
end

function test_area(mesh, sbp, eqn, opts)

  #----------------------------------------------------------------------------
  # test computeNormal
#  println("testing computeNormal")
  nout = mesh.dim*mesh.numNodesPerFace*mesh.numBoundaryFaces
  nin = mesh.dim*mesh.dim*mesh.numNodesPerFace*mesh.numBoundaryFaces
  jac = zeros(nout, nin)

  dxidx = rand(mesh.dim, mesh.dim, mesh.numNodesPerFace, mesh.numBoundaryFaces)

  nrm0 = Utils.computeNormal(mesh, eqn, mesh.bndryfaces, dxidx)
  pert = 1e-6
  for i=1:nin
    dxidx[i] += pert
    nrm = Utils.computeNormal(mesh, eqn, mesh.bndryfaces, dxidx)

    for j=1:nout
      jac[j, i] = (nrm[j] - nrm0[j])/pert
    end

    dxidx[i] -= pert
  end

  # reverse mode
  jac2 = zeros(jac)
  nrm_bar = zeros(nrm0)
  dxidx_bar = zeros(dxidx)
  for i=1:nout
    nrm_bar[i] = 1
    fill!(dxidx_bar, 0.0)
    Utils.computeNormal_rev(mesh, eqn, mesh.bndryfaces, nrm_bar, dxidx, dxidx_bar)

    for j=1:nin
      jac2[i, j] = dxidx_bar[j]
    end

    nrm_bar[i] = 0
  end

  for i=1:nin
    @test isapprox( norm(jac[:, i] - jac2[:, i])/size(jac, 1), 0.0) atol=1e-5
  end

  #----------------------------------------------------------------------------
  # test computeVolumeContribution
#  println("testing computeVolumecontribution")
  vol = calcVolumeContribution!(mesh, eqn, [1])
  @test isapprox( vol, 8) atol=1e-12

  # test computeVolumeContrib
  nout = 1
  nin = mesh.dim*mesh.numNodesPerFace*mesh.numBoundaryFaces
  jac = zeros(nout, nin)

  v0 = calcVolumeContribution!(mesh, eqn, [1])
  pert = 1e-6
  for i=1:nin
    mesh.nrm_bndry[i] += pert
    v_i = calcVolumeContribution!(mesh, eqn, [1])

    jac[1, i] = (v_i - v0)/pert

    mesh.nrm_bndry[i] -= pert
  end

  # reverse mode
  fill!(mesh.nrm_bndry_bar, 0.0)
  Utils.calcVolumeContribution_rev!(mesh, eqn, [1], 1.0)

  @test isapprox( norm(jac - reshape(mesh.nrm_bndry_bar, 1, nin))/size(jac, 2), 0.0) atol=1e-5

  #----------------------------------------------------------------------------
  # test calcProjectedAreaContribution
#  println("testing calcProjectedAreaContribution")

  proj_area = calcProjectedAreaContribution!(mesh, eqn, [1], 1)
  @test isapprox( proj_area, 4) atol=1e-12

  nout = 1
  nin = mesh.dim*mesh.numNodesPerFace*mesh.numBoundaryFaces
  jac = zeros(nout, nin)

  v0 = calcProjectedAreaContribution!(mesh, eqn, [1], 1)
  pert = 1e-6
  for i=1:nin
    mesh.nrm_bndry[i] += pert
    v_i = calcProjectedAreaContribution!(mesh, eqn, [1], 1)

    jac[1, i] = (v_i - v0)/pert

    mesh.nrm_bndry[i] -= pert
  end

  # reverse mode
  fill!(mesh.nrm_bndry_bar, 0.0)
  Utils.calcProjectedAreaContribution_rev!(mesh, eqn, [1], 1, 1.0)

  diffnorm = 0.0
  cnt = 0
  idx = 1
  for i=1:mesh.numBoundaryFaces
    for j=1:mesh.numNodesPerFace
      for k=1:mesh.dim
        if abs(mesh.nrm_bndry[k, j, i]) > 1e-8  # abs is not differentiable at 0
          diffnorm += abs(jac[1, idx] - mesh.nrm_bndry_bar[k, j, i])
          if abs(jac[1, idx] - mesh.nrm_bndry_bar[k, j, i]) > 1e-5
          end

          cnt += 1
        end
        idx += 1
      end
    end
  end

  diffnorm = sqrt(diffnorm)/cnt
  @test isapprox( diffnorm, 0.0) atol=1e-5

  # check the other method
#  println("testing other method")
#  nrm0 = Utils.computeNormal(mesh, eqn, mesh.bndryfaces,)
  nrm0 = copy(mesh.nrm_bndry)
  coords = mesh.coords_bndry

  nrm = copy(nrm0)
  jac = zeros(1, length(nrm0))
  for i=1:length(nrm)
    nrm[i] += pert
    v_i = calcProjectedAreaContribution!(mesh.sbpface, nrm, 1)
    jac[1, i] = (v_i - v0)/pert
    nrm[i] -= pert
  end

  nrm_bar = zeros(nrm)
  calcProjectedAreaContribution_rev!(mesh.sbpface, nrm, nrm_bar, 1, 1.0)

  diffnorm = norm(vec(jac) - vec(nrm_bar))/length(jac)
  diffnorm = 0.0
  cnt = 0
  for i=1:length(jac)
    diff = jac[i] - nrm_bar[i]

    if abs(nrm[i]) > 1e-8
      cnt += 1
      diffnorm += diff*diff
    end
  end

  diffnorm = sqrt(diffnorm)/length(jac)
  @test isapprox( diffnorm, 0.0) atol=1e-5

  return nothing
end

function test_norms(mesh, sbp, eqn, opts)

  u = ones(mesh.numDof)
  v = ones(mesh.numDof)

  # mesh is [0, 4] x [0, 4]
  @test isapprox( calcNorm(eqn, u), sqrt(mesh.numDofPerNode*16)) atol=1e-13

  # test linear
  u_arr = zeros(mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)
  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      x = mesh.coords[1, j, i]
      y = mesh.coords[2, j, i]
      dof = mesh.dofs[1, j, i]
      u[dof] = x + 2*y + 1
    end
  end

#  @test isapprox( calcNorm(eqn, u), 112) atol=1e-13

  # check inner product
  @test isapprox( calcL2InnerProduct(eqn, u, v), 112) atol=1e-13
  @test isapprox( calcL2InnerProduct(eqn, v, u), 112) atol=1e-13

  return nothing
end

function test_utils2(mesh, sbp, eqn, opts)

  @testset "----- Testing Utils 2D -----" begin

    test_calcBCNormal(mesh, eqn)
    test_norms(mesh, sbp, eqn, opts)
    test_identityarray()
    test_sparsematrix()
    test_binning()

    #TODO: test other things too
  end

end

function test_utils3(mesh, sbp, eqn, opts)

  @testset "----- Testing Utils 3D -----" begin

    test_calcBCNormal(mesh, eqn)
    test_area(mesh, sbp, eqn, opts)
  end

end

function test_identityarray()

  @testset "IdentityArray" begin
    @test size(IdentityArray{Int}(1)) == (1,)
    @test size(IdentityArray{Int}(1, 2)) == (1, 2)
    @test size(IdentityArray{Int}(1, 2, 3)) == (1, 2, 3)
    @test size(IdentityArray{Int}(1, 2, 3, 4)) == (1, 2, 3, 4)
    @test size(IdentityArray{Int}(1, 2, 3, 4, 5)) == (1, 2, 3, 4, 5)
    @test size(IdentityArray{Int}(1, 2, 3, 4, 5, 6)) == (1, 2, 3, 4, 5, 6)
    @test size(IdentityArray{Int}(1, 2, 3, 4, 5, 6, 7)) == (1, 2, 3, 4, 5, 6, 7)

    obj = IdentityArray{Int}(5)
    for i=1:5
      @test obj[i] == i
    end

    @test_throws ErrorException (obj[1] = 1)

    obj = IdentityArray{Int}(3, 4)
    for i=1:size(obj, 2)
      for j=1:size(obj, 1)
        @test obj[j, i] == j
      end
    end

    # test linear indexing
    @test obj[4] == 1
    @test_throws ErrorException (obj[1, 2] = 1)

    obj = IdentityArray{Int32}(3, 4)
    for i=1:size(obj, 2)
      for j=1:size(obj, 1)
        @test obj[j, i] == j
        @test typeof(obj[j, i]) == Int32  # make sure the return type matches eltype
      end
    end

    @test isbits(typeof(obj))  # test that this can be stack-allocated

    # also test FullFace
    @test isbits(FullFace{Float64})

    # test sview
    obj1 = sview(obj, :, 2)
    @inferred sview(obj, :, 2)
    @test size(obj1) == (3,)
    for i=1:3
      @test obj1[i] == i
    end

    obj2 = sview(obj, 2:3, 2)
    @inferred sview(obj, 2:3, 2)
    @test size(obj2) == (2,)
    for i=1:2
      @test obj2[i] == i + 1
    end

  end  # testset

end

function test_sparsematrix()

  @testset "SparseMatrixCSC" begin
    # test SparseMatrixCSC
    mat_dense = [1 3 0 0; 2 4 6 0; 0 5 7 0; 0 0 0 8]
    sparse_bnds = [1 1 2 4; 2 3 3 4]

    mat = SparseMatrixCSC(sparse_bnds, Float64)
    @time mat = SparseMatrixCSC(sparse_bnds, Float64)
    mat2 = sparse(mat_dense)

    @test ( mat.colptr )== mat2.colptr
    @test ( mat.rowval )== mat2.rowval
    @test ( mat2[1,1] )== 1
    @test ( mat2[2,1] )== 2
    @test ( mat2[1,2] )== 3
    @test ( mat2[2,2] )== 4
    @test ( mat2[3,2] )== 5
    @test ( mat2[2,3] )== 6
    @test ( mat2[3,3] )== 7
    @test ( mat2[4,4] )== 8

    mat2[3,2] = 9
    @test ( mat2[3,2] )== 9

    mat2[4,4] = 10
    @test ( mat2[4,4] )== 10

    # functions for checking sparsity
    cnt, out_of_bounds = checkSparseColumns(mat_dense, sparse_bnds, 1e-14)
    @test ( cnt )== 0

    cnt, out_of_bounds = checkSparseRows(mat_dense, sparse_bnds, 1e-14)
    @test ( cnt )== 0

    # modify sparse_bnds to verify the functions detect out of bounds entries
    sparse_bnds2 = copy(sparse_bnds)
    sparse_bnds2[2, 1] = 1
    cnt, out_of_bounds = checkSparseColumns(mat_dense, sparse_bnds2, 1e-14)
    @test ( cnt )== 1
    @test ( out_of_bounds )== [true, false, false, false]

    cnt, out_of_bounds = checkSparseRows(mat_dense, sparse_bnds2, 1e-14)
    @test ( cnt )== 1
    @test ( out_of_bounds )== [true, false, false, false]

    cnt = findLarge(mat_dense, 7.5)
    @test ( cnt )== 1

    # test fastfind
    A = rand(Int, 999)
    sort!(A)
    for i=1:length(A)
      idx = Jacobian.fastfind(A, A[i])
      @test ( idx )== i
    end

    A = rand(Int, 1000)
    sort!(A)
    for i=1:length(A)
      idx = Jacobian.fastfind(A, A[i])
      @test ( idx )== i
    end

    A = collect(1:10)
    idx = Jacobian.fastfind(A, 11)
    @test ( idx )== 0

    A = collect(1:2)
    idx = Jacobian.fastfind(A, 2)
    @test ( idx )== 2

  end

  return nothing
end

function test_binning()

  @testset "UniformBinning" begin
    obj = UniformBinning{Int, Float64}(1, 2, 3)

    @test calcBin(obj, 1.1) == 1
    @test calcBin(obj, 1.4) == 2
    @test calcBin(obj, 1.9) == 3
    @test calcBin(obj, 1) == 1
    @test calcBin(obj, 2) == 3

    push!(obj, 1, 1.1)
    push!(obj, 1, 1.2)

    push!(obj, 1, 1.4)
    push!(obj, 1, 1.4)
    push!(obj, 1, 1.5)

    push!(obj, 1, 1.7)
    push!(obj, 1, 1.8)
    push!(obj, 1, 1.9)
    push!(obj, 1, 2.0)

    bin_keycount, val_sum = sumBins(obj)

    @test bin_keycount[1] == 2
    @test bin_keycount[2] == 3
    @test bin_keycount[3] == 4

    @test abs(val_sum[1] - (1.1 + 1.2)) < 1e-15
    @test abs(val_sum[2] - (1.4 + 1.4 + 1.5)) < 1e-15
    @test abs(val_sum[3] - (1.7 + 1.8 + 1.9 + 2.0)) < 1e-15
  end

  return nothing
end

add_func2!(AdvectionTests, test_utils3, test_3d_inputfile, [TAG_REVERSEMODE, TAG_SHORTTEST])
# it doesn't matter what mesh is loaded, so reuse the test_lowlevel one
add_func2!(AdvectionTests, test_utils2, test_lowlevel_inputfile, [TAG_REVERSEMODE, TAG_SHORTTEST, TAG_TMP])


text:
using Images, FileIO, Printf

absdiff(a::RGB{T}, b::RGB{T}) where T = sum(abs(col(a) - col(b)) for col in (red, green, blue))
function pctdiff(A::Matrix{Color{T}}, B::Matrix{Color{T}}) where T
    size(A) != size(B) && throw(ArgumentError("images must be same-size"))
    s = zero(T)
    for (a, b) in zip(A, B)
        s += absdiff(a, b)
    end
    return 100s / 3prod(size(A))
end

img50  = load("data/lenna50.jpg") |> Matrix{RGB{Float64}};
img100 = load("data/lenna100.jpg") |> Matrix{RGB{Float64}};

d = pctdiff(img50, img100)
@printf("Percentage difference: %.4f%%\n", d)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using DataFrames
using CSV
using StringDistances

df = CSV.read("/Users/euclid/projects/ecolumix/epa/src/analysis/parent_co_level_data.csv")

parents_arr = df[:, 2]

parents_arr[ismissing.(parents_arr)] .= ""

findall("3M", parents_arr, Levenshtein(); min_score = 0.8)

function candidates(s, itr, sdist, min_score)
    idxs = findall(s, itr, sdist, min_score = min_score)
    itr[idxs]
end

candidates("3M", parents_arr, Winkler(Jaro()), 0.99)           # upper bound
candidates("3M", parents_arr, Winkler(Jaro()), 0.75)           # lower bound

candidates("A O SMITH", parents_arr, Winkler(Jaro()), 0.86)    # upper bound
candidates("A O SMITH", parents_arr, Winkler(Jaro()), 0.77)    # lower bound

candidates("A O SMITH", parents_arr, Jaccard(2), 0.99)         # upper bound
candidates("A O SMITH", parents_arr, Jaccard(2), 0.31)         # lower bound

candidates("A O SMITH", parents_arr, Overlap(2), 0.99)         # upper bound
candidates("A O SMITH", parents_arr, Overlap(2), 0.63)         # lower bound

candidates("EXXON MOBIL", parents_arr, Winkler(Jaro()), 0.98)  # upper bound
candidates("EXXON MOBIL", parents_arr, Winkler(Jaro()), 0.75)  # lower bound

candidates("SAINT-GOBAIN", parents_arr, Winkler(Jaro()), 0.87)    # upper bound
candidates("SAINT-GOBAIN", parents_arr, Winkler(Jaro()), 0.77)    # lower bound

candidates("SAINT-GOBAIN", parents_arr, Overlap(2), 0.77)         # upper bound
candidates("SAINT-GOBAIN", parents_arr, Overlap(2), 0.51)         # lower bound

###
# NOTE: Oldcastle is owned by CRH Americas, both own other
# so-called parent companies in TRI (e.g., Pike Industries).
###
candidates("OLDCASTLE", parents_arr, Overlap(2), 0.87)    # upper bound
candidates("OLDCASTLE", parents_arr, Overlap(2), 0.63)    # lower bound


text:
@testset "Sinc interpolation based on FFT" begin

    function test_interpolation_sum_fft(N_low, N)
	    x_min = 0.0
	    x_max = 16π

	    xs_low = range(x_min, x_max, length=N_low+1)[1:N_low]
	    xs_high = range(x_min, x_max, length=N)[1:end-1]
	    f(x) = sin(0.5*x) + cos(x) + cos(2 * x) + sin(0.25*x)
	    arr_low = f.(xs_low)
	    arr_high = f.(xs_high)

	    xs_interp = range(x_min, x_max, length=N+1)[1:N]
	    arr_interp = resample(arr_low, N)

	    xs_interp_s = range(x_min, x_max, length=N+1)[1:N]
	    arr_interp_s = FFTResampling.sinc_interpolate_sum(arr_low, N)

        @test ≈(arr_interp[2*N ÷10: N*8÷10], arr_high[2* N ÷10: N*8÷10], rtol=0.05)
        @test ≈(arr_high[2*N ÷10: N*8÷10], arr_interp_s[2*N ÷10: N*8÷10], rtol=0.05)
    end

    test_interpolation_sum_fft(128, 1000)
    test_interpolation_sum_fft(129, 1000)
    test_interpolation_sum_fft(120, 1531)
    test_interpolation_sum_fft(121, 1211)
end

@testset "Downsampling based on frequency cutting" begin
    function test_resample(N_low, N)
	    x_min = 0.0
	    x_max = 16π

	    xs_low = range(x_min, x_max, length=N_low+1)[1:N_low]
	    f(x) = sin(0.5*x) + cos(x) + cos(2 * x) + sin(0.25*x)
	    arr_low = f.(xs_low)

	    xs_interp = range(x_min, x_max, length=N+1)[1:N]
	    arr_interp = resample(arr_low, N)

	    xs_interp_s = range(x_min, x_max, length=N+1)[1:N]
	    arr_interp_s = FFTResampling.sinc_interpolate_sum(arr_low, N)

        arr_ds = FFTResampling.resample(arr_interp, N_low)
        @test ≈(arr_ds, arr_low)
        @test eltype(arr_low) === eltype(arr_ds)
        @test eltype(arr_interp) === eltype(arr_ds)
    end

    test_resample(128, 1000)
    test_resample(128, 1232)
    test_resample(128, 255)
    test_resample(253, 254)
    test_resample(253, 1001)
    test_resample(99, 100101)

end

@testset "Check if resample and interpolate return real result real input" begin
    function test_real(s, s_n, s_n2)
        x = randn(s)
        y = resample(x, s_n, take_real=false)
        y2 = resample(x, s_n2, take_real=false)
        @test all(imag.(y) .< 5e-15)
        @test all(imag.(y2) .< 5e-14)
    end

    test_real((12, 13, 14, 15), (12, 12, 12, 12), (20, 21, 22, 23))
    test_real((12, 13, 14, 15), (12, 13, 13, 13), (20, 22, 22, 23))
    test_real((12, 13, 14, 15), (12, 12, 13, 14), (21, 23, 21, 23))

    test_real((7,7), (6,6), (8, 8))
    test_real(100, 99, 101)
    test_real(101, 89, 101)
    test_real(100, 98, 120)
    test_real(101, 90, 120)

end

@testset "FFT sinc resample and upsample together in 2D" begin

    function test_2D(in_s, out_s)
        x = range(-10.0, 10.0, length=in_s[1] + 1)[1:end-1]
        y = range(-10.0, 10.0, length=in_s[2] + 1)[1:end-1]'
	    arr = abs.(x) .+ abs.(y) .+ sinc.(sqrt.(x .^2 .+ y .^2))
	    arr_interp = resample(arr[1:end, 1:end], out_s);
	    arr_ds = resample(arr_interp, in_s)
        @test arr_ds ≈ arr
    end

    test_2D((128, 128), (150, 150))
    test_2D((128, 128), (151, 151))
    test_2D((129, 129), (150, 150))
    test_2D((129, 129), (151, 151))

    test_2D((150, 128), (151, 150))
    test_2D((128, 128), (151, 153))
    test_2D((129, 128), (150, 153))
    test_2D((129, 128), (129, 153))

    x = range(-10.0, 10.0, length=129)[1:end-1]
    x2 = range(-10.0, 10.0, length=130)[1:end-1]
    x_exact = range(-10.0, 10.0, length=2049)[1:end-1]
    y = x'
    y2 = x2'
    y_exact = x_exact'
    arr = abs.(x) .+ abs.(y) .+sinc.(sqrt.(x .^2 .+ y .^2))
    arr2 = abs.(x) .+ abs.(y) .+sinc.(sqrt.(x .^2 .+ y .^2))
    arr_exact = abs.(x_exact) .+ abs.(y_exact) .+ sinc.(sqrt.(x_exact .^2 .+ y_exact .^2))
    arr_interp = resample(arr[1:end, 1:end], (131, 131));
    arr_interp2 = resample(arr[1:end, 1:end], (512, 512));
    arr_interp3 = resample(arr[1:end, 1:end], (1024, 1024));
    arr_ds = resample(arr_interp, (128, 128))
    arr_ds2 = resample(arr_interp, (128, 128))
    arr_ds23 = resample(arr_interp2, (512, 512))
    arr_ds3 = resample(arr_interp, (128, 128))

    @test ≈(arr_ds3, arr)
    @test ≈(arr_ds2, arr)
    @test ≈(arr_ds, arr)
    @test ≈(arr_ds23, arr_interp2)

end

@testset "FFT sinc resample and upsample together in 2D for a complex signal" begin

    function test_2D(in_s, out_s)
    	x = range(-10.0, 10.0, length=in_s[1] + 1)[1:end-1]
    	y = range(-10.0, 10.0, length=in_s[2] + 1)[1:end-1]'
    	f(x, y) = 1im * (abs(x) + abs(y) + sinc(sqrt(x ^2 + y ^2)))
    	f2(x, y) =  abs(x) + abs(y) + sinc(sqrt((x - 5) ^2 + (y - 5)^2))

    	arr = f.(x, y) .+ f2.(x, y)
    	arr_interp = resample(arr[1:end, 1:end], out_s);
    	arr_ds = resample(arr_interp, in_s)

        @test eltype(arr) === eltype(arr_ds)
        @test eltype(arr_interp) === eltype(arr_ds)
        @test imag(arr) ≈ imag(arr_ds)
        @test real(arr) ≈ real(arr_ds)
    end

    test_2D((128, 128), (150, 150))
    test_2D((128, 128), (151, 151))
    test_2D((129, 129), (150, 150))
    test_2D((129, 129), (151, 151))

    test_2D((150, 128), (151, 150))
    test_2D((128, 128), (151, 153))
    test_2D((129, 128), (150, 153))
    test_2D((129, 128), (129, 153))
end

@testset "FFT sinc resample and upsample together in 2D for a purely imaginary signal" begin
    function test_2D(in_s, out_s)
    	x = range(-10.0, 10.0, length=in_s[1] + 1)[1:end-1]
    	y = range(-10.0, 10.0, length=in_s[2] + 1)[1:end-1]'
    	f(x, y) = 1im * (abs(x) + abs(y) + sinc(sqrt(x ^2 + y ^2)))

    	arr = f.(x, y)
    	arr_interp = resample(arr[1:end, 1:end], out_s);
    	arr_ds = resample(arr_interp, in_s)

        @test imag(arr) ≈ imag(arr_ds)
        @test all(real(arr_ds) .< 1e-14)
        @test all(real(arr_interp) .< 1e-14)
    end

    test_2D((128, 128), (150, 150))
    test_2D((128, 128), (151, 151))
    test_2D((129, 129), (150, 150))
    test_2D((129, 129), (151, 151))

    test_2D((150, 128), (151, 150))
    test_2D((128, 128), (151, 153))
    test_2D((129, 128), (150, 153))
    test_2D((129, 128), (129, 153))

end

return


text:
# Functions related to the counting bound. This is similar
# to the original version, but focusses on the case of finding
# k-cliques in graphs with 2k vertices.
# It also uses several approximations:
# - various approximations of binomial coefficients
# - approximates (1-a)^b
# The bound, again, seems to asymptote at a small positive number.
# Possibly deprecated; I may be moving to only computing
# log-transformed coefficients.

module CountingBound2

export countingBound, approxNumMaximalCliques1, writeCounts

"""
Approximation of binomial(n, k), when n >> k.
From Wikipedia.
"""
function approxNChooseK(n, k)
  (n/k - 0.5)^k * exp(k) / sqrt(2 * pi * k)
end

"""
Approximation of binomial(2n, n), for n large.
From Wikipedia. (Possibly not used).
"""
function approx2KChooseK(n)
	2 ^ (2*n) / sqrt(pi * n)
end

"""
Counting bound, based on Shannon's argument.
  m: number of edges
  w: number of 'wires' -- that is, log2(number of functions),
    which is the number of bits needed to specify a function
  Returns: average number of NAND gates (with unbounded fan-in)
    required to compute any of those functions.
    (This may not be an integer).
"""
function countingBound(m, w)
  m = BigFloat(m)
  w = BigFloat(w)
  b = m - 0.5
  # the "-1" here is because this is the average, not the max.
  sqrt(2*w + b*b) - b - 1
end

"""
	Bound for finding k-cliques in n-vertex graphs, based on
covering a v-vertex graph with k-hypercliques.
	k: the size of the cliques to find
	n: number of vertices on which we're finding cliques
	v: the size of the larger graph
	Returns: minimal number of NAND gates need to find k-cliques
in an n-vertex graph (based on covering v-vertex graphs
with hypercliques).
"""
function cliqueNandBound(k, n, v)
	# force use of higher-precision arithmetic
	k = BigInt(k)
	n = BigInt(n)
	v = BigInt(v)

	# number of gates needed (for the larger problem)
	numGates1 = countingBound(binomial(v, 2), approxNChooseK(v, k))

	# (log of) probability a k-edge is not covered by an n-clique
	logP = - approxNChooseK(v-k, n-k) / (2^(binomial(n, k)-1))

	# number of k-edges not covered
	numNotCovered = exp(logP) * approxNChooseK(v, k)

	# number of expected n-vertex k-regular hypercliques
	numCliques = approxNChooseK(v, n) / (2^binomial(n, k))

	# to get the bound, we consider the number of gates covered by
	# n-cliques, and then divide by how many of those are expected
	numGates = (numGates1 - numNotCovered) / numCliques

end

"""
  Writes counts for some values of k, r, and n
  k: the value of k
  maxN: the maximum value of n
  outputDir: directory in which to write output files
  Side effects: writes files A and b.
    (The choices of r are a bit arbitrary).
"""
function writeCounts(k, maxN, outputDir)
  # values of n to use, from k to maxN
  # (for now, spaced at powers of 2)
  nList = [BigInt(2^i) for i in 1:trunc(log(2,maxN))]
  print(nList)
  print("\n")
  print(typeof(nList[1]))
  print("\n")
  # ??? for some reason, this was giving a syntax error
  # nList = [n for n in nList if n>=k && n<=maxN]
  # therefore, using filter()...
  nList = filter(n -> n>=k && n<=maxN, nList)

  # FIXME should create output directory
  # first, write the bound b
  of = open(outputDir * "/b_k=" * string(k) * "_maxN=" * string(maxN) * ".csv", "w")
  write(of, "k,n,bound,logBound\n")
  for n in nList
    bound = countingBound(binomial(k, 2), binomial(n, k))
    write(of, string(k) * "," * string(n) * "," * string(bound) * ","
      * string(log(bound)) * "\n")
  end
  close(of)

  # then, write the coefficients
  of = open(outputDir * "/A_k=" * string(k) * "_maxN=" * string(maxN) * ".csv", "w")
  write(of, "k,r,n,A,logA\n")
  for n = nList
    # FIXME what should the bound on r be?
    for r = k:min(n, 2*k)
      print("k=" * string(k) * " r=" * string(r) * " n=" * string(n) * "\n")
      bound = countingBound(binomial(k, 2), binomial(n, k))
      A = approxNumMaximalCliques1(k, r, n)
      logA = logApproxNumMaximalCliques1(k, r, n)
      write(of, string(k) * "," * string(r) * "," * string(n)
        * "," * string(A)
        * "," * string(logA) * "\n")
    end
  end
  close(of)

end

end

max_stars_repo_name: joshtburdick/misc

max_issues_repo_name: joshtburdick/misc

max_forks_repo_name: joshtburdick/misc


text:
@testset "Monotones" begin
    @test_throws MethodError relative_entropy([], [])
    @test_throws ErrorException relative_entropy([1.0], [0.5, 0.5])
    @test_throws ErrorException relative_entropy([0.5, 0.5], [1.0])
    @test_throws ErrorException relative_entropy([0.5], [1.0])
    @test_throws ErrorException relative_entropy([1.0], [0.5])

    @test relative_entropy([1.0, 0.0], [0.5, 0.5]) ≈ log(2.0)
    @test relative_entropy([0.5, 0.5], [0.5, 0.5]) ≈ 0.0
    @test isinf(relative_entropy([0.5, 0.5], [1.0, 0.0]))
    @test relative_entropy([0.25, 0.75], [0.5, 0.5]) ≈ 0.75*log(3) - log(2)
end


text:
using HDF5,JLD,KUparser,KUnet,Base.Test
require("mycat.jl")

if !isdefined(:corpus)
    info("Loading data")
    @date @load "yatbazp.tst.jld4"
end
@show ncpu = 10
@show nbatch = 10
@show ndeps = length(deprel)
@show ft = Flist.acl11eager
@show pt = ArcEager13
p0 = pt(1,ndeps)
s1 = corpus[1]

if !isdefined(:net)
    @date net = loadnet("nets/zn11beam4.net4")
    # info("Let's train a reasonable model")
    # h0 = 20000
    # net = [Mmul(h0), Bias(), Relu(),
    #        Mmul(p0.nmove), Bias(), Logp(),
    #        LogpLoss()]
    # @date (p,x,y)=oparse(pt, corpus, ndeps, ft)
    # # @date train(net, x, y; loss=logploss)
    # @date train(net, x, y)
end

@show nbeam = 10
@date (q2,x2,y2)=bparse(pt, corpus[2:2], ndeps, ft, net, nbeam; xy=true)
@show map(size, (q2,x2,y2))
@show evalparse(q2, corpus[2:2])

if false

info("Use gparser as reference for beam=1")
@date gxy=(gp,gx,gy)=gparse(pt, corpus, ndeps, ft, net; xy=true)
@show map(size, gxy)
@show evalparse(gp, corpus)

@show nbeam = 1

info("nbeam=1 Single sentence")
@date p1=bparse(pt, [s1], ndeps, ft, net, nbeam)
@test @show isequal(p1[1],gp[1])
@date (p2,x2,y2)=bparse(pt, [s1], ndeps, ft, net, nbeam; xy=true)
I2 = 1:size(x2,2)
@test @show isequal(p2[1],gp[1])
@test @show isequal(x2,gx[:,I2])
@test @show isequal(y2,gy[:,I2])

info("nbeam=1 Multiple sentences")
@date p3=bparse(pt, corpus, ndeps, ft, net, nbeam)
@test @show isequal(p3,gp)
@date pxy4=bparse(pt, corpus, ndeps, ft, net, nbeam; xy=true)
@show map(size, pxy4)
@show map(size, gxy)
@test @show isequal(pxy4[1],gxy[1])

info("nbeam=1 Batch processing")
@date p5=bparse(pt, corpus, ndeps, ft, net, nbeam, nbatch)
@test @show isequal(p5,gp)
@date pxy6=bparse(pt, corpus, ndeps, ft, net, nbeam, nbatch; xy=true)
@test @show pxyequal(pxy6, pxy4)

info("nbeam=1 Multi-cpu processing")
@date p7=bparse(pt, corpus, ndeps, ft, net, nbeam, nbatch, ncpu)
@test @show isequal(mycat(p7),gp)
@date pxy8=bparse(pt, corpus, ndeps, ft, net, nbeam, nbatch, ncpu; xy=true)
@test @show pxyequal(pxycat(pxy8), pxy4)

@show nbeam = 10

info("nbeam=10 Single sentence")
@date q1=bparse(pt, [s1], ndeps, ft, net, nbeam)
@show evalparse(q1, [s1])
@date (q2,x2,y2)=bparse(pt, [s1], ndeps, ft, net, nbeam; xy=true)
@show map(size, (q2,x2,y2))
@test @show isequal(q2,q1)

info("nbeam=10 Multiple sentences")
@date q3=bparse(pt, corpus, ndeps, ft, net, nbeam)
@show evalparse(q3, corpus)
@test @show isequal(q3[1],q1[1])
@date qxy4=(q4,x4,y4)=bparse(pt, corpus, ndeps, ft, net, nbeam; xy=true)
@show map(size, qxy4)
@test @show isequal(q4, q3)
I2 = 1:size(x2,2)
@test @show isequal(x4[:,I2], x2)
@test @show isequal(y4[:,I2], y2)

info("nbeam=10 Batch processing")
@date q5=bparse(pt, corpus, ndeps, ft, net, nbeam, nbatch)
@test @show isequal(q5,q3)
@date qxy6=bparse(pt, corpus, ndeps, ft, net, nbeam, nbatch; xy=true)
@show map(size, qxy6)
@test @show pxyequal(qxy6, qxy4)

info("nbeam=10 Multi-cpu processing")
@date q7=bparse(pt, corpus, ndeps, ft, net, nbeam, nbatch, ncpu)
@test @show isequal(mycat(q7),qxy6[1])
@date qxy8=bparse(pt, corpus, ndeps, ft, net, nbeam, nbatch, ncpu; xy=true)
@test @show pxyequal(pxycat(qxy8), qxy6)

end


text:
import Oceananigans.TimeSteppers: calculate_tendencies!

using Oceananigans: fields
using Oceananigans.Utils: work_layout

"""
    calculate_tendencies!(model::IncompressibleModel)

Calculate the interior and boundary contributions to tendency terms without the
contribution from non-hydrostatic pressure.
"""
function calculate_tendencies!(model::HydrostaticFreeSurfaceModel)

    # Calculate contributions to momentum and tracer tendencies from fluxes and volume terms in the
    # interior of the domain
    calculate_hydrostatic_free_surface_interior_tendency_contributions!(model.timestepper.Gⁿ,
                                                                        model.architecture,
                                                                        model.grid,
                                                                        model.advection,
                                                                        model.coriolis,
                                                                        model.buoyancy,
                                                                        model.closure,
                                                                        model.velocities,
                                                                        model.free_surface,
                                                                        model.tracers,
                                                                        model.pressure.pHY′,
                                                                        model.diffusivities,
                                                                        model.forcing,
                                                                        model.clock)

    #=
    # Calculate contributions to momentum and tracer tendencies from user-prescribed fluxes across the
    # boundaries of the domain
    calculate_boundary_tendency_contributions!(model.timestepper.Gⁿ,
                                               model.architecture,
                                               model.velocities,
                                               model.tracers,
                                               model.clock,
                                               fields(model))
    =#

    return nothing
end

function calculate_hydrostatic_momentum_tendencies!(tendencies, velocities, arch, grid, advection, coriolis, closure,
                                                    free_surface, tracers, diffusivities, hydrostatic_pressure_anomaly,
                                                    forcings, clock, barrier)

    calculate_Gu_kernel! = calculate_hydrostatic_free_surface_Gu!(device(arch), work_layout(grid, :xyz)...)
    calculate_Gv_kernel! = calculate_hydrostatic_free_surface_Gv!(device(arch), work_layout(grid, :xyz)...)
    calculate_Gη_kernel! = calculate_hydrostatic_free_surface_Gη!(device(arch), work_layout(grid, :xy)...)

    Gu_event = calculate_Gu_kernel!(tendencies.u, grid, advection.momentum, coriolis, closure,
                                    velocities, free_surface, tracers, diffusivities, hydrostatic_pressure_anomaly,
                                    forcings, clock; dependencies = barrier)

    Gv_event = calculate_Gv_kernel!(tendencies.v, grid, advection.momentum, coriolis, closure,
                                    velocities, free_surface, tracers, diffusivities, hydrostatic_pressure_anomaly,
                                    forcings, clock; dependencies = barrier)

    Gη_event = calculate_Gη_kernel!(tendencies.η, grid, velocities, free_surface, tracers,
                                    forcings, clock; dependencies = barrier)

    events = [Gu_event, Gv_event, Gη_event]

    return events
end

""" Store previous value of the source term and calculate current source term. """
function calculate_hydrostatic_free_surface_interior_tendency_contributions!(tendencies,
                                                                             arch,
                                                                             grid,
                                                                             advection,
                                                                             coriolis,
                                                                             buoyancy,
                                                                             closure,
                                                                             velocities,
                                                                             free_surface,
                                                                             tracers,
                                                                             hydrostatic_pressure_anomaly,
                                                                             diffusivities,
                                                                             forcings,
                                                                             clock)

    calculate_Gc_kernel! = calculate_hydrostatic_free_surface_Gc!(device(arch), work_layout(grid, :xyz)...)

    barrier = Event(device(arch))

    events = calculate_hydrostatic_momentum_tendencies!(tendencies, velocities, arch, grid, advection, coriolis, closure,
                                                        free_surface, tracers, diffusivities, hydrostatic_pressure_anomaly,
                                                        forcings, clock, barrier)

    for (tracer_index, tracer_name) in enumerate(propertynames(tracers))
        @inbounds c_tendency = tendencies[tracer_name]
        @inbounds c_advection = advection[tracer_name]
        @inbounds forcing = forcings[tracer_name]

        Gc_event = calculate_Gc_kernel!(c_tendency, grid, Val(tracer_index),
                                        c_advection, closure, buoyancy,
                                        velocities, free_surface, tracers, diffusivities,
                                        forcing, clock; dependencies = barrier)

        push!(events, Gc_event)
    end

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end

#####
##### Tendency calculators for u, v
#####

""" Calculate the right-hand-side of the u-velocity equation. """
@kernel function calculate_hydrostatic_free_surface_Gu!(Gu, grid, args...)
    i, j, k = @index(Global, NTuple)
    @inbounds Gu[i, j, k] = hydrostatic_free_surface_u_velocity_tendency(i, j, k, grid, args...)
end

""" Calculate the right-hand-side of the v-velocity equation. """
@kernel function calculate_hydrostatic_free_surface_Gv!(Gv, grid, args...)
    i, j, k = @index(Global, NTuple)
    @inbounds Gv[i, j, k] = hydrostatic_free_surface_v_velocity_tendency(i, j, k, grid, args...)
end

#####
##### Tendency calculators for tracers
#####

""" Calculate the right-hand-side of the tracer advection-diffusion equation. """
@kernel function calculate_hydrostatic_free_surface_Gc!(Gc, grid, args...)
    i, j, k = @index(Global, NTuple)
    @inbounds Gc[i, j, k] = hydrostatic_free_surface_tracer_tendency(i, j, k, grid, args...)
end

#####
##### Tendency calculators for an explicit free surface
#####

""" Calculate the right-hand-side of the free surface displacement (η) equation. """
@kernel function calculate_hydrostatic_free_surface_Gη!(Gη, grid, args...)
    i, j = @index(Global, NTuple)
    @inbounds Gη[i, j, 1] = free_surface_tendency(i, j, grid, args...)
end


text:
using ChainRulesCore
using ChainRulesOverloadGeneration

using Test

@testset "ChainRulesCore" begin
    @testset "demos" begin
        include("demos/forwarddiffzero.jl")
        include("demos/reversediffzero.jl")
    end

    # Do this after demos run, so that the simple demo code doesn't have to handle
    # anything weird we define for testing purposes
    include("ruleset_loading.jl")
end


text:
using Statistics

if Sys.isunix()
    ENV["JULIA_EDITOR"] = "vim"
elseif Sys.iswindows()
    ENV["JULIA_EDITOR"] = "C:/Program\\ Files/vim/vim81/gvim.exe"
end

"Change the current terminal background from light to dark and vice versa."
function tbg()
    # The shell won't see these changes but it's still useful for us.
    if ENV["SOLARIZED_THEME"] == "dark"
        ENV["SOLARIZED_THEME"] = "light"
    else
        ENV["SOLARIZED_THEME"] = "dark"
    end

    gnome_term_profiles = "/org/gnome/terminal/legacy/profiles:/"
    first_profile = chomp(read(pipeline(`dconf list $gnome_term_profiles`, `head -n 1`), String))
    if ENV["SOLARIZED_THEME"] == "dark"
        bg_color = "'rgb(0,43,54)'"
        fg_color = "'rgb(131,148,150)'"
    else
        bg_color = "'rgb(253,246,227)'"
        fg_color = "'rgb(101,123,131)'"
    end
    run(`dconf write "$gnome_term_profiles$(first_profile)background-color" "$bg_color"`)
    run(`dconf write "$gnome_term_profiles$(first_profile)foreground-color" "$fg_color"`)

    #=
    if !isempty(ENV["ZSH_THEME"])
        run(`source $(ENV["ZSH"])/oh-my-zsh.sh`)
    end

    # FZF
    _gen_fzf_default_opts
    [ -f ~/.fzf.zsh ] && source ~/.fzf.zsh
    =#

    nothing
end

max_stars_repo_name: janEbert/dotfiles

max_issues_repo_name: janEbert/dotfiles

max_forks_repo_name: janEbert/dotfiles


text:
# Canonical form of multivariate normal

### Generic types

immutable MvNormalCanon{T<:Real,P<:AbstractPDMat,V<:Union{Vector,ZeroVector}} <: AbstractMvNormal
    μ::V    # the mean vector
    h::V    # potential vector, i.e. inv(Σ) * μ
    J::P    # precision matrix, i.e. inv(Σ)
end

typealias FullNormalCanon MvNormalCanon{Float64, PDMat{Float64,Matrix{Float64}},Vector{Float64}}
typealias DiagNormalCanon MvNormalCanon{Float64,PDiagMat{Float64,Vector{Float64}},Vector{Float64}}
typealias IsoNormalCanon  MvNormalCanon{Float64,ScalMat{Float64},Vector{Float64}}

typealias ZeroMeanFullNormalCanon MvNormalCanon{Float64,PDMat{Float64,Matrix{Float64}},ZeroVector{Float64}}
typealias ZeroMeanDiagNormalCanon MvNormalCanon{Float64,PDiagMat{Float64,Vector{Float64}},ZeroVector{Float64}}
typealias ZeroMeanIsoNormalCanon  MvNormalCanon{Float64,ScalMat{Float64},ZeroVector{Float64}}

### Constructors

function MvNormalCanon{T<:Real}(μ::Vector{T}, h::Vector{T}, J::AbstractPDMat{T})
    length(μ) == length(h) == dim(J) || throw(DimensionMismatch("Inconsistent argument dimensions"))
    MvNormalCanon{T,typeof(J),typeof(μ)}(μ, h, J)
end

function MvNormalCanon{T<:Real, P<:AbstractPDMat}(μ::Vector{T}, h::Vector{T}, J::P)
    R = promote_type(T, eltype(J))
    MvNormalCanon(Vector{R}(μ), promote_eltype(h, J)...)
end

MvNormalCanon{T<:Real, S<:Real, P<:AbstractPDMat}(μ::Vector{T}, h::Vector{S}, J::P) = MvNormalCanon(promote_eltype(μ, h)..., J)

function MvNormalCanon{P<:AbstractPDMat}(J::P)
    z = ZeroVector(eltype(J), dim(J))
    MvNormalCanon{eltype(J),P,ZeroVector{eltype(J)}}(z, z, J)
end

function MvNormalCanon{T<:Real, P<:AbstractPDMat}(h::Vector{T}, J::P)
    length(h) == dim(J) || throw(DimensionMismatch("Inconsistent argument dimensions"))
    hh, JJ = promote_eltype(h, J)
    MvNormalCanon{eltype(hh),typeof(JJ),typeof(hh)}(JJ \ hh, hh, JJ)
end

MvNormalCanon{T<:Real}(h::Vector{T}, J::Matrix{T}) = MvNormalCanon(h, PDMat(J))
MvNormalCanon{T<:Real}(h::Vector{T}, prec::Vector{T}) = MvNormalCanon(h, PDiagMat(prec))
MvNormalCanon{T<:Real}(h::Vector{T}, prec::T) = MvNormalCanon(h, ScalMat(length(h), prec))

MvNormalCanon{T<:Real, S<:Real}(h::Vector{T}, J::VecOrMat{S}) = MvNormalCanon(promote_eltype(h, J)...)
MvNormalCanon{T<:Real, S<:Real}(h::Vector{T}, prec::S) = MvNormalCanon(promote_eltype(h, prec)...)

MvNormalCanon(J::Matrix) = MvNormalCanon(PDMat(J))
MvNormalCanon(prec::Vector) = MvNormalCanon(PDiagMat(prec))
MvNormalCanon(d::Int, prec) = MvNormalCanon(ScalMat(d, prec))

### Show

distrname(d::IsoNormalCanon) = "IsoNormalCanon"
distrname(d::DiagNormalCanon) = "DiagNormalCanon"
distrname(d::FullNormalCanon) = "FullNormalCanon"

distrname(d::ZeroMeanIsoNormalCanon) = "ZeroMeanIsoNormalCanon"
distrname(d::ZeroMeanDiagNormalCanon) = "ZeroMeanDiagormalCanon"
distrname(d::ZeroMeanFullNormalCanon) = "ZeroMeanFullNormalCanon"

### Conversion
function convert{T<:Real}(::Type{MvNormalCanon{T}}, d::MvNormalCanon)
    MvNormalCanon(convert_eltype(T, d.μ), convert_eltype(T, d.h), convert_eltype(T, d.J))
end
function convert{T<:Real,V<:Union{Vector, ZeroVector}}(::Type{MvNormalCanon{T}}, μ::V, h::V, J::AbstractPDMat)
    MvNormalCanon(convert_eltype(T, μ), convert_eltype(T, h), convert_eltype(T, J))
end

### conversion between conventional form and canonical form

meanform(d::MvNormalCanon) = MvNormal(d.μ, inv(d.J))
# meanform{C, T<:Real}(d::MvNormalCanon{T,C,Vector{T}}) = MvNormal(d.μ, inv(d.J))
# meanform{C, T<:Real}(d::MvNormalCanon{T,C,ZeroVector{T}}) = MvNormal(inv(d.J))

canonform{C, T<:Real}(d::MvNormal{T,C,Vector{T}}) = (J = inv(d.Σ); MvNormalCanon(d.μ, J * d.μ, J))
canonform{C, T<:Real}(d::MvNormal{T,C,ZeroVector{T}}) = MvNormalCanon(inv(d.Σ))

### Basic statistics

length(d::MvNormalCanon) = length(d.μ)
mean(d::MvNormalCanon) = convert(Vector{eltype(d.μ)}, d.μ)
params(d::MvNormalCanon) = (d.μ, d.h, d.J)
@inline partype{T<:Real}(d::MvNormalCanon{T}) = T

var(d::MvNormalCanon) = diag(inv(d.J))
cov(d::MvNormalCanon) = full(inv(d.J))
invcov(d::MvNormalCanon) = full(d.J)
logdetcov(d::MvNormalCanon) = -logdet(d.J)

### Evaluation

sqmahal(d::MvNormalCanon, x::AbstractVector) = quad(d.J, x - d.μ)
sqmahal!(r::AbstractVector, d::MvNormalCanon, x::AbstractMatrix) = quad!(r, d.J, x .- d.μ)

# Sampling (for GenericMvNormal)

unwhiten_winv!(J::AbstractPDMat, x::AbstractVecOrMat) = unwhiten!(inv(J), x)
unwhiten_winv!(J::PDiagMat, x::AbstractVecOrMat) = whiten!(J, x)
unwhiten_winv!(J::ScalMat, x::AbstractVecOrMat) = whiten!(J, x)

_rand!(d::MvNormalCanon, x::AbstractMatrix) = add!(unwhiten_winv!(d.J, randn!(x)), d.μ)


text:
# fetch the field into the low order bits
# all other bits are zeroed

get_sign(x::T) where T = isolate_sign(x) >> sign_offset(T)
get_exponent(x::T) where T = isolate_exponent(x) >> exponent_offset(T)
get_significand(x::T) where T = isolate_significand(x) >> significand_offset(T)
get_sign_exponent(x::T) where T = isolate_sign_exponent(x) >> exponent_offset(T)
get_exponent_significand(x::T) where T = isolate_exponent_significand(x) >> significand_offset(T)

# prepared field is in the low order bits, shift field into its place
# all other bits are zeroed

set_sign(x::T) where T = (x & sign_low_mask(T)) << sign_offset(T)
set_exponent(x::T) where T = (x & exponent_low_mask(T)) << exponent_offset(T)
set_significand(x::T) where T = (x & significand_low_mask(T)) << significand_offset(T)
set_sign_exponent(x::T) where T = (x & sign_exponent_low_mask(T)) << exponent_offset(T)
set_exponent_significand(x::T) where T = (x & exponent_significand_low_mask(T)) << exponent_significand_offset(T)

# isolate the field in place
# (yields the field value, other bits are zeroed)

isolate_sign(x::T) where T = unsigned(x) & sign_mask(T)
isolate_exponent(x::T) where T = unsigned(x) & exponent_mask(T)
isolate_significand(x::T) where T = unsigned(x) & significand_mask(T)
isolate_sign_exponent(x::T) where T = unsigned(x) & sign_exponent_mask(T)
isolate_exponent_significand(x::T) where T = unsigned(x) & sign_filter(T)

# clear the field in place
# (yields the other fields, this field is zeroed)

clear_sign(x::T) where T = unsigned(x) & sign_filter(T)
clear_exponent(x::T) where T = unsigned(x) & exponent_filter(T)
clear_significand(x::T) where T = unsigned(x) & significand_filter(T)
clear_sign_exponent(x::T) where T = unsigned(x) & sign_exponent_filter(T)
clear_exponent_significand(x::T) where T = unsigned(x) & exponent_significand_filter(T)

# field offsets (shift by)

sign_offset(::Type{T}) where T = bitwidth(T) - one(signed(T))
exponent_offset(::Type{T}) where T = sign_offset(T) - exponent_bits(T)
significand_offset(::Type{T}) where T = zero(signed(T))
sign_exponent_offset(::Type{T}) where T = exponent_offset(T)
exponent_significand_offset(::Type{T}) where T = significand_offset(T)

# field filters

sign_filter(::Type{T}) where T = ~zero(unsigned(T)) >>> 1
exponent_filter(::Type{T}) where T = sign_exponent_filter(T) | sign_mask(T)
significand_filter(::Type{T}) where T = ~sign_exponent_filter(T)
sign_exponent_filter(::Type{T}) where T = ~zero(unsigned(T)) >>> (exponent_bits(T) + 1)
exponent_significand_filter(::Type{T}) where T = ~sign_filter(T)

# field masks
# exponent_mask, significand_mask exist in Base
sign_mask(::Type{T}) where T = ~sign_filter(T)
sign_exponent_mask(::Type{T}) where T = ~sign_exponent_filter(T)
exponent_significand_mask(::Type{T}) where T = ~exponent_significand_mask(T)

# place plain masks in the least significant bits

sign_low_mask(::Type{T}) where T = sign_mask(T) >> sign_offset(T)
exponent_low_mask(::Type{T}) where T = exponent_mask(T) >> exponent_offset(T)
significand_low_mask(::Type{T}) where T = significand_mask(T) >> significand_offset(T)
sign_exponent_low_mask(::Type{T}) where T = sign_exponent_mask(T) >> exponent_offset(T)
exponent_significand_low_mask(::Type{T}) where T = exponent_significand_mask(T) >> significand_offset(T)


text:
type_pretty(x::UnionAll) = type_pretty(x.body)
type_pretty(x::DataType) = x.name


text:
@testset "Grid diameter identify" begin
    l = 12
    p = 66

    img = zeros(1, 500)
    for i = 20:length(img)
        if (mod(i, p) < l)
            img[i] = 1
        end
    end

    img1 = img + rand(1:500, 1, 500) / 2500

    img2 = zeros(1, 500)
    for i = 20:length(img)
        img2[i] = (img[i-4] + img[i-3] + img[i-2] + img[i-1] + img[i]) / 5
    end

    img3 = img2 + rand(1:500, 1, 500) / 2500

    # sharp grids
    stript, space = grid_resolution(img)
    @test stript ≈ l atol = 1
    @test space + stript ≈ p atol = 1

    # sharp grids with random
    stript, space = grid_resolution(img1)
    @test stript ≈ l atol = 1
    @test space + stript ≈ p atol = 1

    # moving average grids window 5
    stript, space = grid_resolution(img2)
    @test stript ≈ l atol = 1
    @test space + stript ≈ p atol = 1

    # moving average grids window 5 with random
    stript, space = grid_resolution(img3)
    @test stript ≈ l atol = 1
    @test space + stript ≈ p atol = 1
end


text:
"""
    Bounds.MultiEllipsoid([T=Float64], ndims)
    Bounds.MultiEllipsoid(::AbstractVector{Ellipsoid})

Use multiple [`Ellipsoid`](@ref)s in an optimal clustering to bound prior space. For more details about the bounding algorithm, see the extended help (`??Bounds.MultiEllipsoid`)
"""
struct MultiEllipsoid{T} <: AbstractBoundingSpace{T}
    ellipsoids::Vector{Ellipsoid{T}}
end
MultiEllipsoid(ndims::Integer) = MultiEllipsoid(Float64, ndims)
MultiEllipsoid(T::Type, ndims::Integer) = MultiEllipsoid([Ellipsoid(T, ndims)])

Base.eltype(::MultiEllipsoid{T}) where {T} = T
Base.broadcastable(me::MultiEllipsoid) = (me,)
Base.length(me::MultiEllipsoid) = length(me.ellipsoids)
Base.ndims(me::MultiEllipsoid) = ndims(me.ellipsoids[1])

volume(me::MultiEllipsoid) = sum(volume, me.ellipsoids)

function scale!(me::MultiEllipsoid, factor)
    scale!.(me.ellipsoids, factor)
    return me
end

function fit(::Type{<:MultiEllipsoid}, x::AbstractMatrix; pointvol = 0)
    parent = fit(Ellipsoid, x, pointvol = pointvol)
    ells = fit(MultiEllipsoid, x, parent, pointvol = pointvol)
    return MultiEllipsoid(ells)
end

function fit(::Type{<:MultiEllipsoid}, x::AbstractMatrix, parent::Ellipsoid; pointvol = 0)
    ndim, npoints = size(x)

    # Clustering will fail with fewer than k=2 points
    npoints ≤ 2 && return [parent]

    p1, p2 = endpoints(parent)
    starting_points = hcat(p1, p2)
    R = kmeans!(x, starting_points; maxiter = 10)
    labels = assignments(R)
    x1 = x[:, labels .== 1]
    x2 = x[:, labels .== 2]

    # if either cluster has fewer than ndim points, it is ill-defined
    if size(x1, 2) < 2ndim || size(x2, 2) < 2ndim
        return [parent]
    end

    # Getting bounding ellipsoid for each cluster
    ell1, ell2 = fit.(Ellipsoid, (x1, x2), pointvol = pointvol)

    # If total volume decreased by over half, recurse
    if volume(ell1) + volume(ell2) < 0.5volume(parent)
        return vcat(fit(MultiEllipsoid, x1, ell1, pointvol = pointvol),
                    fit(MultiEllipsoid, x2, ell2, pointvol = pointvol))
    end

    # Otherwise see if total volume is much larger than expected
    # and split into more than 2 clusters
    if volume(parent) > 2npoints * pointvol
        out = vcat(fit(MultiEllipsoid, x1, ell1, pointvol = pointvol),
                    fit(MultiEllipsoid, x2, ell2, pointvol = pointvol))
        sum(volume, out) < 0.5volume(parent) && return out
    end

    # Otherwise, return single bounding ellipse
    return [parent]
end

Base.in(x::AbstractVector, me::MultiEllipsoid) = any(ell->x ∈ ell, me.ellipsoids)

function Base.rand(rng::AbstractRNG, me::MultiEllipsoid)
    length(me) == 1 && return rand(rng, me.ellipsoids[1])

    vols = volume.(me.ellipsoids)
    weights = vols ./ sum(vols)
    local x

    while true
        # Select random ellipsoid
        idx = rand(rng, Categorical(weights))
        ell = me.ellipsoids[idx]

        # Select point
        x = rand(rng, ell)

        # How many ellipsoids is the sample in
        n = count(ell -> x ∈ ell, me.ellipsoids)

        # Only accept with probability 1/n
        (n == 1 || rand(rng) < 1 / n) && break
    end

    return x
end

"""
Returns a random live point and a bounding ellipsoid, since MultiEllipsoid doesn't have
valid transformation axes
"""
function rand_live(rng::AbstractRNG, me::MultiEllipsoid, us)
    idx = rand(rng, Base.axes(us, 2))
    u = us[:, idx]

    # find which Ellipsoid/s it overlaps with
    idxs = findall(ell -> u ∈ ell, me.ellipsoids)
    # TODO if point isn't bounded, update bounds
    if isempty(idxs)
        @warn "You have an unbounded point, somehow. Future PR will re-fit bounds. For now, using a workaround"
        return rand_live(rng, me, us)
    end

    # pick random encompassing ellipsoid
    idx = rand(rng, idxs)

    return u, me.ellipsoids[idx]
end


text:
using AdjacentFloats
using Test

@testset "AdjacentFloats" begin

    for (f, g) in ( (prev_float, prevfloat),
                    (next_float, nextfloat) )

        @testset "Tests for $f" begin

            for x in (1.0, 3.0, sqrt(2.0), sqrt(0.25),
                      1e-10, 1e10, 1e300, 1e-300, Inf)

                @test f(x) == g(x)
                @test f(-x) == g(-x)
            end
        end
    end

    @testset "NaN" begin
        @test isnan(prev_float(NaN))
        @test isnan(next_float(NaN))
    end
end