import sys
import numpy as np
import scipy
import matplotlib as mpl
import numba

print("Python:", sys.version)
print("NumPy:", np.__version__)
print("SciPy:", scipy.__version__)
print("Matplotlib:", mpl.__version__)
print("Numba:", numba.__version__)
print("Installation appears to be working.")

Python: 3.11.14 | packaged by conda-forge | (main, Oct 22 2025, 22:56:31) [Clang 19.1.7 ]
NumPy: 2.3.5
SciPy: 1.17.0
Matplotlib: 3.10.8
Numba: 0.63.1
Installation appears to be working.

import numpy as np

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)) # 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)

import matplotlib.pyplot as plt
%matplotlib inline

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

plt.imshow(data.T, origin='lower', extent=ext)
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
t0 = time.perf_counter()
data = Mandelbrot([-2,1,-1,1], 1000, 1000)
t1 = time.perf_counter()
print('wall time:', t1-t0, 's')

wall time: 15.850245874957182 s

import numpy as np
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.

# The first call includes JIT compilation, so use a small call to compile first.
Mandelbrot2(np.array([-2,1,-1,1]), 10, 10)

t0 = time.perf_counter()
data = Mandelbrot2(np.array([-2,1,-1,1]), 1000, 1000)
t1 = time.perf_counter()
print('wall time after compilation:', t1-t0, 's')

wall time after compilation: 0.7151672090403736 s

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

@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 = data.shape # the 2D array is already allocated; obtain its dimensions
    for i in prange(Nx):
        for j in range(Ny):    # note that we used prange instead of range for parallel loop.
                                # 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.

# Compile the parallel function before timing it.
Mandelbrot3(np.zeros((10,10)), np.array([-2,1,-1,1]), 10)

data = np.zeros((1000,1000))
t0 = time.perf_counter()
Mandelbrot3(data, np.array([-2,1,-1,1]), 1000)
t1 = time.perf_counter()
print('parallel wall time after compilation:', t1-t0, 's')

parallel wall time after compilation: 0.22528633300680667 s

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

%%file mandc.cc

#include <iostream> // for printing to stdout
#include <complex>  // for complex numbers
#include <ctime>    // for measuring time
#include <vector>   // for using array/vector
#include <omp.h>    // openMP and measuing time with openMP
using namespace std;// everything from standard template library is available withouth prefix std
// std::cout and std::endl can be called cout, endl, etc...

int Mandelb(const complex<double>& z0, int max_steps)
{
  complex<double> z=0; // specify type where we start using it. (still need to specify type not like python with auto-typing)
  for (int i=0; i<max_steps; i++){ // specify int type inside the loop, so it is local to the loop. useful for reducing bugs.
    if (abs(z)>2.) return i;
    z = z*z + z0;
  }
  return max_steps;
}

int main()
{
  const int Nx = 1000; // constants are named cost rather than parameter
  const int Ny = 1000;
  int max_steps = 1000; // allowed to change later.
  double ext[]={-2,1,-1,1}; // this specifies fixed array of four numbers and initializes it.

  //double mand[Nx][Ny];
  
  vector<int> mand(Nx*Ny);  // allocating one dimensional array (vector) of size Nx*Ny
  // multidimensional dynamic arrays are not standard in C++. One can use pointers to allocate/deallocate memory,
  // and write one's own class interface for that. Or one has to use extension of C++ (blitz++ is excellent).
  // Unfortunately the standard C++ still does not support standard class for that.
  clock_t startTimec = clock();  // cpu time at the start 
  double start = omp_get_wtime(); // wall time at the start

  #pragma omp parallel for
  for (int i=0; i<Nx; i++){
    for (int j=0; j<Ny; j++){
      double x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.); // x in the interval ext[0]...ext[1]
      double y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.); // y in the interval ext[2]...ext[3]
      mand[i*Ny+j] = Mandelb(complex<double>(x,y), max_steps); // storing values in 2D array using 1D array
      //mand[i*Ny+j] = Mandelb(x+y*complex<double>(0.,1.), max_steps); // storing values in 2D array using 1D array
    }
  }

  clock_t endTimec = clock();
  double diffc = double(endTimec-startTimec)/CLOCKS_PER_SEC; // how to get seconds from cpu time
  double diff = omp_get_wtime()-start;                       // openMP time is already in seconds
  
  clog<<"clock time : "<<diffc<<"s"<<" with wall time="<<diff<<"s "<<endl; // printout of time
  for (int i=0; i<Nx; i++){
    for (int j=0; j<Ny; j++){
      double x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.);
      double y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.);
      cout<<x<<" "<<y<<" "<< 1./mand[i*Ny+j] << endl; // prinout of mandelbrot set
    }
  }
}

