From 3fe6a9aecc06638df4c7ff496a63f34ed779663c Mon Sep 17 00:00:00 2001
From: Yu-Xiang Wang <yuxiangw@users.noreply.github.com>
Date: Fri, 23 Jun 2023 09:35:15 -0700
Subject: [PATCH 01/10] removing specified bounds for brents

---
 autodp/converter.py                           | 14 ++--
 example/example_amplification_by_sampling.py  |  6 +-
 tutorials/tutorial_PATE_with_autoDP.ipynb     | 10 +--
 tutorials/tutorial_calibrator.ipynb           |  8 +--
 .../tutorial_fdp_of_basic_mechanisms.ipynb    |  2 +-
 tutorials/tutorial_new_api.ipynb              | 65 ++++++++++---------
 tutorials/tutorial_online_query_release.ipynb |  4 +-
 .../tutorial_private_deep_learning.ipynb      |  6 +-
 8 files changed, 58 insertions(+), 57 deletions(-)

diff --git a/autodp/converter.py b/autodp/converter.py
index 46621f7..9430dbb 100644
--- a/autodp/converter.py
+++ b/autodp/converter.py
@@ -92,7 +92,7 @@ def fun(x):  # the input the RDP's \alpha
                         result = utils.stable_logsumexp_two(term_1 - np.log(x)- np.log(delta),0)
                         return min(result*1.0/(x - 1), bbghs)
 
-                results = minimize_scalar(fun, method='Brent', bracket=(1, 2), bounds=[1, 100000])
+                results = minimize_scalar(fun, method='Brent', bracket=(1, 2))#, bounds=[1, 100000])
                 if results.success:
                    # print('delta', delta,'eps under rdp', results.fun)
                     return results.fun
@@ -146,7 +146,7 @@ def fun(x):  # the input the RDP's \alpha
                     else:
                         return np.log(1 / delta) / (x - 1) + rdp(x)
 
-            results = minimize_scalar(fun, method='Brent', bracket=(1,2), bounds=[1, alpha_max])
+            results = minimize_scalar(fun, method='Brent', bracket=(1,2))#, bounds=[1, alpha_max])
             if results.success:
                 return results.fun
             else:
@@ -250,7 +250,7 @@ def fun(alpha):
                 return -single_fdp(x)
 
         # This will use brent to start with 1,2.
-        results = minimize_scalar(fun, bracket=(0.5, 2), bounds=(0.5, alpha_max))
+        results = minimize_scalar(fun, bracket=(0.5, 2))#, bounds=(0.5, alpha_max))
         if results.success:
             return -results.fun
         else:
@@ -527,7 +527,7 @@ def fun(alpha):
                 return log_one_minus_fdp_alpha(logx)
 
         # This will use brent to start with 1,2.
-        results = minimize_scalar(fun, bracket=(0.5, 2), bounds=(0.5, alpha_max))
+        results = minimize_scalar(fun, bracket=(0.5, 2))#, bounds=(0.5, alpha_max))
         if results.success:
             return [results.fun, results.x]
         else:
@@ -720,9 +720,9 @@ def normal_equation_loglogx(loglogx):
             bound1 = np.log(-tmp - tmp**2 / 2 - tmp**3 / 6)
         else:
             bound1 = np.log(1-np.exp(fun1(np.log(1-delta))))
-        #results = minimize_scalar(normal_equation, bounds=[-np.inf,0], bracket=[-1,-2])
-        results = minimize_scalar(normal_equation, method="Bounded", bounds=[bound1,0],
-                                  options={'xatol': 1e-10, 'maxiter': 500, 'disp': 0})
+        results = minimize_scalar(normal_equation, bracket=[-1,-2])
+        #results = minimize_scalar(normal_equation, method="Bounded", bounds=[bound1,0],
+        #                          options={'xatol': 1e-10, 'maxiter': 500, 'disp': 0})
         if results.success:
             if abs(results.fun) > 1e-4 and abs(results.x)>1e-10:
                 # This means that we hit xatol (x is close to 0, but
diff --git a/example/example_amplification_by_sampling.py b/example/example_amplification_by_sampling.py
index 96f7319..09c269f 100644
--- a/example/example_amplification_by_sampling.py
+++ b/example/example_amplification_by_sampling.py
@@ -38,9 +38,9 @@
 
 
 # Now let's do subsampling. First we need to use replace-one version of the base mechanisms.
-gm1.replace_one = True
-gm2.replace_one = True
-SVT.replace_one = True
+gm1.neighboring = "replace_one"
+gm2.neighboring = "replace_one"
+SVT.neighboring = "replace_one"
 
 composed_subsampled_mech = compose([subsample(gm1,prob),
                                     subsample(gm2,prob),
diff --git a/tutorials/tutorial_PATE_with_autoDP.ipynb b/tutorials/tutorial_PATE_with_autoDP.ipynb
index 0b575be..843176d 100644
--- a/tutorials/tutorial_PATE_with_autoDP.ipynb
+++ b/tutorials/tutorial_PATE_with_autoDP.ipynb
@@ -54,7 +54,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "2.2540846502197414 1e-06\n",
+      "2.2540846502197396 1e-06\n",
       "4.886554117462211 1e-06\n"
      ]
     },
@@ -62,7 +62,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2555: RuntimeWarning: invalid value encountered in double_scalars\n",
+      "/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py:2884: RuntimeWarning: invalid value encountered in scalar divide\n",
       "  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom\n"
      ]
     }
@@ -111,8 +111,8 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "PATE_eps=2 {'sigma': 22.30476293897311} 1.999999977711552\n",
-      "PATE_eps=0.5 {'sigma': 80.57618519751232} 0.49999999740790046\n"
+      "PATE_eps=2 {'sigma': 22.304762938973138} 1.9999999777115467\n",
+      "PATE_eps=0.5 {'sigma': 80.57618519751239} 0.49999999740789647\n"
      ]
     }
    ],
@@ -177,7 +177,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.8.16"
   }
  },
  "nbformat": 4,
diff --git a/tutorials/tutorial_calibrator.ipynb b/tutorials/tutorial_calibrator.ipynb
index 19d691b..810e2c5 100644
--- a/tutorials/tutorial_calibrator.ipynb
+++ b/tutorials/tutorial_calibrator.ipynb
@@ -33,8 +33,8 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "GM {'sigma': 36.304688756269286} 0.10000000492560832\n",
-      "Ana_GM {'sigma': 36.304691899114694} 0.09999999565548537\n"
+      "GM {'sigma': 36.30468875626822} 0.10000000492561108\n",
+      "Ana_GM {'sigma': 36.304691899114694} 0.09999999565548348\n"
      ]
     }
    ],
@@ -80,7 +80,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2555: RuntimeWarning: invalid value encountered in double_scalars\n",
+      "/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py:2884: RuntimeWarning: invalid value encountered in scalar divide\n",
       "  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom\n"
      ]
     }
@@ -307,7 +307,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.8.16"
   }
  },
  "nbformat": 4,
diff --git a/tutorials/tutorial_fdp_of_basic_mechanisms.ipynb b/tutorials/tutorial_fdp_of_basic_mechanisms.ipynb
index d9de142..d6e3cee 100644
--- a/tutorials/tutorial_fdp_of_basic_mechanisms.ipynb
+++ b/tutorials/tutorial_fdp_of_basic_mechanisms.ipynb
@@ -321,7 +321,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.8.16"
   }
  },
  "nbformat": 4,
diff --git a/tutorials/tutorial_new_api.ipynb b/tutorials/tutorial_new_api.ipynb
index 8e24478..60ec276 100644
--- a/tutorials/tutorial_new_api.ipynb
+++ b/tutorials/tutorial_new_api.ipynb
@@ -62,7 +62,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2555: RuntimeWarning: invalid value encountered in double_scalars\n",
+      "/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py:2884: RuntimeWarning: invalid value encountered in scalar divide\n",
       "  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom\n"
      ]
     }
@@ -135,14 +135,14 @@
       "Generic composition: epsilon(delta) =  2.1309868424169824 , at delta =  1e-06\n",
       "Gaussian composition:  epsilon(delta) =  1.984273919801572 , at delta =  1e-06\n",
       "Generic composition: epsilon(delta) =  1.6484258167240666 , at delta =  0.0001\n",
-      "Gaussian composition:  epsilon(delta) =  1.4867384204500818 , at delta =  0.0001\n"
+      "Gaussian composition:  epsilon(delta) =  1.4867384204500813 , at delta =  0.0001\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2555: RuntimeWarning: invalid value encountered in double_scalars\n",
+      "/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py:2884: RuntimeWarning: invalid value encountered in scalar divide\n",
       "  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom\n"
      ]
     }
