def test_disulfide(self):
     """Test handling of disulfide bonds"""
     mlib = self.get_all_libraries()
     bsep = mdt.features.AtomBondSeparation(mlib,
                                     bins=mdt.uniform_bins(20, 0, 1.0))
     bsep_ss = mdt.features.AtomBondSeparation(mlib,
                                     bins=mdt.uniform_bins(20, 0, 1.0),
                                     disulfide=True)
     env = self.get_environ()
     mdl = modeller.model(env)
     mdl.build_sequence('CC')
     # When SG-SG distance is small enough, an extra bond
     # (separation feature = 1) should be detected, but only with
     # disulfide=True
     for (dist, num) in [(2.6, 11.0), (2.4, 12.0)]:
         sg1 = mdl.residues[0].atoms['SG']
         sg2 = mdl.residues[1].atoms['SG']
         sg1.x = sg1.y = sg1.z = 0.
         sg2.x = sg2.y = 0.
         sg2.z = dist
         a = modeller.alignment(env)
         a.append_model(mdl, atom_files='test', align_codes='test')
         m = mdt.Table(mlib, features=bsep)
         m.add_alignment(a, residue_span_range=(-999,0,0,999))
         self.assertEqual(m[1], 11.0)
         m2 = mdt.Table(mlib, features=bsep_ss)
         m2.add_alignment(a, residue_span_range=(-999,0,0,999))
         self.assertEqual(m2[1], num)

 def test_feature_sidechain_biso(self):
     """Check average sidechain Biso feature"""
     env = self.get_environ()
     mlib = self.get_mdt_library()
     self.assertRaises(ValueError, mdt.features.SidechainBiso, mlib, bins=mdt.uniform_bins(5, 0, 10), protein=3)
     sidechain_biso = mdt.features.SidechainBiso(mlib, bins=mdt.uniform_bins(5, 0, 10))
     mdl = modeller.model(env)
     mdl.build_sequence("A")
     aln = modeller.alignment(env)
     aln.append_model(mdl, align_codes="test")
     s = aln[0]
     # Mainchain atom Biso should be ignored:
     for mainchain in ("N:1", "C:1", "O:1", "OXT:1", "CA:1"):
         s.atoms[mainchain].biso = 1000
     for (biso, bin) in (
         (22, 2),
         (32, 3),  # Map regular values to bins
         (0, -1),  # Zero Biso should be "undefined"
         (1, 3),
     ):  # Biso < 2 is multiplied by 4pi^2
         s.atoms["CB:1"].biso = biso
         m = mdt.Table(mlib, features=sidechain_biso)
         m.add_alignment(aln)
         self.assertEqual(m.shape, (6,))
         self.assertEqual(m.sum(), 1)
         self.assertEqual(m[bin], 1)

 def test_feature_resind_diff(self):
     """Test the residue index difference feature"""
     env = self.get_environ()
     mlib = self.get_mdt_library()
     diff = mdt.features.ResidueIndexDifference(mlib, bins=mdt.uniform_bins(21, -10, 1))
     absdiff = mdt.features.ResidueIndexDifference(mlib, absolute=True, bins=mdt.uniform_bins(21, -10, 1))
     aln = modeller.alignment(env, file="test/data/alignment.ali", align_codes="5fd1")
     m1 = mdt.Table(mlib, features=diff)
     m2 = mdt.Table(mlib, features=absdiff)
     self.assertEqual(m1.symmetric, False)
     self.assertEqual(m2.symmetric, True)
     m1.add_alignment(aln, residue_span_range=(-999, -2, 2, 999))
     m2.add_alignment(aln, residue_span_range=(-999, -2, 2, 999))
     self.assertEqual(m1.sum(), 10920)
     self.assertEqual(m2.sum(), 5460)
     # span range should result in 0, +/- 1 bins being zero:
     for m in (m1, m2):
         self.assertEqual(m[9], 0.0)
         self.assertEqual(m[10], 0.0)
         self.assertEqual(m[11], 0.0)
     # Non-absolute feature should have other bins symmetrically distributed:
     for i in range(9):
         self.assertEqual(m1[i], m[-2 - i])
     # Absolute feature should have no negative values:
     for i in range(9):
         self.assertEqual(m2[i], 0.0)

 def test_feature_hbond(self):
     """Check hydrogen bond features"""
     env = self.get_environ()
     mlib = self.get_mdt_library()
     mlib.hbond_classes.read("data/atmcls-hbda.lib")
     donor = mdt.features.HydrogenBondDonor(mlib, mdt.uniform_bins(7, 1.0, 1.0))
     accep = mdt.features.HydrogenBondAcceptor(mlib, mdt.uniform_bins(7, 1.0, 1.0))
     totchg = mdt.features.HydrogenBondCharge(mlib, mdt.uniform_bins(9, 1.0, 1.0))
     satisf = mdt.features.HydrogenBondSatisfaction(mlib, mdt.uniform_bins(100, 0.0, 10.0))
     self.assertRaises(mdt.MDTError, mlib.hbond_classes.read, "data/atmcls-hbda.lib")
     m = mdt.Table(mlib, features=donor)
     m2 = mdt.Table(mlib, features=accep)
     m3 = mdt.Table(mlib, features=satisf)
     m4 = mdt.Table(mlib, features=totchg)
     aln = modeller.alignment(env, file="test/data/alignment.ali")
     m.add_alignment(aln)
     m2.add_alignment(aln)
     m3.add_alignment(aln)
     m4.add_alignment(aln)
     self.assertInTolerance(m[0], 295.0, 0.0005)
     self.assertInTolerance(m[1], 139.0, 0.0005)
     self.assertEqual(m[-1], 349.0)
     self.assertInTolerance(m2[0], 236.0, 0.0005)
     self.assertInTolerance(m2[1], 223.0, 0.0005)
     self.assertEqual(m2[-1], 168.0)
     self.assertInTolerance(m3[0], 1.0, 0.0005)
     self.assertInTolerance(m3[1], 0.0, 0.0005)
     self.assertEqual(m3[-1], 0.0)
     self.assertInTolerance(m4[0], 78.0, 0.0005)
     self.assertInTolerance(m4[1], 24.0, 0.0005)
     self.assertEqual(m4[-1], 739.0)
     # Exercise writing of hbond information to HDF5 files:
     for t in (m, m2, m3, m4):
         t.write_hdf5("test.hdf5")
         os.unlink("test.hdf5")

 def test_feature_triplet_type(self):
     """Check triplet type features"""
     env = self.get_environ()
     mlib = self.get_mdt_library()
     mlib.tuple_classes.read("data/trpcls.lib")
     tuple_type = mdt.features.TupleType(mlib)
     tuple_type2 = mdt.features.TupleType(mlib, pos2=True)
     tuple_dist = mdt.features.TupleDistance(mlib, bins=mdt.uniform_bins(9, 2.0, 0.2))
     tuple_angle1 = mdt.features.TupleAngle1(mlib, bins=mdt.uniform_bins(6, 0, 30.0))
     tuple_dihed1 = mdt.features.TupleDihedral1(mlib, bins=mdt.uniform_bins(6, -180, 60.0))
     tuple_dihed2 = mdt.features.TupleDihedral2(mlib, bins=mdt.uniform_bins(6, -180, 60.0))
     tuple_dihed3 = mdt.features.TupleDihedral3(mlib, bins=mdt.uniform_bins(6, -180, 60.0))
     self.assertRaises(mdt.MDTError, mlib.tuple_classes.read, "data/trpcls.lib")
     m1 = mdt.Table(mlib, features=tuple_type)
     m2 = mdt.Table(mlib, features=tuple_type2)
     m3 = mdt.Table(mlib, features=tuple_dist)
     m4 = mdt.Table(mlib, features=tuple_angle1)
     m5 = mdt.Table(mlib, features=tuple_dihed1)
     m6 = mdt.Table(mlib, features=tuple_dihed2)
     m7 = mdt.Table(mlib, features=tuple_dihed3)
     aln = modeller.alignment(env, file="test/data/tiny.ali")
     for m in (m1, m2, m3, m4, m5, m6, m7):
         m.add_alignment(aln, residue_span_range=(-9999, 0, 0, 9999))
     self.assertInTolerance(m1[0], 1.0, 0.0005)
     self.assertInTolerance(m1[1], 0.0, 0.0005)
     self.assertInTolerance(m1[2], 1.0, 0.0005)
     self.assertEqual(m1.shape, (236,))
     self.assertEqual(m1[-1], 0.0)
     self.assertInTolerance(m2[0], 60.0, 0.0005)
     self.assertInTolerance(m2[1], 0.0, 0.0005)
     self.assertInTolerance(m2[2], 60.0, 0.0005)
     self.assertEqual(m2.shape, (236,))
     self.assertEqual(m2[-1], 0.0)
     self.assertInTolerance(m3[0], 0.0, 0.0005)
     self.assertInTolerance(m3[1], 82.0, 0.0005)
     self.assertInTolerance(m3[2], 226.0, 0.0005)
     self.assertEqual(m3.shape, (10,))
     self.assertInTolerance(m3[-1], 3018.0, 0.0005)
     self.assertInTolerance(m4[0], 479.0, 0.0005)
     self.assertInTolerance(m4[1], 806.0, 0.0005)
     self.assertInTolerance(m4[2], 471.0, 0.0005)
     self.assertEqual(m4.shape, (7,))
     self.assertEqual(m4[-1], 0.0)
     self.assertInTolerance(m5[0], 556.0, 0.0005)
     self.assertInTolerance(m5[1], 642.0, 0.0005)
     self.assertInTolerance(m5[2], 470.0, 6.0005)
     self.assertEqual(m5.shape, (7,))
     self.assertInTolerance(m5[-1], 180.0, 0.0005)
     self.assertInTolerance(m6[0], 661.0, 0.0005)
     self.assertInTolerance(m6[1], 520.0, 0.0005)
     self.assertInTolerance(m6[2], 545.0, 6.0005)
     self.assertEqual(m6.shape, (7,))
     self.assertInTolerance(m6[-1], 112.0, 0.0005)
     self.assertInTolerance(m7[0], 661.0, 0.0005)
     self.assertInTolerance(m7[1], 520.0, 0.0005)
     self.assertInTolerance(m7[2], 545.0, 6.0005)
     self.assertEqual(m7.shape, (7,))
     self.assertInTolerance(m7[-1], 112.0, 0.0005)

 def test_feature_combination(self):
     """Check that invalid feature combinations are rejected"""
     mlib = self.get_mdt_library()
     atmdist = mdt.features.AtomDistance(mlib,
                                         bins=mdt.uniform_bins(60, 0, 0.5))
     resdist = mdt.features.ResidueDistance(mlib, protein=1,
                                            bins=mdt.uniform_bins(7, 0, 2.0))
     self.assertRaises(ValueError, self.get_test_mdt, mlib,
                       features=(resdist,atmdist))

    def test_bond_span_range(self):
        """Test bond_span_range argument"""
        env = self.get_environ()
        mdl = model(env)
        mdl.build_sequence('A')
        aln = alignment(env)
        aln.append_model(mdl, align_codes='test')
        mlib = self.get_mdt_library()
        mlib.bond_classes.read('data/bndgrp.lib')
        dist = mdt.features.AtomDistance(mlib,
                                         bins=mdt.uniform_bins(60, 0, 0.5))

        # Only 4 direct chemical bonds (N-CA, CA-CB, CA-C, C-O) in ALA; note
        # that bond library does not include OXT so C-OXT interaction is
        # excluded
        m = mdt.Table(mlib, features=dist)
        m.add_alignment(aln, bond_span_range=(1,1),
                        residue_span_range=(0,0,0,0))
        self.assertEqual(m.sample_size, 4.0)

        # Only 2 dihedrals (N-CA-C-O, O-C-CA-CB)
        m = mdt.Table(mlib, features=dist)
        m.add_alignment(aln, bond_span_range=(3,3),
                        residue_span_range=(0,0,0,0))
        self.assertEqual(m.sample_size, 2.0)

        # 4 bonds, 4 angles and 2 dihedrals: 10 in total
        m = mdt.Table(mlib, features=dist)
        m.add_alignment(aln, bond_span_range=(1,3),
                        residue_span_range=(0,0,0,0))
        self.assertEqual(m.sample_size, 10.0)

        # Check for bonds between residues (just the N-C bond here)
        mdl = model(env)
        mdl.build_sequence('AA')
        aln = alignment(env)
        aln.append_model(mdl, align_codes='test')

        # Force a non-symmetric scan (to check handling of bond separation
        # regardless of which order atom indices are in)
        diff = mdt.features.ResidueIndexDifference(mlib,
                                      bins=mdt.uniform_bins(21, -10, 1))
        m = mdt.Table(mlib, features=(dist,diff))
        m.add_alignment(aln, bond_span_range=(0,1),
                        residue_span_range=(-10,-1,1,10))
        self.assertEqual(m.sample_size, 2.0)

        # Bonds never span chains
        mdl = model(env)
        mdl.build_sequence('A/A')
        aln = alignment(env)
        aln.append_model(mdl, align_codes='test')
        m = mdt.Table(mlib, features=dist)
        m.add_alignment(aln, bond_span_range=(0,99999),
                        residue_span_range=(-10,-1,1,10))
        self.assertEqual(m.sample_size, 0.0)

 def test_feature_alpha_content(self):
     """Check alpha content feature"""
     env = self.get_environ()
     mlib = self.get_mdt_library()
     self.assertRaises(ValueError, mdt.features.AlphaContent, mlib, bins=mdt.uniform_bins(10, 0, 0.1), protein=3)
     alpha = mdt.features.AlphaContent(mlib, bins=mdt.uniform_bins(10, 0, 0.1))
     for (alnfile, bin) in (("tiny.ali", 0), ("alignment.ali", 5)):
         m = mdt.Table(mlib, features=alpha)
         a = modeller.alignment(env, file=os.path.join("test", "data", alnfile))
         m.add_alignment(a)
         self.assertEqual(m.shape, (11,))
         self.assertEqual(m.sum(), 1)
         self.assertEqual(m[bin], 1)

    def test_write_2dsplinelib(self):
        """Test the write_2dsplinelib function"""
        mlib = self.get_mdt_library()

        restype = mdt.features.ResidueType(mlib)
        phi = mdt.features.PhiDihedral(mlib,
                                       bins=mdt.uniform_bins(18, -180.0, 40.0))
        psi = mdt.features.PsiDihedral(mlib,
                                       bins=mdt.uniform_bins(18, -180.0, 40.0))
        m = self.get_test_mdt(mlib, [restype, phi, psi])
        mdt.write_2dsplinelib(open('test.out', 'w'), m)

        # Make sure that valid Python code was produced
        code = compile(open('test.out').read(), 'test.out', 'exec')
        os.unlink('test.out')

 def test_abstract(self):
     """Should not be able to instantiate abstract features"""
     env = self.get_environ()
     mlib = self.get_mdt_library()
     bins = mdt.uniform_bins(10, 0, 1.0)
     for feat in [
         mdt.features.Protein,
         mdt.features.ProteinPair,
         mdt.features.Residue,
         mdt.features.ResiduePair,
         mdt.features.AlignedResidue,
         mdt.features.AlignedResiduePair,
         mdt.features.Atom,
         mdt.features.AtomPair,
         mdt.features.Tuple,
         mdt.features.TuplePair,
         mdt.features.ChemicalBond,
     ]:
         self.assertRaises(TypeError, feat, mlib, bins)
     for feat in [
         mdt.features.ResidueFixedBins,
         mdt.features.AtomFixedBins,
         mdt.features.TupleFixedBins,
         mdt.features.ChemicalBondFixedBins,
     ]:
         self.assertRaises(TypeError, feat, mlib)

    def test_integrate(self):
        """Make sure MDT integration works"""
        env = self.get_environ()
        mlib = self.get_mdt_library()
        restyp0 = mdt.features.ResidueType(mlib, protein=0)
        restyp1 = mdt.features.ResidueType(mlib, protein=1)
        chi1 = mdt.features.Chi1Dihedral(mlib,
                                         mdt.uniform_bins(36, -180, 10))
        m1 = mdt.Table(mlib, features=restyp0)
        m2 = mdt.Table(mlib, features=restyp1)
        mboth = mdt.Table(mlib, features=(restyp0,restyp1))
        seq1 = "AFVVTDNCIK"
        seq2 = "DCVEVCPVDC"
        aln = alignment(env)
        aln.append_sequence(seq1)
        aln.append_sequence(seq2)
        for m in (m1, m2, mboth):
            m.add_alignment(aln)

        # Number of features must be correct:
        for features in ((), (restyp0,restyp1), (restyp0,restyp1,chi1)):
            self.assertRaises(ValueError, mboth.integrate, features=features)
        # Features must exist in input MDT:
        self.assertRaises(ValueError, mboth.integrate, features=chi1)
        m1int = mboth.integrate(features=restyp0)
        self.assertMDTsEqual(m1, m1int)
        m2int = mboth.integrate(features=restyp1)
        self.assertMDTsEqual(m2, m2int)

    def test_tuple_pair_bond_span_range(self):
        """Test bond_span_range with tuple pair scan"""
        env = self.get_environ()
        mdl = model(env)
        mdl.build_sequence('A')
        aln = alignment(env)
        aln.append_model(mdl, align_codes='test')
        mlib = self.get_mdt_library()
        mlib.bond_classes.read('data/bndgrp.lib')
        mlib.tuple_classes.read('data/trpcls.lib')
        typ = mdt.features.TupleType(mlib)
        typ2 = mdt.features.TupleType(mlib, pos2=True)
        dist = mdt.features.TupleDistance(mlib,
                                          bins=mdt.uniform_bins(9, 2.0, 0.2))

        m = mdt.Table(mlib, features=dist)
        m.add_alignment(aln, residue_span_range=(0,0,0,0))
        self.assertEqual(m.sample_size, 10.0)

        m = mdt.Table(mlib, features=dist)
        m.add_alignment(aln, bond_span_range=(1,1),
                        residue_span_range=(0,0,0,0))
        # Bond span should restrict interactions to 6
        # (C:CA:CB-CA:C:O, CA:C:O-N:CA:C, CA:C:O-N:CA:CB, and the reverse)
        self.assertEqual(m.sample_size, 6.0)

 def test_chain_span_range(self):
     """Test chain_span_range argument"""
     env = self.get_environ()
     mdl = model(env)
     mdl.build_sequence('G/G')
     aln = alignment(env)
     aln.append_model(mdl, align_codes='test')
     mlib = self.get_mdt_library()
     mlib.tuple_classes.read('data/dblcls.lib')
     tuple_dist = mdt.features.TupleDistance(mlib,
                                       bins=mdt.uniform_bins(49, 2.0, 0.2))
     # All chain differences should be out of range, so this MDT should
     # end up empty:
     m = mdt.Table(mlib, features=tuple_dist)
     m.add_alignment(aln, chain_span_range=(-999, -999, 999, 999),
                     residue_span_range=(-999, 0, 0, 999))
     self.assertEqual(m.sum(), 0.0)
     # Default chain separation should allow intra-chain interactions, so
     # should yield more (56) than only allowing inter-chain
     # interactions (32)
     m = mdt.Table(mlib, features=tuple_dist)
     m.add_alignment(aln, residue_span_range=(-999, 0, 0, 999))
     self.assertEqual(m.sum(), 56.0)
     m = mdt.Table(mlib, features=tuple_dist)
     m.add_alignment(aln, chain_span_range=(-999, -1, 1, 999),
                     residue_span_range=(-999, 0, 0, 999))
     self.assertEqual(m.sum(), 32.0)

    def test_add(self):
        """Check adding MDTs"""
        mlib = self.get_mdt_library()
        bins = mdt.uniform_bins(3, -1.0, 1.5)
        xray0 = mdt.features.XRayResolution(mlib, bins, protein=0)
        xray1 = mdt.features.XRayResolution(mlib, bins, protein=1)
        m1 = mdt.Table(mlib, features=xray0)
        for (n, val) in enumerate((1,2,3,4)):
            m1[n] = val
        m2 = mdt.Table(mlib, features=xray0)
        for (n, val) in enumerate((10,20,30,40)):
            m2[n] = val
        m3 = m1 + m2
        m1 += m2
        self.assertMDTsEqual(m1, m3)
        for (n, val) in enumerate((11,22,33,44)):
            self.assertEqual(m3[n], val)
        # Cannot add if numbers of features are different
        badmdt = mdt.Table(mlib, features=(xray0, xray1))
        self.assertRaises(ValueError, m1.__add__, badmdt)

        # Cannot add if feature types are different
        badmdt = mdt.Table(mlib, features=xray1)
        self.assertRaises(ValueError, m1.__add__, badmdt)

        # Cannot add if starts are different
        badmdt = m2.reshape(features=xray0, offset=1, shape=0)
        self.assertRaises(ValueError, m1.__add__, badmdt)

        # Cannot add if nbins are different
        badmdt = m2.reshape(features=xray0, offset=0, shape=-1)
        self.assertRaises(ValueError, m1.__add__, badmdt)

    def test_bond_span_range_disulfide(self):
        """Test bond_span_range argument with disulfides"""
        env = self.get_environ()
        mdl = model(env)
        mdl.read('1HEL.pdb')
        aln = alignment(env)
        aln.append_model(mdl, align_codes='test')
        mlib = self.get_mdt_library()
        mlib.bond_classes.read('data/bndgrp.lib')
        dist = mdt.features.AtomDistance(mlib,
                                         bins=mdt.uniform_bins(60, 0, 0.5))
        # Four disulfide bond in this structure
        m = mdt.Table(mlib, features=dist)
        m.add_alignment(aln, bond_span_range=(1,1),
                        residue_span_range=(-9999,0,0,9999))

        m2 = mdt.Table(mlib, features=dist)
        m2.add_alignment(aln, bond_span_range=(1,1),
                        residue_span_range=(-9999,0,0,9999), disulfide=True)
        self.assertEqual(m2.sample_size-m.sample_size, 4.0)

        m = mdt.Table(mlib, features=dist)
        m.add_alignment(aln, bond_span_range=(3,3),
                        residue_span_range=(-9999,0,0,9999))

        m2 = mdt.Table(mlib, features=dist)
        m2.add_alignment(aln, bond_span_range=(3,3),
                        residue_span_range=(-9999,0,0,9999), disulfide=True)
        self.assertEqual(m2.sample_size-m.sample_size, 12.0)

    def test_feature_residue_distance_difference(self):
        """Check residue-residue distance difference feature"""
        env = self.get_environ()
        mlib = self.get_mdt_library()
        ddist = mdt.features.ResidueDistanceDifference(mlib, bins=mdt.uniform_bins(20, -10, 1))
        aln = modeller.alignment(env, file="test/data/struc-struc.ali")
        m = mdt.Table(mlib, features=ddist)
        m.add_alignment(aln)
        self.assertEqual(m[9], 20)
        self.assertEqual(m[10], 20)
        self.assertEqual(sum([b for b in m]), 40)
        self.assertEqual(m[-1], 0)

        # Undefined (-999) coordinates in either structure should put
        # features in the undefined bin
        oldx = aln[0].residues[0].atoms["CA"].x
        aln[0].residues[0].atoms["CA"].x = -999
        m = mdt.Table(mlib, features=ddist)
        m.add_alignment(aln)
        self.assertEqual(m[-1], 16)

        aln[0].residues[0].atoms["CA"].x = oldx
        aln[1].residues[0].atoms["CA"].x = -999
        m = mdt.Table(mlib, features=ddist)
        m.add_alignment(aln)
        self.assertEqual(m[-1], 16)

 def test_symmetric(self):
     """Test symmetric/asymmetric residue pair features"""
     env = self.get_environ()
     mlib = self.get_mdt_library()
     bins = mdt.uniform_bins(10, 0, 1.0)
     dist = mdt.features.ResidueDistance(mlib, bins)
     avresacc = mdt.features.AverageResidueAccessibility(mlib, bins)
     avndif = mdt.features.AverageNeighborhoodDifference(mlib, bins)
     diff = mdt.features.ResidueIndexDifference(mlib, bins)
     ddist = mdt.features.ResidueDistanceDifference(mlib, bins)
     avgapdist = mdt.features.AverageGapDistance(mlib, bins)
     sym_features = (avndif, avresacc, avgapdist)
     asym_features = (dist, ddist, diff)
     for a in sym_features:
         m = mdt.Table(mlib, features=a)
         self.assertEqual(m.symmetric, True)
     for a in asym_features:
         m = mdt.Table(mlib, features=a)
         self.assertEqual(m.symmetric, False)
     # Combinations are only symmetric if all features are symmetric:
     for (a, b, symm) in (
         (sym_features[0], sym_features[1], True),
         (asym_features[0], asym_features[1], False),
         (sym_features[0], asym_features[0], False),
     ):
         m = mdt.Table(mlib, features=(a, b))
         self.assertEqual(m.symmetric, symm)

 def test_feature_dihedral_type(self):
     """Check dihedral type features"""
     env = self.get_environ()
     mlib = self.get_mdt_library()
     mlib.dihedral_classes.read("data/impgrp.lib")
     dihedtype = mdt.features.DihedralType(mlib)
     dihedral = mdt.features.Dihedral(mlib, bins=mdt.uniform_bins(288, -180, 1.25))
     self.assertRaises(mdt.MDTError, mlib.dihedral_classes.read, "data/impgrp.lib")
     m = mdt.Table(mlib, features=dihedtype)
     m2 = mdt.Table(mlib, features=dihedral)
     aln = modeller.alignment(env, file="test/data/alignment.ali")
     m.add_alignment(aln)
     m2.add_alignment(aln)
     self.assertInTolerance(m[0], 7.0, 0.0005)
     self.assertInTolerance(m[2], 9.0, 0.0005)
     self.assertInTolerance(m[4], 11.0, 0.0005)
     self.assertEqual(m.shape, (79,))
     self.assertEqual(m[-1], 0.0)
     self.assertInTolerance(m2[143], 60.0, 1.0005)
     self.assertInTolerance(m2[144], 53.0, 1.0005)
     self.assertInTolerance(m2[145], 24.0, 0.0005)
     self.assertEqual(m2.shape, (289,))
     self.assertEqual(m2[-1], 0.0)
     # Exercise writing of dihedral class information to HDF5 files:
     m.write_hdf5("test.hdf5")
     os.unlink("test.hdf5")

    def test_external(self):
        """Check bond separation feature between residues"""
        mlib = self.get_all_libraries()
        attyp1 = mdt.features.AtomType(mlib)
        attyp2 = mdt.features.AtomType(mlib, pos2=True)
        bsep = mdt.features.AtomBondSeparation(mlib,
                                        bins=mdt.uniform_bins(20, 0, 1.0))
        m = self.build_mdt_from_sequence(mlib, [attyp1, attyp2, bsep],
                        'ARN', residue_span_range=(-999,-1,1,999))
        atom_types = {}
        for n, b in enumerate(m.features[0].bins):
            atom_types[b.symbol] = n

        def assertBondSep(at1, at2, numbond, sep):
            bins = [b for b in m[atom_types[at1]][atom_types[at2]]]
            self.assertEqual(sum(bins), numbond)
            self.assertEqual(bins[sep], numbond)

        # Check known bond separations between adjacent residues
        assertBondSep('AN', 'RN', numbond=1, sep=3)
        assertBondSep('AC', 'RN', numbond=1, sep=1)

        # Check known bond separations between non-adjacent residues
        assertBondSep('AN', 'NN', numbond=1, sep=6)
        assertBondSep('AN', 'NC', numbond=1, sep=8)
        assertBondSep('ACB', 'NCG', numbond=1, sep=9)

    def test_feature_iresol(self):
        """Check resolution features"""
        env = self.get_environ()
        mlib = self.get_mdt_library()
        bins = mdt.uniform_bins(3, -1.0, 1.5)
        xray0 = mdt.features.XRayResolution(mlib, bins, protein=0)
        xray0_nmr = mdt.features.XRayResolution(mlib, bins, protein=0, nmr=1.0)
        xray1 = mdt.features.XRayResolution(mlib, bins, protein=1)
        xray2 = mdt.features.XRayResolution(mlib, bins, protein=2)
        # Check valid range for protein argument
        for p in (-1, 3):
            self.assertRaises(ValueError, mdt.features.XRayResolution, mlib, bins, protein=p)
        m = self.get_test_mdt(mlib, features=xray0)
        m2 = self.get_test_mdt(mlib, features=xray1)
        self.assertEqual(m.shape, (4,))
        self.assertEqual([b for b in m], [0.0, 1.0, 1.0, 0.0])
        self.assertMDTDataEqual(m, m2)

        for (code, feat, bin) in (
            ("bin0", xray0, 0),
            ("bin0", xray0_nmr, 1),
            ("bin1", xray0, 1),
            ("bin2", xray0, 2),
            ("undef1", xray0, 3),
            ("undef2", xray0, 3),
        ):
            m = mdt.Table(mlib, features=feat)
            aln = modeller.alignment(env, file="test/data/resol.ali", align_codes=code)
            m.add_alignment(aln)
            self.assertEqual(m[bin], 1.0)