function [constants] = analwheel(indata, varargin);

%% [constants] = analwheel(indata, varargin);
%%
%%  Examples  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  [constants] = analwheel(indata);
%%  [constants] = analwheel(indata, 'tau', 17, 'recfig');
%%  [constants] = analwheel(indata, 'zero');
%%  
%%  Input and output arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%   [argID,    (string)  See below. The values which are unambiguous can 
%%    value]    (varies)  be given without the preceeding argID.
%%
%% Valid argument IDs and corresponding values. Values which are
%% unambiguous (marked with '*') can be given without the preceeding argID.
%%
%%  Embedded Coordiante Parameters  %%%%%%%%%%%%%%%%%%%%%%%%%%
%%      Parameters controlling how embedded coordinate calculations are done
%%   tauType    (*string)  type of delay coordinate estimation: auto, ami (default)
%%  'tau'       (string )  specify what value of tau to use for Poincare map
%%  'lyap'      (*string)  do Lyapunov experiment
%%
%%  Authors  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Erin Miller   (emiller@csss.santafe.edu)
%%  Nicky Grigg   (grigg@csss.santafe.edu)
%%  Aaron Clauset (clauset@csss.santafe.edu)
%%
%%  Santa Fe Institute, Complex Systems Summer School 2003
%%  date created : 06.25.03
%%  last modified: 06.26.03
%%  MATLAB file (not necessarily OCTAVE compatible)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% path(path,'/amd/poptart/root/research/aaron/m_files/');

 %% defaults
p = struct('tau'     , 20, ...             % default delay coordinate
           'd_e'     , 3,  ...             % default number of embeddings
           'lambda'  , 0,  ...             % lyapunov exponent
           'D_o'     , 0,  ...             % capacity (fractal) dimension
           'delay'   , 0,  ...             % length of transient period
           'numsteps', 1000,  ...
           'maxsteps', 0);

 
 f = struct('tauType', 'ami', ...          % type of tau-estimation heuristic
            'd_e'    , 0,     ...          % use user specified d_e
            'recPlot', 0,     ...          % display recurrence plot
            'theta'  , 0,     ...          % analyze theta rather than omega
            'tauPlot', 1,     ...          % display tau-estimation plot
            'recfig' , 0,     ...          % record all images to disk
            'movie'  , 'off', ...          % make movie (.pngs) of attractor
            'lyap'   , 'off');             % do lyapunov experiment
 
 % parse command-line parameters
i=1; 
while i<=length(varargin), 
  argok = 1; 
  if ischar(varargin{i}), 
    switch varargin{i},
		% ODE parameters
     case 'tau',           p.tau         = varargin{i+1}; i = i + 1; f.tauType = 'off';
     case 'd_e',           p.d_e         = varargin{i+1}; i = i + 1; f.d_e = 1;
     case 'theta',         f.theta       = 1;
     case {'ami','auto'},  f.tauType     = varargin{i};
     case 'recPlot',       f.recPlot     = 1;
     case 'tauPlot',       f.tauPlot     = 1;
     case 'recfig',        f.recfig      = 1;
     case 'movie',         f.movie       = varargin{i+1}; i = i + 1;
     case 'lyap',          f.lyap        = 'on';
     otherwise argok=0; 
    end
  end
  if ~argok, 
    disp(['(Analyze Waterwheel) Ignoring invalid argument #' num2str(i+1)]); 
  end
  i = i+1; 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %% Delay Coordinate Embedding Analysis
 % ignore first 10% of run; transients
 p.delay = 240;
 fprintf('\nIgnoring first %6.0f of %8.0f seconds (transients).\n',p.delay, indata.t(end));

 % choose tau offset-value
 if (~strcmp(f.tauType,'off'))      % use tau-approximation heuristic
	fprintf('Calculating appropriate delay coordinate (tau) using');
    switch f.tauType
	   case 'ami'
         fprintf(' average mutual information.\n');
         dt        = 0.1; % sampling frequency (for interpolating to uniform time step)
         numbins   = 100; % number of bins for probabilty distributions
         [omega_uniform p.tau] = estimateTauAMI(p.delay, dt, numbins, indata, f.tauPlot, f.theta);
     case 'auto',
         fprintf(' the autocorrelation function.\n');
         if (f.theta==0)
           a = 0.25*find(ceil(autocorr(indata.omega(p.delay:end)) )==0);
         else 
           a = 0.25*find(ceil(autocorr(indata.theta(p.delay:end)) )==0);
         end;
         p.tau = floor(a(1));  % arbitrary scaling constants - does poorly otherwise
    end;
 end;
 fprintf('Delay coordinate   - tau    : %6.2f\n', p.tau);
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% False Nearest Neighbors algorithm
 
 fprintf('Embedded dimension - d_E    : %1.0f\n', p.d_e);

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Unfold attractor into p.d_e dimensional space using p.tau
 n        = floor(length(omega_uniform)/p.tau)*p.tau;
 embedded = zeros(n-(p.d_e-1)*p.tau,p.d_e);
 for k=1:p.d_e
  embedded(:,k) = omega_uniform((k-1)*p.tau+1:n-(p.d_e-k)*p.tau);
 end

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% plot Poincare Sections via delay coordinate embedding
 fprintf('Drawing Poincare sections.\n');
 figure(4); clf;
 if (p.d_e >= 3)     plot3(embedded(:,1),embedded(:,2),embedded(:,3));
 elseif (p.d_e == 2) plot (embedded(:,1),embedded(:,2)); end;
 title(strcat('Delay Coordinate Embedding; d_{E} = ', int2str(p.d_e), '; \tau = ', num2str(p.tau*dt), ' s'));
 if (f.theta==0)
  xlabel('\omega_{t}');
  ylabel(strcat('\omega_{t+\tau}'));
  zlabel(strcat('\omega_{t+2\tau}'));
