clc
%clear

    %% CASE 

% Matching: Simple Preferential attachment pj*pk

% Updating LINEAR: u_j = <c> - c_j
    % where c = clustering, and <.> is the global average

% Initial Distribution Actions: p ~ (0,1)
% Initial Distribution degree pro/conteo commitment: z ~ (-1,1)
%%  CASE 
% Initialization

niter = 50;  % number of iterations
% in = 0.1; fn = 0.9; st = 0.1; %thresholds
% s = in:st:fn;
% case_st = length(s); %number of thresholds
% T = 30; %period = 1:T;
% n = 1000; % num of people local/global
% p_min = 1e-3; %minimum probab

% DISTRIBUTIONS
p_U_N_G = [1 0 0]; %select distribution probability
z_U_N_G = [1 0 0]; %select distribution individual motivation
G2 = 1; %1=active only fo zG2
   % UNIFORM: (b-a)* rand(n,1) + a;      with range [a,b]
   % GAUSS:   normrnd(0.5,0.15,[n,1]);   a(a<0) = p_min; a(a >1) = 1-p_min;
   % GAMMA:   a = gamrnd(2,2,[n,1]);     a = a/max(a);

% NOTE: both p and z are drawn from one of that distributions so we have 12
% alternative scenarios [pU, pN, pG]*[zU, zN, zG, zG2]
% where zG is skewed toward -1 while zG2 is skewed toward 1
pU_zG2 = zeros(n, T, niter, case_st); %update by (ac- ic): reduce pAE of agent with clustering > avg clust
pU_zG2_T = zeros(n, niter,case_st); %probability when converge (time = T)

%% iterations
for cs = 1:length(s)% threshold
    for k = 1:niter        %initial condition for each iteration
         
        %% initial  Distribution
        if sum(p_U_N_G) > 1 || sum(z_U_N_G) > 1
            disp('error distrib');
           break
        end
        
   % probability
        if p_U_N_G(1) == 1 %Uniform db = 'pU';
            a = rand(n,1); 
                a(a==0) = p_min;
            
        elseif p_U_N_G(2) == 1 %Normal db = 'pN';
            a = normrnd(0.5,0.15,[n,1]);   
                a(a<0) = p_min; a(a >1) = 1-p_min;
            
        elseif p_U_N_G(3) == 1 %Gamma db = 'pG';
            a = gamrnd(2,2,[n,1]);         
                a = a/max(a);            
        end
   pU_zG2(:, 1, k, cs) = a; % initial probab of Agency p_j
        
   % propensity z 
        if z_U_N_G(1) == 1 %Uniform db = 'zU';
            z = 2*rand(n,1) - 1;                  
            
        elseif z_U_N_G(2) == 1 %Normal  db = 'zN';
            z = normrnd(0,0.33,[n,1]);   
                z(z< -1) = -1; z(z >1) = 1;                hist(z)

           
        elseif z_U_N_G(3) == 1 %Gamma db = 'zG';
            z = gamrnd(2,2,[n,1]);         
                z = 2*z/(max(z))-1;               
                %histogram(z)
                if G2 == 1 %    db = 'zG2';
                    z = -z;                 
                end
            
        end
   pU_zG2(:, 1, k, cs) = a; % initial probab of Agency p_j    
        
   % Initialising Adjency matrix
        Adjency = zeros(n,n, T);
        
        %% time steps
        for t = 1:T
            if t > 1 % update prob only from period 2
                clear cc a1 a2 a3 gc icg ac ica deg
                [~,ac,ica] = clust_coeff(Adjency(:,:, t-1), 2); %flag = 2 undirected
               
            % Updating probab        
            p = pU_zG2(:,t-1, k,cs);         
            uc =  exp(z.*(ac - ica) ); %if individual >(<) average then uc <(>) 0
                    if isnan(uc) %== 1                        
                        uc = 1; 
                    end               
                aa = p.*uc ; % update probab
                aa(aa>1) = 1; aa(aa<0) = p_min;
            pU_zG2(:, t, k,cs) = aa;
            clear p 
            
                if min(uc) == 1%abs(min(uc)) <1e-5 && abs(max(uc)) <1e-5 % stop updating
                    pU_zG2_T(:, k,cs) = pU_zG2(:, t-1, k,cs);
                    % filling remaing years 
                    pU_zG2(:, t+1:T, k,cs) = repmat(aa,1, T-t);  
                    break
                end
            clear aa p aad u uc
            end %terminate if-t
            
            %% Matching: preferential attachment with threshold s
                     adj = pU_zG2(:,t, k,cs)*pU_zG2(:,t, k,cs)'; %col_array*row_array
                     adj(adj >= s(cs) ) = 1; %create link
                     adj(adj < s(cs) ) = 0;%not create a link
                     adj = adj - diag(diag(adj)); %zero on diagonal
            Adjency(:,:,t) = adj; %adjency matrix   
            clear u uc adj
            if t == T
                pU_zG2_T(:, k, cs) = pU_zG2(:, T, k,cs);
            end
        end % time
    end % k niter
end %cases s threshold

clear A k cc r c t a1 a2 a3 gc ac as ic row col a aa aad cs Adjency
clear deg ica in fn st u uc z q zz G2 fc ccc db nome pc_f pc pcc 
%nome = char(strcat('Case_z0_n',num2str(n),'_iter_',num2str(niter)) );
%save(nome)