Mercurial > forge

changeset 10583:cfc5f2dd5b4d octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
splines: replace deprecated calls to usage() by print_usage() and do not mix tabs with spaces identation
author carandraug
date Tue, 24 Jul 2012 15:32:47 +0000
parents 9ae794b00c25
children d7080cb246c8
files main/splines/inst/catmullrom.m main/splines/inst/csape.m main/splines/inst/csaps_sel.m main/splines/inst/fnder.m main/splines/inst/fnplt.m main/splines/inst/fnval.m
diffstat 6 files changed, 65 insertions(+), 68 deletions(-) [+]
line wrap: on
line diff
--- a/main/splines/inst/catmullrom.m	Tue Jul 24 15:21:36 2012 +0000
+++ b/main/splines/inst/catmullrom.m	Tue Jul 24 15:32:47 2012 +0000
@@ -53,7 +53,7 @@
     
   for ii = 1:4
     coeff(:,ii) =  ((h00(ii)*p0 + h10(ii)*h.*m0 +...
-		     h01(ii)*p1 + h11(ii)*h.*m1 )./h.^(4-ii))' ;
+                     h01(ii)*p1 + h11(ii)*h.*m1 )./h.^(4-ii))' ;
   end
 
   pp = mkpp (x, coeff);
--- a/main/splines/inst/csape.m	Tue Jul 24 15:21:36 2012 +0000
+++ b/main/splines/inst/csape.m	Tue Jul 24 15:32:47 2012 +0000
@@ -97,10 +97,10 @@
     end
 
     c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * valc(1) 
-	      - c(2,:) * h(1)) / (2 * h(1)); 
+              - c(2,:) * h(1)) / (2 * h(1)); 
     c(n,:) = - (3 / h(n - 1) * (a(n,:) - a(n - 1,:)) - 3 * valc(2) 
 
-		+ c(n - 1,:) * h(n - 1)) / (2 * h(n - 1));
+                + c(n - 1,:) * h(n - 1)) / (2 * h(n - 1));
     b(1:n - 1,:) = diff (a) ./ h(1:n - 1, idx)\
       - h(1:n - 1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:));
     d = diff (c) ./ (3 * h(1:n - 1, idx));
@@ -156,9 +156,9 @@
       dg(1) -=  gamma;
       dg(end) -= h(1) * h(1) / gamma; 
       z = spdiags([[e(:);0],dg,[0;e(:)]],[-1,0,1],n-1,n-1) \ ...
-	  [[gamma; zeros(n-3,1); h(1)],g];
+          [[gamma; zeros(n-3,1); h(1)],g];
       fact = (z(1,2:end) + h(1) * z(end,2:end) / gamma) / ...
-	  (1.0 + z(1,1) + h(1) * z(end,1) / gamma);
+          (1.0 + z(1,1) + h(1) * z(end,1) / gamma);
 
       c(2:n,idx) = z(:,2:end) - z(:,1) * fact;
     endif
@@ -176,8 +176,8 @@
     g(1,:) = 3 / (h(1) + h(2)) * (a(3,:) - a(2,:)\
           - h(2) / h(1) * (a(2,:) - a(1,:)));
     g(n - 2,:) = 3 / (h(n - 1) + h(n - 2)) *\
- 	(h(n - 2) / h(n - 1) * (a(n,:) - a(n - 1,:)) -\
-	 (a(n - 1,:) - a(n - 2,:)));
+        (h(n - 2) / h(n - 1) * (a(n,:) - a(n - 1,:)) -\
+         (a(n - 1,:) - a(n - 2,:)));
 
     if (n > 4)
 
@@ -201,7 +201,7 @@
       c(2:n - 1,:) = spdiags([[ldg(:);0],dg,[0;udg(:)]],[-1,0,1],n-2,n-2) \ g;
       
     else # n == 3
-	    
+            
       dg= [h(1) + 2 * h(2)];
       c(2:n - 1,:) = g/dg(1);
 
@@ -232,7 +232,7 @@
 %!assert (ppval(csape(x',y'),x'), y', 10*eps);
 %!assert (ppval(csape(x',y'),x), y, 10*eps);
 %!assert (ppval(csape(x,[y;y]),x), \
-%!	  [ppval(csape(x,y),x);ppval(csape(x,y),x)], 10*eps)
+%!        [ppval(csape(x,y),x);ppval(csape(x,y),x)], 10*eps)
 
 %!test cond='complete';
 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
