00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00031 #ifndef SCYTHE_ALGORITHM_H
00032 #define SCYTHE_ALGORITHM_H
00033
00034 #include <cmath>
00035 #include <functional>
00036 #include <algorithm>
00037
00038 #ifdef SCYTHE_COMPILE_DIRECT
00039 #include "defs.h"
00040 #include "matrix.h"
00041 #include "matrix_random_access_iterator.h"
00042 #else
00043 #include "scythestat/defs.h"
00044 #include "scythestat/matrix.h"
00045 #include "scythestat/matrix_random_access_iterator.h"
00046 #endif
00047
00048 namespace scythe {
00049 namespace {
00050 typedef unsigned int uint;
00051 }
00052
00053
00054 template <typename T_type, matrix_order ORDER, matrix_style STYLE>
00055 class Matrix;
00056
00062 template <typename T>
00063 struct exponentiate : std::binary_function<T, T, T>
00064 {
00065 T operator() (T base, T exp) const
00066 {
00067 return std::pow(base, exp);
00068 }
00069 };
00070
00076 template <typename T>
00077 struct ax_plus_b : std::binary_function<T,T,T>
00078 {
00079 T a_;
00080 ax_plus_b (T a) : a_ (a) {}
00081 T operator() (T x, T b) const
00082 {
00083 return (a_ * x + b);
00084 }
00085 };
00086
00101 template <typename T, matrix_order O, matrix_style S, class FUNCTOR>
00102 void
00103 for_each_ij_set (Matrix<T,O,S>& M, FUNCTOR func)
00104 {
00105 if (O == Col) {
00106 for (uint j = 0; j < M.cols(); ++j)
00107 for (uint i = 0; i < M.rows(); ++i)
00108 M(i, j) = func(i, j);
00109 } else {
00110 for (uint i = 0; i < M.cols(); ++i)
00111 for (uint j = 0; j < M.rows(); ++j)
00112 M(i, j) = func(i, j);
00113 }
00114 }
00115
00127 template <matrix_order ORDER1, matrix_order ORDER2,
00128 typename T, typename S, matrix_order SO, matrix_style SS,
00129 matrix_order DO, matrix_style DS>
00130 void
00131 copy(const Matrix<T,SO,SS>& source, Matrix<S,DO,DS>& dest)
00132 {
00133 std::copy(source.template begin_f<ORDER1>(),
00134 source.template end_f<ORDER1>(),
00135 dest.template begin_f<ORDER2>());
00136 }
00137
00153 template <matrix_order ORDER1, matrix_order ORDER2,
00154 typename T, matrix_order SO, matrix_style SS,
00155 matrix_order DO, matrix_style DS>
00156 void
00157 copy_recycle (const Matrix<T,SO,SS>& source, Matrix<T,DO,DS>& dest)
00158 {
00159 if (source.size() == dest.size()) {
00160 copy<ORDER1,ORDER2> (source, dest);
00161 } else if (source.size() > dest.size()) {
00162 const_matrix_random_access_iterator<T,ORDER1,SO,SS> s_iter
00163 = source.template begin<ORDER1>();
00164 std::copy(s_iter, s_iter + dest.size(),
00165 dest.template begin_f<ORDER2>());
00166 } else {
00167 const_matrix_random_access_iterator<T,ORDER1,SO,SS> s_begin
00168 = source.template begin<ORDER1> ();
00169 matrix_random_access_iterator<T,ORDER2,DO,DS> d_iter
00170 = dest.template begin<ORDER2>();
00171 matrix_random_access_iterator<T,ORDER2,DO,DS> d_end
00172 = dest.template end<ORDER2>();
00173 while (d_iter != d_end) {
00174 unsigned int span = std::min(source.size(),
00175 (unsigned int) (d_end - d_iter));
00176 d_iter = std::copy(s_begin, s_begin + span, d_iter);
00177 }
00178 }
00179 }
00180
00189 template <class T>
00190 inline T sgn (const T & x)
00191 {
00192 if (x > (T) 0)
00193 return (T) 1;
00194 else if (x < (T) 0)
00195 return (T) -1;
00196 else
00197 return (T) 0;
00198 }
00199
00200 }
00201
00202 #endif