1 /* 2 * Scythe Statistical Library 3 * Copyright (C) 2000-2002 Andrew D. Martin and Kevin M. Quinn; 4 * 2002-present Andrew D. Martin, Kevin M. Quinn, and Daniel 5 * Pemstein. All Rights Reserved. 6 * 7 * This program is free software; you can redistribute it and/or modify 8 * under the terms of the GNU General Public License as published by 9 * Free Software Foundation; either version 2 of the License, or (at 10 * your option) any later version. See the text files COPYING 11 * and LICENSE, distributed with this source code, for further 12 * information. 13 * -------------------------------------------------------------------- 14 * scythestat/algorithm.h 15 */ 16 17 /*! \file algorithm.h 18 * 19 * \brief Generic algorithms for Scythe objects. 20 * 21 * This file provides implementations of a few algorithms that operate 22 * on Scythe objects and also contains the definitions of a handful of 23 * useful function objects. These functions and functors are primarily 24 * intended for use within the library. We add algorithms to this 25 * header as need arises and do not currently attempt to provide a 26 * comprehensive set of generic algorithms for working with Scythe 27 * matrices. 28 * 29 */ 30 31 #ifndef SCYTHE_ALGORITHM_H 32 #define SCYTHE_ALGORITHM_H 33 34 #include <cmath> 35 #include <functional> 36 #include <algorithm> 37 38 #ifdef SCYTHE_COMPILE_DIRECT 39 #include "defs.h" 40 #include "matrix.h" 41 #include "matrix_random_access_iterator.h" 42 #else 43 #include "scythestat/defs.h" 44 #include "scythestat/matrix.h" 45 #include "scythestat/matrix_random_access_iterator.h" 46 #endif 47 48 // These are just goofy 49 50 #ifdef SCYTHE_RPACK 51 #undef DO 52 #undef DS 53 #undef SO 54 #undef SS 55 #endif 56 57 namespace scythe { 58 namespace { 59 typedef unsigned int uint; 60 } 61 62 /* Matrix forward declaration */ 63 template <typename T_type, matrix_order ORDER, matrix_style STYLE> 64 class Matrix; 65 66 /*! \brief A Functor encapsulating exponentiation. 67 * 68 * This function object wraps exponentiation operations for use in 69 * generic algorithms. 70 */ 71 template <typename T> 72 struct exponentiate : std::binary_function<T, T, T> 73 { operatorexponentiate74 T operator() (T base, T exp) const 75 { 76 return std::pow(base, exp); 77 } 78 }; 79 80 /*! \brief A Functor encapsulating \f$ax+b\f$. 81 * 82 * This function object wraps the operation \f$ax+b\f$ for use in 83 * generic algorithms, where a is some constant. 84 */ 85 template <typename T> 86 struct ax_plus_b : std::binary_function<T,T,T> 87 { 88 T a_; ax_plus_bax_plus_b89 ax_plus_b (T a) : a_ (a) {} operatorax_plus_b90 T operator() (T x, T b) const 91 { 92 return (a_ * x + b); 93 } 94 }; 95 96 /*! \brief Iterate through a Matrix in order. 97 * 98 * This function iterates through a Matrix, \a M, in order, 99 * setting each element in the Matrix to the result of an invocation 100 * of the function object, \a func. The () operator of \a func 101 * should take two unsigned integer parameters (i - the row offset 102 * into \a M; j - the column offset into \a M) and return a result 103 * of type T. 104 * 105 * \param M The Matrix to iterate over. 106 * \param func The functor to execute on each iteration. 107 * 108 */ 109 110 template <typename T, matrix_order O, matrix_style S, class FUNCTOR> 111 void for_each_ij_set(Matrix<T,O,S> & M,FUNCTOR func)112 for_each_ij_set (Matrix<T,O,S>& M, FUNCTOR func) 113 { 114 if (O == Col) { 115 for (uint j = 0; j < M.cols(); ++j) 116 for (uint i = 0; i < M.rows(); ++i) 117 M(i, j) = func(i, j); 118 } else { 119 for (uint i = 0; i < M.cols(); ++i) 120 for (uint j = 0; j < M.rows(); ++j) 121 M(i, j) = func(i, j); 122 } 123 } 124 125 /*! \brief Copy the contents of one Matrix into another. 126 * 127 * This function copies the contents of one Matrix into 128 * another, traversing each Matrix in the order specified by the 129 * template terms ORDER1 and ORDER2. This function requires an 130 * explicit template call that specifies ORDER1 and ORDER2. 131 * 132 * \param source The Matrix to copy. 133 * \param dest The Matrix to copy into. 134 */ 135 136 template <matrix_order ORDER1, matrix_order ORDER2, 137 typename T, typename S, matrix_order SO, matrix_style SS, 138 matrix_order DO, matrix_style DS> 139 void copy(const Matrix<T,SO,SS> & source,Matrix<S,DO,DS> & dest)140 copy(const Matrix<T,SO,SS>& source, Matrix<S,DO,DS>& dest) 141 { 142 std::copy(source.template begin_f<ORDER1>(), 143 source.template end_f<ORDER1>(), 144 dest.template begin_f<ORDER2>()); 145 } 146 147 /*! \brief Copy the contents of one Matrix into another. 148 * 149 * This function copies the contents of one Matrix into 150 * another, traversing each Matrix in the order specified by the 151 * template terms ORDER1 and ORDER2. If \a source is larger than \a 152 * dest, the function only copies as many elements from \a source as 153 * will fit in \a dest. On the other hand, if \a source is smaller 154 * than \a dest, the function will start over at the beginning of 155 * \a source, recycling the contents of \a source as many times as 156 * necessary to fill \a dest. This function requires an explicit 157 * template call that specifies ORDER1 and ORDER2. 158 * 159 * \param source The Matrix to copy. 160 * \param dest The Matrix to copy into. 161 */ 162 template <matrix_order ORDER1, matrix_order ORDER2, 163 typename T, matrix_order SO, matrix_style SS, 164 matrix_order DO, matrix_style DS> 165 void copy_recycle(const Matrix<T,SO,SS> & source,Matrix<T,DO,DS> & dest)166 copy_recycle (const Matrix<T,SO,SS>& source, Matrix<T,DO,DS>& dest) 167 { 168 if (source.size() == dest.size()) { 169 copy<ORDER1,ORDER2> (source, dest); 170 } else if (source.size() > dest.size()) { 171 const_matrix_random_access_iterator<T,ORDER1,SO,SS> s_iter 172 = source.template begin<ORDER1>(); 173 std::copy(s_iter, s_iter + dest.size(), 174 dest.template begin_f<ORDER2>()); 175 } else { 176 const_matrix_random_access_iterator<T,ORDER1,SO,SS> s_begin 177 = source.template begin<ORDER1> (); 178 matrix_random_access_iterator<T,ORDER2,DO,DS> d_iter 179 = dest.template begin<ORDER2>(); 180 matrix_random_access_iterator<T,ORDER2,DO,DS> d_end 181 = dest.template end<ORDER2>(); 182 while (d_iter != d_end) { 183 unsigned int span = std::min(source.size(), 184 (unsigned int) (d_end - d_iter)); 185 d_iter = std::copy(s_begin, s_begin + span, d_iter); 186 } 187 } 188 } 189 190 /*! \brief Determine the sign of a number. 191 * 192 * This function compares \a x to (T) 0, returning (T) 1 if \a x is 193 * greater than zero, (T) -1 if \a x is less than zero, and (T) 0 194 * otherwise. 195 * 196 * \param x The value to check. 197 */ 198 template <class T> sgn(const T & x)199 inline T sgn (const T & x) 200 { 201 if (x > (T) 0) 202 return (T) 1; 203 else if (x < (T) 0) 204 return (T) -1; 205 else 206 return (T) 0; 207 } 208 209 } // end namespace scythe 210 211 #endif /* SCYTHE_ALGORITHM_H */ 212