function [X,S]=TriangulationGsedumi(P,E,Eo) % Test a framework with configuration/position P, graph % specified by edge set E, and an objective function specified by by % edges in Eo. The distance of each edge in E will be set as a % constraint; and the sum of edge distances in Eo will be maximized. % E could represent a triangulation graph, and Eo represents edge set % of all virtual edges as in Shamsi, Taheri, Zhu and Ye 2011. % % coded on 12/29/2011. % % Input: % P: d x n matrix where jth column is the position of jth point in R^d % E: 2 x |E| matrix where each column specifies two point-idexes of % an edge which Euclidean distance is set as a constraint. % Eo: 2 x |Eo| matrix where each column specifies two point-idexes of % an edge to be maximized in the sum of all edges in Eo. % % Output: % X: d' x n configuration matrix satisfying the edge distance % constraint, where the first point in X is the origin. % Thus, the second column of X is the positiion of the second point in % P, and so on. % % d'=d if and only if the framework is uniquely localizable with the % added objective function. % % S: a max-rank PSD (dual) stress matrix approximately orthogonal to X % % Solver: sedumi 1.1 by Jos Stern % % Check the dimension of the configuration and the number of edges % [d,n]=size(P); n=n-1; [m,mtotal]=size(E); toler = 1.e-8; % % Compute edge distances specified by E and set up the SDP inputs % b=[]; A=[]; for k =1:mtotal i = E(1,k); j = E(2,k); a=sparse(zeros(n,1)); if (i ~= 1), a(i-1)=1;end; if (j ~= 1), a(j-1)=-1;end; A=[A vec(a*a')]; dd= P(:,i)-P(:,j); rr= dd'*dd; b=[b; rr]; end; % % Set up the SDP objective using edges specified in Eo % [m,objn]=size(Eo); C=vec(sparse(n,n)); for k =1:objn i = Eo(1,k); j = Eo(2,k); a=sparse(zeros(n,1)); if (i ~= 1), a(i-1)=1;end; if (j ~= 1), a(j-1)=-1;end; C=C-vec(a*a'); end; % % Set up algorithm parameter to call the SDP solver sedumi % K.s=[n]; pars.eps=toler; pars.stepdif=0; [Y,y,info] = sedumi(A,b,C,K,pars); % % Prepare the output solutions % Y = mat(Y); X = chol(Y); r = rank(Y,sqrt(toler)); X =[zeros(r,1) X(1:r,:)]; S = mat(C); for k =1:mtotal S=S-y(k)*mat(A(:,k)); end; % % % end of the function