Logo Search packages:      
Sourcecode: octave2.0 version File versions

data.cc

/*

Copyright (C) 1996 John W. Eaton

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 2, 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, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

*/

/*

The function builtin_pwd adapted from a similar function from GNU
Bash, the Bourne Again SHell, copyright (C) 1987, 1989, 1991 Free
Software Foundation, Inc.

*/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <cfloat>
#include <cmath>

#include <string>

#include "lo-ieee.h"
#include "str-vec.h"

#include "defun.h"
#include "error.h"
#include "gripes.h"
#include "help.h"
#include "oct-map.h"
#include "ov.h"
#include "variables.h"
#include "oct-obj.h"
#include "utils.h"

#ifndef MIN
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#endif

#ifndef ABS
#define ABS(x) (((x) < 0) ? (-x) : (x))
#endif

// Should expressions like ones (-1, 5) result in an empty matrix or
// an error?  A positive value means yes.  A negative value means
// yes, but print a warning message.  Zero means it should be
// considered an error.
static int Vtreat_neg_dim_as_zero;

DEFUN (all, args, ,
  "all (X): are all elements of X nonzero?")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 1 && args(0).is_defined ())
    retval = args(0).all ();
  else
    print_usage ("all");

  return retval;
}

DEFUN (any, args, ,
  "any (X): are any elements of X nonzero?")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 1 && args(0).is_defined ())
    retval = args(0).any ();
  else
    print_usage ("any");

  return retval;
}

// These mapping functions may also be useful in other places, eh?

typedef double (*d_dd_fcn) (double, double);

static Matrix
map_d_m (d_dd_fcn f, double x, const Matrix& y)
{
  int nr = y.rows ();
  int nc = y.columns ();

  Matrix retval (nr, nc);

  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      retval (i, j) = f (x, y (i, j));

  return retval;
}

static Matrix
map_m_d (d_dd_fcn f, const Matrix& x, double y)
{
  int nr = x.rows ();
  int nc = x.columns ();

  Matrix retval (nr, nc);

  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      retval (i, j) = f (x (i, j), y);

  return retval;
}

static Matrix
map_m_m (d_dd_fcn f, const Matrix& x, const Matrix& y)
{
  int x_nr = x.rows ();
  int x_nc = x.columns ();

  int y_nr = y.rows ();
  int y_nc = y.columns ();

  assert (x_nr == y_nr && x_nc == y_nc);

  Matrix retval (x_nr, x_nc);

  for (int j = 0; j < x_nc; j++)
    for (int i = 0; i < x_nr; i++)
      retval (i, j) = f (x (i, j), y (i, j));

  return retval;
}

DEFUN (atan2, args, ,
  "atan2 (Y, X): atan (Y / X) in range -pi to pi")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
    {
      octave_value arg_y = args(0);
      octave_value arg_x = args(1);

      int y_nr = arg_y.rows ();
      int y_nc = arg_y.columns ();

      int x_nr = arg_x.rows ();
      int x_nc = arg_x.columns ();

      int arg_y_empty = empty_arg ("atan2", y_nr, y_nc);
      int arg_x_empty = empty_arg ("atan2", x_nr, x_nc);

      if (arg_y_empty > 0 && arg_x_empty > 0)
      return Matrix ();
      else if (arg_y_empty || arg_x_empty)
      return retval;

      int y_is_scalar = (y_nr == 1 && y_nc == 1);
      int x_is_scalar = (x_nr == 1 && x_nc == 1);

      if (y_is_scalar && x_is_scalar)
      {
        double y = arg_y.double_value ();

        if (! error_state)
          {
            double x = arg_x.double_value ();

            if (! error_state)
            retval = atan2 (y, x);
          }
      }
      else if (y_is_scalar)
      {
        double y = arg_y.double_value ();

        if (! error_state)
          {
            Matrix x = arg_x.matrix_value ();

            if (! error_state)
            retval = map_d_m (atan2, y, x);
          }
      }
      else if (x_is_scalar)
      {
        Matrix y = arg_y.matrix_value ();

        if (! error_state)
          {
            double x = arg_x.double_value ();

            if (! error_state)
            retval = map_m_d (atan2, y, x);
          }
      }
      else if (y_nr == x_nr && y_nc == x_nc)
      {
        Matrix y = arg_y.matrix_value ();

        if (! error_state)
          {
            Matrix x = arg_x.matrix_value ();

            if (! error_state)
            retval = map_m_m (atan2, y, x);
          }
      }
      else
      error ("atan2: nonconformant matrices");
    }
  else
    print_usage ("atan2");

  return retval;
}

