[OpenSCAD] Barnsley Fern recursive

[Eric Buijs]

It's probably a crazy idea but I've created a simple script for the Barnsley
Fern fractal. The script works but it brings my PC, an 2011 iMac, to it's
knees (OpenSCAD takes over 10Gb of memory with the script below and I want
to increase the number of objects even further). I was wondering if someone
has any suggestions to improve/optimize the script and make it less memory
hungry. Or is OpenSCAD not suitable for this kind of work. Thanks.

m1 = [[0,0],[0,0.16]];

c1 = [0,0];

m2 = [[0.85,0.04],[-0.04,0.85]];

c2 = [0, 0.16];

m3 = [[0.2,-0.26],[0.23,0.22]];

c3 = [0,1.6];

m4 = [[-0.15,0.28],[0.26,0.24]];

c4 = [0,0.44];

function f1(p) = m1 * p + c1;

function f2(p) = m2 * p + c2;

function f3(p) = m3 * p + c3;

function f4(p) = m4 * p + c4;

module plotCircle(p) {
translate(100 * p) circle(r=0.5,$fn=50);
}

module BarnsleyFern(p,index) {
r = rands(0,1,1)[0];
plotCircle(p);
if (r <= 0.01) {
BarnsleyFern(f1(p), index-1);
}
else if (r > 0.01 && r <=0.86) {
BarnsleyFern(f2(p),index-1);
}
else if (r > 0.86 && r <=0.94) {
BarnsleyFern(f3(p),index-1);
}
else if (r > 0.94 && r < 0.99) {
BarnsleyFern(f4(p),index-1);
}
}


p = [0,0];
BarnsleyFern(p,20000);







--
Sent from: http://forum.openscad.org/

[Colin Smart]

If you are just drawing, consider using a different tool, such as processing.org

It avoids the extra effort of maintaining an ‘object’.

If you want an object, then I would generally look to do the numerical
processing elsewhere and import the shapes for rendering.

Others may have better ideas.

Colin

---------------------------------
Colin Smart (07968 049660)

Post by Eric Buijs
It's probably a crazy idea but I've created a simple script for the Barnsley
Fern fractal. The script works but it brings my PC, an 2011 iMac, to it's
knees (OpenSCAD takes over 10Gb of memory with the script below and I want
to increase the number of objects even further). I was wondering if someone
has any suggestions to improve/optimize the script and make it less memory
hungry. Or is OpenSCAD not suitable for this kind of work. Thanks.
m1 = [[0,0],[0,0.16]];
c1 = [0,0];
m2 = [[0.85,0.04],[-0.04,0.85]];
c2 = [0, 0.16];
m3 = [[0.2,-0.26],[0.23,0.22]];
c3 = [0,1.6];
m4 = [[-0.15,0.28],[0.26,0.24]];
c4 = [0,0.44];
function f1(p) = m1 * p + c1;
function f2(p) = m2 * p + c2;
function f3(p) = m3 * p + c3;
function f4(p) = m4 * p + c4;
module plotCircle(p) {
translate(100 * p) circle(r=0.5,$fn=50);
}
module BarnsleyFern(p,index) {
r = rands(0,1,1)[0];
plotCircle(p);
if (r <= 0.01) {
BarnsleyFern(f1(p), index-1);
}
else if (r > 0.01 && r <=0.86) {
BarnsleyFern(f2(p),index-1);
}
else if (r > 0.86 && r <=0.94) {
BarnsleyFern(f3(p),index-1);
}
else if (r > 0.94 && r < 0.99) {
BarnsleyFern(f4(p),index-1);
}
}
p = [0,0];
BarnsleyFern(p,20000);
--
Sent from: http://forum.openscad.org/
_______________________________________________
OpenSCAD mailing list
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org

[Eric Buijs]

Colin, thanks for your reply. I use P5.js (which is similar to Processing but
uses Javascript instead of Java) from time to time but it's just that I like
to 'torture' OpenSCAD to see where it's limits are ;-)



--
Sent from: http://forum.openscad.org/

[Torsten Paul]

That exact script takes less than a second for me and memory does
not grow over 350MB when running multiple times (using the latest
nightly build on Linux).

ciao,
Torsten.

[Mark Schafer]

yet another good reason for a binary somewhere - even if its not a
release...
I'm on Windows.

Post by Torsten Paul
That exact script takes less than a second for me and memory does
not grow over 350MB when running multiple times (using the latest
nightly build on Linux).
ciao,
  Torsten.

[Torsten Paul]

yet another good reason for a binary somewhere - even if > its not a release... I'm on Windows.

What's wrong with the one from the official download page?

http://www.openscad.org/downloads.html#snapshots

ciao,
Torsten.

[Mark Schafer]

I'm old - I forgot :(

Post by Torsten Paul

