feat(scripts): polynomial and rational approximations

.gitignore

@@ -230,3 +230,6 @@ blocks.db
 # Ignore e2e test artifacts (which clould leak information if commited)
 .ash_history
 .bash_history
+
+# Python
+**/venv/*

scripts/approximations/README.md

@@ -0,0 +1,48 @@
+# Math Operations Approximations
+
+## Context
+
+This is a script to approximate a mathematical operation using polynomial
+and rational approximations.
+
+This script does the following:
+
+1. Computes polynomial and rational approximations of a given function (e^x by default), 
+returning the coefficients.
+
+2. Computes (x,y) coordinates for every kind of approximation given the same x coordinates.
+Plots the results for rough comparison.
+
+
+The following are the resources used to write the script:
+- <https://xn--2-umb.com/22/approximation>
+- <https://sites.tufts.edu/atasissa/files/2019/09/remez.pdf>
+
+## Configuration
+
+There are several parameters that can be changed on the needs basis at the
+top of `main.py`.
+
+Some of the parameters include the function we are approximating, the [x_start, x_end] range of
+the approximation, and the number of terms to be used. For the full parameter list, see `main.py`.
+
+## Usage
+
+Assuming that you are in the root of the repository:
+
+```bash
+# Create a virtual environment.
+python3 -m venv scripts/approximations/venv
+
+# Start the environment
+source scripts/approximations/venv/bin/activate
+
+# Install dependencies in the virtual environment.
+pip install -r scripts/approximations/requirements.txt
+
+# Run the script.
+python3 scripts/approximations/main.py
+
+# Run tests
+python3 scripts/approximations/rational_test.py
+```

scripts/approximations/approximations.py

@@ -0,0 +1,80 @@
+import numpy as np
+from typing import Tuple
+
+import rational
+import polynomial
+import chebyshev
+
+def equispaced_poly_approx(fn, x_start: int, x_end: int, num_terms: int) -> list[np.ndarray]:
+    """ Returns the coefficients for an equispaced polynomial between x_start and x_end with num_terms terms.
+
+    The return value is a list of num_terms polynomial coefficients needed to get the returned y coordinates from returned x coordinates.
+    """
+    # Compute equispaced coordinates.
+    equispaced_nodes_x = np.linspace(x_start, x_end, num_terms)
+    y_nodes = fn(equispaced_nodes_x)
+
+    # Construct a system of linear equations.
+    vandermonde_matrix = polynomial.construct_vandermonde_matrix(equispaced_nodes_x.tolist())
+
+    # Solve the matrix to get the coefficients used in the final approximation polynomial.
+    coef = np.linalg.solve(np.array(vandermonde_matrix), y_nodes)
+
+    return coef
+
+def chebyshev_poly_approx(fn, x_start: int, x_end: int, num_terms: int) -> np.ndarray:
+    """ Returns the coefficients for a polynomial constructed from Chebyshev nodes between x_start and x_end with num_terms terms.
+
+    Equation for Chebyshev nodes:
+    x_i = (x_start + x_end)/2 + (x_end - x_start)/2 * cos((2*i + 1)/(2*num_terms) * pi)
+
+    The return value is a list of num_terms polynomial coefficients needed to get the returned y coordinates from returned x coordinates.
+    """
+    # Compute Chebyshev coordinates.
+    x_estimated, y_estimated = chebyshev.get_nodes(fn, x_start, x_end, num_terms)
+
+    # Construct a system of linear equations.
+    vandermonde_matrix = polynomial.construct_vandermonde_matrix(x_estimated)
+
+    # Solve the matrix to get the coefficients used in the final approximation polynomial.
+    coef = np.linalg.solve(np.array(vandermonde_matrix), y_estimated)
+
+    return coef
+
+def chebyshev_rational_approx(fn, x_start: int, x_end: int, num_terms_numerator: int) -> Tuple[np.array, np.array]:
+    """ Returns a rational approximation between x_start and x_end with num_terms terms
+    using Chebyshev nodes.
+
+    Equation for Chebyshev nodes:
+    x_i = (x_start + x_end)/2 + (x_end - x_start)/2 * cos((2*i + 1)/(2*num_terms) * pi)
+
+    We approximate h(x) = p(x) / q(x)
+
+    Assume num_terms_numerator = 3.
+    Then, we have
+
+    h(x) = (p_0 + p_1 x + p_2 x^2) / (1 + q_1 x + q_2 x^2) 
+
+    The return value is a list with a 2 items where:
+    - item 1: num_terms equispaced x coordinates between x_start and x_end
+    - item 2: num_terms y coordinates for the equispaced x coordinates
+    """
+    # Compute Chebyshev coordinates.
+    x_chebyshev, y_chebyshev = chebyshev.get_nodes(fn, x_start, x_end, num_terms_numerator * 2 - 1)
+
+    # Construct a system of linear equations.
+    vandermonde_matrix = rational.construct_rational_eval_matrix(x_chebyshev, y_chebyshev)
+
+    # Solve the matrix to get the coefficients used in the final approximation polynomial.
+    coef = np.linalg.solve(np.array(vandermonde_matrix), y_chebyshev)
+
+    # first num_terms_numerator values are the numerator coefficients
+    # next num_terms_numerator - 1 values are the denominator coefficients
+    coef_numerator = coef[:num_terms_numerator]
+    coef_denominator = coef[num_terms_numerator:]
+
+    # h(x) = (p_0 + p_1 x + p_2 x^2) / (1 + q_1 x + q_2 x^2)
+    # Therefore, we insert 1 here.
+    coef_denominator = np.insert(coef_denominator, 0, 1, axis=0)
+
+    return [coef_numerator, coef_denominator]

scripts/approximations/chebyshev.py

@@ -0,0 +1,25 @@
+
+import math
+from typing import Tuple
+
+def get_nodes(fn, x_start: int, x_end: int, num_terms: int ) -> Tuple[list, list]:
+    """ Returns Chebyshev nodes between x_start and x_end with num_terms terms
+    and the given function fn.
+
+    Equation for Chebyshev nodes:
+    x_i = (x_start + x_end)/2 + (x_end - x_start)/2 * cos((2*i + 1)/(2*num_terms) * pi)
+
+    The first returned value is a list of x coordinates for the Chebyshev nodes.
+    The second returned value is a list of y coordinates for the Chebyshev nodes.
+    """
+    x_estimated = []
+    y_estimated = []
+
+    for i in range (num_terms):
+        x_i = (x_start + x_end) / 2 + (x_end - x_start) / 2 * math.cos((2*i + 1) * math.pi / (2 * num_terms))
+        y_i = fn(x_i)
+
+        x_estimated.append(x_i)
+        y_estimated.append(y_i)
+
+    return x_estimated, y_estimated

scripts/approximations/main.py

@@ -0,0 +1,122 @@
+import math
+import numpy as np
+import matplotlib.pyplot as plt
+
+import polynomial
+import rational
+import approximations
+
+# This script does the following:
+# - Computes polynomial and rational approximations of a given function (e^x by default).
+# - Computes (x,y) coordinates for every approximation given the same x coordinates.
+# - Plots the results for rough comparison.
+# The following are the resources used to write the script:
+# https://xn--2-umb.com/22/approximation/
+# https://sites.tufts.edu/atasissa/files/2019/09/remez.pdf
+def main():
+
+    ##############################
+    # 1. Configuration Parameters
+
+    # start of the interval to calculate the approximation on
+    x_start = 0
+    # end of the interval to calculate the approximation on
+    x_end = 1
+    # number of terms in the polynomial approximation / numerator of the rational approximation.
+    num_terms_approximation = 3
+
+    # number of (x,y) coordinates used to plot the resulting approximation.
+    num_points_plot = 10
+    # number of (x,y) coordinates used to plot the
+    # actual function that is evenly spaced on the X-axis.
+    num_points_plot_accurate = 50000
+    # function to approximate
+    approximated_fn = lambda x: math.e**x
+
+    # flag controlling whether to plot each approximation.
+    # Plots if true.
+    shouldPlotApproximations = False
+
+    # flag controlling whether to compute max error for each approximation
+    # given the equally spaced x coordinates.
+    # Computes if true.
+    shouldComputeErrorDelta = True
+
+    #####################
+    # 2. Approximations 
+
+    # 2.1. Equispaced Polynomial Approximation
+    coefficients_equispaced_poly = approximations.equispaced_poly_approx(approximated_fn, x_start, x_end, num_terms_approximation)
+
+    # 2.2. Chebyshev Polynomial Approximation
+    coefficients_chebyshev_poly = approximations.chebyshev_poly_approx(approximated_fn, x_start, x_end, num_terms_approximation)
+
+    # 2.3. Chebyshev Rational Approximation
+    numerator_coefficients_chebyshev_rational, denominator_coefficients_chebyshev_rational = approximations.chebyshev_rational_approx(approximated_fn, x_start, x_end, num_terms_approximation)
+
+    # 2.4. Actual With Large Number of Coordinates (evenly spaced on the X-axis)
+    x_accurate = np.linspace(x_start, x_end, num_points_plot_accurate)
+
+    #######################################
+    # 3. Compute (x,y) Coordinates To Plot
+
+    # Equispaced x coordinates to be used for plotting every approximation.
+    plot_nodes_x = np.linspace(x_start, x_end, num_points_plot)
+
+    # 3.1 Equispaced Polynomial Approximation
+    plot_nodes_y_eqispaced_poly = polynomial.evaluate(plot_nodes_x, coefficients_equispaced_poly)
+
+    # 3.2 Chebyshev Polynomial Approximation
+    plot_nodes_y_chebyshev_poly = polynomial.evaluate(plot_nodes_x, coefficients_chebyshev_poly)
+
+    # 3.3 Chebyshev Rational Approximation
+    y_chebyshev_rational = rational.evaluate(plot_nodes_x, numerator_coefficients_chebyshev_rational.tolist(), denominator_coefficients_chebyshev_rational.tolist())
+
+    # 3.4 Actual With Large Number of Coordinate (evenly spaced on the X-axis)
+    y_accurate = approximated_fn(x_accurate)
+
+    #############################
+    # 4. Plot Every Approximation
+
+    if shouldPlotApproximations:
+        # 4.1 Equispaced Polynomial Approximation
+        plt.plot(plot_nodes_x, plot_nodes_y_eqispaced_poly, label="Equispaced Poly")
+
+        # 4.2 Chebyshev Polynomial Approximation
+        plt.plot(plot_nodes_x,plot_nodes_y_chebyshev_poly, label="Chebyshev Poly")
+
+        # 4.3 Chebyshev Rational Approximation
+        plt.plot(plot_nodes_x,y_chebyshev_rational, label="Chebyshev Rational")
+
+        # 4.4 Actual With Large Number of Coordinates (evenly spaced on the X-axis)
+        plt.plot(x_accurate,y_accurate, label=F"Actual - {num_points_plot_accurate} evenly spaced points")
+
+        plt.legend(loc="upper left")
+        plt.grid(True)
+        plt.title(f"Appproximation of e^x on [{x_start}, {x_end}] with {num_terms_approximation} terms")
+        plt.show()
+
+    #############################
+    # 5. Compute Errors
+
+    if shouldComputeErrorDelta:
+        print(F"\n\nMax Error on [{x_start}, {x_end}]")
+        print(F"{num_points_plot} coordinates equally spaced on the X axis")
+        print(F"{num_terms_approximation} polynomial terms and ({num_terms_approximation}, {num_terms_approximation}) rational terms used for approximations.\n\n")
+
+        plot_nodes_y_actual = approximated_fn(plot_nodes_x)
+
+        # 5.1 Equispaced Polynomial Approximation
+        delta_eqispaced_poly = np.abs(plot_nodes_y_eqispaced_poly - plot_nodes_y_actual)
+        print(F"Equispaced Poly: {np.amax(delta_eqispaced_poly)}")
+
+        # 5.2 Chebyshev Polynomial Approximation
+        delta_chebyshev_poly = np.abs(plot_nodes_y_chebyshev_poly - plot_nodes_y_actual)
+        print(F"Chebyshev Poly: {np.amax(delta_chebyshev_poly)}")
+
+        # 5.3 Chebyshev Rational Approximation
+        delta_chebyshev_rational = np.abs(y_chebyshev_rational - plot_nodes_y_actual)
+        print(F"Chebyshev Rational: {np.amax(delta_chebyshev_rational)}")
+
+if __name__ == "__main__":
+    main()

scripts/approximations/polynomial.py

@@ -0,0 +1,41 @@
+
+
+import numpy as np
+
+def construct_vandermonde_matrix(x_list: list) -> list[list]:
+    """ Constructs a Vandermonde matrix for a polynomial approximation.
+    from the list of x values given.
+
+    len(x_list) * len(x_list)
+    x_list is the list of all x values to construct the matrix from.
+
+    V = [1 x_0 x_0^2 ... x_0^{n-1}]
+    [1 x_1 x_2^1 ... x_1^{n-1}]
+    ...
+    [1 x_n x_n^2 ... x_n^{n-1}]
+
+    Vandermonde matrix is a matrix with the terms of a geometric progression in each row.
+    We use Vandermonde matrix to convert the system of linear equations into a linear algebra problem
+    that we can solve.
+    """
+    num_terms = len(x_list)
+
+    matrix = []
+
+    for i in range(num_terms):
+        row = []
+        for j in range(num_terms):
+            row.append(x_list[i]**j)
+        matrix.append(row)
+
+    return matrix
+
+def evaluate(x: np.ndarray, coeffs: np.ndarray) -> np.ndarray:
+    """ Evaluates the polynomial. Given a list of x coordinates and a list of coefficients, returns a list of
+    y coordinates, one for each x coordinate. The coefficients must be in ascending order.
+    """
+    o = len(coeffs)
+    y = 0
+    for i in range(o):
+        y += coeffs[i]*x**i
+    return y

scripts/approximations/rational.py

@@ -0,0 +1,47 @@
+import numpy as np
+
+def construct_rational_eval_matrix(x_list: list, y_list: list) -> list[list]:
+    """ Constructs a matrix to use for computing coefficients for a rational approximation
+        from the list of x and y values given.
+        len(x_list) * len(x_list)
+        x_list is the list of all x values to construct the matrix from.
+        V = [1 x_0 x_0^2 ... x_0^{n-1} - y_0*x_0 - y_0*x_0^2 ... - y_0*x_0^{n-1}]
+            [1 x_1 x_2^1 ... x_1^{n-1} - y_1*x_1 - y_1*x_1^2 ... - y_1*x_1^{n-1}]
+            ...
+            [1 x_n x_n^2 ... x_n^{n-1} - y_n*x_n - y_n*x_n^2 ... - y_n*x_n^{n-1}]
+    """
+    num_terms = (len(x_list) + 1) // 2
+
+    matrix = []
+
+    for i in range(num_terms * 2 - 1):
+        row = []
+        for j in range(num_terms):
+            row.append(x_list[i]**j)
+
+        for j in range(num_terms):
+            # denominator terms
+            if j > 0:
+                row.append(-1 * x_list[i]**j * y_list[i])
+
+        matrix.append(row)
+
+    return matrix
+
+def evaluate(x: np.ndarray, coefficients_numerator: list, coefficients_denominator: list) -> np.ndarray:
+    """ Evaluates the rational function. Assume rational h(x) = p(x) / q(x)
+    Given a list of x coordinates, a list of coefficients of p(x) - coefficients_numerator, and a list of
+    coefficients of q(x) - coefficients_denominator, returns a list of y coordinates, one for each x coordinate.
+    """
+    num_numerator_coeffs = len(coefficients_numerator)
+    num_denominator_coeffs = len(coefficients_denominator)
+
+    p = 0
+    for i in range(num_numerator_coeffs):
+        p += coefficients_numerator[i]*x**i
+
+    q = 0
+    for i in range(num_denominator_coeffs):
+        q += coefficients_denominator[i]*x**i
+
+    return p / q