algorithm.h

Go to the documentation of this file.

/* 
 * Scythe Statistical Library
 * Copyright (C) 2000-2002 Andrew D. Martin and Kevin M. Quinn;
 * 2002-present Andrew D. Martin, Kevin M. Quinn, and Daniel
 * Pemstein.  All Rights Reserved.
 *
 * This program is free software; you can redistribute it and/or modify
 * under the terms of the GNU General Public License as published by
 * Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.  See the text files COPYING
 * and LICENSE, distributed with this source code, for further
 * information.
 * --------------------------------------------------------------------
 * scythestat/algorithm.h
 */

#ifndef SCYTHE_ALGORITHM_H
#define SCYTHE_ALGORITHM_H

#include <cmath>
#include <functional>
#include <algorithm>

#ifdef SCYTHE_COMPILE_DIRECT
#include "defs.h"
#include "matrix.h"
#include "matrix_random_access_iterator.h"
#else
#include "scythestat/defs.h"
#include "scythestat/matrix.h"
#include "scythestat/matrix_random_access_iterator.h"
#endif

// These are just goofy

#ifdef SCYTHE_RPACK
#undef DO
#undef DS
#undef SO
#undef SS
#endif

namespace scythe {
  namespace {
    typedef unsigned int uint;
  }

  /* Matrix forward declaration */
  template <typename T_type, matrix_order ORDER, matrix_style STYLE>
  class Matrix;

  template <typename T>
  struct exponentiate : std::binary_function<T, T, T>
  {
    T operator() (T base, T exp) const
    {
      return std::pow(base, exp);
    }
  };

  template <typename T>
  struct ax_plus_b : std::binary_function<T,T,T>
  {
    T a_;
    ax_plus_b (T a) : a_ (a) {}
    T operator() (T x, T b) const
    {
      return (a_ * x + b);
    }
  };

  template <typename T, matrix_order O, matrix_style S, class FUNCTOR>
  void 
  for_each_ij_set (Matrix<T,O,S>& M, FUNCTOR func)
  {
    if (O == Col) {
      for (uint j = 0; j < M.cols(); ++j)
        for (uint i = 0; i < M.rows(); ++i)
          M(i, j) = func(i, j);
    } else {
      for (uint i = 0; i < M.cols(); ++i)
        for (uint j = 0; j < M.rows(); ++j)
          M(i, j) = func(i, j);
    }
  }

  template <matrix_order ORDER1, matrix_order ORDER2,
            typename T, typename S, matrix_order SO, matrix_style SS,
            matrix_order DO, matrix_style DS>
  void 
  copy(const Matrix<T,SO,SS>& source, Matrix<S,DO,DS>& dest)
  {
    std::copy(source.template begin_f<ORDER1>(), 
              source.template end_f<ORDER1>(),
              dest.template begin_f<ORDER2>());
  }

  template <matrix_order ORDER1, matrix_order ORDER2,
            typename T, matrix_order SO, matrix_style SS,
            matrix_order DO, matrix_style DS>
  void 
  copy_recycle (const Matrix<T,SO,SS>& source, Matrix<T,DO,DS>& dest)
  {
    if (source.size() == dest.size()) {
      copy<ORDER1,ORDER2> (source, dest);
    } else if (source.size() > dest.size()) {
      const_matrix_random_access_iterator<T,ORDER1,SO,SS> s_iter 
        = source.template begin<ORDER1>();
      std::copy(s_iter, s_iter + dest.size(),
                dest.template begin_f<ORDER2>());
    } else {
      const_matrix_random_access_iterator<T,ORDER1,SO,SS> s_begin
        = source.template begin<ORDER1> ();
      matrix_random_access_iterator<T,ORDER2,DO,DS> d_iter 
        = dest.template begin<ORDER2>();
      matrix_random_access_iterator<T,ORDER2,DO,DS> d_end
        = dest.template end<ORDER2>();
      while (d_iter != d_end) {
        unsigned int span = std::min(source.size(), 
            (unsigned int) (d_end - d_iter));
        d_iter = std::copy(s_begin, s_begin + span, d_iter);
      }
    }
  }

  template <class T>
  inline T sgn (const T & x)
  {
    if (x > (T) 0)
      return (T) 1;
    else if (x < (T) 0)
      return (T) -1;
    else
      return (T) 0;
  }

}  // end namespace scythe

#endif /* SCYTHE_ALGORITHM_H */