Fortran-exponents - #3655

Something like this would be nice to have for #3430, since FITS ASCII tables use FORTRAN-inspired string formatting.

If it's already using the c parser that should be fairly straightforward; if it isn't it would certainly benefit from the speedup as well. I'll have a look when I have time; I have some python-level code to check for the various exponent styles as well.

I think the automatic mode (where all of the reasonable options, EDQedq, are supported) may make some sense because it would lessen the need for the user option. I think it's worth implementing and benchmarking to see how the performance compares.

Any sense of how common is the triple-digit exponent thing is? I've never come across it (not that that means anything).

The triple-digit exponent could happen if you computed for example blackbody fluxes in float64 as a function of frequency, since some values would be really small.

Or chemical equilibrium fractions of compounds that don't really exist under the given conditions...
I actually have lots of data where those appear frequently, but they are small enough to be handled at the Python level.
The performance hit certainly would be tolerable, but kind of defeating the original purpose of the fast converter... Don't think it can be easily avoided even when not using the automatic option since the extra conditional(s) are in there in a rather time-critical loop.

I see -- so the triple-digit simplification happens automatically? That's what I was getting at -- whether it was an automatic behavior or something a Fortran programmer would opt-in to. Sounds like the former.

As for the potential performance hit: Let's do the most straightforward implementation, and benchmark it to get a sense of what the hit is. If we can't deal with the hit, a solution would be to have multiple implementations of xstrod, one like it is now, and one with automatic behavior + triple-digit support. When the user specifies "automatic mode" we move a function pointer to point to the latter. Therefore the performance hit only goes into effect when automatic mode is on. But we may not need to do all that -- it may be that other parts of the reader dominate over xstrtod anyway.

Yes, the programmer can avoid it with an explicit format statement like D20.12D3, but by default the letter is just silently dropped.

OK, I've added the auto-detection with a bit of extra optimisation (thought the compiler would actually be able to abstract from the nested ifs...); for a real-world test case with some 2.3 Mio. lines, 26 columns 21 of which are read, maybe 20-40% of the original speed-up of use_fast_converter is lost (i.e. 7.02 s vs. 7.4 s have now increased to 7.12-7.2 s), but that's close to the timing uncertainty IMO. In fact the fastest runs have now occurred with the auto-detection turned on.

I guess the way to add it to the standard parser would be to pre-process the entire string before passing it to strtod, similar to what's suggested here:

fgets(buf, sizeof buf, str);
for (p = buf; p = strchr(p, 'D'); *p = 'E') {} 
for (p = buf; p = strchr(p, 'Q'); *p = 'E') {}

(not tested nor seriously reviewed due to lack of C expertise!)

@dhomeier - thanks, this is looking pretty good from a top-level look! I'll try to go through in detail in the next couple days. Perhaps @mdmueller can take a look as well.

I think just having this in the fast parser is reasonable. Since the pure Python parser really just relies on numpy it would be a significant change to add this capability to the slow readers.

@taldcroft - I rather had the not-so-fast C parser in mind anyway; for the Python parser it would mainly amount to providing a convenience function to specify converters like
(n, lambda s: float(s.replace(b'D',b'e'))) for n in names
but then, if users specified some converters of their own, the trouble might come in...

Reviewing the build failures reported on AppVeyor and Travis, the former fails with the parallel version of test_fortran_reader, not surprising of course since the parallel reader is not supported under Windows. Just wondering if I had missed something in the setup of the parametrize decorator? The only test that specifically checks for the Windows case is _read(), yet none of the other tests seem to have reported failures either.

Travis build fails due to a Sphinx warning:
/home/travis/build/astropy/astropy/build/lib.linux-x86_64-2.7/astropy/io/ascii/ui.py:docstring of astropy.io.ascii.read:98: WARNING: Block quote ends without a blank line; unexpected unindent.
I did not remove any blank lines, don't know though if the 2nd-level indentation is illegal in a docstring.
It looks OK in help(ascii.read), and seemed more structured to me this way...

