We replaced the FORTRAN implementation for BEA with code from package…

… TSP. ME is now calculated using C code.

NEWS.md

@@ -5,6 +5,7 @@
 
 ## Changes
 - We replaced the FORTRAN implementation for BEA with code from package TSP.
+- ME is now calculated using C code.
 
 ## Bug Fixes
 - Added two missing package anchors to palette man page.

R/criterion.matrix.R

@@ -38,17 +38,32 @@ criterion.table <- criterion.matrix
 criterion_ME <- function(x, order = NULL, ...) {
   # ... unused
 
-  if (any(x < 0)) {
-    warning("Bond energy (ME) is only defined for nonnegative matrices. Returning NA.")
+  if (!is.matrix(x))
+    stop("Argument 'x' must be a matrix.")
+  if (!is.double(x))
+    mode(x) <- "double"
+
+  if (any(x < 0) || any(is.infinite(x)) || any(is.na(x))) {
+    warning("Bond energy (ME) is only defined for nonnegative finite matrices. Returning NA.")
     return(NA_real_)
   }
 
-  ### TODO: This could be done with less memory utilization in C
-  if (!is.null(order)) x <- permute(x, order)
-  .5 * sum(x * (rbind(0, x[-nrow(x), , drop = FALSE]) +
-    rbind(x[-1L, , drop = FALSE], 0) +
-    cbind(0, x[, -ncol(x), drop = FALSE]) +
-    cbind(x[, -1L , drop = FALSE], 0)))
+  if (is.null(order)) {
+    rows <- seq(dim(x)[1])
+    cols <- seq(dim(x)[2])
+  } else{
+    rows <- get_order(order, 1)
+    cols <- get_order(order, 2)
+  }
+
+  .Call("measure_of_effectiveness", x, rows, cols)
+
+  ### The R version needs lots of memory
+  #if (!is.null(order)) x <- permute(x, order)
+  #.5 * sum(x * (rbind(0, x[-nrow(x), , drop = FALSE]) +
+  #  rbind(x[-1L, , drop = FALSE], 0) +
+  #  cbind(0, x[, -ncol(x), drop = FALSE]) +
+  #  cbind(x[, -1L , drop = FALSE], 0)))
 }
 
 
@@ -67,8 +82,8 @@ criterion_ME <- function(x, order = NULL, ...) {
     mode(x) <- "double"
 
   if (is.null(order)) {
-    rows <- as.integer(1:dim(x)[1])
-    cols <- as.integer(1:dim(x)[2])
+    rows <- seq(dim(x)[1])
+    cols <- seq(dim(x)[2])
   } else{
     rows <- get_order(order, 1)
     cols <- get_order(order, 2)
@@ -85,6 +100,7 @@ criterion_ME <- function(x, order = NULL, ...) {
 criterion_stress_moore <-
   function(x, order, ...)
     .stress(x, order, "moore")
+
 criterion_stress_neumann <-
   function(x, order, ...)
     .stress(x, order, "neumann")

src/criterion.c

@@ -371,3 +371,49 @@ SEXP bar(SEXP R_dist, SEXP R_order, SEXP R_b) {
 
   return R_out;
 }
+
+
+/*
+ * Measure of effectiveness ME (McCormick et al, 1972)
+ */
+SEXP measure_of_effectiveness(SEXP R_mat, SEXP R_order_row, SEXP R_order_col) {
+
+  int *o_row = INTEGER(R_order_row);
+  int *o_col = INTEGER(R_order_col);
+  double *x = REAL(R_mat);
+  int nrow = INTEGER(getAttrib(R_mat, install("dim")))[0];
+  int ncol = INTEGER(getAttrib(R_mat, install("dim")))[1];
+
+  SEXP R_out;
+  double m = 0;
+  double s;
+  int i, j, ii, jj;
+
+  if (nrow != LENGTH(R_order_row) || ncol!= LENGTH(R_order_col))
+    error("dimenstions of matrix and order do not match!");
+
+  for (i = 0; i < nrow; ++i) {
+    for (j = 0; j < ncol; ++j) {
+      ii = o_row[i] - 1;
+      jj = o_col[j] - 1;
+
+      s = 0;
+      if(i > 0)  s += x[M_POS(nrow, o_row[i - 1] - 1, jj)];
+      if(i < (nrow - 1))  s += x[M_POS(nrow, o_row[i + 1] - 1, jj)];
+      if(j > 0)  s += x[M_POS(nrow, ii, o_col[j - 1] - 1)];
+      if(j < (ncol - 1))  s += x[M_POS(nrow, ii, o_col[j + 1] - 1)];
+
+      m += x[M_POS(nrow, ii, jj)] * s;
+    }
+  }
+
+  m = .5 * m;
+
+  // create R object
+  PROTECT(R_out = NEW_NUMERIC(1));
+  REAL(R_out)[0] = m;
+  UNPROTECT(1);
+
+  return R_out;
+}
+

src/dist.c

@@ -1,6 +1,6 @@
 /*
  * seriation - Infrastructure for seriation
- * Copyrigth (C) 2011 Michael Hahsler, Christian Buchta and Kurt Hornik
+ * Copyright (C) 2011 Michael Hahsler, Christian Buchta and Kurt Hornik
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by

src/init.c

@@ -10,7 +10,7 @@ extern SEXP gradient(SEXP, SEXP, SEXP);
 extern SEXP inertia_criterion(SEXP, SEXP);
 extern SEXP lazy_path_length(SEXP, SEXP);
 extern SEXP least_squares_criterion(SEXP, SEXP);
-//extern SEXP order_greedy(SEXP);
+extern SEXP measure_of_effectiveness(SEXP, SEXP, SEXP);
 extern SEXP order_length(SEXP, SEXP);
 extern SEXP order_optimal(SEXP, SEXP);
 extern SEXP pathdist_floyd(SEXP);
@@ -25,9 +25,6 @@ extern void permNext(void *, void *);
 extern void F77_NAME(arsa)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
 extern void F77_NAME(bburcg)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
 extern void F77_NAME(bbwrcg)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
-//extern void F77_NAME(cbea)(void *, void *, void *, void *, void *, void *, void *);
-//extern void F77_NAME(rbea)(void *, void *, void *, void *, void *, void *, void *);
-extern void F77_NAME(energy)(void *, void *, void *, void *);
 
 static const R_CallMethodDef CallEntries[] = {
     {"ar",                      (DL_FUNC) &ar,                      3},
@@ -36,7 +33,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"inertia_criterion",       (DL_FUNC) &inertia_criterion,       2},
     {"lazy_path_length",        (DL_FUNC) &lazy_path_length,        2},
     {"least_squares_criterion", (DL_FUNC) &least_squares_criterion, 2},
-//    {"order_greedy",            (DL_FUNC) &order_greedy,            1},
+    {"measure_of_effectiveness",(DL_FUNC) &measure_of_effectiveness,3},
     {"order_length",            (DL_FUNC) &order_length,            2},
     {"order_optimal",           (DL_FUNC) &order_optimal,           2},
     {"pathdist_floyd",          (DL_FUNC) &pathdist_floyd,          1},
@@ -57,9 +54,6 @@ static const R_FortranMethodDef FortranEntries[] = {
     {"arsa",   (DL_FUNC) &F77_NAME(arsa),   15},
     {"bburcg", (DL_FUNC) &F77_NAME(bburcg), 10},
     {"bbwrcg", (DL_FUNC) &F77_NAME(bbwrcg), 10},
-//  {"cbea",   (DL_FUNC) &F77_NAME(cbea),    7},
-//  {"rbea",   (DL_FUNC) &F77_NAME(rbea),    7},
-//  {"energy", (DL_FUNC) &F77_NAME(energy),  4},
     {NULL, NULL, 0}
 };
 

src/lt.h

@@ -13,7 +13,9 @@
 #endif
 
 
-/* M_POS to access matrix column-major order by i and j index (starts with 1) */
+/* M_POS to access matrix column-major order by i and j index (starts with 1)
+ * n is the number of rows
+ */
 
 #ifndef M_POS
 #define M_POS(n, i, j) ((i)+(R_xlen_t)(n)*(j))

src/pathdist.c

@@ -1,6 +1,6 @@
 /*
  * seriation - Infrastructure for seriation
- * Copyrigth (C) 2011 Michael Hahsler, Christian Buchta and Kurt Hornik
+ * Copyright (C) 2011 Michael Hahsler, Christian Buchta and Kurt Hornik
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by

src/stress.c

@@ -1,6 +1,6 @@
 /*
  * seriation - Infrastructure for seriation
- * Copyrigth (C) 2011 Michael Hahsler, Christian Buchta and Kurt Hornik
+ * Copyright (C) 2011 Michael Hahsler, Christian Buchta and Kurt Hornik
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by