Thread Subject:

Help in vectorizing code

For a Monte Carlo simulator of compound Poisson noise, I need to add a
variable number of columns of each row of a matrix. So far, I have not
been able to vectorize the for loop and any help will be appreciated.

Here is an example:
>> d = magic(5)

d =

    17 24 1 8 15
    23 5 7 14 16
     4 6 13 20 22
    10 12 19 21 3
    11 18 25 2 9

>> idxs = round(linspace(1,3,5))

idxs =

     1 2 2 3 3

 for k=1:5
    totals(k) = sum(d(k,1:idxs(k)),2);
 end

>> totals'

ans =

    17 28 10 41 54

TIA

On Apr 27, 10:39=A0am, Bob <ralva...@spambob.net> wrote:
> For a Monte Carlo simulator of compound Poisson noise, I need to add a
> variable number of columns of each row of a matrix. So far, I have not
> been able to vectorize the for loop and any help will be appreciated.
>
> Here is an example:
>
> >> d =3D magic(5)
>
> d =3D
>
> =A0 =A0 17 =A0 =A024 =A0 =A0 1 =A0 =A0 8 =A0 =A015
> =A0 =A0 23 =A0 =A0 5 =A0 =A0 7 =A0 =A014 =A0 =A016
> =A0 =A0 =A04 =A0 =A0 6 =A0 =A013 =A0 =A020 =A0 =A022
> =A0 =A0 10 =A0 =A012 =A0 =A019 =A0 =A021 =A0 =A0 3
> =A0 =A0 11 =A0 =A018 =A0 =A025 =A0 =A0 2 =A0 =A0 9
>
> >> idxs =3D round(linspace(1,3,5))
>
> idxs =3D
>
> =A0 =A0 =A01 =A0 =A0 2 =A0 =A0 2 =A0 =A0 3 =A0 =A0 3
>
> =A0for k=3D1:5
> =A0 =A0 totals(k) =3D sum(d(k,1:idxs(k)),2);
> =A0end
>
> >> totals'
>
> ans =3D
>
> =A0 =A0 17 =A0 =A028 =A0 =A010 =A0 =A041 =A0 =A054
>
> TIA

I'm not real sure how to do it using the idxs variable you create. If
you have a way of making it more of a logical matrix, then it's easy:

idxs =3D [1 0 0 0 0;1 1 0 0 0;1 1 0 0 0;1 1 1 0 0;1 1 1 0 0];
sum(d.*idxs,2)

On Apr 27, 12:42=A0pm, jrenfree <jrenf...@gmail.com> wrote:
> I'm not real sure how to do it using the idxs variable you create. =A0If
> you have a way of making it more of a logical matrix, then it's easy:
>
> idxs =3D [1 0 0 0 0;1 1 0 0 0;1 1 0 0 0;1 1 1 0 0;1 1 1 0 0];
> sum(d.*idxs,2)

