import sys, platform
print(sys.executable)
print(platform.machine())

/Users/haule/miniconda3-clean/envs/work311/bin/python
arm64

import numpy as np # because arrays are defined in numpy
# it is desired to use "import numpy as np" to not overwrite python functions with the same name.

def Mand(z0, max_steps):
    z = 0j  # no need to specify type.  
    # To initialize to complex number, just assign 0j==i*0
    for itr in range(max_steps):
        if abs(z)>2:
            return itr
        z = z*z + z0
    return max_steps

def Mandelbrot(ext, Nxy, max_steps):
    """
    ext[4]    -- array of 4 values [min_x,max_x,min_y,max_y]
    Nxy       -- int number of points in x and y direction
    max_steps -- how many steps we will try at most before we conclude the point is in the set
    """
    data = np.zeros((Nxy,Nxy), dtype=np.int32) # initialize a 2D dynamic array
    for i in range(Nxy):
        for j in range(Nxy):
            x = ext[0] + (ext[1]-ext[0])*i/(Nxy-1.)
            y = ext[2] + (ext[3]-ext[2])*j/(Nxy-1.)
            # creating complex number of the fly
            data[i,j] = Mand(x + y*1j, max_steps)  
    return data
# data now contains integers. 
# MandelbrotSet has value 1000, and points not in the set have value <1000.

data = Mandelbrot([-2,1,-1,1], 500, 1000)

print(np.shape(data))
print(data[0,0],data[-1,-1])

(500, 500)
1 2

import matplotlib.pyplot as plt #import pylab as plt    # plotting library
from matplotlib import cm
# next line is non-python line, jupyter notebook directive
%matplotlib inline     

ext=[-2,1,-1,1]
# pylab's function for displaying 2D image

plt.imshow(data.T,extent=ext,origin='lower',cmap=cm.coolwarm)
plt.colorbar()
plt.show()

help(plt.imshow)

Help on function imshow in module matplotlib.pyplot:

