Friday, July 17, 2015

Quickiebrot II

Here are some of those images saved by the script in Quickiebrot I.

CAUTION!  Flickering Images, don't watch if you are sensitive to things like flash photography.

Otherwise, SCROLL DOWN and get dizzy!

Or, better yet, go get my python script and make a better one! (I made this one horrible to encourage you to try it yourself.)




















Quickiebrot I

The story of the discovery of the Mandelbrot set is pretty amazing. It starts with "gee, what's the square root of negative one" and continues with "gee, what happens if I do this?" There is so much out there on the subject - go look!

I like the "Colors of Infinity" interview with Arthur C. Clarke (all over YouTube).

I also like this cool post by Ken Shirriff: 12-minute Mandelbrot: fractals on a 50 year old IBM 1401 mainframe.

This is 2014 but looks like vintage output! https://plus.google.com/photos/+KenShirriff/albums/6116650756716708097/6116651206478534962?pid=6116651206478534962&oid=106338564628446721517

See these YouTube links for really good explanations and how-to's for Julia sets and the Mandelbrot set.

Here is a quickie mandelbrot set imager that you can click-to-move, and zoom. It's just bare bones, but it gets you started. Note: It saves images, you might want to disable that!

Also, the standard slider in matplotlib.widgets is horizontal. I found a modified script to make it vertical and used it instead. That's why the script appears to be so long.

The full python script is given on this page. It's not fancy or optimal, but it gets you started.





Saturday, July 4, 2015

Relax Laplace! I

In a region with no free charge, the electrostatic potential obeys the Laplace equation:

`grad^2 phi = 0`

Jacobi iteration is a simple numerical technique to find solutions to problems like this. For example, if you have a 2D equal-spaced cartesian mesh of points, the numerical gradient

`(del^2 phi)/(del x^2) + (del^2 phi)/(del y^2)` can be approximated as:

(phi[i+1,j] + phi[i-1,j] - 2*phi[i,j]) + (phi[i,j+1] + phi[i,j-1] - 2*phi[i,j])

If that is set to zero, then one can iterate the expression:

next_phi[i,j] = 0.25 * (phi[i+1,j] + phi[i-1,j] +  phi[i,j+1] + phi[i,j-1])

It is slow, but here's a basic example. Scroll down to see cool ways to make it much faster!


Shown as a contour using mostly Matplotlib defaults:



1
2
CS = ax.contour(phi, origin = 'lower')
ax.clabel(CS, inline=1, fontsize=11)

The script to set up the problem (phido_me, and sphere) and then iterate is quite short. Instead of looping over all the indices in the array as you would in a compiled language like C or FORTRAN, in Python we try to just use existing NumPy methods like np.roll in line 29.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
import numpy as np
import matplotlib.pyplot as plt

nx, ny = 101, 101

phi   = np.zeros((ny, nx), dtype = 'float')
do_me = np.ones_like(phi, dtype='bool')

x0, y0, r0 = 40, 65, 12

x = np.arange(nx, dtype = 'float')[None,:]
y = np.arange(ny, dtype = 'float')[:,None]
rsq = (x-x0)**2 + (y-y0)**2

circle = rsq <= r0**2

phi[circle] = 1.0
do_me[circle] = False

do_me[0,:], do_me[-1,:], do_me[:,0], do_me[:,-1] = False, False, False, False

n, nper = 100, 100
phi_hold = np.zeros((n+1, ny, nx))
phi_hold[0] = phi

for i in range(n):
        
    for j in range(nper):
        phi2 = 0.25*(np.roll(phi,  1, axis=0) +
                     np.roll(phi, -1, axis=0) +
                     np.roll(phi,  1, axis=1) +
                     np.roll(phi, -1, axis=1) )

        phi[do_me] = phi2[do_me]

    phi_hold[i+1] = phi

change = phi_hold[1:] - phi_hold[:-1]

cmin, cmax = [], []
for thing in change:
    cmin.append(thing[do_me].min())
    cmax.append(thing[do_me].max())

The actual calculation (lines 26 to 36) takes about 2.2 seconds on my laptop. That's 100*100*101*101 multiplications plus three times that many sums, among other things, or about 6ns per flop. Since a lap-top-CPU-flop is more like about 1 or 1.5ns, there's definitely some room for improvement. (That GPU is just sitting there hungry for work, and could probably do it much faster!)

There's an upcoming post on lap-top-CPU-flops with NumPy; stay tuned!

 So, let's make the problem a little harder by making it 3D, Instead of ~10,000 mesh points we'll have ~1,000,000. At 8 bytes per point, that's beyond the cache, so it's going to be slower still.


and two slices below the sphere:


This took 525 seconds. That's 100*100*101*101*101 multiplications plus five summations, or about 9ns per flop. It's creeping up there! ,And we're still changing at about 3E-04 per 100 iterations. The next Relax Laplace! post will cover some ways to bring that back down.

The 3D version of the script for reference:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
nx, ny, nz = 101, 101, 101

phi   = np.zeros((nz, ny, nx), dtype = 'float')
do_me = np.ones_like(phi, dtype='bool')

x0, y0, z0, r0 = 40, 65, 50, 12

x = np.arange(nx, dtype = 'float')[None,None,:]
y = np.arange(ny, dtype = 'float')[None,:,None]
z = np.arange(ny, dtype = 'float')[:,None,None]
rsq = (x-x0)**2 + (y-y0)**2 + (z-z0)**2

sphere = rsq <= r0**2

phi[sphere] = 1.0
do_me[sphere] = False

do_me[0,:,:], do_me[-1,:,:] = False, False
do_me[:,0,:], do_me[:,-1,:] = False, False
do_me[:,:,0], do_me[:,:,-1] = False, False

n, nper = 100, 100
phi_hold = np.zeros((n+1, nz, ny, nx))
phi_hold[0] = phi

one_sixth = (1./6.)

for i in range(n):
        
    for j in range(nper):
        phi2 = one_sixth*(np.roll(phi,  1, axis=0) +
                          np.roll(phi, -1, axis=0) +
                          np.roll(phi,  1, axis=1) +
                          np.roll(phi, -1, axis=1) +
                          np.roll(phi,  1, axis=2) +
                          np.roll(phi, -1, axis=2) )

        phi[do_me] = phi2[do_me]

    print i

    phi_hold[i+1] = phi






Friday, July 3, 2015

"Don't use jet..."

Wikipedia
Jet is the worst colormap, except for all the others.

No, Winston Churchill was not talking about colormaps. He had more important things on his mind.

Starting from Jake VanderPlas' "How bad is your colormap" and this haiku tweet you can consider the matplotlib colormaps (and references therein) and the issues of luminance. For example, in the "jet" plot (default for imshow in matplotlib, and MatLab by the way) you can see artifacts. The bright blue line is meaningless. It just comes from a poor cmap definition.

But wait, there's more! Don't miss Climate Lab Book's End of the Rainbow, If We Assume's Planck vs WMAP Cosmic Ray Background OMG and this recent Scrap rainbow colour scales article in Nature!

These are some plots of the simulation in my "Relax Laplace!" post (in the works).

I also promise to add a post with some GOOD alternative maps.

Scroll down to the second (Black and White) version to understand the artifacts better.  Keep scrolling to see other gradient patterns.

The Seaborn project is worth reading about also.



...and...


1
2
3
4
5
x     = np.linspace(-1, 1, 200)
X, Y  = np.meshgrid(x, x)
theta = np.arctan2(Y, X)

theta[theta<0] += 2. * np.pi