clear all;
clc;
format long

a=0;                %Initial value of the parameter%
cnt=1;
num_iter=12000;		%Number of points per time series Dp%	
M= 200;			    %Maximum number of bins% 
frac=1/10;		    %Range for the number of bins%
inc= 0.000001;		%Increment used to get the next parameter value%	


while(a<=4)         %While loop. It indicates the maximum value of parameter% 
disp('==========================================================================');

%----------------------------------------------------------------------
%               Generation of time series
%----------------------------------------------------------------------


x1(1)=0.8;          %Initial value of variable x%
i=1;
b(cnt)=a;

    while(i<num_iter)   %generation of time-series%
        x1(i+1) = a* x1(i)*(1-x1(i));
       
        if((x1(i+1)<=-Inf) ||(x1(i+1)==NaN) )
  	        i=i+1;
	        break;
      	end    
         
        if((x1(i+1)>=Inf)|| (x1(i+1)==NaN))
            i=i+1;
	        break;
      	end    
       
       	i=i+1;
    end

    y1=x1;   
    s = size(x1);   
    g(cnt)=s(2);  
   
    x=x1;
    y=y1;

    
    %----------------------------------------------------------------------
    %               Step 1 of the distribution process
    %----------------------------------------------------------------------
    
    cellsiz = g(cnt)/M;
    floor_cell= ceil(cellsiz);

    B =sort(x);
    Z(1) = B(1);
    
    kount=2;
    z_kount=2;
    while(kount<M)
	           
        if((floor_cell*(kount-1))>0 && (floor_cell*(kount-1))<=g(cnt))
            Z(z_kount)=B(floor_cell*(kount-1));
            z_kount=z_kount+1;
        else
            kount=kount+1;
            z_kount=z_kount+1;
            break;
        end   
        kount=kount+1;
  	end    
   
    Z(z_kount)=B(g(cnt));

   	i=1;
    j=1;
    Z1(1)=Z(1);
    while(i<=(z_kount-1))
        if(Z(i)==Z(i+1))
            
        else
            Z1(j+1)=Z(i+1);   
            j=j+1;
        end
       	i=i+1;
   	 end
     Z1(j)=Z(z_kount);
     Z_last=Z(z_kount);
   

    %----------------------------------------------------------------------
    %               Step 2 of the distribution process
    %----------------------------------------------------------------------



  	flag=0;

  	while(flag==0)
        N2=histc(y,Z1);
    	thr = size(N2);
    	thra = thr(2);
    	N1 = N2(1:1:thra-1);
 
    	max_N1 = max(N1);
    	sizN = size(N1);
    	sizNp = sizN(2);
   		i=1;    
    	flagin=0;

     	while((i<sizN(2)) && (sizNp>1) )
            if(N1(i)<(max_N1*frac)) 
                clear Z;
            	Z =Z1;
            	clear Z1;
            	j=1;
            
            	while(j<=i)
             	    Z1(j) = Z(j);
                	j=j+1;
            	end    
            
            	while(j<sizNp)
                    Z1(j)=Z(j+1);
                    j=j+1;
            	end
            
            	Z1(j) = Z_last;
            	flagin=1;    
            	break;
            end
        	i=i+1;
    	end
    
    	if(flagin==0)
     	    flag=1;
    	end    
    end
 
   
    snde=size(Z1);
    snde1=snde(2);
    Z1(snde1)=Z_last;
    clear Z;

    
    if(snde1==3)
        
        if( (N1(2))<((N1(1))*frac) )
            Z = Z1;
            clear Z1;
            
            Z1(1)=Z(1);
            Z1(2)=Z(3);
     
        end
    end
    
    if(snde1>3)
        Z = Z1;
        clear Z1;
        Z1=Z(1:1:(snde1-2));
        Z1(snde1-1)= Z_last;
    end
    
    
    N2=histc(y,Z1);
    thr = size(N2);
    thra = thr(2);
    N1 = N2(1:1:thra-1);
  
    bins = size(N1);
    binn=bins(2);
    bin(cnt)= binn;
    binn


    %----------------------------------------------------------------------
    %               Calculation of auto mutual information if the number of
    %               bins>1.
    %----------------------------------------------------------------------
    
    
    if(bin(cnt)>1)
   	    D1 = N1/g(cnt);
        
    	sizd=size(D1);
    	sizd1= sizd(2);
    
    	county =1;
    	while(county<=sizd1)
            county=county+1;
    	end 

    	county =1;
    	while(county<=sizd1)
        	county=county+1;
    	end 
   
   		sdfg = size(Z1);
   		sdfg1 = sdfg(2);
    
   		county =1;
    	while(county<=sdfg1)
        	county=county+1;
    	end 
   

   	
    	count=1;
    	Ixyv=0;


    	sdmp = size(Z1);
    	sdmp1= sdmp(2);
         
  		while(count<g(cnt))
       		flag=0;
        	j=1;
       		found=0;
        
        	while(flag==0 && j<sdmp1)
           		if(Z1(j)<=y(count)<Z1(j+1))    
                	found=j;
                	flag=1;
            	end
            	j=j+1;

        	end
  
        	if((found>0) && (found<=sdmp1))
            	Py=D1(found);
        	else 
                found=sdfg1-1;
                Py=D1(found);
            end    
        
           	flag1=0;
            j1=1;
            found1=0;

            while(flag1==0 && j1<sdmp1)

                if(Z1(j1)<=x(count)<Z1(j1+1))   
                    found1=j1;
                    flag1=1;
                end
                j1=j1+1;
            end 

            if((found1>0) && (found1<=sdmp1))
                Px=D1(found1);
            else 
                found1=sdfg1-1 ;
                Px=D1(found1);                               
           end    
            
        
           	if((found==found1)  && (found>0) && (found1>0) && (Px>0) && (Py>0) && (found1<=sizd1) )
                Pxy = D1(found);
                Ixyv = Ixyv + abs(Pxy*(log2(Pxy/Px/Py)));
            end     
         
            count = count+1;
        end
    
        
        Ixy(cnt)= Ixyv;
        Ixyv
        cnt=cnt+1;
    
    end

    %----------------------------------------------------------------------
    %               Special case. If the number of bins is 1.
    %----------------------------------------------------------------------
    
    
    if(bin(cnt)==1)
        Ixy(cnt)=-2; %Temporary value of Ixy so that it can be replaced easily later on%
        disp('Ixyv = -2');  
    end   

    %----------------------------------------------------------------------
    
 
    clear x1 y1 p1 p x y  N N1 Z Z1 D1 D temp Ixyv  binn found found1 count ano_count Px Py Pxy flag j1 j sizr sizt s i cellsiz floor_cell B;
    a=a+inc     % getting the nest value of parameter %

end


    %----------------------------------------------------------------------
    %               Special case. If the number of bins is 1.
    %               Replacing the maximum value of Ixy for the
    %               case when number of bins came put to be zero
    %--------------------------------------------------------------------


sinfo = size(Ixy);
sinfo1= sinfo(2)
i=1;
Imax= max(Ixy)


while(i<=sinfo1)
       	if(Ixy(i)>-2)
       		Ixym(i)=Ixy(i);
   	end
    
  	if(Ixy(i)==(-2))
        		Ixym(i)=Imax;
   	end
      
    	i=i+1;
end    

plot(b,Ixym)