imshow(X: 'ArrayLike | PIL.Image.Image', cmap: 'str | Colormap | None' = None, norm: 'str | Normalize | None' = None, *, aspect: "Literal['equal', 'auto'] | float | None" = None, interpolation: 'str | None' = None, alpha: 'float | ArrayLike | None' = None, vmin: 'float | None' = None, vmax: 'float | None' = None, colorizer: 'Colorizer | None' = None, origin: "Literal['upper', 'lower'] | None" = None, extent: 'tuple[float, float, float, float] | None' = None, interpolation_stage: "Literal['data', 'rgba', 'auto'] | None" = None, filternorm: 'bool' = True, filterrad: 'float' = 4.0, resample: 'bool | None' = None, url: 'str | None' = None, data=None, **kwargs) -> 'AxesImage'
    Display data as an image, i.e., on a 2D regular raster.
    
    The input may either be actual RGB(A) data, or 2D scalar data, which
    will be rendered as a pseudocolor image. For displaying a grayscale
    image, set up the colormapping using the parameters
    ``cmap='gray', vmin=0, vmax=255``.
    
    The number of pixels used to render an image is set by the Axes size
    and the figure *dpi*. This can lead to aliasing artifacts when
    the image is resampled, because the displayed image size will usually
    not match the size of *X* (see
    :doc:`/gallery/images_contours_and_fields/image_antialiasing`).
    The resampling can be controlled via the *interpolation* parameter
    and/or :rc:`image.interpolation`.
    
    Parameters
    ----------
    X : array-like or PIL image
        The image data. Supported array shapes are:
    
        - (M, N): an image with scalar data. The values are mapped to
          colors using normalization and a colormap. See parameters *norm*,
          *cmap*, *vmin*, *vmax*.
        - (M, N, 3): an image with RGB values (0-1 float or 0-255 int).
        - (M, N, 4): an image with RGBA values (0-1 float or 0-255 int),
          i.e. including transparency.
    
        The first two dimensions (M, N) define the rows and columns of
        the image.
    
        Out-of-range RGB(A) values are clipped.
    
    cmap : str or `~matplotlib.colors.Colormap`, default: :rc:`image.cmap`
        The Colormap instance or registered colormap name used to map scalar data
        to colors.
    
        This parameter is ignored if *X* is RGB(A).
    
    norm : str or `~matplotlib.colors.Normalize`, optional
        The normalization method used to scale scalar data to the [0, 1] range
        before mapping to colors using *cmap*. By default, a linear scaling is
        used, mapping the lowest value to 0 and the highest to 1.
    
        If given, this can be one of the following:
    
        - An instance of `.Normalize` or one of its subclasses
          (see :ref:`colormapnorms`).
        - A scale name, i.e. one of "linear", "log", "symlog", "logit", etc.  For a
          list of available scales, call `matplotlib.scale.get_scale_names()`.
          In that case, a suitable `.Normalize` subclass is dynamically generated
          and instantiated.
    
        This parameter is ignored if *X* is RGB(A).
    
    vmin, vmax : float, optional
        When using scalar data and no explicit *norm*, *vmin* and *vmax* define
        the data range that the colormap covers. By default, the colormap covers
        the complete value range of the supplied data. It is an error to use
        *vmin*/*vmax* when a *norm* instance is given (but using a `str` *norm*
        name together with *vmin*/*vmax* is acceptable).
    
        This parameter is ignored if *X* is RGB(A).
    
    colorizer : `~matplotlib.colorizer.Colorizer` or None, default: None
        The Colorizer object used to map color to data. If None, a Colorizer
        object is created from a *norm* and *cmap*.
    
        This parameter is ignored if *X* is RGB(A).
    
    aspect : {'equal', 'auto'} or float or None, default: None
        The aspect ratio of the Axes.  This parameter is particularly
        relevant for images since it determines whether data pixels are
        square.
    
        This parameter is a shortcut for explicitly calling
        `.Axes.set_aspect`. See there for further details.
    
        - 'equal': Ensures an aspect ratio of 1. Pixels will be square
          (unless pixel sizes are explicitly made non-square in data
          coordinates using *extent*).
        - 'auto': The Axes is kept fixed and the aspect is adjusted so
          that the data fit in the Axes. In general, this will result in
          non-square pixels.
    
        Normally, None (the default) means to use :rc:`image.aspect`.  However, if
        the image uses a transform that does not contain the axes data transform,
        then None means to not modify the axes aspect at all (in that case, directly
        call `.Axes.set_aspect` if desired).
    
    interpolation : str, default: :rc:`image.interpolation`
        The interpolation method used.
    
        Supported values are 'none', 'auto', 'nearest', 'bilinear',
        'bicubic', 'spline16', 'spline36', 'hanning', 'hamming', 'hermite',
        'kaiser', 'quadric', 'catrom', 'gaussian', 'bessel', 'mitchell',
        'sinc', 'lanczos', 'blackman'.
    
        The data *X* is resampled to the pixel size of the image on the
        figure canvas, using the interpolation method to either up- or
        downsample the data.
    
        If *interpolation* is 'none', then for the ps, pdf, and svg
        backends no down- or upsampling occurs, and the image data is
        passed to the backend as a native image.  Note that different ps,
        pdf, and svg viewers may display these raw pixels differently. On
        other backends, 'none' is the same as 'nearest'.
    
        If *interpolation* is the default 'auto', then 'nearest'
        interpolation is used if the image is upsampled by more than a
        factor of three (i.e. the number of display pixels is at least
        three times the size of the data array).  If the upsampling rate is
        smaller than 3, or the image is downsampled, then 'hanning'
        interpolation is used to act as an anti-aliasing filter, unless the
        image happens to be upsampled by exactly a factor of two or one.
    
        See
        :doc:`/gallery/images_contours_and_fields/interpolation_methods`
        for an overview of the supported interpolation methods, and
        :doc:`/gallery/images_contours_and_fields/image_antialiasing` for
        a discussion of image antialiasing.
    
        Some interpolation methods require an additional radius parameter,
        which can be set by *filterrad*. Additionally, the antigrain image
        resize filter is controlled by the parameter *filternorm*.
    
    interpolation_stage : {'auto', 'data', 'rgba'}, default: 'auto'
        Supported values:
    
        - 'data': Interpolation is carried out on the data provided by the user
          This is useful if interpolating between pixels during upsampling.
        - 'rgba': The interpolation is carried out in RGBA-space after the
          color-mapping has been applied. This is useful if downsampling and
          combining pixels visually.
        - 'auto': Select a suitable interpolation stage automatically. This uses
          'rgba' when downsampling, or upsampling at a rate less than 3, and
          'data' when upsampling at a higher rate.
    
        See :doc:`/gallery/images_contours_and_fields/image_antialiasing` for
        a discussion of image antialiasing.
    
    alpha : float or array-like, optional
        The alpha blending value, between 0 (transparent) and 1 (opaque).
        If *alpha* is an array, the alpha blending values are applied pixel
        by pixel, and *alpha* must have the same shape as *X*.
    
    origin : {'upper', 'lower'}, default: :rc:`image.origin`
        Place the [0, 0] index of the array in the upper left or lower
        left corner of the Axes. The convention (the default) 'upper' is
        typically used for matrices and images.
    
        Note that the vertical axis points upward for 'lower'
        but downward for 'upper'.
    
        See the :ref:`imshow_extent` tutorial for
        examples and a more detailed description.
    
    extent : floats (left, right, bottom, top), optional
        The bounding box in data coordinates that the image will fill.
        These values may be unitful and match the units of the Axes.
        The image is stretched individually along x and y to fill the box.
    
        The default extent is determined by the following conditions.
        Pixels have unit size in data coordinates. Their centers are on
        integer coordinates, and their center coordinates range from 0 to
        columns-1 horizontally and from 0 to rows-1 vertically.
    
        Note that the direction of the vertical axis and thus the default
        values for top and bottom depend on *origin*:
    
        - For ``origin == 'upper'`` the default is
          ``(-0.5, numcols-0.5, numrows-0.5, -0.5)``.
        - For ``origin == 'lower'`` the default is
          ``(-0.5, numcols-0.5, -0.5, numrows-0.5)``.
    
        See the :ref:`imshow_extent` tutorial for
        examples and a more detailed description.
    
    filternorm : bool, default: True
        A parameter for the antigrain image resize filter (see the
        antigrain documentation).  If *filternorm* is set, the filter
        normalizes integer values and corrects the rounding errors. It
        doesn't do anything with the source floating point values, it
        corrects only integers according to the rule of 1.0 which means
        that any sum of pixel weights must be equal to 1.0.  So, the
        filter function must produce a graph of the proper shape.
    
    filterrad : float > 0, default: 4.0
        The filter radius for filters that have a radius parameter, i.e.
        when interpolation is one of: 'sinc', 'lanczos' or 'blackman'.
    
    resample : bool, default: :rc:`image.resample`
        When *True*, use a full resampling method.  When *False*, only
        resample when the output image is larger than the input image.
    
    url : str, optional
        Set the url of the created `.AxesImage`. See `.Artist.set_url`.
    
    Returns
    -------
    `~matplotlib.image.AxesImage`
    
    Other Parameters
    ----------------
    data : indexable object, optional
        If given, all parameters also accept a string ``s``, which is
        interpreted as ``data[s]`` if ``s`` is a key in ``data``.
    
    **kwargs : `~matplotlib.artist.Artist` properties
        These parameters are passed on to the constructor of the
        `.AxesImage` artist.
    
    See Also
    --------
    matshow : Plot a matrix or an array as an image.
    
    Notes
    -----
    
    .. note::
    
        This is the :ref:`pyplot wrapper <pyplot_interface>` for `.axes.Axes.imshow`.
    
    Unless *extent* is used, pixel centers will be located at integer
    coordinates. In other words: the origin will coincide with the center
    of pixel (0, 0).
    
    There are two common representations for RGB images with an alpha
    channel:
    
    -   Straight (unassociated) alpha: R, G, and B channels represent the
        color of the pixel, disregarding its opacity.
    -   Premultiplied (associated) alpha: R, G, and B channels represent
        the color of the pixel, adjusted for its opacity by multiplication.
    
    `~matplotlib.pyplot.imshow` expects RGB images adopting the straight
    (unassociated) alpha representation.

