← Back to branch summary

~ubuntu-branches/debian/sid/octave-tisean/sid

« back to all changes in this revision

Viewing changes to inst/lyap_r.m

  • Committer: Package Import Robot
  • Author(s): Rafael Laboissiere
  • Date: 2017-08-14 12:53:47 UTC
  • Revision ID: package-import@ubuntu.com-20170814125347-ju5owr4dggr53a2n
Tags: upstream-0.2.3
Import upstream version 0.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
inst/lyap_r.m
 
1
## Copyright (C) 1996-2015 Piotr Held
 
2
##
 
3
## This file is part of Octave.
 
4
##
 
5
## Octave is free software; you can redistribute it and/or
 
6
## modify it under the terms of the GNU General Public
 
7
## License as published by the Free Software Foundation;
 
8
## either version 3 of the License, or (at your option) any
 
9
## later version.
 
10
##
 
11
## Octave is distributed in the hope that it will be useful,
 
12
## but WITHOUT ANY WARRANTY; without even the implied
 
13
## warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 
14
## PURPOSE.  See the GNU General Public License for more
 
15
## details.
 
16
##
 
17
## You should have received a copy of the GNU General Public
 
18
## License along with Octave; see the file COPYING.  If not,
 
19
## see <http://www.gnu.org/licenses/>.
 
20
 
 
21
## -*- texinfo -*-
 
22
## @deftypefn{Function File} {@var{output} =} lyap_r (@var{X})
 
23
## @deftypefnx{Function File} {@var{output} =} lyap_r (@var{X}, @var{paramName}, @var{paramValue}, @dots{})
 
24
##
 
25
## Estimates the largest Lyapunov exponent of a given scalar data set using
 
26
## the algorithm described by Resentein et al. on the TISEAN refernce page:
 
27
##
 
28
## http://www.mpipks-dresden.mpg.de/~tisean/Tisean_3.0.1/docs/chaospaper/citation.html
 
29
##
 
30
## @strong{Input}
 
31
##
 
32
## @table @var
 
33
## @item X
 
34
## Must be realvector. The output will be alligned with the input.
 
35
## @end table
 
36
##
 
37
## @strong{Parameters}
 
38
##
 
39
## @table @var
 
40
## @item m
 
41
## Embedding dimension to use [default = 2].
 
42
## @item d
 
43
## Delay used [default = 1].
 
44
## @item t
 
45
## Window around the reference point which should be omitted [default = 0].
 
46
## @item r
 
47
## Minimum length scale for the neighborhood search [default = 1e-3].
 
48
## @item s
 
49
## Number of iterations in time [default = 10].
 
50
## @end table
 
51
##
 
52
## @strong{Switch}
 
53
##
 
54
## @table @var
 
55
## @item verbose
 
56
## Gives information about the current epsilon while performing computation.
 
57
## @end table
 
58
##
 
59
## @strong{Output}
 
60
##
 
61
## Alligned with input. If input was a column vector than output contains two
 
62
## columns. The first contains the iteration number and
 
63
## the second contains the logarithm of the stretching factor for that
 
64
## iteration.
 
65
##
 
66
## @seealso{demo lyap_r, lyap_k, lyap_spec}
 
67
##
 
68
## @strong{Algorithms}
 
69
##
 
70
## The algorithms for this functions have been taken from the TISEAN package.
 
71
## @end deftypefn
 
72
 
 
73
## Author: Piotr Held <pjheld@gmail.com>.
 
74
## This function is based on lyap_r of TISEAN 3.0.1 
 
75
## https://github.com/heggus/Tisean"
 
76
 
 
77
function output = lyap_r (X, varargin)
 
78
 
 
79
  # Initial input validation
 
80
  if (nargin < 1)
 
81
    print_usage;
 
82
  endif
 
83
 
 
84
  if ((isvector (X) == false) || (isreal(X) == false))
 
