Beefy Boxes and Bandwidth Generously Provided by pair Networks
P is for Practical
 
PerlMonks  

Fractal structure: PDL, memory, time

by FFRANK (Beadle)
on Jul 19, 2007 at 21:01 UTC ( #627622=perlquestion: print w/ replies, xml ) Need Help??
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

Comment on Fractal structure: PDL, memory, time
Download Code
Re: Fractal structure: PDL, memory, time
by BrowserUk (Pope) 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.
Re: Fractal structure: PDL, memory, time
by BrowserUk (Pope) 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.
      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

        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.

        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.

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: perlquestion [id://627622]
Approved by Corion
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others drinking their drinks and smoking their pipes about the Monastery: (9)
As of 2014-12-26 11:01 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    Is guessing a good strategy for surviving in the IT business?





    Results (171 votes), past polls