%# P562 HW 4.3 Invariant Measures

global MU

%# ------------------------------------------------

%# Declare functions
function init()
  %# init()
  %# Sets up initial variables and settings
  global MU

  MU = 0.2; % default mu value

  hold on     % we'll keep holding until we plot something new
  gset nokey  % no legend

%  gset terminal "postscript eps color"
  gset out "/dev/null"

endfunction

function y=logistic(x,mu)
  %# y = logistic(x,[mu])
  %#
  %# Calculates the logistic map y = 4*x*mu*(1-x).
  %# If mu is not specified, the global variable "MU" is used.

  %# override the default if asked for
  if (nargin<2) global MU; mu=MU; end

  %# calculate the value of the function and return
  y = 4 .* mu .* x .* (1-x);
endfunction

function plot_logistic(mu)
  %# plot_logistic([mu])
  %#
  %# Plots y=x and y=f(x) on the same axis

  if (nargin<1) global MU; mu=MU; end
  
  %# Set up x-coordinates
  x = [0:0.001:1]';

  %# Clear graph (take off holding), plot, then reset hold
  hold off
  plot(x,x,x,logistic(x,mu));
  hold on
endfunction

function X=iterate(cycles,x0,mu)
  %# X = iterate(cycles,x0,[mu])
  %#
  %# Iterates the logistic map for a given number of cycles,
  %# using supplied initial value and mu
  %#
  %# Outputs a list of all values

  if (nargin<3) global MU; mu=MU; end

  %# Allocate an empty list
  X = zeros(cycles,1);
  X(1) = x0;

  %# Begin iteration
  for i = 2:cycles
    X(i) = logistic(X(i-1),mu);
  end
endfunction

function plot_cycles(X)
  %# plot_cycles(X)
  %#
  %# Plots a history of iteration.
  %# Best if overlayed on top of plot_logistic
  %# using the hold on / hold off commands.

  0; % stop doc comment

  %# First convert X into a format we can use.
  %# For each value in X, we will have two points
  %# on the plot - a horizontal move and a vertical
  %# move.

  %# Allocate empty vectors
  xx = yy = zeros(length(X)*2,1);

  %# Set up first few entries
  xx(1) = xx(2) = X(1);
  yy(1) = yy(2) = X(1);

  %# Now cycle through the rest of X
  for i = 2:length(X)
    xx(2*i-1) = X(i-1);
    xx(2*i) = X(i);
    yy(2*i-1) = X(i);
    yy(2*i) = X(i);
  end

  %# Finally plot
  plot(xx,yy)
endfunction

function rho = invariant_density(x,cycles,mu)
  %# rho = invariant_density(x,cycles,mu)
  %#
  %# Given mu, calculates the invariant density

  if (nargin<3) global MU; mu=MU; end

  %# First pick a "typical" point
  x0 = 0.1;

  %# Iterate x0 N_transient times
  for i=1:1000
    x0 = logistic(x0,mu);
  end

  %# Now call 'iterate' to get a list of values
  X = iterate(cycles,x0,mu);

  %# Use the built-in histogram function
  rho = hist(X,x,length(x)); % normalize to 1
endfunction

function compare_invariant_density(x,rho)
  %# compare_invariant_density(x,rho)
  %#
  %# Given the output to invariant_density,
  %# plots a histogram and the theoretical curve.

  hold off % Clear the plot

  %# Plot the bar graph
  bar(x,rho); 

  %# Overlay the theory graph
  theory = 1./(pi*sqrt(x.*(1-x))); % set up theory
  theory(find(theory>1.1*max(rho)))=1.1*max(rho); % crop theory curve
  hold on
  plot(x,theory,'b');
endfunction

function [points_m,points_x] = bifurcate(mu,cycles,transient)
  %# bifurcation(mu,cycles,transient)
  %#
  %# Takes a range mu and generates a bifurcation diagram
  %# for all the values in the range, with 'cycles' points
  %# each

  hold off % Clear plot

  %# clear output matrices
  points_x = points_m = [];
  
  for i=1:length(mu)
    m = mu(i); % stupid octave won't do for m=mu
    x0 = 0.1;
    x0 = iterate(transient,x0,m)(end);
    XX = iterate(cycles,x0,m); % better be vertical
    points_x = [points_x; XX];
    points_m = [points_m; repmat(m,size(XX))];
  end

  plot(points_m,points_x,'.');
endfunction

function save_graph(name);
  gset("out",name);
  replot;
  gset("out","/dev/null");
endfunction  

%# Run initialization function
init()


%# Make graphs (a)

global GO
if GO

  MU = 0.2;
  plot_logistic
  hold on
  plot_cycles(iterate(100,0.1));
  hold off
  title("mu = 0.2");
  xlabel("x"); ylabel("f(x)");
  save_graph("p04-3_a1.eps")
  
  MU = 0.4;
  plot_logistic
  hold on
  plot_cycles(iterate(100,0.1));
  hold off
  title("mu = 0.4");
  save_graph("p04-3_a2.eps")
  
  MU = 0.6;
  plot_logistic
  hold on
  plot_cycles(iterate(100,0.1));
  hold off
  title("mu = 0.6");
  save_graph("p04-3_a3.eps")
  
  %# Find cycle
  
  MU = 0.8;
  plot_logistic
  hold on
  x0 = iterate(200,0.1)(end);
  plot_cycles(iterate(100,x0));
  hold off
  title("mu = 0.8");
  save_graph("p04-3_a4.eps")

  MU = 0.88;
  plot_logistic
  hold on
  x0 = iterate(200,0.1)(end);
  plot_cycles(iterate(100,x0));
  hold off
  title("mu = 0.88");
  save_graph("p04-3_a4.eps")
  
  %# Show chaos
  
  MU = 1;
  plot_logistic
  hold on
  plot_cycles(iterate(100,0.1));
  hold off
  title("mu = 1.0");
  save_graph("p04-3_a5.eps");
  
  %# Invariant density (b)
  
  MU = 1;
  x = [0:0.01:1]';
  rho = invariant_density(x,250000);
  compare_invariant_density(x,rho);
  title("Invariant density, mu = 1.0");
  ylabel("rho(x)");
  save_graph("p04-3_b1.eps");
  
  %# Cusps
  
  MU = 0.9;
  x = [0:0.0025:1]';
  rho = invariant_density(x,250000);
  X = iterate(25,0.5);
  X = sort(X);
  Y = rho(floor(X*length(rho))+1); % Stupid octave doesn't have interp
  
  bar(x,rho);
  hold on
  plot(X,Y+0.5,'bx'); % plot the first 25 iterates just above the cusps
  hold off
  gset("xr",'[0.32:0.9]')
  title("Invariant density, mu = 0.9");
  save_graph("p04-3_c1.eps");

  %# Bifurcation

  mu = [0.8:0.00025:1.0]';
  cycles = 500;
  transient = 100;
  [RR,XX] = bifurcate(mu,cycles,transient);
  gset("xr",'[0.8:1.0]')
  title("Bifurcation plot");
  xlabel("mu"); ylabel("x");
  save_graph("p04-3_e1.eps");

end % (if GO)