@@ -2,29 +2,37 @@
Loading
2 2
from itertools import product
3 3
import collections
4 4
5 +
import copy
5 6
import numpy as np
6 -
from scipy.stats import multivariate_normal as mvn
7 7
from scipy.stats import norm
8 8
import scipy.spatial as spatial
9 9
import matplotlib.pyplot as plt
10 -
10 +
from persim import images_kernels
11 +
from persim import images_weights
12 +
from joblib import Parallel, delayed
13 +
from deprecated.sphinx import deprecated
11 14
from sklearn.base import TransformerMixin
15 +
from matplotlib.collections import LineCollection
16 +
from scipy.stats import multivariate_normal as mvn
12 17
13 -
__all__ = ["PersImage"]
18 +
__all__ = ["PersImage", "PersistenceImager"]
14 19
20 +
@deprecated(
21 +
    reason="""Replaced with the class :class:`persim.PersistenceImager`.""",
22 +
    version="0.1.5",
23 +
)
15 24
class PersImage(TransformerMixin):
16 25
    """ Initialize a persistence image generator.
17 -
18 -
    Parameters
19 -
    -----------
20 -
26 +
    
27 +
Parameters
28 +
----------
29 +
    
21 30
    pixels : pair of ints like (int, int)
22 31
        Tuple representing number of pixels in return image along x and y axis.
23 32
    spread : float
24 -
        Standard deviation of gaussian kernel
33 +
        Standard deviation of gaussian kernel.
25 34
    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 -
        ::
35 +
        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::
28 36
        
29 37
            {
30 38
                "maxBD": float,
@@ -38,12 +46,6 @@
Loading
38 46
        TODO: Implement this feature.
39 47
        Determine which type of weighting function used, or pass in custom weighting function.
40 48
        Currently only implements linear weighting.
41 -
42 -
43 -
    Usage
44 -
    ------
45 -
46 -
47 49
    """
