function [ best_dist ] = GraphEditDistance( g1, g2, num_missing_nodes)
%GraphEditDistance - calculate the distance between two graphs with missing
%nodes
%   Detailed explanation goes here


%simulated annealing parameters
c = 0.9;
start_temp = 1000;
t = start_temp;

round = 0;
n = size(g1,1);
m = size(g2,1);

%get the maximal size of the two graphs
N = max(n,m);

%if the number of missing nodes is not given, perform the search on all the
%nodes
if nargin < 3
    num_missing_nodes = N;
end

max_iterations = N;

init_distance = 0;

%make the two graphs the same size by padding zeros
if n < N
      %g1 is smaller, make it NxN
    g1(N,N) = 0;
    %the penalty is N-n because we added this number of nodes
    init_distance = N - n;
end
if m < N
      %g2 is smaller, make it NxN and add a penalty of N-m
    g2(N,N) = 0;
    init_distance = N - m;
end

%%%distance = sum(sum(abs(g1 - g2))) / 2;

best_dist = Inf;
for restart = 1 : N*N/4     % sigal 24.12.12 - add 1/4 factor as no improvement seen after N*N/4

      %create a matching vector between the two matrices
      %the first K nodes are known nodes and therefore they match (where     
      %K = N-num_missing_nodes
      %the vector will look like this: 
      %[1,2,3,...,K,  ***random permutation of the numbers between K+1 and N***]
    matching= [1:N-num_missing_nodes, randperm(num_missing_nodes) + N - num_missing_nodes];
    try
          %transform g2 according to the permutation, simply by switching
          %rows and columns according to the matching vector
      perm_graph = g2(matching, matching);
    catch ME1
          x = 90;
    end

     num_failures = 0;
     iter = 0;
    while num_failures < N %&& iter < max_iterations
        iter = iter+1;
        i = randi(num_missing_nodes,1) + N - num_missing_nodes;
        j = randi(num_missing_nodes,1) + N - num_missing_nodes;

        %should we swap i and j?

        row_i = perm_graph(i,:);
        row_j = perm_graph(j,:);
        
        %calculate the current distance due to rows i and j
        curr_dist =     sum(abs(g1(i,:) - row_i)) + sum(abs(g1(j,:) - row_j));
        
        temp = row_i(i);
        row_i(i) = row_i(j);
        row_i(j) = temp;
        
        temp = row_j(i);
        row_j(i) = row_j(j);
        row_j(j) = temp;
        
        %calculate the distance due to rows i and j after swapping them
        new_dist    =     sum(abs(g1(i,:) - row_j)) + sum(abs(g1(j,:) - row_i));

        %if swapping causes an improvement, or randomly according to
        %simulated annealing:
        if new_dist < curr_dist || (rand(1) < exp( - (new_dist - curr_dist) / (c*t) ) && new_dist ~= curr_dist)


            %temp = exp( - (new_dist - curr_dist) / (c*t) )
            %swap the "matching" indexes
            temp = matching(i);
            matching(i) = matching(j);
            matching(j) = temp;

            %swap rows i and j
            temp = perm_graph(i,:);
            perm_graph(i,:) = perm_graph(j,:);
            perm_graph(j,:) = temp;

            %swap columns i and j
            temp = perm_graph(:,i);
            perm_graph(:,i) = perm_graph(:,j);
            perm_graph(:,j) = temp;

            num_failures = 0;    

            round = round+1;

            t = t*c;

        else
            num_failures = num_failures + 1;
        end
    end

    distance = sum(sum(abs(g1 - perm_graph)));

    if distance < best_dist
          best_dist = distance;
          fprintf('restart %d, best_dist %d, N %d\n', full(restart), full(best_dist), full(N));
    end
    

end
%divide by 2 due to symmetry
best_dist = best_dist / 2 + init_distance;

end