octave_packages/general-1.3.1/unvech.m

   1 ## Copyright (C) 2006 Michael Creel <michael.creel@uab.es>
   2 ## Copyright (C) 2009 Jaroslav Hajek <highegg@gmail.com>
   3 ## Copyright (C) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   4 ## Copyright (C) 2011 Carnë Draug <carandraug+dev@gmail.com>
   5 ##
   6 ## This program is free software; you can redistribute it and/or modify it under
   7 ## the terms of the GNU General Public License as published by the Free Software
   8 ## Foundation; either version 3 of the License, or (at your option) any later
   9 ## version.
  10 ##
  11 ## This program is distributed in the hope that it will be useful, but WITHOUT
  12 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  14 ## details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License along with
  17 ## this program; if not, see <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {@var{m} =} unvech (@var{v}, @var{scale})
  21 ## Performs the reverse of @code{vech} on the vector @var{v}.
  22 ##
  23 ## Given a Nx1 array @var{v} describing the lower triangular part of a
  24 ## matrix (as obtained from @code{vech}), it returns the full matrix.
  25 ##
  26 ## The upper triangular part of the matrix will be multiplied by @var{scale} such
  27 ## that 1 and -1 can be used for symmetric and antisymmetric matrix respectively.
  28 ## @var{scale} must be a scalar and defaults to 1.
  29 ##
  30 ## @seealso{vech, ind2sub}
  31 ## @end deftypefn
  32
  33 ## TODO remove subfunction ind2sub_tril after new release of octave that will have
  34 ## it builtin standard ind2sub
  35
  36 function M = unvech (v, scale = 1)
  37
  38   if ( nargin < 1 || nargin > 2 )
  39     print_usage;
  40   elseif ( !ismatrix (v) && any (size (v) != 1) )
  41     error ("V must be a row or column matrix")
  42   elseif ( !isnumeric (scale) || !isscalar (scale) )
  43     error ("SCALE must be a scalar")
  44   endif
  45
  46   N      = length (v);
  47   dim    = (sqrt ( 1 + 8*N ) - 1)/2;
  48   [r, c] = ind2sub_tril (dim, 1:N);   # replace with core ind2sub after octave 3.6
  49   M      = accumarray ([r; c].', v);
  50   M     += scale * tril (M, -1).';
  51
  52 endfunction
  53
  54 function [r c] = ind2sub_tril(N,idx)
  55
  56   endofrow = 0.5*(1:N) .* (2*N:-1:N + 1);
  57   c = lookup(endofrow, idx-1)+1;
  58
  59   r = N - endofrow(c) + idx ;
  60
  61 endfunction
  62
  63 %!assert(unvech([1;0;0;1;0;1]), full(eye(3,3)) );
  64
  65 %!test %symmetric
  66 %! dim = 10;
  67 %! A = tril( floor ( 5*(2*rand(dim)-1) ) );
  68 %! A += A.';
  69 %! M = vech(A);
  70 %! M = unvech(M, 1);
  71 %! assert (A, M);
  72
  73 %!test %antisymmetric
  74 %! dim = 10;
  75 %! A = tril( floor ( 5*(2*rand(dim)-1) ) );
  76 %! A -= A.';
  77 %! M = vech(A);
  78 %! M = unvech(M, -1);
  79 %! assert (A, M);