/* vim: set expandtab shiftwidth=2 softtabstop=2 tw=70: */

#include <Rcpp.h>
using namespace Rcpp;

// Cross-reference work:
// 1. update ../src/registerDynamicSymbol.c with an item for this
// 2. main code should use the autogenerated wrapper in ../R/RcppExports.R
//
// [[Rcpp::export]]
List do_gradient(NumericMatrix m, NumericVector x, NumericVector y)
{
  int nx = x.size();
  int ny = y.size();
  if (ny < 3)
    ::Rf_error("cannot handle case with ny < 3 (FIXME)");
  if (nx < 3)
    ::Rf_error("cannot handle case with nx < 3 (FIXME)");
  if (m.nrow() != nx)
    ::Rf_error("matrix has %d rows, but length(x) is %d", m.nrow(), nx);
  if (m.ncol() != ny)
    ::Rf_error("matrix has %d cols, but length(y) is %d", m.ncol(), ny);
  NumericMatrix gx(nx, ny);
  for (int j = 0; j < ny; j++) {
    gx(     0, j) = (m(     1, j) - m(     0, j)) / (x[   1] - x[   0]);
    gx(nx - 1, j) = (m(nx - 1, j) - m(nx - 2, j)) / (x[nx-1] - x[nx-2]);
  }
  for (int i = 1; i < nx - 1; i++) {
    for (int j = 0; j < ny; j++) {
      gx(i, j) = (m(i+1, j) - m(i-1,j)) / (x[i+1] - x[i-1]);
    }
  }
  NumericMatrix gy(nx, ny);
  for (int i = 0; i < nx; i++) {
    gy(i,      0) = (m(i,      1) - m(i,      0)) / (y[   1] - y[   0]);
    gy(i, ny - 1) = (m(i, ny - 1) - m(i, ny - 2)) / (y[ny-1] - y[ny-2]);
  }
  for (int i = 0; i < nx; i++) {
    for (int j = 1; j < ny - 1; j++) {
      gy(i, j) = (m(i, j+1) - m(i,j-1)) / (y[j+1] - y[j-1]);
    }
  }
  return(List::create(Named("gx")=gx, Named("gy")=gy));
}

