|
FFRANK has asked for the wisdom of the Perl Monks concerning the following question:
Dear Monks,
I was suggested yesterday to try PDL, here's the built of the fractal structure (diagonal 4 out of 4x4 plan +=1):
#!/usr/bin/perl -w
use strict;
use PDL;
use Algorithm::Loops qw(NestedLoops);
my $element = pdl ([1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]);
my $tmp = pdl ($element);
my @depth = qw (2 3 4 5 6 7 8);
foreach my $depth (@depth) {
my $same = pdl ($tmp);
my $fractal = zeroes 4**$depth,4**$depth;
NestedLoops([([0..3])x2],sub{
my ($x,$y) = (@_);
my ($a,$b) = ($x*(4**($depth-1)),$y*(4**($depth-1)));
my ($c,$d) = ($a+((4**($depth-1))-1),$b+((4**($depth-1))-1));
my $square = slice $fractal, "$a:$c,$b:$d";
$square += $same;
},);
for my $coord (0..3) {
my $a = $coord * (4**($depth-1));
my $b = $a + ((4**($depth-1))-1);
my $square = slice $fractal, "$a:$b,$a:$b";
$square += 1;
}
$tmp = pdl ($fractal);
$outFractal = $tmp;
}
At level 6 and more: Out of memory! - but that was fast. So I see that building the structure isn't optimal memorywise, even with PDL.
Frank
|
Re: Fractal structure: PDL, memory, time
by BrowserUk (Patriarch) on Jul 19, 2007 at 22:13 UTC
|
The problem is, that although the counts at each position will be a maximum of depth, meaning that you can get away with using a single byte per position upto depth=255, the size of the matrix grows (hyper?) geometrically. That means that even using a single byte per position, with no further overhead, by the time you reach depth=8, you already require 4GB of storage.
By the time you reached level 16, the storage requirement is 16e18 bytes, which according to this table is roughly 3 times the 'All words ever spoken by human beings':
depth: X x Y required m
+emory:
1: 4 x 4 == 16 _b
2: 16 x 16 == 256 _b
3: 64 x 64 == 4 Kb
4: 256 x 256 == 64 Kb
5: 1024 x 1024 == 1 Mb
6: 4096 x 4096 == 16 Mb
7: 16384 x 16384 == 256 Mb
8: 65536 x 65536 == 4 Gb
9: 262144 x 262144 == 64 Gb
10: 1048576 x 1048576 == 1 Tb
11: 4194304 x 4194304 == 16 Tb
12: 16777216 x 16777216 == 256 Tb
13: 67108864 x 67108864 == 4 Pb
14: 268435456 x 268435456 == 64 Pb
15: 1073741824 x 1073741824 == 1 Eb
16: 4294967296 x 4294967296 == 16 Eb
17: 17179869184 x 17179869184 == 256 Eb
18: 68719476736 x 68719476736 == 4 Zb
19: 274877906944 x 274877906944 == 64 Zb
20: 1099511627776 x 1099511627776 == 1 Yb
21: 4398046511104 x 4398046511104 == 16 Yb
22: 17592186044416 x 17592186044416 == 256 Yb
23: 70368744177664 x 70368744177664 == 4 2^90b
24: 281474976710656 x 281474976710656 == 64 2^90b
25: 1125899906842624 x 1125899906842624 == 1 2^100b
26: 4503599627370496 x 4503599627370496 == 16 2^100b
27: 18014398509481984 x 18014398509481984 == 256 2^100b
28: 72057594037927936 x 72057594037927936 == 4 2^110b
29: 288230376151711740 x 288230376151711740 == 64 2^110b
30: 1152921504606847000 x 1152921504606847000 == 1 2^120b
31: 4611686018427387900 x 4611686018427387900 == 16 2^120b
32: 18446744073709552000 x 18446744073709552000 == 256 2^120b
But if you think about it, even if you could calculate the data at say depth 8, then you would need a screen 70 feet square to view it.
As with all fractals, as you 'zoom' in, you only want to calculate the data required to produce some small subset of the total picture: say 1280x 1024. Because you want the ordered positions in each row, you would have to calculate a 'slice' of the data, but still far less than doing the whole thing.
At depth 1, the value for any given X,Y position is simply: $p = X == Y ? $depth : $depth - 1;
At level 2, the value for any given X,Y position is slightly more complicated:
if( X == Y ) {
return depth;
}
elsif ( ( X,Y ) is in a 4x4 block on the major diagonal ) {
return $depth - 1;
}
else {
## Move the 4x4 block to the origin and treat as for depth -1
X %= 4;
Y %= 4;
depth--;
recurse;
}
This is a recursive data structure, and calculating the value for any given (X,Y,depth) combinaton can be done recursively, quickly and without requiring a lot of memory. The following routine is close (but no cigar!):
#! perl -slw
use strict;
$|++;
our $DEPTH ||= 3;
sub fract2D {
use integer;
my( $x, $y, $depth ) = @_;
if( $x < 16 and $y < ) { ## A block on the major diagonal
return $x == $y ? $depth : $depth -1;
}
elsif( $depth > 1 ) {
return fract2D( $x % 4**$depth, $y % 4**$depth, $depth -1 );
}
else {
return fract2D( $x % 4**$depth, $y % 4**$depth, $depth );
}
}
for my $depth ( $DEPTH ) {
for my $y ( 0 .. 4**$depth-1 ) {
my $c = 0;
for my $x ( 0 .. 4 ** $depth -1 ) {
printf +( ( ++$c % 4 ) == 0 ? " %2d " : " %2d" ),
fract2D( $x, $y, $depth );
}
print +( $y % 4 ) == 3 ? "\n" : '';
}
print '';
}
But I seem to have left my tuits in my other head, and its washday.
Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
"Science is about questioning the status quo. Questioning authority".
In the absence of evidence, opinion is indistinguishable from prejudice.
| [reply]
[d/l]
[select] |
Re: Fractal structure: PDL, memory, time
by BrowserUk (Patriarch) on Jul 22, 2007 at 21:38 UTC
|
Try this. It's correct at levels 1, 2, & 3, and Ithink it is correct for all levels, but it gets harder to discern.
I found this really tough to get right, whether through recursion or iteration. I'd given up and moved on before the pattern dawned on me. I think a mathematician may be able to reduce this to a formula, maybe involving some kind of polynomial, but it's beyond my abilities.
#! perl -slw
use strict;
$|++;
our $DEPTH ||= 3;
sub f{
use integer;
my( $x, $y, $depth ) = @_;
my $bsize = 4**( $depth-1 );
{
if( $x == $y ) {
return $depth;
}
elsif( $x/4 == $y/4 ) {
return $depth-1;
}
elsif( $x/$bsize == $y/$bsize ) {
$x %= 4;
$y %= 4;
--$depth;
redo;
}
else {
$x %= $bsize;
$y %= $bsize;
--$depth;
redo;
}
}
}
for my $depth ( $DEPTH ) {
for my $y ( 0 .. 4**$depth-1 ) {
my $c = 0;
for my $x ( 0 .. 4 ** $depth -1 ) {
printf +( ( ++$c % 4 ) == 0 ? " %2d " : " %2d" ),
f( $x, $y, $depth );
}
print +( $y % 4 ) == 3 ? "\n" : '';
}
print '';
}
__END__
Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
"Science is about questioning the status quo. Questioning authority".
In the absence of evidence, opinion is indistinguishable from prejudice.
| [reply]
[d/l] |
|
|
Hi BrowserUK,
I tried also to get your code working, and could'nt. I verified level 4 after your update and it does not yet match - although close ! (if you compare w/ output from PDL fractal code).
I'll meet a mathematicien and if I find something better, I'll post in this thread.
I'm also looking at simpler solutions using nDeep arrays.
Thanks very much,
Frank
| [reply] |
|
|
Try these. fc() is just an Inline C version of f().
I'm pretty sure that they are now correct, but you should verify them yourself.
#! perl -slw
use strict;
use Inline C => << '__C__', NAME => "_627288"; # => CLEAN_AFTER_BUILD
+=> 0;
#include <stdio.h>
int fc( double x, double y, int d ) {
double bsize = pow( 4, ( d-1 ) );
if( x == y )
return d;
else if( (int)( x / 4 ) == (int)( y / 4 ) )
return d-1;
else if( (int)( x / bsize ) == (int)( y / bsize ) )
return 1 + fc( fmod( x, bsize ), fmod( y, bsize ), d-1 );
else
return fc( fmod( x, bsize ), fmod( y, bsize ), d-1 );
}
__C__
sub f{
my( $x, $y, $d ) = @_;
my $bsize = 4**($d-1);
if( $x == $y ) {
return $d
}
elsif( int( $x / 4 ) == int( $y / 4 ) ) {
return $d-1;
}
elsif( int( $x / $bsize ) == int( $y / $bsize ) ) {
return 1 + f( $x % $bsize, $y % $bsize, $d-1 );
}
else{
return f( $x % $bsize, $y % $bsize, $d-1 );
}
}
Using this I have produced plots of 1000x1000 chunks upto level 31 successfully. At that point the coordinates get so large that they blow the limits of doubles. I'm also fairly sure that you are loosing precision long before level 31. At level 16 everything appears to check out.
I don't think any single set of 20 lines of code has ever (in the best part of 30 years), given me as much trouble to wrap my brain around as this!
Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
"Science is about questioning the status quo. Questioning authority".
In the absence of evidence, opinion is indistinguishable from prejudice.
| [reply]
[d/l]
[select] |
|
|
|
|
| |
|
I verified level 4 after your update and it does not yet match - although close !
Gah. I See where this is headed. I add another ifelse clause to get level 4 right, and then another for level 5...
This has to be possible recursively! It's amazing how something so complex evolves out of such a simple process. It's also amazing how something with so many symmetries and essentially only one dissymmetry could be so hard to generate.
Part of the problem is the shear size the things grow to so rapidly. I wrote a simple png plotter that produces a png of a 1000x1000 chunk and allows you to 'move' around and zoom in--very slowly. But as the level increases, it becomes impossible to keep track of where you are, so it doesn't help much.
Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
"Science is about questioning the status quo. Questioning authority".
In the absence of evidence, opinion is indistinguishable from prejudice.
| [reply] |
|
|
