Simplify the version bump process. ### Description

Simplify the version bump process.

- Use substitution hook in docs to update CDN URL.
  - In `dev` or `prerelease` versions, it will point to `dev` url
  - In `stable` versions, it will point it `vX.Y.Z` url
- Introduce a version bump script `tools/bump_version.py`.


~- Use `bump2version` to bump Python versions (continuing from #2460)~
~- Use `npm version` to bump JS versions~
~- Update `run_docker` script to use latest deployed pre-built docker image when `--pre-built` option is specified.~

### Checklists

- [x] Add new / update outdated documentation

# setup: import numpy ; a = numpy.array([ [i/10., i/10., i/20.] for i in range(44440)], dtype=numpy.double)
# run: hyantes(0, 0, 90, 90, 1, 100, 80, 80, a)

# pythran export hyantes(float, float, float, float, float, float, int,
# int, float[][])
import numpy as np


def hyantes(xmin, ymin, xmax, ymax, step, range_, range_x, range_y, t):
    X, Y = t.shape
    pt = np.zeros((X, Y))
    for i in range(X):
        for j in range(Y):
            for k in t:
                tmp = 6368.0 * np.arccos(
                    np.cos(xmin + step * i)
                    * np.cos(k[0])
                    * np.cos((ymin + step * j) - k[1])
                    + np.sin(xmin + step * i) * np.sin(k[0])
                )
                if tmp < range_:
                    pt[i, j] += k[2] / (1 + tmp)
    return pt

# setup: import numpy ; a = numpy.array([ [i/10., i/10., i/20.] for i in range(44440)], dtype=numpy.double)
# run: hyantes(0, 0, 90, 90, 1, 100, 80, 80, a)

# pythran export hyantes(float, float, float, float, float, float, int,
# int, float[][])
import numpy as np


def hyantes(xmin, ymin, xmax, ymax, step, range_, range_x, range_y, t):
    X, Y = t.shape
    pt = np.zeros((X, Y))
    for i in range(X):
        for j in range(Y):
            for k in t:
                tmp = 6368.0 * np.arccos(
                    np.cos(xmin + step * i)
                    * np.cos(k[0])
                    * np.cos((ymin + step * j) - k[1])
                    + np.sin(xmin + step * i) * np.sin(k[0])
                )
                if tmp < range_:
                    pt[i, j] += k[2] / (1 + tmp)
    return pt
