Computing and plotting the Mandelbrot set in Excel …

… with Python, Numba and pyxll.

As promised in the previous post, this post will look at alternative procedures for calculating and plotting the Mandelbrot set with Python based code.

The spreadsheet presented here and the associated Python code may be downloaded from:
https://interactiveds.com.au/software/FastMandelbrot.zip

Note that pyxll is required to link the Python code to Excel.

Seven alternative Python routines were used to produce the images shown below, using code based on data available at the Github site: https://gist.github.com/jfpuget/60e07a82dece69b011bb

Input data for the two images, and times for the 7 routines to generate the data are shown in the screenshot below:

The second image is zoomed in by a factor of 100,000, which requires many more iterations to define the image, resulting in execution times 20-50 times longer.

The first function uses plain Python, which works very slowly, but illustrates the simplicity of the computation used to generate the complex images, with repeated looping through the procedure:

# 0 All python
def mandelbrot(z,maxiter):
    c = z
    for n in range(maxiter):
        if abs(z) > 2:
            return n
        z = z*z + c
    return maxiter

@xl_func
@xl_arg('width', 'int')
@xl_arg('height', 'int')
@xl_arg('maxiter', 'int')
def mandelbrot_set0(xmin,xmax,ymin,ymax,width,height,maxiter):
    r1 = np.linspace(xmin, xmax, width)
    r2 = np.linspace(ymin, ymax, height)
    res = (r1,r2,[mandelbrot(complex(r, i),maxiter) for r in r1 for i in r2])
    return res

The next two functions replace lists with Numpy arrays, but the result is even slower:

# 1 Relace lists with Numpy arrays
@xl_func
@xl_arg('width', 'int')
@xl_arg('height', 'int')
@xl_arg('maxiter', 'int')
def mandelbrot_set1(xmin,xmax,ymin,ymax,width,height,maxiter):
    r1 = np.linspace(xmin, xmax, width)
    r2 = np.linspace(ymin, ymax, height)
    n3 = np.empty((width,height))
    for i in range(width):
        for j in range(height):
            n3[i,j] = mandelbrot(r1[i] + 1j*r2[j],maxiter)
    res = (r1,r2,n3)
    return res

See the download file for the second function using Numpy arrays (which is even slower).

Using the Numba jit compiler greatly increases the speed of the code, and is very simple to use, just requiring the addition of the @njit decorator to the start of the code:

# 3 Compile code with Numba
@njit
def mandelbrot3(c,maxiter):
    z = c
    for n in range(maxiter):
        if abs(z) > 2:
            return n
        z = z*z + c
    return 0

@xl_func
@xl_arg('width', 'int')
@xl_arg('height', 'int')
@xl_arg('maxiter', 'int')
@njit
def mandelbrot_set3(xmin,xmax,ymin,ymax,width,height,maxiter):
    r1 = np.linspace(xmin, xmax, width)
    r2 = np.linspace(ymin, ymax, height)
    n3 = np.empty((width,height))
    for i in range(width):
        for j in range(height):
            n3[i,j] = mandelbrot3(r1[i] + 1j*r2[j],maxiter)
    res =  (r1,r2,n3)
    return res

Removing the square root operation from the code resulted in a further doubling of the speed, but replacing the complex number operations with two floats made almost no difference (see download file for details). Using the Numba fastmath function (function 7) gave a small increase in performance, but the @guvectorize() decorator resulted in a ten times increase in speed; a total speed improvement of about 1000 times for the second image:

# 6 Use guvectorize
@njit(int32(complex128, int64))
def mandelbrot6(c,maxiter):
    nreal = 0
    real = 0
    imag = 0
    for n in range(maxiter):
        nreal = real*real - imag*imag + c.real
        imag = 2* real*imag + c.imag
        real = nreal;
        if real * real + imag * imag > 4.0:
            return n
    return 0

@guvectorize([(complex128[:], int64[:], int64[:])], '(n),()->(n)',target='parallel')
def mandelbrot_numpy(c, maxit, output):
    maxiter = maxit[0]
    for i in range(c.shape[0]):
        output[i] = mandelbrot6(c[i],maxiter)
    

@xl_func
@xl_arg('width', 'int')
@xl_arg('height', 'int')
@xl_arg('maxiter', 'int')        
def mandelbrot_set6(xmin,xmax,ymin,ymax,width,height,maxiter):
    r1 = np.linspace(xmin, xmax, width, dtype=np.float64)
    r2 = np.linspace(ymin, ymax, height, dtype=np.float64)
    c = r1 + r2[:,None]*1j
    n3 = mandelbrot_numpy(c,maxiter)
    res =  (r1,r2,n3.T)
    return res

