% [e,q] = seglsq(X,c) Segmented least squares
%
% This solves the segmented least squares problem using bottom-up dynamical 
% programming. That is, it splits a 2D dataset into segments and fits a line
% to each segment, so that the segmentation/fit minimises the objective
% function E + c.L, where 1<=L<=n is the number of segments, E is the sum of
% the sums of squared errors within each segment and c>0 is a given constant
% (larger c encourages fewer lines). Each segment consists of a consecutive
% sequence of points, and each point belongs to one and only one segment.
%
% This implementation runs in O(n³) because the computation of the errors for
% each segment e(i,j) is done in a naive loop with runtime O(n³). It can be
% done in O(n²) by precomputing various statistics.
%
% In:
%   X: Nx2 dataset of row vectors.
%   c: nonnegative scalar, tradeoff between fit accuracy and number of lines.
% Out:
%   e: objective function value.
%   q: indices of the beginning of the line segments.

% Copyright (c) 2013 by Max Vladymyrov and Miguel A. Carreira-Perpinan

function [e,q] = seglsq(X,c)

n = size(X,1); r = zeros(n+1,1); s = zeros(n+1,1);
r(1) = 0;

% Dynamic programming loop (with computation of fit errors)
for j=1:n
 e = zeros(j,1);
 for i=1:j-1
   e(i) = lsqerr(X(i:j,:));
 end  
 [r(j+1),s(j+1)] = min(e+c+r(1:j));
end
e = r(n+1);

% Get points at which new segments start
q(1) = s(n+1); k = 1;
while q(k)>1
 q(k+1) = s(q(k)); k = k+1;
end
q = flipdim(q(1:k),2);


% Line parameters (a,b) and sum of squared errors for a dataset A

function [e,a,b] = lsqerr(A)

n = size(A,1);

x1 = sum(A(:,1)); y1 = sum(A(:,2)); 
x2 = sum(A(:,1).^2); xy1 = sum(A(:,1).*A(:,2));

a = (n*xy1 - x1*y1)/(n*x2 - x1.^2); b = (y1 - a*x1)/n;
e = sum((A(:,2) - a*A(:,1) - b).^2);

