diff --git a/src/ACEfit.jl b/src/ACEfit.jl
index 86c6841..acd3e14 100644
--- a/src/ACEfit.jl
+++ b/src/ACEfit.jl
@@ -4,6 +4,7 @@ include("bayesianlinear.jl")
 include("data.jl")
 include("assemble.jl")
 include("solvers.jl")
+include("asp.jl")
 include("fit.jl")
 
 end
diff --git a/src/asp.jl b/src/asp.jl
new file mode 100644
index 0000000..3de00e9
--- /dev/null
+++ b/src/asp.jl
@@ -0,0 +1,157 @@
+
+@doc raw"""
+
+`ASP` : Active Set Pursuit solver
+
+Solves the lasso optimization problem. 
+```math 
+\max_{y} \left( b^T y - \frac{1}{2} λ y^T y \right)
+```
+subject to
+```math
+   \|A^T y\|_{\infty} \leq 1.
+```
+
+### Constructor Keyword arguments 
+```julia
+ACEfit.ASP(; P = I, select = (:byerror, 1.0), 
+            params...)
+``` 
+
+* `select` : Selection criterion for the final solution (required) 
+    * `:final` : final solution (largest computed basis)
+    * `(:byerror, q)` : solution with error within `q` times the minimum error 
+       along the path; if training error is used and `q == 1.0`, then this is 
+       equivalent to to `:final`.
+    * `(:bysize, n)` : best solution with at most `n` non-zero features; if 
+       training error is used, then it will be the solution with exactly `n` 
+       non-zero features. 
+* `P = I` : prior / regularizer (optional)
+
+The remaining kwarguments to `ASP` are parameters for the ASP homotopy solver. 
+
+* `actMax` : Maximum number of active constraints.
+* `min_lambda` : Minimum value for `λ`.  (defaults to 0)
+* `loglevel` : Logging level.
+* `itnMax` : Maximum number of iterations.
+
+### Extended syntax for `solve` 
+
+```julia
+solve(solver::ASP, A, y, Aval=A, yval=y)
+```
+* `A` : `m`-by-`n` design matrix. (required)
+* `b` : `m`-vector. (required)
+* `Aval = nothing` : `p`-by-`n` validation matrix (only for `validate` mode).
+* `bval = nothing` : `p`- validation vector (only for `validate` mode).
+
+If independent `Aval` and `yval` are provided (instead of detaults `A, y`), 
+then the solver will use this separate validation set instead of the training
+set to select the best solution along the model path. 
+"""
+struct ASP
+    P
+    select
+    params
+end
+
+function ASP(; P = I, select, mode=:train, params...)
+    return ASP(P, select, params)
+end
+
+function solve(solver::ASP, A, y, Aval=A, yval=y)
+    # Apply preconditioning
+    AP = A / solver.P
+    
+    tracer = asp_homotopy(AP, y; solver.params...)
+    q = length(tracer)
+    new_tracer = Vector{NamedTuple{(:solution, :λ), Tuple{Any, Any}}}(undef, q)
+
+    for i in 1:q
+        new_tracer[i] = (solution = solver.P \ tracer[i][1], λ = tracer[i][2])
+    end
+
+    xs, in = select_solution(new_tracer, solver, Aval, yval)
+
+   #  println("done.")
+    return Dict(   "C" => xs, 
+                "path" => new_tracer, 
+                "nnzs" => length((new_tracer[in][:solution]).nzind) )
+end
+
+
+function select_solution(tracer, solver, A, y)
+    if solver.select == :final
+      criterion = :final 
+    else
+      criterion, p = solver.select
+    end
+    
+    if criterion == :final
+        return tracer[end][:solution], length(tracer)
+    end 
+
+    if criterion == :byerror
+        maxind = length(tracer)
+        threshold = p 
+    elseif criterion == :bysize
+        maxind = findfirst(t -> length((t[:solution]).nzind) > p, 
+                         tracer) - 1
+        threshold = 1.0                          
+    else 
+        error("Unknown selection criterion: $criterion")
+    end
+
+    errors = [ norm(A * t[:solution] - y) for t in tracer[1:maxind] ]
+    min_error = minimum(errors)
+    for (i, error) in enumerate(errors)
+        if error <= threshold * min_error
+            return tracer[i][:solution], i
+        end
+    end
+
+    error("selection failed for unknown reasons; please file an issue with a MWE to reproduce this error.")
+end
+
+
+#=
+function select_smart(tracer, solver, Aval, yval)
+
+    best_metric = Inf
+    best_iteration = 0
+    validation_metric = 0
+    q = length(tracer)
+    errors = [norm(Aval * t[:solution] - yval) for t in tracer]
+    nnzss = [(t[:solution]).nzind for t in tracer]
+    best_iteration = argmin(errors)
+    validation_metric = errors[best_iteration]
+    validation_end = norm(Aval * tracer[end][:solution] - yval)
+
+    if validation_end < validation_metric #make sure to check the last one too in case q<<100
+        best_iteration = q
+    end
+
+    criterion, threshold = solver.select
+    
+    if criterion == :val
+        return tracer[best_iteration][:solution], best_iteration
+
+    elseif criterion == :byerror
+        for (i, error) in enumerate(errors)
+            if error <= threshold * validation_metric
+                return tracer[i][:solution], i
+            end
+        end
+
+    elseif criterion == :bysize
+        first_index = findfirst(sublist -> threshold in sublist, nnzss)
+        relevant_errors = errors[1:first_index - 1] 
+        min_error = minimum(relevant_errors)
+        min_error_index = findfirst(==(min_error), relevant_errors)
+        return tracer[min_error_index][:solution], min_error_index
+
+    else
+        @error("Unknown selection criterion: $criterion")
+    end
+end
+=#
\ No newline at end of file
diff --git a/src/solvers.jl b/src/solvers.jl
index fc35934..9867271 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -196,95 +196,3 @@ function solve(solver::TruncatedSVD, A, y)
     return Dict{String, Any}("C" => solver.P \ θP)
 end
 
