import numpy as np
from scipy.stats import norm, laplace
import _pickle as pickle
from math import exp
import matplotlib
from collections import Counter
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import scipy.integrate as integrate

class calculate_accountant(object):

    def __init__(self, q, epoch_num, sigma_cl_1, sigma_cl_2, sigma_gd, sigma_dyn_clip, delta):
        self.l_max = 32

        # Sampling prob. for SGD
        self.q = q
        # epoch num for SGD
        self.epoch_num = epoch_num

        # noise for k-means cluster size
        self.sigma_cl_1 = sigma_cl_1
        # noise for k-means cluster centers
        self.sigma_cl_2 = sigma_cl_2
        # noise for Gradient descent
        self.sigma_gd = sigma_gd
        # noise for Adaptive clipping in each GD step
        self.sigma_dyn_clip = sigma_dyn_clip

        self.delta = delta
        # k-means iterations
        self.T_1 = 10
        # SGD iterations
        self.T_2 = epoch_num * int(1.0/q)
        # Total iterations
        self.T = self.T_1 + self.T_2

        # mu_0 = sigma* N(0,1)
        # mu_1 = sigma* N(0,1) + 1
        # mu = (1-q)* mu_0 + q*mu_1 = sigma*N(0,1) + q
        # Probability is not incorporated when the the expected values are computed



    def compute_gauss(self,p_sigma, l, u, d=20, p_sens=1.0):
        e_1 = lambda x: norm.pdf(x,0,p_sens*p_sigma)* (norm.pdf(x,0,p_sens*p_sigma) / ((1-self.q)*norm.pdf(x,0,p_sens*p_sigma) +
            self.q*norm.pdf(x,p_sens,p_sens*p_sigma)))**(u*l)
        e_2 = lambda x: ((1-self.q)*norm.pdf(x,0,p_sens*p_sigma) + self.q*norm.pdf(x,p_sens,p_sens*p_sigma))*\
                (((1-self.q)*norm.pdf(x,0,p_sens*p_sigma) + self.q*norm.pdf(x,p_sens,p_sens*p_sigma)) /
                        norm.pdf(x,0,p_sens*p_sigma))**(u*l)
        E_1, _ = integrate.quad(e_1,-d,d)
        E_1 = np.nan_to_num(E_1)
        E_2, _ = integrate.quad(e_2,-d,d)
        E_2 = np.nan_to_num(E_2)
        return np.log(np.maximum(np.abs(E_1), np.abs(E_2))+0.00000001)

    def compute_alpha(self, p_sigma, q, d):
        #print("Sigma:", p_sigma)
        #print("Precomputing...")
        alpha_values = []
        for l in range(1,self.l_max+1):
            e_1 = lambda x: norm.pdf(x,0,p_sigma)* (norm.pdf(x,0,p_sigma) / ((1-q)*norm.pdf(x,0,p_sigma) +
                q*norm.pdf(x,1.0,p_sigma)))**l
            e_2 = lambda x: ((1-q)*norm.pdf(x,0,p_sigma) + q*norm.pdf(x,1.0,p_sigma))*\
                (((1-q)*norm.pdf(x,0,p_sigma) + q*norm.pdf(x,1.0,p_sigma))/(norm.pdf(x,0,p_sigma)))**l

            E_1, _ = integrate.quad(e_1,-d,d)
            E_2, _ = integrate.quad(e_2,-d,d)
            E_1 = np.nan_to_num(E_1)
            E_2 = np.nan_to_num(E_2)
            val = np.log(max([abs(E_1), abs(E_2)])+0.00000000001)
            #print("l = ", l, " alpha = ", val)
            alpha_values.append(val)
        #print("done.")
        return alpha_values

    def get_privacy(self):
        ## Alpha computation


        # Clustering
        alpha_values_1 = self.compute_alpha(self.sigma_cl_1, 1.0, 100)
        alpha_values_2 = self.compute_alpha(self.sigma_cl_2, 1.0, 500)


        # DP SGD
        alpha_values_3 = []

        for l in range(1, self.l_max+1):
            tmp = []
            for j in range(1, 10):
                u1 = j/10.0
                u2 = (10-j) / 10.0
                alpha_value1 = self.compute_gauss(self.sigma_dyn_clip, l, 1.0/u1, d=50, p_sens=2.0**0.5)
                alpha_value2 = self.compute_gauss(self.sigma_gd, l, 1.0/u2, d=20)
                tmp.append(u1 * alpha_value1 + u2 * alpha_value2)
            val = min(tmp)
            #print("l = ", l, " alpha = ", val)
            alpha_values_3.append(val)

        ## Epsilon computation

        eps_values = [ ]
        cnt = Counter()
        for i in range(self.T):
            epsilon_values = []
            for l in range(1, self.l_max+1):
                # Alpha of Clustering
                a_1 = min(i+1, self.T_1)*(alpha_values_1[l-1] + alpha_values_2[l-1])
                # Alpha of GD
                a_2 = max(0, i+1 - self.T_1)*alpha_values_3[l-1]

                eps_t = (a_1 + a_2 - np.log(self.delta)) / float(l)
                epsilon_values.append(eps_t)

            idx = np.argmin(epsilon_values)
            cnt[idx]+=1
            eps_values.append(epsilon_values[idx])


        print(cnt)

        print("Epsilon after clustering:", eps_values[self.T_1])
        print("Epsilon after GD:", eps_values[self.T-1])

        eps_values = eps_values[self.T_1:]
        eps_values = np.array(eps_values)[np.arange(int(1.0/self.q)-1,self.T_2,int(1.0/self.q))]

        for i in range(self.epoch_num):
            print("Epoch:", i+1, "Eps:", eps_values[i])
        return eps_values


