Mercurial > octave
annotate scripts/polynomial/private/__splinefit__.m @ 19596:0e1f5a750d00
maint: Periodic merge of gui-release to default.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 20 Jan 2015 10:24:46 -0500 |
| parents | be8a12acb20a 446c46af4b42 |
| children | 83792dd9bcc1 |
| rev | line source |
|---|---|
|
14525
624dcb5e978f
maint: note reason __splinefit__.m is private instead of subfunction
John W. Eaton <jwe@octave.org>
parents:
14509
diff
changeset
|
1 ## This function is private because it is maintained by Jonas Lundgren |
|
624dcb5e978f
maint: note reason __splinefit__.m is private instead of subfunction
John W. Eaton <jwe@octave.org>
parents:
14509
diff
changeset
|
2 ## separtely from Octave. Please do not reformat to match Octave coding |
|
624dcb5e978f
maint: note reason __splinefit__.m is private instead of subfunction
John W. Eaton <jwe@octave.org>
parents:
14509
diff
changeset
|
3 ## conventions as that would make it harder to integrate upstream |
|
624dcb5e978f
maint: note reason __splinefit__.m is private instead of subfunction
John W. Eaton <jwe@octave.org>
parents:
14509
diff
changeset
|
4 ## changes. |
|
624dcb5e978f
maint: note reason __splinefit__.m is private instead of subfunction
John W. Eaton <jwe@octave.org>
parents:
14509
diff
changeset
|
5 |
| 14509 | 6 % Copyright (c) 2010, Jonas Lundgren |
| 7 % All rights reserved. | |
|
19593
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
8 % |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
9 % Redistribution and use in source and binary forms, with or without |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
10 % modification, are permitted provided that the following conditions are |
| 14509 | 11 % met: |
|
19593
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
12 % |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
13 % * Redistributions of source code must retain the above copyright |
| 14509 | 14 % notice, this list of conditions and the following disclaimer. |
|
19593
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
15 % * Redistributions in binary form must reproduce the above copyright |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
16 % notice, this list of conditions and the following disclaimer in |
| 14509 | 17 % the documentation and/or other materials provided with the distribution |
|
19593
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
18 % |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
19 % THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
20 % AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
21 % IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
22 % ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
23 % LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
24 % CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
25 % SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
26 % INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
27 % CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) |
|
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17199
diff
changeset
|
28 % ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
| 14509 | 29 % POSSIBILITY OF SUCH DAMAGE. |
| 30 function pp = __splinefit__(varargin) | |
| 31 %SPLINEFIT Fit a spline to noisy data. | |
| 32 % PP = SPLINEFIT(X,Y,BREAKS) fits a piecewise cubic spline with breaks | |
| 33 % (knots) BREAKS to the noisy data (X,Y). X is a vector and Y is a vector | |
| 34 % or an ND array. If Y is an ND array, then X(j) and Y(:,...,:,j) are | |
| 35 % matched. Use PPVAL to evaluate PP. | |
| 36 % | |
| 37 % PP = SPLINEFIT(X,Y,P) where P is a positive integer interpolates the | |
| 38 % breaks linearly from the sorted locations of X. P is the number of | |
| 39 % spline pieces and P+1 is the number of breaks. | |
| 40 % | |
| 41 % OPTIONAL INPUT | |
| 42 % Argument places 4 to 8 are reserved for optional input. | |
| 43 % These optional arguments can be given in any order: | |
| 44 % | |
| 45 % PP = SPLINEFIT(...,'p') applies periodic boundary conditions to | |
| 46 % the spline. The period length is MAX(BREAKS)-MIN(BREAKS). | |
| 47 % | |
| 48 % PP = SPLINEFIT(...,'r') uses robust fitting to reduce the influence | |
| 49 % from outlying data points. Three iterations of weighted least squares | |
| 50 % are performed. Weights are computed from previous residuals. | |
| 51 % | |
| 52 % PP = SPLINEFIT(...,BETA), where 0 < BETA < 1, sets the robust fitting | |
| 53 % parameter BETA and activates robust fitting ('r' can be omitted). | |
| 54 % Default is BETA = 1/2. BETA close to 0 gives all data equal weighting. | |
| 55 % Increase BETA to reduce the influence from outlying data. BETA close | |
| 56 % to 1 may cause instability or rank deficiency. | |
| 57 % | |
| 58 % PP = SPLINEFIT(...,N) sets the spline order to N. Default is a cubic | |
| 59 % spline with order N = 4. A spline with P pieces has P+N-1 degrees of | |
| 60 % freedom. With periodic boundary conditions the degrees of freedom are | |
| 61 % reduced to P. | |
| 62 % | |
| 63 % PP = SPLINEFIT(...,CON) applies linear constraints to the spline. | |
| 64 % CON is a structure with fields 'xc', 'yc' and 'cc': | |
| 65 % 'xc', x-locations (vector) | |
| 66 % 'yc', y-values (vector or ND array) | |
| 67 % 'cc', coefficients (matrix). | |
| 68 % | |
| 69 % Constraints are linear combinations of derivatives of order 0 to N-2 | |
| 70 % according to | |
| 71 % | |
| 72 % cc(1,j)*y(x) + cc(2,j)*y'(x) + ... = yc(:,...,:,j), x = xc(j). | |
| 73 % | |
| 74 % The maximum number of rows for 'cc' is N-1. If omitted or empty 'cc' | |
| 75 % defaults to a single row of ones. Default for 'yc' is a zero array. | |
| 76 % | |
| 77 % EXAMPLES | |
| 78 % | |
| 79 % % Noisy data | |
| 80 % x = linspace(0,2*pi,100); | |
| 81 % y = sin(x) + 0.1*randn(size(x)); | |
| 82 % % Breaks | |
| 83 % breaks = [0:5,2*pi]; | |
| 84 % | |
| 85 % % Fit a spline of order 5 | |
| 86 % pp = splinefit(x,y,breaks,5); | |
| 87 % | |
| 88 % % Fit a spline of order 3 with periodic boundary conditions | |
| 89 % pp = splinefit(x,y,breaks,3,'p'); | |
| 90 % | |
| 91 % % Constraints: y(0) = 0, y'(0) = 1 and y(3) + y"(3) = 0 | |
| 92 % xc = [0 0 3]; | |
| 93 % yc = [0 1 0]; | |
| 94 % cc = [1 0 1; 0 1 0; 0 0 1]; | |
| 95 % con = struct('xc',xc,'yc',yc,'cc',cc); | |
| 96 % | |
| 97 % % Fit a cubic spline with 8 pieces and constraints | |
| 98 % pp = splinefit(x,y,8,con); | |
| 99 % | |
| 100 % % Fit a spline of order 6 with constraints and periodicity | |
| 101 % pp = splinefit(x,y,breaks,con,6,'p'); | |
| 102 % | |
| 103 % See also SPLINE, PPVAL, PPDIFF, PPINT | |
| 104 | |
| 105 % Author: Jonas Lundgren <splinefit@gmail.com> 2010 | |
| 106 | |
| 107 % 2009-05-06 Original SPLINEFIT. | |
| 108 % 2010-06-23 New version of SPLINEFIT based on B-splines. | |
| 109 % 2010-09-01 Robust fitting scheme added. | |
| 110 % 2010-09-01 Support for data containing NaNs. | |
| 111 % 2011-07-01 Robust fitting parameter added. | |
| 112 | |
| 113 % Check number of arguments | |
|
19068
be8a12acb20a
Deprecate nargchk in favor of narginchk.
Rik <rik@octave.org>
parents:
17199
diff
changeset
|
114 narginchk(3,7); |
| 14509 | 115 |
| 116 % Check arguments | |
| 117 [x,y,dim,breaks,n,periodic,beta,constr] = arguments(varargin{:}); | |
| 118 | |
| 119 % Evaluate B-splines | |
| 120 base = splinebase(breaks,n); | |
| 121 pieces = base.pieces; | |
| 122 A = ppval(base,x); | |
| 123 | |
| 124 % Bin data | |
| 125 [junk,ibin] = histc(x,[-inf,breaks(2:end-1),inf]); %#ok | |
| 126 | |
| 127 % Sparse system matrix | |
| 128 mx = numel(x); | |
| 129 ii = [ibin; ones(n-1,mx)]; | |
| 130 ii = cumsum(ii,1); | |
| 131 jj = repmat(1:mx,n,1); | |
| 132 if periodic | |
| 133 ii = mod(ii-1,pieces) + 1; | |
| 134 A = sparse(ii,jj,A,pieces,mx); | |
| 135 else | |
| 136 A = sparse(ii,jj,A,pieces+n-1,mx); | |
| 137 end | |
| 138 | |
| 139 % Don't use the sparse solver for small problems | |
| 140 if pieces < 20*n/log(1.7*n) | |
| 141 A = full(A); | |
| 142 end | |
| 143 | |
| 144 % Solve | |
| 145 if isempty(constr) | |
| 146 % Solve Min norm(u*A-y) | |
| 147 u = lsqsolve(A,y,beta); | |
| 148 else | |
| 149 % Evaluate constraints | |
| 150 B = evalcon(base,constr,periodic); | |
| 151 % Solve constraints | |
| 152 [Z,u0] = solvecon(B,constr); | |
| 153 % Solve Min norm(u*A-y), subject to u*B = yc | |
| 154 y = y - u0*A; | |
| 155 A = Z*A; | |
| 156 v = lsqsolve(A,y,beta); | |
| 157 u = u0 + v*Z; | |
| 158 end | |
| 159 | |
| 160 % Periodic expansion of solution | |
| 161 if periodic | |
| 162 jj = mod(0:pieces+n-2,pieces) + 1; | |
| 163 u = u(:,jj); | |
| 164 end | |
| 165 | |
| 166 % Compute polynomial coefficients | |
| 167 ii = [repmat(1:pieces,1,n); ones(n-1,n*pieces)]; | |
| 168 ii = cumsum(ii,1); | |
| 169 jj = repmat(1:n*pieces,n,1); | |
| 170 C = sparse(ii,jj,base.coefs,pieces+n-1,n*pieces); | |
| 171 coefs = u*C; | |
| 172 coefs = reshape(coefs,[],n); | |
| 173 | |
| 174 % Make piecewise polynomial | |
| 175 pp = mkpp(breaks,coefs,dim); | |
| 176 | |
| 177 | |
| 178 %-------------------------------------------------------------------------- | |
| 179 function [x,y,dim,breaks,n,periodic,beta,constr] = arguments(varargin) | |
| 180 %ARGUMENTS Lengthy input checking | |
| 181 % x Noisy data x-locations (1 x mx) | |
| 182 % y Noisy data y-values (prod(dim) x mx) | |
| 183 % dim Leading dimensions of y | |
| 184 % breaks Breaks (1 x (pieces+1)) | |
| 185 % n Spline order | |
| 186 % periodic True if periodic boundary conditions | |
| 187 % beta Robust fitting parameter, no robust fitting if beta = 0 | |
| 188 % constr Constraint structure | |
| 189 % constr.xc x-locations (1 x nx) | |
| 190 % constr.yc y-values (prod(dim) x nx) | |
| 191 % constr.cc Coefficients (?? x nx) | |
| 192 | |
| 193 % Reshape x-data | |
| 194 x = varargin{1}; | |
| 195 mx = numel(x); | |
| 196 x = reshape(x,1,mx); | |
| 197 | |
| 198 % Remove trailing singleton dimensions from y | |
| 199 y = varargin{2}; | |
| 200 dim = size(y); | |
| 201 while numel(dim) > 1 && dim(end) == 1 | |
| 202 dim(end) = []; | |
| 203 end | |
| 204 my = dim(end); | |
| 205 | |
| 206 % Leading dimensions of y | |
| 207 if numel(dim) > 1 | |
| 208 dim(end) = []; | |
| 209 else | |
| 210 dim = 1; | |
| 211 end | |
| 212 | |
| 213 % Reshape y-data | |
| 214 pdim = prod(dim); | |
| 215 y = reshape(y,pdim,my); | |
| 216 | |
| 217 % Check data size | |
| 218 if mx ~= my | |
| 219 mess = 'Last dimension of array y must equal length of vector x.'; | |
| 220 error('arguments:datasize',mess) | |
| 221 end | |
| 222 | |
| 223 % Treat NaNs in x-data | |
| 224 inan = find(isnan(x)); | |
| 225 if ~isempty(inan) | |
| 226 x(inan) = []; | |
| 227 y(:,inan) = []; | |
| 228 mess = 'All data points with NaN as x-location will be ignored.'; | |
| 229 warning('arguments:nanx',mess) | |
| 230 end | |
| 231 | |
| 232 % Treat NaNs in y-data | |
| 233 inan = find(any(isnan(y),1)); | |
| 234 if ~isempty(inan) | |
| 235 x(inan) = []; | |
| 236 y(:,inan) = []; | |
| 237 mess = 'All data points with NaN in their y-value will be ignored.'; | |
| 238 warning('arguments:nany',mess) | |
| 239 end | |
| 240 | |
| 241 % Check number of data points | |
| 242 mx = numel(x); | |
| 243 if mx == 0 | |
| 244 error('arguments:nodata','There must be at least one data point.') | |
| 245 end | |
| 246 | |
| 247 % Sort data | |
| 248 if any(diff(x) < 0) | |
| 249 [x,isort] = sort(x); | |
| 250 y = y(:,isort); | |
| 251 end | |
| 252 | |
| 253 % Breaks | |
| 254 if isscalar(varargin{3}) | |
| 255 % Number of pieces | |
| 256 p = varargin{3}; | |
| 257 if ~isreal(p) || ~isfinite(p) || p < 1 || fix(p) < p | |
| 258 mess = 'Argument #3 must be a vector or a positive integer.'; | |
| 259 error('arguments:breaks1',mess) | |
| 260 end | |
| 261 if x(1) < x(end) | |
| 262 % Interpolate breaks linearly from x-data | |
| 263 dx = diff(x); | |
| 264 ibreaks = linspace(1,mx,p+1); | |
| 265 [junk,ibin] = histc(ibreaks,[0,2:mx-1,mx+1]); %#ok | |
| 266 breaks = x(ibin) + dx(ibin).*(ibreaks-ibin); | |
| 267 else | |
| 268 breaks = x(1) + linspace(0,1,p+1); | |
| 269 end | |
| 270 else | |
| 271 % Vector of breaks | |
| 272 breaks = reshape(varargin{3},1,[]); | |
| 273 if isempty(breaks) || min(breaks) == max(breaks) | |
| 274 mess = 'At least two unique breaks are required.'; | |
| 275 error('arguments:breaks2',mess); | |
| 276 end | |
| 277 end | |
| 278 | |
| 279 % Unique breaks | |
| 280 if any(diff(breaks) <= 0) | |
| 281 breaks = unique(breaks); | |
| 282 end | |
| 283 | |
| 284 % Optional input defaults | |
| 285 n = 4; % Cubic splines | |
| 286 periodic = false; % No periodic boundaries | |
| 287 robust = false; % No robust fitting scheme | |
| 288 beta = 0.5; % Robust fitting parameter | |
| 289 constr = []; % No constraints | |
| 290 | |
| 291 % Loop over optional arguments | |
| 292 for k = 4:nargin | |
| 293 a = varargin{k}; | |
| 294 if ischar(a) && isscalar(a) && lower(a) == 'p' | |
| 295 % Periodic conditions | |
| 296 periodic = true; | |
| 297 elseif ischar(a) && isscalar(a) && lower(a) == 'r' | |
| 298 % Robust fitting scheme | |
| 299 robust = true; | |
| 300 elseif isreal(a) && isscalar(a) && isfinite(a) && a > 0 && a < 1 | |
| 301 % Robust fitting parameter | |
| 302 beta = a; | |
| 303 robust = true; | |
| 304 elseif isreal(a) && isscalar(a) && isfinite(a) && a > 0 && fix(a) == a | |
| 305 % Spline order | |
| 306 n = a; | |
| 307 elseif isstruct(a) && isscalar(a) | |
| 308 % Constraint structure | |
| 309 constr = a; | |
| 310 else | |
| 311 error('arguments:nonsense','Failed to interpret argument #%d.',k) | |
| 312 end | |
| 313 end | |
| 314 | |
| 315 % No robust fitting | |
| 316 if ~robust | |
| 317 beta = 0; | |
| 318 end | |
| 319 | |
| 320 % Check exterior data | |
| 321 h = diff(breaks); | |
| 322 xlim1 = breaks(1) - 0.01*h(1); | |
| 323 xlim2 = breaks(end) + 0.01*h(end); | |
| 324 if x(1) < xlim1 || x(end) > xlim2 | |
| 325 if periodic | |
| 326 % Move data inside domain | |
| 327 P = breaks(end) - breaks(1); | |
| 328 x = mod(x-breaks(1),P) + breaks(1); | |
| 329 % Sort | |
| 330 [x,isort] = sort(x); | |
| 331 y = y(:,isort); | |
| 332 else | |
| 333 mess = 'Some data points are outside the spline domain.'; | |
| 334 warning('arguments:exteriordata',mess) | |
| 335 end | |
| 336 end | |
| 337 | |
| 338 % Return | |
| 339 if isempty(constr) | |
| 340 return | |
| 341 end | |
| 342 | |
| 343 % Unpack constraints | |
| 344 xc = []; | |
| 345 yc = []; | |
| 346 cc = []; | |
| 347 names = fieldnames(constr); | |
| 348 for k = 1:numel(names) | |
| 349 switch names{k} | |
| 350 case {'xc'} | |
| 351 xc = constr.xc; | |
| 352 case {'yc'} | |
| 353 yc = constr.yc; | |
| 354 case {'cc'} | |
| 355 cc = constr.cc; | |
| 356 otherwise | |
| 357 mess = 'Unknown field ''%s'' in constraint structure.'; | |
| 358 warning('arguments:unknownfield',mess,names{k}) | |
| 359 end | |
| 360 end | |
| 361 | |
| 362 % Check xc | |
| 363 if isempty(xc) | |
| 364 mess = 'Constraints contains no x-locations.'; | |
| 365 error('arguments:emptyxc',mess) | |
| 366 else | |
| 367 nx = numel(xc); | |
| 368 xc = reshape(xc,1,nx); | |
| 369 end | |
| 370 | |
| 371 % Check yc | |
| 372 if isempty(yc) | |
| 373 % Zero array | |
| 374 yc = zeros(pdim,nx); | |
| 375 elseif numel(yc) == 1 | |
| 376 % Constant array | |
| 377 yc = zeros(pdim,nx) + yc; | |
| 378 elseif numel(yc) ~= pdim*nx | |
| 379 % Malformed array | |
| 380 error('arguments:ycsize','Cannot reshape yc to size %dx%d.',pdim,nx) | |
| 381 else | |
| 382 % Reshape array | |
| 383 yc = reshape(yc,pdim,nx); | |
| 384 end | |
| 385 | |
| 386 % Check cc | |
| 387 if isempty(cc) | |
| 388 cc = ones(size(xc)); | |
| 389 elseif numel(size(cc)) ~= 2 | |
| 17199 | 390 error('arguments:ccsize1','Constraint coefficients cc must be 2-D.') |
| 14509 | 391 elseif size(cc,2) ~= nx |
| 392 mess = 'Last dimension of cc must equal length of xc.'; | |
| 393 error('arguments:ccsize2',mess) | |
| 394 end | |
| 395 | |
| 396 % Check high order derivatives | |
| 397 if size(cc,1) >= n | |
| 398 if any(any(cc(n:end,:))) | |
| 399 mess = 'Constraints involve derivatives of order %d or larger.'; | |
| 400 error('arguments:difforder',mess,n-1) | |
| 401 end | |
| 402 cc = cc(1:n-1,:); | |
| 403 end | |
| 404 | |
| 405 % Check exterior constraints | |
| 406 if min(xc) < xlim1 || max(xc) > xlim2 | |
| 407 if periodic | |
| 408 % Move constraints inside domain | |
| 409 P = breaks(end) - breaks(1); | |
| 410 xc = mod(xc-breaks(1),P) + breaks(1); | |
| 411 else | |
| 412 mess = 'Some constraints are outside the spline domain.'; | |
| 413 warning('arguments:exteriorconstr',mess) | |
| 414 end | |
| 415 end | |
| 416 | |
| 417 % Pack constraints | |
| 418 constr = struct('xc',xc,'yc',yc,'cc',cc); | |
| 419 | |
| 420 | |
| 421 %-------------------------------------------------------------------------- | |
| 422 function pp = splinebase(breaks,n) | |
| 423 %SPLINEBASE Generate B-spline base PP of order N for breaks BREAKS | |
| 424 | |
| 425 breaks = breaks(:); % Breaks | |
| 426 breaks0 = breaks'; % Initial breaks | |
| 427 h = diff(breaks); % Spacing | |
| 428 pieces = numel(h); % Number of pieces | |
| 429 deg = n - 1; % Polynomial degree | |
| 430 | |
| 431 % Extend breaks periodically | |
| 432 if deg > 0 | |
| 433 if deg <= pieces | |
| 434 hcopy = h; | |
| 435 else | |
| 436 hcopy = repmat(h,ceil(deg/pieces),1); | |
| 437 end | |
| 438 % to the left | |
| 439 hl = hcopy(end:-1:end-deg+1); | |
| 440 bl = breaks(1) - cumsum(hl); | |
| 441 % and to the right | |
| 442 hr = hcopy(1:deg); | |
| 443 br = breaks(end) + cumsum(hr); | |
| 444 % Add breaks | |
| 445 breaks = [bl(deg:-1:1); breaks; br]; | |
| 446 h = diff(breaks); | |
| 447 pieces = numel(h); | |
| 448 end | |
| 449 | |
| 450 % Initiate polynomial coefficients | |
| 451 coefs = zeros(n*pieces,n); | |
| 452 coefs(1:n:end,1) = 1; | |
| 453 | |
| 454 % Expand h | |
| 455 ii = [1:pieces; ones(deg,pieces)]; | |
| 456 ii = cumsum(ii,1); | |
| 457 ii = min(ii,pieces); | |
| 458 H = h(ii(:)); | |
| 459 | |
| 460 % Recursive generation of B-splines | |
| 461 for k = 2:n | |
| 462 % Antiderivatives of splines | |
| 463 for j = 1:k-1 | |
| 464 coefs(:,j) = coefs(:,j).*H/(k-j); | |
| 465 end | |
| 466 Q = sum(coefs,2); | |
| 467 Q = reshape(Q,n,pieces); | |
| 468 Q = cumsum(Q,1); | |
| 469 c0 = [zeros(1,pieces); Q(1:deg,:)]; | |
| 470 coefs(:,k) = c0(:); | |
| 471 % Normalize antiderivatives by max value | |
| 472 fmax = repmat(Q(n,:),n,1); | |
| 473 fmax = fmax(:); | |
| 474 for j = 1:k | |
| 475 coefs(:,j) = coefs(:,j)./fmax; | |
| 476 end | |
| 477 % Diff of adjacent antiderivatives | |
| 478 coefs(1:end-deg,1:k) = coefs(1:end-deg,1:k) - coefs(n:end,1:k); | |
| 479 coefs(1:n:end,k) = 0; | |
| 480 end | |
| 481 | |
| 482 % Scale coefficients | |
| 483 scale = ones(size(H)); | |
| 484 for k = 1:n-1 | |
| 485 scale = scale./H; | |
| 486 coefs(:,n-k) = scale.*coefs(:,n-k); | |
| 487 end | |
| 488 | |
| 489 % Reduce number of pieces | |
| 490 pieces = pieces - 2*deg; | |
| 491 | |
| 492 % Sort coefficients by interval number | |
| 493 ii = [n*(1:pieces); deg*ones(deg,pieces)]; | |
| 494 ii = cumsum(ii,1); | |
| 495 coefs = coefs(ii(:),:); | |
| 496 | |
| 497 % Make piecewise polynomial | |
| 498 pp = mkpp(breaks0,coefs,n); | |
| 499 | |
| 500 | |
| 501 %-------------------------------------------------------------------------- | |
| 502 function B = evalcon(base,constr,periodic) | |
| 503 %EVALCON Evaluate linear constraints | |
| 504 | |
| 505 % Unpack structures | |
| 506 breaks = base.breaks; | |
| 507 pieces = base.pieces; | |
| 508 n = base.order; | |
| 509 xc = constr.xc; | |
| 510 cc = constr.cc; | |
| 511 | |
| 512 % Bin data | |
| 513 [junk,ibin] = histc(xc,[-inf,breaks(2:end-1),inf]); %#ok | |
| 514 | |
| 515 % Evaluate constraints | |
| 516 nx = numel(xc); | |
| 517 B0 = zeros(n,nx); | |
| 518 for k = 1:size(cc,1) | |
| 519 if any(cc(k,:)) | |
| 520 B0 = B0 + repmat(cc(k,:),n,1).*ppval(base,xc); | |
| 521 end | |
| 522 % Differentiate base | |
| 523 coefs = base.coefs(:,1:n-k); | |
| 524 for j = 1:n-k-1 | |
| 525 coefs(:,j) = (n-k-j+1)*coefs(:,j); | |
| 526 end | |
| 527 base.coefs = coefs; | |
| 528 base.order = n-k; | |
| 529 end | |
| 530 | |
| 531 % Sparse output | |
| 532 ii = [ibin; ones(n-1,nx)]; | |
| 533 ii = cumsum(ii,1); | |
| 534 jj = repmat(1:nx,n,1); | |
| 535 if periodic | |
| 536 ii = mod(ii-1,pieces) + 1; | |
| 537 B = sparse(ii,jj,B0,pieces,nx); | |
| 538 else | |
| 539 B = sparse(ii,jj,B0,pieces+n-1,nx); | |
| 540 end | |
| 541 | |
| 542 | |
| 543 %-------------------------------------------------------------------------- | |
| 544 function [Z,u0] = solvecon(B,constr) | |
| 545 %SOLVECON Find a particular solution u0 and null space Z (Z*B = 0) | |
| 546 % for constraint equation u*B = yc. | |
| 547 | |
| 548 yc = constr.yc; | |
| 549 tol = 1000*eps; | |
| 550 | |
| 551 % Remove blank rows | |
| 552 ii = any(B,2); | |
| 553 B2 = full(B(ii,:)); | |
| 554 | |
| 555 % Null space of B2 | |
| 556 if isempty(B2) | |
| 557 Z2 = []; | |
| 558 else | |
| 559 % QR decomposition with column permutation | |
| 560 [Q,R,dummy] = qr(B2); %#ok | |
| 561 R = abs(R); | |
| 562 jj = all(R < R(1)*tol, 2); | |
| 563 Z2 = Q(:,jj)'; | |
| 564 end | |
| 565 | |
| 566 % Sizes | |
| 567 [m,ncon] = size(B); | |
| 568 m2 = size(B2,1); | |
| 569 nz = size(Z2,1); | |
| 570 | |
| 571 % Sparse null space of B | |
| 572 Z = sparse(nz+1:nz+m-m2,find(~ii),1,nz+m-m2,m); | |
| 573 Z(1:nz,ii) = Z2; | |
| 574 | |
| 575 % Warning rank deficient | |
| 576 if nz + ncon > m2 | |
| 577 mess = 'Rank deficient constraints, rank = %d.'; | |
| 578 warning('solvecon:deficient',mess,m2-nz); | |
| 579 end | |
| 580 | |
| 581 % Particular solution | |
| 582 u0 = zeros(size(yc,1),m); | |
| 583 if any(yc(:)) | |
| 584 % Non-homogeneous case | |
| 585 u0(:,ii) = yc/B2; | |
| 586 % Check solution | |
| 587 if norm(u0*B - yc,'fro') > norm(yc,'fro')*tol | |
| 588 mess = 'Inconsistent constraints. No solution within tolerance.'; | |
| 589 error('solvecon:inconsistent',mess) | |
| 590 end | |
| 591 end | |
| 592 | |
| 593 | |
| 594 %-------------------------------------------------------------------------- | |
| 595 function u = lsqsolve(A,y,beta) | |
| 596 %LSQSOLVE Solve Min norm(u*A-y) | |
| 597 | |
| 598 % Avoid sparse-complex limitations | |
| 599 if issparse(A) && ~isreal(y) | |
| 600 A = full(A); | |
| 601 end | |
| 602 | |
| 603 % Solution | |
| 604 u = y/A; | |
| 605 | |
| 606 % Robust fitting | |
| 607 if beta > 0 | |
| 608 [m,n] = size(y); | |
| 609 alpha = 0.5*beta/(1-beta)/m; | |
| 610 for k = 1:3 | |
| 611 % Residual | |
| 612 r = u*A - y; | |
| 613 rr = r.*conj(r); | |
| 614 rrmean = sum(rr,2)/n; | |
| 615 rrmean(~rrmean) = 1; | |
| 616 rrhat = (alpha./rrmean)'*rr; | |
| 617 % Weights | |
| 618 w = exp(-rrhat); | |
| 619 spw = spdiags(w',0,n,n); | |
| 620 % Solve weighted problem | |
| 621 u = (y*spw)/(A*spw); | |
| 622 end | |
| 623 end | |
| 624 |