Repository: alchemistry/alchemlyb

Issue:
Gromacs parser cannot handle incomplete lines with -
For an incomplete line like

`1 -100 -10 -`

The last `-` will be giving problems as it cannot be cast as a float during 
https://github.com/xiki-tempula/alchemlyb/blob/093ad87eca736f724c9a4b4d6843043ea9861f69/src/alchemlyb/parsing/gmx.py#L308

where dtype=np.float64 would enforce it to be a float and `-` will give

```
File "/Users/XXX/opt/miniconda3/envs/mda_dev/lib/python3.7/site-packages/pandas/io/parsers.py", line 610, in read_csv
    return _read(filepath_or_buffer, kwds)
  File "/Users/XXX/opt/miniconda3/envs/mda_dev/lib/python3.7/site-packages/pandas/io/parsers.py", line 468, in _read
    return parser.read(nrows)
  File "/Users/XXX/opt/miniconda3/envs/mda_dev/lib/python3.7/site-packages/pandas/io/parsers.py", line 1057, in read
    index, columns, col_dict = self._engine.read(nrows)
  File "/Users/XXX/opt/miniconda3/envs/mda_dev/lib/python3.7/site-packages/pandas/io/parsers.py", line 2061, in read
    data = self._reader.read(nrows)
  File "pandas/_libs/parsers.pyx", line 756, in pandas._libs.parsers.TextReader.read
  File "pandas/_libs/parsers.pyx", line 771, in pandas._libs.parsers.TextReader._read_low_memory
  File "pandas/_libs/parsers.pyx", line 850, in pandas._libs.parsers.TextReader._read_rows
  File "pandas/_libs/parsers.pyx", line 982, in pandas._libs.parsers.TextReader._convert_column_data
  File "pandas/_libs/parsers.pyx", line 1056, in pandas._libs.parsers.TextReader._convert_tokens
ValueError: could not convert string to float: '-'
```


Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/alchemistry/alchemlyb/pull/183?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry) Report
> Merging [#183](https://codecov.io/gh/alchemistry/alchemlyb/pull/183?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry) (278557a) into [master](https://codecov.io/gh/alchemistry/alchemlyb/commit/0f80259a173c7ebad22dc43f076e450f83bb3d18?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry) (0f80259) will **decrease** coverage by `0.15%`.
> The diff coverage is `100.00%`.

```diff
@@            Coverage Diff             @@
##           master     #183      +/-   ##
==========================================
- Coverage   98.16%   98.00%   -0.16%     
==========================================
  Files          22       22              
  Lines        1305     1303       -2     
  Branches      286      279       -7     
==========================================
- Hits         1281     1277       -4     
  Misses          5        5              
- Partials       19       21       +2     
```


| [Impacted Files](https://codecov.io/gh/alchemistry/alchemlyb/pull/183?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry) | Coverage Δ | |
|---|---|---|
| [src/alchemlyb/parsing/\_\_init\_\_.py](https://codecov.io/gh/alchemistry/alchemlyb/pull/183/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry#diff-c3JjL2FsY2hlbWx5Yi9wYXJzaW5nL19faW5pdF9fLnB5) | `100.00% <100.00%> (ø)` | |
| [src/alchemlyb/parsing/gmx.py](https://codecov.io/gh/alchemistry/alchemlyb/pull/183/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry#diff-c3JjL2FsY2hlbWx5Yi9wYXJzaW5nL2dteC5weQ==) | `97.57% <100.00%> (-0.55%)` | :arrow_down: |
| [src/alchemlyb/parsing/gomc.py](https://codecov.io/gh/alchemistry/alchemlyb/pull/183/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry#diff-c3JjL2FsY2hlbWx5Yi9wYXJzaW5nL2dvbWMucHk=) | `92.07% <0.00%> (-1.13%)` | :arrow_down: |
| [src/alchemlyb/estimators/mbar\_.py](https://codecov.io/gh/alchemistry/alchemlyb/pull/183/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry#diff-c3JjL2FsY2hlbWx5Yi9lc3RpbWF0b3JzL21iYXJfLnB5) | `98.21% <0.00%> (-0.04%)` | :arrow_down: |
| [src/alchemlyb/parsing/amber.py](https://codecov.io/gh/alchemistry/alchemlyb/pull/183/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry#diff-c3JjL2FsY2hlbWx5Yi9wYXJzaW5nL2FtYmVyLnB5) | `97.46% <0.00%> (-0.03%)` | :arrow_down: |
| [src/alchemlyb/parsing/namd.py](https://codecov.io/gh/alchemistry/alchemlyb/pull/183/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry#diff-c3JjL2FsY2hlbWx5Yi9wYXJzaW5nL25hbWQucHk=) | `99.16% <0.00%> (-0.01%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/alchemistry/alchemlyb/pull/183?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/alchemistry/alchemlyb/pull/183?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry). Last update [0f80259...278557a](https://codecov.io/gh/alchemistry/alchemlyb/pull/183?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alchemistry).

xiki-tempula: @orbeckst Thanks for the review. I made the changes. I wonder how should I bump the version? sorry, I'm not too familiar with the versioner.
orbeckst: Just change the entry in the CHANGELOG.

The *actual* version change happens on release anyway. 
xiki-tempula: @orbeckst Thanks for the review. I have changed the default to True. And add versionchange
It seems that filtering is actually faster but this might just be random fluctuation.

```
dataset = alchemtest.gmx.load_benzene()
time_list = []
for i in range(100):
    start = time.time()
    alchemlyb.concat([gmx.extract_u_nk(filename, T=300, filter=True) for filename in dataset['data']['Coulomb']])
    time_list.append(time.time()-start)
print(np.mean(time_list))
```
gives 0.14853424549102784
Changing the filter=False gives 0.152081116437912



Patch:
diff --git a/CHANGES b/CHANGES
index 824dcc5..94d700f 100644
--- a/CHANGES
+++ b/CHANGES
@@ -13,8 +13,19 @@ The rules for this file:
   * release numbers follow "Semantic Versioning" https://semver.org
 
 ------------------------------------------------------------------------------
-??/??/2022 xiki-tempula
-  * 0.6.1
+
+*/*/2022 xiki-tempula
+
+  * 0.7.0
+
+Changes
+  - gmx parser now defaults to dropping NaN and corrupted lines (filter=True) 
+    (#171, PR #183)
+
+Enhancements
+  - Add filter function to gmx.extract to make it more robust (PR #183): can filter 
+    incomplete/corrupted lines (#126, #171) with filter=True.
+    
 
 Fixes
   - added AutoMBAR to convergence analysis (#189)
diff --git a/src/alchemlyb/parsing/__init__.py b/src/alchemlyb/parsing/__init__.py
index eae423c..60b0f34 100644
--- a/src/alchemlyb/parsing/__init__.py
+++ b/src/alchemlyb/parsing/__init__.py
@@ -7,8 +7,8 @@ def _init_attrs(func):
     and the energy unit is initiated as dataframe.attrs['energy_unit'] = 'kT'.
     '''
     @wraps(func)
-    def wrapper(outfile, T):
-        dataframe = func(outfile, T)
+    def wrapper(outfile, T, *args, **kwargs):
+        dataframe = func(outfile, T, *args, **kwargs)
         if dataframe is not None:
             dataframe.attrs['temperature'] = T
             dataframe.attrs['energy_unit'] = 'kT'
diff --git a/src/alchemlyb/parsing/gmx.py b/src/alchemlyb/parsing/gmx.py
index 3e3678b..ab56e5d 100644
--- a/src/alchemlyb/parsing/gmx.py
+++ b/src/alchemlyb/parsing/gmx.py
@@ -11,7 +11,7 @@ from ..postprocessors.units import R_kJmol
 k_b = R_kJmol
 
 @_init_attrs
-def extract_u_nk(xvg, T):
+def extract_u_nk(xvg, T, filter=True):
     r"""Return reduced potentials `u_nk` from a Hamiltonian differences XVG file.
 
     Parameters
@@ -20,6 +20,10 @@ def extract_u_nk(xvg, T):
         Path to XVG file to extract data from.
     T : float
         Temperature in Kelvin the simulations sampled.
+    filter : bool
+        Filter out the lines that cannot be parsed.
+        Such as rows with incorrect number of Columns and incorrectly
+        formatted numbers (e.g. 123.45.67, nan or -).
 
     Returns
     -------
@@ -49,6 +53,10 @@ def extract_u_nk(xvg, T):
         This leads to slightly different results for GROMACS input compared to
         previous versions of alchemlyb.
 
+    .. versionchanged:: 0.7.0
+        The keyword filter is implemented to ignore the line that cannot be
+        parsed and is turned on by default.
+
     """
 
     h_col_match = r"\xD\f{}H \xl\f{}"
@@ -59,7 +67,7 @@ def extract_u_nk(xvg, T):
     state, lambdas, statevec = _extract_state(xvg)
 
     # extract a DataFrame from XVG data
-    df = _extract_dataframe(xvg)
+    df = _extract_dataframe(xvg, filter=filter)
 
     times = df[df.columns[0]]
 
@@ -128,7 +136,7 @@ def extract_u_nk(xvg, T):
     return u_k
 
 @_init_attrs
-def extract_dHdl(xvg, T):
+def extract_dHdl(xvg, T, filter=True):
     r"""Return gradients `dH/dl` from a Hamiltonian differences XVG file.
 
     Parameters
@@ -137,6 +145,10 @@ def extract_dHdl(xvg, T):
         Path to XVG file to extract data from.
     T : float
         Temperature in Kelvin the simulations sampled.
+    filter : bool
+        Filter out the lines that cannot be parsed.
+        Such as rows with incorrect number of Columns and incorrectly
+        formatted numbers (e.g. 123.45.67, nan or -).
 
     Returns
     -------
@@ -165,6 +177,10 @@ def extract_dHdl(xvg, T):
         This leads to slightly different results for GROMACS input compared to
         previous versions of alchemlyb.
 
+    .. versionchanged:: 0.7.0
+        The keyword filter is implemented to ignore the line that cannot be
+        parsed and is turned on by default.
+
     """
     beta = 1/(k_b * T)
 
@@ -172,7 +188,7 @@ def extract_dHdl(xvg, T):
     state, lambdas, statevec = _extract_state(xvg, headers)
 
     # extract a DataFrame from XVG data
-    df = _extract_dataframe(xvg, headers)
+    df = _extract_dataframe(xvg, headers, filter=filter)
 
     times = df[df.columns[0]]
 
@@ -276,7 +292,7 @@ def _extract_legend(xvg):
     return state_legend
 
 
-def _extract_dataframe(xvg, headers=None):
+def _extract_dataframe(xvg, headers=None, filter=filter):
     """Extract a DataFrame from XVG data using Pandas `read_csv()`.
 
     pd.read_csv() shows the same behavior building pandas Dataframe with better
@@ -289,6 +305,10 @@ def _extract_dataframe(xvg, headers=None):
     headers: dict
        headers dictionary to search header information, reduced I/O by
        reusing if it is already parsed. Direct access by key name
+    filter : bool
+        Filter out the lines that cannot be parsed.
+        Such as rows with incorrect number of Columns and incorrectly
+        formatted numbers (e.g. 123.45.67, nan or -).
 
     """
     if headers is None:
@@ -304,9 +324,23 @@ def _extract_dataframe(xvg, headers=None):
             for i, col, in enumerate(cols)]
 
     header_cnt = len(headers['_raw_lines'])
-    df = pd.read_csv(xvg, sep=r"\s+", header=None, skiprows=header_cnt,
-            na_filter=True, memory_map=True, names=cols, dtype=np.float64,
-            float_precision='high')
+    if not filter:
+        # assumes clean files
+        df = pd.read_csv(xvg, sep=r"\s+", header=None, skiprows=header_cnt,
+                         na_filter=True, memory_map=True, names=cols,
+                         dtype=np.float64,
+                         float_precision='high')
+    else:
+        df = pd.read_csv(xvg, sep=r"\s+", header=None, skiprows=header_cnt,
+                memory_map=True, on_bad_lines='skip')
+        # If names=cols is passed to read_csv, rows with more than the
+        # designated columns will be truncated and used instead of discarded.
+        df.rename(columns={i: name for i, name in enumerate(cols)}, inplace=True)
+        # If dtype=np.float64 and float_precision='high' are passed to read_csv,
+        # 12.345.56 and - cannot be read.
+        df = df.apply(pd.to_numeric, errors='coerce')
+        # drop duplicate
+        df.dropna(inplace=True)
 
     # drop duplicated columns (see PR #86 https://github.com/alchemistry/alchemlyb/pull/86/)
     df = df[df.columns[~df.columns.str.endswith("[duplicated]")]]


Test patch:
diff --git a/src/alchemlyb/tests/parsing/test_gmx.py b/src/alchemlyb/tests/parsing/test_gmx.py
index 063aafe..0c3a087 100644
--- a/src/alchemlyb/tests/parsing/test_gmx.py
+++ b/src/alchemlyb/tests/parsing/test_gmx.py
@@ -2,6 +2,9 @@
 
 """
 
+import bz2
+import pytest
+
 from alchemlyb.parsing.gmx import extract_dHdl, extract_u_nk
 from alchemtest.gmx import load_benzene
 from alchemtest.gmx import load_expanded_ensemble_case_1, load_expanded_ensemble_case_2, load_expanded_ensemble_case_3
@@ -48,7 +51,7 @@ def test_u_nk_case1():
 
     for leg in dataset['data']:
         for filename in dataset['data'][leg]:
-            u_nk = extract_u_nk(filename, T=300)
+            u_nk = extract_u_nk(filename, T=300, filter=False)
 
             assert u_nk.index.names == ['time', 'fep-lambda', 'coul-lambda', 'vdw-lambda', 'restraint-lambda']
 
@@ -62,7 +65,7 @@ def test_dHdl_case1():
 
     for leg in dataset['data']:
         for filename in dataset['data'][leg]:
-            dHdl = extract_dHdl(filename, T=300)
+            dHdl = extract_dHdl(filename, T=300, filter=False)
 
             assert dHdl.index.names == ['time', 'fep-lambda', 'coul-lambda', 'vdw-lambda', 'restraint-lambda']
             assert dHdl.shape == (50001, 4)
@@ -75,7 +78,7 @@ def test_u_nk_case2():
 
     for leg in dataset['data']:
         for filename in dataset['data'][leg]:
-            u_nk = extract_u_nk(filename, T=300)
+            u_nk = extract_u_nk(filename, T=300, filter=False)
 
             assert u_nk.index.names == ['time', 'fep-lambda', 'coul-lambda', 'vdw-lambda', 'restraint-lambda']
 
@@ -89,7 +92,7 @@ def test_u_nk_case3():
 
     for leg in dataset['data']:
         for filename in dataset['data'][leg]:
-            u_nk = extract_u_nk(filename, T=300)
+            u_nk = extract_u_nk(filename, T=300, filter=False)
 
             assert u_nk.index.names == ['time', 'fep-lambda', 'coul-lambda', 'vdw-lambda', 'restraint-lambda']
 
@@ -103,7 +106,7 @@ def test_dHdl_case3():
 
     for leg in dataset['data']:
         for filename in dataset['data'][leg]:
-            dHdl = extract_dHdl(filename, T=300)
+            dHdl = extract_dHdl(filename, T=300, filter=False)
 
             assert dHdl.index.names == ['time', 'fep-lambda', 'coul-lambda', 'vdw-lambda', 'restraint-lambda']
             assert dHdl.shape == (2500, 4)
@@ -195,3 +198,58 @@ def test_extract_dHdl_unit():
     dhdl = extract_dHdl(dataset['data']['Coulomb'][0], 310)
     assert dhdl.attrs['temperature'] == 310
     assert dhdl.attrs['energy_unit'] == 'kT'
+
+class TestRobustGMX():
+    '''Test dropping the row that is wrong in different way'''
+    @staticmethod
+    @pytest.fixture(scope='class')
+    def data():
+        dhdl = extract_dHdl(load_benzene()['data']['Coulomb'][0], 310)
+        with bz2.open(load_benzene()['data']['Coulomb'][0], "rt") as bz_file:
+            text = bz_file.read()
+        return text, len(dhdl)
+
+    def test_sanity(self, data, tmp_path):
+        '''Test if the test routine is working.'''
+        text, length = data
+        new_text = tmp_path / 'text.xvg'
+        new_text.write_text(text)
+        dhdl = extract_dHdl(new_text, 310)
+        assert len(dhdl) == length
+
+    def test_truncated_row(self, data, tmp_path):
+        '''Test the case where the last row has been truncated.'''
+        text, length = data
+        new_text = tmp_path / 'text.xvg'
+        new_text.write_text(text + '40010.0 27.0\n')
+        dhdl = extract_dHdl(new_text, 310, filter=True)
+        assert len(dhdl) == length
+
+    def test_truncated_number(self, data, tmp_path):
+        '''Test the case where the last row has been truncated and a - has
+        been left.'''
+        text, length = data
+        new_text = tmp_path / 'text.xvg'
+        new_text.write_text(text + '40010.0 27.0 -\n')
+        dhdl = extract_dHdl(new_text, 310, filter=True)
+        assert len(dhdl) == length
+
+    def test_weirdnumber(self, data, tmp_path):
+        '''Test the case where the last number has been appended a weird
+        number.'''
+        text, length = data
+        new_text = tmp_path / 'text.xvg'
+        # Note the 27.040010.0 which is the sum of 27.0 and 40010.0
+        new_text.write_text(text + '40010.0 27.040010.0 27.0 0.0 6.7 13.5 20.2 27.0 0.7 27.0 0.0 6.7 '
+                       '13.5 20.2 27.0 0.7\n')
+        dhdl = extract_dHdl(new_text, 310, filter=True)
+        assert len(dhdl) == length
+
+    def test_too_many_cols(self, data, tmp_path):
+        '''Test the case where the row has too many columns.'''
+        text, length = data
+        new_text = tmp_path / 'text.xvg'
+        new_text.write_text(text +
+                       '40010.0 27.0 0.0 6.7 13.5 20.2 27.0 0.7 27.0 0.0 6.7 13.5 20.2 27.0 0.7\n')
+        dhdl = extract_dHdl(new_text, 310, filter=True)
+        assert len(dhdl) == length
diff --git a/src/alchemlyb/tests/test_fep_estimators.py b/src/alchemlyb/tests/test_fep_estimators.py
index e4bf2a6..9b96b5f 100644
--- a/src/alchemlyb/tests/test_fep_estimators.py
+++ b/src/alchemlyb/tests/test_fep_estimators.py
@@ -36,7 +36,7 @@ def gmx_benzene_vdw_u_nk():
 def gmx_expanded_ensemble_case_1():
     dataset = alchemtest.gmx.load_expanded_ensemble_case_1()
 
-    u_nk = alchemlyb.concat([gmx.extract_u_nk(filename, T=300)
+    u_nk = alchemlyb.concat([gmx.extract_u_nk(filename, T=300, filter=False)
                       for filename in dataset['data']['AllStates']])
 
     return u_nk
@@ -44,7 +44,7 @@ def gmx_expanded_ensemble_case_1():
 def gmx_expanded_ensemble_case_2():
     dataset = alchemtest.gmx.load_expanded_ensemble_case_2()
 
-    u_nk = alchemlyb.concat([gmx.extract_u_nk(filename, T=300)
+    u_nk = alchemlyb.concat([gmx.extract_u_nk(filename, T=300, filter=False)
                       for filename in dataset['data']['AllStates']])
 
     return u_nk
@@ -52,7 +52,7 @@ def gmx_expanded_ensemble_case_2():
 def gmx_expanded_ensemble_case_3():
     dataset = alchemtest.gmx.load_expanded_ensemble_case_3()
 
-    u_nk = alchemlyb.concat([gmx.extract_u_nk(filename, T=300)
+    u_nk = alchemlyb.concat([gmx.extract_u_nk(filename, T=300, filter=False)
                       for filename in dataset['data']['AllStates']])
 
     return u_nk
diff --git a/src/alchemlyb/tests/test_ti_estimators.py b/src/alchemlyb/tests/test_ti_estimators.py
index fcadaed..d5d036d 100644
--- a/src/alchemlyb/tests/test_ti_estimators.py
+++ b/src/alchemlyb/tests/test_ti_estimators.py
@@ -34,7 +34,7 @@ def gmx_benzene_vdw_dHdl():
 def gmx_expanded_ensemble_case_1_dHdl():
     dataset = alchemtest.gmx.load_expanded_ensemble_case_1()
 
-    dHdl = alchemlyb.concat([gmx.extract_dHdl(filename, T=300)
+    dHdl = alchemlyb.concat([gmx.extract_dHdl(filename, T=300, filter=False)
                       for filename in dataset['data']['AllStates']])
 
     return dHdl
@@ -42,7 +42,7 @@ def gmx_expanded_ensemble_case_1_dHdl():
 def gmx_expanded_ensemble_case_2_dHdl():
     dataset = alchemtest.gmx.load_expanded_ensemble_case_2()
 
-    dHdl = alchemlyb.concat([gmx.extract_dHdl(filename, T=300)
+    dHdl = alchemlyb.concat([gmx.extract_dHdl(filename, T=300, filter=False)
                       for filename in dataset['data']['AllStates']])
 
     return dHdl
@@ -50,7 +50,7 @@ def gmx_expanded_ensemble_case_2_dHdl():
 def gmx_expanded_ensemble_case_3_dHdl():
     dataset = alchemtest.gmx.load_expanded_ensemble_case_3()
 
-    dHdl = alchemlyb.concat([gmx.extract_dHdl(filename, T=300)
+    dHdl = alchemlyb.concat([gmx.extract_dHdl(filename, T=300, filter=False)
                       for filename in dataset['data']['AllStates']])
 
     return dHdl