Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Segmentation violation due to MEX code
Date: Tue, 7 Jul 2009 13:15:06 +0000 (UTC)
Organization: Imperial College London
Lines: 120
Message-ID: <h2vhoq$oot$1@fred.mathworks.com>
References: <h2tkfd$l84$1@fred.mathworks.com> <5b6d0b68-60e9-41be-a22f-ea755e2db20d@r33g2000yqn.googlegroups.com> <52a78667-52f6-4603-bfca-4cb04888b75a@i6g2000yqj.googlegroups.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: webapp-03-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1246972506 25373 172.30.248.38 (7 Jul 2009 13:15:06 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 7 Jul 2009 13:15:06 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1570235
Xref: news.mathworks.com comp.soft-sys.matlab:553397


Hi again!


Sorry I didn't answer before. I have to admit I have just realized all the experiments you have been posting. Thanks a lot!

The truth is that I posted the topic first through http://mathforum.org/kb/forum.jspa?forumID=80, and then I found the http://www.mathworks.de/matlabcentral/newsreader link and posted again, without knowing that I already had a copy of the first post in the latter forum.


Anyway, the codes to compare are a little different to what you posted, since I changed something that drastically increased the C performance. Here it goes:

#include "mex.h"

/*
 * Calculate_I_Values.c
 * Used to speed up the I mapping calculations 
 *
 * S = Calculate_I_Values(exposures.', (x(:, chan)).', double(y.'));
 *
 * This is a MEX-file for MATLAB.
*/


/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
    /* variable declarations here */
    double *e, *x, *y, *I;
    int eRows, eCols, xRows, xCols, yRows, yCols, index, ySize;
    int m, i, j;
    double num[256], den[256];

    e = mxGetPr(prhs[0]);
    eRows = mxGetN(prhs[0]);
    eCols = mxGetM(prhs[0]);

    x = mxGetPr(prhs[1]);
    xRows = mxGetN(prhs[1]);
    xCols = mxGetM(prhs[1]);

    y = mxGetPr(prhs[2]);
    yRows = mxGetN(prhs[2]);
    yCols = mxGetM(prhs[2]);
    ySize = mxGetNumberOfElements(prhs[2]);


    /* Allocate output memory */
    plhs[0] = mxCreateDoubleMatrix(1, 256, mxREAL);
    I = mxGetPr(plhs[0]);


    /* code here */
    for (m = 0; m <= 255; m++)
        num[m] = den[m] = 0;
    for (i = 0; i < yRows; i++)
        for (j = 0; j < yCols; j++)
        {
            index = i * yCols + j;
            if (index <= ySize)
                m = (int) y[index];
            num[m] += e[j] * x[i];
            den[m]++;
        }
    for (m = 0; m <= 255; m++)    
        if (den[m] == 0) I[m] = num[m];
        else I[m] = num[m] / den[m];
}

Note that now I directly use m as the index, rather than looping for all its range and performing if comparisons.

The code to beat was:
function S = calculate_I(exposures, x, y)

exposuresMat = repmat(exposures, [size(y,1) 1]);
xMat = repmat(x, [1 length(exposures)]);   
v = xMat .* exposuresMat;

% Slower
% S = accumarray(y(:)+1,v(:), [], @mean);

yp1 = y(:)+1;
S = accumarray(yp1,v(:))./accumarray(yp1,1);
S(isnan(S)) = 0;


I have checked the comparison running the functions only once, and measuring with tic and toc. Probably a more precise result can be obtained running them 10 times and averaging.

My results:
tic; S = Calculate_I_Values(exposures.', (x(:, chan)).', double(y.')); toc
Elapsed time is 0.599098 seconds.

tic; S = calculate_I(exposures, x(:, chan), y); toc;
Elapsed time is 6.367185 seconds.
>> tic; S = calculate_I(exposures, x(:, chan), y); toc;
Elapsed time is 1.204878 seconds.
>> tic; S = calculate_I(exposures, x(:, chan), y); toc;
Elapsed time is 1.039880 seconds.

The Matlab function needs to warm up! Jaja.


The arrays:
>> size(exposures)
ans =
     1    16
>> size(x(:, chan))
ans =
      361862           1
>> size(y)
ans =
      361862          16


The returned S values are the same for both functions, but with different length. Yo can check that visually with [S1; [S2; 0].'].', being S1 the answer for the MEX code and S2 the accumarray solution.

I still have to check why I get 1 value less with the second approach, but I am doing a complete different thing today!


Thanks a lot again
Jose