Scheduled service maintenance on November 22


On Friday, November 22, 2024, between 06:00 CET and 18:00 CET, GIN services will undergo planned maintenance. Extended service interruptions should be expected. We will try to keep downtimes to a minimum, but recommend that users avoid critical tasks, large data uploads, or DOI requests during this time.

We apologize for any inconvenience.

npermutek.m 3.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131
  1. function [Matrix,Index] = npermutek(N,K)
  2. %NPERMUTEK Permutation of elements with replacement/repetition.
  3. % MAT = NPERMUTEK(N,K) returns all possible samplings of length K from
  4. % vector N of type: ordered sample with replacement.
  5. % MAT has size (length(N)^K)-by-K, where K must be a scalar.
  6. % [MAT, IDX] = NPERMUTEK(N,K) also returns IDX such that MAT = N(IDX).
  7. % N may be of class: single, double, or char. If N is single or double,
  8. % both MAT and IDX will be of the same class.
  9. %
  10. % For N = 1:M, for some integer M>1, all(MAT(:)==IDX(:)), so there is no
  11. % benefit to calling NPERMUTEK with two output arguments.
  12. %
  13. % Examples:
  14. % MAT = npermutek([2 4 5],2)
  15. %
  16. % MAT =
  17. %
  18. % 2 2
  19. % 2 4
  20. % 2 5
  21. % 4 2
  22. % 4 4
  23. % 4 5
  24. % 5 2
  25. % 5 4
  26. % 5 5
  27. %
  28. % NPERMUTEK also works on characters.
  29. %
  30. % MAT = npermutek(['a' 'b' 'c'],2)
  31. % MAT =
  32. %
  33. % aa
  34. % ab
  35. % ac
  36. % ba
  37. % bb
  38. % bc
  39. % ca
  40. % cb
  41. % cc
  42. %
  43. % See also perms, nchoosek
  44. %
  45. % Also on the web:
  46. % http://mathworld.wolfram.com/BallPicking.html
  47. % See the section on Enumerative combinatorics below:
  48. % http://en.wikipedia.org/wiki/Permutations_and_combinations
  49. % Author: Matt Fig
  50. % Contact: popkenai@yahoo.com
  51. if nargin ~= 2
  52. error('NPERMUTEK requires two arguments. See help.')
  53. end
  54. if isempty(N) || K == 0,
  55. Matrix = [];
  56. Index = Matrix;
  57. return
  58. elseif floor(K) ~= K || K<0 || ~isreal(K) || numel(K)~=1
  59. error('Second argument should be a real positive integer. See help.')
  60. end
  61. LN = numel(N); % Used in calculating the Matrix and Index.
  62. if K==1
  63. Matrix = N(:); % This one is easy to calculate.
  64. Index = (1:LN).';
  65. return
  66. elseif LN==1
  67. Index = ones(K,1);
  68. Matrix = N(1,Index);
  69. return
  70. end
  71. CLS = class(N);
  72. if ischar(N)
  73. CLS = 'double'; % We will deal with this at the end.
  74. flg = 1;
  75. N = double(N);
  76. else
  77. flg = 0;
  78. end
  79. L = LN^K; % This is the number of rows the outputs will have.
  80. Matrix = zeros(L,K,CLS); % Preallocation.
  81. D = diff(N(1:LN)); % Use this for cumsumming later.
  82. LD = length(D); % See comment on LN.
  83. VL = [-sum(D) D].'; % These values will be put into Matrix.
  84. % Now start building the matrix.
  85. TMP = VL(:,ones(L/LN,1,CLS)); % Instead of repmatting.
  86. Matrix(:,K) = TMP(:); % We don't need to do two these in loop.
  87. Matrix(1:LN^(K-1):L,1) = VL; % The first column is the simplest.
  88. if nargout==1
  89. % Here we only have to build Matrix the rest of the way.
  90. for ii = 2:K-1
  91. ROWS = 1:LN^(ii-1):L; % Indices into the rows for this col.
  92. TMP = VL(:,ones(length(ROWS)/(LD+1),1,CLS)); % Match dimension.
  93. Matrix(ROWS,K-ii+1) = TMP(:); % Build it up, insert values.
  94. end
  95. else
  96. % Here we have to finish Matrix and build Index.
  97. Index = zeros(L,K,CLS); % Preallocation.
  98. VL2 = ones(size(VL),CLS); % Follow the logic in VL above.
  99. VL2(1) = 1-LN; % These are the drops for cumsum.
  100. TMP2 = VL2(:,ones(L/LN,1,CLS)); % Instead of repmatting.
  101. Index(:,K) = TMP2(:); % We don't need to do two these in loop.
  102. Index(1:LN^(K-1):L,1) = 1;
  103. for ii = 2:K-1
  104. ROWS = 1:LN^(ii-1):L; % Indices into the rows for this col.
  105. F = ones(length(ROWS)/(LD+1),1,CLS); % Don't do it twice!
  106. TMP = VL(:,F); % Match dimensions.
  107. TMP2 = VL2(:,F);
  108. Matrix(ROWS,K-ii+1) = TMP(:); % Build them up, insert values.
  109. Index(ROWS,K-ii+1) = TMP2(:);
  110. end
  111. Index(1,:) = 1; % The first row must be 1 for proper cumsumming.
  112. Index = cumsum(Index); % This is the time hog.
  113. end
  114. Matrix(1,:) = N(1); % For proper cumsumming.
  115. Matrix = cumsum(Matrix); % This is the time hog.
  116. if flg
  117. Matrix = char(Matrix); % char was implicitly cast to double above.
  118. end