Overwriting mandc.cc

!{ls -l mandc.cc}

-rw-r--r--  1 haule  staff  2656 Aug 30 01:04 mandc.cc

from pathlib import Path

# Homebrew uses versioned names for GNU C++, for example g++-15.
candidates = [p for root in ('/opt/homebrew/bin', '/usr/local/bin')
              for p in Path(root).glob('g++-[0-9]*')]

if not candidates:
    raise FileNotFoundError('GNU C++ was not found. On macOS, install it with: brew install gcc')
    
gxx = str(max(candidates, key=lambda p: int(p.name.rsplit('-', 1)[1])))
print('Using', gxx)

cmd=f'"{gxx}" -fopenmp -O3 -o mandc mandc.cc'

!{cmd}

Using /opt/homebrew/bin/g++-15

!{time mandc > mand.dat}

clock time : 0.948198s with wall time=0.218726s 
mandc > mand.dat  1.47s user 1.08s system 114% cpu 2.220 total

def PlotFile(filename):
    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib.colors import LogNorm

    X,Y,Z = np.loadtxt(filename,unpack=True)
    ext=[X[0],X[-1],Y[0],Y[-1]]
    N = int(np.sqrt(len(Z)))
    dat = Z.reshape(N,N)
    print(ext)
    plt.imshow(dat.T, extent=ext, cmap=plt.cm.coolwarm, norm=LogNorm())
    plt.show()

PlotFile('mand.dat')

[np.float64(-2.0), np.float64(1.0), np.float64(-1.0), np.float64(1.0)]

%%file mandf.f90

INTEGER Function Mandelb(z0, max_steps)
  IMPLICIT NONE ! Every variable needs to be declared. It is very prudent to use that.
  COMPLEX*16, intent(in) :: z0
  INTEGER, intent(in)    :: max_steps
  ! locals
  COMPLEX*16 :: z
  INTEGER    :: i
  z=0.
  do i=1,max_steps
     if (abs(z)>2.) then ! |f(z)|>2 we know it will explode for large n, hence not part of the set.
        Mandelb = i-1      ! return how many iterations it took to know this is not part of the set.
        return 
     end if
     z = z*z + z0
  end do
  Mandelb = max_steps
  return
END Function Mandelb

program mand
  use omp_lib
  IMPLICIT NONE
  ! external function
  INTEGER :: Mandelb ! Need to declare the external function
  ! locals
  INTEGER :: i, j
  REAL*8  :: x, y
  COMPLEX*16  :: z0
  INTEGER, parameter :: Nx = 1000
  INTEGER, parameter :: Ny = 1000
  INTEGER, parameter :: max_steps = 1000
  REAL*8  :: ext(4) = (/-2., 1., -1., 1./) ! The limits of plotting
  REAL*8  :: mande(Nx,Ny)
  REAL*8  :: start, finish, startw, finishw
  
  call cpu_time(start)
  startw  = OMP_get_wtime()
  
  !$OMP PARALLEL DO  PRIVATE(j,x,y,z0)
  do i=1,Nx
     do j=1,Ny
        x = ext(1) + (ext(2)-ext(1))*(i-1.)/(Nx-1.) ! x \in [ext(1)...ext(2)]
        y = ext(3) + (ext(4)-ext(3))*(j-1.)/(Ny-1.) ! y \in [ext(3)...ext(4)]
        z0 = dcmplx(x,y)
        mande(i,j) = Mandelb(z0, max_steps)
     enddo
  enddo
  !$OMP END PARALLEL DO

  finishw = OMP_get_wtime()
  call cpu_time(finish)
  WRITE(0, '("clock time : ",f6.3,"s  wall time=",f6.3,"s")') finish-start, finishw-startw

  do i=1,Nx
     do j=1,Ny
        x = ext(1) + (ext(2)-ext(1))*(i-1.)/(Nx-1.)
        y = ext(3) + (ext(4)-ext(3))*(j-1.)/(Ny-1.)
        print *, x, y, 1./mande(i,j)
     enddo
  enddo
