1########################################################################
2##
3## Copyright (C) 2015-2021 The Octave Project Developers
4##
5## See the file COPYRIGHT.md in the top-level directory of this
6## distribution or <https://octave.org/copyright/>.
7##
8## This file is part of Octave.
9##
10## Octave is free software: you can redistribute it and/or modify it
11## under the terms of the GNU General Public License as published by
12## the Free Software Foundation, either version 3 of the License, or
13## (at your option) any later version.
14##
15## Octave is distributed in the hope that it will be useful, but
16## WITHOUT ANY WARRANTY; without even the implied warranty of
17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18## GNU General Public License for more details.
19##
20## You should have received a copy of the GNU General Public License
21## along with Octave; see the file COPYING.  If not, see
22## <https://www.gnu.org/licenses/>.
23##
24########################################################################
25
26## -*- texinfo -*-
27## @deftypefn  {} {@var{xxx} =} repelem (@var{x}, @var{R})
28## @deftypefnx {} {@var{xxx} =} repelem (@var{x}, @var{R_1}, @dots{}, @var{R_n})
29## Construct an array of repeated elements from @var{x} and repeat
30## instructions @var{R_1}, @dots{}.
31##
32## @var{x} must be a scalar, vector, or N-dimensional array.
33##
34## A repeat instruction @var{R_j} must either be a scalar or a vector.  If the
35## instruction is a scalar then each component of @var{x} in dimension @var{j}
36## is repeated @var{R_j} times.  If the instruction is a vector then it must
37## have the same number of elements as the corresponding dimension @var{j} of
38## @var{x}.  In this case, the @var{k}th component of dimension @var{j} is
39## repeated @code{@var{R_j}(@var{k})} times.
40##
41## If @var{x} is a scalar or vector then @code{repelem} may be called with just
42## a single repeat instruction @var{R} and @code{repelem} will return a vector
43## with the same orientation as the input.
44##
45## If @var{x} is a matrix then at least two @var{R_j}s must be specified.
46##
47## Note: Using @code{repelem} with a vector @var{x} and a vector for @var{R_j}
48## is equivalent to Run Length Decoding.
49##
50## Examples:
51##
52## @example
53## @group
54## A = [1 2 3 4 5];
55## B = [2 1 0 1 2];
56## repelem (A, B)
57##   @result{}   1   1   2   4   5   5
58## @end group
59## @end example
60##
61## @example
62## @group
63## A = magic (3)
64##   @result{} A =
65##        8   1   6
66##        3   5   7
67##        4   9   2
68## B1 = [1 2 3];
69## B2 = 2;
70## repelem (A, B1, B2)
71##   @result{}     8   8   1   1   6   6
72##          3   3   5   5   7   7
73##          3   3   5   5   7   7
74##          4   4   9   9   2   2
75##          4   4   9   9   2   2
76##          4   4   9   9   2   2
77## @end group
78## @end example
79##
80## More @var{R_j} may be specified than the number of dimensions of @var{x}.
81## Any excess @var{R_j} must be scalars (because @var{x}'s size in those
82## dimensions is only 1), and @var{x} will be replicated in those dimensions
83## accordingly.
84##
85## @example
86## @group
87## A = [1 2 3 4 5];
88## B1 = 2;
89## B2 = [2 1 3 0 2];
90## B3 = 3;
91## repelem (A, B1, B2, B3)
92##   @result{}    ans(:,:,1) =
93##            1   1   2   3   3   3   5   5
94##            1   1   2   3   3   3   5   5
95##
96##         ans(:,:,2) =
97##
98##            1   1   2   3   3   3   5   5
99##            1   1   2   3   3   3   5   5
100##
101##         ans(:,:,3) =
102##            1   1   2   3   3   3   5   5
103##            1   1   2   3   3   3   5   5
104## @end group
105## @end example
106##
107## @var{R_j} must be specified in order.  A placeholder of 1 may be used for
108## dimensions which do not need replication.
109##
110## @example
111## @group
112## repelem ([-1, 0; 0, 1], 1, 2, 1, 2)
113##   @result{}  ans(:,:,1,1) =
114##         -1  -1   0   0
115##          0   0   1   1
116##
117##       ans(:,:,1,2) =
118##         -1  -1   0   0
119##          0   0   1   1
120## @end group
121## @end example
122##
123## If fewer @var{R_j} are given than the number of dimensions in @var{x},
124## @code{repelem} will assume @var{R_j} is 1 for those dimensions.
125##
126## @example
127## A = cat (3, [-1 0; 0 1], [-1 0; 0 1])
128##   @result{}  ans(:,:,1) =
129##         -1   0
130##          0   1
131##
132##       ans(:,:,2) =
133##         -1   0
134##          0   1
135##
136## repelem (A,2,3)
137##   @result{}  ans(:,:,1) =
138##         -1  -1  -1   0   0   0
139##         -1  -1  -1   0   0   0
140##          0   0   0   1   1   1
141##          0   0   0   1   1   1
142##
143##       ans(:,:,2) =
144##         -1  -1  -1   0   0   0
145##         -1  -1  -1   0   0   0
146##          0   0   0   1   1   1
147##          0   0   0   1   1   1
148## @end example
149##
150## @code{repelem} preserves the class of @var{x}, and works with strings,
151## cell arrays, NA, and NAN inputs.  If any @var{R_j} is 0 the output will
152## be an empty array.
153##
154## @example
155## @group
156## repelem ("Octave", 2, 3)
157##   @result{}    OOOccctttaaavvveee
158##         OOOccctttaaavvveee
159##
160## repelem ([1 2 3; 1 2 3], 2, 0)
161##   @result{}    [](4x0)
162## @end group
163## @end example
164##
165## @seealso{cat, kron, repmat}
166## @end deftypefn
167
168## Author: Markus Bergholz <markuman@gmail.com>
169## Author: Nicholas R. Jankowski <jankowskin@asme.org>
170
171## As a U.S. government employee, Nicholas R. Jankowski makes no claim
172## of copyright.
173
174## The prepareIdx routine is Copyright (C) 2015 Peter John Acklam
175## <pjacklam@gmail.com>, used with permission.
176
177function retval = repelem (x, varargin)
178
179  if (nargin <= 1)
180    print_usage ();
181
182  elseif (nargin == 2)
183
184    R = varargin{1};
185
186    if (isscalar (R))
187
188      if (! isvector (x))
189        error (["repelem: %dD Array requires %d or more input " ...
190                "arguments, but only %d given"], ...
191               ndims (x), ndims (x) + 1, nargin);
192      endif
193
194      if (iscolumn (x))
195        ## element values repeated R times in a col vector
196        retval = x.'(ones (R, 1), :)(:);
197      else
198        ## element values repeated R times in a row vector
199        retval = x(ones (R, 1), :)(:).';
200      endif
201
202    elseif (isvector (x) && isvector (R))
203
204      ## vector x with vector repeat.
205      if (numel (R) != numel (x))
206        error (["repelem: R1 must either be scalar or have the same " ...
207                "number of elements as the vector to be replicated"]);
208      endif
209
210      ## Basic run-length decoding in function prepareIdx returns
211      ## idx2 as a row vector of element indices in the right positions.
212      idx2 = prepareIdx (R);
213      ## Fill with element values, direction matches element.
214      retval = x(idx2);
215
216    else # catch any arrays passed to x or varargin with nargin==2
217      error (["repelem: when called with only two inputs they must be " ...
218              "either scalars or vectors, not %s and %s."],
219             typeinfo (x), typeinfo (R));
220    endif
221
222  elseif (nargin == 3)  # special optimized case for 2-D (matrices)
223
224    ## Input Validation
225    xsz = size (x);
226    vector_r = ! (cellfun (@numel, varargin) == 1);
227
228    ## 1. Check that all varargin are either scalars or vectors, not arrays.
229    ##    isvector returns true for scalars so one test captures both inputs.
230    if (! (isvector (varargin{1}) && (isvector (varargin{2}))))
231      error ("repelem: R1 and R2 must be scalars or vectors");
232
233    ## 2. check that any repeat vectors have the right length.
234    elseif (any (cellfun (@numel, varargin(vector_r)) != xsz(vector_r)))
235      error (["repelem: R_j vectors must have the same number of elements " ...
236              "as the size of dimension j of X"]);
237    endif
238
239    ## Create index arrays to pass to element.
240    ## (It is no slower to call prepareIdx than to check and do scalars
241    ## directly.)
242    idx1 = prepareIdx (varargin{1}, xsz(1));
243    idx2 = prepareIdx (varargin{2}, xsz(2));
244
245    if (issparse (x))
246      retval = x(idx1, idx2);
247    else
248      ## The ":" at the end takes care of any x dimensions > 2.
249      retval = x(idx1, idx2, :);
250    endif
251
252  else  # (nargin > 3)
253
254    ## Input Validation
255    xsz = size (x);
256    n_xdims = numel (xsz);
257    vector_r = ! (cellfun (@numel, varargin) == 1);
258
259    ## 1. Check that all repeats are scalars or vectors
260    ##    (isvector gives true for scalars);
261    if (! all (cellfun (@isvector, varargin(vector_r))))
262      error ("repelem: R_j must all be scalars or vectors");
263
264    ## 2. Catch any vectors thrown at trailing singletons,
265    ##    which should only have scalars;
266    elseif (find (vector_r, 1, "last") > n_xdims)
267      error ("repelem: R_j for trailing singleton dimensions must be scalar");
268
269    ## 3. Check that the ones that are vectors have the right length.
270    elseif (any (cellfun (@numel, varargin(vector_r)) != xsz(vector_r)))
271      error (["repelem: R_j vectors must have the same number of elements " ...
272              "as the size of dimension j of X"]);
273
274    endif
275
276    n_rpts = nargin - 1;
277    dims_with_vectors_and_scalars = min (n_xdims, n_rpts);
278
279    ## Preallocate idx which will contain index array to be put into element.
280    idx = cell (1, n_rpts);
281
282    ## Use prepareIdx() to fill indices for dimensions that could be
283    ## a scalar or a vector.
284    for i = 1 : dims_with_vectors_and_scalars
285      idx(i) = prepareIdx (varargin{i}, xsz(i));
286    endfor
287
288    ## If there are more varargin inputs than x dimensions, then input tests
289    ## have verified that they are just scalars, so add [1 1 1 1 1 ... 1] to
290    ## those dims to perform concatenation along those dims.
291    if (n_rpts > n_xdims)
292      for i = n_xdims + (1 : (n_rpts - n_xdims))
293        idx(i) = ones (1, varargin{i});
294      endfor
295    endif
296
297    ## Use completed idx to specify repetition of x values in all dimensions.
298    ## The trailing ":" will take care of cases where n_xdims > n_rpts.
299    retval = x(idx{:}, :);
300
301  endif
302
303endfunction
304
305## Return a row vector of indices prepared for replicating.
306function idx = prepareIdx (v, n)
307
308  if (isscalar (v))
309    ## will always return row vector
310    idx = [1:n](ones (v, 1), :)(:).';
311
312  else
313    ## This works for a row or column vector.
314
315    ## Get ending position for each element item.
316    idx_temp = cumsum (v);
317
318    ## Set starting position of each element to 1.
319    idx(idx_temp + 1) = 1;
320
321    ## Set starting position of each element to 1.
322    idx(1) = 1;
323
324    ## Row vector with proper length for output
325    idx = idx(1:idx_temp(end));
326
327    ## with prepared index
328    idx = (find (v != 0))(cumsum (idx));
329
330  endif
331
332endfunction
333
334
335## tests for help examples
336%!assert (repelem ([1 2 3 4 5], [2 1 0 1 2]), [1 1 2 4 5 5])
337%!assert (repelem (magic(3), [1 2 3],2), ...
338%!  [8 8 1 1 6 6;3 3 5 5 7 7;3 3 5 5 7 7;4 4 9 9 2 2;4 4 9 9 2 2;4 4 9 9 2 2])
339%!assert (repelem ([1 2 3 4 5],2,[2 1 3 0 2],3),repmat([1 1 2 3 3 3 5 5],2,1,3))
340%!assert (repelem ([-1 0;0 1],1,2,1,2), repmat([-1 -1 0 0; 0 0 1 1],1,1,1,2))
341%!assert (repelem (cat(3,[-1 0 ; 0 1],[-1 0 ; 0 1]),2,3), ...
342%!  repmat([-1 -1 -1 0 0 0;-1 -1 -1 0 0 0;0 0 0 1 1 1;0 0 0 1 1 1],1,1,2))
343%!assert (repelem ("Octave", 2,3), ["OOOccctttaaavvveee";"OOOccctttaaavvveee"])
344
345## test complex vectors are not Hermitian conjugated
346%!assert (repelem ([i, -i], 2), [i, i, -i, -i])
347%!assert (repelem ([i; -i], 2), [i; i; -i; -i])
348
349## nargin == 2 tests
350%!assert (repelem ([-1 0 1], 2), [-1 -1 0 0 1 1])
351%!assert (repelem ([-1 0 1]', 2), [-1; -1; 0; 0; 1; 1;])
352%!assert (repelem ([-1 0 1], [1 2 1]), [-1 0 0 1])
353%!assert (repelem ([-1 0 1]', [1 2 1]), [-1; 0; 0; 1])
354%!assert (repelem ([1 2 3 4 5]', [2 1 0 1 2]), [1 1 2 4 5 5]')
355
356## nargin == 3 tests
357%!assert (repelem ([1 0;0 -1], 2, 3),
358%!       [1 1 1 0 0 0;1 1 1 0 0 0;0 0 0 -1 -1 -1;0 0 0 -1 -1 -1])
359%!assert (repelem ([1 0; 0 -1], 1,[3 2]), [1 1 1 0 0;0 0 0 -1 -1])
360%!assert (repelem ([1 0; 0 -1], 2,[3 2]),
361%!        [1 1 1 0 0;1 1 1 0 0;0 0 0 -1 -1;0 0 0 -1 -1])
362%!assert (repelem (cat(3,[1 0; 0 -1],[1 0;0 -1]), 1,[3 2]),
363%!        repmat([1 1 1 0 0 ; 0 0 0 -1 -1],1,1,2))
364%!assert (repelem ([1 0; 0 -1], [3 2], 1), [1 0;1 0;1 0;0 -1;0 -1])
365%!assert (repelem ([1 0; 0 -1], [3 2], 2),
366%!        [1 1 0 0;1 1 0 0;1 1 0 0;0 0 -1 -1;0 0 -1 -1])
367%!assert (repelem ([1 0; 0 -1], [2 3] ,[3 2]),
368%!        [1 1 1 0 0;1 1 1 0 0;0 0 0 -1 -1;0 0 0 -1 -1;0 0 0 -1 -1])
369%!assert (repelem (cat(3,[1 1 1 0;0 1 0 0],[1 1 1 1;0 0 0 1],[1 0 0 1;1 1 0 1]),
370%!                2, 3),
371%!        cat(3,[1 1 1 1 1 1 1 1 1 0 0 0
372%!               1 1 1 1 1 1 1 1 1 0 0 0
373%!               0 0 0 1 1 1 0 0 0 0 0 0
374%!               0 0 0 1 1 1 0 0 0 0 0 0],
375%!              [1 1 1 1 1 1 1 1 1 1 1 1
376%!               1 1 1 1 1 1 1 1 1 1 1 1
377%!               0 0 0 0 0 0 0 0 0 1 1 1
378%!               0 0 0 0 0 0 0 0 0 1 1 1],
379%!              [1 1 1 0 0 0 0 0 0 1 1 1
380%!               1 1 1 0 0 0 0 0 0 1 1 1
381%!               1 1 1 1 1 1 0 0 0 1 1 1
382%!               1 1 1 1 1 1 0 0 0 1 1 1]))
383%!assert (repelem (cat(3,[1 1 1 0;0 1 0 0],[1 1 1 1;0 0 0 1],[1 0 0 1;1 1 0 1]),
384%!                2, [3 3 3 3]), ...
385%!        cat(3,[1 1 1 1 1 1 1 1 1 0 0 0
386%!               1 1 1 1 1 1 1 1 1 0 0 0
387%!               0 0 0 1 1 1 0 0 0 0 0 0
388%!               0 0 0 1 1 1 0 0 0 0 0 0], ...
389%!              [1 1 1 1 1 1 1 1 1 1 1 1
390%!               1 1 1 1 1 1 1 1 1 1 1 1
391%!               0 0 0 0 0 0 0 0 0 1 1 1
392%!               0 0 0 0 0 0 0 0 0 1 1 1], ...
393%!              [1 1 1 0 0 0 0 0 0 1 1 1
394%!               1 1 1 0 0 0 0 0 0 1 1 1
395%!               1 1 1 1 1 1 0 0 0 1 1 1
396%!               1 1 1 1 1 1 0 0 0 1 1 1]));
397%!assert (repelem ([1 2 3 4 5], 2,[2 1 2 0 2]), [1 1 2 3 3 5 5;1 1 2 3 3 5 5])
398%
399## nargin > 3 tests
400%!assert (repelem ([1 0;0 -1], 2, 3, 4), ...
401%!        cat(3,[1 1 1 0 0 0;1 1 1 0 0 0;0 0 0 -1 -1 -1;0 0 0 -1 -1 -1], ...
402%!              [1 1 1 0 0 0;1 1 1 0 0 0;0 0 0 -1 -1 -1;0 0 0 -1 -1 -1], ...
403%!              [1 1 1 0 0 0;1 1 1 0 0 0;0 0 0 -1 -1 -1;0 0 0 -1 -1 -1], ...
404%!              [1 1 1 0 0 0;1 1 1 0 0 0;0 0 0 -1 -1 -1;0 0 0 -1 -1 -1]))
405%!assert (repelem (repmat([-1 0;0 1],1,1,2,3),2,2,2), ...
406%!        repmat([-1 -1 0 0;-1 -1 0 0;0 0 1 1; 0 0 1 1],1,1,4,3))
407%!assert (repelem (repmat([-1 0;0 1],1,1,2,3),[2 2],[2 2],2), ...
408%!        repmat([-1 -1 0 0;-1 -1 0 0;0 0 1 1; 0 0 1 1],1,1,4,3))
409%!assert (repelem (repmat([-1 0;0 1],1,1,2,3),2,2,2,2,2), ...
410%!        repmat([-1 -1 0 0;-1 -1 0 0;0 0 1 1; 0 0 1 1],1,1,4,6,2))
411%!assert (repelem ([1,0,-1;-1,0,1],[2 3],[2 3 4],2), ...
412%!        cat(3,[ 1  1 0 0 0 -1 -1 -1 -1
413%!                1  1 0 0 0 -1 -1 -1 -1
414%!               -1 -1 0 0 0  1  1  1  1
415%!               -1 -1 0 0 0  1  1  1  1
416%!               -1 -1 0 0 0  1  1  1  1], ...
417%!              [ 1  1 0 0 0 -1 -1 -1 -1
418%!                1  1 0 0 0 -1 -1 -1 -1
419%!               -1 -1 0 0 0  1  1  1  1
420%!               -1 -1 0 0 0  1  1  1  1
421%!               -1 -1 0 0 0  1  1  1  1]));
422%!assert (repelem ([1 2 3;4 5 6],[0 2],2,2), repmat([4 4 5 5 6 6],2,1,2))
423
424## test with structures
425%!test
426%! a(2).x = 1;
427%! a(2).y = 2;
428%! a(1).x = 3;
429%! a(1).y = 4;
430%! b = repelem (a, 2, [1 3]);
431%! assert (size (b) == [2, 4]);
432%! assert ([b.y], [4 4 2 2 2 2 2 2]);
433
434## test with cell arrays
435%!test
436%! assert (repelem ({-1 0 1},  2), {-1 -1 0 0 1 1});
437%! assert (repelem ({-1 0 1}', 2), {-1; -1; 0; 0; 1; 1;});
438%! assert (repelem ({1 0;0 -1}, 2, 3),
439%!         {1 1 1 0 0 0;1 1 1 0 0 0;0 0 0 -1 -1 -1;0 0 0 -1 -1 -1});
440
441%!test <*54275>
442%! assert (repelem (11:13, [1 3 0]), [11 12 12 12]);
443
444%!test <*59705>
445%! xs = sparse (magic (3));
446%! assert (repelem (xs, 1, 2), ...
447%!         sparse ([8,8,1,1,6,6; 3,3,5,5,7,7; 4,4,9,9,2,2]));
448
449## nargin <= 1 error tests
450%!error (repelem ());
451%!error (repelem ([]));
452%!error (repelem (5));
453%!error (repelem (5,[]));
454%!error (repelem ([1 2 3 3 2 1]));
455%!error (repelem ([1 2 3; 3 2 1]));
456
457## nargin == 2 error tests
458%!error (repelem ([1 2 3; 3 2 1],[]));
459%!error (repelem ([1 2 3; 3 2 1],2));
460%!error (repelem ([1 2 3; 3 2 1],2));
461%!error (repelem ([1 2 3; 3 2 1],[1 2 3]));
462%!error (repelem ([1 2 3; 3 2 1],[1 2 3]'));
463%!error (repelem ([1 2 3; 3 2 1],[1 2 2 1]));
464%!error (repelem ([1 2 3; 3 2 1],[1 2 3;4 5 6]));
465%!error (repelem ([1 2 3 4 5],[1 2 3 4 5;1 2 3 4 5]));
466
467## nargin == 3 error tests
468%!error (repelem ([1 2 3; 3 2 1], 1, [1 2;1 2]));
469%!error (repelem ([1 2 3; 3 2 1], 1, [1 2]));
470%!error (repelem ([1 2 3; 3 2 1], 2, []));
471%!error (repelem ([1 2 3; 3 2 1], [1 2 3], [1 2 3]));
472%!error (repelem ([1 2 3; 3 2 1], [1 2 3], [1 2 3 4]));
473%!error (repelem ([1 2 3; 3 2 1], [1 2], [1 2 3 4]));
474
475## nargin > 3 error tests
476%!error (repelem ([1 2 3; 3 2 1], 1, [1 2;1 2],1,2,3));
477%!error (repelem ([1 2 3; 3 2 1], [],1,2,3));
478%!error (repelem ([1 2 3; 3 2 1], [1 2], [1 2 3],1,2,[1 2;1 2]));
479%!error (repelem ([1 2 3; 3 2 1], [1 2 3], [1 2 3],1,2));
480%!error (repelem ([1 2 3; 3 2 1], [1 2], [1 2 3 4],1,2));
481