df86f981936f545ef4879a9e432fa2de04e35f77
<mask token>


class calculate_accountant(object):

    def __init__(self, q, epoch_num, sigma_cl_1, sigma_cl_2, sigma_gd,
        sigma_dyn_clip, delta):
        self.l_max = 32
        self.q = q
        self.epoch_num = epoch_num
        self.sigma_cl_1 = sigma_cl_1
        self.sigma_cl_2 = sigma_cl_2
        self.sigma_gd = sigma_gd
        self.sigma_dyn_clip = sigma_dyn_clip
        self.delta = delta
        self.T_1 = 10
        self.T_2 = epoch_num * int(1.0 / q)
        self.T = self.T_1 + self.T_2

    def compute_gauss(self, p_sigma, l, u, d=20, p_sens=1.0):
        e_1 = lambda x: norm.pdf(x, 0, p_sens * p_sigma) * (norm.pdf(x, 0, 
            p_sens * p_sigma) / ((1 - self.q) * norm.pdf(x, 0, p_sens *
            p_sigma) + self.q * norm.pdf(x, p_sens, p_sens * p_sigma))) ** (u *
            l)
        e_2 = lambda x: ((1 - self.q) * norm.pdf(x, 0, p_sens * p_sigma) + 
            self.q * norm.pdf(x, p_sens, p_sens * p_sigma)) * (((1 - self.q
            ) * norm.pdf(x, 0, p_sens * p_sigma) + self.q * norm.pdf(x,
            p_sens, p_sens * p_sigma)) / norm.pdf(x, 0, p_sens * p_sigma)) ** (
            u * l)
        E_1, _ = integrate.quad(e_1, -d, d)
        E_1 = np.nan_to_num(E_1)
        E_2, _ = integrate.quad(e_2, -d, d)
        E_2 = np.nan_to_num(E_2)
        return np.log(np.maximum(np.abs(E_1), np.abs(E_2)) + 1e-08)
    <mask token>

    def get_privacy(self):
        alpha_values_1 = self.compute_alpha(self.sigma_cl_1, 1.0, 100)
        alpha_values_2 = self.compute_alpha(self.sigma_cl_2, 1.0, 500)
        alpha_values_3 = []
        for l in range(1, self.l_max + 1):
            tmp = []
            for j in range(1, 10):
                u1 = j / 10.0
                u2 = (10 - j) / 10.0
                alpha_value1 = self.compute_gauss(self.sigma_dyn_clip, l, 
                    1.0 / u1, d=50, p_sens=2.0 ** 0.5)
                alpha_value2 = self.compute_gauss(self.sigma_gd, l, 1.0 /
                    u2, d=20)
                tmp.append(u1 * alpha_value1 + u2 * alpha_value2)
            val = min(tmp)
            alpha_values_3.append(val)
        eps_values = []
        cnt = Counter()
        for i in range(self.T):
            epsilon_values = []
            for l in range(1, self.l_max + 1):
                a_1 = min(i + 1, self.T_1) * (alpha_values_1[l - 1] +
                    alpha_values_2[l - 1])
                a_2 = max(0, i + 1 - self.T_1) * alpha_values_3[l - 1]
                eps_t = (a_1 + a_2 - np.log(self.delta)) / float(l)
                epsilon_values.append(eps_t)
            idx = np.argmin(epsilon_values)
            cnt[idx] += 1
            eps_values.append(epsilon_values[idx])
        print(cnt)
        print('Epsilon after clustering:', eps_values[self.T_1])
        print('Epsilon after GD:', eps_values[self.T - 1])
        eps_values = eps_values[self.T_1:]
        eps_values = np.array(eps_values)[np.arange(int(1.0 / self.q) - 1,
            self.T_2, int(1.0 / self.q))]
        for i in range(self.epoch_num):
            print('Epoch:', i + 1, 'Eps:', eps_values[i])
        return eps_values