85
    error ('Octave:invalid-input-arg', "X must be a realvector");
 
86
  endif
 
87
 
 
88
  # Default parameters
 
89
  embdim  = 2;
 
90
  delay   = 1;
 
91
  mindist = 0;
 
92
  eps0    = 1e-3;
 
93
  steps   = 10;
 
94
 
 
95
  #### Parse the input
 
96
  p = inputParser ();
 
97
  p.FunctionName = "lyap_r";
 
98
 
 
99
  isPositiveIntScalar    = @(x) isreal(x) && isscalar (x) && ...
 
100
                             (x > 0) && (x-round(x) == 0);
 
101
  isNonNegativeIntScalar = @(x) isPositiveIntScalar (x) || (length (x) == 1 ...
 
102
                                                            &&(x == 0));
 
103
  isPositiveScalar       = @(x) isreal(x) && isscalar (x) && (x > 0);
 
104
 
 
105
  p.addParamValue ("m", embdim, isPositiveIntScalar);
 
106
  p.addParamValue ("d", delay, isPositiveIntScalar);
 
107
  p.addParamValue ("t", mindist, isNonNegativeIntScalar);
 
108
  p.addParamValue ("r", eps0, isPositiveScalar);
 
109
  p.addParamValue ("s", steps, isPositiveIntScalar);
 
110
  p.addSwitch ("verbose");
 
111
 
 
112
  p.parse (varargin{:});
 
113
 
 
114
  # Assign inputs
 
115
  embdim  = p.Results.m;
 
116
  delay   = p.Results.d;
 
117
  mindist = p.Results.t;
 
118
  eps0    = p.Results.r;
 
119
  epsset  = !ismember ('r', p.UsingDefaults);
 
120
  steps   = p.Results.s;
 
121
  verbose = p.Results.verbose;
 
122
 
 
123
  # Correct X to always have more rows than columns
 
124
  trnspsd = false;
 
125
  if (rows (X) < columns (X))
 
126
    X = X.';
 
127
    trnspsd = true;
 
128
  endif
 
129
 
 
130
  output = __lyap_r__ (X, embdim, delay, mindist, eps0, epsset, steps,...
 
131
                       verbose);
 
132
 
 
133
  if (trnspsd)
 
134
    output = output.';
 
135
  endif
 
136
 
 
137
endfunction
 
138
 
 
139
%!demo
 
140
%! idx = (1:2500).';
 
141
%! in = sin (idx ./ 360) + cos (idx ./ 180);
 
142
%! mmax = 15;
 
143
%!
 
144
%! cla reset
 
145
%! hold on
 
146
%! for i=2:mmax
 
147
%!   res = lyap_r (in, 'm', i, 'd', 6, 's',400,'t',200);
 
148
%!   plot (res(:,1),res(:,2),'r');
 
149
%! endfor
 
150
%! axis tight
 
151
%! xlabel ("t [flow samples]");
 
152
%! ylabel ("S(eps, embed, t)");
 
153
%! hold off
 
154
%!###############################################################
 
155
 
 
156
 
 
157
%!test
 
158
%! lyap_r_res =  [0 -2.983802;1 -2.980538;2 -2.962341;3 -2.931719;4 -2.891934;5 -2.846183;6 -2.797121;7 -2.74671;8 -2.69629;9 -2.646711;10 -2.598477];
 
159
%! in = sin((1:1000).'./360);
 
160
%! res = lyap_r (in, 'm',4 ,'d',6,'s',10,'t',100);
 
161
%! assert (res, lyap_r_res, -1e-6);
 
162
 
 
163
%% Check if transposed output works correctly
 
164
%!test
 
165
%! res1 = lyap_r(1:100);
 
166
%! res2 = lyap_r((1:100).');
 
167
%! assert(res1.',res2);
 
168
 
 
169
%!error <ranges> lyap_r (1)
 
170
%% Check if program does not run forever
 
171
%!error <too large> lyap_r (1:12)