relaxation.m

best routing protocol

function P = relaxation(X,d,op,stepsize)
%function P = relaxation(X,d,stepsize)
%function [P,P1] = relaxation(X,d,stepsize)
% force-based relaxation algorithm (in [Pierce95])
% embedding from a distance matrix X to a d-dimensional space
% stepsize specifies the "small distance" to move in relaxation algorithm
% MDS initialization if op==1, random if op==2

  n = size(X,1);
  P = zeros(n,d);

  if op==1	% use MDS to initialize embedding
    %[Y,e] = cmdscale(X);	% initial embedding in a p>=d-dimensional space
    %P = Y(:,1:d);	% keep the first d dimensions
    [P e] = mds(X, d);	% using mds function from Yun
  elseif op==2	% random initialization
    maxd = max(max(X));	% max distance
    for i = 1:d
      P(:,i) = maxd*(rand([n 1])-0.5);	% uniform random from [-maxd/2, maxd/2]
    end
  end
  %X1 = squareform(pdist(P));	% new distance matrix from embedding P
  X1 = euclidean_distance(P,P);	% new distance matrix from embedding P

  iter = 0;

  lasterror = Inf;
  error = sum(sum((X1 - X).^2))/2;
  while iter < 100 & error < lasterror	% there is some progress
    F = cell(d,1);
    G = (X1-X)./(X1+eye(n));	% avoid dividing by 0

    for k = 1:d
      Rrep = repmat(P(:,k)',n,1);	% row replication
      Crep = repmat(P(:,k),1,n);	% column replication
      F{k} = (Rrep - Crep).*G;	% force matrix at the k-th dimension
    end

    if stepsize <= 0
      % stepsize adaptive to total error (not realistic, just for temporary use)
      stepsize = sqrt(error)/100;
    end

    for k = 1:d
      P(:,k) = P(:,k) + stepsize*sum(F{k},2);
    end

    %X1 = squareform(pdist(P));	% new distance matrix from embedding P
    X1 = euclidean_distance(P,P);	% new distance matrix from embedding P

    lasterror = error;
    error = sum(sum((X1 - X).^2))/2;
    iter = iter + 1;
%{
%    if iter == 1
%      P1 = P;
%    end
%}
  end

%  fprintf(1,'%d iterations\n',iter);