The idxs are actually samples of a poisson random variable. I did try
to make a mask array using poly2mask from the image processing
toolbox, but that is slower than doing the loop :-(

This is one way to make the mask quickly
mask = zeros(size(d));
mask(:,idxs) = 1;
mask = cumsum(mask,1);

Hth
Darren

Sorry, I think the line above
mask(:,idxs) = 1;
will not work as desired. Should probably be instead

mask([1:5],idxs) = 1;

Haven't the chance to check though.

% The data:
d = magic(5)
idxs = round(linspace(1,3,5));

% The engine.
F = zeros(size(d));
F(:,1) = 1;
F(sub2ind(size(d),1:5,idxs+1)) = -1;
sum(cumsum(F,2).*d,2)

Note: The For loop might be faster still!

Bob Alvarez <ralvarez@spambob.net> wrote in message <128ed907-11c6-4e7a-a67a-736ececef060@s1g2000prd.googlegroups.com>...
> The idxs are actually samples of a poisson random variable. I did try
> to make a mask array using poly2mask from the image processing
> toolbox, but that is slower than doing the loop :-(

  I don't think vectorizing this procedure in this way is a very realistic aim in trying to simulate a truly compound poisson process. There would have to be provision for an unlimited possible number of independent random variables to be summed which would mean a more or less unlimited number of columns in the desired matrix if you were using a fixed d matrix as in your example with a magic square.

  Suppose your independent random variables being summed are the results of the 'randn' output. I would think you would want to proceed something like this:

 R = poissrnd(lambda,m); % <-Or whatever your simulated poisson source
 t = zeros(1,m);
 for k = 1:m
  t(k) = sum(randn(1,R(k)));
 end

This requires only as many samples from 'randn' as are called for by the sum of the counts in R.

  I can conceive of a way of vectorizing this latter method using differences among cumsums of the 'randn' outputs and indexed by the counts in R, but I foresee accuracy difficulties with that approach for very long sequences. Without cumsumming I can see no effective way of eliminating the for-loops.

Roger Stafford

On Apr 27, 8:23=A0pm, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> =A0 I don't think vectorizing this procedure in this way is a very realis=
tic aim in trying to simulate a truly compound poisson process. =A0There wo=
uld have to be provision for an unlimited possible number of independent ra=
ndom variables to be summed which would mean a more or less unlimited numbe=
r of columns in the desired matrix if you were using a fixed d matrix as in=
 your example with a magic square.
>
For every set of trials there will be a maximum number of summed
terms, which determines the size of the matrix. I agree that in
principle the maximum can be arbitrarily large although the
probability of that is infinitesimally small. Using a matrix is
certainly more wasteful but is feasible for small values (say of the
order or hundreds) of the Poisson parameter. Probably larger values
are practically indistinguishable from normal.

Thanks to all for the ideas.

I wrote a mex function called sum2idx.c, which I list below. It is not
'bullet proof' but it gets the job done if used carefully.

Bob

sum2idx.c
-------------------------
/*
rowsums =3D sum2idx(M,idx)
sums rows of matrix M from 1:idx
where idx is a vector of length =3D # of rows of M

   REA 4/26/09
*/

#define V4_COMPAT
#include <matrix.h> /* Matlab matrices */
#include <mex.h>

#include <stddef.h> /* NULL */

#define notDblMtx(it) (!mxIsNumeric(it) || !mxIsDouble(it) ||
mxIsSparse(it) || mxIsComplex(it))

void mexFunction(int nlhs, /* Num return vals on lhs */
mxArray *plhs[], /* Matrices on lhs */
int nrhs, /* Num args on rhs */
const mxArray *prhs[] /* Matrices on rhs */
)
  {
  register double temp, sumval;
  register double *mtx,*idx,*outArray;
  register int i, k,idxval,nrows,ncols,nrows_idx,ncols_idx;
  mxArray *arg;

  if (nrhs !=3D 2) mexErrMsgTxt("sum2idx: requires 2 arguments.");

  /* ARG 1: MATRIX */
  arg =3D prhs[0];
  if notDblMtx(arg) mexErrMsgTxt("sum2idx: M arg must be a real non-
sparse matrix.");
  mtx =3D mxGetPr(arg);
  nrows =3D (int) mxGetM(arg);
  ncols =3D (int) mxGetN(arg);
  if ((nrows*ncols) =3D=3D 0) mexErrMsgTxt("sum2idx: M must be a non-empty
matrix.");

  /* ARG 2: idx vector */
  arg =3D prhs[1];
  if notDblMtx(arg) mexErrMsgTxt("sum2idx: idx arg must be a real non-
sparse 1D vector.");
  idx =3D mxGetPr(arg);
  nrows_idx =3D (int) mxGetM(arg);
  ncols_idx =3D (int) mxGetN(arg);
  /* assume the calling program forces idx to be a column vector*/
  if(ncols_idx!=3D1) mexErrMsgTxt("sum2idx: idx must be 1 D column
vector");
  if(nrows_idx!=3Dnrows) mexErrMsgTxt("sum2idx: idx must be same length
as number of rows of M");

  /* do the sum */
  plhs[0] =3D (mxArray *) mxCreateDoubleMatrix((int)nrows,1,mxREAL);
  if (plhs[0] =3D=3D NULL) mexErrMsgTxt("sum2idx: Error allocating result
matrix");
  //Get a pointer to the data space in our newly allocated memory
  outArray =3D mxGetPr(plhs[0]);
  for (i=3D0; i<nrows; i++){
      outArray[i] =3D 0.0;
      idxval =3D (int)idx[i];
      if(idxval>ncols) mexErrMsgTxt("sum2idx: idx must be <=3D number of
columns of M");
      if(idxval>0){ //do the sum
         for(k=3D0;k<idxval;k++) outArray[i] +=3D mtx[(k*ncols)+i];
      }
  }

  return;

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gt5sqm$7mc$1@fred.mathworks.com>...
> Suppose your independent random variables being summed are the results of the 'randn' output. I would think you would want to proceed something like this:
>
> R = poissrnd(lambda,m); % <-Or whatever your simulated poisson source
> t = zeros(1,m);
> for k = 1:m
> t(k) = sum(randn(1,R(k)));
> end
>
> This requires only as many samples from 'randn' as are called for by the sum of the counts in R.
>
> I can conceive of a way of vectorizing this latter method using differences among cumsums of the 'randn' outputs and indexed by the counts in R, but I foresee accuracy difficulties with that approach for very long sequences. Without cumsumming I can see no effective way of eliminating the for-loops.
---------
  In case it is of interest to you, here is a vectorization of the previous for-loop which I mentioned I could "conceive" of. Again, to generate m random compound poisson values using 'randn' as the independent random variables to be summed, do this:

 R = poissrnd(lambda,[1,m]); % <- Or your simulated poisson source
 p = cumsum([1,R]);
 r = randn(1,p(m+1)-1)+1; % <- Or some other ind. random variables
 s = [0,cumsum(r)];
 t = s(p(2:m+1))-s(p(1:m));

The row vector t will contain the desired compound poisson random values. You can compare the timing on this with the simple for-loop version.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gt6dqg$i52$1@fred.mathworks.com>...
> R = poissrnd(lambda,[1,m]); % <- Or your simulated poisson source
> p = cumsum([1,R]);
> r = randn(1,p(m+1)-1)+1; % <- Or some other ind. random variables
> s = [0,cumsum(r)];
> t = s(p(2:m+1))-s(p(1:m));

  Please remove that +1 from the third line. To be compatible with my other statements it should not be there. It had to do with something I was temporarily testing. The third line should read:

 r = randn(1,p(m+1)-1); % <- Or some other ind. random variables

Some updates:

1. There was a bug in the C code I posted. A new version is at the
end of this post.
2. I timed Roger's vectorized code, the for loop and the Mex code.
Results:
Elapsed time is 0.037361 seconds. % vectorized
Elapsed time is 0.019795 seconds. % for loop
Elapsed time is 0.027539 seconds. % Mex code
Doh! Apparently the JIT compiler is so good it beats the other 2
methods. The test matlab code is:

---
%% code from Roger Stafford of news group
lambda = 1000;
m=1000;
 R = poissrnd(lambda,[1,m]); % <- Or your simulated poisson source
 tic
 p = cumsum([1,R]);
 r = randn(1,p(m+1)-1); % <- Or some other ind. random variables
 s = [0,cumsum(r)];
 t = s(p(2:m+1))-s(p(1:m));
toc
clear p r s

%% do with for loop
tic
t2 = zeros(1,m);
for k = 1:m
     t2(k) = sum(randn(1,R(k)));
end
toc

%% do with mex code
tic
d = randn(m,max(R)+1);
t3 = sum2idx(d,R(:));
toc

----- sum2idx.c (updated) -----
/*
rowsums = sum2idx(M,idx)
sums rows of matrix M from 1:idx
where idx is a vector of length = # of rows of M
   REA 4/28/09
*/

#define V4_COMPAT
#include <matrix.h> /* Matlab matrices */
#include <mex.h>

#include <stddef.h> /* NULL */

#define notDblMtx(it) (!mxIsNumeric(it) || !mxIsDouble(it) ||
mxIsSparse(it) || mxIsComplex(it))

void mexFunction(int nlhs, /* Num return vals on lhs */
mxArray *plhs[], /* Matrices on lhs */
int nrhs, /* Num args on rhs */
const mxArray *prhs[] /* Matrices on rhs */
)
  {
  double temp, sumval;
  double *mtx,*idx,*outArray;
  int kcol, krow,idxval,nrows,ncols,nrows_idx,ncols_idx;
  mxArray *arg;

  if (nrhs != 2) mexErrMsgTxt("sum2idx: requires 2 arguments.");

  /* ARG 1: MATRIX */
  arg = prhs[0];
  if notDblMtx(arg) mexErrMsgTxt("sum2idx: M arg must be a real non-
sparse matrix.");
  mtx = mxGetPr(arg);
  nrows = (int) mxGetM(arg);
  ncols = (int) mxGetN(arg);
  if ((nrows*ncols) == 0) mexErrMsgTxt("sum2idx: M must be a non-empty
matrix.");

  /* ARG 2: idx vector */
  arg = prhs[1];
  if notDblMtx(arg) mexErrMsgTxt("sum2idx: idx arg must be a real non-
sparse 1D vector.");
  idx = mxGetPr(arg);
  nrows_idx = (int) mxGetM(arg);
  ncols_idx = (int) mxGetN(arg);
  /* assume the calling program forces idx to be a column vector*/
  if(ncols_idx!=1) mexErrMsgTxt("sum2idx: idx must be 1 D column
vector");
  if(nrows_idx!=nrows) mexErrMsgTxt("sum2idx: idx must be same length
as number of rows of M");

  /* do the sum */
  plhs[0] = (mxArray *) mxCreateDoubleMatrix((int)nrows,1,mxREAL);
  if (plhs[0] == NULL) mexErrMsgTxt("sum2idx: Error allocating result
matrix");
  //Get a pointer to the data space in our newly allocated memory
  outArray = mxGetPr(plhs[0]);
  for(krow=0;krow<nrows;krow++){
   sumval=0.0;
        idxval = (int)idx[krow];
        if(idxval>ncols) mexErrMsgTxt("sum2idx: idx must be <= number
of columns of M");
        if(idxval>0){ //do the sum
            for(kcol=0;kcol<idxval;kcol++){
                sumval += mtx[(kcol*nrows)+krow];
            }
        }
        outArray[krow] = sumval;
  }


  return;
}

Bob <ralvarez@spambob.net> wrote in message <1e6f5074-edce-417a-b0a5-b62db0fb033f@w31g2000prd.googlegroups.com>...
> For a Monte Carlo simulator of compound Poisson noise, I need to add a
> variable number of columns of each row of a matrix. So far, I have not
> been able to vectorize the for loop and any help will be appreciated.
>
> Here is an example:
> >> d = magic(5)
>
> d =
>
> 17 24 1 8 15
> 23 5 7 14 16
> 4 6 13 20 22
> 10 12 19 21 3
> 11 18 25 2 9
>
> >> idxs = round(linspace(1,3,5))
>
> idxs =
>
> 1 2 2 3 3
>
> for k=1:5
> totals(k) = sum(d(k,1:idxs(k)),2);
> end
>
> >> totals'
>
> ans =
>
> 17 28 10 41 54
>
> TIA

What about simple indexing after applying cumsum:

% data
d = magic(5)
idxs = round(linspace(1,3,5))

% engine
d2 = cumsum(d,2) ;
totals = d2(sub2ind(size(d2),1:numel(idxs),idxs))

hth
Jos

On Apr 28, 11:13 pm, "Jos " <#10...@fileexchange.com> wrote:
> Bob <ralva...@spambob.net> wrote in message <1e6f5074-edce-417a-b0a5-b62d=
b0fb0...@w31g2000prd.googlegroups.com>...
> > For a Monte Carlo simulator of compound Poisson noise, I need to add a
> > variable number of columns of each row of a matrix. So far, I have not
> > been able to vectorize the for loop and any help will be appreciated.
>
> > Here is an example:
> > >> d = magic(5)
>
> > d =
>
> >     17    24     1     8    15
> >     23     5     7    14    16
> >      4     6    13    20    22
> >     10    12    19    21     3
> >     11    18    25     2     9
>
> > >> idxs = round(linspace(1,3,5))
>
> > idxs =
>
> >      1     2     2     3     3
>
> >  for k=1:5
> >     totals(k) = sum(d(k,1:idxs(k)),2);
> >  end
>
> > >> totals'
>
> > ans =
>
> >     17    28    10    41    54
>
> > TIA
>
> What about simple indexing after applying cumsum:
>
> % data
> d = magic(5)
> idxs = round(linspace(1,3,5))
>
> % engine
> d2 = cumsum(d,2) ;
> totals = d2(sub2ind(size(d2),1:numel(idxs),idxs))
>
> hth
> Jos

Good idea. Timed the following code:
%% code from Jos
tic
d = cumsum(randn(m,ceil(max(R))),2);
s = d(sub2ind(size(d),1:numel(R),R));
toc
 result:
Elapsed time is 0.027088 seconds. % Jos
Elapsed time is 0.019429 seconds. % for loop
The for loop seems to win because it does not have to compute a full
matrix, just the specified length random variable at each trial.

On Apr 28, 8:59 am, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> "Roger Stafford" <ellieandrogerxy...@mindspring.com.invalid> wrote in mes=
sage <gt6dqg$i5...@fred.mathworks.com>...
> >  R = poissrnd(lambda,[1,m]); % <- Or your simulated poisson source
> >  p = cumsum([1,R]);
> >  r = randn(1,p(m+1)-1)+1; % <- Or some other ind. random variables
> >  s = [0,cumsum(r)];
> >  t = s(p(2:m+1))-s(p(1:m));
>
>   Please remove that +1 from the third line.  To be compatible with m=
y other statements it should not be there.  It had to do with something I=
 was temporarily testing.  The third line should read:
>
>  r = randn(1,p(m+1)-1); % <- Or some other ind. random variables

Very interesting. Thanks for your effort.

I tested the code with the following and got some interesting results:
-------
%% code from Roger Stafford of news group
tic
lambda = 1000;
m=1000;
 R = poissrnd(lambda,[1,m]); % <- Or your simulated poisson source
 p = cumsum([1,R]);
 r = randn(1,p(m+1)-1); % <- Or some other ind. random variables
 s = [0,cumsum(r)];
 t = s(p(2:m+1))-s(p(1:m));
toc

%% do with for loop
tic
R2 = poissrnd(lambda,[1,m]); % generate these for timing but do not
use them so sums are the same
t2 = zeros(1,m);
for k = 1:m
     t2(k) = sum(randn(1,R(k)));
end
toc
------------------
RESULTS:
Elapsed time is 0.042971 seconds.
Elapsed time is 0.023447 seconds.
The for loop is faster!!
Apparently the JIT compiler is able to handle the for loop code
effectively.

Anyway, thanks to all who answered.