diff --git a/Project.toml b/Project.toml
index f6d9b62b..8454142f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,7 +25,6 @@ StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
 
 [compat]
 ADTypes = "1.20.0"
-Aqua = "0.8.12"
 CommonSolve = "0.2.6"
 ConcreteStructs = "0.2.3"
 DifferentiationInterface = "0.7.13"
@@ -41,18 +40,19 @@ NonlinearSolveBase = "2.25"
 NonlinearSolveFirstOrder = "2"
 PrecompileTools = "1"
 RecursiveArrayTools = "3.33.0, 4.0"
+SafeTestsets = "0.1"
 SciMLBase = "2.90, 3.1"
+SciMLTesting = "1"
 Setfield = "1.1.2"
 StatsAPI = "1.8.0"
 Test = "1.10"
-TestItemRunner = "1.1.0"
 julia = "1.10"
 
 [extras]
-Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["Test", "TestItemRunner", "Aqua", "ExplicitImports"]
+test = ["Test", "SafeTestsets", "SciMLTesting", "ExplicitImports"]
diff --git a/test/qa.jl b/test/explicit_imports.jl
similarity index 58%
rename from test/qa.jl
rename to test/explicit_imports.jl
index d550d168..ff17a15a 100644
--- a/test/qa.jl
+++ b/test/explicit_imports.jl
@@ -1,12 +1,9 @@
-@testitem "Aqua tests" tags = [:qa, :nopre] begin
-    using Aqua
-    Aqua.test_all(CurveFit)
-end
-
-@testitem "Explicit Imports" tags = [:qa] begin
-    using ExplicitImports
-    using StatsAPI
+using CurveFit
+using Test
+using ExplicitImports
+using StatsAPI
 
+@testset "Explicit Imports" begin
     @test check_no_implicit_imports(CurveFit) === nothing
     @test check_no_stale_explicit_imports(CurveFit) === nothing
     @test check_no_self_qualified_accesses(CurveFit) === nothing