import time            # timeing
t0 = time.time()
data = Mandelbrot([-2,1,-1,1], 1000, 1000)
t1 = time.time()
print ('clock time: ',t1-t0,'s')

clock time:  16.060781002044678 s

#from numpy import * # because arrays are defined in numpy
from numba import njit  # This is the new line with numba

@njit   # this is an alias for @jit(nopython=True)
def Mand(z0, max_steps):
    z = 0j  # no need to specify type. 
    # To initialize to complex number, just assign 0j==i*0
    for itr in range(max_steps):
        if abs(z)>2:
            return itr
        z = z*z + z0
    return max_steps

@njit
def Mandelbrot2(ext, Nxy, max_steps):
    """
    ext[4]    -- array of 4 values [min_x,max_x,min_y,max_y]
    Nxy       -- int number of points in x and y direction
    max_steps -- how many steps we will try at most before we conclude the point is in the set
    """
    data = np.zeros((Nxy,Nxy)) # initialize a 2D dynamic array
    for i in range(Nxy):
        for j in range(Nxy):
            x = ext[0] + (ext[1]-ext[0])*i/(Nxy-1.)
            y = ext[2] + (ext[3]-ext[2])*j/(Nxy-1.)
            # creating complex number of the fly
            data[i,j] = Mand(x + y*1j, max_steps)  
    return data
