add breakpoint rank estimation method and more tests

Author: njericha

.gitignore

@@ -1,3 +1,4 @@
 Manifest.toml
 docs/build
 .vscode
+Breakpoint_calc.m

src/Core/curvaturetools.jl

@@ -141,7 +141,11 @@ end
 
 Approximates the signed curvature of a function given evenly spaced samples.
 
-Uses [`d_dx`](@ref) and [`d2_dx2`](@ref) to approximate the first two derivatives.
+# Possible `method`s
+- `:finite_differences`: Approximates first and second derivative with 3rd order finite differences. See [`d_dx`](@ref) and [`d2_dx2`](@ref).
+- `:splines`: Curvature of a third order spline. See [`d_dx_and_d2_dx2_spline`](@ref).
+- `:circles`: Inverse radius of a circle through rolling three points. See [`circle_curvature`](@ref).
+- `:breakpoints`: WARNING does not compute a value that approximates the curvature of a continuous function. Computes the inverse least-squares error of `f.(eachindex(y); z)` and `y` for all `z in eachindex(y)` where `f(x; z) = a + b(min(x, z) - z) + c(max(x, z) - z)`. Useful if `y` looks like two lines. See [`breakpoint_curvature`](@ref).
 """
 function curvature(y::AbstractVector{<:Real}; method=:finite_differences, kwargs...)
     if method == :finite_differences
@@ -153,6 +157,8 @@ function curvature(y::AbstractVector{<:Real}; method=:finite_differences, kwargs
         return @. dy2_dx2 / (1 + dy_dx^2)^1.5
     elseif method == :circles
         return circle_curvature(y; h=1)
+    elseif method == :breakpoints
+        return breakpoint_curvature(y)
     else
         throw(ArgumentError("method $method not implemented"))
     end
@@ -166,6 +172,13 @@ Approximates the signed curvature of a function, scaled to the unit box ``[0,1]^
 Assumes the function is 1 at 0 and (after x dimension is scaled) 0 at 1.
 
 See [`curvature`](@ref).
+
+
+# Possible `method`s
+- `:finite_differences`: Approximates first and second derivative with 3rd order finite differences. See [`d_dx`](@ref) and [`d2_dx2`](@ref).
+- `:splines`: Curvature of a third order spline. See [`d_dx_and_d2_dx2_spline`](@ref).
+- `:circles`: Inverse radius of a circle through rolling three points. See [`circle_curvature`](@ref).
+- `:breakpoints`: WARNING does not compute a value that approximates the curvature of a continuous function. Computes the inverse least-squares error of `f.(eachindex(y); z)` and `y` for all `z in eachindex(y)` where `f(x; z) = a + b(min(x, z) - z) + c(max(x, z) - z)`. Useful if `y` looks like two lines. See [`breakpoint_curvature`](@ref).
 """
 function standard_curvature(y::AbstractVector{<:Real}; method=:finite_differences, kwargs...)
     Δx = 1/length(y)