Having generated the Mandelbrot data, a graph can be generated in Excel using Matplotlib and the pyxll plot function:

 Plot image
@xl_func
@xl_arg('width', 'int')
@xl_arg('height', 'int')
@xl_arg('maxiter', 'int')
def mandelbrot_image(xmin,xmax,ymin,ymax,width=3,height=3,maxiter=80,cmap='hot', normval=0.3):
    dpi = 72
    img_width = dpi * width
    img_height = dpi * height
    x,y,z = mandelbrot_set6(xmin,xmax,ymin,ymax,img_width,img_height,maxiter)
    
    fig, ax = plt.subplots(figsize=(width, height),dpi=72)
    ticks = np.arange(0,img_width,3*dpi)
    x_ticks = np.round(xmin + (xmax-xmin)*ticks/img_width,4)
    plt.xticks(ticks, x_ticks)
    y_ticks = np.round(ymin + (ymax-ymin)*ticks/img_width, 4)
    plt.yticks(ticks, y_ticks)
    
    norm = colors.PowerNorm(normval)
    ax.imshow(z.T,cmap=cmap,origin='lower',norm=norm) 
    pyxll.plot(fig)
    return "Mandelbrot_set6"

The appearance of the image can be varied with the two chart properties “cmap” and colors.PowerNorm. There are 166 named color maps, which have been listed on the spreadsheet using the lengthy code below:

@xl_func()
def colormaps():
        return plt.colormaps()

The appearance of any colormap can be seen in the chart simply by entering the index number in cell B18. Similarly the PowerNorm value (between 0 and 1) may be entered in cell B19.

The chart is created by entering the Mandelbrot_Image function in any cell, as shown in the first two images of this post. More details of the plotting function, and examples of the colormaps, will be given in the next post.

Posted in Charts, Charts, Excel, Link to Python, Maths, Newton, NumPy and SciPy, PyXLL, UDFs | Tagged , , , , , , , | 1 Comment

Animations from VBA

A convenient way to create an animation in Excel is to create on-sheet formulas or user defined functions to generate the required data, then use VBA to iterate through the range of input values so that the chart or image will update for each iteration. A problem with this approach is that often the image will not redraw, while the code is running.

Searching for a solution to this problem, there were surprisingly few discussions on this problem, and most of those I did find I couldn’t get to work, but eventually I found a solution using ActiveWindow.SmallScroll, which has the twin benefits of being very simple, and obvious how it works. Example VBA code is shown below.

Edit 2 Dec. 2021: Modifications to the code for generating the images resulted in very short cycle times, so I added some code to “sleep” between cycles, so the cycle time was at least 1 second. Sleep is a Windows function, so has to be “declared” at the top of the code module. I have also added code to time each iteration, so the sleep period is reduced as the time for generating the image increases. An added benefit of this code is that it has removed the white screen flashing between each iteration, so the resulting video is much smoother.

#If VBA7 Then ' Excel 2010 or later
    Public Declare PtrSafe Sub Sleep Lib "kernel32" (ByVal Milliseconds As LongPtr)
#Else ' Excel 2007 or earlier
    Public Declare PtrSafe Sub Sleep Lib "kernel32" (ByVal Milliseconds As Long)
#End If

Sub Zoomin()
Dim Startfact As Double, inc As Double, steps As Long, Fact As Double
Worksheets("Sheet4").Activate
With Application.ActiveSheet
    Startfact = Range("Start").Value
    inc = Range("Factor").Value
    steps = Range("steps").Value
   
    Fact = Startfact
    For i = 1 To steps
        STime = Timer()
        Range("Scale").Value = Fact
        Range("Iteration") = i
        ActiveWindow.SmallScroll down:=1
        ActiveWindow.SmallScroll up:=1
        Fact = Fact * inc
        TimeDiff = (Timer() - STime) * 1000
        If TimeDiff < 950 Then
            Sleep (1000 - TimeDiff)
        Else
            Sleep (50)
        End If
    Next i
End With
End Sub

An example animation generated with this routine (plus a Python user defined function to generate the images) is shown below. More details of the code for generating the images will be provided in a later post:

Zoom in x2 24 times

With the faster image generation I have added a second video zooming in 48 times, with a final magnification of 2.8E+14 times, which is about as far as you can get with the resolution allowed with 64 bit float values:

Zoom in x2 48 times
Posted in Animation, Excel, Link to Python, Maths, Newton, PyXLL, UDFs, VBA | Tagged , , , , , | 1 Comment

Tanh-Sinh Quadrature v. 5.03

