Infeasible solution for EMD but inputs are in the simplex?

Describe the bug
A call to the EMD solver is raising a UserWarning about the problem being infeasible (Check that a and b are in the simplex). I dug into the code and this seems to be raised only when a value in one of the two inputs is less than zero. However, I am printing the sum, max, and min of the inputs and don't see any odd behavior (sum to one, max values not greater than 1, min values not less than zero; inputs are the same size; type is np.float64).

Could this be an issue with underflow?

To Reproduce
Steps to reproduce the behavior: using python3.7:

import ot
import torch
import random
import numpy as np

# Get distributions
random_dist = torch.softmax(torch.tensor([random.random() * 1000 for _ in range(625)]), dim=0).numpy()
random_dist2 = torch.softmax(torch.tensor([random.random() * 1000 for _ in range(625)]), dim=0).numpy()
print('--- First dist ---')
print('Shape: %s' % random_dist.shape)
print('Sum: %s' % np.sum(random_dist))
print('Max: %s' % random_dist.max())
print('Min: %s' % random_dist.min())
print('--- Second dist ---')
print('Shape: %s' % random_dist2.shape)
print('Sum: %s' % np.sum(random_dist2))
print('Max: %s' % random_dist2.max())
print('Min: %s' % random_dist2.min())

# Create distance matrix
coords = np.reshape(np.asarray(np.meshgrid(np.linspace(0, 24, 25), np.linspace(0, 24, 25)))[[1, 0], :, :].transpose((1, 2, 0)), [-1, 2])
tiled = np.tile(coords[:, np.newaxis, :], [1, coords.shape[0], 1])
ttiled = tiled.transpose((1, 0, 2))
distances = np.array(np.linalg.norm(tiled - ttiled, axis=2), order='C')

# Compute EMD
ot.emd(random_dist, random_dist2, distances)

Example output:

--- First dist ---
Shape: 625
Sum: 1.0
Max: 0.42820615
Min: 0.0
--- Second dist ---
Shape: 625
Sum: 1.0
Max: 0.96281725
Min: 0.0
/home/.../venv/lib/python3.7/site-packages/ot/lp/__init__.py:113: UserWarning: Problem infeasible. Check that a and b are in the simplex
  result_code_string = check_result(result_code)

Expected behavior
I would expect that this error is not thrown as all of the values in both distributions are nonnegative.

Screenshots
N/A, minimal working example above

Desktop (please complete the following information):

OS: I can reproduce it on Linux (Ubuntu 16.04.1) and Mac OS (10.14.6)
Python version: 3.7.0 on Linux; 3.7.5 on Mac
How was POT installed: pip and conda (two different machines)

Output of the following code snippet:

import platform; print(platform.platform())
import sys; print("Python", sys.version)
import numpy; print("NumPy", numpy.__version__)
import scipy; print("SciPy", scipy.__version__)
import ot; print("POT", ot.__version__)

Linux:

Linux-4.15.0-65-generic-x86_64-with-Ubuntu-16.04-xenial
Python 3.7.0 (default, Aug  6 2018, 20:07:46) 
[GCC 5.4.0 20160609]
NumPy 1.17.3
SciPy 1.4.1
POT 0.6.0

Mac:

Darwin-18.7.0-x86_64-i386-64bit
Python 3.7.5 (default, Oct 25 2019, 10:52:18) 
[Clang 4.0.1 (tags/RELEASE_401/final)]
NumPy 1.18.1
SciPy 1.4.1
POT 0.6.0

Additional context
The scale for random.random() can be as low as 10 and the bug will happen. If the scale is 1, it can find a solution. Sounds like a numerical underflow problem to me.

Thanks!