function [MIR_MATRIX,MIR_LINK_NORM,LINKS]=NETWORK_INFERENCE_USING_MIR(DATA0)
%{Funtion created to perform a network inference using Mutual Information Rate.
% INPUTS: DATA

%}


GRAPHIC = 1;      %GRAPHIC
if GRAPHIC==1
fid1=figure;
end
[MIR_LINK_NORM,N]=size(DATA0);
if nargin<3
   
%---------------------------Prealocation
data=DATA0;
[MIR_LINK_NORM,N]=size(data);
%-------------- Selection of the biggest amount of partitions
max_partition=20;
min_partition=10;

MIR_grid=zeros(max_partition,N*(N-1)/2);
Calc_MIR_GRID_sum=zeros(1,N*(N-1)/2);
MIR_MATRIX=zeros(N*(N-1)/2,3);
rr=0;
%--------------Calculation

rN=0;

for grid=min_partition:max_partition
rN=rN+1;
r=1;
rrr=0;
for I=1:N-1
    %grid
   % I
    for J=I+1:N
    %    J
        rrr=rrr+1;
        r=r+1;
        DATA0=zeros(MIR_LINK_NORM,2);
        DATA0(:,1)=data(:,I);
        DATA0(:,2)=data(:,J);
        [Is] = MI_from_DIFF_PARATITIONS(DATA0,grid,grid);
        [~,T]=Transitive_matrix(DATA0,grid);
         Is=Is/T;
        IsR=Is;
        %end
        if grid==max_partition
            rr=rr+1;
            MIR_MATRIX(rr,1)=rr;
            MIR_MATRIX(rr,2)=I;
            MIR_MATRIX(rr,3)=J;
        end
        LINKS(rrr,1)=I;
        LINKS(rrr,2)=J;
        LINKS(rrr,3)=rrr;

        MIR_grid(grid,1)=grid;
MIR_grid(grid,r)=IsR;
    end
end
Calc_MIR_GRID=MIR_grid(grid,2:end);
min_MIR_GRID=min(MIR_grid(grid,2:end));
max_MIR_GRID=max(MIR_grid(grid,2:end));


Calc_MIR_GRID_sum=Calc_MIR_GRID_sum+(Calc_MIR_GRID-min_MIR_GRID)./(max_MIR_GRID-min_MIR_GRID);

Calc_MIR_GRID_sum_pos=Calc_MIR_GRID_sum;



%---------Final PLOT
%%{
if GRAPHIC==1
    bar(Calc_MIR_GRID_sum_pos,'g');
    %set(gca,'XLim',[0 1+N*(N-1)/2])
    %set(gca,'fontsize',40)
    %title(['$MIR$'],'interpret','latex','fontsize',45)
    inp=sprintf('MIR\n GRID = %s',num2str(grid));
    title(inp)
    hold off
    drawnow
end
end
%}
end
IP_nodes(:,1)=MIR_MATRIX(1:end,2);
IP_nodes(:,2)=MIR_MATRIX(1:end,3);
Calc_MIR_GRID_sum=Calc_MIR_GRID_sum/max(Calc_MIR_GRID_sum);
IP_nodes(:,3)=Calc_MIR_GRID_sum;

MIR_MATRIX=zeros(size(DATA0,2));
for k=1:length(IP_nodes)
   MIR_MATRIX(IP_nodes(k,1),IP_nodes(k,2))=IP_nodes(k,3);
   MIR_MATRIX(IP_nodes(k,2),IP_nodes(k,1))=IP_nodes(k,3);
end

end

function [Is] = MI_from_DIFF_PARATITIONS(DATA0,grid_X,grid_Y)
%----------------PREALOCATIONS
[~,~,ppc] = POG_PPC(DATA0,grid_X);
Pij(ppc(:,1))=ppc(:,2)/sum(ppc(:,2)); % JOIN PROBABILITY

if length(grid_X)~=1 && length(grid_Y)~=1
    grid_X=length(grid_X)-1;
    grid_Y=length(grid_Y)-1;
end

