diff --git a/Manifest.toml b/Manifest.toml
index 189abedd..be002f01 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -57,6 +57,11 @@ git-tree-sha1 = "9a11d428dcdc425072af4aea19ab1e8c3e01c032"
 uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 version = "1.3.0"
 
+[[DataAPI]]
+git-tree-sha1 = "8903f0219d3472543fc4b2f5ebaf675a07f817c0"
+uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+version = "1.0.1"
+
 [[DataStructures]]
 deps = ["InteractiveUtils", "OrderedCollections", "Random", "Serialization", "Test"]
 git-tree-sha1 = "ca971f03e146cf144a9e2f2ce59674f5bf0e8038"
@@ -89,15 +94,15 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
 [[Distributions]]
 deps = ["LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
-git-tree-sha1 = "4ae0ff1a3a556383f7b1b004f5fb964b26ac96c9"
+git-tree-sha1 = "baaf9e165ba8a2d11fb4fb3511782ee070ee3694"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.20.0"
+version = "0.21.1"
 
 [[FFTW]]
-deps = ["AbstractFFTs", "BinaryProvider", "Compat", "Conda", "Libdl", "LinearAlgebra", "Reexport", "Test"]
-git-tree-sha1 = "29cda58afbf62f35b1a094882ad6c745a47b2eaa"
+deps = ["AbstractFFTs", "BinaryProvider", "Conda", "Libdl", "LinearAlgebra", "Reexport", "Test"]
+git-tree-sha1 = "03f8776fbdae28c20c0d1d2ae4e090cd1dfcd247"
 uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-version = "0.2.4"
+version = "1.0.0"
 
 [[ForwardDiff]]
 deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "InteractiveUtils", "LinearAlgebra", "NaNMath", "Random", "SparseArrays", "SpecialFunctions", "StaticArrays", "Test"]
@@ -139,10 +144,9 @@ deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
 [[Missings]]
-deps = ["SparseArrays", "Test"]
-git-tree-sha1 = "f0719736664b4358aa9ec173077d4285775f8007"
+git-tree-sha1 = "29858ce6c8ae629cf2d733bffa329619a1c843d0"
 uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
