text:
################################################################################
# VisualizationOptions
################################################################################

@with_kw mutable struct RoadwayVisualizationOptions{T}
    # Options
	"Thickness of the object."
	ϵ::T=0.05

	"Line width."
	line_width::T=0.007

	"Bound width."
	bound_width::T=0.02

	"Bound height."
	bound_height::T=0.05

	"Road color."
	road_color::RGB=RGB(27/255, 50/255, 81/255)

	"Line color."
	line_color::RGB=RGB(255/255, 206/255, 68/255)

	"Road color."
	bound_color::RGB=RGB(75/255, 129/255, 194/255)
end

################################################################################
# set_roadway!
################################################################################

function set_roadway!(vis::Visualizer, opts::HighwayRoadwayOptions{T};
    vis_opts::RoadwayVisualizationOptions=RoadwayVisualizationOptions(), k) where {T}
    clean!(vis)
    ll = opts.lane_length
    lw = opts.lane_width

	ϵ = vis_opts.ϵ
	liw = vis_opts.line_width
	bw = vis_opts.bound_width
	bh = vis_opts.bound_height

    road = Rect3D(Vec(0.0, -lw, -ϵ), Vec(ll, 2lw, ϵ))
    line = Rect3D(Vec(0.0, -liw/2, 0.0), Vec(ll, liw, liw))
	left_bound = Rect3D(Vec(0.0, lw, -ϵ), Vec(ll, bw, ϵ+bh))
    right_bound = Rect3D(Vec(0.0, -lw-bw, -ϵ), Vec(ll, bw, ϵ+bh))

    setobject!(vis["env/roadway$k"]["road"], road, MeshPhongMaterial(color=vis_opts.road_color))
	setobject!(vis["env/roadway$k"]["line"], line, MeshPhongMaterial(color=colorant"yellow"))
	setobject!(vis["env/roadway$k"]["left_bound"], left_bound, MeshPhongMaterial(color=vis_opts.bound_color))
	setobject!(vis["env/roadway$k"]["right_bound"], right_bound, MeshPhongMaterial(color=vis_opts.bound_color))
    return nothing
end

function set_roadway!(vis::Visualizer, opts::MergingRoadwayOptions{T};
    vis_opts::RoadwayVisualizationOptions=RoadwayVisualizationOptions()) where {T}
    clean!(vis)
    ll = opts.lane_length
	lw = opts.lane_width
	mp = opts.merging_point
	θ = opts.merging_angle

	ϵ = vis_opts.ϵ
	liw = vis_opts.line_width
	bw = vis_opts.bound_width
	bh = vis_opts.bound_height

	road = Rect3D(Vec(0.0, -lw, -ϵ), Vec(ll, 2lw, ϵ))
	line = Rect3D(Vec(0.0, -liw/2, 0.0), Vec(ll, liw, liw))
	left_bound = Rect3D(Vec(0.0, lw, -ϵ), Vec(ll, bw, ϵ+bh))
	right_bound = Rect3D(Vec(mp+lw*tan(θ), -lw-bw, -ϵ), Vec(ll-mp-lw*tan(θ), bw, ϵ+bh))
	merging_road = Rect3D(Vec(0., 0., -ϵ), Vec(mp, lw/cos(θ), ϵ))
	merging_right_bound = Rect3D(Vec(0., 0., -ϵ), Vec(mp, bw, ϵ+bh))

    setobject!(vis["env/roadway"]["road"], road, MeshPhongMaterial(color=vis_opts.road_color))
	setobject!(vis["env/roadway"]["line"], line, MeshPhongMaterial(color=colorant"yellow"))
	setobject!(vis["env/roadway"]["left_bound"], left_bound, MeshPhongMaterial(color=vis_opts.bound_color))
	setobject!(vis["env/roadway"]["right_bound"], right_bound, MeshPhongMaterial(color=vis_opts.bound_color))
	setobject!(vis["env/roadway"]["merging_road"], merging_road, MeshPhongMaterial(color=vis_opts.road_color))
	setobject!(vis["env/roadway"]["merging_right_bound"], merging_right_bound, MeshPhongMaterial(color=vis_opts.bound_color))

	tr_road = compose(Translation(mp, 0., 0.), LinearMap(AngleAxis(θ-pi, 0, 0, 1)))
	tr_bound = compose(Translation(mp+lw*tan(θ), -lw, 0.), LinearMap(AngleAxis(θ-pi, 0, 0, 1)))
	settransform!(vis["env/roadway"]["merging_road"], tr_road)
	settransform!(vis["env/roadway"]["merging_right_bound"], tr_bound)
	return nothing
end

function set_roadway!(vis::Visualizer, opts::FourIntersectionRoadwayOptions{T};
	vis_opts::RoadwayVisualizationOptions=RoadwayVisualizationOptions()) where {T}
	clean!(vis)
	ll = opts.lane_length
    lw = opts.lane_width

	ϵ = vis_opts.ϵ
	liw = vis_opts.line_width
	bw = vis_opts.bound_width
	bh = vis_opts.bound_height

	road_horizontal = Rect3D(Vec(-ll/2, -lw, -ϵ), Vec(ll, 2lw, ϵ))
	line_horizontal_1 = Rect3D(Vec(-ll/2, -liw/2, 0.0), Vec(ll/2 - lw, liw, liw))
	line_horizontal_2 = Rect3D(Vec(lw, -liw/2, 0.0), Vec(ll/2 - lw, liw, liw))

	road_vertical = Rect3D(Vec(-lw, -ll/2, -ϵ), Vec(2lw, ll, ϵ))
	line_vertical_1 = Rect3D(Vec(-liw/2, -ll/2, 0.0), Vec(liw, ll/2 - lw, liw))
	line_vertical_2 = Rect3D(Vec(-liw/2, lw, 0.0), Vec(liw, ll/2 - lw, liw))

	horizontal_left_bound_w = Rect3D(Vec(-ll/2, lw, -ϵ), Vec(ll/2 - lw, bw, ϵ+bh))
    horizontal_right_bound_w = Rect3D(Vec(-ll/2, -lw-bw, -ϵ), Vec(ll/2 - lw, bw, ϵ+bh))

	horizontal_left_bound_e = Rect3D(Vec(lw, lw, -ϵ), Vec(ll/2 - lw, bw, ϵ+bh))
    horizontal_right_bound_e = Rect3D(Vec(lw, -lw-bw, -ϵ), Vec(ll/2 - lw, bw, ϵ+bh))

	# horizontal_left_bound = Rect3D(Vec(-ll/2, lw, -ϵ), Vec(ll, bw, ϵ+bh))
    # horizontal_right_bound = Rect3D(Vec(-ll/2, -lw-bw, -ϵ), Vec(ll, bw, ϵ+bh))

	vert_left_bound_n = Rect3D(Vec(lw, -ll/2, -ϵ), Vec(bw, ll/2 -lw, ϵ+bh))
	vert_right_bound_n = Rect3D(Vec(-lw-bw, -ll/2, -ϵ), Vec(bw, ll/2 - lw, ϵ+bh))

	vert_left_bound_s = Rect3D(Vec(lw, lw, -ϵ), Vec(bw, ll/2 - lw, ϵ+bh))
	vert_right_bound_s = Rect3D(Vec(-lw-bw, lw, -ϵ), Vec(bw, ll/2 - lw, ϵ+bh))

	stop_line_west = Rect3D(Vec(-lw, -lw, 0.0), Vec(liw, lw, liw))
	stop_line_east = Rect3D(Vec(lw, 0, 0.0), Vec(liw, lw, liw))
	stop_line_north = Rect3D(Vec(0, -lw, 0.0), Vec(lw, liw, liw))
	stop_line_south = Rect3D(Vec(-lw, lw, 0.0), Vec(lw, liw, liw))

    setobject!(vis["env/roadway"]["road_horizontal"], road_horizontal, MeshPhongMaterial(color=vis_opts.road_color))
	setobject!(vis["env/roadway"]["line_horizontal_1"], line_horizontal_1, MeshPhongMaterial(color=colorant"yellow"))
	setobject!(vis["env/roadway"]["line_horizontal_2"], line_horizontal_2, MeshPhongMaterial(color=colorant"yellow"))
	setobject!(vis["env/roadway"]["horizontal_left_bound_w"], horizontal_left_bound_w, MeshPhongMaterial(color=vis_opts.bound_color))
	setobject!(vis["env/roadway"]["horizontal_right_bound_w"], horizontal_right_bound_w, MeshPhongMaterial(color=vis_opts.bound_color))
	setobject!(vis["env/roadway"]["horizontal_left_bound_e"], horizontal_left_bound_e, MeshPhongMaterial(color=vis_opts.bound_color))
	setobject!(vis["env/roadway"]["horizontal_right_bound_e"], horizontal_right_bound_e, MeshPhongMaterial(color=vis_opts.bound_color))

	setobject!(vis["env/roadway"]["road_vertical"], road_vertical, MeshPhongMaterial(color=vis_opts.road_color))
	setobject!(vis["env/roadway"]["line_vertical_1"], line_vertical_1, MeshPhongMaterial(color=colorant"yellow"))
	setobject!(vis["env/roadway"]["line_vertical_2"], line_vertical_2, MeshPhongMaterial(color=colorant"yellow"))
	setobject!(vis["env/roadway"]["vertical_left_bound_n"], vert_left_bound_n, MeshPhongMaterial(color=vis_opts.bound_color))
	setobject!(vis["env/roadway"]["vertical_right_bound_n"], vert_right_bound_n, MeshPhongMaterial(color=vis_opts.bound_color))

	setobject!(vis["env/roadway"]["vertical_left_bound_s"], vert_left_bound_s, MeshPhongMaterial(color=vis_opts.bound_color))
	setobject!(vis["env/roadway"]["vertical_right_bound_s"], vert_right_bound_s, MeshPhongMaterial(color=vis_opts.bound_color))

	setobject!(vis["env/roadway"]["stop_line_west"], stop_line_west, MeshPhongMaterial(color=colorant"white"))
	setobject!(vis["env/roadway"]["stop_line_east"], stop_line_east, MeshPhongMaterial(color=colorant"white"))
	setobject!(vis["env/roadway"]["stop_line_north"], stop_line_north, MeshPhongMaterial(color=colorant"white"))
	setobject!(vis["env/roadway"]["stop_line_south"], stop_line_south, MeshPhongMaterial(color=colorant"white"))
    return nothing
end


text:
##
using WavePropBase
using WavePropBase.Geometry
using WavePropBase.Interpolation
using WavePropBase.Trees
using StaticArrays
using Test

# recursively check that all point in a cluster tree are in the bounding box
function test_cluster_tree(clt)
    bbox = clt.container
    for iloc in clt.index_range
        x     = clt._elements[iloc]
        x ∈ bbox || (return false)
    end
    if !isroot(clt)
        clt ∈ clt.parent.children || (return false)
    end
    if !isleaf(clt)
        for child in clt.children
            test_cluster_tree(child) || (return false)
        end
    end
    return true
end