@@ -265,11 +265,11 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2555: RuntimeWarning: invalid value encountered in double_scalars\n",
+      "/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py:2884: RuntimeWarning: invalid value encountered in scalar divide\n",
       "  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom\n",
-      "/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2149: RuntimeWarning: invalid value encountered in double_scalars\n",
+      "/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py:2417: RuntimeWarning: invalid value encountered in scalar multiply\n",
       "  tmp2 = (x - v) * (fx - fw)\n",
-      "/Users/yuxiangw/Documents/bitbucket/autodp/autodp/autodp_core.py:233: RuntimeWarning: divide by zero encountered in log\n",
+      "/Users/yuxiangw/Documents/bitbucket/autodp/autodp/autodp_core.py:272: RuntimeWarning: divide by zero encountered in log\n",
       "  fdp = lambda x: 1 - np.exp(fun1(np.log(x)))\n"
      ]
     },
@@ -318,15 +318,17 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2555: RuntimeWarning: invalid value encountered in double_scalars\n",
-      "  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom\n"
+      "/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py:2884: RuntimeWarning: invalid value encountered in scalar divide\n",
+      "  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom\n",
+      "/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py:2417: RuntimeWarning: invalid value encountered in scalar multiply\n",
+      "  tmp2 = (x - v) * (fx - fw)\n"
      ]
     },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "[5.298525912188081, 4.728386984943314, 4.7283856025045115, 4.377178095681224, 4.377178095627265, 1.4142111312027956]\n"
+      "[5.298525912188081, 4.728386984943314, 4.7283856417872805, 4.377178095681228, 4.377178095762625, 1.4142111312027956]\n"
      ]
     },
     {
@@ -406,7 +408,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -436,9 +438,9 @@
     "\n",
     "\n",
     "# Now let's do subsampling. First we need to use replace-one version of the base mechanisms.\n",
-    "gm1.replace_one = True\n",
-    "gm2.replace_one = True\n",
-    "SVT.replace_one = True\n",
+    "gm1.neighboring = \"replace_one\"\n",
+    "gm2.neighboring = \"replace_one\"\n",
+    "SVT.neighboring = \"replace_one\"\n",
     "\n",
     "composed_subsampled_mech = compose([subsample(gm1,prob),\n",
     "                                    subsample(gm2,prob),\n",
@@ -453,17 +455,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 14,
    "metadata": {},
    "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2555: RuntimeWarning: invalid value encountered in double_scalars\n",
-      "  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom\n"
-     ]
-    },
     {
      "name": "stdout",
      "output_type": "stream",
@@ -482,7 +476,7 @@
       "epsilon(delta) =  0.48950668976148426 , at delta =  0.0001\n",
       "---------------------------------------------------\n",
       "Mechanism name is \" Compose:{Subsample:GM1: 30, Subsample:GM2: 50, Subsample:SVT: 10} \"\n",
-      "Parameters are:  {'Subsample:GM1:sigma': 5.0, 'Subsample:GM1:PoissonSample': 0.1, 'Subsample:GM1:Subsample': 0.1, 'Subsample:GM2:sigma': 8.0, 'Subsample:GM2:PoissonSample': 0.1, 'Subsample:GM2:Subsample': 0.1, 'Subsample:SVT:eps': 0.1, 'Subsample:SVT:PoissonSample': 0.1, 'Subsample:SVT:Subsample': 0.1}\n",
+      "Parameters are:  {'Subsample:GM1:sigma': 5.0, 'Subsample:GM1:Subsample': 0.1, 'Subsample:GM2:sigma': 8.0, 'Subsample:GM2:Subsample': 0.1, 'Subsample:SVT:eps': 0.1, 'Subsample:SVT:Subsample': 0.1}\n",
       "epsilon(delta) =  3.161674646617909 , at delta =  1e-06\n",
       "epsilon(delta) =  2.2643681643522973 , at delta =  0.0001\n",
       "------- If qualified for the improved bounds --------\n",
@@ -549,14 +543,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 28,
+   "execution_count": 15,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "GM1 {'sigma': 2.904058395116701} 1.4999997486031456\n",
+      "GM1 {'sigma': 2.904058395116711} 1.4999997486031442\n",
       "Laplace {'eps': 1.500001569629506} 1.500001346499671\n"
      ]
     }
@@ -593,7 +587,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 16,
    "metadata": {},
    "outputs": [
     {
@@ -602,7 +596,7 @@
        "2.0677358464697515"
       ]
      },
-     "execution_count": 14,
+     "execution_count": 16,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -624,7 +618,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 30,
+   "execution_count": 17,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -650,7 +644,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 31,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -686,7 +680,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 32,
+   "execution_count": 19,
    "metadata": {},
    "outputs": [
     {
@@ -748,7 +742,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 33,
+   "execution_count": 20,
    "metadata": {},
    "outputs": [
     {
@@ -769,6 +763,13 @@
     "print(online_ngd.get_approxDP(delta))"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -793,7 +794,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.8.16"
   }
  },
  "nbformat": 4,
diff --git a/tutorials/tutorial_online_query_release.ipynb b/tutorials/tutorial_online_query_release.ipynb
index b647abb..027b76b 100644
--- a/tutorials/tutorial_online_query_release.ipynb
+++ b/tutorials/tutorial_online_query_release.ipynb
@@ -57,7 +57,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2555: RuntimeWarning: invalid value encountered in double_scalars\n",
+      "/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py:2884: RuntimeWarning: invalid value encountered in scalar divide\n",
       "  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom\n"
      ]
     },
@@ -212,7 +212,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.8.16"
   }
  },
  "nbformat": 4,
diff --git a/tutorials/tutorial_private_deep_learning.ipynb b/tutorials/tutorial_private_deep_learning.ipynb
index 11f9dde..5c858ae 100644
--- a/tutorials/tutorial_private_deep_learning.ipynb
+++ b/tutorials/tutorial_private_deep_learning.ipynb
@@ -42,7 +42,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2555: RuntimeWarning: invalid value encountered in double_scalars\n",
+      "/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py:2884: RuntimeWarning: invalid value encountered in scalar divide\n",
       "  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom\n"
      ]
     },
@@ -186,7 +186,7 @@
     {
      "data": {
       "text/plain": [
-       "10.997151214220652"
+       "10.99715121422065"
       ]
      },
      "execution_count": 3,
@@ -383,7 +383,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.8.16"
   }
  },
  "nbformat": 4,

From e67b8ebbf0045b2ef3dd2747dd76c988f2d62287 Mon Sep 17 00:00:00 2001
From: Yu-Xiang Wang <yuxiangw@users.noreply.github.com>
Date: Fri, 23 Jun 2023 10:17:05 -0700
Subject: [PATCH 02/10] fixing find_logx from delta by reintroducing the
 carefully chosen bounds

---
 autodp/converter.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/autodp/converter.py b/autodp/converter.py
index 9430dbb..3c096f5 100644
--- a/autodp/converter.py
+++ b/autodp/converter.py
@@ -720,9 +720,9 @@ def normal_equation_loglogx(loglogx):
             bound1 = np.log(-tmp - tmp**2 / 2 - tmp**3 / 6)
         else:
             bound1 = np.log(1-np.exp(fun1(np.log(1-delta))))
-        results = minimize_scalar(normal_equation, bracket=[-1,-2])
-        #results = minimize_scalar(normal_equation, method="Bounded", bounds=[bound1,0],
-        #                          options={'xatol': 1e-10, 'maxiter': 500, 'disp': 0})
+        #results = minimize_scalar(normal_equation, bracket=[-1,-2])
+        results = minimize_scalar(normal_equation, method="Bounded", bounds=[bound1,0],
+                                  options={'xatol': 1e-10, 'maxiter': 500, 'disp': 0})
         if results.success:
             if abs(results.fun) > 1e-4 and abs(results.x)>1e-10:
                 # This means that we hit xatol (x is close to 0, but

From 827eb2bac1df26d33905af57af0a7922d54fb636 Mon Sep 17 00:00:00 2001
From: Yu-Xiang Wang <yuxiangw@users.noreply.github.com>
Date: Fri, 22 Sep 2023 01:06:38 -0700
Subject: [PATCH 03/10] adding the closed-form rdp2delta implementation

---
 autodp/autodp_core.py |  5 ++-
 autodp/converter.py   | 94 +++++++++++++++++++++++++------------------
 2 files changed, 59 insertions(+), 40 deletions(-)

