Commit 52b9ec70 authored by Frank C Langbein's avatar Frank C Langbein

Matlab spinnet release and cleanup.

parent 9cf74091
This diff is collapsed.
function [F,dF] = dexpma(M,dM,v)
% [F,dF] = dexpma(M,dM,v)
%
% Gradient helper for eal_static_st
N = length(M);
AM = [M zeros(N);dM M];
if exist('v','var')
PSI = expv(AM,[v;v]);
F = PSI(1:N);
dF = PSI(N+1:end);
else
PSI = expm(AM);
F = PSI(1:N,1:N);
dF = PSI(N+1:2*N,1:N);
end
end
\ No newline at end of file
This diff is collapsed.
function plot_graph(obj, type, scale,violations)
% obj . plot_graph (type, scale, violations) - Plot distances as graph
%
% Requires presence of graphviz software.
%
% type - Mode of graph construction: 'fdp' or 'neato'
% scale - Graph scale
% violations - Only print triangle violations
%
S = obj.D >= obj.eps;
D = (obj.D * S + 1) * scale;
filename = tempname ();
fp = fopen ([filename '.dot'], 'w');
fprintf (fp, 'graph %s%d {\n clusterrank=local\n', class(obj.src), obj.N);
fprintf (fp, ' %d [ color=green, style=filled ];\n', 1:obj.N);
if violations
[x,T] = obj.check_triangle_inequality ();
[A B C] = ind2sub(size(T), find (T < -obj.eps));
X = sort([A B C],2);
n = length(A);
colors = floor(varycolor(n) * 255);
for I = 1:n
a = X(I,1);
b = X(I,2);
c = X(I,3);
ok = 1;
for J = 1:I-1
if a == X(J,1) && b == X(J,2) && c == X(J,3)
ok = 0;
break;
end
end
if ok
if S(a,b)
fprintf (fp, ' %d -- %d [ len = %f, color="#%02X%02X%02X" ];\n', a, b, D(a,b), colors(I,1), colors(I,2), colors(I,3));
else
fprintf (fp, ' %d -- %d [ len = %f, color=orange ];\n', a, b, D(a,b));
end
if S(b,c)
fprintf (fp, ' %d -- %d [ len = %f, color="#%02X%02X%02X" ];\n', b, c, D(b,c), colors(I,1), colors(I,2), colors(I,3));
else
fprintf (fp, ' %d -- %d [ len = %f, color=orange ];\n', b, c, D(b,c));
end
if S(c,a)
fprintf (fp, ' %d -- %d [ len = %f, color="#%02X%02X%02X" ];\n', c, a, D(c,a), colors(I,1), colors(I,2), colors(I,3));
else
fprintf (fp, ' %d -- %d [ len = %f, color=orange ];\n', c, a, D(c,a));
end
end
end
col_eq = 'invis';
col_neq = 'invis';
else
col_eq = 'red';
col_neq = 'blue';
end
for k = 1:obj.N
for l = k+1:obj.N
if S(k,l)
fprintf (fp, ' %d -- %d [ len = %f, color=%s ];\n', k, l, D(k,l), col_neq);
else
fprintf (fp, ' %d -- %d [ len = %f, color=%s ];\n', k, l, D(k,l), col_eq);
end
end
end
fprintf (fp, '}\n');
fclose(fp);
system([type ' -T png ' filename '.dot -o ' filename '.png']);
I = imread ([filename '.png'], 'png');
clf
imshow (I);
delete ([filename '.dot']);
delete ([filename '.png']);
pause(0.01);
end
\ No newline at end of file
function [I,L,N] = BellmanFord(w,s)
% [I,L,N] = Bellmanford(w,s)
%
% This function computes the minimum cost paths from a source node s in
% a graph defined by its link cost matrix w. This function uses the
% Bellman-Ford algorithm.
%
% w(i,j) is the cost from node i to node j, w(i,i)=0, and w(i,j)=Inf if
% there is no direct link from node i to node j.
%
% L(h,v) is the minimum communication cost from source node s to
% destination node v, under the constraint that the path from source
% node s to destination node v has no more than (h-1) links. The last
% nonvanishing row of L contains the relevant minimum costs.
%
% On return, L(v) provides minimum cost from source node s to
% destination node v. I is the unidirectional incidence matrix of the
% spanning tree; I(i,j)=1 if there is a direct path from node i to
% node j in a spanning tree. Otherwise, I(i,j)=0.
%
% N is the next node; specifically, N(v) is the next node from source
% node s in the minimum path from source node s to destination node v.
%
% Initialization
Numnode = size(w,1);
I = zeros(Numnode,Numnode);
N = zeros(1,Numnode);
N(1,s) = s;
Lh = zeros(Numnode,Numnode);
Lh(1,:) = Inf;
Lh(:,s) = 0;
% Update
for k = 2:Numnode
for v = 1:Numnode
if v ~= s
[Lh(k,v),vpre] = min(Lh(k-1,:)+(w(:,v))');
if v ~= vpre
I(:,v) = 0;
I(vpre,v) = 1;
end
if vpre == s
N(1,v) = v;
else
N(1,v) = N(1,vpre);
end
end
end
if all((Lh(k,:)-Lh(k-1,:))==0)
L = Lh(k,:);
break
end
end
end
\ No newline at end of file
function [I,L,N] = Dijkstra(w,s)
% [I,L,N] = Dijkstra(w,s)
%
% This function computes the minimum cost paths from a source node s in
% a graph defined by its link cost matrix w. This function uses
% Dijkstra's algorithm.
%
% w(i,j) is the cost from node i to node j, w(i,i)=0, and w(i,j)=Inf if
% there is no direct link from node i to node j.
%
% On return, L(v) provides minimum cost from source node s to
% destination node v. I is the unidirectional incidence matrix of the
% spanning tree; I(i,j)=1 if there is a direct path from node i to node
% j in a spanning tree. Otherwise, I(i,j)=0.
%
% N is the next node; specifically, N(v) is the next node from source node s
% in the minimum path from source node s to destination node v.
%
% Initialization
Numnode = size(w,1);
V = 1:Numnode;
L = w(s,:);
I = zeros(Numnode,Numnode);
N = zeros(1,Numnode);
N(1,s) = s;
T = [s];
% Get Next Node
for k = 2:Numnode
VminusT = setdiff(V,T);
Lj = L(VminusT);
[Lx,index] = min(Lj);
vnew = VminusT(index);
% Find previous node
for i=1:Numnode
if ismember(i,T) == 1
if L(i)+w(i,vnew) == L(vnew)
vpre = i;
if i==s
N(1,vnew) = vnew;
else
N(1,vnew) = N(1,vpre);
end
break
end
end
end
% Update least-cost paths
I(vpre,vnew) = 1;
T = union(T,vnew);
for v = VminusT
if L(v) > L(vnew)+w(vnew,v)
L(v) = L(vnew)+w(vnew,v);
end
end
end
end
\ No newline at end of file
function DD = fix_triangle_violations_isoceles (D, eps, mode, show)
% DD = fix_triangle_violations_isoceles (D, E, mode, show)
%
% Fix triangle violations in distance matrix D by taking
% 'mode' (mean, max) of two largest violation distances
% iteratively,
%
% D - Distance matrix
% eps - Precision
% mode - Mode of correcting distances (@(x) mean(x), @(x) max(x))
% show - If true, verbose mode
% DD - Resulting distance matrix
%
if show ~= 0
D0 = D;
end
DD = D;
N = size(D,1);
T = zeros(N);
while 1,
% Check if distances in D fulfill triangle inequality
for k=1:N
for l=1:N
for j=1:N
T(k,l,j) = D(k,l)+D(l,j)-D(j,k);
end
end
end
x = -min(min(min(T)));
if show ~= 0
count = 0;
display(sprintf('Max triangle inequality violation: %d', max(x,0)));
end
if x < eps
if show ~= 0
display(sprintf('Maximum distance change: %f', max(max(abs(D0-D)))));
end
break
end;
[A B C] = ind2sub(size(T), find (T < -eps));
X = sort([A B C],2);
for I = 1:size(X,1)
a = X(I,1);
b = X(I,2);
c = X(I,3);
ok = 1;
for J = 1:I-1
if a == X(J,1) && b == X(J,2) && c == X(J,3)
ok = 0;
break;
end
end
if ok == 1
% [d,i] = min([D(a,b) D(b,c) D(c,a)]); % preserves smallest edge
[d,i] = min(abs([D(b,c)-D(c,a), D(a,b)-D(c,a), D(a,b)-D(b,c)])); % preserves edge with most different length
if i == 1 % a,b
avg = sqrt (mode([D(b,c),D(c,a)])^2 + D(a,b)^2/2);
DD(b,c) = avg;
DD(c,b) = avg;
DD(c,a) = avg;
DD(a,c) = avg;
elseif i == 2 % b,c
avg = sqrt (mode([D(a,b),D(c,a)])^2 + D(b,c)^2/2);
DD(a,b) = avg;
DD(b,a) = avg;
DD(c,a) = avg;
DD(a,c) = avg;
else % c,a
avg = sqrt (mode([D(a,b),D(b,c)])^2 + D(c,a)^2/2);
DD(a,b) = avg;
DD(b,a) = avg;
DD(b,c) = avg;
DD(c,b) = avg;
end
if show ~= 0
count = count + 1;
e = min([T(a,b,c), T(a,c,b),T(b,a,c), T(b,c,a),T(c,a,b), T(c,b,a)]);
display(sprintf('%3d: triangle %d %d %d / %f %f %f : %f', count, a, b, c, D(a,b), D(b,c), D(a,c), e));
display(sprintf(' => %f %f %f', DD(a,b), DD(b,c), DD(a,c)));
end
end
end
D = DD;
end
end
\ No newline at end of file
function DD = fix_triangle_violations_mincon (D, mode, prob, alg, show)
% DD = fix_triangle_violations_mincon (D, mode, sqp, show)
%
% Fix triangle violations in distance matrix D using linear programming.
%
% D - Distance matrix
% mode - Determines norm for optimisation:
% mode > 0 - L_mode norm
% mode < 0 - L_mode norm with multiplictive distortion
% mode = [N,1] - L_N norm with scaling factor
% alg - Determines algorithm for constraint optmisation:
% alg = 0 trust-region reflective
% alg = 1 interior point
% alg = 2 sqp
% prob - If true, minimise probability change, not distance change
% show - If true, verbose mode
% DD - Resulting distance matrix
%
n = size(D,1);
S = (n * n - n) / 2;
Xorig = [];
for c = 1:n-1;
Xorig = [Xorig; D(c+1:n,c)];
end
X0 = ones (S, 1);
if size(mode,2) == 2
X0 = [X0; 1];
S = S + 1;
end
C = nchoosek (1:n, 3);
cc = size(C,1);
A = zeros (cc * 3,S);
c = 0;
for t = 1:cc
% D(b,a) + D(c,b) - D(c,a) > 0
% D(b,a) - D(c,b) + D(c,a) > 0
% -D(b,a) + D(c,b) + D(c,a) > 0
A(c + [1:3],[get_ind(C(t,2),C(t,1),n) get_ind(C(t,3),C(t,2),n) get_ind(C(t,3),C(t,1),n)]) = [ -1 -1 1; -1 1 -1; 1 -1 -1];
c = c + 3;
end
B = zeros(c,1);
LB = zeros(S, 1);
UB = Inf(S, 1);
if alg == 0
alg='trust-region-reflective';
elseif alg == 1
alg = 'interior-point';
else
alg = 'sqp';
end
if show == 0
dp = 'off';
else
dp = 'Iter';
end
options = optimset('Algorithm', alg, 'Display', dp);
if prob == 0
if size(mode,2) == 2
nn = length(X0);
if mode(1) >= 0
mode = mode(1);
func = @(x)distance_error_scaling (x);
else
mode = -mode(1);
func = @(x)distance_error_scaling_mult (x);
end
elseif mode >= 0
func = @(x)distance_error (x);
else
mode = -mode;
func = @(x)distance_error_mult (x);
end
else
if size(mode,2) == 2
nn = length(X0);
if mode(1) >= 0
mode = mode(1);
func = @(x)prob_error_scaling (x);
else
mode = -mode(1);
func = @(x)prob_error_scaling_mult (x);
end
elseif mode >= 0
func = @(x)prob_error (x);
else
mode = -mode;
func = @(x)prob_error_mult (x);
end
XorigExp = exp(-Xorig);
end
[X,FVAL,EXITFLAG,OptDetails] = fmincon(func, X0, A, B, [], [], LB, UB, [], options);
if show ~= 0
display(OptDetails);
end
DD = zeros (n, n);
p = 1;
l = n-1;
for c = 1:n-1
DD(c+1:n,c) = X(p:p+l-1,1);
p = p + l;
l = l - 1;
end
DD = DD + DD';
if show ~= 0
display(sprintf('Maximum distance change: %f', max(max(abs(D-DD)))));
display(sprintf('Norm of change to transition probability matrix: %f',norm(exp(-D)-exp(-DD))));
end
function e = distance_error (x)
e = double(norm (Xorig - x, mode));
end
function e = distance_error_mult (x)
e = norm ((Xorig - x)./Xorig, mode);
end
function e = distance_error_scaling (x)
e = norm (Xorig * x(nn) - x(1:nn-1), mode);
end
function e = distance_error_scaling_mult (x)
e = norm ((Xorig * x(nn) - x(1:nn-1))./Xorig, mode);
end
function e = prob_error (x)
e = norm (XorigExp - exp(-x), mode);
end
function e = prob_error_mult (x)
e = norm ((XorigExp - exp(-x))./XorigExp, mode);
end
function e = prob_error_scaling (x)
e = norm (XorigExp * x(nn) - exp(-x(1:nn-1)), mode);
end
function e = prob_error_scaling_mult (x)
e = norm ((XorigExp * x(nn) - exp(-x(1:nn-1)))./XorigExp, mode);
end
function p = get_ind(r, c, n)
% r > c !
p = sum(n-c+1:n-1) + (r - c);
end
end
\ No newline at end of file
function DD = fix_triangle_violations_minuncon (D, E, mode, show)
% DD = fix_triangle_violations_minuncon (D, E, mode, show)
%
% Fix triangle violations in distance matrix D using optimisation.
%
% D - Distance matrix
% E - Precision
% mode - Determines norm for optimisatioN:
% mode > 0 - L_mode norm
% mode < 0 - L_mode norm with multiplictive distortion
% mode = [N,1] - L_N norm with scaling factor
% show - If true, verbose mode
% DD - Resulting distance matrix
%
n = size(D,1);
S = (n * n - n) / 2;
Xorig = [];
for c = 1:n-1;
Xorig = [Xorig; D(c+1:n,c)];
end
%X0 = ones (S, 1);
X0 = double(Xorig);
if size(mode,2) == 2
X0 = [X0; 1];
S = S + 1;
end
C = nchoosek (1:n, 3);
cc = size(C,1);
A = zeros (cc * 3,S);
c = 0;
for t = 1:cc
% D(b,a) + D(c,b) - D(c,a) > 0
% D(b,a) - D(c,b) + D(c,a) > 0
% -D(b,a) + D(c,b) + D(c,a) > 0
A(c + [1:3],[get_ind(C(t,2),C(t,1),n) get_ind(C(t,3),C(t,2),n) get_ind(C(t,3),C(t,1),n)]) = [ -1 -1 1; -1 1 -1; 1 -1 -1];
c = c + 3;
end
if show == 0
dp = 'off';
else
dp = 'Iter';
end
if size(mode,2) == 2
nn = length(X0);
if mode(1) >= 0
mode = mode(1);
func = @(x)distance_error_scaling (x);
else
mode = -mode(1);
func = @(x)distance_error_scaling_mult (x);
end
elseif mode >= 0
func = @(x)distance_error (x);
else
mode = -mode;
func = @(x)distance_error_mult (x);
end
options = optimset('Display', dp);
[X,FVAL,EXITFLAG,OptDetails] = fminunc(@(x)err_func(x), X0, options);
if show ~= 0
display(OptDetails);
end
DD = zeros (n, n);
p = 1;
l = n-1;
for c = 1:n-1;
DD(c+1:n,c) = X(p:p+l-1,1);
p = p + l;
l = l - 1;
end
DD = DD + DD';
if show~= 0
display(sprintf('Maximum distance change: %f', max(max(abs(D-DD)))));
end
function e = err_func (x)
tv = A * x;
tv(tv>=E) = 0;
e = double(func(x) + norm(tv, mode));
end