octave-tisean 0.2.3-3 / usr / share / octave / packages / tisean-0.2.3 / surrogates.m

## Copyright (C) 1996-2015 Piotr Held
##
## This file is part of Octave.
##
## Octave 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.
##
## Octave 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 Octave; see the file COPYING.  If not,
## see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn{Function File} {[@var{surro_data}, @var{pars}] =} surrogates (@var{S})
## @deftypefnx{Function File} {[@var{surro_data}, @var{pars}] =} surrogates (@var{S}, @var{paramName}, @var{paramValue}, @dots{})
##
## Generates multivariate surrogate data (implements the iterative Fourier scheme).
## Surrogate data is generated from a dataset with the aim of testing whether the 
## dataset was generated by a given process (null hypothesis). The Fourier scheme 
## assumes that the dataset is the output of a Gaussian linear stochastic process.
## Surrogate data is generally used to test the null hypothesis.
##
## @strong{Input}
##
## @table @var
## @item S
## This function always assumes that each time series is along the longer 
## dimension of matrix @var{S}. It also assumes that every dimension 
## (counting along the shorter dimension) of @var{S} is considered a 
## component of the time series. It's length must be factorizable by only
## 2, 3 and 5. If not the largest submatrix that fulfills this requirement will be
## used. The function @code{endtoend} can be used to determine what is the best
## submatrix for the data and then sending only that submatrix to this program.
## Padding with zeros is @strong{not} and option.
## @end table
##
## @strong {Parameters}
##
## @table @var
## @item n
## Sets the number of surrogates to be calculated. Determines the form of
## the output (see Output section) [default = 1].
## @item i
## The maximum number of permutations. Value '0' yields random permutations or if
## switch @var{exact} is set an unrescaled FFT surrogate. Value '1' is a surrogate
## close to the result of the AAFT procedure, but not quite the same. Value '-1'
## means the program will perform iterations until there is no change between them
## [default = -1].
## @item seed
## Set the seed for the random generator [default = use default seed].
## @end table
##
## @strong {Switch}
##
## @table @var
## @item exact
## This switch makes the spectrum of the output exact rather than a distribution.
## @end table
##
## @strong{Outputs}
##
## @table @var
## @item surro_data
## If parameter @code{n == 1} then this is a matrix that holds the surrogate data.
## If parameter @code{n > 1} then it is @var{n} x 1 cell array of matrixes with
## the data. In both cases the matrixes themselves are alligned with the input.
## @item pars
## This is a matrix of size @var{n} x 2 (if the input components were column
## vectors, otherwise transposed). The first column contains the number of
## iteration it took to generate the @var{i}-th surrogate, whereas the second
## column is the relative discrepency for the @var{i}-th surrogate.
## @end table
##
## @seealso{demo surrogates, endtoend}
##
## @strong{Algorithms}
##
## The algorithms for this functions have been taken from the TISEAN package.
## @end deftypefn

## Author: Piotr Held <pjheld@gmail.com>.
## This function is based on surrogates of TISEAN 3.0.1 
## https://github.com/heggus/Tisean"

function [surro_data, pars] = surrogates (S,varargin)

  if (nargin < 1 || nargout > 2)
    print_usage;
  endif

  if ((ismatrix (S) == false) || (isreal(S) == false) || ...
       (isreal(S) == false))
    error ('Octave:invalid-input-arg', "S must be a realmatrix");
  endif

  # Default values
  nsur           = 1;
  max_iterations = -1; # means until no change
  seed           = 0;

  #### Parse the input
  p = inputParser ();
  p.FunctionName = "surrogates";

  isPositiveIntScalar   = @(x) isreal(x) && isscalar (x) ...
                               && (x > 0) && (x-round(x) == 0);
  isValidIterationValue = @(x) isPositiveIntScalar (x) ...
                               || (isscalar(x) && ((x == 0) || (x == -1)));
  isNonNegativeScalar   = @(x) isreal(x) && isscalar (x) ...
                               && ((x > 0) || (x == 0));

  p.addParamValue ("n", nsur, isPositiveIntScalar);
  p.addParamValue ("i", max_iterations, isValidIterationValue);
  p.addParamValue ("seed", seed, isNonNegativeScalar);
  p.addSwitch ("exact");

  p.parse (varargin{:});

  # Assign input
  nsur           = p.Results.n;
  max_iterations = p.Results.i;
  seed           = p.Results.seed;
  ispec          = p.Results.exact;

  # Check if nsur is not large
  if (nsur > 1000)
    warning ("Octave:tisean", ["Parameter 'n' is larger than 1000, this ", ...
                               "function migth execute a long time and take ",...
                               "up a lot of memory"])
  endif

  # Correct S to always have more rows than columns
  trnspsd = false;
  if (rows (S) < columns (S))
    S = S.';
    trnspsd = true;
  endif

  # Check if the length of 'S' can be factorized by only 2, 3 and 5
  original_length = length (S);
  while (max (factor (length (S))) > 5)
   S(end,:) = [];
  endwhile

  if (original_length > length (S))
    warning ("Octave:tisean",...
             ["The length of 'S' was not factorizable by 2, 3 and 5. ", ...
              "Using only first %d so that it's length would fulfill this ",...
              "requirement"], length (S));
  endif

  # Compute the output
  [surro_data, pars] = __surrogates__ (S, nsur, max_iterations, ispec, seed);

  # If the input was transposed allign output with it
  if (trnspsd)
    pars       = pars.';
    surro_data = cellfun (@transpose, surro_data, 'UniformOutput', false);
  endif

  # If surro_data is a single cell, then change it to be a matrix
  if (isequal (size (surro_data), [1,1]))
    surro_data = surro_data {1};
  endif

endfunction

%!demo
%% 'x' will be a stationary Gaussian linear stochastic process 
%! x = zeros (2000,1);
%! for i = 2:2000
%!   x(i) = 0.7*x(i-1) +  (-6 + sum (rand ([size(1), 12]), 3));
%! endfor
%!
%! # 'spike' is the process above measured s_n (x_n) = x_n^3.
%! spike = x.^3;
%!
%! # Plot the data
%! subplot (2,1,1)
%! plot (spike,'g');
%! axis tight
%! title ("spike")
%! subplot (2,1,2)
%! plot (surrogates(spike),'b');
%! axis tight
%! title ("surrogates")
%!###############################################################

%!shared s,p
%! [s,p] = surrogates (henon (1000).', 'i', 50, 'n', 15);
%! p = p.';

%% Check if parameter i (max iterations works properly)
%!test
%! expected(1:15) = 50;
%! assert (p(:,1), expected.')

%% Check if the relative discrepancy remains similar
%!assert (std (p(:,2)), 0, 5e-4)

%% Check if cell was properly created and transposed
%!assert(iscell (s) && isequal (size (s), [15, 1]))
%!assert(rows (s{1}) < columns (s{1}))

%% Check if shortening data is discovered
%% Warnings are promoted to errors to avoid further computation
%!error <factorizable> warning ("error", "Octave:tisean"); surrogates (henon (11));

%% Check if shortening the data works as expected
%!test
%! warning ("off")
%! expected = surrogates (henon (100), 'seed', 0.25);
%! result   = surrogates (henon (102), 'seed', 0.25);
%! assert (result, expected)

%% Check if checking parameter 'n' works as expected
%!error <long> warning ("error", "Octave:tisean"); surrogates ([1 2], 'n', 1001);