function GrauLeguiaPRE2019example(N,resets,seed,epsilon,rho,noise,homo,varargin)
% Inferring directed networks using a rank-based connectivity measure.
% PRE, 99.1 (2019) 012319.
% For further information, please refer to it, and to GrauLeguiaPRE2019Readme.pdf. 
%How you call it
% N: number of nodes
% resets: number of dynamical resets which we average.
% seed: Seed for the generation of Adjacency matrix
% epsilon: coupling strength
% rho: link density of the Adjacency matrix
%Homo: if 1 homogeneous Lorenz with b=28, if 0 heterogeneous Lorenz with b\in[28,48]if nargin < 1
%   load dataPaper
  % Parameters as in the paper [1]:
  
if nargin==0
  N = 16; 
  rho = 0.1;  
  resets = 1; 
  seed = 1; 
  epsilon =0.1;
  noise=0;
  homo=1;
end

%nš of resets
llavor=resets;
%seed value 
llavor2=seed;




%Matrix we will save of L, Croscorrelation,stand deviation of
%signals

Lmeantot=zeros(N,N,llavor);

%Llavor is number of resets we want to make for the inference 

for seed=1:llavor
    
% for all nš of llavors resets, we keep A the same.
rng(llavor2);

pd = makedist('Normal','mu',epsilon,'sigma',0);
problink=rho;

%Random Adj matrix

A=rand(N,N);
A=A<problink;
A(1:N+1:end)=0;
r = random(pd,N*N,1);
A=double(A);
nlinks=1;
%define adjacency matrix with weigths
eps=zeros(N,N);
for row=1:N
    for colum=1:N
        if A(row,colum) == 1
            eps(row,colum)=r(nlinks);
            nlinks=nlinks+1;
        end
    end
end
%parameters and integration of the Lorentz with a given network
d=3*N;

rng(llavor2);
%Lorenz rho coefficient
if homo==0
R=28 + (48-28).*rand(N,1);
elseif homo==1
R=28.*ones(N,1);
end
rng(seed);
%Initial conditions
x0=rand(1,3*N);
timestep = 0.005;
signallength = 40576;

[signals, times] = IntegratorNetnoise(0,@LorenzLorenzODENetdifR,signallength,timestep,x0,R,eps,noise);

ds = 6;
signals = signals(16001:end,:);%discart transients
signals = signals(1:ds:end,:); %downsample
contai=0;

% select only the x's signals from the x,y,z signals of Lorentz
datax=zeros(length(signals),N);
for i=0:3:3*N-1
    contai=contai+1;
    datax(:,contai) = signals(:,1+i);
end
%Compute the connectivity of the time series
[Lmatrix]= LMultiBig(datax(:,1:N),5,5,5,15);







%saves
Lmeantot(:,:,seed)=Lmatrix;

end
% mean of L for all dynamical resets
Lmeantotfinal=mean(Lmeantot,3);
Av=reshape(A,[],1);
Arv=reshape(Lmeantotfinal,[],1);
Av(1:N+1:end)=[];
Arv(1:N+1:end)=[];    
%Compute ROC curve and Area under ROC curve for each seed
[xpair,ypair,tpair,Aurocpair]=perfcurve(Av,Arv,'1');

%Possible plots
%Plot Matrix of L mean values
% figure;imagesc(Lmeantotfinal);
% %Plot the adjacency matrix
% figure;imagesc(A);
% %plot some signals
% figure;plot3(signals(:,1),signals(:,2),signals(:,3));
% %plot ROC curve
% figure;plot(xpair,ypair,0:0.1:1,0:0.1:1);

end