Following some recent comments Graeme Dennes has released the latest version the Tanh-Sinh Quadrature spreadsheet with some corrections to the test function documentation.

The new file is located (as before) at:

Tanh-Sinh Quadrature

For more information on the last major release see Numerical Integration With Tanh-Sinh Quadrature v 5.0. For more background information and numerous examples search this site for Tanh-Sinh, or select Numerical Integration from the categories drop down.

As always, if you have any questions or comments, please leave a comment below.

Posted in Excel, Maths, Newton, Numerical integration, UDFs, VBA | Tagged , , , | 5 Comments

Some new old songs

Every so often I check YouTube to see what new old music has been posted, and today I found:

Haitian Fight Song, here played by Danny Thompson in Norway in 1968

and Bert Jansch and Danny Thompson playing Thames Lighterman (great pictures as well as great music):

Bert Jansch composition. Recorded at the BBC for John Peel’s Night Ride, broadcast Dec 18, 1968. Apologies for poor sound quality: taken from 50 year old 1/4″ 4-track mono tape running at 3 3/4 ips. The title refers to Pentangle roadie Bobby Cadman, whose previous occupation had been Thames Lighterman. Bert’s song “One For Jo” is also about him, addressed to Bobby’s wife.

Posted in Bach | Tagged , , , | Leave a comment

Using Matplotlib from Excel with pyxll

The pyxll documentation has many examples of plotting in Excel using Matplotlib and other packages, but I find the multiple options confusing and hard to follow, so this post works through the examples in the Matplotlib Users Guide tutorial. The sample spreadsheet and python code can be downloaded from:

Matplotlib Tute.zip

Note that all the examples below were taken from the manual for Release 2.0.2. The current release is 3.4.3, and includes several additional examples, and considerably more background information.

The first example plots a simple line graph using hard coded data. I have added a modified version that will transfer the data from a selected row on the spreadsheet.

import matplotlib.pyplot as plt
from matplotlib import colors

import pyxll
from pyxll import xl_func, xl_arg, xl_return
from pyxll import  plot

@xl_func
def MPLTute_1():
    # Create the figure
    fig = plt.subplots()[0]
    # Plot the data
    plt.plot([1,2,3,4])
    plt.plot([0,2,4,6])
    plt.ylabel('some numbers')
    # Display the figure in Excel
    plot(fig)

@xl_func
@xl_arg('data','float[]')
def MPLTute_1a(data):
    fig, ax = plt.subplots()
    ax.plot(data)
    ax.set(ylabel='some numbers')
    plot(fig)

The first function follows the code in the tutorial as closely as possible, with the addition of a second plotted line. The second function, in addition to reading the data to be plotted from the spreadsheet, follows the coding used in the pyxll examples more closely.

The next example plots XY data, read from the spreadsheet as a list of lists. The string ‘ro’ plots the line as red circles for each data point (see the manual for details):

@xl_func
@xl_arg('data','float[][]')
def MPLTute_2(data):
    fig, ax = plt.subplots()
    ax.plot(data[0], data[1], 'ro')
    ax.set_xbound(0, 6)
    ax.set_ybound(0, 20)
    ax.set(ylabel='some numbers')
    plot(fig)

Function 3 plots 3 lines, each with a different format string. The lines are hard coded functions of the range “t”. The lines are created with a single .plot, as a sequence of 3 sets of X values, Y values, line format string:

@xl_func
def MPLTute_3():
    t = np.arange(0., 5., 0.2)
    fig, ax = plt.subplots()
    ax.plot(t, t, 'r--', t, t**2, 'bs', t, t**3, 'g^')
    plot(fig)

Function 4 plots one of the lines from the previous example, but instead of using a single string to format the line, a dictionary is passed from the spreadsheet, allowing multiple format properties to be specified:

@xl_func
@xl_arg('props','dict<str, var>')
def MPLTute_4(props):
    x = np.arange(0., 5., 0.2)
    y = x**2
    fig, ax = plt.subplots()
    ax.plot(x, y, **props)
    plot(fig)

The spreadsheet includes a list of available properties, but see the manual for full details:

Function 5 plots a histogram, using the plt.hist method:

@xl_func
def MPLTute_5():
    np.random.seed(19680801)
    mu, sigma = 100, 15
    x = mu + sigma * np.random.randn(10000)
    
    fig, ax = plt.subplots()
    # the histogram of the data
    # Changed from the tutorial example: "n, bins, patches = plt.hist(x, 50, normed=1, facecolor='g', alpha=0.75)" 
    # because use of normed now raises an error.
    n, bins, patches = plt.hist(x, 50, density=1, facecolor='g', alpha=0.75)
 
    plt.xlabel('Smarts')
    plt.ylabel('Probability')
    plt.title('Histogram of IQ')
    plt.text(60, .025, r'$\mu=100,\ \sigma=15$')
    plt.axis([40, 160, 0, 0.03])
    plt.grid(True)
    
    plot(fig)