@testset "ClusterTree" begin
    @testset "1d" begin
        points    = SVector.([4,3,1,2,5,-1.0])
        splitter  = GeometricSplitter(nmax=1)
        clt = ClusterTree(points,splitter)
        @test sortperm(points) == clt.loc2glob
        splitter  = GeometricMinimalSplitter(nmax=1)
        clt = ClusterTree(points,splitter)
        @test sortperm(points) == clt.loc2glob
        splitter  = CardinalitySplitter(nmax=1)
        clt = ClusterTree(points,splitter)
        @test sortperm(points) == clt.loc2glob
        splitter  = PrincipalComponentSplitter(nmax=1)
        clt = ClusterTree(points,splitter)
        @test sortperm(points) == clt.loc2glob
        splitter  = DyadicSplitter(nmax=1)
        clt = ClusterTree(points,splitter)
        @test sortperm(points) == clt.loc2glob
    end

    @testset "2d" begin
        points    = rand(SVector{2,Float64},1000)
        splitter  = GeometricSplitter(nmax=1)
        clt = ClusterTree(points,splitter)
        @test test_cluster_tree(clt)
        splitter  = GeometricMinimalSplitter(nmax=32)
        clt = ClusterTree(points,splitter)
        @test test_cluster_tree(clt)
        splitter   = CardinalitySplitter(nmax=32)
        clt = ClusterTree(points,splitter)
        @test test_cluster_tree(clt)
        splitter   = PrincipalComponentSplitter(nmax=32)
        clt = ClusterTree(points,splitter)
        @test test_cluster_tree(clt)
        splitter   = DyadicSplitter(nmax=32)
        clt = ClusterTree(points,splitter)
        @test test_cluster_tree(clt)
    end

    @testset "3d" begin
        points = rand(SVector{3,Float64},1000)
        splitter  = GeometricSplitter(nmax=32)
        clt = ClusterTree(points,splitter)
        @test test_cluster_tree(clt)
        splitter   = GeometricMinimalSplitter(nmax=32)
        clt = ClusterTree(points,splitter)
        @test test_cluster_tree(clt)
        splitter  = CardinalitySplitter(nmax=32)
        clt = ClusterTree(points,splitter)
        @test test_cluster_tree(clt)
        splitter   = PrincipalComponentSplitter(nmax=32)
        clt = ClusterTree(points,splitter)
        @test test_cluster_tree(clt)
        splitter   = DyadicSplitter(nmax=32)
        clt = ClusterTree(points,splitter)
        @test test_cluster_tree(clt)
    end
end

ext: jl

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
export NN,getNN,initTheta
# using TimerOutputs
# to = TimerOutput()
"""
NN Neural Network block

 Y_k+1 = layer{k}(theta{k},Y_k)
"""
mutable struct NN{T, TQ <: Union{Array{T,2},UniformScaling{Int}}} <: AbstractMeganetElement{T}
    layers   ::Array{AbstractMeganetElement{T}, 1} # layers of Neural Network, cell array
    outTimes ::Array{Int,1}
    Q        :: TQ
end

function getNN(layers::Array{AbstractMeganetElement{T}},outTimes=eye(Int,length(layers))[:,end],Q=I) where {T <: Number}
	nt   = length(layers)
    nout = nFeatOut(layers[1])

    for k=2:nt
        if nFeatIn(layers[k]) != nout
            error("Dim. of input features of block $k does not match dim. of output features of block $(k-1)");
        end
        nout = nFeatOut(layers[k])
    end
	return NN(layers,outTimes,Q);
end

import Base.display
function display(this::NN)
    println("-- Neural Network --")
    println("nLayers: \t $(length(this.layers))")
    println("nFeatIn: \t $(nFeatIn(this))")
    println("nFeatOut: \t $(nFeatOut(this))")
    println("nTheta: \t $(nTheta(this))")
end

# ---------- counting thetas, input and output features -----
function nTheta(this::NN)
    n::Int = 0;
    for k=1:length(this.layers)
        n += nTheta(this.layers[k]);
    end
    return n
end
nFeatIn(this::NN)   = nFeatIn(this.layers[1])
nFeatOut(this::NN)::Int = nFeatOut(this.layers[end])

function nDataOut(this::NN)
    n=0;
    for k=1:length(this.layers)
        n = n+this.outTimes[k]* nFeatOut(this.layers[k]);
    end
end

function initTheta(this::NN{T}) where {T <: Number}
    theta = zeros(T,0)
    for k=1:length(this.layers)
        theta = [theta; vec(initTheta(this.layers[k]))]
    end
    return convert(Array{T},theta)
end

# --------- forward problem ----------
function apply(this::NN{T},theta::Array{T},Y::Array{T,2},tmp,doDerivative=true) where {T<:Number}

    nex = div(length(Y),nFeatIn(this))::Int
    nt = length(this.layers)

    if isempty(tmp) #TODO Will have to make sure size of Y doesnt change
        tmp = Array{Any}(nt+1,2)
    end

    if doDerivative
        if isassigned(tmp,1,1)
            #tmp[1,1] .= Y This does not work, need to hack like below :)
            tmp11 = tmp[1,1]
            tmp11 .= Y
        else
            tmp[1,1] = copy(Y)
        end
    end

    Ydata::Array{T,2} = zeros(T,0,nex)
    cnt = 0
    for i=1:nt
        ni = nTheta(this.layers[i])::Int
        if !isassigned(tmp,i,2)
            tmp[i,2] = Array{Any}(0)
        end
        Yd::Array{T,2}, Y, tmp[i,2] = apply(this.layers[i],theta[cnt+(1:ni)],Y,tmp[i,2],doDerivative)

        if this.outTimes[i]==1
            Ydata = [Ydata; this.Q*Yd]
        end
        if doDerivative
            if isassigned(tmp,i+1,1)
                tmp1 = tmp[i+1,1]
                tmp1 .= Y
            else
                tmp[i+1,1] = copy(Y)
            end
        end
        cnt = cnt + ni
    end
    return Ydata,Y,tmp
end

# -------- Jacobian matvecs --------
function JYmv(this::NN{T},dY::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}
    nex = div(length(Y),nFeatIn(this))
    nt = length(this.layers)
    cnt = 0
    dYdata = zeros(T,0,nex)
    for i=1:nt
        ni = nTheta(this.layers[i])
        dY = JYmv(this.layers[i],dY,theta[cnt+(1:ni)],tmp[i,1],tmp[i,2])[2]
        if this.outTimes[i]==1
            dYdata = [dYdata; this.Q*dY]
        end
        cnt = cnt+ni
    end
    return dYdata, dY
end

function  Jmv(this::NN{T},dtheta::Array{T},dY::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}
    nex = div(length(Y),nFeatIn(this))
    nt = length(this.layers);
    if isempty(dY)
        dY = 0*Y
    end

    dYdata = zeros(T,0,nex)
    cnt = 0
    for i=1:nt
        ni = nTheta(this.layers[i])
        dY = Jmv(this.layers[i],dtheta[cnt+(1:ni)],dY,theta[cnt+(1:ni)],
                tmp[i,1],tmp[i,2])[2]
        if this.outTimes[i]==1
            dYdata = [dYdata; this.Q*dY]
        end
        cnt = cnt+ni
    end
    return dYdata,dY
end

# -------- Jacobian' matvecs --------
function JYTmv(this::NN{T},Wdata::Array{T},W::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}

    nex = div(length(Y),nFeatIn(this));
    if !isempty(Wdata)
        Wdata = reshape(Wdata,:,nex);
    end
    if isempty(W)
        W = zero(T)
    elseif length(W)>1
        W     = reshape(W,:,nex)
    end
    nt = length(this.layers)

    cnt = 0; cnt2 = 0;
    for i=nt:-1:1
        ni = nTheta(this.layers[i])
        if this.outTimes[i]==1
            nn = nFeatOut(this.layers[i])
            W = W + this.Q'*Wdata[end-cnt2-nn+1:end-cnt2,:]
            cnt2 = cnt2 + nn
        end
        W  = JYTmv(this.layers[i], W,(T)[],theta[end-cnt-ni+1:end-cnt],
                    tmp[i,1],tmp[i,2])
        cnt = cnt+ni
    end
    return W
end

function JthetaTmv(this::NN{T},Wdata::Array{T},W::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}
	return JTmv(this,Wdata,W,theta,Y,tmp)[1]; # TODO: Why calculating both, Can be more efficient?
end

function JTmv(this::NN{T},Wdata::Array{T},Win::Array{T},theta::Array{T},Y::Array{T},tmp)::Tuple{Array{T,1},Array{T,1}} where {T <: Number}
    # WOW THIS IS HACKED BIG TIME. Need to find a way to type stabalize W (not ez)
    #TODO: Make this type stable - Some internals are not stable
    nex = div(length(Y),nFeatIn(this)::Int)

    if size(Wdata,1)>0
        Wdata = reshape(Wdata,:,nex)
    end

    if length(Win)==0
        W = zeros(T,nFeatOut(this),nex)
    else
        W = reshape(Win,:,nex)
    end
    dtheta = zero(T)*theta
    nt = length(this.layers)

    cnt = 0; cnt2 = 0
    for i=nt:-1:1
        if this.outTimes[i]==1
            nn = nFeatOut(this.layers[i])::Int
            W += this.Q'*Wdata[end-cnt2-nn+1:end-cnt2,:]
            cnt2 = cnt2 + nn
        end
        ni     = nTheta(this.layers[i])::Int

        dmbi,W = JTmv(this.layers[i],W,zeros(T,0),theta[end-cnt-ni+1:end-cnt],tmp[i,1],tmp[i,2])
        dtheta[end-cnt-ni+1:end-cnt]  = dmbi
        cnt += ni
    end

    return  vec(dtheta), vec(W)
end


text:
#=

    This is a port of the SPPL model clickgraph.pynb
    """
    ClickGraph is an information-retrieval model used to compute the posterior probability that a pair
    of items A and B are similar, conditioned on a user's access patterns to these items. This model is
    used, for example, to measure the quality of a search engine, where A and B represent two URLs that
    are returned to a given search query. A search engine is said to be good if the probability that
    URLs A and B are similar is high.

    The probabilistic program below implements the ClickGraph model, for five pairs of URLs and one
    user's click patterns.
    """

    This model:
        parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
            Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

         p_similar    0.6806    0.2866     0.0091    0.0320    88.3928    1.0026       16.3690

    * click_url_A_post (cf [1,1,0,1,0])
            parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
               Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

  click_url_A_post[1]    0.7950    0.4039     0.0128    0.0229   233.1889    0.9992       43.1831
  click_url_A_post[2]    0.7400    0.4389     0.0139    0.0355   274.5914    1.0121       50.8503
  click_url_A_post[3]    0.2370    0.4255     0.0135    0.0197   470.6780    0.9999       87.1626
  click_url_A_post[4]    0.7170    0.4507     0.0143    0.0235   376.7671    0.9990       69.7717
  click_url_A_post[5]    0.2600    0.4389     0.0139    0.0178   425.0068    0.9997       78.7050

    * click_url_B_post (cf [1,1,0,1,0])
           parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
               Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

  click_url_B_post[1]    0.6670    0.4715     0.0149    0.0226   267.1792    0.9990       44.5596
  click_url_B_post[2]    0.7630    0.4255     0.0135    0.0198   402.3125    0.9994       67.0968
  click_url_B_post[3]    0.2880    0.4531     0.0143    0.0180   410.6459    1.0010       68.4866
  click_url_B_post[4]    0.7230    0.4477     0.0142    0.0118   479.3376    0.9994       79.9429
  click_url_B_post[5]    0.2740    0.4462     0.0141    0.0207   513.0703    0.9991       85.5688

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function clickgraph_sppl(click_url_A=missing, click_url_B=missing)
    if click_url_A === missing
        n = 5
        click_url_A = tzeros(n)
        click_url_B = tzeros(n)
    else
        n = length(click_url_A)
    end

    p_similar ~ Beta(1, 1)   # Prior probabilty URLs A and B are similar. (in [0,1])

    # Five file pairs

    similar         = tzeros(n)    # Are url_A[i] and url_B[i] similar? (binary)
    p_click_url_A   = tzeros(n)    # Probability user clicks url_A[i]. (in [0,1])
    p_click_url_B   = tzeros(n)    # Probability user clicks url_B[i]. (in [0,1])
    # click_url_A     = tzeros(n)    # Did user click url_A[i]? (binary)
    # click_url_B     = tzeros(n)    # Did user click url_B[i]? (binary)

    for i in 1:n
        similar[i] ~ Bernoulli(p_similar)

        if similar[i] == true
            p_click_url_A[i] ~ Uniform(0,1)
            p_click_url_B[i] ~ Dirac(p_click_url_A[i])
        else
            p_click_url_A[i] ~ Uniform(0,1)
            p_click_url_B[i] ~ Uniform(0,1)
        end

        click_url_A[i] ~ Bernoulli(p_click_url_A[i])
        click_url_B[i] ~ Bernoulli(p_click_url_B[i])
    end

    click_url_A_post = tzeros(n)
    click_url_B_post = tzeros(n)
    for i in 1:n
        click_url_A_post[i] ~ Bernoulli(p_click_url_A[i])
        click_url_B_post[i] ~ Bernoulli(p_click_url_B[i])
    end