else
  xlabel('\theta_{t}');
  ylabel(strcat('\theta_{t+\tau}'));
  zlabel(strcat('\theta_{t+2\tau}'));
end;
  if (f.recfig>0) print -depsc -f4 malkuswheel_poincare.eps; end;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Recurrence Plot
 if (f.recPlot>0) %  not sure if this plot works
  fprintf('Drawing Recurrence Plot. This may take a while.\n');
  figure(5); clf;
  
  omega_chunk=omega_uniform(1:250);
  
  ij_range=[1:length(omega_chunk)];
  matrix=zeros(length(ij_range));
  
  for ii=ij_range(1:end-1)
      matrix(ii,:)=(omega_chunk(ii) - omega_chunk)';
  end
  
  pcolor(matrix)
  title('Recurrence Plot');
  ylabel('\omega_{j}');
  xlabel('\omega_{i}'); 
  shading flat
 
%   for i=1:size(omega_uniform(p.delay:end))
%    for j=i:size(omega_uniform(p.delay:end))
%      matrix(i,j) = omega_uniform(p.delay+i) - omega_uniform(p.delay+j);
%    end
%   end

%   ij_range=[1:length(omega_uniform(p.delay:end-p.delay))];
% 
%   size(ij_range)
%   for ii=ij_range
%      %omega_uniform(p.delay+ij_range)
%      matrix(ii,ij_range)=omega_uniform(p.delay+ii) - omega_uniform(p.delay+ij_range);
%   end
%   
%   %matrix = triu(matrix) + tril(transpose(matrix));
%   pcolor([matrix zeros(size(matrix,1)); zeros(size(matrix,1)+1)]);
%   title('Recurrence Plot');
%   ylabel('\omega_{j}');
%   xlabel('\omega_{i}'); 
%   if (f.recfig>0) print -depsc -f4 malkuswheel_recurrence.eps; end;
 end;

 %%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Lyapunov Experiment
 if (strcmp(f.lyap, 'on'));
  delta = 0.001;
  rdata_1 = malkuswheel('lyapunov', indata,     0, 't', 15, 'nofigs');
  rdata_2 = malkuswheel('lyapunov', indata, delta, 't', 15, 'nofigs');
  figure(6);

  %% plot angular velocity as a function of time
  subplot(2,1,1);
  if (f.theta==0)
   plot(rdata_1.t/60,rdata_1.omega,'b'); hold on;
   plot(rdata_2.t/60,rdata_2.omega,'r'); hold off;
   ylabel('omega (rad/s)');
  else
   plot(rdata_1.t/60,rdata_1.theta(:,1),'b'); hold on;
   plot(rdata_2.t/60,rdata_2.theta(:,1),'r'); hold off;
   ylabel('theta (rad)');
  end;
  xlabel('time (minutes)');
