% FINDSYNERGY find synergy between inputs
%    [z,s,m] = findsynergy(X,Y)
%       X : input array
%       Y : output array
%       z : columns from input array showing synergy
%       s : synergy ratio
%       m : total area
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [z,s,m] = findsynergy(X,Y)

if size(X,2) ~= size(Y,2)
   error('X and Y are not compatible.');
end

% Check for uniqueness
N = size(X,2);
len_x = size(X,1);
len_y = size(Y,1);

a = [];
z = [];
s = [];
m = [];

% Maximum order
max_order = 0;

% Subtract minimum from Y
%Y = Y - min(min(Y))*ones(size(Y));

for i = 1:N
   x = X(:,i);
   [I,J,V] = find(x~=0);
   num_inputs = length(I);
   % must have 2 or more non-zero inputs
   max_order = max(max_order,num_inputs);
   if num_inputs >= 2
      % initialize storage for inidivual x and y columns
      xi = [];
      yi = [];
      % find individual components
      for j = I'
         % construct individual vector
         xj = ones(len_x,1)*0;
         xj(j) = x(j);
         % look for it
         for k = 1:N
            if X(:,k) == xj
               xi = cat(2,xi,xj);
               yi = cat(2,yi,Y(:,k));
               break;
            end
         end
      end
      % Calculate synergy ratio
      z = cat(2,z,x);
      num = trapz(Y(:,i)+repmat(eps,length(Y(:,i)),1));
      den = sum(trapz(yi+repmat(eps,size(yi,1),size(yi,2))));
      % If areas are too small, set synergy to 0
      m = cat(1,m,num);
      s = cat(1,s,num/den);
   end
end

% Plot synergy heat mat if max_order == 2
if max_order == 2
   figure
   load('cmap_synergy.mat')
   for i = 1:len_x
      for j = (i+1):len_x
         % Pairwise synergy
         zij = z([i,j],:);
         sij = s';
         aij = a';
         [I,J,V] = find(zij==-1);
         zij(:,J) = [];
         sij(:,J) = [];
         aij(:,J) = [];
         for ii = 1:size(zij,2)
            if isnan(sij(ii))
               rectangle('Position',...
                         [i+0.2*zij(1,ii),j+0.2*zij(2,ii),0.15,0.15],...
                         'FaceColor','y');
            else
               c = ceil(sij(ii)*32);
               if c < 1
                  c = 1;
               end
               if c > 64
                  c = 64;
               end
               rectangle('Position',...
                         [i+0.2*zij(1,ii),j+0.2*zij(2,ii),0.15,0.15],...
                         'FaceColor',cmap_synergy(c,:),...
                         'EdgeColor','w');
%                     text('Position',...
%                          [i+0.2*zij(1,ii),j+0.2*zij(2,ii)],...
%                          'String',num2str(aij(ii),'%0.1f'),...
%                          'FontSize',6);
                      
            end
         end
      end
   end
end
         
% Sort by synergy
%[s,i] = sort(s,'descend');
%z = z(:,i);