Mercurial > octave-nkf

comparison scripts/general/cplxpair.m @ 20508:4c2e76cbdc7d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
cplxpair.m: Use tolerance that depends on Z to evaluate pairing (bug #45810). * cplxpair.m: Update docstring to explain how tolerance input is calculated. Validate TOL input is a positive scalar. Use standard code to find the first singleton dimension. When determing whether a conjugate pair exists, use a tolerance of TOL*eps (abs (z(i))). Use Octave coding convestions in BIST tests. Add BIST tests for tolerance input. Add BIST tests for input validation.
author Rik <rik@octave.org>
date Tue, 01 Sep 2015 08:07:57 -0700
parents 7503499a252b
children
comparison
equal deleted inserted replaced
20507:2d415c68213f 20508:4c2e76cbdc7d
18 18
19 ## -*- texinfo -*- 19 ## -*- texinfo -*-
20 ## @deftypefn {Function File} {} cplxpair (@var{z}) 20 ## @deftypefn {Function File} {} cplxpair (@var{z})
21 ## @deftypefnx {Function File} {} cplxpair (@var{z}, @var{tol}) 21 ## @deftypefnx {Function File} {} cplxpair (@var{z}, @var{tol})
22 ## @deftypefnx {Function File} {} cplxpair (@var{z}, @var{tol}, @var{dim}) 22 ## @deftypefnx {Function File} {} cplxpair (@var{z}, @var{tol}, @var{dim})
23 ## Sort the numbers @var{z} into complex conjugate pairs ordered by 23 ## Sort the numbers @var{z} into complex conjugate pairs ordered by increasing
24 ## increasing real part. 24 ## real part.
25 ## 25 ##
26 ## The negative imaginary complex numbers are placed first within each pair. 26 ## The negative imaginary complex numbers are placed first within each pair.
27 ## All real numbers (those with 27 ## All real numbers (those with
28 ## @code{abs (imag (@var{z}) / @var{z}) < @var{tol}}) are placed after the 28 ## @code{abs (imag (@var{z}) / @var{z}) < @var{tol}}) are placed after the
29 ## complex pairs. 29 ## complex pairs.
30 ## 30 ##
31 ## If @var{tol} is unspecified the default value is 100*@code{eps}. 31 ## @var{tol} is a weighting factor which determines the tolerance of matching.
32 ## The default value is 100 and the resulting tolerance for a given complex
33 ## pair is @code{100 * eps (abs (@var{z}(i))}.
32 ## 34 ##
33 ## By default the complex pairs are sorted along the first non-singleton 35 ## By default the complex pairs are sorted along the first non-singleton
34 ## dimension of @var{z}. If @var{dim} is specified, then the complex pairs are 36 ## dimension of @var{z}. If @var{dim} is specified, then the complex pairs are
35 ## sorted along this dimension. 37 ## sorted along this dimension.
36 ## 38 ##
44 ## cplxpair (exp(2i*pi*[0:4]'/5)) == exp(2i*pi*[3; 2; 4; 1; 0]/5) 46 ## cplxpair (exp(2i*pi*[0:4]'/5)) == exp(2i*pi*[3; 2; 4; 1; 0]/5)
45 ## @end smallexample 47 ## @end smallexample
46 ## @end deftypefn 48 ## @end deftypefn
47 49
48 ## FIXME: subsort returned pairs by imaginary magnitude 50 ## FIXME: subsort returned pairs by imaginary magnitude
49 ## FIXME: Why doesn't exp (2i*pi*[0:4]'/5) produce exact conjugates. Does 51 ## FIXME: Why doesn't exp (2i*pi*[0:4]'/5) produce exact conjugates? Does
50 ## FIXME: it in Matlab? The reason is that complex pairs are supposed 52 ## FIXME: it in Matlab? The reason is that complex pairs are supposed
51 ## FIXME: to be exact conjugates, and not rely on a tolerance test. 53 ## FIXME: to be exact conjugates, and not rely on a tolerance test.
52 54
53 ## 2006-05-12 David Bateman - Modified for NDArrays 55 ## 2006-05-12 David Bateman - Modified for NDArrays
54 56
56 58
57 if (nargin < 1 || nargin > 3) 59 if (nargin < 1 || nargin > 3)
58 print_usage (); 60 print_usage ();
59 endif 61 endif
60 62
61 if (length (z) == 0) 63 if (isempty (z))
62 y = zeros (size (z)); 64 y = zeros (size (z));
63 return; 65 return;
64 endif 66 endif
65 67
66 if (nargin < 2 || isempty (tol)) 68 if (nargin < 2 || isempty (tol))
67 if (isa (z, "single")) 69 tol = 100;
68 tol = 100 * eps("single"); 70 elseif (! isscalar (tol) || tol < 0)
69 else 71 error ("cplxpair: TOL must be a positive scalar number")
70 tol = 100*eps; 72 endif
73
74 nd = ndims (z);
75 if (nargin < 3)
76 ## Find the first singleton dimension.
77 sz = size (z);
78 (dim = find (sz > 1, 1)) || (dim = 1);
79 else
80 dim = floor (dim);
81 if (dim < 1 || dim > nd)
82 error ("cplxpair: invalid dimension DIM");
71 endif 83 endif
72 endif 84 endif
73 85
74 nd = ndims (z); 86 ## Move dimension to treat to first position, and convert to a 2-D matrix.
75 orig_dims = size (z);
76 if (nargin < 3)
77 ## Find the first singleton dimension.
78 dim = 0;
79 while (dim < nd && orig_dims(dim+1) == 1)
80 dim++;
81 endwhile
82 dim++;
83 if (dim > nd)
84 dim = 1;
85 endif
86 else
87 dim = floor (dim);
88 if (dim < 1 || dim > nd)
89 error ("cplxpair: invalid dimension along which to sort");
90 endif
91 endif
92
93 ## Move dimension to treat first, and convert to a 2-D matrix.
94 perm = [dim:nd, 1:dim-1]; 87 perm = [dim:nd, 1:dim-1];
95 z = permute (z, perm); 88 z = permute (z, perm);
96 sz = size (z); 89 sz = size (z);
97 n = sz(1); 90 n = sz(1);
98 m = prod (sz) / n; 91 m = prod (sz) / n;
101 ## Sort the sequence in terms of increasing real values. 94 ## Sort the sequence in terms of increasing real values.
102 [q, idx] = sort (real (z), 1); 95 [q, idx] = sort (real (z), 1);
103 z = z(idx + n * ones (n, 1) * [0:m-1]); 96 z = z(idx + n * ones (n, 1) * [0:m-1]);
104 97
105 ## Put the purely real values at the end of the returned list. 98 ## Put the purely real values at the end of the returned list.
106 cls = "double"; 99 cls = ifelse (isa (z, "single"), "single", "double");
107 if (isa (z, "single")) 100 [idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin (cls)) ...
108 cls = "single"; 101 < tol*eps (abs (z)));
109 endif
110 [idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin (cls)) < tol);
111 q = sparse (idxi, idxj, 1, n, m); 102 q = sparse (idxi, idxj, 1, n, m);
112 nr = sum (q, 1); 103 nr = sum (q, 1);
113 [q, idx] = sort (q, 1); 104 [q, idx] = sort (q, 1);
114 z = z(idx); 105 z = z(idx);
115 y = z; 106 y = z;
116 107
117 ## For each remaining z, place the value and its conjugate at the 108 ## For each remaining z, place the value and its conjugate at the start of
118 ## start of the returned list, and remove them from further 109 ## the returned list, and remove them from further consideration.
119 ## consideration.
120 for j = 1:m 110 for j = 1:m
121 p = n - nr(j); 111 p = n - nr(j);
122 for i = 1:2:p 112 for i = 1:2:p
123 if (i+1 > p) 113 if (i+1 > p)
124 error ("cplxpair: could not pair all complex numbers"); 114 error ("cplxpair: could not pair all complex numbers");
125 endif 115 endif
126 [v, idx] = min (abs (z(i+1:p) - conj (z(i)))); 116 [v, idx] = min (abs (z(i+1:p) - conj (z(i))));
127 if (v > tol) 117 if (v > tol*eps (abs (z(i))))
128 error ("cplxpair: could not pair all complex numbers"); 118 error ("cplxpair: could not pair all complex numbers");
129 endif 119 endif
130 if (imag (z(i)) < 0) 120 if (imag (z(i)) < 0)
131 y([i, i+1]) = z([i, idx+i]); 121 y([i, i+1]) = z([i, idx+i]);
132 else 122 else
157 %!shared z 147 %!shared z
158 %! z = exp (2i*pi*[4; 3; 5; 2; 6; 1; 0]/7); 148 %! z = exp (2i*pi*[4; 3; 5; 2; 6; 1; 0]/7);
159 %!assert (cplxpair (z(randperm (7))), z) 149 %!assert (cplxpair (z(randperm (7))), z)
160 %!assert (cplxpair (z(randperm (7))), z) 150 %!assert (cplxpair (z(randperm (7))), z)
161 %!assert (cplxpair (z(randperm (7))), z) 151 %!assert (cplxpair (z(randperm (7))), z)
162 %!assert (cplxpair ([z(randperm(7)),z(randperm(7))]), [z,z]) 152 %!assert (cplxpair ([z(randperm (7)), z(randperm (7))]), [z,z])
163 %!assert (cplxpair ([z(randperm(7)),z(randperm(7))],[],1), [z,z]) 153 %!assert (cplxpair ([z(randperm (7)), z(randperm (7))],[],1), [z,z])
164 %!assert (cplxpair ([z(randperm(7)).';z(randperm(7)).'],[],2), [z.';z.']) 154 %!assert (cplxpair ([z(randperm (7)).'; z(randperm (7)).'],[],2), [z.';z.'])
165 155
166 ## tolerance test 156 ## Test tolerance
167 %!assert (cplxpair ([1i, -1i, 1+(1i*eps)],2*eps), [-1i, 1i, 1+(1i*eps)]) 157 %!assert (cplxpair ([2000 * (1+eps) + 4j; 2000 * (1-eps) - 4j]), ...
158 %! [(2000 - 4j); (2000 + 4j)], 100*eps(200))
159 %!error <could not pair> cplxpair ([2000 * (1+eps) + 4j; 2000 * (1-eps) - 4j], 0);
168 160
161 %!error <could not pair> cplxpair ([2e6 + j; 2e6 - j; 1e-9 * (1 + j); 1e-9 * (1 - 2j)]);
162
163 ## Test input validation
164 %!error cplxpair ()
165 %!error cplxpair (1,2,3,4)
166 %!error <TOL must be .* positive> cplxpair (1, -1)
167 %!error <TOL must be .* scalar number> cplxpair (1, ones (2,2))
168 %!error <invalid dimension DIM> cplxpair (1, [], 3)
169