yet another good reason for a binary somewhere - even if > its not a
release... I'm on Windows.

What's wrong with the one from the official download page?
http://www.openscad.org/downloads.html#snapshots
ciao,
  Torsten.
_______________________________________________
OpenSCAD mailing list
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org

[Eric Buijs]

Thanks Thorsten, I will try the latest development snapshot.



--
Sent from: http://forum.openscad.org/

[Eric Buijs]

You're right Torsten, I've tried the latest experimental build for OSX and it
makes a huge difference in memory. I've updated the code because I hadn't
been very thoughtful about the base case (see below for updated code). One
more thing. The recursive loop refuses to go beyond 1300 (or so). After that
I get the message: ERROR: Recursion detected calling module 'BarnsleyFern'
Any thoughts on that. Thanks guys.

m1 = [[0,0],[0,0.16]];

c1 = [0,0];

m2 = [[0.85,0.04],[-0.04,0.85]];

c2 = [0, 0.16];

m3 = [[0.2,-0.26],[0.23,0.22]];

c3 = [0,1.6];

m4 = [[-0.15,0.28],[0.26,0.24]];

c4 = [0,0.44];

function f1(p) = m1 * p + c1;

function f2(p) = m2 * p + c2;

function f3(p) = m3 * p + c3;

function f4(p) = m4 * p + c4;

module plotCircle(p) {
translate(100 * p) circle(r=0.5,$fn=50);
}

module BarnsleyFern(p,index) {
r = rands(0,1,1)[0];
plotCircle(p);
//echo(r,index);
if (r <= 0.01) {
BarnsleyFern(f1(p), index-1);
}
else if (r > 0.01 && r <=0.86) {
BarnsleyFern(f2(p),index-1);
}
else if (r > 0.86 && r <=0.94) {
BarnsleyFern(f3(p),index-1);
}
else if (index > 0) {
BarnsleyFern(f4(p),index-1);
}
}


p = [0,0];
BarnsleyFern(p,1400);



--
Sent from: http://forum.openscad.org/

[Torsten Paul]

One more thing. The recursive loop refuses to go beyond 1300
(or so). After that I get the message: ERROR: Recursion detected
calling module 'BarnsleyFern' Any thoughts on that. Thanks guys.

That means it's running out of stack space, while not perfect,
it's trying to catch this issue before it just crashes.

I'm not sure if MacOS allows for changing the stack size available
to applications, but if it does, OpenSCAD tries to pick up that
limit.

ulimit -s 65532
/Applications/OpenSCAD.app/Contents/MacOS/OpenSCAD

ciao,
Torsten.

[Eric Buijs]

Thanks again but unfortunately I'm not able to increase the recursive loop
with this (on OSX).



--
Sent from: http://forum.openscad.org/

[Ronaldo]

Eric,

The probabilities in your code does not agree with the definition found in
Wikipedia <https://en.wikipedia.org/wiki/Barnsley_fern> . Besides, a full
recursion stop rule is missing in the recursive module so regardless the
memory allocation size you will always get a stack overflow with probability
one. Anyway, as this fractal may be generated without recursion that would
be a preferable path to go. Finally, you would need a very large number of
circles to get a reasonable picture of the fractal and this may surpass the
limits of the OpenSCAD node tree.



--
Sent from: http://forum.openscad.org/

[doug moen]

Ronaldo said: " Anyway, as this fractal may be generated without recursion
that would
be a preferable path to go."

Yes, but how do you code it without recursion in OpenSCAD?
The Barnsley Fern is an Iterated Function System, where the result at
iteration i+1
depends on the result of iteration i.

The only language feature that supports this kind of iteration, without
using recursion,
is the C-style for loop, which is only available in list comprehensions,
not in modules.
For this, you need a development snapshot with the "C-style for" option
enabled
in Preferences.

So, I would try generating all of the coordinates into a list, using a
C-style for loop
inside a list comprehension. And I'd use a simpler polygon for each point,
to further
reduce resource consumption: like a triangle, square, or hexagon.

Post by Ronaldo
Eric,
The probabilities in your code does not agree with the definition found in
Wikipedia <https://en.wikipedia.org/wiki/Barnsley_fern> . Besides, a full
recursion stop rule is missing in the recursive module so regardless the
memory allocation size you will always get a stack overflow with probability
one. Anyway, as this fractal may be generated without recursion that would
be a preferable path to go. Finally, you would need a very large number of
circles to get a reasonable picture of the fractal and this may surpass the
limits of the OpenSCAD node tree.
--
Sent from: http://forum.openscad.org/
_______________________________________________
OpenSCAD mailing list
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org

[Ronaldo]

Post by doug moen
The only language feature that supports this kind of iteration, without
using recursion,
is the C-style for loop, which is only available in list comprehensions,
not in modules.
For this, you need a development snapshot with the "C-style for" option
enabled
in Preferences.

