NeurIPS 2023 -  TopP&R: Robust Support Estimation Approach for Evaluating Fidelity and Diversity in Generative Models Official Code
import setuptools

with open("README.md", "r") as fh:
    long_description = fh.read()

setuptools.setup(
    name="top_pr",
    version="0.2.1",
    author="Pumjun Kim",
    description="TopP&R: Robust Support Estimation Approach "
                "for Evaluating Fidelity and Diversity in Generative Models.",
    long_description=long_description,
    long_description_content_type="text/markdown",
    url="https://github.com/LAIT-CVLab/TopPR",
    packages=setuptools.find_packages(),
    keywords=['toppr', 'evaluation metric', 'topological metric', 'precision and recall'],
    classifiers=[
        "Programming Language :: Python :: 3",
        "License :: OSI Approved :: MIT License",
        "Operating System :: OS Independent",
    ],
    install_requires=[
        'torch',
        'torchvision',
        'numpy',
        'scikit-learn',
        'scipy',
        'tqdm',
        'matplotlib',
        'gudhi'
    ],
)
setup.py
LAIT-CVLab/TopPR
from .top_pr import compute_top_pr
top_pr/__init__.py
LAIT-CVLab/TopPR
import numpy as np
import torch
import pickle
from torchvision import datasets, transforms

