## Copyright (C) 2016 Chris Adams <octave@chrisada.co.uk>
## 
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*- 
## @deftypefn {} {@var{retval} =} statelevels (@var{input1}, @var{input2})
##
## @seealso{}
## @end deftypefn

#Uses histogram method.
#1 Take a histogram
#2 find ilow and ihigh, min and max non zero hist values
#3 make two subhist
#4 take mean of each for lower and upper

function [S,H,B] = statelevels (A, varargin)
  
  #Setup arguments based on number of given arguments
  #arg1 is bins, arg2 is method
  #setup default values which can be later overwritten
  nBins = 100;
  method = 1;

  nargin = length(varargin);
  if(nargin>0)
    if(!isnumeric(varargin{1}))
      error('statelevels expects a number for argument 2');
    end
    nBins = varargin{1};
  end
  
  if(nargin>1)
    if(strcmp(varargin{2}, 'mode'))
      method = 1;
    elseif(strcmp(varargin{2}, 'mean'))
      method = 0;  
    else
      error('Unexpected argument for method');
    end
  end
    
  #Get size first for various uses
  sz = size(A);
  
  #Row column vector checking
  #Are any of the columns of size one?
  if(sz(1) != 1 && sz(2) != 1)
    error('statelevels expects a vector');
  end
  
  #Switch to row if it's a clumn vector
  if(sz(2) == 1)
    X = A';
  else
    X = A;
  end
  
  #Generate a histogram bin centres that are bounded by signal bounds
  upper = max(X);
  lower = min(X);
  
  #Generate hist vector based on the bounds in the histogram
  histVec = lower + ((1:nBins) - 0.5)' * (upper - lower) / nBins;
 
  #Get a histogram
  #Generate indexes for the histogram
  idx = ceil(nBins * (X-lower)/(upper-lower));  
  #Round zeros up to 1
  idx(idx==0) = 1;
  #Pick only indexes in the range
  idx = idx(idx>=1 & idx<=nBins);
  #Initialise the histogram with zeros
  histogram = zeros(nBins, 1);
  #For every valid index increment the value if it hits
  for i=1:size(idx,2)
    histogram(idx(i)) = histogram(idx(i)) + 1;
  end
  
  #Find the bounds, lowest non zero and highest
  for i=1:nBins
    if(histogram(i) > 0)
      iLow = i;
      break;
    end
  end
  for i=nBins:-1:1
    if(histogram(i) > 0)
      iHigh = i;
      break;
    end
  end

  #Define index pairs that define the lower and upper histogram slices.
  #We'll need these later
  lLow  = iLow;
  lHigh = iLow + floor((iHigh - iLow)/2);
  uLow  = iLow + floor((iHigh - iLow)/2);
  uHigh = iHigh;

  lHist = histogram(iLow:lHigh, 1);
  uHist = histogram(uLow:uHigh, 1);

  #What's the amplitude, size ratio (see ref)
  dy = (upper-lower) / size(histogram,1);
  
  #If mode (default)
  if (method==1)        
    [~, iMax] = max(lHist);
    [~, iMin] = max(uHist);
    S1 = lower + dy * (lLow + iMax(1) - 1.5);
    S2 = lower + dy * (uLow + iMin(1) - 1.5);
  elseif (method==0)
    S1 = lower + dy * dot((lLow:lHigh)-0.5, lHist) / sum(lHist);
    S2 = lower + dy * dot((uLow:uHigh)-0.5, uHist) / sum(uHist);
  end

  S = [S1 S2];
  
  H = histVec;
  
  B = histogram;
  
endfunction