diff --git a/test/expsumfit.jl b/test/expsumfit.jl
index c5cc1969..deee19fe 100644
--- a/test/expsumfit.jl
+++ b/test/expsumfit.jl
@@ -1,105 +1,7 @@
-@testitem "ExpSumFit: Integration rules" tags = [:expsumfit] begin
-    # Trapezoidal rule
-    @test CurveFit.__calc_integral_rules(Rational{Int64}, 1:1; m = 1) == [1 // 2 1 // 2]
+using CurveFit
+using Test
 
-    # Simpson's first (1/3) rule
-    @test CurveFit.__calc_integral_rules(Rational{Int64}, 1:1; m = 2) ==
-        [1 // 3 4 // 3 1 // 3]
-    @test CurveFit.__calc_integral_rules(Rational{Int64}, 2:2; m = 2) ==
-        [2 // 3 4 // 3 0 // 1]
-    @test CurveFit.__calc_integral_rules(Rational{Int64}, 3:3; m = 2) ==
-        [3 // 5 4 // 5 -1 // 15]
-
-    # Simpson's second (3/8) rule
-    @test CurveFit.__calc_integral_rules(Rational{Int64}, 1:1; m = 3) ==
-        [3 // 8 9 // 8 9 // 8 3 // 8]
-    @test CurveFit.__calc_integral_rules(Rational{Int64}, 2:2; m = 3) ==
-        [39 // 40 27 // 10 27 // 40 3 // 20]
-end
-
-@testitem "ExpSumFit: Cumulative integrals" tags = [:expsumfit] begin
-    n = 4
-
-    x = collect(0:0.4:10)
-    y = @. 1 + sin(x)
-    cumints_analytic = @. [(x + 1 - cos(x)) (1 / 2 * x^2 + x - sin(x)) (
-        1 / 6 * x^3 +
-            x^2 / 2 +
-            cos(x) - 1
-    ) (
-        1 /
-            24 *
-            x^4 +
-            x^3 /
-            6 +
-            sin(x) -
-            x
-    )]
-
-    @testset "m = 1 & n = 4" begin
-        coeffs = CurveFit.__calc_integral_rules(Float64, 1:4; m = 1)
-
-        len = length(x)
-        nY, mY = 1 + (len - 1) ÷ 1, 2 * n + 1
-
-        Y = Matrix{Float64}(undef, nY, mY)
-        S = Matrix{Float64}(undef, nY, 4 - 1)
-
-        CurveFit.__cumulative_integrals!(Y, S, coeffs, x, y, 4, 1)
-
-        @test cumints_analytic[:, 1] ≈ Y[:, 1] rtol = 3.0e-3
-        @test cumints_analytic[:, 2] ≈ Y[:, 2] rtol = 3.0e-3
-        @test cumints_analytic[:, 3] ≈ Y[:, 3] rtol = 4.0e-3
-        @test cumints_analytic[:, 4] ≈ Y[:, 4] rtol = 4.0e-3
-    end
-
-    @testset "m = 2 & n = 4" begin
-        coeffs = CurveFit.__calc_integral_rules(Float64, 1:4; m = 2)
-
-        len = length(x)
-        nY, mY = 1 + (len - 1) ÷ 2, 2 * n + 1
-
-        Y = Matrix{Float64}(undef, nY, mY)
-        S = Matrix{Float64}(undef, nY, 4 - 1)
-
-        CurveFit.__cumulative_integrals!(Y, S, coeffs, x, y, 4, 2)
-
-        @test cumints_analytic[1:2:end, 1] ≈ Y[:, 1] rtol = 3.0e-5
-        @test cumints_analytic[1:2:end, 2] ≈ Y[:, 2] rtol = 4.0e-5
-        @test cumints_analytic[1:2:end, 3] ≈ Y[:, 3] rtol = 5.0e-5
-        @test cumints_analytic[1:2:end, 4] ≈ Y[:, 4] rtol = 6.0e-5
-    end
-
-    x = collect(range(0, stop = 2, length = 13))
-    y = @. exp(x)
-    cumints_analytic = @. [(exp(x) - 1) (exp(x) - 1 - x) (exp(x) - 1 - x - x^2 / 2) (
-        exp(x) -
-            1 - x -
-            x^2 /
-            2 -
-            x^3 /
-            6
-    )]
-
-    @testset "m = 2 & n = 4" begin
-        coeffs = CurveFit.__calc_integral_rules(Float64, 1:4; m = 2)
-
-        len = length(x)
-        nY, mY = 1 + (len - 1) ÷ 2, 2 * n + 1
-
-        Y = Matrix{Float64}(undef, nY, mY)
-        S = Matrix{Float64}(undef, nY, 4 - 1)
-
-        CurveFit.__cumulative_integrals!(Y, S, coeffs, x, y, 4, 2)
-
-        @test cumints_analytic[1:2:end, 1] ≈ Y[:, 1] rtol = 5.0e-6
-        @test cumints_analytic[1:2:end, 2] ≈ Y[:, 2] rtol = 3.0e-5
-        @test cumints_analytic[1:2:end, 3] ≈ Y[:, 3] rtol = 3.0e-5
-        @test cumints_analytic[1:2:end, 4] ≈ Y[:, 4] rtol = 4.0e-5
-    end
-end
-
-@testitem "ExpSumFit" tags = [:expsumfit] begin
+@testset "ExpSumFit" begin
     x = collect(0.02:0.02:1.5)
     y = @. 5 * exp(0.5 * x) + 4 * exp(-3 * x) + 2 * exp(-2 * x) - 3 * exp(0.15 * x)
 
diff --git a/test/expsumfit_cumulative_integrals.jl b/test/expsumfit_cumulative_integrals.jl
new file mode 100644
index 00000000..83cde84d
--- /dev/null
+++ b/test/expsumfit_cumulative_integrals.jl
@@ -0,0 +1,84 @@
+using CurveFit
+using Test
+
+@testset "ExpSumFit: Cumulative integrals" begin
+    n = 4
+
+    x = collect(0:0.4:10)
+    y = @. 1 + sin(x)
+    cumints_analytic = @. [(x + 1 - cos(x)) (1 / 2 * x^2 + x - sin(x)) (
+        1 / 6 * x^3 +
+            x^2 / 2 +
+            cos(x) - 1
+    ) (
+        1 /
+            24 *
+            x^4 +
+            x^3 /
+            6 +
+            sin(x) -
+            x
+    )]
+
+    @testset "m = 1 & n = 4" begin
+        coeffs = CurveFit.__calc_integral_rules(Float64, 1:4; m = 1)
+
+        len = length(x)
+        nY, mY = 1 + (len - 1) ÷ 1, 2 * n + 1
+
+        Y = Matrix{Float64}(undef, nY, mY)
+        S = Matrix{Float64}(undef, nY, 4 - 1)
+
+        CurveFit.__cumulative_integrals!(Y, S, coeffs, x, y, 4, 1)
+
+        @test cumints_analytic[:, 1] ≈ Y[:, 1] rtol = 3.0e-3
+        @test cumints_analytic[:, 2] ≈ Y[:, 2] rtol = 3.0e-3
+        @test cumints_analytic[:, 3] ≈ Y[:, 3] rtol = 4.0e-3
+        @test cumints_analytic[:, 4] ≈ Y[:, 4] rtol = 4.0e-3
+    end
+
+    @testset "m = 2 & n = 4" begin
+        coeffs = CurveFit.__calc_integral_rules(Float64, 1:4; m = 2)
+
+        len = length(x)
+        nY, mY = 1 + (len - 1) ÷ 2, 2 * n + 1
+
+        Y = Matrix{Float64}(undef, nY, mY)
+        S = Matrix{Float64}(undef, nY, 4 - 1)
+
+        CurveFit.__cumulative_integrals!(Y, S, coeffs, x, y, 4, 2)
+
+        @test cumints_analytic[1:2:end, 1] ≈ Y[:, 1] rtol = 3.0e-5
+        @test cumints_analytic[1:2:end, 2] ≈ Y[:, 2] rtol = 4.0e-5
+        @test cumints_analytic[1:2:end, 3] ≈ Y[:, 3] rtol = 5.0e-5
+        @test cumints_analytic[1:2:end, 4] ≈ Y[:, 4] rtol = 6.0e-5
+    end
+
+    x = collect(range(0, stop = 2, length = 13))
+    y = @. exp(x)
+    cumints_analytic = @. [(exp(x) - 1) (exp(x) - 1 - x) (exp(x) - 1 - x - x^2 / 2) (
+        exp(x) -
+            1 - x -
+            x^2 /
+            2 -
+            x^3 /
+            6
+    )]
+
+    @testset "m = 2 & n = 4" begin
+        coeffs = CurveFit.__calc_integral_rules(Float64, 1:4; m = 2)
+
+        len = length(x)
+        nY, mY = 1 + (len - 1) ÷ 2, 2 * n + 1
+
+        Y = Matrix{Float64}(undef, nY, mY)
+        S = Matrix{Float64}(undef, nY, 4 - 1)
+
+        CurveFit.__cumulative_integrals!(Y, S, coeffs, x, y, 4, 2)
+
+        @test cumints_analytic[1:2:end, 1] ≈ Y[:, 1] rtol = 5.0e-6
+        @test cumints_analytic[1:2:end, 2] ≈ Y[:, 2] rtol = 3.0e-5
+        @test cumints_analytic[1:2:end, 3] ≈ Y[:, 3] rtol = 3.0e-5
+        @test cumints_analytic[1:2:end, 4] ≈ Y[:, 4] rtol = 4.0e-5
+    end
+end
diff --git a/test/expsumfit_integration_rules.jl b/test/expsumfit_integration_rules.jl
new file mode 100644
index 00000000..a1091d1c
--- /dev/null
+++ b/test/expsumfit_integration_rules.jl
@@ -0,0 +1,21 @@
+using CurveFit
+using Test
+
+@testset "ExpSumFit: Integration rules" begin
+    # Trapezoidal rule
+    @test CurveFit.__calc_integral_rules(Rational{Int64}, 1:1; m = 1) == [1 // 2 1 // 2]
+
+    # Simpson's first (1/3) rule
+    @test CurveFit.__calc_integral_rules(Rational{Int64}, 1:1; m = 2) ==
+        [1 // 3 4 // 3 1 // 3]
+    @test CurveFit.__calc_integral_rules(Rational{Int64}, 2:2; m = 2) ==
+        [2 // 3 4 // 3 0 // 1]
+    @test CurveFit.__calc_integral_rules(Rational{Int64}, 3:3; m = 2) ==
+        [3 // 5 4 // 5 -1 // 15]
+
+    # Simpson's second (3/8) rule
+    @test CurveFit.__calc_integral_rules(Rational{Int64}, 1:1; m = 3) ==
+        [3 // 8 9 // 8 9 // 8 3 // 8]
+    @test CurveFit.__calc_integral_rules(Rational{Int64}, 2:2; m = 3) ==
+        [39 // 40 27 // 10 27 // 40 3 // 20]
+end
diff --git a/test/generic.jl b/test/generic.jl
deleted file mode 100644
index a96e8304..00000000
--- a/test/generic.jl
+++ /dev/null
@@ -1,39 +0,0 @@
-@testitem "CurveFitSolution" tags = [:generic] begin
-    using SciMLBase
-
-    x = 1:10
-    fn(a, x) = @. 1.0 + 2.0 * x + a[1]
-    y = fn([1.0], x)
-    linear_sol = solve(CurveFitProblem(x, y), LinearCurveFitAlgorithm())
-    nonlinear_sol = solve(NonlinearCurveFitProblem(fn, x, y, [1.0]))
-
-    @test SciMLBase.successful_retcode(linear_sol)
-
-    # Smoke test
-    @test contains(repr(MIME"text/plain"(), linear_sol), "residuals mean:")
-    @test contains(repr(MIME"text/plain"(), nonlinear_sol), "residuals mean:")
-end
-
-@testitem "GenericNonlinearCurveFitCache show" tags = [:generic] begin
-    using CurveFit
-    using SciMLBase
-    using NonlinearSolveFirstOrder: LevenbergMarquardt, GaussNewton, TrustRegion
-
-    x = collect(1.0:10.0)
-    fn(a, x) = @. a[1] + a[2] * x
-    y = fn([1.0, 2.0], x)
-    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5], x, y)
-
-    # Smoke tests with various algorithms
-    cache = init(prob)
-    @test_nowarn repr(MIME"text/plain"(), cache)
-
-    cache = init(prob, LevenbergMarquardt())
-    @test_nowarn repr(MIME"text/plain"(), cache)
-
-    # Smoke test with a solved problem
-    cache = init(prob)
-    sol = solve!(cache)
-    @test SciMLBase.successful_retcode(sol)
-    @test_nowarn repr(MIME"text/plain"(), cache)
-end
diff --git a/test/generic_cache_show.jl b/test/generic_cache_show.jl
new file mode 100644
index 00000000..efdaa451
--- /dev/null
+++ b/test/generic_cache_show.jl
@@ -0,0 +1,24 @@
+using CurveFit
+using Test
+using SciMLBase
+using NonlinearSolveFirstOrder: LevenbergMarquardt, GaussNewton, TrustRegion
+
+@testset "GenericNonlinearCurveFitCache show" begin
+    x = collect(1.0:10.0)
+    fn(a, x) = @. a[1] + a[2] * x
+    y = fn([1.0, 2.0], x)
+    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5], x, y)
+
+    # Smoke tests with various algorithms
+    cache = init(prob)
+    @test_nowarn repr(MIME"text/plain"(), cache)
+
+    cache = init(prob, LevenbergMarquardt())
+    @test_nowarn repr(MIME"text/plain"(), cache)
+
+    # Smoke test with a solved problem
+    cache = init(prob)
+    sol = solve!(cache)
+    @test SciMLBase.successful_retcode(sol)
+    @test_nowarn repr(MIME"text/plain"(), cache)
+end
diff --git a/test/generic_curvefitsolution.jl b/test/generic_curvefitsolution.jl
new file mode 100644
index 00000000..6bb2e681
--- /dev/null
+++ b/test/generic_curvefitsolution.jl
@@ -0,0 +1,17 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "CurveFitSolution" begin
+    x = 1:10
+    fn(a, x) = @. 1.0 + 2.0 * x + a[1]
+    y = fn([1.0], x)
+    linear_sol = solve(CurveFitProblem(x, y), LinearCurveFitAlgorithm())
+    nonlinear_sol = solve(NonlinearCurveFitProblem(fn, x, y, [1.0]))
+
+    @test SciMLBase.successful_retcode(linear_sol)
+
+    # Smoke test
+    @test contains(repr(MIME"text/plain"(), linear_sol), "residuals mean:")
+    @test contains(repr(MIME"text/plain"(), nonlinear_sol), "residuals mean:")
+end
diff --git a/test/king_law.jl b/test/king_law.jl
new file mode 100644
index 00000000..2f7f24e8
--- /dev/null
+++ b/test/king_law.jl
@@ -0,0 +1,23 @@
+using CurveFit
+using Test
+
+@testset "King's Law" begin
+    U = range(1, stop = 20, length = 20)
+    A = 5.0
+    B = 1.5
+    n = 0.5
+    E = sqrt.(A .+ B * U .^ n)
+
+    fn(E) = ((E .^ 2 .- A) / B) .^ (1 ./ n)
+
+    prob = CurveFitProblem(E, U)
+    sol = solve(prob, KingCurveFitAlgorithm())
+
+    @testset for val in range(minimum(E), stop = maximum(E), length = 10)
+        @test sol(val) ≈ fn(val)
+    end
+
+    # Sigma not supported
+    prob_sigma = CurveFitProblem(E, U; sigma = ones(length(U)))
+    @test_throws AssertionError solve(prob_sigma, KingCurveFitAlgorithm())
+end
diff --git a/test/king.jl b/test/king_modified.jl
similarity index 60%
rename from test/king.jl
rename to test/king_modified.jl
index eaf4b72d..918c237e 100644
--- a/test/king.jl
+++ b/test/king_modified.jl
@@ -1,27 +1,9 @@
-@testitem "King's Law" tags = [:king] begin
-    U = range(1, stop = 20, length = 20)
-    A = 5.0
-    B = 1.5
-    n = 0.5
-    E = sqrt.(A .+ B * U .^ n)
-
-    fn(E) = ((E .^ 2 .- A) / B) .^ (1 ./ n)
-
-    prob = CurveFitProblem(E, U)
-    sol = solve(prob, KingCurveFitAlgorithm())
-
-    @testset for val in range(minimum(E), stop = maximum(E), length = 10)
-        @test sol(val) ≈ fn(val)
-    end
-
-    # Sigma not supported
-    prob_sigma = CurveFitProblem(E, U; sigma = ones(length(U)))
-    @test_throws AssertionError solve(prob_sigma, KingCurveFitAlgorithm())
-end
-
-@testitem "Modified King's Law" tags = [:king] begin
-    using NonlinearSolveFirstOrder
+using CurveFit
+using Test
+using SciMLBase
+using NonlinearSolveFirstOrder
 
+@testset "Modified King's Law" begin
     U = range(1, stop = 20, length = 20)
     A = 5.0
     B = 1.5
diff --git a/test/lin_choice.jl b/test/lin_choice.jl
index 6cee5569..910db820 100644
--- a/test/lin_choice.jl
+++ b/test/lin_choice.jl
@@ -1,9 +1,10 @@
-@testitem "Linear solver choice testing" tags = [:lin_choice] begin
-    @testset "LUFactorization" begin
-        using CurveFit
-        using LinearSolve
-        using NonlinearSolveFirstOrder
+using CurveFit
+using Test
+using LinearSolve
+using NonlinearSolveFirstOrder
 
+@testset "Linear solver choice testing" begin
+    @testset "LUFactorization" begin
         X = collect(1.0:10.0)
         θ_true = [3.0, 2.0]
 
@@ -23,10 +24,6 @@
     end
 
     @testset "QRFactorization" begin
-        using CurveFit
-        using LinearSolve
-        using NonlinearSolveFirstOrder
-
         X = collect(1.0:10.0)
         θ_true = [3.0, 2.0]
 
@@ -46,10 +43,6 @@
     end
 
     @testset "CholeskyFactorization" begin
-        using CurveFit
-        using LinearSolve
-        using NonlinearSolveFirstOrder
-
         X = collect(1.0:10.0)
         θ_true = [3.0, 2.0]
 
diff --git a/test/linfit.jl b/test/linfit.jl
deleted file mode 100644
index aa716f71..00000000
--- a/test/linfit.jl
+++ /dev/null
@@ -1,161 +0,0 @@
-@testitem "Linear Fit" tags = [:linfit] begin
-    x = range(1, stop = 10, length = 10)
-
-    fn(x) = 1.0 + 2.0 * x
-    y = fn.(x)
-
-    prob = CurveFitProblem(x, y)
-    sol = solve(prob, LinearCurveFitAlgorithm())
-
-    @test sol.u[1] ≈ 2.0
-    @test sol.u[2] ≈ 1.0
-
-    @testset for val in (0.0, 1.5, 4.5, 10.0)
-        @test sol(val) ≈ fn(val)
-    end
-end
-
-@testitem "Log Fit" tags = [:linfit] begin
-    x = range(1, stop = 10, length = 10)
-
-    fn(x) = 1.0 + 2.0 * log(x)
-    y = fn.(x)
-
-    prob = CurveFitProblem(x, y)
-    sol = solve(prob, LogCurveFitAlgorithm())
-
-    @test sol.u[1] ≈ 2.0
-    @test sol.u[2] ≈ 1.0
-
-    @testset for val in (0.0, 1.5, 4.5, 10.0)
-        @test sol(val) ≈ fn(val)
-    end
-end
-
-@testitem "Power Fit" tags = [:linfit] begin
-    x = range(1, stop = 10, length = 10)
-
-    fn(x) = 2.0 * x .^ 0.8
-    y = fn(x)
-
-    prob = CurveFitProblem(x, y)
-    sol = solve(prob, PowerCurveFitAlgorithm())
-
-    @test sol.u[1] ≈ 0.8
-    @test sol.u[2] ≈ log(2.0)
-
-    @testset for val in (0.0, 1.5, 4.5, 10.0)
-        @test sol(val) ≈ fn(val)
-    end
-
-    # Sigma not supported
-    prob_sigma = CurveFitProblem(x, y; sigma = ones(length(y)))
-    @test_throws AssertionError solve(prob_sigma, PowerCurveFitAlgorithm())
-end
-
-@testitem "Exp Fit" tags = [:linfit] begin
-    x = range(1, stop = 10, length = 10)
-
-    fn(x) = 2.0 * exp.(0.8 * x)
-    y = fn(x)
-
-    prob = CurveFitProblem(x, y)
-    sol = solve(prob, ExpCurveFitAlgorithm())
-
-    @test sol.u[1] ≈ 0.8
-    @test sol.u[2] ≈ log(2.0)
-
-    @testset for val in (0.0, 1.5, 4.5, 10.0)
-        @test sol(val) ≈ fn(val)
-    end
-
-    # Sigma not supported
-    prob_sigma = CurveFitProblem(x, y; sigma = ones(length(y)))
-    @test_throws AssertionError solve(prob_sigma, ExpCurveFitAlgorithm())
-end
-
-@testitem "Polynomial Fit" tags = [:linfit] begin
-    using LinearSolve
-
-    x = range(1, stop = 10, length = 10)
-
-    fn(x) = 1.0 + 2.0 * x + 3.0 * x^2 + 0.5 * x^3
-    y = fn.(x)
-
-    prob = CurveFitProblem(x, y)
-    sol = solve(prob, PolynomialFitAlgorithm(degree = 4))
-
-    @test sol.u[1] ≈ 1.0
-    @test sol.u[2] ≈ 2.0
-    @test sol.u[3] ≈ 3.0
-    @test sol.u[4] ≈ 0.5
-    @test sol.u[5] ≈ 0.0 atol = 1.0e-8
-
-    @testset for val in (0.0, 1.5, 4.5, 10.0)
-        @test sol(val) ≈ fn(val)
-    end
-
-    @testset "ill-conditioned" begin
-        true_coeffs = [80.0, -5.0e-18, -7.0e-20, -1.0e-36]
-        x1 = 1.0e10 .* (0:0.1:5)
-        y1 = evalpoly.(x1, (true_coeffs,))
-
-        prob = CurveFitProblem(x1, y1)
-        sol = solve(prob, PolynomialFitAlgorithm(3, QRFactorization()))
-
-        @test sol.u[1] ≈ true_coeffs[1] rtol = 1.0e-5
-        @test sol.u[2] ≈ true_coeffs[2] rtol = 1.0e-5
-        @test sol.u[3] ≈ true_coeffs[3] rtol = 1.0e-5
-        @test sol.u[4] ≈ true_coeffs[4] rtol = 1.0e-5
-
-        @testset for val in (0.0, 1.5, 4.5, 10.0)
-            @test sol(val) ≈ evalpoly(val, true_coeffs)
-        end
-
-        sol = solve(
-            prob,
-            PolynomialFitAlgorithm(3);
-            assumptions = OperatorAssumptions(
-                false; condition = OperatorCondition.VeryIllConditioned
-            )
-        )
-
-        @test sol.u[1] ≈ true_coeffs[1] rtol = 1.0e-5
-        @test sol.u[2] ≈ true_coeffs[2] rtol = 1.0e-5
-        @test sol.u[3] ≈ true_coeffs[3] rtol = 1.0e-5
-        @test sol.u[4] ≈ true_coeffs[4] rtol = 1.0e-5
-
-        @testset for val in (0.0, 1.5, 4.5, 10.0)
-            @test sol(val) ≈ evalpoly(val, true_coeffs)
-        end
-    end
-end
-
-@testitem "Linear Weighted Least Squares" tags = [:linfit] begin
-    using SciMLBase
-
-    x = 1:10
-    a0 = [2.0, 1.0]
-    fn(x) = a0[1] * x + a0[2]
-    y = fn.(x)
-
-    # Add a large outlier
-    y[5] += 20.0
-
-    # Fit without sigma - outlier should affect result
-    prob = CurveFitProblem(x, y)
-    sol_unweighted = solve(prob, LinearCurveFitAlgorithm())
-
-    # Fit with a high sigma on the outlier
-    sigma = ones(length(y))
-    sigma[5] = 100.0
-    prob = CurveFitProblem(x, y; sigma)
-    sol_weighted = solve(prob, LinearCurveFitAlgorithm())
-
-    # Weighted fit should (hopefully) be closer to the true parameters
-    err_no_weight = maximum(abs.(sol_unweighted.u .- a0))
-    err_weighted = maximum(abs.(sol_weighted.u .- a0))
-
-    @test err_weighted < err_no_weight
-    @test SciMLBase.successful_retcode(sol_weighted)
-end
diff --git a/test/linfit_exp.jl b/test/linfit_exp.jl
new file mode 100644
index 00000000..5662df89
--- /dev/null
+++ b/test/linfit_exp.jl
@@ -0,0 +1,23 @@
+using CurveFit
+using Test
+
+@testset "Exp Fit" begin
+    x = range(1, stop = 10, length = 10)
+
+    fn(x) = 2.0 * exp.(0.8 * x)
+    y = fn(x)
+
+    prob = CurveFitProblem(x, y)
+    sol = solve(prob, ExpCurveFitAlgorithm())
+
+    @test sol.u[1] ≈ 0.8
+    @test sol.u[2] ≈ log(2.0)
+
+    @testset for val in (0.0, 1.5, 4.5, 10.0)
+        @test sol(val) ≈ fn(val)
+    end
+
+    # Sigma not supported
+    prob_sigma = CurveFitProblem(x, y; sigma = ones(length(y)))
+    @test_throws AssertionError solve(prob_sigma, ExpCurveFitAlgorithm())
+end
diff --git a/test/linfit_linear.jl b/test/linfit_linear.jl
new file mode 100644
index 00000000..4f8cb1cc
--- /dev/null
+++ b/test/linfit_linear.jl
@@ -0,0 +1,19 @@
+using CurveFit
+using Test
+
+@testset "Linear Fit" begin
+    x = range(1, stop = 10, length = 10)
+
+    fn(x) = 1.0 + 2.0 * x
+    y = fn.(x)
+
+    prob = CurveFitProblem(x, y)
+    sol = solve(prob, LinearCurveFitAlgorithm())
+
+    @test sol.u[1] ≈ 2.0
+    @test sol.u[2] ≈ 1.0
+
+    @testset for val in (0.0, 1.5, 4.5, 10.0)
+        @test sol(val) ≈ fn(val)
+    end
+end
diff --git a/test/linfit_log.jl b/test/linfit_log.jl
new file mode 100644
index 00000000..704a24f1
--- /dev/null
+++ b/test/linfit_log.jl
@@ -0,0 +1,19 @@
+using CurveFit
+using Test
+
+@testset "Log Fit" begin
+    x = range(1, stop = 10, length = 10)
+
+    fn(x) = 1.0 + 2.0 * log(x)
+    y = fn.(x)
+
+    prob = CurveFitProblem(x, y)
+    sol = solve(prob, LogCurveFitAlgorithm())
+
+    @test sol.u[1] ≈ 2.0
+    @test sol.u[2] ≈ 1.0
+
+    @testset for val in (0.0, 1.5, 4.5, 10.0)
+        @test sol(val) ≈ fn(val)
+    end
+end
diff --git a/test/linfit_polynomial.jl b/test/linfit_polynomial.jl
new file mode 100644
index 00000000..31cf99b0
--- /dev/null
+++ b/test/linfit_polynomial.jl
@@ -0,0 +1,58 @@
+using CurveFit
+using Test
+using LinearSolve
+
+@testset "Polynomial Fit" begin
+    x = range(1, stop = 10, length = 10)
+
+    fn(x) = 1.0 + 2.0 * x + 3.0 * x^2 + 0.5 * x^3
+    y = fn.(x)
+
+    prob = CurveFitProblem(x, y)
+    sol = solve(prob, PolynomialFitAlgorithm(degree = 4))
+
+    @test sol.u[1] ≈ 1.0
+    @test sol.u[2] ≈ 2.0
+    @test sol.u[3] ≈ 3.0
+    @test sol.u[4] ≈ 0.5
+    @test sol.u[5] ≈ 0.0 atol = 1.0e-8
+
+    @testset for val in (0.0, 1.5, 4.5, 10.0)
+        @test sol(val) ≈ fn(val)
+    end
+
+    @testset "ill-conditioned" begin
+        true_coeffs = [80.0, -5.0e-18, -7.0e-20, -1.0e-36]
+        x1 = 1.0e10 .* (0:0.1:5)
+        y1 = evalpoly.(x1, (true_coeffs,))
+
+        prob = CurveFitProblem(x1, y1)
+        sol = solve(prob, PolynomialFitAlgorithm(3, QRFactorization()))
+
+        @test sol.u[1] ≈ true_coeffs[1] rtol = 1.0e-5
+        @test sol.u[2] ≈ true_coeffs[2] rtol = 1.0e-5
+        @test sol.u[3] ≈ true_coeffs[3] rtol = 1.0e-5
+        @test sol.u[4] ≈ true_coeffs[4] rtol = 1.0e-5
+
+        @testset for val in (0.0, 1.5, 4.5, 10.0)
+            @test sol(val) ≈ evalpoly(val, true_coeffs)
+        end
+
+        sol = solve(
+            prob,
+            PolynomialFitAlgorithm(3);
+            assumptions = OperatorAssumptions(
+                false; condition = OperatorCondition.VeryIllConditioned
+            )
+        )
+
+        @test sol.u[1] ≈ true_coeffs[1] rtol = 1.0e-5
+        @test sol.u[2] ≈ true_coeffs[2] rtol = 1.0e-5
+        @test sol.u[3] ≈ true_coeffs[3] rtol = 1.0e-5
+        @test sol.u[4] ≈ true_coeffs[4] rtol = 1.0e-5
+
+        @testset for val in (0.0, 1.5, 4.5, 10.0)
+            @test sol(val) ≈ evalpoly(val, true_coeffs)
+        end
+    end
+end
diff --git a/test/linfit_power.jl b/test/linfit_power.jl
new file mode 100644
index 00000000..5134c81b
--- /dev/null
+++ b/test/linfit_power.jl
@@ -0,0 +1,23 @@
+using CurveFit
+using Test
+
+@testset "Power Fit" begin
+    x = range(1, stop = 10, length = 10)
+
+    fn(x) = 2.0 * x .^ 0.8
+    y = fn(x)
+
+    prob = CurveFitProblem(x, y)
+    sol = solve(prob, PowerCurveFitAlgorithm())
+
+    @test sol.u[1] ≈ 0.8
+    @test sol.u[2] ≈ log(2.0)
+
+    @testset for val in (0.0, 1.5, 4.5, 10.0)
+        @test sol(val) ≈ fn(val)
+    end
+
+    # Sigma not supported
+    prob_sigma = CurveFitProblem(x, y; sigma = ones(length(y)))
+    @test_throws AssertionError solve(prob_sigma, PowerCurveFitAlgorithm())
+end
diff --git a/test/linfit_weighted.jl b/test/linfit_weighted.jl
new file mode 100644
index 00000000..94d8c5c0
--- /dev/null
+++ b/test/linfit_weighted.jl
@@ -0,0 +1,30 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "Linear Weighted Least Squares" begin
+    x = 1:10
+    a0 = [2.0, 1.0]
+    fn(x) = a0[1] * x + a0[2]
+    y = fn.(x)
+
+    # Add a large outlier
+    y[5] += 20.0
+
+    # Fit without sigma - outlier should affect result
+    prob = CurveFitProblem(x, y)
+    sol_unweighted = solve(prob, LinearCurveFitAlgorithm())
+
+    # Fit with a high sigma on the outlier
+    sigma = ones(length(y))
+    sigma[5] = 100.0
+    prob = CurveFitProblem(x, y; sigma)
+    sol_weighted = solve(prob, LinearCurveFitAlgorithm())
+
+    # Weighted fit should (hopefully) be closer to the true parameters
+    err_no_weight = maximum(abs.(sol_unweighted.u .- a0))
+    err_weighted = maximum(abs.(sol_weighted.u .- a0))
+
+    @test err_weighted < err_no_weight
+    @test SciMLBase.successful_retcode(sol_weighted)
+end
diff --git a/test/nonlinfit.jl b/test/nonlinfit.jl
deleted file mode 100644
index aa66ddc1..00000000
--- a/test/nonlinfit.jl
+++ /dev/null
@@ -1,391 +0,0 @@
-@testitem "Nonlinear Least Squares: Linear Problem 1" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    x = 1.0:10.0
-    a0 = [3.0, 2.0, 1.0]
-
-    fn(a, x) = @. a[1] + a[2] * x + a[3] * x^2
-    y = fn(a0, x)
-
-    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], x, y)
-    sol = solve(prob)
-
-    @test sol.u ≈ a0
-    @test SciMLBase.successful_retcode(sol.retcode)
-
-    @testset for val in (0.0, 1.5, 4.5, 10.0)
-        @test sol(val) ≈ fn(a0, val)
-    end
-end
-
-@testitem "Nonlinear Least Squares: Nonlinear Problem 1" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    x = 1.0:10.0
-    a0 = [3.0, 2.0, 0.7]
-
-    fn(a, x) = @. a[1] + a[2] * x^a[3]
-    y = fn(a0, x)
-
-    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], x, y)
-    sol = solve(prob)
-
-    @test sol.u ≈ a0
-    @test SciMLBase.successful_retcode(sol.retcode)
-
-    @testset for val in (0.0, 1.5, 4.5, 10.0)
-        @test sol(val) ≈ fn(a0, val)
-    end
-end
-
-@testitem "Nonlinear Least Squares: Linear Problem 2" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    x = 1.0:10.0
-    a0 = [3.0, 2.0, 1.0]
-
-    fn(a, x) = @. a[1] + a[2] * x[:, 1] + a[3] * x[:, 1]^2 - x[:, 2]
-    P = length(x)
-    X = zeros(P, 2)
-    for i in 1:P
-        X[i, 1] = x[i]
-        X[i, 2] = fn(a0, [x[i] 0])[1]
-    end
-
-    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], X)
-    sol = solve(prob)
-
-    @test sol.u ≈ a0 atol = 1.0e-7
-    @test SciMLBase.successful_retcode(sol.retcode)
-
-    @testset for val in (0.0, 1.5, 4.5, 10.0)
-        @test sol([val 0.0])[1] ≈ fn(a0, [val 0.0])[1] atol = 1.0e-7
-    end
-end
-
-@testitem "Nonlinear Least Squares: Nonlinear Problem 2" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    x = 1.0:10.0
-    a0 = [3.0, 2.0, 0.7]
-
-    fn(a, x) = @. a[1] + a[2] * x[:, 1]^a[3] - x[:, 2]
-    P = length(x)
-    X = zeros(P, 2)
-    for i in 1:P
-        X[i, 1] = x[i]
-        X[i, 2] = fn(a0, [x[i] 0])[1]
-    end
-
-    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], X)
-    sol = solve(prob)
-
-    @test sol.u ≈ a0 atol = 1.0e-7
-    @test SciMLBase.successful_retcode(sol.retcode)
-
-    @testset for val in (0.0, 1.5, 4.5, 10.0)
-        @test sol([val 0.0])[1] ≈ fn(a0, [val 0.0])[1] atol = 1.0e-7
-    end
-end
-
-@testitem "Nonlinear Least Squares: reinit!()" tags = [:nonlinfit] begin
-    using CurveFit
-    using SciMLBase
-    using NonlinearSolveBase: NonlinearSolveBase
-
-    # Create an initial problem with a cache
-    x = 1.0:10.0
-    a0 = [3.0, 2.0, 0.7]
-
-    fn(a, x) = @. a[1] + a[2] * x^a[3]
-    y = fn(a0, x)
-
-    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], x, y)
-    cache = CurveFit.init(prob)
-    @test solve!(cache).u ≈ a0 atol = 1.0e-7
-
-    # reinit!() the cache with different parameters and recheck the solve
-    a0 = [4.0, 5.0, 0.2]
-    x = 11.0:20.0
-    y = fn(a0, x)
-
-    CurveFit.reinit!(cache; u0 = [1.0, 1.0, 1.0], x, y)
-    @test solve!(cache).u ≈ a0 atol = 1.0e-7
-
-    # Repeat with an in-place model: NonlinearSolve wraps in-place Float64
-    # problems in `AutoSpecializeCallable`, which `reinit!` must unwrap.
-    fn(resid, a, x) = @. resid = a[1] + a[2] * x^a[3]
-    cache = CurveFit.init(NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], x, y))
-    @test cache.cache.prob.f.f isa NonlinearSolveBase.AutoSpecializeCallable
-    CurveFit.reinit!(cache; u0 = [1.0, 1.0, 1.0], x, y)
-    @test solve!(cache).u ≈ a0 atol = 1.0e-7
-end
-
-@testitem "Nonlinear Weighted Least Squares" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    fn(a, x) = @. a[1] + a[2] * x
-    x = collect(1.0:10.0)
-    a0 = [1.0, 2.0]
-    y = fn(a0, x)
-
-    # Add a large outlier
-    y[5] += 20.0
-
-    # Fit without sigma - outlier should affect result
-    prob_no_weight = NonlinearCurveFitProblem(fn, [0.5, 0.5], x, y)
-    sol_no_weight = solve(prob_no_weight)
-
-    # Fit with a high sigma on the outlier
-    sigma = ones(length(y))
-    sigma[5] = 100.0
-    prob_weighted = NonlinearCurveFitProblem(fn, [0.5, 0.5], x, y, sigma)
-    sol_weighted = solve(prob_weighted)
-
-    # Weighted fit should (hopefully) be closer to the true parameters
-    err_no_weight = maximum(abs.(sol_no_weight.u .- a0))
-    err_weighted = maximum(abs.(sol_weighted.u .- a0))
-
-    @test err_weighted < err_no_weight
-    @test SciMLBase.successful_retcode(sol_weighted)
-end
-
-@testitem "Gauss-Newton curve fitting: Linear problem" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    U = 0.5:0.5:10
-    a0 = [2.0, 1.0, 0.35]
-    E = @. sqrt(a0[1] + a0[2] * U^a0[3])
-
-    X = hcat(E, U)
-    fn(a, x) = @. a[1] + a[2] * x[:, 2]^a[3] - x[:, 1]^2
-
-    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], X)
-    sol = solve(prob)
-
-    @test sol.u ≈ a0 atol = 1.0e-7
-    @test SciMLBase.successful_retcode(sol.retcode)
-
-    @testset for val in range(minimum(E), stop = maximum(E), length = 10)
-        @test sol([val 0.0])[1] ≈ fn(a0, [val 0.0])[1] atol = 1.0e-7
-    end
-end
-
-# Regression test for https://github.com/SciML/CurveFit.jl/issues/69
-@testitem "Issue #69: y ~ a/x with noisy data" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    # Original issue: fitting y ~ a/x failed with NaN when data had noise
-    x = [1.0, 2.0, 3.0, 4.0, 5.0]
-
-    # Test case 1: Exact data (should work perfectly)
-    fn(a, x) = @. a[1] / x
-    y_exact = [1.0, 0.5, 1 / 3, 0.25, 0.2]
-    prob1 = NonlinearCurveFitProblem(fn, [0.1], x, y_exact)
-    sol1 = solve(prob1)
-
-    @test sol1.u[1] ≈ 1.0 atol = 1.0e-6
-    @test SciMLBase.successful_retcode(sol1.retcode)
-    @test !isnan(sol1.u[1])
-
-    # Test case 2: Data with tiny noise (previously failed with NaN)
-    y_noisy = [1.0, 0.5, 1 / 3 + 0.00000001, 0.25, 0.2]
-    prob2 = NonlinearCurveFitProblem(fn, [0.1], x, y_noisy)
-    sol2 = solve(prob2)
-
-    @test sol2.u[1] ≈ 1.0 atol = 1.0e-5
-    @test SciMLBase.successful_retcode(sol2.retcode)
-    @test !isnan(sol2.u[1])
-
-    # Test case 3: y ~ ax exact data
-    fn2(a, x) = @. a[1] * x
-    y2_exact = [1.0, 2.0, 3.0, 4.0, 5.0]
-    prob3 = NonlinearCurveFitProblem(fn2, [0.1], x, y2_exact)
-    sol3 = solve(prob3)
-
-    @test sol3.u[1] ≈ 1.0 atol = 1.0e-6
-    @test SciMLBase.successful_retcode(sol3.retcode)
-    @test !isnan(sol3.u[1])
-
-    # Test case 4: y ~ ax with tiny noise (previously failed with NaN)
-    y2_noisy = [1.0, 2.0, 3.000001, 4.0, 5.0]
-    prob4 = NonlinearCurveFitProblem(fn2, [0.1], x, y2_noisy)
-    sol4 = solve(prob4)
-
-    @test sol4.u[1] ≈ 1.0 atol = 1.0e-5
-    @test SciMLBase.successful_retcode(sol4.retcode)
-    @test !isnan(sol4.u[1])
-end
-
-@testitem "Issue #69: Robustness with larger noise" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    x = [1.0, 2.0, 3.0, 4.0, 5.0]
-
-    # y ~ a/x with larger noise
-    fn(a, x) = @. a[1] / x
-    y_noisy = [1.0 + 0.01, 0.5 - 0.01, 1 / 3 + 0.01, 0.25, 0.2 - 0.005]
-    prob = NonlinearCurveFitProblem(fn, [0.1], x, y_noisy)
-    sol = solve(prob)
-
-    @test sol.u[1] ≈ 1.0 atol = 0.1  # Looser tolerance for noisy data
-    @test SciMLBase.successful_retcode(sol.retcode)
-    @test !isnan(sol.u[1])
-end
-
-# Tests for ScalarModel - Issue #46
-@testitem "ScalarModel: Basic polynomial fitting" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    x = 1.0:10.0
-    a0 = [3.0, 2.0, 1.0]
-
-    # Scalar function (no @.)
-    fn_scalar(a, x) = a[1] + a[2] * x + a[3] * x^2
-
-    # Generate y data by broadcasting the scalar function
-    y = fn_scalar.(Ref(a0), x)
-
-    # Use ScalarModel wrapper
-    prob = NonlinearCurveFitProblem(ScalarModel(fn_scalar), [0.5, 0.5, 0.5], x, y)
-    @test SciMLBase.isinplace(prob.nlfunc)
-    sol = solve(prob)
-    @test sol(x) isa Vector
-
-    @test sol.u ≈ a0
-    @test SciMLBase.successful_retcode(sol.retcode)
-
-    # Test single-point evaluation
-    @testset for val in (0.0, 1.5, 4.5, 10.0)
-        @test sol(val) ≈ fn_scalar(a0, val)
-    end
-end
-
-@testitem "ScalarModel: Nonlinear power function" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    x = 1.0:10.0
-    a0 = [3.0, 2.0, 0.7]
-
-    # Scalar function (no @.)
-    fn_scalar(a, x) = a[1] + a[2] * x^a[3]
-
-    # Generate y data
-    y = fn_scalar.(Ref(a0), x)
-
-    # Use ScalarModel wrapper
-    prob = NonlinearCurveFitProblem(ScalarModel(fn_scalar), [0.5, 0.5, 0.5], x, y)
-    sol = solve(prob)
-
-    @test sol.u ≈ a0
-    @test SciMLBase.successful_retcode(sol.retcode)
-
-    # Test single-point evaluation
-    @testset for val in (0.0, 1.5, 4.5, 10.0)
-        @test sol(val) ≈ fn_scalar(a0, val)
-    end
-end
-
-@testitem "ScalarModel: Exponential decay (LsqFit-style migration)" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    # This test demonstrates migration from LsqFit.jl style
-    # LsqFit: model(x, p) = p[1] * exp(-x * p[2])
-    # CurveFit: model(p, x) = p[1] * exp(-x * p[2])
-
-    x = collect(range(0, stop = 10, length = 20))
-    true_params = [2.5, 0.3]
-
-    # Scalar model function (parameter order: params first, then x)
-    model(p, x) = p[1] * exp(-x * p[2])
-
-    # Generate y data
-    y = model.(Ref(true_params), x)
-
-    # Use ScalarModel wrapper
-    prob = NonlinearCurveFitProblem(ScalarModel(model), [1.0, 0.1], x, y)
-    sol = solve(prob)
-
-    @test sol.u ≈ true_params atol = 1.0e-6
-    @test SciMLBase.successful_retcode(sol.retcode)
-
-    # Verify predictions
-    @test sol(5.0) ≈ model(true_params, 5.0)
-end
-
-@testitem "ScalarModel: Reciprocal function" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    x = [1.0, 2.0, 3.0, 4.0, 5.0]
-    a0 = [2.0]
-
-    # Scalar function: y = a/x
-    fn_scalar(a, x) = a[1] / x
-    y = fn_scalar.(Ref(a0), x)
-
-    prob = NonlinearCurveFitProblem(ScalarModel(fn_scalar), [0.5], x, y)
-    sol = solve(prob)
-
-    @test sol.u[1] ≈ a0[1] atol = 1.0e-6
-    @test SciMLBase.successful_retcode(sol.retcode)
-end
-
-@testitem "ScalarModel: Equivalence with vectorized @. form" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    # Verify that ScalarModel produces the same results as @. form
-
-    x = 1.0:10.0
-    a0 = [3.0, 2.0, 0.7]
-    u0 = [0.5, 0.5, 0.5]
-
-    # Vectorized form (standard CurveFit style)
-    fn_vec(a, x) = @. a[1] + a[2] * x^a[3]
-    y = fn_vec(a0, x)
-
-    prob_vec = NonlinearCurveFitProblem(fn_vec, u0, x, y)
-    sol_vec = solve(prob_vec)
-
-    # Scalar form with ScalarModel
-    fn_scalar(a, x) = a[1] + a[2] * x^a[3]
-
-    prob_scalar = NonlinearCurveFitProblem(ScalarModel(fn_scalar), u0, x, y)
-    sol_scalar = solve(prob_scalar)
-
-    # Both should give the same result
-    @test sol_vec.u ≈ sol_scalar.u
-    @test SciMLBase.successful_retcode(sol_vec.retcode)
-    @test SciMLBase.successful_retcode(sol_scalar.retcode)
-
-    # Both should evaluate the same at any point
-    @testset for val in (0.0, 1.5, 4.5, 10.0)
-        @test sol_vec(val) ≈ sol_scalar(val)
-    end
-end
-
-@testitem "Nonlinear Least Squares: Bounds constrain solution" tags = [:nonlinfit] begin
-    using SciMLBase
-
-    # True params are [3.0, 1.0] but we constrain p[1] within [0.0, 2.0]
-    fn(a, x) = @. a[1] * exp(a[2] * x)
-    x = collect(range(0, 2, length = 20))
-    y = fn.(Ref([3, 1]), x)
-
-    lb = [0.0, -Inf]
-    ub = [2.0, Inf]
-
-    prob = NonlinearCurveFitProblem(fn, [1.0, 0.5], x, y; lb, ub)
-
-    # Verify bounds are stored on the CurveFitProblem
-    @test prob.lb == lb
-    @test prob.ub == ub
-
-    # Test solve() path
-    sol = solve(prob)
-    @test sol.u[1] <= 2.0
-
-    # Test init+solve!() path
-    cache = CurveFit.init(prob)
-    sol2 = solve!(cache)
-    @test sol2.u[1] <= 2.0
-end
diff --git a/test/nonlinfit_bounds.jl b/test/nonlinfit_bounds.jl
new file mode 100644
index 00000000..cd3ca781
--- /dev/null
+++ b/test/nonlinfit_bounds.jl
@@ -0,0 +1,28 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "Nonlinear Least Squares: Bounds constrain solution" begin
+    # True params are [3.0, 1.0] but we constrain p[1] within [0.0, 2.0]
+    fn(a, x) = @. a[1] * exp(a[2] * x)
+    x = collect(range(0, 2, length = 20))
+    y = fn.(Ref([3, 1]), x)
+
+    lb = [0.0, -Inf]
+    ub = [2.0, Inf]
+
+    prob = NonlinearCurveFitProblem(fn, [1.0, 0.5], x, y; lb, ub)
+
+    # Verify bounds are stored on the CurveFitProblem
+    @test prob.lb == lb
+    @test prob.ub == ub
+
+    # Test solve() path
+    sol = solve(prob)
+    @test sol.u[1] <= 2.0
+
+    # Test init+solve!() path
+    cache = CurveFit.init(prob)
+    sol2 = solve!(cache)
+    @test sol2.u[1] <= 2.0
+end
diff --git a/test/nonlinfit_gauss_newton.jl b/test/nonlinfit_gauss_newton.jl
new file mode 100644
index 00000000..8456387a
--- /dev/null
+++ b/test/nonlinfit_gauss_newton.jl
@@ -0,0 +1,22 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "Gauss-Newton curve fitting: Linear problem" begin
+    U = 0.5:0.5:10
+    a0 = [2.0, 1.0, 0.35]
+    E = @. sqrt(a0[1] + a0[2] * U^a0[3])
+
+    X = hcat(E, U)
+    fn(a, x) = @. a[1] + a[2] * x[:, 2]^a[3] - x[:, 1]^2
+
+    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], X)
+    sol = solve(prob)
+
+    @test sol.u ≈ a0 atol = 1.0e-7
+    @test SciMLBase.successful_retcode(sol.retcode)
+
+    @testset for val in range(minimum(E), stop = maximum(E), length = 10)
+        @test sol([val 0.0])[1] ≈ fn(a0, [val 0.0])[1] atol = 1.0e-7
+    end
+end
diff --git a/test/nonlinfit_issue69_noisy.jl b/test/nonlinfit_issue69_noisy.jl
new file mode 100644
index 00000000..456b12f9
--- /dev/null
+++ b/test/nonlinfit_issue69_noisy.jl
@@ -0,0 +1,47 @@
+using CurveFit
+using Test
+using SciMLBase
+
+# Regression test for https://github.com/SciML/CurveFit.jl/issues/69
+@testset "Issue #69: y ~ a/x with noisy data" begin
+    # Original issue: fitting y ~ a/x failed with NaN when data had noise
+    x = [1.0, 2.0, 3.0, 4.0, 5.0]
+
+    # Test case 1: Exact data (should work perfectly)
+    fn(a, x) = @. a[1] / x
+    y_exact = [1.0, 0.5, 1 / 3, 0.25, 0.2]
+    prob1 = NonlinearCurveFitProblem(fn, [0.1], x, y_exact)
+    sol1 = solve(prob1)
+
+    @test sol1.u[1] ≈ 1.0 atol = 1.0e-6
+    @test SciMLBase.successful_retcode(sol1.retcode)
+    @test !isnan(sol1.u[1])
+
+    # Test case 2: Data with tiny noise (previously failed with NaN)
+    y_noisy = [1.0, 0.5, 1 / 3 + 0.00000001, 0.25, 0.2]
+    prob2 = NonlinearCurveFitProblem(fn, [0.1], x, y_noisy)
+    sol2 = solve(prob2)
+
+    @test sol2.u[1] ≈ 1.0 atol = 1.0e-5
+    @test SciMLBase.successful_retcode(sol2.retcode)
+    @test !isnan(sol2.u[1])
+
+    # Test case 3: y ~ ax exact data
+    fn2(a, x) = @. a[1] * x
+    y2_exact = [1.0, 2.0, 3.0, 4.0, 5.0]
+    prob3 = NonlinearCurveFitProblem(fn2, [0.1], x, y2_exact)
+    sol3 = solve(prob3)
+
+    @test sol3.u[1] ≈ 1.0 atol = 1.0e-6
+    @test SciMLBase.successful_retcode(sol3.retcode)
+    @test !isnan(sol3.u[1])
+
+    # Test case 4: y ~ ax with tiny noise (previously failed with NaN)
+    y2_noisy = [1.0, 2.0, 3.000001, 4.0, 5.0]
+    prob4 = NonlinearCurveFitProblem(fn2, [0.1], x, y2_noisy)
+    sol4 = solve(prob4)
+
+    @test sol4.u[1] ≈ 1.0 atol = 1.0e-5
+    @test SciMLBase.successful_retcode(sol4.retcode)
+    @test !isnan(sol4.u[1])
+end
diff --git a/test/nonlinfit_issue69_robust.jl b/test/nonlinfit_issue69_robust.jl
new file mode 100644
index 00000000..f242b343
--- /dev/null
+++ b/test/nonlinfit_issue69_robust.jl
@@ -0,0 +1,17 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "Issue #69: Robustness with larger noise" begin
+    x = [1.0, 2.0, 3.0, 4.0, 5.0]
+
+    # y ~ a/x with larger noise
+    fn(a, x) = @. a[1] / x
+    y_noisy = [1.0 + 0.01, 0.5 - 0.01, 1 / 3 + 0.01, 0.25, 0.2 - 0.005]
+    prob = NonlinearCurveFitProblem(fn, [0.1], x, y_noisy)
+    sol = solve(prob)
+
+    @test sol.u[1] ≈ 1.0 atol = 0.1  # Looser tolerance for noisy data
+    @test SciMLBase.successful_retcode(sol.retcode)
+    @test !isnan(sol.u[1])
+end
diff --git a/test/nonlinfit_linear_problem_1.jl b/test/nonlinfit_linear_problem_1.jl
new file mode 100644
index 00000000..eaa7cfc0
--- /dev/null
+++ b/test/nonlinfit_linear_problem_1.jl
@@ -0,0 +1,21 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "Nonlinear Least Squares: Linear Problem 1" begin
+    x = 1.0:10.0
+    a0 = [3.0, 2.0, 1.0]
+
+    fn(a, x) = @. a[1] + a[2] * x + a[3] * x^2
+    y = fn(a0, x)
+
+    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], x, y)
+    sol = solve(prob)
+
+    @test sol.u ≈ a0
+    @test SciMLBase.successful_retcode(sol.retcode)
+
+    @testset for val in (0.0, 1.5, 4.5, 10.0)
+        @test sol(val) ≈ fn(a0, val)
+    end
+end
diff --git a/test/nonlinfit_linear_problem_2.jl b/test/nonlinfit_linear_problem_2.jl
new file mode 100644
index 00000000..02d91e22
--- /dev/null
+++ b/test/nonlinfit_linear_problem_2.jl
@@ -0,0 +1,26 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "Nonlinear Least Squares: Linear Problem 2" begin
+    x = 1.0:10.0
+    a0 = [3.0, 2.0, 1.0]
+
+    fn(a, x) = @. a[1] + a[2] * x[:, 1] + a[3] * x[:, 1]^2 - x[:, 2]
+    P = length(x)
+    X = zeros(P, 2)
+    for i in 1:P
+        X[i, 1] = x[i]
+        X[i, 2] = fn(a0, [x[i] 0])[1]
+    end
+
+    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], X)
+    sol = solve(prob)
+
+    @test sol.u ≈ a0 atol = 1.0e-7
+    @test SciMLBase.successful_retcode(sol.retcode)
+
+    @testset for val in (0.0, 1.5, 4.5, 10.0)
+        @test sol([val 0.0])[1] ≈ fn(a0, [val 0.0])[1] atol = 1.0e-7
+    end
+end
diff --git a/test/nonlinfit_nonlinear_problem_1.jl b/test/nonlinfit_nonlinear_problem_1.jl
new file mode 100644
index 00000000..7bf4fad1
--- /dev/null
+++ b/test/nonlinfit_nonlinear_problem_1.jl
@@ -0,0 +1,21 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "Nonlinear Least Squares: Nonlinear Problem 1" begin
+    x = 1.0:10.0
+    a0 = [3.0, 2.0, 0.7]
+
+    fn(a, x) = @. a[1] + a[2] * x^a[3]
+    y = fn(a0, x)
+
+    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], x, y)
+    sol = solve(prob)
+
+    @test sol.u ≈ a0
+    @test SciMLBase.successful_retcode(sol.retcode)
+
+    @testset for val in (0.0, 1.5, 4.5, 10.0)
+        @test sol(val) ≈ fn(a0, val)
+    end
+end
diff --git a/test/nonlinfit_nonlinear_problem_2.jl b/test/nonlinfit_nonlinear_problem_2.jl
new file mode 100644
index 00000000..5fef2117
--- /dev/null
+++ b/test/nonlinfit_nonlinear_problem_2.jl
@@ -0,0 +1,26 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "Nonlinear Least Squares: Nonlinear Problem 2" begin
+    x = 1.0:10.0
+    a0 = [3.0, 2.0, 0.7]
+
+    fn(a, x) = @. a[1] + a[2] * x[:, 1]^a[3] - x[:, 2]
+    P = length(x)
+    X = zeros(P, 2)
+    for i in 1:P
+        X[i, 1] = x[i]
+        X[i, 2] = fn(a0, [x[i] 0])[1]
+    end
+
+    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], X)
+    sol = solve(prob)
+
+    @test sol.u ≈ a0 atol = 1.0e-7
+    @test SciMLBase.successful_retcode(sol.retcode)
+
+    @testset for val in (0.0, 1.5, 4.5, 10.0)
+        @test sol([val 0.0])[1] ≈ fn(a0, [val 0.0])[1] atol = 1.0e-7
+    end
+end
diff --git a/test/nonlinfit_reinit.jl b/test/nonlinfit_reinit.jl
new file mode 100644
index 00000000..1f4acf35
--- /dev/null
+++ b/test/nonlinfit_reinit.jl
@@ -0,0 +1,33 @@
+using CurveFit
+using Test
+using SciMLBase
+using NonlinearSolveBase: NonlinearSolveBase
+
+@testset "Nonlinear Least Squares: reinit!()" begin
+    # Create an initial problem with a cache
+    x = 1.0:10.0
+    a0 = [3.0, 2.0, 0.7]
+
+    fn(a, x) = @. a[1] + a[2] * x^a[3]
+    y = fn(a0, x)
+
+    prob = NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], x, y)
+    cache = CurveFit.init(prob)
+    @test solve!(cache).u ≈ a0 atol = 1.0e-7
+
+    # reinit!() the cache with different parameters and recheck the solve
+    a0 = [4.0, 5.0, 0.2]
+    x = 11.0:20.0
+    y = fn(a0, x)
+
+    CurveFit.reinit!(cache; u0 = [1.0, 1.0, 1.0], x, y)
+    @test solve!(cache).u ≈ a0 atol = 1.0e-7
+
+    # Repeat with an in-place model: NonlinearSolve wraps in-place Float64
+    # problems in `AutoSpecializeCallable`, which `reinit!` must unwrap.
+    fn(resid, a, x) = @. resid = a[1] + a[2] * x^a[3]
+    cache = CurveFit.init(NonlinearCurveFitProblem(fn, [0.5, 0.5, 0.5], x, y))
+    @test cache.cache.prob.f.f isa NonlinearSolveBase.AutoSpecializeCallable
+    CurveFit.reinit!(cache; u0 = [1.0, 1.0, 1.0], x, y)
+    @test solve!(cache).u ≈ a0 atol = 1.0e-7
+end
diff --git a/test/nonlinfit_scalarmodel_equivalence.jl b/test/nonlinfit_scalarmodel_equivalence.jl
new file mode 100644
index 00000000..75e3f377
--- /dev/null
+++ b/test/nonlinfit_scalarmodel_equivalence.jl
@@ -0,0 +1,34 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "ScalarModel: Equivalence with vectorized @. form" begin
+    # Verify that ScalarModel produces the same results as @. form
+
+    x = 1.0:10.0
+    a0 = [3.0, 2.0, 0.7]
+    u0 = [0.5, 0.5, 0.5]
+
+    # Vectorized form (standard CurveFit style)
+    fn_vec(a, x) = @. a[1] + a[2] * x^a[3]
+    y = fn_vec(a0, x)
+
+    prob_vec = NonlinearCurveFitProblem(fn_vec, u0, x, y)
+    sol_vec = solve(prob_vec)
+
+    # Scalar form with ScalarModel
+    fn_scalar(a, x) = a[1] + a[2] * x^a[3]
+
+    prob_scalar = NonlinearCurveFitProblem(ScalarModel(fn_scalar), u0, x, y)
+    sol_scalar = solve(prob_scalar)
+
+    # Both should give the same result
+    @test sol_vec.u ≈ sol_scalar.u
+    @test SciMLBase.successful_retcode(sol_vec.retcode)
+    @test SciMLBase.successful_retcode(sol_scalar.retcode)
+
+    # Both should evaluate the same at any point
+    @testset for val in (0.0, 1.5, 4.5, 10.0)
+        @test sol_vec(val) ≈ sol_scalar(val)
+    end
+end
diff --git a/test/nonlinfit_scalarmodel_exp_decay.jl b/test/nonlinfit_scalarmodel_exp_decay.jl
new file mode 100644
index 00000000..20014d5f
--- /dev/null
+++ b/test/nonlinfit_scalarmodel_exp_decay.jl
@@ -0,0 +1,28 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "ScalarModel: Exponential decay (LsqFit-style migration)" begin
+    # This test demonstrates migration from LsqFit.jl style
+    # LsqFit: model(x, p) = p[1] * exp(-x * p[2])
+    # CurveFit: model(p, x) = p[1] * exp(-x * p[2])
+
+    x = collect(range(0, stop = 10, length = 20))
+    true_params = [2.5, 0.3]
+
+    # Scalar model function (parameter order: params first, then x)
+    model(p, x) = p[1] * exp(-x * p[2])
+
+    # Generate y data
+    y = model.(Ref(true_params), x)
+
+    # Use ScalarModel wrapper
+    prob = NonlinearCurveFitProblem(ScalarModel(model), [1.0, 0.1], x, y)
+    sol = solve(prob)
+
+    @test sol.u ≈ true_params atol = 1.0e-6
+    @test SciMLBase.successful_retcode(sol.retcode)
+
+    # Verify predictions
+    @test sol(5.0) ≈ model(true_params, 5.0)
+end
diff --git a/test/nonlinfit_scalarmodel_polynomial.jl b/test/nonlinfit_scalarmodel_polynomial.jl
new file mode 100644
index 00000000..71b0a452
--- /dev/null
+++ b/test/nonlinfit_scalarmodel_polynomial.jl
@@ -0,0 +1,29 @@
+using CurveFit
+using Test
+using SciMLBase
+
+# Tests for ScalarModel - Issue #46
+@testset "ScalarModel: Basic polynomial fitting" begin
+    x = 1.0:10.0
+    a0 = [3.0, 2.0, 1.0]
+
+    # Scalar function (no @.)
+    fn_scalar(a, x) = a[1] + a[2] * x + a[3] * x^2
+
+    # Generate y data by broadcasting the scalar function
+    y = fn_scalar.(Ref(a0), x)
+
+    # Use ScalarModel wrapper
+    prob = NonlinearCurveFitProblem(ScalarModel(fn_scalar), [0.5, 0.5, 0.5], x, y)
+    @test SciMLBase.isinplace(prob.nlfunc)
+    sol = solve(prob)
+    @test sol(x) isa Vector
+
+    @test sol.u ≈ a0
+    @test SciMLBase.successful_retcode(sol.retcode)
+
+    # Test single-point evaluation
+    @testset for val in (0.0, 1.5, 4.5, 10.0)
+        @test sol(val) ≈ fn_scalar(a0, val)
+    end
+end
diff --git a/test/nonlinfit_scalarmodel_power.jl b/test/nonlinfit_scalarmodel_power.jl
new file mode 100644
index 00000000..8ae09538
--- /dev/null
+++ b/test/nonlinfit_scalarmodel_power.jl
@@ -0,0 +1,26 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "ScalarModel: Nonlinear power function" begin
+    x = 1.0:10.0
+    a0 = [3.0, 2.0, 0.7]
+
+    # Scalar function (no @.)
+    fn_scalar(a, x) = a[1] + a[2] * x^a[3]
+
+    # Generate y data
+    y = fn_scalar.(Ref(a0), x)
+
+    # Use ScalarModel wrapper
+    prob = NonlinearCurveFitProblem(ScalarModel(fn_scalar), [0.5, 0.5, 0.5], x, y)
+    sol = solve(prob)
+
+    @test sol.u ≈ a0
+    @test SciMLBase.successful_retcode(sol.retcode)
+
+    # Test single-point evaluation
+    @testset for val in (0.0, 1.5, 4.5, 10.0)
+        @test sol(val) ≈ fn_scalar(a0, val)
+    end
+end
diff --git a/test/nonlinfit_scalarmodel_reciprocal.jl b/test/nonlinfit_scalarmodel_reciprocal.jl
new file mode 100644
index 00000000..4a10afa7
--- /dev/null
+++ b/test/nonlinfit_scalarmodel_reciprocal.jl
@@ -0,0 +1,18 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "ScalarModel: Reciprocal function" begin
+    x = [1.0, 2.0, 3.0, 4.0, 5.0]
+    a0 = [2.0]
+
+    # Scalar function: y = a/x
+    fn_scalar(a, x) = a[1] / x
+    y = fn_scalar.(Ref(a0), x)
+
+    prob = NonlinearCurveFitProblem(ScalarModel(fn_scalar), [0.5], x, y)
+    sol = solve(prob)
+
+    @test sol.u[1] ≈ a0[1] atol = 1.0e-6
+    @test SciMLBase.successful_retcode(sol.retcode)
+end
diff --git a/test/nonlinfit_weighted.jl b/test/nonlinfit_weighted.jl
new file mode 100644
index 00000000..c49fdc8b
--- /dev/null
+++ b/test/nonlinfit_weighted.jl
@@ -0,0 +1,30 @@
+using CurveFit
+using Test
+using SciMLBase
+
+@testset "Nonlinear Weighted Least Squares" begin
+    fn(a, x) = @. a[1] + a[2] * x
+    x = collect(1.0:10.0)
+    a0 = [1.0, 2.0]
+    y = fn(a0, x)
+
+    # Add a large outlier
+    y[5] += 20.0
+
+    # Fit without sigma - outlier should affect result
+    prob_no_weight = NonlinearCurveFitProblem(fn, [0.5, 0.5], x, y)
+    sol_no_weight = solve(prob_no_weight)
+
+    # Fit with a high sigma on the outlier
+    sigma = ones(length(y))
+    sigma[5] = 100.0
+    prob_weighted = NonlinearCurveFitProblem(fn, [0.5, 0.5], x, y, sigma)
+    sol_weighted = solve(prob_weighted)
+
+    # Weighted fit should (hopefully) be closer to the true parameters
+    err_no_weight = maximum(abs.(sol_no_weight.u .- a0))
+    err_weighted = maximum(abs.(sol_weighted.u .- a0))
+
+    @test err_weighted < err_no_weight
+    @test SciMLBase.successful_retcode(sol_weighted)
+end
diff --git a/test/nopre/Project.toml b/test/nopre/Project.toml
new file mode 100644
index 00000000..a6376b5b
--- /dev/null
+++ b/test/nopre/Project.toml
@@ -0,0 +1,16 @@
+[deps]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+CurveFit = "5a033b19-8c74-5913-a970-47c3779ef25c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[sources]
+CurveFit = {path = "../.."}
+
+[compat]
+Aqua = "0.8.12"
+SafeTestsets = "0.1"
+SciMLTesting = "1"
+Test = "1.10"
+julia = "1.10"
diff --git a/test/nopre/aqua.jl b/test/nopre/aqua.jl
new file mode 100644
index 00000000..4f309bff
--- /dev/null
+++ b/test/nopre/aqua.jl
@@ -0,0 +1,7 @@
+using CurveFit
+using Test
+using Aqua
+
+@testset "Aqua tests" begin
+    Aqua.test_all(CurveFit)
+end
diff --git a/test/rationalfit_linear.jl b/test/rationalfit_linear.jl
new file mode 100644
index 00000000..b8ebe943
--- /dev/null
+++ b/test/rationalfit_linear.jl
@@ -0,0 +1,18 @@
+using CurveFit
+using Test
+using LinearSolve, LinearAlgebra
+
+@testset "Linear Rational fit" begin
+    x = range(1, stop = 10, length = 10)
+    r = CurveFit.RationalPolynomial([1.0, 0.0, -2.0], [1.0, 2.0, 3.0])
+    y = r.(x)
+
+    prob = CurveFitProblem(x, y)
+    sol = solve(prob, RationalPolynomialFitAlgorithm(2, 3, QRFactorization(ColumnNorm())))
+
+    @test sol.u ≈ [1.0, 0.0, -2.0, 2.0, 3.0, 0.0] atol = 1.0e-8
+
+    @testset for val in (0.0, 1.5, 4.5, 10.0)
+        @test sol(val) ≈ r(val) atol = 1.0e-8
+    end
+end
diff --git a/test/rationalfit.jl b/test/rationalfit_nonlinear.jl
similarity index 63%
rename from test/rationalfit.jl
rename to test/rationalfit_nonlinear.jl
index e3cabb5b..3f915c27 100644
--- a/test/rationalfit.jl
+++ b/test/rationalfit_nonlinear.jl
@@ -1,23 +1,9 @@
-@testitem "Linear Rational fit" tags = [:rationalfit] begin
-    using LinearSolve, LinearAlgebra
-
-    x = range(1, stop = 10, length = 10)
-    r = CurveFit.RationalPolynomial([1.0, 0.0, -2.0], [1.0, 2.0, 3.0])
-    y = r.(x)
-
-    prob = CurveFitProblem(x, y)
-    sol = solve(prob, RationalPolynomialFitAlgorithm(2, 3, QRFactorization(ColumnNorm())))
-
-    @test sol.u ≈ [1.0, 0.0, -2.0, 2.0, 3.0, 0.0] atol = 1.0e-8
-
-    @testset for val in (0.0, 1.5, 4.5, 10.0)
-        @test sol(val) ≈ r(val) atol = 1.0e-8
-    end
-end
-
-@testitem "Nonlinear Rational fit" tags = [:rationalfit] begin
-    using NonlinearSolveFirstOrder
+using CurveFit
+using Test
+using SciMLBase
+using NonlinearSolveFirstOrder
 