############################################################
############ Gaussian Mode Drop Data Generator #############
############################################################
def gaussian_mode_drop(method, ratio = 1, isnumpy = False):
    # generate mixture of gaussian distribution with multi modes (7 modes)
    num_samples = 1000
    feature_dim = 64
    features = np.random.normal(loc=0.0, scale=1.0, size=[num_samples, feature_dim])
    features = np.vstack((features, np.random.normal(loc=30.0, scale=1.0, size=[num_samples, feature_dim])))
    features = np.vstack((features, np.random.normal(loc=20.0, scale=1.0, size=[num_samples, feature_dim])))
    features = np.vstack((features, np.random.normal(loc=10.0, scale=1.0, size=[num_samples, feature_dim])))
    features = np.vstack((features, np.random.normal(loc=-10.0, scale=1.0, size=[num_samples, feature_dim])))
    features = np.vstack((features, np.random.normal(loc=-20.0, scale=1.0, size=[num_samples, feature_dim])))
    features = np.vstack((features, np.random.normal(loc=-30.0, scale=1.0, size=[num_samples, feature_dim])))

    # assign class
    Class = np.zeros(len(features))
    Class[0:1000] = 0
    Class[1000:1000*2] = 1
    Class[1000*2:1000*3] = 2
    Class[1000*3:1000*4] = 3
    Class[1000*4:1000*5] = 4
    Class[1000*5:1000*6] = 5
    Class[1000*6:1000*7] = 6

    circle0 = np.array([])
    circle1 = np.array([])
    circle2 = np.array([])
    circle3 = np.array([])
    circle4 = np.array([])
    circle5 = np.array([]) 
    circle6 = np.array([])

    for iloop in range(len(Class)):
        if (Class[iloop] == 0):
            circle0 = np.append(circle0, iloop)
        elif (Class[iloop] == 1):
            circle1 = np.append(circle1, iloop)
        elif (Class[iloop] == 2):
            circle2 = np.append(circle2, iloop)
        elif (Class[iloop] == 3):
            circle3 = np.append(circle3, iloop)
        elif (Class[iloop] == 4):
            circle4 = np.append(circle4, iloop)
        elif (Class[iloop] == 5):
            circle5 = np.append(circle5, iloop)
        elif (Class[iloop] == 6):
            circle6 = np.append(circle6, iloop)
        
    circle0 = list(map(int, circle0))
    circle1 = list(map(int, circle1))
    circle2 = list(map(int, circle2))
    circle3 = list(map(int, circle3))
    circle4 = list(map(int, circle4))
    circle5 = list(map(int, circle5))
    circle6 = list(map(int, circle6))

    # sequential Drop
    if (method == 'sequential'):
        if (ratio == 0):
            return features
        elif (ratio == 1):
            temp_data = features[circle0,:]
            temp_data = np.vstack((temp_data, features[circle1,:]))
            temp_data = np.vstack((temp_data, features[circle2,:]))
            temp_data = np.vstack((temp_data, features[circle3,:]))
            temp_data = np.vstack((temp_data, features[circle4,:]))
            temp_data = np.vstack((temp_data, features[circle5,:]))
            temp = temp_data[np.random.choice(range(len(temp_data)), len(circle6), replace = True),:]
            temp = temp + np.random.normal(loc = 0, scale = 0.1, size = [len(circle6), feature_dim])
            temp_data = np.vstack((temp_data, temp))
            return temp_data

        elif (ratio == 2):
            temp_data = features[circle0,:]
            temp_data = np.vstack((temp_data, features[circle1,:]))
            temp_data = np.vstack((temp_data, features[circle2,:]))
            temp_data = np.vstack((temp_data, features[circle3,:]))
            temp_data = np.vstack((temp_data, features[circle4,:]))
            temp = temp_data[np.random.choice(range(len(temp_data)), len(circle5) + len(circle6), replace = True),:]
            temp = temp + np.random.normal(loc = 0, scale = 0.1, size = [len(circle5) + len(circle6), feature_dim])
            temp_data = np.vstack((temp_data, temp))
            return temp_data

        elif (ratio == 3):
            temp_data = features[circle0,:]
            temp_data = np.vstack((temp_data, features[circle1,:]))
            temp_data = np.vstack((temp_data, features[circle2,:]))
            temp_data = np.vstack((temp_data, features[circle3,:]))
            temp = temp_data[np.random.choice(range(len(temp_data)), len(circle4) + len(circle5) + len(circle6), replace = True),:]
            temp = temp + np.random.normal(loc = 0, scale = 0.1, size = [len(circle4) + len(circle5) + len(circle6), feature_dim])
            temp_data = np.vstack((temp_data, temp))
            return temp_data

        elif (ratio == 4):
            temp_data = features[circle0,:]
            temp_data = np.vstack((temp_data, features[circle1,:]))
            temp_data = np.vstack((temp_data, features[circle2,:]))
            temp = temp_data[np.random.choice(range(len(temp_data)), len(circle3) + len(circle4) + len(circle5) + len(circle6), replace = True),:]
            temp = temp + np.random.normal(loc = 0, scale = 0.1, size = [len(circle3) + len(circle4) + len(circle5) + len(circle6), feature_dim])
            temp_data = np.vstack((temp_data, temp))
            return temp_data
        
        elif (ratio == 5):
            temp_data = features[circle0,:]
            temp_data = np.vstack((temp_data, features[circle1,:]))
            temp = temp_data[np.random.choice(range(len(temp_data)), len(circle2) + len(circle3) + len(circle4) + len(circle5) + len(circle6), replace = True),:]
            temp = temp + np.random.normal(loc = 0, scale = 0.1, size = [len(circle2) + len(circle3) + len(circle4) + len(circle5) + len(circle6), feature_dim])
            temp_data = np.vstack((temp_data, temp))
            return temp_data

        elif (ratio == 6):
            temp_data = features[circle0,:]
            temp = temp_data[np.random.choice(range(len(temp_data)), len(circle1) + len(circle2) + len(circle3) + len(circle4) + len(circle5) + len(circle6), replace = True),:]
            temp = temp + np.random.normal(loc = 0, scale = 0.1, size = [len(circle1) + len(circle2) + len(circle3) + len(circle4) + len(circle5) + len(circle6), feature_dim])
            temp_data = np.vstack((temp_data, temp))
            return temp_data

    # simultaneous drop        
    if (method == 'simultaneous'):
        temp_data = features[circle0,:]
        sample_data = features[circle1,:]
        sample_data = np.vstack((sample_data, features[circle2,:]))
        sample_data = np.vstack((sample_data, features[circle3,:]))
        sample_data = np.vstack((sample_data, features[circle4,:]))
        sample_data = np.vstack((sample_data, features[circle5,:]))
        sample_data = np.vstack((sample_data, features[circle6,:]))

        temp = sample_data[np.random.choice(range(len(sample_data)), int(round(len(sample_data)*(1-ratio))), replace = False),:]
        temp2 = temp_data[np.random.choice(range(len(features[circle0,:])), round(len(sample_data)*ratio), replace = True),:]
        temp2 = temp2 + np.random.normal(loc = 0, scale = 0.1, size = [round(len(sample_data)*ratio), feature_dim])
        temp_data = np.vstack((temp_data, temp))
        temp_data = np.vstack((temp_data, temp2))
        return temp_data
