CoCalc Shared Filesresearch / biomath / polytope_plots / rna_poly.py
Authors: Fidel Barrera-Cruz, William A. Stein
Views : 11
1from sage.all import *
2from sage.geometry.polyhedron.parent import Polyhedra
3from sage.geometry.polyhedron.backend_ppl import Polyhedron_QQ_ppl
4from sage.rings.rational_field import QQ
5from sage.geometry.fan import Fan, Cone_of_fan
6from collections import namedtuple
7from pickle import UnpicklingError
8import os.path
9
10pickle_version = 2
11
12class 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
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