sgvm_extract_perm

SGVM_EXTRACT_PERM Map matrix A to a permutation matrix (WARINING INCOMPLETE)
   PERM = SGVM_EXTRACT_PERM(A, OPT)

0001 function perm = sgvm_extract_perm(A, opt)
0002 %SGVM_EXTRACT_PERM Map matrix A to a permutation matrix (WARINING INCOMPLETE)
0003 %   PERM = SGVM_EXTRACT_PERM(A, OPT)
0004 
0005 %   SynGrid
0006 %   Copyright (c) 2018, Power Systems Engineering Research Center (PSERC)
0007 %   by Eran Schweitzer, Arizona State University
0008 %
0009 %   This file is part of SynGrid.
0010 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0011 
0012 n       = size(A, 1);
0013 [r,c,v] = find(A); % find non-zero entries
0014 
0015 % remove very small entries
0016 mask = v > 1e-3;
0017 r = r(mask); c = c(mask); v = v(mask);
0018 
0019 prob = struct();
0020 prob.A  = [];
0021 prob.l  = [];
0022 prob.u  = [];
0023 
0024 total_vars = length(r);
0025 
0026 % ---- sum of elements in a single row = 1 -----
0027 A  = sparse(r.', 1:total_vars, 1, n , total_vars);
0028 lb = ones(n, 1);
0029 ub = ones(n, 1);
0030 prob = update_Albub(prob, A, lb, ub);
0031 
0032 % ---- sum of elements in a single column = 1 -----
0033 A  = sparse(c.', 1:total_vars, 1, n , total_vars);
0034 lb = ones(n, 1);
0035 ub = ones(n, 1);
0036 prob = update_Albub(prob, A, lb, ub);
0037 
0038 % ---- maximize v'*u ---> minimize -v'*u
0039 prob.c = -v;
0040 
0041 % ---- initialize with all variables = 1
0042 prob.x0 = ones(total_vars, 1);
0043 
0044 % --- upper and lower bound
0045 prob.xmin = zeros(total_vars,1);
0046 prob.xmax = ones(total_vars,1);
0047 
0048 % --- variable type is binary
0049 prob.vtype = 'B';
0050 
0051 % --- solve
0052 prob.opt.verbose = opt.verbose;
0053 [x, f, eflag, output, lambda] = miqps_matpower(prob);
0054 if eflag
0055     perm = zeros(n, 1);
0056     mask = abs(x) > 0.1; % make sure solution is a logical array
0057     perm(r(mask)) = c(mask);
0058     
0059     if any(perm == 0)
0060         error('sgvm_extract_perm: some rows were not assigned in the permutation.')
0061     end
0062     if ~all(sort(perm) == (1:n).')
0063         error('sgvm_extract_perm: some columns were not assigned in the permutation.')
0064     end
0065 else
0066     error('sgvm_extract_perm: the MILP failed to converge.')
0067 end
0068 
0069 return
0070 %% old greedy approach
0071 perm = zeros(size(A,1),1);
0072 d = [];
0073 available = true(size(A,1),1);
0074 [~, idx] = sort(A(:), 'descend');
0075 ptr = 1;
0076 while any(perm == 0)
0077     [r, c] = ind2sub(size(A), idx(ptr));
0078     if (perm(r) == 0) && available(c)
0079         perm(r) = c;
0080         available(c) = false;
0081     else
0082         d = append_d(d, r, c, A(r,c));
0083     end
0084     ptr = ptr + 1;
0085 end
0086 
0087 function d = append_d(d, r, c, v)
0088 if isempty(d)
0089     d = [r, c, v];
0090 else
0091     d = [d; r, c, v];
0092 end
0093 
0094 function prob = update_Albub(prob, A, lb, ub)
0095 % simple utility function for updating the constraint and bounds arrays.
0096 
0097 prob.A  = [prob.A; A];
0098 prob.l  = [prob.l; lb];
0099 prob.u  = [prob.u; ub];