top_pr/mode_drop.py
LAIT-CVLab/TopPR
import matplotlib.pyplot as plot
import numpy as np
import random
import os
from sklearn.metrics.pairwise import euclidean_distances
from scipy.spatial import distance
import torch
import gudhi
from tqdm import tqdm

def Fix_seed(seed: int = 42):
    random.seed(seed)
    np.random.seed(seed)
    os.environ["PYTHONHASHSEED"] = str(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)
    torch.backends.cudnn.deterministic = True
    torch.backends.cudnn.benchmark = False 
SEED = 42

###########################################
#   Automatic Grid Search
###########################################
def set_grid(data):
    import numpy as np

    # find min max
    dim = len(data[0])
    mins = np.array([])
    maxs = np.array([])
    for dims in range(dim):
        mins = np.append(mins, min(data[:,dims]))
        maxs = np.append(maxs, max(data[:,dims]))
    
    # set grid
    # 2 dimensional data
    if (len(mins) == 2):
        xval = np.linspace(mins[0], maxs[0], 1000)
        yval = np.linspace(mins[1], maxs[1], 1000)
        positions = np.array([[u,v] for u in xval for v in yval])
    # 3 dimensional data
    elif (len(mins) == 3):
        xval = np.linspace(mins[0], maxs[0], 100)
        yval = np.linspace(mins[1], maxs[1], 100)
        zval = np.linspace(mins[2], maxs[2], 100)
        positions = np.array([[u,v,k] for u in xval for v in yval for k in zval])
    
    return positions

###########################################
#   KDE with Epanechinikov Kernel
###########################################
def compact_KDE(data, position, h, kernel = "cosine"):
    # compact kernel options = {"epanechinikov", "cosine"}
    p_hat = np.array([])
    dist = euclidean_distances(position, data)

    # Epanechinikov kernel
    if (kernel == "epanechinikov"):
        for iloop in range(len(dist)):
            sample_score = dist[iloop][np.where(dist[iloop] ** 2 <= (h**2))]
            p_hat = np.append(p_hat, 
                        (1 / len(data)) * ((3 / (4 * h)) ** len(data[0])) * ((len(sample_score)) - 
                        np.sum(sample_score / (h ** 2))))
        return p_hat
    
    # Cosine kernel
    elif (kernel == "cosine"):
        for iloop in range(len(dist)):
            sample_score = dist[iloop][np.where(dist[iloop] ** 2 <= (h**2))]
            p_hat = np.append(p_hat, 
                        (1 / len(data)) * ((np.pi / (4 * h)) ** len(data[0])) * 
                        np.sum(np.cos((np.pi / 2) * (sample_score / h))))
        return p_hat

###########################################
# Confidence Band
###########################################
def confband_est(data, h, alpha = 0.1, kernel = "cosine", p_val = True, repeat = 10):
    if not isinstance(data, np.ndarray):
        data = np.asarray(data)
    
    if (len(data[0]) <= 3):
        positions = set_grid(data)
    else:
        positions = data
    
    # p_hat
    p_hat = compact_KDE(data, positions, h, kernel = kernel)

    # p_tilde
    theta_star = np.array([])
    for iloop in range(repeat):
        data_bs = data[np.random.choice(np.arange(len(data)), size = len(data), replace = True)]
        p_tilde = compact_KDE(data_bs, positions, h, kernel = kernel)
        
        # theta
        theta_star = np.append(theta_star, np.sqrt(len(data))*np.max(np.abs(p_hat-p_tilde)))
    
    # q_alpha
    if (len(theta_star) == 0):
        q_alpha = 0
    else:
        q_alpha = np.quantile(theta_star, 1 - alpha)
    
    # confidence band
    if (p_val == True):
        return q_alpha/np.sqrt(len(data)), p_hat
    else:
        return q_alpha/np.sqrt(len(data))

