1 | function [enums,Matrices] = mpt_enumerate_binary(Matrices) |
---|
2 | |
---|
3 | binary_var_index = Matrices.binary_var_index; |
---|
4 | notbinary_var_index = setdiff(1:Matrices.nu,binary_var_index); |
---|
5 | |
---|
6 | % Detect and extract pure binary equalities. These are used |
---|
7 | % to detect SOS constraints, and for pruning |
---|
8 | nbin = length(binary_var_index); |
---|
9 | only_binary = ~any(Matrices.Aeq(:,notbinary_var_index),2); |
---|
10 | Aeq_bin = Matrices.Aeq(find(only_binary),binary_var_index); |
---|
11 | %Beq_bin = Beq(find(only_binary),binary_var_index); |
---|
12 | beq_bin = Matrices.beq(find(only_binary),:); |
---|
13 | |
---|
14 | % Keep the rest |
---|
15 | Matrices.Aeq = Matrices.Aeq(find(~only_binary),:); |
---|
16 | Matrices.Beq = Matrices.Beq(find(~only_binary),:); |
---|
17 | Matrices.beq = Matrices.beq(find(~only_binary),:); |
---|
18 | |
---|
19 | % Dummy pre-solve to avoid problems |
---|
20 | if length(beq_bin)>0 |
---|
21 | [ii,jj,kk] = unique([Aeq_bin beq_bin],'rows'); |
---|
22 | Aeq_bin = Aeq_bin(jj,:); |
---|
23 | beq_bin = beq_bin(jj,:); |
---|
24 | end |
---|
25 | |
---|
26 | % Normalize... |
---|
27 | neg_b = find(beq_bin < 0); |
---|
28 | Aeq_bin(neg_b,:) = -Aeq_bin(neg_b,:); |
---|
29 | beq_bin(neg_b,:) = -beq_bin(neg_b,:); |
---|
30 | |
---|
31 | % Detect and extract simple binary LP constraints |
---|
32 | only_binary_lp = find(~any([Matrices.G(:,notbinary_var_index) Matrices.E],2)); |
---|
33 | Abin = Matrices.G(only_binary_lp,binary_var_index); |
---|
34 | bbin = Matrices.W(only_binary_lp); |
---|
35 | |
---|
36 | % Remove pure binary constraints |
---|
37 | Matrices.G(only_binary_lp,:) = []; |
---|
38 | Matrices.W(only_binary_lp,:) = []; |
---|
39 | Matrices.E(only_binary_lp,:) = []; |
---|
40 | |
---|
41 | % if ~isempty(Matrices.binary_var_index) |
---|
42 | % fixed_up = find(Matrices.lb(Matrices.binary_var_index) == 1); |
---|
43 | % fixed_down = find(Matrices.ub(Matrices.binary_var_index) == 0); |
---|
44 | % fixed = [fixed_up fixed_down]; |
---|
45 | % if ~isempty(fixed) |
---|
46 | % % Ok, some of the binary are fixed. |
---|
47 | % % Do not enumerate over these |
---|
48 | % % Fix them and correct the constraints that we use |
---|
49 | % % for enumeration and pruning |
---|
50 | % bbin = bbin - Abin(:,fixed)*Matrices.lb(Matrices.binary_var_index(fixed)); |
---|
51 | % Abin(:,fixed) = []; |
---|
52 | % beq_bin = beq_bin - Aeq_bin(:,fixed)*Matrices.lb(Matrices.binary_var_index(fixed)); |
---|
53 | % Aeq_bin(:,fixed) = []; |
---|
54 | % Matrices.binary_var_index(fixed) = []; |
---|
55 | % binary_var_index = Matrices.binary_var_index; |
---|
56 | % notbinary_var_index = setdiff(1:Matrices.nu,binary_var_index); |
---|
57 | % end |
---|
58 | % end |
---|
59 | |
---|
60 | % Detect groups with constraints sum(d_i) == d_j |
---|
61 | SOS = []; |
---|
62 | variables_in_sos = []; |
---|
63 | for i = 1:size(Aeq_bin,1) |
---|
64 | if beq_bin(i) == 0 |
---|
65 | [ix,jx,sx] = find(Aeq_bin(i,:)); |
---|
66 | if all(abs(sx) == 1) |
---|
67 | j = find(sx == -1); |
---|
68 | if length(j) == 1 & length(jx)>2 |
---|
69 | % Aha, we have sum binary(i) == binary(j) |
---|
70 | j = jx(j); |
---|
71 | jx = setdiff(jx,j); |
---|
72 | this_sos = sparse(1:length(jx),jx,1,length(jx),length(binary_var_index)); |
---|
73 | this_sos(1:end,j)= 1; |
---|
74 | this_sos(end+1,1)=0; |
---|
75 | if isempty(SOS) |
---|
76 | SOS = this_sos; |
---|
77 | else |
---|
78 | % new_sos = []; |
---|
79 | SOS = kron(ones(size(SOS,1),1),this_sos) | kron(SOS,ones(size(this_sos,1),1)); |
---|
80 | % for k = 1:size(SOS,1) |
---|
81 | % for r = 1:size(this_sos,1) |
---|
82 | % new_sos = [new_sos;this_sos(r,:) | SOS(k,:)]; |
---|
83 | % end |
---|
84 | % end |
---|
85 | % if ~isequal(new_sos,new_sos2) |
---|
86 | % error |
---|
87 | % end |
---|
88 | % SOS = new_sos; |
---|
89 | end |
---|
90 | variables_in_sos = [variables_in_sos jx j]; |
---|
91 | end |
---|
92 | end |
---|
93 | end |
---|
94 | end |
---|
95 | |
---|
96 | % Detect groups with constraints sum(d_i) == 1 |
---|
97 | for i = 1:size(Aeq_bin,1) |
---|
98 | if beq_bin(i) == 1 |
---|
99 | [ix,jx,sx] = find(Aeq_bin(i,:)); |
---|
100 | if all(sx == 1) |
---|
101 | % Aha, we have sum binary(jx) == 1 |
---|
102 | this_sos = sparse(1:length(jx),jx,1,length(jx),length(binary_var_index)); |
---|
103 | if isempty(SOS) |
---|
104 | SOS = this_sos; |
---|
105 | else |
---|
106 | % new_sos = []; |
---|
107 | SOS = kron(ones(size(SOS,1),1),this_sos) | kron(SOS,ones(size(this_sos,1),1)); |
---|
108 | % for k = 1:size(SOS,1) |
---|
109 | % for r = 1:size(this_sos,1) |
---|
110 | % new_sos = [new_sos;this_sos(r,:) | SOS(k,:)]; |
---|
111 | % end |
---|
112 | % end |
---|
113 | % SOS = new_sos; |
---|
114 | end |
---|
115 | variables_in_sos = [variables_in_sos jx]; |
---|
116 | end |
---|
117 | end |
---|
118 | end |
---|
119 | |
---|
120 | variables_not_in_sos = setdiff(1:length(binary_var_index),variables_in_sos); |
---|
121 | if ~isempty(variables_not_in_sos) |
---|
122 | % we need to add some more binaries |
---|
123 | n_left = length(variables_not_in_sos); |
---|
124 | all_perms = dec2decbin(0:2^n_left-1,n_left); |
---|
125 | [ix,jx,sx] = find(all_perms);jx = variables_not_in_sos(jx); |
---|
126 | this_sos = sparse(ix,jx,sx,size(all_perms,1),length(binary_var_index)); |
---|
127 | if isempty(SOS) |
---|
128 | SOS = this_sos; |
---|
129 | else |
---|
130 | new_sos = []; |
---|
131 | for k = 1:size(SOS,1) |
---|
132 | for r = 1:size(this_sos,1) |
---|
133 | new_sos = [new_sos;this_sos(r,:) | SOS(k,:)]; |
---|
134 | end |
---|
135 | end |
---|
136 | SOS = new_sos; |
---|
137 | end |
---|
138 | end |
---|
139 | |
---|
140 | enums = unique(SOS,'rows')'; |
---|
141 | |
---|
142 | % FIX :FULL, F**N Linux 6.5 repmat bug |
---|
143 | feasible = find(~any(Aeq_bin*enums-repmat(full(beq_bin),1,size(enums,2)),1)); |
---|
144 | enums = enums(:,feasible); |
---|
145 | |
---|
146 | % Remove binary vertices violating these |
---|
147 | if ~isempty(Abin) |
---|
148 | remove = find(any(Abin*enums - repmat(bbin,1,size(enums,2)) >0,1)); |
---|
149 | % keep = setdiff(1:size(enums,2),remove); |
---|
150 | % mapper = find(ismember(1:size(enums,2),keep)); |
---|
151 | % keep(remove)=0;keep2 = 1:size(enums,2);keep2(find(keep)) = find(keep); |
---|
152 | % map = find(keep); |
---|
153 | % remove = find(ismember(j,remove)); |
---|
154 | % i(remove) = []; |
---|
155 | % j(remove) = []; |
---|
156 | % k(remove) = []; |
---|
157 | % enums = sparse(i, |
---|
158 | % |
---|
159 | % [i,j,k] = find(enums); |
---|
160 | enums(:,remove)=[]; |
---|
161 | end |
---|