Thread Subject: Feeding array subsets into mexCallMATLAB

Hi,

I'm a bit of a newbie to the world of C/Mex. I want to convolve some data from a Mex file and consequently will use mexCallMATLAB to do so. Given that z = conv(x,y), x is a subset of data that I'm processing in a big loop and I want y to be a column of an input matrix fed to the Mex function. So I set up my inputs to the convolution:
mxArray *conv_in[2];

the y matrix is obtained from the input:
double *sd_in;
sd_in = mxGetPr(prhs[1]);

which is a 2-D matrix from which I want to use one column. I guess the assignment would be something like:
conv_in[1] = sd_in[<index denoting some subset>]; /* at least that's the MATLAB way of doing it */

and then call the convolution:
mexCallMATLAB(1,conv_out,2,conv_in,"conv");

But I have no idea how to obtain that subset. I guess I would need to write the appropriate data to another variable and assign that to conv_in[1], but am not sure of the most efficient way of doing this. Any thoughts would be greatly appreciated!

(Of course, the same principal will need to be applied to the x data)

Many thanks,

Chris

"Christopher Hummersone" <wheely_chairs@hotmail.com> wrote in message <h7oq7t$t64$1@fred.mathworks.com>...
>
> I want to convolve some data from a Mex file and consequently will use mexCallMATLAB to do so. Given that z = conv(x,y), x is a subset of data that I'm processing in a big loop and I want y to be a column of an input matrix fed to the Mex function. So I set up my inputs to the convolution:
> mxArray *conv_in[2];
>
> the y matrix is obtained from the input:
> double *sd_in;
> sd_in = mxGetPr(prhs[1]);
>
> which is a 2-D matrix from which I want to use one column. I guess the assignment would be something like:
> conv_in[1] = sd_in[<index denoting some subset>]; /* at least that's the MATLAB way of doing it */

No, what you have above is assigning a double array subset directly to an mxArray * variable. Basically, you have to do what you are apparently trying to avoid ... make a variable that has a copy of the data first and then pass that to conv. So you need to make the copy of the column either on the MATLAB side or the mex side, and then pass that to conv via your mexCallMATLAB call.

> and then call the convolution:
> mexCallMATLAB(1,conv_out,2,conv_in,"conv");
>
> But I have no idea how to obtain that subset. I guess I would need to write the appropriate data to another variable and assign that to conv_in[1], but am not sure of the most efficient way of doing this. Any thoughts would be greatly appreciated!

Do it on the MATLAB side since it will be much easier. Writing all of the code to create the column copy could be done in the mex routine, of course, but it won't save you any time so why bother doing it there? But maybe you have other reasons for doing this step in the mex routine itself. If so, just use mxCreateDoubleMatrix to create a new mxArray variable. Then use mxGetPr to get the address of the real double data, and copy your desired input column data into it. Then use that result in your mexCallMATLAB call. After you are done, use mxDestroyArray on your temporary column variable.

James Tursa

Hey James,

Thanks for your reply. I'm doing everything in mex simply because the MATLAB function may be called millions of times so wanted to maximise efficiency. It seems to be doing pretty much what I want although I am getting regular bus errors. In my attempts to debug I have been compiling with three switches: -v -g -argcheck. However, when I use the switches the code compiles and runs fine, if I don't, it compiles but I get bus errors. Any thoughts? I can post the code if necessary (I'm trying to avoid it as it's 230 lines). Btw I'm running MATLAB 7.4 on an intel Mac.

Thanks a lot,

Chris

On 9 Sep, 13:08, "Christopher Hummersone" <wheely_cha...@hotmail.com>
wrote:
> Hey James,
>
> Thanks for your reply. I'm doing everything in mex simply because the MATLAB function may be called millions of times so wanted to maximise efficiency.

Then do *everything* on the MEX side of things. The overhead
of passing data in and out of the matlab engine will easily
outweigh most of the benefits you want from going MEX.

> It seems to be doing pretty much what I want although I am getting regular bus errors. In my attempts to debug I have been compiling with three switches: -v -g -argcheck. However, when I use the switches the code compiles and runs fine, if I don't, it compiles but I get bus errors.  Any thoughts?

You have a rougue pointer somewere - a pointer that tries
to acces some area of memory it is not supposed to. In the
latetr case the program tries to access un-allocated memory.
You have the same error when you run in debug mode, but
coincidence has it the program accesses accessible memory in
that case. Of course, it's a part of memory the program is not
supposed to touch, so you are in for undefined behaviour.