-
-@doc raw"""
-`struct ASP` : Active Set Pursuit sparse solver
-    solves the following optimization problem using the homotopy approach:
-
-    ```math 
-    \max_{y} \left( b^T y - \frac{1}{2} λ y^T y \right)
-    ```
-        subject to
-        
-    ```math
-        \|A^T y\|_{\infty} \leq 1.
-    ```
-
-    * Input
-    * `A` : `m`-by-`n` explicit matrix or linear operator.
-    * `b` : `m`-vector.
-
-    * Solver parameters
-    * `min_lambda` : Minimum value for `λ`. Defaults to zero if not provided.
-    * `loglevel` : Logging level.
-    * `itnMax` : Maximum number of iterations.
-    * `actMax` : Maximum number of active constraints.
-
-    Constructor
-    ```julia
-    ACEfit.ASP(; P = I, select, params)
-    ``` 
-    where 
-    - `P` : right-preconditioner / tychonov operator
-    - `select`: Selection mode for the final solution.
-    - `(:byerror, th)`: Selects the smallest active set fit within a factor `th` of the smallest fit error.
-    - `(:final, nothing)`: Returns the final iterate.
-    - `params`: The solver parameters, passed as named arguments.
-"""
-struct ASP
-    P::Any
-    select::Tuple
-    params::NamedTuple
-end
-
-function ASP(; P = I, select, params...)
-    params_tuple = NamedTuple(params)
-    return ASP(P, select, params_tuple)
-end
-
-function solve(solver::ASP, A, y)
-    # Apply preconditioning
-    AP = A / solver.P
-    
-    tracer = asp_homotopy(AP, y; solver.params[1]...)
-    
-    new_tracer = Vector{NamedTuple{(:solution, :λ), Tuple{Any, Any}}}(undef, length(tracer))
-
-    for i in 1:length(tracer)
-        new_tracer[i] = (solution = solver.P \ tracer[i][1], λ = tracer[i][2])
-    end
-
-    # Select the final solution based on the criterion
-    xs, in = select_solution(new_tracer, solver, A, y)
-    
-    println("done.")
-    return Dict("C" => xs, "path" => new_tracer, "nnzs" => length((tracer[in][1]).nzind) )
-end
-
-function select_solution(tracer, solver, A, y)
-    criterion, threshold = solver.select
-    
-    if criterion == :final
-        return tracer[end][1], length(tracer)
-
-    elseif criterion == :byerror
-        errors = [norm(A * t[1] - y) for t in tracer]
-        min_error = minimum(errors)
-        
-        # Find the solution with the smallest error within the threshold
-        for (i, error) in enumerate(errors)
-            if error <= threshold * min_error
-                return tracer[i][1], i
-            end
-        end
-    elseif criterion == :bysize
-        for i in 1:length(tracer)
-            if length((tracer[i][1]).nzind) == threshold
-                return tracer[i][1], i
-            end
-        end
-    else
-        @error("Unknown selection criterion: $criterion")
-    end
-end
-
diff --git a/test/Project.toml b/test/Project.toml
index 2284da9..53ae33c 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,4 +4,5 @@ MLJ = "add582a8-e3ab-11e8-2d5e-e98b27df1bc7"
 MLJLinearModels = "6ee0df7b-362f-4a72-a706-9e79364fb692"
 MLJScikitLearnInterface = "5ae90465-5518-4432-b9d2-8a1def2f0cab"
 PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/test_linearsolvers.jl b/test/test_linearsolvers.jl