%  axis([0 rdata_1.t(end)/60 1.1*min(rdata_1.omega) 1.1*max(rdata_1.omega)]);
  legend('x_o','x_o+\Delta');

  %% plot last 2 minutes of omega-vs-time data
  if (rdata_1.t(end)/60>=3)
   subplot(2,1,2);
   a = find(rdata_1.t<180);
   if (f.theta==0)
    plot( rdata_1.t(1:a(end)),rdata_1.omega(1:a(end)),'b'); hold on;
    plot( rdata_2.t(1:a(end)),rdata_2.omega(1:a(end)),'r'); hold off;
    ylabel('omega (rad/s)');
   else
    plot( rdata_1.t(1:a(end)),rdata_1.theta(1:a(end)),'b'); hold on;
    plot( rdata_2.t(1:a(end)),rdata_2.theta(1:a(end)),'r'); hold off;
    ylabel('theta (rad)');
   end;
    xlabel('first 3 minutes');
%   axis([0 newdata.t(a(end)) 1.1*min(newdata.omega) 1.1*max(newdata.omega)]);
  legend('x_o','x_{o+\Delta}');
  end;

  %% plot overlay poincare sections
  figure(7); clf;
  [romega_1 rtau_1] = estimateTauAMI(0, dt, numbins, rdata_1, 0, f.theta);
  [romega_2 rtau_2] = estimateTauAMI(0, dt, numbins, rdata_2, 0, f.theta);
  rn_1              = floor(length(romega_1)/rtau_1)*rtau_1;
  rn_2              = floor(length(romega_2)/rtau_2)*rtau_2;
  rembed_1          = zeros(rn_1-(p.d_e-1)*rtau_1,p.d_e);
  rembed_2          = zeros(rn_2-(p.d_e-1)*rtau_2,p.d_e);
  for k=1:p.d_e
   rembed_1(:,k) = romega_1((k-1)*rtau_1+1:rn_1-(p.d_e-k)*rtau_1);
   rembed_2(:,k) = romega_2((k-1)*rtau_2+1:rn_2-(p.d_e-k)*rtau_2);
  end
  if (p.d_e >= 3)
      plot3(rembed_1(:,1),rembed_1(:,2),rembed_1(:,3),'b'); hold on;
      plot3(rembed_2(:,1),rembed_2(:,2),rembed_2(:,3),'r'); hold off;
  elseif (p.d_e == 2)
      plot (rembed_1(:,1),rembed_1(:,2),'b'); hold on;
      plot (rembed_2(:,1),rembed_2(:,2),'r'); hold off;
  end;
 
  title(strcat('Delay Coordinate Embedding; d_{E} = ', int2str(p.d_e), '; \tau_1 = ', int2str(rtau_1), '; \tau_2 = ', int2str(rtau_2)));
  if (f.theta==0)
   xlabel('\omega_{t}');
   ylabel(strcat('\omega_{t+\tau}'));
   zlabel(strcat('\omega_{t+2\tau}'));
  else
   xlabel('\theta_{t}');
   ylabel(strcat('\theta_{t+\tau}'));
   zlabel(strcat('\theta_{t+2\tau}'));
  end;
  legend('x_o','x_{o+\Delta}');

  fprintf('Lyanpunov exponent - lambda : %6.4f (not working yet)\n', p.lambda);
 end;
 fprintf('Capacity dimension - D_o    : %6.4f (not working yet)\n', p.D_o);

 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Make Movie of .png files of delta-separated trajectories
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 if ~(strcmp(f.movie,'off'))
  figure(9); clf;
  p.maxstep = 300; ofs = 120;
  step = ceil(size(indata.omega(p.delay:end-p.tau),1)/p.numsteps);
  % for i=1:step:size(indata.omega(p.delay:end-p.tau),1)
