IsotopeQuantity test fails with datetime assert

[cosama]

All recent test have been failing. The error message is:

E       AssertionError: assert datetime.datetime(2021, 1, 4, 23, 42, 22, 196186) == datetime.datetime(2021, 1, 4, 23, 42, 34, 811280)
E        +  where datetime.datetime(2021, 1, 4, 23, 42, 22, 196186) = <bound method IsotopeQuantity.time_when of <becquerel.tools.isotope_qty.IsotopeQuantity object at 0x7efd1c4023d0>>(**{'uci': 10.047000000000002})
E        +    where <bound method IsotopeQuantity.time_when of <becquerel.tools.isotope_qty.IsotopeQuantity object at 0x7efd1c4023d0>> = <becquerel.tools.isotope_qty.IsotopeQuantity object at 0x7efd1c4023d0>.time_when
E        +  and   datetime.datetime(2021, 1, 4, 23, 42, 34, 811280) = <becquerel.tools.isotope_qty.IsotopeQuantity object at 0x7efd1c4023d0>.ref_date

Anybody has an idea what goes on here. It looks like the python2 test still work, so maybe a version issue with one of the libraries?

[markbandstra]

Weird, it looks like the 2.7 and 3.6 tests are passing but the 3.7 and 3.8 tests are failing. The full failure message for the most recent test of main is:

___________________ test_isotopequantity_time_when[K-40-uci] ___________________
641
642iq = <becquerel.tools.isotope_qty.IsotopeQuantity object at 0x7f46b09c3850>
643kw = 'uci'
644
645    @pytest.mark.parametrize('kw', ('atoms', 'g', 'bq', 'uci'))
646    def test_isotopequantity_time_when(iq, kw):
647        """Test IsotopeQuantity.time_when()"""
648    
649        ref_qty = getattr(iq, kw + '_at')(iq.ref_date)
650    
651        kwarg = {kw: ref_qty}
652>       assert iq.time_when(**kwarg) == iq.ref_date
653E       AssertionError: assert datetime.datetime(2020, 12, 30, 21, 20, 55, 549746) == datetime.datetime(2020, 12, 30, 21, 21, 8, 164840)
654E        +  where datetime.datetime(2020, 12, 30, 21, 20, 55, 549746) = <bound method IsotopeQuantity.time_when of <becquerel.tools.isotope_qty.IsotopeQuantity object at 0x7f46b09c3850>>(**{'uci': 10.047000000000002})
655E        +    where <bound method IsotopeQuantity.time_when of <becquerel.tools.isotope_qty.IsotopeQuantity object at 0x7f46b09c3850>> = <becquerel.tools.isotope_qty.IsotopeQuantity object at 0x7f46b09c3850>.time_when
656E        +  and   datetime.datetime(2020, 12, 30, 21, 21, 8, 164840) = <becquerel.tools.isotope_qty.IsotopeQuantity object at 0x7f46b09c3850>.ref_date
657
658tests/isotope_qty_test.py:201: AssertionError

Which is a failure in test_isotopequantity_time_when. I am able to replicate this locally with 3.8. I'll look into it, although it's weird that this test is failing after being in the code since May 2017.

[markbandstra]

Very strange error. Of these isotopes that the test is run against:

    'Bi-212',
    'Cs-137',
    'K-40',
    'Cs-134',
    'N-13',
    'Rn-220',
    'Tc-99m',
    'Tl-208',
    'Th-232',

and these kwargs: ('atoms', 'g', 'bq', 'uci'), only that single test (K-40 with uci) is failing. No idea yet but I suspect a subtle floating point issue.

[markbandstra]

I printed out some debugging information and was able to confirm a weird floating point issue.

Here is what I printed out for the code as-is:

isotope_qty_test: iq: radioisotope: K-40
isotope_qty_test: iq: date:         2021-01-04 17:33:28.697044
isotope_qty_test: iq: kwargs:       {'uci': 10.047}
isotope_qty_test: test_isotopequantity_time_when: iq:                    371739.00000000006 Bq of K-40 (at 2021-01-04 17:33:28.697044)
isotope_qty_test: test_isotopequantity_time_when: ref_qty:               10.047000000000002
isotope_qty_test: test_isotopequantity_time_when: kwarg:                 {'uci': 10.047000000000002}
IsotopeQuantity.time_when: ref_atoms: 2.1119730744881784e+22
IsotopeQuantity.time_when: target:    2.1119730744881788e+22
IsotopeQuantity.time_when: dt:        -12.615093572023143
isotope_qty_test: test_isotopequantity_time_when: iq.time_when(**kwarg): 2021-01-04 17:33:16.081950
isotope_qty_test: test_isotopequantity_time_when: ref_date:              2021-01-04 17:33:28.697044
IsotopeQuantity.time_when: ref_atoms: 2.1119730744881784e+22
IsotopeQuantity.time_when: target:    2.1119730744881788e+22
IsotopeQuantity.time_when: dt:        -12.615093572023143

(The error occurs because that last line for dt is not 0.0.)

Then I increased the activity in uCi by a factor of 10 from 10.047 to 100.47 in the test fixture and the error went away!

isotope_qty_test: iq: radioisotope: K-40
isotope_qty_test: iq: date:         2021-01-04 17:32:17.773989
isotope_qty_test: iq: kwargs:       {'uci': 100.47}
isotope_qty_test: test_isotopequantity_time_when: iq:                    3717390.0 Bq of K-40 (at 2021-01-04 17:32:17.773989)
isotope_qty_test: test_isotopequantity_time_when: ref_qty:               100.47
isotope_qty_test: test_isotopequantity_time_when: kwarg:                 {'uci': 100.47}
IsotopeQuantity.time_when: ref_atoms: 2.111973074488178e+23
IsotopeQuantity.time_when: target:    2.111973074488178e+23
IsotopeQuantity.time_when: dt:        -0.0
isotope_qty_test: test_isotopequantity_time_when: iq.time_when(**kwarg): 2021-01-04 17:32:17.773989
isotope_qty_test: test_isotopequantity_time_when: ref_date:              2021-01-04 17:32:17.773989
IsotopeQuantity.time_when: ref_atoms: 2.111973074488178e+23
IsotopeQuantity.time_when: target:    2.111973074488178e+23
IsotopeQuantity.time_when: dt:        -0.0

(In the last line, now dt is effectively 0.0, so the test passes.)

Here's a minimum working example of the failure:

import datetime
import numpy as np
import becquerel as bq

# original instantiation
iq0 = bq.IsotopeQuantity('K-40', date=datetime.datetime.now(), uci=10.047)
# multiply activity by 10
iq1 = bq.IsotopeQuantity('K-40', date=datetime.datetime.now(), uci=100.47)

print(iq0.uci_at(iq0.ref_date))
print(iq1.uci_at(iq1.ref_date))
assert iq0.uci_at(iq0.ref_date) == iq1.uci_at(iq1.ref_date)

10.047000000000002
100.47
AssertionError

I did some more testing and found the place were nearby uci values result in slightly different results, presumably due to rounding errors:

import datetime
import numpy as np
import becquerel as bq

iq0 = bq.IsotopeQuantity('K-40', date=datetime.datetime.now(), uci=11.050504791)
iq1 = bq.IsotopeQuantity('K-40', date=datetime.datetime.now(), uci=11.050504790)

print(iq0.uci_at(iq0.ref_date))
print(iq1.uci_at(iq1.ref_date))

11.050504791
11.050504790000002

Any thoughts on what this could be or how to address it?

[cosama]

We just need to use np.close in the tests. There is nothing else really we can do about that. Weird though it didn't show up earlier.

[cosama]

It is an issue with a package, not sure which one though. I run all test, they passed, then I updated all the required packages and it failed.

[markbandstra]

Seems like if we use numpy.isclose we will need to extract the seconds and set a tolerance. Here is an actual minimum working example:

import datetime
import numpy as np
import becquerel as bq

iso = 'K-40'
dt0 = datetime.datetime.now()

iq1 = bq.IsotopeQuantity(iso, date=dt0, uci=10.047)
iq2 = bq.IsotopeQuantity(iso, date=dt0, uci=100.47)

dt1 = iq1.time_when(uci=iq1.uci_at(dt0))
dt2 = iq2.time_when(uci=iq2.uci_at(dt0))

print(dt0)
print(dt1)
print(dt2)

diff1 = dt1 - dt0
diff2 = dt2 - dt0

print(diff1.total_seconds())
print(diff2.total_seconds())

print(np.isclose(diff1.total_seconds(), 0, atol=15.))
print(np.isclose(diff2.total_seconds(), 0, atol=15.))

2021-01-04 20:09:11.317735
2021-01-04 20:08:58.702641
2021-01-04 20:09:11.317735
-12.615094
0.0
True
True

[cosama]

@markbandstra I am looking into this right now as well. Quite confusing to be honest, still a bit puzzled where it arises. I think it is from

iq = bq.IsotopeQuantity("K-40", date=datetime.datetime.now(), uci=10.047)
a = 10.047
a, a*bq.tools.isotope_qty.UCI_TO_BQ/iq.decay_const*iq.decay_const/bq.tools.isotope_qty.UCI_TO_BQ

results in:

(10.047, 10.047000000000002)

but not sure why this didn't show up earlier.

[cosama]

One way to fix this would be to calculate all the possible quantities during initialization and enforce the quantity used for initialization to be unchanged. The __init__ would have to look something like this (not tested and not sure if there is a more elegant way):

    def __init__(self, isotope, date=None, stability=1e18, **kwargs):

        self._init_isotope(isotope, stability)
        self._init_date(date)
        self._ref_quantities = {}

        # list that defines the order of insertion
        order = ["atoms", "bq", "uci", "q"]

        # dictionary with functions that define how to calculate all quantities in a circular manner
        conversions = dict(
            atoms=lambda : self._ref_quantities["g"] * N_AV / self.isotope.A,
            bq=lambda : self._ref_quantities["atoms"] * self.decay_const,
            uci=lambda : self._ref_quantities["bq"] / UCI_TO_BQ,
            g=lambda : self._ref_quantities["uci"] * UCI_TO_BQ /
                       self.decay_const / N_AV * self.isotope.A
        )

        # this will rotate the order so that the provided kwarg is at [0]
        while order[0] not in kwargs:
            order.append(order.pop(0))

        # the kwarg provided will be the base, all other quantities will be calculated from it
        first = order.pop(0)
        self._ref_quantities[first] = kwargs[first]
        for i in order:
            self._ref_quantities[i] = conversions[i]()

Then we use self._ref_quantities[quantity], everywhere.

We could also just define all the permutations:

if "atoms" in kwargs:
    order = ["atoms", "bq", "uci", "q"]
elif "bq" in kwargs:
    order = ["bq", "uci", "q", "atoms"]
elif "uci" in kwargs:
    order = ["uci", "q", "atoms", "bq"]
elif "q" in kwargs:
    order = ["q", "atoms", "bq", "uci"]
else:
    raise ValueError("Unknown input")

[cosama]

Alternatively we could also define the relative tolerance with respect to the half-live or the inverse of the decay constant:

np.isclose(diff1.total_seconds()*iq.decay_const/np.log(2), 0, atol=1e-9.)

[jccurtis]

Regarding permutations above, some DRY could be done with just rolling the list: https://stackoverflow.com/questions/2150108/efficient-way-to-rotate-a-list-in-python

[cosama]

@jccurtis Yeah, my first example above was rolling the list:

        # this will rotate the order so that the provided kwarg is at [0]
        while order[0] not in kwargs:
            order.append(order.pop(0))

Just wasn't sure what would be more clear.

If you guys think this makes sense I am happy to give it a shot.