
![Super fast Python: Numba](super-fast-python-numba/super-fast-python-numba.jpg)

# Super fast Python (Part-5): Numba

This is the fifth and last post in the series on Python performance and Optimization. The series points out the utilization of inbuilt libraries, low-level code conversions, and other Python implementations to speed-up Python. The other posts included in this series are 

- (Part-1): [Why Python is slow?](https://santhalakshminarayana.github.io/blog/super-fast-python-why-python-slow)
- (Part-2): [Good practices to write fast Python code](https://santhalakshminarayana.github.io/blog/super-fast-python-good-practices)
- (Part-3): [Multi-processing in Python](https://santhalakshminarayana.github.io/blog/super-fast-python-multi-processing)
- (Part-4): [Use Cython to get speed as fast as C](https://santhalakshminarayana.github.io/blog/super-fast-python-cython)
- (Part-5): Use Numba to speed up Python Functions (this post)

In the last post about [Cython to speed-up Python code](https://santhalakshminarayana.github.io/blog/super-fast-python-cython), we discussed writing Python code in C-style, compiling that code separately into an object file, and using that generated file as an import directly into Python. But, not all people would feel comfortable writing C-style code or even some might not know C at all. So, to deal with such cases and get the performant efficient code to speed-up Python, one can use [Numba](https://numba.pydata.org/) instead. Numba translates Python code to machine code that executes almost as fast as C/C++ if optimized correctly.

## What is Numba?

Numba is a JIT (just-in-time) compiler that takes Python byte code and compiles it into machine code directly using [LLVM](https://llvm.org/) compiling mechanism. JIT is a type of interpreter that compiles frequently called code into machine code and caches that generated machine code to be used later for faster execution type. Here, Numba also takes Python code and generated machine code which the Python interpreter calls directly instead of interpreting and converting to machine code each time. Numba works best for numerical calculations, Array and Numpy operations, and loops. With Numba, we can write vectorized operations and parallelized loops to run on either CPU or GPU. 

Numba decorators are one of the many ways to invoke the JIT compilation. Numba provides different decorators to compile code in different modes and types, the common decorators used in Numba are:
- @jit - invoke JIT compilation for the provided function
- @njit - @jit decorator with enabling strict no-python mode
- @vectorize - convert normal functions into Numpy like **ufuncs**
- @guvectorize - generalized **ufuncs** for higher dimensional arrays
- @stencil - make a function behave as a kernel for a **stencil** like operation

Numba also provides different options to pass for some of these decorators to configure the JIT compilation behavior
- nopython
- parallel
- cache
- nogil
- fastmath
- boundscheck
- error_model
- cuda

## Numba @jit

**@jit** decorator takes the Python function that needs to machine code compiled. When we make a call to the function we provided to **@jit**, upon the first time calling, Numba compiles the function, caches the machine code, and this machine code is directly used for the execution. As compilation takes time, the first-time call to the function gives some latency. But, for consecutive function calls in the same runtime, just the cached machine code is used instead of re-compiling every time. 

Let's consider the following simple function *solve_expression* as an example. *solve_expression* takes some arguments, checks some conditions, and calculates the final polynomial expression. 

```python
def solve_expression(x, a, b, c, d):
    A, B, C, D = a, b, c, d
    if a > 10.1:
        A = 2 * a
    if 2.6 <= b < 8.3:
        B = b - 1/b
    if c > 4.5:
        C = 4
    if d < 9.0:
        D = d ** 2
    
    return A*(x**3) + B*(x**2) + C*(x) + D
```

Now, use the **@jit** decorator to compile this function into machine code as 

```python
from numba import jit

@jit
def solve_expression(x, a, b, c, d):
    A, B, C, D = a, b, c, d
    if a > 10.1:
        A = 2 * a
    if 2.6 <= b < 8.3:
        B = b - 1/b
    if c > 4.5:
        C = 4
    if d < 9.0:
        D = d ** 2
    
    return A*(x**3) + B*(x**2) + C*(x) + D
```

In the code snippet, we have imported the **jit** function decorator and decorated *solve_expression* with it.

```python
x, a, b, c, d = 2, 13, 1.2, 4, 7

res = solve_expression(x, a, b, c, d)
``` 

> Inspect the Intermediate Representation (IR) of the function using solve_expression.inspect_types()

As Numba **@jit** defers the JIT compilation until it encounters the first call to the function, the function call with arguments *solve_expression(x, a, b, c, d)* takes some time for executions. But, function calls later at this point will be fast. 

Now compare the speeds of the normal Python function and Numba JIT decorated function. Using *solve_expression.py_func()*, we can invoke the normal python function of this JIT decorated function.

```python
%% timeit
res = solve_expression.py_func(x, a, b, c, d)

'''Output
912 ns ± 3.47 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
'''
```

The normal Python function takes approx. 900 nanoseconds. 

```python
%%timeit
res = solve_expression(x, a, b, c, d)

'''Output
277 ns ± 1.36 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
'''
```

And the Numba version takes approx. 280 nanoseconds. The Numba version is 3x times faster than the pure Python function.

### Compilation options

For **@jit** decorator, we can pass multiple options to configure the compilation behavior

- **nopython**: if True, enables no-python mode making code execution without Python interpreter interference
- **nogil**: if True, releases the GIL (only when *nopython=True*) inside the compiled function, useful for concurrent execution such as threads
- **cache**: if True, store the compiled code in local storage and use this code whenever the function is called instead of re-compiling for every runtime
- **parallel**: if True, enables automatic parallelization
- **fastmath**: if True, uses faster math operations but less safe floating-point operations

Apart from these, there are some more options available. You can check all options at [@jit reference](https://numba.readthedocs.io/en/stable/reference/jit-compilation.html#jit-functions).

## Lazy and Eager compilation

Note that, as Python function arguments can take any time of arguments, Numba compiles the function for that specific type of the argument passed. If a new type of argument is passed to the function while calling, Numba re-compiles the code for that specific type.

### Lazy compilation

The above *solve_expression()* function takes any type of argument. If we provide a function like this Numba calculates the optimization steps to be done based on the argument types provided at the first function call. In this way, Numba infers the argument types and compiles the specific version of the same function for different types. 

Ex: if we change the argument types like this

```python
# previous values, x, a, b, c, d = 2, 13, 1.2, 4, 7
# new values
x, a, b, c, d = 3.9, 12, 5, 9.1, 14

res = solve_expression(x, a, b, c, d)
```

As the previous compiled function expects types for *x=int, a=int, b=float, c=int, d=int*, and in the latest function call we have changed some argument types. So, Numba re-compiles the function for new argument types. This mode of the compilation of code is called lazy compilation because Numba compiles for specific argument types only if it encounters them.

### Eager compilation

For function overloading, we can specify function signatures with argument types and return types in a list with the least significant precision at the top. [Numba types](https://numba.readthedocs.io/en/stable/reference/types.html) follow Numpy convention types with different precision levels.

```python
@jit(['int32(int32, int32)',
      'i4(int32, int64)', 
      '(f4, f8)',
      'f8(f4, f4)'])
def func1(a, b):
    return a + b
```

In the above function, we passed function signatures as a list of strings. The syntax is return type is specified first and argument types are specified after. It is allowed to have no return type specified. Numba will infer the return type automatically and use that specification. Calling the function with argument types not provided in the list raises an error.

## @njit or @jit(nopython=True)

Numba **@jit** operates in two modes *nopython* and *object*. In **nopython** mode, no interference of the Python interpreter is required and execution is very fast compared to normal mode. **object** mode is the same as calling function without **@jit**. 

Normally with the **@jit** decorator, Numba tries to compile in **nopython** mode. If any part of the code cannot be compiled due to the presence of code that is not supported by Numba like a heterogenous dictionary, some string methods, etc, then Numba fallbacks to **object** or normal Python mode for compilation. But still, it will improve the performance when loops are involved. If there is no code to optimize, **@jit** in object mode runs slower than the normal Python version as Numba compilation involves several function call overheads.

We can enable the strict *nopython* mode by passing the option to **@jit** as *@jit(nopython=True)*. Numba also provides a separate decorator for this option. *@njit* decorator is an alias for *@jit(nopython=True)*.

```python
@jit(nopython=True)
def f(a, b):
	return a + b

# or

@njit
def f(a, b):
	return a + b
```

With **nopython** mode, if any code is present that requires object mode compilation, Numba will raise an error. The primary goal is to write functions that can be implemented in strict no-python mode.

---

## Numba and Numpy

Numba is best for Numpy arrays and supports some [Numpy features](https://numba.readthedocs.io/en/stable/reference/numpysupported.html) in no-python mode.

```python
@njit(['int64[:](int64[:, :], int64[:, :])'])
def f(a, b):
    c = np.empty(a.shape[0], dtype='int64')
    
    for i in range(a.shape[0]):
        c[i] = a[i].sum() * b[i].sum()
    
    return c
```

The function takes 2 2D Numpy int64 arrays as arguments with the return type as an array of float64. The function calculates the sum of the product of the sum of each row of *a* and *b*.

```python
x, y, n = 1000, 1000, 1_000_000
l, h = 0, 100
a = np.random.randint(l, h, n).reshape(x, y)
b = np.random.randint(l, h, n).reshape(x, y)
```

If we compare both normal Numpy calculation of the above function and Numba compiled code, 

```python
%%timeit
# normal Python calculation
res = a.sum(axis=1) * b.sum(axis=1)

'''Output:
1.4 ms ± 7.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
'''
```

```python
%%timeit 
# Numba compiled function
res = f(a, b)

'''Output:
1.19 ms ± 15.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
'''
```

Numba compiled version takes **1.19 ms** which is slightly faster than normal calculation with Numpy **1.4ms**.

## @vectorize decorator

One of the main reasons for Numpy speed is that most of the Numpy functions are [ufunc](https://numpy.org/doc/stable/user/basics.ufuncs.html)s (universal functions) that are vectorized and implemented inside compiled layer of Numpy and hence the speed. We might run into a situation where we cannot find any existing Numpy functions for use and write a workaround by combining Numpy functions into a single operation. This new operation is not optimized and we might lose the speed that Numpy provides due to custom loops. To solve this problem, we can make any function as **ufunc** that comes with vectorization and speed. 

With [@vectorize](https://numba.readthedocs.io/en/stable/user/vectorize.html), Numba provides functionality to create custom vectorized universal functions. The universal function takes scalar values and returns a scalar value and these functions are applied over any Numpy arrays where array values are passed as single scalar values and an outside loop is automatically generated.

```python
@vectorize(['float64(int64, int64)'])
def v_expr(a, b):
    return (a**2) + (a*3) + (b/2) + 10
```

Here, we created a simple universal function that takes two scalar values and returns a calculated expression. 

```python
n = 1_000_000
l, h = 0, 100
a = np.random.randint(l, h, n)
b = np.random.randint(l, h, n)
```

The speed comparison of the Numpy expression and vectorized ufunc function gives 

```python
%%timeit
# general way of calculating Numpy expression
res = (a**2) + (a*3) + (b/2) + 10

'''Output:
8.66 ms ± 71.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
'''
```

```python
%%timeit
# calculate expression with ufunc
res = v_expr(a, b)

'''Output:
2.17 ms ± 96.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
'''
```

the optimized function **v_expr** takes **2.17 ms** which is 4x faster than the normal calculation with Numpy.

---

In this blog, we discussed basic Numba decorators and compared Numba compiled functions with normal Python functions. Explore other features of Numba like 
- [automatic parallelization](https://numba.readthedocs.io/en/stable/user/parallel.html)
- [loop parallelization with prange](https://numba.readthedocs.io/en/stable/user/parallel.html#explicit-parallel-loops)
- [stencil kernel implementation for arrays](https://numba.readthedocs.io/en/stable/user/stencil.html)
- [ahead-of-time compilation](https://numba.readthedocs.io/en/stable/user/pycc.html)

Remember that not every function can be passed to Numba as there are some limitations like 
- Numba only supports a subset of Python code like classes, multi-dimensional dictionaries, etc, which are not supported yet. 
- As object mode compilation takes more time than the normal Python mode in some cases, it is better to check the speed of both normal and compiled code execution speed.
- Support for external libraries like Pandas is not supported.
- As Numba re-implements some Numpy APIs, there may be different behavior expected.

---

### References
- [https://github.com/ContinuumIO/gtc2020-numba](https://github.com/ContinuumIO/gtc2020-numba)
- [https://www.nvidia.com/en-us/glossary/data-science/numba/](https://www.nvidia.com/en-us/glossary/data-science/numba/)
- [https://www.chrisvoncsefalvay.com/2019/03/23/jit-fast/](https://www.chrisvoncsefalvay.com/2019/03/23/jit-fast/)
- [https://towardsdatascience.com/numba-weapon-of-mass-optimization-43cdeb76c7da](https://towardsdatascience.com/numba-weapon-of-mass-optimization-43cdeb76c7da)
- [https://www.infoworld.com/article/3622013/speed-up-your-python-with-numba.html](https://www.infoworld.com/article/3622013/speed-up-your-python-with-numba.html)
- [https://coderzcolumn.com/tutorials/python/numba](https://coderzcolumn.com/tutorials/python/numba)
- [https://towardsdatascience.com/numpy-ufuncs-the-magic-behind-vectorized-functions-8cc3ba56aa2c](https://towardsdatascience.com/numpy-ufuncs-the-magic-behind-vectorized-functions-8cc3ba56aa2c)


![Super fast Python: Cython](super-fast-python-cython/super-fast-python-cython.jpg)

# Super fast Python (Part-4): Cython

This is the fourth post in the series on Python performance and Optimization. The series points out the utilization of inbuilt libraries, low-level code conversions, and other Python implementations to speed-up Python. The other posts included in this series are 

- (Part-1): [Why Python is slow?](https://santhalakshminarayana.github.io/blog/super-fast-python-why-python-slow)
- (Part-2): [Good practices to write fast Python code](https://santhalakshminarayana.github.io/blog/super-fast-python-good-practices)
- (Part-3): [Multi-processing in Python](https://santhalakshminarayana.github.io/blog/super-fast-python-multi-processing)
- (Part-4): Use Cython to get speed as fast as C (this post)
- (Part-5): [Use Numba to speed up Python Functions and Numeric calculations](https://santhalakshminarayana.github.io/blog/super-fast-python-numba)

In the last post, we discussed **multiprocessing** to optimize Python code by utilizing parallel computing with multiple CPU cores. Multi-processing is useful when we can split certain parts of code into parallel tasks and then execute them parallelly. 

Imagine when we don't have any parallelizable code and it is taking huge time in Python compared to **C/C++**, how to optimize the code then? One thing we can do is to convert the Python code into a low-level programming language like C/C++ and [embed that C/C++ code in Python](https://docs.python.org/3/extending/extending.html). Understanding and writing C/C++ code is hard for someone who is not familiar with C/C++. Fortunately, we can get the performance speed as fast as C/C++ by writing the Python code in [Cython](https://cython.org/) which is a superset of Python that provides functionality to write C-Extensions for Python.

## How Cython can improve Speed?

In our previous discussion on [Why Python is slow?](https://santhalakshminarayana.github.io/blog/super-fast-python-why-python-slow), we learned that the major speed issues arise in Python due to

- interpretation of generated bytecode
- dynamic types and their management

With Cython, we can take care of the above problems with both compiled code instead of interpretation and static typing instead of dynamic typing. 

Cython translates the Python code into C extension code and compiles the C code into an object code that can be imported directly into Python. Also, we can make some changes like adding static types that improve the execution speed drastically over dynamic typing. Cython also supports the usage of C/C++ libraries and functions that are super-fast compared to Python libraries and functions.

### Install Cython

[Install the Cython](https://cython.readthedocs.io/en/stable/src/quickstart/install.html) with Pip as follows

```bash
pip install Cython
```

> To compile the C code into an object file, we need a C compiler. Ubuntu comes with **gcc** by default. For other platforms like Windows, install the C compiler if not installed previously.

---

## Cython usage

There are multiple ways to use Cython like building manually,  using in Jupyter as an extension, or importing directly like a Python module without compilation using *pyximport*.

### Build a Cython module manually

We write the Cython code in a file with the extension *.pyx* instead of the normal Python extension *.py*.

The compiler translates '.pyx' Cython file into a '.c' C file and then compiles that C file to a sharable object file '.so' (or '.pyd' on Windows). We tell the compiler those build instructions and compilation options by writing a [setup.py](https://cython.readthedocs.io/en/stable/src/userguide/source_files_and_compilation.html#basic-setup-py) file.

### Cython compilation

- Translate the *.pyx* source code to a *.c* file with additional wrappers of Python extension code.
- Compile the *.c* file with a C compiler to a platform-specific shared object file *.so* that can be imported directly into Python.

```python:calculate_y_cython.pyx
def fx(x, a, b, c):
    d = (a + b) * c
    return a*x + b*(x**2) + c*(x**3) + d

def y(x, n, a, b, c):
    af = fx(x, a, b, c)**2
    k = abs(n//2 - x)
    bf = fx(k, a, b, c)
    return (af-bf)/(n-1 + 1e-12)

def calculate_y_cython(n):
    ys = []
    a, b, c = 2, 5, -4
    for i in range(n):
        ys.append(y(i, n, a, b, c))
    
    return ys
```

The above file *calculate_y.pyx* contains Cython code that looks the same as Python without any optimizations. Now to compile the above file, write a *setup.py* file as following

```python:setup.py
from distutils.core import setup
from Cython.Build import cythonize

setup(
    name='Calculate Expression Y',
    ext_modules=cythonize(
        "calculate_y_cython.pyx", 
        compiler_directives={"language_level": "3"},
    )
)
```

In the above *setup.py*, for the **setup()**, we have passed the optional name to our Cython module for *name* and the path to the *.pyx* file for the *ext_module* parameter.

Build the sharable object file '.so' using the following command in the terminal

```bash
python setup.py build_ext --inplace
or
python setup.py build_ext --inplace --quiet
```

This will generate a translated '.c' C file and a compiled '.so' file.

We can import the above compiled *calculate_y_cython* module directly into Python runtime.

```python:main.py
from time import perf_counter
from calculate_y_cython import calculate_y_cython

def fx(x, a, b, c):
    d = (a + b) * c
    return a*x + b*(x**2) + c*(x**3) + d

def y(x, n, a, b, c):
    af = fx(x, a, b, c)**2
    k = abs(n//2 - x)
    bf = fx(k, a, b, c)
    return (af-bf)/(n-1 + 1e-12)

def calculate_y_py(n):
    ys = []
    a, b, c = 2, 5, -4
    for i in range(n):
        ys.append(y(i, n, a, b, c))
    
    return ys

def main():
    n = 100000
    # Python implementation
    atime = perf_counter()
    res = calculate_y_py(n)
    print(f'Python Time: {perf_counter()-atime:.2}')
    
    atime = perf_counter()
    res = calculate_y_cython(n)
    print(f'Cython Time: {perf_counter()-atime:.2}')
    
if __name__=="__main__":
    main()

"""Output:
Python Time: 0.19
Cython Time: 0.16
"""
```

> At the time of writing, the latest Python version is Python 3.11 which is incredibly faster (30-60% in some cases) than earlier versions. So, to understand the Cython potential, I'm running the scripts in Python 3.8.10 on my old system with 8GB Intel(R) i5-8250U 1.60GHz CPU on HP Laptop 15-da0xxx.

In the above *main.py*, we have imported the *calculate_y_cython* module at line 2. If we look at the output to check the time taken for normal Python implementation and Cython version, they are **0.19** and **0.16** seconds respectively. The time difference is very low and we didn't gain much from Cython because we haven't done any optimization steps like 

- static typing
- limit calling Python's libraries
- reducing Python's PyObject usage

### Cython Annotations

Cython provides an easy way to check where we can optimize our Cython code by using Cython annotate feature. With the following command (*-a* denotes annotate and *-3* denotes **language_level** which is Python3), cython generates an HTML file that we can open in the browser to check for optimizable code. Another option is to pass the **annotate=True** parameter to *cythonize()* function call in **ext_modules** in *setup*.

```bash
cython -a calculate_y_cython.pyx -3
```

If we open the generated HTML file in the browser, it will look like this

![Cython Annotation:=:80](super-fast-python-cython/cython-annotations.png)

The more yellow lines the more interaction with the Python interpreter. Our goal should be converting as many yellow lines as to white lines that denote pure Cython code. We discuss the optimization part in a later section.

### Cython as an extension in Jupyter

Cython can be imported and used in Jupyter directly as an extension without the need for any additional build/compilation steps. 

First load the Cython extension using *%load_ext cython*, and then, for the cell that is to be Cythonized, use the magic command *%%cython* at the top of that cell as shown in the following image.

![Cython Jupyter:=:40](super-fast-python-cython/cython-jupyter.png)

We can show the interactive Cython annotations in Jupyter just like we have generated the HTML file above. To show annotations, pass annotate option to *%%cython -a* magic command.

### Import '.pyx' using pyximport

While developing or debugging, for each change in the '.pyx' file, running *setup.py* is a repetitive task and cumbersome. Instead, we can dynamically import '.pyx' to Python directly using [pyximport](https://cython.readthedocs.io/en/stable/src/userguide/source_files_and_compilation.html#pyximport) without any external build and compilation. *pyximport* takes care of compiling and building in the background without calling *cythonize()* internally. So, while importing the '.pyx' file, it will take some time to be imported as a regular Python module.

```python
import pyximport
pyximport.install(language_level=3)

from calculate_y_cython import calculate_y_cython
```

Though it is easy to work with Cython using **pyximport**, there are some [limitations with pyximport](https://cython.readthedocs.io/en/stable/src/userguide/source_files_and_compilation.html#limitations) and it is also not flexible as normal setup.

> It is not recommended to use **pyximport** while distributing Python packages and modules.

---

## Optimize Cython

In the previous section, we have seen in the Cython annotations HTML file that many areas in the code need to be optimized for more speed. There are several ways to improve speed like

- define static types
- use C-libraries and functions
- utilize OpenMP for parallel computing

### Define static types to variables

Since Python is a dynamic typing language, we can define C-like data types for Python variables in Cython. To declare C variables, prefix the **cdef** keyword to the variable declaration which is the same as declaring variables in C.

The syntax for variable declaration is

```python
cdef type variable_name = initilization_value {optional}
```

The **type** can be any of the acceptable C data types.

### Define function definitions

Python functions are defined using **def** and [C functions in Cython](https://cython.readthedocs.io/en/latest/src/userguide/language_basics.html#python-functions-vs-c-functions) are defined with the keyword **cdef**. The difference between **def** and **cdef** is that with the former declaration, it can be called from anywhere both local and external modules where as **cdef** functions are only module level. Also, Cython wraps the **def** function that is defined inside '.pyx' into a Python object. There is another declaration called **cpdef** which behaves as **def** when called from outside the module and behaves as **cdef** inside the module, so it is faster inside the same module function call.

Cython also provides support for writing function definitions just like C. We can define parameter types and return types. 

```python
cdef(or cpdef) type function_name(type parameter1, ...)
```

If no type is specified, the parameters and return values are treated as Python objects which need Python interpretation and they are slow.

Now, make some changes to *calculate_y_cython.pyx* with static type declarations and function definitions as,

> [Call the C functions in Cython](https://cython.readthedocs.io/en/stable/src/tutorial/external.html) inplace of Python functions that reduces the Python interaction.

```python:calculate_y_cython.pyx
ctypedef long int li
ctypedef long long int lli

cdef lli fx(li x, int a, int b, int c):
    cdef int d = (a + b) * c
    return a*x + b*(x**2) + c*(x**3) + d

cdef double y(li x, li n, int a, int b, int c):
    cdef:
        lli af = fx(x, a, b, c)**2
        li k = abs(n//2 - x)
        lli bf = fx(k, a, b, c)
    return (af-bf)/(n-1 + 1e-12)

cpdef calculate_y_cython(li n):
    ys = []
    cdef:
        int a = 2, b = 5, c = -4
        li i = 0
    for i in range(n):
        ys.append(y(i, n, a, b, c))
    
    return ys
```

> speed can be improved further by disabling bound checking(@cython.boundscheck(False)) and negative indexing(@cython.wraparound(False)) [compiler directive instructions](https://cython.readthedocs.io/en/stable/src/userguide/source_files_and_compilation.html#compiler-directives).

Check the annotations for the Cython code, most of the yellow lines are now changed to white and only list operations are dark yellow because lists are Python objects. In a later section, we discuss optimizing lists using C arrays.

![Optimized Cython Annotation:=:70](super-fast-python-cython/optimized_cython.jpg)

If we check the time for the optimized Cython code,

```python
%%timeit -n 100
n = 100000
res = calculate_y_cython(n)

"""Output:
2.04 ms ± 166 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
"""
```

it takes **0.002** seconds which is approximately **100x** times faster than the normal Python version that takes **0.19** seconds. By adding only static type and some tweaks, we made Python 100x faster. The powers of Cython are not limited to only static typing. We can further speed up the code with Numpy.

---

## 1000x times faster with MemoryView and Numpy

In the Cython annotation snippet above, we see that the list operations of **ys** are still yellow as *list()* is a python object and Cython cannot optimize Python objects interaction. 

To solve this, we can convert the list type to [memoryview arrays](https://cython.readthedocs.io/en/stable/src/tutorial/array.html), and operations like *append()* to indexing just like looping in C.

```python:calculate_y_cython.pyx
cimport cython
import numpy as np
cimport numpy as np

ctypedef long int li
ctypedef long long int lli

@cython.boundscheck(False)  # Deactivate bounds checking
@cython.wraparound(False)   # Deactivate negative indexing
cdef lli fx(li x, int a, int b, int c):
    cdef int d = (a + b) * c
    return a*x + b*(x**2) + c*(x**3) + d

@cython.boundscheck(False)
@cython.wraparound(False)
cdef double y(li x, li n, int a, int b, int c):
    cdef:
        lli af = fx(x, a, b, c)**2
        li k = abs(n//2 - x)
        lli bf = fx(k, a, b, c)
    return (af-bf)/(n-1 + 1e-12)

@cython.boundscheck(False)
@cython.wraparound(False)
cpdef calculate_y_cython(li n):
    cdef double[:] ys = np.empty(n, dtype=np.float64)
    cdef:
        int a = 2, b = 5, c = -4
        li i = 0
    for i in range(n):
        ys[i] = y(i, n, a, b, c)
    return ys
```

Before building the object code, we need to change [build instructions to link Numpy](https://cython.readthedocs.io/en/stable/src/tutorial/numpy.html) as a dependency like following

```python:setup.py
from distutils.core import setup
from distutils.extension import Extension
from Cython.Build import cythonize
import numpy

ext_modules = [
    Extension("calculate_y_function",
              sources=["calculate_y_cython.pyx"],
              libraries=["m"],
              compiler_directives={"language_level": "3"},
              )
]

setup(name="calculate Y function",
      ext_modules=cythonize(ext_modules),
      include_dirs=[numpy.get_include()])
```

The parameter *include_dirs* includes external libraries like Numpy here.

```python
%%timeit -n 1000
n = 100000
res = calculate_y_cython(n)

"""Output:
208 µs ± 11.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
"""
```

The fully optimized version of Cython runs at **0.0002** seconds which is approximately **1000x** faster than the normal Python code. 

The usage of Numpy and MemoryViews in Cython needs separate discussion and I will write about Classes, C-math, Numpy, MemoryViews, and OpenMP in the Advanced Cython series later.

---

Cython is a very powerful extension that we can use to speed up Python code. Sometimes it may not be possible to convert Python objects straight away like dictionaries (use C++ maps), so it's better to use Cython for repetitive tasks like loops, general functions with simple statements (like declaration and usage only), and math operations.

Learn more about Cython by referencing

- [Speeding up basic object operations in Cython
](http://blog.behnel.de/posts/tuning-basic-object-operations-in-cython.html)
- [Faster Python made easier with Cython’s pure Python mode](https://www.infoworld.com/article/3648539/faster-python-made-easier-with-cythons-pure-python-mode.html)
- [An Introduction to Just Enough Cython to be Useful](https://www.peterbaumgartner.com/blog/intro-to-just-enough-cython-to-be-useful/)
- [Cython notes and tips](https://nicolas-hug.com/blog/cython_notes)
- [(Github Issue) Improve cython interface with a more user-friendly compiling interface](https://github.com/cython/cython/issues/3974)


![Super fast Python: Good Practices](super-fast-python-good-practices/super-fast-python-good-practices.jpg)

# Super fast Python (Part-2): Good practices

In the earlier post on [why Python is slow?](https://santhalakshminarayana.github.io/blog/super-fast-python-why-python-slow), we discussed slowness is in Python due to its internal design of some essential components like GIL, dynamic typing, and interpretation.

In this blog, we discuss some good practices to speed up Python incredibly faster.

This is the second post in the series on Python performance and Optimization. The series points out the utilization of inbuilt libraries, low-level code conversions, and other Python implementations to speed-up Python. The other posts included in this series are 

- (Part-1): [Why Python is slow?](https://santhalakshminarayana.github.io/blog/super-fast-python-why-python-slow)
- (Part-2): Good practices to write fast Python code (this post)
- (Part-3): [Multi-processing in Python](https://santhalakshminarayana.github.io/blog/super-fast-python-multi-processing)
- (Part-4): [Use Cython to get speed as fast as C](https://santhalakshminarayana.github.io/blog/super-fast-python-cython)
- (Part-5): [Use Numba to speed up Python Functions and Numeric calculations](https://santhalakshminarayana.github.io/blog/super-fast-python-numba)

The following section describes the various good practices one can use to make Python super speed (up to 30% or more) without any external support like PyPy, Cython, Numpy, etc.

## Python good practices for super fast code

### Use built-in data structures and libraries

As Python data types are implemented directly in **C**, using the built-in types like list, map, and trees, compared to custom types we define, really helps the program to run faster.

Also, use built-in libraries for common algorithms like counting the duplicates, summing all list elements, finding the maximum element, etc, because these are all already written in **C** and compiled which makes these functions run faster than custom functions we write.

```python
from random import randint

rand_nums = [randint(1, 100) for _ in range(100000)]
```
Create 100000 random numbers between 1 and 100.

```python
%%timeit
cc = 0
for i in range(len(rand_nums)):
    cc += rand_nums[i]
```
```python:output
6.02 ms ± 94 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
```

Manually summing up the all numbers over running a loop takes approximately 6ms.

```python
%%timeit
cc = sum(rand_nums)
```
```python:output
332 µs ± 15 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
```

Using built-in *sum()* takes only 0.3ms approximately.

### Local Variables vs Global Variables

When we call a function (a routine), the system pauses the code execution at the call site in the current routine (say **main()**) where the call has been made and places the called function at the top of the call stack. Imagine if we have defined numerous global variables and made multiple function calls. The system has to make sure that all these global variables should be available for any routine placed in the call stack at all times. The system has to provide a lookup mechanism for both local and global variables for each routine. And with global variables, this lookup mechanism may take some time than local variables.

### Import the sub-modules and functions directly

When importing any module to use its sub-modules, classes, or functions, import them directly instead of importing just the module. When accessing objects using *.*, it triggers dictionary lookup using *\_\_getattribute\_\_*. If we call the object multiple times using *.*, that may increase the program time.

```python
# instead of this
import abc
def_obj = abc.Def()

# do this
from abc import Def
def_obj = Def()
```

### Limit the usage of '.'

Speaking of the lookup time with the module's object references, the same can be applied to the referencing of properties and functions of an object(both custom and in-built). 

```python
%%timeit
ll = []
for i in range(len(rand_nums)):
    ll.append(rand_nums[i])
```
```python:output
6.67 ms ± 84.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
```

Adding the list elements by appending with *list.append()* takes more time than following code because the function is assigned to a variable (functions are first-class citizens in Python) and used inside the loop. This simple practice avoids referencing the functions with **'.'** too often and finally limits the need for dictionary lookup.

```python
%%timeit
ll = []
ll_append = ll.append
for i in range(len(rand_nums)):
    ll_append(rand_nums[i])
```
```python:output
5.45 ms ± 116 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
```

### Avoid writing functions unnecessary

It's good to have code separability by using functions for each independent task. But, as functions in Python are relatively more expensive than C/C++ due to boxing and unboxing dynamic variables and other factors, limit writing functions for unnecessary cases like one-liners. 

### Don't wrap lambdas around functions

Overuse or misuse of lambdas is not a good practice. It's common to wrap functions inside lambdas which do the same thing without wrapping.

Consider the following two functions for sorting a list based on absolute values.

```python
def fun_sort_with_lambda(l):
    return sorted(l, key=lambda x: abs(x))

def fun_sort_without_lambda(l):
    return sorted(l, key=abs)
```

If we look at the CPython bytecode for the above functions with lambda expression passed as a key and with *abs* function object as a key,

```python
>>> from dis import dis
>>> dis(fun_sort_with_lambda)
  2           0 LOAD_GLOBAL              0 (sorted)
              2 LOAD_FAST                0 (l)
              4 LOAD_CONST               1 (<code object <lambda> at 0x7fc51b3a19d0, file "<ipython-input-62-c4147c242c71>", line 2>)
              6 LOAD_CONST               2 ('fun_sort_with_lambda.<locals>.<lambda>')
              8 MAKE_FUNCTION            0
             10 LOAD_CONST               3 (('key',))
             12 CALL_FUNCTION_KW         2
             14 RETURN_VALUE

Disassembly of <code object <lambda> at 0x7fc51b3a19d0, file "<ipython-input-62-c4147c242c71>", line 2>:
  2           0 LOAD_GLOBAL              0 (abs)
              2 LOAD_FAST                0 (x)
              4 CALL_FUNCTION            1
              6 RETURN_VALUE

>>> dis(fun_sort_without_lambda)
  2           0 LOAD_GLOBAL              0 (sorted)
              2 LOAD_FAST                0 (l)
              4 LOAD_GLOBAL              1 (abs)
              6 LOAD_CONST               1 (('key',))
              8 CALL_FUNCTION_KW         2
             10 RETURN_VALUE
```

for the function *fun_sort_with_lambda*, there is an additional function has been generated for lambda. We can avoid this function generation without using lambda as we can see in function *fun_sort_without_lambda*. 

### List comprehension is fast

When operating over lists like data structures, list comprehension is faster than traditional methods like looping, functional programming, etc.

```python
%%timeit
rand_nums = []
for _ in range(1000):
    rand_nums.append(randint(1, 100))
```
```output
603 µs ± 15.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
```

The list comprehension version of the above code snippet runs faster.

```python:output
%%timeit
rand_nums = [randint(1, 100) for _ in range(1000)]
```
```python:output
565 µs ± 11 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
```

---

The optimization practices are not limited to the above approaches. We can check how much time it is taking for each line by using libraries like [CProfile](https://docs.python.org/3/library/profile.html) and making changes to run faster. In the next blog, we discuss how to improve Python computing efficiency using multiprocessing.

---

### References
- [Python Speed](https://wiki.python.org/moin/PythonSpeed)
- [Python Optimization](https://aglowiditsolutions.com/blog/python-optimization/)
- [Making Python Programs Blazingly Fast](https://martinheinz.dev/blog/13)


![Super fast Python: Why Python is Slow?](super-fast-python-why-python-slow/super-fast-python-why-python-slow.jpg)

# Super fast Python (Part-1): Why Python is Slow?

Python is an interpreted, high-level, and dynamically typed programming language. Developers prefer Python because of its easy-to-learn, flexibility, fast development, easy-to-code, readability, and many other development choices. From nowhere in the 2000s to the most used programming language right now, Python has come a long way with help of a strong community. But, there has been always a discussion on Python's choice for computation-intensive tasks like Machine Learning as Python is generally slow compared to widely used languages like C/C++ and Java. Even with the development of computation-efficient libraries and packages like Numpy, still Python is slow for general usage. 

Python's core development team has been working to make Python as fast as C/C++. They set a goal to make each Python release significantly faster than the earlier release. The current [Python 3.11 is up to 10-60 percent faster than Python 3.10](https://docs.python.org/3.11/whatsnew/3.11.html#faster-cpython). Let's hope we will reach a state where Python is at least at the same speed level as Java if not C++. 

One common practice to tackle the speed issues in production is to upgrade the hardware resources or upscale the cloud infrastructure that increases the project budget. As Python's core development team is trying to improve Python, it's up to us to leverage the core libraries and code practices to make code faster at the developer's end and eventually use fewer resources and budget.

This is the first post in the series on Python performance and Optimization. The series points out the utilization of inbuilt libraries, low-level code conversions, and other Python implementations to speed-up Python. The other posts included in this series are 

- (Part-1): Why Python is slow? (this post)
- (Part-2): [Good practices to write fast Python code](https://santhalakshminarayana.github.io/blog/super-fast-python-good-practices)
- (Part-3): [Multi-processing in Python](https://santhalakshminarayana.github.io/blog/super-fast-python-multi-processing)
- (Part-4): [Use Cython to get speed as fast as C](https://santhalakshminarayana.github.io/blog/super-fast-python-cython)
- (Part-5): [Use Numba to speed up Python Functions and Numeric calculations](https://santhalakshminarayana.github.io/blog/super-fast-python-numba)

---

## Why Python is slow?

Before looking into how we can optimize the Python code, we should look at first why Python is slow. 

Some of the reasons for Python's slowness is due to its design of core details like how code executes, type-inference, and memory management. 

### Python implementation by interpretation

Python is a programming language that talks about syntax and rules to write programs. Executing the code written is done by the [programming language implementation](https://en.wikipedia.org/wiki/Programming_language_implementation). This can be of two categories - Compilation and Interpretation. **CPython** is a Python implementation written in **C** that applies an interpretation approach to execute the python code written. CPython is the default runtime and reference implementation of Python and there are other runtimes like PyPy, Cython, Jython, etc., that take different execution approaches for different use cases. CPython is both an interpreter (widely represented) and a compiler as it complies the Python code to Python bytecode and then interprets it for the specific platform using Python Virtual Machine (PVM).

CPython uses Global Interpreter Lock (GIL) on each CPython interpreter process. This means, within a single process, only a single thread processes the Python bytecode. We will see later why this behavior is both good and bad.

As interpretation is usually slow compared to compilation it is understood that Python is slow, but maybe up to 2x-3x times slower than C/C++. But it is not true, Java, which also interprets the code, is still many times faster than Python and this asks questions on what are other factors for slowness in Python. 

Changing the runtime of Python from interpretation to compilation like Cython and embedded C-code, or JIT with PyPy, we can improve the speed of Python many times. We will talk about Cython, and how to use it with Python, in a later article.

### Comes the dynamically typing comes the problem 

One of the beautiful features of Python is dynamic typing and many people like the way it is. But with dynamic typing, there is an additional burden on the interpreter to keep track of the type of variables that makes less scope for optimization. As Java is statically typed, the interpreter can optimize the bytecode generation and can interpret it faster than Python. 

```python
a = 1 # a as int
b = a * 2 # b as int 

a = 'python ' # a as string
# c as a string with operations
# of string and int
c = a * b 
```

In the above snippet, the interpreter has to keep track of the *type* of **a** from top to bottom. If not dynamic typing is supported in Python, we have to declare types for every variable and the operation **c = a * b** is not possible then. So, with flexible support from Python, there is some overhead with the interpreter too.

### That is an object, this too

In Python, everything is an object. Even functions too. If one can remember how objects are referenced in C++ (yes, the pointers), apply the same concept in Python for everything including built-in primitive (not exactly) types which are objects too and are specially taken care of. 

Assigning memory to the objects is done by creating actual memory to hold an object and a reference is given to the variables that point to the real object in the memory. Every time the variable changes its value, instead of changing the value in the memory, a new object is created and the variable is given the new location of the object that holds the new value. As CPython is implemented in C, objects created in Python are called PyObjects which are struct in C that refers to all Python objects.

```python
>>> a = 10
>>> print(id(a))
9801536

>>> a = 12
>>> print(id(a))
9801600
``` 

We can see the address of **a** changes every time we change/assign the value because of how Python manages variables and their values. In C/C++, for variables, the addressing is done by creating an actual location for the variable instead of the value, and if the value changes, it just overwrites the previous value.

Python has to create new objects, keep their references, delete unused objects, and repeat the cycle. This continuous cycle of object creation and deletion whenever variable values changes make runtime slow.

The major overhead that needs to address is Python's way of handling collections like Lists. In *C*, we create the arrays with fixed memory and the address of that array is fixed, and continuous memory allocation happens for all the elements starting from the fixed memory point. But in Python, for each value, an [object is created](http://jakevdp.github.io/blog/2014/05/09/why-python-is-slow/) anywhere in the memory (not continuous), and the list holds only references to these objects. Working on these arbitrary memory locations make things complex and eventually takes more time.

```python
>>> l = [1, 2, 3, 4]
>>> for i in l: 
        print(id(i))

OUTPUT:
-------

9801248
9801280
9801312
9801344
```

In the above snippet, we can observe that the addresses for each element in the list are not the same or continuous.

One of the reasons why Numpy is faster is because it creates fixed memory of the array elements. The below snippet prints out the memory address of each element in the array and they are the same, meaning Numpy create an array with continuous memory address and they are easy to operate.

```python
>>> a = np.array([1, 2, 3, 4])
>>> for i in a: 
        print(i.__array_interface__['data'][0]) 

OUTPUT:
-------

37936064
37936064
37936064
37936064
```


### GIL, good and bad:=:gil-good-and-bad

In **C**, there is no inbuilt support for garbage collection. One has to manually de-allocate/free the memory. Unlike C, Python does garbage collection by using a mechanism called reference count. For every PyObject, Python keeps track of a count of how many references are pointing to the current object. If no reference is pointing to the current object, then, Python frees the memory by deleting the unused object. 

Earlier we talked about GIL that, in the current process, GIL locks the process to have only a single thread to run the Python interpreter. This does not mean we cannot use more threads, but it is hard and takes more time compared to multi-threading in other languages because CPython creates a separate environment for each thread. This limits multi-threading in Python to use only in situations where threads wait for external processes to complete like I/O operations, network requests, etc. 

When multiple threads access the same resource, it is hard to keep track of the reference count of objects. GIL is required in CPython to ensure thread-safe as it allows only a single thread to process. As Python is a general-purpose language, in most cases there is no need to implement multi-threading, and there is no problem with GIL.

But when there is a need to optimize CPU-intensive tasks, one must use multi-processing to utilize all cores and reduce the overall speed. We will talk about different use cases for multi-threading and multi-processing in the later part of this series.

---

Now, we have learned that there are a lot of things that contribute to the slowness of Python compared to other languages. We cannot change the design of Python but we can change how we write the code by using good libraries, code snippets, and changing implementations. 

In the upcoming posts, we discuss how we can use both internal and external libraries like *multiprocessing*, *Cython*, and *Numba* to make Python blazingly faster.


![Whiteboard image enhancement in Python](whiteboard-enhance/whiteboard-image-enhancement.jpg)

# Whiteboard image enhancement using OpenCV

Whiteboard images generally contain less contrast and low brightness as they would be captured in mobile under normal room light conditions. Enhancing whiteboard images makes text readable and gives an image with high contrast and brightness. 

We will apply different image-processing techniques to enhance whiteboard images using OpenCV in Python. From this [whiteboard-cleaner](https://gist.github.com/lelandbatey/8677901) gist that enhances whiteboard images using [ImageMagick](https://imagemagick.org/), we will implement those ImageMagick methods in Python.

In that script, the following ImageMagick functions were used to enhance whiteboard images

```shell
-morphology Convolve DoG:15,100,0 -negate -normalize -blur 0x1 -channel RBG -level 60%,91%,0.1
``` 
Above command applies image enhancing functions in order 
- **-morphology Convolve DoG:15, 100, 0**: Difference of Gaussian (DoG) with kernel_radius=15, sigma1=100, and sigma2=0
- **-negate**: Negative of image
- **-normalize**: Contrast stretch image with black=0.15% and white=0.05%
- **-blur 0x1**: Gaussian blur with sigma=1
- **-level 60%,91%,0.1: Stretch image with black=60% and white=91%, and Gamma correction by gamma=0.1

As I found some difficulty for exactly converting the ImageMagick C code to Python, I have changed the order and parameters that would give close results.

---

We will apply series of image-processing methods and effects to enhance whiteboard images in the following order
- Difference of Gaussian (DoG)
- Negative effect
- Contrast Stretching
- Gaussian blur
- Gamma correction
- Color balance

You can find the full code in my Github repository [whiteboard-image-enhance](https://github.com/santhalakshminarayana/whiteboard-image-enhance)

### Import packages and read image

```python
import cv2
import numpy as np

img = cv2.imread('input.jpg')
```

![Whiteboard image input:=:70:=:Input Whiteboard image](whiteboard-enhance/whiteboard-image.jpg)

### Difference of Gaussian (DoG)

[Difference of Gaussians (DoG)](https://en.wikipedia.org/wiki/Difference_of_Gaussians) is the difference of two Gaussian kernel convoluted images. DoG image is obtained by subtracting two Gaussian blurred images with different kernel radius and variance.

Normally $I_{dog}$ (DoG of image) is calculated by subtracting $I_{g1}$ and $I_{g2}$ which are convoluted images with two different Gaussian kernels. But ImageMagick applies convolution after subtracting and scaling two gaussian kernels.

So, we will first subtract two different Gaussian kernels, scale and normalize the dog-kernel to the zero-summing kernel (sum of all elements ~ 0.0)  and then apply convolution.

```python
def normalize_kernel(kernel, k_width, k_height, scaling_factor = 1.0):
    '''Zero-summing normalize kernel'''
    
    K_EPS = 1.0e-12
    # positive and negative sum of kernel values
    pos_range, neg_range = 0, 0
    for i in range(k_width * k_height):
        if abs(kernel[i]) < K_EPS:
            kernel[i] = 0.0
        if kernel[i] < 0:
            neg_range += kernel[i]
        else:
            pos_range += kernel[i]
    
    # scaling factor for positive and negative range
    pos_scale, neg_scale = pos_range, -neg_range
    if abs(pos_range) >= K_EPS:
        pos_scale = pos_range
    else:
        pos_sacle = 1.0
    if abs(neg_range) >= K_EPS:
        neg_scale = 1.0
    else:
        neg_scale = -neg_range
        
    pos_scale = scaling_factor / pos_scale
    neg_scale = scaling_factor / neg_scale
    
    # scale kernel values for zero-summing kernel
    for i in range(k_width * k_height):
        if (not np.nan == kernel[i]):
            kernel[i] *= pos_scale if kernel[i] >= 0 else neg_scale
            
    return kernel

def dog(img, k_size, sigma_1, sigma_2):
    '''Difference of Gaussian by subtracting kernel 1 and kernel 2'''
    
    k_width = k_height = k_size
    x = y = (k_width - 1) // 2
    kernel = np.zeros(k_width * k_height)
    
    # first gaussian kernal
    if sigma_1 > 0:
        co_1 = 1 / (2 * sigma_1 * sigma_1)
        co_2 = 1 / (2 * np.pi * sigma_1 * sigma_1)
        i = 0
        for v in range(-y, y + 1):
            for u in range(-x, x + 1):
                kernel[i] = np.exp(-(u*u + v*v) * co_1) * co_2
                i += 1
    # unity kernel
    else:
        kernel[x + y * k_width] = 1.0
    
    # subtract second gaussian from kernel
    if sigma_2 > 0:
        co_1 = 1 / (2 * sigma_2 * sigma_2)
        co_2 = 1 / (2 * np.pi * sigma_2 * sigma_2)
        i = 0
        for v in range(-y, y + 1):
            for u in range(-x, x + 1):
                kernel[i] -= np.exp(-(u*u + v*v) * co_1) * co_2
                i += 1
    # unity kernel
    else:
        kernel[x + y * k_width] -= 1.0
    
    # zero-normalize scling kernel with scaling factor 1.0
    norm_kernel = normalize_kernel(kernel, k_width, k_height, scaling_factor = 1.0)
    
    # apply filter with norm_kernel
    return cv2.filter2D(img, -1, norm_kernel.reshape(k_width, k_height))
``` 

Get Difference of Gaussian (DoG) for image by calling *dog()* function with **radius = 15, sigma_1 = 100, sigma_2 = 0**. Here **radius = 15** for both kernels, for first kernel, **sigma = 100**, and for second kernel, **sigma = 0**. Kernel with **sigma = 0** creates unity kernel means convolution with this kernel gives same image. 

```python
dog_img = dog(img, 15, 100, 0)
```

After applying DoG, the resultant image looks like

![Difference of Gaussian (DoG) Image:=:70:=:Difference of Gaussian (DoG) image](whiteboard-enhance/difference-of-gaussian.jpg)

---

### Negative Image

For $I_{dog}$, get a negative image which is just an inversion of colors (255 - image).

```python
def negate(img):
    '''Negative of image'''
    
    return cv2.bitwise_not(img)
```

```python
negative_img = negate(dog_img)
```

The result of the inversion image is
![Negative image:=:70:=:Negative image](whiteboard-enhance/negative.jpg)

Image content is not much visible as we inverted an image whose most of the pixels are black. So, to improve the contrast, we apply contrast-stretch enhancement for the negative image.

---

### Contrast Stretching

Contrast stretching of an image is the same as histogram equalization but we cap some percentage of pixel values to black (0) and white (255).

```python
def get_black_white_indices(hist, tot_count, black_count, white_count):
    '''Blacking and Whiting out indices same as color balance'''

    black_ind = 0
    white_ind = 255
    co = 0
    for i in range(len(hist)):
        co += hist[i]
        if co > black_count:
            black_ind = i
            break
            
    co = 0
    for i in range(len(hist) - 1, -1, -1):
        co += hist[i]
        if co > (tot_count - white_count):
            white_ind = i
            break
    
    return [black_ind, white_ind]

def contrast_stretch(img, black_point, white_point):
    '''Contrast stretch image with black and white cap'''
    
    tot_count = img.shape[0] * img.shape[1]
    black_count = tot_count * black_point / 100
    white_count= tot_count * white_point / 100
    ch_hists = []
    # calculate histogram for each channel
    for ch in cv2.split(img):
        ch_hists.append(cv2.calcHist([ch], [0], None, [256], (0, 256)).flatten().tolist())
    
    # get black and white percentage indices
    black_white_indices = []
    for hist in ch_hists:
        black_white_indices.append(get_black_white_indices(hist, tot_count, black_count, white_count))
        
    stretch_map = np.zeros((3, 256), dtype = 'uint8')
    
    # stretch histogram 
    for curr_ch in range(len(black_white_indices)):
        black_ind, white_ind = black_white_indices[curr_ch]
        for i in range(stretch_map.shape[1]):
            if i < black_ind:
                stretch_map[curr_ch][i] = 0
            else:
                if i > white_ind:
                    stretch_map[curr_ch][i] = 255
                else:
                    if (white_ind - black_ind) > 0:
                        stretch_map[curr_ch][i] = round((i - black_ind) / (white_ind - black_ind)) * 255
                    else:
                        stretch_map[curr_ch][i] = 0
    
    # stretch image
    ch_stretch = []
    for i, ch in enumerate(cv2.split(img)):
        ch_stretch.append(cv2.LUT(ch, stretch_map[i]))
        
    return cv2.merge(ch_stretch)
```

For each image channel, calculate cummulative histogram sum, and then cap pixels based on **black_point = 2** and **white_point = 99.5** percentage.

```python
contrast_stretch_img = contrast_stretch(negative_img, 2, 99.5)
``` 

Negative image after contrast stretching is

![Contrast stretch image:=:70:=:Contrast strech image](whiteboard-enhance/contrast-stretch.jpg)

---

### Gaussin Blur & Gamma Correction

Contrast stretching image contains noise, so blur the image with Gaussian kernel. As Gaussian distribution kernel can be linearly separable, we apply convolution with the same 1D-kernel along the x-axis and y-axis for performance (negligible for small kernels and low-res images).

```python
def fast_gaussian_blur(img, ksize, sigma):
    '''Gussian blur using linear separable property of Gaussian distribution'''
    
    kernel_1d = cv2.getGaussianKernel(ksize, sigma)
    return cv2.sepFilter2D(img, -1, kernel_1d, kernel_1d)

def gamma(img, gamma_value):
    '''Gamma correction of image'''
    
    i_gamma = 1 / gamma_value
    lut = np.array([((i / 255) ** i_gamma) * 255 for i in np.arange(0, 256)], dtype = 'uint8')
    return cv2.LUT(img, lut)
```

Apply Gaussian blur with **kernel_size = 3** and **sigma = 1**.

```python
blur_img = fast_gaussian_blur(contrast_stretch_img, 3, 1)
```

Blurred image after noise suppression is

![Gaussian blur image:=:70:=:Blurred image](whiteboard-enhance/gaussian-blur.jpg)

Now apply Gamma correction to enhance the blurred image with **gamma_value = 1.1**.

```python
gamma_img = gamma(blur_img, 1.1)
```

Blurred image Gamma corrected looks like

![Gamma correction:=:70:=:Gamma corrected image](whiteboard-enhance/gamma-correction.jpg)

---

### Color Balance

Color balance of an image is same as contrast-stretching method above but they are different in implementation. Above contrast-stretching is an implementation based on [ImageMagick-ContrastStretchImage()](https://imagemagick.org/api/MagickCore/enhance_8c.html), and [color balance](https://gist.github.com/DavidYKay/9dad6c4ab0d8d7dbf3dc) is based on [Simplest Color Balance](https://www.ipol.im/pub/art/2011/llmps-scb/article.pdf).

```python
def color_balance(img, low_per, high_per):
    '''Contrast stretch image by histogram equilization with black and white cap'''
    
    tot_pix = img.shape[1] * img.shape[0]
    # no.of pixels to black-out and white-out
    low_count = tot_pix * low_per / 100
    high_count = tot_pix * (100 - high_per) / 100
    
    cs_img = []
    # for each channel, apply contrast-stretch
    for ch in cv2.split(img):
        # cummulative histogram sum of channel
        cum_hist_sum = np.cumsum(cv2.calcHist([ch], [0], None, [256], (0, 256)))

        # find indices for blacking and whiting out pixels
        li, hi = np.searchsorted(cum_hist_sum, (low_count, high_count))
        if (li == hi):
            cs_img.append(ch)
            continue
        # lut with min-max normalization for [0-255] bins
        lut = np.array([0 if i < li 
                        else (255 if i > hi else round((i - li) / (hi - li) * 255)) 
                        for i in np.arange(0, 256)], dtype = 'uint8')
        # constrast-stretch channel
        cs_ch = cv2.LUT(ch, lut)
        cs_img.append(cs_ch)
        
    return cv2.merge(cs_img)
```

Enhance image by passing Gamma corrected image to **color_balance()** with parameters **low_per = 2** and **high_per = 1**.

```python
color_balanced_img = color_balance(gamma_img, 2, 1)
```

The final enhanced whiteboard image is
![Whiteboard image enhance:=:70:=:Whiteboard image enhanced](whiteboard-enhance/whiteboard-image-enhanced.jpg)

You can find out full code at my Github repository file [whiteboard_image_enhance.py](https://github.com/santhalakshminarayana/whiteboard-image-enhance/blob/main/whiteboard_image_enhance.py)

---

## Results

![img_2:=:100](whiteboard-enhance/img_2.jpg)

![img_3:=:100](whiteboard-enhance/img_3.jpg)

![img_4:=:100](whiteboard-enhance/img_4.jpg)


![Color Properties, Models and Spaces](color-theory/color-theory.jpg)

# Color Theory

> Before reading color theory, it is recommended to know about color science that explains how color is defined, represented, and quantified in my previous article about [color science](https://santhalakshminarayana.github.io/blog/color-science)

The Color theory deals with multiple color terminology, models, and schemes that we use in our daily day-to-day life. The color theory describes relationships between different colors and color schemes. It acts like guidelines for better marketing and design of a product. Color theory is a set of rules that combines art and science. The color theory defines the logical structure of color to create color palettes, color schemes, aesthetically pleasing color designs, and color psychology.

From an e-commerce website to an offline store, still photography to a motion picture, interior design to an art gallery, we can see the usage of color theory to attract, please, or astonish the customers.

In fields like image processing, digital photography, graphic design, and game design, color plays a vital role and one must understand basic color theory concepts like color properties, color models, and color spaces before learning other concepts. 


## Color Wheel

A color wheel is a set of colors distributed over a circular disc with some rules that signify the meaning of colors. A color wheel can have numerous colors, and generally, a color wheel is a combination of 12 basic colors from primary, secondary, and teritary colors.

Hue, saturation and lightness are three elements used to define a color.

## Hue

Hue is the dominant color among similar colors in a range of the visible spectrum. Hues are the basic colors to which adding white and black produce multiple colors. Hues are generally classified into major colors we see in nature. They are Red, Yellow, Green, Blue, Violet, and Orange. 

**Primary colors** - Primary hues are colors that cannot be produced by mixing other hues. In painting world **red**, **yellow**, and **blue** are primary hues. In digital world **red**, **green**, and **blue** are primary colors.

**Secondary colors** - Mixing two primary hues creates secondary hues:
- **Green** - (Blue + Yellow)
- **Orange** - (Yellow + Red)
- **Purple** - (Red + Blue)

**Tertiary colors** - In the color wheel, tertiary hues are mixtures of adjacent primary and secondary hues. Thus, six territory hues are produced.
- **Chartreuse** - (Yellow + Green)
- **Teal** - (Green + Blue)
- **Violet** - (Blue + Purple)
- **Magenta** - (Purple + Red)
- **Vermillion** - (Red + Orange)
- **Amber** - (Orange + Yellow)

![Color Wheel - Primary, Secondary, and Tertiary Hues:=:50:=:Color wheel](color-theory/color-wheel.jpg)

## Saturation/Chroma

Saturation/Chroma is the purity of a hue. It defines the color intensity. A pure color is 100% saturated or fully saturated. Mixing gray (white + black) color to a pure hue decreases the intensity and finally reaches a 0% saturation level. A color can be desaturated by mixing with its complementary color. A hue at 0% saturation is a gray color. Saturation defines a range of color intensities from pure color (100%) to gray (0%). 

![Saturation:=:70](color-theory/saturation.jpg)

### Tint

When white color is mixed with a color, the resultant is the Tint of that color. Adding white color makes color lighter than the original. The amount of white added defines the tint of color.

![Tints:=:70](color-theory/tints.jpg)

### Shade

A shade of a color is the blackness of the color after being mixed with black color. The shade of a color is darker than the original color.

![Shades:=:70](color-theory/shades.jpg)

### Tone 

A Tone is mixing both white and black colors to a color. The resultant color is grayer than the original.

![Tones:=:70](color-theory/tones.jpg)

Although Chroma and Saturation speak the same, Chroma is an absolute term that defines color value a range from gray (0%) to pure hue (100%), and saturation defines the brilliance of a color relative to gray.

## Lightness/Luminance

Lightness/Luminance define how light (100%) or dark (0%) a color is. All pure hues have a lightness value of 50%. Lightness is measured relative to the brightness of white color. In color perception, lightness tells how much light is reflected from a surface. 

![Lightness:=:70](color-theory/lightness.jpg)

---

## Color properties and color terms

### Chromatic & Achromatic colors

Achromatic colors are colors that have lightness but no hue or saturation. They are black, white and gray. Mixing complementary colors produces achromatic colors.

Chromatic colors are the ones that have any amount of hue and saturation. The presence of saturation defines chromatic colors. All other colors except achromatic colors are chromatic colors.

### Chromatic value/Chromatic intensity/Chromaticity

Chromatic value/intensity is the measure of how light or dark a color is. Depending on how light or dark a hue is, Tints and Shades of color are created. Chromaticity is similar to saturation as they both speak about how much color is chromatic (a hue with some saturation). Chromaticity is also defined as color independent of lightness (only hue and saturation).


### Luminance

Luminance is the intensity of light emitted, reflected, or passes through per unit area in a given direction.

### Brightness

Brightness is a relative term that describes how much light is shining from something. It is the average lightness of an object. It is not a color term but a perception of our eyes created by color's lightness. It is not a quantitative term but can be scaled (%). It is quite similar to luminance in human color vision. [Luminance](http://www.workwithcolor.com/color-luminance-2233.htm) is the brightness of a color perceived by our eyes. A high-intensity color looks more bright. The green color looks more bright than the blue color. A hue with more saturation is brighter than a desaturated hue. A color with the same luminance appears in different brightness depending on the surrounding colors.

### Contrast

Contrast is the difference between the luminance of two or more colors placed adjacent or overlapped. High contrast colors appear more bright than low contrast. Generally, contrast is used in the context of foreground and background colors. When equal luminance colors are placed on foreground and background, they look similar and there is no difference observed. Depending on the position of colors in the color wheel, contrast levels vary. Colors next to each other have low contrast and opposite colors have high contrast.

### Color contrast ratio

Color contrast ratio is the quantitative term to define the contrast between two colors. They are measured based on the relative luminance of different colors. The [Web Content Accessibility Guidelines (WCGA)](https://www.w3.org/WAI/WCAG22/Understanding/contrast-minimum.html) published guidelines for calculating contrast ratio of color and these are considered as standard rules in web design.

### Luma

Luma describes the presence of achromatic signal (white, gray, and black) intensity in a color. Luma describes how bright a color is. The difference between luminance and luma is, luminance is derived from linear RGB values, whereas luma is derived from non-linear (gamma-corrected) RGB values. Luminance is based on standard [relative luminance ratios of RGB](https://santhalakshminarayana.github.io/blog/color-science#relative-luminance), and luma is based on relative gamma-corrected RGB ratios.

### Metamerism (Color metameres)

Metamerism occurs when two colors look the same under one lighting condition but, different when light source/conditions change. This usually happens when we buy clothes in a store like a shirt appears light blue in the store but when we brought it in sunlight the shirt looks pale blue. 

### Color Temperature

The color temperature of the light source is the temperature of black-body radiation when black-body is heated at a constant temperature. As the black body observes all spectrum wavelengths, when heated up at a constant temperature, it emits certain wavelengths which are perceived as color. Color temperature is measured in Kelvin (&degK). At different color temperatures, the body is seen as different colors. 

![Color Temperature:=:60:=:Colors at different temperatures emitted by a black-body](color-theory/color-temperature.jpg)

The lower the temperature, the lower the color appears warm, and the higher the temperature, the cooler the color appears.

### Color Rendering Index (CRI)

Color Rendering Index (CRI) is a measurement of how light affects the appearance of a color. It defines how accurately colors can be distinguished under a light source. Light sources with different color temperatures illuminate objects as different colors and sometimes we cannot distinguish object colors in extremely warm or cool light sources. CRI indicates how well a light source can reproduce colors with the same temperature for natural light such as the sun. It is a quantitative term measured on a scale of 0-100. A score of 90 or more is considered excellent and less than 80 is considered poor.

Color temperature and Color Rendering Index (CRI) are two vital terms that affect the selection of light sources and conditions. Art galleries, museums, photographers, and product displays stores choose light sources based on these terms. 

---

## Color Harmony or Color schemes

Color harmony is the theory of selecting/combining colors to create a visually pleasing color scheme. Color harmony creates aesthetically rich color combinations that maintain harmony and engages users with the product. The color wheel is the basic logical structure for creating color schemes.

### Monochromatic

The monochromatic color scheme contains single color with variations of tints, tones, and shades.

![Monochromatic Harmony:=:25:=:Monochromatic Harmony](color-theory/monochromatic.jpg)

### Analogous

Analogous color schemes are created by grouping the main color with adjacent colors to it in the color wheel.

![Analogous Harmony:=:25:=:Analogous Harmony](color-theory/analogous.jpg)

### Complementary

Complementary colors are two colors selected in a way that colors are opposite to each other in the color wheel.

![Complementary Harmony:=:25:=:Complementary Harmony](color-theory/complementary.jpg)

### Split Complementary

In the split complementary scheme, the main color is selected and the other two colors are adjacent colors to the complementary color of the main color.

![Split Complementary Harmony:=:25:=:Split Complementary Harmony](color-theory/split-complementary.jpg)

### Triad

Triad is a group of three colors that are equally distant from each other in the color wheel and form a triangle.

![Triad Harmony:=:25:=:Triad Harmony](color-theory/triad.jpg)

### Square

The square color scheme has four colors that are the same distance to each other forming a square or rhombus in the color. 

![Square Harmony:=:25:=:Square Harmony](color-theory/square.jpg)

Other harmonies include diad, diad complementary, rectangle, and polygon schemes.


## Color models

Color models are the organization of colors and their mixing (additive or subtractive) based on human color perception. A color model defines color and its properties in a mathematical way.

### RGB

In the RGB color model, Red, Green, and Blue are three primary colors that are mixed to produce all other colors by [additive color mixing](https://santhalakshminarayana.github.io/blog/color-science#additve-subtractive-color-mixing). Display systems and digital photography stores color values based on RGB model.

Based on the depth (no. of bits allocated to store) of a color, RGB values are scaled in a certain range. For an 8-bit system, the range is 0-255, and when all colors are mixed at minimum level then it designates black color and at maximum level designates white color. If RGB colors are mapped to 3-dimensional space with unit vectors X, Y, and Z-axis as R, G, and B colors, then the vector space forms a cube with all possible colors in it formed by the combination of RGB at different values.

![RGB color model:=:35](color-theory/rgb-cube.jpeg)

With the RGB model, one cannot easily identify important color terms like lightness, contrast, brightness, tones, tint, and shades.

### HSV (HSB)

Hue, Saturation, and Value (or Brightness) (HSV or HSB) is a cylindrical color model that is based on hue, saturation, and brightness described above. 

Hue is arranged with all hues taking up a certain range in a 0-360 degrees circle. Starting red at 0&deg, green at 120&deg, blue at 240&deg, and again wrapping back to red at 360&deg.

Saturation has a range from 0%-100%, starting from 0% (gray) at the center of the base to 100% (pure color/hue) at the circumference. Saturation controls tints and shades

Value/Brightness also has a range of 0%-100%, with 0 (black) at the bottom of the base to 100% (no black) at the top.

![HSV:=:35](color-theory/hsv.jpeg)

If brightness is at 0% then the color is black irrespective of hue and saturation. But for white color, the brightness value is 100% and saturation is kept at 0% while hue can be anything.

### HSL

Hue, Saturation, and Lightness (HSL) is a bi-cone (double) model while H and S are the same as HSB. 

![HSL Cone:=:35](color-theory/hsl-cone.jpeg)

Lightness has a range from 0%-100%, starting from 0% (black) at bottom of the cone to 100%(white) at top of the cone. All pure hues have a lightness value of 50%.

![HSL:=:35](color-theory/hsl.jpeg)

Color is white if lightness is 100% irrespective of hue and saturation. Same for black color with lightness at 0% irrespective of hue and saturation values. 

## Color spaces

A color space is based on the color model that maps colors to a set of colors like sRGB and Adobe RGB and is reproduced by display systems.

### L\*a\*b* (CIELAB) color space 

L\*a\*b* color space is a device-independent color space in which color is expressed in three components:
- **L***- Lightness
- **a***- Green-Red
- **b***- Blue-Yellow

It was defined by CIE and also called CIELAB (Lab without * is color space defined by Hunter Lab).

In L\*a\*b\* color space, L indicates lightness, and a\* and b\* indicates chromaticity coordinates. Unlike RGB and HSV, L\*a\*b* represents colors in the sphere where L\*, a\*, and b\* are orthogonal axis to each other. CIELAB is a perceptually uniform color space i.e. a standard color system that reflects the color representation close to human vision color perception. 

![Lab color space:=:35:=:Lab sphere$$$https://sensing.konicaminolta.asia/what-is-cie-1976-lab-color-space/](color-theory/3d-lab.jpg)

L\*a\*b* is based on [opponent process theory](https://santhalakshminarayana.github.io/blog/color-science#opponent-process-theory) where red-green and blue-yellow form opponent pairs. 

L\*-axis is a positive axis with coordinates black at 0 to white at 100. 
a\*-axis is present along horizontal from left to right in the range of -128 to 127, +a\* is the red axis, and -a\* is the green axis.
b\*-axis runs along top to down in the range of -128 to 127, +b\* axis is the yellow axis, and -b\* axis is the blue axis.

![Lab color space:=:40:=:Lab color space$$$https://www.xrite.com/blog/lab-color-space](color-theory/lab-color-space.jpg)

At the center, colors are achromatic colors, and saturation increases moving towards the circumference.

### YUV and YCbCr

[YUV](https://en.wikipedia.org/wiki/YUV) and [YCbCr](https://en.wikipedia.org/wiki/YCbCr) are a family of color spaces used to encode color data. In both models, the Y component refers to luminance but practically Y'UV and Y'CbCr are used. Y' component is luma which is scaled in the range of 0-100. The other two components are chrominance components derived by the difference of Y' on blue and red values.

YUV, historically used for transmitting video data for analog devices. In the olden days, as there was a need to support both black and white display and color displays, a conversion has be to applied for RGB data to transmit signals without disturbing black-and-white channels. As the black-and-white image is the same as the luminance channel, YUV was developed to pass chrominance channels along with luminance. Later, YCbCr was developed to store, compress and encode digital photographs. YUV is for analog TVs and YCbCr is for digital TVs. 

Using [chroma subsampling](https://www.matrox.com/en/video/media/guides-articles/introduction-color-spaces-video), image data can be transmitted and stored at a low-bit rate for chrominance channels UV and CbCr, because human eyes are tolerant for perceptually equal colors. This made transmitting or storing image data in [Y'UV or Y'CbCr](https://docs.microsoft.com/en-us/windows/win32/medfound/about-yuv-video) requires less space: a high bit-rate luma channel and low bit-rate chrominance channels.

Conversion from **RGB** to **Y'UV** is computed as follows:
$$
Y' = W_{R}R + W_{G}G + W_{B}B
$$
$$
Y' = 0.299R + 0.587G + 0.114B
$$
$$
U = U_{max}\frac{B - Y'}{1 - W_B}
$$
$$
V = V_{max}\frac{R - Y'}{1 - W_R}
$$

Where $W_R$, $W_G$, and $W_B$ are relative ratios of non-linear gamma-corrected RGB values. $U_{max}$ and $V_{max}$ are maximum values of the numerical range selected. For PAL and NTSC, these values were $U_{max} = 0.436$ and $V_{max} = 0.615$.

Similarly, **RGB** to **Y'CbCr** is computed as
$$
Y' = K_{R}R + K_{G}G + K_{B}B
$$
$$
Y' = 0.299R + 0.587G + 0.114B
$$
$$
C_b = \frac{1}{2}\frac{B - Y'}{1 - K_B}
$$
$$
C_r = \frac{1}{2}\frac{R - Y'}{1 - K_R}
$$

Where $K_R$, $K_G$, and $K_B$ are relative ratios of non-linear gamma-corrected RGB values and satisfies $K_R + K_G + K_R = 1$.

---

## Device-dependent and device-independent color spaces

**Device-dependent** color models mean they depend on the display system subset of colors or reference color space. This model explains the physical device output rather than human color vision. Coordinates used to display colors in device-dependent models change when the display system changes its light source. RGB, CMY, HSV, and HSL are device-dependent color models.

**Device-independent** color spaces are universal references and coordinates used to specify colors will produce the same color in all conditions because they define color output based on human color vision. These are used to convert device-dependent color spaces across display devices. CIELAB, CIELUV, YCbCr, and YUV are device-independent color spaces.