% for i=1:step:p.maxstep
  for i=1:p.maxstep
   plot(indata.omega(p.delay:end-p.tau), indata.omega(p.delay+p.tau:end),'c:'); hold on;
   switch movieFlag
   case 'points'
    plot(indata.omega(i+ofs), indata.omega(p.tau+i+ofs),'b.'); 
    plot(newdata.omega(i+ofs), newdata.omega(p.tau+i+ofs),'r.'); hold off

   case 'delay'
    plot(indata.omega(p.delay+i), indata.omega(p.delay+p.tau+i),'b.'); 
    plot(newdata.omega(p.delay+i), newdata.omega(p.delay+p.tau+i),'r.'); hold off
  
   case 'stream'
    plot(indata.omega(p.delay+i), indata.omega(p.delay+p.tau+i),'b.', 'EraseMode','background'); 
    plot(newdata.omega(p.delay+i), newdata.omega(p.delay+p.tau+i),'r.','EraseMode','background'); hold off
   end;
   axis([-2 2 -2 2]);
   set(gca,'xticklabel',[], 'yticklabel',[], 'xtick', [], 'ytick', []);
   legend('attractor','x_o','x_o+\Delta',4);
   if (i<10)
    filename = strcat('malkuswheel_attractor_000',int2str(i),'.png');
   elseif (i<100)
    filename = strcat('malkuswheel_attractor_00',int2str(i),'.png');
   elseif (i<1000)
    filename = strcat('malkuswheel_attractor_0',int2str(i),'.png');
   else
    filename = strcat('malkuswheel_attractor_',int2str(i),'.png');
   end;
   print('-dpng', '-f9', filename);
    fprintf('%s\n',filename);
   drawnow;
  end;
 end;
 
 % calculate capacity dimension (fractal dimension)
 p.D_o = 0;
 
 % calculate Lyapunov exponent
 p.lambda = 0;

 constants = p;
 fprintf('\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% algorithm for calculating the autocorrelation
function [auto] = autocorr(omega)
 if (size(omega,2)==1) omega = transpose(omega); end;
 for i=1:ceil(size(omega,2)/2)
  a = corrcoef([omega(i:end) omega(1:i-1)],omega);
  auto(i) = a(1,2);
  if (mod(i,floor(size(omega,2)/20))==0) fprintf('\t%3.0f%% done.\n',100*i/(floor(size(omega,2)/2))); end;
 end
 auto = [auto auto(ceil(size(omega,2)/2):-1:1)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% algorithm for finding average mutual information (AMI) of vector
function av_I=find_av_I(omega_uniform,T_lag,omega_bins,numbins)

 % individual probabilities
 omega_delay = [omega_uniform(1:end-T_lag), omega_uniform(T_lag+1:end)]; 
 omega_hist  = hist(omega_delay,omega_bins); 
 omega_prob  = omega_hist/sum(omega_hist(:,1),1);

 % tabulate joint distribution
 binomega   = 1+floor((numbins-1).*(omega_uniform-min(omega_uniform))/(max(omega_uniform)-min(omega_uniform)));
 joint_dist = zeros(numbins);                  % initialize empty joint distribution matrix
 for ii=1:size(binomega,1)-T_lag               % count co-occurrences
     joint_dist(binomega(ii),binomega(ii+T_lag)) = joint_dist(binomega(ii),binomega(ii+T_lag)) +1;
 end
 joint_dist=joint_dist/(sum(sum(joint_dist))); % normalize joint distribution

 % calculate the AMI
 warning off
 prob_ratio = joint_dist./(omega_prob(:,1)*omega_prob(:,2)');
 warning on
 prob_ratio(find(joint_dist==0)) = 1;
 prob_ratio(isinf(prob_ratio))   = 1;
 mutual_information              = log2(prob_ratio);

 % return result
 av_I=sum(sum(joint_dist.*mutual_information));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% AMI wrapper algorithm
function [omega_uniform, tau] = estimateTauAMI(delay, dt, numbins, data, tauPlot, thetaFlag)
 t_uniform = [delay:dt:max(data.t)]'; %uniform time vector for interpolation

 % spline interpolate omega to uniform distribution of timesteps; set-up bins for AMI
 if (thetaFlag==0)
   omega_uniform = interp1(data.t, data.omega, t_uniform, 'spline');
 else
%   omega_uniform = interp1(data.t, data.theta(:,1), t_uniform, 'spline');
   omega_uniform = [diff(interp1(data.t, sin(data.theta(:,1)), t_uniform, 'spline'))/dt;0];
   
 end;
 omega_bins    = linspace(min(omega_uniform), max(omega_uniform), numbins)';

 % initialize AMI minimum search; then calculate AMIs until minimum found (greedy)
 av_I(1) = find_av_I(omega_uniform,1, omega_bins, numbins);
 av_I(2) = find_av_I(omega_uniform,2, omega_bins, numbins);
 t       = 2;
 while (av_I(t-1) > av_I(t))
    t       = t+1;
    av_I(t) = find_av_I(omega_uniform, t,omega_bins, numbins);
 end;
 % stopped loop when first minimum found
 if (tauPlot > 0) figure(10); plot([1:t],av_I,'o'); end;
 tau = t-1;

