1function v = allVL1(n, L1, L1ops, MaxNbSol) 2% All integer permutations with sum criteria 3% 4% function v=allVL1(n, L1); OR 5% v=allVL1(n, L1, L1opt); 6% v=allVL1(n, L1, L1opt, MaxNbSol); 7% 8% INPUT 9% n: length of the vector 10% L1: target L1 norm 11% L1ops: optional string ('==' or '<=' or '<') 12% default value is '==' 13% MaxNbSol: integer, returns at most MaxNbSol permutations. 14% When MaxNbSol is NaN, allVL1 returns the total number of all possible 15% permutations, which is useful to check the feasibility before getting 16% the permutations. 17% OUTPUT: 18% v: (m x n) array such as: sum(v,2) == L1, 19% (or <= or < depending on L1ops) 20% all elements of v is naturel numbers {0,1,...} 21% v contains all (=m) possible combinations 22% v is sorted by sum (L1 norm), then by dictionnary sorting criteria 23% class(v) is same as class(L1) 24% Algorithm: 25% Recursive 26% Remark: 27% allVL1(n,L1-n)+1 for natural numbers defined as {1,2,...} 28% Example: 29% This function can be used to generate all orders of all 30% multivariable polynomials of degree p in R^n: 31% Order = allVL1(n, p) 32% Author: Bruno Luong (brunoluong@yahoo.com) 33% Original, 30/nov/2007 34% Version 1.1, 30/apr/2008: Add H1 line as suggested by John D'Errico 35% 1.2, 17/may/2009: Possibility to get the number of permutations 36% alone (set fourth parameter MaxNbSol to NaN) 37% 1.3, 16/Sep/2009: Correct bug for number of solution 38% 1.4, 18/Dec/2010: + non-recursive engine 39 40% Retrieved from https://www.mathworks.com/matlabcentral/fileexchange/17818-all-permutations-of-integers-with-sum-criteria 41% ========================================================================= 42% Copyright (C) 2007-2010 Bruno Luong <brunoluong@yahoo.com> 43% Copyright (C) 2020 Dynare Team 44% 45% This file is part of Dynare. 46% 47% Dynare is free software: you can redistribute it and/or modify 48% it under the terms of the GNU General Public License as published by 49% the Free Software Foundation, either version 3 of the License, or 50% (at your option) any later version. 51% 52% Dynare is distributed in the hope that it will be useful, 53% but WITHOUT ANY WARRANTY; without even the implied warranty of 54% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 55% GNU General Public License for more details. 56% 57% You should have received a copy of the GNU General Public License 58% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 59% ========================================================================= 60global MaxCounter; 61 62if nargin<3 || isempty(L1ops) 63 L1ops = '=='; 64end 65 66n = floor(n); % make sure n is integer 67 68if n<1 69 v = []; 70 return 71end 72 73if nargin<4 || isempty(MaxNbSol) 74 MaxCounter = Inf; 75else 76 MaxCounter = MaxNbSol; 77end 78Counter(0); 79 80switch L1ops 81 case {'==' '='}, 82 if isnan(MaxCounter) 83 % return the number of solutions 84 v = nchoosek(n+L1-1,L1); % nchoosek(n+L1-1,n-1) 85 else 86 v = allVL1eq(n, L1); 87 end 88 case '<=', % call allVL1eq for various sum targets 89 if isnan(MaxCounter) 90 % return the number of solutions 91 %v = nchoosek(n+L1,L1)*factorial(n-L1); BUG <- 16/Sep/2009: 92 v = 0; 93 for j=0:L1 94 v = v + nchoosek(n+j-1,j); 95 end 96 % See pascal's 11th identity, the sum doesn't seem to 97 % simplify to a fix formula 98 else 99 v = cell2mat(arrayfun(@(j) allVL1eq(n, j), (0:L1)', ... 100 'UniformOutput', false)); 101 end 102 case '<', 103 v = allVL1(n, L1-1, '<=', MaxCounter); 104 otherwise 105 error('allVL1: unknown L1ops') 106end 107 108end % allVL1 109 110%% 111function v = allVL1eq(n, L1) 112 113global MaxCounter; 114 115n = feval(class(L1),n); 116s = n+L1; 117sd = double(n)+double(L1); 118notoverflowed = double(s)==sd; 119if isinf(MaxCounter) && notoverflowed 120 v = allVL1nonrecurs(n, L1); 121else 122 v = allVL1recurs(n, L1); 123end 124 125end % allVL1eq 126 127%% Recursive engine 128function v = allVL1recurs(n, L1, head) 129% function v=allVL1eq(n, L1); 130% INPUT 131% n: length of the vector 132% L1: desired L1 norm 133% head: optional parameter to by concatenate in the first column 134% of the output 135% OUTPUT: 136% if head is not defined 137% v: (m x n) array such as sum(v,2)==L1 138% all elements of v is naturel numbers {0,1,...} 139% v contains all (=m) possible combinations 140% v is (dictionnary) sorted 141% Algorithm: 142% Recursive 143 144global MaxCounter; 145 146if n==1 147 if Counter < MaxCounter 148 v = L1; 149 else 150 v = zeros(0,1,class(L1)); 151 end 152else % recursive call 153 v = cell2mat(arrayfun(@(j) allVL1recurs(n-1, L1-j, j), (0:L1)', ... 154 'UniformOutput', false)); 155end 156 157if nargin>=3 % add a head column 158 v = [head+zeros(size(v,1),1,class(head)) v]; 159end 160 161end % allVL1recurs 162 163%% 164function res=Counter(newval) 165 persistent counter; 166 if nargin>=1 167 counter = newval; 168 res = counter; 169 else 170 res = counter; 171 counter = counter+1; 172 end 173end % Counter 174 175%% Non-recursive engine 176function v = allVL1nonrecurs(n, L1) 177% function v=allVL1eq(n, L1); 178% INPUT 179% n: length of the vector 180% L1: desired L1 norm 181% OUTPUT: 182% if head is not defined 183% v: (m x n) array such as sum(v,2)==L1 184% all elements of v is naturel numbers {0,1,...} 185% v contains all (=m) possible combinations 186% v is (dictionnary) sorted 187% Algorithm: 188% NonRecursive 189 190% Chose (n-1) the splitting points of the array [0:(n+L1)] 191s = nchoosek(1:n+L1-1,n-1); 192m = size(s,1); 193 194s1 = zeros(m,1,class(L1)); 195s2 = (n+L1)+s1; 196 197v = diff([s1 s s2],1,2); % m x n 198v = v-1; 199 200end % allVL1nonrecurs 201