Divided by zero encountered when fitting a RL model

I think I have a hypothesis. So in my likelihood function, there is a line:
utility_ = reward_*split_
All three are aesara tensors but reward_ may include NaN. Is it true that when I simply call the likelihood function, this line works fine because just like in numpy, NaN times anything will still give you NaN. However, pymc requires gradient computation from it, in which case NaN may break it. Is it true? If so, is there a NaN safe way to multiply two tensors?