From 5748362c372a133d9aa05ddd4fd78ff59e5e3415 Mon Sep 17 00:00:00 2001 From: stefan-aws Date: Thu, 14 Mar 2024 00:05:00 +0000 Subject: [PATCH 01/12] skeleton --- build/py/run_benchmarks.sh | 3 + docs/py/Benchmarks.py | 59 +++++++ docs/py/IBM.py | 306 +++++++++++++++++++++++++++++++++++++ 3 files changed, 368 insertions(+) create mode 100644 build/py/run_benchmarks.sh create mode 100644 docs/py/Benchmarks.py create mode 100644 docs/py/IBM.py diff --git a/build/py/run_benchmarks.sh b/build/py/run_benchmarks.sh new file mode 100644 index 00000000..0cadb2d9 --- /dev/null +++ b/build/py/run_benchmarks.sh @@ -0,0 +1,3 @@ +#1/bin/bash + +PYTHONPATH=.:build/py/DafnyVMC-py:docs/py/IBM.py python3 docs/py/Benchmarks.py \ No newline at end of file diff --git a/docs/py/Benchmarks.py b/docs/py/Benchmarks.py new file mode 100644 index 00000000..60260eba --- /dev/null +++ b/docs/py/Benchmarks.py @@ -0,0 +1,59 @@ +import timeit +import secrets +import numpy +import matplotlib.pyplot as plt +from decimal import Decimal +import DafnyVMC +from IBM import sample_dgauss + +sigma_range = 1.0 +step = 0.001 +sigmas = numpy.arange(0.000001, sigma_range, step).tolist() + +vmc_mean = [] +vmc_std = [] +ibm_mean = [] +ibm_std = [] + +fig,ax1 = plt.subplots() + +r = DafnyVMC.Random() + +for sigma in sigmas: + vmc = [] + ibm = [] + for i in range(1100): + num, denom = Decimal(sigma).as_integer_ratio() + start_time = timeit.default_timer() + r.DiscreteGaussianSample(num, denom) + elapsed = timeit.default_timer() - start_time + vmc.append(elapsed) + + for i in range(1100): + start_time = timeit.default_timer() + sample_dgauss(sigma ** 2, rng=secrets.SystemRandom()) + elapsed = timeit.default_timer() - start_time + ibm.append(elapsed) + + vmc = numpy.array(vmc[-1000:]) + ibm = numpy.array(ibm[-1000:]) + + vmc_mean.append(vmc.mean()*1000.0) + vmc_std.append(vmc.std()*1000.0) + ibm_mean.append(ibm.mean()*1000.0) + ibm_std.append(ibm.std()*1000.0) + +ax1.plot(sigmas, vmc_mean, color='green', linewidth=1.0, label='VMC') +ax1.fill_between(sigmas, numpy.array(vmc_mean)-0.5*numpy.array(vmc_std), numpy.array(vmc_mean)+0.5*numpy.array(vmc_std), + alpha=0.2, facecolor='k', + linewidth=2, linestyle='dashdot', antialiased=True) + +ax1.plot(sigmas, ibm_mean, color='red', linewidth=1.0, label='IBM') +ax1.fill_between(sigmas, numpy.array(ibm_mean)-0.5*numpy.array(ibm_std), numpy.array(ibm_mean)+0.5*numpy.array(ibm_std), + alpha=0.2, facecolor='y', + linewidth=2, linestyle='dashdot', antialiased=True) + +ax1.set_xlabel("Sigma") +ax1.set_ylabel("Sampling Time (ms)") +plt.legend(loc = 'best') +plt.savefig('Benchmarks.pdf') \ No newline at end of file diff --git a/docs/py/IBM.py b/docs/py/IBM.py new file mode 100644 index 00000000..c141b42e --- /dev/null +++ b/docs/py/IBM.py @@ -0,0 +1,306 @@ +# Module: IBM + +# Implementation of exact discrete gaussian distribution sampler +# See https://arxiv.org/abs/2004.00010 +# - Thomas Steinke dgauss@thomas-steinke.net 2020 + +import random #Default random number generator, +#random.SecureRandom() provides high-quality randomness from /dev/urandom or similar +from fractions import Fraction #we will work with rational numbers + +#sample uniformly from range(m) +#all randomness comes from calling this +def sample_uniform(m,rng): + assert isinstance(m,int) #python 3 + #assert isinstance(m,(int,long)) #python 2 + assert m>0 + return rng.randrange(m) + +#sample from a Bernoulli(p) distribution +#assumes p is a rational number in [0,1] +def sample_bernoulli(p,rng): + assert isinstance(p,Fraction) + assert 0 <= p <= 1 + m=sample_uniform(p.denominator,rng) + if m < p.numerator: + return 1 + else: + return 0 + +#sample from a Bernoulli(exp(-x)) distribution +#assumes x is a rational number in [0,1] +def sample_bernoulli_exp1(x,rng): + assert isinstance(x,Fraction) + assert 0 <= x <= 1 + k=1 + while True: + if sample_bernoulli(x/k,rng)==1: + k=k+1 + else: + break + return k%2 + +#sample from a Bernoulli(exp(-x)) distribution +#assumes x is a rational number >=0 +def sample_bernoulli_exp(x,rng): + assert isinstance(x,Fraction) + assert x >= 0 + #Sample floor(x) independent Bernoulli(exp(-1)) + #If all are 1, return Bernoulli(exp(-(x-floor(x)))) + while x>1: + if sample_bernoulli_exp1(Fraction(1,1),rng)==1: + x=x-1 + else: + return 0 + return sample_bernoulli_exp1(x,rng) + +#sample from a geometric(1-exp(-x)) distribution +#assumes x is a rational number >= 0 +def sample_geometric_exp_slow(x,rng): + assert isinstance(x,Fraction) + assert x >= 0 + k=0 + while True: + if sample_bernoulli_exp(x,rng)==1: + k=k+1 + else: + return k + +#sample from a geometric(1-exp(-x)) distribution +#assumes x >= 0 rational +def sample_geometric_exp_fast(x,rng): + assert isinstance(x,Fraction) + if x==0: return 0 #degenerate case + assert x>0 + + t=x.denominator + while True: + u=sample_uniform(t,rng) + b=sample_bernoulli_exp(Fraction(u,t),rng) + if b==1: + break + v=sample_geometric_exp_slow(Fraction(1,1),rng) + value = v*t+u + return value//x.numerator + +#sample from a discrete Laplace(scale) distribution +#Returns integer x with Pr[x] = exp(-abs(x)/scale)*(exp(1/scale)-1)/(exp(1/scale)+1) +#casts scale to Fraction +#assumes scale>=0 +def sample_dlaplace(scale,rng=None): + if rng is None: + rng = random.SystemRandom() + scale = Fraction(scale) + assert scale >= 0 + while True: + sign=sample_bernoulli(Fraction(1,2),rng) + magnitude=sample_geometric_exp_fast(1/scale,rng) + if sign==1 and magnitude==0: continue + return magnitude*(1-2*sign) + +#compute floor(sqrt(x)) exactly +#only requires comparisons between x and integer +def floorsqrt(x): + assert x >= 0 + #a,b integers + a=0 #maintain a^2<=x + b=1 #maintain b^2>x + while b*b <= x: + b=2*b #double to get upper bound + #now do binary search + while a+1=0 +def sample_dgauss(sigma2,rng=None): + if rng is None: + rng = random.SystemRandom() + sigma2=Fraction(sigma2) + if sigma2==0: return 0 #degenerate case + assert sigma2 > 0 + t = floorsqrt(sigma2)+1 + while True: + candidate = sample_dlaplace(t,rng=rng) + bias=((abs(candidate)-sigma2/t)**2)/(2*sigma2) + if sample_bernoulli_exp(bias,rng)==1: + return candidate + +######################################################################### +#DONE That's it! Now some utilities + +import math #need this, code below is no longer exact + +#Compute the normalizing constant of the discrete gaussian +#i.e. sum_{x in Z} exp(-x^2/2sigma2) +#By Poisson summation formula, this is equivalent to +# sqrt{2*pi*sigma2}*sum_{y in Z} exp(-2*pi^2*sigma2*y^2) +#For small sigma2 the former converges faster +#For large sigma2, the latter converges faster +#crossover at sigma2=1/2*pi +#For intermediate sigma2, this code will compute both and check +def normalizing_constant(sigma2): + original=None + poisson=None + if sigma2<=1: + original = 0 + x=1000 #summation stops at exp(-x^2/2sigma2)<=exp(-500,000) + while x>0: + original = original + math.exp(-x*x/(2.0*sigma2)) + x = x - 1 #sum from small to large for improved accuracy + original = 2*original + 1 #symmetrize and add x=0 + if sigma2*100 >= 1: + poisson = 0 + y = 1000 #summation stops at exp(-y^2*2*pi^2*sigma2)<=exp(-190,000) + while y>0: + poisson = poisson + math.exp(-math.pi*math.pi*sigma2*2*y*y) + y = y - 1 #sum from small to large + poisson = math.sqrt(2*math.pi*sigma2)*(1+2*poisson) + if poisson is None: return original + if original is None: return poisson + #if we have computed both, check equality + scale = max(1,math.sqrt(2*math.pi*sigma2)) #tight-ish lower bound on constant + assert -1e-15*scale <= original-poisson <= 1e-15*scale + #10^-15 is about as much precision as we can expect from double precision floating point numbers + #64-bit float has 56-bit mantissa 10^-15 ~= 2^-50 + return (original+poisson)/2 + +#compute the variance of discrete gaussian +#mean is zero, thus: +#var = sum_{x in Z} x^2*exp(-x^2/(2*sigma2)) / normalizing_constant(sigma2) +#By Poisson summation formula, we have equivalent expression: +# variance(sigma2) = sigma2 * (1 - 4*pi^2*sigma2*variance(1/(4*pi^2*sigma2)) ) +#See lemma 20 https://arxiv.org/pdf/2004.00010v3.pdf#page=17 +#alternative expression converges faster when sigma2 is large +#crossover point (in terms of convergence) is sigma2=1/(2*pi) +#for intermediate values of sigma2, we compute both expressions and check +def variance(sigma2): + original=None + poisson=None + if sigma2<=1: #compute primary expression + original=0 + x = 1000 #summation stops at exp(-x^2/2sigma2)<=exp(-500,000) + while x>0: #sum from small to large for improved accuracy + original = original + x*x*math.exp(-x*x/(2.0*sigma2)) + x=x-1 + original = 2*original/normalizing_constant(sigma2) + if sigma2*100>=1: + poisson=0 #we will compute sum_{y in Z} y^2 * exp(-2*pi^2*sigma2*y^2) + y=1000 #summation stops at exp(-y^2*2*pi^2*sigma2)<=exp(-190,000) + while y>0: #sum from small to large + poisson = poisson + y*y*math.exp(-y*y*2*sigma2*math.pi*math.pi) + y=y-1 + poisson = 2*poisson/normalizing_constant(1/(4*sigma2*math.pi*math.pi)) + #next convert from variance(1/(4*pi^2*sigma2)) to variance(sigma2) + poisson = sigma2*(1-4*sigma2*poisson*math.pi*math.pi) + if original is None: return poisson + if poisson is None: return original + #if we have computed both check equality + assert -1e-15*sigma2 <= original-poisson <= 1e-15*sigma2 + return (original+poisson)/2 + +######################################################################### +#DONE Now some basic testing code + +import matplotlib.pyplot as plt #only needed for testing +import time #only needed for testing + +#This generates n samples from sample_dgauss(sigma2) +#It times this and releases statistics +#produces a histogram plot if plot==True +#if plot==None it will only produce a histogram if it's not too large +#can save image instead of displaying by specifying a path e.g., save="plot.png" +def plot_histogram(sigma2,n,save=None,plot=None): + #generate samples + before=time.time() + samples = [sample_dgauss(sigma2) for i in range(n)] + after=time.time() + print("generated "+str(n)+" samples in "+str(after-before)+" seconds ("+str(n/(after-before))+" samples per second) for sigma^2="+str(sigma2)) + #now process + samples.sort() + values=[] + counts=[] + counter=None + prev=None + for sample in samples: + if prev is None: #initializing + prev=sample + counter=1 + elif sample==prev: #still same element + counter=counter+1 + else: + #add prev to histogram + values.append(prev) + counts.append(counter) + #start counting + prev=sample + counter=1 + #add final value + values.append(prev) + counts.append(counter) + + #print & sum + sum=0 + sumsquared=0 + kl=0 #compute KL divergence betwen empirical distribution and true distribution + norm_const=normalizing_constant(sigma2) + true_var=variance(sigma2) + for i in range(len(values)): + if len(values)<=100: #don't print too much + print(str(values[i])+":\t"+str(counts[i])) + sum = sum + values[i]*counts[i] + sumsquared = sumsquared + values[i]*values[i]*counts[i] + kl = kl + counts[i]*(math.log(counts[i]*norm_const/n)+values[i]*values[i]/(2.0*sigma2)) + mean = Fraction(sum,n) + var=Fraction(sumsquared,n) + kl=kl/n + print("mean="+str(float(mean))+" (true=0)") + print("variance="+str(float(var))+" (true="+str(true_var)+")") + print("KL(empirical||true)="+str(kl)) # https://en.wikipedia.org/wiki/G-test + assert kl>0 #kl divergence always >=0 and ==0 iff empirical==true, which is impossible + #now plot + if plot is None: + plot = (len(values)<=1000) #don't plot if huge + if not plot: return + ideal_counts = [n*math.exp(-x*x/(2.0*sigma2))/norm_const for x in values] + plt.bar(values, counts) + plt.plot(values, ideal_counts,'r') + plt.title("Histogram of samples from discrete Gaussian\nsigma^2="+str(sigma2)+" n="+str(n)) + if save is None: + plt.show() + else: + plt.savefig(save) + plt.clf() + +if __name__ == '__main__': + print("This is the discrete Gaussian sampler") + print("See the paper https://arxiv.org/abs/2004.00010") + print("Now running some basic testing code") + print("Start by calculating normalizing constant and variance for different values") + #some test code for normalizing_constant and variance functions + for sigma2 in [0.1**100,0.1**6,0.001,0.01,0.03,0.05,0.08,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,10,100,10**6,10**20,10**100]: + #internal asserts do some testing when 0.01<=sigma2<=1 + c=normalizing_constant(sigma2) + v=variance(sigma2) + #print + print("sigma^2="+str(sigma2) + ":\tnorm_const=" + str(c) + "=sqrt{2*pi}*sigma*" + str(c/math.sqrt(2*math.pi*sigma2)) + "\tvar=" + str(v)) + #print a few samples + #for i in range(20): print sample_dgauss(1) + #plot histogram and statistics + #includes timing + print("Now run the sampler") + print("Start with very large sigma^2=10^100 -- for timing purposes only") + plot_histogram(10**100,100000,plot=False) #large var, this will just be for timing + print("Now sigma^2=10 -- will display a histogram") + plot_histogram(10,100000) #small var, this will produce plot From c9a51dcc6fd0791dc73cb11bbf1250f2e75e426d Mon Sep 17 00:00:00 2001 From: stefan-aws Date: Thu, 14 Mar 2024 17:55:29 +0000 Subject: [PATCH 02/12] fix --- docs/py/Benchmarks.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/docs/py/Benchmarks.py b/docs/py/Benchmarks.py index 60260eba..abc67b6b 100644 --- a/docs/py/Benchmarks.py +++ b/docs/py/Benchmarks.py @@ -5,6 +5,7 @@ from decimal import Decimal import DafnyVMC from IBM import sample_dgauss +from datetime import datetime sigma_range = 1.0 step = 0.001 @@ -22,16 +23,18 @@ for sigma in sigmas: vmc = [] ibm = [] + sigma_num, sigma_denom = Decimal(sigma).as_integer_ratio() + sigma_squared = sigma ** 2 + for i in range(1100): - num, denom = Decimal(sigma).as_integer_ratio() start_time = timeit.default_timer() - r.DiscreteGaussianSample(num, denom) + r.DiscreteGaussianSample(sigma_num, sigma_denom) elapsed = timeit.default_timer() - start_time vmc.append(elapsed) for i in range(1100): start_time = timeit.default_timer() - sample_dgauss(sigma ** 2, rng=secrets.SystemRandom()) + sample_dgauss(sigma_squared, rng=secrets.SystemRandom()) elapsed = timeit.default_timer() - start_time ibm.append(elapsed) @@ -56,4 +59,6 @@ ax1.set_xlabel("Sigma") ax1.set_ylabel("Sampling Time (ms)") plt.legend(loc = 'best') -plt.savefig('Benchmarks.pdf') \ No newline at end of file +now = datetime.now() +filename = 'Benchmarks' + now.strftime("%H_%M_%S") + '.pdf' +plt.savefig(filename) \ No newline at end of file From d481afabf757bb05036163216368af65f4b0d5af Mon Sep 17 00:00:00 2001 From: stefan-aws Date: Thu, 14 Mar 2024 17:56:18 +0000 Subject: [PATCH 03/12] range --- docs/py/Benchmarks.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/py/Benchmarks.py b/docs/py/Benchmarks.py index abc67b6b..a0685035 100644 --- a/docs/py/Benchmarks.py +++ b/docs/py/Benchmarks.py @@ -7,8 +7,8 @@ from IBM import sample_dgauss from datetime import datetime -sigma_range = 1.0 -step = 0.001 +sigma_range = 100 +step = 1 sigmas = numpy.arange(0.000001, sigma_range, step).tolist() vmc_mean = [] From d5cbb4ddbba7f1bdf036b8dc02312c23cf0118d4 Mon Sep 17 00:00:00 2001 From: stefan-aws Date: Fri, 15 Mar 2024 12:03:32 +0000 Subject: [PATCH 04/12] submodules --- .gitmodules | 6 + .../benchmarks.py} | 40 ++- .../Benchmarks/differential-privacy-library | 1 + .../discrete-gaussian-differential-privacy | 1 + docs/py/IBM.py | 306 ------------------ 5 files changed, 40 insertions(+), 314 deletions(-) create mode 100644 .gitmodules rename docs/py/{Benchmarks.py => Benchmarks/benchmarks.py} (62%) create mode 160000 docs/py/Benchmarks/differential-privacy-library create mode 160000 docs/py/Benchmarks/discrete-gaussian-differential-privacy delete mode 100644 docs/py/IBM.py diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 00000000..12b39b25 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,6 @@ +[submodule "docs/py/Benchmarks/discrete-gaussian-differential-privacy"] + path = docs/py/Benchmarks/discrete-gaussian-differential-privacy + url = https://github.com/IBM/discrete-gaussian-differential-privacy/ +[submodule "docs/py/Benchmarks/differential-privacy-library"] + path = docs/py/Benchmarks/differential-privacy-library + url = https://github.com/IBM/differential-privacy-library/ diff --git a/docs/py/Benchmarks.py b/docs/py/Benchmarks/benchmarks.py similarity index 62% rename from docs/py/Benchmarks.py rename to docs/py/Benchmarks/benchmarks.py index a0685035..074196d2 100644 --- a/docs/py/Benchmarks.py +++ b/docs/py/Benchmarks/benchmarks.py @@ -4,25 +4,33 @@ import matplotlib.pyplot as plt from decimal import Decimal import DafnyVMC -from IBM import sample_dgauss +from diffprivlib.mechanisms import GaussianDiscrete +import discretegauss from datetime import datetime - -sigma_range = 100 -step = 1 -sigmas = numpy.arange(0.000001, sigma_range, step).tolist() +import tqdm vmc_mean = [] vmc_std = [] ibm_mean = [] ibm_std = [] +ibm2_mean = [] +ibm2_std = [] fig,ax1 = plt.subplots() +rng = secrets.SystemRandom() r = DafnyVMC.Random() -for sigma in sigmas: +sigmas = [] +for epsilon_times_100 in tqdm.tqdm(range(1, 500, 2)): vmc = [] ibm = [] + ibm2= [] + + g = GaussianDiscrete(epsilon=0.01 * epsilon_times_100, delta=0.00001) + sigma = g._scale + sigmas += [sigma] + sigma_num, sigma_denom = Decimal(sigma).as_integer_ratio() sigma_squared = sigma ** 2 @@ -34,28 +42,44 @@ for i in range(1100): start_time = timeit.default_timer() - sample_dgauss(sigma_squared, rng=secrets.SystemRandom()) + discretegauss.sample_dgauss(sigma_squared, rng) elapsed = timeit.default_timer() - start_time ibm.append(elapsed) + for i in range(1100): + start_time = timeit.default_timer() + g.randomise(0) + elapsed = timeit.default_timer() - start_time + ibm2.append(elapsed) + vmc = numpy.array(vmc[-1000:]) ibm = numpy.array(ibm[-1000:]) + ibm2 = numpy.array(ibm2[-1000:]) vmc_mean.append(vmc.mean()*1000.0) vmc_std.append(vmc.std()*1000.0) ibm_mean.append(ibm.mean()*1000.0) ibm_std.append(ibm.std()*1000.0) + ibm2_mean.append(ibm2.mean()*1000.0) + ibm2_std.append(ibm2.std()*1000.0) + +print(sigmas) ax1.plot(sigmas, vmc_mean, color='green', linewidth=1.0, label='VMC') ax1.fill_between(sigmas, numpy.array(vmc_mean)-0.5*numpy.array(vmc_std), numpy.array(vmc_mean)+0.5*numpy.array(vmc_std), alpha=0.2, facecolor='k', linewidth=2, linestyle='dashdot', antialiased=True) -ax1.plot(sigmas, ibm_mean, color='red', linewidth=1.0, label='IBM') +ax1.plot(sigmas, ibm_mean, color='red', linewidth=1.0, label='IBM-DPL') ax1.fill_between(sigmas, numpy.array(ibm_mean)-0.5*numpy.array(ibm_std), numpy.array(ibm_mean)+0.5*numpy.array(ibm_std), alpha=0.2, facecolor='y', linewidth=2, linestyle='dashdot', antialiased=True) +ax1.plot(sigmas, ibm2_mean, color='purple', linewidth=1.0, label='IBM-DGDP') +ax1.fill_between(sigmas, numpy.array(ibm2_mean)-0.5*numpy.array(ibm2_std), numpy.array(ibm2_mean)+0.5*numpy.array(ibm2_std), + alpha=0.2, facecolor='y', + linewidth=2, linestyle='dashdot', antialiased=True) + ax1.set_xlabel("Sigma") ax1.set_ylabel("Sampling Time (ms)") plt.legend(loc = 'best') diff --git a/docs/py/Benchmarks/differential-privacy-library b/docs/py/Benchmarks/differential-privacy-library new file mode 160000 index 00000000..2ec58659 --- /dev/null +++ b/docs/py/Benchmarks/differential-privacy-library @@ -0,0 +1 @@ +Subproject commit 2ec58659f8a1795a02f73d1fb0506477b3345a55 diff --git a/docs/py/Benchmarks/discrete-gaussian-differential-privacy b/docs/py/Benchmarks/discrete-gaussian-differential-privacy new file mode 160000 index 00000000..cb190d2a --- /dev/null +++ b/docs/py/Benchmarks/discrete-gaussian-differential-privacy @@ -0,0 +1 @@ +Subproject commit cb190d2a990a78eff6e21159203bc888e095f01b diff --git a/docs/py/IBM.py b/docs/py/IBM.py deleted file mode 100644 index c141b42e..00000000 --- a/docs/py/IBM.py +++ /dev/null @@ -1,306 +0,0 @@ -# Module: IBM - -# Implementation of exact discrete gaussian distribution sampler -# See https://arxiv.org/abs/2004.00010 -# - Thomas Steinke dgauss@thomas-steinke.net 2020 - -import random #Default random number generator, -#random.SecureRandom() provides high-quality randomness from /dev/urandom or similar -from fractions import Fraction #we will work with rational numbers - -#sample uniformly from range(m) -#all randomness comes from calling this -def sample_uniform(m,rng): - assert isinstance(m,int) #python 3 - #assert isinstance(m,(int,long)) #python 2 - assert m>0 - return rng.randrange(m) - -#sample from a Bernoulli(p) distribution -#assumes p is a rational number in [0,1] -def sample_bernoulli(p,rng): - assert isinstance(p,Fraction) - assert 0 <= p <= 1 - m=sample_uniform(p.denominator,rng) - if m < p.numerator: - return 1 - else: - return 0 - -#sample from a Bernoulli(exp(-x)) distribution -#assumes x is a rational number in [0,1] -def sample_bernoulli_exp1(x,rng): - assert isinstance(x,Fraction) - assert 0 <= x <= 1 - k=1 - while True: - if sample_bernoulli(x/k,rng)==1: - k=k+1 - else: - break - return k%2 - -#sample from a Bernoulli(exp(-x)) distribution -#assumes x is a rational number >=0 -def sample_bernoulli_exp(x,rng): - assert isinstance(x,Fraction) - assert x >= 0 - #Sample floor(x) independent Bernoulli(exp(-1)) - #If all are 1, return Bernoulli(exp(-(x-floor(x)))) - while x>1: - if sample_bernoulli_exp1(Fraction(1,1),rng)==1: - x=x-1 - else: - return 0 - return sample_bernoulli_exp1(x,rng) - -#sample from a geometric(1-exp(-x)) distribution -#assumes x is a rational number >= 0 -def sample_geometric_exp_slow(x,rng): - assert isinstance(x,Fraction) - assert x >= 0 - k=0 - while True: - if sample_bernoulli_exp(x,rng)==1: - k=k+1 - else: - return k - -#sample from a geometric(1-exp(-x)) distribution -#assumes x >= 0 rational -def sample_geometric_exp_fast(x,rng): - assert isinstance(x,Fraction) - if x==0: return 0 #degenerate case - assert x>0 - - t=x.denominator - while True: - u=sample_uniform(t,rng) - b=sample_bernoulli_exp(Fraction(u,t),rng) - if b==1: - break - v=sample_geometric_exp_slow(Fraction(1,1),rng) - value = v*t+u - return value//x.numerator - -#sample from a discrete Laplace(scale) distribution -#Returns integer x with Pr[x] = exp(-abs(x)/scale)*(exp(1/scale)-1)/(exp(1/scale)+1) -#casts scale to Fraction -#assumes scale>=0 -def sample_dlaplace(scale,rng=None): - if rng is None: - rng = random.SystemRandom() - scale = Fraction(scale) - assert scale >= 0 - while True: - sign=sample_bernoulli(Fraction(1,2),rng) - magnitude=sample_geometric_exp_fast(1/scale,rng) - if sign==1 and magnitude==0: continue - return magnitude*(1-2*sign) - -#compute floor(sqrt(x)) exactly -#only requires comparisons between x and integer -def floorsqrt(x): - assert x >= 0 - #a,b integers - a=0 #maintain a^2<=x - b=1 #maintain b^2>x - while b*b <= x: - b=2*b #double to get upper bound - #now do binary search - while a+1=0 -def sample_dgauss(sigma2,rng=None): - if rng is None: - rng = random.SystemRandom() - sigma2=Fraction(sigma2) - if sigma2==0: return 0 #degenerate case - assert sigma2 > 0 - t = floorsqrt(sigma2)+1 - while True: - candidate = sample_dlaplace(t,rng=rng) - bias=((abs(candidate)-sigma2/t)**2)/(2*sigma2) - if sample_bernoulli_exp(bias,rng)==1: - return candidate - -######################################################################### -#DONE That's it! Now some utilities - -import math #need this, code below is no longer exact - -#Compute the normalizing constant of the discrete gaussian -#i.e. sum_{x in Z} exp(-x^2/2sigma2) -#By Poisson summation formula, this is equivalent to -# sqrt{2*pi*sigma2}*sum_{y in Z} exp(-2*pi^2*sigma2*y^2) -#For small sigma2 the former converges faster -#For large sigma2, the latter converges faster -#crossover at sigma2=1/2*pi -#For intermediate sigma2, this code will compute both and check -def normalizing_constant(sigma2): - original=None - poisson=None - if sigma2<=1: - original = 0 - x=1000 #summation stops at exp(-x^2/2sigma2)<=exp(-500,000) - while x>0: - original = original + math.exp(-x*x/(2.0*sigma2)) - x = x - 1 #sum from small to large for improved accuracy - original = 2*original + 1 #symmetrize and add x=0 - if sigma2*100 >= 1: - poisson = 0 - y = 1000 #summation stops at exp(-y^2*2*pi^2*sigma2)<=exp(-190,000) - while y>0: - poisson = poisson + math.exp(-math.pi*math.pi*sigma2*2*y*y) - y = y - 1 #sum from small to large - poisson = math.sqrt(2*math.pi*sigma2)*(1+2*poisson) - if poisson is None: return original - if original is None: return poisson - #if we have computed both, check equality - scale = max(1,math.sqrt(2*math.pi*sigma2)) #tight-ish lower bound on constant - assert -1e-15*scale <= original-poisson <= 1e-15*scale - #10^-15 is about as much precision as we can expect from double precision floating point numbers - #64-bit float has 56-bit mantissa 10^-15 ~= 2^-50 - return (original+poisson)/2 - -#compute the variance of discrete gaussian -#mean is zero, thus: -#var = sum_{x in Z} x^2*exp(-x^2/(2*sigma2)) / normalizing_constant(sigma2) -#By Poisson summation formula, we have equivalent expression: -# variance(sigma2) = sigma2 * (1 - 4*pi^2*sigma2*variance(1/(4*pi^2*sigma2)) ) -#See lemma 20 https://arxiv.org/pdf/2004.00010v3.pdf#page=17 -#alternative expression converges faster when sigma2 is large -#crossover point (in terms of convergence) is sigma2=1/(2*pi) -#for intermediate values of sigma2, we compute both expressions and check -def variance(sigma2): - original=None - poisson=None - if sigma2<=1: #compute primary expression - original=0 - x = 1000 #summation stops at exp(-x^2/2sigma2)<=exp(-500,000) - while x>0: #sum from small to large for improved accuracy - original = original + x*x*math.exp(-x*x/(2.0*sigma2)) - x=x-1 - original = 2*original/normalizing_constant(sigma2) - if sigma2*100>=1: - poisson=0 #we will compute sum_{y in Z} y^2 * exp(-2*pi^2*sigma2*y^2) - y=1000 #summation stops at exp(-y^2*2*pi^2*sigma2)<=exp(-190,000) - while y>0: #sum from small to large - poisson = poisson + y*y*math.exp(-y*y*2*sigma2*math.pi*math.pi) - y=y-1 - poisson = 2*poisson/normalizing_constant(1/(4*sigma2*math.pi*math.pi)) - #next convert from variance(1/(4*pi^2*sigma2)) to variance(sigma2) - poisson = sigma2*(1-4*sigma2*poisson*math.pi*math.pi) - if original is None: return poisson - if poisson is None: return original - #if we have computed both check equality - assert -1e-15*sigma2 <= original-poisson <= 1e-15*sigma2 - return (original+poisson)/2 - -######################################################################### -#DONE Now some basic testing code - -import matplotlib.pyplot as plt #only needed for testing -import time #only needed for testing - -#This generates n samples from sample_dgauss(sigma2) -#It times this and releases statistics -#produces a histogram plot if plot==True -#if plot==None it will only produce a histogram if it's not too large -#can save image instead of displaying by specifying a path e.g., save="plot.png" -def plot_histogram(sigma2,n,save=None,plot=None): - #generate samples - before=time.time() - samples = [sample_dgauss(sigma2) for i in range(n)] - after=time.time() - print("generated "+str(n)+" samples in "+str(after-before)+" seconds ("+str(n/(after-before))+" samples per second) for sigma^2="+str(sigma2)) - #now process - samples.sort() - values=[] - counts=[] - counter=None - prev=None - for sample in samples: - if prev is None: #initializing - prev=sample - counter=1 - elif sample==prev: #still same element - counter=counter+1 - else: - #add prev to histogram - values.append(prev) - counts.append(counter) - #start counting - prev=sample - counter=1 - #add final value - values.append(prev) - counts.append(counter) - - #print & sum - sum=0 - sumsquared=0 - kl=0 #compute KL divergence betwen empirical distribution and true distribution - norm_const=normalizing_constant(sigma2) - true_var=variance(sigma2) - for i in range(len(values)): - if len(values)<=100: #don't print too much - print(str(values[i])+":\t"+str(counts[i])) - sum = sum + values[i]*counts[i] - sumsquared = sumsquared + values[i]*values[i]*counts[i] - kl = kl + counts[i]*(math.log(counts[i]*norm_const/n)+values[i]*values[i]/(2.0*sigma2)) - mean = Fraction(sum,n) - var=Fraction(sumsquared,n) - kl=kl/n - print("mean="+str(float(mean))+" (true=0)") - print("variance="+str(float(var))+" (true="+str(true_var)+")") - print("KL(empirical||true)="+str(kl)) # https://en.wikipedia.org/wiki/G-test - assert kl>0 #kl divergence always >=0 and ==0 iff empirical==true, which is impossible - #now plot - if plot is None: - plot = (len(values)<=1000) #don't plot if huge - if not plot: return - ideal_counts = [n*math.exp(-x*x/(2.0*sigma2))/norm_const for x in values] - plt.bar(values, counts) - plt.plot(values, ideal_counts,'r') - plt.title("Histogram of samples from discrete Gaussian\nsigma^2="+str(sigma2)+" n="+str(n)) - if save is None: - plt.show() - else: - plt.savefig(save) - plt.clf() - -if __name__ == '__main__': - print("This is the discrete Gaussian sampler") - print("See the paper https://arxiv.org/abs/2004.00010") - print("Now running some basic testing code") - print("Start by calculating normalizing constant and variance for different values") - #some test code for normalizing_constant and variance functions - for sigma2 in [0.1**100,0.1**6,0.001,0.01,0.03,0.05,0.08,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,10,100,10**6,10**20,10**100]: - #internal asserts do some testing when 0.01<=sigma2<=1 - c=normalizing_constant(sigma2) - v=variance(sigma2) - #print - print("sigma^2="+str(sigma2) + ":\tnorm_const=" + str(c) + "=sqrt{2*pi}*sigma*" + str(c/math.sqrt(2*math.pi*sigma2)) + "\tvar=" + str(v)) - #print a few samples - #for i in range(20): print sample_dgauss(1) - #plot histogram and statistics - #includes timing - print("Now run the sampler") - print("Start with very large sigma^2=10^100 -- for timing purposes only") - plot_histogram(10**100,100000,plot=False) #large var, this will just be for timing - print("Now sigma^2=10 -- will display a histogram") - plot_histogram(10,100000) #small var, this will produce plot From c7fcfc6517691c74c2d1748fcad5c445f55cd89e Mon Sep 17 00:00:00 2001 From: stefan-aws Date: Fri, 15 Mar 2024 12:16:32 +0000 Subject: [PATCH 05/12] permissions --- build/py/run_benchmarks.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) mode change 100644 => 100755 build/py/run_benchmarks.sh diff --git a/build/py/run_benchmarks.sh b/build/py/run_benchmarks.sh old mode 100644 new mode 100755 index 0cadb2d9..791c90b1 --- a/build/py/run_benchmarks.sh +++ b/build/py/run_benchmarks.sh @@ -1,3 +1,3 @@ #1/bin/bash -PYTHONPATH=.:build/py/DafnyVMC-py:docs/py/IBM.py python3 docs/py/Benchmarks.py \ No newline at end of file +PYTHONPATH=.:build/py/DafnyVMC-py:docs/py/benchmarks/differential-privacy-library:docs/py/benchmarks/discrete-gaussian-differential-privacy python3 docs/py/benchmarks/benchmarks.py \ No newline at end of file From 28ee4242153fee7be250474081a4c9f8974ec369 Mon Sep 17 00:00:00 2001 From: stefan-aws Date: Fri, 15 Mar 2024 12:21:42 +0000 Subject: [PATCH 06/12] renaming --- build/py/{run_benchmarks.sh => run_gaussian_benchmarks.sh} | 2 +- docs/py/Benchmarks/{benchmarks.py => gaussian_benchmarks.py} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename build/py/{run_benchmarks.sh => run_gaussian_benchmarks.sh} (69%) rename docs/py/Benchmarks/{benchmarks.py => gaussian_benchmarks.py} (100%) diff --git a/build/py/run_benchmarks.sh b/build/py/run_gaussian_benchmarks.sh similarity index 69% rename from build/py/run_benchmarks.sh rename to build/py/run_gaussian_benchmarks.sh index 791c90b1..f5ad6af7 100755 --- a/build/py/run_benchmarks.sh +++ b/build/py/run_gaussian_benchmarks.sh @@ -1,3 +1,3 @@ #1/bin/bash -PYTHONPATH=.:build/py/DafnyVMC-py:docs/py/benchmarks/differential-privacy-library:docs/py/benchmarks/discrete-gaussian-differential-privacy python3 docs/py/benchmarks/benchmarks.py \ No newline at end of file +PYTHONPATH=.:build/py/DafnyVMC-py:docs/py/benchmarks/differential-privacy-library:docs/py/benchmarks/discrete-gaussian-differential-privacy python3 docs/py/benchmarks/gaussian_benchmarks.py \ No newline at end of file diff --git a/docs/py/Benchmarks/benchmarks.py b/docs/py/Benchmarks/gaussian_benchmarks.py similarity index 100% rename from docs/py/Benchmarks/benchmarks.py rename to docs/py/Benchmarks/gaussian_benchmarks.py From 1cae5349f2e54d0413094d2359e965055167ebf7 Mon Sep 17 00:00:00 2001 From: Stefan Zetzsche <120379523+stefan-aws@users.noreply.github.com> Date: Fri, 15 Mar 2024 13:19:37 +0000 Subject: [PATCH 07/12] Update docs/py/Benchmarks/gaussian_benchmarks.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Tancrède Lepoint --- docs/py/Benchmarks/gaussian_benchmarks.py | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/py/Benchmarks/gaussian_benchmarks.py b/docs/py/Benchmarks/gaussian_benchmarks.py index 074196d2..493ba7df 100644 --- a/docs/py/Benchmarks/gaussian_benchmarks.py +++ b/docs/py/Benchmarks/gaussian_benchmarks.py @@ -48,6 +48,7 @@ for i in range(1100): start_time = timeit.default_timer() + # The sampler is not directly accessible, so we call `.randomize(0)` instead, as it adds a noise drawn according to a discrete Gaussian to `0`. g.randomise(0) elapsed = timeit.default_timer() - start_time ibm2.append(elapsed) From 732b97f76910a87a477fa50619bfd55a6d4d2fd8 Mon Sep 17 00:00:00 2001 From: Stefan Zetzsche <120379523+stefan-aws@users.noreply.github.com> Date: Fri, 15 Mar 2024 13:19:45 +0000 Subject: [PATCH 08/12] Update docs/py/Benchmarks/gaussian_benchmarks.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Tancrède Lepoint --- docs/py/Benchmarks/gaussian_benchmarks.py | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/py/Benchmarks/gaussian_benchmarks.py b/docs/py/Benchmarks/gaussian_benchmarks.py index 493ba7df..7cd7e59a 100644 --- a/docs/py/Benchmarks/gaussian_benchmarks.py +++ b/docs/py/Benchmarks/gaussian_benchmarks.py @@ -64,7 +64,6 @@ ibm2_mean.append(ibm2.mean()*1000.0) ibm2_std.append(ibm2.std()*1000.0) -print(sigmas) ax1.plot(sigmas, vmc_mean, color='green', linewidth=1.0, label='VMC') ax1.fill_between(sigmas, numpy.array(vmc_mean)-0.5*numpy.array(vmc_std), numpy.array(vmc_mean)+0.5*numpy.array(vmc_std), From b3fb37207b064a0b119417e8e6c5af20010d2a52 Mon Sep 17 00:00:00 2001 From: Stefan Zetzsche <120379523+stefan-aws@users.noreply.github.com> Date: Fri, 15 Mar 2024 13:20:05 +0000 Subject: [PATCH 09/12] Update docs/py/Benchmarks/gaussian_benchmarks.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Tancrède Lepoint --- docs/py/Benchmarks/gaussian_benchmarks.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/py/Benchmarks/gaussian_benchmarks.py b/docs/py/Benchmarks/gaussian_benchmarks.py index 7cd7e59a..bd29dd20 100644 --- a/docs/py/Benchmarks/gaussian_benchmarks.py +++ b/docs/py/Benchmarks/gaussian_benchmarks.py @@ -27,6 +27,8 @@ ibm = [] ibm2= [] + # The GaussianDiscrete class does not expose the sampler directly, and needs to be instantiated with `(epsilon, delta)`. + # We access its `_scale` member to get the values `sigma`'s needed by `DafnyVMC` and `discretegauss`. g = GaussianDiscrete(epsilon=0.01 * epsilon_times_100, delta=0.00001) sigma = g._scale sigmas += [sigma] From cc51a2661fbbcdc102f2013993d34a22bfe4d9bf Mon Sep 17 00:00:00 2001 From: stefan-aws Date: Fri, 15 Mar 2024 13:27:19 +0000 Subject: [PATCH 10/12] new variable naming, fix labeling, randomize -> randomise --- docs/py/Benchmarks/gaussian_benchmarks.py | 38 +++++++++++------------ 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/docs/py/Benchmarks/gaussian_benchmarks.py b/docs/py/Benchmarks/gaussian_benchmarks.py index bd29dd20..ec6b9207 100644 --- a/docs/py/Benchmarks/gaussian_benchmarks.py +++ b/docs/py/Benchmarks/gaussian_benchmarks.py @@ -11,10 +11,10 @@ vmc_mean = [] vmc_std = [] -ibm_mean = [] -ibm_std = [] -ibm2_mean = [] -ibm2_std = [] +ibm_dgdp_mean = [] +ibm_dgdp_std = [] +ibm_dpl_mean = [] +ibm_dpl_std = [] fig,ax1 = plt.subplots() @@ -24,8 +24,8 @@ sigmas = [] for epsilon_times_100 in tqdm.tqdm(range(1, 500, 2)): vmc = [] - ibm = [] - ibm2= [] + ibm_dgdp = [] + ibm_dpl = [] # The GaussianDiscrete class does not expose the sampler directly, and needs to be instantiated with `(epsilon, delta)`. # We access its `_scale` member to get the values `sigma`'s needed by `DafnyVMC` and `discretegauss`. @@ -46,25 +46,25 @@ start_time = timeit.default_timer() discretegauss.sample_dgauss(sigma_squared, rng) elapsed = timeit.default_timer() - start_time - ibm.append(elapsed) + ibm_dgdp.append(elapsed) for i in range(1100): start_time = timeit.default_timer() - # The sampler is not directly accessible, so we call `.randomize(0)` instead, as it adds a noise drawn according to a discrete Gaussian to `0`. + # The sampler is not directly accessible, so we call `.randomise(0)` instead, as it adds a noise drawn according to a discrete Gaussian to `0`. g.randomise(0) elapsed = timeit.default_timer() - start_time - ibm2.append(elapsed) + ibm_dpl.append(elapsed) vmc = numpy.array(vmc[-1000:]) - ibm = numpy.array(ibm[-1000:]) - ibm2 = numpy.array(ibm2[-1000:]) + ibm_dgdp = numpy.array(ibm_dgdp[-1000:]) + ibm_dpl = numpy.array(ibm_dpl[-1000:]) vmc_mean.append(vmc.mean()*1000.0) vmc_std.append(vmc.std()*1000.0) - ibm_mean.append(ibm.mean()*1000.0) - ibm_std.append(ibm.std()*1000.0) - ibm2_mean.append(ibm2.mean()*1000.0) - ibm2_std.append(ibm2.std()*1000.0) + ibm_dgdp_mean.append(ibm_dgdp.mean()*1000.0) + ibm_dgdp_std.append(ibm_dgdp.std()*1000.0) + ibm_dpl_mean.append(ibm_dpl.mean()*1000.0) + ibm_dpl_std.append(ibm_dpl.std()*1000.0) ax1.plot(sigmas, vmc_mean, color='green', linewidth=1.0, label='VMC') @@ -72,13 +72,13 @@ alpha=0.2, facecolor='k', linewidth=2, linestyle='dashdot', antialiased=True) -ax1.plot(sigmas, ibm_mean, color='red', linewidth=1.0, label='IBM-DPL') -ax1.fill_between(sigmas, numpy.array(ibm_mean)-0.5*numpy.array(ibm_std), numpy.array(ibm_mean)+0.5*numpy.array(ibm_std), +ax1.plot(sigmas, ibm_dgdp_mean, color='red', linewidth=1.0, label='IBM-DGDP') +ax1.fill_between(sigmas, numpy.array(ibm_dgdp_mean)-0.5*numpy.array(ibm_dgdp_std), numpy.array(ibm_dgdp_mean)+0.5*numpy.array(ibm_dgdp_std), alpha=0.2, facecolor='y', linewidth=2, linestyle='dashdot', antialiased=True) -ax1.plot(sigmas, ibm2_mean, color='purple', linewidth=1.0, label='IBM-DGDP') -ax1.fill_between(sigmas, numpy.array(ibm2_mean)-0.5*numpy.array(ibm2_std), numpy.array(ibm2_mean)+0.5*numpy.array(ibm2_std), +ax1.plot(sigmas, ibm_dpl_mean, color='purple', linewidth=1.0, label='IBM-DPL') +ax1.fill_between(sigmas, numpy.array(ibm_dpl_mean)-0.5*numpy.array(ibm_dpl_std), numpy.array(ibm_dpl_mean)+0.5*numpy.array(ibm_dpl_std), alpha=0.2, facecolor='y', linewidth=2, linestyle='dashdot', antialiased=True) From b2fdb646c47b0112c02e7d3b7a73f97b512df70f Mon Sep 17 00:00:00 2001 From: stefan-aws Date: Fri, 15 Mar 2024 16:23:53 +0000 Subject: [PATCH 11/12] diagrams --- docs/py/Benchmarks/gaussian_benchmarks.py | 2 +- docs/py/Benchmarks/gaussian_diagrams.py | 45 +++++++++++++++++++++++ 2 files changed, 46 insertions(+), 1 deletion(-) create mode 100644 docs/py/Benchmarks/gaussian_diagrams.py diff --git a/docs/py/Benchmarks/gaussian_benchmarks.py b/docs/py/Benchmarks/gaussian_benchmarks.py index ec6b9207..40bc8b7f 100644 --- a/docs/py/Benchmarks/gaussian_benchmarks.py +++ b/docs/py/Benchmarks/gaussian_benchmarks.py @@ -86,5 +86,5 @@ ax1.set_ylabel("Sampling Time (ms)") plt.legend(loc = 'best') now = datetime.now() -filename = 'Benchmarks' + now.strftime("%H_%M_%S") + '.pdf' +filename = 'GaussianBenchmarks' + now.strftime("%H%M%S") + '.pdf' plt.savefig(filename) \ No newline at end of file diff --git a/docs/py/Benchmarks/gaussian_diagrams.py b/docs/py/Benchmarks/gaussian_diagrams.py new file mode 100644 index 00000000..2e460e72 --- /dev/null +++ b/docs/py/Benchmarks/gaussian_diagrams.py @@ -0,0 +1,45 @@ +import matplotlib.pyplot as plt +import secrets +from decimal import Decimal +from datetime import datetime +import DafnyVMC +import discretegauss +from diffprivlib.mechanisms import GaussianDiscrete + +fig, axs = plt.subplots(8, 3, figsize=(20, 20)) + +rng = secrets.SystemRandom() +r = DafnyVMC.Random() + +for i in range(8): + vmc_data = [] + ibm_dgdp_data = [] + ibm_dpl_data = [] + + epsilon_times_100 = 1 + (i**2)*2.5 + g = GaussianDiscrete(epsilon=0.01 * epsilon_times_100, delta=0.00001) + + sigma = g._scale + sigma_squared = sigma ** 2 + sigma_num, sigma_denom = Decimal(sigma).as_integer_ratio() + + title_vmc = 'VMC, Sigma = ' + str(sigma) + title_ibm_dgdp = 'IBM-DGDP, Sigma = ' + str(sigma) + title_ibm_dpl = 'IBM-DPL, Sigma = ' + str(sigma) + + for _ in range(100000): + vmc_data.append(r.DiscreteGaussianSample(sigma_num, sigma_denom)) + ibm_dgdp_data.append(discretegauss.sample_dgauss(sigma_squared, rng)) + ibm_dpl_data.append(g.randomise(0)) + + axs[i, 0].hist(vmc_data, color='lightgreen', ec='black', bins=50) + axs[i, 0].set_title(title_vmc) + axs[i, 1].hist(ibm_dgdp_data, color='lightgreen', ec='black', bins=50) + axs[i, 1].set_title(title_ibm_dgdp) + axs[i, 2].hist(ibm_dpl_data, color='lightgreen', ec='black', bins=50) + axs[i, 2].set_title(title_ibm_dpl) + +now = datetime.now() +filename = 'GaussianDiagrams' + now.strftime("%H%M%S") + '.pdf' +plt.subplots_adjust(wspace=0.4, hspace=0.4) +plt.savefig(filename) \ No newline at end of file From e0c7f16f1454c7d501db44db31a8fd1cddc4e3b4 Mon Sep 17 00:00:00 2001 From: stefan-aws Date: Fri, 15 Mar 2024 16:41:38 +0000 Subject: [PATCH 12/12] missed file --- build/py/run_gaussian_diagrams.sh | 3 +++ 1 file changed, 3 insertions(+) create mode 100755 build/py/run_gaussian_diagrams.sh diff --git a/build/py/run_gaussian_diagrams.sh b/build/py/run_gaussian_diagrams.sh new file mode 100755 index 00000000..0a86feee --- /dev/null +++ b/build/py/run_gaussian_diagrams.sh @@ -0,0 +1,3 @@ +#1/bin/bash + +PYTHONPATH=.:build/py/DafnyVMC-py:docs/py/benchmarks/differential-privacy-library:docs/py/benchmarks/discrete-gaussian-differential-privacy python3 docs/py/benchmarks/gaussian_diagrams.py \ No newline at end of file