You are right, Doug, and the following code takes that path. The nested for
in the module BarnsleyFern was needed to circumvent the OpenSCAD for limit
of at most 9999 elements. Anyway, this code is valid to at most 9999x9999
iterations.

m1 = [[0,0],[0,0.16]];
c1 = [0,0];
m2 = [[0.85,0.04],[-0.04,0.85]];
c2 = [0, 1.6];
m3 = [[0.2,-0.26],[0.23,0.22]];
c3 = [0, 1.6];
m4 = [[-0.15,0.28],[0.26,0.24]];
c4 = [0, 0.44];

function f1(p) = m1 * p + c1;
function f2(p) = m2 * p + c2;
function f3(p) = m3 * p + c3;
function f4(p) = m4 * p + c4;

module dot(p) translate(100 * p) square();

module BarnsleyFern(index) {
q = BarnsleyFern([0,0],index);
echo(len=len(q));
for(i=[0:9999:len(q)-1]) {
r = len(q)-i;
for(j=[0:min(r,9999)-1])
dot(q[i+j]);
}
}

function BarnsleyFern(p,index) =
[ for(i=0, q=p; i<index; i=i+1,
r = rands(0,1,1)[0],
q =
r &lt;= 0.01 ?
f1(q) :
r > 0.01 && r <=0.86?
f2(q) :
r > 0.86 && r <=0.93?
f3(q) :
f4(q) )
q ];


BarnsleyFern(200000);

The code generated the following image in 25sec:

<Loading Image...>

That is a detail look:

<Loading Image...>




--
Sent from: http://forum.openscad.org/

[doug moen]

Nice work.

Post by Ronaldo

Post by doug moen
The only language feature that supports this kind of iteration, without
using recursion,
is the C-style for loop, which is only available in list comprehensions,
not in modules.
For this, you need a development snapshot with the "C-style for" option
enabled
in Preferences.

You are right, Doug, and the following code takes that path. The nested for
in the module BarnsleyFern was needed to circumvent the OpenSCAD for limit
of at most 9999 elements. Anyway, this code is valid to at most 9999x9999
iterations.
m1 = [[0,0],[0,0.16]];
c1 = [0,0];
m2 = [[0.85,0.04],[-0.04,0.85]];
c2 = [0, 1.6];
m3 = [[0.2,-0.26],[0.23,0.22]];
c3 = [0, 1.6];
m4 = [[-0.15,0.28],[0.26,0.24]];
c4 = [0, 0.44];
function f1(p) = m1 * p + c1;
function f2(p) = m2 * p + c2;
function f3(p) = m3 * p + c3;
function f4(p) = m4 * p + c4;
module dot(p) translate(100 * p) square();
module BarnsleyFern(index) {
q = BarnsleyFern([0,0],index);
echo(len=len(q));
for(i=[0:9999:len(q)-1]) {
r = len(q)-i;
for(j=[0:min(r,9999)-1])
dot(q[i+j]);
}
}
function BarnsleyFern(p,index) =
[ for(i=0, q=p; i<index; i=i+1,
r = rands(0,1,1)[0],
q =
r &lt;= 0.01 ?
r > 0.01 && r <=0.86?
r > 0.86 && r <=0.93?
f4(q) )
q ];
BarnsleyFern(200000);
<http://forum.openscad.org/file/t1275/BarnsleyFern_200000.png>
<http://forum.openscad.org/file/t1275/BarnsleyFern_200000_detail.png>
--
Sent from: http://forum.openscad.org/
_______________________________________________
OpenSCAD mailing list
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org

[Eric Buijs]

Ronaldo, first I take my hat of for that code. It works perfectly (once I
translated the &lt; to smaller than, <). I wasn't even aware of the C-style
for option in OpenSCAD and therefore discarded the possibility to use
iteration. The Barnsley Fern looks pretty neat taking only 1 minute and 18
seconds on my old computer. Thank you all guys.



--
Sent from: http://forum.openscad.org/

[Ronaldo]

The C-style for loop available for list comprehension is not the only way to
address the Barnsley Fern fractal drawing because its building process is
tail recursive. As OpenSCAD is able to eliminate recursion of some forms of
tail recursive functions (not modules), it is in principle possible to have
an iteration written in a recursive way. Find bellow an alternative
recursion to the function and module defined in my previous post.

module BarnsleyFern(index) {
q = BarnsleyFern(index);
for(i=[0:9999:len(q)-1]) {
r = len(q)-i;
for(j=[0:min(r,9999)-1])
dot(q[i+j]);
}
}

function nextP(p) =
let(r = rands(0,1,1)[0])
r <= 0.01 ?
f1(p) :
r > 0.01 && r <=0.86?
f2(p) :
r > 0.86 && r <=0.93?
f3(p) :
f4(p) ;

