I think that the `NaN`s originated not in the computation of `softmax`, but more probably in the gradient of the expression, and scaling the softmax after its computation would solve the problem, much like changing the value of the argument of softmax currently addresses the issue.
The `rescale_softmax` function would look like:
``` Python
def rescale_softmax(sm, min_val):
    n_classes = sm.shape[-1]
    return sm * (1 - n_classes * min_val) + min_val
```
Then you can use, in the base case, `rescale_softmax(plushmax(x), min_val)`, or, in the case where `eps > 0` (if that even makes sense to use that):
```Python
n_classes = x.shape[-1]
pm = plushmax(x, eps)
mv = max(min_val, (min_val - pm.min()) / (1 - n_classes * pm.min()))
out = rescale_softmax(pm, min_val=mv)
```