Rune

Thanks a lot for that Rune. I'm only just starting with C/Mex so it's probably a stupid coding error, but of course my lack of experience makes tracking it down 10 times harder. I've included it below, maybe you'll be able to spot a silly mistake just at a glance!

//begin C code

#include "math.h"
#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    //function martin_xcorr(hc_L,hc_R,frameCount,numchans,frame_length,maxlag,window,inhib,noverlap,maxdeg,sd,warpers)

// declarations
double *hc_L,*hc_R,*frameCount_in,*frame_length_in,*maxlag_in,*window,*inhib,*noverlap_in,*ccg_out,\
*maxdeg_in,*sd,*warpers_in,*ccg_azi_out,*sccg_out,*tau,*ccg_warp_in,*warp_data,*warp_result,\
*conv_data_in,*conv_in_gauss,*conv_result,maxlag_temp,*ccg_raw,*ccg_r,*gaussian;


int numchans,frameCount,frame_length,maxlag,noverlap,numsamples,offset,sample,sample_late,\
lengthCCG,leftn,rightn,ccg_index,ccg_r_index,ccg_raw_index,\
firstFrame,maxdeg,newLengthCCG;

    mwSize ccg_r_dims[3],ccg_raw_dims[3],ccg_dims[3],sccg_dims[3];
    
    mwIndex i,j,n,m,o;
    
mxArray *warp_out[1],*warp_in[3],*conv_out[1],*conv_in[2],*ccg_raw_mx[1],*ccg_r_mx[1],*gaussian_mx[1];

// inputs
hc_L = mxGetPr(prhs[0]);
hc_R = mxGetPr(prhs[1]);
frameCount_in = mxGetPr(prhs[2]);
frame_length_in = mxGetPr(prhs[3]);
maxlag_in = mxGetPr(prhs[4]);
window = mxGetPr(prhs[5]);
inhib = mxGetPr(prhs[6]);
noverlap_in = mxGetPr(prhs[7]);
maxdeg_in = mxGetPr(prhs[8]);
sd = mxGetPr(prhs[9]);
warpers_in = mxGetPr(prhs[10]);
numsamples = mxGetM(prhs[0]);
numchans = mxGetN(prhs[0]);

// derive variables
frameCount = (int) *frameCount_in;
frame_length = (int) *frame_length_in;
maxlag = (int) *maxlag_in;
noverlap = (int) *noverlap_in;
maxdeg = (int) *maxdeg_in;
offset = (int) noverlap*frame_length;
lengthCCG = (int) (2*maxlag)+1;
newLengthCCG = (int) (2*maxdeg)+1;
firstFrame = (int) 1;

// matrix dimensions
ccg_r_dims[0] = frame_length;
ccg_r_dims[1] = lengthCCG;
ccg_r_dims[2] = numchans;
ccg_raw_dims[0] = (noverlap+1)*frame_length;
    ccg_raw_dims[1] = lengthCCG;
    ccg_raw_dims[2] = numchans;
ccg_dims[0] = numchans;
ccg_dims[1] = lengthCCG;
ccg_dims[2] = frameCount;
sccg_dims[0] = numchans;
sccg_dims[1] = newLengthCCG;
sccg_dims[2] = frameCount;

// temp variables
ccg_raw_mx[0] = mxCreateNumericArray((mwSize)3,ccg_raw_dims,mxDOUBLE_CLASS,mxREAL);
ccg_raw = mxGetPr(ccg_raw_mx[0]);
ccg_r_mx[0] = mxCreateNumericArray((mwSize)3,ccg_r_dims,mxDOUBLE_CLASS,mxREAL);
ccg_r = mxGetPr(ccg_r_mx[0]);
gaussian_mx[0] = mxCreateDoubleMatrix((mwSize)newLengthCCG,(mwSize)numchans,mxREAL);
gaussian = mxGetPr(gaussian_mx[0]);

// outputs
plhs[0] = mxCreateNumericArray((mwSize)3,ccg_dims,mxDOUBLE_CLASS,mxREAL);
plhs[1] = mxCreateNumericArray((mwSize)3,sccg_dims,mxDOUBLE_CLASS,mxREAL);
plhs[2] = mxCreateNumericArray((mwSize)3,sccg_dims,mxDOUBLE_CLASS,mxREAL);
    //mexPrintf("frameCount: %d \n",frameCount);
    
