function [d,s,domset,g,threshold,varsin,varsout,varstotal,edges,density,A,cor,time,b] = cormds(threshold,myset,vars)

%This function uses complex network theory to engage a matrix of variables
%and observations, construct a graph from the adjacency matrix and finally 
%estimate the minimum dominating set.  

% Created by Sarantis Georgios and Theophilos Papadimitriou 
% 01/04/2013

%Input: 
%1. Threshold under which the nodes will NOT be connected
%2. myset is the matrix of variables under examination. The rows must contain
%the observations and the columns must contain the variables. 
%3. vars is a cell array containing the variable names

%Output:
%1. d is the number of nodes in the dominating set.
%2. domset is a cell array containing the nodes of the dominating set. The
%first column describes the vertexes as lebeled by the graph. The third
%column describes the vertexes by their initial number.
%3. s is a structure containinng the possible sub-networks that might occur after
%the thresholding proccess. 
%4. g is the graph of the dominating set. 
%5. threshold returns the selected threshold. 
%6. varsin are the variables which remain for the calculation of the
%dominating set. 
%7. varsout are the variables which are excluded from the dominating set
%process since they are isolated (zero-neighbors) and have to be otherwise studied.
%8. varstotal is a cell array which contains the whole variable set, numbered for ease
%9. A is the Adjacency matrix calculated using Pearson cross-correlations

%The function first computes the various networks derived after the
%thresholding. Afterwards the individual nodes are removed and the 
%dominating set is computed for the remaining nodes.

%Calculation of the cross-correlations
b=corrcoef(myset);
cor=b;

%delete values under a certain threshold
b(b <= threshold) = 0;
b(b > threshold) = 1;

A=b;

%Find subnetworks (start)
a=b;
N = length(b);
end_of_net = 1;
current_start = 1;
current_net = 1;
sub_nets = zeros(N,1);
total_nodes = 0;

while (end_of_net)
    all_points_found = 1;
    current_nodes = find(a(current_start,:));
    while (all_points_found)
        card_nodes = length(current_nodes);
        % I'd like to find an elegant way to union many vectors at once
        for i=1:card_nodes
            current_nodes = union(current_nodes,find(a(current_nodes(i),:)));
        end
        if card_nodes == length(current_nodes)
            all_points_found = 0;
        end
    end
    sub_nets(current_nodes) = current_net;
    total_nodes = total_nodes + card_nodes;
    current_net = current_net+1;
    
    if total_nodes<N
        %reinitialization
        current_start = min(find(sub_nets<1));
    else
        end_of_net = 0;
    end
    
    %reinitialization
    current_start = min(find(sub_nets<1));
end
ind = max(sub_nets);

for i = 1:ind
s(i).Indices= find(sub_nets == i);
end
%Find subnetworks (end)

%assign zero values to the diagonal
b=b-eye(length(b));

%Calculate number of edges and network density
n = size(b,1);
edges=sum(b(:))/2;
density=edges/(n*(n-1)/2);

%remove individual nodes (nodes with zero neighbours-row sum=0). 
c=b;
sumb = sum(c);
coorsumb = (sumb == 0);
c(coorsumb,:) = [];
c(:,coorsumb) = [];
old_position = 1:n;
new_position = old_position(~coorsumb);
b=c;

%Create 2 cell arrays containing the deleted and the remaining variables
%for the calculation of the dominating set
k=1;j=1;
varsout=[];
for i=1:length(A)
    if coorsumb(i)== 0
        varsin{1,j}=i;
        varsin{2,j}=vars{i};
        j=j+1;
    elseif coorsumb(i)==1 
        varsout{1,k}=i;
        varsout{2,k}=vars{i};
        k=k+1;
    end
end
varsin=varsin';
varsout=varsout';

%Assign numbers to the original data set
varnum=1:length(vars);
varnum=num2cell(varnum);
varstotal=[varnum; vars];
varstotal=varstotal';

%create an empty graph
g=graph;

%ATTENTION! Prior to the next function the matgraph toolbox must be 
%dowloaded and set to path!
set_matrix(g,b);

%Identify the dominating set. 
tic;
[d,S]=dom(g);
time=toc;
set1=num2cell(S);
domset=[set1 varsin(S,2) varsin(S,1)];

%Label the nodes according to the initial variables' array
lab=varsin(:,2)';
for i=1:length(lab)
    label(g,i,lab(i));
end

%Draw the network. The dominating set is shown with different colour.
notS = setdiff(1:nv(g), S);
p = partition({S,notS});
clf;
ldraw(g)
hold on
cdraw(g,p)

end