index 4535ba9..234be82 100644
--- a/test/test_linearsolvers.jl
+++ b/test/test_linearsolvers.jl
@@ -1,9 +1,14 @@
 
 using ACEfit
-using LinearAlgebra
+using LinearAlgebra, Random, Test 
+using Random
 using PythonCall
 
+##
+
 @info("Test Solver on overdetermined system")
+
+Random.seed!(1234)
 Nobs = 10_000
 Nfeat = 100
 A1 = randn(Nobs, Nfeat) / sqrt(Nobs)
@@ -11,9 +16,12 @@ U, S1, V = svd(A1)
 S = 1e-4 .+ ((S1 .- S1[end]) / (S1[1] - S1[end])).^2
 A = U * Diagonal(S) * V'
 c_ref = randn(Nfeat)
-y = A * c_ref + 1e-3 * randn(Nobs) / sqrt(Nobs)
+epsn = 1e-5 
+y = A * c_ref + epsn * randn(Nobs) / sqrt(Nobs)
 P = Diagonal(1.0 .+ rand(Nfeat))
 
+##
+
 @info(" ... QR")
 solver = ACEfit.QR()
 results = ACEfit.solve(solver, A, y)
@@ -21,6 +29,10 @@ C = results["C"]
 @show norm(A * C - y)
 @show norm(C)
 @show norm(C - c_ref)
+@test norm(A * C - y) < 10 * epsn 
+@test norm(C - c_ref) < 100 * epsn 
+
+##
 
 @info(" ... regularised QR, λ = 1e-5")
 solver = ACEfit.QR(lambda = 1e-5, P = P)
@@ -29,6 +41,11 @@ C = results["C"]
 @show norm(A * C - y)
 @show norm(C)
 @show norm(C - c_ref)
+@test norm(A * C - y) < 10 * epsn 
+@test norm(C - c_ref) < 1000 * epsn 
+
+
+##
 
 @info(" ... regularised QR, λ = 1e-2")
 solver = ACEfit.QR(lambda = 1e-2, P = P)
@@ -38,6 +55,11 @@ C = results["C"]
 @show norm(C)
 @show norm(C - c_ref)
 
+@test norm(A * C - y) < 1
+@test norm(C - c_ref) < 10
+
+##
+
 @info(" ... RRQR, rtol = 1e-15")
 solver = ACEfit.RRQR(rtol = 1e-15, P = P)
 results = ACEfit.solve(solver, A, y)
@@ -46,6 +68,10 @@ C = results["C"]
 @show norm(C)
 @show norm(C - c_ref)
 
+@test norm(A * C - y) < 10 * epsn 
+@test norm(C - c_ref) < 100 * epsn 
+
+##
 
 @info(" ... RRQR, rtol = 1e-5")
 solver = ACEfit.RRQR(rtol = 1e-5, P = P)
@@ -54,6 +80,10 @@ C = results["C"]
 @show norm(A * C - y)
 @show norm(C)
 @show norm(C - c_ref)
+@test norm(A * C - y) < 10 * epsn 
+@test norm(C - c_ref) < 100 * epsn 
+
+##
 
 @info(" ... RRQR, rtol = 1e-3")
 solver = ACEfit.RRQR(rtol = 1e-3, P = P)
@@ -63,6 +93,11 @@ C = results["C"]
 @show norm(C)
 @show norm(C - c_ref)
 
+@test norm(A * C - y) < 1
+@test norm(C - c_ref) < 1
+
+##
+
 @info(" ... LSQR")
 solver = ACEfit.LSQR(damp = 0, atol = 1e-6)
 results = ACEfit.solve(solver, A, y)
@@ -71,6 +106,11 @@ C = results["C"]
 @show norm(C)
 @show norm(C - c_ref)
 
+@test norm(A * C - y) < 10 * epsn 
+@test norm(C - c_ref) < 1
+
+##
+
 @info(" ... SKLEARN_BRR")
 solver = ACEfit.SKLEARN_BRR()
 results = ACEfit.solve(solver, A, y)