###########################################
# BandWidth estimator
###########################################
def bandwidth_est(data, bandwidth_list = [], confidence_band = False, kernel = 'cosine', alpha = 0.1, Plot = False):
    # non-compact kernel list = {'gaussian','exponential'} | compact kernel list = {'tophat','epanechnikov','linear','cosine'}
    if not isinstance(data, np.ndarray):
        data = np.asarray(data)

    # estimate bandwidth candidates with "balloon estimator" (variable-bandwidth estimation)
    if (len(bandwidth_list) == 0):
        dist = distance.cdist(data, data, metric = "euclidean")
        dist = dist[:-50,]
        for iloop in range(len(dist)):
            if (iloop == 0):
                balloon_est = np.array(sorted(dist[iloop,(iloop+1):])[50-1])
            else:
                balloon_est = np.append(balloon_est, sorted(dist[iloop,(iloop+1):])[50-1])
        balloon_est = sorted(balloon_est)
        bandwidth_list = balloon_est[int(len(balloon_est)*0.05)-1] # top 5% estimated bandwidth
        bandwidth_list = np.append(bandwidth_list, balloon_est[int(len(balloon_est)*0.2)-1]) # top 20% estimated bandwidth
        bandwidth_list = np.append(bandwidth_list, balloon_est[int(len(balloon_est)*0.35)-1]) # top 35% estimated bandwidth
        bandwidth_list = np.append(bandwidth_list, balloon_est[int(len(balloon_est)*0.5)-1]) # median estimated bandwidth
        bandwidth_list = np.append(bandwidth_list, balloon_est[int(len(balloon_est)*0.65)-1]) # top 65% estimated bandwidth
        bandwidth_list = np.append(bandwidth_list, balloon_est[int(len(balloon_est)*0.8)-1]) # top 80% estimated bandwidth
        bandwidth_list = np.append(bandwidth_list, balloon_est[int(len(balloon_est)*0.95)-1]) # top 95% estimated bandwidth

    # estimate bandwidth
    n_h0 = np.array([])
    s_h0 = np.array([])
    cn_list = np.array([])
    for h in tqdm(bandwidth_list):
        # confidence band & p_hat
        cn = confband_est(data, h, alpha = alpha, kernel = kernel, p_val = False)
        cn_list = np.append(cn_list, cn)

        grid = set_grid(data)
        # filtration
        p_hat = compact_KDE(data, grid, h, kernel = kernel)
        PD = gudhi.CubicalComplex(dimensions = [round(len(grid)**(1/grid.shape[1])),round(len(grid)**(1/grid.shape[1]))],
                                top_dimensional_cells = -p_hat).persistence()
        
        # measure life length of all homology
        l_h0 = np.array([])
        for iloop in range(len(PD)):
            if (PD[iloop][0] == 0):
                if (np.abs(PD[iloop][1][1]-PD[iloop][1][0]) != float("inf")):
                    l_h0 = np.append(l_h0, np.abs(PD[iloop][1][1]-PD[iloop][1][0]))
        
        # N(h)
        n_h0 = np.append(n_h0, sum(l_h0 > cn)+1)
        
        # S(h)
        S_h0 = l_h0 - cn
        s_h0 = np.append(s_h0, sum(list(filter(lambda S_h0 : S_h0 > 0, S_h0))))
        print('bandwidth: ',h,', N_0(h): ',n_h0[-1],', S_0(h): ',s_h0[-1],', cn: ',cn)
    
    if (Plot == True): 
        fig = plot.figure(figsize = (10,2))
        for i in range(1,3):
            axes = fig.add_subplot(1,2,i)
            if (i == 1):
                axes.set_title(r"N(h) for $H_0$",fontsize = 15)
                plot.plot(bandwidth_list, n_h0, color = [133/255, 185/255, 190/255], linestyle = '-', linewidth = 4)
                plot.scatter(bandwidth_list, n_h0, color = [133/255, 185/255, 190/255], s=100)
                
            elif (i == 2):
                axes.set_title(r"S(h) for $H_0$",fontsize = 15)
                plot.plot(bandwidth_list, s_h0, color = [255/255, 110/255, 97/255], linestyle = '-', linewidth = 4)
                plot.scatter(bandwidth_list, s_h0, color = [255/255, 110/255, 97/255], s=100)

    try:
        if (sum(s_h0 == max(s_h0)) == 1):
            if (confidence_band == True):
                return bandwidth_list[s_h0.tolist().index(max(s_h0))], cn_list[s_h0.tolist().index(max(s_h0))]
            elif (confidence_band == False):
                return bandwidth_list[s_h0.tolist().index(max(s_h0))]
        else:
            return bandwidth_list[0]
    except Exception as e:
        print(e)
        raise SystemExit

