Недавно я проводил вычислительные испытания с mgcv
GAM. Некоторые оригинальные функции изменены, а некоторые добавлены. Чтобы не нарушать совместимость, для каждой функции, которую я хочу изменить, я создаю новую версию с .zheyuan
surfix в имени функции. Например, с функцией Sl.fit
для подгонки методом наименьших квадратов со штрафом у меня будет Sl.fit.zheyuan
. Я бы просто собрал все функции R, которые я пишу, в отдельный сценарий R "zheyuan.R". Добавив этот файл в каталог R
источника пакета mgcv_1.8-17
и скомпилировав этот измененный пакет в локальный путь, я мог загрузить его для тестирования.
У меня нет проблем с добавлением подпрограмм R, но не при добавлении подпрограмм C. При установке измененного пакета ошибок не возникает, но когда я вызываю функцию-оболочку R моей добавленной подпрограммы C, я получаю сообщение об ошибке, как в заголовке моего вопроса. Если вас интересует мой случай, вы можете выполнить следующие шаги, чтобы воспроизвести такую ошибку.
Шаг 1. Загрузите последний исходный код пакета
Загрузите версию 1.8-17 по указанной выше ссылке. Такая ссылка исчезнет, когда новый mgcv
релиз будет опубликован на CRAN. Но вы всегда можете перейти на mgcv
страницу CRAN, чтобы загрузить последний релиз.
Давайте untar
источник. Сначала удалите файл MD5
, чтобы мы не получали раздражающее предупреждение MD5 при компиляции измененной версии. В дальнейшем мы добавим новый материал в каталог R
и каталог src
.
Шаг 2. Создайте сценарий R
Рассмотрим следующую функцию-оболочку R:
RX <- function (R, X) {
X <- X + 0
.Call("C_mgcv_RX", R, X)
X
}
Создайте сценарий R "zheyuan.R" для размещения этой функции и добавьте его в mgcv/R
.
Шаг 3. Добавление подпрограмм C
Подпрограммы C, касающиеся матричных вычислений, обычно находятся в src/mat.c
. Итак, давайте добавим новую функцию в конец этого скрипта:
void mgcv_RX (SEXP R, SEXP X) {
int nrowX = nrows(X);
int ncolX = ncols(X);
double one = 1.0;
F77_CALL(dtrmm)("l", "u", "n", "n", &nrowX, &ncolX, &one, REAL(R), &nrowX, REAL(X), &nrowX);
}
Это простая процедура умножения верхней треугольной матрицы R
на прямоугольную матрицу X
. Выходная матрица перезапишет X
. Для этого будет вызван уровень 3 BLAS dtrmm
. Нам не нужно беспокоиться о файлах заголовков или связывании среды выполнения с библиотекой BLAS. Заголовки доступны в mat.c
, а связывание с BLAS управляется R.
Шаг 4. Зарегистрируйте процедуру C
Вышеупомянутого недостаточно. Каждая процедура C в mgcv
появится в трех местах. Например, давайте попробуем выполнить поиск в собственной подпрограмме C:
grep mgcv_RPPt mgcv/src/*
# mgcv/src/init.c: {"mgcv_RPPt",(DL_FUNC)&mgcv_RPPt,3},
# mgcv/src/mat.c:void mgcv_RPPt(SEXP a,SEXP r, SEXP NT) {
# mgcv/src/mgcv.h:void mgcv_RPPt(SEXP a,SEXP r, SEXP NT);
Нам также необходимо добавить файл заголовка mgcv.h
и зарегистрировать эту процедуру C в init.c
.
Добавим
void mgcv_RX (SEXP R, SEXP X);
до конца mgcv.h
и внутри init.c
выполните:
R_CallMethodDef CallMethods[] = {
{"mgcv_pmmult2", (DL_FUNC) &mgcv_pmmult2,5},
{"mgcv_Rpiqr", (DL_FUNC) &mgcv_Rpiqr,5},
{"mgcv_tmm",(DL_FUNC)&mgcv_tmm,5},
{"mgcv_Rpbsi",(DL_FUNC)&mgcv_Rpbsi,2},
{"mgcv_RPPt",(DL_FUNC)&mgcv_RPPt,3},
{"mgcv_Rpchol",(DL_FUNC)&mgcv_Rpchol,4},
{"mgcv_Rpforwardsolve",(DL_FUNC)&mgcv_Rpforwardsolve,3},
{"mgcv_Rpcross",(DL_FUNC)&mgcv_Rpcross,3},
{"mgcv_RX",(DL_FUNC)&mgcv_RX,2}, // we add this line
{NULL, NULL, 0}
};
Шаг 5. Выполните и загрузите
tar
измененная папка mgcv
на mgcv.tar.gz
.
Откройте новый чистый сеанс R (возможно, вам понадобится R --vanilla
для запуска). Затем укажите путь к локальной библиотеке и запустите:
path <- getwd() ## let's just use current working directory
## make sure you move "mgcv.tar.gz" into current working path
install.packages("mgcv.tar.gz", repos = NULL, lib = path)
library(mgcv, lib.loc = path)
Шаг 6: проверьте и получите ошибку
R <- matrix(runif(25), 5)
R[lower.tri(R)] <- 0
X <- matrix(runif(25), 5)
mgcv:::RX(R, X) ## function is not exported, so use `mgcv:::` to find it
# Error in .Call("C_mgcv_RX", R, X) :
# "C_mgcv_RX" not resolved from current namespace (mgcv)
Может ли кто-нибудь объяснить, почему и как это решить?