function CUTOFF = ot_cutoff(RUNS,plt,tt)
%
%function CUTOFF = ot_cutoff(RUNS,plt,tt)
%
% Overturn Toolbox
% David S. Ciochetto
% Dalhousie University Physical Oceanography
% 2005 Mar 24 (Last Modified)
%
% This function is to calculate the recommended runlength cutoff value
% given in Galbraith and Kelley (1996) (GK96). Runs which are <= cutoff
% are argued to be indistinguishable from random noise. (This work has
% shown that this statistical property applies to uniformly distributed
% random noise and not to normally distributed noise. Noise with normal
% distributions can generate random runlengths greater than 2^-n. This
% is a caution only as it is suspected that the xi test will remove the
% suspect overturns with runlengths above this cutoff.)
% 
% This function takes a vector of runlengths as input and returns the
% cutoff value recommended by GK96. The vector of runlengths may come from
% a single or multiple profiles that have been collected. Runlengths can
% be calculated in the finction runlength.m If that is the case, send the
% output from RL of runlength.m as input RUNS.
% 
% To calculate the cutoff, this function compares the histogram of the
% runs sent to it to the theoretical probability density function (PDF)
% outlined in Galbraith and Kelley (1996)... "For a random, un-correlated
% time series with equal numbers of positive and negative values, the PDF
% of run length is P(n) = 2^-n; where P(n) is the probability of a run of
% length n." They recommend setting a cutoff runlength value as the first
% runlength which has greater probability than twice the theoretical PDF.
%
% INPUTS:
%   RUNS  =  the runlength data. May be a vector or matrix. The runs do
%        not need ot be sorted or ordered in any manner. Runlengths from
%        all profiles included in this input data will be used.
%   plt (optional) = a string indicating the users desire to
%        generate a plot of the results. ['n'(default)/'bar'/'semi']
%        'bar' - plots the histogram on linear axes.
%        'semi - plots the histogram on a semilog axis.
%   tt (optional) = a string sent by the user to title the check plot.
%
% OUTPUT:
%   CUTOFF = The integer cutoff value below which runs are considered
%        to be indistinguishable from random noise according to Galbraith
%        and Kelley (1996).
%
% NOTES:
%
% REFERENCES:
%   Galbraith, P. S. and Kelley, D. E., 1996 "Identifying overturns in CTD
%   profiles," Journal of Atmospheric	and Oceanic Technology, 13, 688-702.
%
% USES:
%
% SEE ALSO:
%	  runlength.m,
%

% Check if data is being sent
if ~exist('RUNS')
  disp(' ')
  disp('ot_cutoff.m has been called but no data was sent.')
  disp('Try ot_cutoff(RL) where RL is a vector of runlengths,')
  disp(' which may have been calculated in runlength.m')
  error('type help ot_cutoff.m for more information.')
end

% Set default values
if ~exist('plt')
  plt = 'n';
end

if ~exist('tt')
  tt ='Analysis of Runlength Data';
end

% Check if the user sent a matrix rather than a vector and
% convert the data to a single vector if required.
if size(RUNS,2)>1
  RUNS=reshape(RUNS,prod(size(RUNS)),1);
end

% Convert to a column if not.
RUNS=RUNS(:);

% Make an axis of possible integer runlenghts
N=[1:max(RUNS)];

% Put runs into a histogram with integer bins based on integers set by N
H=hist(RUNS,N);

% Normalize the results to the total area, make it a probability from 0 - 1
y = H/sum(H);

% Theoretical random uniform noise profile is 2.^(-N)
% Recommended threshold value is 2*2.^(-N)

% To compare the profiles make the values of y with no runs Nan
y(y == 0) = NaN;

% Compare the distribution of the runs to the threshold, take the first
% value >= the threshold, runs of this length can be distinguished from
% uniform random noise from the various sources.
CUTOFF = N( min( find( (y - 2*2.^(-N)) >= 0 ) ) );

% NOTE: unique(RUNS) will give an axis of integers with run lenghts if needed.

switch plt
case 'bar',
  figure
  bar(N,y);
  set(gca,'FontSize',12);
  title(tt);
  ylabel('Probability of run of length n, P(n)  [0-1]');
  xlabel('n (Solid Line with dots is P(n) = 2^{-n}, dotted line is GK96 threshold value.)');
  hold on; 
  % Add the theoretical distribution in yellow + cutoff in red
  plot(N,2.^(-N),'y.-',N,2*2.^(-N),'r:');
  hold off;
  %L = line([CUTOFF CUTOFF], get(gca,'YLim'),'LineStyle',':');
  L = line([CUTOFF CUTOFF], [10*2.^(-CUTOFF) 0.2],'Marker','v','MarkerFaceColor',[1 1 1]*0.4);
case 'semi',
  figure
  semilogy(N,y,'kx-',N,2.^(-N),'y.-',N,2*2.^(-N),'r:');
  set(gca,'FontSize',12);
  title(tt);
  ylabel('Probability of run of length n, P(n)  [0-1]');
  xlabel('n (Solid Line with dots is P(n) = 2^{-n}, dotted line is GK96 threshold value.)');
  L = line([CUTOFF CUTOFF], [10*2.^(-CUTOFF) 0.2],'Marker','v','MarkerFaceColor',[1 1 1]*0.4);
end