# data now contains integers. 
# MandelbrotSet has value 1000, and points not in the set have value <1000.

import time            # timeing
t0 = time.time()
data = Mandelbrot2(np.array([-2,1,-1,1]), 1000, 1000)
t1 = time.time()
print ('clock time: ',t1-t0,'s')

clock time:  0.756770133972168 s

import numpy as np
from numba import njit, prange  # This is the new line with numba

@njit   # this is an alias for @jit(nopython=True)
def Mand(z0, max_steps):
    z = 0j  # no need to specify type. 
    # To initialize to complex number, just assign 0j==i*0
    for itr in range(max_steps):
        if abs(z)>2:
            return itr
        z = z*z + z0
    return max_steps

@njit(parallel=True)
def Mandelbrot3(data, ext, max_steps):
    """
    ext[4]    -- array of 4 values [min_x,max_x,min_y,max_y]
    Nxy       -- int number of points in x and y direction
    max_steps -- how many steps we will try at most before we conclude the point is in the set
    """
    Nx,Ny = np.shape(data) # 2D array should be already allocated we get its size
    for i in prange(Nx):
        for j in range(Ny):    # note that we used prange instead of range.
                                # this switches off parallelization of this loop, so that
                                # only the outside loop over i is parallelized.
            x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.)
            y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.)
            # creating complex number of the fly
            data[i,j] = Mand(x + y*1j, max_steps)  
# data now contains integers. 
# MandelbrotSet has value 1000, and points not in the set have value <1000.

import time            # timeing
data = np.zeros((1000,1000))
t0 = time.time()
Mandelbrot3(data, np.array([-2,1,-1,1]), 1000)
t1 = time.time()
print ('clock time: ',t1-t0,'s')

clock time:  0.2283637523651123 s

import matplotlib.cm as cm

plt.imshow(np.log(data.T), extent=[-2,1,-1,1], cmap=cm.coolwarm, origin='lower')
plt.colorbar()
plt.show()

CPP_code="""
#include "pybind11/pybind11.h"
#include "pybind11/numpy.h"
#include "pybind11/stl.h"
#include <cstdint>

namespace py = pybind11;
using namespace std;

void mand(py::array_t<double>& data, int Nx, int Ny, int max_steps, const vector<int>& ext)
{
  auto dat = data.mutable_unchecked<2>();
  #pragma omp parallel for
  for (int i=0; i<Nx; i++){
    for (int j=0; j<Ny; j++){
      dat(j,i) = max_steps;
      double x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.);
      double y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.);
      complex<double> z0(x,y);
      complex<double> z=0;
      for (int itr=0; itr<max_steps; itr++){
        if (norm(z)>4.){
          dat(j,i) = itr;
          break;
        }
        z = z*z + z0;
      }
    }
  }
}

PYBIND11_MODULE(imanc,m){
  m.doc() = "pybind11 wrap for mandelbrot";
  m.def("mand", &mand);
}
"""