The comment noting the change from the code in the tutorial refers to the document for Release 2.02. The current tutorial (Release 3.4.3) has been corrected.

Function 6 illustrates the addition of text and graphics to the graph:

@xl_func
def MPLTute_6():
    fig, ax = plt.subplots()
    t = np.arange(0.0, 5.0, 0.01)
    s = np.cos(2*np.pi*t)
    line, = plt.plot(t, s, lw=2)
    plt.annotate('local max', xy=(2, 1), xytext=(3, 1.5),
    arrowprops=dict(facecolor='black', shrink=0.05),)
    plt.ylim(-2,2)

    plot(fig)

Finally Functions 7 and 7a illustrate the use of different axis types, and returning multiple graphs from a single function:

@xl_func
@xl_arg('scaletype', 'int')
def MPLTute_7(scaletype):
    fig, ax = plt.subplots()
    from matplotlib.ticker import NullFormatter # useful for `logit` scale
    # Fixing random state for reproducibility
    np.random.seed(19680801)
    # make up some data in the interval ]0, 1[
    y = np.random.normal(loc=0.5, scale=0.4, size=1000)
    y = y[(y > 0) & (y < 1)]
    y.sort()
    x = np.arange(len(y))
    # plot with selected axes scale
    if scaletype == 1:
        # plt.figure(1)
        # # linear
        # plt.subplot(221)
        plt.plot(x, y)
        plt.yscale('linear')
        plt.title('linear')
        plt.grid(True)
        # log
    elif scaletype == 2:
        # plt.subplot(222)
        plt.plot(x, y)
        plt.yscale('log')
        plt.title('log')
        plt.grid(True)
        # symmetric log
    elif scaletype == 3:
        # plt.subplot(223)
        plt.plot(x, y - y.mean())
        plt.yscale('symlog', linthreshy=0.01)
        plt.title('symlog')
        plt.grid(True)
        # logit
    elif scaletype == 4:
        # plt.subplot(224)
        plt.plot(x, y)
        plt.yscale('logit')
        plt.title('logit')
        plt.grid(True)
    else:
        return 'scaletype must be between 1 and 4'
    # Format the minor tick labels of the y-axis into empty strings with
    # `NullFormatter`, to avoid cumbering the axis with too many labels.
    plt.gca().yaxis.set_minor_formatter(NullFormatter())
    # Adjust the subplot layout, because the logit one may take more space
    # than usual, due to y-tick labels like "1 - 10^{-3}"
    # Not required, Excel version returns only 1 chart
    # plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=0.25,
    # wspace=0.35)

    plot(fig)


@xl_func
def MPLTute_7a():
    # Multiplot version
    fig, ax = plt.subplots()
    from matplotlib.ticker import NullFormatter # useful for `logit` scale
    # Fixing random state for reproducibility
    np.random.seed(19680801)
    # make up some data in the interval ]0, 1[
    y = np.random.normal(loc=0.5, scale=0.4, size=1000)
    y = y[(y > 0) & (y < 1)]
    y.sort()
    x = np.arange(len(y))
    # plot with various axes scales
    
    fig = plt.figure(1)
    # # linear
    plt.subplot(221)
    plt.plot(x, y)
    plt.yscale('linear')
    plt.title('linear')
    plt.grid(True)
    # log
    plt.subplot(222)
    plt.plot(x, y)
    plt.yscale('log')
    plt.title('log')
    plt.grid(True)
    # symmetric log
    plt.subplot(223)
    plt.plot(x, y - y.mean())
    plt.yscale('symlog', linthreshy=0.01)
    plt.title('symlog')
    plt.grid(True)
    # logit
    plt.subplot(224)
    plt.plot(x, y)
    plt.yscale('logit')
    plt.title('logit')
    plt.grid(True)

    # Format the minor tick labels of the y-axis into empty strings with
    # `NullFormatter`, to avoid cumbering the axis with too many labels.
    plt.gca().yaxis.set_minor_formatter(NullFormatter())
    # Adjust the subplot layout, because the logit one may take more space
    # than usual, due to y-tick labels like "1 - 10^{-3}"
    plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=0.25,
    wspace=0.35)

    plot(fig)

Posted in Charts, Charts, Drawing, Excel, Link to Python, PyXLL, UDFs | Tagged , , , , , , | 1 Comment