octave_packages/m/general/cplxpair.m

   1 ## Copyright (C) 2000-2012 Paul Kienzle
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {} cplxpair (@var{z})
  21 ## @deftypefnx {Function File} {} cplxpair (@var{z}, @var{tol})
  22 ## @deftypefnx {Function File} {} cplxpair (@var{z}, @var{tol}, @var{dim})
  23 ## Sort the numbers @var{z} into complex conjugate pairs ordered by
  24 ## increasing real part.  Place the negative imaginary complex number
  25 ## first within each pair.  Place all the real numbers (those with
  26 ## @code{abs (imag (@var{z}) / @var{z}) < @var{tol})}) after the
  27 ## complex pairs.
  28 ##
  29 ## If @var{tol} is unspecified the default value is 100*@code{eps}.
  30 ##
  31 ## By default the complex pairs are sorted along the first non-singleton
  32 ## dimension of @var{z}.  If @var{dim} is specified, then the complex
  33 ## pairs are sorted along this dimension.
  34 ##
  35 ## Signal an error if some complex numbers could not be paired.  Signal an
  36 ## error if all complex numbers are not exact conjugates (to within
  37 ## @var{tol}).  Note that there is no defined order for pairs with identical
  38 ## real parts but differing imaginary parts.
  39 ## @c Set example in small font to prevent overfull line
  40 ##
  41 ## @smallexample
  42 ## cplxpair (exp(2i*pi*[0:4]'/5)) == exp(2i*pi*[3; 2; 4; 1; 0]/5)
  43 ## @end smallexample
  44 ## @end deftypefn
  45
  46 ## FIXME: subsort returned pairs by imaginary magnitude
  47 ## FIXME: Why doesn't exp(2i*pi*[0:4]'/5) produce exact conjugates. Does
  48 ## FIXME:    it in Matlab?  The reason is that complex pairs are supposed
  49 ## FIXME:    to be exact conjugates, and not rely on a tolerance test.
  50
  51 ## 2006-05-12 David Bateman - Modified for NDArrays
  52
  53 function y = cplxpair (z, tol, dim)
  54
  55   if nargin < 1 || nargin > 3
  56     print_usage ();
  57   endif
  58
  59   if (length (z) == 0)
  60     y = zeros (size (z));
  61     return;
  62   endif
  63
  64   if (nargin < 2 || isempty (tol))
  65     if (isa (z, "single"))
  66       tol = 100 * eps("single");
  67     else
  68       tol = 100*eps;
  69     endif
  70   endif
  71
  72   nd = ndims (z);
  73   orig_dims = size (z);
  74   if (nargin < 3)
  75     ## Find the first singleton dimension.
  76     dim = 0;
  77     while (dim < nd && orig_dims(dim+1) == 1)
  78       dim++;
  79     endwhile
  80     dim++;
  81     if (dim > nd)
  82       dim = 1;
  83     endif
  84   else
  85     dim = floor(dim);
  86     if (dim < 1 || dim > nd)
  87       error ("cplxpair: invalid dimension along which to sort");
  88     endif
  89   endif
  90
  91   ## Move dimension to treat first, and convert to a 2-D matrix.
  92   perm = [dim:nd, 1:dim-1];
  93   z = permute (z, perm);
  94   sz = size (z);
  95   n = sz (1);
  96   m = prod (sz) / n;
  97   z = reshape (z, n, m);
  98
  99   ## Sort the sequence in terms of increasing real values.
 100   [q, idx] = sort (real (z), 1);
 101   z = z(idx + n * ones (n, 1) * [0:m-1]);
 102
 103   ## Put the purely real values at the end of the returned list.
 104   cls = "double";
 105   if (isa (z, "single"))
 106     cls = "single";
 107   endif
 108   [idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin(cls)) < tol);
 109   q = sparse (idxi, idxj, 1, n, m);
 110   nr = sum (q, 1);
 111   [q, idx] = sort (q, 1);
 112   z = z(idx);
 113   y = z;
 114
 115   ## For each remaining z, place the value and its conjugate at the
 116   ## start of the returned list, and remove them from further
 117   ## consideration.
 118   for j = 1:m
 119     p = n - nr(j);
 120     for i = 1:2:p
 121       if (i+1 > p)
 122         error ("cplxpair: could not pair all complex numbers");
 123       endif
 124       [v, idx] = min (abs (z(i+1:p) - conj (z(i))));
 125       if (v > tol)
 126         error ("cplxpair: could not pair all complex numbers");
 127       endif
 128       if (imag (z(i)) < 0)
 129         y([i, i+1]) = z([i, idx+i]);
 130       else
 131         y([i, i+1]) = z([idx+i, i]);
 132       endif
 133       z(idx+i) = z(i+1);
 134     endfor
 135   endfor
 136
 137   ## Reshape the output matrix.
 138   y = ipermute (reshape (y, sz), perm);
 139
 140 endfunction
 141
 142 %!demo
 143 %! [ cplxpair(exp(2i*pi*[0:4]'/5)), exp(2i*pi*[3; 2; 4; 1; 0]/5) ]
 144
 145 %!assert (isempty(cplxpair([])));
 146 %!assert (cplxpair(1), 1)
 147 %!assert (cplxpair([1+1i, 1-1i]), [1-1i, 1+1i])
 148 %!assert (cplxpair([1+1i, 1+1i, 1, 1-1i, 1-1i, 2]), \
 149 %!        [1-1i, 1+1i, 1-1i, 1+1i, 1, 2])
 150 %!assert (cplxpair([1+1i; 1+1i; 1; 1-1i; 1-1i; 2]), \
 151 %!        [1-1i; 1+1i; 1-1i; 1+1i; 1; 2])
 152 %!assert (cplxpair([0, 1, 2]), [0, 1, 2]);
 153
 154 %!shared z
 155 %! z=exp(2i*pi*[4; 3; 5; 2; 6; 1; 0]/7);
 156 %!assert (cplxpair(z(randperm(7))), z);
 157 %!assert (cplxpair(z(randperm(7))), z);
 158 %!assert (cplxpair(z(randperm(7))), z);
 159 %!assert (cplxpair([z(randperm(7)),z(randperm(7))]),[z,z])
 160 %!assert (cplxpair([z(randperm(7)),z(randperm(7))],[],1),[z,z])
 161 %!assert (cplxpair([z(randperm(7)).';z(randperm(7)).'],[],2),[z.';z.'])
 162
 163 %!## tolerance test
 164 %!assert (cplxpair([1i, -1i, 1+(1i*eps)],2*eps), [-1i, 1i, 1+(1i*eps)]);