octave-tisean 0.2.3-3 / usr / share / octave / packages / tisean-0.2.3 / c1.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{output} =} c1 (@var{S})
## @deftypefnx{Function File} {@var{output} =} c1 (@var{S}, @var{paramName}, @var{paramValue}, @dots{})
##
## Computers curves for the fixed mass computation of information dimension
## (mentioned in TISEAN 3.0.1 documentation).
##
## A logarithmic range of masses between 1/N and 1 is realised by varying the
## neighbour order k as well as the subsequence length n. For a given mass
## k/n, n is chosen as small is possible as long as k is not smaller than the
## value specified by parameter @var{k} .
##
## You will probably use the auxiliary functions c2d or c2t to process the
## output further. The formula used for the Gaussian kernel correlation sum
## does not apply to the information dimension. 
##
## @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.
## @end table
##
## @strong{Parameters}
##
## @table @var
## @item mindim
## The minimum embedding dimension [default = 1].
## @item maxdim
## The maximum embedding dimension [default = 10].
## @item d
## The delay used [default = 1].
## @item t
## Minimum time separation [default = 0].
## @item n
## The number of reference points. That number of points are selected at
## random from all time indices [default = 100].
## @item res
## Resolution, values per octave [default = 2].
## @item i
## Seed for the random numbers [use default seed].
## @item k
## Maximum number of neighbors [default = 100].
## @end table
##
## @strong{Switch}
##
## @table @var
## @item verbose
## Display information about current mass during execution.
## @end table
##
## @strong{Output}
##
## The output is a @var{maxdim} - @var{mindim} + 1 x 1 struct array with the
## following fields:
## @table @var
## @item dim
## The embedding dimension of the struct.
## @item c1
## A matrix with two collumns that contain the following data:
## @enumerate
## @item
## radius
## @item
## 'mass'
## @end enumerate
## @end table
##
## @seealso{demo c1, c2d, c2t}
##
## @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 c1 of TISEAN 3.0.1 
## https://github.com/heggus/Tisean"

function output = c1 (S, varargin)

  if (nargin < 1)
    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
  mindim     = 1;
  maxdim     = 10;
  delay      = 1;
  tmin       = 0;
  cmin       = 100;
  resolution = 2;
  seed       = 0;
  kmax       = 100;

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

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

  p.addParamValue ("mindim", mindim, isPositiveIntScalar);
  p.addParamValue ("maxdim", maxdim, isPositiveIntScalar);
  p.addParamValue ("d", delay, isPositiveIntScalar);
  p.addParamValue ("t", tmin, isNonNegativeIntScalar);
  p.addParamValue ("n", cmin, isPositiveIntScalar);
  p.addParamValue ("res", resolution, isPositiveIntScalar);
  p.addParamValue ("i", seed, isNonNegativeScalar);
  p.addParamValue ("k", kmax, isPositiveIntScalar);
  p.addSwitch ("verbose");

  p.parse (varargin{:});

  # Assign input
  mindim     = p.Results.mindim;
  maxdim     = p.Results.maxdim;
  delay      = p.Results.d;
  tmin       = p.Results.t;
  cmin       = p.Results.n; 
  resolution = p.Results.res;
  seed       = p.Results.i;
  kmax       = p.Results.k;
  verbose    = p.Results.verbose;

  if (mindim > maxdim)
    warning ("Octave:tisean", ["Parameter 'mindim' is greater than ", ...
                               "'maxdim', setting 'mindim' = 'maxdim'"]);
    mindim = maxdim;
  endif

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

  output = __c1__ (S, mindim, maxdim, delay, tmin, cmin, resolution, seed,
                   kmax, verbose);

endfunction

%!demo
%! res   = c1 (henon (5000)(:,1), 'd', 1, 'maxdim', 6, 't',50, 'n', 500);
%! slope = c2d (res, 2);
%!
%! do_plot_c1  = @(x) semilogx (x{1}(:,1),x{1}(:,2),'g');
%! hold on
%! arrayfun (do_plot_c1, {slope.d});
%! plot ([5e-4 1],[1.2 1.2])
%! hold off
%! axis tight
%! ylim ([0 3]);
%! xlabel ("Epsilon")
%! ylabel ("Local slopes");
%! title ("Information dimension")
%!###############################################################