DEFUN (cumprod, args, ,
  "cumprod (X): cumulative products")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 1)
    {
      octave_value arg = args(0);

      if (arg.is_real_type ())
      {
        Matrix tmp = arg.matrix_value ();

        if (! error_state)
          retval(0) = tmp.cumprod ();
      }
      else if (arg.is_complex_type ())
      {
        ComplexMatrix tmp = arg.complex_matrix_value ();

        if (! error_state)
          retval(0) = tmp.cumprod ();
      }
      else
      {
        gripe_wrong_type_arg ("cumprod", arg);
        return retval;
      }
    }
  else
    print_usage ("cumprod");

  return retval;
}

DEFUN (cumsum, args, ,
  "cumsum (X): cumulative sums")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 1)
    {
      octave_value arg = args(0);

      if (arg.is_real_type ())
      {
        Matrix tmp = arg.matrix_value ();

        if (! error_state)
          retval(0) = tmp.cumsum ();
      }
      else if (arg.is_complex_type ())
      {
        ComplexMatrix tmp = arg.complex_matrix_value ();

        if (! error_state)
          retval(0) = tmp.cumsum ();
      }
      else
      {
        gripe_wrong_type_arg ("cumsum", arg);
        return retval;
      }
    }
  else
    print_usage ("cumsum");

  return retval;
}

static octave_value
make_diag (const Matrix& v, int k)
{
  int nr = v.rows ();
  int nc = v.columns ();
  assert (nc == 1 || nr == 1);

  octave_value retval;

  int roff = 0;
  int coff = 0;
  if (k > 0)
    {
      roff = 0;
      coff = k;
    }
  else if (k < 0)
    {
      roff = -k;
      coff = 0;
    }

  if (nr == 1)
    {
      int n = nc + ABS (k);
      Matrix m (n, n, 0.0);
      for (int i = 0; i < nc; i++)
      m (i+roff, i+coff) = v (0, i);
      retval = octave_value (m);
    }
  else
    {
      int n = nr + ABS (k);
      Matrix m (n, n, 0.0);
      for (int i = 0; i < nr; i++)
      m (i+roff, i+coff) = v (i, 0);
      retval = octave_value (m);
    }

  return retval;
}

static octave_value
make_diag (const ComplexMatrix& v, int k)
{
  int nr = v.rows ();
  int nc = v.columns ();
  assert (nc == 1 || nr == 1);

  octave_value retval;

  int roff = 0;
  int coff = 0;
  if (k > 0)
    {
      roff = 0;
      coff = k;
    }
  else if (k < 0)
    {
      roff = -k;
      coff = 0;
    }

  if (nr == 1)
    {
      int n = nc + ABS (k);
      ComplexMatrix m (n, n, 0.0);
      for (int i = 0; i < nc; i++)
      m (i+roff, i+coff) = v (0, i);
      retval = octave_value (m);
    }
  else
    {
      int n = nr + ABS (k);
      ComplexMatrix m (n, n, 0.0);
      for (int i = 0; i < nr; i++)
      m (i+roff, i+coff) = v (i, 0);
      retval = octave_value (m);
    }

  return retval;
}

