Loss of precision with `+=`
Opened this issue · 11 comments
It looks like the implementation of the += operator for function_trees leads to precision loss.
The following code (taken from ReMRChem) gives two very different results for rho and n which ideally should be identical.
def setup(self):
rho = vp.FunctionTree(self.mra)
rho.setZero()
rholist = []
for i in range(0, len(self.Psi)):
dens = self.Psi[i].overlap_density(self.Psi[i], self.prec)
print("density i = ", i)
print(dens)
rho += dens.real
rholist.append(dens.real)
#rho.crop(self.prec)
n = rholist[0] + rholist[1]
self.potential = (4.0*np.pi) * self.poisson(n).crop(self.prec)
Commit of ReMRChem showing the issue: 2b6052d94ea4f68715f840b62f9cb7
Suggestions:
- reimplement
__iadd__to mirror__add__ - keep
__add__written with pybind11, remove__iadd__from the pybind11 binding set, re-implement it in pure Python to use__add__and then dynamically attach it to the FunctionTree class.
Number 2 is incorrect: you want to keep __add__ written with pybind11, remove __iadd__ from the pybind11 binding set, re-implement it in pure Python to use __add__ and then dynamically attach it to the FunctionTree class.
Updated your suggestion above
Hmm, I like the dynamically add solution. I guess we need to do something like this then:
https://pybind11.readthedocs.io/en/stable/classes.html?highlight=dynamically%20add#dynamic-attributes
Where would you suggest putting the python code @robertodr ?
A third alternative would be to write:
while (refine_grid(*out, *inp)); as Stig mentions in Zulip. That would be correcting the current implementation.
It's true that Stig's suggestion (probably) solves the issue, though using two different implementations for the same things it's not good.1
As to how to implement dynamic attributes in practice. You'd have to add a Python file under src/vampyr, say _function_tree.py. The drawback of this approach is that you need to repeat the same code for D=1,2,3 🤦🏻
Footnotes
-
And I realize that I should have spotted this earlier during code review! ↩
@robertodr can you write an example of how to use it for a method in f_tree method in 1D?
I can do it for the other methods and dims once I have the blueprint
I tried to do it with a "hello_world" kind of example, but didn't get it to work properly.
On vacation until 5th July 🥥🌴
There is a fix for addition on master now. All other overloaded operators should be double-checked (and in case re-implemented) before closing this issue.
I can try, but then I need to compile from scratch. Not sure I'm able to 😉