<mask token>


class calculate_accountant(object):

    def __init__(self, q, epoch_num, sigma_cl_1, sigma_cl_2, sigma_gd,
        sigma_dyn_clip, delta):
        self.l_max = 32
        self.q = q
        self.epoch_num = epoch_num
        self.sigma_cl_1 = sigma_cl_1
        self.sigma_cl_2 = sigma_cl_2
        self.sigma_gd = sigma_gd
        self.sigma_dyn_clip = sigma_dyn_clip
        self.delta = delta
        self.T_1 = 10
        self.T_2 = epoch_num * int(1.0 / q)
        self.T = self.T_1 + self.T_2

    def compute_gauss(self, p_sigma, l, u, d=20, p_sens=1.0):
        e_1 = lambda x: norm.pdf(x, 0, p_sens * p_sigma) * (norm.pdf(x, 0, 
            p_sens * p_sigma) / ((1 - self.q) * norm.pdf(x, 0, p_sens *
            p_sigma) + self.q * norm.pdf(x, p_sens, p_sens * p_sigma))) ** (u *
            l)
        e_2 = lambda x: ((1 - self.q) * norm.pdf(x, 0, p_sens * p_sigma) + 
            self.q * norm.pdf(x, p_sens, p_sens * p_sigma)) * (((1 - self.q
            ) * norm.pdf(x, 0, p_sens * p_sigma) + self.q * norm.pdf(x,
            p_sens, p_sens * p_sigma)) / norm.pdf(x, 0, p_sens * p_sigma)) ** (
            u * l)
        E_1, _ = integrate.quad(e_1, -d, d)
        E_1 = np.nan_to_num(E_1)
        E_2, _ = integrate.quad(e_2, -d, d)
        E_2 = np.nan_to_num(E_2)
        return np.log(np.maximum(np.abs(E_1), np.abs(E_2)) + 1e-08)

    def compute_alpha(self, p_sigma, q, d):
        alpha_values = []
        for l in range(1, self.l_max + 1):
            e_1 = lambda x: norm.pdf(x, 0, p_sigma) * (norm.pdf(x, 0,
                p_sigma) / ((1 - q) * norm.pdf(x, 0, p_sigma) + q * norm.
                pdf(x, 1.0, p_sigma))) ** l
            e_2 = lambda x: ((1 - q) * norm.pdf(x, 0, p_sigma) + q * norm.
                pdf(x, 1.0, p_sigma)) * (((1 - q) * norm.pdf(x, 0, p_sigma) +
                q * norm.pdf(x, 1.0, p_sigma)) / norm.pdf(x, 0, p_sigma)) ** l
            E_1, _ = integrate.quad(e_1, -d, d)
            E_2, _ = integrate.quad(e_2, -d, d)
            E_1 = np.nan_to_num(E_1)
            E_2 = np.nan_to_num(E_2)
            val = np.log(max([abs(E_1), abs(E_2)]) + 1e-11)
            alpha_values.append(val)
        return alpha_values

    def get_privacy(self):
        alpha_values_1 = self.compute_alpha(self.sigma_cl_1, 1.0, 100)
        alpha_values_2 = self.compute_alpha(self.sigma_cl_2, 1.0, 500)
        alpha_values_3 = []
        for l in range(1, self.l_max + 1):
            tmp = []
            for j in range(1, 10):
                u1 = j / 10.0
                u2 = (10 - j) / 10.0
                alpha_value1 = self.compute_gauss(self.sigma_dyn_clip, l, 
                    1.0 / u1, d=50, p_sens=2.0 ** 0.5)
                alpha_value2 = self.compute_gauss(self.sigma_gd, l, 1.0 /
                    u2, d=20)
                tmp.append(u1 * alpha_value1 + u2 * alpha_value2)
            val = min(tmp)
            alpha_values_3.append(val)
        eps_values = []
        cnt = Counter()
        for i in range(self.T):
            epsilon_values = []
            for l in range(1, self.l_max + 1):
                a_1 = min(i + 1, self.T_1) * (alpha_values_1[l - 1] +
                    alpha_values_2[l - 1])
                a_2 = max(0, i + 1 - self.T_1) * alpha_values_3[l - 1]
                eps_t = (a_1 + a_2 - np.log(self.delta)) / float(l)
                epsilon_values.append(eps_t)
            idx = np.argmin(epsilon_values)
            cnt[idx] += 1
            eps_values.append(epsilon_values[idx])
        print(cnt)
        print('Epsilon after clustering:', eps_values[self.T_1])
        print('Epsilon after GD:', eps_values[self.T - 1])
        eps_values = eps_values[self.T_1:]
        eps_values = np.array(eps_values)[np.arange(int(1.0 / self.q) - 1,
            self.T_2, int(1.0 / self.q))]
        for i in range(self.epoch_num):
            print('Epoch:', i + 1, 'Eps:', eps_values[i])
        return eps_values

