From da35a9d973e6ce282822a45873b8492d4ef97ccf Mon Sep 17 00:00:00 2001 From: Michael Hahsler Date: Tue, 14 May 2024 14:36:08 -0500 Subject: [PATCH] We replaced the FORTRAN implementation for BEA with code from package TSP. ME is now calculated using C code. --- NEWS.md | 1 + R/criterion.matrix.R | 36 ++++++++++++++++++++++++---------- src/criterion.c | 46 ++++++++++++++++++++++++++++++++++++++++++++ src/dist.c | 2 +- src/init.c | 10 ++-------- src/lt.h | 4 +++- src/pathdist.c | 2 +- src/stress.c | 2 +- 8 files changed, 81 insertions(+), 22 deletions(-) diff --git a/NEWS.md b/NEWS.md index 60ae0fa..344a071 100644 --- a/NEWS.md +++ b/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. diff --git a/R/criterion.matrix.R b/R/criterion.matrix.R index 2db9f30..88ec3c5 100644 --- a/R/criterion.matrix.R +++ b/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") diff --git a/src/criterion.c b/src/criterion.c index 9a173e4..bdb3093 100644 --- a/src/criterion.c +++ b/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; +} + diff --git a/src/dist.c b/src/dist.c index ea897f6..f973ec4 100644 --- a/src/dist.c +++ b/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 diff --git a/src/init.c b/src/init.c index 74d2011..6d1e938 100644 --- a/src/init.c +++ b/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} }; diff --git a/src/lt.h b/src/lt.h index abc9ae9..a47d63a 100644 --- a/src/lt.h +++ b/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)) diff --git a/src/pathdist.c b/src/pathdist.c index 0552d49..ad2b213 100644 --- a/src/pathdist.c +++ b/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 diff --git a/src/stress.c b/src/stress.c index 1b378a2..e6c6287 100644 --- a/src/stress.c +++ b/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