Path: news.mathworks.com!newsfeed-00.mathworks.com!oleane.net!oleane!nerim.net!proxad.net!feeder1-2.proxad.net!74.125.46.134.MISMATCH!postnews.google.com!y33g2000prg.googlegroups.com!not-for-mail
From: Bob Alvarez <ralvarez@spambob.net>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Help in vectorizing code
Date: Mon, 27 Apr 2009 21:28:34 -0700 (PDT)
Organization: http://groups.google.com
Lines: 96
Message-ID: <e950fb95-5fa0-420c-a1ce-1ce64a360ebe@y33g2000prg.googlegroups.com>
References: <1e6f5074-edce-417a-b0a5-b62db0fb033f@w31g2000prd.googlegroups.com> 
	<128ed907-11c6-4e7a-a67a-736ececef060@s1g2000prd.googlegroups.com> 
	<gt5sqm$7mc$1@fred.mathworks.com>
NNTP-Posting-Host: 216.103.212.87
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: quoted-printable
X-Trace: posting.google.com 1240892914 16307 127.0.0.1 (28 Apr 2009 04:28:34 GMT)
X-Complaints-To: groups-abuse@google.com
NNTP-Posting-Date: Tue, 28 Apr 2009 04:28:34 +0000 (UTC)
Complaints-To: groups-abuse@google.com
Injection-Info: y33g2000prg.googlegroups.com; posting-host=216.103.212.87; 
	posting-account=pNrldAoAAAB5GlVBDJqvQa-_ca5yFHCO
User-Agent: G2/1.0
X-HTTP-UserAgent: Mozilla/5.0 (Windows; U; Windows NT 5.1; en-US; rv:1.9.0.9) 
	Gecko/2009040821 Firefox/3.0.9 (.NET CLR 3.5.30729),gzip(gfe),gzip(gfe)
Xref: news.mathworks.com comp.soft-sys.matlab:535868


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;