% 
% logistic map 
% the idea is based on 
% R. C. Hilborn, "Chaos and Nonlinear Dynamics - an introduction for
% scientists and engineers", second edition, Oxford University Press, 2000.
% page 19-25.

% logistic map given as: x(n+1) = A * x(n) * ( 1 - x(n))
% defines the mathematical model to describe the growth of biological
% populations. (refer May, 1976 article)
% usually, x values is normalized between 0<= x(n) <= 1 and 0 <= A <= 4.
% In whole equation x(n+1) = f(x(n)), f() is called as iterated map
% function or just map function.

% one more concept worth noting is: any value of x which satifying the
% condition : x` = f(x`) is called as fixed point of the iterated map.
% for our example, x` = ZERO.  We could have one more value of 
% x` = 1 - (1/A) For A < 1, as x cannot be negative (it is population,
% hence), x` = ZERO is the only fixed point of iterate map
% For A > 1, x` has two values  betweeb the range 1 and 0.

% ========================================================================
% The matlab code for Chaos Theory - 
%                 Nonlinear Statistical Modeling of Speech
% date : July 14, 2005 Thursday.
% IES: http://www.cavs.msstate.edu/hse/ies/projects/nsf_nonlinear
%
% ========================================================================
% ========================================================================
%         
%       PART ONE : LOGISTIC MAP
%        Concept of 'what are Logistic Maps' [figure 1]
%        behavior and estimation to some initial value [figures 2 ands 3] 
%
% ========================================================================

% setting up the environment
% clearing all the variables, clearing the command window, closing all the
% figures (plots)
%
clear all;
clc;
close all;

% x is varied between 0 and 1
%
x = 0:0.001:1;

% The iterate gain factor is initialized to A = 0.1
%
A = 0.1;

% taking the number of elements in x, 
% to know how many times the logistic map will iterate
%
len = size(x,2);

