0

I am trying to stretch mesh points using density functions. For example, given the following point distribtion (uniformally distributed):

enter image description here

The code below will change the distribution to something like this: enter image description here

import numpy as np
from scipy.integrate import quad
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
from scipy.special import erf

# the x interval limits
a = 0.0
b = 3.0

# The density function, normalize it
_density_func = lambda x: 5*x
density_func = lambda x: _density_func(x) / quad(_density_func, a, b)[0]*(b-a)

# The xarr and yarr
Npts = 50
xarr = np.linspace(a, b, Npts)
yarr = np.zeros_like(xarr)

# Calculate the primitive function F = integral of density_func, from a, to t normalized by the int(density_func, a, b)
# F = np.vectorize(lambda t: quad(density_func, a, t)[0] / quad(density_func, a, b)[0])
F = np.vectorize(lambda t: quad(density_func, a, t)[0])

# if debug is true, print some info
debug = True
if debug:
    print('The range of xarr is: [', a, b, ']')
    print('The range of F(xarr) is: [', F(xarr).min(), F(xarr).max(), ']')

# Calculate the new x distribution of points using the inverse function of F.
# Use the trick of interpolation to calculate the inverse function,i.e: interp1d(y, x)(x)
xnew = interp1d(F(xarr), xarr)(xarr) 

# plot the result
plt.scatter(xnew, yarr, facecolors='none', edgecolors='black')
plt.show()

When I run this script, I get the following error:

The range of xarr is: [ 0.0 3.0 ]
The range of F(xarr) is: [ 0.0 2.9999999999999996 ]
Traceback (most recent call last):
  File "meshDensity.py", line 38, in <module>
    xnew = interp1d(F(xarr), xarr)(xarr)
  File "C:\Users\navaro\Anaconda3\lib\site-packages\scipy\interpolate\polyint.py", line 79, in __call__
    y = self._evaluate(x)
  File "C:\Users\navaro\Anaconda3\lib\site-packages\scipy\interpolate\interpolate.py", line 664, in _evaluate
    below_bounds, above_bounds = self._check_bounds(x_new)
  File "C:\Users\navaro\Anaconda3\lib\site-packages\scipy\interpolate\interpolate.py", line 696, in _check_bounds
    raise ValueError("A value in x_new is above the interpolation "
ValueError: A value in x_new is above the interpolation range.

As you can see, the problem is the right limit of F(xarr) is 2.9999999999999996 instead of the exact value 3.0.

Could you please suggest any solution to this rounding error problem? I appreciate any help.

Edit: a temporary solution is to use mpmath.quad function with mpmath.mp.dps = 20 but this makes the script relatively slow.

Community
  • 1
  • 1

2 Answers2

0

This kind of problem is inherent to using floating point numbers on a computer: you have to cut off the decimals somewhere

As you can see, the error is 4.440892098500626e-16, which is probably around the floating point accuracy on your system

The mathematical approach is probably correct, but I would suggest handling this through software, i.e. do you really need 16 decimal places? If you put a couple round(num,6) in there, you would make your code more robust to these kinds of issues

More specifically, your code runs if I replace your density func with:

density_func = lambda x: round(_density_func(x) / quad(_density_func, a, b)[0]*(b-a),6)

EDIT: funnily enough, I ran into the same issue just now: 3.0-2.9999999999999996 = 4.440892098500626e-16

where we both know that it should have been 0.0000000000000004

Funky things happen when you're dealing with numbers at the edge of your floating point accuracy

imochoa
  • 736
  • 4
  • 8
  • Thank you for your suggestion. unfortunately I couldn't understand how could this solve the problem. –  Oct 29 '19 at 15:26
  • `round` is not a panacea, because the rounded value may still not be perfectly representable in binary. See https://stackoverflow.com/q/588004/5987 – Mark Ransom Oct 29 '19 at 15:28
  • @IamNotaMathematician I updated the post to mark the place where the rounding would go and now it runs on my computer. But take a minute to understand what's going on because this kind of issue is pretty common – imochoa Oct 29 '19 at 15:36
  • @imochoa: Sorry, I don't think this is a solution because my script is a part of a larger program that needs to deal with meshes with high precision. –  Oct 29 '19 at 15:38
  • @IamNotaMathematician well the sad truth is that floating point operations are not exact, your accuracy is going to be limited by the datatype you're using (see https://stackoverflow.com/questions/38477908/smallest-positive-float64-number) so, you can take 10 decimal places instead of 6 but if you're dealing with smaller numbers than that.... consider changing units, working on log scale or something a long those lines – imochoa Oct 29 '19 at 15:48
  • @imochoa: I will use `mpmath.quad` instead of `scipy.integrate.quad` –  Oct 29 '19 at 15:49
0

I solved my problem using arbitrary precision arithmetics module, mpmath.

import numpy as np
import mpmath as mp
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

mp.mp.dps = 18
# the x interval limits
a = 0.0
b = 3.0

# The density function, normalize it
_density_func = lambda x: 5*x
density_func = lambda x: _density_func(x) / mp.quad(_density_func, [a, b])*(b-a)

# The xarr and yarr
Npts = 50
xarr = np.linspace(a, b, Npts)
yarr = np.zeros_like(xarr)

# Calculate the primitive function F = integral of density_func, from a, to t.
F = np.vectorize(lambda t: mp.quad(density_func, [a, t]))

# if debug is true, print some info
debug = True
if debug:
    print('The range of xarr is: [', a, b, ']')
    print('The range of F(xarr) is: [', F(xarr).min(), F(xarr).max(), ']')

# Calculate the new x distribution of points using the inverse function of F.
# Use the trick of interpolation to calculate the inverse function,i.e: interp1d(y, x)(x)
xnew = interp1d(F(xarr), xarr)(xarr) 

# plot the result
plt.scatter(xnew, yarr, facecolors='none', edgecolors='black')
plt.show()

After running the script, I get:

The range of xarr is: [ 0.0 3.0 ]
The range of F(xarr) is: [ 0.0 3.0 ]

enter image description here