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

/* -*- coding: utf-8 -*- */

/* Copyright (C) 1996-2015 Piotr Held <pjheld@gmail.com>

 *

 * 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/>.

 */

/********************************************************************/

/********************************************************************/

#define HELPTEXT "\

-*- texinfo -*-\n\

@deftypefn{Function File} {@var{cleaned} =} lazy (@var{X}, @var{m}, @var{rv})\n\

@deftypefnx{Function File} {[@var{cleaned}, @var{diff}] =} lazy (@var{X}, @var{m}, @var{rv})\n\

@deftypefnx{Function File} {@dots{} =} lazy (@var{X}, @var{m}, @var{rv}, @var{imax})\n\

\n\

Performs simple nonlinear noise reduction\n\

\n\

@strong {Inputs}\n\

\n\

@table @var\n\

@item X\n\

Must be realvector. If it is a row vector then the output will be row vectors as well.\n\

@item m\n\

Embedding dimension. Must be postive integer.\n\

@item rv\n\

If @var{rv} > 0 then it is equal to the absolute radius of the neighbourhoods. If @var{rv} < 0 then its opposite (-@var{rv}) is equal to the fraction of standard deviation used. It cannot be equal 0.\n\

@item imax\n\

The number of iterations [default = 1].\n\

@end table\n\

\n\

@strong {Output}\n\

\n\

@table @var\n\

@item cleaned\n\

Vector containing the cleaned data.\n\

@item diff\n\

Difference between the clean and noisy data. \n\

@end table\n\

\n\

See the demo for example of how lazy works. \n\

\n\

@strong{Algorithm}@*\n\

Uses TISEAN package lazy\n\

@end deftypefn"

/******************************************************************************/

/******************************************************************************/

/* Author: Piotr Held <pjheld@gmail.com>.

 * This function is based on lazy of TISEAN 3.0.1 

 * https://github.com/heggus/Tisean"

 */

 

#include <octave/oct.h>

#include <octave/f77-fcn.h>

 

#define DEFAULT_IMAX 1

 

// In order to avoid clobbered warnings transposed is initialized globally.

bool transposed;

 

extern "C"

{

  F77_RET_T

  F77_FUNC (ts_lazy, TS_LAZY)

            (const int& m, const double& rv,

             const int& imax, const int& lines_read,

             double* in_out1, double* in_out2);

}

 

 

DEFUN_DLD (lazy, args, nargout, HELPTEXT)

{

  octave_value_list retval;

  int nargin = args.length ();

 

 

  if ((nargin != 4) && (nargin != 3))

    {

      print_usage ();

    }

  else if (nargout > 2)

    {

      error_with_id ("Octave:invalid-fun-call", \

                     "Lazy only produces two outputs");

    }

  else

    {

    // Assigning inputs

      Matrix in_out1 = args(0).matrix_value();

      int m          = args(1).int_value();

      double rv      = args(2).double_value();

      int imax       = DEFAULT_IMAX;

 

      if (nargin == 4)

        imax = args(3).int_value();

 

// --- DATA VALIDATION ---

 

    // Checking if matrix is vector.

 

      int rows = in_out1.rows();

      int cols = in_out1.cols();

 

      if (((rows != 1) && (cols != 1)) || (cols == 0) || (rows == 0))

        error_with_id ("Octave:invalid-input-arg",\

                       "input X must be a vector");

 

    // Checking parameters

      if (m < 1)

        error_with_id ("Octave:invalid-input-arg",\

                       "Embedding dimension (M) must be a positive integer");

 

      if (rv == 0.)

        error_with_id ("Octave:invalid-input-arg",\

                       "Set either radious of neighbourhoods"

                       " or fraction of std deviation");

 

      if (imax < 1)

        error_with_id ("Octave:invalid-input-arg",\

                       "Number of iterations (IMAX) must be a positive "

                       "integer");

 

      if (! error_state)

        {

        // If vector is in 1 row: transpose (we will transpose the output to fit)

            transposed = 0;

 

          if ((rows == 1) && (cols > 1))

            {

              transposed = 1;

              in_out1    = in_out1.transpose();

            }

 

          int lines_read = in_out1.numel();

          NDArray in_out2 (Matrix (lines_read, 1));

 

          F77_XFCN (ts_lazy, TS_LAZY,

                    (m, rv, imax, lines_read,

                      in_out1.fortran_vec(), in_out2.fortran_vec()));

 

          // Transpose the output to resemble the input

          if (transposed)

          {

            in_out1 = in_out1.transpose();

            in_out2 = in_out2.transpose();

          }

 

          retval(0) = in_out1;

          retval(1) = in_out2;

        }

    }

  return retval;

}

 

/*

%!demo

%! hen    = henon (10000);

%! "The following line is equvalent to 'addnoise -v0.02 hen' from TISEAN";

%! hen    = hen + std (hen) * 0.02 .* (-6 + sum (rand ([size(hen), 12]), 3));

%! hendel = delay (hen(:,1));

%! henlaz = lazy (hen(:,1),7,-0.06,3);

%! henlaz = delay (henlaz);

%!

%! subplot (2,3,1)

%! plot (hendel(:,1), hendel(:,2), 'b.','markersize', 3);

%! title ("Noisy data");

%! pbaspect ([1 1 1]);

%! axis tight

%! axis off

%!

%! subplot (2,3,4)

%! plot (henlaz(:,1), henlaz(:,2),'r.','markersize', 3);

%! title ("Clean data");

%! pbaspect ([1 1 1]);

%! axis tight

%! axis off

%!

%! subplot (2,3,[2 3 5 6])

%! plot (hendel(:,1), hendel(:,2), 'b.','markersize', 3,...

%!       henlaz(:,1), henlaz(:,2),'r.','markersize', 3);

%! legend ("Noisy", "Clean");

%! title ("Superimposed data");

%! axis tight

 

%!###############################################################

 

%!fail("lazy([(1:10);(1:10)],7,-0.06)");

 

%!fail("lazy((1:10),0,0.04)");

 

%!fail("[a,b,c] = lazy((1:10),1,0.05)");

 

%!test

%! "In is generated from Octave using 'in = 1 + 0.5 * rand(10,1);'";

%! in = [1.47007925526322;1.168775342017635;1.10943000146922; 1.174293926353764;  1.075741574572656; 1.373465364407417; 1.089417388489702; 1.403669883669071;1.452726826806777;  1.016960990335037];

%! "res was generated using 'lazy -m1 -v0.06 in.dat' from TISEAN 'lazy'";

%! res = [1.47007930, 0.00000000; 1.17153454, -2.75921822E-03; 1.10942996, 0.00000000; 1.17153454, 2.75933743E-03; 1.07574153, 0.00000000; 1.37346542, 0.00000000; 1.08941734, 0.00000000; 1.40366983, 0.00000000; 1.45272684, 0.00000000; 1.01696098, 0.00000000];

%! [al,bl] = lazy(in, 1, -0.06);

%! assert([al,bl],res,1e-6);

*/