diff --git a/.gitignore b/.gitignore
index b6dc83d..adafdf5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,3 +6,4 @@ experiments
 *.mat
 *.csv
 *.pdf
+Manifest.toml
diff --git a/Project.toml b/Project.toml
new file mode 100644
index 0000000..cf63e32
--- /dev/null
+++ b/Project.toml
@@ -0,0 +1,15 @@
+name = "SumProductNetworks"
+uuid = "5f6e642e-680c-5a9c-a175-7c23ed4da89e"
+version = "0.1.4"
+
+[deps]
+AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+
+[compat]
+julia = ">= 1.0"
diff --git a/src/SumProductNetworks.jl b/src/SumProductNetworks.jl
index 2119def..a94e30c 100644
--- a/src/SumProductNetworks.jl
+++ b/src/SumProductNetworks.jl
@@ -20,17 +20,22 @@ using SparseArrays
 using Random
 using Printf
 
+# include custom distributions
+include("distributions.jl")
+
 # include general implementations
 include("nodes.jl")
 include("nodeFunctions.jl")
 include("networkFunctions.jl")
-
-# include approach specific implementations
-include("bmiTest.jl")
-include("structureUtilities.jl")
+include("regiongraphs.jl")
 
 # include utilities
 include("utilityFunctions.jl")
 include("io.jl")
 
+# include approach specific implementations
+include("bmiTest.jl")
+include("structureUtilities.jl")
+include("ratspn.jl")
+
 end # module
