[RFC]: add `math/base/tools/normhermitepolyf`
kgryte opened this issue · 9 comments
Description
This RFC proposes adding the package @stdlib/math/base/tools/normhermitepolyf for evaluating a normalized Hermite polynomial using single-precision floating-point arithmetic.
Package: @stdlib/math/base/tools/normhermitepolyf
Alias: normhermitepolyf
Related Issues
None.
Questions
No.
Other
This RFC has prior art. Namely,
@stdlib/math/base/tools/normhermitepoly: uses double-precision floating-point arithmetic@stdlib/math/base/tools/evalpolyf: emulates single-precision floating-point arithmetic usingfloat64ToFloat32.
This RFC proposes effectively copying the normhermitepoly package and modifying the implementation, tests, and examples (as in evalpolyf) to use float64ToFloat32 for emulating single-precision arithmetic.
Checklist
- I have read and understood the Code of Conduct.
- Searched for existing issues and pull requests.
- The issue name begins with
RFC:.
Hi @kgryte , can I work on this? I will try it. Thank you!
@Daniel777y As @xaman27x is working on #2018, feel free to work on this issue. I'll go ahead and assign you.
Hi @kgryte , I've implemented all the files but didn't pass the test (about 5% not passed). I guess I get stuck somewhere and need help.
I am not quite sure if it's because of the precision of test cases or of my implementation. Do I need to convert those cases to np.float32 in python?
Here's my implementation:
function normhermitepolyf( n, x ) {
var y1;
var y2;
var y3;
var i;
x = float64ToFloat32( x );
if ( isnan( n ) || isnanf( x ) || n < 0 || !isint( n ) ) {
return NaN;
}
if ( n === 0 ) {
// `x` is completely canceled from the expression:
return 1.0;
}
if ( n === 1 ) {
return x;
}
y2 = 1.0;
y3 = 0.0;
for ( i = n; i > 1; i-- ) {
y1 = float64ToFloat32( float64ToFloat32( x * y2 ) - float64ToFloat32( i * y3 ) );
y3 = y2;
y2 = y1;
}
return float64ToFloat32( float64ToFloat32( x * y2 ) - y3 );
}Here're the fixtures generator:
"""Generate fixtures."""
import os
import json
import numpy as np
from scipy.special import eval_hermitenorm
# Get the file path:
FILE = os.path.realpath(__file__)
# Extract the directory in which this file resides:
DIR = os.path.dirname(FILE)
def gen(n, x, name):
"""Generate fixture data and write to file.
# Arguments
* `n`: degree(s)
* `x`: domain
* `name::str`: output filename
# Examples
``` python
python> n = 1
python> x = linspace(-1000, 1000, 2001)
python> gen(n, x, './data.json')
```
"""
y = eval_hermitenorm(n, x)
if isinstance(n, np.ndarray):
data = {
"n": n.tolist(),
"x": x.tolist(),
"expected": y.tolist()
}
else:
data = {
"n": n,
"x": x.tolist(),
"expected": y.tolist()
}
# Based on the script directory, create an output filepath:
filepath = os.path.join(DIR, name)
# Write the data to the output filepath as JSON:
with open(filepath, "w", encoding="utf-8") as outfile:
json.dump(data, outfile)
def main():
"""Generate fixture data."""
# Random values across `n` and `x`:
n = np.random.randint(1, 100, 1000)
x = np.float32(np.random.random(1000)*100.0)
gen(n, x, "random2.json")
# Medium negative:
x = np.linspace(-709.78, -1.0, 1000, dtype=np.float32)
gen(1, x, "medium_negative_1.json")
gen(2, x, "medium_negative_2.json")
gen(5, x, "medium_negative_5.json")
# Medium positive:
x = np.linspace(1.0, 709.78, 1000, dtype=np.float32)
gen(1, x, "medium_positive_1.json")
gen(2, x, "medium_positive_2.json")
gen(5, x, "medium_positive_5.json")
# Small positive:
x = np.linspace(2.0**-51, 1.0, 1000, dtype=np.float32)
gen(1, x, "small_positive_1.json")
gen(2, x, "small_positive_2.json")
gen(5, x, "small_positive_5.json")
# Small negative:
x = np.linspace(-1.0, -2.0**-51, 1000, dtype=np.float32)
gen(1, x, "small_negative_1.json")
gen(2, x, "small_negative_2.json")
gen(5, x, "small_negative_5.json")
# Tiny values:
x = np.linspace(-2.0**-51, 2.0**-51, 1000, dtype=np.float32)
gen(1, x, "tiny_1.json")
gen(2, x, "tiny_2.json")
gen(5, x, "tiny_5.json")
if __name__ == "__main__":
main()Here's a unit of tests:
tape( 'the function accurately computes `He1(x)` for small positive numbers', function test( t ) {
var expected;
var delta;
var tol;
var x;
var v;
var i;
var n;
var e;
n = smallPositive1.n;
x = smallPositive1.x;
expected = smallPositive1.expected;
for ( i = 0; i < x.length; i++ ) {
v = normhermitepoly( n, x[ i ] );
// I also convert the expected value to float32 here
e = float64ToFloat32( expected[ i ] );
delta = abs( v - e );
tol = EPS * abs( e );
t.ok( delta <= tol, 'within tolerance. n: ' + n + '. x: ' + x[ i ] + '. Value: ' + v + '. Expected: ' + expected[ i ] + '. Delta: ' + delta + '. Tolerance: ' + tol + '.' );
}
t.end();
});If I need to include more information or there're other examples for reference, please let me know.
Thank you!
@Daniel777y You can tweak the tolerance level a bit. E.g., 2.0 * EPS. That is fine.
@Daniel777y Your implementation looks sound. It's possible that SciPy upcasts inputs to float64, resulting in higher intermediate precision. As long as the accuracy is not significantly far off, we can live with some minor divergence.
Thank you for response @kgryte ! After tweaking the tolerance level, it's better.
But because of float32, some test cases would be large (or small) enough to be Infinity (or -Infinity). The calculation contains subtraction, so the return value can be NaN. Do I need to make sure that it will not return NaN?
Here're some potential ideas to pass tests:
- Don't convert middle variables, only convert the input and output to float32. Then if the return value and expected value are both
Infinity, accept it. - Adjust the test functions, try to figure out the whether the
NaNis due toInfinitysubtraction, and then accept those with expected value ofInfinity. - Simply not include large and small cases temporarily.
But for now I am not sure if these methods can accurately calculate the results.
I also walked through the packages evalpolyf and evalrationalf. But the range of the test cases are kinda small, so they might not include the Infinity and NaN cases.
Or maybe there're other good ideas to deal with this problem, so I am looking forward to your suggestions.