% Delivered to the Government in Accordance with Contract FA8650-10-C-7080
% 
%  (c) 2011 Massachusetts Institute of Technology
% Developed with the sponsorship of the Defense Advanced Research Projects Agency (DARPA).
%  
% Permission is hereby granted, free of charge, to any person obtaining a copy of this data, including any software or models in source or binary form, as well as any drawings, specifications, and documentation (collectively "the Data"), to deal in the Data without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Data, and to permit persons to whom the Data is furnished to do so, subject to the following conditions:
%  
% The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Data.
%  
% THE DATA IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS, SPONSORS, DEVELOPERS, CONTRIBUTORS, OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE DATA OR THE USE OR OTHER DEALINGS IN THE DATA.
% 



function [C,normalized_C,E,normalized_E,total_entropy,normalized_entropy,X,A,entropy_matrix,rank_A,num_edges,num_subgraphs,r,p] = dynamical_complexity(filename,conf_level,A_init);
% Function Call: 
%
% [C,normalized_C,E,normalized_E,total_entropy,normalized_entropy,X,A,entropy_matrix,rank_A,num_edges,num_subgraphs,r,p] = dynamical_complexity(filename,conf_level,A_init);
%
% This function computes the dynamical density in both, absolute and normalized forms
% dynamical complexity = E(A)*entropy where A = dependency structure matrix for functional attributes considered
%
%
% Input:- csv file containing observed values of (n) functional attributes over (m) states of the system. 
%         States could also be replaced by time-steps, depending on the application.
%
%       - confidence level for the significance test on pairwise correlations among functional attributes (default: 0.95)
%       - Intial A matrix (binary 0/1 matrix),expressing known interactions among functional responses (optional argument)
%
%
% Output: C = dynamical complexity
%         E = sume of singular values of dependency structure matrix A
%         r = correlation matrix (n*n)
%         p = p-value matrix corresponding to correlation matrix (n*n) -
%             each terms gives the p-value of corresponding correlation coeff.
%         A = interaction matrix (binary) of functional responses
%         entropy_matrix = image entropy of pairwise scatter plots of functional responses (each entry in the matrix is normalised)
%         total_entropy = sum of all elements of the entropy_matrix
%         normalized_entropy = total_entropy/(n*(n-1)) where n is the number of functional variables
%         rank_A = rank of interaction matrix A.
%         num_subgraphs = #connected sugbraphs as expressed through interaction matrix A
%         num_edges = #interactions as expressed through interaction matrix A
%
%
% Examples:
%           [C,normalized_C,E,normalized_E,total_entropy,normalized_entropy,X,A,entropy_matrix] = dynamical_complexity('3_Spool_Turbofan.csv',0.95);close all
%
%           [C,normalized_C,E,normalized_E,total_entropy,normalized_entropy,X,A,entropy_matrix] = dynamical_complexity('Turbojet.csv',0.8);close all
%
%
numeric=csvread(filename);
[m n]=size(numeric);
if(m<n)
    X_base=numeric';
    num_vars=m;
    num_data_pnts = n;
else
    X_base=numeric;
    num_vars=n;
    num_data_pnts = m;
end
%
% Should we normalize the X matrix?
a=min(X_base,[],1);b=max(X_base,[],1);  %min and max of each colum in data matrix..
for i=1:num_data_pnts
    for j=1:num_vars
        X(i,j)=(X_base(i,j)-a(j))/(b(j)-a(j));
        if(isnan(X(i,j)))
            X(i,j)=0;
        end
    end
end
%Build the dynamic connectivity matrix based on the correlation structure
%observed in the data matrix X..
if(nargin<2)
    conf_level=0.95;     %specify confidence level for no correlation..
end
%
if(nargin==3)
    [r,s]=size(A_init);
    if(r==num_vars)
        A=A_init;
    else
        disp(['Inappropriate initial matrix size...reverting to empty initial matrix']);
        A=zeros(num_vars,num_vars);
    end
else
    A=zeros(num_vars,num_vars);
end
%
[r,p]=corrcoef(X);
[i,j]=find(p<(1-conf_level));
el=[i j];
k=length(i);
%Build dynamic/behavioral connectivity matrix..
for i=1:k
    tmp=el(i,:);
    A(tmp(1),tmp(2))=1;
end
%
if(size(A,1)==2)
    A=[0 1;1 0];
end
%
rank_A=rank(A);
E=sum(svd(A));
num_edges=sum(sum(A));
num_subgraphs=num_conn_comp(A);
%
%
ent=0;
n=num_vars;
entropy_matrix=zeros(n,n);
for i=1:n-1
    for j=i+1:n
        x=X(:,i);
        y=X(:,j);
        plot(x,y,'r.');xlabel('x');ylabel('y');
        grid on
        saveas(gcf,'scatter.jpg');
        I = imread('scatter.jpg');
        a=size(I);
        b=a(1:2);
        k=min(b);
        delete('scatter.jpg');
        max_ent=log2(k);
        ent = entropy(I)/max_ent;
        entropy_matrix(i,j)=ent;
        %
        plot(y,x,'r.');xlabel('y');ylabel('x');
        grid on
        saveas(gcf,'scatter1.jpg');
        I1 = imread('scatter1.jpg');
        a1=size(I1);
        b1=a1(1:2);
        k1=min(b1);
        delete('scatter1.jpg');
        max_ent1=log2(k1);
        ent1 = entropy(I1)/max_ent1;
        entropy_matrix(j,i)=ent1;
    end
end
close all
% keyboard
%
%correction for no correlation cases...
[i,j]=find(A==0);
list=[i j];
count=size(list,1);
for i=1:count
    tmp1=list(i,1);
    tmp2=list(i,2);
    entropy_matrix(tmp1,tmp2)=0;
end
%
total_entropy=sum(sum(entropy_matrix));
normalized_entropy=total_entropy/(n*(n-1));
%
%dynamical complexity = structure*entropy..
C = E*total_entropy;
normalized_E=E/num_vars;
normalized_C=normalized_E*normalized_entropy;
%
end

