Mercurial > hg > octave-nkf

changeset 6721:01036667884a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-06-14 06:56:41 by dbateman]
author dbateman
date Thu, 14 Jun 2007 06:56:42 +0000
parents fa2f5d4e55db
children 5b09d433171c
files doc/ChangeLog doc/interpreter/Makefile.in doc/interpreter/interp.txi doc/interpreter/interpimages.m scripts/ChangeLog scripts/general/__splinen__.m scripts/general/interp1.m scripts/general/interpn.m scripts/polynomial/mkpp.m
diffstat 9 files changed, 171 insertions(+), 18 deletions(-) [+]
line wrap: on
line diff
--- a/doc/ChangeLog
+++ b/doc/ChangeLog
@@ -1,3 +1,12 @@
+2007-06-14  David Bateman  <dbateman@free.fr>
+
+	* interpreter/Makefile.in (SCRIPT_SORCES): add interimages.m
+	(INTERPIMAGES*): New variables. Add targets for them
+	(IMAGES_EPS, IMAGES_PDF, IMAGES_PNG): Add the INTERPIMAGES.
+	* interpreter/interpimages.m: New function
+	* interpreter/interp.txi: Add text about second derivation of
+	splines	and add figures.
+
 2007-06-12  David Bateman  <dbateman@free.fr>
 
 	* interpreter/numbers.txi: Document that 64-bit arithmetic is
--- a/doc/interpreter/Makefile.in
+++ b/doc/interpreter/Makefile.in
@@ -18,7 +18,7 @@
 INSTALL_PROGRAM = @INSTALL_PROGRAM@
 INSTALL_DATA = @INSTALL_DATA@
 
-SCRIPT_SOURCES = sparseimages.m
+SCRIPT_SOURCES = sparseimages.m interpimages.m
 
 EXAMPLE_FILES_NODIR = \
   addtwomatrices.cc \
@@ -43,16 +43,20 @@
 
 EXAMPLE_FILES = $(addprefix $(top_srcdir)/examples/, $(EXAMPLE_FILES_NODIR))
 
+INTERPIMAGES = interpft interpn interpderiv
+INTERPIMAGES_EPS = $(addsuffix .eps, $(INTERPIMAGES))
+INTERPIMAGES_PDF = $(addsuffix .pdf, $(INTERPIMAGES))
+INTERPIMAGES_PNG = $(addsuffix .png, $(INTERPIMAGES))
+
 SPARSEIMAGES_1 = gplot grid spmatrix spchol spcholperm
-
 SPARSEIMAGES_EPS = $(addsuffix .eps, $(SPARSEIMAGES_1))
 SPARSEIMAGES_PDF = $(addsuffix .pdf, $(SPARSEIMAGES_1))
 SPARSEIMAGES_PNG = $(addsuffix .png, $(SPARSEIMAGES_1))
 SPARSEIMAGES_TXT = $(addsuffix .txt, $(SPARSEIMAGES_1))
 
-IMAGES_EPS = $(SPARSEIMAGES_EPS)
-IMAGES_PDF = $(SPARSEIMAGES_PDF)
-IMAGES_PNG = $(SPARSEIMAGES_PNG)
+IMAGES_EPS = $(SPARSEIMAGES_EPS) $(INTERPIMAGES_EPS)
+IMAGES_PDF = $(SPARSEIMAGES_PDF) $(INTERPIMAGES_PDF)
+IMAGES_PNG = $(SPARSEIMAGES_PNG) $(INTERPIMAGES_PNG)
 IMAGES_TXT = $(SPARSEIMAGES_TXT)
 
 HTML_IMAGES_PNG = $(addprefix HTML/, $(IMAGES_PNG))
@@ -207,6 +211,9 @@
     --eval "$(notdir $(basename $<)) ('$(notdir $(basename $@))', '$(patsubst .%,%, $(suffix $@))'); sleep (1);"
 endef
 
+$(INTERPIMAGES_EPS) $(INTERPIMAGES_PNG) $(INTERPIMAGES_TXT): interpimages.m
+	$(run-octave)
+
 $(SPARSEIMAGES_EPS) $(SPARSEIMAGES_PNG) $(SPARSEIMAGES_TXT): sparseimages.m
 	$(run-octave)
 
