diff --git a/advancedUse/experiment/POPrandom/solvePOParrow.m b/advancedUse/experiment/POPrandom/solvePOParrow.m index 8539764..eed1fc3 100644 --- a/advancedUse/experiment/POPrandom/solvePOParrow.m +++ b/advancedUse/experiment/POPrandom/solvePOParrow.m @@ -54,29 +54,14 @@ function solvePOParrow(printFileName,degree,a,b,c,l,isBin,addComplement,solver) % end if nargin == 4 - solver = 'BP'; + solver = 'BP'; end % [objPoly, I01, Icomp]=genPOParrow2(d,a,b,c,l,isBin,addComplement) - [objPoly, I01, Icomp] = genPOParrow(degree,a,b,c,l,isBin,addComplement); - - params.solver = solver; - params.newBPswitch = false; - - %% to DNN - eqPolySys = []; - params.relaxOrder = ceil(degree / 2); - [Q0vec, H0vec, H1vec, PSDcone, polyCone, UbdIX] ... - = BBPOPtoCOP(objPoly, eqPolySys, I01, Icomp, params); + [objPoly, I01, Icomp] = genPOParrow(degree,a,b,c,l,isBin,addComplement); + + [sol, info, params] = solvePOPpipeline(objPoly, I01, Icomp, degree, solver); - %% solve - params.UbdIX = UbdIX; - params.UbdObjVal = evalPoly(objPoly, zeros(n, 1)); - params.delta1 = 1.0e-4; - params.delta2 = 0; - params.delta = 0; - [sol, info] = solveCOP(Q0vec, H0vec, H1vec, PSDcone, polyCone, params); - %% print filename = printFileName; fid = fopen(filename, 'a'); diff --git a/advancedUse/experiment/POPrandom/solvePOPchordal.m b/advancedUse/experiment/POPrandom/solvePOPchordal.m index 170028d..c1f6368 100644 --- a/advancedUse/experiment/POPrandom/solvePOPchordal.m +++ b/advancedUse/experiment/POPrandom/solvePOPchordal.m @@ -13,40 +13,14 @@ function solvePOPchordal(printFileName,degree,nDim,radiorange,isBin,addComplemen % addComplement : 1 for adding random complementarity constraints % solver : 'BP, 'sdpnalplus', 'sedumi' -[objPoly, I01, Icomp] = genPOPchordal(degree,nDim,radiorange,isBin,addComplement); - -% Icomp = []; -% if addComplement -% Icomp = complementarityConstraint(2, 2, clique, randSeed); -% end - if nargin == 4 - solver = 'BP'; + solver = 'BP'; end - - if strcmp(solver, 'newBP') - params.solver = 'BP'; - params.newBPswitch = true; - else - params.solver = solver; - params.newBPswitch = false; - end - - %% to DNN - eqPolySys = []; - params.relaxOrder = ceil(degree / 2); - [Q0vec, H0vec, H1vec, PSDcone, polyCone, UbdIX] ... - = BBPOPtoCOP(objPoly, eqPolySys, I01, Icomp, params); - %% solve - params.UbdIX = UbdIX; - params.UbdObjVal = evalPoly(objPoly, zeros(nDim, 1)); - params.delta1 = 1.0e-4; - params.delta2 = 0; - params.delta = 0; - [sol, info] = solveCOP(Q0vec, H0vec, H1vec, PSDcone, polyCone, params); - - + [objPoly, I01, Icomp] = genPOPchordal(degree,nDim,radiorange,isBin,addComplement); + + [sol, info, params] = solvePOPpipeline(objPoly, I01, Icomp, degree, solver); + %% print fid = fopen(printFileName, 'a'); % fprintf(fid,'degree,n,rRange,isBin,addComp,rOder,solver,LBv,sec,iter(APGR),termcode,iterBP\n'); diff --git a/advancedUse/experiment/POPrandom/solvePOPdense.m b/advancedUse/experiment/POPrandom/solvePOPdense.m index c3baf35..3234f4a 100644 --- a/advancedUse/experiment/POPrandom/solvePOPdense.m +++ b/advancedUse/experiment/POPrandom/solvePOPdense.m @@ -10,35 +10,13 @@ function solvePOPdense(printFileName,degree,nDim,isBin,addComplement,solver) % solver : 'BP, 'sdpnalplus', 'sedumi' if nargin == 4 - solver = 'BP'; + solver = 'BP'; end - - [objPoly,I01,Icomp] = genPOPdense(degree,nDim,isBin,addComplement); - if strcmp(solver, 'newBP') - params.solver = 'BP'; - params.newBPswitch = true; - else - params.solver = solver; - params.newBPswitch = false; - end - %% to DNN - eqPolySys = []; - params.relaxOrder = ceil(degree / 2); - [Q0vec, H0vec, H1vec, PSDcone, polyCone, UbdIX] ... - = BBPOPtoCOP(objPoly, eqPolySys, I01, Icomp, params); - %% solve - params.delta1 = 1.0e-4; - params.delta2 = 0; - params.delta = 0; - params.UbdIX = UbdIX; - params.UbdObjVal = evalPoly(objPoly, zeros(nDim, 1)); - % params.deltaRelBP = 1.0e-4; - % params.deltaAbsBP = 0; - % params.deltaDifRelBP = 0; - [sol, info] = solveCOP(Q0vec, H0vec, H1vec, PSDcone, polyCone, params); - - + [objPoly, I01, Icomp] = genPOPdense(degree,nDim,isBin,addComplement); + + [sol, info, params] = solvePOPpipeline(objPoly, I01, Icomp, degree, solver); + %% print fid = fopen(printFileName, 'a'); % fprintf(fid,'degree,n,isBin,addComp,rOrder,solver,LBv,sec,iter,termcode\n'); diff --git a/advancedUse/experiment/POPrandom/solvePOPpipeline.m b/advancedUse/experiment/POPrandom/solvePOPpipeline.m new file mode 100644 index 0000000..c9e06b1 --- /dev/null +++ b/advancedUse/experiment/POPrandom/solvePOPpipeline.m @@ -0,0 +1,21 @@ +function [sol, info, params] = solvePOPpipeline(objPoly, I01, Icomp, degree, solver) +%SOLVEPOPPIPELINE Shared solver pipeline: BBPOPtoCOP -> solveCOP. +% Common logic extracted from solvePOParrow, solvePOPchordal, solvePOPdense. + if strcmp(solver, 'newBP') + params.solver = 'BP'; + params.newBPswitch = true; + else + params.solver = solver; + params.newBPswitch = false; + end + eqPolySys = []; + params.relaxOrder = ceil(degree / 2); + [Q0vec, H0vec, H1vec, PSDcone, polyCone, UbdIX] ... + = BBPOPtoCOP(objPoly, eqPolySys, I01, Icomp, params); + params.UbdIX = UbdIX; + params.UbdObjVal = evalPoly(objPoly, zeros(objPoly.dimVar, 1)); + params.delta1 = 1.0e-4; + params.delta2 = 0; + params.delta = 0; + [sol, info] = solveCOP(Q0vec, H0vec, H1vec, PSDcone, polyCone, params); +end diff --git a/instances/BIQ/biqread.m b/instances/BIQ/biqread.m index 3da74b8..b63fc48 100644 --- a/instances/BIQ/biqread.m +++ b/instances/BIQ/biqread.m @@ -6,14 +6,7 @@ function [Q] = biqread(instance) - fid = fopen(instance,'r'); - if (fid == -1); error('file cannot be opened'); end - - [datavec,count] = fscanf(fid,'%c'); - fclose(fid); - linefeeds = findstr(datavec,char(10)); - datavec(linefeeds) = blanks(length(linefeeds)); - datavec = sscanf(datavec,'%f'); + datavec = readInstanceFile(instance); n = datavec(1); nz = datavec(2); tmp = datavec(3:end); diff --git a/instances/POPrandom/findChordalCliques.m b/instances/POPrandom/findChordalCliques.m new file mode 100644 index 0000000..04248b4 --- /dev/null +++ b/instances/POPrandom/findChordalCliques.m @@ -0,0 +1,40 @@ +function [clique, R] = findChordalCliques(adjMat) +% FINDCHORDALCLIQUES Find maximal cliques via chordal extension (Cholesky). +% Input: +% adjMat : symmetric adjacency matrix of a sparse graph +% Output: +% clique.Set : logical matrix; each row is a maximal clique (0-1 encoding) +% clique.NoC : number of maximal cliques +% clique.maxC : maximum clique size +% clique.minC : minimum clique size +% R : upper triangular adjacency of the chordal extension + + nDim = size(adjMat,1); + adjMat = spones(adjMat); + adjMat = adjMat + adjMat'; + adjMat = adjMat + 10*nDim*speye(nDim); + + I = symamd(adjMat); + [R, p] = chol(adjMat(I,I)); + if (p > 0) + error('Correlative sparsity matrix is not positive definite.'); + end + + R = spones(R); + [~, orig_idx] = sort(I); + remainIdx = 1; + for i = 2:nDim + checkSet = R(i,i:nDim); + noOfone = length(find(checkSet)); + cliqueResult = R(1:i-1,i:nDim)*checkSet'; + yesno = find(cliqueResult == noOfone); + if ~any(yesno) + remainIdx = [remainIdx; i]; + end + end + clique.Set = logical(R(remainIdx, orig_idx)); + clique.NoC = size(clique.Set,1); + sumClique = sum(clique.Set,2); + clique.maxC = max(sumClique); + clique.minC = min(sumClique); +end diff --git a/instances/POPrandom/genPOPchordal.m b/instances/POPrandom/genPOPchordal.m index 5026d18..fcfd998 100644 --- a/instances/POPrandom/genPOPchordal.m +++ b/instances/POPrandom/genPOPchordal.m @@ -8,9 +8,9 @@ %% radiorange : radiorange which controls the sparsity of the chordal graph generated %% for example, radiorange = 0.1 %% isBin : 1 for the binary constraint, 0 for the box constraint -%% addComplement : 1 for adding random complementarity constraints +%% addComplement : 1 for adding random complementarity constraints - relaxOrder = ceil(degree/2); + relaxOrder = ceil(degree/2); params.delta1 = 1.0e-4; params.delta = 0; @@ -18,32 +18,16 @@ I01 = isBin .* true(1,nDim); rng('default'); rng(randSeed); - [adjMat] = generateSparseGraph(radiorange,nDim,randSeed); - adjMat = spones(adjMat); - [clique,R] = genCliques(adjMat); - % clique: information on the maximal cliques, this variable is the - % structure with the following fields. - % - % clique.Set: The maximal clieques, which are represented by 0-1. - % e.g. The number of variable is 3 and the maximal - % cliques are C1={1,2} and C2={2,3}. - % - % Then, clique.Set(1,:) = [1,1,0];(= C1) - % clique.Set(2,:) = [0,1,1];(= C2) - % clique.NoC: The number of cliques - % clique.maxC: The maximum size in the maximum cliques - % clique.minC: The minimum size in the maximum cliques - % R: The upper triangular adjacency matrix of a chordal graph - - %%%%%%%%%% + [adjMat] = generateSparseGraph(radiorange,nDim,randSeed); + adjMat = spones(adjMat); + [clique,R] = findChordalCliques(adjMat); + Icomp = []; if addComplement Icomp = complementarityConstraint(2, 2, clique, randSeed); end - %%%%%%%%% - + if 1 - %spy(R); [m, n] = size(clique.Set); A = sparse(n, n); for ii = 1:m @@ -51,17 +35,10 @@ end spy(A); end - - mDim = size(clique.Set,1); -% supportsT = []; + + mDim = size(clique.Set,1); supports = cell(mDim, 1); - for i=1:mDim -% idx=find(clique.Set(i,:)); -% cDim = length(idx); -% oneSupports = enumMonosUptoDeg(cDim,degree); -% sDim = size(oneSupports,1); -% [I,J,V] = find(oneSupports); -% supportsT = [supportsT,sparse(idx(J),I,V,nDim,sDim,length(V))]; + for i = 1:mDim idx = clique.Set(i,:); cDim = sum(idx); oneSupports = enumMonosUptoDeg(cDim, degree); @@ -70,10 +47,8 @@ supports{i}(:, idx) = oneSupports; end supports = cell2mat(supports); -% objPoly.supports = unique(supportsT','rows','first','legacy'); supports = fastUnique(supports); objPoly.supports = sortrows(supports); - %objPoly.supports = unique(supports, 'rows', 'first', 'legacy'); objPoly.noTerms = size(objPoly.supports,1); objPoly.coef = 2*rand(objPoly.noTerms,1) - ones(objPoly.noTerms,1); objPoly.typeCone = 1; @@ -81,448 +56,5 @@ objPoly.dimVar = nDim; objPoly.degree = degree; x0 = zeros(nDim,1); - objPoly.UbdObjVal = evalPoly(objPoly,x0); % evalPolynomials(objPoly,x0); -end - -function [clique,R] = genCliques(adjMat) - % - % adjMat: symmetric adjacency matrix of a sparse graph - % - % clique: information on the maximal cliques, this variable is the - % structure with the following fields. - % - % clique.Set: The maximal clieques, which are represented by 0-1. - % e.g. The number of variable is 3 and the maximal - % cliques are C1={1,2} and C2={2,3}. - % - % Then, clique.Set(1,:) = [1,1,0];(= C1) - % clique.Set(2,:) = [0,1,1];(= C2) - % clique.NoC: The number of cliques - % clique.maxC: The maximum size in the maximum cliques - % clique.minC: The minimum size in the maximum cliques - - nDim = size(adjMat,1); - %% make rmat symmetric - adjMat = spones(adjMat); - adjMat = adjMat+adjMat'; - %% make rmat diagonally dominant - adjMat = adjMat + 10*nDim*speye(nDim); - %% - %% Chordal Extension by Cholesky decomposition - %% - %% minimum degree ordering - I = symamd(adjMat); - %% cholesky decomposition - [R,p] = chol(adjMat(I,I)); - if (p > 0) - error('Correlative sparsity matrix is not positive definite.'); - end - - %% - %% Step3 - %% Finding the maxmal clieques - %% - - %% put 1 for nonzero element of R - R = spones(R); - % Cliques = spones(R); - [value,orig_idx] = sort(I); - remainIdx = 1; - for i=2:nDim - checkSet = R(i,i:nDim); - one = find(checkSet); - noOfone = length(one); - cliqueResult = R(1:i-1,i:nDim)*checkSet'; - yesno = find(cliqueResult == noOfone); - %% - %% Remove the set included into other set. - %% - if ~any(yesno) - remainIdx = [remainIdx;i]; - end - end - clique.Set = logical(R(remainIdx,orig_idx)); - %% - %%Clique Information - %% - clique.NoC = size(clique.Set,1); - sumClique = sum(clique.Set,2); - clique.maxC = max(sumClique); - clique.minC = min(sumClique); - -end - -function [distanceMatrix0]=generateSparseGraph(radiorange,noOfSensors,randSeed) - % A MATLAB program for generating a sensor network localization problem - % Sunyoung Kim, Masakazu Kojima^* and Hayato Waki - % July 28, 2008 - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Input - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % sDim : the dimension of the space in which sensors and anchors are - % located; sDim is either 2 or 3. - sDim=2; - % When sDim = 2, sensors and anchors are located in [0,1]^2. - % When sDim = 3, sensors and anchors are located in [0,1]^3. - % noisyFac : noisy factor - % = 0 --- no noise - noisyFac=0; - % = \sigma > 0 --- noise with the distribution N(0,\sigma). - % radiorange : radio range; - % If \|x_p - \x_q\| <= radio range, a distance (with noise) is given between x_p and x_q. - % If \|x_p - \x_q\| > radio range, no distance is given between x_p and x_q. - % Here all sensors are placed in [0,1]^{sDim} randomely. - % noOfSensors : the number of sensors. - % anchorType : parameters on how we locate anchors. - anchorType=10; - % anchorType = 0 - % all anchors are placed at grid points on the boundary and interior - % of [0,1]^{sDim} - % anchorType = 1 - % all anchors are placed at grid points in the interior of - % [0,1]^{sDim} - % anchorType = 2 - % all anchors are placed randomly in [0,1]^{sDim} - % anchorType = 3 - % sDim + 1 anchors on the origin and the coordinate axis - % anchorType = 4 - % sDim + 1 anchors near the center - % anchorType = 5 - % noOfAnchors are placed randomly + - % anchors are placed at the corners of [0,1]^2 or [0,1]^3 - % anchorType = 10 - % no anchor - % noOfAnchors : the number of anchors - noOfAnchors=0; - % randSeed : a random seed number for a random distribution of sensors and - % anchors when anchorType = 2 - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Output - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % xMatrix0 : an sDim x (noOfSensors + noOfAnchors) matrix to represent - % locations of sensors and anchors, where noOfSensors sensors - % are placed in the first noOfSensors columns, and - % noOfAnchors are place in the last noOfAnchors columns. - % distanceMatrix0 : a noOfSensors x (noOfSensors + noOfAnchors) upper - % triangular matrix to represent distances from sensors - % to sensors and anchors, where the (p,q)th element - % d_{pq} denotes the distance from the pth sensor to the qth - % sensor or anchor. When noisyFac = 0, - % d_{pq} = norm(xMatrix0(:,p) - xMatrix0(:,q)) - % if norm(xMatrix0(:,p) - xMatrix0(:,q)) <= radiorange, - % d_{pq} = 0 - % if norm(xMatrix0(:,p) - xMatrix0(:,q)) > radiorange, - % When \sigma = noisyFac > 0, - % d_{pq} = max(1+epsilon,0.1)*(norm(xMatrix0(:,p) - xMatrix0(:,q)) - % if norm(xMatrix0(:,p) - xMatrix0(:,q)) <= radiorange, - % d_{pq} = 0 - % if norm(xMatrix0(:,p) - xMatrix0(:,q)) > radiorange, - % Here epsilon is a random number chosen from N(0,\sigma). - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - xbd = 1.0; - - % Checking anchorType and noOfAnchors ---> - if ((anchorType < 0) || (anchorType > 5)) && (anchorType ~= 10) - error('anchorType should be 0, 1, 2, 3, 4 and 10.'); - elseif mod(anchorType,1) > 0 - error('anchorType should be 0, 1, 2, 3, 4 and 10.'); - elseif (anchorType == 0) || (anchorType == 1) - if noOfAnchors < power(2,sDim) - error('When anchorType = 0 or 1, noOfAnchors has to be n^{sDim} for some integer n >= 2.'); - else - gridForAnchors = power(noOfAnchors,1/sDim); - if mod(gridForAnchors,1) > 0 - error('When anchorType = 0 or 1, noOfAnchors has to be n^{sDim} for some integer n >= 2.'); - end - end - elseif (anchorType == 2) - if noOfAnchors < sDim+1 - error('noOfAnchors should be not less than sDim+1 when anchorType = 2.'); - end - elseif (anchorType == 3) || (anchorType == 4) - if noOfAnchors ~= sDim+1 - error('noOfAnchors should be sDim+1 when anchorType = 3 or 4.'); - end - elseif anchorType == 5 - if noOfAnchors < 0 - error('noOfAnchors should be a nonnegative integer.'); - end - elseif anchorType == 10 - if noOfAnchors > 0 - error('noOfAnchors should be 0 when anchorType = 10.'); - end - end - % <--- Checking anchorType and noOfAnchors - - startingTime = tic; % cputime; - - rng('default'); - rng(randSeed+1000); - - - % placing anchors on grids ---> - if (anchorType == 0) || (anchorType == 1) - gridForAnchors = power(noOfAnchors,1/sDim); - if noOfAnchors == 1 - noOfAnchors = 1; - anchorMatrix = ones(sDim,1)*0.5; - xMatrix0 = [xbd*rand(sDim,noOfSensors),anchorMatrix]; - else % noOfAnchors > 1 - xMatrix0 = xbd*rand(sDim,noOfSensors+noOfAnchors); - anchorMatrix = sparse(sDim,noOfAnchors); - h0 = xbd/(gridForAnchors+2); - if anchorType == 0 - h = xbd / (gridForAnchors - 1); - origin = zeros(sDim,1); - else - h = (xbd-h0) / (gridForAnchors-1); - origin = (h0/2) * ones(sDim,1); - end - pointer = 0; - if sDim == 2 - for i = 1:gridForAnchors - for j=1:gridForAnchors - pointer = pointer+1; - anchorMatrix(1,pointer) = (i-1)*h + origin(1,1); - anchorMatrix(2,pointer) = (j-1)*h + origin(2,1); - end - end - elseif sDim == 3 - for i = 1:gridForAnchors - for j=1:gridForAnchors - for k=1:gridForAnchors - pointer = pointer+1; - anchorMatrix(1,pointer) = (i-1)*h + origin(1,1); - anchorMatrix(2,pointer) = (j-1)*h + origin(2,1); - anchorMatrix(3,pointer) = (k-1)*h + origin(3,1); - end - end - end - end - xMatrix0 = [xMatrix0(:,1:noOfSensors), anchorMatrix]; - end - % <--- placing anchors on grids - % placing anchors randomly ---> - elseif anchorType == 2 - xMatrix0 = xbd*rand(sDim,noOfSensors+noOfAnchors); - % <--- placing anchors randomly - % placing anchors on the coordinate axes ---> - elseif anchorType == 3 - noOfAnchors = sDim + 1; - xMatrix0 = xbd*rand(sDim,noOfSensors+noOfAnchors); - anchorMatrix = sparse(sDim,noOfAnchors); - pointer = 1; - anchorMatrix(:,pointer) = zeros(sDim,1); - for i=1:sDim - pointer = pointer+1; - anchorMatrix(:,pointer) = zeros(sDim,1); - anchorMatrix(i,pointer) = xbd*0.5; - end - xMatrix0 = [xMatrix0(:,1:noOfSensors), anchorMatrix]; - % <--- placing anchors on the coordinate axes - % placing anchors in a small box near the center ---> - elseif anchorType == 4 - noOfAnchors = sDim + 1; - xMatrix0 = xbd*rand(sDim,noOfSensors+noOfAnchors); - pointer = 1; - anchorMatrix(:,pointer) = 0.5*xbd*ones(sDim,1); - for i=1:sDim - pointer = pointer+1; - anchorMatrix(:,pointer) = 0.5*xbd*ones(sDim,1); - anchorMatrix(i,pointer) = 0.6*xbd; - end - xMatrix0 = [xMatrix0(:,1:noOfSensors), anchorMatrix]; - elseif anchorType == 5 - xMatrix0 = xbd*rand(sDim,noOfSensors); - if sDim == 2 - for i=0:1 - for j=0:1 - xMatrix0 = [xMatrix0,xbd*[i;j]]; - end - end - elseif sDim == 3 - for i=0:1 - for j=0:1 - for k=0:1 - xMatrix0 = [xMatrix0,xbd*[i;j;k]]; - end - end - end - end - xMatrix0 = [xMatrix0,xbd*rand(sDim,noOfAnchors)]; - if sDim == 2 - noOfAnchors = noOfAnchors + 4; - elseif sDim == 3 - noOfAnchors = noOfAnchors + 8; - end - elseif anchorType == 10 - noOfAnchors = 0; - xMatrix0 = xbd*rand(sDim,noOfSensors); - end - % <--- placing anchors in a small box near the center - - randSeed = randSeed + 13; - - [distanceMatrix0] = computeDistance2(noOfSensors,xMatrix0,radiorange,noisyFac,randSeed); - - fprintf('## elapsed time for generating a sensor network problem = %8.2f\n',toc(startingTime)); - - if 0 == 1 - sensorDistMat = distanceMatrix0(:,1:noOfSensors)+distanceMatrix0(:,1:noOfSensors)'+(noOfSensors+1)*speye(noOfSensors,noOfSensors); - permutation = symamd(sensorDistMat); - UMat = chol(sensorDistMat(permutation,permutation)); - figure(1); - spy(UMat + UMat'); - end + objPoly.UbdObjVal = evalPoly(objPoly,x0); end - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% function [distanceMatrix0] = computeDistance(noOfSensors,xMatrix0,radiorange,noisyFac,randSeed) -% -% sDim = size(xMatrix0,1); -% randn('seed',randSeed); -% noOfAnchors = size(xMatrix0,2) - noOfSensors; -% distanceMatrix0 = sparse(noOfSensors,noOfSensors+noOfAnchors); -% -% startingTime = tic; -% -% for r = noOfSensors+1:noOfSensors+noOfAnchors -% for p=1:noOfSensors -% d0 = norm(xMatrix0(:,p)-xMatrix0(:,r)); -% if d0 <= radiorange -% rate = max([1+randn(1,1)*noisyFac,0.1]); -% distanceMatrix0(p,r) = d0*rate; -% end -% end -% end -% -% fprintf('##0 %6.1f\n',toc(startingTime)); -% startingTime = tic; -% -% idx = [1:noOfSensors]; -% zeroPone = repmat(0.1,1,noOfSensors+noOfAnchors); -% for q = idx -% tmpIdx = repmat(q,1,q-1); -% tmpMat = xMatrix0(:,1:q-1)-xMatrix0(:,tmpIdx); -% tmpMat = tmpMat.*tmpMat; -% tmpMat = sum(tmpMat,1); -% tmpMat = sqrt(tmpMat); -% I = find(tmpMat <= radiorange); -% if ~isempty(I) -% s = length(I); -% tmp = [1+randn(1,s)*noisyFac;zeroPone(1,1:s)]; -% tmpMat(1,I) = tmpMat(1,I).* max(tmp,[],1); -% distanceMatrix0(I,q) = tmpMat(1,I)'; -% end -% end -% -% fprintf('##1 %6.1f\n',toc(startingTime)); -% -% return - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -function [distanceMatrix0] ... - = computeDistance2(noOfSensors,xMatrix0,radiorange,noisyFac,randSeed) - noOfAnchors = size(xMatrix0,2) - noOfSensors; - - if noisyFac < 1.0e-10; - % startingTime = tic; - rowIdx = []; - colIdx = []; - value = []; - for r=noOfSensors+1:noOfSensors+noOfAnchors - tmpIdx = repmat(r,1,noOfSensors); - tmpMat = xMatrix0(:,1:noOfSensors)-xMatrix0(:,tmpIdx); - normVector = sqrt((sum(tmpMat .* tmpMat,1))); - rowAdd = find(normVector < radiorange); - lenRowAdd = length(rowAdd); - if lenRowAdd > 0 - rowIdx = [rowIdx,rowAdd]; - colIdx = [colIdx,repmat(r-noOfSensors,1,lenRowAdd)]; - % valueAdd = normVector(rowAdd) .* (ones(1,lenRowAdd)+randn(1,s)*noisyFac); - value = [value,normVector(rowAdd)]; - end - end - senToAnchDistMat = sparse(rowIdx,colIdx,value,noOfSensors,noOfAnchors); - % fprintf('##0 %6.1f\n',toc(startingTime)); - % startingTime = tic; - distanceMatrix0 = sparse(noOfSensors,noOfSensors); - idx = [1:noOfSensors]; - for q = idx - tmpIdx = repmat(q,1,q-1); - tmpMat = xMatrix0(:,1:q-1)-xMatrix0(:,tmpIdx); - tmpMat = tmpMat.*tmpMat; - tmpMat = sum(tmpMat,1); - tmpMat = sqrt(tmpMat); - I = find(tmpMat <= radiorange); - if ~isempty(I) - % s = length(I); - % tmp = [1+randn(1,s)*noisyFac;zeroPone(1,1:s)]; - % tmpMat(1,I) = tmpMat(1,I).* max(tmp,[],1); - distanceMatrix0(I,q) = tmpMat(1,I)'; - end - end - % fprintf('##1 %6.1f\n',toc(startingTime)); - else - randn('seed',randSeed); - zeroPone = repmat(0.1,1,noOfSensors+noOfAnchors); - % startingTime = tic; - rowIdx = []; - colIdx = []; - value = []; - for r=noOfSensors+1:noOfSensors+noOfAnchors - tmpIdx = repmat(r,1,noOfSensors); - tmpMat = xMatrix0(:,1:noOfSensors)-xMatrix0(:,tmpIdx); - normVector = sqrt((sum(tmpMat .* tmpMat,1))); - rowAdd = find(normVector < radiorange); - lenRowAdd = length(rowAdd); - if lenRowAdd > 0 - rowIdx = [rowIdx,rowAdd]; - colIdx = [colIdx,repmat(r-noOfSensors,1,lenRowAdd)]; - rateVector = max([ones(1,lenRowAdd)+randn(1,lenRowAdd)*noisyFac;zeroPone(1,1:lenRowAdd)],[],1); - valueAdd = normVector(rowAdd) .* rateVector; - value = [value,valueAdd]; - end - end - senToAnchDistMat = sparse(rowIdx,colIdx,value,noOfSensors,noOfAnchors); - % fprintf('##0 %6.1f\n',toc(startingTime)); - % startingTime = tic; - distanceMatrix0 = sparse(noOfSensors,noOfSensors); - idx = [1:noOfSensors]; - for q = idx - tmpIdx = repmat(q,1,q-1); - tmpMat = xMatrix0(:,1:q-1)-xMatrix0(:,tmpIdx); - tmpMat = tmpMat.*tmpMat; - tmpMat = sum(tmpMat,1); - tmpMat = sqrt(tmpMat); - I = find(tmpMat <= radiorange); - if ~isempty(I) - s = length(I); - tmp = [1+randn(1,s)*noisyFac;zeroPone(1,1:s)]; - tmpMat(1,I) = tmpMat(1,I).* max(tmp,[],1); - distanceMatrix0(I,q) = tmpMat(1,I)'; - end - end - % fprintf('##1 %6.1f\n',toc(startingTime)); - end - - distanceMatrix0 = [distanceMatrix0, senToAnchDistMat]; - % rowIdx = []; - % colIdx = []; - % value = []; - % - % for p=1:noOfSensors-1 - % tmpIdx = repmat(p,1,noOfSensors-p); - % tmpMat = xMatrix0(:,p+1:noOfSensors)-xMatrix0(:,tmpIdx); - % normVector = sqrt((sum(tmpMat .* tmpMat,1))); - % rowAdd = find(normVector < radiorange); - % if ~isempty(rowAdd) - % rowIdx = [rowIdx,rowAdd+p]; - % colIdx = [colIdx,repmat(p,1,length(rowAdd))]; - % value = [value,normVector(rowAdd)]; - % end - % end - % distanceMatrix0= [sparse(rowIdx,colIdx,value,noOfSensors,noOfSensors)',senToAnchDistMat]; -end - diff --git a/instances/POPrandom/generateSparseGraph.m b/instances/POPrandom/generateSparseGraph.m new file mode 100644 index 0000000..4e8d4b5 --- /dev/null +++ b/instances/POPrandom/generateSparseGraph.m @@ -0,0 +1,143 @@ +function [distanceMatrix0] = generateSparseGraph(radiorange, noOfSensors, randSeed) +% GENERATESPARSEGRAPH Generate a sparse sensor network graph. +% Returns a distance matrix for a randomly placed sensor network where +% only sensors within radiorange of each other are connected. +% +% Input: +% radiorange : radio range threshold for connectivity +% noOfSensors : number of sensors +% randSeed : random seed +% +% Originally part of a sensor network localization problem generator by +% Sunyoung Kim, Masakazu Kojima and Hayato Waki, July 28, 2008. + + sDim = 2; + noisyFac = 0; + anchorType = 10; + noOfAnchors = 0; + + % Validate anchorType and noOfAnchors + if ((anchorType < 0) || (anchorType > 5)) && (anchorType ~= 10) + error('anchorType should be 0, 1, 2, 3, 4 and 10.'); + elseif mod(anchorType,1) > 0 + error('anchorType should be 0, 1, 2, 3, 4 and 10.'); + elseif (anchorType == 0) || (anchorType == 1) + if noOfAnchors < power(2,sDim) + error('When anchorType = 0 or 1, noOfAnchors has to be n^{sDim} for some integer n >= 2.'); + else + gridForAnchors = power(noOfAnchors,1/sDim); + if mod(gridForAnchors,1) > 0 + error('When anchorType = 0 or 1, noOfAnchors has to be n^{sDim} for some integer n >= 2.'); + end + end + elseif (anchorType == 2) + if noOfAnchors < sDim+1 + error('noOfAnchors should be not less than sDim+1 when anchorType = 2.'); + end + elseif (anchorType == 3) || (anchorType == 4) + if noOfAnchors ~= sDim+1 + error('noOfAnchors should be sDim+1 when anchorType = 3 or 4.'); + end + elseif anchorType == 5 + if noOfAnchors < 0 + error('noOfAnchors should be a nonnegative integer.'); + end + elseif anchorType == 10 + if noOfAnchors > 0 + error('noOfAnchors should be 0 when anchorType = 10.'); + end + end + + startingTime = tic; + + rng('default'); + rng(randSeed+1000); + + xbd = 1.0; + + if anchorType == 10 + noOfAnchors = 0; + xMatrix0 = xbd*rand(sDim,noOfSensors); + end + + randSeed = randSeed + 13; + + [distanceMatrix0] = computeDistance2(noOfSensors,xMatrix0,radiorange,noisyFac,randSeed); + + fprintf('## elapsed time for generating a sensor network problem = %8.2f\n',toc(startingTime)); +end + +function [distanceMatrix0] = computeDistance2(noOfSensors, xMatrix0, radiorange, noisyFac, randSeed) + noOfAnchors = size(xMatrix0,2) - noOfSensors; + + if noisyFac < 1.0e-10 + rowIdx = []; + colIdx = []; + value = []; + for r = noOfSensors+1:noOfSensors+noOfAnchors + tmpIdx = repmat(r,1,noOfSensors); + tmpMat = xMatrix0(:,1:noOfSensors)-xMatrix0(:,tmpIdx); + normVector = sqrt((sum(tmpMat .* tmpMat,1))); + rowAdd = find(normVector < radiorange); + lenRowAdd = length(rowAdd); + if lenRowAdd > 0 + rowIdx = [rowIdx,rowAdd]; + colIdx = [colIdx,repmat(r-noOfSensors,1,lenRowAdd)]; + value = [value,normVector(rowAdd)]; + end + end + senToAnchDistMat = sparse(rowIdx,colIdx,value,noOfSensors,noOfAnchors); + distanceMatrix0 = sparse(noOfSensors,noOfSensors); + idx = [1:noOfSensors]; + for q = idx + tmpIdx = repmat(q,1,q-1); + tmpMat = xMatrix0(:,1:q-1)-xMatrix0(:,tmpIdx); + tmpMat = tmpMat.*tmpMat; + tmpMat = sum(tmpMat,1); + tmpMat = sqrt(tmpMat); + I = find(tmpMat <= radiorange); + if ~isempty(I) + distanceMatrix0(I,q) = tmpMat(1,I)'; + end + end + else + randn('seed',randSeed); + zeroPone = repmat(0.1,1,noOfSensors+noOfAnchors); + rowIdx = []; + colIdx = []; + value = []; + for r = noOfSensors+1:noOfSensors+noOfAnchors + tmpIdx = repmat(r,1,noOfSensors); + tmpMat = xMatrix0(:,1:noOfSensors)-xMatrix0(:,tmpIdx); + normVector = sqrt((sum(tmpMat .* tmpMat,1))); + rowAdd = find(normVector < radiorange); + lenRowAdd = length(rowAdd); + if lenRowAdd > 0 + rowIdx = [rowIdx,rowAdd]; + colIdx = [colIdx,repmat(r-noOfSensors,1,lenRowAdd)]; + rateVector = max([ones(1,lenRowAdd)+randn(1,lenRowAdd)*noisyFac;zeroPone(1,1:lenRowAdd)],[],1); + valueAdd = normVector(rowAdd) .* rateVector; + value = [value,valueAdd]; + end + end + senToAnchDistMat = sparse(rowIdx,colIdx,value,noOfSensors,noOfAnchors); + distanceMatrix0 = sparse(noOfSensors,noOfSensors); + idx = [1:noOfSensors]; + for q = idx + tmpIdx = repmat(q,1,q-1); + tmpMat = xMatrix0(:,1:q-1)-xMatrix0(:,tmpIdx); + tmpMat = tmpMat.*tmpMat; + tmpMat = sum(tmpMat,1); + tmpMat = sqrt(tmpMat); + I = find(tmpMat <= radiorange); + if ~isempty(I) + s = length(I); + tmp = [1+randn(1,s)*noisyFac;zeroPone(1,1:s)]; + tmpMat(1,I) = tmpMat(1,I).* max(tmp,[],1); + distanceMatrix0(I,q) = tmpMat(1,I)'; + end + end + end + + distanceMatrix0 = [distanceMatrix0, senToAnchDistMat]; +end diff --git a/instances/QAP/qapread.m b/instances/QAP/qapread.m index ad6bc8c..d705d3b 100644 --- a/instances/QAP/qapread.m +++ b/instances/QAP/qapread.m @@ -5,15 +5,7 @@ filename = [instance '.dat']; - fid = fopen(filename,'r'); - if (fid == -1); error('file cannot be opened'); end - - [datavec, ~] = fscanf(fid, '%c'); - fclose('all'); - - linefeeds = strfind(datavec, char(10)); - datavec(linefeeds) = blanks(length(linefeeds)); - datavec = sscanf(datavec, '%f'); + datavec = readInstanceFile(filename); nn = datavec(1); n2 = nn * nn; diff --git a/instances/QAP/qapreadBP.m b/instances/QAP/qapreadBP.m index 6182619..6d2e77a 100644 --- a/instances/QAP/qapreadBP.m +++ b/instances/QAP/qapreadBP.m @@ -16,15 +16,7 @@ filename = [instance '.dat']; - fid = fopen(filename,'r'); - if (fid == -1); error('file cannot be opened'); end - - [datavec, ~] = fscanf(fid, '%c'); - fclose('all'); - - linefeeds = strfind(datavec, char(10)); - datavec(linefeeds) = blanks(length(linefeeds)); - datavec = sscanf(datavec, '%f'); + datavec = readInstanceFile(filename); nn = datavec(1); n2 = nn * nn; diff --git a/instances/QAP/qapreadnal.m b/instances/QAP/qapreadnal.m index 75ff978..cbf9dde 100644 --- a/instances/QAP/qapreadnal.m +++ b/instances/QAP/qapreadnal.m @@ -11,14 +11,7 @@ fname = [instance '.dat']; - fid = fopen(fname,'r'); - if (fid == -1); error('file cannot be opened'); end - - [datavec,count] = fscanf(fid,'%c'); - fclose('all'); - linefeeds = findstr(datavec,char(10)); - datavec(linefeeds) = blanks(length(linefeeds)); - datavec = sscanf(datavec,'%f'); + datavec = readInstanceFile(fname); n = datavec(1); n2 = n*n; aa = datavec(2:n2+1); diff --git a/relaxation/BBPOPtoCOP.m b/relaxation/BBPOPtoCOP.m index 7579361..6e0a2e8 100644 --- a/relaxation/BBPOPtoCOP.m +++ b/relaxation/BBPOPtoCOP.m @@ -13,7 +13,7 @@ % % Note: % The input ineqPolySys (which must be a SOS function) -% and the corresponding output Qpvecs are used +% and the corresponding output Qpvecs are used % for the Lagrangian-DNN relaxation. @@ -33,7 +33,7 @@ end if params.printyesDNN >= 1 - fprintf('\n----------------------------------------------'); + fprintf('\n----------------------------------------------'); fprintf('\n BBPOPtoCOP'); fprintf('\n----------------------------------------------\n'); end @@ -106,7 +106,7 @@ fprintf('Generating coefficient matrices H0, Q0, Qp...'); stime = tic; end -H0vec = genH0vec(polyCone.varStructure);%sparse(1,1,1,length(Q0vec),1); +H0vec = genH0vec(polyCone.varStructure); Q0vec = genQ0vec(polyCone.varStructure, objPoly); [~, Qpvecs] = genH1vec(polyCone.varStructure, ineqPolySys, params); if params.printyesDNN >= 2 @@ -133,17 +133,13 @@ %% function [H1vec, Qpvecs] = genH1vec(varStructure, ineqPolySys, params) - %idxSupports = varStructure.idxSupports; - %sizeblk = cellfun(@(x) size(x, 1), idxSupports); BMat = varStructure.BMat; sumBMat = sum(BMat, 1); mm = length(ineqPolySys); Qpvecs = zeros(size(BMat, 1), mm); - %Qpvecs = cell(1, mm); for ii = 1:mm avecoef = ineqPolySys{ii}.coef ./ sumBMat(varStructure.sysic{ii})'; Qpvecs(:, ii) = BMat(:, varStructure.sysic{ii}) * avecoef; - %Qpvecs{ii} = BMat(:, varStructure.sysic{ii}) * sparse(avecoef); end if mm == 0 H1vec = sparse(size(BMat, 1), 1); @@ -152,218 +148,11 @@ end end -%% -function UbdIX = estimateUbdIX(polyCone, Icomp, params) - sizeblk = polyCone.varStructure.sizeblk; - if params.LagDNN - UbdIX = sum((0.5 * (sizeblk - 1)) + 1); - return - end - - eta = max(sum(Icomp, 1)); - - vS = polyCone.varStructure; - isDiag = vS.rowIdx == vS.colIdx; - numAppearInDiag = accumarray(vS.momentic(isDiag), 1, [length(vS.ia), 1]); - cardiB = sum(numAppearInDiag .* ~polyCone.eq0Var); - if eta == 0 - UbdIX = cardiB; %sum(numAppearInDiag .* ~eq0Var); - return - end - - - isNNZinDiag = logical(numAppearInDiag) & ~polyCone.eq0Var; - numAppearNNZDiag = numAppearInDiag(isNNZinDiag); - diagSupports = vS.supports(isNNZinDiag, :); - c = @(S) sum(any(diagSupports(:, S), 2) .* numAppearNNZDiag); - f = @(S) sum(any(Icomp(:, S), 2)); - - [~, val] = GSC(c, f, size(Icomp, 2)); - val = full(val); - - Heta = sum(1 ./ (1:eta)); - dd = full(max(sum(spones(diagSupports), 2))); - LB = ceil(val / (dd * Heta)); - UbdIX = floor(cardiB - LB); - if params.printyesDNN >= 2 - fprintf('cardiB: %d, UbdIX:%d, LB:%d, val:%d, dd:%d, Heta:%g\n', cardiB, UbdIX, LB, val, dd, Heta); - end -end - -%% -function [newObjPoly, newIneqPolySys, newIdxSupports, newI01, newIcomp, xIdx] ... - = addSlack(objPoly, ineqPolySys, idxSupports, I01, Icomp) - %ADDSLACK Add slack variable - % - - numXVar = length(I01); - - [uniqueIdx, ~, icIdx] = fastUnique(cell2mat(idxSupports)); - originalBlkSize = cellfun(@(x) size(x, 1), idxSupports); - icIdx = mat2cell(icIdx, originalBlkSize, 1); - numSlk = size(uniqueIdx, 1) - 1; - numSlkInBlk = cellfun(@(x) length(x), icIdx) - 1; - - %% poly - newObjPoly.coef = objPoly.coef; - newObjPoly.supports = [objPoly.supports, sparse(size(objPoly.supports, 1), numSlk)]; - - mm = length(ineqPolySys); - newIneqPolySys = cell(mm + numSlk, 1); - for ii = 1:mm - newIneqPolySys{ii}.coef = ineqPolySys{ii}.coef; - newIneqPolySys{ii}.supports = ... - [ineqPolySys{ii}.supports, sparse(size(ineqPolySys{ii}.supports, 1), numSlk)]; - end - %% addeq - % is it ok? - for ii = 1:numSlk - newIneqPolySys{ii + mm}.coef = [-1; 1; 1]; - newIneqPolySys{ii + mm}.supports = [sparse(1, numXVar + numSlk); %const - uniqueIdx(ii + 1, :), sparse(1, numSlk); % x^alpha - sparse(1, numXVar + ii, 1, 1, numXVar + numSlk)]; % slk_alpha - end - - %% idxSupports - newIdxSupports = cell(length(idxSupports), 1); - for ii = 1:length(icIdx) - newIdxSupports{ii} = ... - [idxSupports{ii}, sparse(originalBlkSize(ii), numSlk); ... - sparse(numSlkInBlk(ii), numXVar), ... - sparse(1:numSlkInBlk(ii), icIdx{ii}(2:end) - 1, 1, numSlkInBlk(ii), numSlk)]; - end - %newBlkSize = cellfun(@(x) size(x, 1), newIdxSupports); - - %% I01 - slkI01 = ~any(uniqueIdx(2:end, ~I01), 2)'; % binary = there is no non-binary variable - newI01 = [I01, slkI01]; - - %% Icomp - forcomp = find(slkI01); - newIcomp = logical([Icomp, sparse(size(Icomp, 1), numSlk);... - uniqueIdx(forcomp + 1, :),... - sparse(1:length(forcomp), forcomp, 1, length(forcomp), numSlk)]); - - %% xIdx - xIdx = [true(1, numXVar), false(1, numSlk)]; -end - %% function PSDcone = genPSDcone(varStructure) PSDcone = varStructure.sizeblk; end -%% -function [polyCone, ineqPolySys] = genPolyCone(objPoly, ineqPolySys, I01, Icomp, params) - numVarX = size(objPoly.supports, 2); - eqSupports = cellfun(@(poly) poly.supports, ineqPolySys, 'UniformOutput', false); - allSupports = [objPoly.supports; Icomp; cell2mat(eqSupports)]; - cliques = genCliques(allSupports, params); % mathcal{C}_k - -% if params.sparseSW && params.cliqueMIP -% cliques = cliqueMIP(cliques, 100); -% end - - sparseIdxSupports = genIdxSupports(cliques, params.relaxOrder, I01); % mathcal{A}_k - sparseEigTime = sum(cellfun(@(x) size(x, 1), sparseIdxSupports).^3); - - denseIdxSupports = genIdxSupports(true(1, numVarX), params.relaxOrder, I01); - denseEigTime = size(denseIdxSupports{1}, 1)^3; - if denseEigTime < sparseEigTime - idxSupports = denseIdxSupports; - else - idxSupports = sparseIdxSupports; - end - - if params.LagDNN - [objPoly, ineqPolySys, idxSupports, I01, Icomp, xIdx] ... - = addSlack(objPoly, ineqPolySys, idxSupports, I01, Icomp); - else - xIdx = true(1, numVarX); - end - - lowTriMomentMat = genLowTriMomentMat(idxSupports, I01); % blk,row,col,support - varStructure = findCommonMonos(lowTriMomentMat, objPoly, ineqPolySys, numVarX); - varStructure.sizeblk = countSizeblk(idxSupports); - varStructure.BMat = genBasisMat(varStructure); - varStructure.idxSupports = idxSupports; - varStructure.xIdx = xIdx; - - [chain,eq0Var,lbd,ubd] ... - = genEqIneq(varStructure.supports, I01, Icomp, params.DNNSW); - if params.findArborescence - chain = findLongChain(varStructure.supports, eq0Var, chain); - chain = findArborescence(varStructure.supports, eq0Var, chain); - elseif params.findLongChain - chain = findLongChain(varStructure.supports, eq0Var, chain); - end - polyCone.chain = chain; - polyCone.singleton = ~any(chain, 2); - polyCone.eq0Var = eq0Var; - polyCone.lbd = lbd; - polyCone.ubd = ubd; %??? - polyCone.nonneg = params.DNNSW; %for sdp solver - polyCone.varStructure = varStructure; - polyCone.LagDNN = params.LagDNN; - - % for proj - [permIdx, ~] = find(chain); - cardinality = full(sum(varStructure.BMat, 1))'; - polyCone.varStructure.cardinality = cardinality; - polyCone.cardinalityPerm = cardinality(permIdx)'; - - polyCone.BMatPerm = varStructure.BMat(:, permIdx); - - lengthEachChain = full(sum(chain, 1)); - polyCone.lengthEachChain = lengthEachChain; - polyCone.I = cumsum([0, lengthEachChain]); - polyCone.topIdx = polyCone.I + 1; - polyCone.tailIdx = lengthEachChain; - polyCone.permIdx = permIdx; -end - -%% -function sizeblk = countSizeblk(idxSupports) - numblk = length(idxSupports); - sizeblk = zeros(1, numblk); - for ii = 1:numblk - sizeblk(ii) = size(idxSupports{ii}, 1); - end -end - -%% -function varStructure = findCommonMonos(momentMat, objPoly, ineqPolySys, numVarX) - - numVar = size(objPoly.supports, 2); % with slack - eqSupports = cellfun(@(poly) poly.supports, ineqPolySys, 'UniformOutput', false); - [varStructure.supports, varStructure.ia, ic] =... - fastUnique([momentMat.supports; objPoly.supports; cell2mat(eqSupports); ... - speye(numVarX), sparse(numVarX, numVar - numVarX)]); - - numMoment = size(momentMat.supports, 1); - numObj = size(objPoly.supports, 1); - numSys = cellfun(@(supp) size(supp, 1), eqSupports, 'UniformOutput', false); - - numBlk = [numMoment; numObj; cell2mat(numSys); numVarX]; - icBlk = mat2cell(ic, numBlk, 1); - - %varStructure.momentic = ic(1:numMoment); - %varStructure.objic = ic(numMoment + (1:numObj)); - %varStructure.xic = ic(numMoment + numObj + (1:numVarX)); - varStructure.momentic = icBlk{1}; - varStructure.objic = icBlk{2}; - varStructure.sysic = cell(length(numSys), 1); - for ii = 1:length(numSys) - varStructure.sysic{ii} = icBlk{ii+2}; - end - varStructure.xic = icBlk{end}; - varStructure.rowIdx = momentMat.rowIdx; - varStructure.colIdx = momentMat.colIdx; - varStructure.blkIdx = momentMat.blkIdx; - varStructure.computeMomentIdx = momentMat.computeMomentIdx; - -end - %% function Q0vec = genQ0vec(varStructure, objPoly) sumBMat = sum(varStructure.BMat, 1); @@ -372,8 +161,7 @@ end %% -% function H0vec = genH0vec(varStructure) bb = varStructure.BMat(:, 1); H0vec = (1 / sum(bb)) * bb; -end \ No newline at end of file +end diff --git a/relaxation/addSlack.m b/relaxation/addSlack.m new file mode 100644 index 0000000..27c8461 --- /dev/null +++ b/relaxation/addSlack.m @@ -0,0 +1,53 @@ +function [newObjPoly, newIneqPolySys, newIdxSupports, newI01, newIcomp, xIdx] ... + = addSlack(objPoly, ineqPolySys, idxSupports, I01, Icomp) +%ADDSLACK Add slack variable for Lagrangian-DNN relaxation. + + numXVar = length(I01); + + [uniqueIdx, ~, icIdx] = fastUnique(cell2mat(idxSupports)); + originalBlkSize = cellfun(@(x) size(x, 1), idxSupports); + icIdx = mat2cell(icIdx, originalBlkSize, 1); + numSlk = size(uniqueIdx, 1) - 1; + numSlkInBlk = cellfun(@(x) length(x), icIdx) - 1; + + %% poly + newObjPoly.coef = objPoly.coef; + newObjPoly.supports = [objPoly.supports, sparse(size(objPoly.supports, 1), numSlk)]; + + mm = length(ineqPolySys); + newIneqPolySys = cell(mm + numSlk, 1); + for ii = 1:mm + newIneqPolySys{ii}.coef = ineqPolySys{ii}.coef; + newIneqPolySys{ii}.supports = ... + [ineqPolySys{ii}.supports, sparse(size(ineqPolySys{ii}.supports, 1), numSlk)]; + end + %% addeq + for ii = 1:numSlk + newIneqPolySys{ii + mm}.coef = [-1; 1; 1]; + newIneqPolySys{ii + mm}.supports = [sparse(1, numXVar + numSlk); %const + uniqueIdx(ii + 1, :), sparse(1, numSlk); % x^alpha + sparse(1, numXVar + ii, 1, 1, numXVar + numSlk)]; % slk_alpha + end + + %% idxSupports + newIdxSupports = cell(length(idxSupports), 1); + for ii = 1:length(icIdx) + newIdxSupports{ii} = ... + [idxSupports{ii}, sparse(originalBlkSize(ii), numSlk); ... + sparse(numSlkInBlk(ii), numXVar), ... + sparse(1:numSlkInBlk(ii), icIdx{ii}(2:end) - 1, 1, numSlkInBlk(ii), numSlk)]; + end + + %% I01 + slkI01 = ~any(uniqueIdx(2:end, ~I01), 2)'; + newI01 = [I01, slkI01]; + + %% Icomp + forcomp = find(slkI01); + newIcomp = logical([Icomp, sparse(size(Icomp, 1), numSlk); ... + uniqueIdx(forcomp + 1, :), ... + sparse(1:length(forcomp), forcomp, 1, length(forcomp), numSlk)]); + + %% xIdx + xIdx = [true(1, numXVar), false(1, numSlk)]; +end diff --git a/relaxation/findCommonMonos.m b/relaxation/findCommonMonos.m new file mode 100644 index 0000000..e2ce15e --- /dev/null +++ b/relaxation/findCommonMonos.m @@ -0,0 +1,28 @@ +function varStructure = findCommonMonos(momentMat, objPoly, ineqPolySys, numVarX) +%FINDCOMMONMONOS Collect all unique monomials across moment matrix and polynomials. + + numVar = size(objPoly.supports, 2); % with slack + eqSupports = cellfun(@(poly) poly.supports, ineqPolySys, 'UniformOutput', false); + [varStructure.supports, varStructure.ia, ic] = ... + fastUnique([momentMat.supports; objPoly.supports; cell2mat(eqSupports); ... + speye(numVarX), sparse(numVarX, numVar - numVarX)]); + + numMoment = size(momentMat.supports, 1); + numObj = size(objPoly.supports, 1); + numSys = cellfun(@(supp) size(supp, 1), eqSupports, 'UniformOutput', false); + + numBlk = [numMoment; numObj; cell2mat(numSys); numVarX]; + icBlk = mat2cell(ic, numBlk, 1); + + varStructure.momentic = icBlk{1}; + varStructure.objic = icBlk{2}; + varStructure.sysic = cell(length(numSys), 1); + for ii = 1:length(numSys) + varStructure.sysic{ii} = icBlk{ii+2}; + end + varStructure.xic = icBlk{end}; + varStructure.rowIdx = momentMat.rowIdx; + varStructure.colIdx = momentMat.colIdx; + varStructure.blkIdx = momentMat.blkIdx; + varStructure.computeMomentIdx = momentMat.computeMomentIdx; +end diff --git a/relaxation/genPolyCone.m b/relaxation/genPolyCone.m new file mode 100644 index 0000000..f09ada1 --- /dev/null +++ b/relaxation/genPolyCone.m @@ -0,0 +1,72 @@ +function [polyCone, ineqPolySys] = genPolyCone(objPoly, ineqPolySys, I01, Icomp, params) +%GENPOLYCONE Generate the polynomial cone (L) for the DNN relaxation. + + numVarX = size(objPoly.supports, 2); + eqSupports = cellfun(@(poly) poly.supports, ineqPolySys, 'UniformOutput', false); + allSupports = [objPoly.supports; Icomp; cell2mat(eqSupports)]; + cliques = genCliques(allSupports, params); + + sparseIdxSupports = genIdxSupports(cliques, params.relaxOrder, I01); + sparseEigTime = sum(cellfun(@(x) size(x, 1), sparseIdxSupports).^3); + + denseIdxSupports = genIdxSupports(true(1, numVarX), params.relaxOrder, I01); + denseEigTime = size(denseIdxSupports{1}, 1)^3; + if denseEigTime < sparseEigTime + idxSupports = denseIdxSupports; + else + idxSupports = sparseIdxSupports; + end + + if params.LagDNN + [objPoly, ineqPolySys, idxSupports, I01, Icomp, xIdx] ... + = addSlack(objPoly, ineqPolySys, idxSupports, I01, Icomp); + else + xIdx = true(1, numVarX); + end + + lowTriMomentMat = genLowTriMomentMat(idxSupports, I01); + varStructure = findCommonMonos(lowTriMomentMat, objPoly, ineqPolySys, numVarX); + varStructure.sizeblk = countSizeblk(idxSupports); + varStructure.BMat = genBasisMat(varStructure); + varStructure.idxSupports = idxSupports; + varStructure.xIdx = xIdx; + + [chain, eq0Var, lbd, ubd] ... + = genEqIneq(varStructure.supports, I01, Icomp, params.DNNSW); + if params.findArborescence + chain = findLongChain(varStructure.supports, eq0Var, chain); + chain = findArborescence(varStructure.supports, eq0Var, chain); + elseif params.findLongChain + chain = findLongChain(varStructure.supports, eq0Var, chain); + end + polyCone.chain = chain; + polyCone.singleton = ~any(chain, 2); + polyCone.eq0Var = eq0Var; + polyCone.lbd = lbd; + polyCone.ubd = ubd; + polyCone.nonneg = params.DNNSW; + polyCone.varStructure = varStructure; + polyCone.LagDNN = params.LagDNN; + + [permIdx, ~] = find(chain); + cardinality = full(sum(varStructure.BMat, 1))'; + polyCone.varStructure.cardinality = cardinality; + polyCone.cardinalityPerm = cardinality(permIdx)'; + + polyCone.BMatPerm = varStructure.BMat(:, permIdx); + + lengthEachChain = full(sum(chain, 1)); + polyCone.lengthEachChain = lengthEachChain; + polyCone.I = cumsum([0, lengthEachChain]); + polyCone.topIdx = polyCone.I + 1; + polyCone.tailIdx = lengthEachChain; + polyCone.permIdx = permIdx; +end + +function sizeblk = countSizeblk(idxSupports) + numblk = length(idxSupports); + sizeblk = zeros(1, numblk); + for ii = 1:numblk + sizeblk(ii) = size(idxSupports{ii}, 1); + end +end diff --git a/util/readInstanceFile.m b/util/readInstanceFile.m new file mode 100644 index 0000000..c8ac035 --- /dev/null +++ b/util/readInstanceFile.m @@ -0,0 +1,13 @@ +function datavec = readInstanceFile(filename) +%READINSTANCEFILE Read a space/newline-delimited numeric instance file. +% Returns a column vector of all numeric values in the file. + fid = fopen(filename, 'r'); + if fid == -1 + error('Cannot open file: %s', filename); + end + [rawdata, ~] = fscanf(fid, '%c'); + fclose(fid); + linefeeds = strfind(rawdata, char(10)); + rawdata(linefeeds) = blanks(length(linefeeds)); + datavec = sscanf(rawdata, '%f'); +end