Loss of precision with `+=`

[ilfreddy]

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

[ilfreddy]

Suggestions:

1. reimplement __iadd__ to mirror __add__
2. 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.

[robertodr]

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.

[ilfreddy]

Updated your suggestion above

[bjorgve]

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.

[robertodr]

It's true that Stig's suggestion (probably) solves the issue, though using two different implementations for the same things it's not good.¹

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

1. And I realize that I should have spotted this earlier during code review! ↩

[bjorgve]

@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.

[robertodr]

On vacation until 5th July 🥥🌴

[robertodr]

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.

[robertodr]

@ilfreddy can you check if this is still an issue?

[robertodr]

@ilfreddy can you check if this is still an issue? If it isn't, please close.