Mandelbrot shallow tour

[#1]

JSL for the video.

picpath = ""; // "D:\Chales\pic\m";
//http://www.fractalforums.com/mandelbrot-and-julia-set/good-coordinates-for-mandelbrot-zooms/
//
//Points on the boundary of the Mandelbrot set generally have the most interesting orbits.  The easiest boundary points to compute are:
//
//* the spike along the negative real axis
//* the boundary of the main cardioid: r = (1 - cos(theta))/2, x = r*cos(theta)+0.25, y = r*sin(theta)
//* the boundary of the period 2 disk: r = 0.25, x = r*cos(theta)-1, y = r*sin(theta)
eps = 1.6; // slightly outside
deep = 10;
steps = 8000;
angle = 3;
scale = 2; // <<<<<<< video made with 20, very slow
xsize = 192 * scale;
ysize = 108 * scale;
diag = Sqrt( xsize ^ 2 + ysize ^ 2 );
mat = J( ysize, xsize, . );
New Window( "mandelbrot", bb = Border Box( Text Box( "wait" ) ) );
Wait( .2 );
num = 1e6;
For( t = 0, t < steps, t += 1,
	If( eps > 1e-9, eps *= .997 );
	deep += .3;
	angle += eps / 50;
	radius = (1 - Cos( angle )) / 2 + eps;
	x = radius * Cos( angle ) + 0.25;
	y = radius * Sin( angle );
	top = y - eps;
	bottom = y + eps;
	left = x - eps * xsize / ysize;
	right = x + eps * xsize / ysize;
	d = Floor( Min( 1e5, deep ) );
	Parallel Assign( {d = d}, mat[r( top, bottom ), c( left, right )] = Mandelbrot( d, 4, c, r ) );
	mat = mat / d;
	red = mat ^ .5;
	green = mat ^ 9;
	both = Loc( red :* green == 1 );
	red[both] = .1;
	green[both] = .4;
	blue = 1 - mat ^ .2;
	blue[both] = .1;
	img = New Image( "rgb", {red, green, blue} );
	(bb << child) << delete;
	bb << append( img );
	num += 1;
	If( picpath != "", img << saveimage( picpath || Char( num ) || ".png" ) );
	bb << updatewindow;
	Wait( 0 );
);

[David_Burnham]

I started out creating GKS programs in the 80s.  It would hours of microVAX compute time to render a fractal.  Thanks for sharing.

[Craige_Hales]

 Welcome, and thanks! The video took hours to render as well, but it is 4K resolution, and 1000s of frames.

[MathStatChem]

I have been having some fun with this script, doing some modifications and creating the image series to convert into a video to make my own mandelbrot tours.  

Question:  The Parallel Assign() function -- I'm looking for syntax reference on it?  I don't see it in the JMP scripting books, but there is an entry in the JMP Scripting Guide under Help.  It's not clear to me what you are doing with the Parallel Assign function in your script.  

Parallel Assign( {d = d}, mat[r( top, bottom ), c( left, right )] = Mandelbrot( d, 4, c, r ) );

what does the r(top,bottom) and c(left, right) do for you?  top, bottom, left, and right are determined just prior in the script to determine the coordinate boundaries of the frame you are creating.  I assume it is doing some sort of interpolation between top, bottom and left,right to get the y and x coordinate to use?  Is there a general description of this syntax somewhere?  

[Craige_Hales]

It looks like I never documented that, I'll see what I can do. 

Yes, it does a linear interpolation of 1..nrows into top..bottom. (similar for c, ncols, left, right).

If you omit the (top,bottom) bit, the r value in the expression will be 1..nrows.

The {d = d} part needs explaining too. Each thread gets its own copy of everything it needs to reference; the left side variable is the worker thread's name and the right side is the main thread's name. On a 4-cpu machine that assignment is done four times. If it happens to be a huge matrix you are copying in you might not see a performance improvement.

Mandelbrot needs no other data and scales nicely. Most real problems need to read from another array and maybe do some bookkeeping as well.

If you want to peek under the covers, try this:

m=J(10,10,0);
x=0;
parallelassign({x=x++},m[r,c]=x);
show(x,m);

My 8-cpu machine produces this:

x = 8;
m =
[ 0 1 2 3 4 5 6 7 0 1,
2 3 4 5 6 7 0 1 2 3,
4 5 6 7 0 1 2 3 4 5,
6 7 0 1 2 3 4 5 6 7,
0 1 2 3 4 5 6 7 0 1,
2 3 4 5 6 7 0 1 2 3,
4 5 6 7 0 1 2 3 4 5,
6 7 0 1 2 3 4 5 6 7,
0 1 2 3 4 5 6 7 0 1,
2 3 4 5 6 7 0 1 2 3];

You should not write code that depends on that pattern, of course.

Also notice the right hand side of the matrix assignment can be more complicated.

m=J(400,600,0);
parallelassign({},m[r(-1,1),c(-1,1)]=( 
    z=sqrt(r^2+c^2); 
    z=if(z<.707,.707,z)-.5; 
    1-z; // last item in ; sequence is result
));
m=m/max(m); // for heat color
newwindow("example",newimage(heatcolor(m,"spectral")));

"z" only exists inside each of the threads, a separate copy for each.

That flat area, dark red, was created by the if statement

And yes, you can use for statements. You should not try to do any form of I/O except writing to the log. Don't make display boxes, don't update display boxes, don't read or write disk files. If it doesn't crash, it will still be wrong.

You can use globals, but they use a mutex to serialize access and the code may be no faster than unthreaded code.

m=J(3,5,0);
global:words="";
parallelassign({},m[r(-1,1),c(-1,1)]=( 
    z=sqrt(r^2+c^2); 
    show(r+0,c+0,r+c,z); // r and c are special, use an expression
    global:words||=char(z);
    z=if(z<.707,.707,z)-.5; 
    1-z; // last item in ; sequence is result
));
show(words);

global:words is built in an unpredictable order, determined by how fast different threads run.

JSL is complicated (internally) and so is threading. The parallelassign function tries to keep everything sane (those copies, for example). If you find a way to work around something it is trying to keep you from doing, you may be damaging internal data structures when two threads access them at the same time, and that will cause JMP to crash at some point. Save often.

Oh...one final note about 1D matrices: use two subscripts! There are some confusing situations if you have a 1xN row vector and try to use a single subscript. Use R,C and ignore the constant R==1 value.

[MathStatChem]

Thanks for the detailed info on Parallel Assign()!  My brain is swimming with ideas on how I could use it.  

[MathStatChem]

Check out my video!  https://youtu.be/o9M5GZSKA6o

[Craige_Hales]

Cool, the first one is the hardest. I've settled in on a slightly different path using Blender (which usese ffmpeg internally) rather than ffmpeg directly. I give up some control over the ffmpeg parameters (which I didn't want to spend time learning) and get a pretty cool video sequence editor in exchange (which has its own significant learning curve).