Contour plots with Excel and Matplotlib – 2 

The procedures discussed in the previous post required the contour data to be arranged on a regular rectangular grid, with the data points listed in a 2D array. As far as I know there is no built in alternative in Excel, but Matplotlib provides functions that allow contour plots to be generated from data at random locations. Examples of two alternative methods are given at: Contour plot of irregularly spaced data

I have modified the code to plot the data in Excel, and to return two separate plots:

@xl_func
def contour_examples(out=1):
    
    np.random.seed(19680801)
    npts = 200
    ngridx = 100
    ngridy = 200
    x = np.random.uniform(-2, 2, npts)
    y = np.random.uniform(-2, 2, npts)
    z = x * np.exp(-x**2 - y**2)

    # fig, (ax1, ax2) = plt.subplots(nrows=2)
    if out == 1:
        fig, ax1 = plt.subplots()

        # -----------------------
        # Interpolation on a grid
        # -----------------------
        # A contour plot of irregularly spaced data coordinates
        # via interpolation on a grid.

        # Create grid values first.
        xi = np.linspace(-2.1, 2.1, ngridx)
        yi = np.linspace(-2.1, 2.1, ngridy)

        # Linearly interpolate the data (x, y) on a grid defined by (xi, yi).
        triang = tri.Triangulation(x, y)
        interpolator = tri.LinearTriInterpolator(triang, z)
        Xi, Yi = np.meshgrid(xi, yi)
        zi = interpolator(Xi, Yi)

        # Note that scipy.interpolate provides means to interpolate data on a grid
        # as well. The following would be an alternative to the four lines above:
        # from scipy.interpolate import griddata
        # zi = griddata((x, y), z, (xi[None, :], yi[:, None]), method='linear')

        ax1.contour(xi, yi, zi, levels=14, linewidths=0.5, colors='k')
        cntr1 = ax1.contourf(xi, yi, zi, levels=14, cmap="RdBu_r")

        fig.colorbar(cntr1, ax=ax1)
        ax1.plot(x, y, 'ko', ms=3)
        ax1.set(xlim=(-2, 2), ylim=(-2, 2))
        ax1.set_title('grid and contour (%d points, %d grid points)' %
                    (npts, ngridx * ngridy))
        pyxll.plot(fig)
        return 'Contour example1'
    else:

    # ----------
    # Tricontour
    # ----------
    # Directly supply the unordered, irregularly spaced coordinates
    # to tricontour.
        fig, ax2 = plt.subplots()

        ax2.tricontour(x, y, z, levels=14, linewidths=0.5, colors='k')
        cntr2 = ax2.tricontourf(x, y, z, levels=14, cmap="RdBu_r")

        fig.colorbar(cntr2, ax=ax2)
        ax2.plot(x, y, 'ko', ms=3)
        ax2.set(xlim=(-2, 2), ylim=(-2, 2))
        ax2.set_title('tricontour (%d points)' % npts)

        plt.subplots_adjust(hspace=0.5)
        # plt.show()
        pyxll.plot(fig)
        return 'Contour example2'

The resulting plots are identical to those from the original code, except that the axes are not plotted to scale:

Interpolate on a grid
Use Matplotlib tricontour

The tricontour option has been used with the data from the simple FEA analysis used in the previous post, allowing the input data to include all the mid-side nodes. The code has also been modified to allow plotting to scale, and use of a Matplotlib colormap. The resulting plot is shown alongside the original Strand7 plot:

In the revised code the data is input as three 1D arrays, the plot width and height may be specified, or set height = 0 to plot to scale, and the number of contour levels, and any colormap may be specified:

@xl_func
@xl_arg('x', 'numpy_array', ndim = 1)
@xl_arg('y', 'numpy_array', ndim = 1)
@xl_arg('z', 'numpy_array', ndim = 1)
@xl_arg('margins', 'numpy_array', ndim = 1)
@xl_arg('levels', 'int')
@xl_arg('title', 'str')
def contour_plot(x, y, z, width = 10., height = 0., ms = 0, title = '', margins = [], levels = 8, cmap = 'RdBu_r'):
    # ----------
    # Tricontour
    # ----------
    # Directly supply the unordered, irregularly spaced coordinates
    # to tricontour.
    xmin = np.min(x)
    xmax = np.max(x)
    ymin = np.min(y)
    ymax = np.max(y)
    if height == 0: height = (width-1) * (ymax-ymin)/(xmax-xmin)
    fig, ax = plt.subplots(figsize=(width, height), dpi=72)
    if margins == []:
        plt.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.1)
    else:
        left = 0.1 if np.isnan(margins[0]) else margins[0] 
        right = 0.9 if np.isnan(margins[1]) else margins[1]
        bottom = 0.1  if np.isnan(margins[2]) else margins[2]
        top = 0.9 if np.isnan(margins[3]) else margins[3]
        plt.subplots_adjust(left=left, right=right, bottom=bottom, top=top)
   
    ax.tricontour(x, y, z, levels=levels, linewidths=0.0, colors='k')
    cntr2 = ax.tricontourf(x, y, z, levels=levels, cmap=cmap)

    fig.colorbar(cntr2, ax=ax)
    ax.plot(x, y, 'ko', ms=ms)
    xmin = np.min(x)
    xmax = np.max(x)
    ymin = np.min(y)
    ymax = np.max(y)
    ax.set(xlim=(xmin, xmax), ylim=(ymin, ymax))
    ax.set_title(title)

    pyxll.plot(fig)
    return 'Contour plot'