end

click_url_A = [1,1,0,1,0]
click_url_B = [1,1,0,1,0]
model = clickgraph_sppl(click_url_A, click_url_B)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5),  1_000)
# chns = sample(model, PG(5),  MCMCThreads(), 10_000, 4)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns)
# display(plot(chns))

display(group(chns,:click_url_A_post))
display(group(chns,:click_url_B_post))

println("\nPredict:")
mpred = clickgraph_sppl(missing,missing)
pred = predict(mpred, chns)

display(pred)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#####
##### Functions for non-rotating models
#####

@inline x_f_cross_U(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U) where FT = zero(FT)
@inline y_f_cross_U(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U) where FT = zero(FT)
@inline z_f_cross_U(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U) where FT = zero(FT)


text:
################################################################################
# Matrix containing the number of observations and the condition
# number per pixel

@doc raw"""
    mutable struct NobsMatrixElement

This structure encodes the condition matrix for a pixel in the map. It
is essentially matrix M_p in Eq. (10) of KurkiSuonio2009, but with
steroids, as it implements the following fields:

- `m`: the 3×3 matrix in Eq. (10)
- `invm`: the 3×3 inverse of `m`
- `invcond`: the inverse of the condition number *rcond* for `m`; this
  field is useful to check whether the samples in a TOD and the attack
  angles are sufficient to constrain a unique solution for I, Q, and U
  for that pixel.
- `neglected`: a Boolean flag telling whether the pixel was skipped or
  not during map-making (this usually happens if `invcond` is too
  small)

A typical use case is an array of `NobsMatrixElement` objects in an
array that get updated via [`update_nobs_matrix!`](@ref), like in the
following example:

```julia
NPIX = 10   # Number of pixels in the map
nobs_matrix = [NobsMatrixElement() for i in 1:NPIX]

# The variables psi_angle, sigma_squared, pixidx, and flagged
# must have been defined somewhere else; they contain the TOD
update_nobs_matrix!(nobs_matrix, psi_angle, sigma_squared, pixidx, flagged)
```

"""
mutable struct NobsMatrixElement{T <: Real}
    # We do not define "m" as symmetric, as we want to update its
    # fields one by one; however, we only consider the upper part
    m::Array{T, 2}
    invm::Symmetric{T, Array{T, 2}}
    invcond::T

    # If "true", do not consider this pixel in the solution, as it was
    # not covered enough by the scanning strategy
    neglected::Bool

    NobsMatrixElement{T}() where {T <: Real} = new(
        Symmetric(zeros(T, 3, 3)),
        Symmetric(zeros(T, 3, 3)),
        0.0,
        true,
    )
end

function Base.show(io::IO, nobsmatr::NobsMatrixElement{T}) where {T <: Real}
    print(io, "NobsMatrixElement(invcond => $(nobsmatr.invcond), neglected => $(nobsmatr.neglected))")
end

function Base.show(io::IO, ::MIME"text/plain", nobsmatr::NobsMatrixElement{T}) where {T <: Real}
    println(
        io,
        @sprintf(
            """NobsMatrixElement:
Inverse condition number.......... %e
Neglected......................... %s
Matrix:
    %.5f   %.5f   %.5f
    %.5f   %.5f   %.5f
    %.5f   %.5f   %.5f
Inverse matrix:
    %.5f   %.5f   %.5f
    %.5f   %.5f   %.5f
    %.5f   %.5f   %.5f
    """,
            nobsmatr.invcond,
            nobsmatr.neglected,
            nobsmatr.m[1, 1], nobsmatr.m[1, 2], nobsmatr.m[1, 3],
            nobsmatr.m[2, 1], nobsmatr.m[2, 2], nobsmatr.m[2, 3],
            nobsmatr.m[3, 1], nobsmatr.m[3, 2], nobsmatr.m[3, 3],
            nobsmatr.invm[1, 1], nobsmatr.invm[1, 2], nobsmatr.invm[1, 3],
            nobsmatr.invm[2, 1], nobsmatr.invm[2, 2], nobsmatr.invm[2, 3],
            nobsmatr.invm[3, 1], nobsmatr.invm[3, 2], nobsmatr.invm[3, 3],
        )
    )
end

@doc raw"""
    update_nobs!(nobs::NobsMatrixElement{T}; threshold = 1e-7)

This function makes sure that all the elements in `nobs` are
synchronized. It should be called whenever the field `nobs.m` (matrix
M_p in Eq. 10 of KurkiSuonio2009) has changed.

"""
function update_nobs!(nobs::NobsMatrixElement{T}; threshold = 1e-7) where {T <: Real}
    c = cond(nobs.m)

    if isfinite(c)
        nobs.invm = inv(Symmetric(nobs.m))
        nobs.invcond = 1 / c
        nobs.neglected = nobs.invcond < threshold
    else
        nobs.invcond = 0.0
        nobs.neglected = true
    end
end

@doc raw"""
    update_nobs_matrix!(nobs_matrix::Vector{NobsMatrixElement{T}}, psi_angle, sigma_squared, pixidx, flagged)

Apply Eq. (10) of KurkiSuonio2009 iteratively on the samples of a TOD
to update matrices ``M_p`` in `nobs_matrix`. The meaning of the
parameters is the following:

- `nobs_matrix` is the structure that gets updated by this function
- `psi_angle` is an array of `N` elements, containing the polarization
  angles (in radians)
- `sigma_squared` is an array of `N` elements, each containing the
  value of σ^2 for the samples in the TOD
- `pixidx` is an array of `N` elements, containing the pixel index
  observed by the TOD
- `flagged` is a Boolean array of `N` elements; `true` means that the
  sample in the TOD should not be used to produce the map. This can be
  used to produce jackknives and to neglect moving objects in the TOD

"""
function update_nobs_matrix!(
    nobs_matrix::Vector{NobsMatrixElement{T}},
    psi_angle,
    sigma_squared,
    pixidx,
    flagged;
    comm = nothing,
) where {T <: Real}

    @assert length(psi_angle) == length(sigma_squared)
    @assert length(psi_angle) == length(pixidx)
    @assert length(psi_angle) == length(flagged)

    for (idx, curpix, curpsi, curflagged) in zip(
        1:length(pixidx),
        pixidx,
        psi_angle,
        flagged,
    )
        curflagged && continue

        constant = 1 / sigma_squared[idx]

        sin_term, cos_term = sincos(2 * curpsi)
        nobs_matrix[curpix].m[1, 1] += constant
        nobs_matrix[curpix].m[1, 2] += cos_term * constant
        nobs_matrix[curpix].m[1, 3] += sin_term * constant

        nobs_matrix[curpix].m[2, 2] += cos_term^2 * constant
        nobs_matrix[curpix].m[2, 3] += cos_term * sin_term * constant

        nobs_matrix[curpix].m[3, 3] += sin_term^2 * constant
    end

    if use_mpi() && comm != nothing
        for curpix in eachindex(nobs_matrix)
            # TODO: repeated calls to MPI_AllReduce are tremendously slow!
            MPI.Allreduce!(MPI.IN_PLACE, nobs_matrix[curpix].m, +, comm)
        end
    end
end


