Thread Subject:

Getting rid of for loop

Hi,
can someone please help me get rid of the for loop in this code?

I tried but got stuck on:
allZ = rand(2,N);
allZ(allZ<=0.025)=fernM1*z;


CODE -------------------

%
% Graphics demo
%
figure(1) % Establish a graphic window - use the same on each time
clf % Clear it and set the axis
axis([-5 5 0 10]);
hold on % Points accumulate in the window as an animation

% Initial values
z = [0;1];

%
% 4 Matrices are key to the IFS Fractal description
%
fernM1 = [0 0 ; 0 0.16];
fernM2 = [0.2 -0.26 ; 0.23 0.22];
fernM3 = [-0.15 0.28 ; 0.26 0.24];
fernM4 = [0.85 0.04 ; -0.04 0.85];

% How many points are we making this time? Allocate the array
N = 50000;
allZ = zeros(2, N);
for kk = 1:N
    dd = rand();
    if dd <= 0.025
        z = fernM1 * z;
    elseif (dd > 0.025) && (dd <= 0.125)
        z = fernM2 * z + [0 ; 1.6];
    elseif (dd > 0.125) && (dd <= 0.225)
        z = fernM3 * z + [0 ; 0.44];
    else
        z = fernM4 * z + [0 ; 1.6];
    end
    allZ(:, kk) = z;
end

% plot the rest
plot(allZ(1,:), allZ(2,:), '.g','MarkerSize', 4)
display('The fern is done')

END CODE -------------------

"Peder Schmedling" <peder.schmedling@sensonor.no> wrote in message <i7q8nk$lg7$1@fred.mathworks.com>...
> Hi,
> can someone please help me get rid of the for loop in this code?
>
> I tried but got stuck on:
> allZ = rand(2,N);
> allZ(allZ<=0.025)=fernM1*z;
>
>
> CODE -------------------
>
> %
> % Graphics demo
> %
> figure(1) % Establish a graphic window - use the same on each time
> clf % Clear it and set the axis
> axis([-5 5 0 10]);
> hold on % Points accumulate in the window as an animation
>
> % Initial values
> z = [0;1];
>
> %
> % 4 Matrices are key to the IFS Fractal description
> %
> fernM1 = [0 0 ; 0 0.16];
> fernM2 = [0.2 -0.26 ; 0.23 0.22];
> fernM3 = [-0.15 0.28 ; 0.26 0.24];
> fernM4 = [0.85 0.04 ; -0.04 0.85];
>
> % How many points are we making this time? Allocate the array
> N = 50000;
> allZ = zeros(2, N);
> for kk = 1:N
> dd = rand();
> if dd <= 0.025
> z = fernM1 * z;
> elseif (dd > 0.025) && (dd <= 0.125)
> z = fernM2 * z + [0 ; 1.6];
> elseif (dd > 0.125) && (dd <= 0.225)
> z = fernM3 * z + [0 ; 0.44];
> else
> z = fernM4 * z + [0 ; 1.6];
> end
> allZ(:, kk) = z;
> end
>
> % plot the rest
> plot(allZ(1,:), allZ(2,:), '.g','MarkerSize', 4)
> display('The fern is done')
>
> END CODE -------------------

A hint:

N = 50000;
allZ = zeros(2, N);
dd = rand(N,1);

idx = dd <= 0.025;
allZ(:,idx) = repmat(fernM1 * z,1,sum(idx));

idx = (dd > 0.025 & dd <= 0.125);
allZ(:,idx) = repmat(fernM2 * z + [0 ; 1.6],1,sum(idx));

And so on

Oleg