@@ -79,6 +119,8 @@ C = results["C"]
 @show norm(C)
 @show norm(C - c_ref)
 
+##
+
 @info(" ... SKLEARN_ARD")
 solver = ACEfit.SKLEARN_ARD()
 results = ACEfit.solve(solver, A, y)
@@ -87,6 +129,8 @@ C = results["C"]
 @show norm(C)
 @show norm(C - c_ref)
 
+##
+
 @info(" ... BLR")
 solver = ACEfit.BLR()
 results = ACEfit.solve(solver, A, y)
@@ -95,6 +139,11 @@ C = results["C"]
 @show norm(C)
 @show norm(C - c_ref)
 
+@test norm(A * C - y) < 10 * epsn 
+@test norm(C - c_ref) < 1
+
+##
+
 @info(" ... TruncatedSVD(; rtol = 1e-5)")
 solver = ACEfit.TruncatedSVD(; rtol = 1e-5)
 results = ACEfit.solve(solver, A, y)
@@ -103,40 +152,61 @@ C = results["C"]
 @show norm(C)
 @show norm(C - c_ref)
 
-@info(" ... TruncatedSVD(; rtol = 1e-4)")
-solver = ACEfit.TruncatedSVD(; rtol=1e-4)
-results = ACEfit.solve(solver, A, y)
-C = results["C"]
-@show norm(A * C - y)
-@show norm(C)
-@show norm(C - c_ref)
+@test norm(A * C - y) < 10 * epsn 
+@test norm(C - c_ref) < 100 * epsn 
 
-@info(" ... ASP_homotopy selected by error")
-solver = ACEfit.ASP(P = P, select = (:byerror,1.5), params = (loglevel=0, traceFlag=true))
-results = ACEfit.solve(solver, A, y)
-C = results["C"]
-full_path = results["path"]
-@show results["nnzs"]
-@show norm(A * C - y)
-@show norm(C)
-@show norm(C - c_ref)
+##
 
-@info(" ... ASP_homotopy selected by size")
-solver = ACEfit.ASP(P = P, select = (:bysize,50), params = (loglevel=0, traceFlag=true))
+@info(" ... TruncatedSVD(; rtol = 1e-4)")
+solver = ACEfit.TruncatedSVD(; rtol=1e-4)
 results = ACEfit.solve(solver, A, y)
 C = results["C"]
-full_path = results["path"]
-@show results["nnzs"]
 @show norm(A * C - y)
 @show norm(C)
 @show norm(C - c_ref)
 
-@info(" ... ASP_homotopy final solution")
-solver = ACEfit.ASP(P = P, select = (:final,nothing), params = (loglevel=0, traceFlag=true))
-results = ACEfit.solve(solver, A, y)
-C = results["C"]
-full_path = results["path"]
-@show results["nnzs"]
-@show norm(A * C - y)
-@show norm(C)
-@show norm(C - c_ref)
+@test norm(A * C - y) < 10 * epsn 
+@test norm(C - c_ref) < 1
+
+##
+
+@info(" ... ASP")
+shuffled_indices = shuffle(1:length(y))
+train_indices = shuffled_indices[1:round(Int, 0.85 * length(y))]
+val_indices = shuffled_indices[round(Int, 0.85 * length(y)) + 1:end]
+At = A[train_indices,:]
+Av = A[val_indices,:]
+yt = y[train_indices]
+yv = y[val_indices]
+
+for (select, tolr, tolc) in [ (:final, 10*epsn, 1), 
+                             ( (:byerror,1.3), 10*epsn, 1), 
+                            ( (:bysize,73), 1, 10) ]
+    @show select 
+    local solver, results, C 
+    solver = ACEfit.ASP(P=I, select = select, loglevel=0, traceFlag=true)
+    # without validation 
+    results = ACEfit.solve(solver, A, y)
+    C = results["C"]
+    full_path = results["path"]
+    @show results["nnzs"]
+    @show norm(A * C - y)
+    @show norm(C)
+    @show norm(C - c_ref)
+
+    @test norm(A * C - y) < tolr
+    @test norm(C - c_ref) < tolc
+    
+
+    # with validation 
+    results = ACEfit.solve(solver, At, yt, Av, yv)
+    C = results["C"]
+    full_path = results["path"]
+    @show results["nnzs"]
+    @show norm(Av * C - yv)
+    @show norm(C)
+    @show norm(C - c_ref)
+
+    @test norm(Av * C - yv) < tolr
+    @test norm(C - c_ref) < tolc
+end