<mask token>
matplotlib.use('Agg')
<mask token>


class calculate_accountant(object):

    def __init__(self, q, epoch_num, sigma_cl_1, sigma_cl_2, sigma_gd,
        sigma_dyn_clip, delta):
        self.l_max = 32
        self.q = q
        self.epoch_num = epoch_num
        self.sigma_cl_1 = sigma_cl_1
        self.sigma_cl_2 = sigma_cl_2
        self.sigma_gd = sigma_gd
        self.sigma_dyn_clip = sigma_dyn_clip
        self.delta = delta
        self.T_1 = 10
        self.T_2 = epoch_num * int(1.0 / q)
        self.T = self.T_1 + self.T_2

    def compute_gauss(self, p_sigma, l, u, d=20, p_sens=1.0):
        e_1 = lambda x: norm.pdf(x, 0, p_sens * p_sigma) * (norm.pdf(x, 0, 
            p_sens * p_sigma) / ((1 - self.q) * norm.pdf(x, 0, p_sens *
            p_sigma) + self.q * norm.pdf(x, p_sens, p_sens * p_sigma))) ** (u *
            l)
        e_2 = lambda x: ((1 - self.q) * norm.pdf(x, 0, p_sens * p_sigma) + 
            self.q * norm.pdf(x, p_sens, p_sens * p_sigma)) * (((1 - self.q
            ) * norm.pdf(x, 0, p_sens * p_sigma) + self.q * norm.pdf(x,
            p_sens, p_sens * p_sigma)) / norm.pdf(x, 0, p_sens * p_sigma)) ** (
            u * l)
        E_1, _ = integrate.quad(e_1, -d, d)
        E_1 = np.nan_to_num(E_1)
        E_2, _ = integrate.quad(e_2, -d, d)
        E_2 = np.nan_to_num(E_2)
        return np.log(np.maximum(np.abs(E_1), np.abs(E_2)) + 1e-08)

    def compute_alpha(self, p_sigma, q, d):
        alpha_values = []
        for l in range(1, self.l_max + 1):
            e_1 = lambda x: norm.pdf(x, 0, p_sigma) * (norm.pdf(x, 0,
                p_sigma) / ((1 - q) * norm.pdf(x, 0, p_sigma) + q * norm.
                pdf(x, 1.0, p_sigma))) ** l
            e_2 = lambda x: ((1 - q) * norm.pdf(x, 0, p_sigma) + q * norm.
                pdf(x, 1.0, p_sigma)) * (((1 - q) * norm.pdf(x, 0, p_sigma) +
                q * norm.pdf(x, 1.0, p_sigma)) / norm.pdf(x, 0, p_sigma)) ** l
            E_1, _ = integrate.quad(e_1, -d, d)
            E_2, _ = integrate.quad(e_2, -d, d)
            E_1 = np.nan_to_num(E_1)
            E_2 = np.nan_to_num(E_2)
            val = np.log(max([abs(E_1), abs(E_2)]) + 1e-11)
            alpha_values.append(val)
        return alpha_values

    def get_privacy(self):
        alpha_values_1 = self.compute_alpha(self.sigma_cl_1, 1.0, 100)
        alpha_values_2 = self.compute_alpha(self.sigma_cl_2, 1.0, 500)
        alpha_values_3 = []
        for l in range(1, self.l_max + 1):
            tmp = []
            for j in range(1, 10):
                u1 = j / 10.0
                u2 = (10 - j) / 10.0
                alpha_value1 = self.compute_gauss(self.sigma_dyn_clip, l, 
                    1.0 / u1, d=50, p_sens=2.0 ** 0.5)
                alpha_value2 = self.compute_gauss(self.sigma_gd, l, 1.0 /
                    u2, d=20)
                tmp.append(u1 * alpha_value1 + u2 * alpha_value2)
            val = min(tmp)
            alpha_values_3.append(val)
        eps_values = []
        cnt = Counter()
        for i in range(self.T):
            epsilon_values = []
            for l in range(1, self.l_max + 1):
                a_1 = min(i + 1, self.T_1) * (alpha_values_1[l - 1] +
                    alpha_values_2[l - 1])
                a_2 = max(0, i + 1 - self.T_1) * alpha_values_3[l - 1]
                eps_t = (a_1 + a_2 - np.log(self.delta)) / float(l)
                epsilon_values.append(eps_t)
            idx = np.argmin(epsilon_values)
            cnt[idx] += 1
            eps_values.append(epsilon_values[idx])
        print(cnt)
        print('Epsilon after clustering:', eps_values[self.T_1])
        print('Epsilon after GD:', eps_values[self.T - 1])
        eps_values = eps_values[self.T_1:]
        eps_values = np.array(eps_values)[np.arange(int(1.0 / self.q) - 1,
            self.T_2, int(1.0 / self.q))]
        for i in range(self.epoch_num):
            print('Epoch:', i + 1, 'Eps:', eps_values[i])
        return eps_values