"Peder Schmedling" <peder.schmedling@sensonor.no> wrote in message <i7q8nk$lg7$1@fred.mathworks.com>...
> Hi,
> can someone please help me get rid of the for loop in this code?
>
> I tried but got stuck on:
> allZ = rand(2,N);
> allZ(allZ<=0.025)=fernM1*z;
>
>
> CODE -------------------
>
> %
> % Graphics demo
> %
> figure(1) % Establish a graphic window - use the same on each time
> clf % Clear it and set the axis
> axis([-5 5 0 10]);
> hold on % Points accumulate in the window as an animation
>
> % Initial values
> z = [0;1];
>
> %
> % 4 Matrices are key to the IFS Fractal description
> %
> fernM1 = [0 0 ; 0 0.16];
> fernM2 = [0.2 -0.26 ; 0.23 0.22];
> fernM3 = [-0.15 0.28 ; 0.26 0.24];
> fernM4 = [0.85 0.04 ; -0.04 0.85];
>
> % How many points are we making this time? Allocate the array
> N = 50000;
> allZ = zeros(2, N);
> for kk = 1:N
> dd = rand();
> if dd <= 0.025
> z = fernM1 * z;
> elseif (dd > 0.025) && (dd <= 0.125)
> z = fernM2 * z + [0 ; 1.6];
> elseif (dd > 0.125) && (dd <= 0.225)
> z = fernM3 * z + [0 ; 0.44];
> else
> z = fernM4 * z + [0 ; 1.6];
> end
> allZ(:, kk) = z;
> end
>
> % plot the rest
> plot(allZ(1,:), allZ(2,:), '.g','MarkerSize', 4)
> display('The fern is done')
>
> END CODE -------------------

Like this?

z = [0;1];
N = 50000;
allZ = zeros(2,N);
dd = rand(1,N);
allZ(:,dd < 0.025) = fernM1*z;
allZ(:,(dd > 0.025) & (dd <= 0.125)) = fernM2*z+[0;1.6];
allZ(:,(dd > 0.125) & (dd <= 0.225)) = fernM3*z+[0;0.44];
allZ(:,(dd > 0.225)) = fernM4*z+[0;1.6];
%Not tested in ML

"Oleg Komarov" <oleg.komarovRemove.this@hotmail.it> wrote in message
> A hint:
>
> N = 50000;
> allZ = zeros(2, N);
> dd = rand(N,1);
>
> idx = dd <= 0.025;
> allZ(:,idx) = repmat(fernM1 * z,1,sum(idx));
>
> idx = (dd > 0.025 & dd <= 0.125);
> allZ(:,idx) = repmat(fernM2 * z + [0 ; 1.6],1,sum(idx));
>
> And so on
>
> Oleg

I don't think this will work because z itself is changing in the loop.

"Matt Fig" <spamanon@yahoo.com> wrote in message <i7qi9r$l4q$1@fred.mathworks.com>...
> "Oleg Komarov" <oleg.komarovRemove.this@hotmail.it> wrote in message
> > A hint:
> >
> > N = 50000;
> > allZ = zeros(2, N);
> > dd = rand(N,1);
> >
> > idx = dd <= 0.025;
> > allZ(:,idx) = repmat(fernM1 * z,1,sum(idx));
> >
> > idx = (dd > 0.025 & dd <= 0.125);
> > allZ(:,idx) = repmat(fernM2 * z + [0 ; 1.6],1,sum(idx));
> >
> > And so on
> >
> > Oleg
>
> I don't think this will work because z itself is changing in the loop.

Thats what I concluded with also (using my limited knowledge) but I was hoping that there were some gurus out there who could prove me wrong..

Oleg: I tried your code and it's kind of up the alley I was going. But it fails, probably due to what Matt is pointing out..

> > I don't think this will work because z itself is changing in the loop.
>
> Thats what I concluded with also (using my limited knowledge) but I was hoping that there were some gurus out there who could prove me wrong..
>
> Oleg: I tried your code and it's kind of up the alley I was going. But it fails, probably due to what Matt is pointing out..

That is true and my example will not work for the same reason.

Dear Peder,

> can someone please help me get rid of the for loop in this code?

Do you want to get rid of the loop of accelerate the function?

For the later case:
1. Move the RAND call outside the loop:
  d = rand(1, N);
  ...
  dd = d(kk);

2. For the most likely interval, the most IF checks are processed. Re-order them:
  if dd > 0.225 % Catch 77% of cases at first
    z = fernM4 * z + [0 ; 1.6];
  elseif (dd > 0.125) % Not needed: && (dd <= 0.225)
     z = fernM3 * z + [0 ; 0.44];
  elseif (dd > 0.025) % Not needed: && (dd <= 0.125)
      z = fernM2 * z + [0 ; 1.6];
  else % if dd <= 0.025
     z = fernM1 * z;
  end

Think: In the real life, you'd never check if d<0.125, if you have already sure, that d is not >=0.125 !