For the user API (the key in fast_reader dict), what about calling this exponent_style with following possible values:

• E (default), D, Q (case-insensitive): used as-is
• fortran (case-insensitive): same as the current A which allows any fortran-compatible style.

So not really any change to the C code, just a more obvious front-end interface.

About the Sphinx warnings, you'll just need to play around until you get something that builds locally without warnings. I'm often surprised by what Sphinx doesn't like. I do like your idea of making it visually parseable.

This change will also benefit from updating a couple of places in the narrative docs so people actually know about. First is in the page about the fast reader:
http://astropy.readthedocs.org/en/latest/io/ascii/fast_ascii_io.html

Next, I would suggest a new section called "Fortran-style exponents" just after the Converters section and at the same section header level as Converters:
http://astropy.readthedocs.org/en/latest/io/ascii/read.html#converters

Humm, running build_sphinx on my Mac returns a couple of warnings
build succeeded, 2028 warnings.
I'll try to watch out for those related to the new docstrings; I could see now that the html doc for astropy.io.ascii.ui.read.html did not render the 2nd indentation level right and apparently needed the blank line before that (seems to be the Block quote ends without a blank line; unexpected unindent. warning from Sphinx). Now the interactive help has an extra line, but probably can live with that.

I'll upload those doc fixes after changing the API to your suggestions. Any objections against silently accepting any substring of Fortran for the automatic mode?

2028 warnings? That's not normal...

Don't know if the sphinx installation is broken - it's not terribly old, 1.2.3, but the python 3.4 variant does not work at all... Most of the warnings seem to be of this kind:

/Users/derek/git/astropy/docs/api/astropy.wcs.WCS.rst:40: WARNING: py:obj reference target not found: int
/Users/derek/git/astropy/docs/api/astropy.wcs.WCS.rst:76: WARNING: py:obj reference target not found: None
/Users/derek/git/astropy/docs/api/astropy.wcs.WCS.rst:82: WARNING: py:obj reference target not found: True
/Users/derek/git/astropy/docs/api/astropy.wcs.WCS.rst:89: WARNING: py:obj reference target not found: True
/Users/derek/git/astropy/docs/api/astropy.wcs.WCS.rst:115: WARNING: py:obj reference target not found: exceptions.ValueError

That seems like errors that would happen if there was an issue with the internet connection. If you re-run the build, do you still get that many warnings?

1633 - the connection should be fast enough, not sure if the firewall might cause any problems.

Near the beginning of the doc build it should say something like:

loading intersphinx inventory from http://docs.python.org/objects.inv...
loading intersphinx inventory from http://matplotlib.org/objects.inv...
loading intersphinx inventory from http://docs.scipy.org/doc/scipy/reference/objects.inv...
loading intersphinx inventory from http://docs.h5py.org/en/latest/objects.inv...
loading intersphinx inventory from http://docs.scipy.org/doc/numpy/objects.inv...

