home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Geek Gadgets 1
/
ADE-1.bin
/
ade-dist
/
octave-1.1.1p1-src.tgz
/
tar.out
/
fsf
/
octave
/
src
/
data.cc
< prev
next >
Wrap
C/C++ Source or Header
|
1996-09-28
|
19KB
|
1,040 lines
// data.cc -*- C++ -*-
/*
Copyright (C) 1992, 1993, 1994, 1995 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, 675 Mass Ave, Cambridge, MA 02139, 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 "tree-const.h"
#include "user-prefs.h"
#include "help.h"
#include "utils.h"
#include "error.h"
#include "gripes.h"
#include "defun.h"
#ifndef MIN
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#endif
#ifndef ABS
#define ABS(x) (((x) < 0) ? (-x) : (x))
#endif
DEFUN ("all", Fall, Sall, 1, 1,
"all (X): are all elements of X nonzero?")
{
Octave_object retval;
int nargin = args.length ();
if (nargin == 1 && args(0).is_defined ())
retval = args(0).all ();
else
print_usage ("all");
return retval;
}
DEFUN ("any", Fany, Sany, 1, 1,
"any (X): are any elements of X nonzero?")
{
Octave_object 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_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.elem (i, j) = f (x, y.elem (i, j));
return retval;
}
static Matrix
map (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.elem (i, j) = f (x.elem (i, j), y);
return retval;
}
static Matrix
map (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.elem (i, j) = f (x.elem (i, j), y.elem (i, j));
return retval;
}
DEFUN ("atan2", Fatan2, Satan2, 2, 1,
"atan2 (Y, X): atan (Y / X) in range -pi to pi")
{
Octave_object retval;
int nargin = args.length ();
if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
{
tree_constant arg_y = args(0);
tree_constant 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 (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 (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 (atan2, y, x);
}
}
else
error ("atan2: nonconformant matrices");
}
else
print_usage ("atan2");
return retval;
}
DEFUN ("cumprod", Fcumprod, Scumprod, 1, 1,
"cumprod (X): cumulative products")
{
Octave_object retval;
int nargin = args.length ();
if (nargin == 1)
{
tree_constant 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", Fcumsum, Scumsum, 1, 1,
"cumsum (X): cumulative sums")
{
Octave_object retval;
int nargin = args.length ();
if (nargin == 1)
{
tree_constant 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 tree_constant
make_diag (const Matrix& v, int k)
{
int nr = v.rows ();
int nc = v.columns ();
assert (nc == 1 || nr == 1);
tree_constant 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.elem (i+roff, i+coff) = v.elem (0, i);
retval = tree_constant (m);
}
else
{
int n = nr + ABS (k);
Matrix m (n, n, 0.0);
for (int i = 0; i < nr; i++)
m.elem (i+roff, i+coff) = v.elem (i, 0);
retval = tree_constant (m);
}
return retval;
}
static tree_constant
make_diag (const ComplexMatrix& v, int k)
{
int nr = v.rows ();
int nc = v.columns ();
assert (nc == 1 || nr == 1);
tree_constant 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.elem (i+roff, i+coff) = v.elem (0, i);
retval = tree_constant (m);
}
else
{
int n = nr + ABS (k);
ComplexMatrix m (n, n, 0.0);
for (int i = 0; i < nr; i++)
m.elem (i+roff, i+coff) = v.elem (i, 0);
retval = tree_constant (m);
}
return retval;
}
static tree_constant
make_diag (const tree_constant& arg)
{
tree_constant 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 tree_constant
make_diag (const tree_constant& a, const tree_constant& b)
{
tree_constant retval;
double tmp = b.double_value ();
if (error_state)
{
error ("diag: invalid second argument");
return retval;
}
int k = NINT (tmp);
int n = ABS (k) + 1;
if (a.is_real_type ())
{
if (a.is_scalar_type ())
{
double d = a.double_value ();
if (k == 0)
retval = d;
else if (k > 0)
{
Matrix m (n, n, 0.0);
m.elem (0, k) = d;
retval = m;
}
else if (k < 0)
{
Matrix m (n, n, 0.0);
m.elem (-k, 0) = d;
retval = m;
}
}
else if (a.is_matrix_type ())
{
Matrix m = a.matrix_value ();
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
gripe_wrong_type_arg ("diag", a);
}
else if (a.is_complex_type ())
{
if (a.is_scalar_type ())
{
Complex c = a.complex_value ();
if (k == 0)
retval = c;
else if (k > 0)
{
ComplexMatrix m (n, n, 0.0);
m.elem (0, k) = c;
retval = m;
}
else if (k < 0)
{
ComplexMatrix m (n, n, 0.0);
m.elem (-k, 0) = c;
retval = m;
}
}
else if (a.is_matrix_type ())
{
ComplexMatrix cm = a.complex_matrix_value ();
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);
}
else
gripe_wrong_type_arg ("diag", a);
return retval;
}
DEFUN ("diag", Fdiag, Sdiag, 2, 1,
"diag (X [,k]): form/extract diagonals")
{
Octave_object 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", Fprod, Sprod, 1, 1,
"prod (X): products")
{
Octave_object retval;
int nargin = args.length ();
if (nargin == 1)
{
tree_constant 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", Fsize, Ssize, 2, 1,
"[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_object 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.elem (0, 0) = nr;
m.elem (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", Fsum, Ssum, 1, 1,
"sum (X): sum of elements")
{
Octave_object retval;
int nargin = args.length ();
if (nargin == 1)
{
tree_constant 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", Fsumsq, Ssumsq, 1, 1,
"sumsq (X): sum of squares of elements")
{
Octave_object retval;
int nargin = args.length ();
if (nargin == 1)
{
tree_constant 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_struct", Fis_struct, Sis_struct, 1, 1,
"is_struct (x): return nonzero if x is a structure")
{
Octave_object retval;
int nargin = args.length ();
if (nargin == 1)
{
tree_constant arg = args(0);
if (arg.is_map ())
retval = 1.0;
else
retval = 0.0;
}
else
print_usage ("is_struct");
return retval;
}
static void
check_dimensions (int& nr, int& nc, const char *warnfor)
{
if (nr < 0 || nc < 0)
{
if (user_pref.treat_neg_dim_as_zero)
{
nr = (nr < 0) ? 0 : nr;
nc = (nc < 0) ? 0 : nc;
if (user_pref.treat_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 tree_constant& 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.elem (0));
nc = NINT (v.elem (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 tree_constant& a, const tree_constant& b,
const char *warn_for, int& nr, int& nc)
{
nr = NINT (a.double_value ());
nc = 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 tree_constant
fill_matrix (const tree_constant& a, double val, const char *warn_for)
{
int nr, nc;
get_dimensions (a, warn_for, nr, nc);
if (error_state)
return tree_constant ();
Matrix m (nr, nc, val);
return m;
}
static tree_constant
fill_matrix (const tree_constant& a, const tree_constant& 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 tree_constant ();
Matrix m (nr, nc, val);
return m;
}
DEFUN ("ones", Fones, Sones, 2, 1,
"ones (N), ones (N, M), ones (X): create a matrix of all ones")
{
Octave_object 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", Fzeros, Szeros, 2, 1,
"zeros (N), zeros (N, M), zeros (X): create a matrix of all zeros")
{
Octave_object 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 tree_constant
identity_matrix (const tree_constant& a)
{
int nr, nc;
get_dimensions (a, "eye", nr, nc); // May set error_state.
if (error_state)
return tree_constant ();
Matrix m (nr, nc, 0.0);
if (nr > 0 && nc > 0)
{
int n = MIN (nr, nc);
for (int i = 0; i < n; i++)
m.elem (i, i) = 1.0;
}
return m;
}
static tree_constant
identity_matrix (const tree_constant& a, const tree_constant& b)
{
int nr, nc;
get_dimensions (a, b, "eye", nr, nc); // May set error_state.
if (error_state)
return tree_constant ();
Matrix m (nr, nc, 0.0);
if (nr > 0 && nc > 0)
{
int n = MIN (nr, nc);
for (int i = 0; i < n; i++)
m.elem (i, i) = 1.0;
}
return m;
}
DEFUN ("eye", Feye, Seye, 2, 1,
"eye (N), eye (N, M), eye (X): create an identity matrix")
{
Octave_object 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", Flinspace, Slinspace, 2, 1,
"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_object retval;
int nargin = args.length ();
int npoints = 100;
if (nargin == 3)
{
double n = args(2).double_value ();
if (! error_state)
npoints = NINT (n);
}
else
print_usage ("linspace");
if (! error_state)
{
if (npoints > 1)
{
tree_constant arg_1 = args(0);
tree_constant 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) = tree_constant (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) = tree_constant (rv, 0);
}
}
}
else
error ("linspace: npoints must be greater than 2");
}
return retval;
}
/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; page-delimiter: "^/\\*" ***
;;; End: ***
*/