@@ -183,6 +196,8 @@ function standard_curvature(y::AbstractVector{<:Real}; method=:finite_difference
         return @. dy2_dx2 / (1 + dy_dx^2)^1.5
     elseif method == :circles
         return circle_curvature(y / max(1,maximum(y)); h=Δx)
+    elseif method == :breakpoints
+        return breakpoint_curvature(y) # best breakpoint unaffected by scaling and stretching
     else
         throw(ArgumentError("method $method not implemented"))
     end
@@ -264,3 +279,51 @@ function signed_circle_curvature((a,f),(b,g),(c,h))
     sign = g > (f+h)/2 ? -1 : 1
     return sign / r
 end
+
+"""
+    breakpoint_model_coefficients(xs, ys, breakpoint)
+
+Least squares fit data ``(x_i, y_i)``
+
+``\\min_{a,b,c} 0.5\\sum_{i} (f(x_i; a,b,c) - y_i)^2``
+
+with the model
+
+``f(x; a,b,c) = a + b(\\min(x, z) - x) + c(\\max(x, z) - x)``
+
+for some fixed ``z``.
+"""
+function breakpoint_model_coefficients(xs, ys, z)
+    n = length(xs)
+    @assert n == length(ys)
+    M = hcat(ones(n), (min.(xs, z) .- z), (max.(xs, z) .- z))
+    a, b, c = M \ ys
+    return a, b, c
+end
+
+breakpoint_model(a, b, c, z) = x -> a + b*(min(x, z) - z) + c*(max(x, z) - z)
+
+function breakpoint_error(xs, ys, z)
+    a, b, c = breakpoint_model_coefficients(xs, ys, z)
+    f = breakpoint_model(a, b, c, z)
+    return norm2(@. f(xs) - ys)
+    # equivalent to sum(((x, y),) -> (f(x) - y)^2, zip(xs, ys))
+end
+
+best_breakpoint(xs, ys; breakpoints=xs) = argmin(z -> breakpoint_error(xs, ys, z), breakpoints)
+
+"""
+    breakpoint_curvature(y)
+
+This is a hacked way to fit the data `y` with a breakpoint model,
+which can be called by `k = standard_curvature(...; model=:breakpoints)`
+
+This lets us call `argmax(k)` to get the breakpoint that minimizes the model error.
+
+See [`breakpoint_model_coefficients`](@ref).
+"""
+function breakpoint_curvature(y)
+    x = eachindex(y)
+    errors = [breakpoint_error(x, y, z) for z in x]
+    return 1 ./ errors
+end

src/Core/rankdetection.jl

@@ -7,7 +7,7 @@ Selects the rank that maximizes the standard curvature of the Relative Error (as
 
 # Keywords
 - `online_rank_estimation`: `false`. Set to `true` to stop testing larger ranks after the first peak in curvature
-- `curvature_method`: `:splines`. Can also pick `:finite_differences` (faster but less accurate) or `circles` (fastest and smallest memory but more sensitive to results from `factorize`)
+- `curvature_method`: `:splines`. Can also pick `:finite_differences` (faster but less accurate) or `circles` (fastest and smallest memory but more sensitive to results from `factorize`). Set to `:breakpoints` to pick the rank `R` that minimizes least-squares error in the model `f(r) = a + b(min(r, R) - R) + c(max(r, R) - R)` and the errors.
 - `model`: `Tucker1`. Only rank detection with `Tucker1` and `CPDecomposition` is currently implemented
 - `max_rank`: `max_possible_rank(Y, model)`. Test ranks from `1` up to `max_rank`. Defaults to largest possible rank under the model
 - `rank`: `nothing`. If a rank is passed, rank detection is ignored and `factorize(Y; kwargs...)` is called

test/runtests.jl

@@ -894,6 +894,42 @@ end
 
     check_slice = 1:10
     @test isapprox(decomposition[check_slice], Y[check_slice]; rtol=0.01) # should be within 1% error
+
+    N = 33 # 1 plus a power of 2
+    R = 5
+    D = 2
+
+    # More constrained problem
+    matrices = [abs_randn(N, R) for _ in 1:D]
+    l1scale_cols!.(matrices)
+    Ydecomp = CPDecomposition(Tuple(matrices))#abs_randn
+    @assert all(check.(simplex_cols!, factors(Ydecomp)))
+    Y = array(Ydecomp)
+
+    scaleB_rescaleA! = ConstraintUpdate(0, l1scale_1slices! ∘ nonnegative!;
+        whats_rescaled=(x -> eachcol(factor(x, 1)))
+    )
+    nonnegativeB! = ConstraintUpdate(0, nonnegative!)
+    nonnegativeA! = ConstraintUpdate(1, nonnegative!)
+    #[l1scale_1slices! ∘ nonnegative!, nonnegative!]
+
+    options = (
+        rank=3,
+        momentum=true,
+        model=Tucker1,
+        tolerance=(1e-5),
+        converged=(GradientNNCone),
+        do_subblock_updates=false,
+        constrain_init=true,
+        constraints=[scaleB_rescaleA!, nonnegativeA!],
+        stats=[Iteration, ObjectiveValue, GradientNNCone, RelativeError],
+        maxiter=200
+    )
+
+    decomposition, stats, kwargs = multiscale_factorize(Y; options...)
+
+    check_slice = 1:10
+    @test isapprox(decomposition[check_slice], Y[check_slice]; rtol=0.02) # should be within 2% error
     end
 end
 
@@ -927,6 +963,13 @@ end
     @test MAPE(k_circles, k_true) < 0.03 # 3%
     @test MAPE(k_finite_differences, k_true) < 0.03 # 3%
 
+    # Break point method
+    ys = [13,10,5,4.5,4,3.6,3,2.5]; xs = [0,1,2,3,4,5,6,7]; z=2
+    a,b,c = BlockTensorFactorization.Core.breakpoint_model_coefficients(xs, ys, z)
+
+    @test isapprox([a,b,c], [5.2,-4.1,-0.5], rtol=0.01)
+
+    # Rank detect
     T = Tucker1((10, 10, 10), 3)
     Y = array(T)
     decomposition, stats, kwargs, final_rel_errors = rank_detect_factorize(Y; model=Tucker1, curvature_method=:splines)
@@ -938,10 +981,17 @@ end
     decomposition, stats, kwargs, final_rel_errors = rank_detect_factorize(Y; model=Tucker1, curvature_method=:finite_differences)
     @test kwargs[:rank] == 3
 
+    decomposition, stats, kwargs, final_rel_errors = rank_detect_factorize(Y; model=Tucker1, curvature_method=:breakpoints)
+    @test kwargs[:rank] == 3
+
     T = CPDecomposition((10, 11, 12), 4)
     Y = array(T)
-    decomposition, stats, kwargs, final_rel_errors = rank_detect_factorize(Y; model=CPDecomposition, curvature_method=:splines, online_rank_estimation=true)
-    @test kwargs[:rank] == 4
+    V = zeros(Int, 5)
+    for i in 1:5
+        decomposition, stats, kwargs, final_rel_errors = rank_detect_factorize(Y; model=CPDecomposition, curvature_method=:splines, online_rank_estimation=true, tolerance=0.01)
+        V[i] = kwargs[:rank]
+    end
+    @test count(x -> x == 4, V) ≥ 3 # should predict 4 most of the time
 end
 
 end