But this saves just 10% computing time on my Matlab2009a for the computations - the drawing needs much more time at all.

Kind regards, Jan

I think Matlab does not have recursion built-in function beside FILTER, which cannot process varying coefficients like yours. The only enhancement I can think of if get rid of calling RAND within the loop:

% Initial values
z = [0;1];

%
% 4 Matrices are key to the IFS Fractal description
%
fernM1 = [0 0 ; 0 0.16];
fernM2 = [0.2 -0.26 ; 0.23 0.22];
fernM3 = [-0.15 0.28 ; 0.26 0.24];
fernM4 = [0.85 0.04 ; -0.04 0.85];

N = 50000;
allZ = zeros(2, N);
dd = rand(1,N);
[~, where] = histc(dd, [0 0/025 0.125 0.225 1]);
fern = cat(3,fernM1,fernM2,fernM3,fernM4);
a = cat(2, [0; 0], [0 ; 1.6], [0 ; 0.44], [0 ; 1.6]);
for kk = 1:N
    w = where(kk);
    z = fern(:,:,w)*z + a(:,w);
    allZ(:,kk) = z;
end

% Bruno

Sorry for that typo at the line:

[~, where] = histc(dd, [0 0.025 0.125 0.225 1]);

% Bruno

Dear Bruno,

> for kk = 1:N
> w = where(kk);
> z = fern(:,:,w)*z + a(:,w);
> allZ(:,kk) = z;
> end

I've tried exactly this, but was slower than with the IF method on Matlab 2009a and a single core CPU. Perhaps it is a benefit to omit the "+ a(:, w)" if it is zero.
Is it faster on your machine?

BTW. I'm impressed by the resulting picture anytime I see it again.

Kind regards, Jan

On Sep 28, 1:23 am, "Jan Simon" <matlab.THIS_Y...@nMINUSsimon.de>
wrote:
> Dear Bruno,
>
> > for kk = 1:N
> >     w = where(kk);
> >     z = fern(:,:,w)*z + a(:,w);
> >     allZ(:,kk) = z;
> > end
>
> I've tried exactly this, but was slower than with the IF method on Matlab 2009a and a single core CPU. Perhaps it is a benefit to omit the "+ a(:, w)" if it is zero.
> Is it faster on your machine?
>
> BTW. I'm impressed by the resulting picture anytime I see it again.
>
> Kind regards, Jan

Yes, me too :-)
This should be the Barnsley fern: http://en.wikipedia.org/wiki/Barnsley's_fern

Maybe one of you guys can come up with a "Matlab way" of computing the
fern with N number of data points when yuo see the actual model on the
wiki page..

Anyways thank you all for your input so far, I'm learning by the hour!

"Jan Simon" <matlab.THIS_YEAR@nMINUSsimon.de> wrote in message <i7r90q$7sn$1@fred.mathworks.com>...
> Dear Bruno,
>
> > for kk = 1:N
> > w = where(kk);
> > z = fern(:,:,w)*z + a(:,w);
> > allZ(:,kk) = z;
> > end
>
> I've tried exactly this, but was slower than with the IF method on Matlab 2009a and a single core CPU. Perhaps it is a benefit to omit the "+ a(:, w)" if it is zero.
> Is it faster on your machine?
>

Yes it's 30% faster on my machine. I think the extra time is to copy the slice of array fern at each loop, not because of "+ a(:, w)". Here is an even faster way on my machine because no copy is necessary:

% Initial values
z = [0;1];

%
% 4 Matrices are key to the IFS Fractal description
%
fernM1 = [0 0 ; 0 0.16];
fernM2 = [0.2 -0.26 ; 0.23 0.22];
fernM3 = [-0.15 0.28 ; 0.26 0.24];
fernM4 = [0.85 0.04 ; -0.04 0.85];

tic
N = 50000;
allZ = zeros(2, N);
dd = rand(1,N);
[~, where] = histc(dd, [0 0.025 0.125 0.225 1]);
fern = { fernM1,fernM2,fernM3,fernM4 };
a = cat(2, [0; 0], [0 ; 1.6], [0 ; 0.44], [0 ; 1.6]);
for kk = 1:N
    w = where(kk);
    z = fern{w}*z + a(:,w);
    allZ(:,kk) = z;
end
toc

% Bruno