text:
using JuliaMaps
using Test

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


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libbtf"], :libbtf),
    LibraryProduct(prefix, String["libsundials_sunlinsolspfgmr"], :libsundials_sunlinsolspfgmr),
    LibraryProduct(prefix, String["libsundials_ida"], :libsundials_ida),
    LibraryProduct(prefix, String["libsundials_cvode"], :libsundials_cvode),
    LibraryProduct(prefix, String["libsundials_cvodes"], :libsundials_cvodes),
    LibraryProduct(prefix, String["libcolamd"], :libcolamd),
    LibraryProduct(prefix, String["libsundials_sunmatrixdense"], :libsundials_sunmatrixdense),
    LibraryProduct(prefix, String["libsundials_sunlinsolspbcgs"], :libsundials_sunlinsolspbcgs),
    LibraryProduct(prefix, String["libsundials_idas"], :libsundials_idas),
    LibraryProduct(prefix, String["libsundials_nvecserial"], :libsundials_nvecserial),
    LibraryProduct(prefix, String["libsundials_sunlinsoldense"], :libsundials_sunlinsoldense),
    LibraryProduct(prefix, String["libsundials_sunlinsolspgmr"], :libsundials_sunlinsolspgmr),
    LibraryProduct(prefix, String["libsundials_sunlinsolpcg"], :libsundials_sunlinsolpcg),
    LibraryProduct(prefix, String["libsundials_sunlinsolsptfqmr"], :libsundials_sunlinsolsptfqmr),
    LibraryProduct(prefix, String["libsundials_sunlinsolklu"], :libsundials_sunlinsolklu),
    LibraryProduct(prefix, String["libsundials_sunmatrixsparse"], :libsundials_sunmatrixsparse),
    LibraryProduct(prefix, String["libsundials_sunlinsolband"], :libsundials_sunlinsolband),
    LibraryProduct(prefix, String["libsundials_sunmatrixband"], :libsundials_sunmatrixband),
    LibraryProduct(prefix, String["libsundials_kinsol"], :libsundials_kinsol),
    LibraryProduct(prefix, String["libsundials_arkode"], :libsundials_arkode),
    LibraryProduct(prefix, String["libklu"], :libklu),
    LibraryProduct(prefix, String["libsuitesparseconfig"], :libsuitesparseconfig),
    LibraryProduct(prefix, String["libamd"], :libamd),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaDiffEq/SundialsBuilder/releases/download/sundials-3.1.1.3"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/Sundials.aarch64-linux-gnu.tar.gz", "f555ce54c06e5b9f4323fe05e0491d067bc097924795f9a42c83f3980424cf10"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/Sundials.arm-linux-gnueabihf.tar.gz", "e5642827eea06b469aed69bc4ce1d485905c17aa645af411f8bc2cc9475f5092"),
    Linux(:i686, :glibc) => ("$bin_prefix/Sundials.i686-linux-gnu.tar.gz", "b180b28986a076525d938aab0a5ebf87147bee8c6b12dc1f9634f2c1d576cdde"),
    Windows(:i686) => ("$bin_prefix/Sundials.i686-w64-mingw32.tar.gz", "c457fe74ebe6d2fdd74b6ea8346fc23ac74ac04e867ce987ae32d03ccf7beac8"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/Sundials.powerpc64le-linux-gnu.tar.gz", "7307849df35003f25511686183ecb20628ba1f46dac1e86b2fdc4ecbe88fc602"),
    MacOS(:x86_64) => ("$bin_prefix/Sundials.x86_64-apple-darwin14.tar.gz", "8099e1391fffc69cc88a8ff3a12ea783d07663e3b2f13902b5b5155527616c51"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/Sundials.x86_64-linux-gnu.tar.gz", "a8393cc1fbbfd468c174234b8d6aad1a82cd3bd76f597c956da5eeea9e30b9f6"),
    Windows(:x86_64) => ("$bin_prefix/Sundials.x86_64-w64-mingw32.tar.gz", "d0f8aa079121a74b8315541fe05a208f4581f4c9fe6208d81447eb540a6199d0"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
"""
    static_slice(A, sizes)
    static_slice(A, code)

Slice array `A` into an `AbstractArray{StaticArray}`,
with the slices always belonging to the first few indices, i.e. `code = (:,:,...,*,*)`.
Prefers to get `sizes = Size(size(A)[1:ncolons])` for stability,
but if given code instead it will do what it can.

Typically produces `reshape(reinterpret(SArray...`; copy this to get `Array{SArray}`,
or call `static_slice(A, sizes, false)` to omit the reshape.
"""
@inline function static_slice(A::AbstractArray{T,N}, sizes::Size{tup}) where {T,N,tup}
    IN = length(tup)
    for d in 1:IN
        size(A,d) == tup[d] || throw(DimensionMismatch("cannot slice array of size $(size(A)) using Size$tup"))
        first(axes(A,d)) == 1 || throw(ArgumentError("cannot creat static slices with offset indices"))
    end
    IT = SArray{Tuple{tup...}, T, IN, prod(tup)}
    finalaxes = axes(A)[1+IN:end]
    reshape(reinterpret(IT, vec(A)), finalaxes)
end

function static_slice(A::AbstractArray, code::Tuple)
    iscodesorted(code) || error("expected to make slices of left-most dimensions")
    tup = ntuple(d -> size(A,d), countcolons(code))
    static_slice(A, Size(tup))
end

"""
    static_glue(A)

Glues the output of `static_slice` back into one array, again with `code = (:,:,...,*,*)`.
For `MArray` slices, which can't be reinterpreted, this reverts to `red_glue`.
"""
@inline function static_glue(A::AbstractArray{IT}) where {IT<:SArray}
    finalaxes = (axes(IT)..., axes(A)...)
    reshape(reinterpret(eltype(IT), A), finalaxes)
end

function static_glue(A::AbstractArray{IT}) where {IT<:MArray}
    LazyStack.stack(A)
end

# This is now used only by maybestaticsizes...
iscodesorted(code::Tuple{}) = true
iscodesorted(code::Tuple{Any}) = true
iscodesorted(code::Tuple) = !(code[1]===(*) && code[2]===(:)) && iscodesorted(code[2:end])


text:
using StaticNumbers
export remainder_DGModel

import ..BalanceLaws:
    vars_state,
    init_state_prognostic!,
    init_state_auxiliary!,
    flux_first_order!,
    flux_second_order!,
    compute_gradient_flux!,
    compute_gradient_argument!,
    source!,
    transform_post_gradient_laplacian!,
    wavespeed,
    boundary_conditions,
    boundary_state!,
    update_auxiliary_state!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

"""
    RemBL(
        main::BalanceLaw,
        subcomponents::Tuple,
        maindir::Direction,
        subsdir::Tuple,
    )

Balance law for holding remainder model information. Direction is put here since
direction is so intertwined with the DGModel_kernels, that it is easier to hande
this in this container.
"""
struct RemBL{M, S, MD, SD} <: BalanceLaw
    main::M
    subs::S
    maindir::MD
    subsdir::SD
end

"""
    rembl_has_subs_direction(
        direction::Direction,
        rem_balance_law::RemBL,
    )

Query whether the `rem_balance_law` has any subcomponent balance laws operating
in the direction `direction`
"""
@generated function rembl_has_subs_direction(
    ::Dir,
    ::RemBL{MainBL, SubsBL, MainDir, SubsDir},
) where {Dir <: Direction, MainBL, SubsBL, MainDir, SubsDir <: Tuple}
    if Dir in SubsDir.types
        return :(true)
    else
        return :(false)
    end
end

"""
    remainder_DGModel(
        maindg::DGModel,
        subsdg::NTuple{NumModels, DGModel};
        numerical_flux_first_order,
        numerical_flux_second_order,
        numerical_flux_gradient,
        state_auxiliary,
        state_gradient_flux,
        states_higher_order,
        diffusion_direction,
        modeldata,
    )

Constructs a `DGModel` from the `maindg` model and the tuple of `subsdg` models.
The concept of a remainder model is that it computes the contribution of the
model after subtracting all of the subcomponents.

By default the numerical fluxes are set to be a tuple of the main models
numerical flux and the splitting is done at the PDE level (e.g., the remainder
model is calculated prior to discretization). If instead a tuple of numerical
fluxes is passed in the main numerical flux is evaluated first and then the
subcomponent numerical fluxes are subtracted off. This is discretely different
(for the Rusanov / local Lax-Friedrichs flux) than defining a numerical flux for
the remainder of the physics model.

The other parameters are set to the value in the `maindg` component, mainly the
data and arrays are aliased to the `maindg` values.
"""
function remainder_DGModel(
    maindg::DGModel,
    subsdg::NTuple{NumModels, DGModel};
    numerical_flux_first_order = maindg.numerical_flux_first_order,
    numerical_flux_second_order = maindg.numerical_flux_second_order,
    numerical_flux_gradient = maindg.numerical_flux_gradient,
    state_auxiliary = maindg.state_auxiliary,
    state_gradient_flux = maindg.state_gradient_flux,
    states_higher_order = maindg.states_higher_order,
    diffusion_direction = maindg.diffusion_direction,
    modeldata = maindg.modeldata,
) where {NumModels}
    FT = eltype(state_auxiliary)

    # If any of these asserts fail, the remainder model will need to be extended
    # to allow for it; see `flux_first_order!` and `source!` below.
    for subdg in subsdg
        @assert number_states(subdg.balance_law, Prognostic()) <=
                number_states(maindg.balance_law, Prognostic())

        @assert number_states(subdg.balance_law, Auxiliary()) ==
                number_states(maindg.balance_law, Auxiliary())

        @assert number_states(subdg.balance_law, Gradient()) == 0
        @assert number_states(subdg.balance_law, GradientFlux()) == 0

        @assert number_states(subdg.balance_law, GradientLaplacian()) == 0
        @assert number_states(subdg.balance_law, Hyperdiffusive()) == 0

        # Do not currenlty support nested remainder models
        # For this to work the way directions and numerical fluxes are handled
        # would need to be updated.
        @assert !(subdg.balance_law isa RemBL)

        @assert number_states(subdg.balance_law, UpwardIntegrals()) == 0
        @assert number_states(subdg.balance_law, DownwardIntegrals()) == 0

        # The remainder model requires that the subcomponent direction be
        # included in the main model directions
        @assert (
            maindg.direction isa EveryDirection ||
            maindg.direction === subdg.direction
        )
    end
    balance_law = RemBL(
        maindg.balance_law,
        ntuple(i -> subsdg[i].balance_law, length(subsdg)),
        maindg.direction,
        ntuple(i -> subsdg[i].direction, length(subsdg)),
    )

    DGModel(
        balance_law,
        maindg.grid,
        numerical_flux_first_order,
        numerical_flux_second_order,
        numerical_flux_gradient,
        state_auxiliary,
        state_gradient_flux,
        states_higher_order,
        maindg.direction,
        diffusion_direction,
        modeldata,
    )
end

# Inherit most of the functionality from the main model
vars_state(rem_balance_law::RemBL, st::AbstractStateType, FT) =
    vars_state(rem_balance_law.main, st, FT)

update_auxiliary_state!(dg::DGModel, rem_balance_law::RemBL, args...) =
    update_auxiliary_state!(dg, rem_balance_law.main, args...)

nodal_update_auxiliary_state!(dg::DGModel, rem_balance_law::RemBL, args...) =
    nodal_update_auxiliary_state!(dg, rem_balance_law.main, args...)

update_auxiliary_state_gradient!(dg::DGModel, rem_balance_law::RemBL, args...) =
    update_auxiliary_state_gradient!(dg, rem_balance_law.main, args...)

integral_load_auxiliary_state!(rem_balance_law::RemBL, args...) =
    integral_load_auxiliary_state!(rem_balance_law.main, args...)

integral_set_auxiliary_state!(rem_balance_law::RemBL, args...) =
    integral_set_auxiliary_state!(rem_balance_law.main, args...)

reverse_integral_load_auxiliary_state!(rem_balance_law::RemBL, args...) =
    reverse_integral_load_auxiliary_state!(rem_balance_law.main, args...)

reverse_integral_set_auxiliary_state!(rem_balance_law::RemBL, args...) =
    reverse_integral_set_auxiliary_state!(rem_balance_law.main, args...)

transform_post_gradient_laplacian!(rem_balance_law::RemBL, args...) =
    transform_post_gradient_laplacian!(rem_balance_law.main, args...)

flux_second_order!(rem_balance_law::RemBL, args...) =
    flux_second_order!(rem_balance_law.main, args...)

compute_gradient_argument!(rem_balance_law::RemBL, args...) =
    compute_gradient_argument!(rem_balance_law.main, args...)

compute_gradient_flux!(rem_balance_law::RemBL, args...) =
    compute_gradient_flux!(rem_balance_law.main, args...)

boundary_conditions(rem_balance_law::RemBL) =
    boundary_conditions(rem_balance_law.main)

boundary_state!(nf, bc, rem_balance_law::RemBL, args...) =
    boundary_state!(nf, bc, rem_balance_law.main, args...)

init_state_auxiliary!(rem_balance_law::RemBL, args...) =
    init_state_auxiliary!(rem_balance_law.main, args...)

nodal_init_state_auxiliary!(rem_balance_law::RemBL, args...) =
    nodal_init_state_auxiliary!(rem_balance_law.main, args...)

init_state_prognostic!(rem_balance_law::RemBL, args...) =
    init_state_prognostic!(rem_balance_law.main, args...)

"""
    function flux_first_order!(
        rem_balance_law::RemBL,
        flux::Grad,
        state::Vars,
        aux::Vars,
        t::Real,
        directions,
    )

Evaluate the remainder flux by first evaluating the main flux and subtracting
the subcomponent fluxes.

Only models which have directions that are included in the `directions` tuple
are evaluated. When these models are evaluated the models underlying `direction`
is passed (not the original `directions` argument).
"""
function flux_first_order!(
    rem_balance_law::RemBL,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    ::Dirs,
) where {NumDirs, Dirs <: NTuple{NumDirs, Direction}}
    m = parent(flux)
    if rem_balance_law.maindir isa Union{Dirs.types...}
        flux_first_order!(
            rem_balance_law.main,
            flux,
            state,
            aux,
            t,
            (rem_balance_law.maindir,),
        )
    end

    flux_s = similar(flux)
    m_s = parent(flux_s)

    # Force the loop to unroll to get type stability on the GPU
    @inbounds ntuple(Val(length(rem_balance_law.subs))) do k
        Base.@_inline_meta
        if rem_balance_law.subsdir[k] isa Union{Dirs.types...}
            sub = rem_balance_law.subs[k]
            fill!(m_s, -zero(eltype(m_s)))
            flux_first_order!(
                sub,
                flux_s,
                state,
                aux,
                t,
                (rem_balance_law.subsdir[k],),
            )
            m .-= m_s
        end
    end
    nothing
end

"""
    function source!(
        rem_balance_law::RemBL,
        source::Vars,
        state::Vars,
        diffusive::Vars,
        aux::Vars,
        t::Real,
        directions,
    )

Evaluate the remainder source by first evaluating the main source and subtracting
the subcomponent sources.

Only models which have directions that are included in the `directions` tuple
are evaluated. When these models are evaluated the models underlying `direction`
is passed (not the original `directions` argument).
"""
function source!(
    rem_balance_law::RemBL,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    ::Dirs,
) where {NumDirs, Dirs <: NTuple{NumDirs, Direction}}
    m = parent(source)
    if EveryDirection() isa Union{Dirs.types...} ||
       rem_balance_law.maindir isa EveryDirection ||
       rem_balance_law.maindir isa Union{Dirs.types...}
        source!(
            rem_balance_law.main,
            source,
            state,
            diffusive,
            aux,
            t,
            (rem_balance_law.maindir,),
        )
    end

    source_s = similar(source)
    m_s = parent(source_s)

    # Force the loop to unroll to get type stability on the GPU
    ntuple(Val(length(rem_balance_law.subs))) do k
        Base.@_inline_meta
        @inbounds if EveryDirection() isa Union{Dirs.types...} ||
                     rem_balance_law.subsdir[k] isa EveryDirection ||
                     rem_balance_law.subsdir[k] isa Union{Dirs.types...}
            sub = rem_balance_law.subs[k]
            fill!(m_s, -zero(eltype(m_s)))
            source!(
                sub,
                source_s,
                state,
                diffusive,
                aux,
                t,
                (rem_balance_law.subsdir[k],),
            )
            m .-= m_s
        end
    end
    nothing
end

"""
    function wavespeed(
        rem_balance_law::RemBL,
        args...,
    )

The wavespeed for a remainder model is defined to be the difference of the
wavespeed of the main model and the sum of the subcomponents.

Note: Defining the wavespeed in this manner can result in a smaller value than
the actually wavespeed of the remainder physics model depending on the
composition of the models.
"""
function wavespeed(
    rem_balance_law::RemBL,
    nM,
    state::Vars,
    aux::Vars,
    t::Real,
    dir::Dirs,
) where {ND, Dirs <: NTuple{ND, Direction}}
    FT = eltype(state)

    ws = fill(
        -zero(FT),
        MVector{number_states(rem_balance_law.main, Prognostic()), FT},
    )
    rs = fill(
        -zero(FT),
        MVector{number_states(rem_balance_law.main, Prognostic()), FT},
    )

    # Compute the main components wavespeed
    if rem_balance_law.maindir isa Union{Dirs.types...}
        ws .= wavespeed(
            rem_balance_law.main,
            nM,
            state,
            aux,
            t,
            (rem_balance_law.maindir,),
        )
    end

    # Compute the sub components wavespeed
    for (sub, subdir) in zip(rem_balance_law.subs, rem_balance_law.subsdir)
        @inbounds if subdir isa Union{Dirs.types...}
            num_state = static(number_states(sub, Prognostic()))
            rs[static(1):num_state] .+=
                wavespeed(sub, nM, state, aux, t, (subdir,))
        end
    end

    ws .-= rs

    return ws
end

# Here the fluxes are pirated to handle the case of tuples of fluxes
import ..DGMethods.NumericalFluxes:
    NumericalFluxFirstOrder,
    numerical_flux_first_order!,
    numerical_boundary_flux_first_order!,
    normal_boundary_flux_second_order!

"""
    function numerical_flux_first_order!(
        numerical_fluxes::Tuple{
            NumericalFluxFirstOrder,
            NTuple{NumSubFluxes, NumericalFluxFirstOrder},
        },
        rem_balance_law::RemBL,
        fluxᵀn::Vars{S},
        normal_vector::SVector,
        state_prognostic⁻::Vars{S},
        state_auxiliary⁻::Vars{A},
        state_prognostic⁺::Vars{S},
        state_auxiliary⁺::Vars{A},
        t,
        directions,
    )

When the `numerical_fluxes` are a tuple and the balance law is a remainder
balance law the main components numerical flux is evaluated then all the
subcomponent numerical fluxes are evaluated and subtracted.

Only models which have directions that are included in the `directions` tuple
are evaluated. When these models are evaluated the models underlying `direction`
is passed (not the original `directions` argument).
"""
function numerical_flux_first_order!(
    numerical_fluxes::Tuple{
        NumericalFluxFirstOrder,
        NTuple{NumSubFluxes, NumericalFluxFirstOrder},
    },
    rem_balance_law::RemBL,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    ::Dirs,
) where {NumSubFluxes, S, A, Dirs <: NTuple{2, Direction}}
    # Call the numerical flux for the main model
    if rem_balance_law.maindir isa Union{Dirs.types...}
        @inbounds numerical_flux_first_order!(
            numerical_fluxes[1],
            rem_balance_law.main,
            fluxᵀn,
            normal_vector,
            state_prognostic⁻,
            state_auxiliary⁻,
            state_prognostic⁺,
            state_auxiliary⁺,
            t,
            (rem_balance_law.maindir,),
        )
    end

    # Create put the sub model fluxes
    a_fluxᵀn = parent(fluxᵀn)

    # Force the loop to unroll to get type stability on the GPU
    ntuple(Val(length(rem_balance_law.subs))) do k
        Base.@_inline_meta
        @inbounds if rem_balance_law.subsdir[k] isa Union{Dirs.types...}
            sub = rem_balance_law.subs[k]
            nf = numerical_fluxes[2][k]

            FT = eltype(a_fluxᵀn)
            num_state_prognostic = number_states(sub, Prognostic())

            a_sub_fluxᵀn = MVector{num_state_prognostic, FT}(undef)
            a_sub_state_prognostic⁻ = MVector{num_state_prognostic, FT}(undef)
            a_sub_state_prognostic⁺ = MVector{num_state_prognostic, FT}(undef)

            state_rng = static(1):static(number_states(sub, Prognostic()))
            a_sub_fluxᵀn .= a_fluxᵀn[state_rng]
            a_sub_state_prognostic⁻ .= parent(state_prognostic⁻)[state_rng]
            a_sub_state_prognostic⁺ .= parent(state_prognostic⁺)[state_rng]

            # compute this submodels flux
            fill!(a_sub_fluxᵀn, -zero(eltype(a_sub_fluxᵀn)))
            numerical_flux_first_order!(
                nf,
                sub,
                Vars{vars_state(sub, Prognostic(), FT)}(a_sub_fluxᵀn),
                normal_vector,
                Vars{vars_state(sub, Prognostic(), FT)}(
                    a_sub_state_prognostic⁻,
                ),
                state_auxiliary⁻,
                Vars{vars_state(sub, Prognostic(), FT)}(
                    a_sub_state_prognostic⁺,
                ),
                state_auxiliary⁺,
                t,
                (rem_balance_law.subsdir[k],),
            )

            # Subtract off this sub models flux
            a_fluxᵀn[state_rng] .-= a_sub_fluxᵀn
        end
    end
end

"""
    function numerical_boundary_flux_first_order!(
        numerical_fluxes::Tuple{
            NumericalFluxFirstOrder,
            NTuple{NumSubFluxes, NumericalFluxFirstOrder},
        },
        bc,
        rem_balance_law::RemBL,
        fluxᵀn::Vars{S},
        normal_vector::SVector,
        state_prognostic⁻::Vars{S},
        state_auxiliary⁻::Vars{A},
        state_prognostic⁺::Vars{S},
        state_auxiliary⁺::Vars{A},
        t,
        directions,
        args...,
    )

When the `numerical_fluxes` are a tuple and the balance law is a remainder
balance law the main components numerical flux is evaluated then all the
subcomponent numerical fluxes are evaluated and subtracted.

Only models which have directions that are included in the `directions` tuple
are evaluated. When these models are evaluated the models underlying `direction`
is passed (not the original `directions` argument).
"""
function numerical_boundary_flux_first_order!(
    numerical_fluxes::Tuple{
        NumericalFluxFirstOrder,
        NTuple{NumSubFluxes, NumericalFluxFirstOrder},
    },
    bc,
    rem_balance_law::RemBL,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    ::Dirs,
    state_prognostic1⁻::Vars{S},
    state_auxiliary1⁻::Vars{A},
) where {NumSubFluxes, S, A, Dirs <: NTuple{2, Direction}}
    # Since the fluxes are allowed to modified these we need backups so they can
    # be reset as we go
    a_state_prognostic⁺ = parent(state_prognostic⁺)
    a_state_auxiliary⁺ = parent(state_auxiliary⁺)

    a_back_state_prognostic⁺ = copy(a_state_prognostic⁺)
    a_back_state_auxiliary⁺ = copy(a_state_auxiliary⁺)

    # Call the numerical flux for the main model
    if rem_balance_law.maindir isa Union{Dirs.types...}
        @inbounds numerical_boundary_flux_first_order!(
            numerical_fluxes[1],
            bc,
            rem_balance_law.main,
            fluxᵀn,
            normal_vector,
            state_prognostic⁻,
            state_auxiliary⁻,
            state_prognostic⁺,
            state_auxiliary⁺,
            t,
            (rem_balance_law.maindir,),
            state_prognostic1⁻,
            state_auxiliary1⁻,
        )
    end

    # Create put the sub model fluxes
    a_fluxᵀn = parent(fluxᵀn)

    # Force the loop to unroll to get type stability on the GPU
    ntuple(Val(length(rem_balance_law.subs))) do k
        Base.@_inline_meta
        @inbounds if rem_balance_law.subsdir[k] isa Union{Dirs.types...}
            sub = rem_balance_law.subs[k]
            nf = numerical_fluxes[2][k]

            # reset the plus-side data
            a_state_prognostic⁺ .= a_back_state_prognostic⁺
            a_state_auxiliary⁺ .= a_back_state_auxiliary⁺

            FT = eltype(a_fluxᵀn)
            num_state_prognostic = number_states(sub, Prognostic())

            a_sub_fluxᵀn = MVector{num_state_prognostic, FT}(undef)
            a_sub_state_prognostic⁻ = MVector{num_state_prognostic, FT}(undef)
            a_sub_state_prognostic⁺ = MVector{num_state_prognostic, FT}(undef)
            a_sub_state_prognostic1⁻ = MVector{num_state_prognostic, FT}(undef)

            state_rng = static(1):static(number_states(sub, Prognostic()))
            a_sub_fluxᵀn .= a_fluxᵀn[state_rng]
            a_sub_state_prognostic⁻ .= parent(state_prognostic⁻)[state_rng]
            a_sub_state_prognostic⁺ .= parent(state_prognostic⁺)[state_rng]
            a_sub_state_prognostic1⁻ .= parent(state_prognostic1⁻)[state_rng]

            # compute this submodels flux
            fill!(a_sub_fluxᵀn, -zero(eltype(a_sub_fluxᵀn)))
            numerical_boundary_flux_first_order!(
                nf,
                bc,
                sub,
                Vars{vars_state(sub, Prognostic(), FT)}(a_sub_fluxᵀn),
                normal_vector,
                Vars{vars_state(sub, Prognostic(), FT)}(
                    a_sub_state_prognostic⁻,
                ),
                state_auxiliary⁻,
                Vars{vars_state(sub, Prognostic(), FT)}(
                    a_sub_state_prognostic⁺,
                ),
                state_auxiliary⁺,
                t,
                (rem_balance_law.subsdir[k],),
                Vars{vars_state(sub, Prognostic(), FT)}(
                    a_sub_state_prognostic1⁻,
                ),
                state_auxiliary1⁻,
            )

            # Subtract off this sub models flux
            a_fluxᵀn[state_rng] .-= a_sub_fluxᵀn
        end
    end
end

"""
    normal_boundary_flux_second_order!(nf, rem_balance_law::RemBL, args...)

Currently the main models `normal_boundary_flux_second_order!` is called. If the
subcomponents models have second order terms this would need to be updated.
"""
normal_boundary_flux_second_order!(
    nf,
    bc,
    rem_balance_law::RemBL,
    fluxᵀn::Vars{S},
    args...,
) where {S} = normal_boundary_flux_second_order!(
    nf,
    bc,
    rem_balance_law.main,
    fluxᵀn,
    args...,
)

max_forks_repo_name: jm-c/CLIMA


text:
using Base: Indices, tail, PermutedDimsArrays.PermutedDimsArray

struct SliceData{transpose,N,Axs}
    signals::NTuple{N,Signal{Int}}
    axs::Axs
end

function Base.show(io::IO, sd::SliceData{transpose,N}) where {transpose,N}
    println(io, "SliceData{transpose=$transpose}:")
    for i = 1:N
        println(io, "  ", axisname(sd.axs[i]), ": ", value(sd.signals[i]))
    end
end
axisname(ax::Axis) = axisnames(ax)[1]

"""
    SliceData{transpose::Bool}(signals::NTuple{N,Signal{Int}}, axes::NTuple{N,Axes})

Specify slice information for a (possibly) multidimensional
image. `signals` hold the currently-selected slices for the selected
`axes`, all of which are effectively "orthogonal" to the plane in the
viewer.
"""
SliceData{transpose}(signals::NTuple{N,Signal{Int}}, axs::NTuple{N,Axis}) where {transpose,N} =
    SliceData{transpose,N,typeof(axs)}(signals, axs)
SliceData{transpose}() where {transpose} = SliceData{transpose}((), ())

Base.isempty(sd::SliceData{transpose,N}) where {transpose,N} = N == 0
Base.length(sd::SliceData{transpose,N}) where {transpose,N} = N

"""
    roi(A) -> zr::Signal(ZoomRegion), slicedata::SliceData
    roi(A, dims=(1,2)) -> zr::Signal(ZoomRegion), slicedata::SliceData
    roi(A, (:namey, :namex)) -> zr::Signal(ZoomRegion), slicedata::SliceData

Create the initial "region of interest" for viewing `A`. For
multidimensional objects, optionally select two dimensions (the first
two, by default) for slicing. The outputs `zr` and `slicedata`
describe the within-view and player-controlled axes, respectively.

See also: [`slice2d`](@ref).
"""
roi(A) = roi(A, (1,2))

roi(A, dims) = roi(axes(A), dims)
roi(A, axs::Tuple{Symbol,Symbol}) = roi(AxisArrays.axes(A), axs)

function roi(inds::Indices, dims::Dims{2})
    dims[1] != dims[2] || error("entries in dims must be distinct, got ", dims)
    zr = ZoomRegion(inds[[dims...]])
    sigs, axs = [], []
    for i = 1:length(inds)
        if !(i ∈ dims)
            ind = inds[i]
            push!(sigs, Signal(first(ind)))
            push!(axs, Axis{i}(ind))
        end
    end
    Signal(zr), SliceData{dims[2] < dims[1]}((sigs...,), (axs...,))
end

function roi(axs::NTuple{N,Axis}, axnames::Tuple{Symbol,Symbol}) where N
    axnames[1] != axnames[2] || error("entries in axnames must be distinct, got ", axnames)
    names = axisnames(axs...)
    dims = indexin([axnames...], [names...])
    inds = map(v->axes(v, 1), axisvalues(axs...))
    zr = ZoomRegion(inds[[dims...]])
    sigs, axs = [], []
    for (i, n) in enumerate(names)
        if !(n ∈ axnames)
            ind = inds[i]
            push!(sigs, Signal(first(ind)))
            push!(axs, Axis{n}(ind))
        end
    end
    Signal(zr), SliceData{dims[2] < dims[1]}((sigs...,), (axs...,))
end

"""
    slice2d(A, zr, sd) -> A2

Create a two-dimensional slice `A2` using the current ZoomRegion `zr`
and SliceData `sd`.
"""
function slice2d(img, zr::ZoomRegion, sd::SliceData{false})
    slice2d(img, makeroi(zr, false), makeslices(sd)...)
end
function slice2d(img, zr::ZoomRegion, sd::SliceData{true})
    transposedview(slice2d(img, makeroi(zr, true), makeslices(sd)...))
end

function makeroi(zr::ZoomRegion, transpose::Bool)
    rngs = map(UnitRange{Int}, (zr.currentview.y, zr.currentview.x))
    transpose ? (rngs[2], rngs[1]) : rngs
end

makeslices(sd::SliceData) = makeslices(sd.axs, sd.signals)
makeslices(axs::NTuple{N,Axis}, sigs::NTuple{N,Signal}) where {N} =
    map((ax,s) -> ax(value(s)), axs, sigs)

function slice2d(img::AbstractArray, roi, slices::Axis...)
    inds = sliceinds(img, roi, slices...)
    view(img, inds...)
end
function slice2d(img, roi, slices::Axis...)
    inds = sliceinds(img, roi, slices...)
    img[inds...]
end

# Infer whether we're using positional or named axes. This allows us
# to use positional labels even with AxisArrays.
abstract type TagType end
struct Positional <: TagType end
struct Named <: TagType end

TagType(::Type{Axis{name}}) where {name} = isa(name, Integer) ? Positional() : Named()
TagType(::Type{Axis{name,T}}) where {name,T} = TagType(Axis{name})
TagType(ax::Axis) = TagType(typeof(ax))

@inline tagtype(axs::Axis...) = check_same(map(TagType, axs)...)
@inline check_same(tt1::T, tt2::T, tts...) where {T<:TagType} = check_same(tt2, tts...)
@inline check_same(tt1::T, tts...) where {T<:TagType} = tt1
check_same(tt1::T) where {T<:TagType} = tt1
check_same() = Positional()

@noinline check_same(tt1::TagType, tt2::TagType, tts...) = error("must use either positional or named")

"""
    inds = sliceinds(img, roi, slices...)
    inds = sliceinds(axs, roi, slices...)

Return an indices-tuple `inds` that selects the region-of-interest `roi` at particular `slices`
along orthogonal coordinates. Use `AxisArrays.Axis` to indicate out-of-order axes, either
by dimension number (e.g., `Axis{3}`) or by name (`Axis{:z}`).
If you're using dimension numbers, `slices` must list them in increasing order.

# Examples:
```julia
julia> using AxisArrays: Axis

# Mimic a 4d array
julia> axs = (1:1080, 1:1920, 1:20, 1:1000)
(1:1080, 1:1920, 1:20, 1:1000)

julia> ImageView.sliceinds(axs, (6:10, 1:5), Axis{3}(7), Axis{4}(15))
(6:10, 1:5, 7, 15)

julia> ImageView.sliceinds(axs, (6:10, 1:5), Axis{2}(7), Axis{4}(15))
(6:10, 7, 1:5, 15)

# Mimic a 4d AxisArray
julia> axs = (Axis{:y}(1:1080), Axis{:x}(1:1920), Axis{:z}(1:20), Axis{:t}(1:1000))
(Axis{:y,UnitRange{Int64}}(1:1080), Axis{:x,UnitRange{Int64}}(1:1920), Axis{:z,UnitRange{Int64}}(1:20), Axis{:t,UnitRange{Int64}}(1:1000))

julia> ImageView.sliceinds(axs, (6:10, 1:5), Axis{:t}(15), Axis{:z}(7))   # out-of-order OK
(6:10, 1:5, 7, 15)

julia> ImageView.sliceinds(axs, (6:10, 1:5), Axis{:y}(15), Axis{:z}(7))
(15, 6:10, 7, 1:5)
```
"""
sliceinds(img, zoomranges, slices...) =
    sliceinds_t(tagtype(slices...), img, zoomranges, slices...)
sliceinds_t(::Positional, img, zoomranges, slices...) =
    sliceinds(axes(img), zoomranges, slices...)
sliceinds_t(::Named, img, zoomranges, slices...) =
    sliceinds(AxisArrays.axes(img), zoomranges, slices...)

# Positional
# Here we insist the axes are supplied in increasing order
sliceinds(axs::Tuple{}, zoomranges::Tuple{}) = ()
@inline sliceinds(inds::Indices, zoomranges, slices...) =
    _sliceinds((false,), inds, zoomranges, slices...)  # seed so that d-matching works below

@inline _sliceinds(out, inds::NTuple{N,Any}, ::Tuple{}) where {N} = tail(out)::NTuple{N,Any}
@inline _sliceinds(out, inds, zoomranges) =
    _sliceinds((out..., zoomranges[1]), inds, tail(zoomranges))
@inline _sliceinds(out::NTuple{d,Any}, inds, zoomranges, slice1::Axis{d}, slices...) where {d} =
    _sliceinds((out..., slice1.val), inds, zoomranges, slices...)
@inline _sliceinds(out, inds, zoomranges, slice1::Axis{d}, slices...) where {d} =
    _sliceinds((out..., zoomranges[1]), inds, tail(zoomranges), slice1, slices...)

# Named
# Here we allow any ordering of axes
@inline function sliceinds(axs::NTuple{N,Axis}, zoomranges, slices...) where N
    ind, newzoomranges, newslices = pickaxis((), axs[1], zoomranges, slices...)
    # ind, newzoomranges, newslices = pickaxis(axs[1], zoomranges, slices...)
    (ind, sliceinds(tail(axs), newzoomranges, newslices...)...)
end

@inline pickaxis(out, ax::Axis, zoomranges, slice1::Axis, slices...) =
    _pickaxis(out, AxisArrays.samesym(ax, slice1), ax, zoomranges, slice1, slices...)
@inline _pickaxis(out, ::Val{true}, ax, zoomranges, slice1, slices...) =
    slice1.val, zoomranges, (out..., slices...)
@inline _pickaxis(out, ::Val{false}, ax, zoomranges, slice1, slices...) =
    pickaxis((out..., slice1), ax, zoomranges, slices...)
@inline pickaxis(out, ax::Axis, zoomranges) =
    zoomranges[1], tail(zoomranges), out

# @generated function pickaxis(ax::Axis{name}, zoomranges, slices...) where name
#     idx = findfirst(x->axisnames(x)[1] == name, slices)
#     if idx == 0
#         return quote
#             zoomranges[1], tail(zoomranges), slices
#         end
#     end
#     newslices_exprs = [:(slices[$i]) for i in setdiff(1:length(slices), idx)]
#     quote
#         slices[$idx].val, zoomranges, ($(newslices_exprs...),)
#     end
# end

# For type stability we don't want to use permuteddimsview
transposedview(A::AbstractMatrix) =
    PermutedDimsArray{eltype(A),2,(2,1),(2,1),typeof(A)}(A)

function transposedview(A::AxisArray{T,2}) where T
    axs = AxisArrays.axes(A)
    AxisArray(transposedview(A.data), (axs[2], axs[1]))
end


text:
"""
    SummarizedExperiment(comm::CommunityProfile)

Converts a CommunityProfile to a SummarizedExperiment object.

# Arguments
- `comm::CommunityProfile`: the CommunityProfile to be converted.
"""
function SummarizedExperiment(comm::CommunityProfile)

    counts = abundances(comm)
    assays = OrderedDict{String, AbstractArray}("comm" => counts)

    rowdata = convert_rowdata(comm)

    coldata = DataFrame()
    symlist = [i for i in keys(metadata(comm)[1])]

    for sym in symlist

        coldata[!, sym] = [i[sym] for i in metadata(comm)]

    end

    rename!(coldata, Dict(:sample => "name"))

    # new(assays, rowdata, coldata, Dict{String,Any}())

    SummarizedExperiment(assays, rowdata, coldata)

end

"""
    CommunityProfile(se::SummarizedExperiment)

Converts a SummarizedExperiment object to a CommunityProfile.

# Arguments
- `se::SummarizedExperiment`: the SummarizedExperiment object to be converted.
"""
function CommunityProfile(se::SummarizedExperiment)

    samps = MicrobiomeSample.(se.coldata.name)

    taxa = Vector{Taxon}(undef, size(se.rowdata, 1))
    tax_ranks = ["strain", "subspecies", "species", "genus", "phamily", "order", "class", "phylum", "kingdom"]

    for (idx, row) in enumerate(eachrow(se.rowdata))

        not_missing_taxa = names(row)[collect(.~ ismissing.(values(row)))]

        for taxon in tax_ranks

            if sum(taxon .== not_missing_taxa) > 0

                taxa[idx] = Taxon(row[Symbol(taxon)], Symbol(taxon))

            end

        end

    end

    mat = assay(se, 1)

    comm = CommunityProfile(mat, taxa, samps)

    for key in Symbol.(names(se.coldata)[2:end])

        for (idx, row) in enumerate(eachrow(se.coldata))

            set!(samples(comm)[idx], key, row[key])

        end

    end

    comm

end

function convert_rowdata(comm::CommunityProfile{<:Real, Taxon, MicrobiomeSample})

    rowdata = DataFrame(
        name = featurenames(comm)
    )

    taxa = unique([feature.rank for feature in features(comm)])

    for taxon in taxa

        rowdata[!, taxon] = Vector{Union{Missing, String}}(missing, size(rowdata, 1))

    end

    for (idx, feature) in enumerate(features(comm))

        rowdata[idx, feature.rank] = feature.name

    end

    rowdata

end

function convert_rowdata(comm::CommunityProfile{<:Real, GeneFunction, MicrobiomeSample})

    rowdata = DataFrame(
        name = featurenames(comm)
    )

    taxa = unique([feature.taxon.rank for feature in features(comm)])

    for taxon in taxa

        rowdata[!, taxon] = Vector{Union{Missing, String}}(missing, size(rowdata, 1))

    end

    for (idx, feature) in enumerate(features(comm))

        rowdata[idx, feature.taxon.rank] = feature.taxon.name

    end

    rowdata

end

using Microbiome
# generate array with sample data
samps = MicrobiomeSample.(["s$i" for i in 1:10]);

# generate array with feature data
taxa = [[Taxon("s$i", :species) for i in 1:10]; [Taxon("g$i", :genus) for i in 1:10]];

# generate matrix with random entries
mat = rand(20, 10);

# create cp instance
comm = CommunityProfile(mat, taxa, samps)

# add some metadata about origin
for (sample, value) in zip(1:length(samples(comm)), ["origin$i" for i in 1:length(samples(comm))])

    set!(samples(comm)[sample], :origin, value)

end

gene_functions = GeneFunction.(["f$i" for i in 1:20], taxa)

comm = CommunityProfile(mat, gene_functions, samps)
se = SummarizedExperiment(comm)
comm


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

using Markdown
using InteractiveUtils

# ╔═╡ 60ccea97-08b4-417e-81e0-5110ce532d51
using PlutoUI

# ╔═╡ c9d6131c-d577-4523-9f1e-6f4543dd5f44
using Measures

# ╔═╡ 5115e0e4-0e37-4787-9b41-2a18605dc457
using LaTeXStrings

# ╔═╡ 504ac24f-5e3b-48fc-9ea7-b684274022f8
md"""
# Basic Julia types for data visualization
"""

# ╔═╡ 934c4b0c-913b-44f0-a880-50141b853165
md"""
## Functions
"""

# ╔═╡ b5c95cbd-1fe8-4796-92d0-3b4ae309c8d9
function addition(x, y)
    x + y
end

# ╔═╡ 4df5192a-c576-480f-ab1f-bbb73273a8cf
subtraction(x, y) = x - y

# ╔═╡ c27f7c14-c6bb-43f6-b231-37d66fd9dd81
md"""
### Higher-order and anonymous functions
"""

# ╔═╡ 2ed06105-6ca5-4cdb-b176-b4eedf24db43
sum(abs, [-1, 1])

# ╔═╡ 48dafa99-d4ff-4ec3-b58d-5a06ae49db26
sum(x -> abs(x), [-1, 1])

# ╔═╡ 4e4d3536-2360-4790-9962-d761ea565af6
sum([-1, 1]) do x
    abs(x)
end

# ╔═╡ ef345764-ffb7-4e59-a99a-30cbe3f383d6
md"""
## Julia types
"""

# ╔═╡ fa6157ea-49e0-4b48-b83c-fcfa334cc3a6
typeof("Hello World")

# ╔═╡ 4b5de8d6-c778-42e3-9a91-e63ff235913a
typeof(2)

# ╔═╡ c1ee4f85-b9f0-4ab3-acbc-50573c02a459
numbers = 1:5

# ╔═╡ dabf264b-93d0-492f-be5a-47cc3efb9b0a
Dump(numbers)

# ╔═╡ 8fdef782-3d8f-4560-9054-371bf49a705d
numbers.start

# ╔═╡ b99c8424-1da3-4fde-8d9a-f698c1511d03
numbers.stop

# ╔═╡ 49d1b1c3-45e4-41ef-91ed-562ee220ac5f
md"""
### Multiple dispatch
"""

# ╔═╡ 9699f41d-6f92-4528-9e38-a94c8ee72f25
concatenate(a::String, b::String) = a * b

# ╔═╡ 47b84567-81ad-4ddd-aebb-5cbaab2a1543
concatenate(a::Int, b::Int) = parse(Int, string(a) * string(b))

# ╔═╡ 67f9a128-d276-4455-827f-66cd58030517
concatenate("Hello", "World")

# ╔═╡ 0f77e7ed-1b7e-4004-b8cc-2061292ada5e
concatenate(1, 2)

# ╔═╡ a4b1caa7-bb57-49bc-bec0-6600a6d70951
methods(concatenate)

# ╔═╡ a1ac508b-9503-479f-8b16-3039c0d16c9b
concatenate(string(1), string(2))

# ╔═╡ 3a0be158-bd80-4837-a040-dc11470e3c58
convert(Float64, 1)

# ╔═╡ 9f915932-d7f3-47a9-a85c-2182531e8c46
Float64(1)

# ╔═╡ 0abeb1d5-1413-4dd4-8f20-9d95b866fa63
string(1)

# ╔═╡ 5634427a-e14d-4c37-86a4-6da6ca656c12
parse(Int, "1")

# ╔═╡ 72fc678c-06c7-43d7-9716-b91a1425aa89
md"""
### Representing numerical values
"""

# ╔═╡ 57757638-993e-4fbe-b0fc-faf0cc34e995
md"""
`Int`
"""

# ╔═╡ ec4f8362-39ca-4cf6-af8b-857e152d82b7
-1

# ╔═╡ 082a675a-d7c5-4c33-a5d7-1e83570f16fd
5

# ╔═╡ 6947633b-6896-4bc6-929c-50a68002a94d
10_000

# ╔═╡ 61371d87-acb9-448a-a875-f5745224c5df
md"`Float64`"

# ╔═╡ b5b0311e-9c16-4b9a-8091-43838b7feb8f
1.0

# ╔═╡ 588fd62f-509c-4b99-8e23-9b60e14d5602
-2e3

# ╔═╡ c2d8adc5-1ef2-4c46-bfb2-9ea30ca12e89
13.5E10

# ╔═╡ c54a18cd-eddc-4e20-bb7e-1b33e832b422
.1

# ╔═╡ c03b8c69-64e3-4223-8765-2818ea52f5a2
1.

# ╔═╡ 38e70ac5-fa8d-4ca8-9c71-c104a48b1e87
NaN

# ╔═╡ 05a088e1-60b6-47df-bf61-fddf548563e0
md"Units"

# ╔═╡ 1bcbfd07-375a-4fcb-913b-cd159fce54e7
10.5cm

# ╔═╡ 9f9573af-b65a-4246-aefa-64cc48ea9dcd
md"""
### Representing text
"""

# ╔═╡ ac0abe2d-6fe0-4a11-b212-0ca1ecc68ae2
"Single line String"

# ╔═╡ cc3a3304-c9b5-45e5-8209-9c94bfc611dd
"""
Multiline
String
"""

# ╔═╡ 0cc52da8-303d-45fd-9390-7b203cd15773
'a'

# ╔═╡ dfd31ce8-4c00-4ebe-a5f9-4d50b4c2e3cd
md"**Markdown** string"

# ╔═╡ 4f79eb5e-3cd6-42a0-b557-5a5821a9428b
html"HTML string"

# ╔═╡ c91a721d-1cb5-4839-8566-842c996e61a1
L"LaTeX string"

# ╔═╡ daf19907-5675-4430-886d-a69f3b71aed9
md"`Symbol`"

# ╔═╡ c71d8f39-4911-450d-bb97-8b93a61693a7
:var1

# ╔═╡ 3e72cf06-2247-4a71-bbaa-c1d420ed83e4
Symbol("1var")

# ╔═╡ ca7472e6-b08b-410b-a921-cccdd5fbfc31
md"""
### Collections
"""

# ╔═╡ 1d1b47c9-429c-41ec-9e09-38bddd2e3f7e
md"`Tuple`"

# ╔═╡ 8c9c54d7-403c-48a7-b63e-12c7fdd5659c
("x", 0)

# ╔═╡ 69631bb0-f279-4bf2-86d0-0f21bc6d3314
(1,)

# ╔═╡ 7e07074d-4aa3-459c-9c90-24718815d865
md"`Array`"

# ╔═╡ 86fab7a8-6d48-42dc-b963-df227162fb39
[1]

# ╔═╡ ff173169-77aa-4e6b-988d-0d1f7bd15a02
[1, 2]

# ╔═╡ 96e44dec-8ed7-4b87-b0de-c069d22e54f7
[1 2; 3 4]

# ╔═╡ 7bcf656c-e397-4544-a2c1-7c5a082a42d6
[ 1 2
  3 4 ]

# ╔═╡ bcd44701-7992-45fc-be4d-e4a13cfbcc63
zeros(4, 2)

# ╔═╡ 5a93658c-5807-46ef-8faf-108347a38bb4
zeros(10)

# ╔═╡ d66723b9-7db7-4943-9a34-b55a930fa207
collection = collect(1:5)

# ╔═╡ db71a959-3ef0-4da5-bd27-deaa4c47837b
md"Indexing"

# ╔═╡ 666df65b-3679-42ec-877e-4947b35f97a5
collection[1]

# ╔═╡ 3348ed56-0d26-407e-90b1-ba1a82e9f3d3
md"Array comprehension"

# ╔═╡ bff3d2f1-8255-4fc6-ac6f-43285d5d6989
[sqrt(x) for x in 1:5]

# ╔═╡ 5ea04cfd-7d8c-450e-9335-65f9a13887ae
[sqrt(x) for x in 1:10 if x % 2 != 0]

# ╔═╡ d8f5d3ae-5175-43e4-8279-77024d69856f
md"Broadcasting"

# ╔═╡ e60f8bdc-b825-4a69-aec0-89ea01e590fe
sqrt.(collection) .+ 1

# ╔═╡ 86deb0c1-f291-4bd9-861d-4355fe2dff81
md"""
---
"""

# ╔═╡ 29a2068e-47f7-4327-b52c-457f521a970f
PlutoUI.TableOfContents()

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
LaTeXStrings = "~1.2.1"
Measures = "~0.3.1"
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Measures]]
git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╟─504ac24f-5e3b-48fc-9ea7-b684274022f8
# ╟─934c4b0c-913b-44f0-a880-50141b853165
# ╠═b5c95cbd-1fe8-4796-92d0-3b4ae309c8d9
# ╠═4df5192a-c576-480f-ab1f-bbb73273a8cf
# ╟─c27f7c14-c6bb-43f6-b231-37d66fd9dd81
# ╠═2ed06105-6ca5-4cdb-b176-b4eedf24db43
# ╠═48dafa99-d4ff-4ec3-b58d-5a06ae49db26
# ╠═4e4d3536-2360-4790-9962-d761ea565af6
# ╟─ef345764-ffb7-4e59-a99a-30cbe3f383d6
# ╠═fa6157ea-49e0-4b48-b83c-fcfa334cc3a6
# ╠═4b5de8d6-c778-42e3-9a91-e63ff235913a
# ╠═60ccea97-08b4-417e-81e0-5110ce532d51
# ╠═c1ee4f85-b9f0-4ab3-acbc-50573c02a459
# ╠═dabf264b-93d0-492f-be5a-47cc3efb9b0a
# ╠═8fdef782-3d8f-4560-9054-371bf49a705d
# ╠═b99c8424-1da3-4fde-8d9a-f698c1511d03
# ╟─49d1b1c3-45e4-41ef-91ed-562ee220ac5f
# ╠═9699f41d-6f92-4528-9e38-a94c8ee72f25
# ╠═67f9a128-d276-4455-827f-66cd58030517
# ╠═47b84567-81ad-4ddd-aebb-5cbaab2a1543
# ╠═0f77e7ed-1b7e-4004-b8cc-2061292ada5e
# ╠═a4b1caa7-bb57-49bc-bec0-6600a6d70951
# ╠═a1ac508b-9503-479f-8b16-3039c0d16c9b
# ╠═3a0be158-bd80-4837-a040-dc11470e3c58
# ╠═9f915932-d7f3-47a9-a85c-2182531e8c46
# ╠═0abeb1d5-1413-4dd4-8f20-9d95b866fa63
# ╠═5634427a-e14d-4c37-86a4-6da6ca656c12
# ╟─72fc678c-06c7-43d7-9716-b91a1425aa89
# ╟─57757638-993e-4fbe-b0fc-faf0cc34e995
# ╠═ec4f8362-39ca-4cf6-af8b-857e152d82b7
# ╠═082a675a-d7c5-4c33-a5d7-1e83570f16fd
# ╠═6947633b-6896-4bc6-929c-50a68002a94d
# ╟─61371d87-acb9-448a-a875-f5745224c5df
# ╠═b5b0311e-9c16-4b9a-8091-43838b7feb8f
# ╠═588fd62f-509c-4b99-8e23-9b60e14d5602
# ╠═c2d8adc5-1ef2-4c46-bfb2-9ea30ca12e89
# ╠═c54a18cd-eddc-4e20-bb7e-1b33e832b422
# ╠═c03b8c69-64e3-4223-8765-2818ea52f5a2
# ╠═38e70ac5-fa8d-4ca8-9c71-c104a48b1e87
# ╟─05a088e1-60b6-47df-bf61-fddf548563e0
# ╠═c9d6131c-d577-4523-9f1e-6f4543dd5f44
# ╠═1bcbfd07-375a-4fcb-913b-cd159fce54e7
# ╟─9f9573af-b65a-4246-aefa-64cc48ea9dcd
# ╠═ac0abe2d-6fe0-4a11-b212-0ca1ecc68ae2
# ╠═cc3a3304-c9b5-45e5-8209-9c94bfc611dd
# ╠═0cc52da8-303d-45fd-9390-7b203cd15773
# ╠═5115e0e4-0e37-4787-9b41-2a18605dc457
# ╠═dfd31ce8-4c00-4ebe-a5f9-4d50b4c2e3cd
# ╠═4f79eb5e-3cd6-42a0-b557-5a5821a9428b
# ╠═c91a721d-1cb5-4839-8566-842c996e61a1
# ╟─daf19907-5675-4430-886d-a69f3b71aed9
# ╠═c71d8f39-4911-450d-bb97-8b93a61693a7
# ╠═3e72cf06-2247-4a71-bbaa-c1d420ed83e4
# ╟─ca7472e6-b08b-410b-a921-cccdd5fbfc31
# ╟─1d1b47c9-429c-41ec-9e09-38bddd2e3f7e
# ╠═8c9c54d7-403c-48a7-b63e-12c7fdd5659c
# ╠═69631bb0-f279-4bf2-86d0-0f21bc6d3314
# ╟─7e07074d-4aa3-459c-9c90-24718815d865
# ╠═86fab7a8-6d48-42dc-b963-df227162fb39
# ╠═ff173169-77aa-4e6b-988d-0d1f7bd15a02
# ╠═96e44dec-8ed7-4b87-b0de-c069d22e54f7
# ╠═7bcf656c-e397-4544-a2c1-7c5a082a42d6
# ╠═bcd44701-7992-45fc-be4d-e4a13cfbcc63
# ╠═5a93658c-5807-46ef-8faf-108347a38bb4
# ╠═d66723b9-7db7-4943-9a34-b55a930fa207
# ╟─db71a959-3ef0-4da5-bd27-deaa4c47837b
# ╠═666df65b-3679-42ec-877e-4947b35f97a5
# ╟─3348ed56-0d26-407e-90b1-ba1a82e9f3d3
# ╠═bff3d2f1-8255-4fc6-ac6f-43285d5d6989
# ╠═5ea04cfd-7d8c-450e-9335-65f9a13887ae
# ╟─d8f5d3ae-5175-43e4-8279-77024d69856f
# ╠═e60f8bdc-b825-4a69-aec0-89ea01e590fe
# ╟─86deb0c1-f291-4bd9-861d-4355fe2dff81
# ╠═29a2068e-47f7-4327-b52c-457f521a970f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
export LangMatlab

# Data structure for the language.
"""
    LangMatlab()

Code generation for the Matlab language.
"""
struct LangMatlab end

# Language specific operations
comment(::LangMatlab, s) = "% $s"

slotname(::LangMatlab, i) = "memslots{$i}"

function assign_coeff(::LangMatlab, v, i)
    if imag(v) == 0
        ("coeff$i", "coeff$i = $(v)")
    else
        ("coeff$i", "coeff$i = $(real(v)) + 1i*$(imag(v))")
    end
end

# Code generation.
function function_definition(
    lang::LangMatlab,
    graph,
    T,
    funname,
    precomputed_nodes,
)
    code = init_code(lang)
    input_variables = join(precomputed_nodes, ", ")
    push_code!(
        code,
        "function output = $funname($input_variables)",
        ind_lvl = 0,
    )
    return code
end

function function_init(lang::LangMatlab, T, mem, graph, precomputed_nodes)
    code = init_code(lang)
    push_code!(code, "n = size(A,1);")
    push_code!(code, "I = eye(n,n);")
    return code
end

function init_mem(lang::LangMatlab, max_nof_nodes, precomputed_nodes)
    mem = CodeMem(max_nof_nodes, i -> slotname(lang, i))
    return mem
end

function function_end(lang::LangMatlab, graph, mem)
    code = init_code(lang)
    push_code!(code, "output = $(graph.outputs[end]);")
    push_code!(code, "end", ind_lvl = 0)
    return code
end

function execute_operation!(lang::LangMatlab, T, graph, node, dealloc_list, mem)
    op = graph.operations[node]
    parent1 = graph.parents[node][1]
    parent2 = graph.parents[node][2]

    code = init_code(lang)
    push_comment!(code, "Computing $node with operation: $op")
    # Needs to be updated
    if op == :mult
        push_code!(code, "$node = $parent1 * $parent2;")
    elseif op == :ldiv
        push_code!(code, "$node = $parent1 \\ $parent2;")
    elseif op == :lincomb
        (coeff1, coeff1_code) = assign_coeff(lang, graph.coeffs[node][1], 1)
        push_code!(code, "$coeff1_code;")
        (coeff2, coeff2_code) = assign_coeff(lang, graph.coeffs[node][2], 2)
        push_code!(code, "$coeff2_code;")
        push_code!(code, "$node = $coeff1*$parent1 + $coeff2*$parent2;")
    end
    return (code, "$node")
end


text:
"""
    sum_ap(first_term, diff, num_terms)

Finds sum of a arithmetic progression series

# Input parameters

- first_term : first term of the series
- diff       : common difference between consecutive terms
- num_terms  : number of terms in the series till which we count sum

# Example

```julia
sum_ap(1, 1, 10)    # returns 55.0
sum_ap(1, 10, 100)  # returns 49600.0
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function sum_ap(first_term, diff, num_terms)
    #formula for sum of the ap series
    sum = (num_terms / 2) * (2 * first_term + (num_terms - 1) * diff)
    return sum
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia