octave_packages/linear-algebra-2.2.0/ndcovlt.m

   1 ## Copyright (C) 2010 VZLU Prague, a.s., Czech Republic
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn{Function File} {@var{y} =} ndcovlt (@var{x}, @var{t1}, @var{t2}, @dots{})
  18 ## Computes an n-dimensional covariant linear transform of an n-d tensor, given a
  19 ## transformation matrix for each dimension. The number of columns of each transformation
  20 ## matrix must match the corresponding extent of @var{x}, and the number of rows determines
  21 ## the corresponding extent of @var{y}. For example:
  22 ##
  23 ## @example
  24 ##   size (@var{x}, 2) == columns (@var{t2})
  25 ##   size (@var{y}, 2) == rows (@var{t2})
  26 ## @end example
  27 ##
  28 ## The element @code{@var{y}(i1, i2, @dots{})} is defined as a sum of
  29 ##
  30 ## @example
  31 ##   @var{x}(j1, j2, @dots{}) * @var{t1}(i1, j1) * @var{t2}(i2, j2) * @dots{}
  32 ## @end example
  33 ##
  34 ## over all j1, j2, @dots{}. For two dimensions, this reduces to
  35 ## @example
  36 ##   @var{y} = @var{t1} * @var{x} * @var{t2}.'
  37 ## @end example
  38 ##
  39 ## [] passed as a transformation matrix is converted to identity matrix for
  40 ## the corresponding dimension.
  41 ##
  42 ## @end deftypefn
  43
  44 ## Author: Jaroslav Hajek <highegg@gmail.com>
  45
  46 function y = ndcovlt (x, varargin)
  47   nd = max (ndims (x), nargin - 1);
  48   varargin = resize (varargin, 1, nd);
  49
  50   # check dimensions
  51   for i = 1:nd
  52     ti = varargin{i};
  53     if (isnumeric (ti) && ndims (ti) == 2)
  54       [r, c] = size (ti);
  55       if (r + c == 0)
  56         varargin{i} = eye (size (x, i));
  57       elseif (c != size (x, i))
  58         error ("ndcovt: dimension mismatch for x-th transformation matrix");
  59       endif
  60     else
  61       error ("ndcovt: transformation matrices must be numeric 2d matrices");
  62     endif
  63   endfor
  64
  65   if (isempty (x))
  66     szy = cellfun (@rows, varargin);
  67     y = zeros (szy);
  68     return
  69   endif
  70
  71   ldp = [2:nd, 1];
  72   ## First transformation.
  73   y = ldtrans (x, varargin{1});
  74
  75   ## Always shift one dimension.
  76   for i = 2:nd-1
  77     y = ldtrans (permute (y, ldp), varargin{i});
  78   endfor
  79
  80   ## Permute to normal order now to save one permutation.
  81   if (nd > 2)
  82     y = ipermute (y, [nd-1:nd, 1:nd-2]);
  83   endif
  84
  85   ## Now multiply from the right.
  86   szy = size (y);
  87   szy(end+1:nd-1) = 1;
  88   m = varargin{nd};
  89   szy(nd) = rows (m);
  90   y = reshape (y, [], size (y, nd));
  91   y = reshape (y * m.', szy);
  92
  93 endfunction
  94
  95 function y = ldtrans (x, m)
  96   sz = size (x);
  97   sz(1) = rows (m);
  98   y = reshape (m * x(:,:), sz);
  99 endfunction