function [X,S]=TriangulationGdsdp(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: dsdp 5.8 by Benson and Ye % % Check the dimension of the configuration and the number of edges % [d,n]=size(P); n=n-1; AC{1,1}='SDP'; AC{1,2}=n; [m,mtotal]=size(E); toler = 1.e-6; % % Compute edge distances specified by E and set up the SDP constraints % b=[]; A=[]; for k =1:mtotal i = E(1,k); j = E(2,k); a=sparse(n,1); if (i ~= 1), a(i-1)=1;end; if (j ~= 1), a(j-1)=-1;end; A=[A dvec(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=dvec(sparse(n,n)); for k =1:objn i = Eo(1,k); j = Eo(2,k); a=sparse(n,1); if (i ~= 1), a(i-1)=1;end; if (j ~= 1), a(j-1)=-1;end; C=C-dvec(a*a'); end; % % Set up algorithm parameter to call the SDP solver dsdp % AC{1,3}=[A C]; %OPTIONS.gaptol=toler; [STAT,y,X] = dsdp(b,AC); % % Prepare the output solutions % Y = dmat(X{1}); X = chol(Y); r = rank(Y,sqrt(toler)); X =[zeros(r,1) X(1:r,:)]; S = dmat(C); for k =1:mtotal S=S-y(k)*dmat(A(:,k)); end; % % % end of the function