static octave_value
make_diag (const octave_value& arg)
{
  octave_value retval;

  if (arg.is_real_type ())
    {
      Matrix m = arg.matrix_value ();

      if (! error_state)
      {
        int nr = m.rows ();
        int nc = m.columns ();

        if (nr == 0 || nc == 0)
          retval = Matrix ();
        else if (nr == 1 || nc == 1)
          retval = make_diag (m, 0);
        else
          {
            ColumnVector v = m.diag ();
            if (v.capacity () > 0)
            retval = v;
          }
      }
      else
      gripe_wrong_type_arg ("diag", arg);
    }
  else if (arg.is_complex_type ())
    {
      ComplexMatrix cm = arg.complex_matrix_value ();

      if (! error_state)
      {
        int nr = cm.rows ();
        int nc = cm.columns ();

        if (nr == 0 || nc == 0)
          retval = Matrix ();
        else if (nr == 1 || nc == 1)
          retval = make_diag (cm, 0);
        else
          {
            ComplexColumnVector v = cm.diag ();
            if (v.capacity () > 0)
            retval = v;
          }
      }
      else
      gripe_wrong_type_arg ("diag", arg);
    }
  else
    gripe_wrong_type_arg ("diag", arg);

  return retval;
}

static octave_value
make_diag (const octave_value& a, const octave_value& b)
{
  octave_value retval;

  double tmp = b.double_value ();

  if (error_state)
    {
      error ("diag: invalid second argument");      
      return retval;
    }

  int k = NINT (tmp);

  if (a.is_real_type ())
    {
      Matrix m = a.matrix_value ();

      if (! error_state)
      {
        int nr = m.rows ();
        int nc = m.columns ();

        if (nr == 0 || nc == 0)
          retval = Matrix ();
        else if (nr == 1 || nc == 1)
          retval = make_diag (m, k);
        else
          {
            ColumnVector d = m.diag (k);
            retval = d;
          }
      }
    }
  else if (a.is_complex_type ())
    {
      ComplexMatrix cm = a.complex_matrix_value ();

      if (! error_state)
      {
        int nr = cm.rows ();
        int nc = cm.columns ();

        if (nr == 0 || nc == 0)
          retval = Matrix ();
        else if (nr == 1 || nc == 1)
          retval = make_diag (cm, k);
        else
          {
            ComplexColumnVector d = cm.diag (k);
            retval = d;
          }
      }
    }
  else
    gripe_wrong_type_arg ("diag", a);

  return retval;
}

DEFUN (diag, args, ,
  "diag (X [,k]): form/extract diagonals")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 1 && args(0).is_defined ())
    retval = make_diag (args(0));
  else if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
    retval = make_diag (args(0), args(1));
  else
    print_usage ("diag");

  return retval;
}

DEFUN (prod, args, ,
  "prod (X): products")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 1)
    {
      octave_value arg = args(0);

      if (arg.is_real_type ())
      {
        Matrix tmp = arg.matrix_value ();

        if (! error_state)
          retval(0) = tmp.prod ();
      }
      else if (arg.is_complex_type ())
      {
        ComplexMatrix tmp = arg.complex_matrix_value ();

        if (! error_state)
          retval(0) = tmp.prod ();
      }
      else
      {
        gripe_wrong_type_arg ("prod", arg);
        return retval;
      }
    }
  else
    print_usage ("prod");

  return retval;
}

DEFUN (size, args, nargout,
  "[m, n] = size (x): return rows and columns of X\n\
\n\
d = size (x): return number of rows and columns of x as a row vector\n\
\n\
m = size (x, 1): return number of rows in x\n\
m = size (x, 2): return number of columns in x")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 1 && nargout < 3)
    {
      int nr = args(0).rows ();
      int nc = args(0).columns ();

      if (nargout == 0 || nargout == 1)
      {
        Matrix m (1, 2);
        m (0, 0) = nr;
        m (0, 1) = nc;
        retval = m;
      }
      else if (nargout == 2)
      {
        retval(1) = (double) nc;
        retval(0) = (double) nr;
      }
    }
  else if (nargin == 2 && nargout < 2)
    {
      int nd = NINT (args(1).double_value ());

      if (error_state)
      error ("size: expecting scalar as second argument");
      else
      {
        if (nd == 1)
          retval(0) = (double) (args(0).rows ());
        else if (nd == 2)
          retval(0) = (double) (args(0).columns ());
        else
          error ("size: invalid second argument -- expecting 1 or 2");
      }
    }
  else
    print_usage ("size");

  return retval;
}