ccg_out = mxGetPr(plhs[0]);
ccg_azi_out = mxGetPr(plhs[1]);
sccg_out = mxGetPr(plhs[2]);

// initiate tau for warping
warp_in[0] = mxCreateDoubleMatrix((mwSize)lengthCCG,(mwSize)1,mxREAL);
tau = mxGetPr(warp_in[0]);
maxlag_temp = (double) maxlag;
for ( m = 0; m < lengthCCG; m++ ) {
tau[m] = (((double) m - maxlag_temp)/maxlag_temp);
}
// initiate ccg input to warping
warp_in[1] = mxCreateDoubleMatrix((mwSize)lengthCCG,1,mxREAL);
ccg_warp_in = mxGetPr(warp_in[1]);
// initiate warping array
warp_in[2] = mxCreateDoubleMatrix((mwSize)newLengthCCG,1,mxREAL);
warp_data = mxGetPr(warp_in[2]);
warp_result = mxGetPr(warp_out[0]);
    //warp_in[3] = mxCreateString("spline");

// initiate gaussians for skeleton ccg
    
conv_in[0] = mxCreateDoubleMatrix((mwSize)newLengthCCG,(mwSize)1,mxREAL);
conv_in[1] = mxCreateDoubleMatrix((mwSize)newLengthCCG,(mwSize)1,mxREAL);
conv_data_in = mxGetPr(conv_in[0]);
conv_in_gauss = mxGetPr(conv_in[1]);
conv_result = mxGetPr(conv_out[0]);
for ( i = 0; i < numchans; i++ ) {
for ( m = 0; m < newLengthCCG; m++ ) {
gaussian[i*newLengthCCG+m] = exp((-pow(m-(*maxdeg_in),2))/(2*(pow(sd[i],2))));
}
}

// Do the calculations
for ( j = 0; j < frameCount; j++ ) {
for ( i = 0; i < numchans; i++ ) {
if ( firstFrame==1 ) {

for ( n = 0; n < frame_length*(noverlap+1); n++ ) {
sample = i*numsamples+n;
for ( m = 0; m < lengthCCG; m++ ) { // raw cross-correlations
ccg_raw_index = n+((m+(i*lengthCCG))*((noverlap+1)*frame_length));
leftn = fmax(sample-maxlag+m,sample);
rightn = fmax(sample+maxlag-m,sample);
ccg_raw[ccg_raw_index] = hc_L[leftn]*hc_R[rightn];
}
}

for ( n = 0; n < frame_length; n++ ) {
sample = i*numsamples+n;
for ( m = 0; m < lengthCCG; m++ ) { // integrate
ccg_raw_index = n+((m+(i*lengthCCG))*((noverlap+1)*frame_length));
ccg_r_index = n+((m+(i*lengthCCG))*frame_length);
ccg_index = i+((m+(j*lengthCCG))*numchans);
ccg_r[ccg_r_index] = window[0]*ccg_raw[ccg_raw_index];
for ( o = 1; o < noverlap*frame_length; o++ ) {
ccg_r[ccg_r_index] += window[o]*ccg_raw[ccg_raw_index+o];
}
ccg_r[ccg_r_index] /= (frame_length*noverlap); // running ccg
ccg_r[ccg_r_index] *= inhib[sample];
ccg_out[ccg_index] += ccg_r[ccg_r_index]/frame_length; // write to o/p
}
}
}

else { // if not first frame (firstFrame = 0)
for ( n = 0; n < noverlap*frame_length; n++ ) { // move ccgs
for ( m = 0; m < lengthCCG; m++ ) {
ccg_raw_index = n+((m+(i*lengthCCG))*((noverlap+1)*frame_length));
ccg_raw[ccg_raw_index] = ccg_raw[ccg_raw_index+frame_length];
}
}
for ( n = 0; n < frame_length; n++ ) {
sample_late = (i*numsamples)+(j*frame_length)+n+offset;
sample = (i*numsamples)+(j*frame_length)+n;
for ( m = 0; m < lengthCCG; m++ ) {
// indices
ccg_r_index = n+((m+(i*lengthCCG))*frame_length);
ccg_raw_index = n+((m+(i*lengthCCG))*((noverlap+1)*frame_length));
ccg_index = i+((m+(j*lengthCCG))*numchans);
leftn = fmax(sample_late-maxlag+m,sample_late);
rightn = fmax(sample_late+maxlag-m,sample_late);

// calculations
ccg_raw[ccg_raw_index+offset] = hc_L[leftn]*hc_R[rightn];
ccg_r[ccg_r_index] = window[0]*ccg_raw[ccg_raw_index];
for ( o = 1; o < noverlap*frame_length; o++ ) { // integrate
ccg_r[ccg_r_index] += window[o]*ccg_raw[ccg_raw_index+o];
}
ccg_r[ccg_r_index] /= (frame_length*noverlap); // running ccg
ccg_r[ccg_r_index] *= inhib[sample];
ccg_out[ccg_index] += ccg_r[ccg_r_index]/frame_length; // write to o/p
}
}
}
} // end frequency loop
firstFrame = 0;
} // end frame loop
    