A more complex example illustrates a problem with the tricontour function when plotting over an area with concave segments. The image below shows vertical deflections from a Strand7 analysis of a retaining wall with piled foundations:

The Matplotlib plot with the same data inserts an additional triangular area in front of the wall:

Limiting the plot to the region to the right of the wall face avoids this problem, and generates a plot closely matching the Strand7 output:

It seems that it is nor straightforward to create a contour plot on a shape with a concave boundary. The best information I could find was:

However, this may not solve all problems, particular where the problem is defined on an irregular domain. Also, in the case where the domain has one or more concave areas, the delaunay triangulation may result in generating spurious triangles exterior to the domain. In such cases, these rogue triangles have to be removed from the triangulation in order to achieve the correct surface representation. For these situations, the user may have to explicitly include the delaunay triangulation calculation so that these triangles can be removed programmatically. Under these circumstances, the following code could replace the previous plot code:

Response from Graham G at: https://stackoverflow.com/questions/9170838/surface-plots-in-matplotlib

See the Stackoverflow discussion for more details and suggested code.

Posted in Arrays, Charts, Excel, Finite Element Analysis, Link to Python, Newton, NumPy and SciPy, PyXLL, Strand7, UDFs | Tagged , , , , | 1 Comment

Contour plots with Excel and Matplotlib

Contour plots are widely used for a great variety of purposes, but facilities for producing them in Excel are very limited, and those that are provided are very well hidden. In this post I will look at some simple examples from the Python Matplotlib library, linking the input and resulting plots to Excel using pyxll, and compare with Excel’s capabilities for producing contour plots. The following post will look at some of the much more extensive facilities in Matplotlib.

The first example plot (below) is based on the Matplotlib documentation at:
Matplotlib Contourf example
modified to accept input of the X and Y data and the colour bar data from Excel:

The Python code accepts input of the X and Y data and the data required to generate the contour colours and ranges:

@xl_func()
@xl_arg('levels', 'float[]')
@xl_arg('colors', 'numpy_array<var>', ndim = 2)
def contour_plot_ex1(Xvals, Yvals, levels, colors):
    # Use Numpy meshgrid to convert the 1D input lists to 2D grids
    X, Y = np.meshgrid(Xvals, Yvals)
    
    # Calculate the Z array with Numpy functions
    Z1 = np.exp(-X**2 - Y**2)
    Z2 = np.exp(-(X - 1)**2 - (Y - 1)**2)
    Z = (Z1 - Z2) * 2

    # If the first element of colors is a string, extract the first column and convert to a list
    if isinstance(colors[0,0], str):
        colors = colors[:,0].tolist()
    else:
        colors = colors.tolist()
    fig,ax=plt.subplots(1,1)
    # Generate the contour plot with contourf, using contour bands specified in the levels list, and colours specified in colors
    cp = ax.contourf(X, Y, Z, levels=levels, extend='both', colors= colors) 
    
    fig.colorbar(cp) # Add a colorbar to a plot
    ax.set_title('Filled Contours Plot')
    # Send to Excel with pyxll.plot
    pyxll.plot(fig)
    return "Contour Plot Example 1"

Input for the ‘levels’ and ‘colors’ arguments is shown in the screenshot below:

The X and Y data is listed in single column ranges (or in this case the same range for X and Y). ‘levels’ is a single column range listing the values for each contour boundary, and ‘colors’ may be a single column range with color codes, or a three or four column range specifying the color components (see Matplotlib documentation for details).

The contourf function can also be used as an alternative for plotting the Mandelbrot set (although it is much slower in generating the plot than the procedures used in previous posts):

@xl_func()
@xl_arg('width', 'int')
@xl_arg('height', 'int')
@xl_arg('maxiter1', 'int')
@xl_arg('itfactor', 'int')
@xl_arg('levels', 'float[]')
@xl_arg('colors', 'float[][]')
def contour_mandelplot(xc,xrange,yc,yrange, maxiter1, zoomfact, itfactor, width, height, levels, colors):
    dpi = 72
    img_width = (dpi * width)
    img_height = np.int(dpi * height*1.5)
    xrange0 = xrange
    xrange = xrange/zoomfact
    yrange = yrange/zoomfact
    xmin = xc - xrange/2
    xmax = xmin + xrange
    ymin = yc - yrange/2
    ymax = ymin + yrange
    
    if zoomfact > 1:
        maxiter = np.int(np.log(xrange0/xrange)*itfactor+maxiter1)
        if maxiter < maxiter1: maxiter = maxiter1
    else:
        maxiter = maxiter1
    X,Y,Z = mandelbrot_set5(xmin,xmax,ymin,ymax,img_width,img_height,maxiter)
    
    fig,ax=plt.subplots(1,1)
    cp = ax.contourf(X, Y, Z, levels=levels, extend='both', colors= colors) 
    
    ax.set_title('Mandelbrot Contour Plot')
    pyxll.plot(fig)
    return maxiter 

