[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[patch] Re-arrange matrix min/max reductions.
|
From: |
Edward Jason Riedy |
|
Subject: |
[patch] Re-arrange matrix min/max reductions. |
|
Date: |
Thu, 17 Aug 2000 13:45:53 -0700 |
The included patch re-arranges the min/max reductions, but does not
change the NaN semantics (yet). This means {row,column}_{min,max}
have different semantics than the rest of the min/max reductions
scattered through Octave.
It does, however, change the complex matrix min/max semantics. If
a matrix is complex, it is complex everywhere after this patch.
Octave previously would apply real min/max comparisons if a column
or row is entirely real. A real min over [-5, 4] returns -5, but
a complex one returns 4 (magnitude).
I can see this leading to confusion, as someone doesn't necessarily
know if an intermediate result just happens to have an all-real row
or column.
I'm not sure how to hook mx-reductions.h into the build system,
or even it it needs to be. It's meant to be entirely internal.
There are places with heavy macros used to define helper functions
and members. Equivalent template-based solutions are much longer
and more convoluted, or at least the ones I tried.
If this one's ok, I'll change the interpreter routines in minmax.cc
in my next burst of spare time.
Jason
liboctave/ChangeLog:
* mx-reductions.h: New file. Helps define matrix reductions.
* dMatrix.h (Matrix::row_min, Matrix::row_max, Matrix::column_min,
Matrix::column_max): add bool ignore_nan parameter (default false).
* dMatrix.cc: include mx-reductions.h, define reduction helpers
(Matrix::row_min): Replace with mx-reductions.h:DEF_COMP_OP call.
(Matrix::row_max): Ditto.
(Matrix::column_min): Ditto.
(Matrix::column_max): Ditto.
* cMatrix.hh: Changes analagous to those in dMatrix.h.
* cMatrix.cc: Changes analagous to those in dMatrix.cc.
* lo-mappers.h: Add xmin_nan, xmax_nan declarations.
* lo-mappers.cc (xmin_nan): Added double and Complex as copy of
old xmin.
(xmax_nan): Ditto.
(xmin): Replaced xisnan(x) with xisnan(y) to prefer returning
numbers rather than NaNs.
(xmax): Ditto.
* oct-cmplx.h: Include <functional>
(ComplexLess): std::binary_function implementing less-than over
complex numbers by magnitude.
(ComplexGreater): Analagous, for greater-than.
src/ChangeLog:
* minmax.cc (SM_COMP): Macro added to define scalar-matrix
comparisons.
(MS_COMP): Same for matrix-scalar.
(MM_COMP): Same for matrix-matrix.
(min): All definitions replaced by calls to the above macros.
(max): Ditto.
diff -X ../diff-ignore-patterns -urP ../CVS-srcs/octave/liboctave/CMatrix.cc
liboctave/CMatrix.cc
--- ../CVS-srcs/octave/liboctave/CMatrix.cc Mon Feb 7 20:35:41 2000
+++ liboctave/CMatrix.cc Thu Aug 17 11:45:45 2000
@@ -54,6 +54,7 @@
#include "mx-cm-s.h"
#include "mx-inlines.cc"
#include "oct-cmplx.h"
+#include "mx-reductions.h"
// Fortran functions we call.
@@ -2616,253 +2617,16 @@
return retval;
}
-ComplexColumnVector
-ComplexMatrix::row_min (void) const
-{
- Array<int> index;
- return row_min (index);
-}
-
-ComplexColumnVector
-ComplexMatrix::row_min (Array<int>& index) const
-{
- ComplexColumnVector result;
-
- int nr = rows ();
- int nc = cols ();
-
- if (nr > 0 && nc > 0)
- {
- result.resize (nr);
- index.resize (nr);
-
- for (int i = 0; i < nr; i++)
- {
- int idx_j = 0;
-
- Complex tmp_min = elem (i, idx_j);
-
- bool real_only = row_is_real_only (i);
-
- double abs_min = real_only ? real (tmp_min) : abs (tmp_min);
-
- if (xisnan (tmp_min))
- idx_j = -1;
- else
- {
- for (int j = 1; j < nc; j++)
- {
- Complex tmp = elem (i, j);
-
- double abs_tmp = real_only ? real (tmp) : abs (tmp);
-
- if (xisnan (tmp))
- {
- idx_j = -1;
- break;
- }
- else if (abs_tmp < abs_min)
- {
- idx_j = j;
- tmp_min = tmp;
- abs_min = abs_tmp;
- }
- }
-
- result.elem (i) = (idx_j < 0) ? Complex_NaN_result : tmp_min;
- index.elem (i) = idx_j;
- }
- }
- }
-
- return result;
-}
-
-ComplexColumnVector
-ComplexMatrix::row_max (void) const
-{
- Array<int> index;
- return row_max (index);
-}
-
-ComplexColumnVector
-ComplexMatrix::row_max (Array<int>& index) const
-{
- ComplexColumnVector result;
-
- int nr = rows ();
- int nc = cols ();
-
- if (nr > 0 && nc > 0)
- {
- result.resize (nr);
- index.resize (nr);
-
- for (int i = 0; i < nr; i++)
- {
- int idx_j = 0;
-
- Complex tmp_max = elem (i, idx_j);
-
- bool real_only = row_is_real_only (i);
-
- double abs_max = real_only ? real (tmp_max) : abs (tmp_max);
-
- if (xisnan (tmp_max))
- idx_j = -1;
- else
- {
- for (int j = 1; j < nc; j++)
- {
- Complex tmp = elem (i, j);
-
- double abs_tmp = real_only ? real (tmp) : abs (tmp);
-
- if (xisnan (tmp))
- {
- idx_j = -1;
- break;
- }
- else if (abs_tmp > abs_max)
- {
- idx_j = j;
- tmp_max = tmp;
- abs_max = abs_tmp;
- }
- }
-
- result.elem (i) = (idx_j < 0) ? Complex_NaN_result : tmp_max;
- index.elem (i) = idx_j;
- }
- }
- }
-
- return result;
-}
-
-ComplexRowVector
-ComplexMatrix::column_min (void) const
-{
- Array<int> index;
- return column_min (index);
-}
-
-ComplexRowVector
-ComplexMatrix::column_min (Array<int>& index) const
-{
- ComplexRowVector result;
-
- int nr = rows ();
- int nc = cols ();
-
- if (nr > 0 && nc > 0)
- {
- result.resize (nc);
- index.resize (nc);
-
- for (int j = 0; j < nc; j++)
- {
- int idx_i = 0;
-
- Complex tmp_min = elem (idx_i, j);
-
- bool real_only = column_is_real_only (j);
-
- double abs_min = real_only ? real (tmp_min) : abs (tmp_min);
-
- if (xisnan (tmp_min))
- idx_i = -1;
- else
- {
- for (int i = 1; i < nr; i++)
- {
- Complex tmp = elem (i, j);
-
- double abs_tmp = real_only ? real (tmp) : abs (tmp);
-
- if (xisnan (tmp))
- {
- idx_i = -1;
- break;
- }
- else if (abs_tmp < abs_min)
- {
- idx_i = i;
- tmp_min = tmp;
- abs_min = abs_tmp;
- }
- }
-
- result.elem (j) = (idx_i < 0) ? Complex_NaN_result : tmp_min;
- index.elem (j) = idx_i;
- }
- }
- }
-
- return result;
-}
-
-ComplexRowVector
-ComplexMatrix::column_max (void) const
-{
- Array<int> index;
- return column_max (index);
-}
-
-ComplexRowVector
-ComplexMatrix::column_max (Array<int>& index) const
-{
- ComplexRowVector result;
-
- int nr = rows ();
- int nc = cols ();
-
- if (nr > 0 && nc > 0)
- {
- result.resize (nc);
- index.resize (nc);
-
- for (int j = 0; j < nc; j++)
- {
- int idx_i = 0;
-
- Complex tmp_max = elem (idx_i, j);
-
- bool real_only = column_is_real_only (j);
-
- double abs_max = real_only ? real (tmp_max) : abs (tmp_max);
-
- if (xisnan (tmp_max))
- idx_i = -1;
- else
- {
- for (int i = 1; i < nr; i++)
- {
- Complex tmp = elem (i, j);
-
- double abs_tmp = real_only ? real (tmp) : abs (tmp);
-
- if (xisnan (tmp))
- {
- idx_i = -1;
- break;
- }
- else if (abs_tmp > abs_max)
- {
- idx_i = i;
- tmp_max = tmp;
- abs_max = abs_tmp;
- }
- }
-
- result.elem (j) = (idx_i < 0) ? Complex_NaN_result : tmp_max;
- index.elem (j) = idx_i;
- }
- }
- }
-
- return result;
-}
+// Use definition macros from mx-reductions.h
+DEF_REDUCTION_HELPER(row, ComplexMatrix, ComplexColumn, nr)
+DEF_REDUCTION_HELPER(column, ComplexMatrix, ComplexRow, nc)
+
+// Use definition macros from mx-reductions.h
+// Defines row_min, row_max, column_min, column_max
+DEF_COMP_OP(row, ComplexMatrix, ComplexColumn, Complex, min, ComplexLess)
+DEF_COMP_OP(row, ComplexMatrix, ComplexColumn, Complex, max, ComplexGreater)
+DEF_COMP_OP(column, ComplexMatrix, ComplexRow, Complex, min, ComplexLess)
+DEF_COMP_OP(column, ComplexMatrix, ComplexRow, Complex, max, ComplexGreater)
// i/o
diff -X ../diff-ignore-patterns -urP ../CVS-srcs/octave/liboctave/CMatrix.h
liboctave/CMatrix.h
--- ../CVS-srcs/octave/liboctave/CMatrix.h Mon Feb 7 20:35:41 2000
+++ liboctave/CMatrix.h Thu Aug 17 11:38:52 2000
@@ -256,17 +256,21 @@
bool row_is_real_only (int) const;
bool column_is_real_only (int) const;
- ComplexColumnVector row_min (void) const;
- ComplexColumnVector row_max (void) const;
+ ComplexColumnVector row_min (bool ignore_nan = false) const;
+ ComplexColumnVector row_max (bool ignore_nan = false) const;
- ComplexColumnVector row_min (Array<int>& index) const;
- ComplexColumnVector row_max (Array<int>& index) const;
+ ComplexColumnVector row_min (Array<int>& index, bool ignore_nan = false)
+ const;
+ ComplexColumnVector row_max (Array<int>& index, bool ignore_nan = false)
+ const;
- ComplexRowVector column_min (void) const;
- ComplexRowVector column_max (void) const;
+ ComplexRowVector column_min (bool ignore_nan = false) const;
+ ComplexRowVector column_max (bool ignore_nan = false) const;
- ComplexRowVector column_min (Array<int>& index) const;
- ComplexRowVector column_max (Array<int>& index) const;
+ ComplexRowVector column_min (Array<int>& index, bool ignore_nan = false)
+ const;
+ ComplexRowVector column_max (Array<int>& index, bool ignore_nan = false)
+ const;
// i/o
diff -X ../diff-ignore-patterns -urP ../CVS-srcs/octave/liboctave/dMatrix.cc
liboctave/dMatrix.cc
--- ../CVS-srcs/octave/liboctave/dMatrix.cc Mon Feb 7 20:35:44 2000
+++ liboctave/dMatrix.cc Thu Aug 17 11:45:30 2000
@@ -31,6 +31,7 @@
#include <cfloat>
+#include <functional>
#include <iostream>
#include "byte-swap.h"
@@ -49,6 +50,8 @@
#include "mx-dm-m.h"
#include "mx-inlines.cc"
#include "oct-cmplx.h"
+#include "mx-reductions.h"
+
// Fortran functions we call.
@@ -2173,225 +2176,16 @@
return d;
}
-ColumnVector
-Matrix::row_min (void) const
-{
- Array<int> index;
- return row_min (index);
-}
-
-ColumnVector
-Matrix::row_min (Array<int>& index) const
-{
- ColumnVector result;
-
- int nr = rows ();
- int nc = cols ();
-
- if (nr > 0 && nc > 0)
- {
- result.resize (nr);
- index.resize (nr);
-
- for (int i = 0; i < nr; i++)
- {
- int idx_j = 0;
-
- double tmp_min = elem (i, idx_j);
-
- if (xisnan (tmp_min))
- idx_j = -1;
- else
- {
- for (int j = 1; j < nc; j++)
- {
- double tmp = elem (i, j);
-
- if (xisnan (tmp))
- {
- idx_j = -1;
- break;
- }
- else if (tmp < tmp_min)
- {
- idx_j = j;
- tmp_min = tmp;
- }
- }
- }
-
- result.elem (i) = (idx_j < 0) ? octave_NaN : tmp_min;
- index.elem (i) = idx_j;
- }
- }
-
- return result;
-}
-
-ColumnVector
-Matrix::row_max (void) const
-{
- Array<int> index;
- return row_max (index);
-}
-
-ColumnVector
-Matrix::row_max (Array<int>& index) const
-{
- ColumnVector result;
-
- int nr = rows ();
- int nc = cols ();
-
- if (nr > 0 && nc > 0)
- {
- result.resize (nr);
- index.resize (nr);
-
- for (int i = 0; i < nr; i++)
- {
- int idx_j = 0;
-
- double tmp_max = elem (i, idx_j);
-
- if (xisnan (tmp_max))
- idx_j = -1;
- else
- {
- for (int j = 1; j < nc; j++)
- {
- double tmp = elem (i, j);
-
- if (xisnan (tmp))
- {
- idx_j = -1;
- break;
- }
- else if (tmp > tmp_max)
- {
- idx_j = j;
- tmp_max = tmp;
- }
- }
- }
-
- result.elem (i) = (idx_j < 0) ? octave_NaN : tmp_max;
- index.elem (i) = idx_j;
- }
- }
-
- return result;
-}
-
-RowVector
-Matrix::column_min (void) const
-{
- Array<int> index;
- return column_min (index);
-}
-
-RowVector
-Matrix::column_min (Array<int>& index) const
-{
- RowVector result;
-
- int nr = rows ();
- int nc = cols ();
-
- if (nr > 0 && nc > 0)
- {
- result.resize (nc);
- index.resize (nc);
-
- for (int j = 0; j < nc; j++)
- {
- int idx_i = 0;
-
- double tmp_min = elem (idx_i, j);
-
- if (xisnan (tmp_min))
- idx_i = -1;
- else
- {
- for (int i = 1; i < nr; i++)
- {
- double tmp = elem (i, j);
-
- if (xisnan (tmp))
- {
- idx_i = -1;
- break;
- }
- else if (tmp < tmp_min)
- {
- idx_i = i;
- tmp_min = tmp;
- }
- }
- }
-
- result.elem (j) = (idx_i < 0) ? octave_NaN : tmp_min;
- index.elem (j) = idx_i;
- }
- }
-
- return result;
-}
-
-RowVector
-Matrix::column_max (void) const
-{
- Array<int> index;
- return column_max (index);
-}
-
-RowVector
-Matrix::column_max (Array<int>& index) const
-{
- RowVector result;
-
- int nr = rows ();
- int nc = cols ();
-
- if (nr > 0 && nc > 0)
- {
- result.resize (nc);
- index.resize (nc);
-
- for (int j = 0; j < nc; j++)
- {
- int idx_i = 0;
-
- double tmp_max = elem (idx_i, j);
-
- if (xisnan (tmp_max))
- idx_i = -1;
- else
- {
- for (int i = 1; i < nr; i++)
- {
- double tmp = elem (i, j);
-
- if (xisnan (tmp))
- {
- idx_i = -1;
- break;
- }
- else if (tmp > tmp_max)
- {
- idx_i = i;
- tmp_max = tmp;
- }
- }
- }
-
- result.elem (j) = (idx_i < 0) ? octave_NaN : tmp_max;
- index.elem (j) = idx_i;
- }
- }
-
- return result;
-}
+// Use definition macros from mx-reductions.h
+DEF_REDUCTION_HELPER(row, Matrix, Column, nr)
+DEF_REDUCTION_HELPER(column, Matrix, Row, nc)
+
+// Use definition macros from mx-reductions.h
+// Defines row_min, row_max, column_min, column_max
+DEF_COMP_OP(row, Matrix, Column, double, min, std::less<double>)
+DEF_COMP_OP(row, Matrix, Column, double, max, std::greater<double>)
+DEF_COMP_OP(column, Matrix, Row, double, min, std::less<double>)
+DEF_COMP_OP(column, Matrix, Row, double, max, std::greater<double>)
std::ostream&
operator << (std::ostream& os, const Matrix& a)
diff -X ../diff-ignore-patterns -urP ../CVS-srcs/octave/liboctave/dMatrix.h
liboctave/dMatrix.h
--- ../CVS-srcs/octave/liboctave/dMatrix.h Fri Jun 30 02:30:45 2000
+++ liboctave/dMatrix.h Thu Aug 17 10:17:04 2000
@@ -210,17 +210,17 @@
ColumnVector diag (void) const;
ColumnVector diag (int k) const;
- ColumnVector row_min (void) const;
- ColumnVector row_max (void) const;
+ ColumnVector row_min (bool ignore_nan = false) const;
+ ColumnVector row_max (bool ignore_nan = false) const;
- ColumnVector row_min (Array<int>& index) const;
- ColumnVector row_max (Array<int>& index) const;
+ ColumnVector row_min (Array<int>& index, bool ignore_nan = false) const;
+ ColumnVector row_max (Array<int>& index, bool ignore_nan = false) const;
- RowVector column_min (void) const;
- RowVector column_max (void) const;
+ RowVector column_min (bool ignore_nan = false) const;
+ RowVector column_max (bool ignore_nan = false) const;
- RowVector column_min (Array<int>& index) const;
- RowVector column_max (Array<int>& index) const;
+ RowVector column_min (Array<int>& index, bool ignore_nan = false) const;
+ RowVector column_max (Array<int>& index, bool ignore_nan = false) const;
// i/o
Binary files ../CVS-srcs/octave/liboctave/liboctave.so and
liboctave/liboctave.so differ
Binary files ../CVS-srcs/octave/liboctave/liboctave.so.2.1.31 and
liboctave/liboctave.so.2.1.31 differ
diff -X ../diff-ignore-patterns -urP ../CVS-srcs/octave/liboctave/lo-mappers.cc
liboctave/lo-mappers.cc
--- ../CVS-srcs/octave/liboctave/lo-mappers.cc Thu Sep 9 22:17:01 1999
+++ liboctave/lo-mappers.cc Wed Aug 16 16:15:49 2000
@@ -185,12 +185,24 @@
double
xmin (double x, double y)
{
- return x < y ? x : (xisnan (x) ? x : y);
+ return x < y ? x : (xisnan (y) ? x : y);
}
double
xmax (double x, double y)
{
+ return x > y ? x : (xisnan (y) ? x : y);
+}
+
+double
+xmin_nan (double x, double y)
+{
+ return x < y ? x : (xisnan (x) ? x : y);
+}
+
+double
+xmax_nan (double x, double y)
+{
return x > y ? x : (xisnan (x) ? x : y);
}
@@ -315,11 +327,23 @@
Complex
xmin (const Complex& x, const Complex& y)
{
- return abs (x) < abs (y) ? x : (xisnan (x) ? x : y);
+ return abs (x) < abs (y) ? x : (xisnan (y) ? x : y);
}
Complex
xmax (const Complex& x, const Complex& y)
+{
+ return abs (x) > abs (y) ? x : (xisnan (y) ? x : y);
+}
+
+Complex
+xmin_nan (const Complex& x, const Complex& y)
+{
+ return abs (x) < abs (y) ? x : (xisnan (x) ? x : y);
+}
+
+Complex
+xmax_nan (const Complex& x, const Complex& y)
{
return abs (x) > abs (y) ? x : (xisnan (x) ? x : y);
}
diff -X ../diff-ignore-patterns -urP ../CVS-srcs/octave/liboctave/lo-mappers.h
liboctave/lo-mappers.h
--- ../CVS-srcs/octave/liboctave/lo-mappers.h Mon Jul 12 20:35:32 1999
+++ liboctave/lo-mappers.h Wed Aug 16 16:16:14 2000
@@ -41,6 +41,8 @@
extern double xmin (double x, double y);
extern double xmax (double x, double y);
+extern double xmin_nan (double x, double y);
+extern double xmax_nan (double x, double y);
extern Complex acos (const Complex& x);
extern Complex acosh (const Complex& x);
@@ -63,6 +65,8 @@
extern Complex xmin (const Complex& x, const Complex& y);
extern Complex xmax (const Complex& x, const Complex& y);
+extern Complex xmin_nan (const Complex& x, const Complex& y);
+extern Complex xmax_nan (const Complex& x, const Complex& y);
#endif
diff -X ../diff-ignore-patterns -urP
../CVS-srcs/octave/liboctave/mx-reductions.h liboctave/mx-reductions.h
--- ../CVS-srcs/octave/liboctave/mx-reductions.h Wed Dec 31 16:00:00 1969
+++ liboctave/mx-reductions.h Thu Aug 17 13:15:17 2000
@@ -0,0 +1,246 @@
+/*
+
+Copyright (C) 2000 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 template functions and macros are to help define reduction member
+ functions like Matrix::row_min.
+*/
+
+#if !defined (octave_mx_reductions_h)
+#define octave_mx_reductions_h 1
+
+#include "lo-mappers.h"
+#include "Array.h"
+#include "MArray2.h"
+
+/*
+ The do_over_{rows,columns} template functions wrap the iteration
+ order. The ActionObj should define a 1-d reduction operation with
+ three function members:
+ * begin (value)
+ begin a row or column reduction (implied col/row inner_index 0)
+ * apply (inner_index, value)
+ accumulate value at col/row inner_index into the reduction
+ * end (outer_index)
+ finalize the reduction over row/col outer_index
+*/
+
+template <typename T, typename ActionObj>
+inline void
+do_over_rows (const MArray2<T>& A, ActionObj& obj)
+{
+ const int nc = A.cols();
+ const int nr = A.rows();
+
+ for (int i = 0; i < nr; ++i)
+ {
+ obj.begin(A.elem(i, 0));
+ for (int j = 1; j < nc; ++j)
+ {
+ obj.apply(j, A.elem(i, j));
+ }
+ obj.end(i);
+ }
+}
+
+template <typename T, typename ActionObj>
+inline void
+do_over_columns (const MArray2<T>& A, ActionObj& obj)
+{
+ const int nc = A.cols();
+ const int nr = A.rows();
+
+ for (int j = 0; j < nc; ++j)
+ {
+ obj.begin(A.elem(0, j));
+ for (int i = 1; i < nr; ++i)
+ {
+ obj.apply(i, A.elem(i, j));
+ }
+ obj.end(j);
+ }
+}
+
+/*
+ PredReducer satisfies the requirements for ActionObj above. It
+ reduces by applying the predicate from template param. Pred.
+
+ PredReducer_NaN is similar, but NaNs are treated specially. A
+ NaN value will induce an index of -1 and a value of NaN rather
+ than simply falling through the predicate.
+
+ In large matrices with fairly many Infs (or NaNs), it may be
+ faster to cause the do_over_* reduction loops to terminate
+ early. The ActionObj could either return a coded value or
+ throw an exception for that effect. Both choices render the
+ code even uglier, and so are avoided here.
+*/
+
+template <typename T, typename Pred>
+class PredReducer
+{
+public:
+ inline PredReducer (Array<T>& stored_val, Array<int>& stored_index,
+ Pred pred = Pred())
+ : stored_val(stored_val), stored_index(stored_index), pred(pred)
+ { }
+
+ inline void begin (T val)
+ {
+ tmp_val = val;
+ tmp_idx = 0;
+ }
+
+ inline void apply (int idx, T val)
+ {
+ if (pred(val, tmp_val))
+ {
+ tmp_val = val;
+ tmp_idx = idx;
+ }
+ }
+
+ inline void end (int outer_idx)
+ {
+ stored_val.elem(outer_idx) = tmp_val;
+ stored_index.elem(outer_idx) = tmp_idx;
+ }
+
+private:
+
+ Array<T>& stored_val;
+ Array<int>& stored_index;
+ Pred pred;
+
+ T tmp_val;
+ int tmp_idx;
+};
+
+template <typename T, typename Pred>
+class PredReducer_NaN
+{
+public:
+ inline PredReducer_NaN (Array<T>& stored_val, Array<int>& stored_index,
+ Pred pred = Pred())
+ : stored_val(stored_val), stored_index(stored_index), pred(pred)
+ { }
+
+ inline void begin (T val)
+ {
+ tmp_val = val;
+ if (xisnan (val))
+ {
+ tmp_idx = -1;
+ }
+ else
+ {
+ tmp_idx = 0;
+ }
+ }
+
+ inline void apply (int idx, T val)
+ {
+ if (xisnan (val))
+ {
+ tmp_val = val;
+ tmp_idx = -1;
+ return;
+ }
+
+ // NaN should fail pred...
+ if (pred (val, tmp_val))
+ {
+ tmp_val = val;
+ tmp_idx = idx;
+ }
+ }
+
+ inline void end (int outer_idx)
+ {
+ stored_val.elem(outer_idx) = tmp_val;
+ stored_index.elem(outer_idx) = tmp_idx;
+ }
+
+private:
+
+ Array<T>& stored_val;
+ Array<int>& stored_index;
+ Pred pred;
+
+ T tmp_val;
+ int tmp_idx;
+};
+
+/*
+ DEF_REDUCTION_HELPER generates a simple, non-member helper function
+ that applies a predicate reducer to a matrix.
+
+ DEF_COMP_OP uses the DEF_REDUCTION_HELPER functions to define a
+ matrix type's min- and max-reduction members.
+
+ The non-uniformity between Matrix and ComplexMatrix seems to be
+ captured more easily through macros than templates.
+*/
+
+#define DEF_REDUCTION_HELPER(Side, MatType, VectorType, Size) \
+/* GENERATED BY DEF_REDUCTION_HELPER */ \
+template <typename PredReducer> \
+static inline VectorType ## Vector \
+Side ## _reduction (const MatType& A, Array<int>& index) \
+{ \
+ VectorType ## Vector result; \
+ const int nr = A.rows (); \
+ const int nc = A.cols (); \
+ \
+ if (nr > 0 && nc > 0) \
+ { \
+ result.resize (Size); \
+ index.resize (Size); \
+ \
+ PredReducer pred (result, index); \
+ do_over_ ## Side ## s (A, pred); \
+ } \
+ \
+ return result; \
+}
+
+#define DEF_COMP_OP(Side, MatType, VectorType, Type, Op, Comp) \
+/* GENERATED BY DEF_COMP_OP */ \
+VectorType ## Vector \
+MatType:: ## Side ## _ ## Op (bool ignore_nan) const \
+{ \
+ Array<int> index; \
+ return Side ## _ ## Op (index, ignore_nan); \
+} \
+ \
+VectorType ## Vector \
+MatType:: ## Side ## _ ## Op (Array<int>& index, bool ignore_nan) const \
+{ \
+ if (ignore_nan) \
+ return Side ## _reduction< PredReducer<Type, Comp > > \
+ (*this, index); \
+ else \
+ return Side ## _reduction< PredReducer_NaN<Type, Comp > > \
+ (*this, index); \
+}
+
+#endif
diff -X ../diff-ignore-patterns -urP ../CVS-srcs/octave/liboctave/oct-cmplx.h
liboctave/oct-cmplx.h
--- ../CVS-srcs/octave/liboctave/oct-cmplx.h Tue Feb 1 02:07:22 2000
+++ liboctave/oct-cmplx.h Wed Aug 16 15:32:02 2000
@@ -24,8 +24,29 @@
#define octave_oct_cmplx_h 1
#include <complex>
+#include <functional>
typedef std::complex<double> Complex;
+
+struct ComplexLess
+ : public std::binary_function<Complex, Complex, bool>
+{
+ bool operator() (const Complex& x, const Complex& y) const
+ {
+ return abs(x) < abs(y);
+ }
+};
+
+struct ComplexGreater
+ : public std::binary_function<Complex, Complex, bool>
+{
+ bool operator() (const Complex& x, const Complex& y) const
+ {
+ return abs(x) > abs(y);
+ }
+};
+
+
#endif
diff -X ../diff-ignore-patterns -ur
../CVS-srcs/octave/src/DLD-FUNCTIONS/minmax.cc src/DLD-FUNCTIONS/minmax.cc
--- ../CVS-srcs/octave/src/DLD-FUNCTIONS/minmax.cc Tue Apr 11 12:02:05 2000
+++ src/DLD-FUNCTIONS/minmax.cc Thu Aug 17 12:55:35 2000
@@ -34,245 +34,100 @@
#include "gripes.h"
#include "oct-obj.h"
-// XXX FIXME XXX -- it would be nice to share code among the min/max
-// functions below.
+// Scalar-matrix comparisons
-static Matrix
-min (double d, const Matrix& m)
-{
- int nr = m.rows ();
- int nc = m.columns ();
-
- Matrix result (nr, nc);
-
- for (int j = 0; j < nc; j++)
- for (int i = 0; i < nr; i++)
- result (i, j) = xmin (d, m (i, j));
+#define SM_COMP(ScalarType, MatrixType, Op) \
+static MatrixType \
+Op (const ScalarType d, const MatrixType& m, const bool ignore_nan = false) \
+{ \
+ const int nr = m.rows (); \
+ const int nc = m.columns (); \
+ \
+ MatrixType result (nr, nc); \
+ \
+ for (int j = 0; j < nc; j++) \
+ for (int i = 0; i < nr; i++) \
+ if (ignore_nan) \
+ result (i, j) = x ## Op (d, m (i, j)); \
+ else \
+ result (i, j) = x ## Op ## _nan (d, m (i, j)); \
+ \
+ return result; \
+}
+
+SM_COMP(double, Matrix, min)
+SM_COMP(double, Matrix, max)
+
+SM_COMP(Complex, ComplexMatrix, min)
+SM_COMP(Complex, ComplexMatrix, max)
+
+#undef SM_COMP
+
+// Matrix-scalar comparisons
+
+#define MS_COMP(ScalarType, MatrixType, Op) \
+static MatrixType \
+Op (const MatrixType& m, const ScalarType d, const bool ignore_nan = false) \
+{ \
+ const int nr = m.rows (); \
+ const int nc = m.columns (); \
+ \
+ MatrixType result (nr, nc); \
+ \
+ for (int j = 0; j < nc; j++) \
+ for (int i = 0; i < nr; i++) \
+ if (ignore_nan) \
+ result (i, j) = x ## Op (m (i, j), d); \
+ else \
+ result (i, j) = x ## Op ## _nan (m (i, j), d); \
+ \
+ return result; \
+}
+
+MS_COMP(double, Matrix, min)
+MS_COMP(double, Matrix, max)
+
+MS_COMP(Complex, ComplexMatrix, min)
+MS_COMP(Complex, ComplexMatrix, max)
+
+#undef MS_COMP
+
+// Matrix-matrix comparisons
+
+#define MM_COMP(Op, MatType) \
+static MatType \
+Op (const MatType& a, const MatType& b, const bool ignore_nan = false) \
+{ \
+ const int nr = a.rows(); \
+ int nc = a.columns (); \
+ if (nr != b.rows () || nc != b.columns ()) \
+ { \
+ error ("two-arg min expecting args of same size"); \
+ return MatType (); \
+ } \
+ \
+ MatType result (nr, nc); \
+ \
+ for (int j = 0; j < nc; j++) \
+ for (int i = 0; i < nr; i++) \
+ if (ignore_nan) \
+ result (i, j) = x ## Op (a (i, j), b (i, j)); \
+ else \
+ result (i, j) = x ## Op ## _nan (a (i, j), b (i, j)); \
+ \
+ return result; \
+}
+
+MM_COMP(min, Matrix)
+MM_COMP(max, Matrix)
- return result;
-}
-
-static Matrix
-min (const Matrix& m, double d)
-{
- int nr = m.rows ();
- int nc = m.columns ();
-
- Matrix result (nr, nc);
-
- for (int j = 0; j < nc; j++)
- for (int i = 0; i < nr; i++)
- result (i, j) = xmin (m (i, j), d);
-
- return result;
-}
-
-static ComplexMatrix
-min (const Complex& c, const ComplexMatrix& m)
-{
- int nr = m.rows ();
- int nc = m.columns ();
-
- ComplexMatrix result (nr, nc);
-
- for (int j = 0; j < nc; j++)
- for (int i = 0; i < nr; i++)
- result (i, j) = xmin (c, m (i, j));
-
- return result;
-}
-
-static ComplexMatrix
-min (const ComplexMatrix& m, const Complex& c)
-{
- int nr = m.rows ();
- int nc = m.columns ();
-
- ComplexMatrix result (nr, nc);
-
- for (int j = 0; j < nc; j++)
- for (int i = 0; i < nr; i++)
- result (i, j) = xmin (m (i, j), c);
-
- return result;
-}
-
-static Matrix
-min (const Matrix& a, const Matrix& b)
-{
- int nr = a.rows ();
- int nc = a.columns ();
- if (nr != b.rows () || nc != b.columns ())
- {
- error ("two-arg min expecting args of same size");
- return Matrix ();
- }
-
- Matrix result (nr, nc);
-
- for (int j = 0; j < nc; j++)
- for (int i = 0; i < nr; i++)
- result (i, j) = xmin (a (i, j), b (i, j));
-
- return result;
-}
-
-static ComplexMatrix
-min (const ComplexMatrix& a, const ComplexMatrix& b)
-{
- int nr = a.rows ();
- int nc = a.columns ();
- if (nr != b.rows () || nc != b.columns ())
- {
- error ("two-arg min expecting args of same size");
- return ComplexMatrix ();
- }
-
- ComplexMatrix result (nr, nc);
-
- for (int j = 0; j < nc; j++)
- {
- int columns_are_real_only = 1;
- for (int i = 0; i < nr; i++)
- if (imag (a (i, j)) != 0.0 || imag (b (i, j)) != 0.0)
- {
- columns_are_real_only = 0;
- break;
- }
-
- if (columns_are_real_only)
- {
- for (int i = 0; i < nr; i++)
- result (i, j) = xmin (real (a (i, j)), real (b (i, j)));
- }
- else
- {
- for (int i = 0; i < nr; i++)
- result (i, j) = xmin (a (i, j), b (i, j));
- }
- }
-
- return result;
-}
-
-static Matrix
-max (double d, const Matrix& m)
-{
- int nr = m.rows ();
- int nc = m.columns ();
-
- Matrix result (nr, nc);
-
- for (int j = 0; j < nc; j++)
- for (int i = 0; i < nr; i++)
- result (i, j) = xmax (d, m (i, j));
+MM_COMP(min, ComplexMatrix)
+MM_COMP(max, ComplexMatrix)
- return result;
-}
-
-static Matrix
-max (const Matrix& m, double d)
-{
- int nr = m.rows ();
- int nc = m.columns ();
-
- Matrix result (nr, nc);
-
- for (int j = 0; j < nc; j++)
- for (int i = 0; i < nr; i++)
- result (i, j) = xmax (m (i, j), d);
-
- return result;
-}
-
-static ComplexMatrix
-max (const Complex& c, const ComplexMatrix& m)
-{
- int nr = m.rows ();
- int nc = m.columns ();
-
- ComplexMatrix result (nr, nc);
-
- for (int j = 0; j < nc; j++)
- for (int i = 0; i < nr; i++)
- result (i, j) = xmax (c, m (i, j));
-
- return result;
-}
-
-static ComplexMatrix
-max (const ComplexMatrix& m, const Complex& c)
-{
- int nr = m.rows ();
- int nc = m.columns ();
-
- ComplexMatrix result (nr, nc);
-
- for (int j = 0; j < nc; j++)
- for (int i = 0; i < nr; i++)
- result (i, j) = xmax (m (i, j), c);
-
- return result;
-}
-
-static Matrix
-max (const Matrix& a, const Matrix& b)
-{
- int nr = a.rows ();
- int nc = a.columns ();
- if (nr != b.rows () || nc != b.columns ())
- {
- error ("two-arg max expecting args of same size");
- return Matrix ();
- }
-
- Matrix result (nr, nc);
-
- for (int j = 0; j < nc; j++)
- for (int i = 0; i < nr; i++)
- result (i, j) = xmax (a (i, j), b (i, j));
-
- return result;
-}
-
-static ComplexMatrix
-max (const ComplexMatrix& a, const ComplexMatrix& b)
-{
- int nr = a.rows ();
- int nc = a.columns ();
- if (nr != b.rows () || nc != b.columns ())
- {
- error ("two-arg max expecting args of same size");
- return ComplexMatrix ();
- }
-
- ComplexMatrix result (nr, nc);
-
- for (int j = 0; j < nc; j++)
- {
- int columns_are_real_only = 1;
- for (int i = 0; i < nr; i++)
- if (imag (a (i, j)) != 0.0 || imag (b (i, j)) != 0.0)
- {
- columns_are_real_only = 0;
- break;
- }
-
- if (columns_are_real_only)
- {
- for (int i = 0; i < nr; i++)
- result (i, j) = xmax (real (a (i, j)), real (b (i, j)));
- }
- else
- {
- for (int i = 0; i < nr; i++)
- result (i, j) = xmax (a (i, j), b (i, j));
- }
- }
-
- return result;
-}
+#undef MM_COMP
+// interpreter routines
+
DEFUN_DLD (min, args, nargout,
"-*- texinfo -*-\n\
For a vector argument, return the minimum value. For a matrix\n\
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [patch] Re-arrange matrix min/max reductions.,
Edward Jason Riedy <=