-version = "0.4.1"
+version = "0.4.2"
 
 [[Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
@@ -167,9 +171,9 @@ version = "0.9.9"
 
 [[Parsers]]
 deps = ["Dates", "Test"]
-git-tree-sha1 = "db2b35dedab3c0e46dc15996d170af07a5ab91c9"
+git-tree-sha1 = "ef0af6c8601db18c282d092ccbd2f01f3f0cd70b"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "0.3.6"
+version = "0.3.7"
 
 [[Pkg]]
 deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
@@ -187,9 +191,9 @@ uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
 [[QuadGK]]
 deps = ["DataStructures", "LinearAlgebra"]
-git-tree-sha1 = "389fd27b958a33df5234772cc464b663b208273d"
+git-tree-sha1 = "438630b843c210b375b2a246329200c113acc61b"
 uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
-version = "2.0.4"
+version = "2.1.0"
 
 [[REPL]]
 deps = ["InteractiveUtils", "Markdown", "Sockets"]
@@ -251,10 +255,10 @@ deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[StatsBase]]
-deps = ["DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"]
-git-tree-sha1 = "8a0f4b09c7426478ab677245ab2b0b68552143c7"
+deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"]
+git-tree-sha1 = "c53e809e63fe5cf5de13632090bc3520649c9950"
 uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-version = "0.30.0"
+version = "0.32.0"
 
 [[StatsFuns]]
 deps = ["Rmath", "SpecialFunctions", "Test"]
diff --git a/NEWS.md b/NEWS.md
index 966f6f90..ba71c9a8 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,10 @@
 # PolyChaos Release Notes
 
+## Version 0.2.1
+- increased code coverage
+- examples work with new type hierarchy
+- updated dependencies
+
 ## Version 0.2.0
 - new type hierarchy that involves four abstract types:
     - `AbstractMeasure`
diff --git a/Project.toml b/Project.toml
index 160a9c63..69653458 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "PolyChaos"
 uuid = "8d666b04-775d-5f6e-b778-5ac7c70f65a3"
 authors = ["timueh <tillmann.muehlpfordt@kit.edu>"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 AdaptiveRejectionSampling = "c75e803d-635f-53bd-ab7d-544e482d8c75"
diff --git a/src/PolyChaos.jl b/src/PolyChaos.jl
index 1e3ece82..4a967900 100644
--- a/src/PolyChaos.jl
+++ b/src/PolyChaos.jl
@@ -16,14 +16,14 @@ include("typesMeasures.jl")
 include("typesQuad.jl")
 include("typesOrthoPolys.jl")
 include("typesTensor.jl")
-include("multiIndices.jl")
+include("multi_indices.jl")
 include("evaluate.jl")
 include("quadrature_rules.jl")
 include("recurrence_coefficients_monic.jl")
 include("stieltjes.jl")
 include("densities.jl")
 include("auxfuns.jl")
-include("scalarproduct.jl")
+include("scalar_product.jl")
 include("tensor.jl")
 include("show.jl")
 include("polynomial_chaos.jl")
diff --git a/src/densities.jl b/src/densities.jl
index 12ce2195..70b7706c 100644
--- a/src/densities.jl
+++ b/src/densities.jl
@@ -1,8 +1,6 @@
 export  w_legendre,
         build_w_jacobi,
-        build_w_jacobi01,
         w_jacobi,
-        w_jacobi01,
         w_laguerre,
         w_hermite,
         build_w_genhermite,
@@ -17,24 +15,18 @@ export  w_legendre,
         build_w_beta,
         build_w_gamma
 
+_throwError(t) = throw(DomainError(t, "not in support"))
+
 function w_legendre(t)
-    -1. <= t <= 1. ? 1. : throw(error("$t not in support"))
+    -1. <= t <= 1. ? 1. : _throwError(t)
 end
 
 function build_w_jacobi(a,b)
     return t->w_jacobi(t,a,b)
 end
-function w_jacobi(t,a,b)
-    -1. <=t<= 1. ? ((1-t)^a*(1+t)^b) : throw(error("$t not in support"))
-end
-
-function build_w_jacobi01(a,b)
-    @assert a > -1. && b > -1. "Invalid shape parameters"
-    return t->w_jacobi01(t,a,b)
-end
 
-function w_jacobi01(t,a,b)
-    0. <=t<= 1. ? ((1-t)^a*t^b) : throw(error("$t not in support"))
+function w_jacobi(t,a,b)
+    -1. <= t <= 1. ? (1-t)^a*(1+t)^b : _throwError(t)
 end
 
 function w_hermite(t)
@@ -45,7 +37,7 @@ function build_w_genhermite(mu)
     return t->w_genhermite(t,mu)
 end
 
-function w_genhermite(t,μ::Float64)
+function w_genhermite(t,μ)
     abs(t)^(2*μ)*exp(-t^2)
 end
 
@@ -53,51 +45,44 @@ function build_w_genlaguerre(a)
     return t->w_genlaguerre(t,a)
 end
 
-function w_genlaguerre(t,a)
-    t>=0. ? Float64(t)^a*exp(-Float64(t)) : throw(error("$t not in support"))
-end
-
 function w_laguerre(t)
-    t>=0. ? exp(-t) : throw(error("$t not in support"))
+    t >= 0. ? exp(-t) : _throwError(t)
 end
 
 function w_meixner_pollaczek(t,lambda,phi)
-    1/(2*pi)*exp((2*phi-pi)*t)*abs(gamma(lambda+im*t))^2
+    1 / (2pi) * exp((2*phi - pi)*t) * abs(gamma(lambda+im*t))^2
 end
 
 function build_w_meixner_pollaczek(lambda,phi)
     t->w_meixner_pollaczek(t,lambda,phi)
 end
 
-
 ##################################################
 # probability density functions
 function w_gaussian(t)
-    1/(sqrt(2*pi))*exp(-0.5*t^2)
+    1 / (sqrt(2*pi))*exp(-0.5*t^2)
 end
 
 function build_w_beta(α,β)
     return t->w_beta(t,α,β)
 end
+
 function w_beta(t,α,β)
-    t^(α-1)*(1-t)^(β-1)/beta(α,β)
+    -1 <= t <= 1 ? t^(α-1)*(1-t)^(β-1)/beta(α,β) : _throwError(t)
 end
 
-
-
 function build_w_gamma(α)
     return t->w_gamma(t,α)
 end
+
 function w_gamma(t,α)
-    t>=0. ? (1/gamma(α)*Float64(t)^(α-1)*exp(-t)) : (error("$t not in support"))
+    t >= 0. ? (1/gamma(α)*Float64(t)^(α-1)*exp(-t)) : _throwError(t)
 end
 
-
 function w_uniform01(t)
-    0. <=t<= 1. ? 1. : (error("$t not in support"))
+    0. <= t <= 1. ? 1. : _throwError(t)
 end
 
-
 function w_logistic(t)
     0.25*sech(0.5t)^2
 end
diff --git a/src/multiIndices.jl b/src/multiIndices.jl
deleted file mode 100644
index a50aa397..00000000
--- a/src/multiIndices.jl
+++ /dev/null
@@ -1,146 +0,0 @@
-export calculateMultiIndices, findUnivariateIndices, calculateMultiIndices_interaction
-
-function calculateMultiIndices(d::Int, n::Int)
-    # d denotes dimension of random variables/number of sources of uncertainty,
-    # n the maximum degree of multivariate basis
-    # function to calculate indices of multivariate basis following the algorithm
-    # from "Spectral Methods for Uncertainty Quantification; Le Maitre, Knio;2014" p.516-517
-    # No::BigInt = factorial(BigInt(d+n))/(factorial(Bigint(d))*factorial(BigInt(n)));  #number of polynomials of multivariate basis
-  n < 0 && throw(error("maximum degree must be non-negative"))
-  d <= 0 && throw(error("number of uncertainties must be positive"))
-  # catch case n == 0 --> No-d==0
-  n == 0 && return zeros(Int64,1,d)
-  # non-pathological cases begin here
-  No = numberPolynomials(d,n)
-  inds = vcat(zeros(Int64,1,d),Matrix(1I,d,d),zeros(Int64, No-d-1, d));  #initiate index matrix for basis
-  pi=ones(Int64,No,d);
-
-  for k=2:No
-    g = 0;
-    for l=1:d
-      pi[k,l]=sum(pi[k-1,:])-g;
-      g=g+pi[k-1,l];
-    end
-  end
-
-  P=d+1;
-  for k=2:n
-    L = P;
-    for j=1:d, m=L-pi[k,j]+1:L
-        P=P+1;
-        inds[P,:]=inds[m,:];
-        inds[P,j]=inds[P,j]+1;
-    end
-  end
-
-  return inds
-end
-
-function calculateMultiIndices(a::Array{T,1}) where T<:Integer
-  #a contains maximum degrees of univariate bases
-  #function to calculate indices of multivariate basis following the algorithm
-  #from "Spectral Methods for Uncertainty Quantification; Le Maitre, Knio;2014" p.516-517
-  d = length(a)
-  n = maximum(a)
-  No = numberPolynomials(d,n)  #number of polynomials of multivariate basis
-  if d==0
-    print("Input zero in indicesMulti!")
-    return zeros(Int64, 1,d)
-  end
-
-  # catch case n == 0 --> No-d==0
-  if n==0
-    inds=zeros(Int64,1,d)
-    return inds
-  end
-
-  inds = vcat(zeros(Int64,1,d),Matrix(1I,d,d),zeros(Int64, No-d-1, d));  #initiate index matrix for basis
-  pi=ones(Int64,No,d);
-
-  for k=2:No
-    g = 0;
-    for l=1:d
-      pi[k,l]=sum(pi[k-1,:])-g;
-      g=g+pi[k-1,l];
-    end
-  end
-
-  P=d+1;
-  for k=2:n
-    L = P;
-    for j=1:d, m=L-pi[k,j]+1:L
-        P=P+1;
-        inds[P,:]=inds[m,:];
-        inds[P,j]=inds[P,j]+1;
-    end
-  end
-
-  #remove elements of high order
-  #TODO more efficient way to do that?
-  for j=1:length(a)
-    k=size(inds,1)
-    l=1
-    while l<=k
-      if inds[l,j]>a[j]
-        inds=vcat(inds[1:l-1,:],inds[l+1:end,:])
-        k-=1
-        l-=1
-      end
-      l+=1
-    end
-  end
-
-  return inds
-end
-
-
-"""
-    computes the number of polynomials with a multivariate basis
-    `(d+n)!/(d!+n!)`
-"""
-function numberPolynomials(a::Array{T,1}) where T<:Integer
-    d = length(a)
-    n = maximum(a)
-
-    x = max(d,n)
-    y = min(d,n)
-
-    return (prod(UInt128(x+1):UInt128(d+n))/factorial(UInt128(y)))
-end
-
-function numberPolynomials(d::Int64, n::Int64)
-    x = max(d,n)
-    y = min(d,n)
-    return UInt128(prod(UInt128(x+1):UInt128(d+n))/factorial(UInt128(y)))
-end
-
-"""
-    findUnivariateIndices(i::Int64,ind::Matrix{Int64})::Vector{Int64}
-Given the multi-index `ind` this function returns all entries of the multivariate basis
-that correspond to the `i`th univariate basis.
-"""
-function findUnivariateIndices(i::Int64,ind::Matrix{Int64})::Vector{Int64}
-  l,p = size(ind)
-  @assert i<=p "basis is $p-variate, you requested $i-variate"
-  deg = ind[end,end]
-  @assert deg>=0 "invalid degree"
-  col = ind[:,i]
-  myind = zeros(Int64,deg)
-  for deg_ = 1:deg
-      myind[deg_] = findfirst(x->x==deg_,col)
-  end
-  pushfirst!(myind,1)
-end
-
-#################################################################
-function calculateMultiIndices_interaction(nξ::Int,deg::Int,j::Int,p::Int)
-    inds = calculateMultiIndices(nξ,deg)
-    get_interaction(inds,j,p)
-end
-
-function get_interaction(inds::Matrix{Int},j::Int,p::Int)
-    j < 0 && throw(error("interaction order must be non-negative"))
-    j > size(inds,2) && throw(error("interaction order cannot be greater than number of uncertainties"))
-    Iterators.filter(x -> count(!iszero,x) == j && sum(x) == p, eachrow(inds))
-    # alternative -> collect
-end
diff --git a/src/multi_indices.jl b/src/multi_indices.jl
new file mode 100644
index 00000000..1540b583
--- /dev/null
+++ b/src/multi_indices.jl
@@ -0,0 +1,79 @@
+export  calculateMultiIndices,
+        findUnivariateIndices
+
+function calculateMultiIndices(d::Int, n::Int)
+    # d denotes dimension of random variables/number of sources of uncertainty,
+    # n the maximum degree of multivariate basis
+    # function to calculate indices of multivariate basis following the algorithm
+    # from "Spectral Methods for Uncertainty Quantification; Le Maitre, Knio;2014" p.516-517
+    # No::BigInt = factorial(BigInt(d+n))/(factorial(Bigint(d))*factorial(BigInt(n)));  #number of polynomials of multivariate basis
+  n < 0 && throw(DomainError(n, "maximum degree must be non-negative"))
+  d <= 0 && throw(DomainError(d, "number of uncertainties must be positive"))
+  # catch case n == 0 --> No-d==0
+  n == 0 && return zeros(Int64,1,d)
+  # non-pathological cases begin here
+  No = numberPolynomials(d,n)
+  inds = vcat(zeros(Int64,1,d),Matrix(1I,d,d),zeros(Int64, No-d-1, d));  #initiate index matrix for basis
+  pi=ones(Int64,No,d);
+
+  for k=2:No
+    g = 0;
+    for l=1:d
+      pi[k,l]=sum(pi[k-1,:])-g;
+      g=g+pi[k-1,l];
+    end
+  end
+
+  P=d+1;
+  for k=2:n
+    L = P;
+    for j=1:d, m=L-pi[k,j]+1:L
+        P=P+1;
+        inds[P,:]=inds[m,:];
+        inds[P,j]=inds[P,j]+1;
+    end
+  end
+
+  return inds
+end
+
+"""
+    computes the number of polynomials with a multivariate basis
+    `(d+n)!/(d!+n!)`
+"""
+
+function numberPolynomials(d::Int64, n::Int64)
+    x, y = max(d,n), min(d,n)
+    return UInt128(prod(UInt128(x+1):UInt128(d+n)) ÷ factorial(UInt128(y)))
+end
+
+"""
+    findUnivariateIndices(i::Int64,ind::Matrix{Int64})::Vector{Int64}
+Given the multi-index `ind` this function returns all entries of the multivariate basis
+that correspond to the `i`th univariate basis.
+"""
+function findUnivariateIndices(i::Int64,ind::Matrix{Int64})::Vector{Int64}
+  l,p = size(ind)
+  @assert i<=p "basis is $p-variate, you requested $i-variate"
+  deg = ind[end,end]
+  @assert deg>=0 "invalid degree"
+  col = ind[:,i]
+  myind = zeros(Int64,deg)
+  for deg_ = 1:deg
+      myind[deg_] = findfirst(x->x==deg_,col)
+  end
+  pushfirst!(myind,1)
+end
+
+#################################################################
+# function calculateMultiIndices_interaction(nξ::Int,deg::Int,j::Int,p::Int)
+#     inds = calculateMultiIndices(nξ,deg)
+#     get_interaction(inds,j,p)
+# end
+
+# function get_interaction(inds::Matrix{Int},j::Int,p::Int)
+#     j < 0 && throw(error("interaction order must be non-negative"))
+#     j > size(inds,2) && throw(error("interaction order cannot be greater than number of uncertainties"))
+#     Iterators.filter(x -> count(!iszero,x) == j && sum(x) == p, eachrow(inds))
+#     # alternative -> collect
+# end
diff --git a/src/quadrature_rules.jl b/src/quadrature_rules.jl
index 3a366bd5..b1db43af 100644
--- a/src/quadrature_rules.jl
+++ b/src/quadrature_rules.jl
@@ -96,13 +96,13 @@ function golubwelsch(α::Vector{<:Real}, β::Vector{<:Real}, maxiter::Int=30)
     a[idx], (β0*w.^2)[idx]
 end
 
-golubwelsch(op::OrthoPoly) = golubwelsch(op.α,op.β)
+golubwelsch(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) = golubwelsch(op.α,op.β)
 
 """
     gauss(N::Int,α::Vector{<:Real},β::Vector{<:Real})
     gauss(α::Vector{<:Real},β::Vector{<:Real})
-    gauss(N::Int,op::OrthoPoly)
-    gauss(op::OrthoPoly)
+    gauss(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly})
+    gauss(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly})
 Gauss quadrature rule, also known as Golub-Welsch algorithm
 
 `gauss()` generates the `N` Gauss quadrature nodes and weights for a given weight function.
@@ -124,15 +124,18 @@ function gauss(N::Int,α::Vector{<:Real},β::Vector{<:Real})
     @assert N0 >= N  "not enough recurrence coefficients"
     @inbounds golubwelsch(α[1:N],β[1:N])
 end
+
 gauss(α::Vector{<:Real},β::Vector{<:Real}) = gauss(length(α)-1,α,β)
-gauss(N::Int,op::OrthoPoly) = gauss(N::Int,op.α,op.β)
-gauss(op::OrthoPoly) = gauss(op.α,op.β)
+
+gauss(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) = gauss(N::Int,op.α,op.β)
+
+gauss(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) = gauss(op.α,op.β)
 
 """
     radau(N::Int,α::Vector{<:Real},β::Vector{<:Real},end0::Real)
     radau(α::Vector{<:Real},β::Vector{<:Real},end0::Real)
-    radau(N::Int,op::OrthoPoly,end0::Real)
-    radau(op::OrthoPoly,end0::Real)
+    radau(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real)
+    radau(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real)
 Gauss-Radau quadrature rule.
 Given a weight function encoded by the recurrence coefficients `(α,β)`for the associated
 orthogonal polynomials, the function generates the
@@ -164,14 +167,14 @@ function radau(N::Int,α_::Vector{<:Real},β::Vector{<:Real},end0::Real)
     gauss(N+1,α,β)
 end
 radau(α::Vector{<:Real},β::Vector{<:Real},end0::Real) = radau(length(α)-2,α,β,end0)
-radau(N::Int,op::OrthoPoly,end0::Real) = radau(N,op.α,op.β,end0)
-radau(op::OrthoPoly,end0::Real) = radau(op.α,op.β,end0::Real)
+radau(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real) = radau(N,op.α,op.β,end0)
+radau(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real) = radau(op.α,op.β,end0::Real)
 
 """
     lobatto(N::Int,α::Vector{<:Real},β::Vector{<:Real},endl::Real,endr::Real)
     lobatto(α::Vector{<:Real},β::Vector{<:Real},endl::Real,endr::Real)
-    lobatto(N::Int,op::OrthoPoly,endl::Real,endr::Real)
-    lobatto(op::OrthoPoly,endl::Real,endr::Real)
+    lobatto(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real)
+    lobatto(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real)
 Gauss-Lobatto quadrature rule.
 Given a weight function encoded by the recurrence coefficients for the associated
 orthogonal polynomials, the function generates
@@ -209,5 +212,5 @@ function lobatto(N::Int,α_::Vector{<:Real},β_::Vector{<:Real},endl::Real,endr:
     gauss(N+2,α,β)
 end
 lobatto(α::Vector{<:Real},β::Vector{<:Real},endl::Real,endr::Real) = lobatto(length(α)-3,α,β,endl,endr)
-lobatto(N::Int,op::OrthoPoly,endl::Real,endr::Real) = lobatto(N,op.α,op.β,endl,endr)
-lobatto(op::OrthoPoly,endl::Real,endr::Real) = lobatto(op.α,op.β,endl,endr)
+lobatto(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real) = lobatto(N,op.α,op.β,endl,endr)
+lobatto(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real) = lobatto(op.α,op.β,endl,endr)
diff --git a/src/recurrence_coefficients_monic.jl b/src/recurrence_coefficients_monic.jl
index aa0090d6..755b6419 100644
--- a/src/recurrence_coefficients_monic.jl
+++ b/src/recurrence_coefficients_monic.jl
@@ -194,17 +194,49 @@ that are orthogonal on ``(0,1)`` relative to ``w(t) = 1``.
 rm_legendre01(N::Int) = rm_jacobi01(N)
 
 
-function rm_chebyshev1(N::Int)
-    @assert N>=0 "N has to be non-negative"
-    α = zeros(Float64,N)
-    if N == 1
-        return α, [pi]
-    elseif N == 2
-        return α, [pi; 0.5]
+"""
+    rm_meixner_pollaczek(N::Int,lambda::Real,phi::Real)
+    rm_meixner_pollaczek(N::Int,lambda::Real)
+
+ Creates `N` recurrence coefficients for monic
+ Meixner-Pollaczek polynomials with parameters λ and ϕ. These are orthogonal on
+ ``[-\\infty,\\infty]`` relative to the weight function ``w(t)=(2 \\pi)^{-1} \\exp{(2 \\phi-\\pi)t} |\\Gamma(\\lambda+ i t)|^2``.
+
+ The call `rm_meixner_pollaczek(n,lambda)` is the same as `rm_meixner_pollaczek(n,lambda,pi/2)`.
+"""
+function rm_meixner_pollaczek(N::Int,lambda::Real,phi::Real)
+    @assert N>=0 && lambda>0. && phi>0. "parameter(s) out of range"
+    N == 0 && return Array{Float64,1}(undef,0), Array{Float64,1}(undef,0)
+    n=1:N;
+    sinphi=sin(phi); lam2=2*lambda;
+    ab=zeros(Float64,N,2)
+    if sinphi==1
+      ab[:,1]=zeros(N);
     else
-        return α, pushfirst!(pushfirst!(0.25*ones(N-2),0.5),pi)
+        ab[:,1]= -(n .+ (lambda-1))/tan(phi);
     end
+    ab[:,2]= (n .- 1).*(n .+ (lam2-2))/(4*sinphi^2);
+    ab[1,2]=gamma(lam2)/(2*sinphi)^lam2;
+    return ab[:,1], ab[:,2]
 end
+rm_meixner_pollaczek(N::Int,lambda::Real) = rm_meixner_pollaczek(N,lambda,pi/2)
+
+
+###################################################################
+###################################################################
+###################################################################
+
+# function rm_chebyshev1(N::Int)
+#     @assert N>=0 "N has to be non-negative"
+#     α = zeros(Float64,N)
+#     if N == 1
+#         return α, [pi]
+#     elseif N == 2
+#         return α, [pi; 0.5]
+#     else
+#         return α, pushfirst!(pushfirst!(0.25*ones(N-2),0.5),pi)
+#     end
+# end
 #     rm_hahn(N::Int,a::Real,b::Real)
 #     rm_hahn(N::Int,a::Real)
 #     rm_hahn(N::Int)
@@ -248,30 +280,3 @@ end
 # end
 # rm_hahn(N::Int,a::Real) = rm_hahn(N,a,a)
 # rm_hahn(N::Int) = rm_hahn(N,0.,0.)
-
-"""
-    rm_meixner_pollaczek(N::Int,lambda::Real,phi::Real)
-    rm_meixner_pollaczek(N::Int,lambda::Real)
-
- Creates `N` recurrence coefficients for monic
- Meixner-Pollaczek polynomials with parameters λ and ϕ. These are orthogonal on
- ``[-\\infty,\\infty]`` relative to the weight function ``w(t)=(2 \\pi)^{-1} \\exp{(2 \\phi-\\pi)t} |\\Gamma(\\lambda+ i t)|^2``.
-
- The call `rm_meixner_pollaczek(n,lambda)` is the same as `rm_meixner_pollaczek(n,lambda,pi/2)`.
-"""
-function rm_meixner_pollaczek(N::Int,lambda::Real,phi::Real)
-    @assert N>=0 && lambda>0. && phi>0. "parameter(s) out of range"
-    N == 0 && return Array{Float64,1}(undef,0), Array{Float64,1}(undef,0)
-    n=1:N;
-    sinphi=sin(phi); lam2=2*lambda;
-    ab=zeros(Float64,N,2)
-    if sinphi==1
-      ab[:,1]=zeros(N);
-    else
-        ab[:,1]= -(n .+ (lambda-1))/tan(phi);
-    end
-    ab[:,2]= (n .- 1).*(n .+ (lam2-2))/(4*sinphi^2);
-    ab[1,2]=gamma(lam2)/(2*sinphi)^lam2;
-    return ab[:,1], ab[:,2]
-end
-rm_meixner_pollaczek(N::Int,lambda::Real) = rm_meixner_pollaczek(N,lambda,pi/2)
diff --git a/src/scalarproduct.jl b/src/scalar_product.jl
similarity index 100%
rename from src/scalarproduct.jl
rename to src/scalar_product.jl
diff --git a/src/stieltjes.jl b/src/stieltjes.jl
index 9ee1989c..9039d8fe 100644
--- a/src/stieltjes.jl
+++ b/src/stieltjes.jl
@@ -88,31 +88,6 @@ function lanczos(N::Int64,nodes::Vector{Float64},weights::Vector{Float64};remove
     @inbounds p0[Base.OneTo(N)], p1[Base.OneTo(N)]
 end
 
-
-function mcdis_twologistics(n::Int,p1::Vector{Float64},p2::Vector{Float64};Mmax::Int=100,eps0::Real=1e-9)
-    M0 = n
-    Mcap = 0
-    Mi = M0
-    α, β = zeros(n), zeros(n)
-    for i=1:Mmax
-        nai, wai = quadrature_logistic(Mi,p1[1],p1[2])
-        nbi, wbi = quadrature_logistic(Mi,p2[1],p2[2])
-        wai, wbi = 0.5*wai, 0.5*wbi
-        ni, wi   = [nai; nbi], [wai; wbi]
-        α_, β_ = stieltjes(n,ni,wi)
-        err = maximum( abs(β_[k]-β[k])/β_[k] for k=1:n )
-        display("err = $err")
-        if err<=eps0
-            Mcap = i-1
-            return α_, β_
-        else
-            α, β, = α_, β_
-            i==1 ? Mi = M0+1 : Mi = Int(2^(floor((i-1)/5))*n)
-        end
-    end
-    warn("Algorithm did not terminate after $Mmax iterations.")
-end
-
 """
     mcdiscretization(N::Int64,quads::Vector{},discretemeasure::Matrix{Float64}=zeros(0,2);discretization::Function=stieltjes,Nmax::Integer=300,ε::Float64=1e-8,gaussquad::Bool=false)
 This routine returns ``N`` recurrence coefficients of the polynomials that are
diff --git a/test/auxfuns.jl b/test/auxfuns.jl
index 0d967245..f02eddf0 100644
--- a/test/auxfuns.jl
+++ b/test/auxfuns.jl
@@ -1,5 +1,17 @@
 using Test, PolyChaos
 
+@test nw(EmptyQuad()) |> typeof == Array{Float64,2}
+@test nw(EmptyQuad()) |> length == 0
+
+nodes, weights = gauss(rm_legendre01(5)...)
+myQuad = Quad("myQuad", length(nodes), nodes, weights)
+foo(x) = 6*x^5
+res = 1
+@test isapprox(integrate(foo, nodes, weights), res)
+@test isapprox(integrate(foo, myQuad), res)
+@test_throws DomainError integrate(foo, EmptyQuad())
+@test isapprox(integrate(foo, Uniform01OrthoPoly(4)), res)
+
 d, nunc = 5, 10
 op = genLaguerreOrthoPoly(d,1.23,addQuadrature = false)
 opq = genLaguerreOrthoPoly(d,1.23)
diff --git a/test/densities.jl b/test/densities.jl
new file mode 100644
index 00000000..65443464
--- /dev/null
+++ b/test/densities.jl
@@ -0,0 +1,20 @@
+using PolyChaos, Test
+
+outliers = [-1.3, 1.3]
+
+@test_throws DomainError PolyChaos._throwError(3)
+
+densities = [w_legendre,
+			 w_uniform01,
+			 build_w_jacobi(1.2,3.4),
+			 build_w_beta(1.2,3.4)
+]
+
+@testset "Supports of densities" begin
+	for ρ in densities, x in outliers
+		@test_throws DomainError ρ(x)
+	end
+end
+
+@test_throws DomainError w_laguerre(-1)
+@test_throws DomainError build_w_gamma(1)(-1)
\ No newline at end of file
diff --git a/test/evaluate.jl b/test/evaluate.jl
index 506753c6..d14e3308 100644
--- a/test/evaluate.jl
+++ b/test/evaluate.jl
@@ -5,11 +5,23 @@ op = genLaguerreOrthoPoly(deg,1.23)
 
 n = 10
 X = [ exp(3 + rand()), exp.(3 .+ rand(n)) ]
+
 @testset "evaluate univariate basis at point(s)" begin
     for x in X
         @test evaluate(collect(0:deg),x,op.α,op.β) == evaluate(collect(0:deg),x,op)
         for d in 0:deg
-            @test evaluate(d,x,op.α,op.β) == evaluate(d,x,op)
+            if typeof(x) <: Real
+                @test evaluate(d,x,op.α,op.β) == evaluate(d,x,op) == evaluate(d,[x],op.α,op.β)
+            else
+                @test evaluate(d,x,op.α,op.β) == evaluate(d,x,op)
+            end
+        end
+
+        res = hcat(map(i->evaluate(i,x,op.α,op.β),collect(0:op.deg))...)
+        if typeof(x) <: Real
+            @test evaluate(x,op) == evaluate([x], op) == evaluate(collect(0:op.deg),[x],op) == evaluate(collect(0:op.deg),[x],op.α,op.β) == res
+        else
+            @test evaluate(x, op) == evaluate(collect(0:op.deg),x,op) == evaluate(collect(0:op.deg),x,op.α,op.β) == res
         end
     end
 end
@@ -29,6 +41,12 @@ X = [ exp.(3 .+ rand(nop)) , exp.(3 .+ rand(n,nop)) ]
 @testset "evaluate multivariate basis at point(s)" begin
     for x in X
         @test evaluate(mop.ind,x,coeffs(mop)...) == evaluate(mop.ind,x,mop) == evaluate(x,mop)
+        if typeof(x) <: Vector{<:Real}
+            @test evaluate(mop.ind,x,coeffs(mop)...) == evaluate(mop.ind,x,mop) == evaluate(x,mop) == evaluate(mop.ind,reshape(x,1,length(x)),coeffs(mop)...)
+            @test evaluate(x,mop) == evaluate(mop.ind,reshape(x,1,length(x)),mop) == evaluate(mop.ind,reshape(x,1,length(x)),coeffs(mop)...)
+        else
+            @test evaluate(x,mop) == evaluate(mop.ind,x,mop) == evaluate(mop.ind,x,coeffs(mop)...)
+        end
         for d in eachrow(mop.ind[1:3,:])
             @test evaluate(collect(d),x,coeffs(mop)...) == evaluate(collect(d),x,mop)
         end
diff --git a/test/multi_indices.jl b/test/multi_indices.jl
new file mode 100644
index 00000000..f9b36799
--- /dev/null
+++ b/test/multi_indices.jl
@@ -0,0 +1,26 @@
+using Test, PolyChaos
+
+A = [ 0  0  0  0;
+ 1  0  0  0;
+ 0  1  0  0;
+ 0  0  1  0;
+ 0  0  0  1;
+ 2  0  0  0;
+ 1  1  0  0;
+ 1  0  1  0;
+ 1  0  0  1;
+ 0  2  0  0;
+ 0  1  1  0;
+ 0  1  0  1;
+ 0  0  2  0;
+ 0  0  1  1;
+ 0  0  0  2]
+
+@test calculateMultiIndices(4, 2) == A
+@test calculateMultiIndices(3, 0) == zeros(Int64, 1, 3)
+@test_throws DomainError calculateMultiIndices( 3, -1)
+@test_throws DomainError calculateMultiIndices( 0,  1)
+@test_throws DomainError calculateMultiIndices(-1,  1)
+
+@test PolyChaos.numberPolynomials(3,4) == factorial(3+4) ÷ ( factorial(3)*factorial(4) )
+@test findUnivariateIndices(2,A) == [1,3,10]
\ No newline at end of file
diff --git a/test/quadrature_rules.jl b/test/quadrature_rules.jl
index 4af1e706..43aa1299 100644
--- a/test/quadrature_rules.jl
+++ b/test/quadrature_rules.jl
@@ -123,3 +123,18 @@ intervals = [(-5.,5.),(-Inf,Inf),(-Inf,-1.),(1.,Inf),(-3.,5.),(-10.,1.)]
             @test isapprox(exp.(-abs.(nod2))'*wei2-analytic,0;atol=1e-4)
         end
 end
+
+degree = 9
+N = degree + 1
+α, β = rm_legendre01(degree + 1)
+op = Uniform01OrthoPoly(degree)
+endl, endr = 0, 1
+
+myMeas = Measure("myMeasure", t->1, (0, 1), true)
+op_num = OrthoPoly("myOp", degree, myMeas)
+
+@testset "OrthoPoly overloading" begin
+    @test gauss(N-1, α, β) == gauss(α, β) == gauss(N-1, op) == gauss(op)
+    @test radau(N-2, α, β, endr) == radau(α, β, endr) == radau(N-2, op, endr) == radau(op, endr)
+    @test lobatto(N-3, α, β, endl, endr) == lobatto(α, β, endl, endr) == lobatto(N-3, op, endl, endr) == lobatto(op, endl, endr)
+end
\ No newline at end of file
diff --git a/test/recurrence_coefficients.jl b/test/recurrence_coefficients.jl
index 451e5892..001328f0 100644
--- a/test/recurrence_coefficients.jl
+++ b/test/recurrence_coefficients.jl
@@ -47,6 +47,7 @@ end
                     @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
                     αβcom = coeffs(JacobiOrthoPoly(n-1, al, be, addQuadrature=false))
                     @test isapprox(norm(αβref-[αβcom[:,1];αβcom[:,2]],Inf),0.;atol=tol)
+                    @test rm_jacobi(n,al) == rm_jacobi(n,al,al)
                     close(myfile)
             end
     end
@@ -58,6 +59,7 @@ end
                     myfile = open("dataRecCoeffs/jac01$(n)al$(al)be$(be).txt")
                     αβref = parse.(Float64,readlines(myfile))
                     αβcom = rm_jacobi01(n,al,be)
+                    @test rm_jacobi01(n,al) == rm_jacobi01(n,al,al)
                     @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
                     close(myfile)
             end
@@ -89,6 +91,9 @@ meipol = 0.1:0.2:2
                 αβcom = coeffs(MeixnerPollaczekOrthoPoly(n-1, lambda, phi, addQuadrature=false))
                 @test isapprox(norm(αβref-[αβcom[:,1];αβcom[:,2]],Inf),0.;atol=tol)
                 close(myfile)
+                @test rm_meixner_pollaczek(n,lambda) == rm_meixner_pollaczek(n,lambda,pi/2)
+                @test rm_meixner_pollaczek(n,lambda,pi/2)[1] == zeros(n)
             end
+
     end
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 0d2d2a76..cb45805b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,4 +7,6 @@ include("quadrature_rules.jl")
 include("show.jl")
 include("types.jl")
 include("auxfuns.jl")
+include("densities.jl")
 include("polynomial_chaos.jl")
+include("multi_indices.jl")
\ No newline at end of file
diff --git a/todos.md b/todos.md
deleted file mode 100644
index 597f7a69..00000000
--- a/todos.md
+++ /dev/null
@@ -1,16 +0,0 @@
-# Todos for PolyChaos.jl
-
-## High priority
-  - documentation:
-    - make conversion to Latex work
-    - examples:
-      - quadrature rules (fejer, fejer2, clenshaw-curtis, gauss, radau, lobatto)
-      - ~~discretization procedures (stieltjes, lanczos), see `test/discretization.jl`~~
-      - `quadgp()`
-      - optimizing Rosenbrock/quadratic function
-      - DC-OPF with stochastic uncertainties
-
-
-## Low priority
-  - orthonormal polynomials
-  - arbitrary polynomials