@@ -240,7 +240,7 @@
 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
 %!assert (ppval(csape(x,[y;y],cond),x), \
-%!	  [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
+%!        [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
 
 %!test cond='variational';
 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
@@ -248,7 +248,7 @@
 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
 %!assert (ppval(csape(x,[y;y],cond),x), \
-%!	  [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
+%!        [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
 
 %!test cond='second';
 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
@@ -256,7 +256,7 @@
 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
 %!assert (ppval(csape(x,[y;y],cond),x), \
-%!	  [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
+%!        [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
 
 %!test cond='periodic';
 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
@@ -264,7 +264,7 @@
 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
 %!assert (ppval(csape(x,[y;y],cond),x), \
-%!	  [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
+%!        [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
 
 %!test cond='not-a-knot';
 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
@@ -272,4 +272,4 @@
 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
 %!assert (ppval(csape(x,[y;y],cond),x), \
-%!	  [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
+%!        [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
--- a/main/splines/inst/csaps_sel.m	Tue Jul 24 15:21:36 2012 +0000
+++ b/main/splines/inst/csaps_sel.m	Tue Jul 24 15:32:47 2012 +0000
@@ -85,7 +85,7 @@
 
   if isempty(crit)
     crit = 'aicc';
-  end  
+  end
 
   h = diff(x);
 
@@ -144,87 +144,85 @@
 
 
 function H = influence_matrix(p, U, D, n, w) #returns influence matrix for given p
-	H = speye(n) - U * diag(D ./ (D + (p / (6*(1-p))))) * U';
-	H = diag(1 ./ sqrt(w)) * H * diag(sqrt(w)); #rescale to original units	
+  H = speye(n) - U * diag(D ./ (D + (p / (6*(1-p))))) * U';
+  H = diag(1 ./ sqrt(w)) * H * diag(sqrt(w)); #rescale to original units	
 endfunction	
 
 function [MSR, Ht] = penalty_terms(H, y, w)
-	MSR = mean(w .* (y - (H*y)) .^ 2); #mean square residual
-	Ht = trace(H); #effective number of fitted parameters
+  MSR = mean(w .* (y - (H*y)) .^ 2); #mean square residual
+  Ht = trace(H); #effective number of fitted parameters
 endfunction
 
 function Hd = influence_matrix_diag_chol(p, QT, R, y, w, n)
-	#LDL factorization of 6*(1-p)*QT*diag(1 ./ w)*QT' + p*R
-	U = chol(6*(1-p)*QT*diag(1 ./ w)*QT' + p*R, 'upper');
-	d = 1 ./ diag(U);
-	U = diag(d)*U; 
-	d = d .^ 2;
-	#5 central bands in the inverse of 6*(1-p)*QT*diag(1 ./ w)*QT' + p*R
-	Binv = banded_matrix_inverse(d, U, 2);
-	Hd = diag(speye(n) - (6*(1-p))*diag(1 ./ w)*QT'*Binv*QT);	
+  #LDL factorization of 6*(1-p)*QT*diag(1 ./ w)*QT' + p*R
+  U = chol(6*(1-p)*QT*diag(1 ./ w)*QT' + p*R, 'upper');
+  d = 1 ./ diag(U);
+  U = diag(d)*U; 
+  d = d .^ 2;
+  #5 central bands in the inverse of 6*(1-p)*QT*diag(1 ./ w)*QT' + p*R
+  Binv = banded_matrix_inverse(d, U, 2);
+  Hd = diag(speye(n) - (6*(1-p))*diag(1 ./ w)*QT'*Binv*QT);	
 endfunction
 
 function [MSR, Ht] = penalty_terms_chol(p, QT, R, y, w, n)
-	#LDL factorization of 6*(1-p)*QT*diag(1 ./ w)*QT' + p*R
-	U = chol(6*(1-p)*QT*diag(1 ./ w)*QT' + p*R, 'upper');
-	d = 1 ./ diag(U);
-	U = diag(d)*U; 
-	d = d .^ 2;
-	Binv = banded_matrix_inverse(d, U, 2); #5 central bands in the inverse of 6*(1-p)*QT*diag(1 ./ w)*QT' + p*R
-	Ht = 2 + trace(speye(n-2) - (6*(1-p))*QT*diag(1 ./ w)*QT'*Binv);
-	MSR = mean(w .* ((6*(1-p)*diag(1 ./ w)*QT'*((6*(1-p)*QT*diag(1 ./ w)*QT' + p*R) \ (QT*y)))) .^ 2);
+  #LDL factorization of 6*(1-p)*QT*diag(1 ./ w)*QT' + p*R
+  U = chol(6*(1-p)*QT*diag(1 ./ w)*QT' + p*R, 'upper');
+  d = 1 ./ diag(U);
+  U = diag(d)*U; 
+  d = d .^ 2;
+  Binv = banded_matrix_inverse(d, U, 2); #5 central bands in the inverse of 6*(1-p)*QT*diag(1 ./ w)*QT' + p*R
+  Ht = 2 + trace(speye(n-2) - (6*(1-p))*QT*diag(1 ./ w)*QT'*Binv);
+  MSR = mean(w .* ((6*(1-p)*diag(1 ./ w)*QT'*((6*(1-p)*QT*diag(1 ./ w)*QT' + p*R) \ (QT*y)))) .^ 2);
 endfunction
 
 function J = aicc(MSR, Ht, n)
-	J = mean(log(MSR)(:)) + 2 * (Ht + 1) / max(n - Ht - 2, 0); #hurvich98, taking the average if there are multiple data sets as in woltring86 
+  J = mean(log(MSR)(:)) + 2 * (Ht + 1) / max(n - Ht - 2, 0); #hurvich98, taking the average if there are multiple data sets as in woltring86 
 endfunction
 
 function J = aic(MSR, Ht, n)
-	J = mean(log(MSR)(:)) + 2 * Ht / n;
+  J = mean(log(MSR)(:)) + 2 * Ht / n;
 endfunction
 
 function J = gcv(MSR, Ht, n)
-	J = mean(log(MSR)(:)) - 2 * log(1 - Ht / n);
+  J = mean(log(MSR)(:)) - 2 * log(1 - Ht / n);
 endfunction
 
 function J = msr_bound(MSR, Ht, n)
-	J = mean(MSR(:) - 1) .^ 2;
+  J = mean(MSR(:) - 1) .^ 2;
 endfunction
 
 function J = penalty_compute(p, U, D, y, w, n, crit) #evaluates a user-supplied penalty function crit at given p
-	H = influence_matrix(p, U, D, n, w);
-	[MSR, Ht] = penalty_terms(H, y, w);
-	J = feval(crit, MSR, Ht, n);
-	if ~isfinite(J)
-		J = Inf;
-	endif 	
+  H = influence_matrix(p, U, D, n, w);
+  [MSR, Ht] = penalty_terms(H, y, w);
+  J = feval(crit, MSR, Ht, n);
+  if ~isfinite(J)
+    J = Inf;
+  endif
 endfunction
 
 function J = penalty_compute_chol(p, QT, R, y, w, n, crit) #evaluates a user-supplied penalty function crit at given p
-	[MSR, Ht] = penalty_terms_chol(p, QT, R, y, w, n);
-	J = feval(crit, MSR, Ht, n);
-	if ~isfinite(J)
-		J = Inf;
-	endif 	
+  [MSR, Ht] = penalty_terms_chol(p, QT, R, y, w, n);
+  J = feval(crit, MSR, Ht, n);
+  if ~isfinite(J)
+    J = Inf;
+  endif
 endfunction
 
 function Binv = banded_matrix_inverse(d, U, m) #given a (2m+1)-banded, symmetric n x n matrix B = U'*inv(diag(d))*U, where U is unit upper triangular with bandwidth (m+1), returns Binv, a sparse symmetric matrix containing the central 2m+1 bands of the inverse of B
 #Reference: Hutchinson and de Hoog 1985
-	Binv = diag(d);
-	n = rows(U);
-	for i = n:(-1):1
-		p = min(m, n - i);
-		for l = 1:p
-			for k = 1:p
-				Binv(i, i+l) -= U(i, i+k)*Binv(i + k, i + l);
-			end
-			Binv(i, i) -= U(i, i+l)*Binv(i, i+l);
-		end
-		Binv(i+(1:p), i) = Binv(i, i+(1:p))'; #add the lower triangular elements
-	end
+  Binv = diag(d);
+  n = rows(U);
+  for i = n:(-1):1
+    p = min(m, n - i);
+    for l = 1:p
+      for k = 1:p
+        Binv(i, i+l) -= U(i, i+k)*Binv(i + k, i + l);
+      end
+      Binv(i, i) -= U(i, i+l)*Binv(i, i+l);
+    end
+    Binv(i+(1:p), i) = Binv(i, i+(1:p))'; #add the lower triangular elements
+  end
 endfunction
-			
-
 
 %!shared x,y,ret,p,sigma2,unc_y
 %! x = [0:0.01:1]'; y = sin(x);
@@ -251,4 +249,3 @@
 tic; [ret,p,sigma2,unc_y]=csaps_sel(x,y,x,w); toc
 
 %}
-
--- a/main/splines/inst/fnder.m	Tue Jul 24 15:21:36 2012 +0000
+++ b/main/splines/inst/fnder.m	Tue Jul 24 15:32:47 2012 +0000
@@ -26,7 +26,7 @@
 function dpp = fnder (pp, o)
 
   if (nargin < 1 || nargin > 2)
-    usage ("fnder (pp [, order])");
+    print_usage;
   endif
   if (nargin < 2)
     o = 1;
--- a/main/splines/inst/fnplt.m	Tue Jul 24 15:21:36 2012 +0000
+++ b/main/splines/inst/fnplt.m	Tue Jul 24 15:32:47 2012 +0000
@@ -31,7 +31,7 @@
 function [x, y] = fnplt (pp, plt)
 
   if (nargin < 1 || nargin > 2)
-    usage ("[x, y] = fnplt (pp [, plotstring])");
+    print_usage;
   endif
   if (nargin < 2)
     plt = "r;;";
--- a/main/splines/inst/fnval.m	Tue Jul 24 15:21:36 2012 +0000
+++ b/main/splines/inst/fnval.m	Tue Jul 24 15:32:47 2012 +0000
@@ -9,7 +9,7 @@
     # XXX FIXME XXX ignoring left continuous vs. right continuous option
     if isstruct(a), r=ppval(a,b); else r=ppval(b,a); end
   else
-    usage("r=fnval(pp,x) || r=fnval(x,pp)");
+    print_usage;
   end
 end