diff --git a/src/distributions.jl b/src/distributions.jl
new file mode 100644
index 0000000..d055a86
--- /dev/null
+++ b/src/distributions.jl
@@ -0,0 +1,32 @@
+export NormalInverseGamma
+
+## Gaussian with Normal Inverse Gamma prior ##
+struct NormalInverseGamma{T<:Real} <: Distribution{Multivariate,Continuous}
+    μ::T
+    ν::T
+    a::T
+    b::T
+end
+NormalInverseGamma() = NormalInverseGamma(0.0, 1.0, 1.0, 1.0)
+
+Distributions.length(d::NormalInverseGamma) = 2
+
+function Distributions.logpdf(d::NormalInverseGamma, μ::T, σ²::T) where {T}
+    lp = _invgammalogpdf(d.a, d.b, σ²)
+    lp += normlogpdf(d.μ, sqrt(σ²) / d.ν, μ)
+    return l
+end
+
+@inline _invgammalogpdf(a::T, b::T, x::T) where {T<:Real} = a*log(b)-lgamma(a)-(a+one(T))*log(x)-b/x
+
+function Distributions._rand!(rng::AbstractRNG, d::NormalInverseGamma, out::AbstractArray{T,2}) where {T}
+    for p in eachslice(out, dims=[2])
+        @inbounds begin
+            p[2] = rand(InverseGamma(d.a, d.b))
+            p[1] = rand(Normal(d.μ, sqrt(p[2] / d.ν )))
+        end
+    end
+    return out
+end
+
+@inline Distributions.params(d::NormalInverseGamma) = (d.μ, d.ν, d.a, d.b)
diff --git a/src/ratspn.jl b/src/ratspn.jl
new file mode 100644
index 0000000..7f74906
--- /dev/null
+++ b/src/ratspn.jl
@@ -0,0 +1,103 @@
+export ratspn
+#export templateLeaves, templatePartition, templateRegion
+
+function templateLeaves(likelihoods::AbstractVector{<:Distribution},
+                        priors::AbstractVector{<:Distribution},
+                        N::Int, D::Int, K::Int)
+    scopeVec = zeros(Bool, D)
+    parameters = map(prior -> rand(prior,K), priors)
+
+    return FactorizedDistributionGraphNode(
+                                    gensym("fact"),
+                                    scopeVec,
+                                    likelihoods,
+                                    priors,
+                                    parameters
+                                   )
+end
+
+function templatePartition(likelihoods::AbstractVector,
+                           N::Int, D::Int,
+                           K_sum::Int, K_prod::Int,
+                           J::Int, K::Int,
+                           depth::Int, maxdepth::Int)
+
+    children = if depth == maxdepth
+        map(k -> templateLeaves(alpha_leaf_prior, priors_leaf, likelihoods, sstats, N, D, K), 1:K_prod)
+    else
+        map(k -> templateRegion(alpha_region_prior, alpha_partition_prior, alpha_leaf_prior,
+                                priors_leaf, likelihoods, sstats, N, D, K_sum, K_prod, J, K, depth+1, maxdepth), 1:K_prod)
+    end
+
+    K_ = mapreduce(child -> length(child), *, children)
+    scopeVec = zeros(Bool, D)
+    obsVec = zeros(Bool, N, K_)
+
+    return PartitionGraphNode(
+                              gensym("partition"),
+                              scopeVec,
+                              obsVec,
+                              prior,
+                              children
+                             )
+end
+
+function buildRegion(llhs::AbstractVector, priors::AbstractVector,
+                        N::Int, D::Int,
+                        K_sum::Int, K_prod::Int,
+                        J::Int, K::Int,
+                        depth::Int, maxdepth::Int; root = false)
+
+    K_ = root ? 1 : K_sum
+    children = if depth == maxdepth
+        map(k -> buildLeaves(likelihoods, priors, N, D, K), 1:J)
+    else
+        map(k -> buildPartition(likelihoods, priors, N, D, K_sum, K_prod, J, K, depth+1, maxdepth), 1:J)
+    end
+
+    Ch = sum(length.(children))
+    scopeVec = zeros(Bool, D)
+    logweights = rand(Dirichlet(Ch, 1.0), K_)
+    active = zeros(Bool, size(logweights)...)
+    @assert size(logweights) == (Ch, K_)
+
+    return RegionGraphNode( gensym("region"), scopeVec, logweights, children )
+end
+
+
+function ratspn(x::AbstractMatrix{T};
+               Ksplits::Int=2,
+               Kparts::Int=2,
+               Ksums::Int=2,
+               Kdists::Int=5,
+               maxdepth::Int=2
+              ) where {T<:Real}
+
+    N,D = size(x)
+    isdiscrete = map(d -> all(isinteger, x[:,d]), 1:D)
+
+    K = map(d -> isdiscrete[d] ? length(unique(x[:,d])) : Inf, 1:D)
+
+    llhs = map(d -> isdiscrete[d] ? Categorical(K[d]) : Normal(), 1:D)
+    priors = map(d -> isdiscrete[d] ? Dirichlet(K[d], 1.0) : NormalInverseGamma(), 1:D)
+
+    return ratspn(N,D, llhs, priors, Ksplits, Kparts, Ksums, Kdists, maxdepth)
+end
+
+function ratspn(N::Int, D::Int,
+                llhs::AbstractVector{<:Distribution},
+                priors::AbstractVector{<:Distribution},
+                Ksplits::Int, # Number of partitions under each region
+                Kparts::Int, # Number of sub-regions under each partition
+                Ksums::Int, # Number of sum nodes per region
+                Kdists::Int, # Number of distibutions per terminal region
+                maxdepth::Int # Maximum number of consecutive region-partition pairs
+                )
+
+    # sanity checks
+    @assert length(llhs) == D == length(priors)
+
+    
+
+    #return templateRegion()
+end
diff --git a/src/regiongraphs.jl b/src/regiongraphs.jl
new file mode 100644
index 0000000..3f6a9f3
--- /dev/null
+++ b/src/regiongraphs.jl
@@ -0,0 +1,246 @@
+export draw
+export FactorizedAtomicRegion
+
+abstract type AbstractRegionGraphNode end
+
+function setscope!(n::AbstractRegionGraphNode, dims::Vector{Int})
+    fill!(n.scopeVecs, false)
+    @inbounds n.scopeVecs[dims] .= true
+    empty!(n.scope)
+    append!(n.scope, dims)
+end
+nscope(d::AbstractRegionGraphNode) = sum(d.scopeVecs)
+hasscope(d::AbstractRegionGraphNode) = any(d.scopeVecs)
+scope(d::AbstractRegionGraphNode) = d.scope
+
+function setscope!(n::AbstractRegionGraphNode, i::Int, s::Bool)
+    @inbounds n.scopeVecs[i] = s
+    push!(n.scope, s)
+end
+
+id(d::AbstractRegionGraphNode) = d.id
+
+"""
+    topoligicalorder(root::AbstractRegionGraphNode)
+
+Return topological order of the nodes in the SPN.
+"""
+function topoligicalorder(root::AbstractRegionGraphNode)
+  visitedNodes = Vector{AbstractRegionGraphNode}()
+  visitNode!(root, visitedNodes)
+  return visitedNodes
+end
+
+function visitNode!(node::AbstractRegionGraphNode, visitedNodes)
+  # check if we have already visited this node
+  if !(node in visitedNodes)
+
+    # visit node
+    if haschildren(node)
+        for child in children(node)
+            visitNode!(child, visitedNodes)
+        end
+    end
+    push!(visitedNodes, node)
+  end
+end
+
+"""
+    Region node.
+
+Parameters:
+
+* ids::Symbol                                           Id
+* scopeVecs::Vector{Bool}                               Active dimensions (D)
+* logweights::Matrix{<:Real}                            Log weights of sum nodes (Ch x K)
+* prior::Dirichlet                                      Prior for sum nodes
+* children::Vector{AbstractRegionGraphNode}             Children of region
+
+"""
+struct SumRegion{T<:Real,V<:AbstractRegionGraphNode} <: AbstractRegionGraphNode
+    id::Symbol
+    scopeVecs::Vector{Bool}
+    logweights::Matrix{T}
+    children::Vector{V}
+end
+
+haschildren(n::SumRegion) = !isempty(n.children)
+children(n::SumRegion) = n.children
+
+function updatescope!(n::SumRegion)
+    updatescope!.(children(n))
+    setscope!(n, scope(first(n.children)))
+end
+
+# Multi-threaded LSE
+function _logsumexp(y::AbstractMatrix{T}, lw::AbstractMatrix{T}) where {T<:Real}
+    r = zeros(T,size(lw,2),size(y,2))
+    Threads.@threads for j in 1:size(y,2)
+        for k in 1:size(lw,2)
+            @inbounds begin
+                yi_max = y[1,j] + lw[1,k]
+                for i in 2:size(y,1)
+                    yi_max = max(y[i,j]+lw[i,k], yi_max)
+                end
+                s = zero(T)
+                for i in 1:size(y,1)
+                    s += exp(y[i,j]+lw[i,k] - yi_max)
+                end
+                r[k,j] = log(s) + yi_max
+            end
+        end
+    end
+    return transpose(r)
+end
+
+function _getchildlogpdf(child::AbstractRegionGraphNode, out::AxisArray{V}) where {V<:AbstractMatrix}
+    return out[id(child)]
+end
+
+function _getchildlogpdf(child::PartitionGraphNode, out::AxisArray{V}) where {V<:AbstractMatrix}
+    childids = map(id, children(child))
+    lp_ = out[childids]
+    return reduce(_cross_prod, lp_)
+end
+
+function logpdf(n::SumRegion, x::AbstractMatrix{T}) where {T<:Real}
+    out = zeros(T, size(x,1), lengh(n))
+    logdf!(n, x, out)
+    return out
+end
+
+function logpdf!(n::SumRegion{T,Tc},
+                 x::AbstractMatrix{T},
+                 out::AbstractMatrix{T}) where {T<:Real,Tc<:FactorizedAtomicRegion}
+
+    for k1 = 1:length()
+
+
+
+
+    lp_ = _getchildlogpdf.(children(d), Ref(out))
+    out[id(d)] = _logsumexp(reduce(hcat, lp_)', d.logweights)
+    return out
+end
+
+"""
+    logmllh!(n::RegionGraphNode, x::AbstractMatrix{T}, out::AxisArray{V}, y::AxisArray, d::Int)
+
+Log marginal likelihood of a region node in a region graph for dimension d. (in-place)
+"""
+function logmllh!(n::RegionGraphNode,
+                  x::AbstractMatrix{T},
+                  out::AxisArray{V},
+                  y::AxisArray,
+                  d::Int) where {T<:Real,V<:AbstractVector}
+    K = length(n)
+    lp_ = out[id.(children(n))]
+    lp_ = reduce(vcat, lp_) .* n.active
+    out[id(n)] = vec(sum(lp_, dims = 1))
+    return out
+end
+
+function _cross_prod(x1::AbstractMatrix{T}, x2::AbstractMatrix{T}) where {T<:Real}
+    nx, ny = size(x1,2), size(x2,2)
+    r = zeros(T, size(x1,1), nx*ny)
+    Threads.@threads for j in 1:nx*ny
+        @inbounds begin
+            j1, j2 = Base._ind2sub((nx,ny), j)
+            r[:,j] .= x1[:,j1] .+ x2[:,j2]
+        end
+    end
+    return r
+end
+
+"""
+    Factorized distribution of the scope.
+
+Parameters:
+
+* ids::Symbol                               Id
+* scopeVecs::Vector{Bool}                   Active dimensions (D)
+* likelihoods::Vector{Type{Distribution}}   Likelihood functions (D)
+* priors::Vector{Distribution}              Priors for each dimension (D)
+* parameters::Vector{Matrix}                Parameters for each k≦K for each dimension (D)
+
+"""
+struct FactorizedAtomicRegion <: AbstractRegionGraphNode
+    id::Symbol
+    K::Int
+    scopeVecs::Vector{Bool}
+    scope::Vector{Int}
+    likelihoods::Vector{Function}
+    priors::Vector{Distribution}
+    parameters::Vector{Matrix}
+end
+
+function FactorizedAtomicRegion(likelihoods::Vector{Function},
+                                         parameters::Vector{<:Matrix},
+                                         D::Int;
+                                         priors::Vector{Distribution}=Vector{Distribution}()
+                                        )
+    K = size(first(parameters), 2)
+    for p in parameters
+        @assert size(p,2) == K
+    end
+    svec = falses(D)
+    return FactorizedAtomicRegion(gensym(), K, Int[], svec, likelihoods, priors, parameters)
+end
+
+function draw(priors::Vector{Distribution}, K::Int)
+    return map(prior -> prior isa UnivariateDistribution ? reshape(rand(prior,K),1,:) : rand(prior,K), priors)
+end
+
+@inline haschildren(n::FactorizedAtomicRegion) = false
+@inline updatescope!(n::FactorizedAtomicRegion) = scope(d)
+@inline length(n::FactorizedAtomicRegion) = n.K
+
+function apply!(n::FactorizedAtomicRegion,
+                x::AbstractMatrix{T},
+                out::AbstractVector{V},
+                k::Int
+               ) where {T,V}
+    for d in scope(n)
+        @inbounds begin
+            θ = @view(n.parameters[d][:,k])
+            out[:] += n.likelihoods[d].(@view(x[:,d]), θ...)
+        end
+    end
+end
+
+function _logpdf!(n::FactorizedAtomicRegion,
+                  x::AbstractMatrix{T},
+                  out::AbstractMatrix{V}) where {T<:Real,V<:Real}
+    fill!(out, zero(V))
+    for k in 1:lenght(n)
+        apply!(n, x, @view(out[:,k]), k)
+    end
+    return out
+end
+
+function logpdf(n::FactorizedAtomicRegion, x::AbstractMatrix{T}) where {T<:Real}
+    N = size(x,1)
+    K = length(n)
+    lp_ = zeros(Float64, N, K)
+    return _logpdf!(n, x, lp_)
+end
+
+function _pdf!(n::FactorizedAtomicRegion,
+                  x::AbstractMatrix{T},
+                  out::AbstractMatrix{V}) where {T<:Real,V<:Real}
+    fill!(out, zero(V))
+    ds = scope(n)
+    for d in ds
+        for (θ, lp) in zip(eachslice(n.parameters[d], dims=2), eachslice(out, dims=2))
+            @inbounds lp[:] .*= exp.(n.likelihoods[d].(@view(x[:,d]), θ...))
+        end
+    end
+    return out
+end
+
+function pdf(n::FactorizedAtomicRegion, x::AbstractMatrix{T}) where {T<:Real}
+    N = size(x,1)
+    K = length(n)
+    lp_ = zeros(Float64, N, K)
+    return _pdf!(n, x, lp_)
+end
diff --git a/test/regiongraph_test.jl b/test/regiongraph_test.jl
new file mode 100644
index 0000000..6ed51f5
--- /dev/null
+++ b/test/regiongraph_test.jl
@@ -0,0 +1,61 @@
+using SumProductNetworks
+using StatsFuns
+using BenchmarkTools
+using Test, Random
+
+Random.seed!(1)
+
+# some test data
+x = hcat(randn(100,1), rand(1:10, 100,1), rand(Bool,100,1))
+
+catlogpdf(x::Int, p...) = log(p[x])
+catlogpdf(x::AbstractFloat, p...) = catlogpdf(Int(x), p...)
+bernlogpdf(x::Int, p) = x == 0 ? 1-p : p
+bernlogpdf(x::AbstractFloat, p) = bernlogpdf(Int(x), p)
+
+likelihoods = Function[(x, μ, σ) -> normlogpdf(μ, σ, x), catlogpdf, bernlogpdf]
+priors = Distribution[NormalInverseGamma(), Dirichlet(10, 1.0), Beta()]
+
+@testset "graph nodes" begin
+    @testset "atomic region" begin
+
+        K = 10
+        parameters = draw(priors, K)
+        node = FactorizedAtomicRegion( likelihoods, parameters, 3 )
+        @test node isa FactorizedAtomicRegion
+
+        lp = logpdf(node, x)
+        @test size(lp) == (100, K)
+
+        @test all(lp .== 0.0)
+
+        setscope!(node, [1])
+        lp = logpdf(node, x)
+        lp2 = zero(lp)
+        for (p,l) in zip(eachslice(parameters[1], dims=2), eachslice(lp2, dims=2))
+            l[:] = normlogpdf.(p[1], p[2], x[:,1])
+        end
+        @test all(lp .== lp2)
+
+        setscope!(node, [1,2])
+        @test scope(node) == [1,2]
+        lp = logpdf(node, x)
+        lp2 = zero(lp)
+        for (p,l) in zip(eachslice(parameters[1], dims=2), eachslice(lp2, dims=2))
+            l[:] = normlogpdf.(p[1], p[2], x[:,1])
+        end
+        for (p,l) in zip(eachslice(parameters[2], dims=2), eachslice(lp2, dims=2))
+            l[:] += catlogpdf.(x[:,2], p...)
+        end
+        @test all(lp .== lp2)
+
+    end
+    @testset "sum region" begin
+
+
+
+    end
+end
+
+# build a region graph
+#spn = ratspn(x)