diff --git a/autodp/autodp_core.py b/autodp/autodp_core.py
index b2577a0..c564e9f 100644
--- a/autodp/autodp_core.py
+++ b/autodp/autodp_core.py
@@ -200,7 +200,6 @@ def approx_dp_func(delta1):
         elif type_of_update == 'RDP':
             # function output RDP eps as a function of alpha
             self.RenyiDP = converter.pointwise_minimum(self.RenyiDP, func)
-            self.approx_delta = converter.pointwise_minimum(self.approx_delta, converter.rdp_to_delta(self.RenyiDP))
             if fDP_based_conversion:
 
                 fdp_log, fdp_grad_log = converter.rdp_to_fdp_and_fdp_grad_log(func)
@@ -234,12 +233,16 @@ def approx_dp_func(delta1):
                                                             converter.fdp_fdp_grad_to_approxdp(
                                                                 fdp_log, fdp_grad_log,
                                                                 log_flag=True))
+                self.approx_delta = converter.pointwise_minimum(self.approx_delta,
+                          converter.rdp_to_delta(self.RenyiDP, BBGHS_conversion=BBGHS_conversion))
 
                 # self.approxDP = converter.pointwise_minimum(self.approxDP,
                 #                                            converter.fdp_to_approxdp(self.fDP))
             else:
                 self.approxDP = converter.pointwise_minimum(self.approxDP,
                          converter.rdp_to_approxdp(self.RenyiDP, BBGHS_conversion=BBGHS_conversion))
+                self.approx_delta = converter.pointwise_minimum(self.approx_delta,
+                          converter.rdp_to_delta(self.RenyiDP, BBGHS_conversion=BBGHS_conversion))
                 self.fDP = converter.pointwise_maximum(self.fDP,
                                                        converter.approxdp_func_to_fdp(
                                                            self.approxDP))
diff --git a/autodp/converter.py b/autodp/converter.py
index 3c096f5..bd740f6 100644
--- a/autodp/converter.py
+++ b/autodp/converter.py
@@ -60,7 +60,7 @@ def approxdp(delta):
     return approxdp
 
 
-def rdp_to_delta(rdp):
+def rdp_to_delta(rdp, BBGHS_conversion=True):
     """
     from RDP to delta with a fixed epsilon
     """
@@ -70,48 +70,64 @@ def approx_delta(eps, naive=False):
         """
         approximate delta as a function of epsilon
         """
-
-        def get_eps(delta):
-            if delta == 0:
-                return rdp(np.inf)
+        def get_delta(alpha):
+            if BBGHS_conversion:
+                # Also by Canonne et al., 2020
+                return np.minimum(
+                    (np.exp((alpha-1)*(rdp(alpha)-eps))/alpha)*((1-1/alpha)**(alpha-1)), 1.0)
             else:
-                def fun(x):  # the input the RDP's \alpha
-                    if x <= 1:
-                        return np.inf
-                    else:
-
-                        if naive:
-                            return np.log(1 / delta) / (x - 1) + rdp(x)
-                        bbghs = np.maximum(rdp(x) + np.log((x-1)/x)
-                                          - (np.log(delta) + np.log(x))/(x-1), 0)
-                        """
-                        The following is for optimal conversion
-                        1/(alpha -1 )log(e^{(alpha-1)*rdp -1}/(alpha*delta) +1 )
-                        """
-                        sign, term_1= utils.stable_log_diff_exp((x-1)*rdp(x),0)
-                        result = utils.stable_logsumexp_two(term_1 - np.log(x)- np.log(delta),0)
-                        return min(result*1.0/(x - 1), bbghs)
-
-                results = minimize_scalar(fun, method='Brent', bracket=(1, 2))#, bounds=[1, 100000])
-                if results.success:
-                   # print('delta', delta,'eps under rdp', results.fun)
-                    return results.fun
-                else:
-                    return np.inf
+                # Naive conversion from Mironov
+                return np.minimum(np.exp((alpha-1)*(rdp(alpha)-eps)),1.0)
 
-
-        def err(delta):
-            current_eps = get_eps(delta)
-            #print('current delta', delta, 'eps', current_eps)
-            return abs(eps - current_eps)
-
-        results = minimize_scalar(err, method='bounded', bounds=[0, 0.1],options={'xatol':1e-14})
+        results = minimize_scalar(get_delta, method = 'Brent', bracket=(1,2))
         if results.success:
-            #print('results', results.x)
-            return results.x
+            # print('delta', delta,'eps under rdp', results.fun)
+            return results.fun
         else:
             print('not found')
-            return 1
+            return 1.0
+
+        # def get_eps(delta):
+        #     if delta == 0:
+        #         return rdp(np.inf)
+        #     else:
+        #         def fun(x):  # the input the RDP's \alpha
+        #             if x <= 1:
+        #                 return np.inf
+        #             else:
+        #
+        #                 if naive:
+        #                     return np.log(1 / delta) / (x - 1) + rdp(x)
+        #                 bbghs = np.maximum(rdp(x) + np.log((x-1)/x)
+        #                                   - (np.log(delta) + np.log(x))/(x-1), 0)
+        #                 """
+        #                 The following is for optimal conversion
+        #                 1/(alpha -1 )log(e^{(alpha-1)*rdp -1}/(alpha*delta) +1 )
+        #                 """
+        #                 sign, term_1= utils.stable_log_diff_exp((x-1)*rdp(x),0)
+        #                 result = utils.stable_logsumexp_two(term_1 - np.log(x)- np.log(delta),0)
+        #                 return min(result*1.0/(x - 1), bbghs)
+        #
+        #         results = minimize_scalar(fun, method='Brent', bracket=(1, 2))#, bounds=[1, 100000])
+        #         if results.success:
+        #            # print('delta', delta,'eps under rdp', results.fun)
+        #             return results.fun
+        #         else:
+        #             return np.inf
+        #
+        #
+        # def err(delta):
+        #     current_eps = get_eps(delta)
+        #     #print('current delta', delta, 'eps', current_eps)
+        #     return abs(eps - current_eps)
+        #
+        # results = minimize_scalar(err, method='bounded', bounds=[0, 0.1],options={'xatol':1e-14})
+        # if results.success:
+        #     #print('results', results.x)
+        #     return results.x
+        # else:
+        #     print('not found')
+        #     return 1
 
 
     return approx_delta
@@ -244,7 +260,7 @@ def fdp(x):
 
         def fun(alpha):
             if alpha < 0.5:
-                return np.inf
+                return 0
             else:
                 single_fdp = single_rdp_to_fdp(alpha, rdp(alpha))
                 return -single_fdp(x)

From 12390bd85d201d1e71868559e48f56161274923f Mon Sep 17 00:00:00 2001
From: Yu-Xiang Wang <yuxiangw@users.noreply.github.com>
Date: Fri, 3 Nov 2023 14:57:45 -0700
Subject: [PATCH 04/10] fixing naming convention and its use in GaussianSVT

---
 autodp/dp_bank.py       | 2 +-
 autodp/mechanism_zoo.py | 2 +-
 setup.py                | 4 ++--
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/autodp/dp_bank.py b/autodp/dp_bank.py
index 4fbcbaf..d4a1f1e 100644
--- a/autodp/dp_bank.py
+++ b/autodp/dp_bank.py
@@ -85,7 +85,7 @@ def fun(x):
     else:
         raise RuntimeError(f"Failed to find epsilon: {results.flag}")
 