Pj=sum(reshape(Pij,grid_X,grid_Y),1);              %ROW PROBABILITY
Pi=sum(reshape(Pij,grid_X,grid_Y),2);              %COLUMN PROBABILITY
Pi_and_j.Pi=Pi;
Pi_and_j.Pj=Pj;
        %---------------------EZE MI---------------------
        %Is(r)=sum(sum((Pij.*log(Pij./(Pj*Pi)))));
        
        Is(1)=-1*(nansum((Pi).*log(Pi))+nansum((Pj).*log(Pj))...
            -sum((nansum((Pij).*log(Pij)))));
       %  Is(r)=-1*(sum((nansum((Pij).*log(Pij))))); % JOINT ENTROPY
H.H1=-1*(nansum((Pi).*log(Pi))); 
H.H2=-1*(nansum((Pj).*log(Pj)));
H.Hj=-sum((nansum((Pij).*log(Pij))));

%% correlation
if nargout>1
C=0;
Pij=reshape(Pij,grid_X,grid_Y);
for i=1:grid_X
    for j=1:grid_Y
        if Pij(i,j)<1e-100 || Pj(j)<1e-100 
            continue
        else
        C=C+(Pij(i,j)-Pi(i)*Pj(j))/Pj(j);
        end
    end
end
end

end

function [coordinates,data,ppc] = POG_PPC(DATA0,GRID)
%{
Function POG (points on grid) creates a Grid of boxes and identifies to wich
Box correspond each point. 
-------OUTPUT-------------
coordinates is a (LENGTH(DATA0),4) matrix with:
column 1 = points number (point identification)
column 2 = coordinate i of the column where the points belong
column 3 = coordinate j of the row where the points belong
column 4 = Joint coordinate of cells.
data = the normalization of the data values into an arrange from 0 to 1
------ Variables----------
DATA0 = time series data
GRID  = Grid for Partition

%}
%---------------Normalization----------------
%in this section we are going to normalizate the values of the data serie.
%tic;

NORMALIZATION=1;


[m]=size(DATA0,1);
max_1=max(DATA0(:,1));
max_2=max(DATA0(:,2));
%
min_1=min(DATA0(:,1));
min_2=min(DATA0(:,2));
%
if NORMALIZATION == 1
    DATA0(:,1)=(DATA0(:,1)-min_1)/(max_1-min_1); 
    DATA0(:,2)=(DATA0(:,2)-min_2)/(max_2-min_2);
end
data=DATA0;
DATA0=DATA0*GRID;
%---------------End Normalization-------------
%
%---------------Box creation -----------------
%
%t=0
%e=linspace(0,1,ne+1);       % the ne+1 is used to obtain a more exact division
%of the space)
coordinates=zeros(m,4);               %vector of locations (vl)
lc=DATA0==0;
DATA0(lc)=0.1;

for i=1:m
    coordinates(i)=i;                             %point identification
    coordinates(i+m)=ceil(DATA0(i));               %coordinate i of the box
    coordinates(i+2*m)=ceil(DATA0(i+m));           %Coordinate j of the box
    coordinates(i+3*m)=coordinates(i+2*m)+(coordinates(i+m)-1)*GRID;  %Coordinate 1D of the box
end
%-----------PPB
mx=max(coordinates(:,4));
k=zeros(mx,1);
ppc=zeros(mx,2);
for i=1:m
    k(coordinates(i,4))=k(coordinates(i,4))+1;
end
c=(1:mx)';
ppc(:,1)=c;
ppc(:,2)=k;
%a=ppb(:,2)==0;
%ppb(a,:)=[];
mlg=size(ppc,1);
lg=zeros(mlg,1);
lg=(ppc(:,2)>0);
nob=sum(lg);
lgm=ppc(:,2)==0;
ppb1=ppc;
ppb1(lgm,:)=[];
mnp=mean(ppb1(:,2));
%nob=nob-2*ne;
end


function [C_A,T,GSP]=Transitive_matrix(DATA0,GRID)

[coordinates] = POG_PPC(DATA0,GRID);

C_A=full(sparse(coordinates(1:end-1,4),coordinates(2:end,4),1));

GSP=graphallshortestpaths(sparse(C_A));

l=isinf(GSP);
B=GSP;
B(l)=0;
MB=max(B,[],2);
l=MB==0;
T=mean(MB(~l));
end