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