-def eps_generalized_gaussian(x, sigma, delta,k, c, c_tilde):
+def get_eps_gaussian_svt(x, sigma, delta, k, c, c_tilde):
     """
     submodule for generalized SVT with Gaussian noise
     we want to partition c into [c/c'] parts, each part using (k choose c')
diff --git a/autodp/mechanism_zoo.py b/autodp/mechanism_zoo.py
index f0f1725..e03a61f 100644
--- a/autodp/mechanism_zoo.py
+++ b/autodp/mechanism_zoo.py
@@ -431,7 +431,7 @@ def __init__(self, params=None,approxDP_off=False, name='StageWiseMechanism'):
         self.delta0 = 0
 
         if not approxDP_off:  # Direct implementation of approxDP
-            new_approxdp = lambda x: dp_bank.get_generalized_gaussian(params, x)
+            new_approxdp = lambda x: dp_bank.get_gaussian_svt(params, x)
             self.propagate_updates(new_approxdp, 'approxDP_func')
 
 
diff --git a/setup.py b/setup.py
index a264b39..83d35b1 100644
--- a/setup.py
+++ b/setup.py
@@ -18,14 +18,14 @@ def _parse_requirements(path):
 
 setup(
    name='autodp',
-   version='0.2.1b',
+   version='0.2.3',
    description='Automating Differential Privacy Computation',
    license="Apache",
    long_description="The package helps researchers and developers to correctly use advanced methods in differential privacy and obtain provable DP guarantees. The core of the package is an analytical moments accountant that keeps track of Renyi Differential Privacy in analytical forms.",
    author='Yu-Xiang Wang',
    author_email='yuxiangw@cs.ucsb.edu',
    url='https://github.com/yuxiangw/autodp',
-   download_url = 'https://github.com/yuxiangw/autodp/archive/refs/tags/v0.2.1b.tar.gz',   
+   download_url = 'https://github.com/yuxiangw/autodp/archive/refs/tags/v0.2.3.tar.gz',   
    keywords = ['Differential Privacy','Moments Accountant','Renyi Differential Privacy'],
    packages=['autodp'],  #same as name
    install_requires=[_parse_requirements('requirements.txt')], #external packages as dependencies

From 0abd50c52ff5062d16ac75b2a54628928af90d49 Mon Sep 17 00:00:00 2001
From: Yu-Xiang Wang <yuxiangw@users.noreply.github.com>
Date: Fri, 3 Nov 2023 15:03:49 -0700
Subject: [PATCH 05/10] fixing function names

---
 autodp/mechanism_zoo.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/autodp/mechanism_zoo.py b/autodp/mechanism_zoo.py
index e03a61f..af008eb 100644
--- a/autodp/mechanism_zoo.py
+++ b/autodp/mechanism_zoo.py
@@ -431,7 +431,7 @@ def __init__(self, params=None,approxDP_off=False, name='StageWiseMechanism'):
         self.delta0 = 0
 
         if not approxDP_off:  # Direct implementation of approxDP
-            new_approxdp = lambda x: dp_bank.get_gaussian_svt(params, x)
+            new_approxdp = lambda x: dp_bank.get_eps_gaussian_svt(params, x)
             self.propagate_updates(new_approxdp, 'approxDP_func')
 
 

From 92daefd42db9b032bb58877593a2c19968b1ff7d Mon Sep 17 00:00:00 2001
From: Yuqing <yuqingzhu@ucsb.edu>
Date: Sun, 5 Nov 2023 23:43:16 -0800
Subject: [PATCH 06/10] fix issue in pull request rdp_dp

---
 autodp/dp_bank.py                             | 137 +++++-----
 autodp/mechanism_zoo.py                       |  11 +-
 .../tutorial_sparse_vector_technique.ipynb    | 233 ++++++++++++++++++
 3 files changed, 319 insertions(+), 62 deletions(-)
 create mode 100644 tutorials/tutorial_sparse_vector_technique.ipynb

diff --git a/autodp/dp_bank.py b/autodp/dp_bank.py
index d4a1f1e..bbfb778 100644
--- a/autodp/dp_bank.py
+++ b/autodp/dp_bank.py
@@ -21,8 +21,8 @@ def get_eps_rdp(func, delta):
     :param delta:
     :return: The corresponding epsilon
     """
-    assert(delta >= 0)
-    acct = rdp_acct.anaRDPacct(m=10,m_max=10)
+    assert (delta >= 0)
+    acct = rdp_acct.anaRDPacct(m=10, m_max=10)
     acct.compose_mechanism(func)
     return acct.get_eps(delta)
 
@@ -34,38 +34,38 @@ def get_eps_rdp_subsampled(func, delta, prob):
     :param delta:
     :return: The corresponding epsilon
     """
-    assert(delta >= 0)
-    assert(prob >=0)
-    if prob==0:
+    assert (delta >= 0)
+    assert (prob >= 0)
+    if prob == 0:
         return 0
     elif prob == 1:
-        return get_eps_rdp(func,delta)
+        return get_eps_rdp(func, delta)
     else:
         acct = rdp_acct.anaRDPacct()
-        acct.compose_subsampled_mechanism(func,prob)
+        acct.compose_subsampled_mechanism(func, prob)
         return acct.get_eps(delta)
 
 
 # Get the eps and delta for a single Gaussian mechanism
 def get_eps_gaussian(sigma, delta):
     """ This function calculates the eps for Gaussian Mech given sigma and delta"""
-    assert(delta >= 0)
-    func = lambda x: rdp_bank.RDP_gaussian({'sigma':sigma},x)
-    return get_eps_rdp(func,delta)
+    assert (delta >= 0)
+    func = lambda x: rdp_bank.RDP_gaussian({'sigma': sigma}, x)
+    return get_eps_rdp(func, delta)
 
 
-def get_logdelta_ana_gaussian(sigma,eps):
+def get_logdelta_ana_gaussian(sigma, eps):
     """ This function calculates the delta parameter for analytical gaussian mechanism given eps"""
-    assert(eps>=0)
+    assert (eps >= 0)
     s, mag = utils.stable_log_diff_exp(norm.logcdf(0.5 / sigma - eps * sigma),
-                                       eps + norm.logcdf(-0.5/sigma - eps * sigma))
+                                       eps + norm.logcdf(-0.5 / sigma - eps * sigma))
     return mag
 
 
 def get_eps_ana_gaussian(sigma, delta):
     """ This function calculates the gaussian mechanism given sigma and delta using analytical GM"""
     # Basically inverting the above function by solving a nonlinear equation
-    assert(delta >=0 and delta <=1)
+    assert (delta >= 0 and delta <= 1)
 
     if delta == 0:
         return np.inf
@@ -78,71 +78,94 @@ def fun(x):
         else:
             return get_logdelta_ana_gaussian(sigma, x) - np.log(delta)
 
-    eps_upperbound = 1/2/sigma**2+1/sigma*np.sqrt(2*np.log(1/delta))
-    results = root_scalar(fun,bracket=[0, eps_upperbound])
+    eps_upperbound = 1 / 2 / sigma ** 2 + 1 / sigma * np.sqrt(2 * np.log(1 / delta))
+    results = root_scalar(fun, bracket=[0, eps_upperbound])
     if results.converged:
         return results.root
     else:
         raise RuntimeError(f"Failed to find epsilon: {results.flag}")
 
-def get_eps_gaussian_svt(x, sigma, delta, k, c, c_tilde):
+
+# def get_eps_gaussian_svt(x, sigma, delta, k, c, c_tilde):
+def get_eps_gaussian_svt(params, delta):
     """
-    submodule for generalized SVT with Gaussian noise
+    Implements Theorem 12 (stage-wise length-capped Gaussian SVT) in Zhu et.al., NeurIPS-20.
+
     we want to partition c into [c/c'] parts, each part using (k choose c')
     need to check whether (k choose c') > log(1/delta')
