error_setting_debugpoint

  % Newman fast community finding algorithm
  % Source: "Fast algorithm for detecting community structure in networks", Mark Newman
  % Input: adjacency matrix
  % Output: group (cluster) formation over time, modularity metric for each cluster breakdown
  % Other functions used: numedges.m
  % Originally: June 6, 2007, GB, last modified May 4, 2011

  function [groups_hist,Q]=newman_comm_fast(adj)

      profile on
      n = size(adj,1);   % number of vertices
      isSymmetricGraph = issymmetric(adj);

      groups_hist = cell(1, n);
      Q = zeros(1,n);
      groups_hist_entry = cell(1, n);
      groups = cell(1, n);

      edgeSum = adj;

      for i = 1:n
          groups{i} = i; % each node starts out as one community
          groups_hist_entry = i;
      end
      groups_hist{1}{i} = groups_hist_entry;

      nedges = numedges(adj); % compute the total number of edges

      Q(1) = QfnOptimized(groups,adj, nedges);

      for s = 1:n-1 % all possible iterations
          dQ=zeros(length(groups));

          for i = 1:length(groups)
              group_i=groups{i};
              for j = i+1:length(groups)
                  group_j=groups{j};

                  dQ(i,j)=0; % default
                  % check if connected

                  connected=0;
                  if edgeSum(i, j) + edgeSum(j, i) > 0
                      connected = 1;
                  end

                  if connected && ~(isempty(group_i)) && ~(isempty(group_j))
                      dQ(i,j) = deltaQOptimized( edgeSum, i, j, nedges, isSymmetricGraph);
                  end
              end
          end

          % pick max nonzero dQ
          dQ(dQ==0)=-Inf;
          [~,minij]=max(dQ(:));

          [mini,minj] = ind2sub([size(dQ,1),size(dQ,1)],minij);   % i and j corresponding to min dQ
          %update the edge sums for the group merge.

          if isSymmetricGraph
              edgeSum(mini, mini) = edgeSum(mini, mini) + edgeSum(minj, minj) + edgeSum(minj, mini);
              for l=1:length(groups)
                  if l ~= mini && l ~= minj
                      edgeSum(mini, l) = edgeSum(mini, l) + edgeSum(minj, l);
                      edgeSum(l, mini) = edgeSum(mini, l);
                  end
              end
          else
              edgeSum(mini, mini) = edgeSum(mini, mini) + edgeSum(minj, minj) + edgeSum(mini, minj) + edgeSum(minj, mini);
              for l=1:length(groups)
                  if l ~= mini && l ~= minj
                      edgeSum(mini, l) = edgeSum(mini, l) + edgeSum(minj, l);
                      edgeSum(l, mini) = edgeSum(l, mini) + edgeSum(l, minj);
                  end
              end
          end
          edgeSum(length(groups), mini) = edgeSum(length(groups), mini) + edgeSum(length(groups), minj);
          edgeSum(mini, length(groups)) = edgeSum(mini, length(groups)) + edgeSum(minj, length(groups));

          edgeSum(minj,:) = edgeSum(length(groups),:);
          edgeSum(:,minj) = edgeSum(:, length(groups));
          edgeSum(minj, minj) = edgeSum(length(groups), length(groups));

          %reduce the size of array matrix
          edgeSum = edgeSum(1:length(groups)-1, 1:length(groups)-1);

          %merge the groups.
          groups{mini} = sort([groups{mini} groups{minj}]);       % merge i and j into ith location.
          groups{minj} = groups{length(groups)};                  % swap minj with last group
          groups = groups(1:length(groups)-1);

          groups_hist{s+1} = groups; % save current snapshot
          Q(s+1) = Qfn(groups,adj);

      end % end of all iterations

      num_modules = zeros(1, length(groups_hist));
      for g=1:length(groups_hist)
          num_modules(g) = length(groups_hist{g});
      end

      profile viewer

      set(gcf,'Color',[1 1 1],'Colormap',jet);
      plot(num_modules,Q,'ko-')
      xlabel('number of modules / communities')
      ylabel('modularity metric, Q')
  end

  function dQ=deltaQOptimized(edgeSum, i, j, nedges, isSymmetricGraph)
  % Optimized version of deltaQ. For details see deltaQ.

      e_ii = edgeSum(i, i)/nedges;

      e_jj = edgeSum(j, j)/nedges;
      if isSymmetricGraph
          crossCommunityEdge = edgeSum(j, j) + edgeSum(j, j) + edgeSum(i, j);
      else
          crossCommunityEdge = edgeSum(j, j) + edgeSum(j, j) + edgeSum(i, j) + edgeSum(j, i);
      end
      e_ij = crossCommunityEdge / nedges - e_ii - e_jj;

      a_i = sum(edgeSum(i, :))/nedges-e_ii;
      a_j = sum(edgeSum(j, :)) / nedges - e_jj;
      dQ = 2*(e_ij-a_i*a_j);
  end

  function Q = QfnOptimized(modules,adj, nedges)
  % Optimized verison of Qfn. For details see Qfn

      Q = 0;
      for m=1:length(modules)
          e_mm=numedges(adj(modules{m},modules{m}))/nedges;
          a_m=sum(sum(adj(modules{m},:)))/nedges - e_mm; % counting e_mm only once

          Q = Q + (e_mm - a_m^2);
      end
  end