Friday, August 30, 2024

Fractal Dimension of Objects in Python

Target audience: Beginner
Estimated reading time: 5'

Are you finding it challenging to configure a convolutional neural network (CNN) for modeling 3D objects?

The complexity of a 3D object can pose a considerable challenge when tuning the parameters of a 3D convolutional neural network. Fractal analysis [ref 1] offers a way to measure the complexity of key features, volumes, and boundaries within the object, providing valuable insights that can help data scientists fine-tune their models for better performance.

Table of contents
Follow me on LinkedIn

What you will learn: How to evaluate the complexity of a 3-dimension object using fractal dimension index.

Notes

  • This article is a follow up on Fractal Dimension of Images in Python
  • Environments: Python  3.11,  Matplotlib 3.9.
  • Source code is available at  Github.com/patnicolas/Data_Exploration/fractal
  • To enhance the readability of the algorithm implementations, we have omitted non-essential code elements like error checking, comments, exceptions, validation of class and method arguments, scoping qualifiers, and import statements.

Introduction

As described in a previous article, Fractal Dimension of Images in Python - Overview , a fractal dimension is a measure used to describe the complexity of fractal patterns or sets by quantifying the ratio of change in detail relative to the change in scale [ref 2].
Among the various approaches to estimate the fractal dimension, from variation, structure function methods, root mean square and R/S analysis, we selected the box counting method because of its simplicity and visualization capability.

Point cloud
To evaluate our method for calculating the fractal dimension index of an object, it's necessary to simulate or represent a 3D object. This is achieved by generating a cluster of random data points across the x, y, and z axes, commonly referred to as a point cloud.

Box counting method
The box counting method [ref 3] is described in a previous article, Box-counting method

If N is the number of measurements units for our counting boxes or cubesand eps the related scaling factor, the fractal dimension index is computed as\[ D=- \displaystyle \lim_{\epsilon  \to 0} \frac{log(N)}{log(\epsilon)} \simeq - \frac{log(N)}{log(eps)} \ \ with \ N=r^3  \]We will use the height of the square box r as our measurement unit for images.
The fractal dimension index varies from 2 for very simple object to 3 for objects with complex pattern.

Implementation

Let's define a class, FractalDimObject, that encapsulates the calculation of fractal dimension and the creation of wrapping boxes. 

This class provides two constructors:
  1. The default constructor, __init__, which accepts a 3D array xyz and a threshold value near 1.
  2. An alternative constructor, build, which generates a 3D array of shape (size, size, size) to simulate a 3-dimension object.