-    k is the maximam number of queries to answer for each chunk
-    x is log delta for each chunk, it needs to be negative
+    k is the maximum number of queries to answer for each chunk
+    x is delta for each chunk, it needs to be negative
     :param x:
     :param sigma:
     :param delta:
     :return:
     """
-    acct = dp_acct.DP_acct()
-    per_delta = np.exp(x) # per_delta for each c' chunk
-    coeff = comb(k,c_tilde)
-    assert per_delta < 1.0/(coeff)
-    #compute the eps per step with 1/(sigma_1**2) + sqrt(2/simga_1**2 *(log k + log(1/epr_delta)))
-    # compose eps for each chunk
-    while c:
-        if c>= c_tilde:
-            c = c - c_tilde
-        else:
-            # the remaining part c // c_tilde
-            c = 0
-            c_tilde = c
-            coeff = comb(k, c_tilde)
-        per_eps = (1.0+c_tilde) / (2*sigma ** 2) + np.sqrt((1.0 +c_tilde) / (2*sigma ** 2) * (np.log(coeff) - x))
-        acct.update_DPlosses(per_eps, per_delta)
-
-    compose_eps = acct.get_eps(delta)
-    return compose_eps
-
-
-def get_eps_laplace(b,delta):
-    assert(delta >= 0)
-    func = lambda x: rdp_bank.RDP_laplace({'b':b},x)
-    return get_eps_rdp(func,delta)
+    # {'sigma': params['sigma'], 'k': params['k'], 'c': params['c']}
+    sigma = params['sigma']
+    k = params.get('k', 1)
+    c = params.get('c', 1)
+    c_tilde = params.get('c_tilde', int(np.sqrt(c)))  # default c_tilde is sqrt(c)
+    coeff = comb(k, c_tilde)
+
+    # the function below returns the composed epsilon if each chunk is assigned a e^x budget of delta (the delta'
+    # used in Theorem 12).
+    def search_per_delta_each_chunk(x, sigma, c, c_tilde, k, coeff):
+        acct = dp_acct.DP_acct()
+        per_delta = np.exp(x)  # per_delta for each c' chunk
+        assert per_delta < 1.0 / (coeff)
+        # compute the eps per step with 1/(sigma_1**2) + sqrt(2/simga_1**2 *(log k + log(1/epr_delta)))
+        # compose eps for each chunk
+        while c:
+            if c >= c_tilde:
+                c = c - c_tilde
+            else:
+                # the remaining part c // c_tilde
+                c = 0
+                c_tilde = c
+                coeff = comb(k, c_tilde)
+            per_eps = (1.0 + c_tilde) / (2 * sigma ** 2) + np.sqrt(
+                (1.0 + c_tilde) / (2 * sigma ** 2) * (np.log(coeff) - x))
+            acct.update_DPlosses(per_eps, per_delta)
+        compose_eps = acct.get_eps(delta)
+        return compose_eps  # composed epsilon over c/c' chunks.
+
+    fun = lambda x: search_per_delta_each_chunk(x, sigma, c, c_tilde, k, coeff)
+    # the bound of per chunk delta (delta').
+    const = 10
+    right_bound = min(delta / (c_tilde + 1), 1.0 / coeff)
+    left_bound =delta / (c_tilde * 10)
+    # ensures the left bound is smaller than the right bound
+    while 10 * left_bound > right_bound:
+        const = const * 10
+        left_bound = left_bound * 0.1
+    results = minimize_scalar(fun, method='Bounded', bounds=[np.log(left_bound), np.log(right_bound)], options={'disp': False})
+    return results.fun
+
+
+def get_eps_laplace(b, delta):
+    assert (delta >= 0)
+    func = lambda x: rdp_bank.RDP_laplace({'b': b}, x)
+    return get_eps_rdp(func, delta)
 
 
-def get_eps_randresp(p,delta):
-    assert(delta >= 0)
-    func = lambda x: rdp_bank.RDP_randresponse({'p':p},x)
+def get_eps_randresp(p, delta):
+    assert (delta >= 0)
+    func = lambda x: rdp_bank.RDP_randresponse({'p': p}, x)
     return get_eps_rdp(func, delta)
 
 
+def get_eps_randresp_optimal(p, delta):
+    assert (delta >= 0)
+    if p < 0.5:
+        p = 1 - p
 
-def get_eps_randresp_optimal(p,delta):
-    assert(delta >= 0)
-    if p<0.5:
-        p = 1-p
-
-    if p==0 or p==1:
+    if p == 0 or p == 1:
         return np.inf
 
-    eps_max = np.log(p) - np.log(1-p)
+    eps_max = np.log(p) - np.log(1 - p)
     if delta == 0:
         return eps_max
-    elif delta >= 2*p - 1:
+    elif delta >= 2 * p - 1:
         return 0.0
     else:
-        return np.log(p-delta) - np.log(1 - p)
+        return np.log(p - delta) - np.log(1 - p)
diff --git a/autodp/mechanism_zoo.py b/autodp/mechanism_zoo.py
index af008eb..0882ca3 100644
--- a/autodp/mechanism_zoo.py
+++ b/autodp/mechanism_zoo.py
@@ -379,7 +379,7 @@ class GaussianSVT_Mechanism(Mechanism):
     The Gaussian-based SVT mechanism is described in
     https://papers.nips.cc/paper/2020/file/e9bf14a419d77534105016f5ec122d62-Paper.pdf
 
-    The mechanism takes the parameter k and c.k is the maximum length before
+    The mechanism takes the parameter k and c. k is the maximum length before
     the algorithm stops. c is the cut-off parameter.  Setting rdp_c_1 = True implies that
     we use RDP-based Gaussian-SVT with c=1, else c>1.
 
@@ -401,15 +401,16 @@ def __init__(self,params,name='GaussianSVT', rdp_c_1=True):
 class LaplaceSVT_Mechanism(Mechanism):
     """
     Laplace SVT (c>=1) with a RDP description.