import numpy as np
from scipy.stats import norm, laplace
import _pickle as pickle
from math import exp
import matplotlib
from collections import Counter
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import scipy.integrate as integrate


class calculate_accountant(object):

    def __init__(self, q, epoch_num, sigma_cl_1, sigma_cl_2, sigma_gd,
        sigma_dyn_clip, delta):
        self.l_max = 32
        self.q = q
        self.epoch_num = epoch_num
        self.sigma_cl_1 = sigma_cl_1
        self.sigma_cl_2 = sigma_cl_2
        self.sigma_gd = sigma_gd
        self.sigma_dyn_clip = sigma_dyn_clip
        self.delta = delta
        self.T_1 = 10
        self.T_2 = epoch_num * int(1.0 / q)
        self.T = self.T_1 + self.T_2

    def compute_gauss(self, p_sigma, l, u, d=20, p_sens=1.0):
        e_1 = lambda x: norm.pdf(x, 0, p_sens * p_sigma) * (norm.pdf(x, 0, 
            p_sens * p_sigma) / ((1 - self.q) * norm.pdf(x, 0, p_sens *
            p_sigma) + self.q * norm.pdf(x, p_sens, p_sens * p_sigma))) ** (u *
            l)
        e_2 = lambda x: ((1 - self.q) * norm.pdf(x, 0, p_sens * p_sigma) + 
            self.q * norm.pdf(x, p_sens, p_sens * p_sigma)) * (((1 - self.q
            ) * norm.pdf(x, 0, p_sens * p_sigma) + self.q * norm.pdf(x,
            p_sens, p_sens * p_sigma)) / norm.pdf(x, 0, p_sens * p_sigma)) ** (
            u * l)
        E_1, _ = integrate.quad(e_1, -d, d)
        E_1 = np.nan_to_num(E_1)
        E_2, _ = integrate.quad(e_2, -d, d)
        E_2 = np.nan_to_num(E_2)
        return np.log(np.maximum(np.abs(E_1), np.abs(E_2)) + 1e-08)

    def compute_alpha(self, p_sigma, q, d):
        alpha_values = []
        for l in range(1, self.l_max + 1):
            e_1 = lambda x: norm.pdf(x, 0, p_sigma) * (norm.pdf(x, 0,
                p_sigma) / ((1 - q) * norm.pdf(x, 0, p_sigma) + q * norm.
                pdf(x, 1.0, p_sigma))) ** l
            e_2 = lambda x: ((1 - q) * norm.pdf(x, 0, p_sigma) + q * norm.
                pdf(x, 1.0, p_sigma)) * (((1 - q) * norm.pdf(x, 0, p_sigma) +
                q * norm.pdf(x, 1.0, p_sigma)) / norm.pdf(x, 0, p_sigma)) ** l
            E_1, _ = integrate.quad(e_1, -d, d)
            E_2, _ = integrate.quad(e_2, -d, d)
            E_1 = np.nan_to_num(E_1)
            E_2 = np.nan_to_num(E_2)
            val = np.log(max([abs(E_1), abs(E_2)]) + 1e-11)
            alpha_values.append(val)
        return alpha_values

    def get_privacy(self):
        alpha_values_1 = self.compute_alpha(self.sigma_cl_1, 1.0, 100)
        alpha_values_2 = self.compute_alpha(self.sigma_cl_2, 1.0, 500)
        alpha_values_3 = []
        for l in range(1, self.l_max + 1):
            tmp = []
            for j in range(1, 10):
                u1 = j / 10.0
                u2 = (10 - j) / 10.0
                alpha_value1 = self.compute_gauss(self.sigma_dyn_clip, l, 
                    1.0 / u1, d=50, p_sens=2.0 ** 0.5)
                alpha_value2 = self.compute_gauss(self.sigma_gd, l, 1.0 /
                    u2, d=20)
                tmp.append(u1 * alpha_value1 + u2 * alpha_value2)
            val = min(tmp)
            alpha_values_3.append(val)
        eps_values = []
        cnt = Counter()
        for i in range(self.T):
            epsilon_values = []
            for l in range(1, self.l_max + 1):
                a_1 = min(i + 1, self.T_1) * (alpha_values_1[l - 1] +
                    alpha_values_2[l - 1])
                a_2 = max(0, i + 1 - self.T_1) * alpha_values_3[l - 1]
                eps_t = (a_1 + a_2 - np.log(self.delta)) / float(l)
                epsilon_values.append(eps_t)
            idx = np.argmin(epsilon_values)
            cnt[idx] += 1
            eps_values.append(epsilon_values[idx])
        print(cnt)
        print('Epsilon after clustering:', eps_values[self.T_1])
        print('Epsilon after GD:', eps_values[self.T - 1])
        eps_values = eps_values[self.T_1:]
        eps_values = np.array(eps_values)[np.arange(int(1.0 / self.q) - 1,
            self.T_2, int(1.0 / self.q))]
        for i in range(self.epoch_num):
            print('Epoch:', i + 1, 'Eps:', eps_values[i])
        return eps_values