--- a/doc/interpreter/interp.txi
+++ b/doc/interpreter/interp.txi
@@ -15,6 +15,43 @@
 
 @DOCSTRING(interp1)
 
+There are some important differences between the various interpolation
+methods. The 'spline' method enforces that both the first and second
+derivatives of the interpolated values have a continuous derivative,
+whereas the other methods do not. This can be demonstrated by the code
+
+@example
+@group
+t = 0 : 0.3 : pi; dt = t(2)-t(1);
+n = length (t); k = 100; dti = dt*n/k;
+ti = t(1) + [0 : k-1]*dti;
+y = sin (4*t + 0.3) .* cos (3*t - 0.1);
+ddyc = diff(diff(interp1(t,y,ti,'cubic'))./dti)./dti;
+ddys = diff(diff(interp1(t,y,ti,'spline'))./dti)./dti;
+ddyp = diff(diff(interp1(t,y,ti,'pchip'))./dti)./dti;
+plot (ti(2:end-1), ddyc,'g+',ti(2:end-1),ddys,'b*', ...
+      ti(2:end-1),ddyp,'c^');
+legend('cubic','spline','pchip');
+@end group
+@end example
+
+@ifnotinfo
+@noindent
+The result of which can be seen in @ref{fig:interpderiv}.
+
+@float Figure,fig:interpderiv
+@image{interpderiv,8cm}
+@caption{Comparison of second derivative of interpolated values for
+various interpolation methods}
+@end float
+@end ifnotinfo
+
+This means that in general the 'spline' method results in smooth
+results. If the function to be interpolated is in fact smooth, then
+'spline' will give excellent results. However, if the function to be
+evaluated is in some manner discontinuous, then 'cubic' or 'pchip'
+interpolation might give better results.
+
 Fourier interpolation, is a resampling technique where a signal is
 converted to the frequency domain, padded with zeros and then
 reconverted to the time domain.
@@ -40,8 +77,20 @@
 @end group
 @end example
 
+@ifinfo
 which demonstrates the poor behavior of Fourier interpolation for non
 periodic functions.
+@end ifinfo
+@ifnotinfo
+which demonstrates the poor behavior of Fourier interpolation for non
+periodic functions, as can be seen in @ref{fig:interpft}.
+
+@float Figure,fig:interpft
+@image{interpft,8cm}
+@caption{Comparison of @code{interp1} and @code{interpft} for non
+periodic data}
+@end float
+@end ifnotinfo
 
 In additional the support function @code{spline} and @code{lookup} that
 underlie the @code{interp1} function can be called directly.
@@ -81,11 +130,12 @@
 f = @@(x,y,z) x.^2 - y - z.^2;
 [xx, yy, zz] = meshgrid (x, y, z);
 v = f (xx,yy,zz);
-xi = yi = zi = -1:0.5:1;
+xi = yi = zi = -1:0.1:1;
 [xxi, yyi, zzi] = meshgrid (xi, yi, zi);
-vi = interp3(x, y, z, v, xxi, yyi, zzi);
+vi = interp3(x, y, z, v, xxi, yyi, zzi, 'spline');
 [xxi, yyi, zzi] = ndgrid (xi, yi, zi);
-vi2 = interpn(x, y, z, v, xxi, yyi, zzi);
+vi2 = interpn(x, y, z, v, xxi, yyi, zzi, 'spline');
+mesh (yi, zi, squeeze (vi2(1,:,:)));
 @end group
 @end example
 
@@ -93,6 +143,14 @@
 where @code{vi} and @code{vi2} are identical. The reversal of the
 dimensions is treated in the @code{meshgrid} and @code{ndgrid} functions
 respectively.
+@ifnotinfo
+The result of this code can be seen in @ref{fig:interpn}.
+
+@float Figure,fig:interpn
+@image{interpn,8cm}
+@caption{Demonstration of the use of @code{interpn}}
+@end float
+@end ifnotinfo
 
 In additional the support function @code{bicubic} that underlies the
 cubic interpolation of @code{interp2} function can be called directly.
