Вот как это сделать с помощью рекурсии:
permsum <- function(s,m) if (m==1L) matrix(s) else do.call(rbind,lapply(seq_len(s-m+1L),function(x) unname(cbind(x,permsum(s-x,m-1L)))));
res <- permsum(100L,4L);
head(res);
## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 97
## [2,] 1 1 2 96
## [3,] 1 1 3 95
## [4,] 1 1 4 94
## [5,] 1 1 5 93
## [6,] 1 1 6 92
tail(res);
## [,1] [,2] [,3] [,4]
## [156844,] 95 2 2 1
## [156845,] 95 3 1 1
## [156846,] 96 1 1 2
## [156847,] 96 1 2 1
## [156848,] 96 2 1 1
## [156849,] 97 1 1 1
Вы можете разделить на 100, чтобы получить дроби, а не целые числа:
head(res)/100;
## [,1] [,2] [,3] [,4]
## [1,] 0.01 0.01 0.01 0.97
## [2,] 0.01 0.01 0.02 0.96
## [3,] 0.01 0.01 0.03 0.95
## [4,] 0.01 0.01 0.04 0.94
## [5,] 0.01 0.01 0.05 0.93
## [6,] 0.01 0.01 0.06 0.92
Объяснение
Сначала давайте определим входы:
s Это целевое значение, к которому должна суммироваться каждая строка в выходной матрице.m Это количество столбцов, которые необходимо создать в выходной матрице.
Более эффективно и надежно вычислять результат, используя целочисленную арифметику, а не арифметику с плавающей запятой, поэтому я разработал свое решение для работы только с целыми числами. Следовательно, s является скалярным целым числом, представляющим целевую целочисленную сумму.
Теперь давайте рассмотрим последовательность, сгенерированную seq_len() для небазового случая:
seq_len(s-m+1L)
Это генерирует последовательность от 1 до максимально возможного значения, которое может быть частью суммы до s с оставшимися m столбцами. Например, подумайте о случае s=100,m=4: максимальное число, которое мы можем использовать, — это 97, участвующее в сумме 97+1+1+1. Каждый оставшийся столбец уменьшает максимально возможное значение на 1, поэтому мы должны вычесть m из s при вычислении длины последовательности.
Каждый элемент сгенерированной последовательности следует рассматривать как один из возможных «выборов» слагаемого при суммировании.
do.call(rbind,lapply(seq_len(s-m+1L),function(x) ...))
Затем для каждого выбора мы должны выполнить рекурсию. Мы можем использовать lapply() для этого.
Чтобы перейти вперед, лямбда сделает один рекурсивный вызов permsum(), а затем cbind() вернет значение с текущим выбором. Это создаст матрицу всегда одинаковой ширины для этого уровня рекурсии. Следовательно, вызов lapply() вернет список матриц одинаковой ширины. Затем мы должны связать их вместе, поэтому здесь мы должны использовать трюк do.call(rbind,...).
unname(cbind(x,permsum(s-x,m-1L)))
Тело лямбды довольно простое; мы cbind() текущий выбор x с возвращаемым значением рекурсивного вызова, завершая суммирование для этой подматрицы. К сожалению, мы должны вызывать unname(), иначе каждый столбец, который в конечном итоге будет установлен из аргумента x, будет иметь имя столбца x.
Самая важная деталь здесь — выбор аргументов для рекурсивного вызова. Во-первых, поскольку лямбда-аргумент x был только что выбран во время текущей рекурсивной оценки, мы должны вычесть его из s, чтобы получить новую цель суммирования, за достижение которой будет отвечать предстоящий рекурсивный вызов. Следовательно, первый аргумент становится s-x. Во-вторых, поскольку выбор x занимает один столбец, мы должны вычесть 1 из m, так что рекурсивный вызов будет иметь на один столбец меньше выходной матрицы.
if (m==1L) matrix(s) else ...
Наконец, давайте рассмотрим базовый случай. При каждом вычислении рекурсивной функции мы должны проверять, достигло ли m 1, и в этом случае мы можем просто вернуть требуемую сумму s.
Расхождение с плавающей запятой
Я изучил несоответствие между моими результатами и результатами psidom. Например:
library(data.table);
bgoldst <- function(s,m) permsum(s,m)/s;
psidom <- function(ss,m) { raw <- do.call(data.table::CJ,rep(list(ss),m)); raw[rowSums(raw)==1,]; };
## helper function to sort a matrix by columns
smp <- function(m) m[do.call(order,as.data.frame(m)),];
s <- 100L; m <- 3L; ss <- seq_len(s-1L)/s;
x <- smp(bgoldst(s,m));
y <- smp(unname(as.matrix(psidom(ss,m))));
nrow(x);
## [1] 4851
nrow(y);
## [1] 4809
Таким образом, между нашими двумя результатами есть расхождение в 42 строки. Я решил попытаться найти, какие именно перестановки были пропущены с помощью следующей строки кода. По сути, он сравнивает каждый элемент двух матриц и выводит результат сравнения в виде логической матрицы. Мы можем просмотреть прокрутку вниз, чтобы найти первую отличающуюся строку. Ниже приведен фрагмент вывода:
x==do.call(rbind,c(list(y),rep(list(NA),nrow(x)-nrow(y))));
## [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE
## [4,] TRUE TRUE TRUE
## [5,] TRUE TRUE TRUE
##
## ... snip ...
##
## [24,] TRUE TRUE TRUE
## [25,] TRUE TRUE TRUE
## [26,] TRUE TRUE TRUE
## [27,] TRUE TRUE TRUE
## [28,] TRUE TRUE TRUE
## [29,] TRUE FALSE FALSE
## [30,] TRUE FALSE FALSE
## [31,] TRUE FALSE FALSE
## [32,] TRUE FALSE FALSE
## [33,] TRUE FALSE FALSE
##
## ... snip ...
Итак, в строке 29 у нас есть первое несоответствие. Вот окно вокруг этой строки в каждой матрице перестановок:
win <- 27:31;
x[win,]; y[win,];
## [,1] [,2] [,3]
## [1,] 0.01 0.27 0.72
## [2,] 0.01 0.28 0.71
## [3,] 0.01 0.29 0.70 (missing from y)
## [4,] 0.01 0.30 0.69 (missing from y)
## [5,] 0.01 0.31 0.68
## [,1] [,2] [,3]
## [1,] 0.01 0.27 0.72
## [2,] 0.01 0.28 0.71
## [3,] 0.01 0.31 0.68
## [4,] 0.01 0.32 0.67
## [5,] 0.01 0.33 0.66
Интересно, что отсутствующие перестановки обычно дают в сумме ровно 1, когда вы вычисляете сумму вручную. Сначала я подумал, что это функция CJ() data.table делает что-то странное с числами с плавающей запятой, но дальнейшее тестирование, похоже, указывает на то, что это делает rowSums():
0.01+0.29+0.70==1;
## [1] TRUE
ss[1L]+ss[29L]+ss[70L]==1;
## [1] TRUE
rowSums(CJ(0.01,0.29,0.70))==1; ## looks like CJ()'s fault, but wait...
## [1] FALSE
cj <- CJ(0.01,0.29,0.70);
cj$V1+cj$V2+cj$V3==1; ## not CJ()'s fault
## [1] TRUE
rowSums(matrix(c(0.01,0.29,0.70),1L,byrow=T))==1; ## rowSums()'s fault
## [1] FALSE
Мы можем обойти эту rowSums() причуду, применив ручной (и несколько произвольный) допуск в сравнении с плавающей запятой. Для этого нам нужно взять абсолютную разницу, а затем выполнить сравнение меньше чем с допуском:
abs(rowSums(CJ(0.01,0.29,0.70))-1)<1e-10;
## [1] TRUE
Следовательно:
psidom2 <- function(ss,m) { raw <- do.call(data.table::CJ,rep(list(ss),m)); raw[abs(rowSums(raw)-1)<1e-10,]; };
y <- smp(unname(as.matrix(psidom2(ss,m))));
nrow(y);
## [1] 4851
identical(x,y);
## [1] TRUE
Комбинации
Спасибо Джозефу Вуду за указание на то, что это действительно перестановки. Первоначально я назвал свою функцию combsum(), но переименовал ее в permsum(), чтобы отразить это откровение. И, как предположил Джозеф, можно модифицировать алгоритм для создания комбинаций, что можно сделать следующим образом, теперь соответствуя имени combsum():
combsum <- function(s,m,l=s) if (m==1L) matrix(s) else do.call(rbind,lapply(seq((s+m-1L)%/%m,min(l,s-m+1L)),function(x) unname(cbind(x,combsum(s-x,m-1L,x)))));
res <- combsum(100L,4L);
head(res);
## [,1] [,2] [,3] [,4]
## [1,] 25 25 25 25
## [2,] 26 25 25 24
## [3,] 26 26 24 24
## [4,] 26 26 25 23
## [5,] 26 26 26 22
## [6,] 27 25 24 24
tail(res);
## [,1] [,2] [,3] [,4]
## [7148,] 94 3 2 1
## [7149,] 94 4 1 1
## [7150,] 95 2 2 1
## [7151,] 95 3 1 1
## [7152,] 96 2 1 1
## [7153,] 97 1 1 1
Для этого потребовалось 3 изменения.
Во-первых, я добавил новый параметр l, что означает «лимит». По сути, чтобы гарантировать, что каждая рекурсия генерирует уникальную комбинацию, я заставляю каждый выбор быть меньше или равным любому предыдущему выбору в текущей комбинации. Это потребовало принятия текущего верхнего предела в качестве параметра l. В вызове верхнего уровня l можно просто установить по умолчанию на s, что в любом случае слишком велико для случаев, когда m>1, но это не проблема, поскольку это всего лишь один из двух верхних пределов, которые будут применяться во время генерации последовательности.
Вторым изменением, конечно же, было передать последний выбор x в качестве аргумента l при выполнении рекурсивного вызова в лямбда-выражении lapply().
Последнее изменение самое сложное. Последовательность выбора теперь должна быть вычислена следующим образом:
seq((s+m-1L)%/%m,min(l,s-m+1L))
Нижний предел должен был быть повышен с 1, используемого в permsum(), до самого низкого возможного выбора, который все еще допускал бы комбинацию убывающей величины. Наименьший возможный выбор, конечно, зависит от того, сколько столбцов еще предстоит создать; чем больше столбцов, тем больше «места» мы должны оставить для будущих выборов. Формула состоит в том, чтобы целочисленное деление s на m, но мы также должны эффективно «округлить», поэтому я добавляю m-1L перед делением. Я также рассматривал возможность деления с плавающей запятой, а затем вызова as.integer(ceiling(...)), но я думаю, что полностью целочисленный подход намного лучше.
Например, рассмотрим случай s=10,m=3. Чтобы получить сумму 10 с оставшимися 3 столбцами, мы не можем сделать выборку меньше 4, потому что тогда у нас не будет достаточного количества, чтобы произвести 10 без восхождения по комбинации. В этом случае формула делит 12 на 3, чтобы получить 4.
Верхний предел можно вычислить по той же формуле, что и в permsum(), за исключением того, что мы должны также применить текущий предел l, используя вызов min().
Я убедился, что моя новая функция combsum() ведет себя идентично функции IntegerPartitionsOfLength() Джозефа для многих случайных тестовых случаев со следующим кодом:
## helper function to sort a matrix within each row and then by columns
smc <- function(m) smp(t(apply(m,1L,sort)));
## test loop
for (i in seq_len(1000L)) {
repeat {
s <- sample(1:100,1L);
m <- sample(2:5,1L);
if (s>=m) break;
};
x <- combsum(s,m);
y <- IntegerPartitionsOfLength(s,m);
cat(paste0(s,',',m,'\n'));
if (!identical(smc(x),smc(y))) stop('bad.');
};
Бенчмаркинг
Общий автономный тестовый код:
library(microbenchmark);
library(data.table);
library(partitions);
library(gtools);
permsum <- function(s,m) if (m==1L) matrix(s) else do.call(rbind,lapply(seq_len(s-m+1L),function(x) unname(cbind(x,permsum(s-x,m-1L)))));
combsum <- function(s,m,l=s) if (m==1L) matrix(s) else do.call(rbind,lapply(seq((s+m-1L)%/%m,min(l,s-m+1L)),function(x) unname(cbind(x,combsum(s-x,m-1L,x)))));
IntegerPartitionsOfLength <- function(n, Lim, combsOnly = TRUE) { a <- 0L:n; k <- 2L; a[2L] <- n; MyParts <- vector("list", length=P(n)); count <- 0L; while (!(k==1L) && k <= Lim + 1L) { x <- a[k-1L]+1L; y <- a[k]-1L; k <- k-1L; while (x<=y && k <= Lim) {a[k] <- x; y <- y-x; k <- k+1L}; a[k] <- x+y; if (k==Lim) { count <- count+1L; MyParts[[count]] <- a[1L:k]; }; }; MyParts <- MyParts[1:count]; if (combsOnly) {do.call(rbind, MyParts)} else {MyParts}; };
GetDecimalReps <- function(s,m) { myPerms <- permutations(m,m); lim <- nrow(myPerms); intParts <- IntegerPartitionsOfLength(s,m,FALSE); do.call(rbind, lapply(intParts, function(x) { unique(t(sapply(1L:lim, function(y) x[myPerms[y, ]]))); })); };
smp <- function(m) m[do.call(order,as.data.frame(m)),];
smc <- function(m) smp(t(apply(m,1L,sort)));
bgoldst.perm <- function(s,m) permsum(s,m)/s;
psidom2 <- function(ss,m) { raw <- do.call(data.table::CJ,rep(list(ss),m)); raw[abs(rowSums(raw)-1)<1e-10,]; };
joseph.perm <- function(s,m) GetDecimalReps(s,m)/s;
bgoldst.comb <- function(s,m) combsum(s,m)/s;
joseph.comb <- function(s,m) IntegerPartitionsOfLength(s,m)/s;
Перестановки
## small scale
s <- 10L; m <- 3L; ss <- seq_len(s-1L)/s;
ex <- smp(bgoldst.perm(s,m));
identical(ex,smp(unname(as.matrix(psidom2(ss,m)))));
## [1] TRUE
identical(ex,smp(joseph.perm(s,m)));
## [1] TRUE
microbenchmark(bgoldst.perm(s,m),psidom2(ss,m),joseph.perm(s,m));
## Unit: microseconds
## expr min lq mean median uq max neval
## bgoldst.perm(s, m) 347.254 389.5920 469.1011 420.383 478.7575 1869.697 100
## psidom2(ss, m) 702.206 830.5015 1007.5111 907.265 1038.3405 2618.089 100
## joseph.perm(s, m) 1225.225 1392.8640 1722.0070 1506.833 1860.0745 4411.234 100
## large scale
s <- 100L; m <- 4L; ss <- seq_len(s-1L)/s;
ex <- smp(bgoldst.perm(s,m));
identical(ex,smp(unname(as.matrix(psidom2(ss,m)))));
## [1] TRUE
identical(ex,smp(joseph.perm(s,m)));
## [1] TRUE
microbenchmark(bgoldst.perm(s,m),psidom2(ss,m),joseph.perm(s,m),times=5L);
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst.perm(s, m) 1.286856 1.304177 1.426376 1.374411 1.399850 1.766585 5
## psidom2(ss, m) 6.673545 7.046951 7.416161 7.115375 7.629177 8.615757 5
## joseph.perm(s, m) 5.299452 10.499891 13.769363 12.680607 15.107748 25.259117 5
## very large scale
s <- 100L; m <- 5L; ss <- seq_len(s-1L)/s;
ex <- smp(bgoldst.perm(s,m));
identical(ex,smp(unname(as.matrix(psidom2(ss,m)))));
## Error: cannot allocate vector of size 70.9 Gb
identical(ex,smp(joseph.perm(s,m)));
## [1] TRUE
microbenchmark(bgoldst.perm(s,m),joseph.perm(s,m),times=1L);
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst.perm(s, m) 28.58359 28.58359 28.58359 28.58359 28.58359 28.58359 1
## joseph.perm(s, m) 50.51965 50.51965 50.51965 50.51965 50.51965 50.51965 1
Комбинации
## small-scale
s <- 10L; m <- 3L;
ex <- smc(bgoldst.comb(s,m));
identical(ex,smc(joseph.comb(s,m)));
## [1] TRUE
microbenchmark(bgoldst.comb(s,m),joseph.comb(s,m));
## Unit: microseconds
## expr min lq mean median uq max neval
## bgoldst.comb(s, m) 161.225 179.6145 205.0898 187.3120 199.5005 1310.328 100
## joseph.comb(s, m) 172.344 191.8025 204.5681 197.7895 205.2735 437.489 100
## large-scale
s <- 100L; m <- 4L;
ex <- smc(bgoldst.comb(s,m));
identical(ex,smc(joseph.comb(s,m)));
## [1] TRUE
microbenchmark(bgoldst.comb(s,m),joseph.comb(s,m),times=5L);
## Unit: milliseconds
## expr min lq mean median uq max neval
## bgoldst.comb(s, m) 409.0708 485.9739 556.4792 591.4774 627.419 668.4548 5
## joseph.comb(s, m) 2164.2134 3315.0138 3317.9725 3540.6240 3713.732 3856.2793 5
## very large scale
s <- 100L; m <- 6L;
ex <- smc(bgoldst.comb(s,m));
identical(ex,smc(joseph.comb(s,m)));
## [1] TRUE
microbenchmark(bgoldst.comb(s,m),joseph.comb(s,m),times=1L);
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst.comb(s, m) 2.498588 2.498588 2.498588 2.498588 2.498588 2.498588 1
## joseph.comb(s, m) 12.344261 12.344261 12.344261 12.344261 12.344261 12.344261 1
person
bgoldst
schedule
07.06.2016
combinationsилиpermutations. Судя по всему, кажется, что вы ищете «перестановки». См. это. - person Joseph Wood schedule 07.06.2016