mxDestroyArray(*warp_out);
    mxDestroyArray(*warp_in);
    mxDestroyArray(*conv_out);
    mxDestroyArray(*conv_in);
    mxDestroyArray(*ccg_raw_mx);
    mxDestroyArray(*ccg_r_mx);
    mxDestroyArray(*gaussian_mx);
return;
}

//end C code

On 9 Sep, 14:47, "Christopher Hummersone" <wheely_cha...@hotmail.com>
wrote:
> Thanks a lot for that Rune. I'm only just starting with C/Mex so it's probably a stupid coding error, but of course my lack of experience makes tracking it down 10 times harder. I've included it below, maybe you'll be able to spot a silly mistake just at a glance!

Hah! No one (and I mean _no_one_) spot those kinds of rouge
pointers at a glance. It takes time and effort to learn
how to program C, for precisely the reasons that causes
the problems you have.

There is no way to proceed other than meticulously work your
way through the code, checking every single details yourself.
Yeah, I know, that takes insane amounts of time, but that's
how the world works. If you want to reap the benefits of
programming C, this is the rite of initiation you will have
to go through.

Rune

> Hah! No one (and I mean _no_one_) spot those kinds of rouge
> pointers at a glance. It takes time and effort to learn
> how to program C, for precisely the reasons that causes
> the problems you have.

Well that makes me feel slightly better, but scared at the same time!

> There is no way to proceed other than meticulously work your
> way through the code, checking every single details yourself.
> Yeah, I know, that takes insane amounts of time, but that's
> how the world works. If you want to reap the benefits of
> programming C, this is the rite of initiation you will have
> to go through.

Stupid question I know, but what details should I check for?

Thanks a lot for your time :-)

Chris

On 9 Sep, 17:35, "Christopher Hummersone" <wheely_cha...@hotmail.com>
wrote:
> > Hah! No one (and I mean _no_one_) spot those kinds of rouge
> > pointers at a glance. It takes time and effort to learn
> > how to program C, for precisely the reasons that causes
> > the problems you have.
>
> Well that makes me feel slightly better, but scared at the same time!

You should be scared. Once upon a time I spent a couple of
months searching for a bug like the one you are dealing with.
True, the code was a bit bigger, about 20 A4 pages in print,
but still.

> > There is no way to proceed other than meticulously work your
> > way through the code, checking every single details yourself.
> > Yeah, I know, that takes insane amounts of time, but that's
> > how the world works. If you want to reap the benefits of
> > programming C, this is the rite of initiation you will have
> > to go through.
>
> Stupid question I know, but what details should I check for?

Anything and everything that can go wrong. If it was
possible to make a short concise answer to that question,
everybody and their brother would program C in no time, flat.

Rune

> > Stupid question I know, but what details should I check for?
>
> Anything and everything that can go wrong. If it was
> possible to make a short concise answer to that question,
> everybody and their brother would program C in no time, flat.
>
> Rune

Well it didn't actually take that long to find the problem and, as I suspected, it was something stupid (trying to get a pointer to the output of mexCallMATLAB before calling the function - d'oh!)! Oh well, it works beautifully now, and has potentially reduced the total processing time from 12 weeks to 1! That makes all the pain very worthwhile!

Thanks a lot for all the advice :-)

Chris

> > Stupid question I know, but what details should I check for?
>
> Anything and everything that can go wrong. If it was
> possible to make a short concise answer to that question,
> everybody and their brother would program C in no time, flat.
>
> Rune

Well it didn't actually take that long to find the problem and, as I suspected, it was something stupid (trying to get a pointer to the output of mexCallMATLAB before calling the function - d'oh!)! Oh well, it works beautifully now, and has potentially reduced the total processing time from 12 weeks to 1! That makes all the pain very worthwhile!

Thanks a lot for all the advice :-)

Chris