From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: please help optimize this ('find' is too slow)
Date: Wed, 10 Dec 2008 20:35:04 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 43
Message-ID: <ghp95o$27o$>
References: <gh46f4$pg2$> <gh7ror$eno$> <gh7vs8$fgq$> <gh9faq$fb4$> <gh9qmi$l5c$> <gh9ske$6he$> <gha1vj$850$> <ghefbe$ell$> <gher83$sfd$> <ghfk6p$mu6$> <ghi6f9$3um$> <ghkqra$pq9$>
Reply-To: <HIDDEN>
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: 1228941304 2296 (10 Dec 2008 20:35:04 GMT)
NNTP-Posting-Date: Wed, 10 Dec 2008 20:35:04 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: comp.soft-sys.matlab:506191

"Jveer " <> wrote in message <ghkqra$pq9$>...
> .......
> Your new idea sounds very promising. I'm afraid i have no idea how to put that into code though. i reckon as long as the sort is outside the loop it'll be fine. it 0.09s on the sample data i used. profiler results .tiff can be found on the link above. i know it's a lot to ask but any chance you could put that into code?
> .......

  I worked out the details of that scheme I told you about, Jayveer, and include it below.  However, I do so with misgivings, since my impression is that it is very much slower than your method.  You can try it out and see.

  It does as I described with three preliminary sorts followed by a for-loop that handles fifty-three of the cubic cells at a time and therefore makes only some two percent as many trips through the loop.  However the processing during each trip is much more time-consuming.  Only at the point "q = find(r>0)" does it narrow down the length of the arrays it is dealing with to a small number.  So use it at your own risk.  I compared its results with those of yours on some fairly large random data sets and they do agree.

  Note that the routine fails unless the lower limit is less than or equal to the corresponding upper limit for each coordinate in each cell.  The row vectors X, Y, and Z are your P(1,:), P(2,:), and P(3,:) and the row vectors Lx, Ly, Lz, Ux, Uy, and Uz come from your NCx, NCy, and NCz in the obvious way.  The quantity called w refers to the number of bits available in one 'double' format number, and it must not exceed 53, though it may be made less if desired.

  If nothing else, this example can serve to show you things like 1) how to find the inverse of a permutation, 2) how 53 binary results can be packed into a single matlab 'double' number, or 3) how the 'cumsum' function can be used to assist in doing so.

n = length(X);
m = length(Lx);
w = 53;
[t,px] = sort([Lx,Ux,X]); px(px) = 1:2*m+n;
[t,py] = sort([Ly,Uy,Y]); py(py) = 1:2*m+n;
[t,pz] = sort([Lz,Uz,Z]); pz(pz) = 1:2*m+n;
PiC = zeros(1,n);
for k = 0:w:m-1
 u = k+1:min(k+w,m);
 v = 2.^(0:length(u)-1);
 u = [u u+m];
 v = [v -v];
 s = zeros(1,2*m+n);
 s(px(u)) = v;
 s = cumsum(s);
 rx = s(px(2*m+1:2*m+n));
 s = zeros(1,2*m+n);
 s(py(u)) = v;
 s = cumsum(s);
 ry = s(py(2*m+1:2*m+n));
 s = zeros(1,2*m+n);
 s(pz(u)) = v;
 s = cumsum(s);
 rz = s(pz(2*m+1:2*m+n));
 r = bitand(bitand(rx,ry),rz);
 q = find(r>0);
 [f,e] = log2(r(q));
 PiC(q) = e+k;