keflavich / pyradex

Python interface to RADEX
BSD 3-Clause "New" or "Revised" License
19 stars 12 forks source link

Arrays are not almost equal to 5 decimals #31

Closed yangcht closed 4 years ago

yangcht commented 4 years ago

Hi there! I am trying to run pyradex under Python 3.8.3 (with numpy version 1.19.1). However, I encountered an error regarding the decimals. It seems like that there's something related numpy version?

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
~/emcee/emcee_radex.py in <module>
    184 # Note that N_CO is totaly degenerated with deltav, so in
    185 # principle we should fit for N_CO/deltav
--> 186 R = pyradex.Radex(species='co', datapath="radex_moldata",
    187                   density={'oH2':fortho*10.**10.0,'pH2':(1-fortho)*10.**10.0},
    188                   column=10.0**6.0,

~/opt/anaconda3/envs/py38/lib/python3.8/site-packages/pyradex/core.py in __init__(self, collider_densities, density, total_density, temperature, species, column, column_per_bin, tbackground, deltav, abundance, datapath, escapeProbGeom, outfile, logfile, debug, mu, source_area)
    356             self.column_per_bin = column_per_bin
    357         elif column is not None:
--> 358             self.column_per_bin = column
    359         else:
    360             self._locked_parameter = 'density'

~/opt/anaconda3/envs/py38/lib/python3.8/site-packages/pyradex/core.py in column_per_bin(self, col)
    811 
    812             invab = (self.total_density / (self.column / self.length)).decompose().value
--> 813             np.testing.assert_almost_equal(invab, 1/self.abundance, decimal=5)
    814 
    815     @property

    [... skipping hidden 1 frame]

AssertionError: 
Arrays are not almost equal to 5 decimals
 ACTUAL: 3.0856775814671917e+22
 DESIRED: 3.0856775814671912e+22
keflavich commented 4 years ago

huh, that's annoying... those are equal to 5 decimal places... you can try hacking the decimal= to be a smaller value and see if that skips this issue.

yangcht commented 4 years ago

Changing to decimal=1 seems not working. Python 3.7 is working fine though.

keflavich commented 4 years ago

thanks. This is weird. I added this assertion step to make sure we don't end up in a non-self-consistent parameter state, but I don't understand why it's going wrong.

SmirnGreg commented 4 years ago

Hi @keflavich,

I think this is because numpy.testing.assert_almost_equal is abs(desired-actual) < 1.5 * 10**(-decimal) being abolute error not relative error. https://numpy.org/doc/stable/reference/generated/numpy.testing.assert_almost_equal.html with both atol and rtol.

I think you wanted to use https://numpy.org/doc/stable/reference/generated/numpy.allclose.html instead. AssertionErorr should not in principle be raised outside of [unit]tests. I think https://docs.python.org/3/library/exceptions.html?highlight=assertionerror#RuntimeError RuntimeError or UserError would suite the case better. It also might be a good idea to define your own error for this, something like

class PyRadexInconsistencyError(UserError, RuntimeError):
    pass

if you want to explicitly raise an error here.

keflavich commented 4 years ago

You're right, we want a relative tolerance here because these numbers are so huge.

@SmirnGreg you've solved the problem - could you submit it as a PR?