new file mode 100644
--- /dev/null
+++ b/doc/interpreter/interpimages.m
@@ -0,0 +1,45 @@
+function interpimages (nm, typ)
+  bury_output ();
+  if (strcmp (nm, "interpft"))
+    t = 0 : 0.3 : pi; dt = t(2)-t(1);
+    n = length (t); k = 100;
+    ti = t(1) + [0 : k-1]*dt*n/k;
+    y = sin (4*t + 0.3) .* cos (3*t - 0.1);
+    yp = sin (4*ti + 0.3) .* cos (3*ti - 0.1);
+    plot (ti, yp, 'g', ti, interp1(t, y, ti, 'spline'), 'b', ...
+	  ti, interpft (y, k), 'c', t, y, 'r+');
+    legend ('sin(4t+0.3)cos(3t-0.1','spline','interpft','data');
+    print (strcat (nm, ".", typ), strcat ("-d", typ))
+  elseif (strcmp (nm, "interpn"))
+    x = y = z = -1:1;
+    f = @(x,y,z) x.^2 - y - z.^2;
+    [xx, yy, zz] = meshgrid (x, y, z);
+    v = f (xx,yy,zz);
+    xi = yi = zi = -1:0.1:1;
+    [xxi, yyi, zzi] = ndgrid (xi, yi, zi);
+    vi = interpn(x, y, z, v, xxi, yyi, zzi, 'spline');
+    mesh (yi, zi, squeeze (vi(1,:,:)));
+    print (strcat (nm, ".", typ), strcat ("-d", typ))
+  elseif (strcmp (nm, "interpderiv"))
+    t = 0 : 0.3 : pi; dt = t(2)-t(1);
+    n = length (t); k = 100; dti = dt*n/k;
+    ti = t(1) + [0 : k-1]*dti;
+    y = sin (4*t + 0.3) .* cos (3*t - 0.1);
+    ddyc = diff(diff(interp1(t,y,ti,'cubic'))./dti)./dti;
+    ddys = diff(diff(interp1(t,y,ti,'spline'))./dti)./dti;
+    ddyp = diff(diff(interp1(t,y,ti,'pchip'))./dti)./dti;
+    plot (ti(2:end-1), ddyc,'g+',ti(2:end-1),ddys,'b*', ...
+          ti(2:end-1),ddyp,'c^');
+    legend('cubic','spline','pchip');
+    print (strcat (nm, ".", typ), strcat ("-d", typ))
+  endif
+  bury_output ();  
+endfunction
+
+## Use this function before plotting commands and after every call to
+## print since print() resets output to stdout (unfortunately, gnpulot
+## can't pop output as it can the terminal type).
+function bury_output ()
+  f = figure (1);
+  set (f, "visible", "off");
+endfunction
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,13 @@
+2007-06-14  David Bateman  <dbateman@free.fr>
+
+	* general/__splinen__.m: Check also for ND vectors. Fix for N > 2,
+	as permutation of results was incorrect.
+	* general/interp1.m: Add demo on second derivative
+	* general/interpn.m: Convert "y" to vectors for __splinen__
+	call. Add 3D demo.
+
+	* polynomial/mkpp.m: Correction for matrices of 3 or more dimensions.
+
 2007-06-13  John W. Eaton  <jwe@octave.org>
 
 	* miscellaneous/mkoctfile.m: Quote args too.
--- a/scripts/general/__splinen__.m
+++ b/scripts/general/__splinen__.m
@@ -28,14 +28,14 @@
   if (nargin != 5)
     error ("Incorrect number of arguments");
   endif
-
-  if (!iscell (x) || length(x) < ndims(y) || any (! cellfun (@isvector, x)) ||
-      !iscell (xi) || length(xi) < ndims(y) || any (! cellfun (@isvector, xi)))
+  isvec = @(x) prod(size(x)) == length(x);   # ND isvector function
+  if (!iscell (x) || length(x) < ndims(y) || any (! cellfun (isvec, x)) ||
+      !iscell (xi) || length(xi) < ndims(y) || any (! cellfun (isvec, xi)))
     error ("%s: non gridded data or dimensions inconsistent", f);
   endif
   yi = y;
   for i = length(x):-1:1
