PrivacyScanData/truncate/new2/model-change-paper-main/util/rkhs.py

# Python implementation of the Laplace Kernel RKHS data-based norm of Karzand and Nowak
# @author : Kevin Miller

'''
Currently only implemented for binary (as is the case in Karzand and Nowak paper)

So current implementation doesn't scale well, because we are calculating the Gram matrix on the go.
    * Try doing Nystrom on Gram matrix, use this to compare?
    * Could use the full matrix, and show how long it takes then? This will probably do better than our method?

THOUGHT- Should we just phrase my method as a generalization and speed improvement of Nowak and Karzand?
    * Then all the tests their better performance can be put into perspective.
'''

import numpy as np
from scipy.spatial.distance import cdist
from argparse import ArgumentParser

import matplotlib.pyplot as plt
import pickle



class RKHSClassifierOld(object):
    def __init__(self, X, sigma):
        self.X = X
        self.N = self.X.shape[0]
        self.sigma = sigma
        self.f = None

    def calculate_model(self, labeled, y):
        self.labeled = labeled
        self.unlabeled = list(filter(lambda x: x not in self.labeled, range(self.N)))
        self.y = y

        self.K_inv = np.linalg.inv(np.exp(-cdist(self.X[self.labeled, :], self.X[self.labeled,:]) / self.sigma))
        self.f = np.empty(self.N)
        self.f[self.labeled] = self.y
        self.K_ul = np.exp(-cdist(self.X[self.unlabeled, :], self.X[self.labeled,:]) / self.sigma)
        self.f[self.unlabeled] = self.K_ul @ self.K_inv @ np.array(self.y) # Kul @ Kll^{-1} y

    def update_model_old(self, Q, yQ):
        len_lab_old = len(self.labeled)
        if self.f is None:
            print("No previous model calculated, so we will do initial calculation")
            self.calculate_model(Q, yQ)
            return
        self.labeled += Q
        self.y += yQ

        aQ = np.exp(-cdist(self.X[self.labeled,:], self.X[Q, :])/ self.sigma)
        aQ1 = aQ[:len_lab_old,:]
        aQ2 = aQ[-len(Q):,:]
        Z = np.linalg.inv(aQ2 - aQ1.T @ self.K_inv @ aQ1)
        A12 = self.K_inv @ aQ1 @ Z
        A11 = A12 @ aQ1.T
        A11.ravel()[::len_lab_old+1] += 1.
        A11 = A11 @ self.K_inv
        self.K_inv = np.hstack((A11, -A12))
        self.K_inv = np.vstack((self.K_inv, np.hstack((-A12.T, Z))))
        unl_Qi, unl_Q = zip(*list(filter(lambda x: x[1] not in Q, enumerate(self.unlabeled))))
        K_Qu_Q = np.exp(-cdist(self.X[unl_Q,:], self.X[Q,:])/self.sigma)
        self.K_ul = np.hstack((self.K_ul[unl_Qi, :], K_Qu_Q))
        self.unlabeled = list(filter(lambda x: x not in self.labeled, range(self.N)))
        self.f[Q] = np.array(yQ)
        self.f[self.unlabeled] = self.K_ul @ self.K_inv @ np.array(self.y) # Kul @ Kll^{-1} y
        return

    def update_model(self, Q, yQ):
        unl_Q = [self.unlabeled.index(k) for k in Q]
        unl_notQ = list(filter(lambda x: x not in unl_Q, range(len(self.unlabeled))))
        notQ = list(filter(lambda x: x not in Q, self.unlabeled))

        # update f
        SQ_inv = np.linalg.inv(np.exp(-cdist(self.X[Q,:], self.X[Q,:])/self.sigma))
        K_notQ_Q = np.exp(-cdist(self.X[notQ,:], self.X[Q,:])/self.sigma)
        sc = K_notQ_Q - self.K_ul[unl_notQ,:] @ self.K_inv @ self.K_ul[unl_Q,:].T
        self.f[notQ] += sc @ SQ_inv @ (np.array(yQ) - self.f[Q])
        self.f[Q] = yQ

        # update self.K_inv, self.K_ul for future calculations
        # calculate self.K_inv new by block matrix inversion formula
        Mat = self.K_ul[unl_Q,:].T @ SQ_inv @ self.K_ul[unl_Q,:] @ self.K_inv
        Mat.ravel()[::len(self.labeled)+1] += 1.0
        new_Kinv_lower_left = -SQ_inv @ self.K_ul[unl_Q,:] @ self.K_inv
        bottom_new_Kinv = np.hstack((new_Kinv_lower_left, SQ_inv))
        self.K_inv = self.K_inv @ Mat
        self.K_inv = np.hstack((self.K_inv, new_Kinv_lower_left.T))
        self.K_inv = np.vstack((self.K_inv, bottom_new_Kinv))

        # self.K_ul
        self.K_ul = self.K_ul[unl_notQ,:]
        self.K_ul = np.hstack((self.K_ul, K_notQ_Q))

        self.labeled.extend(Q)
        self.unlabeled = notQ

        return
    # def look_ahead_db_norm(self, k):
    #     ak = self.K_ul[self.unlabeled.index(k), :]
    #     Z = 1./(1. - np.inner(ak, self.K_inv @ ak))
    #     A12 = self.K_inv @ ak * Z
    #     akK_inv = self.K_inv @ ak
    #     A11 = self.K_inv + Z * (np.outer(akK_inv, akK_inv))
    #     K_inv_new = np.hstack((A11, -A12[:,np.newaxis]))
    #     K_inv_new = np.vstack((K_inv_new, np.hstack((-A12[np.newaxis,:], np.array([[Z]])))))
    #     f_k = self.f.copy()
    #     f_k[k] = np.sign(self.f[k])# f_k [k] = yk = lowest interpolating label
    #     unl_i, unl_k = zip(*list(filter(lambda x : x[1] != k, enumerate(self.unlabeled))))
    #     K_ku_k = np.exp(-cdist(self.X[unl_k,:], self.X[k,:][np.newaxis,:])/self.sigma)
    #     f_k[list(unl_k)] = np.hstack((self.K_ul[unl_i, :], K_ku_k)) @ K_inv_new \
    #                             @ np.array(self.y + [f_k[k]])
    #     return np.linalg.norm(f_k - self.f)#**2./len(self.unlabeled)
    #
    # def look_ahead_db_norm2(self, k):
    #     modelk = RKHSClassifier(self.X, self.sigma)
    #     modelk.calculate_model(self.labeled[:] + [k], self.y[:] + [np.sign(self.f[k])])
    #
    #     return np.linalg.norm(modelk.f - self.f)

    def look_ahead_db_norms(self, Cand):
        unl_Cand = [self.unlabeled.index(k) for k in Cand]
        KCandl = self.K_ul[unl_Cand,:]
        sc = np.exp(-cdist(self.X[self.unlabeled, :], self.X[Cand,:])/self.sigma) \
                - self.K_ul @ self.K_inv @ KCandl.T
        return (1. - np.absolute(self.f[Cand]))*np.linalg.norm(sc, axis=0)/sc[unl_Cand, range(len(Cand))]

        vals = []
        for ii, i in enumerate(unl_Cand):
            vals.append((1. - np.absolute(self.f[self.unlabeled[i]]))/(1. - B[i,ii]) * np.linalg.norm(B[:,ii]))
        return np.array(vals)

    # def look_ahead_db_norms2(self, Cand):
    #     unl_Cand = list(filter(lambda x: self.unlabeled[x] in Cand, range(len(self.unlabeled))))
    #     sc_submatrix = np.exp(-cdist(self.X[self.unlabeled, :], self.X[Cand,:])/self.sigma) \
    #                         - self.K_ul @ self.K_inv @ self.K_ul[unl_Cand,:].T
    #     return np.abs(1. - self.f[Cand]) *  np.linalg.norm(sc_submatrix, axis=0) / np.abs(sc_submatrix[unl_Cand, range(len(Cand))])

    # def look_ahead_db_norms2(self, Cand):
    #     unl_Cand = [self.unlabeled.index(k) for k in Cand]
    #     sc = np.exp(-cdist(self.X[self.unlabeled, :], self.X[self.unlabeled,:])/self.sigma) \
    #                         - self.K_ul @ self.K_inv @ self.K_ul.T
    #     return (1. - np.absolute(self.f[Cand])) *  np.linalg.norm(sc[:,unl_Cand], axis=0) / np.abs(sc[unl_Cand, unl_Cand])





class RKHSClassifier(object):
    def __init__(self, X, sigma):
        self.N = X.shape[0]
        self.sigma = sigma
        self.f = None
        self.K = np.exp(-cdist(X, X)/self.sigma) # calculate full, dense kernel matrix upfront
        self.modelname = 'rkhs'
        self.nc = 2

    def calculate_model(self, labeled, y):
        self.labeled = labeled
        self.unlabeled = list(filter(lambda x: x not in self.labeled, range(self.N)))
        self.y = y

        self.K_inv = np.linalg.inv(self.K[np.ix_(self.labeled, self.labeled)])
        self.f = np.empty(self.N)
        self.f[self.labeled] = self.y
        self.K_ul = self.K[np.ix_(self.unlabeled, self.labeled)]
        self.f[self.unlabeled] = self.K_ul @ self.K_inv @ np.array(self.y) # Kul @ Kll^{-1} y

    def update_model(self, Q, yQ):
        unl_Q = [self.unlabeled.index(k) for k in Q]
        unl_notQ = list(filter(lambda x: x not in unl_Q, range(len(self.unlabeled))))
        notQ = list(filter(lambda x: x not in Q, self.unlabeled))

        # update f
        SQ_inv = np.linalg.inv(self.K[np.ix_(Q, Q)])
        K_notQ_Q = self.K[np.ix_(notQ, Q)]
        sc = K_notQ_Q - self.K_ul[unl_notQ,:] @ self.K_inv @ self.K_ul[unl_Q,:].T
        self.f[notQ] += sc @ SQ_inv @ (np.array(yQ) - self.f[Q])
        self.f[Q] = yQ

        # update self.K_inv, self.K_ul for future calculations
        # calculate self.K_inv new by block matrix inversion formula
        Mat = self.K_ul[unl_Q,:].T @ SQ_inv @ self.K_ul[unl_Q,:] @ self.K_inv
        Mat.ravel()[::len(self.labeled)+1] += 1.0
        new_Kinv_lower_left = -SQ_inv @ self.K_ul[unl_Q,:] @ self.K_inv
        bottom_new_Kinv = np.hstack((new_Kinv_lower_left, SQ_inv))
        self.K_inv = self.K_inv @ Mat
        self.K_inv = np.hstack((self.K_inv, new_Kinv_lower_left.T))
        self.K_inv = np.vstack((self.K_inv, bottom_new_Kinv))

        # self.K_ul
        self.K_ul = self.K_ul[unl_notQ,:]
        self.K_ul = np.hstack((self.K_ul, K_notQ_Q))

        self.labeled.extend(Q)
        self.unlabeled = notQ

        return


    def look_ahead_db_norms(self, Cand):
        unl_Cand = [self.unlabeled.index(k) for k in Cand]
        KCandl = self.K_ul[unl_Cand,:]
        sc = self.K[np.ix_(self.unlabeled, Cand)] - self.K_ul @ self.K_inv @ KCandl.T
        return (1. - np.absolute(self.f[Cand]))*np.linalg.norm(sc, axis=0)/sc[unl_Cand, range(len(Cand))]

    def return_parts(self, Cand):
        unl_Cand = [self.unlabeled.index(k) for k in Cand]
        KCandl = self.K_ul[unl_Cand,:]
        sc = self.K[np.ix_(self.unlabeled, Cand)] - self.K_ul @ self.K_inv @ KCandl.T
        return (1. - np.absolute(self.f[Cand])), np.linalg.norm(sc, axis=0), 1./sc[unl_Cand, range(len(Cand))]



