Mercurial > octave

comparison scripts/specfun/perms.m @ 31100:05729b0e39e4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
perms.m: Much faster technique to generate unique permutations
author Michael Leitner <michael.leitner@frm2.tum.de>
date Mon, 20 Jun 2022 17:27:39 -0400
parents b949f8a631e2
children 46e15523ca06
comparison
equal deleted inserted replaced
31099:7e6e70dd1c0d 31100:05729b0e39e4
32 ## @code{factorial (@var{n}) * @var{n}}, where @var{n} is the length of 32 ## @code{factorial (@var{n}) * @var{n}}, where @var{n} is the length of
33 ## @var{v}. Any repeated elements are included in the output. 33 ## @var{v}. Any repeated elements are included in the output.
34 ## 34 ##
35 ## If the optional argument @qcode{"unique"} is given, then only unique 35 ## If the optional argument @qcode{"unique"} is given, then only unique
36 ## permutations are returned, using less memory and generally taking less time 36 ## permutations are returned, using less memory and generally taking less time
37 ## than calling @code{unique (perms (@var{v}), "rows")}. In this case, the 37 ## than calling @code{unique (perms (@var{v}), "rows")}.
38 ## results are not guaranteed to be in any particular order.
39 ## 38 ##
40 ## Example 1 39 ## Example 1
41 ## 40 ##
42 ## @example 41 ## @example
43 ## @group 42 ## @group
81 80
82 if (nargin < 1) 81 if (nargin < 1)
83 print_usage (); 82 print_usage ();
84 endif 83 endif
85 84
85 unique_v = 0;
86 if (nargin == 2) 86 if (nargin == 2)
87 if (! strcmpi (opt, "unique")) 87 if (! strcmpi (opt, "unique"))
88 error ('perms: option must be the string "unique".'); 88 error ('perms: option must be the string "unique".');
89 else
90 unique_v = 1;
89 endif 91 endif
90 uniqv = unique (v); 92 endif
91 n = numel (uniqv); 93
92 if (n == 1)
93 P = v;
94 return;
95 elseif (n < numel (v))
96 P = unique_perms (v, uniqv);
97 return;
98 endif
99 endif
100
101 ## If we are here, then we want all permutations, not just unique ones.
102 ## Or the user passed the "unique" option but the input had no repetitions.
103 v = v(:).'; # convert to row vector 94 v = v(:).'; # convert to row vector
104 if (isnumeric (v) || ischar (v)) 95 if (isnumeric (v) || ischar (v))
105 ## Order of output is only dependent on the actual values for 96 ## Order of output is only dependent on the actual values for
106 ## character and numeric arrays. 97 ## character and numeric arrays.
107 v = sort (v, "ascend"); 98 v = sort (v, "ascend");
117 case 2 108 case 2
118 P = [v([2 1]);v]; 109 P = [v([2 1]);v];
119 case 3 110 case 3
120 P = v([3 2 1; 3 1 2; 2 3 1; 2 1 3; 1 3 2; 1 2 3]); 111 P = v([3 2 1; 3 1 2; 2 3 1; 2 1 3; 1 3 2; 1 2 3]);
121 endswitch 112 endswitch
113 if (unique_v)
114 P = unique (P, "rows");
115 endif
122 else 116 else
123 v = v(end:-1:1); 117 if (! unique_v)
124 n-= 1; 118 v = v(end:-1:1);
125 119 n-= 1;
126 idx = zeros (factorial (n), n); 120
127 idx(1:6, n-2:n) = [1, 2, 3;1, 3, 2;2, 1, 3;2, 3, 1;3, 1, 2;3, 2, 1]+(n-3); 121 idx = zeros (factorial (n), n);
128 f = 2; # jump-start for efficiency with medium n 122 idx(1:6, n-2:n) = [1, 2, 3;1, 3, 2;2, 1, 3;2, 3, 1;3, 1, 2;3, 2, 1]+(n-3);
129 for j = 3:n-1 123 f = 2; # jump-start for efficiency with medium n
130 b = 1:n; 124 for j = 3:n-1
131 f *= j; 125 b = 1:n;
132 perm = idx(1:f, n-(j-1):n); 126 f *= j;
133 idx(1:(j+1)*f, n-j) = (n-j:n)(ones (f, 1),:)(:); 127 perm = idx(1:f, n-(j-1):n);
134 for i=0:j 128 idx(1:(j+1)*f, n-j) = (n-j:n)(ones (f, 1),:)(:);
135 b(i+n-j) -= 1; 129 for i=0:j
136 idx((1:f)+i*f, n-(j-1):n) = b(perm); 130 b(i+n-j) -= 1;
131 idx((1:f)+i*f, n-(j-1):n) = b(perm);
132 endfor
137 endfor 133 endfor
138 endfor 134
139 135 n += 1;
140 n += 1; 136 f *= n-1;
141 f *= n-1; 137 P = v(1)(ones (factorial (n), n));
142 P = v(1)(ones (factorial (n), n)); 138 P(:,1) = v(ones (f, 1),:)(:);
143 P(:,1) = v(ones (f, 1),:)(:); 139
144 140 for i = 1:n
145 for i = 1:n 141 b = v([1:i-1 i+1:n]);
146 b = v([1:i-1 i+1:n]); 142 P((1:f)+(i-1)*f, 2:end) = b(idx);
147 P((1:f)+(i-1)*f, 2:end) = b(idx); 143 endfor
148 endfor 144 else # user wants unique permutations
145 [v, i, j] = unique(v);
146 h = accumarray (j, 1)';
147 idx = m_perms (h);
148 P = v (sortrows (idx'));
149 endif
149 endif 150 endif
150 151
151 endfunction 152 endfunction
152 153
153 function P = unique_perms (v, uniqv) 154 function out_perms = m_perms (multis)
154 lenvec = 1:numel (v); 155
155 n_uniq = numel (uniqv); 156 l = numel (multis);
156 157 if (l == 1)
157 ## Calculate frequencies. Slightly faster to not use hist() or histc(). 158 out_perms = uint8 (ones (multis, 1));
158 for i = n_uniq:-1:1 159 else
159 freq(i) = nnz (v == uniqv(i)); 160 p1 = m_perms (multis (1:floor (l/2)));
160 endfor 161 p2 = m_perms (multis (floor(l/2)+1:l)) + max (p1 (:, 1));
161 162 l1 = rows (p1);
162 ## Performance is improved by handling the most frequent elements first. 163 l2 = rows (p2);
163 ## This can change the output order though. Currently we do not promise 164 cp1 = columns (p1);
164 ## the results will be in any specific order. 165 cp2 = columns (p2);
165 [freq, order] = sort (freq, "descend"); 166
166 uniqv = uniqv(order); 167 p = nchoosek (1:l1+l2, l1);
167 168 rp = rows(p);
168 ## Call nchoosek for each choice and do a Cartesian set product, 169
169 ## avoiding collisions. 170 ii = false (l1+l2, rp);
170 idx = nchoosek (lenvec, freq(1)); 171 ii(p + (0:rp - 1)' * (l1 + l2)) = true;
171 for i = 2:n_uniq 172 out_perms = zeros (l1 + l2, cp1 * cp2 * rp, "uint8");
172 this = nchoosek (lenvec, freq(i)); 173 out_perms (repmat ( ii, cp1 * cp2, 1)(:)) = repmat (p1, cp2, rp)(:);
173 174 out_perms (repmat (!ii, cp1 * cp2, 1)(:)) = repmat (p2, 1, cp1 * rp)(:);
174 new = zeros (rows (idx) + 1e5, columns (idx) + columns (this)); 175 endif
175 ## Preallocate 1e5 extra rows.
176 pos = 0;
177 for j = 1:rows (this)
178 tmp = idx;
179 for k = this(j,:)
180 tmp = tmp(all (tmp != k, 2), :);
181 endfor
182
183 r = rows (tmp);
184 c = columns (tmp);
185 if (pos + r > rows (new)) # Allocate more memory on the fly
186 new(pos + r + 1e5, 1) = 0;
187 endif
188
189 new(pos+1:pos+r, 1:c) = tmp;
190 new(pos+1:pos+r, c+1:end) = repmat (this(j,:), r, 1);
191 pos += r;
192 endfor
193 idx = new(1:pos, :);
194 endfor
195 clear new tmp this # Free memory before building P
196
197 idx = (idx-1) * rows (idx) + (1:rows (idx))';
198
199 ## Initialize P to be the same size as idx with the same class as v.
200 P = repmat (v, rows(idx), 1);
201 pos = 0;
202 for i = 1:n_uniq
203 P(idx(:, (pos + 1):(pos + freq(i)))) = uniqv(i);
204 pos += freq(i);
205 endfor
206 176
207 endfunction 177 endfunction
208
209 178
210 %!assert (rows (perms (1:6)), factorial (6)) 179 %!assert (rows (perms (1:6)), factorial (6))
211 %!assert (perms (pi), pi) 180 %!assert (perms (pi), pi)
212 %!assert (perms ([pi, e]), [pi, e; e, pi]) 181 %!assert (perms ([pi, e]), [pi, e; e, pi])
213 %!assert (perms ([1,2,3]), [3,2,1;3,1,2;2,3,1;2,1,3;1,3,2;1,2,3]) 182 %!assert (perms ([1,2,3]), [3,2,1;3,1,2;2,3,1;2,1,3;1,3,2;1,2,3])