-
Notifications
You must be signed in to change notification settings - Fork 2
/
rna_poly.py
173 lines (135 loc) · 6.65 KB
/
rna_poly.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
from sage.all import *
from sage.geometry.polyhedron.parent import Polyhedra
from sage.geometry.polyhedron.backend_ppl import Polyhedron_QQ_ppl
from sage.rings.rational_field import QQ
from sage.geometry.fan import Fan, Cone_of_fan
from collections import namedtuple
from pickle import UnpicklingError
import os.path
pickle_version = 2
class RNAPolytope(Polyhedron_QQ_ppl):
def __init__(self, points):
"""
Construct an RNAPolytope from a collection of points.
Keyword arguments:
points -- An Iterable of points, each of which has (at least) members `structure` (holding arbitrary data about the point) and `vector` (giving the coordinates of the point in a format which can be handled by `Polytope`)
"""
self._normal_fan = None
self._vertex_from_cone = None
parent = Polyhedra(QQ, len(points[0].vector))
vertex_list = [point.vector for point in points]
super(RNAPolytope, self).__init__(parent, [vertex_list, None, None], None)
structure_dict = {tuple(map(QQ, point.vector)): point.structure for point in points}
self._structures = structure_dict
def _build_normal_fan(self):
if self._normal_fan is None:
cones = [[ieq.index() for ieq in vertex.incident()] for vertex in self.vertices()]
rays = [ ieq.A() for ieq in self.inequalities() ]
fan = Fan(cones, rays, check=False, is_complete = True)
self._normal_fan = fan
if self._vertex_from_cone is None:
conedict = {}
for vertex in self.vertices():
cone = self.normal_fan().cone_containing([ieq.A() for ieq in vertex.incident()])
conedict[cone] = vertex
self._vertex_from_cone = conedict
def normal_fan(self):
"""
Return the normal fan of `self`.
"""
if not hasattr(self, "_normal_fan"):
self._build_normal_fan()
return self._normal_fan
def _build_d1_slices(self):
subspace = Polyhedron(eqns=[(-1, 0, 0, 0, 1)]) # Hyperplane d = 1
slices = set()
# Construct the slices
for cone in self.normal_fan().cones(codim=0):
coneslice = cone.polyhedron().intersection(subspace)
# If the slice does not have volume, we're done
if coneslice.dimension() < 3:
continue
# Otherwise, project to Q3
vecprojector = lambda vec: vec[0:3] # Drop the d coordinate
newverts = lambda polyslice: [vecprojector(vert) for vert in polyslice.vertices()] # Projects the defining vertices of a slice into Q3
newrays = lambda polyslice: [vecprojector(ray) for ray in polyslice.rays()] # Projects the defining rays of a slice into Q3
polyprojector = lambda polyslice: Polyhedron(vertices = newverts(polyslice), rays = newrays(polyslice)) # Projects a polytope from Q4 to Q3 using the helper functions we just defined
# We can then apply polyprojector to our slice to move it down to Q3
projslice = polyprojector(coneslice)
# For future use, we'll store the original cone and the slice in Q4 as a member of this object
projslice.original_cone = cone
projslice.original_slice = coneslice
# As well as the original vertex and structure associated to those cones
projslice.original_vertex = self.vertex_from_cone(cone)
projslice.original_structure = self[projslice.original_vertex]
# And, finally, we add this slice to our list
slices.add(projslice)
# Once we've finished the loop, all the slices have been processed, so we store the results
self._d1_slices = slices
def d1_slices(self):
"""
Return the d=1 slices of self as Polyhedra in Q3
"""
if not hasattr(self, "_d1_slices"):
self._build_d1_slices()
return self._d1_slices
def vertex_from_cone(self, cone):
assert cone in self.normal_fan().cones(self.dim())
return self._vertex_from_cone[cone]
def __getitem__(self, key):
try:
processed_key = tuple(map(QQ, key)) # Try to cast the input to a tuple of rationals
return self._structures[processed_key] # Cast the input to a tuple of rationals
except TypeError:
return None # Bad key!
@classmethod
def construct_from_file(cls, polyfile):
"""
Construct an RNAPolytope from a file in polyfile format.
Keyword arguments:
polyfile -- A file representing the polytope structure
"""
filebase = os.path.splitext(polyfile)[0]
# Construct the polytope
try:
# If the polytope is available in a pickle, we should use that
thepoly = load(filebase) # Attempt to load the polytope from a pickle file
# If the pickled polytope is obsolete, however, we need to rebuild it
if cls.poly_is_obsolete(thepoly):
raise UnpicklingError
except (IOError, UnpicklingError, AssertionError):
# In any of these exception cases, the load failed, so we generate the polytope from the points
# Read the point data from the specified file
points = []
with open(polyfile) as f:
for line in f:
newline = line.split('#')[0] # Only use content before the comment symbol
if newline.strip() != '':
elts = newline.split()
structure = elts[1]
multiloops = QQ(elts[2])
unpaired = QQ(elts[3])
branches = QQ(elts[4])
w = QQ(elts[5])
energy = QQ(elts[6])
coords = (multiloops, unpaired, branches, w)
newpoint = namedtuple('RNApoint', ['structure', 'vector', 'energy'])
newpoint.vector = coords
newpoint.structure = structure
newpoint.energy = energy
points.append(newpoint)
thepoly = RNAPolytope(points)
thepoly._pickle_version = pickle_version
thepoly._build_normal_fan()
thepoly._build_d1_slices()
thepoly.dump(filebase, -1)
return thepoly
@classmethod
def poly_is_obsolete(cls, poly):
"""
Test whether the polytope object is obsolete, in which case it should be upgraded
"""
if (not hasattr(poly, "_pickle_version") or poly._pickle_version < pickle_version):
return true
else:
return false