+@testset "Nonlinear Rational fit" begin
     x = range(1, stop = 10, length = 10)
     r = CurveFit.RationalPolynomial([1.0, 0.0, -2.0], [1.0, 2.0, 3.0])
     y = r.(x)
diff --git a/test/runtests.jl b/test/runtests.jl
index 1fa8ed0f..a18a7cc4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,13 +1,2 @@
-using TestItemRunner
-
-# First try to get the group from the argument list, otherwise fall back to the
-# environment variable.
-# Example command for running the linfit tests: Pkg.test(; test_args=["nonlinfit"])
-const GROUP = Symbol(get(ARGS, 1, get(ENV, "GROUP", "All")))
-
-if GROUP == :All
-    # Run all tests except those tagged with :nopre
-    @run_package_tests filter = ti -> !(:nopre in ti.tags)
-else
-    @run_package_tests verbose = true filter = ti -> (GROUP in ti.tags)
-end
+using SciMLTesting
+run_tests()
diff --git a/test/stats.jl b/test/stats.jl
index 9a00b853..fc2d6e8d 100644
--- a/test/stats.jl
+++ b/test/stats.jl
@@ -1,9 +1,11 @@
-@testitem "StatsAPI Integration" tags = [:stats] begin
-    using StatsAPI
-    using NonlinearSolveFirstOrder
-    using LinearAlgebra
-    using LinearSolve
-
+using CurveFit
+using Test
+using StatsAPI
+using NonlinearSolveFirstOrder
+using LinearAlgebra
+using LinearSolve
+
+@testset "StatsAPI Integration" begin
     @testset "Linear Fit" begin
         # y = 2x + 1
         x = [1.0, 2.0, 3.0, 4.0, 5.0]