DEFUN (sum, args, ,
  "sum (X): sum of elements")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 1)
    {
      octave_value arg = args(0);

      if (arg.is_real_type ())
      {
        Matrix tmp = arg.matrix_value ();

        if (! error_state)
          retval(0) = tmp.sum ();
      }
      else if (arg.is_complex_type ())
      {
        ComplexMatrix tmp = arg.complex_matrix_value ();

        if (! error_state)
          retval(0) = tmp.sum ();
      }
      else
      {
        gripe_wrong_type_arg ("sum", arg);
        return retval;
      }
    }
  else
    print_usage ("sum");

  return retval;
}

DEFUN (sumsq, args, ,
  "sumsq (X): sum of squares of elements.\n\
\n\
This function is equivalent to computing\n\
\n\
  sum (X .* conj (X))\n\
\n\
but it uses less memory and avoids calling conj if X is real.")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 1)
    {
      octave_value arg = args(0);

      if (arg.is_real_type ())
      {
        Matrix tmp = arg.matrix_value ();

        if (! error_state)
          retval(0) = tmp.sumsq ();
      }
      else if (arg.is_complex_type ())
      {
        ComplexMatrix tmp = arg.complex_matrix_value ();

        if (! error_state)
          retval(0) = tmp.sumsq ();
      }
      else
      {
        gripe_wrong_type_arg ("sumsq", arg);
        return retval;
      }
    }
  else
    print_usage ("sumsq");

  return retval;
}

DEFUN (is_complex, args, ,
  "is_complex (x): return nonzero if x is a complex numeric object")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).is_complex_type () ? 1.0 : 0.0;
  else
    print_usage ("is_complex");

  return retval;
}

DEFUN (isnumeric, args, ,
  "isnumeric (x): return nonzero if x is a numeric object")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).is_numeric_type () ? 1.0 : 0.0;
  else
    print_usage ("isnumeric");

  return retval;
}

DEFUN (is_struct, args, ,
  "is_struct (x): return nonzero if x is a structure")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 1)
    {
      octave_value arg = args(0);

      if (arg.is_map ())
      retval = 1.0;
      else
      retval = 0.0;
    }
  else
    print_usage ("is_struct");

  return retval;
}

DEFUN (struct_elements, args, ,
  "struct_elements (S)\n\
\n\
Return a list of the names of the elements of the structure S.")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 1)
    {
      if (args (0).is_map ())
      {
        Octave_map m = args(0).map_value ();
        retval(0) = m.make_name_list ();
      }
      else
      gripe_wrong_type_arg ("struct_elements", args (0));
    }
  else
    print_usage ("struct_elements");

  return retval;
}

DEFUN (struct_contains, args, ,
  "struct_contains (S, NAME)\n\
\n\
Return nonzero if S is a structure with element NAME.\n\
S must be a structure and NAME must be a string.")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 2)
    {
      retval = 0.0;

      if (args(0).is_map () && args(1).is_string ())
      {
        string s = args(1).string_value ();
        octave_value tmp = args(0).struct_elt_val (s, true);
        retval = (double) tmp.is_defined ();
      }
      else
      print_usage ("struct_contains");
    }
  else
    print_usage ("struct_contains");

  return retval;
}