%% testing if it works with default parameters
%!test
%% res_c1 was generated by TISEAN 3.0.1
%! res_c1 =  [0.00204866193 0.000562021392;0.00458388403 0.00152773259;0.00732116727 0.00251880521;0.0114847319 0.00351527636;0.0178178754 0.00551304966;0.0229953099 0.00751305372;0.0337692834 0.0115151331;0.0432167873 0.0155181978;0.070036374 0.0225279164;0.0809115991 0.0315354168;0.105735436 0.045548249;0.137048692 0.0635655001;0.198512435 0.0905919373;0.249011219 0.127565667;0.336370319 0.180583045;0.417932093 0.255131334;0.533079803 0.361822665;0.760820627 0.510262668;1.17550445 0.726286411;0.00397582119 0.000562021392;0.00887203496 0.00152773259;0.0129663227 0.00251880521;0.0233632866 0.00451370422;0.0256882776 0.00551304966;0.0377443098 0.00851340499;0.0444159135 0.0115151331;0.0625304207 0.0165190343;0.0792045668 0.0225279164;0.112165064 0.0325363018;0.142160922 0.045548249;0.183529988 0.0645664856;0.247998312 0.0905919373;0.365789354 0.128554597;0.479705334 0.181571603;0.648174882 0.25644654;0.828185916 0.361822665;1.0652355 0.512893081;1.51780188 0.726286411];
%! clear __c1__
%! res = c1 (henon (1000), 'maxdim', 2);
%! res_mat = cell2mat({res.c1}.');
%% row 23 and 25 are excluded because TISEAN data was calculated using floats
%% this program uses doubles
%! good_idx = [1:22,24,26:38];
%! assert (res_mat(good_idx,:), res_c1(good_idx,:),-1e-5);
%% bad_idx are used as the idx of those that were further apart than the rest
%! bad_idx = setdiff (1:length(res_c1),good_idx);
%! assert (res_mat(bad_idx,:), res_c1(bad_idx,:),6e-3);

%% testing if works with other-than-default parameters
%!test
%! res_c1 = [0.0343293324 0.000628733949;0.0685995296 0.00170907611;0.104914345 0.00281779026;0.147767007 0.00393254356;0.187380955 0.00616745465;0.242883027 0.00840486214;0.320845127 0.0117625576;0.424331248 0.0162406135;0.5311203 0.0229629427;0.638983071 0.0319196843;0.77487731 0.0453561395;0.95007056 0.0643920898;1.13252032 0.0912670717;1.33960843 0.129390419;1.51084828 0.182906702;1.65632439 0.259117812;1.78787661 0.365813404;1.93386149 0.518235624;2.13333058 0.731626809;0.0578132086 0.000630145252;0.106152184 0.00171291234;0.163010269 0.0028241151;0.214579925 0.00394137064;0.280628532 0.00618129829;0.353775769 0.00842372701;0.448207438 0.0117889596;0.582583368 0.0162770674;0.693925321 0.0230144858;0.805182397 0.0319913328;0.948441863 0.0454579443;1.12319851 0.0645366237;1.28063464 0.0914719254;1.44319057 0.129558906;1.57120717 0.183243573;1.71812892 0.259117812;1.83539987 0.367163211;1.98069465 0.518235624;2.16835737 0.737046182;0.0653617978 0.000630145252;0.120350979 0.00171291234;0.182360902 0.0028241151;0.22783263 0.00394137064;0.300161123 0.00618129829;0.384833038 0.00842372701;0.497192621 0.0117889596;0.616096795 0.0162770674;0.750326335 0.0230144858;0.908238411 0.0331135131;1.03522205 0.0465802066;1.23226273 0.0656589493;1.37276983 0.0925942659;1.53525198 0.130750626;1.65205204 0.184261546;1.77773058 0.260474414;1.88386154 0.368523091;2.01883292 0.520948827;2.19400668 0.737046182];
%! clear __c1__
%! res = c1 (henon (1000), 'mindim', 8, 'd',2,'i',0.5,'t',50,'n',500);
%% row 1, 2, 22 are excluded because TISEAN data was calculated using floats
%% this program uses doubles
%! good_idx = [3:21,23:57];
%! assert (cell2mat({res.c1}.')(good_idx,:), res_c1(good_idx,:), 1e-5);
%% bad_idx are used as the idx of those that were further apart than the rest
%! bad_idx = setdiff (1:length(res_c1),good_idx);
%! assert (cell2mat({res.c1}.')(bad_idx,:), res_c1(bad_idx,:),-3e-3);

%% Test input validation
%% Promote warnings to error to not execute program
%!error <greater> warning("error", "Octave:tisean"); c1 ([1 2 3], 'mindim', 3, 'maxdim', 2);