void Min::backward(const tensor_type& data_input, const tensor_type& gradient_output)
    {
      gradient_input.resizeLike(data_input);
      gradient_input.setZero();
      
      if (dimension) {
	for (int row = 0; row != data_input.rows(); ++ row)
	  gradient_input.row(row)[indices[row]] = gradient_output.col(0)[row];
      } else {
	for (int col = 0; col != data_input.cols(); ++ col)
	  gradient_input.col(col)[indices[col]] = gradient_output.col(0)[col];
      }
    }

    static void apply( const tensor_type & tensor ,
                       const MatrixValue * const a ,
                       const VectorValue * const x ,
                       VectorValue * const y )
    {

        const size_type nDim = tensor.dimension();

        // Loop over i
        for ( size_type i = 0; i < nDim; ++i) {
            VectorValue ytmp = 0;

            // Loop over k for this i
            const size_type nk = tensor.num_k(i);
            const size_type kBeg = tensor.k_begin(i);
            const size_type kEnd = kBeg + nk;
            for (size_type kEntry = kBeg; kEntry < kEnd; ++kEntry) {
                const size_type k = tensor.k_coord(kEntry);
                const MatrixValue ak = a[k];
                const VectorValue xk = x[k];

                // Loop over j for this i,k
                const size_type nj = tensor.num_j(kEntry);
                const size_type jBeg = tensor.j_begin(kEntry);
                const size_type jEnd = jBeg + nj;
                for (size_type jEntry = jBeg; jEntry < jEnd; ++jEntry) {
                    const size_type j = tensor.j_coord(jEntry);
                    ytmp += tensor.value(jEntry) * ( a[j] * xk + ak * x[j] );
                }
            }

            y[i] += ytmp ;
        }
    }

  KOKKOS_INLINE_FUNCTION
  static void apply( const tensor_type & tensor ,
                     const MatrixValue * const a ,
                     const VectorValue * const x ,
                           VectorValue * const y )
  {
    const size_type nk = tensor.num_k();

    // Loop over k
    for ( size_type k = 0; k < nk; ++k) {
      const MatrixValue ak = a[k];
      const VectorValue xk = x[k];

      // Loop over j for this k
      const size_type nj = tensor.num_j(k);
      const size_type jBeg = tensor.j_begin(k);
      const size_type jEnd = jBeg + nj;
      for (size_type jEntry = jBeg; jEntry < jEnd; ++jEntry) {
        const size_type j = tensor.j_coord(jEntry);
        VectorValue tmp = a[j] * xk + ak * x[j];

        // Loop over i for this k,j
        const size_type ni = tensor.num_i(jEntry);
        const size_type iBeg = tensor.i_begin(jEntry);
        const size_type iEnd = iBeg + ni;
        for (size_type iEntry = iBeg; iEntry < iEnd; ++iEntry) {
          const size_type i = tensor.i_coord(iEntry);
          y[i] += tensor.value(iEntry) * tmp;
        }
      }
    }
  }

void fill_from_file(tensor_type& tens, const std::string& str)
{
   std::ifstream ifs(str);
   for(int k = 0 ; k < tens.extent<2>(); ++k)
      for(int j = 0 ; j < tens.extent<1>(); ++j)
         for(int i = 0 ; i < tens.extent<0>(); ++i)
      {
         ifs >> tens.at(i,j,k);
      }
}

 void operator ()(const tensor_type &delta, const tensor_type &x, tensor_type &x_gradient) {
   CHECK_EQ(x.order(), delta.order());
   // TODO(robertsdionne): Support arbitrary tensor order with an n-dimensional iterator.
   for (auto i = 0; i < x.shape().at(0); ++i) {
     for (auto j = 0; j < x.shape().at(1); ++j) {
       x_gradient.set({i, j}, delta.at({i, j}) * (x.at({i, j}) > F(0)));
     }
   }
 }

    void SoftMax::forward(const tensor_type& data_input)
    {
      // use of redux for computing logsum...
      
      const double infty = - std::numeric_limits<double>::infinity();
      
      double logsum = infty;
      for (difference_type i = 0; i != data_input.rows(); ++ i) {
	const double value = data_input.col(0)[i];
	
	if (logsum == infty)
	  logsum = value;
	else if (value > infty) {
	  if (logsum >= value)
	    logsum = logsum + utils::mathop::log1p(std::exp(value - logsum));
	  else
	    logsum = value  + utils::mathop::log1p(std::exp(logsum - value));
	}
      }
      
      data_output = (data_input.array() - logsum).exp();
    }

