replace sort(), sortperm() with sort!() and sortperm!() for lts, lta, and lms

Author: jbytecode

src/lms.jl

@@ -70,24 +70,22 @@ function lms(X::AbstractMatrix{Float64}, y::AbstractVector{Float64}; iters = not
     bestobjective = Inf
     bestparamaters = Array{Float64}(undef, p)
     bestres = Array{Float64}(undef, n)
-    origres = Array{Float64}(undef, n)
     indices = collect(1:n)
     kindices = collect(p:n)
     betas = Array{Float64}(undef, p)
     res = Array{Float64}(undef, n)
 
     for _ = 1:iters
         try
-            k = rand(kindices, 1)[1]
+            k = rand(kindices)
             sampledindices = sample(indices, k, replace = false)
             betas = X[sampledindices, :] \ y[sampledindices]
-            origres = y .- X * betas
-            res = sort(origres .^ 2.0)
+            res = sort!((y .- X * betas) .^ 2.0)
             m2 = res[h]
             if m2 < bestobjective
-                bestparamaters .= betas
+                bestparamaters = betas
                 bestobjective = m2
-                bestres .= origres
+                bestres = y .- X * betas
             end
         catch e
             @warn e

src/lta.jl

@@ -85,15 +85,14 @@ function lta(X::AbstractMatrix{Float64}, y::AbstractVector{Float64}; exact = fal
         psubsets = collect(combinations(1:n, p))
     else
         iters = p * 3000
-        psubsets = [sample(1:n, p, replace = false) for i = 1:iters]
+        psubsets = [sample(1:n, p, replace = false) for _ = 1:iters]
     end
 
     function lta_cost(subsetindices::Array{Int,1})::Tuple{Float64,Vector{Float64}}
         try
             betas = X[subsetindices, :] \ y[subsetindices]
-            res_abs = abs.(y .- X * betas)
-            ordered_res = sort(res_abs)
-            cost = sum(ordered_res[1:h])
+            sorted_res_abs = sort!(abs.(y .- X * betas))
+            cost = sum(sorted_res_abs[1:h])
             return (cost, betas)
         catch
             return (Inf64, [])

src/lts.jl

@@ -45,15 +45,17 @@ function iterateCSteps(
     subsetindices::Array{Int,1},
     h::Int; eps::Float64 = 0.01, maxiter::Int = 10000
 )
+    n = length(y)
     oldobjective::Float64 = Inf64
     objective::Float64 = Inf64
     iter::Int = 0
+    sortedresindices = Array{Int}(undef, n)
     while iter < maxiter
         tempols = ols(X[subsetindices, :], y[subsetindices])
         res = y - X * coef(tempols)
-        sortedresindices = sortperm(abs.(res))
+        sortperm!(sortedresindices, abs.(res))
         subsetindices = sortedresindices[1:h]
-        objective = sum(sort(res .^ 2.0)[1:h])
+        objective = sum(sort!(res .^ 2.0)[1:h])
         if isapprox(oldobjective, objective, atol=eps)
             break
         end

test.jl

@@ -1,39 +1,9 @@
 using LinRegOutliers
+using BenchmarkTools
 
-import LinRegOutliers.OrdinaryLeastSquares: ols, wls, coef, residuals
-import Distributions: median
-
-#sett = createRegressionSetting(@formula(calls ~ year), phones)
-#sett = createRegressionSetting(@formula(y ~ x1 + x2 + x3), hbk)
-
-function generatedata(
-    n::Int,
-    p::Int,
-    contamination::Float64,
-    direction::Symbol,
-)::Tuple{Array{Float64,2},Array{Float64,1}}
-    if !(direction in [:y, :x])
-        error("- Direction must be :x or :y")
-    end 
-    on = ones(Float64, n)
-    pxvars = p - 1 
-    xvars = randn(n, pxvars)
-    b = [5.0 for i = 1:pxvars]
-    y = 5.0 .+ xvars * b + randn(n)
-    totalcontamination = round(Int, contamination * n)
-    if direction == :y
-        y[1:totalcontamination] .= maximum(y) .+ abs.(rand(totalcontamination) * 5.0)
-    elseif direction == :x
-        xvars[1:totalcontamination, :] .=
-            maximum(xvars) .+ abs.(rand(totalcontamination) * 5.0)
-    end 
-    return (hcat(on, xvars), y)
-end
 
 
 
 sett = createRegressionSetting(@formula(calls ~ year), phones)
-x = designMatrix(sett)
-y = responseVector(sett)
-x, y = generatedata(1000, 25, 0.30, :x)
-result = gwcga(x, y)
+
+@btime lta($sett)