diff --git a/Project.toml b/Project.toml
index dc781c6..04719d8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,18 +4,14 @@ authors = ["Joseph Tindall <joeytindall94@googlemail.com>", "Ryan Levy <rlevy@fl
 version = "0.1.0"
 
 [deps]
-DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
-Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 ITensorNetworks = "2919e153-833c-4bdc-8836-1ea460a35fc7"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-ITensorNetworks = "0.6.0"
+ITensorNetworks = "0.8.1"
diff --git a/src/elementary_functions.jl b/src/elementary_functions.jl
index 1d18ff7..ce4370b 100644
--- a/src/elementary_functions.jl
+++ b/src/elementary_functions.jl
@@ -36,6 +36,7 @@ function exp_itensornetwork(
   s::IndsNetworkMap;
   k=default_k_value(),
   a=default_a_value(),
+  c=default_c_value(),
   dimension::Int=default_dimension(),
 )
   ψ = const_itensornetwork(s)
@@ -45,6 +46,8 @@ function exp_itensornetwork(
     ψ[v] = ITensor(exp(a / Lx) * exp.(k * index_values_to_scalars(s, sind)), inds(ψ[v]))
   end
 
+  ψ[first(dimension_vertices(ψ, dimension))] *= c
+
   return ψ
 end
 
@@ -54,13 +57,11 @@ function cosh_itensornetwork(
   s::IndsNetworkMap;
   k=default_k_value(),
   a=default_a_value(),
+  c=default_c_value(),
   dimension::Int=default_dimension(),
 )
-  ψ1 = exp_itensornetwork(s; a, k, dimension)
-  ψ2 = exp_itensornetwork(s; a=-a, k=-k, dimension)
-
-  ψ1[first(vertices(ψ1))] *= 0.5
-  ψ2[first(vertices(ψ1))] *= 0.5
+  ψ1 = exp_itensornetwork(s; a, k, c=0.5 * c, dimension)
+  ψ2 = exp_itensornetwork(s; a=-a, k=-k, c=0.5 * c, dimension)
 
   return ψ1 + ψ2
 end
@@ -71,13 +72,11 @@ function sinh_itensornetwork(
   s::IndsNetworkMap;
   k=default_k_value(),
   a=default_a_value(),
+  c=default_c_value(),
   dimension::Int=default_dimension(),
 )
-  ψ1 = exp_itensornetwork(s; a, k, dimension)
-  ψ2 = exp_itensornetwork(s; a=-a, k=-k, dimension)
-
-  ψ1[first(vertices(ψ1))] *= 0.5
-  ψ2[first(vertices(ψ1))] *= -0.5
+  ψ1 = exp_itensornetwork(s; a, k, c=0.5 * c, dimension)
+  ψ2 = exp_itensornetwork(s; a=-a, k=-k, c=-0.5 * c, dimension)
 
   return ψ1 + ψ2
 end
@@ -88,16 +87,19 @@ function tanh_itensornetwork(
   s::IndsNetworkMap;
   k=default_k_value(),
   a=default_a_value(),
+  c=default_c_value(),
   nterms::Int=default_nterms(),
   dimension::Int=default_dimension(),
 )
   ψ = const_itensornetwork(s)
   for n in 1:nterms
     ψt = exp_itensornetwork(s; a=-2 * n * a, k=-2 * k * n, dimension)
-    ψt[first(vertices(ψt))] *= 2 * ((-1)^n)
+    ψt[first(dimension_vertices(ψ, dimension))] *= 2 * ((-1)^n)
     ψ = ψ + ψt
   end
 
+  ψ[first(dimension_vertices(ψ, dimension))] *= c
+
   return ψ
 end
 
@@ -107,13 +109,11 @@ function cos_itensornetwork(
   s::IndsNetworkMap;
   k=default_k_value(),
   a=default_a_value(),
+  c=default_c_value(),
   dimension::Int=default_dimension(),
 )
-  ψ1 = exp_itensornetwork(s; a=a * im, k=k * im, dimension)
-  ψ2 = exp_itensornetwork(s; a=-a * im, k=-k * im, dimension)
-
-  ψ1[first(vertices(ψ1))] *= 0.5
-  ψ2[first(vertices(ψ1))] *= 0.5
+  ψ1 = exp_itensornetwork(s; a=a * im, k=k * im, c=0.5 * c, dimension)
+  ψ2 = exp_itensornetwork(s; a=-a * im, k=-k * im, c=0.5 * c, dimension)
 
   return ψ1 + ψ2
 end
@@ -124,13 +124,11 @@ function sin_itensornetwork(
   s::IndsNetworkMap;
   k=default_k_value(),
   a=default_a_value(),
+  c=default_c_value(),
   dimension::Int=default_dimension(),
 )
-  ψ1 = exp_itensornetwork(s; a=a * im, k=k * im, dimension)
-  ψ2 = exp_itensornetwork(s; a=-a * im, k=-k * im, dimension)
-
-  ψ1[first(vertices(ψ1))] *= -0.5 * im
-  ψ2[first(vertices(ψ1))] *= 0.5 * im
+  ψ1 = exp_itensornetwork(s; a=a * im, k=k * im, c=-0.5 * im * c, dimension)
+  ψ2 = exp_itensornetwork(s; a=-a * im, k=-k * im, c=0.5 * im * c, dimension)
 
   return ψ1 + ψ2
 end
@@ -138,9 +136,14 @@ end
 """Build a representation of the function f(x) = sum_{i=0}^{n}coeffs[i+1]*(x)^{i} on the graph structure specified
 by indsnetwork"""
 function polynomial_itensornetwork(
-  s::IndsNetworkMap, coeffs::Vector; dimension::Int=default_dimension()
+  s::IndsNetworkMap,
+  coeffs::Vector;
+  dimension::Int=default_dimension(),
+  k=default_k_value(),
+  c=default_c_value(),
 )
   n = length(coeffs)
+  coeffs = [c * (k^(i - 1)) for (i, c) in enumerate(coeffs)]
   #First treeify the index network (ignore edges that form loops)
   _s = indsnetwork(s)
   g = underlying_graph(_s)
@@ -179,6 +182,8 @@ function polynomial_itensornetwork(
     end
   end
 
+  ψ[first(dim_vertices)] *= c
+
   #Put the transfer tensors in, these are special tensors that
   # go on the digits (sites) that don't correspond to the desired dimension
   for v in setdiff(vertices(ψ), dim_vertices)
diff --git a/src/indexmap.jl b/src/indexmap.jl
index 02bde97..978df9a 100644
--- a/src/indexmap.jl
+++ b/src/indexmap.jl
@@ -58,7 +58,7 @@ dimensions(im::IndexMap, inds::Vector{Index}) = dimension.(inds)
 digit(im::IndexMap, ind::Index) = index_digit(im)[ind]
 digits(im::IndexMap, inds::Vector{Index}) = digit.(inds)
 function index_value_to_scalar(im::IndexMap, ind::Index, value::Int)
-  return value / (dim(ind)^digit(im, ind))
+  return (value) / (dim(ind)^digit(im, ind))
 end
 function index_values_to_scalars(im::IndexMap, ind::Index)
   return [index_value_to_scalar(im, ind, i) for i in 0:(dim(ind) - 1)]
diff --git a/src/indsnetworkmap.jl b/src/indsnetworkmap.jl
index d2e8e10..5216bbf 100644
--- a/src/indsnetworkmap.jl
+++ b/src/indsnetworkmap.jl
@@ -1,7 +1,6 @@
 using Base: Base
 using Graphs: Graphs
 using NamedGraphs: NamedGraphs
-using DataGraphs: DataGraphs
 using ITensors: ITensors
 using ITensorNetworks:
   ITensorNetworks, AbstractIndsNetwork, IndsNetwork, data_graph, underlying_graph
@@ -17,11 +16,11 @@ indexmap(inm::IndsNetworkMap) = inm.indexmap
 indtype(inm::IndsNetworkMap) = indtype(typeof(indsnetwork(inm)))
 indtype(::Type{<:IndsNetworkMap{V,I,IN,IM}}) where {V,I,IN,IM} = I
 ITensorNetworks.data_graph(inm::IndsNetworkMap) = data_graph(indsnetwork(inm))
-function DataGraphs.underlying_graph(inm::IndsNetworkMap)
+function ITensorNetworks.underlying_graph(inm::IndsNetworkMap)
   return underlying_graph(data_graph(indsnetwork(inm)))
 end
 NamedGraphs.vertextype(::Type{<:IndsNetworkMap{V,I,IN,IM}}) where {V,I,IN,IM} = V
-DataGraphs.underlying_graph_type(G::Type{<:IndsNetworkMap}) = NamedGraph{vertextype(G)}
+ITensorNetworks.underlying_graph_type(G::Type{<:IndsNetworkMap}) = NamedGraph{vertextype(G)}
 Graphs.is_directed(::Type{<:IndsNetworkMap}) = false
 
 function Base.copy(inm::IndsNetworkMap)
diff --git a/src/itensornetworkfunction.jl b/src/itensornetworkfunction.jl
index 7e7fcc5..c6c1e4e 100644
--- a/src/itensornetworkfunction.jl
+++ b/src/itensornetworkfunction.jl
@@ -1,8 +1,11 @@
 using Base: Base
-using ITensorNetworks: ITensorNetworks, AbstractITensorNetwork, data_graph, data_graph_type
+using ITensorNetworks:
+  ITensorNetworks, AbstractITensorNetwork, data_graph, data_graph_type, scalar
 using ITensors: ITensor, dim, contract, siteinds, onehot
 using Graphs: Graphs
 
+default_contraction_alg() = "bp"
+
 struct ITensorNetworkFunction{V,TN<:AbstractITensorNetwork{V},INM<:IndsNetworkMap} <:
        AbstractITensorNetwork{V}
   itensornetwork::TN
@@ -70,19 +73,25 @@ function project(fitn::ITensorNetworkFunction, ind_to_ind_value_map)
   return fitn
 end
 
-function calculate_fxyz(fitn::ITensorNetworkFunction, xs::Vector, dimensions::Vector{Int})
+function calculate_fxyz(
+  fitn::ITensorNetworkFunction,
+  xs::Vector,
+  dimensions::Vector{Int};
+  alg=default_contraction_alg(),
+  kwargs...,
+)
   ind_to_ind_value_map = calculate_ind_values(fitn, xs, dimensions)
   fitn_xyz = project(fitn, ind_to_ind_value_map)
-  return contract(fitn_xyz)[]
+  return scalar(itensornetwork(fitn_xyz); alg, kwargs...)
 end
 
-function calculate_fxyz(fitn::ITensorNetworkFunction, xs::Vector)
-  return calculate_fxyz(fitn, xs, [i for i in 1:length(xs)])
+function calculate_fxyz(fitn::ITensorNetworkFunction, xs::Vector; kwargs...)
+  return calculate_fxyz(fitn, xs, [i for i in 1:length(xs)]; kwargs...)
 end
 
-function calculate_fx(fitn::ITensorNetworkFunction, x::Number)
+function calculate_fx(fitn::ITensorNetworkFunction, x::Number; kwargs...)
   @assert dimension(fitn) == 1
-  return calculate_fxyz(fitn, [x], [1])
+  return calculate_fxyz(fitn, [x], [1]; kwargs...)
 end
 
 function ITensorNetworks.truncate(fitn::ITensorNetworkFunction; kwargs...)
diff --git a/test/test_itensorfunction.jl b/test/test_itensorfunction.jl
index 09b0887..85cff30 100644
--- a/test/test_itensorfunction.jl
+++ b/test/test_itensorfunction.jl
@@ -5,9 +5,9 @@ using Graphs: SimpleGraph, uniform_tree
 using NamedGraphs: NamedGraph, named_grid, vertices, named_comb_tree, rename_vertices
 using ITensors: siteinds
 using Dictionaries: Dictionary
-using SplitApplyCombine: group
-using Random: seed!
-using Distributions: Uniform
+using Random: Random
+
+Random.seed!(1234)
 
 @testset "test constructor from ITensorNetwork" begin
   L = 10
@@ -22,7 +22,10 @@ using Distributions: Uniform
   @test dimension_vertices(fψ, 1) == vertices(fψ)
   @test dimension(fψ) == 1
 
-  dim_vertices = collect(values(group(v -> first(v) < Int(0.5 * L), vertices(ψ))))
+  dim_vertices = [
+    filter(v -> first(v) < Int(0.5 * L), vertices(ψ)),
+    filter(v -> first(v) >= Int(0.5 * L), vertices(ψ)),
+  ]
   fψ = ITensorNetworkFunction(ψ, dim_vertices)
   @test union(Set(dimension_vertices(fψ, 1)), Set(dimension_vertices(fψ, 2))) ==
     Set(vertices(fψ))
@@ -56,15 +59,16 @@ end
     @testset "test $name in binary" begin
       Lx, Ly = 2, 3
       g = named_comb_tree((2, 3))
-      a = 1.2
-      k = 0.125
+      a = rand()
+      k = rand()
+      c = rand()
 
       s = continuous_siteinds(g)
 
       x = 0.625
-      ψ_fx = net_func(s; k, a)
+      ψ_fx = net_func(s; k, a, c)
       fx_x = calculate_fx(ψ_fx, x)
-      @test func(k * x + a) ≈ fx_x
+      @test c * func(k * x + a) ≈ fx_x
     end
   end
 
@@ -79,33 +83,34 @@ end
     @testset "test $name in trinary" begin
       Lx, Ly = 2, 3
       g = named_comb_tree((2, 3))
-      a = 1.2
-      k = 0.125
-      b = 3
+      a = rand()
+      k = rand()
+      c = rand()
 
       s = continuous_siteinds(g; base=3)
 
       x = (5.0 / 9.0)
-      ψ_fx = net_func(s; k, a)
+      ψ_fx = net_func(s; k, a, c)
       fx_x = calculate_fx(ψ_fx, x)
-      @test func(k * x + a) ≈ fx_x
+      @test c * func(k * x + a) ≈ fx_x
     end
   end
 
   @testset "test tanh" begin
     L = 10
     g = named_grid((L, 1))
-    a = 1.3
-    k = 0.15
+    a = rand()
+    k = rand()
+    c = rand()
     nterms = 50
 
     s = continuous_siteinds(g)
 
     x = 0.625
-    ψ_fx = tanh_itn(s; k, a, nterms)
+    ψ_fx = tanh_itn(s; k, a, c, nterms)
     fx_x = calculate_fx(ψ_fx, x)
 
-    @test tanh(k * x + a) ≈ fx_x
+    @test c * tanh(k * x + a) ≈ fx_x
   end
 
   @testset "test poly" begin
@@ -114,18 +119,19 @@ end
 
     ###Generate a series of random polynomials on random graphs. Evaluate them at random x values"""
     for deg in degrees
-      seed!(1234 * deg)
       g = NamedGraph(SimpleGraph(uniform_tree(L)))
       g = rename_vertices(g, Dict(zip(vertices(g), [(v, 1) for v in vertices(g)])))
       s = continuous_siteinds(g)
+      k = rand()
+      c = rand()
 
-      coeffs = [rand(Uniform(-2, 2)) for i in 1:(deg + 1)]
+      coeffs = [rand() for i in 1:(deg + 1)]
 
       x = 0.875
-      ψ_fx = poly_itn(s, coeffs)
+      ψ_fx = poly_itn(s, coeffs; k, c)
       fx_x = calculate_fx(ψ_fx, x)
 
-      fx_exact = sum([coeffs[i] * (x^(i - 1)) for i in 1:(deg + 1)])
+      fx_exact = c * sum([coeffs[i] * ((k * x)^(i - 1)) for i in 1:(deg + 1)])
       @test fx_x ≈ fx_exact atol = 1e-4
     end
   end
@@ -162,15 +168,16 @@ end
 
   for (name, net_func, func) in funcs
     @testset "test $name" begin
-      a = 1.2
-      k = 0.125
+      a = rand()
+      k = rand()
+      c = rand()
 
-      ψ_fx = net_func(s; k, a, dimension=1)
-      ψ_fy = net_func(s; k, a, dimension=2)
+      ψ_fx = net_func(s; k, a, c, dimension=1)
+      ψ_fy = net_func(s; k, a, c, dimension=2)
 
       ψ_fxy = ψ_fx + ψ_fy
       fxy_xy = calculate_fxyz(ψ_fxy, [x, y], [1, 2])
-      @test func(k * x + a) + func(k * y + a) ≈ fxy_xy
+      @test c * func(k * x + a) + c * func(k * y + a) ≈ fxy_xy
     end
   end
 
@@ -178,17 +185,18 @@ end
   @testset "test tanh" begin
     L = 3
     g = named_grid((L, 2))
-    a = 1.3
-    k = 0.15
-    nterms = 10
+    a = rand()
+    k = rand()
+    c = rand()
+    nterms = 20
     s = continuous_siteinds(g; map_dimension=2)
 
     x, y = 0.625, 0.875
-    ψ_fx = tanh_itn(s; k, a, nterms, dimension=1)
-    ψ_fy = tanh_itn(s; k, a, nterms, dimension=2)
+    ψ_fx = tanh_itn(s; k, a, c, nterms, dimension=1)
+    ψ_fy = tanh_itn(s; k, a, c, nterms, dimension=2)
 
     ψ_fxy = ψ_fx + ψ_fy
-    fxy_xy = calculate_fxyz(ψ_fxy, [x, y], [1, 2])
-    @test tanh(k * x + a) + tanh(k * y + a) ≈ fxy_xy
+    fxy_xy = calculate_fxyz(ψ_fxy, [x, y], [1, 2]; alg="exact")
+    @test c * tanh(k * x + a) + c * tanh(k * y + a) ≈ fxy_xy
   end
 end
diff --git a/test/test_operators.jl b/test/test_operators.jl
index c03b5b7..9a0c28b 100644
--- a/test/test_operators.jl
+++ b/test/test_operators.jl
@@ -29,7 +29,7 @@ using Dictionaries: Dictionary
 
   xs = [0.0, 0.25, 0.625, 0.875, 1.0 - delta]
   for x in xs
-    @test 1.0 + calculate_fx(ψ_fx, x) ≈ 1.0 + sin(2.0 * pi * x) rtol = 1e-3
+    @test 1.0 + calculate_fx(ψ_fx, x; alg="exact") ≈ 1.0 + sin(2.0 * pi * x) rtol = 1e-3
     @test 1.0 + calculate_fx(ψ_f1x, x) ≈ 1.0 + 2.0 * pi * cos(2.0 * pi * x) rtol = 1e-3
     @test 1.0 + calculate_fx(ψ_f2x, x) ≈ 1.0 + -1.0 * (2.0 * pi)^2 * sin(2.0 * pi * x) rtol =
       1e-3
@@ -235,7 +235,7 @@ end
       @test fx_xplus ≈ calculate_fx(ψ_fx_pshift_pbc, x) atol = 1e-8
       @test fx_xplus ≈ calculate_fx(ψ_fx_pshift_neumann, x) atol = 1e-8
     elseif x == 1.0 - 2.0 * delta || x == 1.0 - delta
-      @test calculate_fx(ψ_fx_pshift_dirichlet, x) ≈ 0.0 atol = 1e-8
+      @test calculate_fx(ψ_fx_pshift_dirichlet, x; alg="exact") ≈ 0.0 atol = 1e-8
       @test calculate_fx(ψ_fx_pshift_pbc, x) ≈ calculate_fx(ψ_fx, x + 2.0 * delta - 1.0) atol =
         1e-8
       @test calculate_fx(ψ_fx_pshift_neumann, x) ≈ calculate_fx(ψ_fx, x) atol = 1e-8