static void
check_dimensions (int& nr, int& nc, const char *warnfor)
{
  if (nr < 0 || nc < 0)
    {
      if (Vtreat_neg_dim_as_zero)
      {
        nr = (nr < 0) ? 0 : nr;
        nc = (nc < 0) ? 0 : nc;

        if (Vtreat_neg_dim_as_zero < 0)
          warning ("%s: converting negative dimension to zero",
                 warnfor);
      }
      else
      error ("%s: can't create a matrix with negative dimensions",
             warnfor);
    }
}

static void
get_dimensions (const octave_value& a, const char *warn_for,
            int& nr, int& nc)
{
  if (a.is_scalar_type ())
    {
      double tmp = a.double_value ();
      nr = nc = NINT (tmp);
    }
  else
    {
      nr = a.rows ();
      nc = a.columns ();

      if ((nr == 1 && nc == 2) || (nr == 2 && nc == 1))
      {
        ColumnVector v = a.vector_value ();

        if (error_state)
          return;

        nr = NINT (v (0));
        nc = NINT (v (1));
      }
      else
      warning ("%s (A): use %s (size (A)) instead", warn_for, warn_for);
    }

  check_dimensions (nr, nc, warn_for); // May set error_state.
}

static void
get_dimensions (const octave_value& a, const octave_value& b,
            const char *warn_for, int& nr, int& nc)
{
  nr = (a.is_empty ()) ? 0 : NINT (a.double_value ());
  nc = (b.is_empty ()) ? 0 : NINT (b.double_value ());

  if (error_state)
    error ("%s: expecting two scalar arguments", warn_for);
  else
    check_dimensions (nr, nc, warn_for); // May set error_state.
}

static octave_value
fill_matrix (const octave_value& a, double val, const char *warn_for)
{
  int nr, nc;
  get_dimensions (a, warn_for, nr, nc);

  if (error_state)
    return  octave_value ();

  Matrix m (nr, nc, val);

  return m;
}

static octave_value
fill_matrix (const octave_value& a, const octave_value& b,
           double val, const char *warn_for)
{
  int nr, nc;
  get_dimensions (a, b, warn_for, nr, nc); // May set error_state.

  if (error_state)
    return octave_value ();

  Matrix m (nr, nc, val);

  return m;
}

DEFUN (ones, args, ,
  "ones (N), ones (N, M), ones (X): create a matrix of all ones")
{
  octave_value_list retval;

  int nargin = args.length ();

  switch (nargin)
    {
    case 0:
      retval = 1.0;
      break;

    case 1:
      retval = fill_matrix (args(0), 1.0, "ones");
      break;

    case 2:
      retval = fill_matrix (args(0), args(1), 1.0, "ones");
      break;

    default:
      print_usage ("ones");
      break;
    }

  return retval;
}

DEFUN (zeros, args, ,
  "zeros (N), zeros (N, M), zeros (X): create a matrix of all zeros")
{
  octave_value_list retval;

  int nargin = args.length ();

  switch (nargin)
    {
    case 0:
      retval = 0.0;
      break;

    case 1:
      retval = fill_matrix (args(0), 0.0, "zeros");
      break;

    case 2:
      retval = fill_matrix (args(0), args(1), 0.0, "zeros");
      break;

    default:
      print_usage ("zeros");
      break;
    }

  return retval;
}

static octave_value
identity_matrix (const octave_value& a)
{
  int nr, nc;
  get_dimensions (a, "eye", nr, nc); // May set error_state.

  if (error_state)
    return octave_value ();

  Matrix m (nr, nc, 0.0);

  if (nr > 0 && nc > 0)
    {
      int n = MIN (nr, nc);
      for (int i = 0; i < n; i++)
      m (i, i) = 1.0;
    }

  return m;
}

static octave_value
identity_matrix (const octave_value& a, const octave_value& b)
{
  int nr, nc;
  get_dimensions (a, b, "eye", nr, nc);  // May set error_state.

  if (error_state)
    return octave_value ();

  Matrix m (nr, nc, 0.0);

  if (nr > 0 && nc > 0)
    {
      int n = MIN (nr, nc);
      for (int i = 0; i < n; i++)
      m (i, i) = 1.0;
    }

  return m;
}