class FractalDimObject(object):
    def __init__(self, xyz: np.array, threshold: float) -> None:
        self.xyz = xyz
        self.threshold = threshold


    @classmethod
    def build(cls, size: int, threshold: float) -> Self:
        _xyz = np.zeros((size, size, size))

        # Create a 3D fractal-like structure such as cube
        for x in range(size):            # Width
            for y in range(size):        # Depth
                for z in range(size):    # Height
                    if (x // 2 + y // 2) % 2 == 0:    # Condition for non-zero values
                        _xyz[x, y, z] = random.gauss(size//2, size)

        return cls(_xyz, threshold)

The alternative constructor generates a test array for evaluation purposes. It starts by initializing the array with values of 0.0, and then assigns random Gaussian-distributed values to a specific subset of the array

The values used in define the 3D object is visualized below.

Fig. 1 Visualization of the point cloud representing 3D object


The __call__ method performs the fractal dimension calculation and tracks the relationship between box counts and sizes in three stages:
  1. It determines the box sizes for non-zero elements in the array.
  2. It counts the number of boxes for each size.
  3. It applies a linear regression to the logarithms of the box sizes and their corresponding counts.
def __call__(self) -> (np.array, List[int], List[int]):
     # Step 1 Extract the sizes of array
     sizes = self.__extract_sizes()
     sizes_list = list(sizes)
     sizes_list.reverse()

     # Step 2 Count the number of boxes of each size
     counts = [self.__count_boxes(int(size)) for size in sizes_list]

     # Step 3 Fit the points to a line log(counts) = a.log(sizes) + b
     coefficients = np.polyfit(np.log(sizes), np.log(counts), 1)
     return -coefficients[0], sizes, counts

The __extract_sizes method is responsible for generating the box sizes, as detailed in the Appendix. The implementation of __count_boxes, which counts the wrapping boxes for a given size, follows a similar approach to the method used in calculating the fractal dimension of images.

def __count_boxes(self, box_size: int) -> int:
     sx, sy, sz = self.xyz.shape
     count = 0
        
      for i in range(0, sz, box_size):
          for j in range(0, sy, box_size):
             for k in range(0, sz, box_size):
                  # Wraps the non-zero values (object) with boxes
                 data = self.xyz[i:i+box_size, j:j+box_size, k:k+box_size]

                 if np.any(data):    # For non-zero values (inside object)
                    count += 1
     return count


Evaluation

Let's compute the fractal dimension of the array representing a 3D object with an initial 3D sampling grid 1024 x 1024 x 1024 

import math

grid_size = 1024      # Grid size 
threshold = 0.92
        
fractal_dim_object = FractalDimObject.build(grid_size, threshold)
coefficient, counts, sizes = fractal_dim_object()
print(coefficient)
      

Output: 2.7456

Finally let's plot the profile of box sizes vs. box counts.


Fig. 2 Plot of box sizes vs box counts size = exp(dim*counts)

The plot reflects the linear regression of log (size) and log (counts).

References


------------------
Patrick Nicolas has over 25 years of experience in software and data engineering, architecture design and end-to-end deployment and support with extensive knowledge in machine learning. 
He has been director of data engineering at Aideo Technologies since 2017 and he is the author of "Scala for Machine Learning", Packt Publishing ISBN 978-1-78712-238-3 and Geometric Learning in Python Newsletter on LinkedIn.


Appendix

Generate the array of box sizes used to compute the number of boxes to cover a given 3D object.

def __extract_sizes(self) -> np.array:
     # Remove values close to 1.0
     filtered = (self.xyz < self.threshold)
        
     # Minimal dimension of box size
     min_dim = min(filtered.shape)

     # Greatest power of 2 less than or equal to p
     n = 2 ** np.floor(np.log(min_dim) / np.log(2))
        
     # Extract the sizes
     size_x: int = int(np.log(n) / np.log(2))

     return np.arange(size_x, 1, -1) * 2

Monday, August 19, 2024

Performance Improvement in Numpy 2.x & Lapack 3.9

Target audience: Beginner
Estimated reading time: 4'
The performance of training geometric learning and complex deep learning models directly affects time to market and development costs. Numpy 2.x, combined with the accelerator and the latest BLAS/LAPACK libraries, can reduce execution time by up to 38%.


Table of contents
Overview
Linear Algebra
      Evaluation
      Implementation
       Evaluation
Follow me on LinkedIn

What you will learn: Enhanced performance in linear algebra computations using Numpy 2.x with the ILP64 accelerator and LAPACK 3.9.1 library

Notes

  • Environments: Python  3.11,  Matplotlib 3.9, Numpy 1.2.6/2.1, LAPACK 3.9.1
  • Source code is available at Github.com/patnicolas/Data_Exploration/numpy2
  • To enhance the readability of the algorithm implementations, we have omitted non-essential code elements like error checking, comments, exceptions, validation of class and method arguments, scoping qualifiers, and import statements.


Overview

Non-linear deep learning and geometric learning depend on geodesics, Riemannian metrics, and exponential maps to train models and predict outcomes. These operations are computationally intensive, making them ideal candidates for performance enhancements in the latest release of Numpy, version 2.1 [ref 1]. 
Linear algebra operations, on the other hand, rely on the BLAS and LAPACK libraries, which can be optimized for specific platforms [ref 2].

This study is performed on a MacOS M2 Maxlaptop, utilizing optimizations of various mathematical libraries like LAPACK and BLAS from the Accelerate framework [ref 3]. Specifically, we use the ILP64 interface with the LAPACK 3.9.1 library alongside Numpy 2.1 [ref 4].

In the following two sections, we will assess the performance enhancements in Numpy 2.1 when utilizing the ILP64 accelerator for the LINPACK library on macOS. The evaluation will cover two scenarios:
- Basic Linear Algebra operations
- Fast Fourier Transform

Linear Algebra

Let's compare the performance of Numpy 2.1 with the LINPACK accelerator 3.9.1 against Numpy 1.2.6 by performing the following computations:
  • Sigmoid function on a Numpy array of varying sizes
  • Summation of two 3D Numpy arrays
  • Multiplication of two 3D Numpy arrays
  • Dot product of two arrays.

Implementation

Let's start by creating keys for our linear algebra operators.

SIGMOID: AnyStr = 'sigmoid'
ADD: AnyStr = 'add'
MUL: AnyStr = 'mul'
DOT: AnyStr = 'dot'


We wraps the basic operations on Numpy array into the class LinearAlgebraEval with two constructors
  • Default, __init__ with input array and shape
  • build method to generate a random 3D array.
The four methods/operations—sigmoid, add, mul, and dot product—are straightforward and self-explanatory.

class LinearAlgebraEval(object):
    def __init__(self, x: np.array, shape: List[int]) -> None:
        self.x = x.reshape(shape)

    @classmethod
    def build(cls, size: int) -> Self:
        _x = np.random.uniform(0.0, 100.0, size)
        return cls(_x, [10, 100, -1])

    def sigmoid(self, a: float) -> np.array:
        return 1.0/(1.0+ np.exp(-self.x * a))

    def add(self, y: np.array) -> np.array:
        return self.x + y.reshape(self.x.shape)

    def mul(self, y: np.array) -> np.array:
        z = 2.0* self.x * y.reshape(self.x.shape)
        return z * z

    def dot(self, y: np.array) -> np.array:
        return np.dot(self.x.reshape(-1), y.reshape(-1))
   

The performance evaluation is conducted by an instance of the LinearAlgebraPerfEval class. The key method __call__ relies on a dictionary with operation as key and list of execution times for these operations.

class LinearAlgebraPerfEval(object):
    def __init__(self, sizes: List[int]) -> None:
        self.sizes = sizes

    def __call__(self) -> Dict[AnyStr, List[float]]:
        _perf = {}
        lin_algebra = LinearAlgebraEval.build(self.sizes[0])
        _perf[SIGMOID] = [LinearAlgebraPerfEval.__performance(lin_algebra, SIGMOID)]
        _perf[ADD] = [LinearAlgebraPerfEval.__performance(lin_algebra, ADD)]
        _perf[MUL] = [LinearAlgebraPerfEval.__performance(lin_algebra, MUL)]
        _perf[DOT] = [ LinearAlgebraPerfEval.__performance(lin_algebra, DOT)]

        for index in range(1, len(self.sizes)):
            # Alternative constructor
            lin_algebra = LinearAlgebraEval.build(self.sizes[index])

            # The method __record executes each of the four operations for
            # a variable size of Numpy array.

            _perf = LinearAlgebraPerfEval.__record(_perf, lin_algebra, SIGMOID)
            _perf = LinearAlgebraPerfEval.__record(_perf, lin_algebra, ADD)
            _perf = LinearAlgebraPerfEval.__record(_perf, lin_algebra, MUL)
            _perf = LinearAlgebraPerfEval.__record(_perf, lin_algebra, DOT)
        return _perf


The private, helper method, __record described in the Appendix, executes and time  each of the four operations for set of Numpy arrays of increasing size. The timing relies on the timeit decorator described in the Appendix.

Evaluation

Let's plot the execution time for these 4 operations with a size of 3D arrays (10, 100, -1) varying between 1million to 60 million values.

Numpy 1.2.6 with OpenBlas


Numpy 2.1 with ILP64- BLAS/LAPACK 3.9.1



Numpy 2.1 and ILP-64 have no impact on the performance of the Sigmoid but reduce the average time to execute and addition, multiplication and dot product on these large arrays by 25 to 40%


Fast Fourier Transform

Let's define a signal as a sum of 4 sinusoidal functions with 4 different frequency modes.\[ f(t)=\sum_{i=1}^{4} sin(2.\pi.f_{i}.t)\]The signal is then sampled in the interval [0, 1] with various number of data points.

Implementation

We encapsulate the generation, sampling, and extraction of the frequency spectrum within the `FFTEval` class. 
The constructor iteratively adds sine functions to implement the desired formula. 
The `compute` method samples the signal over the interval [0, 1] and calls the Numpy `fft.fftfreq` method to generate the frequency distribution. The timing of the execution in the method compute uses the timeit decorator.

class FFTEval(object):
    def __init__(self, frequencies: List[int], samples: int) -> None:
        self.F = 1.0/samples
        self.x = np.arange(0, 1.0, self.F)
        pi_2 = 2*np.pi
        
       # Composed the sinusoidal signal/function
        self.signal = np.sin(pi_2*frequencies[0]*self.x)
        if len(frequencies) > 1:
            for f in frequencies[1:]:
                self.signal += np.sin(pi_2*f*self.x)

    @timeit
    def compute(self) -> (np.array, np.array):
        num_samples = len(self.signal)

        # Select the positive amplitudes and frequencies
        num_half_samples = num_samples//2

        # Invoke the Numpy FFT function to extract frequencies
        freqs = np.fft.fftfreq(num_samples, self.F)
        _freqs = freqs[:num_half_samples]
        _amplitudes = np.abs(np.fft.fft(self.signal)[:num_half_samples]) / num_samples
        return _freqs, _amplitudes

The class FFTPerfEval implements the execution and collection of the execution time for the FFT with various frequency of samples (number of samples in the interval [0, 1])

class FFTPerfEval(PerfEval):
    def __init__(self, sizes: List[int]) -> None:
        super(FFTPerfEval, self).__init__(sizes)

    def __call__(self) -> List[float]:
        durations = []

        # Collect the execution time for each of the number of
        # samples defined in the constructor
        for samples in self.sizes:
            durations.append(FFTPerfEval.__compute(samples))
        return durations

    @staticmethod
    def __compute(sz: int) -> float:
        frequencies = [4, 7, 11, 17]
        fft_eval = FFTEval(frequencies, sz)
        # Time it
        return fft_eval.compute()


Evaluation

Numpy 1.2.6 with OpenBlas



Numpy 2.1 with BLAS/LAPACK 3.9.1



Relative performance improvement


The switch to Numpy 2.1 with ILP-64 on MacOS 14.6.1 shows an average improvement of 33%

References

[1Numpy
[2Lapack


------------------
Patrick Nicolas has over 25 years of experience in software and data engineering, architecture design and end-to-end deployment and support with extensive knowledge in machine learning. 
He has been director of data engineering at Aideo Technologies since 2017 and he is the author of "Scala for Machine Learning", Packt Publishing ISBN 978-1-78712-238-3 
and Geometric Learning in Python Newsletter on LinkedIn.

Appendix

Standard timing decorator in Python

def timeit(func):
    import time

    def wrapper(*args, **kwargs):
        start_time = time.time()
        func(*args, **kwargs)
        end_time = time.time()
        return end_time - start_time
    return wrapper


Methods to compute and record the execution time for basic linear algebra operations. The execution time is recorded through the decorator timer

 @staticmethod
    def __record(perf: Dict[AnyStr, List[float]],
                 lin_algebra: LinearAlgebraEval,
                 op: AnyStr) -> Dict[AnyStr, List[float]]:
        duration = LinearAlgebraPerfEval.__performance(lin_algebra, op)
        lst = perf[op]
        lst.append(duration)
        perf[op] = lst
        return perf

    @timeit
    @staticmethod
    def __performance(lin_algebra: LinearAlgebraEval, op: AnyStr) -> np.array:
        match op:
            case "sigmoid":
                return lin_algebra.sigmoid(8.5)
            case "add":
                return lin_algebra.add(lin_algebra.x)
            case "mul":
                return lin_algebra.mul(lin_algebra.x)
            case "dot":
                return lin_algebra.dot(lin_algebra.