[37] | 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 |
---|