Pytraj nastruct needs to be updated to handle [nxyz] data

[drroe]

Cpptraj nastruct now saves the base pair axis normal (Z) vector with aspect [nxyz], see Amber-MD/cpptraj#1018.

Unfortunately, this breaks pytraj nastruct since it doesn't know how to return vector (x y z ox oy oz) data.

[hainm]

this breaks pytraj

how broken it is?

[drroe]

The nastruct test fails. I'm not sure if the entire command is broken or if its just the specific test. Everything else seems ok.

[drroe]

I'll post the error when I get home.

[hainm]

sure. thanks Dan.

[drroe]

Here is the full output:

$ pytest test_nastruct.py 
============================================ test session starts ============================================
platform linux -- Python 3.8.12, pytest-6.2.5, py-1.10.0, pluggy-1.0.0
rootdir: /home/droe/GitHub/pytraj
collected 3 items                                                                                           

test_nastruct.py .F.                                                                                  [100%]

================================================= FAILURES ==================================================
_________________________________________ TestNastruct.test_nupars __________________________________________

self = <test_analysis.test_nastruct.TestNastruct testMethod=test_nupars>

    def test_nupars(self):
        pdb_fn = fn('Test_NAstruct/adh026.3.pdb')
        traj = pt.iterload(pdb_fn, pdb_fn)
        data = pt.nastruct(traj)
    
        # default
        text = '''
        parm {}
        trajin {}
        nastruct groovecalc 3dna
        '''.format(
            fn('Test_NAstruct/adh026.3.pdb'), fn('Test_NAstruct/adh026.3.pdb'))
    
        state = pt.load_cpptraj_state(text)
        state.run()
    
        for key in ['major', 'minor', 'twist']:
            cpp_data = np.array(
                [x.values for x in state.data if x.aspect == key])
            # need to transpose to get shape=(n_frames, n_pairs)
            cpp_data = cpp_data.T
            aa_eq(data[key][1], cpp_data)
    
        # TODO: assert
>       data._summary(np.mean, indices=None)

test_nastruct.py:37: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
/home/droe/Programs/anaconda3/envs/pytrajtest/lib/python3.8/site-packages/pytraj-2.0.6.dev0-py3.8-linux-x86_64.egg/pytraj/analysis/nucleic_acid_analysis.py:225: in _summary
    values = self[k][1]
/home/droe/Programs/anaconda3/envs/pytrajtest/lib/python3.8/site-packages/pytraj-2.0.6.dev0-py3.8-linux-x86_64.egg/pytraj/analysis/nucleic_acid_analysis.py:177: in __getitem__
    return self.__getattr__(key)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

self = <nupars, keys = ['pucker', 'shear', 'stretch', 'stagger', 'buckle', 'prop', 'open', 'hb', 'bp', 'nxyz', 'shift', 'slide', 'rise', 'tilt', 'roll', 'twist', 'xdisp', 'ydisp', 'hrise', 'incl', 'tip', 'htwist', 'zp', 'major', 'minor']>
aspect = 'nxyz'

    def __getattr__(self, aspect):
        '''self.nxyz, ...
        '''
        # data is a list of DataArray
        data = self._dict[aspect]
        arr = np.empty((len(data), len(data[0])), dtype='f8')
    
        keylist = []
        for idx, arr0 in enumerate(data):
            keylist.append(arr0.key)
>           arr[idx] = arr0.values
E           ValueError: could not broadcast input array from shape (3,3) into shape (3,)

/home/droe/Programs/anaconda3/envs/pytrajtest/lib/python3.8/site-packages/pytraj-2.0.6.dev0-py3.8-linux-x86_64.egg/pytraj/analysis/nucleic_acid_analysis.py:189: ValueError
========================================== short test summary info ==========================================
FAILED test_nastruct.py::TestNastruct::test_nupars - ValueError: could not broadcast input array from shap...
======================================== 1 failed, 2 passed in 0.34s ========================================

My guess is that the underlying issue is that the __getattr__ function in nucleic_acid_anlysis.py chokes because it is trying to convert everything into 1D scalar arrays, i.e. arr = np.empty((len(data), len(data[0])), dtype='f8'), and when it encounters the [nxyz] data which is DataSet_Vector instead of a DataSet_1D type it chokes. Can we add something like if aspect == [nxyz] return ndarray where ndarray has dimension [len, 6], with the 6 values being the vector x y z and origin x y z?

[hainm]

thanks Dan. I will try to update pytraj tomorrow.

[drroe]

@hainm Since this might not be a quick fix I'm going to try to revert this part of the cpptraj PR for now behind an ifdef. I want to try to get this into the Amber release and the deadline is today. Let's keep it on the radar though since I think it's useful info to have.

[hainm]

@drroe it turns out the change in pytraj is quite simple (I only run the affected test)

$ git diff origin/master  -- pytraj/analysis/nucleic_acid_analysis.py
diff --git a/pytraj/analysis/nucleic_acid_analysis.py b/pytraj/analysis/nucleic_acid_analysis.py
index e4a9641ea..4ba6bafe1 100644
--- a/pytraj/analysis/nucleic_acid_analysis.py
+++ b/pytraj/analysis/nucleic_acid_analysis.py
@@ -190,7 +190,7 @@ class nupars(object):
         return keylist, arr.T
 
     def keys(self):
-        return list(self._dict)
+        return [k for k in list(self._dict) if k != 'nxyz']
 
     def __dir__(self):
         '''for autocompletion in ipython

[hainm]

@drroe I've updated pytraj. Please try again.