end program mand

Overwriting mandf.f90

cmd="gfortran -fopenmp -O3 -o mandf mandf.f90"

!{cmd}

!{time mandf > mandf.dat}

clock time :  0.932s  wall time= 0.210s
mandf > mandf.dat  1.81s user 0.05s system 140% cpu 1.325 total

PlotFile('mandf.dat')

[np.float64(-2.0), np.float64(1.0), np.float64(-1.0), np.float64(1.0)]

import imanc
imanc.mand(data, Nx, Ny, max_steps, ext)

%%file imanc.cc   

#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);
}

Overwriting imanc.cc

!{ls -l imanc.cc}

-rw-r--r--  1 haule  staff  828 Aug 30 01:05 imanc.cc

import subprocess
import sys
import sysconfig

includes = subprocess.check_output([sys.executable, '-m', 'pybind11', '--includes'], text=True).strip()
print(f'The includes we need for pybind11: {includes}\n')
extension_suffix = sysconfig.get_config_var('EXT_SUFFIX')
print(f'modules have suffix {extension_suffix}\n')

cmd = (f'"{gxx}" -fopenmp {includes} -undefined dynamic_lookup -O3 -shared -std=c++11 -fPIC imanc.cc -o imanc{extension_suffix}')
print(f'compilation comand: {cmd}')

!{cmd}

The includes we need for pybind11: -I/Users/haule/miniconda3-clean/envs/work311/include/python3.11 -I/Users/haule/miniconda3-clean/envs/work311/lib/python3.11/site-packages/pybind11/include

modules have suffix .cpython-311-darwin.so

compilation comand: "/opt/homebrew/bin/g++-15" -fopenmp -I/Users/haule/miniconda3-clean/envs/work311/include/python3.11 -I/Users/haule/miniconda3-clean/envs/work311/lib/python3.11/site-packages/pybind11/include -undefined dynamic_lookup -O3 -shared -std=c++11 -fPIC imanc.cc -o imanc.cpython-311-darwin.so

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

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

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

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

pybind11: walltime:  0.20748345798347145

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

%%file mandel.f90