% enter into the loop, for different values of A starting with 0.1
% 
while ( A < 4)
  
    % just keep on storing the values 
    % here in, it is important to note that we are checking for the
    % population growth plot over the years, (as every time, the iteration
    % or the logistic map runs on different value of previous year
    % population
    % 
    for i = 1:1:len
        mapValue(i) = A * x(i) * ( 1 - x(i));
    end;
    
    % plot the variation in population for the corresponding values of x
    % at a particular value of iterate gain factor
    %
    plot(x, mapValue);
    
    % setting up the information for the plot
    % setting up the limits for the axis, line thickness, labels,
    % specifying the text at a particular point
    %
    xlim([min(x) max(x)]);
    ylim([0 max(mapValue)+ 0.1]);
    set(findobj(gca,'Type','line','Color',[0 0 1]),'Color','red',...
    'LineWidth',1.2)
    xlabel('0 \leq x \leq 1');
    ylabel('f(x)');
    text(x((len+1)/2), mapValue((len+1)/2)+0.01,...
        ['A = ', num2str(A)],...
        'HorizontalAlignment','center');
    
    % you can decide the change in the iterate growth factor
    %
    A = A + 0.5
    
    % as we want to see all the plots to compare the behavior 
    % refer figure 1.10. pp.20 of the reference
    %
    hold on;
end;
title('Logistic Maps');


% Just getting the values for a seed and time evolution
% The population shifts between two points only
%
% population is normalized between 0 and 1
%
iterate = 20;
yinitial = 0.250;
A = 3.1; % picking some value between 0 and 4 (limits not included)

y(1) = yinitial;

% this is the logistic map function
% y(i) stores the ith year population based on (i-1)th population
%
for i = 2:1:iterate
    y(i) = A * y(i-1) * ( 1 - y(i-1));
end;

% plotting it out
% you will see only few values here as population settles down between two
% values after a couple of iterations
%
figure;
for i = 1:1:iterate-1
    plot(y(i),y(i+1),'--rs','LineWidth', 2, 'MarkerEdgeColor','k',...
        'MarkerFaceColor','g', 'MarkerSize',10);
    hold on;
end;

title(['logistic map - with ystart = ',...
    num2str(yinitial), ' & A = ', num2str(A)]);
xlabel('previous year population');
ylabel('this year population');

% another way of looking at it
% results similar to table 1.1 pp 23 of reference
% 
figure;
plot(y);
xlim([1 iterate]);
title([' Population Evolution (start value = ',...
    num2str(yinitial), ' & A = ', num2str(A), ')']);
xlabel('year index');
ylabel('year population');

% ========================================================================
%         
%       PART TWO : BIFURCATION DIAGRAM
%                Concept of Attractors and basin of attraction
%
% ========================================================================

% part for bifurcation diagram for logistic maps
% The idea is to understand the concept of period-doubling phenomenon
% The basic of attractors and attractors 

% for example, with any value of A < 1, the population terminates to zero.
% Hence, ZERO attracts all the values towards it, for A < 1.
% for 1 < A < 3 : the population reaches a single value and settles down
% at A = 3: we see a sudden change, the trajectory settles down between two
% values. trajectory values oscillates back and forth between two values.
% since, the population fraction returns to same value every two years,
% this is period-2 behavior. 
% Similarly, at value close to A = 3.4498, we see trajectory settles down
% between four values => period-4 behavior.. and so on.

% just for the time being
% to concentrate of this part of the code
%
bifurIterate = 200;
xstart = 0.1;
A = 0.1;
figure;

% idea is to computation by picking up a value of A, then pick up a
% starting point between 0 and 1, iterate the map function 100 times to
% allow the trajectory to appraoch the attractor values and then plotting
% the next 100 values of x. pp. 24 of reference
%
while (A < 4)
    
    % restart with some other value of population, so as to check for the
    % attractors and iterate growth factor A
    % 
    xmap = xstart + 0.1;
    A = A + 0.01;
    
    for i = 2:1:bifurIterate
        % xmap(i) = A * xmap(i-1) * ( 1 - xmap(i-1));
        xmap = A * xmap * (1 - xmap);
        
        if ( i > 100)
            plot(A, xmap, 'color', 'b');   
            hold on;
        end;
    end;
end;

title('Bifurcation Diagram for Logistic map - period doubling phenomenon');
xlabel('iterate growth factor A');
ylabel('steady state value(s) of x');
text(3.02,0.65, ['A = ', num2str(3)], 'color', 'red');

% ========================================================================
%         
%       PART TWO : CHAOTIC BEHAVIOR - 
%                 concept of 'Sensitive dependence to initial condition'
%
% ========================================================================

% the identication of a 'chaotic' behavior
% sensitive dependence to initial condition effect
% even with starting initial conditions which are very close by, 
% [0.4 and 0.4010 and 0.4005] - three conditions checked for
% after about 10th or 12th iterate, you find the difference 
% in the pupolation estimate 
% 
clear all;
clc;
close all;

A = 3.99;
sensIterate = 13;

xstart_1 = 0.4000;
xstart_2 = 0.4010;
xstart_3 = 0.4005;

xmap_1(1) = xstart_1;
xmap_2(1) = xstart_2;
xmap_3(1) = xstart_3;

for i = 2:1:sensIterate
    xmap_1(i) =  A * xmap_1(i-1) * ( 1 - xmap_1(i-1));
    xmap_2(i) =  A * xmap_2(i-1) * ( 1 - xmap_2(i-1));
    xmap_3(i) =  A * xmap_3(i-1) * ( 1 - xmap_3(i-1));
end;

figure;
plot(xmap_1, 'color', 'red');
hold on;
plot(xmap_2, 'color', 'blue');
hold on;
plot(xmap_3, 'color', 'black');
xlabel('year index');
ylabel(' population');

title(['Sensitive dependence to initial condition', ...
    '   [ red = ',  num2str(xstart_1),...
    ' blue = ', num2str(xstart_2),...
    ' black = ', num2str(xstart_3), ' ]']);

% putting the values of population at the specific points on the graph
% for 10th iteration
%
text(10, xmap_1(10),[num2str(xmap_1(10))], 'color', 'red');
text(10, xmap_2(10),[num2str(xmap_2(10))], 'color', 'blue');
text(10, xmap_3(10),[num2str(xmap_3(10))], 'color', 'black');

% for 12th iteration
%
text(12, xmap_1(12),[num2str(xmap_1(12))], 'color', 'red');
text(12, xmap_2(12),[num2str(xmap_2(12))], 'color', 'blue');
text(12, xmap_3(12),[num2str(xmap_3(12))], 'color', 'black');