A more practical example is plotting finite element analysis results, where we will compare Matplotlib alternatives and Excel contour plots with output from the FEA program Strand7. The first FEA example is a very simple 2D plane strain analysis with 6 rectangular plates, restrained horizontally along the two vertical edges, vertically and horizontally along the base, with a downward deflection of 10 mm applied to the top left corner. The Strand7 vertical deflection results are shown below:

Strand7 Y deflection contour output

Each plate in the Strand7 model is defined by 8 nodes (the corners and the midpoint of each side), but the Excel contour plot and the standard contourf input require a rectangular grid, with an equal number of points in each row, and for the Excel plot, equal spacing between each point. The contoured data is also required to be in a 2D grid. To satisfy those requirements the deflections were extracted along the top and bottom edge of each plate as shown below:

The Matplotlib code, modified to accept input of the Z data was:

@xl_func()
@xl_arg('levels', 'float[]')
@xl_arg('colors', 'float[][]')
def contour_plot2(Xvals, Yvals, Z, levels, colors):
    X, Y = np.meshgrid(Xvals, Yvals)
   
    fig,ax=plt.subplots(1,1)
    cp = ax.contourf(X, Y, Z, levels=levels, extend='both', colors= colors) 
    
    fig.colorbar(cp) # Add a colorbar to a plot
    ax.set_title('Filled Contours Plot')
    pyxll.plot(fig)
    return "Contour Plot Example 2"

and the resulting output:

Matplotlib contour plot

The same input data shown above can be used to generate a contour plot in Excel as follows:

  1. Select the XYZ data arranged as a single range (Range I4:N8 in the data screenshot below).
  2. Select the Insert tab, then the ‘Waterfall, Funnel, Stock, Surface, or Radar chart’ group
  3. Click on the 3rd icon from the ‘Surface’ selection (Contour)
Excel contour plot input

The resulting chart, with all default settings, is:

Excel contour plot with default settings

The Matplotlib results are reasonably close to the Strand7 plot, the main differences being:

  • The horizontal and vertical dimensions are not plotted to equal scale
  • The interpolation of the contours is more jagged, because the vertical mid-side nodes were ignored

Dealing with these differences will be discussed in the next post.

The Excel chart with default settings on the other hand is almost useless. As well as the scale issue:

  • The chart is inverted, with Y values increasing from top to bottom
  • There are only 3 contour bands, resulting in a greatly simplified plot
  • Graded shading of the plot, intended to improve the readability of 3D surface charts, is also applied by default to 2D contour plots, making contour bands difficult to interpret, and in this case almost unreadable

These formatting problems can be fixed, but the procedures for doing it are very far from obvious. An excellent article by John Peltier provides detailed instructions, and the Excel 2007 procedure is still applicable to Excel 365.

Contour and Surface Charts in Excel 2007

In addition to the suggestions in the link, the Y axis being listed in reverse order seems to cause plotting problems which can be fixed by:

  • Select the vertical axis
  • Select Axis Options (right hand icon at the top of the dialog)
  • Select ‘Series in reverse order’

The chart then looks like:

Excel contour plot after adjustment of formatting
Posted in Arrays, Charts, Charts, Excel, Finite Element Analysis, Link to Python, NumPy and SciPy, PyXLL, Strand7, UDFs | Tagged , , , , , | 1 Comment

Plotting Mandelbrot in Excel with Matplotlib

The functions used in the previous post generate an array of integers, the number of iterations for the function return value to exceed 2. This post looks in more detail at how this data can be plotted in Excel using the Matplotlib imshow function and pyxll.

Matplotlib.imshow documentation

First the Mandelbrot data is generated using the fastest of the functions described in the previous post, then new Matplotlib fig and ax objects are created:

@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, show_ticks = True):
    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)

Axis “tick” values are created at the required spacing and precision, or alternatively are set not to display:

 if show_ticks:
        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)
    else:
        plt.xticks([])
        plt.yticks([])

The imshow function requires cmap and norm arguments:

cmap str or Colormap, default: rcParams["image.cmap"] (default: 'viridis'): The Colormap instance or registered colormap name used to map scalar data to colors.

norm Normalize, optional: The Normalize instance used to scale scalar data to the [0, 1] range before mapping to colors using cmap. By default, a linear scaling mapping the lowest value to 0 and the highest to 1 is used. 

Norm is set using colors.PowerNorm:

Linearly map a given value to the 0-1 range and then apply a power-law normalization over that range.

ax.imshow is used to generate the image, then pyxll.plot to copy the image to Excel:

    norm = colors.PowerNorm(normval)
    ax.imshow(z.T,cmap=cmap,origin='lower',norm=norm) 
    pyxll.plot(fig)

The animations below step through all 166 of the “cmap” values with three different “norm” values (double click any image to display full screen):

Norm = 0.3:

Norm = 0.2:

Norm = 0.25:

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

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