-    yi = spline (x{i}, yi, xi{i}).';
+    yi = permute (spline (x{i}, yi, xi{i}), [length(x),1:length(x)-1]);
   endfor
 
   [xi{:}] = ndgrid (xi{:});
--- a/scripts/general/interp1.m
+++ b/scripts/general/interp1.m
@@ -345,6 +345,19 @@
 %! %--------------------------------------------------------
 %! % confirm that interpolated function matches the original
 
+%!demo
+%! t = 0 : 0.3 : pi; dt = t(2)-t(1);
+%! n = length (t); k = 100; dti = dt*n/k;
+%! ti = t(1) + [0 : k-1]*dti;
+%! y = sin (4*t + 0.3) .* cos (3*t - 0.1);
+%! ddyc = diff(diff(interp1(t,y,ti,'cubic'))./dti)./dti;
+%! ddys = diff(diff(interp1(t,y,ti,'spline'))./dti)./dti;
+%! ddyp = diff(diff(interp1(t,y,ti,'pchip'))./dti)./dti;
+%! plot (ti(2:end-1), ddyc,'g+',ti(2:end-1),ddys,'b*', ...
+%!       ti(2:end-1),ddyp,'c^');
+%! legend('cubic','spline','pchip');
+%! title("Second derivative of interpolated 'sin (4*t + 0.3) .* cos (3*t - 0.1)'");
+
 ## For each type of interpolated test, confirm that the interpolated
 ## value at the knots match the values at the knots.  Points away
 ## from the knots are requested, but only 'nearest' and 'linear'
--- a/scripts/general/interpn.m
+++ b/scripts/general/interpn.m
@@ -128,22 +128,22 @@
       endif
       idx (1 : nd) = {1};
       idx (i) = ":";
-      x{i} = x{i}(idx{:});
+      x{i} = x{i}(idx{:})(:);
     endfor
     idx (1 : nd) = {1};
     idx (1) = ":";
-    x{1} = x{1}(idx{:});
+    x{1} = x{1}(idx{:})(:);
   endif
 
   if (strcmp (method, "linear") || strcmp (method, "nearest"))
     if (all (cellfun (@isvector, y)))
       [y{:}] = ndgrid (y{:});
     endif
-  elseif (any (! cellfun (@isvector, x)))
+  elseif (any (! cellfun (@isvector, y)))
     for i = 1 : nd
       idx (1 : nd) = {1};
       idx (i) = ":";
-      y{i} = y{i}(idx{:});
+      y{i} = y{i}(idx{:})(:).';
     endfor
   endif
 
@@ -170,7 +170,7 @@
     endfor
     vi(idx) = extrapval;
     vi = reshape (vi, yshape); 
-  elseif (strcmp (method, "spline")) 
+  elseif (strcmp (method, "spline"))
     vi = __splinen__ (x, v, y, extrapval, "interpn");
   elseif (strcmp (method, "cubic")) 
     error ("cubic interpolation not yet implemented");
@@ -216,3 +216,14 @@
 %! [x,y] = meshgrid(x,y); 
 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 
+
+%!demo
+%! x = y = z = -1:1;
+%! f = @(x,y,z) x.^2 - y - z.^2;
+%! [xx, yy, zz] = meshgrid (x, y, z);
+%! v = f (xx,yy,zz);
+%! xi = yi = zi = -1:0.1:1;
+%! [xxi, yyi, zzi] = ndgrid (xi, yi, zi);
+%! vi = interpn(x, y, z, v, xxi, yyi, zzi, 'spline');
+%! mesh (yi, zi, squeeze (vi(1,:,:)));
+
--- a/scripts/polynomial/mkpp.m
+++ b/scripts/polynomial/mkpp.m
@@ -52,7 +52,7 @@
     d = round (rows (P) / pp.n); 
   endif
   pp.d = d;
-  if (pp.n*d != rows (P))
+  if (pp.n * prod (d) != rows (P))
     error ("mkpp: num intervals in x doesn't match num polynomials in P");
   endif
 endfunction