Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / linear-algebra-2.2.0 / ndcovlt.m

diff --git a/octave_packages/linear-algebra-2.2.0/ndcovlt.m b/octave_packages/linear-algebra-2.2.0/ndcovlt.m
new file mode 100644 (file)
index 0000000..e2bab20
--- /dev/null
+++ b/octave_packages/linear-algebra-2.2.0/ndcovlt.m
@@ -0,0 +1,99 @@
+## Copyright (C) 2010 VZLU Prague, a.s., Czech Republic
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn{Function File} {@var{y} =} ndcovlt (@var{x}, @var{t1}, @var{t2}, @dots{})
+## Computes an n-dimensional covariant linear transform of an n-d tensor, given a
+## transformation matrix for each dimension. The number of columns of each transformation
+## matrix must match the corresponding extent of @var{x}, and the number of rows determines
+## the corresponding extent of @var{y}. For example:
+## 
+## @example
+##   size (@var{x}, 2) == columns (@var{t2})
+##   size (@var{y}, 2) == rows (@var{t2})
+## @end example
+##
+## The element @code{@var{y}(i1, i2, @dots{})} is defined as a sum of
+##
+## @example
+##   @var{x}(j1, j2, @dots{}) * @var{t1}(i1, j1) * @var{t2}(i2, j2) * @dots{}
+## @end example
+##
+## over all j1, j2, @dots{}. For two dimensions, this reduces to
+## @example
+##   @var{y} = @var{t1} * @var{x} * @var{t2}.'
+## @end example
+## 
+## [] passed as a transformation matrix is converted to identity matrix for
+## the corresponding dimension.
+##
+## @end deftypefn
+
+## Author: Jaroslav Hajek <highegg@gmail.com>
+
+function y = ndcovlt (x, varargin)
+  nd = max (ndims (x), nargin - 1);
+  varargin = resize (varargin, 1, nd);
+
+  # check dimensions
+  for i = 1:nd
+    ti = varargin{i};
+    if (isnumeric (ti) && ndims (ti) == 2)
+      [r, c] = size (ti);
+      if (r + c == 0)
+        varargin{i} = eye (size (x, i));
+      elseif (c != size (x, i))
+        error ("ndcovt: dimension mismatch for x-th transformation matrix");
+      endif
+    else
+      error ("ndcovt: transformation matrices must be numeric 2d matrices");
+    endif
+  endfor
+
+  if (isempty (x))
+    szy = cellfun (@rows, varargin);
+    y = zeros (szy);
+    return
+  endif
+
+  ldp = [2:nd, 1];
+  ## First transformation.
+  y = ldtrans (x, varargin{1});
+
+  ## Always shift one dimension.
+  for i = 2:nd-1
+    y = ldtrans (permute (y, ldp), varargin{i});
+  endfor
+
+  ## Permute to normal order now to save one permutation.
+  if (nd > 2)
+    y = ipermute (y, [nd-1:nd, 1:nd-2]);
+  endif
+
+  ## Now multiply from the right.
+  szy = size (y);
+  szy(end+1:nd-1) = 1;
+  m = varargin{nd};
+  szy(nd) = rows (m);
+  y = reshape (y, [], size (y, nd));
+  y = reshape (y * m.', szy);
+
+endfunction
+
+function y = ldtrans (x, m)
+  sz = size (x);
+  sz(1) = rows (m);
+  y = reshape (m * x(:,:), sz);
+endfunction