|
| 1 | +function output = cssd(x,y,p,gamma,xx,delta) |
| 2 | +%CSSD Cubic smoothing spline with discontinuities |
| 3 | +% |
| 4 | +% cssd(x, y, p, gamma, xx, delta) computes a cubic smoothing spline with discontinuities for the |
| 5 | +% given data (x,y). The data values may be scalars or vectors. Data points with the |
| 6 | +% same site are replaced by their (weighted) average as in the builtin csaps |
| 7 | +% function. |
| 8 | +% |
| 9 | +% Input |
| 10 | +% x: vector of data sites |
| 11 | +% |
| 12 | +% y: vector of same lenght as x or matrix where y(:,i) is a data vector at site x(i) |
| 13 | +% |
| 14 | +% p: parameter between 0 and 1 that weights the rougness penalty |
| 15 | +% (high values result in smoother curves). Use CSSD_CV for automatic |
| 16 | +% selection. |
| 17 | +% |
| 18 | +% gamma: parameter between 0 and Infinity that weights the discontiuity |
| 19 | +% penalty (high values result in less discontinuities, gamma = Inf returns |
| 20 | +% a classical smoothing spline). Use CSSD_CV for automatic |
| 21 | +% selection. |
| 22 | +% |
| 23 | +% xx: (optional) evaluation points for the result |
| 24 | +% |
| 25 | +% delta: (optional) weights of the data sites. delta may be thought of as the |
| 26 | +% standard deviation of the at site x_i. Should have the same size as x. |
| 27 | +% - Note: The Matlab built in spline function csaps uses a different weight |
| 28 | +% convention (w). delta is related to Matlab's w by w = 1./delta.^2 |
| 29 | +% - Note for vector-valued data: Weights are assumed to be identical over |
| 30 | +% vector-components. (Componentwise weights might be supported in a future version.) |
| 31 | +% |
| 32 | +% Output |
| 33 | +% output = cssd(...) |
| 34 | +% output.pp: ppform of a smoothing spline with discontinuities; if xx is specified, |
| 35 | +% the evaluation of the result at the points xx is returned |
| 36 | +% output.discont: locations of detected discontinuities, the locations are a |
| 37 | +% subset of the midpoints of the data sites x |
| 38 | +% output.interval_cell: a list of discrete indices between two discontinuities |
| 39 | +% output.pp_cell: a list of the cubic splines corresponding to the indices in interval_cell |
| 40 | +% |
| 41 | +% See also CSAPS, CSSD_CV |
| 42 | + |
| 43 | +%%% BEGIN CHECK ARGUMENTS |
| 44 | +if nargin<5, xx = []; end |
| 45 | +if nargin<6, delta = []; end |
| 46 | + |
| 47 | +if isempty(delta), delta = ones(size(x)); end |
| 48 | + |
| 49 | +assert( (0 <= p) && (p <= 1), 'The p parameter must fulfill 0 <= p <= 1') |
| 50 | +assert( 0 <= gamma, 'The gamma parameter must fulfill 0 < gamma') |
| 51 | + |
| 52 | +% Matlab uses the parameter w which is related to delta of De Boor's book by w = 1./delta.^2 |
| 53 | +w = 1./delta.^2; |
| 54 | + |
| 55 | +% checks arguments and creates column vectors (chckxywp is Matlab built in) |
| 56 | +[xi,yi,~,wi] = chckxywp(x,y,2,w,p); |
| 57 | +deltai = sqrt(1./wi); |
| 58 | + |
| 59 | +% Note: from now on we use the xi, yi, wi, deltai versions |
| 60 | +%%% END CHECK ARGUMENTS |
| 61 | + |
| 62 | +[N,D] = size(yi); |
| 63 | + |
| 64 | +% if gamma == Inf (discontinuity has infinite penalty), we may directly |
| 65 | +% compute a classical smoothing spline |
| 66 | +% also, if p == 1, we may straight compute an interpolating spline, no |
| 67 | +% matter how large gamma is (smoothness costs are equal to 0) |
| 68 | +if (gamma == Inf) || (p == 1) |
| 69 | + pp = csaps(xi,yi',p,[],wi); |
| 70 | + discont = []; |
| 71 | + interval_cell = {1:N}; |
| 72 | + pp_cell = {fnxtr(pp,2)}; |
| 73 | +else |
| 74 | + % F stores Bellmann values |
| 75 | + F = zeros(N, 1); |
| 76 | + % partition: stores the optimal partition |
| 77 | + partition = zeros(N, 1); |
| 78 | + |
| 79 | + %%% BEGIN PIECEWISE LINEAR CASE |
| 80 | + if p == 0 % the piecewise linear case |
| 81 | + B = [ones(N,1), xi]./deltai(:); |
| 82 | + rhs = yi./deltai(:); |
| 83 | + % precompute eps_1r for r=1,...,N |
| 84 | + A = [B, rhs]; |
| 85 | + G = planerot(A(1:2,1)); |
| 86 | + A(1:2, :) = G*A(1:2, :); |
| 87 | + eps_1r = 0; |
| 88 | + % loop starts from index three because eps_11 and eps_12 are zero |
| 89 | + for r=3:N |
| 90 | + G = planerot(A([1,r],1)); |
| 91 | + A([1,r],:) = G * A([1,r],:); |
| 92 | + G = planerot(A([2,r],2)); |
| 93 | + A([2,r],2:end) = G * A([2,r],2:end); |
| 94 | + eps_1r = eps_1r + sum(A(r, 3:end).^2); |
| 95 | + % store the eps_1r as the initial Bellman value corresponding to a |
| 96 | + % solution without discontinuities |
| 97 | + F(r) = eps_1r; |
| 98 | + end |
| 99 | + %%% BEGIN MAIN LOOP |
| 100 | + for rb=2:N |
| 101 | + |
| 102 | + % best left bound (blb) initialized with 1 corresponding to interval 1:rb |
| 103 | + % corresponding Bellman value has been set in the precomputation |
| 104 | + blb = 1; |
| 105 | + |
| 106 | + A = [B(1:rb,:), rhs(1:rb,:)]; |
| 107 | + |
| 108 | + eps_lr = 0; |
| 109 | + % the loop is performed in reverse order so that we may use pruning |
| 110 | + for lb = rb-1:-1:2 |
| 111 | + if lb == rb-1 |
| 112 | + G = planerot(A([end,end-1],1)); |
| 113 | + A([end,end-1], :) = G*A([end,end-1], :); |
| 114 | + else |
| 115 | + G = planerot(A([end, lb],1)); |
| 116 | + A([end,lb], :) = G * A([end,lb], :); |
| 117 | + G = planerot(A([end-1, lb],2)); |
| 118 | + A([end-1,lb], 2:end) = G*A([end-1,lb], 2:end); |
| 119 | + eps_lr = eps_lr + sum(A(lb,3:end).^2); |
| 120 | + end |
| 121 | + |
| 122 | + % check if setting a discontinuity between lb-1 and lb gives a |
| 123 | + % better energy |
| 124 | + candidate_value = F(lb-1) + gamma + eps_lr; |
| 125 | + if candidate_value < F(rb) |
| 126 | + F(rb ) = candidate_value; |
| 127 | + blb = lb; |
| 128 | + end |
| 129 | + % store the best left bound corresponding to the right bound rb |
| 130 | + partition( rb ) = blb-1; |
| 131 | + end |
| 132 | + end |
| 133 | + %%% END MAIN LOOP |
| 134 | + |
| 135 | + %%% END PIECEWISE LINEAR CASE |
| 136 | + |
| 137 | + |
| 138 | + %%% BEGIN CSSD CASE |
| 139 | + else % this is the standard case: gamma > 0 and 0 < p < 1 |
| 140 | + %%% BEGIN PRECOMPUTATIONS |
| 141 | + beta = sqrt(1-p); |
| 142 | + alpha = sqrt(p)./deltai; |
| 143 | + d = diff(xi); % xi is sorted ascendingly |
| 144 | + |
| 145 | + % precompute eps_1r for r=1,...,N |
| 146 | + [eps_1r, R, z] = startEpsLR(yi(1:2,:), d(1), alpha(1:2), beta); |
| 147 | + % loop starts from index three because eps_11 and eps_12 are zero |
| 148 | + for r=3:N |
| 149 | + [eps_1r, R, z] = updateEpsLR(eps_1r, R, yi(r,:), d(r-1), z, alpha(r), beta); |
| 150 | + % store the eps_1r as the initial Bellman value corresponding to a |
| 151 | + % solution without discontinuities |
| 152 | + F(r) = eps_1r; |
| 153 | + end |
| 154 | + %%% END PRECOMPUTATIONS |
| 155 | + |
| 156 | + %%% BEGIN MAIN LOOP |
| 157 | + for rb=2:N |
| 158 | + |
| 159 | + % best left bound (blb) initialized with 1 corresponding to interval 1:rb |
| 160 | + % corresponding Bellman value has been set in the precomputation |
| 161 | + blb = 1; |
| 162 | + |
| 163 | + % the loop is performed in reverse order so that we may use pruning |
| 164 | + for lb = rb-1:-1:2 |
| 165 | + if lb == rb-1 |
| 166 | + % get start configuration and store start state in R, z |
| 167 | + [eps_lr, R, z] = startEpsLR(yi([rb,rb-1],:), d(rb-1), alpha([rb,rb-1]), beta); |
| 168 | + else |
| 169 | + % perform fast energy update and store current state in R, z |
| 170 | + [eps_lr, R, z] = updateEpsLR(eps_lr, R, yi(lb,:), d(lb), z, alpha(lb), beta); |
| 171 | + end |
| 172 | + |
| 173 | + % pruning to skip unreachable configurations |
| 174 | + % (if this condition is met the following if-condition cannot never |
| 175 | + % be fulfilled because eps_lr is monote increasing and F >= 0.) |
| 176 | + if (eps_lr + gamma) >= F(rb) |
| 177 | + break |
| 178 | + end |
| 179 | + |
| 180 | + % check if setting a discontinuity between lb-1 and lb gives a |
| 181 | + % better energy |
| 182 | + candidate_value = F(lb-1) + gamma + eps_lr; |
| 183 | + if candidate_value < F(rb) |
| 184 | + F(rb ) = candidate_value; |
| 185 | + blb = lb; |
| 186 | + end |
| 187 | + |
| 188 | + end |
| 189 | + % store the best left bound corresponding to the right bound rb |
| 190 | + partition( rb ) = blb-1; |
| 191 | + end |
| 192 | + %%% END MAIN LOOP |
| 193 | + |
| 194 | + |
| 195 | + end |
| 196 | + %%% END CSSD CASE |
| 197 | + |
| 198 | + %%% BEGIN RECONSTRUCTION |
| 199 | + % the discontinuity locations are coded in the array 'partition'. The |
| 200 | + % vector [partition(rb)+1:rb] gives the indices of between two |
| 201 | + % discontinuity locations. We start from behind with [partition(N)+1:N] and |
| 202 | + % successively compute the preceding intervals. |
| 203 | + rb = N; |
| 204 | + pp_cell = {}; |
| 205 | + interval_cell = {}; |
| 206 | + discont = []; |
| 207 | + upper_discont = xi(end) + 1; |
| 208 | + while rb > 0 |
| 209 | + % partition(rb) stores corresponding optimal left bound lb |
| 210 | + lb = partition(rb); |
| 211 | + if lb == 0 |
| 212 | + lower_discont = xi(1) - 1; |
| 213 | + else |
| 214 | + lower_discont = (xi(lb+1) + xi(lb)) /2; |
| 215 | + end |
| 216 | + interval = (lb+1) : rb; |
| 217 | + interval_cell{end+1} = interval; %#ok<AGROW> (runtime not critical in this part of the algorithm) |
| 218 | + if length(interval) == 1 % this case should happen rarely but may happen e.g. for data of uneven length and low gamma parameter |
| 219 | + ymtx = zeros(4,D); |
| 220 | + ymtx(:,D) = yi(interval, :); |
| 221 | + pp = ppmak([lower_discont, upper_discont], ymtx', D); |
| 222 | + else |
| 223 | + pp = csaps(xi(interval),yi(interval, :)', p, [], wi(interval)); |
| 224 | + pp = linext_pp(pp, lower_discont, upper_discont); |
| 225 | + pp = embed_pptocubic(pp); |
| 226 | + end |
| 227 | + pp_cell{end+1} = pp; %#ok<AGROW> (runtime not critical in this part of the algorithm) |
| 228 | + % continue with next right bound |
| 229 | + rb = lb; |
| 230 | + upper_discont = lower_discont; |
| 231 | + discont(end+1) = lower_discont; %#ok<AGROW> (runtime not critical in this part of the algorithm) |
| 232 | + end |
| 233 | + %%% END RECONSTRUCTION |
| 234 | + |
| 235 | + %%% BEGIN MAKE PP FORM |
| 236 | + pp_cell = flip(pp_cell); % the pp's were computed in reverse order which is fixed here |
| 237 | + interval_cell = flip(interval_cell); |
| 238 | + pp = merge_ppcell(pp_cell); |
| 239 | + %%% END MAKE PP FORM |
| 240 | + |
| 241 | + discont = flip(discont(1:end-1))'; % the discontinuities were computed in reverse order which is fixed here |
| 242 | + |
| 243 | +end |
| 244 | + |
| 245 | + |
| 246 | +%%% BEGIN SET OUTPUT |
| 247 | +output.pp = pp; |
| 248 | +output.discont = discont; |
| 249 | +output.interval_cell = interval_cell; |
| 250 | +output.pp_cell = pp_cell; |
| 251 | +output.discont_idx = zeros(numel(interval_cell)-1, 1); |
| 252 | +for i = 1:numel(output.discont_idx) |
| 253 | + output.discont_idx(i) = output.interval_cell{i}(end); |
| 254 | +end |
| 255 | + |
| 256 | +if isempty(xx) |
| 257 | + output.yy = []; |
| 258 | +else |
| 259 | + output.yy = ppval(pp, xx); |
| 260 | +end |
| 261 | + |
| 262 | +fun_cell = cell(numel(pp_cell),1); |
| 263 | +for i=1:numel(pp_cell) |
| 264 | + fun_cell{i} = @(xx) ppval(pp_cell{i}, xx); |
| 265 | +end |
| 266 | +output.pcw_fun = PcwFunReal([-Inf; discont(:); Inf], fun_cell); |
| 267 | +%%% END SET OUTPUT |
| 268 | + |
| 269 | +end |
| 270 | + |
| 271 | + |
0 commit comments