Mercurial > hg > octave-nkf
annotate scripts/general/cplxpair.m @ 9141:c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Split section into "Exponents and Logarithms" and "Utility Functions"
Use Tex in many more of the doc strings for pretty printing in pdf format.
| author | Rik <rdrider0-list@yahoo.com> |
|---|---|
| date | Mon, 20 Apr 2009 17:16:09 -0700 |
| parents | 1bf0ce0930be |
| children | 5247e89688e1 |
| rev | line source |
|---|---|
| 8920 | 1 ## Copyright (C) 2000, 2006, 2007, 2008, 2009 Paul Kienzle |
| 5820 | 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 | |
| 7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
| 8 ## your option) any later version. | |
| 5820 | 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 | |
| 7016 | 16 ## along with Octave; see the file COPYING. If not, see |
| 17 ## <http://www.gnu.org/licenses/>. | |
| 5820 | 18 |
| 19 ## -*- texinfo -*- | |
|
9141
c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9051
diff
changeset
|
20 ## @deftypefn {Function File} {} cplxpair (@var{z}) |
|
c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9051
diff
changeset
|
21 ## @deftypefnx {Function File} {} cplxpair (@var{z}, @var{tol}) |
|
c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9051
diff
changeset
|
22 ## @deftypefnx {Function File} {} cplxpair (@var{z}, @var{tol}, @var{dim}) |
| 5820 | 23 ## Sort the numbers @var{z} into complex conjugate pairs ordered by |
|
9141
c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9051
diff
changeset
|
24 ## increasing real part. Place the negative imaginary complex number |
|
c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9051
diff
changeset
|
25 ## first within each pair. Place all the real numbers (those with |
|
c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9051
diff
changeset
|
26 ## @code{abs (imag (@var{z}) / @var{z}) < @var{tol})}) after the |
|
c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9051
diff
changeset
|
27 ## complex pairs. |
|
c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9051
diff
changeset
|
28 ## |
|
c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9051
diff
changeset
|
29 ## If @var{tol} is unspecified the default value is 100*@var{eps}. |
| 5820 | 30 ## |
| 31 ## By default the complex pairs are sorted along the first non-singleton | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
32 ## dimension of @var{z}. If @var{dim} is specified, then the complex |
| 5820 | 33 ## pairs are sorted along this dimension. |
| 34 ## | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
35 ## Signal an error if some complex numbers could not be paired. Requires |
|
9141
c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9051
diff
changeset
|
36 ## all complex numbers to be exact conjugates within @var{tol}, or signals an |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
37 ## error. Note that there are no guarantees on the order of the returned |
| 5820 | 38 ## pairs with identical real parts but differing imaginary parts. |
| 39 ## | |
| 7031 | 40 ## @c Using 'smallexample' to make text fit in page when using 'smallbook' |
| 41 ## @smallexample | |
| 5831 | 42 ## cplxpair (exp(2i*pi*[0:4]'/5)) == exp(2i*pi*[3; 2; 4; 1; 0]/5) |
| 7031 | 43 ## @end smallexample |
| 5820 | 44 ## @end deftypefn |
| 45 | |
|
8202
cf59d542f33e
replace all TODOs and XXXs with FIXMEs
Jaroslav Hajek <highegg@gmail.com>
parents:
7795
diff
changeset
|
46 ## FIXME: subsort returned pairs by imaginary magnitude |
|
cf59d542f33e
replace all TODOs and XXXs with FIXMEs
Jaroslav Hajek <highegg@gmail.com>
parents:
7795
diff
changeset
|
47 ## FIXME: Why doesn't exp(2i*pi*[0:4]'/5) produce exact conjugates. Does |
|
cf59d542f33e
replace all TODOs and XXXs with FIXMEs
Jaroslav Hajek <highegg@gmail.com>
parents:
7795
diff
changeset
|
48 ## FIXME: it in Matlab? The reason is that complex pairs are supposed |
|
cf59d542f33e
replace all TODOs and XXXs with FIXMEs
Jaroslav Hajek <highegg@gmail.com>
parents:
7795
diff
changeset
|
49 ## FIXME: to be exact conjugates, and not rely on a tolerance test. |
| 5820 | 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 | |
| 6046 | 56 print_usage (); |
| 5820 | 57 endif |
| 58 | |
| 59 if (length (z) == 0) | |
| 60 y = zeros (size (z)); | |
| 61 return; | |
| 62 endif | |
| 63 | |
| 64 if (nargin < 2 || isempty (tol)) | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
65 if (isa (z, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
66 tol = 100 * eps("single"); |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
67 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
68 tol = 100*eps; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
69 endif |
| 5820 | 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 | |
| 8507 | 91 ## Move dimension to treat first, and convert to a 2-D matrix. |
| 5820 | 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 | |
| 8507 | 99 ## Sort the sequence in terms of increasing real values. |
| 5820 | 100 [q, idx] = sort (real (z), 1); |
| 101 z = z(idx + n * ones (n, 1) * [0:m-1]); | |
| 102 | |
| 8507 | 103 ## Put the purely real values at the end of the returned list. |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
104 cls = "double"; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
105 if (isa (z, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
106 cls = "single"; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
107 endif |
| 8507 | 108 [idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin(cls)) < tol); |
| 5820 | 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); | |
| 8507 | 120 for i = 1:2:p |
| 5820 | 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 | |
| 8507 | 137 ## Reshape the output matrix. |
| 5820 | 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)]); |