import numpy as np
from scipy.stats import norm, laplace
import _pickle as pickle
from math import exp
import matplotlib
from collections import Counter
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import scipy.integrate as integrate

class calculate_accountant(object):

    def __init__(self, q, epoch_num, sigma_cl_1, sigma_cl_2, sigma_gd, sigma_dyn_clip, delta):
        self.l_max = 32

        # Sampling prob. for SGD
        self.q = q
        # epoch num for SGD
        self.epoch_num = epoch_num

        # noise for k-means cluster size
        self.sigma_cl_1 = sigma_cl_1
        # noise for k-means cluster centers
        self.sigma_cl_2 = sigma_cl_2
        # noise for Gradient descent
        self.sigma_gd = sigma_gd
        # noise for Adaptive clipping in each GD step
        self.sigma_dyn_clip = sigma_dyn_clip

        self.delta = delta
        # k-means iterations
        self.T_1 = 10
        # SGD iterations
        self.T_2 = epoch_num * int(1.0/q)
        # Total iterations
        self.T = self.T_1 + self.T_2

        # mu_0 = sigma* N(0,1)
        # mu_1 = sigma* N(0,1) + 1
        # mu = (1-q)* mu_0 + q*mu_1 = sigma*N(0,1) + q
        # Probability is not incorporated when the the expected values are computed



    def compute_gauss(self,p_sigma, l, u, d=20, p_sens=1.0):
        e_1 = lambda x: norm.pdf(x,0,p_sens*p_sigma)* (norm.pdf(x,0,p_sens*p_sigma) / ((1-self.q)*norm.pdf(x,0,p_sens*p_sigma) +
            self.q*norm.pdf(x,p_sens,p_sens*p_sigma)))**(u*l)
        e_2 = lambda x: ((1-self.q)*norm.pdf(x,0,p_sens*p_sigma) + self.q*norm.pdf(x,p_sens,p_sens*p_sigma))*\
                (((1-self.q)*norm.pdf(x,0,p_sens*p_sigma) + self.q*norm.pdf(x,p_sens,p_sens*p_sigma)) /
                        norm.pdf(x,0,p_sens*p_sigma))**(u*l)
        E_1, _ = integrate.quad(e_1,-d,d)
        E_1 = np.nan_to_num(E_1)
        E_2, _ = integrate.quad(e_2,-d,d)
        E_2 = np.nan_to_num(E_2)
        return np.log(np.maximum(np.abs(E_1), np.abs(E_2))+0.00000001)

    def compute_alpha(self, p_sigma, q, d):
        #print("Sigma:", p_sigma)
        #print("Precomputing...")
        alpha_values = []
        for l in range(1,self.l_max+1):
            e_1 = lambda x: norm.pdf(x,0,p_sigma)* (norm.pdf(x,0,p_sigma) / ((1-q)*norm.pdf(x,0,p_sigma) +
                q*norm.pdf(x,1.0,p_sigma)))**l
            e_2 = lambda x: ((1-q)*norm.pdf(x,0,p_sigma) + q*norm.pdf(x,1.0,p_sigma))*\
                (((1-q)*norm.pdf(x,0,p_sigma) + q*norm.pdf(x,1.0,p_sigma))/(norm.pdf(x,0,p_sigma)))**l

            E_1, _ = integrate.quad(e_1,-d,d)
            E_2, _ = integrate.quad(e_2,-d,d)
            E_1 = np.nan_to_num(E_1)
            E_2 = np.nan_to_num(E_2)
            val = np.log(max([abs(E_1), abs(E_2)])+0.00000000001)
            #print("l = ", l, " alpha = ", val)
            alpha_values.append(val)
        #print("done.")
        return alpha_values

    def get_privacy(self):
        ## Alpha computation


        # Clustering
        alpha_values_1 = self.compute_alpha(self.sigma_cl_1, 1.0, 100)
        alpha_values_2 = self.compute_alpha(self.sigma_cl_2, 1.0, 500)


        # DP SGD
        alpha_values_3 = []

        for l in range(1, self.l_max+1):
            tmp = []
            for j in range(1, 10):
                u1 = j/10.0
                u2 = (10-j) / 10.0
                alpha_value1 = self.compute_gauss(self.sigma_dyn_clip, l, 1.0/u1, d=50, p_sens=2.0**0.5)
                alpha_value2 = self.compute_gauss(self.sigma_gd, l, 1.0/u2, d=20)
                tmp.append(u1 * alpha_value1 + u2 * alpha_value2)
            val = min(tmp)
            #print("l = ", l, " alpha = ", val)
            alpha_values_3.append(val)

        ## Epsilon computation

        eps_values = [ ]
        cnt = Counter()
        for i in range(self.T):
            epsilon_values = []
            for l in range(1, self.l_max+1):
                # Alpha of Clustering
                a_1 = min(i+1, self.T_1)*(alpha_values_1[l-1] + alpha_values_2[l-1])
                # Alpha of GD
                a_2 = max(0, i+1 - self.T_1)*alpha_values_3[l-1]

                eps_t = (a_1 + a_2 - np.log(self.delta)) / float(l)
                epsilon_values.append(eps_t)

            idx = np.argmin(epsilon_values)
            cnt[idx]+=1
            eps_values.append(epsilon_values[idx])


        print(cnt)

        print("Epsilon after clustering:", eps_values[self.T_1])
        print("Epsilon after GD:", eps_values[self.T-1])

        eps_values = eps_values[self.T_1:]
        eps_values = np.array(eps_values)[np.arange(int(1.0/self.q)-1,self.T_2,int(1.0/self.q))]

        for i in range(self.epoch_num):
            print("Epoch:", i+1, "Eps:", eps_values[i])
        return eps_values