DEFUN (eye, args, ,
  "eye (N), eye (N, M), eye (X): create an identity matrix")
{
  octave_value_list retval;

  int nargin = args.length ();

  switch (nargin)
    {
    case 0:
      retval = 1.0;
      break;

    case 1:
      retval = identity_matrix (args(0));
      break;

    case 2:
      retval = identity_matrix (args(0), args(1));
      break;

    default:
      print_usage ("eye");
      break;
    }

  return retval;
}

DEFUN (linspace, args, ,
  "usage: linspace (x1, x2, n)\n\
\n\
Return a vector of n equally spaced points between x1 and x2\n\
inclusive.\n\
\n\
If the final argument is omitted, n = 100 is assumed.\n\
\n\
All three arguments must be scalars.\n\
\n\
See also: logspace")
{
  octave_value_list retval;

  int nargin = args.length ();

  int npoints = 100;

  if (nargin != 2 && nargin != 3)
    {
      print_usage ("linspace");
      return retval;
    }

  if (nargin == 3)
    {
      double n = args(2).double_value ();

      if (! error_state)
      npoints = NINT (n);
    }

  if (! error_state)
    {
      if (npoints > 1)
      {
        octave_value arg_1 = args(0);
        octave_value arg_2 = args(1);

        if (arg_1.is_complex_type () || arg_2.is_complex_type ())
          {
            Complex x1 = arg_1.complex_value ();
            Complex x2 = arg_2.complex_value ();

            if (! error_state)
            {
              ComplexRowVector rv = linspace (x1, x2, npoints);

              if (! error_state)
                retval (0) = octave_value (rv, 0);
            }
          }
        else
          {
            double x1 = arg_1.double_value ();
            double x2 = arg_2.double_value ();

            if (! error_state)
            {
              RowVector rv = linspace (x1, x2, npoints);

              if (! error_state)
                retval (0) = octave_value (rv, 0);
            }
          }
      }
      else
      error ("linspace: npoints must be greater than 2");
    }

  return retval;
}

static int
treat_neg_dim_as_zero (void)
{
  Vtreat_neg_dim_as_zero = check_preference ("treat_neg_dim_as_zero");

  return 0;
}

void
symbols_of_data (void)
{
  DEFCONST (I, Complex (0.0, 1.0),
    "sqrt (-1)");

  DEFCONST (Inf, octave_Inf,
    "infinity");

  DEFCONST (J, Complex (0.0, 1.0),
    "sqrt (-1)");

  DEFCONST (NaN, octave_NaN,
    "not a number");

#if defined (M_E)
  double e_val = M_E;
#else
  double e_val = exp (1.0);
#endif

  DEFCONST (e, e_val,
    "exp (1)");

  DEFCONST (eps, DBL_EPSILON,
    "machine precision");

  DEFCONST (i, Complex (0.0, 1.0),
    "sqrt (-1)");

  DEFCONST (inf, octave_Inf,
    "infinity");

  DEFCONST (j, Complex (0.0, 1.0),
    "sqrt (-1)");

  DEFCONST (nan, octave_NaN,
    "not a number");

#if defined (M_PI)
  double pi_val = M_PI;
#else
  double pi_val = 4.0 * atan (1.0);
#endif

  DEFCONST (pi, pi_val,
    "ratio of the circumference of a circle to its diameter");

  DEFCONST (realmax, DBL_MAX,
    "realmax (): return largest representable floating point number");

  DEFCONST (realmin, DBL_MIN,
    "realmin (): return smallest representable floating point number");

  DEFVAR (treat_neg_dim_as_zero, 0.0, 0, treat_neg_dim_as_zero,
    "convert negative dimensions to zero");
}

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/

Generated by  Doxygen 1.6.0   Back to index