function BarnsleyFern(index, bag=[[0,0]]) =
index==0 ?
bag :
BarnsleyFern(index-1, concat( [nextP(bag[0])], bag ));

Although exploring the tail recursion elimination in an elegant solution,
this approach is not so efficient as the previous one. My previous code have
a linear behavior: the running time is approximately proportional to the
parameter index (O(n)). For my surprise, although the above code run as an
iteration, its running time is O(n2). A closer look solves the apparent
mystery: the main operation in the function BarnsleyFern is a concat of a
single element list with a growing list (bag); possibly concat operates by
rewriting the output list from the input lists so the total number of
operations of all concats is O(n2). This approach would benefit from a
better implementation of concat.




--
Sent from: http://forum.openscad.org/

[Eric Buijs]

Ronaldo, a nice comparison of iterative vs tail recursive. The latter took
almost half an hour to finish on my PC, indeed not so efficient. Thanks.



--
Sent from: http://forum.openscad.org/

[Torsten Paul]

Post by Eric Buijs
Ronaldo, a nice comparison of iterative vs tail recursive. The latter took
almost half an hour to finish on my PC, indeed not so efficient. Thanks.

That is not the fault of the tail recursion elimination. This part is
pretty much the same as a normal loop. As Ronaldo mentioned the slowness
comes from the part where it generates the data by copying everything
for each iteration.

ciao,
Torsten.

[Ronaldo Persiano]

Torsten,

I partially disagree. I don't know the intrinsic processes of OpenSCAD but,
from the running times I got, I deduct that the iterative step of
generating one element of a list in the C-style solution is constant time
in contrast with the linear time the *concat *is done in the tail recursion
elimination strategy due to the successive copies. However, I don't see any
way to write a tail recursion solution to generate a list without resorting
to *concat*. So, tail recursion will be generally less competitive than
iterations for list generations.

Post by Torsten Paul

Post by Eric Buijs
Ronaldo, a nice comparison of iterative vs tail recursive. The latter

took

Post by Eric Buijs
almost half an hour to finish on my PC, indeed not so efficient. Thanks.

That is not the fault of the tail recursion elimination. This part is
pretty much the same as a normal loop. As Ronaldo mentioned the slowness
comes from the part where it generates the data by copying everything
for each iteration.
ciao,
Torsten.
_______________________________________________
OpenSCAD mailing list
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org

[Torsten Paul]

However, I don't see any way to write a tail recursion solution
to generate a list without resorting to /concat/.

True, but that still makes it an issue with how concat works and
not with the tail-recursion elimination which *is* running as a
loop.
There is probably a way to prevent all that copying as by
definition there is no modification of existing vectors.

So, tail recursion will be generally less competitive than
iterations for list generations.

Yes, I guess that's a fair point looking from user perspective.

ciao,
Torsten.

[doug moen]

Ronaldo: "This approach would benefit from a
better implementation of concat."

Torsten: "There is probably a way to prevent all that copying as by
definition there is no modification of existing vectors."

Pretty much every functional language has an operation that appends an
element onto the head of a list in constant time. In Lisp, this operation
is called "cons". The result is internally represented as a linked list
that can be traversed from front to back in linear time, same as an array,
but indexing into the middle of the linked list requires linear time
instead of constant time.

So we could consider adding a 'cons(element, list)' operation that returns
a new list in constant time. Internally, there would be two list
representations: the existing representation as a C++ std::vector, and a
new representation (a linked list node or "cons" node) that contains the
first element and a list containing the remaining elements. But both
representations are just lists from the user's perspective: they print the
same and support the same set of operations.

If you google "functional data structures", you will see that cleverer and
more complicated solutions exist. An alternative is to replace std::vector
with a general purpose functional list data structure. And that might speed
up `concat` in Ronaldo's program. By contrast, with my solution, `concat`
behaves the same as before, and returns a flat list represented as a
std::vector. So existing code needs to be changed to use `cons` in order to
get a performance benefit.

But, maybe my solution is easier to implement? I was thinking that
Value::type() would return VECTOR for a linked list, and Value::toVector()
would internally modify a linked list to a vector, modifying the Value to
the new representation, before returning a vector. But some operations
would be modified to detect a linked list and operate efficiently on it
without converting it to a vector. The idea is to limit the amount of code
that needs to be changed.

Post by Torsten Paul

However, I don't see any way to write a tail recursion solution
to generate a list without resorting to /concat/.

True, but that still makes it an issue with how concat works and
not with the tail-recursion elimination which *is* running as a
loop.
There is probably a way to prevent all that copying as by
definition there is no modification of existing vectors.

So, tail recursion will be generally less competitive than
iterations for list generations.

Yes, I guess that's a fair point looking from user perspective.
ciao,
Torsten.
_______________________________________________
OpenSCAD mailing list
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org