void print(const tensor_type& tens, const std::string& str)
{
   std::ofstream of(str);
   of << std::showpos;
   using limit_type = std::numeric_limits<typename tensor_type::value_type>;
   of << std::scientific;
   of << std::setprecision(limit_type::max_digits10);
   for(int k = 0 ; k < tens.extent<2>(); ++k)
      for(int j = 0 ; j < tens.extent<1>(); ++j)
         for(int i = 0 ; i < tens.extent<0>(); ++i)
      {
         of << tens.at(i,j,k) << std::endl;
      }
}

 void operator ()(const tensor_type &x, tensor_type &y) {
   using std::max;
   CHECK_EQ(x.order(), y.order());
   // TODO(robertsdionne): Support arbitrary tensor order with an n-dimensional iterator.
   for (auto i = 0; i < x.shape().at(0); ++i) {
     for (auto j = 0; j < x.shape().at(1); ++j) {
       y.set({i, j}, max(F(0), x.at({i, j})));
     }
   }
 }

    void Min::forward(const tensor_type& data_input)
    {
      if (dimension) {
	// select row and compute column-min
	data_output.resize(data_input.rows(), 1);
	indices.resize(data_input.rows());
	
	for (size_type row = 0; row != data_input.rows(); ++ row) {
	  int col_min = 0;
	  data_output.col(0)[row] = data_input.row(row).minCoeff(&col_min);
	  indices[row] = col_min;
	}
      } else {
	// select column and compute row-min!
	data_output.resize(data_input.cols(), 1);
	indices.resize(data_input.cols());
	
	for (size_type col = 0; col != data_input.cols(); ++ col) {
	  int row_min = 0;
	  data_output.col(0)[col] = data_input.col(col).minCoeff(&row_min);
	  indices[col] = row_min;
	}
      }
    }

 static size_type matrix_size( const tensor_type & tensor )
 {
     return tensor.dimension();
 }

 KOKKOS_INLINE_FUNCTION
 static size_type vector_size( const tensor_type & tensor )
 { return tensor.dimension(); }

 void Exp::forward(const tensor_type& data_input)
 {
   data_output = data_input.array().exp();
 }

 void Exp::backward(const tensor_type& data_input, const tensor_type& gradient_output)
 {
   gradient_input.array() = gradient_output.array() * data_output.array();
 }

 void Log::backward(const tensor_type& data_input, const tensor_type& gradient_output)
 {
   gradient_input = data_input.array() / gradient_output.array();
 }

 static size_type vector_size( const tensor_type & tensor )
 {
     return tensor.dimension();
 }

  void operator ()(const tensor_type &x, tensor_type &y) {
    // y = F ∗ x + b
    for (auto i = 0; i < y.shape().at(0); ++i) {
      for (auto j = 0; j < y.shape().at(1); ++j) {
        for (auto s = 0; s < filter_.shape().at(0); ++s) {
          for (auto t = 0; t < filter_.shape().at(1); ++t) {
            for (auto u = 0; u < y.shape().at(2); ++u) {
              F output_value = F(1) * y.at({i, j, u});
              for (auto v = 0; v < filter_.shape().at(3); ++v) {
                output_value += F(1) * filter_.at({s, t, u, v}) * x.at({i + s, j + t, v});
              }
              y.set({i, j, u}, output_value);
            }
          }
        }
      }
    }

    for (auto i = 0; i < y.shape().at(0); ++i) {
      for (auto j = 0; j < y.shape().at(1); ++j) {
        for (auto k = 0; k < y.shape().at(2); ++k) {
          y.set({i, j, k}, F(1) * y.at({i, j, k}) + F(1) * bias_.at({k}));
        }
      }
    }
  }

 static void rand_fill(tensor_type& tensor) {
   for(std::size_t i = 0ul; i < tensor.size(); ++i)
     tensor[i] = GlobalFixture::world->rand() % 27;
 }

 void SoftMax::backward(const tensor_type& data_input, const tensor_type& gradient_output)
 {
   const double sum = (gradient_output.array() * data_output.array()).sum();
   
   gradient_input = (gradient_output.array() - sum) * data_output.array();
 }

 void Log::forward(const tensor_type& data_input)
 {
   data_output = data_input.array().log();
 }