with open('imanc.cc', mode='w') as file:
    file.write(CPP_code)
#file.close()

cmd = 'ls -ltr'
!{cmd}

total 3528
-rw-r--r--@ 1 haule  staff  588532 Jan  2  2024 Interactive MandelbrotSet.html
-rw-r--r--@ 1 haule  staff  548683 Jan 11  2024 Introduction to Comp Phys 488.html
drwxr-xr-x  6 haule  staff     192 Jan 22  2025 pybind11
-rwxr-xr-x  1 haule  staff    2009 Jan 30  2025 manp_dynamic.py
drwxr-xr-x  3 haule  staff      96 Aug 31 16:21 anaconda_projects
-rw-r--r--  1 haule  staff      72 Jan 19 20:17 Untitled.ipynb
-rw-r--r--  1 haule  staff      72 Jan 19 20:17 Untitled1.ipynb
-rw-r--r--  1 haule  staff      72 Jan 19 20:30 Untitled2.ipynb
-rwxr-xr-x  1 haule  staff  210568 Jan 19 21:30 imanc.so
-rw-r--r--  1 haule  staff  140573 Jan 19 23:27 Interactive MandelbrotSet.ipynb
-rw-r--r--  1 haule  staff  286988 Jan 20 00:45 Introduction to Comp Phys 488.ipynb
-rw-r--r--  1 haule  staff     828 Jan 20 00:46 imanc.cc

#cmd="g++ `python -m pybind11 --includes` -undefined dynamic_lookup -O3 -shared -std=c++11 -fPIC imanc.cc -o imanc.so"
cmd="g++-14 -O3 -Wall -shared -std=c++17 -fopenmp -fPIC `python -m pybind11 --includes` -undefined dynamic_lookup imanc.cc -o imanc.so"

!{cmd}

import imanc
print(help(imanc))

Help on module imanc:

NAME
    imanc - pybind11 wrap for mandelbrot

FUNCTIONS
    mand(...) method of pybind11_builtins.pybind11_detail_function_record_v1_system_libstdcpp_gxx_abi_1xxx_use_cxx11_abi_1 instance
        mand(arg0: typing.Annotated[numpy.typing.ArrayLike, numpy.float64], arg1: typing.SupportsInt, arg2: typing.SupportsInt, arg3: typing.SupportsInt, arg4: collections.abc.Sequence[typing.SupportsInt]) -> None

FILE
    /Users/haule/Teaching/488/Mandelbrot/imanc.so


None

import imanc  # We now import pybind11 module, produced with C++ code

data3 = np.ones((1000,1000))

t0 = time.time()
imanc.mand(data3, 1000, 1000, 1000, [-2,1,-1,1])
t1 = time.time()

print('pybind11: walltime: ', t1-t0)

pybind11: walltime:  0.20749211311340332

import matplotlib.cm as cm

plt.imshow(-np.log(data3), extent=ext, cmap=cm.coolwarm, origin='lower') 
plt.show()

Installation¶

Mandelbrot set¶

First checking the environment¶

Optional: Using C++ with `pybind11` to speed up the code¶

Optional installation of C++. These instructions are MAC specific¶

Return to `pybind11` to speed up the code¶

Popularity of programing languages¶

Installation¶

Mandelbrot set¶

First checking the environment¶

Optional: Using C++ with pybind11 to speed up the code¶

Optional installation of C++. These instructions are MAC specific¶

Return to pybind11 to speed up the code¶

Popularity of programing languages¶

Optional: Using C++ with `pybind11` to speed up the code¶

Return to `pybind11` to speed up the code¶