These are all http requests. Are there any messages about them failing? You may try removing the docs/_build directory before the rebuild in case any of those files are incomplete from an earlier failed request (that's just a shot in the dark, of course).

@mdboom - thanks, you're right; they are failing:

loading intersphinx inventory from http://docs.python.org/objects.inv...
WARNING: intersphinx inventory 'http://docs.python.org/objects.inv' not fetchable due to <class 'urllib2.URLError'>: <urlopen error [SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed (_ssl.c:581)>

This is a "feature" of the latest openssl (at least the one in use on Mac OS since ~January); no idea if it is correct in not trusting certain certificate servers...
Not so easy to work around this, but I'll try to test the doc build on a Linux machine for the next commit.

Just pinging - is there anything left to be done at this point, suggestions for labels (I missed to mark it as ENH, it's obviously io.ascii, no idea about how release-affecting this is)?
I guess some feedback/discussion on the documentation could be useful.

Finally realised that the segfaults on Travis occur on the same cases as known from test_data_out_of_range, therefore marked it xfail know.

I note that there are still three different outcomes of too large exponents:

>>> text = 'a b c d\n10.1E+199 3.14e+313 2048e+306 0.6E-414'
>>> ascii.read(text, format='basic', guess=False, fast_reader=False)
<Table length=1>
    a        b       c       d   
 float64  float64 float64 float64
--------- ------- ------- -------
1.01e+200     inf     inf     0.0
>>> ascii.read(text, format='basic', guess=False, fast_reader={'exponent_style': 'e'})
<Table length=1>
    a         b         c        d    
 float64     str9      str9     str8  
--------- --------- --------- --------
1.01e+200 3.14e+313 2048e+306 0.6E-414
>>> ascii.read(text, format='basic', guess=False, fast_reader={'exponent_style': 'fortran'})
<Table length=1>
    a         b         c        d   
 float64     str9      str9   float64
--------- --------- --------- -------
1.01e+200 3.14e+313 2048e+306     0.0

Should one ultimately try to make the C reader return inf for too large exponents as well for consistency with the Python parser or even allow float128 as return values (I imagine this should be feasible by modifying xstrtod to return exponent and mantissa as separate numbers)?

One more thought, whatever the opinion on reading in extended-precision floats, the above listing shows one difference remaining between the existing fast_reader and the results using exponent_style extension: input with very small exponents is returned as float64 rounded to zero rather than as string type.
Arguably the C parser's way of returning string arrays is, while not very pretty, for the time being still the most useful behaviour, since it at least allows one to repair the data pedestrian style after read-in:

newarr = np.empty(1, dtype=[('a', 'f8'), ('b', 'f16'), ('c', 'f16'), ('d', 'f16')])
newarr ['a'] = readarr['a']
newarr['b'][0] = np.float128(str(readarr['b'][0]))
...

So is there a general feeling it would be desirable to ensure consistency with the existing fast_reader functionality for the case of column 'd' above?

About the too-large exponents, for this PR I would suggest maintaining consistency with Python and returning inf or 0. Not all platforms support float128 and in any case it would be better to have this discussion in a separate context.

Overall it is difficult to look at the tests and get a good understanding of what is being tested, in particular with regard to what parses and what does not. What would help is defining this in terms of a uniform test matrix like:

text_value expected exponent_style
...

You could actually do the same in code by just making a big list of tuples with those same three value expected exponent_style. The point is that the coverage of various format possibilities (correct and incorrect) is crucial here. It will be much more obvious when organized in a single table.

OK, looks like the Travis errors are checked at last. The fast readers now return inf or 0 as float64 for out_of_range values consistent with the Python parser, only the detailed results for some borderline cases may differ:

>>> asc.read('a\n 2.5e-324', format='basic', guess=False, fast_reader=False)['a']
<Column name='a' dtype='float64' length=1>
4.94065645841e-324
>>> asc.read('a\n 2.5e-324', format='basic', guess=False, fast_reader={'use_fast_converter': True})['a']
<Column name='a' dtype='float64' length=1>
4.94065645841e-324
>>> asc.read('a\n 2.5e-324', format='basic', guess=False, fast_reader={'use_fast_converter': False})['a']
<Column name='a' dtype='str8' length=1>
2.5e-324
>>> asc.read('a\n 2.45e-324', format='basic', guess=False, fast_reader={'use_fast_converter': False})['a']
<Column name='a' dtype='float64' length=1>
0.0
>>> asc.read('a\n 2.5e-308', format='basic', guess=False, fast_reader={'use_fast_converter': False})['a']
<Column name='a' dtype='float64' length=1>
2.5e-308

It seems for numbers larger than the smallest representable float64, 2.5e-324, but smaller than 1/ np.finfo(np.float).max, the builtin strtod (at least for the Darwin gcc here) raises errno == ERANGE without returning a value of HUGE_VAL or 0. May need some more investigation, but I am beginning to wonder if we should trigger an OVERFLOW_ERROR in cparser at all, or just return inf or 0 in all those cases.
Will look into the tests structure later; the new fast_reader behaviour changed the outcome of test_convert_overflow as well, returning extremely long integer inputs as valid float64 or inf now.

Ah, I guess strtod raises the ERANGE flag because those numbers cannot be represented with the specified resolution of 1.e-15 (as is obvious above). But if the policy is not to fall back to str as long as possible, any inaccurate float would still be better than an even less accurate 0, right?

I obviously haven't dug into this at the level of detail you have been doing, but I think the guiding principle here should be predictability and consistency. I.e. is there a simple way to predict what will happen for different cases, and is it consistent with behavior of Python and especially Numpy. I'm not super worried about the extreme edge cases like 2.5e-324.

I actually can't think of any case where the ERANGE flag should raise a ValueError, as all parsers now consistently return those numbers as float64. I think it would make sense to raise a Warning however, similar to the OverflowError warning issued by the core.py converter when the range of IntType is exceeded.

BTW long integer conversion is still handled inconsistently between the Python parser, which switches right away to str when the range is exceeded, and the fast readers that try the float converter first. But in PR #2234 the agreement was eventually to fall back to str from int, so maybe that should be the common ground.
In any case, if the sequence of converters should remain int -> float -> str, it would be very desirable to have an option to override this to be able to read a column in string format at all. The try_float option in cparser unfortunately seems only settable from the file header.

To summarise, I see 3 options for handling input values extending the respective numerical range:

1. int overflows are returned as str by the standard reader and promoted to float (which may include inf) by the fast readers; float inputs remain float64.
2. int overflows are returned as str by the standard and fast readers; float inputs remain float64.
3. int overflows are promoted to float (which may include inf) by all readers; float inputs remain float64.

*1 is the present situation, *2 would be the consistent solution keeping the spirit of #2234, so I'd tend to this option. In principle one could think of a 4th option always returning float overflows as str as well, but I don't see a straightforward way to identify overflows in the convert_numpy return, to begin with. I'd rather add some comments on explicitly specifying converters to return an S dtype to the documentation, if users want to preserve the full information.

I would vote for option 2 as well.

OK, it's been some time again before I got back to this; I have implemented option 2 from above
int overflows are returned as str by the standard and fast readers; float inputs remain float64.
plus a Warning for the first overflow in either int or float, as well as possible. Some notes:

1. the actual fallback in mixed integer + float columns depends on order: if an int overflow shows up first, the column will be cast to str, if a float appears before that, it will be float. This same behaviour is already produced by the Python parser, and I don't see an easy way to fix it without reprocessing the entire column in question. I've put in an xfail test however to keep it on the radar.
2. The parser now ignores leading zeros in the input; working on that turned up that the integer parser was called with auto-detection flag set: strtol(str, &tmp, 0), which caused it to interpret leading zeros as indicating an octal number. Depending on the digits this would either return the octal or raise a CONVERSION_ERROR. Given that it is quite unlikely anyone has used this feature, and it is also again inconsistent with the Python parser, I've changed the call to strtol(str, &tmp, 10).

I've changed the invalid exponents tests according to taldcroft's suggestions, but left the multi-row input for some of the other tests, since columns with mixed types are part of the test case (see above).

However I feel rather lost with the new CI system - I cannot see any details at all on what's happening on AppVeyor, and CircleCI only shows this quite unhelpful message:

ImportError: No module named jinja2

Also I don't find any way to investigate the merge conflicts reported for astropy/io/ascii/tests/test_c_reader.py...

Could not find any way to work around the conflicts here, so I've merged the branch locally into the updated master and submitted this as a new PR #5552.