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