diff --git a/src/solve.jl b/src/solve.jl
index 707155a..5495c70 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -549,9 +549,33 @@ function poly_minus_constant_is_greater{T}(Λ::PiecewiseQuadratic{T}, ρ::Quadra
     return true
 end
 
+# Continuous case
+""" `get_transition_costs(g, t, lazy; tol=1e-3)`
+    Return `(ℓ, cost_last)`: the transition cost and end-cost.
+    Slightly unsafe since `tol` is ignored (not needed) on discrete problems
+"""
+function get_transition_costs(g, t, lazy; tol=1e-3)
+    ℓ = lazy ? TransitionCostContinuous{Float64}(g, t, tol) :
+               compute_transition_costs(g, t, tol)
+    # Continouous case, no cost at endpoint
+    cost_last = zero(QuadraticPolynomial{Float64})
+    return ℓ, cost_last
+end
+
+# Discrete case, tol is not used here, but in signature to enable dispatch
+function get_transition_costs(g::AbstractArray, t, lazy; tol=1e-3)
+    ℓ = lazy ? TransitionCostDiscrete{Float64}(g, t) :
+               compute_discrete_transition_costs(g, t)
+    if t[1] != 1 || t[end] != length(g)
+        warn("In fit_pwl_constrained: The supplied grid t only covers the range ($(t[1]),$(t[end])) while the range of indices for g is (1,$(length(g))). No costs will be considered outside the range of t.")
+    end
+    # Discrete case, so cost at endpoint is quadratic
+    cost_last = QuadraticPolynomial(1.0, -2*g[t[end]], g[t[end]]^2)
+    return ℓ, cost_last
+end
 
 
-## Convenience function
+## Public interface
 """
     I, Y, v = fit_pwl_regularized(g::AbstractArray, ζ; t=1:length(g), lazy=true)
     I, Y, v = fit_pwl_regularized(g, t, ζ, tol=1e-3; lazy=true)
@@ -571,19 +595,13 @@ i.e. so that `f[t[I]] .= Y`.
 `tol` specifies the relative tolerance sent to `quadg` kused when calculating the integrals (continuous case).
 """
 function  fit_pwl_regularized(g::AbstractArray, ζ; t=1:length(g), lazy=true)
-    ℓ = lazy ? TransitionCostDiscrete{Float64}(g, t) :
-               compute_discrete_transition_costs(g, t)
-    # Discrete case, so cost at endpoint is quadratic
-    cost_last = QuadraticPolynomial(1.0, -2*g[end], g[end]^2)
+    ℓ, cost_last = get_transition_costs(g, t, lazy)
     fit_pwl_regularized_internal(ℓ, cost_last, ζ)
 end
 
 # Continuous function
 function  fit_pwl_regularized(g, t, ζ, tol=1e-3; lazy=true)
-    ℓ = lazy ? TransitionCostContinuous{Float64}(g, t, tol) :
-               compute_transition_costs(g, t, tol)
-    # Continouous case, no cost at endpoint
-    cost_last = zero(QuadraticPolynomial{Float64})
+    ℓ, cost_last = get_transition_costs(g, t, lazy, tol=tol)
     fit_pwl_regularized_internal(ℓ, cost_last, ζ)
 end
 
@@ -615,26 +633,22 @@ i.e. so that `f[t[I]] .= Y`.
 `tol` specifies the relative tolerance sent to `quadg` kused when calculating the integrals (continuous case).
 """
 function  fit_pwl_constrained(g::AbstractArray, M; t=1:length(g), lazy=false)
-    ℓ = lazy ? TransitionCostDiscrete{Float64}(g, t) :
-               compute_discrete_transition_costs(g, t)
-    cost_last = QuadraticPolynomial(1.0, -2*g[end], g[end]^2)
+    ℓ, cost_last = get_transition_costs(g, t, lazy)
     fit_pwl_constrained_internal(ℓ, cost_last, M)
 end
 
 # Continuous function
 function  fit_pwl_constrained(g, t, M, tol=1e-3; lazy=false)
-    ℓ = lazy ? TransitionCostContinuous{Float64}(g, t, tol) :
-               compute_transition_costs(g, t, tol)
-    cost_last = zero(QuadraticPolynomial{Float64})
+    ℓ, cost_last = get_transition_costs(g, t, lazy, tol=tol)
     fit_pwl_constrained_internal(ℓ, cost_last, M)
 end
 
-function  fit_pwl_constrained_internal(ℓ, cost_last, M)
+function  fit_pwl_constrained_internal{T}(ℓ::AbstractTransitionCost{T}, cost_last, M)
     Λ = pwq_dp_constrained(ℓ, cost_last, M);
 
     Ivec = Vector{Vector{Int}}(M)
-    Yvec = Vector{Vector{Float64}}(M)
-    fvec = Vector{Float64}(M)
+    Yvec = Vector{Vector{T}}(M)
+    fvec = Vector{T}(M)
 
     for m=1:M
         Ivec[m], Yvec[m], fvec[m] = recover_solution(Λ[m, 1], ℓ, cost_last)
diff --git a/src/types/QuadraticPolynomial.jl b/src/types/QuadraticPolynomial.jl
index 6369c77..b01efd9 100644
--- a/src/types/QuadraticPolynomial.jl
+++ b/src/types/QuadraticPolynomial.jl
@@ -20,6 +20,7 @@ QuadraticPolynomial{T}() where {T} = new()
 end
 
 QuadraticPolynomial{T}(a::T,b::T,c::T) = QuadraticPolynomial{T}(a,b,c)
+QuadraticPolynomial(a,b,c) = QuadraticPolynomial(promote(a,b,c)...)
 
 function QuadraticPolynomial(x::Vector)
     QuadraticPolynomial(x[3], x[2], x[1])
diff --git a/test/auxilliary_test_fcns.jl b/test/auxilliary_test_fcns.jl
index 451c054..114a5ee 100644
--- a/test/auxilliary_test_fcns.jl
+++ b/test/auxilliary_test_fcns.jl
@@ -2,25 +2,32 @@
 circle_segment(N) = sin.(acos.(linspace(-1, 1, N)))
 linear_trend(N) = (0:N-1)
 
-# Finds the solution to the cardinality constrained problem for
-# m=1:M segments using brute force
-function brute_force_multi(g, t, M; tol=1e-3)
+# function get_transition_costs(g, t; tol=1e-3)
+#     ℓ = TransitionCostContinuous{Float64}(g, t, tol=tol)
+#     V_N = zero(QuadraticPolynomial{Float64})
+#     return ℓ, V_N
+# end
+# function get_transition_costs(g::AbstractArray, t; tol=1e-3)
+#     ℓ = compute_discrete_transition_costs(g, t)
+#     V_N = QuadraticPolynomial(1.0, -2*g[t[end]], g[t[end]]^2)
+#     return ℓ, V_N
+# end
+
+
+
+""" brute_force_multi(g, M, t=1:length(g); tol=1e-3)
+    Finds the solution to the cardinality constrained problem for
+    m=1:M segments using brute force with possible breakpoints on grid t.
+    `t` has to be supplied if `g` is not `AbstractArray`
+"""
+function brute_force_multi(g, M, t=1:length(g); tol=1e-3)
 
     if length(g) <= M
         warn("Specified M too large. Setting M=length(g)-1.")
         M = length(g)-1
     end
 
-    ℓ = []
-    V_N = []
-
-    if typeof(g) <: AbstractArray
-        ℓ = compute_discrete_transition_costs(g)
-        V_N = QuadraticPolynomial(1.0, -2*g[end], g[end]^2)
-    else
-        ℓ = TransitionCostContinuous{Float64}(g, t, tol=tol)
-        V_N = zero(QuadraticPolynomial{Float64})
-    end
+    ℓ, V_N = DynamicApproximations.get_transition_costs(g, t, true, tol=tol)
 
     I_vec = Vector{Vector{Int64}}(M)
     Y_vec = Vector{Vector{Float64}}(M)
@@ -28,7 +35,7 @@ function brute_force_multi(g, t, M; tol=1e-3)
 
     for m=1:M
         (I_bf, Y_bf, f_bf) = brute_force_search(ℓ, V_N, m)
-        I_vec[m] = I_bf
+        I_vec[m] = t[I_bf]
         Y_vec[m] = Y_bf
         f_vec[m] = f_bf
     end
diff --git a/test/brute_force_search.jl b/test/brute_force_search.jl
index 47a8281..a3f64a4 100644
--- a/test/brute_force_search.jl
+++ b/test/brute_force_search.jl
@@ -16,7 +16,7 @@ for problem_fcn in ["discontinuous1",
     g, ζ_vec, I_sols, f_sols = @eval $(Symbol(problem_fcn))()
 
     M = length(I_sols)
-    I_vec, _, f_vec = brute_force_multi(g, -1, M)
+    I_vec, _, f_vec = brute_force_multi(g, M)
 
     @testset "Data set: $problem_fcn, brute_force_search m=$m" for m in 1:length(f_sols)
         @test f_vec[m] ≈ f_sols[m]  atol=1e-10
diff --git a/test/generate_problem_files.jl b/test/generate_problem_files.jl
index 4e9fcda..79390d3 100644
--- a/test/generate_problem_files.jl
+++ b/test/generate_problem_files.jl
@@ -61,7 +61,7 @@ problem_data[10] = ("super_exponential5",
 for (problem_name, g, M) in problem_data[10:end]
 
 
-    @time I_vec, Y_vec, f_vec = brute_force_multi(g, 1, M)
+    @time I_vec, Y_vec, f_vec = brute_force_multi(g, M)
 
     # It should not be beneficial to use more segments than has been computed ...
     ζ_vec = logspace(log10(f_vec[2]), log10(0.8*(max(f_vec[end-1] - f_vec[end], f_vec[end-1]))), 10)
diff --git a/test/runtests.jl b/test/runtests.jl
index f58a26c..dfa3fc1 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -12,6 +12,7 @@ tests = [
     "exact_tests.jl",
     "fit_pwl_constrained.jl",
     "fit_pwl_regularized.jl",
+    "subsample_test.jl",
     "examples.jl"]
 
 @testset "All Tests" begin
diff --git a/test/subsample_test.jl b/test/subsample_test.jl
new file mode 100644
index 0000000..3748c9d
--- /dev/null
+++ b/test/subsample_test.jl
@@ -0,0 +1,65 @@
+using DynamicApproximations, Base.Test
+
+include(joinpath(Pkg.dir("DynamicApproximations"),"test","auxilliary_test_fcns.jl"))
+
+srand(1234)
+
+# Brute force up to nbr segments:
+M = 7
+# Total number of points:
+N = 40
+# Number of internal gridpoints:
+n = 12
+
+#Function
+g = cumsum(randn(N).^2)
+
+t_grids = Vector{Int}[1:N, [[1;sort!(shuffle(2:N-1)[1:n]);N] for i = 1:10]...]
+@testset "Testing subgrid, regularized, with endpoints, i=$i" for (i,t_grid) in enumerate(t_grids)
+    I_sol, Y_sol, f_sol = brute_force_multi(g, M, t_grid, tol=1e-3)
+
+    I, Y, f = fit_pwl_constrained(g, M, t=t_grid)
+    I_t = getindex.([t_grid], I)
+    @test I_t == I_sol
+    @test Y ≈ Y_sol atol = 1e-10 # Crazy that this is possible
+    @test f ≈ f_sol atol = 1e-10
+
+    ζ_vec = logspace(log10(f_sol[2]), log10(0.8*(max(f_sol[end-1] - f_sol[end], f_sol[end-1]))), 3)
+    for ζ in ζ_vec
+        I2, Y2, f2 = fit_pwl_regularized(g, ζ, t=t_grid)
+        I2_t = t_grid[I2]
+        #Make sure regularized solution is in constrained solution
+        idx = find([I2_t] .== I_t)
+        @test length(idx) == 1
+        if length(idx) == 1
+            @test Y2 ≈ Y_sol[idx[1]] atol = 1e-10
+            @test f2 ≈ f_sol[idx[1]] atol = 1e-10
+        end
+    end
+end
+
+println("Testing sub-grid that is not covering the full range of input, expect warnings.")
+# Test when endpoints are not included
+t_grids = Vector{Int}[1:N, [sort!(shuffle(2:N-1)[1:n]) for i = 1:10]...]
+@testset "Testing subgrid, constrained, without endpoints, i=$i" for (i, t_grid) in enumerate(t_grids)
+    I_sol, Y_sol, f_sol = brute_force_multi(g, M-2, t_grid, tol=1e-3)
+
+    I, Y, f = fit_pwl_constrained(g, M-2, t=t_grid)
+    I_t = getindex.([t_grid], I)
+    @test I_t == I_sol
+    @test Y ≈ Y_sol atol = 1e-10 # Crazy that this is possible
+    @test f ≈ f_sol atol = 1e-10
+
+    ζ_vec = logspace(log10(f_sol[2]), log10(0.8*(max(f_sol[end-1] - f_sol[end], f_sol[end-1]))), 3)
+    for ζ in ζ_vec
+        I2, Y2, f2 = fit_pwl_regularized(g, ζ, t=t_grid)
+        I2_t = t_grid[I2]
+        #Make sure regularized solution is in constrained solution
+        idx = find([I2_t] .== I_t)
+        @test length(idx) == 1
+        if length(idx) == 1
+            @test Y2 ≈ Y_sol[idx[1]] atol = 1e-10
+            @test f2 ≈ f_sol[idx[1]] atol = 1e-10
+        end
+    end
+end