Pointers in mex files (mxCalloc)
in [Matlab]
|
From: James Tursa on 5 Aug 2010 14:39 "Aino" <aino.tietavainen(a)removeThis.helsinki.fi> wrote in message <i3ej3d$lc7$1(a)fred.mathworks.com>... > > So below is my code as it is now. I left the seg_i, seg_j and Cim as they were previously (I tried to declare them as arrays of size M+1 or N-M as "double seg_i[M+1], seg_j[M+1], Cim[N-M];", but this just created new errors). It compiles, but doesn't still work as it should. I believe the problem lies in here: > > seg_i[ii]=cop[i+ii];//TROUBLE, > > and more particularly in the left hand side (so again pointer stuff..). I don't understand why seg_i is not getting the cop[i+ii]-value in ii. Is it so that I am now trying to put double (cop[i+ii]) in to a double (seg_i[ii])? If so, I cannot figure out how to get around it. I know this must seen trivial to you, but I am really stuck here.. Can you provide a small example of inputs and expected outputs? Also, can you provide the m-code that you are trying to functionally duplicate with your mex file? James Tursa
From: Aino on 5 Aug 2010 15:41 "James Tursa" <aclassyguy_with_a_k_not_a_c(a)hotmail.com> wrote in message <i3f0g9$o0f$1(a)fred.mathworks.com>... > "Aino" <aino.tietavainen(a)removeThis.helsinki.fi> wrote in message <i3ej3d$lc7$1(a)fred.mathworks.com>... > > > > So below is my code as it is now. I left the seg_i, seg_j and Cim as they were previously (I tried to declare them as arrays of size M+1 or N-M as "double seg_i[M+1], seg_j[M+1], Cim[N-M];", but this just created new errors). It compiles, but doesn't still work as it should. I believe the problem lies in here: > > > > seg_i[ii]=cop[i+ii];//TROUBLE, > > > > and more particularly in the left hand side (so again pointer stuff..). I don't understand why seg_i is not getting the cop[i+ii]-value in ii. Is it so that I am now trying to put double (cop[i+ii]) in to a double (seg_i[ii])? If so, I cannot figure out how to get around it. I know this must seen trivial to you, but I am really stuck here.. > > Can you provide a small example of inputs and expected outputs? Also, can you provide the m-code that you are trying to functionally duplicate with your mex file? > > > James Tursa Here it comes: t=1:0.1:50; y=sin(t); M=3; r=0.1*std(y); c=0.01; [FuzzSampEnt]=FuzzySampEntropyv2(y,M,r,c) Using this, I got: FuzzSampEnt = 0.1355 Here's the Matlab code (non-parfor version): %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [FuzzSampEnt]=FuzzySampEntropyv2(cop,M,r,c) N=length(cop); expo=(log(log(2^c))/log(r)); CIM=[]; for m=[M,M+1] Cim=NaN(1,N-m); for i=1:N-m count=0; seg_i=cop(i:i+m-1); for j=1:N-m seg_j=cop(j:j+m-1); d=max(abs(seg_i-seg_j)); %if d=0, self match if d>0 count=count+exp(-(d^expo/c)); end end Cim(i)=(count)/(N-m-1); end CIM=cat(2,CIM,sum(Cim)/(N-m)); end FuzzSampEnt=-log(CIM(2)/CIM(1)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Another set of inputs/outputs: M=2; r=0.5*std(y); c=0.1;%same y gives me: FuzzSampEnt = 0.1049 Thanks, Aino
From: James Tursa on 5 Aug 2010 16:48 "Aino" <aino.tietavainen(a)removeThis.helsinki.fi> wrote in message <i3f44v$go5$1(a)fred.mathworks.com>... I haven't looked at your C code in much detail, so there may be some further improvements that can be made. Mostly I was just trying to get your M code to match your C code. Here are the results. The main change was the operator precedence in one of your equations. I like to use the following advice: Except for +, -, *, and /, use parentheses in all expressions. Also, your seg_i and seg_j don't need to be allocated. And I avoid using fabs because the lcc compiler version of this just sucks for speed. James Tursa #include <math.h> #include "mex.h" /*Subroutine to calculate the Fuzzy Sample Entropy*/ void FSE(double *cop, int M, double r, double c, mwSize N, double *FuzzSampEnt) { mwIndex m, i, ii, j, jj, k, l; double expo, d, *seg_i, *seg_j, count, sumCim, *Cim, CIM[2], f; expo=(log(log(pow(2, c)))/log(r)); Cim = mxCalloc(N-M+1, sizeof(double)); for (m=M;m<M+2;m++) { for (i=0; i<N-m; i++) { count=0; seg_i = cop + i; // No need to allocate & copy, just point! for (j=0; j<N-m; j++) { seg_j = cop + j; // No need to allocate & copy, just point! d=0; /*Finding the maximum absolute value of seg_i-seg_j*/ for (k=0; k<m; k++) { f = seg_i[k]-seg_j[k]; if( f < 0 ) f = -f; if (d < f) { d = f; } } /*if d=0, self match*/ if (d>0) { // count += exp(-(pow(d, expo/c))); <-- changed this line count += exp(-(pow(d, expo)/c)); } } Cim[i] = count/(N-m-1); } sumCim = 0; for(l=0;l<N-m;l++){ sumCim += Cim[l]; } if (m==M){ CIM[0]=sumCim/(N-m); } else if (m==M+1){ CIM[1]=sumCim/(N-m); } } *FuzzSampEnt = -log(CIM[1]/ CIM[0]); mxFree(Cim); } /*the gateway function*/ void mexFunction(int nlhs, mxArray *plhs[ ], int nrhs, const mxArray *prhs[ ]) { int M; double *cop, *FuzzSampEnt; double c, r; mwSize N; if(nrhs!=4) { mexErrMsgTxt("Four inputs required (cop,M,r,c)."); } /*Assign pointers*/ cop=mxGetPr(prhs[0]); N=mxGetNumberOfElements(prhs[0]); M=(int)(mxGetScalar(prhs[1])); r=mxGetScalar(prhs[2]); c=mxGetScalar(prhs[3]); plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); FuzzSampEnt = mxGetPr(plhs[0]); /*Call the subroutine */ FSE(cop, M, r, c, N, FuzzSampEnt); return; } Here are the results I get: t=1:0.1:50; y=sin(t); M=3; r=0.1*std(y); c=0.01; FuzzySampEntropyv2(y,M,r,c) FuzzySampEntropyv2_c(y,M,r,c) disp(' '); M=2; r=0.5*std(y); c=0.1;%same y FuzzySampEntropyv2(y,M,r,c) FuzzySampEntropyv2_c(y,M,r,c) >> fuzz ans = 0.1355 ans = 0.1355 ans = 0.1049 ans = 0.1049
From: Aino on 5 Aug 2010 18:00 "James Tursa" <aclassyguy_with_a_k_not_a_c(a)hotmail.com> wrote in message <i3f824$13v$1(a)fred.mathworks.com>... > "Aino" <aino.tietavainen(a)removeThis.helsinki.fi> wrote in message <i3f44v$go5$1(a)fred.mathworks.com>... > > I haven't looked at your C code in much detail, so there may be some further improvements that can be made. Mostly I was just trying to get your M code to match your C code. Here are the results. The main change was the operator precedence in one of your equations. I like to use the following advice: Except for +, -, *, and /, use parentheses in all expressions. Also, your seg_i and seg_j don't need to be allocated. And I avoid using fabs because the lcc compiler version of this just sucks for speed. > > James Tursa > > #include <math.h> > #include "mex.h" > > /*Subroutine to calculate the Fuzzy Sample Entropy*/ > void FSE(double *cop, int M, double r, double c, mwSize N, double *FuzzSampEnt) { > > mwIndex m, i, ii, j, jj, k, l; > double expo, d, *seg_i, *seg_j, count, sumCim, *Cim, CIM[2], f; > > expo=(log(log(pow(2, c)))/log(r)); > > Cim = mxCalloc(N-M+1, sizeof(double)); > > for (m=M;m<M+2;m++) { > > for (i=0; i<N-m; i++) { > count=0; > seg_i = cop + i; // No need to allocate & copy, just point! > > for (j=0; j<N-m; j++) { > seg_j = cop + j; // No need to allocate & copy, just point! > d=0; > /*Finding the maximum absolute value of seg_i-seg_j*/ > for (k=0; k<m; k++) { > f = seg_i[k]-seg_j[k]; > if( f < 0 ) f = -f; > if (d < f) { > d = f; > } > } > > /*if d=0, self match*/ > if (d>0) { > // count += exp(-(pow(d, expo/c))); <-- changed this line > count += exp(-(pow(d, expo)/c)); > } > > } > Cim[i] = count/(N-m-1); > } > > sumCim = 0; > for(l=0;l<N-m;l++){ > sumCim += Cim[l]; > } > > if (m==M){ > CIM[0]=sumCim/(N-m); > } > else if (m==M+1){ > CIM[1]=sumCim/(N-m); > } > > } > > *FuzzSampEnt = -log(CIM[1]/ CIM[0]); > > mxFree(Cim); > > } > > > /*the gateway function*/ > void mexFunction(int nlhs, mxArray *plhs[ ], > int nrhs, const mxArray *prhs[ ]) { > > int M; > double *cop, *FuzzSampEnt; > double c, r; > mwSize N; > > if(nrhs!=4) { > mexErrMsgTxt("Four inputs required (cop,M,r,c)."); > } > > /*Assign pointers*/ > cop=mxGetPr(prhs[0]); > N=mxGetNumberOfElements(prhs[0]); > M=(int)(mxGetScalar(prhs[1])); > r=mxGetScalar(prhs[2]); > c=mxGetScalar(prhs[3]); > > plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); > FuzzSampEnt = mxGetPr(plhs[0]); > > /*Call the subroutine */ > FSE(cop, M, r, c, N, FuzzSampEnt); > > return; > } > > > Here are the results I get: > > t=1:0.1:50; > y=sin(t); > M=3; r=0.1*std(y); c=0.01; > FuzzySampEntropyv2(y,M,r,c) > FuzzySampEntropyv2_c(y,M,r,c) > disp(' '); > M=2; r=0.5*std(y); c=0.1;%same y > FuzzySampEntropyv2(y,M,r,c) > FuzzySampEntropyv2_c(y,M,r,c) > > > >> fuzz > ans = > 0.1355 > ans = > 0.1355 > > ans = > 0.1049 > ans = > 0.1049 Thank you so much, Mr. Tursa! It works! And so easy, seg_i=cop+i. (That syntax is quite different than in Matlab.) I thought that the code requires to know somehow how many elements will be stored in seg_i. And you got rid of the whole ii and jj loops! In my computer the plain version with cop length 5000 dp takes 200 s, the parfor-loop version takes 100 s and this one takes a little more than 13 s! In my university server the numbers are 180 s, 27 s and 6.7s, so this is certainly quite an improvement! Now I can actually start calculating the results. :) With both of the computers. I also bought a C programming book today, maybe after I get the code running I have some time to start reading about how to make these files in the future. The speed of the m-files has been my problem for so many times, it always comes back to haunt me. Thank you so much, you have tough me a lot! Best Regards, Aino :)
From: James Tursa on 5 Aug 2010 18:47
"Aino" <aino.tietavainen(a)removeThis.helsinki.fi> wrote in message <i3fc9m$bg$1(a)fred.mathworks.com>... > > Thank you so much, Mr. Tursa! > > It works! And so easy, seg_i=cop+i. (That syntax is quite different than in Matlab.) I thought that the code requires to know somehow how many elements will be stored in seg_i. And you got rid of the whole ii and jj loops! That's one of the main advantages you can exploit in mex files. Every time you do an array slice operation in MATLAB it forces a data copy of the entire slice. If it is in a loop then this can really drag performance. But using a mex file you can often access the array slice directly without making a copy first and just loop through the data directly (loops are *not* bad in C). James Tursa |