###########################################
#   Top P&R
###########################################
def compute_top_pr(*, real_features, fake_features, alpha = 0.1, kernel = "cosine", random_proj = True, f1_score = True):
    """
    Computing Top Precision and Recall
        Args: 
            real_features (n, d): input real features
            fake_features (n, d): input fake features
            alpha (float): significance level alpha in confidence band estimation (default=0.1)
            kernel (str): kernel for KDE                                          (default='cosine')
            random_proj (bool): If true, it will add linear layer from Pytorch library for random projection. (default=True)
                                However, If the dimension of the feature is lower than 32, even though random_proj is True, random projection will not be turned on.
            f1_score (bool): If True, it caculates f1 score for getting a 1-score evaluation (default=True)
        Returns:
            evaluation score (dict): fidelity, diversity and (opitionally f1 score.)

    """
    # match data format for random projection
    if (torch.is_tensor(real_features) == False):
        real_features = torch.tensor(real_features, dtype = torch.float32)
    if (torch.is_tensor(fake_features) == False):
        fake_features = torch.tensor(fake_features, dtype = torch.float32)
        
    # random projection
    if ((random_proj == True) and (real_features.size()[1] > 32)):
        Fix_seed(SEED)
        projection = torch.nn.Linear(real_features.size()[1], 32, bias = False).eval()
        torch.manual_seed(99)
        torch.nn.init.xavier_normal_(projection.weight)
        for param in projection.parameters():
            param.requires_grad_(False)
        real_features = projection(real_features)
        fake_features = projection(fake_features)

    # to numpy
    real_features = real_features.detach().cpu().numpy()
    fake_features = fake_features.detach().cpu().numpy()

    # use bandwidth estimator to calculate Top P&R
    if (len(real_features[0]) <= 3):
        bandwidth_r, c_r = bandwidth_est(real_features, bandwidth_list = [], confidence_band = True, alpha = alpha)
        bandwidth_f, c_g = bandwidth_est(fake_features, bandwidth_list = [], confidence_band = True, alpha = alpha)

    # use balloon estimator to calculate Top P&R
    else:
        k = len(real_features[0]) * 5
        dist = distance.cdist(real_features, real_features, metric = "euclidean")
        dist = dist[:-k,]
        for iloop in range(len(dist)):
            if (iloop == 0):
                balloon_est = np.array(sorted(dist[iloop,(iloop+1):])[k-1])
            else:
                balloon_est = np.append(balloon_est, sorted(dist[iloop,(iloop+1):])[k-1])
        balloon_est = sorted(balloon_est)
        bandwidth_r = balloon_est[int(len(balloon_est)*0.5)-1] # median balloon estimated bandwidth

        # estimate bandwidth for fake with "balloon estimator" (variable-bandwidth estimation)
        dist = distance.cdist(fake_features, fake_features, metric = "euclidean")
        dist = dist[:-k,]
        for iloop in range(len(dist)):
            if (iloop == 0):
                balloon_est = np.array(sorted(dist[iloop,(iloop+1):])[k-1])
            else:
                balloon_est = np.append(balloon_est, sorted(dist[iloop,(iloop+1):])[k-1])
        balloon_est = sorted(balloon_est)
        bandwidth_f = balloon_est[int(len(balloon_est)*0.5)-1] # median balloon estimated bandwidth

        # estimation of confidence band and manifold
        c_r, score_rr = confband_est(data = real_features, h = bandwidth_r, alpha = alpha, kernel = kernel, p_val = True)
        c_g, score_gg = confband_est(data = fake_features, h = bandwidth_f, alpha = alpha, kernel = kernel, p_val = True)
    
    # count significant real & fake samples
    num_real = np.sum(score_rr > c_r)
    num_fake = np.sum(score_gg > c_g)

    # count significant fake samples on real manifold
    score_rg = compact_KDE(fake_features, real_features, bandwidth_f, kernel = kernel)
    inter_r = np.sum((score_rr > c_r) * (score_rg > c_g))

    # count significant real samples on fake manifold
    score_gr = compact_KDE(real_features, fake_features, bandwidth_r, kernel = kernel)
    inter_g = np.sum((score_gg > c_g) * (score_gr > c_r))

    # topological precision and recall
    t_precision = inter_g / num_fake
    t_recall = inter_r / num_real

    # top f1-score
    if (f1_score == True):
        if (t_precision > 0.0001 and t_recall > 0.0001):
            F1_score = 2/((1/t_precision) + (1/t_recall))
        else:
            F1_score = 0
        return dict(fidelity = t_precision, diversity = t_recall, Top_F1 = F1_score)
    else:
        return dict(fidelity = t_precision, diversity = t_recall)
top_pr/top_pr.py
LAIT-CVLab/TopPR