48 50
49 51
    def __init__(
@@ -210,3 +212,597 @@
Loading
210 212
        for i, img in enumerate(imgs):
211 213
            ax.imshow(img, cmap=plt.get_cmap("plasma"))
212 214
            ax.axis("off")
215 +
216 +
217 +
class PersistenceImager(TransformerMixin):
218 +
    """Transformer which converts persistence diagrams into persistence images.
219 +
220 +
    Parameters
221 +
    ----------
222 +
    birth_range : pair of floats
223 +
        Range of persistence pair birth values covered by the persistence image (default: (0.0, 1.0)).
224 +
    pers_range : pair of floats
225 +
        Range of persistence pair persistence (death-birth) values covered by the persistence image (default: (0.0, 1.0)).
226 +
    pixel_size : float
227 +
        Dimensions of each square pixel (default: 0.2).
228 +
    weight : callable or str in ['persistence', 'linear_ramp']
229 +
        Function which weights the birth-persistence plane (default: 'persistence').
230 +
    weight_params : dict
231 +
        Arguments needed to specify the weight function (default: {'n': 1.0}).
232 +
    kernel : callable or str in ['gaussian', 'uniform']
233 +
        Cumulative distribution function defining the kernel (default: 'gaussian').
234 +
    kernel_params : dict
235 +
        Arguments needed to specify the kernel function (default: {'sigma': [[1.0, 0.0], [0.0, 1.0]]}).
236 +
237 +
238 +
    Example
239 +
    -------
240 +
    First instantiate a PersistenceImager() object::
241 +
242 +
        > from persim import PersistenceImager
243 +
        > pimgr = PersistenceImager(pixel_size=0.2, birth_range=(0,1))
244 +
245 +
246 +
    Printing a PersistenceImager() object will print its hyperparameters::
247 +
    
248 +
        > print(pimgr)
249 +
        
250 +
        PersistenceImager(birth_range=(0.0, 1.0), pers_range=(0.0, 1.0), pixel_size=0.2, weight=persistence, weight_params={'n': 1.0}, kernel=gaussian, kernel_params={'sigma': [[1.0, 0.0], [0.0, 1.0]]})
251 +
252 +
253 +
    PersistenceImager() attributes can be adjusted at or after instantiation. Updating attributes of a PersistenceImager() object will automatically update all other dependent attributes::
254 +
    
255 +
        > pimgr.pixel_size = 0.1
256 +
        > pimgr.birth_range = (0, 2)
257 +
        > print(pimgr)
258 +
        > print(pimgr.resolution)
259 +
    
260 +
        PersistenceImager(birth_range=(0.0, 2.0), pers_range=(0.0, 1.0), pixel_size=0.1, weight=persistence, weight_params={'n': 1.0}, kernel=gaussian, kernel_params={'sigma': [[1.0, 0.0], [0.0, 1.0]]})
261 +
        (20, 10)
262 +
263 +
264 +
    The `fit()` method can be called on one or more (-,2) numpy.ndarrays to automatically determine the miniumum birth and persistence ranges needed to capture all persistence pairs. The ranges and resolution are automatically adjusted to accomodate the specified pixel size. The option `skew=True` specifies that the diagram is currently in birth-death coordinates and must first be transformed to birth-persistence coordinates::
265 +
    
266 +
        > import numpy as np
267 +
        > pimgr = PersistenceImager(pixel_size=0.5)
268 +
        > pdgms = [np.array([[0.5, 0.8], [0.7, 2.2], [2.5, 4.0]]),
269 +
                   np.array([[0.1, 0.2], [3.1, 3.3], [1.6, 2.9]]),
270 +
                   np.array([[0.2, 1.5], [0.4, 0.6], [0.2, 2.6]])]
271 +
        > pimgr.fit(pdgms, skew=True)
272 +
        > print(pimgr)
273 +
        > print(pimgr.resolution)
274 +
    
275 +
        PersistenceImager(birth_range=(0.1, 3.1), pers_range=(-8.326672684688674e-17, 2.5), pixel_size=0.5, weight=persistence, weight_params={'n': 1.0}, kernel=gaussian, kernel_params={'sigma': [[1.0, 0.0], [0.0, 1.0]]})
276 +
        (6, 5)
277 +
278 +
279 +
    The `transform()` method can then be called on one or more (-,2) numpy.ndarrays to generate persistence images from diagrams. The option `skew=True` specifies that the diagrams are currently in birth-death coordinates and must first be transformed to birth-persistence coordinates::
280 +
281 +
        > pimgs = pimgr.transform(pdgms, skew=True)
282 +
        > pimgs[0]
283 +
    
284 +
        array([[0.03999068, 0.05688393, 0.06672051, 0.06341749, 0.04820814],
285 +
               [0.04506697, 0.06556791, 0.07809764, 0.07495246, 0.05730671],
286 +
               [0.04454486, 0.06674611, 0.08104366, 0.07869919, 0.06058808],
287 +
               [0.04113063, 0.0636504 , 0.07884635, 0.07747833, 0.06005714],
288 +
               [0.03625436, 0.05757744, 0.07242608, 0.07180125, 0.05593626],
289 +
               [0.02922239, 0.04712024, 0.05979033, 0.05956698, 0.04653357]]) 
290 +
291 +
292 +
    Notes
293 +
    -----
294 +
    [1] Adams et. al., "Persistence Images: A Stable Vector Representation of Persistent Homology," Journal of Machine Learning Research, vol. 18, pp. 1-35, 2017. http://www.jmlr.org/papers/volume18/16-337/16-337.pdf
295 +
    """
296 +
    
297 +
    def __init__(self, birth_range=None, pers_range=None, pixel_size=None, weight=None, weight_params=None, kernel=None, kernel_params=None):
298 +
        """ PersistenceImager constructor method
299 +
        """       
300 +
        # set defaults
301 +
        if birth_range is None:
302 +
            birth_range = (0.0, 1.0)
303 +
        if pers_range is None:
304 +
            pers_range = (0.0, 1.0)
305 +
        if pixel_size is None:
306 +
            pixel_size = 0.2
307 +
        if weight is None:
308 +
            weight = images_weights.persistence
309 +
        if kernel is None:
310 +
            kernel = images_kernels.gaussian
311 +
        if weight_params is None:
312 +
            weight_params = {'n': 1.0}
313 +
        if kernel_params is None:
314 +
            kernel_params = {'sigma': [[1.0, 0.0], [0.0, 1.0]]}
315 +
        
316 +
        # validate parameters
317 +
        self._validate_parameters(birth_range=birth_range, pers_range=pers_range, pixel_size=pixel_size, weight=weight, weight_params=weight_params, kernel=kernel, kernel_params=kernel_params)
318 +
        
319 +
        self.weight, self.kernel = self._ensure_callable(weight=weight, kernel=kernel)
320 +
        self.weight_params = weight_params
321 +
        self.kernel_params = kernel_params
322 +
        self._pixel_size = pixel_size
323 +
        self._birth_range = birth_range
324 +
        self._pers_range = pers_range
325 +
        self._width = birth_range[1] - birth_range[0]
326 +
        self._height = pers_range[1] - pers_range[0]
327 +
        self._resolution = (int(self._width / self._pixel_size), int(self._height / self._pixel_size))
328 +
        self._create_mesh()
329 +
        
330 +
    @property
331 +
    def width(self):
332 +
        """
333 +
        Persistence image width.
334 +
        
335 +
        Returns
336 +
        -------
337 +
        width : float
338 +
            The width of the region of the birth-persistence plane covered by the persistence image in birth units.
339 +
        """
340 +
        return self._width
341 +
342 +
    @property
343 +
    def height(self):
344 +
        """
345 +
        Persistence image height.
346 +
        
347 +
        Returns
348 +
        -------
349 +
        height : float
350 +
            The height of the region of the birth-persistence plane covered by the persistence image in persistence units.
351 +
        """
352 +
        return self._height
353 +
354 +
    @property
355 +
    def resolution(self):
356 +
        """
357 +
        Persistence image resolution.
358 +
        
359 +
        Returns
360 +
        -------
361 +
        resolution : pair of ints (width, height)
362 +
            The number of pixels along each dimension of the persistence image, determined by the birth and persistence ranges and the pixel size.
363 +
        """
364 +
        return self._resolution
365 +
366 +
    @property
367 +
    def pixel_size(self):
368 +
        """
369 +
        Persistence image square pixel dimensions.
370 +
        
371 +
        Returns
372 +
        -------
373 +
        pixel_size : float
374 +
            The width (and height) in birth/persistence units of each square pixel in the persistence image. 
375 +
        """
376 +
        return self._pixel_size
377 +
378 +
    @pixel_size.setter
379 +
    def pixel_size(self, val):
380 +
        self._pixel_size = val
381 +
        self._width = int(np.ceil((self.birth_range[1] - self.birth_range[0]) / self.pixel_size)) * self.pixel_size
382 +
        self._height = int(np.ceil((self.pers_range[1] - self.pers_range[0]) / self.pixel_size)) * self.pixel_size
383 +
        self._resolution = (int(self.width / self.pixel_size), int(self.height / self.pixel_size))
384 +
        self._create_mesh()
385 +
386 +
    @property
387 +
    def birth_range(self):
388 +
        """
389 +
        Range of birth values covered by the persistence image.
390 +
        
391 +
        Returns
392 +
        -------
393 +
        birth_range : pair of floats (min. birth, max. birth)
394 +
            The minimum and maximum birth values covered by the persistence image.
395 +
        """
396 +
        return self._birth_range
397 +
398 +
    @birth_range.setter
399 +
    def birth_range(self, val):
400 +
        self._birth_range = val
401 +
        self._width = int(np.ceil((self.birth_range[1] - self.birth_range[0]) / self.pixel_size)) * self._pixel_size
402 +
        self._resolution = (int(self.width / self.pixel_size), int(self.height / self.pixel_size))
403 +
        self._create_mesh()
404 +
405 +
    @property
406 +
    def pers_range(self):
407 +
        """
408 +
        Range of persistence values covered by the persistence image.
409 +
        
410 +
        Returns
411 +
        -------
412 +
        pers_range : pair of floats (min. persistence, max. persistence)
413 +
            The minimum and maximum persistence values covered by the persistence image.
414 +
        """
415 +
        return self._pers_range
416 +
417 +
    @pers_range.setter
418 +
    def pers_range(self, val):
419 +
        self._pers_range = val
420 +
        self._height = int(np.ceil((self.pers_range[1] - self.pers_range[0]) / self.pixel_size)) * self._pixel_size
421 +
        self._resolution = (int(self.width / self.pixel_size), int(self.height / self.pixel_size))
422 +
        self._create_mesh()
423 +
        
424 +
    def __repr__(self):
425 +
        import pprint as pp
426 +
        params = tuple([self.birth_range, self.pers_range, self.pixel_size, self.weight.__name__, pp.pformat(self.weight_params), self.kernel.__name__, pp.pformat(self.kernel_params)])
427 +
        repr_str = 'PersistenceImager(birth_range=%s, pers_range=%s, pixel_size=%s, weight=%s, weight_params=%s, kernel=%s, kernel_params=%s)' % params
428 +
        return repr_str
429 +
    
430 +
    def _ensure_callable(self, weight=None, kernel=None):
431 +
        valid_weights_dict = {'persistence': images_weights.persistence, 'linear_ramp': images_weights.linear_ramp}
432 +
        valid_kernels_dict = {'gaussian': images_kernels.gaussian, 'uniform': images_kernels.uniform}
433 +
        
434 +
        if isinstance(weight, str):
435 +
            weight = valid_weights_dict[weight]
436 +
            
437 +
        if isinstance(kernel, str):
438 +
            kernel = valid_kernels_dict[kernel]
439 +
            
440 +
        return weight, kernel
441 +
    
442 +
    def _validate_parameters(self, birth_range=None, pers_range=None, pixel_size=None, weight=None, weight_params=None, kernel=None, kernel_params=None):
443 +
        valid_weights = ['persistence', 'linear_ramp']
444 +
        valid_kernels = ['gaussian', 'uniform']
445 +
        
446 +
        # validate birth_range
447 +
        if isinstance(birth_range, tuple):
448 +
            if len(birth_range) != 2:
449 +
                raise ValueError("birth_range must be a pair: (min. birth, max. birth)")
450 +
            elif (not isinstance(birth_range[0], (int, float))) or (not isinstance(birth_range[1], (int, float))):
451 +
                raise ValueError("birth_range must be a pair of numbers: (min. birth, max. birth)")
452 +
        else:
453 +
            raise ValueError("birth_range must be a tuple")
454 +
        
455 +
        # validate pers_range
456 +
        if isinstance(pers_range, tuple):
457 +
            if len(pers_range) != 2:
458 +
                raise ValueError("pers_range must be a pair: (min. persistence, max. persistence)")
459 +
            elif (not isinstance(pers_range[0], (int, float))) or (not isinstance(pers_range[1], (int, float))):
460 +
                raise ValueError("pers_range must be a pair of numbers: (min. persistence, max. persistence)")
461 +
        else:
462 +
            raise ValueError("pers_range must be a tuple")
463 +
            
464 +
        # validate pixel_size
465 +
        if not isinstance(pixel_size, (int, float)):
466 +
            raise ValueError("pixel_size must be an int or float")            
467 +
        
468 +
        # validate weight
469 +
        if not callable(weight):
470 +
            if isinstance(weight, str):
471 +
                if weight not in ['persistence', 'linear_ramp']:
472 +
                    raise ValueError("weight must be callable or a str in %s" %valid_weights)
473 +
            else:
474 +
                raise ValueError("weight must be callable or a str")
475 +
        
476 +
        # validate weight_params
477 +
        if not isinstance(weight_params, dict):
478 +
            raise ValueError("weight_params must be a dict")
479 +
       
480 +
        # validate kernel
481 +
        if not callable(kernel):
482 +
            if isinstance(kernel, str):
483 +
                if kernel not in valid_kernels:
484 +
                    raise ValueError("kernel must be callable or a str in %s" %valid_kernels)
485 +
            else:
486 +
                raise ValueError("kernel must be callable or a str")
487 +
        
488 +
        # validate kernel_params
489 +
        if not isinstance(kernel_params, dict):
490 +
            raise ValueError("kernel_params must be a dict")
491 +
            
492 +
    def _create_mesh(self):
493 +
        # padding around specified image ranges as a result of incommensurable ranges and pixel width
494 +
        db = self._width - (self._birth_range[1] - self._birth_range[0])
495 +
        dp = self._height - (self._pers_range[1] - self._pers_range[0])
496 +
497 +
        # adjust image ranges to accommodate incommensurable ranges and pixel width
498 +
        self._birth_range = (self._birth_range[0] - db / 2, self._birth_range[1] + db / 2)
499 +
        self._pers_range = (self._pers_range[0] - dp / 2, self._pers_range[1] + dp / 2)
500 +
        
501 +
        # construct linear spaces defining pixel locations
502 +
        self._bpnts = np.linspace(self._birth_range[0], self._birth_range[1] + self._pixel_size,
503 +
                                           self._resolution[0] + 1, endpoint=False, dtype=np.float64)
504 +
        self._ppnts = np.linspace(self._pers_range[0], self._pers_range[1] + self._pixel_size,
505 +
                                           self._resolution[1] + 1, endpoint=False, dtype=np.float64)
506 +
507 +
    def fit(self, pers_dgms, skew=True):
508 +
        """ Choose persistence image range parameters which minimally enclose all persistence pairs across one or more persistence diagrams.
509 +
        
510 +
        Parameters
511 +
        ----------
512 +
        pers_dgms : one or an iterable of (-,2) numpy.ndarrays
513 +
            Collection of one or more persistence diagrams.
514 +
        skew : boolean 
515 +
            Flag indicating if diagram(s) need to first be converted to birth-persistence coordinates (default: True).
516 +
        """
517 +
        min_birth = np.Inf
518 +
        max_birth = -np.Inf
519 +
        min_pers = np.Inf
520 +
        max_pers = -np.Inf
521 +
        
522 +
        # convert to a list of diagrams if necessary 
523 +
        pers_dgms, singular = self._ensure_iterable(pers_dgms)
524 +
        
525 +
        # loop over diagrams to determine the maximum extent of the pairs contained in the birth-persistence plane
526 +
        for pers_dgm in pers_dgms:
527 +
            pers_dgm = np.copy(pers_dgm)
528 +
            if skew:
529 +
                pers_dgm[:, 1] = pers_dgm[:, 1] - pers_dgm[:, 0]
530 +
531 +
            min_b, min_p = pers_dgm.min(axis=0)
532 +
            max_b, max_p = pers_dgm.max(axis=0)
533 +
534 +
            if min_b < min_birth:
535 +
                min_birth = min_b
536 +
537 +
            if min_p < min_pers:
538 +
                min_pers = min_p
539 +
540 +
            if max_b > max_birth:
541 +
                max_birth = max_b
542 +
543 +
            if max_p > max_pers:
544 +
                max_pers = max_p
545 +
546 +
        self.birth_range = (min_birth, max_birth)
547 +
        self.pers_range = (min_pers, max_pers)
548 +
549 +
    def transform(self, pers_dgms, skew=True, n_jobs=None):
550 +
        """ Transform a persistence diagram or an iterable containing a collection of persistence diagrams into 
551 +
        persistence images.
552 +
        
553 +
        Parameters
554 +
        ----------
555 +
        pers_dgms : one or an iterable of (-,2) numpy.ndarrays
556 +
            Collection of one or more persistence diagrams.
557 +
        skew : boolean 
558 +
            Flag indicating if diagram(s) need to first be converted to birth-persistence coordinates (default: True).
559 +
        n_jobs : int
560 +
            Number of cores to use to transform a collection of persistence diagrams into persistence images (default: None, uses a single core).
561 +
562 +
        Returns
563 +
        -------
564 +
        list
565 +
            Collection of numpy.ndarrays encoding the persistence images in the same order as pers_dgms.
566 +
        """
567 +
        if n_jobs is not None:
568 +
            parallelize = True
569 +
        else:
570 +
            parallelize = False
571 +
            
572 +
        # if diagram is empty, return empty image
573 +
        if len(pers_dgms) == 0:
574 +
            return np.zeros(self.resolution)
575 +
        
576 +
        # convert to a list of diagrams if necessary 
577 +
        pers_dgms, singular = self._ensure_iterable(pers_dgms)
578 +
        
579 +
        if parallelize:
580 +
            pers_imgs = Parallel(n_jobs=n_jobs)(delayed(_transform)(pers_dgm, skew, self.resolution, self.weight, self.weight_params, self.kernel, self.kernel_params, self._bpnts, self._ppnts) for pers_dgm in pers_dgms)
581 +
        else:
582 +
            pers_imgs = [_transform(pers_dgm, skew=skew, resolution=self.resolution, weight=self.weight, weight_params=self.weight_params, kernel=self.kernel, kernel_params=self.kernel_params, _bpnts=self._bpnts, _ppnts=self._ppnts) for pers_dgm in pers_dgms]
583 +
        
584 +
        if singular:
585 +
            pers_imgs = pers_imgs[0]
586 +
        
587 +
        return pers_imgs
588 +
589 +
    def fit_transform(self, pers_dgms, skew=True):
590 +
        """ Choose persistence image range parameters which minimally enclose all persistence pairs across one or more persistence diagrams and transform the persistence diagrams into persistence images.
591 +
        
592 +
        Parameters
593 +
        ----------
594 +
        pers_dgms : one or an iterable of (-,2) numpy.ndarray
595 +
            Collection of one or more persistence diagrams.
596 +
        skew : boolean 
597 +
            Flag indicating if diagram(s) need to first be converted to birth-persistence coordinates (default: True).
598 +
        
599 +
600 +
        Returns
601 +
        -------
602 +
        list
603 +
            Collection of numpy.ndarrays encoding the persistence images in the same order as pers_dgms.
604 +
        """
605 +
        pers_dgms = copy.deepcopy(pers_dgms)
606 +
607 +
        # fit imager parameters
608 +
        self.fit(pers_dgms, skew=skew)
609 +
610 +
        # transform diagrams to images
611 +
        pers_imgs = self.transform(pers_dgms, skew=skew)
612 +
613 +
        return pers_imgs
614 +
 
615 +
    def _ensure_iterable(self, pers_dgms):
616 +
        # if first entry of first entry is not iterable, then diagrams is singular and we need to make it a list of diagrams
617 +
        try:
618 +
            singular = not isinstance(pers_dgms[0][0], collections.Iterable)
619 +
        except IndexError:
620 +
            singular = False
621 +
622 +
        if singular:
623 +
            pers_dgms = [pers_dgms]
624 +
            
625 +
        return pers_dgms, singular
626 +
627 +
    def plot_diagram(self, pers_dgm, skew=True, ax=None, out_file=None):
628 +
        """ Plot a persistence diagram.
629 +
        
630 +
        Parameters
631 +
        ----------
632 +
        pers_dgm : (-,2) numpy.ndarray
633 +
            A persistence diagram.
634 +
        skew : boolean 
635 +
            Flag indicating if diagram needs to first be converted to birth-persistence coordinates (default: True).
636 +
        ax : matplotlib.Axes
637 +
            Instance of a matplotlib.Axes object in which to plot (default: None, generates a new figure)
638 +
        out_file : str
639 +
            Path and file name including extension to save the figure (default: None, figure not saved).
640 +
641 +
        Returns
642 +
        -------
643 +
        matplotlib.Axes
644 +
            The matplotlib.Axes which contains the persistence diagram
645 +
        """
646 +
        pers_dgm = np.copy(pers_dgm)
647 +
648 +
        if skew:
649 +
            pers_dgm[:, 1] = pers_dgm[:, 1] - pers_dgm[:, 0]
650 +
            ylabel = 'persistence'
651 +
        else:
652 +
            ylabel = 'death'
653 +
654 +
        # setup plot range
655 +
        plot_buff_frac = 0.05
656 +
        bmin = np.min((np.min(pers_dgm[:, 0]), np.min(self._bpnts)))
657 +
        bmax = np.max((np.max(pers_dgm[:, 0]), np.max(self._bpnts)))
658 +
        b_plot_buff = (bmax - bmin) * plot_buff_frac
659 +
        bmin -= b_plot_buff
660 +
        bmax += b_plot_buff
661 +
662 +
        pmin = np.min((np.min(pers_dgm[:, 1]), np.min(self._ppnts)))
663 +
        pmax = np.max((np.max(pers_dgm[:, 1]), np.max(self._ppnts)))
664 +
        p_plot_buff = (pmax - pmin) * plot_buff_frac
665 +
        pmin -= p_plot_buff
666 +
        pmax += p_plot_buff
667 +
668 +
        ax = ax or plt.gca()
669 +
        ax.set_xlim(bmin, bmax)
670 +
        ax.set_ylim(pmin, pmax)
671 +
672 +
        # compute reasonable line width for pixel overlay (initially 1/50th of the width of a pixel)
673 +
        linewidth = (1/50 * self.pixel_size) * 72 * plt.gcf().bbox_inches.width * ax.get_position().width / \
674 +
                    np.min((bmax - bmin, pmax - pmin))
675 +
676 +
        # plot the persistence image grid
677 +
        if skew:
678 +
            hlines = np.column_stack(np.broadcast_arrays(self._bpnts[0], self._ppnts, self._bpnts[-1], self._ppnts))
679 +
            vlines = np.column_stack(np.broadcast_arrays(self._bpnts, self._ppnts[0], self._bpnts, self._ppnts[-1]))
680 +
            lines = np.concatenate([hlines, vlines]).reshape(-1, 2, 2)
681 +
            line_collection = LineCollection(lines, color='black', linewidths=linewidth)
682 +
            ax.add_collection(line_collection)       
683 +
684 +
        # plot persistence diagram
685 +
        ax.scatter(pers_dgm[:, 0], pers_dgm[:, 1])
686 +
687 +
        # plot diagonal if necessary
688 +
        if not skew:
689 +
            min_diag = np.min((np.min(ax.get_xlim()), np.min(ax.get_ylim())))
690 +
            max_diag = np.min((np.max(ax.get_xlim()), np.max(ax.get_ylim())))
691 +
            ax.plot([min_diag, max_diag], [min_diag, max_diag])
692 +
693 +
        # fix and label axes
694 +
        ax.set_aspect('equal')
695 +
        ax.set_xlabel('birth')
696 +
        ax.set_ylabel(ylabel)
697 +
698 +
        # optionally save figure
699 +
        if out_file:
700 +
            plt.savefig(out_file, bbox_inches='tight')
701 +
        
702 +
        return ax
703 +
704 +
705 +
    def plot_image(self, pers_img, ax=None, out_file=None):
706 +
        """ Plot a persistence image.
707 +
        
708 +
        Parameters
709 +
        ----------
710 +
        pers_img : (M,N) numpy.ndarray
711 +
            A persistence image, as output by PersistenceImager().transform()
712 +
        ax : matplotlib.Axes
713 +
            Instance of a matplotlib.Axes object in which to plot (default: None, generates a new figure)
714 +
        out_file : str
715 +
            Path and file name including extension to save the figure (default: None, figure not saved).
716 +
717 +
        Returns
718 +
        -------
719 +
        matplotlib.Axes
720 +
            The matplotlib.Axes which contains the persistence image
721 +
        """
722 +
        ax = ax or plt.gca()
723 +
        ax.matshow(pers_img.T, **{'origin': 'lower'})
724 +
725 +
        # fix and label axes
726 +
        ax.set_xlabel('birth')
727 +
        ax.set_ylabel('persistence')
728 +
        ax.get_xaxis().set_ticks([])
729 +
        ax.get_yaxis().set_ticks([])
730 +
731 +
        # optionally save figure
732 +
        if out_file:
733 +
            plt.savefig(out_file, bbox_inches='tight')
734 +
        
735 +
        return ax
736 +
737 +
738 +
def _transform(pers_dgm, skew=True, resolution=None, weight=None, weight_params=None, kernel=None, kernel_params=None, _bpnts=None, _ppnts=None):
739 +
        """ Transform a persistence diagram into a persistence image.
740 +
        
741 +
        Parameters
742 +
        ----------
743 +
        pers_dgm : (-,2) numpy.ndarray
744 +
            A persistence diagrams.
745 +
        skew : boolean 
746 +
            Flag indicating if diagram(s) need to first be converted to birth-persistence coordinates (default: True).
747 +
        resolution : pair of ints
748 +
            The number of pixels along the birth and persistence axes in the persistence image.
749 +
        weight : callable
750 +
            Function which weights the birth-persistence plane.
751 +
        weight_params : dict
752 +
            Arguments needed to specify the weight function.
753 +
        kernel : callable
754 +
            Cumulative distribution function defining the kernel.
755 +
        kernel_params : dict
756 +
            Arguments needed to specify the kernel function.
757 +
        _bpnts : (N,) numpy.ndarray
758 +
            The birth coordinates of the persistence image pixel locations.
759 +
        _ppnts : (M,) numpy.ndarray
760 +
            The persistence coordinates of the persistence image pixel locations.
761 +
            
762 +
        Returns
763 +
        -------
764 +
        numpy.ndarray
765 +
            (M,N) numpy.ndarray encoding the persistence image corresponding to pers_dgm.
766 +
        """
767 +
        pers_dgm = np.copy(pers_dgm)
768 +
        pers_img = np.zeros(resolution)
769 +
        n = pers_dgm.shape[0]
770 +
        general_flag = True
771 +
772 +
        # if necessary convert from birth-death coordinates to birth-persistence coordinates
773 +
        if skew:
774 +
            pers_dgm[:, 1] = pers_dgm[:, 1] - pers_dgm[:, 0]
775 +
776 +
        # compute weights at each persistence pair
777 +
        wts = weight(pers_dgm[:, 0], pers_dgm[:, 1], **weight_params)
778 +
779 +
        # handle the special case of a standard, isotropic Gaussian kernel
780 +
        if kernel == images_kernels.gaussian:
781 +
            general_flag = False
782 +
            sigma = kernel_params['sigma']
783 +
784 +
            # sigma is specified by a single variance
785 +
            if isinstance(sigma, (int, float)):
786 +
                sigma = np.array([[sigma, 0.0], [0.0, sigma]], dtype=np.float64)
787 +
788 +
            if (sigma[0][0] == sigma[1][1] and sigma[0][1] == 0.0):
789 +
                sigma = np.sqrt(sigma[0][0])
790 +
                for i in range(n):
791 +
                    ncdf_b = images_kernels.norm_cdf((_bpnts - pers_dgm[i, 0]) / sigma)
792 +
                    ncdf_p = images_kernels.norm_cdf((_ppnts - pers_dgm[i, 1]) / sigma)
793 +
                    curr_img = ncdf_p[None, :] * ncdf_b[:, None]
794 +
                    pers_img += wts[i]*(curr_img[1:, 1:] - curr_img[:-1, 1:] - curr_img[1:, :-1] + curr_img[:-1, :-1])
795 +
            else:
796 +
                general_flag = True
797 +
798 +
        # handle the general case
799 +
        if general_flag:
800 +
            bb, pp = np.meshgrid(_bpnts, _ppnts, indexing='ij')
801 +
            bb = bb.flatten(order='C')
802 +
            pp = pp.flatten(order='C')
803 +
            for i in range(n):
804 +
                curr_img = np.reshape(kernel(bb, pp, mu=pers_dgm[i, :], **kernel_params),
805 +
                                      (resolution[0]+1, resolution[1]+1), order='C')
806 +
                pers_img += wts[i]*(curr_img[1:, 1:] - curr_img[:-1, 1:] - curr_img[1:, :-1] + curr_img[:-1, :-1])
807 +
808 +
        return pers_img

@@ -0,0 +1,247 @@
Loading
1 +
"""
2 +
Kernel functions for PersistenceImager() transformer:
3 +
4 +
A valid kernel is a Python function of the form 
5 +
6 +
kernel(x, y, mu=(birth, persistence), **kwargs) 
7 +
8 +
defining a cumulative distribution function(CDF) such that kernel(x, y) = P(X <= x, Y <=y), where x and y are numpy arrays of equal length. 
9 +
10 +
The required parameter mu defines the dependance of the kernel on the location of a persistence pair and is usually taken to be the mean of the probability distribution function associated to kernel CDF.
11 +
"""
12 +
import numpy as np
13 +
from scipy.special import erfc
14 +
15 +
def uniform(x, y, mu=None, width=1, height=1):
16 +
    w1 = np.maximum(x - (mu[0] - width/2), 0)
17 +
    h1 = np.maximum(y - (mu[1] - height/2), 0)
18 +
    
19 +
    w = np.minimum(w1, width)
20 +
    h = np.minimum(h1, height)
21 +
22 +
    return w*h / (width*height)
23 +
24 +
25 +
def gaussian(birth, pers, mu=None, sigma=None):
26 +
    """ Optimized bivariate normal cumulative distribution function for computing persistence images using a Gaussian kernel.
27 +
    
28 +
    Parameters
29 +
    ----------
30 +
    birth : (M,) numpy.ndarray
31 +
        Birth coordinate(s) of pixel corners.
32 +
    pers : (N,) numpy.ndarray
33 +
        Persistence coordinates of pixel corners.
34 +
    mu : (2,) numpy.ndarray
35 +
        Coordinates of the distribution mean (birth-persistence pairs).
36 +
    sigma : float or (2,2) numpy.ndarray 
37 +
        Distribution's covariance matrix or the equal variances if the distribution is standard isotropic.
38 +
    
39 +
    Returns
40 +
    -------
41 +
    float
42 +
        Value of joint CDF at (birth, pers), i.e.,  P(X <= birth, Y <= pers).
43 +
    """
44 +
    if mu is None:
45 +
        mu = np.array([0.0, 0.0], dtype=np.float64)
46 +
    if sigma is None:
47 +
        sigma = np.array([[1.0, 0.0], [0.0, 1.0]], dtype=np.float64)
48 +
49 +
    if sigma[0][1] == 0.0:
50 +
        return sbvn_cdf(birth, pers,
51 +
                         mu_x=mu[0], mu_y=mu[1], sigma_x=sigma[0][0], sigma_y=sigma[1][1])
52 +
    else:
53 +
        return bvn_cdf(birth, pers,
54 +
                        mu_x=mu[0], mu_y=mu[1], sigma_xx=sigma[0][0], sigma_yy=sigma[1][1], sigma_xy=sigma[0][1])
55 +
56 +
57 +
def norm_cdf(x):
58 +
    """ Univariate normal cumulative distribution function (CDF) with mean 0.0 and standard deviation 1.0.
59 +
    
60 +
    Parameters
61 +
    ----------
62 +
    x : float
63 +
        Value at which to evaluate the CDF (upper limit).
64 +
    
65 +
    Returns
66 +
    -------
67 +
    float
68 +
        Value of CDF at x, i.e., P(X <= x), for X ~ N[0,1].
69 +
    """
70 +
    return erfc(-x / np.sqrt(2.0)) / 2.0
71 +
72 +
73 +
def sbvn_cdf(x, y, mu_x=0.0, mu_y=0.0, sigma_x=1.0, sigma_y=1.0):
74 +
    """ Standard bivariate normal cumulative distribution function with specified mean and variances.
75 +
    
76 +
    Parameters
77 +
    ----------
78 +
    x : float or numpy.ndarray of floats
79 +
        x-coordinate(s) at which to evaluate the CDF (upper limit).
80 +
    y : float or numpy.ndarray of floats 
81 +
        y-coordinate(s) at which to evaluate the CDF (upper limit).
82 +
    mu_x : float
83 +
        x-coordinate of the mean.
84 +
    mu_y : float
85 +
        y-coordinate of the mean.
86 +
    sigma_x : float
87 +
        Variance in x.
88 +
    sigma_y : float
89 +
        Variance in y.
90 +
    
91 +
    Returns
92 +
    -------
93 +
    float
94 +
        Value of joint CDF at (x, y), i.e.,  P(X <= birth, Y <= pers).
95 +
    """
96 +
    x = (x - mu_x) / np.sqrt(sigma_x)
97 +
    y = (y - mu_y) / np.sqrt(sigma_y)
98 +
    return norm_cdf(x) * norm_cdf(y)
99 +
100 +
101 +
def bvn_cdf(x, y, mu_x=0.0, mu_y=0.0, sigma_xx=1.0, sigma_yy=1.0, sigma_xy=0.0):
102 +
    """ Bivariate normal cumulative distribution function with specified mean and covariance matrix.
103 +
    
104 +
    Parameters
105 +
    ----------
106 +
    x : float or numpy.ndarray of floats
107 +
        x-coordinate(s) at which to evaluate the CDF (upper limit).
108 +
    y : float or numpy.ndarray of floats 
109 +
        y-coordinate(s) at which to evaluate the CDF (upper limit).
110 +
    mu_x : float
111 +
        x-coordinate of the mean.
112 +
    mu_y : float
113 +
        y-coordinate of the mean.
114 +
    sigma_x : float
115 +
        Variance in x.
116 +
    sigma_y : float
117 +
        Variance in y.
118 +
    sigma_xy : float
119 +
        Covariance of x and y.
120 +
        
121 +
    Returns
122 +
    -------
123 +
    float
124 +
        Value of joint CDF at (x, y), i.e.,  P(X <= birth, Y <= pers).
125 +
        
126 +
    Notes
127 +
    -----
128 +
    Based on the Matlab implementations by Thomas H. Jørgensen (http://www.tjeconomics.com/code/) and Alan Genz (http://www.math.wsu.edu/math/faculty/genz/software/matlab/bvnl.m) using the approach described by Drezner and Wesolowsky (https://doi.org/10.1080/00949659008811236).
129 +
    """
130 +
    dh = -(x - mu_x) / np.sqrt(sigma_xx)
131 +
    dk = -(y - mu_y) / np.sqrt(sigma_yy)
132 +
133 +
    hk = np.multiply(dh, dk)
134 +
    r = sigma_xy / np.sqrt(sigma_xx * sigma_yy)
135 +
136 +
    lg, w, x = gauss_legendre_quad(r)
137 +
138 +
    dim1 = np.ones((len(dh),), dtype=np.float64)
139 +
    dim2 = np.ones((lg,), dtype=np.float64)
140 +
    bvn = np.zeros((len(dh),), dtype=np.float64)
141 +
142 +
    if abs(r) < 0.925:
143 +
        hs = (np.multiply(dh, dh) + np.multiply(dk, dk)) / 2.0
144 +
        asr = np.arcsin(r)
145 +
        sn1 = np.sin(asr * (1.0 - x) / 2.0)
146 +
        sn2 = np.sin(asr * (1.0 + x) / 2.0)
147 +
        dim1w = np.outer(dim1, w)
148 +
        hkdim2 = np.outer(hk, dim2)
149 +
        hsdim2 = np.outer(hs, dim2)
150 +
        dim1sn1 = np.outer(dim1, sn1)
151 +
        dim1sn2 = np.outer(dim1, sn2)
152 +
        sn12 = np.multiply(sn1, sn1)
153 +
        sn22 = np.multiply(sn2, sn2)
154 +
        bvn = asr * np.sum(np.multiply(dim1w, np.exp(np.divide(np.multiply(dim1sn1, hkdim2) - hsdim2,
155 +
                                                               (1 - np.outer(dim1, sn12))))) +
156 +
                           np.multiply(dim1w, np.exp(np.divide(np.multiply(dim1sn2, hkdim2) - hsdim2,
157 +
                                                               (1 - np.outer(dim1, sn22))))), axis=1) / (4 * np.pi) \
158 +
              + np.multiply(norm_cdf(-dh), norm_cdf(-dk))
159 +
    else:
160 +
        if r < 0:
161 +
            dk = -dk
162 +
            hk = -hk
163 +
164 +
        if abs(r) < 1:
165 +
            opmr = (1.0 - r) * (1.0 + r)
166 +
            sopmr = np.sqrt(opmr)
167 +
            xmy2 = np.multiply(dh - dk, dh - dk)
168 +
            xmy = np.sqrt(xmy2)
169 +
            rhk8 = (4.0 - hk) / 8.0
170 +
            rhk16 = (12.0 - hk) / 16.0
171 +
            asr = -1.0 * (np.divide(xmy2, opmr) + hk) / 2.0
172 +
173 +
            ind = asr > 100
174 +
            bvn[ind] = sopmr * np.multiply(np.exp(asr[ind]),
175 +
                                           1.0 - np.multiply(np.multiply(rhk8[ind], xmy2[ind] - opmr),
176 +
                                                             (1.0 - np.multiply(rhk16[ind], xmy2[ind]) / 5.0) / 3.0)
177 +
                                           + np.multiply(rhk8[ind], rhk16[ind]) * opmr * opmr / 5.0)
178 +
179 +
            ind = hk > -100
180 +
            ncdfxmyt = np.sqrt(2.0 * np.pi) * norm_cdf(-xmy / sopmr)
181 +
            bvn[ind] = bvn[ind] - np.multiply(np.multiply(np.multiply(np.exp(-hk[ind] / 2.0), ncdfxmyt[ind]), xmy[ind]),
182 +
                                              1.0 - np.multiply(np.multiply(rhk8[ind], xmy2[ind]),
183 +
                                                                (1.0 - np.multiply(rhk16[ind], xmy2[ind]) / 5.0) / 3.0))
184 +
            sopmr = sopmr / 2
185 +
            for ix in [-1, 1]:
186 +
                xs = np.multiply(sopmr + sopmr * ix * x, sopmr + sopmr * ix * x)
187 +
                rs = np.sqrt(1 - xs)
188 +
                xmy2dim2 = np.outer(xmy2, dim2)
189 +
                dim1xs = np.outer(dim1, xs)
190 +
                dim1rs = np.outer(dim1, rs)
191 +
                dim1w = np.outer(dim1, w)
192 +
                rhk16dim2 = np.outer(rhk16, dim2)
193 +
                hkdim2 = np.outer(hk, dim2)
194 +
                asr1 = -1.0 * (np.divide(xmy2dim2, dim1xs) + hkdim2) / 2.0
195 +
196 +
                ind1 = asr1 > -100
197 +
                cdim2 = np.outer(rhk8, dim2)
198 +
                sp1 = 1.0 + np.multiply(np.multiply(cdim2, dim1xs), 1.0 + np.multiply(rhk16dim2, dim1xs))
199 +
                ep1 = np.divide(np.exp(np.divide(-np.multiply(hkdim2, (1.0 - dim1rs)),
200 +
                                                 2.0 * (1.0 + dim1rs))), dim1rs)
201 +
                bvn = bvn + np.sum(np.multiply(np.multiply(np.multiply(sopmr, dim1w), np.exp(np.multiply(asr1, ind1))),
202 +
                                               np.multiply(ep1, ind1) - np.multiply(sp1, ind1)), axis=1)
203 +
            bvn = -bvn / (2.0 * np.pi)
204 +
205 +
        if r > 0:
206 +
            bvn = bvn + norm_cdf(-np.maximum(dh, dk))
207 +
        elif r < 0:
208 +
            bvn = -bvn + np.maximum(0, norm_cdf(-dh) - norm_cdf(-dk))
209 +
210 +
    return bvn
211 +
212 +
213 +
def gauss_legendre_quad(r):
214 +
    """ Return weights and abscissae for the Legendre-Gauss quadrature integral approximation.
215 +
    
216 +
    Parameters
217 +
    ----------
218 +
    r : float
219 +
        Correlation
220 +
    
221 +
    Returns
222 +
    -------
223 +
    tuple
224 +
        Number of points in the Gaussian quadrature rule, quadrature weights, and quadrature points.
225 +
    """
226 +
    if np.abs(r) < 0.3:
227 +
        lg = 3
228 +
        w = np.array([0.1713244923791705, 0.3607615730481384, 0.4679139345726904])
229 +
        x = np.array([0.9324695142031522, 0.6612093864662647, 0.2386191860831970])
230 +
    elif np.abs(r) < 0.75:
231 +
        lg = 6
232 +
        w = np.array([.04717533638651177, 0.1069393259953183, 0.1600783285433464,
233 +
                      0.2031674267230659, 0.2334925365383547, 0.2491470458134029])
234 +
        x = np.array([0.9815606342467191, 0.9041172563704750, 0.7699026741943050,
235 +
                      0.5873179542866171, 0.3678314989981802, 0.1252334085114692])
236 +
    else:
237 +
        lg = 10
238 +
        w = np.array([0.01761400713915212, 0.04060142980038694, 0.06267204833410906,
239 +
                      0.08327674157670475, 0.1019301198172404, 0.1181945319615184,
240 +
                      0.1316886384491766, 0.1420961093183821, 0.1491729864726037,
241 +
                      0.1527533871307259])
242 +
        x = np.array([0.9931285991850949, 0.9639719272779138, 0.9122344282513259,
243 +
                      0.8391169718222188, 0.7463319064601508, 0.6360536807265150,
244 +
                      0.5108670019508271, 0.3737060887154196, 0.2277858511416451,
245 +
                      0.07652652113349733])
246 +
247 +
    return lg, w, x

@@ -122,12 +122,12 @@
Loading
122 122
123 123
    Parameters
124 124
    -----------
125 -
    AG: np.array (|V(G)|×|V(G)|)
126 -
        (Sparse) adjacency matrix of graph G, or an iterable of
125 +
    AG: (N,N) np.array
126 +
        (Sparse) adjacency matrix of graph G with N vertices, or an iterable of
127 127
        adjacency matrices if AH=None.
128 -
    AH: np.array (|V(H)|×|V(H)|)
129 -
        (Sparse) adjacency matrix of graph H, or None.
130 -
    mapping_sample_size_order: np.array (2)
128 +
    AH: (M,M) np.array
129 +
        (Sparse) adjacency matrix of graph H with M vertices, or None.
130 +
    mapping_sample_size_order: (2,) np.array 
131 131
        Parameter that regulates the number of mappings to sample when
132 132
        tightening upper bound of the mGH distance.
133 133
@@ -184,12 +184,12 @@
Loading
184 184
185 185
    Parameters
186 186
    -----------
187 -
    AG: np.array (|V(G)|×|V(G)|)
188 -
        (Sparse) adjacency matrix of simple unweighted graph G.
187 +
    AG: (N,N) np.array
188 +
        (Sparse) adjacency matrix of simple unweighted graph G with N vertices.
189 189
190 190
    Returns
191 191
    --------
192 -
    DG: np.array (|V(G)|×|V(G)|)
192 +
    DG: (N,N)  np.array
193 193
        (Dense) distance matrix of the compact metric space
194 194
        representation of G based on its shortest path lengths.
195 195
    """

@@ -0,0 +1,70 @@
Loading
1 +
"""
2 +
Weight Functions for PersistenceImager() transformer:
3 +
4 +
A valid weight function is a Python function of the form 
5 +
6 +
weight(birth, persistence, **kwargs) 
7 +
8 +
defining a scalar-valued function over the birth-persistence plane, where birth and persistence are assumed to be numpy arrays of equal length. To ensure stability, functions should vanish continuously at the line persistence = 0.
9 +
"""
10 +
import numpy as np
11 +
12 +
def linear_ramp(birth, pers, low=0.0, high=1.0, start=0.0, end=1.0):
13 +
    """ Continuous peicewise linear ramp function which is constant below and above specified input values.
14 +
    
15 +
    Parameters
16 +
    ----------
17 +
    birth : (N,) numpy.ndarray
18 +
        Birth coordinates of a collection of persistence pairs.
19 +
    pers : (N,) numpy.ndarray
20 +
        Persistence coordinates of a collection of persistence pairs.
21 +
    low : float
22 +
        Minimum weight. 
23 +
    high : float
24 +
       Maximum weight.
25 +
    start : float
26 +
        Start persistence value of linear transition from low to high weight.
27 +
    end : float
28 +
        End persistence value of linear transition from low to high weight.
29 +
    
30 +
    Returns
31 +
    -------
32 +
    (N,) numpy.ndarray
33 +
        Weights at the persistence pairs.
34 +
    """
35 +
    try:
36 +
        n = len(birth)
37 +
    except:
38 +
        n = 1
39 +
        birth = [birth]
40 +
        pers = [pers]
41 +
42 +
    w = np.zeros((n,))
43 +
    for i in range(n):
44 +
        if pers[i] < start:
45 +
            w[i] = low
46 +
        elif pers[i] > end:
47 +
            w[i] = high
48 +
        else:
49 +
            w[i] = (pers[i] - start) * (high - low) / (end - start) + low
50 +
51 +
    return w
52 +
53 +
def persistence(birth, pers, n=1.0):
54 +
    """ Continuous monotonic function which weight a persistence pair (b,p) by p^n for some n > 0.
55 +
    
56 +
    Parameters
57 +
    ----------
58 +
    birth : (N,) numpy.ndarray
59 +
        Birth coordinates of a collection of persistence pairs.
60 +
    pers : (N,) numpy.ndarray
61 +
        Persistence coordinates of a collection of persistence pairs.
62 +
    n : positive float
63 +
        Exponent of persistence weighting function.
64 +
    
65 +
    Returns
66 +
    -------
67 +
    (N,) numpy.ndarray
68 +
        Weights at the persistence pairs.
69 +
    """
70 +
    return pers ** n
Files Coverage
persim 75.30%
Project Totals (12 files) 75.30%
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading