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

import numpy as np
import scipy as sp
from scipy.optimize import root
import scipy.linalg as sla
from scipy.stats import norm



def get_acc(u, labels, unlabeled = None):
    u_ = np.sign(u)
    u_[u_ == 0] = 1
    if unlabeled is None:
        corr = sum(1.*(u_ == labels))
        return corr, corr/u.shape[0]
    else:
        corr = sum(1.*(u_[unlabeled] == labels[unlabeled]))
        return corr, corr/len(unlabeled)
def get_acc_multi(u, labels, unlabeled=None):
    """
    Assuming that u and labels are NOT in one-hot encoding. i.e. ith entry of
    u and labels is the integer class in {0,1,2,...num_classes}.
    """
    if unlabeled is None:
        corr = sum(1.*(u == labels))
        return corr, corr/u.shape[0]
    else:
        corr = sum(1.*(u[unlabeled] == labels[unlabeled]))
        return corr, corr/len(unlabeled)


############# Gaussian Regression functions (Full Storage Model)##############

def gr_C(Z, gamma, d, v):
    H_d = len(Z) *[1./gamma**2.]
    vZ = v[Z,:]
    post = sp.sparse.diags(d, format='csr') \
           + vZ.T @ sp.sparse.diags(H_d, format='csr') @ vZ
    return v @ sp.linalg.inv(post) @ v.T

def gr_map(Z, y, gamma, d, v):
    C = gr_C(Z, gamma, d, v)
    return (C[:,Z] @ y)/(gamma * gamma)



########### ProbitNorm Jacobian and Hessian entry Calculations ##########
def jac_calc(uj_, yj_, gamma):
    return -yj_*(norm.pdf(uj_*yj_, scale=gamma)/norm.cdf(uj_*yj_, scale=gamma))

def hess_calc(uj_, yj_, gamma):
     return (norm.pdf(uj_*yj_,scale=gamma)**2 - pdf_deriv(uj_*yj_, gamma)*norm.cdf(uj_*yj_,scale=gamma))/(norm.cdf(uj_*yj_,scale=gamma)**2)



########### (Full Storage Model) Functions #############
def probit_map_dr(Z_, yvals, gamma, Ct):
    """
    Probit MAP estimator, using dimensionality reduction via Representer Theorem.
    *** This uses cdf of normal distribution ***
    """
    Ctp = Ct[np.ix_(Z_,Z_)]
    Jj = len(yvals)

    def f(x):
        vec = [jac_calc(x[j], yj, gamma) for j,yj in enumerate(yvals)]
        if np.any(vec == np.inf):
            print('smally in f')
        return x + Ctp @ vec

    def fprime(x):
        H = Ctp * np.array([hess_calc(x[j],yj, gamma)
                            for j, yj in enumerate(yvals)])
        if np.any(H == np.inf):
            print('smally')
        H[np.diag_indices(Jj)] += 1.0
        return H

    x0 = np.random.rand(Jj)
    x0[np.array(yvals) < 0] *= -1
    res = root(f, x0, jac=fprime)
    #print(np.allclose(0., f(res.x)))
    tmp = sla.inv(Ctp) @ res.x
    return Ct[:, Z_] @ tmp


def Hess_inv(u, Z, y, gamma, Ct):
    Ctp = Ct[np.ix_(Z, Z)]
    H_d = np.zeros(len(y))
    for i, (j, yj) in enumerate(zip(Z, y)):
        H_d[i] = 1./hess_calc(u[j], yj, gamma)
    temp = sp.linalg.inv(sp.sparse.diags(H_d, format='csr') + Ctp)
    return Ct - Ct[:, Z] @ temp @ Ct[Z, :]

def Hess_inv_st(u, Z, y, w, v, gamma, debug=False):
    H_d = np.zeros(len(Z))
    vZ = v[Z,:]
    for i, (j, yj) in enumerate(zip(Z, y)):
        H_d[i] = hess_calc2(u[j], yj, gamma)
    post = sp.sparse.diags(w, format='csr') \
           + vZ.T @ sp.sparse.diags(H_d, format='csr') @ vZ
    return v @ sp.linalg.inv(post) @ v.T




################ Logistic and Hessian entry Calculations##################
def jac_calc2(uj_, yj_, gamma):
    return -yj_*(np.exp(-uj_*yj_/gamma))/(gamma*(1.0 + np.exp(-uj_*yj_/gamma)))

def hess_calc2(uj_, yj_, gamma):
    return np.exp(-uj_*yj_/gamma)/(gamma**2. * (1. + np.exp(-uj_*yj_/gamma))**2.)


########### (Full Storage Model) Functions #############
def probit_map_dr2(Z_, yvals, gamma, Ct):
    """
    Probit MAP estimator, using dimensionality reduction via Representer Theorem.
    *** This uses logistic cdf ***
    """
    Ctp = Ct[np.ix_(Z_,Z_)]
    Jj = len(yvals)

    def f(x):
        vec = [jac_calc2(x[j], yj, gamma)
               for j,yj in enumerate(yvals)]
        return x + Ctp @ vec

    def fprime(x):
        H = Ctp * np.array([hess_calc2(x[j], yj, gamma)
                            for j, yj in enumerate(yvals)])
        if np.any(H == np.inf):
            print('smally')
        H[np.diag_indices(Jj)] += 1.0
        return H

    x0 = np.random.rand(Jj)
    x0[np.array(yvals) < 0] *= -1
    res = root(f, x0, jac=fprime)
    tmp = sla.inv(Ctp) @ res.x
    return Ct[:, Z_] @ tmp


def Hess2_inv(u, Z, y, gamma, Ct):
    Ctp = Ct[np.ix_(Z, Z)]
    H_d = np.zeros(len(y))
    for i, (j, yj) in enumerate(zip(Z, y)):
        H_d[i] = 1./hess_calc2(u[j], yj, gamma)
    temp = sp.linalg.inv(sp.sparse.diags(H_d, format='csr') + Ctp)
    return Ct - Ct[:, Z] @ temp @ Ct[Z, :]


def Hess_inv_st2(u, Z, y, w, v, gamma, debug=False):
    H_d = np.zeros(len(Z))
    vZ = v[Z,:]
    for i, (j, yj) in enumerate(zip(Z, y)):
        H_d[i] = hess_calc2(u[j], yj, gamma)
    post = sp.sparse.diags(w, format='csr') \
           + vZ.T @ sp.sparse.diags(H_d, format='csr') @ vZ
    return v @ sp.linalg.inv(post) @ v.T






##################################################################################
################### Helper Functions for Reduced Model ###########################
##################################################################################

def Hess_inv_st_alpha(alpha, y, d, v_Z, gamma):
    '''
    This method keeps everything in the "alpha"-space, of dimension num_eig x num_eig
    '''
    Dprime = sp.sparse.diags([1./hess_calc(np.inner(v_Z[i,:],alpha), yi, gamma) for i, yi in enumerate(y)], format='csr')
    A = d[:, np.newaxis] * ((v_Z.T @ sp.linalg.inv(Dprime + v_Z @ (d[:, np.newaxis] * (v_Z.T))) @ v_Z) * d)
    return np.diag(d) - A

def Hess_inv_st2_alpha(alpha, y, d, v_Z, gamma):
    '''
    This method keeps everything in the "alpha"-space, of dimension num_eig x num_eig
    '''
    Dprime = sp.sparse.diags([1./hess_calc2(np.inner(v_Z[i,:],alpha), yi, gamma) for i, yi in enumerate(y)], format='csr')
    A = d[:, np.newaxis] * ((v_Z.T @ sp.linalg.inv(Dprime + v_Z @ (d[:, np.newaxis] * (v_Z.T))) @ v_Z) * d)
    return np.diag(d) - A


def probit_map_st_alpha(Z, y,  gamma, w, v):
    n,l = v.shape[1], len(y)
    def f(x):
        vec = np.zeros(l)
        for j, (i,yi) in enumerate(zip(Z,y)):
            vec[j] = - jac_calc(np.inner(v[i,:], x), yi, gamma)
        return w * x  - v[Z,:].T @ vec
    def fprime(x):
        vec = np.zeros(l)
        for j, (i,yi) in enumerate(zip(Z, y)):
            vec[j] = -hess_calc(np.inner(v[i,:], x), yi, gamma)

        H = (-v[Z,:].T * vec) @ v[Z,:]
        H[np.diag_indices(n)] += w
        return H
    x0 = np.random.rand(len(w))
    res = root(f, x0, jac=fprime)
    #print(f"Root Finding is successful: {res.success}")
    return res.x


def probit_map_st2_alpha(Z, y,  gamma, w, v):
    n,l = v.shape[1], len(y)
    def f(x):
        vec = np.zeros(l)
        for j, (i,yi) in enumerate(zip(Z,y)):
            vec[j] = - jac_calc2(np.inner(v[i,:], x), yi, gamma)
        return w * x  - v[Z,:].T @ vec
    def fprime(x):
        vec = np.zeros(l)
        for j, (i,yi) in enumerate(zip(Z, y)):
            vec[j] = -hess_calc2(np.inner(v[i,:], x), yi, gamma)

        H = (-v[Z,:].T * vec) @ v[Z,:]
        H[np.diag_indices(n)] += w
        return H
    x0 = np.random.rand(len(w))
    res = root(f, x0, jac=fprime)
    #print(f"Root Finding is successful: {res.success}")
    return res.x







# class Classifier(object):
#     def __init__(self, name, gamma, tau, v=None, w=None, Ct=None):
#         self.gamma = gamma
#         self.tau = tau
#         self.v = v
#         self.w = w
#         if Ct is not None:
#             self.Ct = Ct
#         else:
#             self.d = (self.tau ** (2.)) * ((self.w + self.tau**2.) ** (-1.))
#             self.Ct = (self.v * self.d) @ self.v.T
#         self.name = name
#         return
#
#     def get_m(self, Z, y):
#         if self.name == "probit":
#             if len(y) <= len(self.w):
#                 return probit_map_dr(Z, y, self.gamma, self.Ct)
#             else:
#                 return probit_map_st(Z, y, self.gamma, 1./self.d, self.v)
#         elif self.name == "probit2":
#             if len(y) <= len(self.w):
#                 return probit_map_dr2(Z, y, self.gamma, self.Ct)
#             else:
#                 return probit_map_st2(Z, y, self.gamma, 1./self.d, self.v)
#         elif self.name == "gr":
#             #return gr_map(Z, y, self.gamma, self.Ct)
#             return gr_map(Z, y, self.gamma, 1./self.d, self.v)
#         else:
#             pass
#
#     def get_C(self, Z, y, m):
#         if self.name in ["probit", "probit-st"]:
#             if len(y) > len(self.w):
#                 return Hess_inv_st(m, Z, y, 1./self.d, self.v, self.gamma)
#             else:
#                 return Hess_inv(m, Z, y, self.gamma, self.Ct)
#         elif self.name in ["probit2", "probit2-st"]:
#             if len(y) > len(self.w):
#                 return Hess_inv_st2(m, Z, y, 1./self.d, self.v, self.gamma)
#             else:
#                 return Hess2_inv(m, Z, y, self.gamma, self.Ct)
#         elif self.name in ["gr"]:
#             #return gr_C(Z, self.gamma, self.Ct)
#             return gr_C(Z, self.gamma, 1./self.d, self.v)
#         else:
#             pass
#
#
# class Classifier_HF(object):
#     def __init__(self, tau, L):
#         self.tau = tau
#         if sps.issparse(L):
#             self.sparse = True
#             self.L = L + self.tau**2. * sps.eye(L.shape[0])
#         else:
#             self.sparse = False
#             self.L = L + self.tau**2. * np.eye(L.shape[0])
#         return
#
#     def get_m(self, Z, y):
#         Zbar = list(filter(lambda x: x not in Z, range(self.L.shape[0])))
#         self.m = -self.get_C(Z, Zbar=Zbar) @ self.L[np.ix_(Zbar, Z)] @ y
#         return self.m
#
#     def get_C(self, Z, Zbar=None):
#         if Zbar is None:
#             Zbar = list(filter(lambda x: x not in Z, range(self.L.shape[0])))
#         if self.sparse:
#             self.C = scipy.sparse.linalg.inv(self.L[np.ix_(Zbar, Zbar)]).toarray() # inverse will be dense anyway
#             return self.C
#         else:
#             self.C = sla.inv(self.L[np.ix_(Zbar, Zbar)])
#             return self.C
#
#
# ###############################################################################
# ############### Gaussian Regression Helper Functions ##########################
#
#
# def get_init_post(C_inv, labeled, gamma2):
#     """
#     calculate the risk of each unlabeled point
#     C_inv: prior inverse (i.e. graph Laplacian L)
#     """
#     N = C_inv.shape[0]
#     #unlabeled = list(filter(lambda x: x not in labeled, range(N)))
#     B_diag = [1 if i in labeled else 0 for i in range(N)]
#     B = sp.sparse.diags(B_diag, format='csr')
#     return sp.linalg.inv(C_inv + B/gamma2)
#
#
# def calc_next_m(m, C, y, lab, k, y_k, gamma2):
#     ck = C[k,:]
#     ckk = ck[k]
#     ip = np.dot(ck[lab], y[lab])
#     val = ((gamma2)*y_k -ip )/(gamma2*(gamma2 + ckk))
#     m_k = m + val*ck
#     return m_k
#
#
# def get_probs(m, sigmoid=False):
#     if sigmoid:
#         return 1./(1. + np.exp(-3.*m))
#     m_probs = m.copy()
#     # simple fix to get probabilities that respect the 0 threshold
#     m_probs[np.where(m_probs >0)] /= 2.*np.max(m_probs)
#     m_probs[np.where(m_probs <0)] /= -2.*np.min(m_probs)
#     m_probs += 0.5
#     return m_probs
#
#
# def EEM_full(k, m, C, y, lab, unlab, m_probs, gamma2):
#     N = C.shape[0]
#     m_at_k = m_probs[k]
#     m_k_p1 = calc_next_m(m, C, y, lab, k, 1., gamma2)
#     m_k_p1 = get_probs(m_k_p1)
#     risk = m_at_k*np.sum([min(m_k_p1[i], 1.- m_k_p1[i]) for i in range(N)])
#     m_k_m1 = calc_next_m(m, C, y, lab, k, -1., gamma2)
#     m_k_m1 = get_probs(m_k_m1)
#     risk += (1.-m_at_k)*np.sum([min(m_k_m1[i], 1.- m_k_m1[i]) for i in range(N)])
#     return risk
#
#
# def EEM_opt_record(m, C, y, labeled, unlabeled, gamma2):
#     m_probs = get_probs(m)
#     N = C.shape[0]
#     risks = [EEM_full(j, m, C, y, labeled, unlabeled, m_probs, gamma2) for j in range(N)]
#     k = np.argmin(risks)
#     return k, risks
#
#
# def Sigma_opt(C, unlabeled, gamma2):
#     sums = np.sum(C[np.ix_(unlabeled,unlabeled)], axis=1)
#     sums = np.asarray(sums).flatten()**2.
#     s_opt = sums/(gamma2 + np.diag(C)[unlabeled])
#     k_max = unlabeled[np.argmax(s_opt)]
#     return k_max
#
#
# def V_opt_record(C, unlabeled, gamma2):
#     ips = np.array([np.inner(C[k,:], C[k,:]) for k in unlabeled]).flatten()
#     v_opt = ips/(gamma2 + np.diag(C)[unlabeled])
#     k_max = unlabeled[np.argmax(v_opt)]
#     return k_max, v_opt
#
#
# def V_opt_record2(u, C, unlabeled, gamma2, lam):
#     ips = np.array([np.inner(C[k,:], C[k,:]) for k in unlabeled]).flatten()
#     v_opt = ips/(gamma2 + np.diag(C)[unlabeled])
#     print(np.max(v_opt), np.min(v_opt))
#     v_opt += lam*(np.max(v_opt) + np.min(v_opt))*0.5*(1./np.absolute(u[unlabeled])) # add term that bias toward decision boundary
#     k_max = unlabeled[np.argmax(v_opt)]
#     return k_max, v_opt
#
#
# def Sigma_opt_record(C, unlabeled, gamma2):
#     sums = np.sum(C[np.ix_(unlabeled,unlabeled)], axis=1)
#     sums = np.asarray(sums).flatten()**2.
#     s_opt = sums/(gamma2 + np.diag(C)[unlabeled])
#     k_max = unlabeled[np.argmax(s_opt)]
#     return k_max, s_opt
#
#
#
#
# ################################################################################
# ################## Plotting Helper Functions ###################################
#
#
# def plot_iter(m, X, labels, labeled, k_next=-1, title=None, subplot=False):
#     '''
#     Assuming labels are +1, -1
#     '''
#     m1 = np.where(m >= 0)[0]
#     m2 = np.where(m < 0)[0]
#
#     #sup1 = list(set(labeled).intersection(set(np.where(labels == 1)[0])))
#     #sup2 = list(set(labeled).intersection(set(np.where(labels == -1)[0])))
#
#     corr1 = list(set(m1).intersection(set(np.where(labels == 1)[0])))
#     incorr1 = list(set(m2).intersection(set(np.where(labels == 1)[0])))
#     corr2 = list(set(m2).intersection(set(np.where(labels == -1)[0])))
#     incorr2 = list(set(m1).intersection(set(np.where(labels == -1)[0])))
#
#     print("\tnum incorrect = %d" % (len(incorr1) + len(incorr2)))
#
#     if k_next >= 0:
#         plt.scatter(X[k_next,0], X[k_next,1], marker= 's', c='y', alpha= 0.7, s=200) # plot the new points to be included
#         plt.title(r'Dataset with Label for %s added' % str(k_next))
#     elif k_next == -1:
#         if not title:
#             plt.title(r'Dataset with Initial Labeling')
#         else:
#             plt.title(title)
#
#     plt.scatter(X[corr1,0], X[corr1,1], marker='x', c='b', alpha=0.2)
#     plt.scatter(X[incorr1,0], X[incorr1,1], marker='x', c='r', alpha=0.2)
#     plt.scatter(X[corr2,0], X[corr2,1], marker='o', c='r',alpha=0.15)
#     plt.scatter(X[incorr2,0], X[incorr2,1], marker='o', c='b',alpha=0.15)
#     # plt.scatter(X[corr1,0], X[corr1,1], marker='x', c='b', alpha=0.8)
#     # plt.scatter(X[incorr1,0], X[incorr1,1], marker='x', c='r', alpha=0.8)
#     # plt.scatter(X[corr2,0], X[corr2,1], marker='o', c='r',alpha=0.8)
#     # plt.scatter(X[incorr2,0], X[incorr2,1], marker='o', c='b',alpha=0.8)
#
#     sup1 = list(set(labeled).intersection(set(np.where(labels == 1)[0])))
#     sup2 = list(set(labeled).intersection(set(np.where(labels == -1)[0])))
#     plt.scatter(X[sup1,0], X[sup1,1], marker='x', c='k', alpha=1.0)
#     plt.scatter(X[sup2,0], X[sup2,1], marker='o', c='k', alpha=1.0)
#     plt.axis('equal')
#     if subplot:
#         return
#     plt.show()
#     return








#
# ################################################################################
# ########## MAP estimators of Probit model (with normal distribution's pdf )
# ################################################################################
#
#
def pdf_deriv(t, gamma):
    return -t*norm.pdf(t, scale = gamma)/(gamma**2)
#
# def jac_calc(uj_, yj_, gamma):
#     return -yj_*(norm.pdf(uj_*yj_, scale=gamma)/norm.cdf(uj_*yj_, scale=gamma))
#
#
# def hess_calc(uj_, yj_, gamma):
#     return (norm.pdf(uj_*yj_,scale=gamma)**2 - pdf_deriv(uj_*yj_, gamma)*norm.cdf(uj_*yj_,scale=gamma))/(norm.cdf(uj_*yj_,scale=gamma)**2)
#
# def Hess(u, y, labeled, Lt, gamma, debug=False):
#     """
#     Assuming matrices are sparse, since L_tau should be relatively sparse,
#         and we are perturbing with diagonal matrix.
#     """
#     H_d = np.zeros(u.shape[0])
#     for j, yj in zip(labeled, y):
#         H_d[j] = hess_calc(u[j], yj, gamma)
#     if debug:
#         print(H_d[np.nonzero(H_d)])
#     # if np.any(H_d == np.inf):
#     #     print('smally')
#     return Lt + sp.sparse.diags(H_d, format='csr')
#
#
# def J(u, y, labeled, Lt, gamma, debug=False):
#     vec = np.zeros(u.shape[0])
#     for j, yj in zip(labeled,y):
#         vec[j] = -yj*norm.pdf(u[j]*yj, gamma)/norm.cdf(u[j]*yj, gamma)
#     if debug:
#         print(vec[np.nonzero(vec)])
#     return Lt @ u + vec
#
#
# def probit_map_dr(Z_, yvals, gamma, Ct):
#     """
#     Probit MAP estimator, using dimensionality reduction via Representer Theorem.
#     *** This uses cdf of normal distribution ***
#     """
#     Ctp = Ct[np.ix_(Z_,Z_)]
#     Jj = len(yvals)
#
#     def f(x):
#         vec = [jac_calc(x[j], yj, gamma) for j,yj in enumerate(yvals)]
#         if np.any(vec == np.inf):
#             print('smally in f')
#         return x + Ctp @ vec
#
#     def fprime(x):
#         H = Ctp * np.array([hess_calc(x[j],yj, gamma)
#                             for j, yj in enumerate(yvals)])
#         if np.any(H == np.inf):
#             print('smally')
#         H[np.diag_indices(Jj)] += 1.0
#         return H
#
#     x0 = np.random.rand(Jj)
#     x0[np.array(yvals) < 0] *= -1
#     res = root(f, x0, jac=fprime)
#     #print(np.allclose(0., f(res.x)))
#     tmp = sla.inv(Ctp) @ res.x
#     return Ct[:, Z_] @ tmp
#
#
#
# ################################################################################
# ################### Probit with logit likelihood
# ################################################################################
#
# ###################### Psi = cdf of logistic #######################
# def log_pdf(t, g):
#     return np.exp(-t/g)/(g*(1. + np.exp(-t/g)**2.))
#
#
# def log_cdf(t, g):
#     return 1.0/(1.0 + np.exp(-t/g))
#
#
# def log_pdf_deriv(t, g):
#     return -np.exp(-t/g)/((g*(1. + np.exp(-t/g)))**2.)
#
# def jac_calc2(uj_, yj_, gamma):
#     return -yj_*(np.exp(-uj_*yj_/gamma))/(gamma*(1.0 + np.exp(-uj_*yj_/gamma)))
#
# def hess_calc2(uj_, yj_, gamma):
#     return np.exp(-uj_*yj_/gamma)/(gamma**2. * (1. + np.exp(-uj_*yj_/gamma))**2.)
#
#
# def Hess2(u, y, labeled, Lt, gamma, debug=False):
#     """
#     Assuming matrices are sparse, since L_tau should be relatively sparse,
#         and we are perturbing with diagonal matrix.
#     """
#     H_d = np.zeros(u.shape[0])
#     for j, yj in zip(labeled, y):
#         H_d[j] = hess_calc2(u[j], yj, gamma)
#     if debug:
#         print(H_d[np.nonzero(H_d)])
#     if np.any(H_d == np.inf):
#         print('smally')
#     return Lt + sp.sparse.diags(H_d, format='csr')
#
#
# def J2(u, y, labeled, Lt, gamma, debug=False):
#     vec = np.zeros(u.shape[0])
#     for j, yj in zip(labeled,y):
#         vec[j] = jac_calc2(u[j], yj, gamma)
#     if debug:
#         print(vec[np.nonzero(vec)])
#     return Lt @ u + vec
#
#
# def probit_map_dr2(Z_, yvals, gamma, Ct):
#     """
#     Probit MAP estimator, using dimensionality reduction via Representer Theorem.
#     *** This uses logistic cdf ***
#     """
#     Ctp = Ct[np.ix_(Z_,Z_)]
#     Jj = len(yvals)
#
#     def f(x):
#         vec = [jac_calc2(x[j], yj, gamma)
#                for j,yj in enumerate(yvals)]
#         return x + Ctp @ vec
#
#     def fprime(x):
#         H = Ctp * np.array([hess_calc2(x[j], yj, gamma)
#                             for j, yj in enumerate(yvals)])
#         if np.any(H == np.inf):
#             print('smally')
#         H[np.diag_indices(Jj)] += 1.0
#         return H
#
#     x0 = np.random.rand(Jj)
#     x0[np.array(yvals) < 0] *= -1
#     res = root(f, x0, jac=fprime)
#     tmp = sla.inv(Ctp) @ res.x
#     return Ct[:, Z_] @ tmp
#
#
#
# ################################################################################
# # Spectral truncation
# ################################################################################
# """
#
# def pdf_deriv(t, gamma):
#     return -t*norm.pdf(t, scale = gamma)/(gamma**2)
#
# def jac_calc(uj_, yj_, gamma):
#     return -yj_*(norm.pdf(uj_*yj_, scale=gamma)/norm.cdf(uj_*yj_, scale=gamma))
#
#
# def hess_calc(uj_, yj_, gamma):
#     return (norm.pdf(uj_*yj_,scale=gamma)**2 - pdf_deriv(uj_*yj_, gamma)*norm.cdf(uj_*yj_,scale=gamma))/(norm.cdf(uj_*yj_,scale=gamma)**2)
#
# """
# def probit_map_st(Z, y, gamma, w, v):
#     N = v.shape[0]
#     n = v.shape[1]
#     def f(x):
#         vec = np.zeros(N)
#         tmp = v @ x
#         for i, yi in zip(Z, y):
#             vec[i] = - jac_calc(tmp[i], yi, gamma)
#             #vec[i] = yi*norm.pdf(tmp[i]*yi, scale=gamma)/norm.cdf(tmp[i]*yi, scale=gamma)
#         return w * x  - v.T @ vec
#     def fprime(x):
#         tmp = v @ x
#         vec = np.zeros(N)
#         for i, yi in zip(Z, y):
#             vec[i] = -hess_calc(tmp[i], yi, gamma)
#             #vec[i] = (pdf_deriv(tmp[i]*yi,gamma)*norm.cdf(tmp[i]*yi,scale=gamma)
#             #             - norm.pdf(tmp[i]*yi,scale=gamma)**2)/ (norm.cdf(tmp[i]*yi,scale=gamma)**2)
#         H = (-v.T * vec) @ v
#         H[np.diag_indices(n)] += w
#         return H
#     x0 = np.random.rand(len(w))
#     res = root(f, x0, jac=fprime)
#     #print(f"Root Finding is successful: {res.success}")
#     return v @ res.x
#
# # OLD Code
# # def Hess_inv_st(u, Z, y, w, v, gamma, debug=False):
# #     """
# #     Assuming matrices are sparse, since L_tau should be relatively sparse,
# #         and we are perturbing with diagonal matrix.
# #     """
# #     H_d = np.zeros(u.shape[0])
# #     for j, yj in zip(Z, y):
# #         H_d[j] = hess_calc(u[j], yj, gamma)
# #     post = sp.sparse.diags(w, format='csr') \
# #            + v.T @ sp.sparse.diags(H_d, format='csr') @ v
# #     return v @ sp.linalg.inv(post) @ v.T
#
# def Hess_inv_st(u, Z, y, w, v, gamma, debug=False):
#     H_d = np.zeros(len(Z))
#     vZ = v[Z,:]
#     for i, (j, yj) in enumerate(zip(Z, y)):
#         H_d[i] = hess_calc2(u[j], yj, gamma)
#     post = sp.sparse.diags(w, format='csr') \
#            + vZ.T @ sp.sparse.diags(H_d, format='csr') @ vZ
#     return v @ sp.linalg.inv(post) @ v.T
#
#
#
# def probit_map_st2(Z, y,  gamma, w, v):
#     N = v.shape[0]
#     n = v.shape[1]
#     def f(x):
#         vec = np.zeros(N)
#         tmp = v @ x
#         for i, yi in zip(Z, y):
#             vec[i] = - jac_calc2(tmp[i], yi, gamma)
#             #vec[i] = yi*norm.pdf(tmp[i]*yi, scale=gamma)/norm.cdf(tmp[i]*yi, scale=gamma)
#         return w * x  - v.T @ vec
#     def fprime(x):
#         tmp = v @ x
#         vec = np.zeros(N)
#         for i, yi in zip(Z, y):
#             vec[i] = -hess_calc2(tmp[i], yi, gamma)
#             #vec[i] = (pdf_deriv(tmp[i]*yi,gamma)*norm.cdf(tmp[i]*yi,scale=gamma)
#             #             - norm.pdf(tmp[i]*yi,scale=gamma)**2)/ (norm.cdf(tmp[i]*yi,scale=gamma)**2)
#         H = (-v.T * vec) @ v
#         H[np.diag_indices(n)] += w
#         return H
#     x0 = np.random.rand(len(w))
#     res = root(f, x0, jac=fprime)
#     #print(f"Root Finding is successful: {res.success}")
#     return v @ res.x
#
# # OLD code
# # def Hess_inv_st2(u, Z, y, w, v, gamma, debug=False):
# #     """
# #     Assuming matrices are sparse, since L_tau should be relatively sparse,
# #         and we are perturbing with diagonal matrix.
# #     """
# #     H_d = np.zeros(u.shape[0])
# #     for j, yj in zip(Z, y):
# #         H_d[j] = hess_calc2(u[j], yj, gamma)
# #     post = sp.sparse.diags(w, format='csr') \
# #            + v.T @ sp.sparse.diags(H_d, format='csr') @ v
# #     return v @ sp.linalg.inv(post) @ v.T
#
# def Hess_inv_st2(u, Z, y, w, v, gamma, debug=False):
#     H_d = np.zeros(len(Z))
#     vZ = v[Z,:]
#     for i, (j, yj) in enumerate(zip(Z, y)):
#         H_d[i] = hess_calc2(u[j], yj, gamma)
#     post = sp.sparse.diags(w, format='csr') \
#            + vZ.T @ sp.sparse.diags(H_d, format='csr') @ vZ
#     return v @ sp.linalg.inv(post) @ v.T
#
#
# def Hess_inv(u, Z, y, gamma, Ct):
#     Ctp = Ct[np.ix_(Z, Z)]
#     H_d = np.zeros(len(y))
#     for i, (j, yj) in enumerate(zip(Z, y)):
#         H_d[i] = 1./hess_calc(u[j], yj, gamma)
#     temp = sp.linalg.inv(sp.sparse.diags(H_d, format='csr') + Ctp)
#     return Ct - Ct[:, Z] @ temp @ Ct[Z, :]
#
# def Hess2_inv(u, Z, y, gamma, Ct):
#     Ctp = Ct[np.ix_(Z, Z)]
#     H_d = np.zeros(len(y))
#     for i, (j, yj) in enumerate(zip(Z, y)):
#         H_d[i] = 1./hess_calc2(u[j], yj, gamma)
#     temp = sp.linalg.inv(sp.sparse.diags(H_d, format='csr') + Ctp)
#     return Ct - Ct[:, Z] @ temp @ Ct[Z, :]
#
# # OLD inefficient
# # def gr_C(Z, gamma, Ct):
# #     Ctp = Ct[np.ix_(Z, Z)]
# #     H_d = np.ones(len(Z)) * (gamma * gamma)
# #     temp = sp.linalg.inv(sp.sparse.diags(H_d, format='csr') + Ctp)
# #     return Ct - Ct[:, Z] @ temp @ Ct[Z, :]
#
# def gr_C(Z, gamma, d, v):
#     H_d = len(Z) *[1./gamma**2.]
#     vZ = v[Z,:]
#     post = sp.sparse.diags(d, format='csr') \
#            + vZ.T @ sp.sparse.diags(H_d, format='csr') @ vZ
#     return v @ sp.linalg.inv(post) @ v.T
#
# def gr_map(Z, y, gamma, d, v):
#     C = gr_C(Z, gamma, d, v)
#     return (C[:,Z] @ y)/(gamma * gamma)