!-----------------------------------------------------
! Produces Mandelbrot plot in the range [-ext[0]:ext[1]]x[ext[2]:ext[3]]
! It uses formula z = z*z + z0 iteratively until
! asb(z) gets bigger than 2 
! (deciding that z0 is not in mandelbrot)
! The value returned is 1/(#-iterations to escape)
!-----------------------------------------------------
SUBROUTINE Mandelb(data, ext, Nx, Ny, max_steps)
  IMPLICIT NONE ! Don't use any implicit names of variables!
  ! Function arguments
  REAL*8, intent(out) :: data(Nx,Ny)
  REAL*8, intent(in)  :: ext(4)              ! [xa,xb,ya,yb]
  INTEGER, intent(in) :: max_steps
  INTEGER, intent(in) :: Nx, Ny
  !f2py integer optional, intent(in)          :: max_steps=1000
  !  !f2py integer intent(hide), depend(data) :: Nx=shape(data,0)  ! it will be hidden automatically
  !  !f2py integer intent(hide), depend(data) :: Ny=shape(data,1)  ! it will be hidden automatically
  ! Local variables
  INTEGER    :: i, j, itt
  COMPLEX*16 :: z0, z
  REAL*8     :: x, y
  data(:,:) = max_steps
  !$OMP PARALLEL DO  PRIVATE(j,x,y,z0,z,itt)
  DO i=1,Nx
     DO j=1,Ny
        x = ext(1) + (ext(2)-ext(1))*(i-1.)/(Nx-1.)
        y = ext(3) + (ext(4)-ext(3))*(j-1.)/(Ny-1.)
        z0 = dcmplx(x,y)
        z=0
        DO itt=1,max_steps
           IF (abs(z)>2.) THEN
              data(i,j) = itt-1 !1./itt         ! result is number of iterations
              EXIT
           ENDIF
           z = z**2 + z0             ! f(z) = z**2+z0 -> z
        ENDDO
        !if (abs(z)<2) data(i,j) = max_steps
     ENDDO
  ENDDO
  !$OMP END PARALLEL DO
  RETURN
END SUBROUTINE Mandelb

Overwriting mandel.f90

!{f2py -c mandel.f90 --f90flags='-fopenmp' -m mandel -lgomp}

/Users/haule/miniconda3-clean/envs/work311/lib/python3.11/site-packages/numpy/f2py/f2py2e.py:737: VisibleDeprecationWarning: 
distutils has been deprecated since NumPy 1.26.x
Use the Meson backend instead, or generate wrappers without -c and use a custom build script
  builder = build_backend(
running build
running config_cc
INFO: unifying config_cc, config, build_clib, build_ext, build commands --compiler options
running config_fc
INFO: unifying config_fc, config, build_clib, build_ext, build commands --fcompiler options
running build_src
INFO: build_src
INFO: building extension "mandel" sources
INFO: f2py options: []
INFO: f2py:> /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11/mandelmodule.c
creating /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11
Reading fortran codes...
	Reading file 'mandel.f90' (format:free)
Post-processing...
	Block: mandel
			Block: mandelb
Applying post-processing hooks...
  character_backward_compatibility_hook
Post-processing (stage 2)...
Building modules...
    Building module "mandel"...
    Generating possibly empty wrappers"
    Maybe empty "mandel-f2pywrappers.f"
        Constructing wrapper function "mandelb"...
          data = mandelb(ext,nx,ny,[max_steps])
    Wrote C/API module "mandel" to file "/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11/mandelmodule.c"
INFO:   adding '/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11/fortranobject.c' to sources.
INFO:   adding '/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11' to include_dirs.
copying /Users/haule/miniconda3-clean/envs/work311/lib/python3.11/site-packages/numpy/f2py/src/fortranobject.c -> /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11
copying /Users/haule/miniconda3-clean/envs/work311/lib/python3.11/site-packages/numpy/f2py/src/fortranobject.h -> /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11
INFO:   adding '/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11/mandel-f2pywrappers.f' to sources.
INFO: build_src: building npy-pkg config files
/Users/haule/miniconda3-clean/envs/work311/lib/python3.11/site-packages/setuptools/_distutils/cmd.py:90: SetuptoolsDeprecationWarning: setup.py install is deprecated.
!!

        ********************************************************************************
        Please avoid running ``setup.py`` directly.
        Instead, use pypa/build, pypa/installer or other
        standards-based tools.

        This deprecation is overdue, please update your project and remove deprecated
        calls to avoid build errors in the future.

        See https://blog.ganssle.io/articles/2021/10/setup-py-deprecated.html for details.
        ********************************************************************************

!!
  self.initialize_options()
running build_ext
INFO: customize Compiler
INFO: customize Compiler using build_ext
INFO: get_default_fcompiler: matching types: '['gnu95', 'nag', 'nagfor', 'absoft', 'ibm', 'intel', 'gnu', 'g95', 'pg']'
INFO: customize Gnu95FCompiler
INFO: Found executable /opt/homebrew/bin/gfortran
INFO: customize Gnu95FCompiler
INFO: customize Gnu95FCompiler using build_ext
INFO: building 'mandel' extension
INFO: compiling C sources
INFO: C compiler: clang -DNDEBUG -fwrapv -O2 -Wall -fPIC -O2 -isystem /Users/haule/miniconda3-clean/envs/work311/include -arch arm64 -fPIC -O2 -isystem /Users/haule/miniconda3-clean/envs/work311/include -arch arm64 -I/usr/local/opt/openjdk@8/include

creating /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11
INFO: compile options: '-DNPY_DISABLE_OPTIMIZATION=1 -I/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11 -I/Users/haule/miniconda3-clean/envs/work311/lib/python3.11/site-packages/numpy/_core/include -I/Users/haule/miniconda3-clean/envs/work311/include/python3.11 -c'
INFO: clang: /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11/mandelmodule.c
INFO: clang: /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11/fortranobject.c
INFO: compiling Fortran sources
INFO: Fortran f77 compiler: /opt/homebrew/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops
Fortran f90 compiler: /opt/homebrew/bin/gfortran -fopenmp -fPIC -O3 -funroll-loops
Fortran fix compiler: /opt/homebrew/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -fopenmp -fPIC -O3 -funroll-loops
INFO: compile options: '-I/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11 -I/Users/haule/miniconda3-clean/envs/work311/lib/python3.11/site-packages/numpy/_core/include -I/Users/haule/miniconda3-clean/envs/work311/include/python3.11 -c'
INFO: gfortran:f90: mandel.f90
INFO: gfortran:f77: /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11/mandel-f2pywrappers.f
INFO: /opt/homebrew/bin/gfortran -Wall -g -Wall -g -undefined dynamic_lookup -bundle /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11/mandelmodule.o /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11/fortranobject.o /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/mandel.o /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq/src.macosx-11.0-arm64-3.11/mandel-f2pywrappers.o -L/opt/homebrew/Cellar/gcc/15.2.0_1/bin/../lib/gcc/current/gcc/aarch64-apple-darwin25/15 -L/opt/homebrew/Cellar/gcc/15.2.0_1/bin/../lib/gcc/current/gcc/aarch64-apple-darwin25/15/../../.. -L/opt/homebrew/Cellar/gcc/15.2.0_1/bin/../lib/gcc/current/gcc/aarch64-apple-darwin25/15/../../.. -lgomp -lgfortran -o ./mandel.cpython-311-darwin.so
ld: warning: ignoring duplicate libraries: '-lgfortran'
ld: warning: building for macOS-11.0, but linking with dylib '/opt/homebrew/opt/gcc/lib/gcc/current/libgomp.1.dylib' which was built for newer version 26.0
ld: warning: building for macOS-11.0, but linking with dylib '/opt/homebrew/opt/gcc/lib/gcc/current/libgfortran.5.dylib' which was built for newer version 26.0
ld: warning: building for macOS-11.0, but linking with dylib '/opt/homebrew/opt/gcc/lib/gcc/current/libquadmath.0.dylib' which was built for newer version 26.0
Removing build directory /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpviez7ocq

import mandel  # importing module created by f2py
import time

ext=[-2,1,-1,1]
Nx = Ny = 1000
max_steps = 1000
tc = time.process_time() # CPU time
tw = time.perf_counter() # wall time
data = mandel.mandelb(ext,Nx,Ny)
print('# wall time:', time.perf_counter()-tw, 's  CPU time:', time.process_time()-tc, 's')

# wall time: 0.21303945896215737 s  CPU time: 0.9238530000000011 s

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

Installation¶

Mandelbrot set¶

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

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

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

Install Homebrew¶

Install gnu-gcc and gfortran¶

(Optional) Install gnuplot and gsl¶

C++ and fortran with openMP implementation of the same algorithm¶

Return to `pybind11` to speed up the code¶

`f2py` as the alternative to `pybind11`¶

What did we learn from the optimization?¶

Popularity of programing languages¶

Installation¶

Mandelbrot set¶

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

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

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

Install Homebrew¶

Install gnu-gcc and gfortran¶

(Optional) Install gnuplot and gsl¶

C++ and fortran with openMP implementation of the same algorithm¶

Return to pybind11 to speed up the code¶

f2py as the alternative to pybind11¶

What did we learn from the optimization?¶

Popularity of programing languages¶

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

Return to `pybind11` to speed up the code¶

`f2py` as the alternative to `pybind11`¶