Как преобразовать рекурсивную функцию в код mex?

У меня есть рекурсивная функция, выбираемая в коде MATLAB следующим образом:

   function nk=choose(n, k)
        if (k == 0)
            nk=1;
        else
            nk=(n * choose(n - 1, k - 1)) / k;
        end
    end

Код используется для вычисления комбинации между n и k. Я хочу ускорить его, используя код mex. Я попытался написать код mex как

double choose(double* n, double* k)
{
   if (k==0) 
        return 1;
   else
        return (n * choose(n - 1, k - 1)) / k;
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *n, *k, *nk;
    int mrows, ncols;
    plhs[0] = mxCreateDoubleMatrix(1,1, mxREAL);
    /* Assign pointers to each input and output. */
    n = mxGetPr(prhs[0]);    
    k = mxGetPr(prhs[1]);
    nk = mxGetPr(plhs[0]);
    /* Call the recursive function. */
    nk=choose(n,k);
}

Однако это не работает. Не могли бы вы помочь мне изменить код mex, который может реализовать приведенный выше код MATLAB? Спасибо


matlab mex
person Jame    schedule 24.08.2016    source источник
comment
Что ж, в таких языках программирования, как C, рекурсия обычно не является предпочтительным инструментом программирования для достижения больших скоростей — каждый уровень рекурсии требует, чтобы вы хранили стек, упорядочивали аргументы, передаваемые вашей функции, вызывали функцию, позволяли функции делать это. вещь и, в конце концов, всплывают значения, восстанавливая стек на каждом этапе... накладные расходы существенны для простых функций. Итак, учитывая, что ваша проблема не диктует, вы делаете это рекурсивно, не делайте этого рекурсивно.   -  person Marcus Müller    schedule 24.08.2016
comment
Я думаю, что моя реализация в коде mex была неправильной. Не могли бы вы помочь мне исправить это   -  person Jame    schedule 24.08.2016
comment
Нет: дело в том, что как бы правильно вы ни реализовали свой мекс, он упускает суть.   -  person Marcus Müller    schedule 24.08.2016
comment
Возможно, он пытается реализовать пример неэффективных методов программирования в качестве примера того, чего следует избегать при программировании.   -  person Sembei Norimaki    schedule 24.08.2016
comment
На самом деле, у него есть проблема XY. Это не работает — тоже не лучшее описание проблемы.   -  person dasdingonesin    schedule 24.08.2016
comment
вы могли бы сделать это, может быть, с matlabFunction(), но, как говорят остальные, не используйте здесь рекурсию.   -  person Ander Biguri    schedule 24.08.2016

Ответы (2)


arrow_upward
1
arrow_downward

Следующий код исправляет вашу реализацию C mex.
Проблема, конечно, не в рекурсии...
В вашем коде вместо значений используются указатели (в C важно использовать указатели только в нужных местах).

Вы можете использовать встроенную функцию Matlab: nchoosek
См.: http://www.mathworks.com/help/matlab/ref/nchoosek.html

Работает следующий код:

//choose.c

#include "mex.h"

double choose(double n, double k)
{
    if (k==0) 
    {
        return 1;
    }
    else
    {
        return (n * choose(n - 1, k - 1)) / k;
    }
}


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *n, *k, *nk;
    int mrows, ncols;
    plhs[0] = mxCreateDoubleMatrix(1,1, mxREAL);
    /* Assign pointers to each input and output. */
    n = mxGetPr(prhs[0]);    
    k = mxGetPr(prhs[1]);
    nk = mxGetPr(plhs[0]);

    /* Call the recursive function. */
    //nk=choose(n,k);
    *nk = choose(*n, *k);
}

Скомпилируйте его в Matlab: mex choose.c

Выполнить:
choose(10,5)
ans =

252

Это не неэффективная реализация...
Я помогаю исправить вашу реализацию, чтобы использовать ее как "неэффективный пример".

Измерьте выполнение реализации rahnema1:
tic;n = 1000000;k = 500000;nk = prod((k+1:n) .* prod((1:n-k).^ (-1/(n-k))) );toc
Прошедшее время составляет 0,022855 секунды.

Измерьте выполнение реализации Choose.mexw64:
tic;n = 1000000;k = 500000;nk = Choose(1000000, 500000);toc
Истекшее время составляет 0,007952 секунды.
(заняло немного меньше времени, чем prod((k+1:n).* prod((1:n-k).^ (-1/(n-k))))).

Измерение рекурсии Matlab, получение ошибки (даже для n=700 и k=500):
ic;n = 700;k = 500;nk = RecursiveFunctionTest(n, k);toc
Достигнут максимальный предел рекурсии 500 . Используйте set(0,'RecursionLimit',N), чтобы изменить ограничение. Имейте в виду, что превышение доступного пространства стека может привести к сбою MATLAB и/или вашего компьютера.

tic;n = 700;k = 400;nk = RecursiveFunctionTest(n, k);toc
Истекшее время составляет 0,005635 секунды. Очень неэффективно...

Измерение встроенной функции Matlab nchoosek:
tic;nchoosek(1000000, 500000);toc
Внимание: результат может быть неточным. Коэффициент больше 9,007199e+15 и точен только до 15 цифр. В nchoosek на 92 Прошедшее время составляет 0,005081 секунды.

Вывод:
Вам нужно реализовать mex-файл C без использования рекурсии и принять меры.

Измерить без рекурсии:

static double factorial(double number) 
{
    int x;
    double fac = 1;

    if (number == 0)
    {
        return 1.0;
    }    

    for (x = 2; x <= (int)number; x++)
    {
        fac = fac * x;
    }

    return fac;
}



double choose(double n, double k)
{
    if (k == 0) 
    {
        return 1.0;
    }
    else
    {
        //n!/((n–k)! k!) 
        return factorial(n)/(factorial(n-k)*factorial(k));
    }
}

tic;choose(1000000, 500000);toc

Прошедшее время составляет 0,003079 секунды.

Быстрее...

person Rotem    schedule 24.08.2016
comment
Отличный ответ. Я хочу использовать код mex, потому что он быстрее, чем файл MATLAB, хотя может быть немного меньше времени. Но когда мы используем его в цикле (for), это действительно полезно. Спасибо - person Jame; 25.08.2016

arrow_upward
1
arrow_downward

нет необходимости использовать биномиальные коэффициенты, которые можно реализовать в Matlab:

function nk=nchoosek2(n, k)
    if n-k > k
        nk = prod((k+1:n) .* prod((1:n-k).^ (-1/(n-k))));
    else
        nk = prod((n-k+1:n) .* prod((1:k).^ (-1/k)) ) ;
    end
end
person rahnema1    schedule 24.08.2016