From octave-maintainers-request at bevo dot che dot wisc dot edu Sun Sep 3 23:25:55 2000 Subject: [PATCH] Add NaN-awareness to matrix min/max operators (2/4) From: Edward Jason Riedy To: octave-maintainers at bevo dot che dot wisc dot edu Date: Sun, 03 Sep 2000 21:25:54 -0700 The next in the series. I put the default OCTAVE_IGNORE_NAN_DEFAULT in mx-defs.h simply because I couldn't see a better place. Jason liboctave/ChangeLog * mx-reductions.h: New file. Helps define matrix reductions. * mx-defs.h (OCTAVE_IGNORE_NAN_DEFAULT): New define for ignore_nan default value. * dMatrix.h (Matrix::row_min, Matrix::row_max, Matrix::column_min, Matrix::column_max): add bool ignore_nan parameter. * 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. * oct-cmplx.h: Include (ComplexLess): std::binary_function implementing less-than over complex numbers by magnitude. (ComplexGreater): Analagous, for greater-than. --- octave.orig/liboctave/CMatrix.cc Mon Feb 7 20:35:41 2000 +++ octave.2/liboctave/CMatrix.cc Sat Sep 2 20:06:29 2000 at @ -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. at @ -2616,253 +2617,16 @@ return retval; } -ComplexColumnVector -ComplexMatrix::row_min (void) const -{ - Array index; - return row_min (index); -} - -ComplexColumnVector -ComplexMatrix::row_min (Array& 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 index; - return row_max (index); -} - -ComplexColumnVector -ComplexMatrix::row_max (Array& 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 index; - return column_min (index); -} - -ComplexRowVector -ComplexMatrix::column_min (Array& 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 index; - return column_max (index); -} - -ComplexRowVector -ComplexMatrix::column_max (Array& 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 --new-file -r -u octave.orig/liboctave/CMatrix.h octave.2/liboctave/CMatrix.h --- octave.orig/liboctave/CMatrix.h Mon Feb 7 20:35:41 2000 +++ octave.2/liboctave/CMatrix.h Sat Sep 2 20:15:09 2000 at @ -256,17 +256,29 @@ 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 = OCTAVE_IGNORE_NAN_DEFAULT) + const; + ComplexColumnVector row_max (bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) + const; - ComplexColumnVector row_min (Array& index) const; - ComplexColumnVector row_max (Array& index) const; + ComplexColumnVector row_min (Array& index, + bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) + const; + ComplexColumnVector row_max (Array& index, + bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) + const; - ComplexRowVector column_min (void) const; - ComplexRowVector column_max (void) const; + ComplexRowVector column_min (bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) + const; + ComplexRowVector column_max (bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) + const; - ComplexRowVector column_min (Array& index) const; - ComplexRowVector column_max (Array& index) const; + ComplexRowVector column_min (Array& index, + bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) + const; + ComplexRowVector column_max (Array& index, + bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) + const; // i/o diff --new-file -r -u octave.orig/liboctave/dMatrix.cc octave.2/liboctave/dMatrix.cc --- octave.orig/liboctave/dMatrix.cc Mon Feb 7 20:35:44 2000 +++ octave.2/liboctave/dMatrix.cc Sat Sep 2 20:06:29 2000 at @ -31,6 +31,7 @@ #include +#include #include #include "byte-swap.h" at @ -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. at @ -2173,225 +2176,16 @@ return d; } -ColumnVector -Matrix::row_min (void) const -{ - Array index; - return row_min (index); -} - -ColumnVector -Matrix::row_min (Array& 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 index; - return row_max (index); -} - -ColumnVector -Matrix::row_max (Array& 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 index; - return column_min (index); -} - -RowVector -Matrix::column_min (Array& 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 index; - return column_max (index); -} - -RowVector -Matrix::column_max (Array& 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) +DEF_COMP_OP(row, Matrix, Column, double, max, std::greater) +DEF_COMP_OP(column, Matrix, Row, double, min, std::less) +DEF_COMP_OP(column, Matrix, Row, double, max, std::greater) std::ostream& operator << (std::ostream& os, const Matrix& a) diff --new-file -r -u octave.orig/liboctave/dMatrix.h octave.2/liboctave/dMatrix.h --- octave.orig/liboctave/dMatrix.h Fri Jun 30 02:30:45 2000 +++ octave.2/liboctave/dMatrix.h Sat Sep 2 20:15:12 2000 at @ -210,17 +210,21 @@ 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 = OCTAVE_IGNORE_NAN_DEFAULT) const; + ColumnVector row_max (bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) const; - ColumnVector row_min (Array& index) const; - ColumnVector row_max (Array& index) const; + ColumnVector row_min (Array& index, + bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) const; + ColumnVector row_max (Array& index, + bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) const; - RowVector column_min (void) const; - RowVector column_max (void) const; + RowVector column_min (bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) const; + RowVector column_max (bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) const; - RowVector column_min (Array& index) const; - RowVector column_max (Array& index) const; + RowVector column_min (Array& index, + bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) const; + RowVector column_max (Array& index, + bool ignore_nan = OCTAVE_IGNORE_NAN_DEFAULT) const; // i/o diff --new-file -r -u octave.orig/liboctave/mx-defs.h octave.2/liboctave/mx-defs.h --- octave.orig/liboctave/mx-defs.h Mon Jan 31 20:06:17 2000 +++ octave.2/liboctave/mx-defs.h Sat Sep 2 20:15:27 2000 at @ -23,6 +23,10 @@ #if !defined (octave_mx_defs_h) #define octave_mx_defs_h 1 +// Behavioral defaults. + +#define OCTAVE_IGNORE_NAN_DEFAULT true + // Classes we declare. class Matrix; diff --new-file -r -u octave.orig/liboctave/mx-reductions.h octave.2/liboctave/mx-reductions.h --- octave.orig/liboctave/mx-reductions.h Wed Dec 31 16:00:00 1969 +++ octave.2/liboctave/mx-reductions.h Sat Sep 2 20:09:08 2000 at @ -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 +inline void +do_over_rows (const MArray2& 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 +inline void +do_over_columns (const MArray2& 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 +class PredReducer +{ +public: + inline PredReducer (Array& stored_val, Array& 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) || (xisnan (tmp_val) && !xisnan (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& stored_val; + Array& stored_index; + Pred pred; + + T tmp_val; + int tmp_idx; +}; + +template +class PredReducer_NaN +{ +public: + inline PredReducer_NaN (Array& stored_val, Array& 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& stored_val; + Array& 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 \ +static inline VectorType ## Vector \ +Side ## _reduction (const MatType& A, Array& 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 index; \ + return Side ## _ ## Op (index, ignore_nan); \ +} \ + \ +VectorType ## Vector \ +MatType:: ## Side ## _ ## Op (Array& index, bool ignore_nan) const \ +{ \ + if (ignore_nan) \ + return Side ## _reduction< PredReducer > \ + (*this, index); \ + else \ + return Side ## _reduction< PredReducer_NaN > \ + (*this, index); \ +} + +#endif diff --new-file -r -u octave.orig/liboctave/oct-cmplx.h octave.2/liboctave/oct-cmplx.h --- octave.orig/liboctave/oct-cmplx.h Tue Feb 1 02:07:22 2000 +++ octave.2/liboctave/oct-cmplx.h Sat Sep 2 20:06:29 2000 at @ -24,8 +24,29 @@ #define octave_oct_cmplx_h 1 #include +#include typedef std::complex Complex; + +struct ComplexLess + : public std::binary_function +{ + bool operator() (const Complex& x, const Complex& y) const + { + return abs(x) < abs(y); + } +}; + +struct ComplexGreater + : public std::binary_function +{ + bool operator() (const Complex& x, const Complex& y) const + { + return abs(x) > abs(y); + } +}; + + #endif