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