Merge branch 'pinyili2' into 'master'

Pinyili2

See merge request !1

.gitignore

 __pycache__
 .backup
-
+*.swp
+.ipynb*
+.nfs*
+.DS_Store
 /build
 *.egg-info

README.md

@@ -15,6 +15,8 @@ Alternatively, the package can be used to initialize simulations using the oxDNA
   * mdanalysis >= 0.18 
   * cadnano >= 2.5.2.1
   * appdirs >= 1.4
+  * pandas >= 1.0
+  * scadnano >= 0.1 (optional)
   * [oxDNA](https://dna.physics.ox.ac.uk/index.php/Main_Page) (optional)
 
 ## Installation
@@ -41,8 +43,8 @@ cd ../../
 
 ## Possibly setup python dependencies (you may want to set up a virtual environement) using commands below
 # conda install "numpy>=1.14" "scipy>=1.1" "appdirs>=1.4"
-# conda install -c conda-forge "mdanalysis>=0.18" 
-# pip3 install "cadnano>=2.5.2.1"
+# conda install -c conda-forge "mdanalysis>=0.18" "pandas>=1.0"
+# pip3 install "cadnano>=2.5.2.1" "scadnano>=0.1" 
 
 git clone https://gitlab.engr.illinois.edu/tbgl/tools/mrdna
 cd mrdna

mrdna/RELEASE-VERSION

-1.0a.dev71
+1.0a.dev74

mrdna/__init__.py

@@ -7,11 +7,11 @@ def _get_username():
     except:
         return None
 
-logging.basicConfig(format='%(name)s: %(levelname)s: %(message)s')
+logging.basicConfig(format='%(asctime)s %(name)s: %(levelname)s: %(message)s')
 logger = logging.getLogger(__name__)
 logger.setLevel(logging.INFO)
 _ch = logging.StreamHandler()
-_ch.setFormatter(logging.Formatter('%(name)s: %(levelname)s: %(message)s'))
+_ch.setFormatter(logging.Formatter('%(asctime)s %(name)s: %(levelname)s: %(message)s'))
 logger.addHandler(_ch)
 logger.propagate = False
 

mrdna/arbdmodel/coords.py

@@ -110,6 +110,52 @@ def quaternion_to_matrix(q):
     return np.array(R)
 
 def quaternion_from_matrix( R ):
+    q = np.empty(4)
+    if R[2,2] < 0:
+        if R[0,0] > R[1,1]:
+            trace = 1.0 + R[0,0] - R[1,1] - R[2,2]
+            s = 2.0 * np.sqrt(trace)
+            if R[1,2] < R[2,1]: s = -s
+            q[0] = (R[1,2] - R[2,1]) / s
+            q[1] = 0.25 * s
+            q[2] = (R[0,1] + R[1,0]) / s
+            q[3] = (R[2,0] + R[0,2]) / s
+            if np.isclose(trace,1) and np.all(np.isclose([x for i,x in enumerate(q) if i != 1],0)):
+                q[1] = 1
+        else:
+            trace = 1.0 - R[0,0] + R[1,1] - R[2,2]
+            s = 2.0 * np.sqrt(trace)
+            if R[2,0] < R[0,2]: s = -s
+            q[0] = (R[2,0] - R[0,2]) / s
+            q[1] = (R[0,1] + R[1,0]) / s
+            q[2] = 0.25 * s
+            q[3] = (R[1,2] + R[2,1]) / s
+            if np.isclose(trace,1) and np.all(np.isclose([x for i,x in enumerate(q) if i != 2],0)):
+                q[2] = 1
+    else:
+        if R[0,0] < -R[1,1]:
+            trace = 1.0 - R[0,0] - R[1,1] + R[2,2]
+            s = 2.0 * np.sqrt(trace)
+            if R[0,1] < R[1,0]: s = -s
+            q[0] = (R[0,1] - R[1,0]) / s
+            q[1] = (R[2,0] + R[0,2]) / s
+            q[2] = (R[1,2] + R[2,1]) / s
+            q[3] = 0.25 * s
+            if np.isclose(trace,1) and np.all(np.isclose([x for i,x in enumerate(q) if i != 3],0)):
+                q[3] = 1
+        else:
+            trace = 1.0 + R[0,0] + R[1,1] + R[2,2]
+            s = 2.0 * np.sqrt(trace)
+            q[0] = 0.25 * s
+            q[1] = (R[1,2] - R[2,1]) / s
+            q[2] = (R[2,0] - R[0,2]) / s
+            q[3] = (R[0,1] - R[1,0]) / s
+            if np.isclose(trace,1) and np.all(np.isclose([x for i,x in enumerate(q) if i != 0],0)):
+                q[0] = 1
+    assert( q[0] >= 0 )
+    return q
+
+def __quaternion_from_matrix__deprecated( R ):
     e1 = R[0]
     e2 = R[1]
     e3 = R[2]
@@ -165,9 +211,13 @@ def quaternion_inverse(q):
 
 def quaternion_exp(q,t):
     assert(len(q) == 4)
+    if q[0] > 1 and np.isclose(q[0],1): q[0] = 1
+    if q[0] < -1 and np.isclose(q[0],-1): q[0] = -1
     omega = np.arccos(q[0])
     v = q[1:]
-    v = v/np.linalg.norm(v)
+    vl = np.linalg.norm(v)
+    if vl != 0:
+        v = v/vl
     qexp = np.empty(4)
     qexp[0] = np.cos(omega*t)
     qexp[1:] = np.sin(omega*t)*v

mrdna/readers/__init__.py

 ## TODO: make module this package conform to a single style for input/output
 
-def read_cadnano(json_file, sequence=None, fill_sequence='T', **model_parameters):
+def read_cadnano_legacy(json_file, sequence=None, fill_sequence='T', **model_parameters):
     from .cadnano_segments import read_json_file
     from .cadnano_segments import read_model as model_from_cadnano_json
 
     data = read_json_file(json_file)
     return model_from_cadnano_json(data, sequence, fill_sequence, **model_parameters)
 
+def read_cadnano(json_file,sequence=None,**model_parameters):
+    from .segmentmodel_from_cadnano import mrdna_model_from_cadnano
+    
+    return mrdna_model_from_cadnano(json_file,seq=sequence,**model_parameters)
+def read_scadnano(sc_file, **model_parameters):
+    from .segmentmodel_from_scadnano import mrdna_model_from_scadnano
+    return mrdna_model_from_scadnano(sc_file, **model_parameters)
+
 def read_vhelix(maya_file, **model_parameters):
     from .polygon_mesh import parse_maya_file, convert_maya_bases_to_segment_model
 
@@ -28,12 +36,16 @@ def read_list(infile,**model_parameters):
     from .segmentmodel_from_lists import model_from_basepair_stack_3prime
     return model_from_basepair_stack_3prime(coords, bp, stack, three_prime)
 
+def read_pandas(df,coordinate_col="r",bp_col="bp",stack_col="stack",three_prime_col="threeprime",
+                        seq_col="seq",orientation_col="orientation"):
+    from .segmentmodel_from_lists import model_from_pandas
+    return model_from_pandas(df,coordinate_col,bp_col,stack_col,three_prime_col,
+                        seq_col,orientation_col)
 
 def read_atomic_pdb(pdb_file, **model_parameters):
     from .segmentmodel_from_pdb import SegmentModelFromPdb
     return SegmentModelFromPdb(pdb_file)
 
-def read_oxdna(coordinate_file, topology_file, idealized_coordinate_file=None, **model_parameters):
-    """ Idealized coordinate file: coordinates for detecting stacks and base pairs, defaults to coordinate_file """ 
+def read_oxdna(coordinate_file, topology_file, virt2nuc=None, **model_parameters):
     from .segmentmodel_from_oxdna import mrdna_model_from_oxdna
-    return mrdna_model_from_oxdna(coordinate_file, topology_file, idealized_coordinate_file, **model_parameters)
+    return mrdna_model_from_oxdna(coordinate_file, topology_file,virt2nuc, **model_parameters)

mrdna/readers/cadnano_segments.py

@@ -9,7 +9,7 @@ from ..arbdmodel.coords import readArbdCoords, readAvgArbdCoords, rotationAboutA
 from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
 from ..model.dna_sequence import m13 as m13seq
 
-
+## Only testing on cadnano2.5
 ## TODO: separate SegmentModel from ArbdModel so multiple parts can be combined
 ## TODO: catch circular strands in "get_5prime" cadnano calls
 ## TODO: handle special motifs
@@ -17,166 +17,7 @@ from ..model.dna_sequence import m13 as m13seq
 ##   - helices that should be stacked across an empty region (crossovers from and end in the helix to another end in the helix)
 ##   - circular constructs
 
-def combineRegionLists(loHi1,loHi2,intersect=False):
-
-    """Combines two lists of (lo,hi) pairs specifying integer
-    regions a single list of regions.  """
-
-    ## Validate input
-    for l in (loHi1,loHi2):
-        ## Assert each region in lists is sorted
-        for pair in l:
-            assert(len(pair) == 2)
-            assert(pair[0] <= pair[1])
-
-    if len(loHi1) == 0:
-        if intersect:
-            return []
-        else:
-            return loHi2
-    if len(loHi2) == 0:
-        if intersect:
-            return []
-        else:
-            return loHi1
-
-    ## Break input into lists of compact regions
-    compactRegions1,compactRegions2 = [[],[]]
-    for compactRegions,loHi in zip(
-            [compactRegions1,compactRegions2],
-            [loHi1,loHi2]):
-        tmp = []
-        lastHi = loHi[0][0]-1
-        for lo,hi in loHi:
-            if lo-1 != lastHi:
-                compactRegions.append(tmp)
-                tmp = []
-            tmp.append((lo,hi))
-            lastHi = hi
-        if len(tmp) > 0:
-            compactRegions.append(tmp)
-
-    ## Build result
-    result = []
-    region = []
-    i,j = [0,0]
-    compactRegions1.append([[1e10]])
-    compactRegions2.append([[1e10]])
-    while i < len(compactRegions1)-1 or j < len(compactRegions2)-1:
-        cr1 = compactRegions1[i]
-        cr2 = compactRegions2[j]
-
-        ## initialize region
-        if len(region) == 0:
-            if cr1[0][0] <= cr2[0][0]:
-                region = cr1
-                i += 1
-                continue
-            else:
-                region = cr2
-                j += 1
-                continue
-
-        if region[-1][-1] >= cr1[0][0]:
-            region = combineCompactRegionLists(region, cr1, intersect=False)
-            i+=1
-        elif region[-1][-1] >= cr2[0][0]:
-            region = combineCompactRegionLists(region, cr2, intersect=False)
-            j+=1
-        else:
-            result.extend(region)
-            region = []
-
-    assert( len(region) > 0 )
-    result.extend(region)
-    result = sorted(result)
-
-    # print("loHi1:",loHi1)
-    # print("loHi2:",loHi2)
-    # print(result,"\n")
-
-    if intersect:
-        lo = max( [loHi1[0][0], loHi2[0][0]] )
-        hi = min( [loHi1[-1][1], loHi2[-1][1]] )
-        result = [r for r in result if r[0] >= lo and r[1] <= hi]
-
-    return result
-
-def combineCompactRegionLists(loHi1,loHi2,intersect=False):
-
-    """Combines two lists of (lo,hi) pairs specifying regions within a
-    compact integer set into a single list of regions.
-
-    examples:
-    loHi1 = [[0,4],[5,7]]
-    loHi2 = [[2,4],[5,9]]
-    out = [(0, 1), (2, 4), (5, 7), (8, 9)]
-
-    loHi1 = [[0,3],[5,7]]
-    loHi2 = [[2,4],[5,9]]
-    out = [(0, 1), (2, 3), (4, 4), (5, 7), (8, 9)]
-    """
-
-    ## Validate input
-    for l in (loHi1,loHi2):
-        ## Assert each region in lists is sorted
-        for pair in l:
-            assert(len(pair) == 2)
-            assert(pair[0] <= pair[1])
-        ## Assert lists are compact
-        for pair1,pair2 in zip(l[::2],l[1::2]):
-            assert(pair1[1]+1 == pair2[0])
-
-    if len(loHi1) == 0:
-        if intersect:
-            return []
-        else:
-            return loHi2
-    if len(loHi2) == 0:
-        if intersect:
-            return []
-        else:
-            return loHi1
-
-    ## Find the ends of the region
-    lo = min( [loHi1[0][0], loHi2[0][0]] )
-    hi = max( [loHi1[-1][1], loHi2[-1][1]] )
-
-    ## Make a list of indices where each region will be split
-    splitAfter = []
-    for l,h in loHi2:
-        if l != lo:
-            splitAfter.append(l-1)
-        if h != hi:
-            splitAfter.append(h)
-
-    for l,h in loHi1:
-        if l != lo:
-            splitAfter.append(l-1)
-        if h != hi:
-            splitAfter.append(h)
-    splitAfter = sorted(list(set(splitAfter)))
-
-    # print("splitAfter:",splitAfter)
-
-    split=[]
-    last = -2
-    for s in splitAfter:
-        split.append(s)
-        last = s
-
-    # print("split:",split)
-    returnList = [(i+1,j) if i != j else (i,j) for i,j in zip([lo-1]+split,split+[hi])]
-
-    if intersect:
-        lo = max( [loHi1[0][0], loHi2[0][0]] )
-        hi = min( [loHi1[-1][1], loHi2[-1][1]] )
-        returnList = [r for r in returnList if r[0] >= lo and r[1] <= hi]
 
-    # print("loHi1:",loHi1)
-    # print("loHi2:",loHi2)
-    # print(returnList,"\n")
-    return returnList
 
 class cadnano_part(SegmentModel):
     def __init__(self, part, 
@@ -700,6 +541,166 @@ def read_model(json_data, sequence=None, fill_sequence='T', **kwargs):
 # pynvml.nvmlShutdown()
 # gpus = [0,1,2]
 # print(gpus)
+def combineRegionLists(loHi1,loHi2,intersect=False):
+
+    """Combines two lists of (lo,hi) pairs specifying integer
+    regions a single list of regions.  """
+
+    ## Validate input
+    for l in (loHi1,loHi2):
+        ## Assert each region in lists is sorted
+        for pair in l:
+            assert(len(pair) == 2)
+            assert(pair[0] <= pair[1])
+
+    if len(loHi1) == 0:
+        if intersect:
+            return []
+        else:
+            return loHi2
+    if len(loHi2) == 0:
+        if intersect:
+            return []
+        else:
+            return loHi1
+
+    ## Break input into lists of compact regions
+    compactRegions1,compactRegions2 = [[],[]]
+    for compactRegions,loHi in zip(
+            [compactRegions1,compactRegions2],
+            [loHi1,loHi2]):
+        tmp = []
+        lastHi = loHi[0][0]-1
+        for lo,hi in loHi:
+            if lo-1 != lastHi:
+                compactRegions.append(tmp)
+                tmp = []
+            tmp.append((lo,hi))
+            lastHi = hi
+        if len(tmp) > 0:
+            compactRegions.append(tmp)
+
+    ## Build result
+    result = []
+    region = []
+    i,j = [0,0]
+    compactRegions1.append([[1e10]])
+    compactRegions2.append([[1e10]])
+    while i < len(compactRegions1)-1 or j < len(compactRegions2)-1:
+        cr1 = compactRegions1[i]
+        cr2 = compactRegions2[j]
+
+        ## initialize region
+        if len(region) == 0:
+            if cr1[0][0] <= cr2[0][0]:
+                region = cr1
+                i += 1
+                continue
+            else:
+                region = cr2
+                j += 1
+                continue
+
+        if region[-1][-1] >= cr1[0][0]:
+            region = combineCompactRegionLists(region, cr1, intersect=False)
+            i+=1
+        elif region[-1][-1] >= cr2[0][0]:
+            region = combineCompactRegionLists(region, cr2, intersect=False)
+            j+=1
+        else:
+            result.extend(region)
+            region = []
+
+    assert( len(region) > 0 )
+    result.extend(region)
+    result = sorted(result)
+
+    # print("loHi1:",loHi1)
+    # print("loHi2:",loHi2)
+    # print(result,"\n")
+
+    if intersect:
+        lo = max( [loHi1[0][0], loHi2[0][0]] )
+        hi = min( [loHi1[-1][1], loHi2[-1][1]] )
+        result = [r for r in result if r[0] >= lo and r[1] <= hi]
+
+    return result
+
+def combineCompactRegionLists(loHi1,loHi2,intersect=False):
+
+    """Combines two lists of (lo,hi) pairs specifying regions within a
+    compact integer set into a single list of regions.
+
+    examples:
+    loHi1 = [[0,4],[5,7]]
+    loHi2 = [[2,4],[5,9]]
+    out = [(0, 1), (2, 4), (5, 7), (8, 9)]
+
+    loHi1 = [[0,3],[5,7]]
+    loHi2 = [[2,4],[5,9]]
+    out = [(0, 1), (2, 3), (4, 4), (5, 7), (8, 9)]
+    """
+
+    ## Validate input
+    for l in (loHi1,loHi2):
+        ## Assert each region in lists is sorted
+        for pair in l:
+            assert(len(pair) == 2)
+            assert(pair[0] <= pair[1])
+        ## Assert lists are compact
+        for pair1,pair2 in zip(l[::2],l[1::2]):
+            assert(pair1[1]+1 == pair2[0])
+
+    if len(loHi1) == 0:
+        if intersect:
+            return []
+        else:
+            return loHi2
+    if len(loHi2) == 0:
+        if intersect:
+            return []
+        else:
+            return loHi1
+
+    ## Find the ends of the region
+    lo = min( [loHi1[0][0], loHi2[0][0]] )
+    hi = max( [loHi1[-1][1], loHi2[-1][1]] )
+
+    ## Make a list of indices where each region will be split
+    splitAfter = []
+    for l,h in loHi2:
+        if l != lo:
+            splitAfter.append(l-1)
+        if h != hi:
+            splitAfter.append(h)
+
+    for l,h in loHi1:
+        if l != lo:
+            splitAfter.append(l-1)
+        if h != hi:
+            splitAfter.append(h)
+    splitAfter = sorted(list(set(splitAfter)))
+
+    # print("splitAfter:",splitAfter)
+
+    split=[]
+    last = -2
+    for s in splitAfter:
+        split.append(s)
+        last = s
+
+    # print("split:",split)
+    returnList = [(i+1,j) if i != j else (i,j) for i,j in zip([lo-1]+split,split+[hi])]
+
+    if intersect:
+        lo = max( [loHi1[0][0], loHi2[0][0]] )
+        hi = min( [loHi1[-1][1], loHi2[-1][1]] )
+        returnList = [r for r in returnList if r[0] >= lo and r[1] <= hi]
+
+    # print("loHi1:",loHi1)
+    # print("loHi2:",loHi2)
+    # print(returnList,"\n")
+    return returnList
 
 if __name__ == '__main__':
     loHi1 = [[0,4],[5,7]]