From 5f55d1beaf379ccd0d5d31274fdd054887c652d6 Mon Sep 17 00:00:00 2001 From: Roman Date: Sat, 10 Dec 2022 15:07:29 -0500 Subject: [PATCH] feat(scripts): polynomial and rational approximations (#3620) * feat(scripts): polynomial and rational approximations * typos * md lint * osmomath things * rename solve -> evaluate * terminology corrections for actual function * correct terminology for rational matrix * Revert "osmomath things" This reverts commit 90f137ed903859460f332c4065ab4c699f5ca8e6. * compute and print max error * increase num_points_plot * polynomial_num_terms = numerator_terms + denominator_terms * add ability to approximate errors * clean up plot errors * clean up * clean up * clean up * clean up * refactores equispaced with higher precision * fix chebyshev poly * refactor chebyshev rational * clean up * updates * updates * updates * plot errors on range * isolate exponent approximation choice with expected precision * add comments in main.py * README updates * update readme and requirements * requirements.txt and gitignore updates * clean ups and docs * readme * update test test_construct_matrix_3_3 * feat(osmomath): BigDec power function with integer exponent, overflow tests at max spot price (#3676) * feat(osmomath): BigDec power function with integer exponent, overflow tests at max spot price * euler's number * nolint * more tests * changelog * feat(osmomath): mutative `PowerIntegerMut` function on `osmomath.BigDec` (#3680) --- .gitignore | 4 + CHANGELOG.md | 3 +- osmomath/decimal.go | 28 ++++ osmomath/decimal_test.go | 198 ++++++++++++++++++++++- osmomath/math.go | 4 + scripts/approximations/README.md | 65 ++++++++ scripts/approximations/approximations.py | 191 ++++++++++++++++++++++ scripts/approximations/chebyshev.py | 29 ++++ scripts/approximations/main.py | 194 ++++++++++++++++++++++ scripts/approximations/polynomial.py | 41 +++++ scripts/approximations/rational.py | 44 +++++ scripts/approximations/rational_test.py | 36 +++++ scripts/approximations/requirements.txt | 16 ++ 13 files changed, 844 insertions(+), 9 deletions(-) create mode 100644 scripts/approximations/README.md create mode 100644 scripts/approximations/approximations.py create mode 100644 scripts/approximations/chebyshev.py create mode 100644 scripts/approximations/main.py create mode 100644 scripts/approximations/polynomial.py create mode 100644 scripts/approximations/rational.py create mode 100644 scripts/approximations/rational_test.py create mode 100644 scripts/approximations/requirements.txt diff --git a/.gitignore b/.gitignore index 705d20fc0c9..f9b128d098d 100644 --- a/.gitignore +++ b/.gitignore @@ -230,3 +230,7 @@ blocks.db # Ignore e2e test artifacts (which clould leak information if commited) .ash_history .bash_history + +# Python +**/venv/* +**/*.pyc diff --git a/CHANGELOG.md b/CHANGELOG.md index 2db6b4efb45..813327ef9da 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -47,6 +47,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * [#3609](https://github.com/osmosis-labs/osmosis/pull/3609) Add Downtime-detection module. * [#2788](https://github.com/osmosis-labs/osmosis/pull/2788) Add logarithm base 2 implementation. * [#3677](https://github.com/osmosis-labs/osmosis/pull/3677) Add methods for cloning and mutative multiplication on osmomath.BigDec. +* [#3676](https://github.com/osmosis-labs/osmosis/pull/3676) implement `PowerInteger` function on `osmomath.BigDec` +* [#3678](https://github.com/osmosis-labs/osmosis/pull/3678) implement mutative `PowerIntegerMut` function on `osmomath.BigDec`. ### Bug fixes @@ -89,7 +91,6 @@ This release includes stableswap, and expands the IBC safety & composability fun - This PR adds v2 queries for doing the correct thing, and giving people time to migrate from v1beta1 queries to v2. - It also changes cosmwasm to only allow the v2 queries, as no contracts on Osmosis mainnet uses the v1beta1 queries. - ### Bug fixes * [#2803](https://github.com/osmosis-labs/osmosis/pull/2803) Fix total pool liquidity CLI query. diff --git a/osmomath/decimal.go b/osmomath/decimal.go index 55dbd399ff9..cba13ca44a7 100644 --- a/osmomath/decimal.go +++ b/osmomath/decimal.go @@ -989,3 +989,31 @@ func (x BigDec) CustomBaseLog(base BigDec) BigDec { return y } + +// PowerInteger takes a given decimal to an integer power +// and returns the result. Non-mutative. Uses square and multiply +// algorithm for performing the calculation. +func (d BigDec) PowerInteger(power uint64) BigDec { + clone := d.Clone() + return clone.PowerIntegerMut(power) +} + +// PowerIntegerMut takes a given decimal to an integer power +// and returns the result. Mutative. Uses square and multiply +// algorithm for performing the calculation. +func (d BigDec) PowerIntegerMut(power uint64) BigDec { + if power == 0 { + return OneDec() + } + tmp := OneDec() + + for i := power; i > 1; { + if i%2 != 0 { + tmp = tmp.MulMut(d) + } + i /= 2 + d = d.MulMut(d) + } + + return d.MulMut(tmp) +} diff --git a/osmomath/decimal_test.go b/osmomath/decimal_test.go index 1b09c9f0ef3..d881f86ba55 100644 --- a/osmomath/decimal_test.go +++ b/osmomath/decimal_test.go @@ -24,6 +24,20 @@ func TestDecimalTestSuite(t *testing.T) { suite.Run(t, new(decimalTestSuite)) } +// assertMutResult given expected value after applying a math operation, a start value, +// mutative and non mutative results with start values, asserts that mutation are only applied +// to the mutative versions. Also, asserts that both results match the expected value. +func (s *decimalTestSuite) assertMutResult(expectedResult, startValue, mutativeResult, nonMutativeResult, mutativeStartValue, nonMutativeStartValue BigDec) { + // assert both results are as expected. + s.Require().Equal(expectedResult, mutativeResult) + s.Require().Equal(expectedResult, nonMutativeResult) + + // assert that mutative method mutated the receiver + s.Require().Equal(mutativeStartValue, expectedResult) + // assert that non-mutative method did not mutate the receiver + s.Require().Equal(nonMutativeStartValue, startValue) +} + func TestDecApproxEq(t *testing.T) { // d1 = 0.55, d2 = 0.6, tol = 0.1 d1 := NewDecWithPrec(55, 2) @@ -1031,6 +1045,139 @@ func (s *decimalTestSuite) TestCustomBaseLog() { } } +func (s *decimalTestSuite) TestPowerInteger() { + var expectedErrTolerance = MustNewDecFromStr("0.000000000000000000000000000000100000") + + tests := map[string]struct { + base BigDec + exponent uint64 + expectedResult BigDec + + expectedToleranceOverwrite BigDec + }{ + "0^2": { + base: ZeroDec(), + exponent: 2, + + expectedResult: ZeroDec(), + }, + "1^2": { + base: OneDec(), + exponent: 2, + + expectedResult: OneDec(), + }, + "4^4": { + base: MustNewDecFromStr("4"), + exponent: 4, + + expectedResult: MustNewDecFromStr("256"), + }, + "5^3": { + base: MustNewDecFromStr("5"), + exponent: 4, + + expectedResult: MustNewDecFromStr("625"), + }, + "e^10": { + base: eulersNumber, + exponent: 10, + + // https://www.wolframalpha.com/input?i=e%5E10+41+digits + expectedResult: MustNewDecFromStr("22026.465794806716516957900645284244366354"), + }, + "geom twap overflow: 2^log_2{max spot price + 1}": { + base: twoBigDec, + // add 1 for simplicity of calculation to isolate overflow. + exponent: uint64(BigDecFromSDKDec(gammtypes.MaxSpotPrice).Add(OneDec()).LogBase2().TruncateInt().Uint64()), + + // https://www.wolframalpha.com/input?i=2%5E%28floor%28+log+base+2+%282%5E128%29%29%29+++39+digits + expectedResult: MustNewDecFromStr("340282366920938463463374607431768211456"), + }, + "geom twap overflow: 2^log_2{max spot price}": { + base: twoBigDec, + exponent: uint64(BigDecFromSDKDec(gammtypes.MaxSpotPrice).LogBase2().TruncateInt().Uint64()), + + // https://www.wolframalpha.com/input?i=2%5E%28floor%28+log+base+2+%282%5E128+-+1%29%29%29+++39+digits + expectedResult: MustNewDecFromStr("170141183460469231731687303715884105728"), + }, + "geom twap overflow: 2^log_2{max spot price / 2 - 2017}": { // 2017 is prime. + base: twoBigDec, + exponent: uint64(BigDecFromSDKDec(gammtypes.MaxSpotPrice.Quo(sdk.NewDec(2)).Sub(sdk.NewDec(2017))).LogBase2().TruncateInt().Uint64()), + + // https://www.wolframalpha.com/input?i=e%5E10+41+digits + expectedResult: MustNewDecFromStr("85070591730234615865843651857942052864"), + }, + + // sdk.Dec test vectors copied from osmosis-labs/cosmos-sdk: + + "1.0 ^ (10) => 1.0": { + base: OneDec(), + exponent: 10, + + expectedResult: OneDec(), + }, + "0.5 ^ 2 => 0.25": { + base: NewDecWithPrec(5, 1), + exponent: 2, + + expectedResult: NewDecWithPrec(25, 2), + }, + "0.2 ^ 2 => 0.04": { + base: NewDecWithPrec(2, 1), + exponent: 2, + + expectedResult: NewDecWithPrec(4, 2), + }, + "3 ^ 3 => 27": { + base: NewBigDec(3), + exponent: 3, + + expectedResult: NewBigDec(27), + }, + "-3 ^ 4 = 81": { + base: NewBigDec(-3), + exponent: 4, + + expectedResult: NewBigDec(81), + }, + "-3 ^ 50 = 717897987691852588770249": { + base: NewBigDec(-3), + exponent: 50, + + expectedResult: MustNewDecFromStr("717897987691852588770249"), + }, + "-3 ^ 51 = -2153693963075557766310747": { + base: NewBigDec(-3), + exponent: 51, + + expectedResult: MustNewDecFromStr("-2153693963075557766310747"), + }, + "1.414213562373095049 ^ 2 = 2": { + base: NewDecWithPrec(1414213562373095049, 18), + exponent: 2, + + expectedResult: NewBigDec(2), + expectedToleranceOverwrite: MustNewDecFromStr("0.0000000000000000006"), + }, + } + + for name, tc := range tests { + tc := tc + s.Run(name, func() { + + tolerance := expectedErrTolerance + if !tc.expectedToleranceOverwrite.IsNil() { + tolerance = tc.expectedToleranceOverwrite + } + + actualResult := tc.base.PowerInteger(tc.exponent) + require.True(DecApproxEq(s.T(), tc.expectedResult, actualResult, tolerance)) + + }) + } +} + func (s *decimalTestSuite) TestClone() { // The value to change the underlying copy's @@ -1099,16 +1246,51 @@ func (s *decimalTestSuite) TestMul_Mutation() { startNonMut := tc.startValue.Clone() resultMut := startMut.MulMut(mulBy) - result := startNonMut.Mul(mulBy) + resultNonMut := startNonMut.Mul(mulBy) + + s.assertMutResult(tc.expectedMulResult, tc.startValue, resultMut, resultNonMut, startMut, startNonMut) + }) + } +} + +// TestMul_Mutation tests that PowerIntegerMut mutates the receiver +// while PowerInteger is not. +func (s *decimalTestSuite) TestPowerInteger_Mutation() { + + exponent := uint64(2) + + tests := map[string]struct { + startValue BigDec + expectedResult BigDec + }{ + "1": { + startValue: OneDec(), + expectedResult: OneDec(), + }, + "-3": { + startValue: MustNewDecFromStr("-3"), + expectedResult: MustNewDecFromStr("9"), + }, + "0": { + startValue: ZeroDec(), + expectedResult: ZeroDec(), + }, + "4": { + startValue: MustNewDecFromStr("4.5"), + expectedResult: MustNewDecFromStr("20.25"), + }, + } + + for name, tc := range tests { + s.Run(name, func() { + + startMut := tc.startValue.Clone() + startNonMut := tc.startValue.Clone() - // assert both results are as expectde. - s.Require().Equal(tc.expectedMulResult, resultMut) - s.Require().Equal(tc.expectedMulResult, result) + resultMut := startMut.PowerIntegerMut(exponent) + resultNonMut := startNonMut.PowerInteger(exponent) - // assert MulMut mutated the receiver - s.Require().Equal(tc.expectedMulResult, startMut) - // assert Mul did not mutate the receiver - s.Require().Equal(tc.startValue, startNonMut) + s.assertMutResult(tc.expectedResult, tc.startValue, resultMut, resultNonMut, startMut, startNonMut) }) } } diff --git a/osmomath/math.go b/osmomath/math.go index d8f341bb605..0e07cac8216 100644 --- a/osmomath/math.go +++ b/osmomath/math.go @@ -18,6 +18,10 @@ var ( one_half sdk.Dec = sdk.MustNewDecFromStr("0.5") one sdk.Dec = sdk.OneDec() two sdk.Dec = sdk.MustNewDecFromStr("2") + + // https://www.wolframalpha.com/input?i=2.718281828459045235360287471352662498&assumption=%22ClashPrefs%22+-%3E+%7B%22Math%22%7D + // nolint: unused + eulersNumber = MustNewDecFromStr("2.718281828459045235360287471352662498") ) // Returns the internal "power precision". diff --git a/scripts/approximations/README.md b/scripts/approximations/README.md new file mode 100644 index 00000000000..ebadb88198a --- /dev/null +++ b/scripts/approximations/README.md @@ -0,0 +1,65 @@ +# Math Operations Approximations + +## Context + +This is a set of scripts to approximate a mathematical operation using polynomial +and rational approximations. + +The `main` script is in its respective function of `main.py`. It does the following: + +1. Computes polynomial and rational approximations of a given function (e^x by default), +returning the coefficients. + +1. Computes (x,y) coordinates for every kind of approximation given the same x coordinates. +Plots the results for rough comparison. + +1. Plots the results for rough comparison. + +2. Computes the max error for every approximation given the same x coordinates. + +3. Computes and plots max errors for every approximation with a varying number of parameters. + +In other words, this script runs various approximation methods, plots their results and deltas +from actual function values. It can be configured to print the maximum error. +The exact behavior is controlled by the global variables at the top of `main.py`. + +The following are the resources used to create this: +- +- + +In `main.py`, there is also an `exponent_approximation_choice` script. + +This is a shorter and simpler version of `main` that isolates the 13-parameter +Chebyshev Rational approximation of e^x. We are planning to use it in production. +Therefore, we need to perform coefficient truncations to 36 decimal points +(the max `osmomath` supported precision). This truncation is applied +to `exponent_approximation_choice` but not `main`. + +## Configuration + +Several parameters 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 and have Sympy installed: + +```bash +# Create a virtual environment. +python3 -m venv ~/approx-venv + +# Start the environment +source ~/approx-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 +``` diff --git a/scripts/approximations/approximations.py b/scripts/approximations/approximations.py new file mode 100644 index 00000000000..19849698088 --- /dev/null +++ b/scripts/approximations/approximations.py @@ -0,0 +1,191 @@ +import sympy as sp + +import rational +import polynomial +import chebyshev + +def linspace(start: sp.Float, end: sp.Float, num_points: int) -> list[sp.Float]: + """ Given [start, end] on the x-axis, creates the desired number of + equispaced points and returns them in increasing order as a list. + """ + if num_points == 1: + return [(end - start) / 2] + + result: list[sp.Float] = [] + for i in range(num_points): + cur_coord = sp.Float(start + i * (end - start) / (num_points - 1), 200) + + if cur_coord is sp.nan: + raise ValueError("cur_coord is nan") + + result.append(cur_coord) + + return result + +def approx_and_eval_all(approximated_fn, num_parameters: int, x_coordinates) -> tuple: + """ Performs all supported approximations of the given function, evaluates + each wih the given x-coordinates. + + Returns y-coordinates as a tuple in the following order: + - Evaluated Equispaced Polynomial + - Evaluated Chebyshev Polynomial + - Evaluated Chebyshev Rational + - True Y-coordinates as determined by Sympy's implementation of the function. + """ + x_start = x_coordinates[0] + x_end = x_coordinates[-1] + + # Equispaced Polynomial Approximation + coefficients_equispaced_poly = equispaced_poly_approx(approximated_fn, x_start, x_end, num_parameters) + y_eqispaced_poly = polynomial.evaluate(x_coordinates, coefficients_equispaced_poly) + + # Chebyshev Polynomial Approximation + coefficients_chebyshev_poly = chebyshev_poly_approx(approximated_fn, x_start, x_end, num_parameters) + y_chebyshev_poly = polynomial.evaluate(x_coordinates, coefficients_chebyshev_poly) + + # Chebyshev Rational Approximation + numerator_coefficients_chebyshev_rational, denominator_coefficients_chebyshev_rational = chebyshev_rational_approx(approximated_fn, x_start, x_end, num_parameters) + y_chebyshev_rational = rational.evaluate(x_coordinates, numerator_coefficients_chebyshev_rational, denominator_coefficients_chebyshev_rational) + + # Actual + y_actual = get_y_actual(approximated_fn, x_coordinates) + + return (y_eqispaced_poly, y_chebyshev_poly, y_chebyshev_rational, y_actual) + +def get_y_actual(approximated_fn, x_coordinates) -> list[sp.Float]: + """ Given x coordinates and a function, returns a list of y coordinates, one for each x coordinate + that represent the true (x,y) coordinates of the function, as calculated by sp. + """ + y_actual = [] + for x in x_coordinates: + y_actual.append(approximated_fn(x)) + return y_actual + +def compute_max_error(y_approximation, y_actual) -> sp.Float: + """ Given an approximated list of y values and actual y values, computes and returns + the maximum error delta between them. + + CONTRACT: + - y_approximation and y_actual must be the same length + - for every i in range(len(y_approximation)), y_approximation[i] and y_actual[i] must correspond to the + same x coordinate + """ + if len(y_approximation) != len(y_actual): + raise ValueError(F"y_approximation ({len(y_approximation)}) and y_actual ({len(y_actual)}) must be the same length.") + + max: sp.Float = None + for i in range(len(y_approximation)): + cur_abs = sp.Abs(y_approximation[i] - y_actual[i]) + if cur_abs is sp.nan: + raise ValueError(F"cur_abs is nan. y_approximation[i] ({y_approximation[i]}) and y_actual[i] ({y_actual[i]})") + if max is None: + max = cur_abs + else: + max = sp.Max(max, cur_abs) + return max + +def compute_error_range(y_approximation: list, y_actual: list) -> list[sp.Float]: + """ Given an approximated list of y values and actual y values, computes and returns + error deltas between them. + + CONTRACT: + - y_approximation and y_actual must be the same length + - for every i in range(len(y_approximation)), y_approximation[i] and y_actual[i] must correspond to the + same x coordinate + """ + result = [] + for i in range(len(y_approximation)): + cur_abs = sp.Abs(y_approximation[i] - y_actual[i]) + if cur_abs is sp.nan: + raise ValueError(F"cur_abs is nan. y_approximation[i] ({y_approximation[i]}) and y_actual[i] ({y_actual[i]})") + result.append(cur_abs) + return result + +def equispaced_poly_approx(fn, x_start: sp.Float, x_end: sp.Float, num_terms: int): + """ 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 = linspace(x_start, x_end, num_terms) + y_nodes = sp.Matrix([fn(x) for x in equispaced_nodes_x]) + + # Construct a system of linear equations. + vandermonde_matrix = polynomial.construct_vandermonde_matrix(equispaced_nodes_x) + + coef = vandermonde_matrix.solve(y_nodes) + + return coef + +def chebyshev_poly_approx(fn, x_start: int, x_end: int, num_terms: int): + """ 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 = vandermonde_matrix.solve(sp.Matrix(y_estimated)) + + return coef + +def chebyshev_rational_approx(fn, x_start: int, x_end: int, num_parameters: int): + """ 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 + """ + if num_parameters % 2 == 0: + # if num_parameters is 6, we want (3, 4) terms + # assume h(x) = p(x) / q(x) + # we always set the first term of q(x) to 1 for ease of calculation. + # so, we want p(x) to have 3 terms and q(x) to have 4 terms. + num_terms_numerator = num_parameters // 2 + num_terms_denominator = num_parameters // 2 + 1 + else: + # if num_parameters is 5, we want (3, 3) terms + # assume h(x) = p(x) / q(x) + # we always set the first term of q(x) to 1 for ease of calculation. + # so, we want p(x) to have 3 terms and q(x) to have 4 terms. + num_terms_numerator = num_parameters // 2 + 1 + num_terms_denominator = num_parameters // 2 + 1 + + # Compute Chebyshev coordinates. + x_chebyshev, y_chebyshev = chebyshev.get_nodes(fn, x_start, x_end, num_parameters) + + # Construct a system of linear equations. + matrix = rational.construct_rational_eval_matrix(x_chebyshev, y_chebyshev, num_terms_numerator, num_terms_denominator) + + # Solve the matrix to get the coefficients used in the final approximation polynomial. + # coef = np.linalg.solve(np.array(matrix), y_chebyshev) + coef = matrix.solve(sp.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 = [1] + coef_denominator + + return [coef_numerator, coef_denominator] diff --git a/scripts/approximations/chebyshev.py b/scripts/approximations/chebyshev.py new file mode 100644 index 00000000000..0c2986d6f19 --- /dev/null +++ b/scripts/approximations/chebyshev.py @@ -0,0 +1,29 @@ +import sympy as sp +from typing import Tuple + +# N.B. We must evaluate x values early. +# Otherwise, solving for the coefficients will take a long time. +_x_eval_precision = 200 +_x_eval_chop = 1e-200 + +def get_nodes(fn, x_start: sp.Float, x_end: sp.Float, 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 * sp.cos((2*sp.Float(i,_x_eval_precision) + 1).evalf(chop=_x_eval_chop) * sp.pi.evalf(chop=_x_eval_chop) / (2 * sp.Float(num_terms, _x_eval_precision)))) + y_i = fn(x_i) + + x_estimated.append(x_i) + y_estimated.append(y_i) + + return x_estimated, y_estimated diff --git a/scripts/approximations/main.py b/scripts/approximations/main.py new file mode 100644 index 00000000000..b91433855c9 --- /dev/null +++ b/scripts/approximations/main.py @@ -0,0 +1,194 @@ +import matplotlib.pyplot as plt +import sympy as sp + +import approximations +import rational + +########################## +# 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 paramters to use for the approximations. +num_parameters = 13 + +# number of paramters to use for plotting error deltas. +num_parameters_errors = 30 + +# number of (x,y) coordinates used to plot the resulting approximation. +num_points_plot = 100000 + +# function to approximate +approximated_fn = lambda x: sp.Pow(sp.E, x) + +# fixed point precision used in Osmosis `osmomath` package. +osmomath_precision = 36 + +# flag controlling whether to plot each approximation. +# Plots if true. +shouldPlotApproximations = True + +# flag controlling whether to compute max error for each approximation +# given the equally spaced x coordinates. +# Computes if true. +shouldComputeErrorDelta = True + +# flag controlling whether to plot errors over a range. +# Currently, does so only for Chebyshev Rational Approximation. +# Computes if true. +shouldPlotErrorRange = True + +# flag controlling whether to plot max error for every approximation +# with a varying number of parameters. This is useful to find the most +# optimal number of parameters to use for each kind of approximation. +# Plots if true. +shouldPlotMaxError = True + +def plot_error_range(x_coordinates, y_approximation, y_actual): + """ Given x coordinates that correspond to approximated y coordinates and actual y coordinates, + compute the deltas between y approximated and y actual and plot them in log scale on y. + """ + error_deltas = approximations.compute_error_range(y_approximation, y_actual) + + plt.semilogy(x_coordinates, error_deltas) + + plt.grid(True) + plt.title(f"Chebyshev Rational e^x Errors on [{x_start}, {x_end}]. {num_parameters} params, {num_points_plot} points") + plt.show() + +# 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. +# - Computes the max error for every approximation given the same x coordinates. +# - Computes and plots max errors for every approximation with a varying number of parameters. +# This script runs various approximation methods, plots their results and deltas +# from actual function values. The script can also be configured to print the maximum error. +# The exact behavior is controlled by the global +# variables at the top of the file. +# The following are the resources used: +# https://xn--2-umb.com/22/approximation/ +# https://sites.tufts.edu/atasissa/files/2019/09/remez.pdf +def main(): + # Equispaced x coordinates to be used for plotting every approximation. + x_coordinates = approximations.linspace(x_start, x_end, num_points_plot) + + if shouldComputeErrorDelta or shouldPlotApproximations or shouldPlotErrorRange: + ############################################### + # Approximation With Given Number of Parameters + y_eqispaced_poly, y_chebyshev_poly, y_chebyshev_rational, y_actual = approximations.approx_and_eval_all(approximated_fn, num_parameters, x_coordinates) + + ################ + # 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_parameters} parameters used\n\n") + + # Equispaced Polynomial Approximation + max_error_equispaced_poly = approximations.compute_max_error(y_eqispaced_poly, y_actual) + print(F"Equispaced Poly: {max_error_equispaced_poly.evalf(chop=1e-100)}") + + # Chebyshev Polynomial Approximation + max_error_chebyshev_poly = approximations.compute_max_error(y_chebyshev_poly, y_actual) + print(F"Chebyshev Poly: {max_error_chebyshev_poly.evalf(chop=1e-100)}") + + # Chebyshev Rational Approximation + max_error_chebyshev_rational = approximations.compute_max_error(y_chebyshev_rational, y_actual) + print(F"Chebyshev Rational: {max_error_chebyshev_rational.evalf(chop=1e-100)}") + + if shouldPlotErrorRange: + plot_error_range(x_coordinates, y_chebyshev_rational, y_actual) + + ############################### + # Plot Every Approximation Kind + if shouldPlotApproximations: + # Equispaced Polynomial Approximation + plt.plot(x_coordinates, y_eqispaced_poly, label="Equispaced Poly") + + # Chebyshev Polynomial Approximation + plt.plot(x_coordinates, y_chebyshev_poly, label="Chebyshev Poly") + + # Chebyshev Rational Approximation + plt.plot(x_coordinates, y_chebyshev_rational, label="Chebyshev Rational") + + # Actual With Large Number of Coordinates (evenly spaced on the X-axis) + plt.plot(x_coordinates, y_actual, label=F"Actual") + + plt.legend(loc="upper left") + plt.grid(True) + plt.title(f"Appproximation of e^x on [{x_start}, {x_end}] with {num_parameters} parameters") + plt.show() + + ##################################################### + # Calculate Errors Given Varying Number of Parameters + if shouldPlotMaxError: + x_axis = [] + + deltas_eqispaced_poly = [] + deltas_chebyshev_poly = [] + deltas_chebyshev_rational = [] + + ################ + # Compute Deltas + # The deltas are taken from actual function values for different number of parameters + # This is needed to find the most optimal number of parameters to use. + for num_parameters_current in range(1, num_parameters_errors + 1): + x_axis.append(int(num_parameters_current)) + y_eqispaced_poly, y_chebyshev_poly, y_chebyshev_rational, y_actual = approximations.approx_and_eval_all(approximated_fn, num_parameters_current, x_coordinates) + + deltas_eqispaced_poly.append(approximations.compute_max_error(y_eqispaced_poly, y_actual)) + deltas_chebyshev_poly.append(approximations.compute_max_error(y_chebyshev_poly, y_actual)) + deltas_chebyshev_rational.append(approximations.compute_max_error(y_chebyshev_rational, y_actual)) + + ################## + # Plot the results + + # Equispaced Polynomial Approximation + plt.semilogy(x_axis, deltas_eqispaced_poly, label="Equispaced Poly") + + # Chebyshev Polynomial Approximation + plt.semilogy(x_axis, deltas_chebyshev_poly, label="Chebyshev Poly") + + # Chebyshev Rational Approximation + plt.semilogy(x_axis, deltas_chebyshev_rational, label="Chebyshev Rational") + + plt.legend(loc="upper left") + plt.grid(True) + plt.title(f"Approximation Errors on [{x_start}, {x_end}]") + plt.gca().invert_yaxis() + plt.xlabel('Number of Parameters') + plt.ylabel(F"-log_10{{ max | f'(x) - f(x) | }} where x is in [{x_start}, {x_end}]") + plt.show() + +# This script isolates the 13-parameter Chebyshev Rational approximation of e^x +# We are planning to use it in production. Therefore, we need to peform coefficient +# truncations to 36 decimal points (the max osmomath supported precision). +def exponent_approximation_choice(): + # Equispaced x coordinates to be used for plotting every approximation. + x_coordinates = approximations.linspace(x_start, x_end, num_points_plot) + x_coordinates = [sp.Float(sp.N(coef, osmomath_precision + 1), osmomath_precision + 1) for coef in x_coordinates] + + # Chebyshev Rational Approximation to get the coefficients. + coef_numerator, coef_denominator = approximations.chebyshev_rational_approx(approximated_fn, x_start, x_end, num_parameters) + + # Truncate the coefficients to osmomath precision. + coef_numerator = [sp.Float(sp.N(coef, osmomath_precision + 1), osmomath_precision + 1) for coef in coef_numerator] + coef_denominator = [sp.Float(sp.N(coef, osmomath_precision + 1), osmomath_precision + 1) for coef in coef_denominator] + + # Evaluate approximation. + y_chebyshev_rational = rational.evaluate(x_coordinates, coef_numerator, coef_denominator) + + # Compute Actual Values + y_actual = approximations.get_y_actual(approximated_fn, x_coordinates) + + plot_error_range(x_coordinates, y_chebyshev_rational, y_actual) + +if __name__ == "__main__": + # Uncomment to run the main script. + #main() + exponent_approximation_choice() diff --git a/scripts/approximations/polynomial.py b/scripts/approximations/polynomial.py new file mode 100644 index 00000000000..5c6af20431a --- /dev/null +++ b/scripts/approximations/polynomial.py @@ -0,0 +1,41 @@ +import sympy as sp + +def construct_vandermonde_matrix(x_list: list[sp.Float]) -> sp.Matrix: + """ 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(sp.Pow(x_list[i], j)) + matrix.append(row) + + return sp.Matrix(matrix) + +def evaluate(x, coeffs): + """ 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. + """ + y = [] + for x_i in x: + y_i = 0 + for i in range(len(coeffs)): + y_i += coeffs[i]*sp.Pow(x_i, i) + y.append(y_i) + return y diff --git a/scripts/approximations/rational.py b/scripts/approximations/rational.py new file mode 100644 index 00000000000..4adea2d9188 --- /dev/null +++ b/scripts/approximations/rational.py @@ -0,0 +1,44 @@ +import sympy as sp + +import polynomial + +def construct_rational_eval_matrix(x_list: list, y_list: list, num_terms_numerator: int, num_terms_denominator) -> sp.Matrix: + """ 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}] + """ + matrix = [] + + for i in range(num_terms_numerator + num_terms_denominator - 1): + row = [] + for j in range(num_terms_numerator): + row.append(sp.Pow(x_list[i], j)) + + for j in range(num_terms_denominator): + # denominator terms + if j > 0: + row.append(-1 * sp.Pow(x_list[i], j) * y_list[i]) + + matrix.append(row) + + return sp.Matrix(matrix) + +def evaluate(x: list, coefficients_numerator: list, coefficients_denominator: list): + """ 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. + """ + p = polynomial.evaluate(x, coefficients_numerator) + q = polynomial.evaluate(x, coefficients_denominator) + + result = [] + + for i in range(len(x)): + result.append(p[i] / q[i]) + + return result diff --git a/scripts/approximations/rational_test.py b/scripts/approximations/rational_test.py new file mode 100644 index 00000000000..099b15d5ac9 --- /dev/null +++ b/scripts/approximations/rational_test.py @@ -0,0 +1,36 @@ +import unittest +import math +import sympy as sp + +import rational + +class TestChebyshevRational(unittest.TestCase): + + # assume (3,3) rational + # h(x) = p(x) / q(x) = (p_0 + p_1 * x + p_2 * x^2) / (1 + q_1 * x + q_2 * x^2) + # p_0 + p_1 * x + p_2 * x^2 - q_1 * y * x - q_2 * y * x^2 = y + # construct 5 x 5 matrix solving for p and q given x_0 to x_4 and y_0 to y_4 + # assume function is y = e**x + def test_construct_matrix_3_3(self): + fn = lambda x: math.e**x + x = [1, 2, 3, 4, 5] + y = list(map(fn, x)) + + coeffs = rational.construct_rational_eval_matrix(x, y, 3, 3) + + # correct matrix size + self.assertEqual(len(x) * len(x), len(coeffs)) + + # first row is correct + for i in range(len(x)): + x_i = x[i] + y_i = y[i] + + expected_row = [sp.Pow(x_i, 0), sp.Pow(x_i, 1), sp.Pow(x_i, 2), -1 * sp.Pow(x_i, 1) * y_i, -1 * sp.Pow(x_i, 2) * y_i] + + actual_row = coeffs.row(i).tolist()[0] + + self.assertEqual(expected_row, actual_row) + +if __name__ == '__main__': + unittest.main() diff --git a/scripts/approximations/requirements.txt b/scripts/approximations/requirements.txt new file mode 100644 index 00000000000..ce1140c3e00 --- /dev/null +++ b/scripts/approximations/requirements.txt @@ -0,0 +1,16 @@ +contourpy==1.0.6 +cycler==0.11.0 +fonttools==4.38.0 +kiwisolver==1.4.4 +matplotlib==3.6.2 +mpmath==1.2.1 +numpy==1.23.5 +packaging==22.0 +Pillow==9.3.0 +pyparsing==3.0.9 +PyQt5==5.15.7 +PyQt5-Qt5==5.15.2 +PyQt5-sip==12.11.0 +python-dateutil==2.8.2 +six==1.16.0 +sympy==1.11.1