1 3
from __future__ import division
2 3
from itertools import product
3 3
import collections
4

5 3
import numpy as np
6 3
from scipy.stats import multivariate_normal as mvn
7 3
from scipy.stats import norm
8 3
import scipy.spatial as spatial
9 3
import matplotlib.pyplot as plt
10

11 3
from sklearn.base import TransformerMixin
12

13 3
__all__ = ["PersImage"]
14

15 3
class PersImage(TransformerMixin):
16
    """ Initialize a persistence image generator.
17

18
    Parameters
19
    -----------
20

21
    pixels : pair of ints like (int, int)
22
        Tuple representing number of pixels in return image along x and y axis.
23
    spread : float
24
        Standard deviation of gaussian kernel
25
    specs : dict
26
        Parameters for shape of image with respect to diagram domain. This is used if you would like images to have a particular range. Shaped like 
27
        ::
28
        
29
            {
30
                "maxBD": float,
31
                "minBD": float
32
            }
33

34
    kernel_type : string or ...
35
        TODO: Implement this feature.
36
        Determine which type of kernel used in the convolution, or pass in custom kernel. Currently only implements Gaussian.
37
    weighting_type : string or ...
38
        TODO: Implement this feature.
39
        Determine which type of weighting function used, or pass in custom weighting function.
40
        Currently only implements linear weighting.
41

42

43
    Usage
44
    ------
45

46

47
    """
48

49 3
    def __init__(
50
        self,
51
        pixels=(20, 20),
52
        spread=None,
53
        specs=None,
54
        kernel_type="gaussian",
55
        weighting_type="linear",
56
        verbose=True,
57
    ):
58

59 3
        self.specs = specs
60 3
        self.kernel_type = kernel_type
61 3
        self.weighting_type = weighting_type
62 3
        self.spread = spread
63 3
        self.nx, self.ny = pixels
64

65 3
        if verbose:
66 3
            print(
67
                'PersImage(pixels={}, spread={}, specs={}, kernel_type="{}", weighting_type="{}")'.format(
68
                    pixels, spread, specs, kernel_type, weighting_type
69
                )
70
            )
71

72 3
    def transform(self, diagrams):
73
        """ Convert diagram or list of diagrams to a persistence image.
74

75
        Parameters
76
        -----------
77

78
        diagrams : list of or singleton diagram, list of pairs. [(birth, death)]
79
            Persistence diagrams to be converted to persistence images. It is assumed they are in (birth, death) format. Can input a list of diagrams or a single diagram.
80

81
        """
82
        # if diagram is empty, return empty image
83 3
        if len(diagrams) == 0:
84 3
            return np.zeros((self.nx, self.ny))
85
        # if first entry of first entry is not iterable, then diagrams is singular and we need to make it a list of diagrams
86 3
        try:
87 3
            singular = not isinstance(diagrams[0][0], collections.Iterable)
88 3
        except IndexError:
89 3
            singular = False
90

91 3
        if singular:
92 3
            diagrams = [diagrams]
93

94 3
        dgs = [np.copy(diagram) for diagram in diagrams]
95 3
        landscapes = [PersImage.to_landscape(dg) for dg in dgs]
96

97 3
        if not self.specs:
98 3
            self.specs = {
99
                "maxBD": np.max([np.max(np.vstack((landscape, np.zeros((1, 2))))) 
100
                                 for landscape in landscapes] + [0]),
101
                "minBD": np.min([np.min(np.vstack((landscape, np.zeros((1, 2))))) 
102
                                 for landscape in landscapes] + [0]),
103
            }
104 3
        imgs = [self._transform(dgm) for dgm in landscapes]
105

106
        # Make sure we return one item.
107 3
        if singular:
108 3
            imgs = imgs[0]
109

110 3
        return imgs
111

112 3
    def _transform(self, landscape):
113
        # Define an NxN grid over our landscape
114 3
        maxBD = self.specs["maxBD"]
115 3
        minBD = min(self.specs["minBD"], 0)  # at least show 0, maybe lower
116

117
        # Same bins in x and y axis
118 3
        dx = maxBD / (self.ny)
119 3
        xs_lower = np.linspace(minBD, maxBD, self.nx)
120 3
        xs_upper = np.linspace(minBD, maxBD, self.nx) + dx
121

122 3
        ys_lower = np.linspace(0, maxBD, self.ny)
123 3
        ys_upper = np.linspace(0, maxBD, self.ny) + dx
124

125 3
        weighting = self.weighting(landscape)
126

127
        # Define zeros
128 3
        img = np.zeros((self.nx, self.ny))
129

130
        # Implement this as a `summed-area table` - it'll be way faster
131 3
        spread = self.spread if self.spread else dx
132 3
        for point in landscape:
133 3
            x_smooth = norm.cdf(xs_upper, point[0], spread) - norm.cdf(
134
                xs_lower, point[0], spread
135
            )
136 3
            y_smooth = norm.cdf(ys_upper, point[1], spread) - norm.cdf(
137
                ys_lower, point[1], spread
138
            )
139 3
            img += np.outer(x_smooth, y_smooth) * weighting(point)
140 3
        img = img.T[::-1]
141 3
        return img
142

143 3
    def weighting(self, landscape=None):
144
        """ Define a weighting function, 
145
                for stability results to hold, the function must be 0 at y=0.    
146
        """
147

148
        # TODO: Implement a logistic function
149
        # TODO: use self.weighting_type to choose function
150

151 3
        if landscape is not None:
152 3
            if len(landscape) > 0:
153 3
                maxy = np.max(landscape[:, 1])
154
            else: 
155 3
                maxy = 1
156

157 3
        def linear(interval):
158
            # linear function of y such that f(0) = 0 and f(max(y)) = 1
159 3
            d = interval[1]
160 3
            return (1 / maxy) * d if landscape is not None else d
161

162 3
        def pw_linear(interval):
163
            """ This is the function defined as w_b(t) in the original PI paper
164

165
                Take b to be maxy/self.ny to effectively zero out the bottom pixel row
166
            """
167

168 0
            t = interval[1]
169 0
            b = maxy / self.ny
170

171 3
            if t <= 0:
172 0
                return 0
173 3
            if 0 < t < b:
174 0
                return t / b
175 3
            if b <= t:
176 0
                return 1
177

178 3
        return linear
179

180 3
    def kernel(self, spread=1):
181
        """ This will return whatever kind of kernel we want to use.
182
            Must have signature (ndarray size NxM, ndarray size 1xM) -> ndarray size Nx1
183
        """
184
        # TODO: use self.kernel_type to choose function
185

186 3
        def gaussian(data, pixel):
187 3
            return mvn.pdf(data, mean=pixel, cov=spread)
188

189 3
        return gaussian
190

191 3
    @staticmethod
192 1
    def to_landscape(diagram):
193
        """ Convert a diagram to a landscape
194
            (b,d) -> (b, d-b)
195
        """
196 3
        diagram[:, 1] -= diagram[:, 0]
197

198 3
        return diagram
199

200 3
    def show(self, imgs, ax=None):
201
        """ Visualize the persistence image
202

203
        """
204

205 0
        ax = ax or plt.gca()
206

207 3
        if type(imgs) is not list:
208 0
            imgs = [imgs]
209

210 3
        for i, img in enumerate(imgs):
211 0
            ax.imshow(img, cmap=plt.get_cmap("plasma"))
212 0
            ax.axis("off")

Read our documentation on viewing source code .

Loading