diff --git a/test/test_groups.toml b/test/test_groups.toml
index efba756b..a3ce918e 100644
--- a/test/test_groups.toml
+++ b/test/test_groups.toml
@@ -1,13 +1,23 @@
 # Root test-group matrix for CurveFit.jl, consumed by the reusable
 # grouped-tests.yml workflow via scripts/compute_affected_sublibraries.jl
-# --root-matrix. runtests.jl reads the GROUP env var and dispatches via
-# TestItemRunner tags: "All" runs everything except :nopre-tagged items,
-# "nopre" runs only :nopre-tagged items.
+# --root-matrix, and by SciMLTesting's run_tests() folder-discovery mode in
+# runtests.jl. Each table key is a GROUP value passed to runtests.jl.
+#
+#   "Core"  -> the top-level test/*.jl files (the functional fits plus the
+#              ExplicitImports check). Runs on 1/lts/pre.
+#   "nopre" -> the test/nopre/ folder (Aqua), which is brittle on prereleases,
+#              so it runs only on 1/lts. `in_all = false` keeps it out of a
+#              local GROUP=All aggregate run.
+#
+# This replaces the former TestItemRunner tag split (GROUP=All ran every
+# @testitem except :nopre; GROUP=nopre ran only the :nopre-tagged Aqua item):
+# the same tests run under the same two (group x version) CI cells.
 
-[All]
+[Core]
 versions = ["1", "lts", "pre"]
 os = ["ubuntu-latest", "macOS-latest", "windows-latest"]
 
 [nopre]
+in_all = false
 versions = ["1", "lts"]
 os = ["ubuntu-latest", "macOS-latest", "windows-latest"]