-    The mechanism takes the parameter k and sigma.
-    k is the maximum length before the algorithm stops
+    b: the noise scale to perturb the threshold and the query.
+    k: the maximum length before the algorithm steps. k could be infinite.
+
     We provide the RDP implementation and pure-DP implementation
     """
     def __init__(self,params,name='GaussianSVT'):
         Mechanism.__init__(self)
         self.name=name
-        self.params={'b':params['b'],'k':params['k'], 'c':params['c']}
-
+        valid_keys = ['b', 'k', 'c']
+        self.params = dict(filter(lambda tuple: tuple[0] in valid_keys, params.items()))
         new_rdp = lambda x: rdp_bank.RDP_svt_laplace(self.params, x)
         self.propagate_updates(new_rdp, 'RDP')
 
diff --git a/tutorials/tutorial_sparse_vector_technique.ipynb b/tutorials/tutorial_sparse_vector_technique.ipynb
new file mode 100644
index 0000000..8962440
--- /dev/null
+++ b/tutorials/tutorial_sparse_vector_technique.ipynb
@@ -0,0 +1,233 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# In this tutorial, we discuss sparse vector techniques.\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Sparse Vector Technique** (SVT) is one of the most fundamental algorithmic tools in differential privacy that\n",
+    "allows the algorithm to screen potentially an unbounded number of adaptively chosen queries while paying a cost of privacy\n",
+    "only for a small number of queries that passes a predefined threshold.\n",
+    "\n",
+    "In this tutorial, we revisit the classic SVT (based on Laplace noise) and recent new variants of SVT algorithms (e.g., SVT with Gaussian noise) from\n",
+    "[Zhu et.al., NeurIPS-20](https://proceedings.neurips.cc/paper/2020/hash/e9bf14a419d77534105016f5ec122d62-Abstract.html)\n",
+    "\n",
+    "\n",
+    "In SVT, the input is a stream of adaptively chosen queries $q_1, ..., q_k, ...$. The queries are provided with a sequence of thresholds\n",
+    "$T_1, ..., T_k, ...$. The goal of SVT is to release a binary vector $(\\perp, \\top)^k$ at every time $k$,\n",
+    "$\\top$ indicates that the correponding query answers $q_i(D)$ is above the threshold $T_i$ and $\\perp$ indivates below.\n",
+    "To release this vector differential privately, the classic SVT first perturbs the threshold with a Laplace noise\n",
+    "$\\rho$. Then each individual query $q_i(D)$ is perturbed with another Laplace noise $\\nu_i$ before comparing against the\n",
+    "perturbed threshold $T_i + \\rho$ to determine the binary decision, until the stopping condition --- the c-th $\\top$ arrives.\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from autodp.mechanism_zoo import ExactGaussianMechanism, PureDP_Mechanism\n",
+    "from autodp.transformer_zoo import Composition\n",
+    "import matplotlib.pyplot as plt\n",
+    "%matplotlib inline "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Let's first define the classic SVT.\n",
+    "\n",
+    "We define a Laplace-based SVT with c=1."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 23,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py:2884: RuntimeWarning: invalid value encountered in double_scalars\n",
+      "  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom\n",
+      "/Users/zhuyuqing/github_proj/YQ-autoDP/autodp/autodp/rdp_bank.py:363: RuntimeWarning: overflow encountered in sinh\n",
+      "  cdp_bound = np.sinh(alpha * tilde_eps) - np.sinh((alpha - 1) * tilde_eps)\n",
+      "/Users/zhuyuqing/github_proj/YQ-autoDP/autodp/autodp/rdp_bank.py:363: RuntimeWarning: invalid value encountered in double_scalars\n",
+      "  cdp_bound = np.sinh(alpha * tilde_eps) - np.sinh((alpha - 1) * tilde_eps)\n"
+     ]
+    }
+   ],
+   "source": [
+    "from autodp.mechanism_zoo import GaussianSVT_Mechanism, LaplaceSVT_Mechanism\n",
+    "import numpy as np\n",
+    "\n",
+    "# b denotes the Laplace noise scale of \\rho (\\nu = 2\\rho by default). k is the maximum length before the algorithm stop. k could\n",
+    "# be infinite in the classic SVT. When k is not required, we provide an improved RDP bound.\n",
+    "#\n",
+    "params = {'b':1,'k':100, 'c':1}\n",
+    "lap_svt = LaplaceSVT_Mechanism(params)\n",
+    "delta = 1e-5\n",
+    "eps = lap_svt.get_eps(delta)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##  We next provide a Gaussian-noise based SVT.\n",
+    "\n",
+    "Gaussian SVT adds Gaussian noise to perturb both the threshold and the query."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# setting rdp_c_1=True implies we use RDP-based Gaussian-SVT with c=1. sigma_nu is the noise scale added to each query and\n",
+    "# sigma is the noise scale added to the threshold.\n",
+    "params = {'rdp_c_1':True, 'sigma':1, 'sigma_nu':1, 'k':100}\n",
+    "gau_svt = GaussianSVT_Mechanism(params)\n",
+    "delta = 1e-5\n",
+    "eps = gau_svt.get_eps(delta)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## How can we compare two SVT variants?\n",
+    "\n",
+    "Given a predetermined privacy budget $(\\epsilon, \\delta)$ and the cut-off c, we compare the length each\n",
+    "SVT-like algorithm can screen before stopping.\n",
+    "\n",
+    "We estimate the length of answered queries with Negative Binomial Distribution. For example, in the case of Gaussian-SVT,\n",
+    "with the threshold noise $z \\sim N(0, \\sigma_1^2)$ and the query noise $\\nu \\sim N(0,\\sigma_2^2)$, denote K is\n",
+    "the number of queries answered when hits c. Notie that $K|\\rho=z$ follows a Negative Binomial\n",
+    "Distribution, $E[K|z]$ can be estimated with $\\frac{cF_\\nu(T+z)}{1-F_\\nu(T+z)}$, where\n",
+    "$F_\\nu$ is the CDF of the noise $\\nu$ and queries are all zeros.\n",
+    "\n",
+    "Let's consider setting the cut-off parameter $c=20$, the threshold T=700 for all queries and $\\delta=1e-10$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "KeyError",
+     "evalue": "'c'",
+     "output_type": "error",
+     "traceback": [
+      "\u001B[0;31m---------------------------------------------------------------------------\u001B[0m",
+      "\u001B[0;31mKeyError\u001B[0m                                  Traceback (most recent call last)",
+      "\u001B[0;32m<ipython-input-25-d5c35c8f8838>\u001B[0m in \u001B[0;36m<module>\u001B[0;34m\u001B[0m\n\u001B[1;32m      6\u001B[0m \u001B[0mgeneral_calibrate\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mgeneralized_eps_delta_calibrator\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m      7\u001B[0m \u001B[0mparams_lap\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0;34m{\u001B[0m\u001B[0;34m'b'\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;32mNone\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0;34m'k'\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;36m100\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;34m'c'\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;36m20\u001B[0m\u001B[0;34m}\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m----> 8\u001B[0;31m \u001B[0mlap_svt_mech\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mgeneral_calibrate\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mLaplaceSVT_Mechanism\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0meps\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mdelta\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;34m[\u001B[0m\u001B[0;36m0\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;36m1000\u001B[0m\u001B[0;34m]\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mparams\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0mparams\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mpara_name\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;34m'b'\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mname\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;34m'lap_SVT'\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m      9\u001B[0m \u001B[0mprint\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mlap_svt_mech\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mname\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mlap_svt_mech\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mparams\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mlap_svt_mech\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mget_approxDP\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mdelta\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m     10\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m~/github_proj/YQ-autoDP/autodp/autodp/autodp_core.py\u001B[0m in \u001B[0;36m__call__\u001B[0;34m(self, *args, **kwargs)\u001B[0m\n\u001B[1;32m    453\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    454\u001B[0m     \u001B[0;32mdef\u001B[0m \u001B[0m__call__\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mself\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;34m*\u001B[0m\u001B[0margs\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;34m**\u001B[0m\u001B[0mkwargs\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m--> 455\u001B[0;31m         \u001B[0;32mreturn\u001B[0m \u001B[0mself\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mcalibrate\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0;34m*\u001B[0m\u001B[0margs\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;34m**\u001B[0m\u001B[0mkwargs\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m",
+      "\u001B[0;32m~/github_proj/YQ-autoDP/autodp/autodp/calibrator_zoo.py\u001B[0m in \u001B[0;36mparam_from_eps_delta\u001B[0;34m(self, mech_class, eps, delta, bounds, params, para_name, name)\u001B[0m\n\u001B[1;32m     78\u001B[0m             \u001B[0;32mreturn\u001B[0m \u001B[0mabs\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0meps\u001B[0m\u001B[0;34m-\u001B[0m\u001B[0mget_eps\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mx\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mdelta\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m     79\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m---> 80\u001B[0;31m         \u001B[0mresults\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mminimize_scalar\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0merr\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mmethod\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;34m'bounded'\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mbounds\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0mbounds\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m     81\u001B[0m         \u001B[0;32mif\u001B[0m \u001B[0mresults\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0msuccess\u001B[0m \u001B[0;32mand\u001B[0m \u001B[0mresults\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mfun\u001B[0m \u001B[0;34m<\u001B[0m \u001B[0;36m1e-3\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m     82\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m/usr/local/lib/python3.8/site-packages/scipy/optimize/_minimize.py\u001B[0m in \u001B[0;36mminimize_scalar\u001B[0;34m(fun, bracket, bounds, args, method, tol, options)\u001B[0m\n\u001B[1;32m    911\u001B[0m             raise ValueError('The `bounds` parameter is mandatory for '\n\u001B[1;32m    912\u001B[0m                              'method `bounded`.')\n\u001B[0;32m--> 913\u001B[0;31m         \u001B[0;32mreturn\u001B[0m \u001B[0m_minimize_scalar_bounded\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mfun\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mbounds\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0margs\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;34m**\u001B[0m\u001B[0moptions\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m    914\u001B[0m     \u001B[0;32melif\u001B[0m \u001B[0mmeth\u001B[0m \u001B[0;34m==\u001B[0m \u001B[0;34m'golden'\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    915\u001B[0m         \u001B[0;32mreturn\u001B[0m \u001B[0m_minimize_scalar_golden\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mfun\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mbracket\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0margs\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;34m**\u001B[0m\u001B[0moptions\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py\u001B[0m in \u001B[0;36m_minimize_scalar_bounded\u001B[0;34m(func, bounds, args, xatol, maxiter, disp, **unknown_options)\u001B[0m\n\u001B[1;32m   2211\u001B[0m     \u001B[0mrat\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0me\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0;36m0.0\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m   2212\u001B[0m     \u001B[0mx\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mxf\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m-> 2213\u001B[0;31m     \u001B[0mfx\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mfunc\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mx\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;34m*\u001B[0m\u001B[0margs\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m   2214\u001B[0m     \u001B[0mnum\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0;36m1\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m   2215\u001B[0m     \u001B[0mfmin_data\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0;34m(\u001B[0m\u001B[0;36m1\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mxf\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mfx\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m~/github_proj/YQ-autoDP/autodp/autodp/calibrator_zoo.py\u001B[0m in \u001B[0;36merr\u001B[0;34m(x)\u001B[0m\n\u001B[1;32m     76\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m     77\u001B[0m         \u001B[0;32mdef\u001B[0m \u001B[0merr\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mx\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m---> 78\u001B[0;31m             \u001B[0;32mreturn\u001B[0m \u001B[0mabs\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0meps\u001B[0m\u001B[0;34m-\u001B[0m\u001B[0mget_eps\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mx\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mdelta\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m     79\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m     80\u001B[0m         \u001B[0mresults\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mminimize_scalar\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0merr\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mmethod\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;34m'bounded'\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mbounds\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0mbounds\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m~/github_proj/YQ-autoDP/autodp/autodp/calibrator_zoo.py\u001B[0m in \u001B[0;36mget_eps\u001B[0;34m(x, delta)\u001B[0m\n\u001B[1;32m     73\u001B[0m                 \u001B[0mparams\u001B[0m\u001B[0;34m[\u001B[0m\u001B[0mpara_name\u001B[0m\u001B[0;34m]\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mx\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m     74\u001B[0m                 \u001B[0mmech\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mmech_class\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mparams\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m---> 75\u001B[0;31m             \u001B[0;32mreturn\u001B[0m \u001B[0mmech\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mget_approxDP\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mdelta\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m     76\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m     77\u001B[0m         \u001B[0;32mdef\u001B[0m \u001B[0merr\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mx\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m~/github_proj/YQ-autoDP/autodp/autodp/autodp_core.py\u001B[0m in \u001B[0;36mget_approxDP\u001B[0;34m(self, delta)\u001B[0m\n\u001B[1;32m    111\u001B[0m     \u001B[0;32mdef\u001B[0m \u001B[0mget_approxDP\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mself\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mdelta\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    112\u001B[0m         \u001B[0;31m# Output eps as a function of delta\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m--> 113\u001B[0;31m         \u001B[0;32mreturn\u001B[0m \u001B[0mself\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mapproxDP\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mdelta\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m    114\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    115\u001B[0m     \u001B[0;32mdef\u001B[0m \u001B[0mget_approx_delta\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mself\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0meps\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m~/github_proj/YQ-autoDP/autodp/autodp/converter.py\u001B[0m in \u001B[0;36mmin_f1_f2\u001B[0;34m(x)\u001B[0m\n\u001B[1;32m   1132\u001B[0m \u001B[0;32mdef\u001B[0m \u001B[0mpointwise_minimum\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mf1\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mf2\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m   1133\u001B[0m     \u001B[0;32mdef\u001B[0m \u001B[0mmin_f1_f2\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mx\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m-> 1134\u001B[0;31m         \u001B[0;32mreturn\u001B[0m \u001B[0mnp\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mminimum\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mf1\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mx\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mf2\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mx\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m   1135\u001B[0m     \u001B[0;32mreturn\u001B[0m \u001B[0mmin_f1_f2\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m   1136\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m~/github_proj/YQ-autoDP/autodp/autodp/converter.py\u001B[0m in \u001B[0;36mapproxdp\u001B[0;34m(delta)\u001B[0m\n\u001B[1;32m    163\u001B[0m                         \u001B[0;32mreturn\u001B[0m \u001B[0mnp\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mlog\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0;36m1\u001B[0m \u001B[0;34m/\u001B[0m \u001B[0mdelta\u001B[0m\u001B[0;34m)\u001B[0m \u001B[0;34m/\u001B[0m \u001B[0;34m(\u001B[0m\u001B[0mx\u001B[0m \u001B[0;34m-\u001B[0m \u001B[0;36m1\u001B[0m\u001B[0;34m)\u001B[0m \u001B[0;34m+\u001B[0m \u001B[0mrdp\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mx\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    164\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m--> 165\u001B[0;31m             \u001B[0mresults\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mminimize_scalar\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mfun\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mmethod\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;34m'Brent'\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mbracket\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0;36m1\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0;36m2\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;31m#, bounds=[1, alpha_max])\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m    166\u001B[0m             \u001B[0;32mif\u001B[0m \u001B[0mresults\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0msuccess\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    167\u001B[0m                 \u001B[0;32mreturn\u001B[0m \u001B[0mresults\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mfun\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m/usr/local/lib/python3.8/site-packages/scipy/optimize/_minimize.py\u001B[0m in \u001B[0;36mminimize_scalar\u001B[0;34m(fun, bracket, bounds, args, method, tol, options)\u001B[0m\n\u001B[1;32m    906\u001B[0m         \u001B[0;32mreturn\u001B[0m \u001B[0mmethod\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mfun\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0margs\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0margs\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mbracket\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0mbracket\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mbounds\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0mbounds\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;34m**\u001B[0m\u001B[0moptions\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    907\u001B[0m     \u001B[0;32melif\u001B[0m \u001B[0mmeth\u001B[0m \u001B[0;34m==\u001B[0m \u001B[0;34m'brent'\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m--> 908\u001B[0;31m         \u001B[0;32mreturn\u001B[0m \u001B[0m_minimize_scalar_brent\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mfun\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mbracket\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0margs\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;34m**\u001B[0m\u001B[0moptions\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m    909\u001B[0m     \u001B[0;32melif\u001B[0m \u001B[0mmeth\u001B[0m \u001B[0;34m==\u001B[0m \u001B[0;34m'bounded'\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    910\u001B[0m         \u001B[0;32mif\u001B[0m \u001B[0mbounds\u001B[0m \u001B[0;32mis\u001B[0m \u001B[0;32mNone\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py\u001B[0m in \u001B[0;36m_minimize_scalar_brent\u001B[0;34m(func, brack, args, xtol, maxiter, disp, **unknown_options)\u001B[0m\n\u001B[1;32m   2604\u001B[0m                   full_output=True, maxiter=maxiter, disp=disp)\n\u001B[1;32m   2605\u001B[0m     \u001B[0mbrent\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mset_bracket\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mbrack\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m-> 2606\u001B[0;31m     \u001B[0mbrent\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0moptimize\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m   2607\u001B[0m     \u001B[0mx\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mfval\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mnit\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mnfev\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mbrent\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mget_result\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mfull_output\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;32mTrue\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m   2608\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py\u001B[0m in \u001B[0;36moptimize\u001B[0;34m(self)\u001B[0m\n\u001B[1;32m   2375\u001B[0m         \u001B[0;31m# set up for optimization\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m   2376\u001B[0m         \u001B[0mfunc\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mself\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mfunc\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m-> 2377\u001B[0;31m         \u001B[0mxa\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mxb\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mxc\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mfa\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mfb\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mfc\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mfuncalls\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mself\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mget_bracket_info\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m   2378\u001B[0m         \u001B[0m_mintol\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mself\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0m_mintol\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m   2379\u001B[0m         \u001B[0m_cg\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mself\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0m_cg\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py\u001B[0m in \u001B[0;36mget_bracket_info\u001B[0;34m(self)\u001B[0m\n\u001B[1;32m   2344\u001B[0m             \u001B[0mxa\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mxb\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mxc\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mfa\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mfb\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mfc\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mfuncalls\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mbracket\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mfunc\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0margs\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0margs\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m   2345\u001B[0m         \u001B[0;32melif\u001B[0m \u001B[0mlen\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mbrack\u001B[0m\u001B[0;34m)\u001B[0m \u001B[0;34m==\u001B[0m \u001B[0;36m2\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m-> 2346\u001B[0;31m             xa, xb, xc, fa, fb, fc, funcalls = bracket(func, xa=brack[0],\n\u001B[0m\u001B[1;32m   2347\u001B[0m                                                        xb=brack[1], args=args)\n\u001B[1;32m   2348\u001B[0m         \u001B[0;32melif\u001B[0m \u001B[0mlen\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mbrack\u001B[0m\u001B[0;34m)\u001B[0m \u001B[0;34m==\u001B[0m \u001B[0;36m3\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py\u001B[0m in \u001B[0;36mbracket\u001B[0;34m(func, xa, xb, args, grow_limit, maxiter)\u001B[0m\n\u001B[1;32m   2866\u001B[0m     \u001B[0m_verysmall_num\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0;36m1e-21\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m   2867\u001B[0m     \u001B[0mfa\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mfunc\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0;34m*\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mxa\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0;34m)\u001B[0m \u001B[0;34m+\u001B[0m \u001B[0margs\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m-> 2868\u001B[0;31m     \u001B[0mfb\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mfunc\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0;34m*\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mxb\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0;34m)\u001B[0m \u001B[0;34m+\u001B[0m \u001B[0margs\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m   2869\u001B[0m     \u001B[0;32mif\u001B[0m \u001B[0;34m(\u001B[0m\u001B[0mfa\u001B[0m \u001B[0;34m<\u001B[0m \u001B[0mfb\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m:\u001B[0m                      \u001B[0;31m# Switch so fa > fb\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m   2870\u001B[0m         \u001B[0mxa\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mxb\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mxb\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mxa\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m~/github_proj/YQ-autoDP/autodp/autodp/converter.py\u001B[0m in \u001B[0;36mfun\u001B[0;34m(x)\u001B[0m\n\u001B[1;32m    158\u001B[0m                 \u001B[0;32melse\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    159\u001B[0m                     \u001B[0;32mif\u001B[0m \u001B[0mBBGHS_conversion\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m--> 160\u001B[0;31m                         return np.maximum(rdp(x) + np.log((x-1)/x)\n\u001B[0m\u001B[1;32m    161\u001B[0m                                           - (np.log(delta) + np.log(x))/(x-1), 0)\n\u001B[1;32m    162\u001B[0m                     \u001B[0;32melse\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m~/github_proj/YQ-autoDP/autodp/autodp/converter.py\u001B[0m in \u001B[0;36mmin_f1_f2\u001B[0;34m(x)\u001B[0m\n\u001B[1;32m   1132\u001B[0m \u001B[0;32mdef\u001B[0m \u001B[0mpointwise_minimum\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mf1\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mf2\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m   1133\u001B[0m     \u001B[0;32mdef\u001B[0m \u001B[0mmin_f1_f2\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mx\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m-> 1134\u001B[0;31m         \u001B[0;32mreturn\u001B[0m \u001B[0mnp\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mminimum\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mf1\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mx\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mf2\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mx\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m   1135\u001B[0m     \u001B[0;32mreturn\u001B[0m \u001B[0mmin_f1_f2\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m   1136\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m~/github_proj/YQ-autoDP/autodp/autodp/mechanism_zoo.py\u001B[0m in \u001B[0;36m<lambda>\u001B[0;34m(x)\u001B[0m\n\u001B[1;32m    411\u001B[0m         \u001B[0mself\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mname\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0mname\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    412\u001B[0m         \u001B[0mvalid_keys\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0;34m[\u001B[0m\u001B[0;34m'b'\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;34m'k'\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;34m'c'\u001B[0m\u001B[0;34m]\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m--> 413\u001B[0;31m         \u001B[0mself\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mparams\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mdict\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mfilter\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0;32mlambda\u001B[0m \u001B[0mtuple\u001B[0m\u001B[0;34m:\u001B[0m \u001B[0mtuple\u001B[0m\u001B[0;34m[\u001B[0m\u001B[0;36m0\u001B[0m\u001B[0;34m]\u001B[0m \u001B[0;32min\u001B[0m \u001B[0mvalid_keys\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mparams\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mitems\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m    414\u001B[0m         \u001B[0mnew_rdp\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0;32mlambda\u001B[0m \u001B[0mx\u001B[0m\u001B[0;34m:\u001B[0m \u001B[0mrdp_bank\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mRDP_svt_laplace\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mself\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mparams\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mx\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    415\u001B[0m         \u001B[0mself\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mpropagate_updates\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mnew_rdp\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;34m'RDP'\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;32m~/github_proj/YQ-autoDP/autodp/autodp/rdp_bank.py\u001B[0m in \u001B[0;36mRDP_svt_laplace\u001B[0;34m(params, alpha)\u001B[0m\n\u001B[1;32m    341\u001B[0m     \u001B[0mb\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mparams\u001B[0m\u001B[0;34m[\u001B[0m\u001B[0;34m'b'\u001B[0m\u001B[0;34m]\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    342\u001B[0m     \u001B[0mk\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mparams\u001B[0m\u001B[0;34m[\u001B[0m\u001B[0;34m'k'\u001B[0m\u001B[0;34m]\u001B[0m  \u001B[0;31m# the algorithm stops either k is achieved or c is achieved\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m--> 343\u001B[0;31m     \u001B[0mc\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mmax\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mparams\u001B[0m\u001B[0;34m[\u001B[0m\u001B[0;34m'c'\u001B[0m\u001B[0;34m]\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0;36m1\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m    344\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    345\u001B[0m     \u001B[0malpha\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0;36m1.0\u001B[0m \u001B[0;34m*\u001B[0m \u001B[0malpha\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
+      "\u001B[0;31mKeyError\u001B[0m: 'c'"
+     ]
+    }
+   ],
+   "source": [
+    "# calibrate noise with a given privacy budget.\n",
+    "from autodp.calibrator_zoo import eps_delta_calibrator,generalized_eps_delta_calibrator\n",
+    "from autodp.mechanism_zoo import LaplaceSVT_Mechanism, GaussianSVT_Mechanism\n",
+    "\n",
+    "eps = 0.1\n",
+    "general_calibrate = generalized_eps_delta_calibrator()\n",
+    "params_lap = {'b':None,'k':100, 'c':20}\n",
+    "lap_svt_mech = general_calibrate(LaplaceSVT_Mechanism, eps, delta, [0, 1000], params=params, para_name='b', name='lap_SVT')\n",
+    "print(lap_svt_mech.name, lap_svt_mech.params, lap_svt_mech.get_approxDP(delta))\n",
+    "\n",
+    "\n",
+    "print(f'in algorithm f eps,delta = ({eps},{delta}) ==> Noise level b=', params['b'])\n",
+    "lambda_rho = params['b']\n",
+    "ret_number = []\n",
+    "lambda_nu = 2*lambda_rho # this is the default choice of lambda_nu used in SVT.\n",
+    "repeat_n = 10000\n",
+    "margin = 700 # the threshold\n",
+    "c = 20\n",
+    "# estimate the empirical mean of the length of the answered queries.\n",
+    "for i in range(repeat_n):\n",
+    "    rho = np.random.laplace(loc=0, scale=lambda_rho)\n",
+    "    fail_prob = 0.5 * np.exp(-(margin + rho) / lambda_nu)\n",
+    "    expect_steps = c * (1 - fail_prob) / fail_prob\n",
+    "    ret_number.append(expect_steps)\n",
+    "ret_laplace = sum(ret_number) * 1.0 / repeat_n\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
\ No newline at end of file

From ebd5a242b2ca7d11c1e7f540bc79748188a76010 Mon Sep 17 00:00:00 2001
From: Yuqing <yuqingzhu@ucsb.edu>
Date: Sun, 5 Nov 2023 23:53:54 -0800
Subject: [PATCH 07/10] fix signature

---
 autodp/transformer_zoo.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/autodp/transformer_zoo.py b/autodp/transformer_zoo.py
index fb0751d..d3c7f9a 100644
--- a/autodp/transformer_zoo.py
+++ b/autodp/transformer_zoo.py
@@ -148,7 +148,7 @@ def __init__(self):
         Composition.__init__(self)
         self.name = 'ComposeGaussian'
 
-    def compose(self, mechanism_list, coeff_list):
+    def compose(self, mechanism_list, coeff_list, RDP_compose_only=True, BBGHS_conversion=True, fDP_based_conversion=True):
         # Make sure that the list contains only Gaussian mechanisms
         for mech in mechanism_list:
             assert(isinstance(mech, mechanism_zoo.GaussianMechanism)

From 270f4523476f1e86c063191c284a0063450c1abb Mon Sep 17 00:00:00 2001
From: Yuqing <yuqingzhu@ucsb.edu>
Date: Sun, 5 Nov 2023 23:59:51 -0800
Subject: [PATCH 08/10] fix again

---
 autodp/transformer_zoo.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/autodp/transformer_zoo.py b/autodp/transformer_zoo.py
index d3c7f9a..7346474 100644
--- a/autodp/transformer_zoo.py
+++ b/autodp/transformer_zoo.py
@@ -141,14 +141,14 @@ def update_params(self, mechanism_list):
 
 
 # composition of only Gaussian mechanisms
-class ComposeGaussian(Composition):
+class ComposeGaussian(Transformer):
     """ CompositionGaussian is a specialized composation function of ONLY Guassian mechanisms
     output a Mechanism that represents the composed mechanism"""
     def __init__(self):
         Composition.__init__(self)
         self.name = 'ComposeGaussian'
 
-    def compose(self, mechanism_list, coeff_list, RDP_compose_only=True, BBGHS_conversion=True, fDP_based_conversion=True):
+    def compose(self, mechanism_list, coeff_list):
         # Make sure that the list contains only Gaussian mechanisms
         for mech in mechanism_list:
             assert(isinstance(mech, mechanism_zoo.GaussianMechanism)

From 81c410f2e9c261eac2e6928f9d7aa441fb2ee102 Mon Sep 17 00:00:00 2001
From: Yuqing <yuqingzhu@ucsb.edu>
Date: Mon, 6 Nov 2023 00:06:39 -0800
Subject: [PATCH 09/10] fix again +1

---
 autodp/transformer_zoo.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/autodp/transformer_zoo.py b/autodp/transformer_zoo.py
index 7346474..a27e8ed 100644
--- a/autodp/transformer_zoo.py
+++ b/autodp/transformer_zoo.py
@@ -141,14 +141,14 @@ def update_params(self, mechanism_list):
 
 
 # composition of only Gaussian mechanisms
-class ComposeGaussian(Transformer):
+class ComposeGaussian(Composition):
     """ CompositionGaussian is a specialized composation function of ONLY Guassian mechanisms
     output a Mechanism that represents the composed mechanism"""
     def __init__(self):
         Composition.__init__(self)
         self.name = 'ComposeGaussian'
 
-    def compose(self, mechanism_list, coeff_list):
+    def compose(self, mechanism_list, coeff_list, RDP_compose_only=True, BBGHS_conversion=True, fDP_based_conversion=False):
         # Make sure that the list contains only Gaussian mechanisms
         for mech in mechanism_list:
             assert(isinstance(mech, mechanism_zoo.GaussianMechanism)

From fae0105c5631e106330df7c9d23195edf0acd458 Mon Sep 17 00:00:00 2001
From: Yuqing <yuqingzhu@ucsb.edu>
Date: Tue, 14 Nov 2023 00:41:45 -0800
Subject: [PATCH 10/10] fix logpq in phi_bank

---
 test/unit_test_sampling_pld.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/unit_test_sampling_pld.py b/test/unit_test_sampling_pld.py
index 71ef81a..3e64b5e 100644
--- a/test/unit_test_sampling_pld.py
+++ b/test/unit_test_sampling_pld.py
@@ -17,7 +17,7 @@ def _pld_sample_remove_only(sigma):
     t_list = [1.1**x for x in range(2, 20)]
     # Phi-function implementation
     gm = GaussianMechanism(sigma, name='GM3', RDP_off=True, approxDP_off=True, phi_off=False)
-    phi_gm = lambda t: gm.log_phi_p(t)
+
     # sampling ratio
     gamma = 0.01
     # direct compute the phi function of subsample gaussian for "removal only" neighboring relationship