From help-request at octave dot org Wed May 18 19:18:19 2005
Subject: Re: for i=1:bignum, x(a(i)) += y(i); end
From: Miroslaw Kwasniak <mirek at zind dot ikem dot pwr dot wroc dot pl>
To: help at octave dot org
Date: Thu, 19 May 2005 02:17:35 +0200

On Wed, May 18, 2005 at 02:39:14PM -0400, Etienne Grossmann wrote:
> 
> 
>   Hi All,
> 
> has anyone written an oct-file for this?
> 
>   for i=1:bignum, x(a(i)) += y(i); end

Not oct, but m-function with about 15x speedup (and 2x slowup for small
bignum:)

[xa,a_unique] = sum_by_index(y,a);
x(a_unique)+=xa;


Timing results (t0-loop, t1-sum_by_index):

    octaveX:29> bignum=100;n_classes=10; sum_x
    t0 = 0.012109 t1 = 0.021936 Delta_max = 0

    octaveX:30> bignum=100;n_classes=100; sum_x
    t0 = 0.012310 t1 = 0.022559 Delta_max = 0

    octaveX:31> bignum=10000;n_classes=100; sum_x
    t0 = 0.75447  t1 = 0.059115 Delta_max = 0

    octaveX:32> bignum=10000;n_classes=10; sum_x
    t0 = 0.75039  t1 = 0.058044 Delta_max = 0

    octaveX:33> bignum=10000;n_classes=1000; sum_x
    t0 = 0.74890  t1 = 0.064619 Delta_max = 0

    octaveX:34> bignum=100000;n_classes=100; sum_x
    t0 = 7.4814   t1 = 0.47814 Delta_max = 0

  
Mirek

%%%%%%%%%%%%%%%% 

function [y,a] = sum_by_index(x,a)

  [ a i ] = sort(a); 
  x = x(i);

% Find location of unique values of a
  i = [ -2*abs(a(1))-1, a(:)' ];
  i = find(diff(i));
  a = a(i);

  lx = length(x);
  li = length(i);

% row indices
  r = ones(1,lx);
  r(i(2:end)) = 1-diff(i);
  r = cumsum(r);

% column indices
  c = zeros(1,lx);
  c(i)= 1:li;
  c(i(2:end)) -= 1:li-1;
  c = cumsum(c);

% Place x columnwise & sum
  R = max(r);
  C = length(a);
  y = zeros(1,R*C);
  y(sub2ind([R C],r,c)) = x;
  y = sum(reshape(y,R,C));

endfunction

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Test script

%bignum=100000
%n_classes=500

x=rand(1,bignum);
a=round(n_classes*x)+1;
x=a+(1:length(a))/100;

y1=y=zeros(1,max(a));

tic
for i=1:length(a), y(a(i)) += x(i); end
t0=toc

tic
[ya,a_unique] = sum_by_index(x,a);
y1(a_unique)+=ya;
t1=toc

Delta_max=max(abs(y-y1))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------