if __name__ == "__main__":
    import sys
    sys.path.append('..')
    import os
    from runs_util import get_acc
    parser = ArgumentParser(description="Read in previous RKHS run and check classifier")
    parser.add_argument("--loc", default='../checker2/db-rkhs-2000-0.1-0.1/rand-top-5-100-1.txt', type=str)
    parser.add_argument("--Xloc", default='../checker2/X_labels.npz', type=str)
    parser.add_argument("--sigma", default='0.1', type=str)
    args = parser.parse_args()
    print(float(args.sigma))
    labeled = []
    with open(args.loc, 'r') as f:
        for i, line in enumerate(f.readlines()):
            # read in init_labeled, and initial accuracy
            line = line.split(',')
            labeled.extend([int(x) for x in line[:-2]])
            if i == 0:
                num_init = len(labeled)

    lab_set = set(labeled)
    print(len(lab_set), len(labeled))
    data = np.load(args.Xloc, allow_pickle=True)
    X, labels = data['X'], data['labels']

    model = RKHSClassifier(X, sigma=float(args.sigma))
    model_new = RKHSClassifierNew(X, sigma=float(args.sigma))


    model.calculate_model(labeled[:20], list(labels[labeled[:20]]))
    model_new.calculate_model(labeled[:20], list(labels[labeled[:20]]))

    assert np.allclose(model.f, model_new.f)

    Cand = list(np.random.choice(model.unlabeled, 50))
    orig_vals = model.look_ahead_db_norms(Cand[:])
    new_vals = model_new.look_ahead_db_norms(Cand[:])

    assert np.allclose(orig_vals, new_vals)

    model.update_model(Cand[:5], list(labels[Cand[:5]]))
    model_new.update_model(Cand[:5], list(labels[Cand[:5]]))

    assert np.allclose(model.f, model_new.f)
    print("passed all tests!")
    # for i in np.arange(0,100,10):
    #     K = num_init+i*5
    #     model = RKHSClassifier(X, sigma=float(args.sigma))
    #     model.calculate_model(labeled[:K], list(labels[labeled[:K]]))
    #     model2 = RKHSClassifier(X, sigma=float(args.sigma))
    #     model2.calculate_model(labeled[:K], list(labels[labeled[:K]]))
    #
    #
    #
    #     assert np.allclose(model.f, model2.f)
    #
    #     Q = list(np.random.choice(model.unlabeled, 5))
    #
    #     model3 = RKHSClassifier(X, sigma=float(args.sigma))
    #     model3.calculate_model(labeled[:K] + Q, list(labels[labeled[:K]]) + list(labels[Q]))
    #
    #     #print(list(labels[Q]))
    #     model.update_model_old(Q, list(labels[Q]))
    #     model2.update_model(Q, list(labels[Q]))
    #
    #     if np.allclose(model.f, model2.f):
    #         print("both models are close in the update")
    #     else:
    #         print("true to old - " + str(np.allclose(model.f, model3.f)))
    #         print("true to new - " + str(np.allclose(model2.f, model3.f)))
    #         print(np.linalg.norm(model.f - model2.f)/np.linalg.norm(model.f))
    #         # print(labels[model.labeled], model.f[model.labeled])
    #         # print(labels[model2.labeled], model.f[model2.labeled])
    #         plt.scatter(range(model.N), model3.f, label='true', marker='^')
    #         plt.scatter(range(model.N), model.f, label='old', marker='x')
    #         plt.scatter(range(model.N), model2.f, label='new', marker='.')
    #
    #         plt.legend()
    #         plt.savefig('./comp-%d.png' % i)
    #         plt.show(0)
    #         plt.close()
    #         print()
    #     #assert np.allclose(model.f, model2.f)
    #
    #     # print(os.path.exists('rkhs-model-0.npz'))
    #     # data = np.load('rkhs-model-%d.npz' % i)
    #     # print(type(data))
    #     # print(list(data.keys()))
    #     # saved_f, saved_lab = data['f'], data['lab']
    #     # print(model.labeled)
    #     # print(saved_lab)
    #     # print(np.allclose(saved_f, model.f))
    #     # print(K, len(model.labeled), get_acc(model.f, labels, unlabeled=model.unlabeled)[1], get_acc(saved_f, labels, unlabeled=model.unlabeled)[1])
    #     # fig, (ax1, ax2) = plt.subplots(1,2)
    #     # ax1.scatter(X[:, 0], X[:, 1], c=labels)
    #     # ax1.set_title("Ground Truth")
    #     # ax2.scatter(X[:,0], X[:,1], c=np.sign(model.f))
    #     # ax2.scatter(X[labeled[:K],0], X[labeled[:K],1], c='k', marker='^')
    #     # ax2.set_title("Calculated -- acc = {:.4f}".format(get_acc(model.f, labels, unlabeled=model.unlabeled)[1]))
    #     # plt.savefig('./check2-rkhs-%d.png' % K)
    #     # plt.show(0)