Repository: optimagic-dev/optimagic (language: python)

Issue:
Allow pytrees in bootstrap
Allow bootstrap outcomes (statistics) to be general pytrees. This closes #318.


Discussion / hints:

Fix patch:
diff --git a/docs/source/how_to_guides/inference/how_to_do_bootstrap_inference.ipynb b/docs/source/how_to_guides/inference/how_to_do_bootstrap_inference.ipynb
index 419ae1cbf..fe9639ea1 100644
--- a/docs/source/how_to_guides/inference/how_to_do_bootstrap_inference.ipynb
+++ b/docs/source/how_to_guides/inference/how_to_do_bootstrap_inference.ipynb
@@ -3,43 +3,40 @@
   {
    "cell_type": "code",
    "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "import pandas as pd\n",
     "import seaborn as sns\n",
     "import statsmodels.api as sm"
-   ],
-   "outputs": [],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "# Bootstrap Tutorial\n",
     "\n",
-    "This notebook contains a tutorial on how to use the bootstrap functionality provided by estimagic. We start with the simplest possible example of calculating standard errors and confidence intervals for an OLS estimator without and clustering. Then we progress to more advanced examples.\n",
-    "\n",
-    "In this example, we will work with the \"exercise\" example dataset taken from the seaborn library.\n",
+    "This notebook contains a tutorial on how to use the bootstrap functionality provided by estimagic. We start with the simplest possible example of calculating standard errors and confidence intervals for an OLS estimator without as well as with clustering. Then we progress to more advanced examples.\n",
     "\n",
-    "The working example will be a linear regression to investigate the effects of exercise time on pulse.\n",
+    "In the example here, we will work with the \"exercise\" example dataset taken from the seaborn library.\n",
     "\n",
+    "The working example will be a linear regression to investigate the effects of exercise time on pulse."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
     "## Prepare the dataset"
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": 2,
-   "source": [
-    "df = sns.load_dataset(\"exercise\", index_col=0)\n",
-    "replacements = {\"1 min\": 1, \"15 min\": 15, \"30 min\": 30}\n",
-    "df = df.replace({\"time\": replacements})\n",
-    "df[\"constant\"] = 1\n",
-    "\n",
-    "df.head()"
-   ],
+   "metadata": {},
    "outputs": [
     {
-     "output_type": "execute_result",
      "data": {
       "text/html": [
        "<div>\n",
@@ -119,305 +116,249 @@
        "</div>"
       ],
       "text/plain": [
-       "   id     diet  pulse  time  kind  constant\n",
-       "0   1  low fat     85     1  rest         1\n",
-       "1   1  low fat     85    15  rest         1\n",
-       "2   1  low fat     88    30  rest         1\n",
-       "3   2  low fat     90     1  rest         1\n",
-       "4   2  low fat     92    15  rest         1"
+       "   id     diet  pulse time  kind  constant\n",
+       "0   1  low fat     85    1  rest         1\n",
+       "1   1  low fat     85   15  rest         1\n",
+       "2   1  low fat     88   30  rest         1\n",
+       "3   2  low fat     90    1  rest         1\n",
+       "4   2  low fat     92   15  rest         1"
       ]
      },
+     "execution_count": 2,
      "metadata": {},
-     "execution_count": 2
+     "output_type": "execute_result"
     }
    ],
-   "metadata": {}
+   "source": [
+    "df = sns.load_dataset(\"exercise\", index_col=0)\n",
+    "replacements = {\"1 min\": 1, \"15 min\": 15, \"30 min\": 30}\n",
+    "df = df.replace({\"time\": replacements})\n",
+    "df[\"constant\"] = 1\n",
+    "\n",
+    "df.head()"
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "## Doing a very simple bootstrap"
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "The first thing we need is a function that calculates the bootstrap outcome, given an empirical or re-sampled dataset. The bootstrap outcome is the quantity for which you want to calculate standard errors and confidence intervals. In most applications those are just parameter estimates.\n",
     "\n",
-    "In our case, we want to regress \"pulse\" on \"time\" and a constant, and the function looks as follows:"
-   ],
-   "metadata": {}
+    "In our case, we want to regress \"pulse\" on \"time\" and a constant. Our outcome function looks as follows:"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "def ols_fit(data):\n",
-    "\n",
     "    y = data[\"pulse\"]\n",
     "    x = data[[\"constant\", \"time\"]]\n",
     "    params = sm.OLS(y, x).fit().params\n",
+    "\n",
     "    return params"
-   ],
-   "outputs": [],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
-    "The outcome function must return a pandas Series.\n",
-    "\n",
-    "Now we can calculate confidence intervals and standard errors."
-   ],
-   "metadata": {}
+    "In general, the user-specified outcome function may return any pytree (e.g. numpy.ndarray, pandas.DataFrame, dict etc.). In the example here, it returns a pandas.Series."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we are ready to calculate confidence intervals and standard errors."
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(constant    90.889924\n",
+       " time         0.143938\n",
+       " dtype: float64,\n",
+       " constant    96.701376\n",
+       " time         0.647134\n",
+       " dtype: float64)"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
     "from estimagic.inference import bootstrap\n",
     "\n",
     "results_without_cluster = bootstrap(data=df, outcome=ols_fit)\n",
-    "results_without_cluster[\"summary\"]"
-   ],
+    "results_without_cluster.ci()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
    "outputs": [
     {
-     "output_type": "execute_result",
      "data": {
-      "text/html": [
-       "<div>\n",
-       "<style scoped>\n",
-       "    .dataframe tbody tr th:only-of-type {\n",
-       "        vertical-align: middle;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe tbody tr th {\n",
-       "        vertical-align: top;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe thead th {\n",
-       "        text-align: right;\n",
-       "    }\n",
-       "</style>\n",
-       "<table border=\"1\" class=\"dataframe\">\n",
-       "  <thead>\n",
-       "    <tr style=\"text-align: right;\">\n",
-       "      <th></th>\n",
-       "      <th>mean</th>\n",
-       "      <th>std</th>\n",
-       "      <th>lower_ci</th>\n",
-       "      <th>upper_ci</th>\n",
-       "    </tr>\n",
-       "  </thead>\n",
-       "  <tbody>\n",
-       "    <tr>\n",
-       "      <th>constant</th>\n",
-       "      <td>93.739117</td>\n",
-       "      <td>1.480086</td>\n",
-       "      <td>90.956284</td>\n",
-       "      <td>96.627072</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>time</th>\n",
-       "      <td>0.391281</td>\n",
-       "      <td>0.125367</td>\n",
-       "      <td>0.152010</td>\n",
-       "      <td>0.650922</td>\n",
-       "    </tr>\n",
-       "  </tbody>\n",
-       "</table>\n",
-       "</div>"
-      ],
       "text/plain": [
-       "               mean       std   lower_ci   upper_ci\n",
-       "constant  93.739117  1.480086  90.956284  96.627072\n",
-       "time       0.391281  0.125367   0.152010   0.650922"
+       "constant    1.491852\n",
+       "time        0.128982\n",
+       "dtype: float64"
       ]
      },
+     "execution_count": 5,
      "metadata": {},
-     "execution_count": 4
+     "output_type": "execute_result"
     }
    ],
-   "metadata": {}
+   "source": [
+    "results_without_cluster.se()"
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
-    "The above function call represents the minimum that a user has to specify, making full use of the default options, such as drawing a 1000 bootstrap draws, using the \"percentile\" bootstrap confidence interval, not making use of parallelization, etc.\n",
+    "The above function call represents the minimum that a user has to specify, making full use of the default options, such as drawing a 1_000 bootstrap draws, using the \"percentile\" bootstrap confidence interval, not making use of parallelization, etc.\n",
     "\n",
-    "If we would for example want to make 10000 draws, while parallelizing on two cores and using a \"bca\" type confidence interval, we would simply call the following."
-   ],
-   "metadata": {}
+    "If we would, for example, want to take 10_000 draws, while parallelizing on two cores, and using a \"bc\" type confidence interval, we would simply call the following:"
+   ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
-   "source": [
-    "results_without_cluster2 = bootstrap(\n",
-    "    data=df, outcome=ols_fit, n_draws=10_000, ci_method=\"bca\", n_cores=2\n",
-    ")\n",
-    "\n",
-    "results_without_cluster2[\"summary\"]"
-   ],
+   "execution_count": 6,
+   "metadata": {},
    "outputs": [
     {
-     "output_type": "execute_result",
      "data": {
-      "text/html": [
-       "<div>\n",
-       "<style scoped>\n",
-       "    .dataframe tbody tr th:only-of-type {\n",
-       "        vertical-align: middle;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe tbody tr th {\n",
-       "        vertical-align: top;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe thead th {\n",
-       "        text-align: right;\n",
-       "    }\n",
-       "</style>\n",
-       "<table border=\"1\" class=\"dataframe\">\n",
-       "  <thead>\n",
-       "    <tr style=\"text-align: right;\">\n",
-       "      <th></th>\n",
-       "      <th>mean</th>\n",
-       "      <th>std</th>\n",
-       "      <th>lower_ci</th>\n",
-       "      <th>upper_ci</th>\n",
-       "    </tr>\n",
-       "  </thead>\n",
-       "  <tbody>\n",
-       "    <tr>\n",
-       "      <th>constant</th>\n",
-       "      <td>93.838738</td>\n",
-       "      <td>1.50108</td>\n",
-       "      <td>91.745650</td>\n",
-       "      <td>96.883387</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>time</th>\n",
-       "      <td>0.382765</td>\n",
-       "      <td>0.12273</td>\n",
-       "      <td>0.141525</td>\n",
-       "      <td>0.560802</td>\n",
-       "    </tr>\n",
-       "  </tbody>\n",
-       "</table>\n",
-       "</div>"
-      ],
       "text/plain": [
-       "               mean      std   lower_ci   upper_ci\n",
-       "constant  93.838738  1.50108  91.745650  96.883387\n",
-       "time       0.382765  0.12273   0.141525   0.560802"
+       "(constant    91.344758\n",
+       " time         0.200908\n",
+       " dtype: float64,\n",
+       " constant    96.282862\n",
+       " time         0.610949\n",
+       " dtype: float64)"
       ]
      },
+     "execution_count": 6,
      "metadata": {},
-     "execution_count": 5
+     "output_type": "execute_result"
     }
    ],
-   "metadata": {}
+   "source": [
+    "results_without_cluster2 = bootstrap(\n",
+    "    data=df, outcome=ols_fit, n_draws=10_000, n_cores=2\n",
+    ")\n",
+    "\n",
+    "results_without_cluster2.ci(ci_method=\"bc\")"
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "## Doing a clustered bootstrap"
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "In the cluster robust variant of the bootstrap, the original dataset is divided into clusters according to the values of some user-specified variable, and then clusters are drawn uniformly with replacement in order to create the different bootstrap samples. \n",
     "\n",
     "In order to use the cluster robust boostrap, we simply specify which variable to cluster by. In the example we are working with, it seems sensible to cluster on individuals, i.e. on the column \"id\" of our dataset."
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
-   "source": [
-    "results_with_cluster = bootstrap(data=df, outcome=ols_fit, cluster_by=\"id\")\n",
-    "\n",
-    "results_with_cluster[\"summary\"]"
-   ],
+   "execution_count": 7,
+   "metadata": {},
    "outputs": [
     {
-     "output_type": "execute_result",
      "data": {
-      "text/html": [
-       "<div>\n",
-       "<style scoped>\n",
-       "    .dataframe tbody tr th:only-of-type {\n",
-       "        vertical-align: middle;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe tbody tr th {\n",
-       "        vertical-align: top;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe thead th {\n",
-       "        text-align: right;\n",
-       "    }\n",
-       "</style>\n",
-       "<table border=\"1\" class=\"dataframe\">\n",
-       "  <thead>\n",
-       "    <tr style=\"text-align: right;\">\n",
-       "      <th></th>\n",
-       "      <th>mean</th>\n",
-       "      <th>std</th>\n",
-       "      <th>lower_ci</th>\n",
-       "      <th>upper_ci</th>\n",
-       "    </tr>\n",
-       "  </thead>\n",
-       "  <tbody>\n",
-       "    <tr>\n",
-       "      <th>constant</th>\n",
-       "      <td>93.787766</td>\n",
-       "      <td>1.177882</td>\n",
-       "      <td>91.538014</td>\n",
-       "      <td>96.098402</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>time</th>\n",
-       "      <td>0.387924</td>\n",
-       "      <td>0.101298</td>\n",
-       "      <td>0.194408</td>\n",
-       "      <td>0.580074</td>\n",
-       "    </tr>\n",
-       "  </tbody>\n",
-       "</table>\n",
-       "</div>"
-      ],
       "text/plain": [
-       "               mean       std   lower_ci   upper_ci\n",
-       "constant  93.787766  1.177882  91.538014  96.098402\n",
-       "time       0.387924  0.101298   0.194408   0.580074"
+       "constant    1.176006\n",
+       "time        0.100272\n",
+       "dtype: float64"
       ]
      },
+     "execution_count": 7,
      "metadata": {},
-     "execution_count": 6
+     "output_type": "execute_result"
     }
    ],
-   "metadata": {}
+   "source": [
+    "results_with_cluster = bootstrap(data=df, outcome=ols_fit, cluster_by=\"id\")\n",
+    "\n",
+    "results_with_cluster.se()"
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
-    "We can see that the estimated standard errors are indeed of a smaller magnitude when we use the cluster robust bootstrap. \n",
+    "We can see that the estimated standard errors are indeed of smaller magnitude when we use the cluster robust bootstrap. \n",
     "\n",
-    "Finally, we can compare the results to full sample regressions using the statsmodels regression function and see that the cluster robust bootstrap gives standard error estimates very close to the cluster robust regression, while the regular bootstrap seems to overestimate the standard errors of both coefficients.\n",
+    "Finally, we can compare our bootstrap results to a regression on the full sample using statsmodels' OLS function.\n",
+    "We see that the cluster robust bootstrap yields standard error estimates very close to the ones of the cluster robust regression, while the regular bootstrap seems to overestimate the standard errors of both coefficients.\n",
     "\n",
-    "Note: We would not expect the asymptotic statsmodels standard errors to be exactly the same as the bootstrapped standard errors.\n"
-   ],
-   "metadata": {}
+    "**Note**: We would not expect the asymptotic statsmodels standard errors to be exactly the same as the bootstrapped standard errors.\n"
+   ]
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "                            OLS Regression Results                            \n",
+      "==============================================================================\n",
+      "Dep. Variable:                  pulse   R-squared:                       0.096\n",
+      "Model:                            OLS   Adj. R-squared:                  0.086\n",
+      "Method:                 Least Squares   F-statistic:                     13.75\n",
+      "Date:                Tue, 21 Jun 2022   Prob (F-statistic):           0.000879\n",
+      "Time:                        11:01:55   Log-Likelihood:                -365.51\n",
+      "No. Observations:                  90   AIC:                             735.0\n",
+      "Df Residuals:                      88   BIC:                             740.0\n",
+      "Df Model:                           1                                         \n",
+      "Covariance Type:              cluster                                         \n",
+      "==============================================================================\n",
+      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
+      "------------------------------------------------------------------------------\n",
+      "constant      93.7611      1.205     77.837      0.000      91.400      96.122\n",
+      "time           0.3873      0.104      3.708      0.000       0.183       0.592\n",
+      "==============================================================================\n",
+      "Omnibus:                       20.828   Durbin-Watson:                   0.827\n",
+      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               26.313\n",
+      "Skew:                           1.173   Prob(JB):                     1.93e-06\n",
+      "Kurtosis:                       4.231   Cond. No.                         31.7\n",
+      "==============================================================================\n",
+      "\n",
+      "Notes:\n",
+      "[1] Standard Errors are robust to cluster correlation (cluster)\n"
+     ]
+    }
+   ],
    "source": [
     "y = df[\"pulse\"]\n",
     "x = df[[\"constant\", \"time\"]]\n",
@@ -425,348 +366,173 @@
     "\n",
     "cluster_robust_ols = sm.OLS(y, x).fit(cov_type=\"cluster\", cov_kwds={\"groups\": df[\"id\"]})\n",
     "\n",
-    "cluster_robust_ols.summary()"
-   ],
-   "outputs": [
-    {
-     "output_type": "execute_result",
-     "data": {
-      "text/html": [
-       "<table class=\"simpletable\">\n",
-       "<caption>OLS Regression Results</caption>\n",
-       "<tr>\n",
-       "  <th>Dep. Variable:</th>          <td>pulse</td>      <th>  R-squared:         </th> <td>   0.096</td>\n",
-       "</tr>\n",
-       "<tr>\n",
-       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.086</td>\n",
-       "</tr>\n",
-       "<tr>\n",
-       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   13.75</td>\n",
-       "</tr>\n",
-       "<tr>\n",
-       "  <th>Date:</th>             <td>Mon, 06 Sep 2021</td> <th>  Prob (F-statistic):</th> <td>0.000879</td>\n",
-       "</tr>\n",
-       "<tr>\n",
-       "  <th>Time:</th>                 <td>11:58:21</td>     <th>  Log-Likelihood:    </th> <td> -365.51</td>\n",
-       "</tr>\n",
-       "<tr>\n",
-       "  <th>No. Observations:</th>      <td>    90</td>      <th>  AIC:               </th> <td>   735.0</td>\n",
-       "</tr>\n",
-       "<tr>\n",
-       "  <th>Df Residuals:</th>          <td>    88</td>      <th>  BIC:               </th> <td>   740.0</td>\n",
-       "</tr>\n",
-       "<tr>\n",
-       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
-       "</tr>\n",
-       "<tr>\n",
-       "  <th>Covariance Type:</th>       <td>cluster</td>     <th>                     </th>     <td> </td>   \n",
-       "</tr>\n",
-       "</table>\n",
-       "<table class=\"simpletable\">\n",
-       "<tr>\n",
-       "      <td></td>        <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
-       "</tr>\n",
-       "<tr>\n",
-       "  <th>constant</th> <td>   93.7611</td> <td>    1.205</td> <td>   77.837</td> <td> 0.000</td> <td>   91.400</td> <td>   96.122</td>\n",
-       "</tr>\n",
-       "<tr>\n",
-       "  <th>time</th>     <td>    0.3873</td> <td>    0.104</td> <td>    3.708</td> <td> 0.000</td> <td>    0.183</td> <td>    0.592</td>\n",
-       "</tr>\n",
-       "</table>\n",
-       "<table class=\"simpletable\">\n",
-       "<tr>\n",
-       "  <th>Omnibus:</th>       <td>20.828</td> <th>  Durbin-Watson:     </th> <td>   0.827</td>\n",
-       "</tr>\n",
-       "<tr>\n",
-       "  <th>Prob(Omnibus):</th> <td> 0.000</td> <th>  Jarque-Bera (JB):  </th> <td>  26.313</td>\n",
-       "</tr>\n",
-       "<tr>\n",
-       "  <th>Skew:</th>          <td> 1.173</td> <th>  Prob(JB):          </th> <td>1.93e-06</td>\n",
-       "</tr>\n",
-       "<tr>\n",
-       "  <th>Kurtosis:</th>      <td> 4.231</td> <th>  Cond. No.          </th> <td>    31.7</td>\n",
-       "</tr>\n",
-       "</table><br/><br/>Notes:<br/>[1] Standard Errors are robust to cluster correlation (cluster)"
-      ],
-      "text/plain": [
-       "<class 'statsmodels.iolib.summary.Summary'>\n",
-       "\"\"\"\n",
-       "                            OLS Regression Results                            \n",
-       "==============================================================================\n",
-       "Dep. Variable:                  pulse   R-squared:                       0.096\n",
-       "Model:                            OLS   Adj. R-squared:                  0.086\n",
-       "Method:                 Least Squares   F-statistic:                     13.75\n",
-       "Date:                Mon, 06 Sep 2021   Prob (F-statistic):           0.000879\n",
-       "Time:                        11:58:21   Log-Likelihood:                -365.51\n",
-       "No. Observations:                  90   AIC:                             735.0\n",
-       "Df Residuals:                      88   BIC:                             740.0\n",
-       "Df Model:                           1                                         \n",
-       "Covariance Type:              cluster                                         \n",
-       "==============================================================================\n",
-       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
-       "------------------------------------------------------------------------------\n",
-       "constant      93.7611      1.205     77.837      0.000      91.400      96.122\n",
-       "time           0.3873      0.104      3.708      0.000       0.183       0.592\n",
-       "==============================================================================\n",
-       "Omnibus:                       20.828   Durbin-Watson:                   0.827\n",
-       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               26.313\n",
-       "Skew:                           1.173   Prob(JB):                     1.93e-06\n",
-       "Kurtosis:                       4.231   Cond. No.                         31.7\n",
-       "==============================================================================\n",
-       "\n",
-       "Notes:\n",
-       "[1] Standard Errors are robust to cluster correlation (cluster)\n",
-       "\"\"\""
-      ]
-     },
-     "metadata": {},
-     "execution_count": 7
-    }
-   ],
-   "metadata": {}
+    "print(cluster_robust_ols.summary())"
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "## Splitting up the process"
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
-    "In many situations the above examples are enough. However, sometimes it is important to split up the bootstrapping process into smaller steps. Examples for such situations are:\n",
+    "In many situations, the above procedure is enough. However, sometimes it may be important to split the bootstrapping process up into smaller steps. Examples for such situations are:\n",
     "\n",
     "- You want to look at the bootstrap estimates\n",
     "- You want to do a bootstrap with a low number of draws first and add more draws later without duplicated calculations\n",
     "- You have more bootstrap outcomes than just the parameters\n",
     "\n",
-    "### Accessing bootsrap outcomes\n",
+    "### Accessing bootstrap outcomes\n",
     "\n",
-    "The bootstrap outcomes are stored in the dictionary you get when calling the bootstrap function. "
-   ],
-   "metadata": {}
+    "The bootstrap outcomes are stored in the results object you get back when calling the bootstrap function. "
+   ]
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
-   "source": [
-    "my_outcomes = bootstrap(data=df, outcome=ols_fit, seed=1234)[\"outcomes\"]\n",
-    "my_outcomes.head()"
-   ],
+   "execution_count": 9,
+   "metadata": {},
    "outputs": [
     {
-     "output_type": "execute_result",
      "data": {
-      "text/html": [
-       "<div>\n",
-       "<style scoped>\n",
-       "    .dataframe tbody tr th:only-of-type {\n",
-       "        vertical-align: middle;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe tbody tr th {\n",
-       "        vertical-align: top;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe thead th {\n",
-       "        text-align: right;\n",
-       "    }\n",
-       "</style>\n",
-       "<table border=\"1\" class=\"dataframe\">\n",
-       "  <thead>\n",
-       "    <tr style=\"text-align: right;\">\n",
-       "      <th></th>\n",
-       "      <th>constant</th>\n",
-       "      <th>time</th>\n",
-       "    </tr>\n",
-       "  </thead>\n",
-       "  <tbody>\n",
-       "    <tr>\n",
-       "      <th>0</th>\n",
-       "      <td>92.377004</td>\n",
-       "      <td>0.465963</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>1</th>\n",
-       "      <td>94.549027</td>\n",
-       "      <td>0.219700</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>2</th>\n",
-       "      <td>94.171886</td>\n",
-       "      <td>0.298749</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>3</th>\n",
-       "      <td>94.087851</td>\n",
-       "      <td>0.346593</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>4</th>\n",
-       "      <td>93.624816</td>\n",
-       "      <td>0.466219</td>\n",
-       "    </tr>\n",
-       "  </tbody>\n",
-       "</table>\n",
-       "</div>"
-      ],
       "text/plain": [
-       "    constant      time\n",
-       "0  92.377004  0.465963\n",
-       "1  94.549027  0.219700\n",
-       "2  94.171886  0.298749\n",
-       "3  94.087851  0.346593\n",
-       "4  93.624816  0.466219"
+       "[constant    92.377004\n",
+       " time         0.465963\n",
+       " dtype: float64,\n",
+       " constant    94.549027\n",
+       " time         0.219700\n",
+       " dtype: float64,\n",
+       " constant    94.171886\n",
+       " time         0.298749\n",
+       " dtype: float64,\n",
+       " constant    94.087851\n",
+       " time         0.346593\n",
+       " dtype: float64,\n",
+       " constant    93.624816\n",
+       " time         0.466219\n",
+       " dtype: float64]"
       ]
      },
+     "execution_count": 9,
      "metadata": {},
-     "execution_count": 8
+     "output_type": "execute_result"
     }
    ],
-   "metadata": {}
+   "source": [
+    "result = bootstrap(data=df, outcome=ols_fit, seed=1234)\n",
+    "my_outcomes = result.outcomes\n",
+    "\n",
+    "my_outcomes[:5]"
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "To further compare the cluster bootstrap to the uniform bootstrap, let's plot the sampling distribution of the parameters on time. We can again see that the standard error is smaller when we cluster on the subject id. "
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
    "source": [
-    "my_outcomes_clustered = bootstrap(data=df, outcome=ols_fit, seed=1234, cluster_by=\"id\")[\n",
-    "    \"outcomes\"\n",
-    "]\n",
-    "\n",
-    "# clustered distribution in blue\n",
-    "sns.histplot(my_outcomes_clustered[\"time\"], kde=True, stat=\"density\", linewidth=0)\n",
-    "\n",
-    "# non-clustered distribution in orange\n",
-    "sns.histplot(my_outcomes[\"time\"], kde=True, stat=\"density\", linewidth=0, color=\"orange\")"
-   ],
+    "result_clustered = bootstrap(data=df, outcome=ols_fit, seed=1234, cluster_by=\"id\")\n",
+    "my_outcomes_clustered = result_clustered.outcomes"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
    "outputs": [
     {
-     "output_type": "execute_result",
-     "data": {
-      "text/plain": [
-       "<AxesSubplot:xlabel='time', ylabel='Density'>"
-      ]
-     },
-     "metadata": {},
-     "execution_count": 9
-    },
-    {
-     "output_type": "display_data",
      "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABCN0lEQVR4nO3dd3zb5bX48c+RZHnvvRKP7A04CWGUvVdbaAudcGkpbentvLfj1/bSdXu7b4G2lNLFLYWWQim0hFH2SALZiTMdx/GIHe9tWZb0/P74KmASD3lo2ef9eukVS3r01Ylj5+hZ5xFjDEoppWYvW7gDUEopFV6aCJRSapbTRKCUUrOcJgKllJrlNBEopdQs5wh3ABOVlZVlSkpKwh2GUkpFlS1btrQaY7JHei7qEkFJSQmbN28OdxhKKRVVROTIaM/p0JBSSs1ymgiUUmqW00SglFKznCYCpZSa5TQRKKXULKeJQCmlZjlNBEopNcsFPRGIiF1EtonIP0Z4TkTkDhGpEpGdInJqsONRSin1dqHoEXwG2DvKc5cB8/23W4BfhiAepZRSwwR1Z7GIFAFXAN8FPj9Ck2uA+4x1Os5GEUkTkXxjTGMw41Iz00+fORBw289dtCCIkSgVXYLdI/hf4D8B3yjPFwJ1w+7X+x9TSikVIkFLBCJyJdBsjNkyVrMRHjvp7EwRuUVENovI5paWlmmLUSmlVHB7BGcCV4tIDfAgcL6I/PGENvVA8bD7RcDREy9kjLnHGFNhjKnIzh6xeJ5SSqlJCtocgTHmK8BXAETkXOCLxpgPntDsMeA2EXkQWAt06fyAGm4i4/5KqckJeRlqEbkVwBhzN/AEcDlQBfQDN4U6HjXz9A566Hd7sImQHOcg1mEPd0hKRbSQJAJjzAvAC/6v7x72uAE+FYoY1MzWO+hhe10nB4/10O3yvO25nORYFuYms6QghbgYTQpKnSjqDqZRajifMWyr7WRDdRs+n6EkK5FVxWkkx8XgM4aOPjc1bf28XNXKpsPtVJSkc8qctHCHrVRE0USgotagx8v63U0caeunLCuRs+dnkZbgPKnd2rJMWnoG2VjdxmuH2jh4rJd3ripkfm5yGKJWKvJorSEVlVxDXh7e0kBdez/nL8rhyhX5IyaB47KTY7lqZQFXrcind9DDlXe+wh9eq8EanVRqdtMegYo6bo+Pv21roL3fzdUrC5ibmRjwa8uyk/hAShz/2nuM/3qskgffqOX8RTk4bCN/JtIdyGo20B6BiirGGJ7e00RL7yCXL8+bUBI4LjHWwdUrCzi9NIO9jT08srWBfrdn/BcqNUNpIlBR5fWadg619PGO+dmUZSVN+joiwtqyTC5flkdLzyAPvlFHa+/gNEaqVPTQRKCiRmPXAJuq21mUl8zKotRpueb83GSuO60InzH8ZXMd1S2903JdpaKJJgIVFYa8Pp6qPEZynINzF2YjMlKZqsnJTYnj+oo5pCc4eXxnI2/UtOsksppVNBGoqLDpcDtdA0NctCQ3KDuFk+IcvOe0IhbkJvHaoTaerGzC4x2taK5SM4smAhXxWnsH2VbbwZL8FIrSE4L2Pg67jUuX5nFGeSYHjvXy0JZ6GrsGgvZ+SkUKTQQqohljeOlgC067jbPmZwX9/USE1SUZXLUin45+N1fc8QovHdDS52pm030EKqIdae+nrn2Ad8zPIj6EdYLKspO4fvUcNh1u4yO/e51PnzePz1y4ALtt9LkJPSFNRSvtEaiI5fMZXqlqJTU+hhVFaSF//4xEJ3//1Flce2oRdzxXxQfv3URzjyvkcSgVbJoIVMR6srKJtl43p5dljPlJPJjinXZ+9J6V/OC6FWyr6+CKO15hw6G2sMSiVLBoIlARyecz/OxfB8lIcLIgAorDvbeimEc/dSbJcQ4+cO9Gfv58FT6fLjFVM4POEaiwGG88/VBLL/uP9XDp0jxs07hnYCoW5aXw2G1n8ZVHdvHDp/bz+uF2fvq+VWQkjl7sTqlooD0CFZG21XaSHOdgfs7ky0gEQ1KsgzuuX8V33rmMDYfauOrOV3Q3sop6QUsEIhInIq+LyA4RqRSRb47Q5lwR6RKR7f7bN4IVj4oezT0uGjoHWFWUhi1McwNjERE+ePpc/vqJdbiGvLz3VxvY39QT7rCUmrRg9ggGgfONMSuBVcClInL6CO1eNsas8t++FcR4VJTYXtdJjF1YWpAS7lDGtKIojYduXYfdJnzoN5voGhgKd0hKTUrQEoGxHO8zx/hvOrumxtQ36OFAUy9L8lOIjYLzhcuyk/i/m9fi9vp4dFsDg0PecIek1IQFdY5AROwish1oBp4xxmwaodk6//DRehFZOsp1bhGRzSKyuaVFd3nOZLsauvAaw8ritHCHErAFucn8+sMVdLuGWF/ZpAXrVNQJaiIwxniNMauAImCNiCw7oclWYK5/+OhO4NFRrnOPMabCGFORnZ0dzJBVGHl9hp31XZRmJZI+xrGTkWh1SQbnLMjmSFs/22o7wx2OUhMSklVDxphO4AXg0hMe7z4+fGSMeQKIEZHgF5RREelwax8DQ15WFE7PWQOhtrwwlfLsRF471KaH3KioEsxVQ9kikub/Oh64ENh3Qps88ReWF5E1/nh02+Ystaexm8RYO3Myg1dhNJhEhPMX5eB02Hh2b7MOEamoEcweQT7wvIjsBN7AmiP4h4jcKiK3+ttcB+wWkR3AHcD1Rn97ZqW+QQ81bX0sykuJmA1kk5HgdHDW/Cyaul3saewOdzhKBSRoO4uNMTuBU0Z4/O5hX98F3BWsGFT02H+sB2NgcV74y0lM1eK8ZHY3dPFqVRvzc5JxOnTfpops+hOqws4Yw57GbnJTYslMig13OFMmIrxjfjYDQ16213WGOxylxqWJQIVdS+8gbb1uluRH9gayichLjaMsK5EttR24dG+BinCaCFTY7W3swS4SEVVGp9PpZZm4PT7tFaiIp9VHVVj5jGF/Uw+l2YnEReBO4omcOnai7ORYSrMS2VnfxWlz04mx6+cuFZn0J1OFVUPHAANDXhZEWJXR6XLanHQGhry6gkhFNE0EKqyqmntx2ISSrMRwhxIUBWlx5KXEsb22U/cVqIiliUCFjc8Yqlp6KclKnLHDJiLCyqJUOgeGqOsYCHc4So1oZv72qajQ2Omi3+2NuMNnptu8nCTiYmzsqu8KdyhKjUgTgQqbg8092G1CSebMHBY6zmG3sTQ/lUOtvfQNesIdjlIn0USgwsIcHxbKTJgVO2+XFqZgDOzTk8xUBJr5v4EqIjV2uegb9DJvhg8LHZee4CQvJY69jd06aawijiYCFRYHm3ux24TSGbpaaCSL85Np63PT0qMlqlVk0USgQs4YQ3VLL8Xp8cQ6Im8TWbAsyE3GLsJeHR5SEUYTgQq5g829dLs8lGXNjmGh4+JirLMWqpp7dXhIRRRNBCrknt3bDDCrhoWOm5eTRO+ghx26lFRFEE0EKuSe3XuMnORYkuJmX6mrsqxEbALrdzeGOxSl3qSJQIVUe5+brbUds7I3ANbwUHF6Aut3NenwkIoYwTyzOE5EXheRHSJSKSLfHKGNiMgdIlIlIjtF5NRgxaMiw/P7mvGZ2TksdNy8nCRq2/u1EJ2KGMHsmw8C5xtjekUkBnhFRNYbYzYOa3MZMN9/Wwv80v+nmgl23n7SQ89tKicnPomrbL/G1jf6SzcmfjpoYYVbeXYSLxxoYf2uJpYWpIY7HKWCemaxAXr9d2P8txP7wtcA9/nbbhSRNBHJN8boAOoM5PYKLx5N5aq5bdgi9Hz60/vunNLrA0lg8U47a0szeGJ3I1+4eAEiEfrNULNGUOcIRMQuItuBZuAZY8ymE5oUAnXD7tf7HzvxOreIyGYR2dzS0hK0eFVwvd6STK/HzgWFneEOZXxBHr+/bHk+1S19HGzuHb+xUkEW1GUbxhgvsEpE0oC/icgyY8zuYU1G+ih00m+gMeYe4B6AiooKnWGLUs83pOK0+Tgzr4ftteGO5u1SvbUscj3GEtcjJPjaicGFDztuSaLHlkebvZxO+1yMTM8GuEuW5PL1R3fzzJ5jM+6IThV9QrJ+zxjTKSIvAJcCwxNBPVA87H4RcDQUManQe7ExldNze4h3+IL6PhMZ3knwtTDHvZF03xEM0GfLpt1RzhAJCB7iTBdp3lqyvftxSyJ1MWtoti8GmVpn+v5NteQkx3L/xiO4PWN/Pz530YIpvZdS4wlaIhCRbGDInwTigQuB75/Q7DHgNhF5EGuSuEvnB2am+j4nVd3x3DAvMob27GaQOUMbyPPsxkOs/z/4JbhtI+x2Nj7SfEcoGtpCuft5cmx7qHJeiMuWPqUYSjITeaOmHdeQNyLPa1azRzB7BPnAH0TEjjUX8RdjzD9E5FYAY8zdwBPA5UAV0A/cFMR4VBi91GitjjknP/w7apO9R5nvfhqn6eOoYyX1MWvwSuzoLxAbnfZSOm0lZHoPUup+kRWuP1PlvJB2x7xJx1GSlcDrNe3Utvfr8JAKq2CuGtoJnDLC43cP+9oAnwpWDCpyvHg0lcLEQcpTXOELwhjyPDsoGXqVQUlhd+y19NrzAn+9CG2OBXTbCljoXs9C95McMWdwNGZy219yU+KIi7FR09qniUCF1ezb469CbsgnvNqUwtUlbYRrpaQYL6VDL5HrqaTNXsYh54V4xTmpaw3ZkqiMfRfz3M8yd+g17GaQupjTmehfzibC3IxEatr6McboMlIVNpoIVNBtbU2k12MP27CQw7hYMPgkqb566h2nTeo/7RMZcXDQeRFedwxFni14xcnRmNMmfJ2SrAT2H+vhWM8geSlxU4pJqcnSRKCC7sWjqTjExxl5oS+pEOfrYNHgP4g1PRx0XkirY9H0XVxsVDvPw+b2MHdoA0MSP+FLzM2wSm3UtPZpIlBho4lABd2Ljamcmt1Lckxwl42eKMVbz8LB9RiEPbHvpMdeMP1vIsIh5wXEDLoodz/PZV2fp8NROqFLPJ90Oq3NnVCWOf3xKRUArT6qgqp5wEFlRyLn5Ie2N5A7tJvFg4/hlkR2xb03OEnAz4id/bGX0mvLYYH7SRJ9zRN6/aq0Fqr7Uul3e4IUoVJj00SggurlEC8btRk3pe7nKRt6gS5bMbvjrmPQlhL09/WJk32xVzIkCSwcXI/DDAT82pVpLRiEuvbAX6PUdNJEoILqxcZUsuKGWJLeH/T3SvC1cF3njeR5Kql3nMa+2CsmvTJoMjwSz/7Yy4gx/cwffApMYENhpYndJNrd1LYH/3uk1Eg0Eaig8frg5aYU3pHfFfRqo+WDz/Ch9qvI8VRywHkxdc51Uy4DMRl9thyqneeS5qtnztCGgF5jE1iaam0s08NqVDhoIlBBs7sjkY7BmKAOCzl9PVzc/RWu7r6NHnsBf0p/mDZHeGvztDgW0+RYTqFnG+me6oBesyy1jd5BD50DQ0GOTqmT6aohFTQvHk1BMJwdhGWjNuNmxcCDrO3/BXGmi40Jn2BTwifxidM6EinMamLOIsnXRLn7WXbackauYTTMitRWAGrb+0lPCN1wllKgiUBNo58+c+Bt99cfTqQssYv9R49N23s4zAALXf9kTf/dpPnqqI05nZcT/5PmmKXT9h7TwYidg85LWOH6M/PcT7Mn9p1jDlXlxA6QEuegrr2flUVpIYtTKdBEoILENeTlYG8a7yo8NPWLGUOWdz9LXQ+zxPUocaabZvsiHkn9NUdizp7yLuFgcdnSqHaew3z3vyj0bKEhZvWobUWgOCOBg829+HwGW6Qe4aZmJE0EKijq2vsxCCvTWid3AeMj37ODeYNPM2/wGdJ8dXiJ4WDsxeyMu4GGmIqITQDDtdoXkmavpXjodbptRfTY80dtOycjgcqj3RzrcZGfOvFdykpNliYCFRQ1bf0k2t2UJ01gotgYknzHyPIeYLnrzyT5WvASQ61zHW84b+FQ7IUM2DKCF3QwiHDYeS7Jribmu59mR9z78MrIpSSK0xMAqGsf0ESgQkoTgZp2xhiOtPexPLUNu4y/HNJm3OR49pDv2Umc6caHnUPO8zkYewmHnefitkV3iWavODnovISlgw9T7n6eA85LR+zNxDvt5CTHUtvez5rSKEt4KqppIlDTrrXXTd+glxWFYw8L2Y2bfM928oZ2EMMg3bYC6mNW024v49WkL4Yo2tDotedSF7OWuUMbyPHupdmxZMR2xRkJbKvtwO3x4XTo6m4VGkH7SRORYhF5XkT2ikiliHxmhDbnikiXiGz3374RrHhU6NS09QGMOT+Q6TnIKtcfKR56nR57Abtir6Uy7t20OBaPfVpYFDvqOJUuWxEl7peI83WM2KY4PR6fgaOdWm5ChU4wewQe4AvGmK0ikgxsEZFnjDF7Tmj3sjHmyiDGoULsSFs/2cmxZDhPXtBvN25K3S+Q7T1Ary2Hfc4r6LPnhiHKMBChynkhK1wPMH/waXbHXYeRt59VXJgWj90m1Lb3U5KVGKZA1WwTtB6BMabRGLPV/3UPsBcoDNb7qcgwOOTlaNcAJZkJJz0X5+tkmeshsrwHqY1Zy67Y62ZPEvBz25I45LyAJNNC8dDGk5532G0UpMZR26F1h1TohGQQUkRKsM4v3jTC0+tEZIeIrBeRyNoVpCbMqpcDJZlv/zSb6D3GMtdDxJh+9sReY62pD0MtoEjQ4SijybGMQs82Ur11Jz0/JyOBtl43fYNallqFRtB/E0UkCXgY+Kwx5sRaA1uBucaYlcCdwKOjXOMWEdksIptbWlqCGq+ampq2fmIdtredtpXsPcrSwUfxSiy74t5Lt70ojBFGhiMxZ9Iv6cxzP3NSyeriDP8yUu0VqBAJaiIQkRisJHC/MeaRE583xnQbY3r9Xz8BxIhI1gjt7jHGVBhjKrKzs4MZspoCYww1bX3MzUh4c2dskreJxYOPMyhJVMa+m0FbapijjAw+ieFg7CU4jIty93MwrOpodnIscQ6blqVWIRPMVUMC/AbYa4z5ySht8vztEJE1/njaghWTCq6WnkH63d43JznjfB0sHnycIUlgT+w14xZem236bVkciTmDDO9hcj2733zcJkJxRgJ17QNallqFRECJQEQeFpErRCY0qHsm8CHg/GHLQy8XkVtF5FZ/m+uA3SKyA7gDuN7oT37UqmmzPsHOzUwg1tfJosF/YLCxJ/YahjQJjKjJsZIO2xxKhl4h03PwzcfnZCTQO+iho1/LUqvgC3T56C+Bm4A7ROQh4PfGmH1jvcAY8wowZjEYY8xdwF0BxqAiXE1bHznJsSTFGK7s+iyxpoc9se8KyVGRUUuEQ7EXsmLgAS7r/gIPpD+EV2LfnCfQ4SEVCgElAmPMv4B/iUgqcAPwjIjUAb8G/miM0Y8ts1xnv5umLherSzM4t/e7zBnaQJXzgjGLrI3l9L47pznCyDUkCRyKvYDFg//gnN7v8Vzy7aTGx5AaH6OJQIVEwEM9IpIJ3Ah8FNgG/Aw4FXgmKJGpqPLSwVYM8O60Z1npeoA34j9Ki2NxuMOKGp32EjbH38xK1wMscj0GQHFGPA0dAwx5Azv7WKnJCnSO4BHgZSABuMoYc7Ux5s/GmE8DOvireGFfM4sTj/Ju7/eoi1nDq4mfD3dIUeeVxM9TH7OaC3u+QabnIHMyEnB7feyo6wx3aGqGC7RHcK8xZokx5nvGmEYAEasgjDGmImjRqajg8xk2HGzgrjk/wCtO1if/6KTSCWp8Rhw8kfwT3JLIld2fpjzNWjfx8sFJnumgVIACTQTfGeGxDdMZiIpeOxu6+HjK3ZTHVPFU8vdnXdmI6dRnz+GJlJ+Q5j3Cla5vkJvi5NUqTQQquMZMBP51/qcB8SJyioic6r+dizVMpBTV2//MjVn/4PXYD3M49txwhxP16p1reTXxcywcXM/H855gW10nPS5dj6GCZ7xVQ5dgTRAXAcM3hfUAXw1STCqaDLZxTudXqaWMjcn/Ee5oZozN8R+lYGgbH+EunorPZWP1aVy0RHtaKjjG7BEYY/5gjDkPuNEYc96w29UjlYxQs0/fq58gxdbF1sKf4RVnuMOZOcTGk8k/oNNWzN1zv0flge3hjkjNYOMNDX3Q/2WJiHz+xFsI4lORrPYhEpse4mfHbuC0VeeHO5oZx21L5rHUX+K0G67u+iQM9YY7JDVDjTdZfLyWcBKQPMJNzVYDx+CNT1DlXcxztpve3Amrpleno5QXs+9krqMG10sfAKN7CtT0G3OOwBjzK/+f3wxNOCpqbPl3zFAPnzj0HS47Q88bCqbS5e/iu/e/yjfk17Drm7BCfx3V9Ap0Q9kPRCRFRGJE5FkRaR02bKRmm/rHoPYv7E7/DAddxVy0JC/cEc1oi/KSecz1HjZwNez+FtQ+FO6Q1AwT6D6Ci/2HylwJ1AMLAF0iMhu5u+CNT0Dacn5+7F3kp8axrFCLygWTiHDuwhw+XfUxTOYZsOEj0L4l3GGpGSTQRBDj//Ny4AFjTHuQ4lGRbvuXwNWE69R7eOFgJxctycV/pIQKovMW5tDaL2wv/Q3EZsOLV0N/Q7jDUjNEoIngcRHZB1QAz4pINuAKXlgqIjW/BFW/goWf5bmWubiGfFyyVIeFQuHsBVnYbcIz1cA5j8NQN7x0DXi0OqmauoASgTHmy8A6oMJfcroPuCaYgakI43XBpo9BYims+Bb/3NlIVpKTtaUZ4Y5sVkiJi6FibjrP7WuG9BVwxp+gfas1TKQridQUTeTEscXA+0Tkw1gni10cnJBURNr1Leg5AGvvod8Xy3P7mrl0WR4Oe1CPvVbDnL8oh31NPTR2DUDRVXDKD6Dur9ZKIqWmINBVQ/8H/Ag4C1jtv41ZdVREikXkeRHZKyKVIvKZEdqIiNwhIlUislNETp3E30EFW8d22PsDKLsR8i7k+X0tDAx5uWJ5Qbgjm1XOW5QDwPP7WqwHFn0Bym6yVhLVPBDGyFS0C/SoygpgyQTPE/YAXzDGbBWRZGCLiDxjjNkzrM1lwHz/bS3WkZhrJ/AeKth8Hth4M8RmwSk/BuCfu46SlRTLGh0WCqn5OUkUpsXz/P5m3r92DojA6ruh9xBsvAmSyiBLf33UxAXar98NTGhW0BjTaIzZ6v+6B9gLnLjz6BrgPmPZCKSJyOTONlTBsf8O6NgKp90BsRn0DXp4bl8zly/Pw27T1UKhJCKcvyiHV6taGfR4rQftTjjrYYgvsCaP++rCG6SKSoEmgixgj4g8JSKPHb8F+iYiUgKcAmw64alCYPhPbj0nJwtE5BYR2Swim1taWgJ9WzVVfXWw6xtQcAXMeQ8Az+1rxjXk44rlmq/D4bxF2fS7vWyqHraCOy4Lzv0HeAfgxau0JpGasECHhm6f7BuISBLwMPBZ/6a0tz09wktOGn4yxtwD3ANQUVExkeEpNRVbP2utSKm40xqGAP65s5Gc5FgqSnRYKBzOKM8iPsbO03uaeMeC7LeeOPIXKLgKjvwJnlwNc9775r/ZuFbcHpRYVfQIdPnoi0ANEOP/+g1g63ivE5EYrCRw/yhlq+uB4mH3i4CjgcSkgqzhCah7BJZ9HZJKAeh2DfH8/mYuX56vw0JhEhdj5/xFOTy5+xhe3wmfiZLnQf4l0LMPjj0XngBVVAqoRyAiHwNuATKAcqzhm7uBC8Z4jQC/AfYaY34ySrPHgNtE5EGsSeKu42ciqzDy9MPm2yBlMT9ruArf0QMA7D7axaDHx4Dby0+fORDmIGee0/vuHPmJnZlvu3tpSgb/7C1n83M/Zm3OCcNAGWvA1QKtr1gT/OkrgxStmkkCHRr6FLAG/xi/MeagiOSM85ozgQ8Bu0Rku/+xrwJz/Ne4G3gCq2xFFdAP3DSR4FWQVH4X+g7DBS+wZsOv3nz46fo15MfFco3910hf+MKb7c4r6CTW7mN9bcbJiUAECi4DdzscfRxiMyCheOQLKeUX6GTxoDHGffyOiDgYYSx/OGPMK8YYMcasMMas8t+eMMbc7U8C+FcLfcoYU26MWW6M2Tz5v4qaFl17Ye8PofTDkHvOmw+3DsaxtyeDs7KOBjz0rIIjKcbHOfldPFmXzomjQwCIHYrfA44Uq1KpRyeP1dgCTQQvishXsQ6xvwh4CHg8eGGpsDAG3vgkOJLglB++7alXWq3NY2dl6RROJLisuIOmASfb2xJHbuCItyaMvQNQ97CWoVBjCjQRfBloAXYBH8ca0vlasIJSYVLzR2h+AVZ9H+LeGvkzxkoEC5PbyYkbCF986k0XFHYSY7OGh0YVnwcFV0JfDRx7NmSxqegT0ByBMcYnIo8CjxpjdCF/NNl5e2DtPANw8C6IL4K++re97nBfCg0DSdxcWhmEANVkpDi9nJ3XzRN16Xz1lLrRh+vSV8JAPbS+BvGFkLokpHGq6DDe4fUiIreLSCuwD9gvIi0i8o3QhKdCpvk5axih4IqT1p+/3FqAQ3ycnqkLuiLJpcUdNPTFjj48dFzepVaCb/i7taJIqROMNzT0WazVP6uNMZnGmAysZZ5nisjngh2cChHXMevEq4wKazhhGI9PeLW1gFPTm0lyeMIUoBrJpcUdxNm9/LU6a+yGNru1M1wcUPcX8A6GJkAVNcZLBB8GbjDGHD7+gDGmGvig/zkV7YyBxifBHgc555309JaOHHo8Ts7Lrg9DcGosKU4vlxV38NiRDFyecZZyxaRA8XUw2AYNj1n/7kr5jZcIYowxrSc+6J8niBmhvYo23fusycSc86yVJid4vrmIDOcAK9JO+jFQEeA9Za30DDl4qj59/MZJpZB7AXTvgbaNwQ9ORY3xEoF7ks+paODzQNPTEJsDGaed9PTRPic7u7I4J7sBrSgRmU7P7aEocZCHxhseOi7rDEheBE3/gn7t5SnLeIlgpYh0j3DrAZaHIkAVRK2vwVCnVZ9GTv5R+Gt1JgbhnGw9JD1S2QSuLW3l1aYUGvqc479ABIquhphk63Qzry4HVuMkAmOM3RiTMsIt2RijQ0PRbKgbWl6BlEXWgSYn8Bl4qDqLpSlt5OregYh2XVkrBuHh6szxGwPY4635gqEeqNf5AjWxM4vVTNL0LOCDvJGPnt54LJm6vjjOy9Hhg0hXnOTmjNxuHqrOGrnkxEgSiiDvQqtS6YFRit2pWUMTwWzU3wBdOyFzHThHnmT8U1U2KTEeVmccC3FwajKuL2+hri+OF46mBv6izNMheQFs+yK0aZmv2UwTwWxjDBz7F9gTIPusEZu0DFirUK4ra8Vp0xo10eCyOR3kxrv53f7cwF8kAoXvhLg8eOW94O4MVngqwmkimG16D/mXi74D7LEjNvlLdTZDPhsfmK+7UKNFjM3wofnNvNyUysGuuMBf6IiHMx+E/lrY9FGdL5ilAj2PQM0Ex3sDMWmQXjFiE6/PGhZal9tNeYqLZt0+EHYbqtsCanfDPAd3VeZzz948fnh6TeBvkH0GrPxv2P4lOPgLWPCpyQWqopb2CGaTrl1WOYnc862yAyN4sTGVhr5YPji/OcTBqanKjPNwfXkrfzucGdhS0uEWfxEKLoetn4f2cU+hVTOMJoLZwuexzrGNy4fUZaM2u/9gDtlxbi4u6gxdbGrafGxxEwD37M0bp+UJxAan/wFis635gqHuIESnIlXQEoGI/FZEmkVk9yjPnysiXSKy3X/TiqbB1P4GDHVZSwZHqVlc3+fkuaOpXF/eSoxNx4qjUWGim+vK2nigKpu63gn2CuKyrPmCvsPWAUU6XzBrBHOO4PfAXcB9Y7R52RhzZRBjUABeF7S8DEnlI24eO+6BqmxE4Pp5OkkczT6zvIFHazL40Y4ifnZm9fgv2Hn72+9nvwNq7oehXkhfNfZrV9w+9vMqKgStR2CMeQloD9b11QS0vmqVEsi9YNQmbq/w50PZnF/QSWGilpGKZvkJQ9y86Bh/P5LJzraEiV8g+2xImAuNT1jVStWMF+45gnUiskNE1ovI0tEaicgtIrJZRDa3tOin1Qnx9EHbJkhdCvH5ozZ7uj6NVleMLhmdIW5d0khG7BD/va144iM8YoPid/vPL/irNb+kZrRwJoKtwFxjzErgTuDR0RoaY+4xxlQYYyqys7NDFd/M0Poa+IYg+5wxm913MIeixEHOye8KUWAqmJJjfHxm2VE2Nqfw/ER2Gx8XkwKF14CryVpyrGa0sCUCY0y3MabX//UTQIyIBFhLVwXE1Qxtb0DqcogbPYFWdsTzenMKH1lwTMtNzyA3zGuhLHmA27fMod8ziV/1lIWQucbqUXYfmP4AVcQIWyIQkTwRa/mKiKzxx6IDktNpzw/AeCBn7N7A7/bnkuDw8t5y3T02kzjthu+traG2N44f7yic3EVyL7JKUDQ8qktKZ7CgrRoSkQeAc4EsEakH/gv/qWbGmLuB64BPiIgHGACuN0bXq02bgSZrl2jaCogdvTxxq8vBYzWZvK+8hVSnN4QBqlBYm9PLB+c389v9uVwxt51Ts/qAwHcrA6wrvhYO3QP1f4OSD414doWKbkFLBMaYG8Z5/i6s5aUqGPZ8H3xuayngGO4/mIPbZ+PGhbqTeKb60qo6nmtI5UubSvnHpZXE2if4eSs2C/Ivh4a/W2dY5Iz9M6Wij9YamokGGqHqbij9EMRmnPT08U+DHp/wu33LWZXWQnNrg9YVmqGSY3x8d80RbnphAT/ZWchXTpnEGRNpK6G3GppfgMQSSJwz3WGqMNI+3kxU+T/WSqFlXx+z2Ya2PDqH4rg070iIAlPhcl5BFx+Y18yv9ubzwtGUiV9ABAquAGca1D+iR1zOMJoIZpr+Bqj6FZTdOOYuYmNgfVMJhfG9rEjVrsBs8PVTa1mU1s8XNpTR7h65BPmY7LFQdK11xGXD41qCYgbRRDDTVH4PjBeWfm3MZgd60zjcl8oleUdGKz2kZpg4h+GuMw/R77Hx86oVgR9rOVxCIeRdAN17oWPLtMeowkMTwUzSVweHfg3l/wZJJWM2Xd84l0T7EGdnHQ1NbCoizEt18c2KWvZ0Z/JI/bzJXSRznVW3qvEp6ByxpqSKMpoIZpLK/wYMLP1/YzY75orn9fY8zs+pI86uS0Znm/eUtXJOdj2PNJSzo3MSezhFoOid1lDRq+8DT/+0x6hCSxPBTNF3BKp/A2U3j7ui45+NJdjFcGm+ThLPRiJwU8keihJ6uatqBW2DEzja8jhHEhS9C7r2wNbPTX+QKqQ0EcwUld+z/lz61TGbtbocvNBcxNlZDWQ4B0MQmIpEsXYfn5u/Da+x8b8HV+HxTWKiKKkclnwJqu6B2oenP0gVMpoIZoK+I1D9Wyj/KCQWj9n0vgM5eIyNKwpqQhObilj58f18vGwXVb1p3F+7cHIXWfFtyKiA12+xVqypqKSJYCao/G9AYMlXxmzWN2TjDwdyOS29mcL4vtDEpiLa2sxjXJZXw5NNJWxom+DxlgC2GDjjfuvwo403gvFNe4wq+HRncbTrOwKHfgvzbhm3N/DgoWy63A6unh/AqVUq6kykftBw75+zn6reVO45tIy5CT0UTPRDQsoCOO2n8PrHYf/PYJHOGUQb7RFEu93ftYqALR27N+DyCr/em8eanG7mJ+uZA+otDpvhM/N3EGPz8b8HVuH2TeK/hfKPQeHVsP3L0Llr+oNUQaWJIJr11kD176xfwoSiMZv+5VA2TQNOPrNM9w2ok2XGuvhE+U7qBpJ5sHbBxC8gAmvvBWc6vPp+a6hIRQ1NBNGs8nhv4MtjNhv0Cr+ozGd1dg9n5PaEKDgVbU5Jb+Xi3COsbyphZ+fopctHFZcNp/8OunbD9rFXr6nIookgWvUehurfW3MD4/YGsqzewPKjWk5CjekDc/dTGN/LLw8tp3soZuIXKLgM5n8K9v8UGp+Z/gBVUGgiiFbHewNLAugN7CmgIruHM3P1hCk1NqfNx23zdtDjcfLbw0smd5FTfggpi2HjR2BQDx2MBsE8oey3wJVAszFm2QjPC/Az4HKgH7jRGLM1WPHMKL3VUP0HmH+rVQRsDH+pzqKx38kP1h7W3oAKSEliD9cWVfGXugW80d7IurIxypDsvH3kx7PPhup74Zl3QPF1jPrDt2KU16uQCmaP4PfApWM8fxkw33+7BfhlEGOZWXZ/F8Q+bm+g3+3hjl0FrM7u4aw87Q2owF2Vf5i5Cd389vASutz2iV8gPh9yzoPuPdClq4giXdASgTHmJaB9jCbXAPcZy0YgTUTygxXPjNG9Hw7/wT83MHZv4LevHKbF5eTLq+q1N6AmxGEzfLx8N91DTr6zdez9KaPKOgMSiuHoE+DWJcuRLJwbygqBumH36/2PNYYnnAi28/a3vq59yJobsDlH75YD7S4Hv3p+ORcV9XBadm+QA1QzUWliN1cW1PBQdRlXz23n7PwJ9irFZhWmq7obGh6Fkg+PPkSkwiqck8Uj/USMeFSGiNwiIptFZHNLS0uQw4pgA0etrnbmOqv64xh+XplPn8fOf66cxPm0SvldW1RFabKL/9o8h0HvJP4Td6ZD3iXQVwNtG6c9PjU9wpkI6oHhfc4iYMTdTsaYe4wxFcaYiuzs7JAEF5GOPQf2eMhaN2az+j4n/3cwh/eUtTI/VTf2qMlz2nz812lHqO6J5959k6hFBJB+CiQvgGPPgqt5egNU0yKcieAx4MNiOR3oMsbosNBoeg9D7yFrNYZ97Prx/7OtCJsYPrtcq0GqqTu3oJtLijq4a3c+DX3OiV9ABAqvAlss1P8NfHoYUqQJWiIQkQeADcBCEakXkZtF5FYRudXf5AmgGqgCfg18MlixRD1jrE9TjhTIWD1m09ebk/hHbSa3LmkiP2EoRAGqme7rp9ViYPITx44kKxm4mqD5hekMTU2DoE0WG2NuGOd5A3wqWO8/o3TthoEGKLwGbKP/k3l98M0tcyhIGOTji5tCGKCayY5XNb26IJ6/1C3gV1sPsSJt5I1i68rGKE2RsgjSVkHrq5A8f9yT9FTo6M7iSOfpg6Z/QVw+pK0cs+lD1VlUdiTylVPqiXdoXXg1va7MP0xeXB+/r1nC0GRONAPIvxRiUqH+UfDqCXmRQhNBpNv7I/B0Q/4lYy6963Lb+eGOIlZn93DlnLG2byg1OTE2w0dK9tLoSuSfjaWTu4g91jr4fqgDmp6e1vjU5GkiiGT99bDn+5CyBBLnjtn0RzsK6XA7+K/TanWptgqaVWmtVKQf49GGMtrdsZO7SOJcyDoTOrZC/ePTG6CaFD2hLJJt/woer5ed7goGxzh9qqo3lT8ezOGmhcdYltEfwgDVbPTBufv4jx1n8afahdw2b+fkLpJzLvRWwesfhazdVglrFTbaI4hUx16Emj+yJeHfGLSljNrMa4R7q5eS7hzk8yt0uagKvty4Aa7Ir+HV1gIO9KRN7iI2h7Xr2N1pHXxvRtxLqkJEewSRyDsIb3wcEkt5Pf5WKvrvHbXpU01zONKfwmfnb2NX3Szeda1C6prCal5sKeT3NYv5zrIN2CYzHBmXCyu/C9v+wzpbo/ym6Q5TBUh7BJGo8ntWcbnVv8Qj8aM2axuM5aG6+axKa2ZNxrEQBqhmuzi7lw/M3c/hvlRebBm7+OGYFn4Ocs6BLZ+xNk2qsNBEEGm69sGe78HcG6DgkjGb3ndkMT6Em0r26gSxCrkzMhtZmNzBg7UL6PNMcnDBZod1f7C+3vAR3XUcJpoIIonxWUNC9gQ49adjNt3Skc3r7Xm8u7CKnLiBEAWo1FtE4CMle+nxOHmkft7kL5Q4FyruhJaXYd9Ppi9AFTBNBJFk/8+g+SU49UcQnztqs36Pg98eXkpxQg9X5NeELj6lTlCa2M15OfU8dWwODQOJU7jQh6H43bDza9AxyZVIatI0EUSKzl2w/ctQeDWU/duYTf9Uu5AOdywfL9uFw6arLVR4va/4ALE2L/fVLJr84h8RWP0rq2z1hg/qruMQ00QQCbwueO0D1i/B2nvH3EFc2ZXBs83FXJFfQ3mSHj+pwi8lZohri6rY2ZXNsw2pk79QXBas/Y31oWjn16cvQDUuTQSRYMf/s3741/52zI01g14b91QvIy+uj+uKDoYwQKXGdnFuLYXxvXx76yQPsDmu8ArrGNa9P7KGSVVIaCIIt/rHrAmy+Z+EwsvHbPqX+vk0DybwsbLdxNq1qJyKHA6b4cNz93KkN27yB9gcd8qPIakcXn0/DOiy6FDQDWWhsvP2kx8bbIVDv4b4AohJH7HN6X1WaYmqnlTWN5ZwYU4tS1I6ghmpUpOyIq2NS4o6uHN3PlfPbaM4yT25C8UkwdkPwdNnwCvvgQueBVvM9Aar3kZ7BOHiHYTaP4M4oPi9Y54z4PEJv6peRrrTxQ1z9ocwSKUm5vaKIzgE/t8bJVOrGpG+ypova3kZtn5husJTo9BEEA7GWPXYB9ug+D3gHHuC7dGGcuoHkrm5dA8JDt1woyJXfsIQX1xZz0uNqTx2JGNqFyt5v7Xz+MCdUH3f9ASoRhTURCAil4rIfhGpEpEvj/D8uSLSJSLb/bdvBDOeiGAMND0FPfsg72JIKhmz+eG+FB49WsaZWUc5NV1rCanI96H5zazM7OVbW+bQOWif2sVO+QHknmdVKW16bnoCVCcJ5pnFduDnwGXAEuAGEVkyQtOXjTGr/LdvBSueiNH6GrRtgsy11m0MLq/w86rlpDjc3FiyN0QBKjU1dht8b00NnW4H39wyxeMobQ44+2FIXgAvvRM6dkxLjOrtgtkjWANUGWOqjTFu4EHgmiC+X+Tr2A7H/gWpyyBv7BPHAH64vYiGgWQ+Xr6bJIceRK+ix5L0AW5bepS/1WSxvjZ9ahdzpsN5T1pDqC9cBr010xKjekswE0EhUDfsfr3/sROtE5EdIrJeRJaOdCERuUVENovI5paWKB0e6doDDY9BYpl1CP04SeC1Y8n8Zn8eF+UeYWVaa4iCVGr63LaskRUZfXz1jbk0D0xxgWJCEZz7JHgG4LmLoK92eoJUQHCXj470P92J6wi2AnONMb0icjnwKDD/pBcZcw9wD0BFRUX01VSo+RPU/dX6YZ4z9gohgG63nS9uKKU02cX75xwIUZBKTa8Ym+En66q54smlfHlTKb855+DJn3923j6xixZfCzV/hCdWwKWbIXkKxe7Um4LZI6gHiofdLwKODm9gjOk2xvT6v34CiBGRrCDGFHrVf4DXPmhVWJz7Qevw7nHcvmUOxwac/GRdNXF2XSWkote8VBdfWlnPc0fT+P2BnKlfMKEISj8CPjc8czZ0Vk79miqoieANYL6IlIqIE7geeGx4AxHJE7E+I4jIGn88ox/OG02Mgb0/gY03Qd6FMPf9YHeO+7K/12TwyOEsPrX0KKdk9YUgUKWC68aFx7iwsIPvbi1mS0vS1C8Ynw+lN1pfP3MWHF0/9WvOckFLBMYYD3Ab8BSwF/iLMaZSRG4VkVv9za4DdovIDuAO4HpjZsDhpT4PbL4Ntn3B6sqe81hAOyOruuL4yuslrM7u4dPLGkMQqFLBZxP48brDFCa6+eQr5bRMdb4AIC4HLn4VEufAC1fArm9b53moSQnqPgJjzBPGmAXGmHJjzHf9j91tjLnb//VdxpilxpiVxpjTjTGvBTOekBjqsZa5HfwFLP5POOvPYI8b92X9HhuffKWceLuPO888RIyWl1YzSKrTyy/PrqLLbefTr5Yz5JuGI/WSyuDiDVDyAdj1DXjxKuivn/p1ZyHdWTydOivhqdXQ+CSsvhtO+T7I+N9iY+Brb8zlYFc8/3tGNXkJulRUzTxL0gf47zVH2Nicwpc2TbEExXGOBFh3H1TcBceeh38sgf136pGXE6SJYLocvh+eWgPuTjj/WZj/8YBfevfePB45nMVnlh/l7Hw9Y0DNXO8ubePzy+t55HAWP9hRND0XFYEFn4IrdkPWOtjy7/D06XD0KaYn28x8mgimaqgHNt5snaqUcRpctg1yzwn45f+sTef724u5ck4b/77s6PgvUCrKfXpZI++f18wv9+Tzu/3TsJLouKQya+PZGfeD6xi8cCk8cyYcfVLnD8ahiWAqml+GJ1bC4d/Dkq9Y5XLj8wN++ZaWRD73WhmnZfXwo3WHsU3DsKlSkU4Evl1xhIuLOvjmlrncu2/087kndfGS98NVB2H1L605gxcug8cXwp4fgitKN6QGmZ5HMBlDPbDzG9Zh80mlcOHLkH3GhC6xsy2Bm15YQH6Cm3veUUWcXbuwKvptqA5s9fe6skzuOvMQn32tjO9snUPfkJ1/X3Z0vA33gbPHwvxboewmqP0rVP0Ktv8n7Pgq5J4Pc66Fondaq48UEm2rNSsqKszmzZvD8+Y7b4fufda6ZU83ZFRA7kUB7Q8YbltrIh9+fgGpTi8PXLBvzAM8Av3FUiqarCvLBMDjgy+/XsJfq7P5yIJjfO3UuuCtmHO1QOcO6N4L7nZrIUf2WVbJl/xLIHXJuKVfopmIbDHGVIz0nPYIAtW1B2ruh94qiMuFOe+xdjlO0KbmJG5+YQGZcUP86YL9FCZO8hQnpWYAhw1+sLaGdKeXX+/LY29HAnedVUVOvGf63ywu29rcmXuB1SOoe9i6bfuCdUsoskrD518CeRdAbOb0xxChtEcwHlcL7Lrd6lqKA3LeYZWPlonVWTcG/u9gDt/aUsycpEHuv2A/+QEsE9UegZqJjvcIhvt7TQZf3lRCstPLT9Yd5qy8EK2gc3dB7yH/rRp8Luvx+ELr7OSkMogvApv/d37F7aGJa5ppj2AyXM2w90dw4OfgG4R5t4I9HhyJE77UgMfG196Yy8OHs7igoJOfnFFNqlPXOSs13DUl7SxMG+CTL5fzwecWcm1pK//vlDoy4oLQOxjOmQoZp1o344OBo1bPv+eQdVRmy0sgMdYu5sRSaN9qHaUZwB6haKE9ghO9LQG4YO4NsOzrkLJw4pUSgecbUvn65rnU98Xy2eUN/PuyoxNaHaQ9AjUTjdQjOM7lEe6sLOBXe/JIcXr55NJGPjCvhXhHGJaAel3QVwO9h6HvMAz6Vx05M6yT0/IusIaakudH/PzCWD2C2ZUIdt4++nOuFmjfZJ2AZLzW4TE574DYyRVDreqK48c7C1lfl8G8lAG+u6aGtTm9gP7nrlQgavuTuK9mMZXdmWTFDXHL4iauK20Nfg9hLEM9Vs/h2LPQ9C/o9x+5klBkJYTjiSGhIHwxjkITwXE7b3/7fWOsLmDbJmt8UOyQuhyyz5x0AtjdnsAvKvNZX5dOrN3Hp5c18rFFTTiHLQ/VRKBU4BxJc/nZ7gJeaUolxubjgsJO3lnSzhm53aSEc4jVGGv1Ud/ht3oM3gHrOWemNZSU4L8509/qMYRpjkHnCE401G0dG9m5w/qHdCRBznnWzuBJzAG0uRw8fiSDh6qzqOxIJDnGwyeXNvJvC4+RGc5PL0rNAKtzevnj+QfY2xHPw4ezeLQmkyfrMrCL4ZSsXtZk97A0o5+l6f0UJw5iD9XQvYi1sig201pKbgy4mvyJoQa69kLHNqutI8mfFIqhbbM1xzDOAVWhNHt6BF4XvPYh6x+mtxowkFgC6adCypK3VgQEwBjY3xXPsw1pPNeQytbWJAzC8ow+ritt5Z2lbWNOBmuPQKnAnTifMOQTtrYm8nJjKi83plDZkYDHWP/7O8RHQaKbokQ3RYmDFCS6yYt3kxs/RG6Cm7z4IdJjPaEZzjfGmlPor7WO1uyvg6FO6zlHorX6MHMtZK6GzDWQMNJJvtNHh4YADv0WNt0MMSmQtsrKyM7xD9XeUN2Gz0DDQBL7utPZ35PO3p502t3xAJQldnFKejNrMo4xJ6F34nEppaZkyCfUDyRT05fMMVcCLYPxb966hmIxJ5yaGyNe0p2D1i3GRYZzkHSni6zYAQri+8iL68dpe2tieqyJ7YkH2/32xOBqBvzv5UiylqzGF/hvheCIf/vrpzCspENDYG0AO/a8tfxrnGVfzQMOKtsTqexI4Lm6uRzoSaPXY+0eTotxsSi5gxVpVaxKayXdORiK6JVSo4ixGUoTuylNPHnfgccndA7F0u6Oo8Nt/dnujqVjyLpf05/Cts5YBn1v/VcoGLJjByiI76Ugvo8aH5SnuChPcU19qDcmxVqIkrrMuu8bsgrkDTRYy1b7G6Bn/7D26dbEc3wBxOWDuyOgD7ATNXsSQUyytTlkmG63neruOKp74jjUHUdlewKVHQm0uN4qGZEf18fq9GMsTOlgUXIHObEDkb5KTCnl57AZsmJdZMW6Rm1jDPR7HbQMxnN0IJGjA0kcdSVydCCRyq5Mnmh8a9g4PXbozaRQnuJiXsoA5SkuiiY7N2GLsVYcDa9S4HXBQOOw5FAPXf6zmWMz4bSfTOKNxhbURCAilwI/A+zAvcaY/znhefE/fznQD9xojNkazJi+taWY3e2JVPfE0ep66/hIuxjmpw5wdn43y9L7WZrRx5L0fnbXabVCpWYyEUh0eEh09FCS2PO253wG5uTmU9UVzyH/B8ZDXfH8qz6NPw++9f+H0+ajNNnFutwebq+onVpA9jirmGVS6VuPefpgoAnK/21q1x5F0BKBiNiBnwMXAfXAGyLymDFmz7BmlwHz/be1wC/9fwZNVXc8Xa4hlid3kpfTT0FcH/nxfeTG9uMYVuzK1wu7dchfqVnNJlDf3EgcsNQJS7MA/8rynqEYGl2JNBzvRQwk4vIGacmSIxGSyyFtWXAuH5SrWtYAVcaYagAReRC4BhieCK4B7vMfWL9RRNJEJN8YE7ST2+8774Cu2lFKTVlyzBDJMZ0sSO5887FpnVgOoWAmgkKgbtj9ek7+tD9Sm0LgbYlARG4BbvHf7RWR/YROFtAawvebimiJNVrihOiJNVrihOiJNQLj/OZoTwQS69zRnghmIhhpSvXEtaqBtMEYcw9wz3QENVEisnm0JVeRJlpijZY4IXpijZY4IXpijZY4YeqxBnMPXj1QPOx+EXDiobyBtFFKKRVEwUwEbwDzRaRURJzA9cBjJ7R5DPiwWE4HuoI5P6CUUupkQRsaMsZ4ROQ24Cms5aO/NcZUisit/ufvBp7AWjpahbV89KZgxTMFYRmSmqRoiTVa4oToiTVa4oToiTVa4oQpxhp1JSaUUkpNr5lzxI5SSqlJ0USglFKznCYCPxG5VET2i0iViHx5hOcXicgGERkUkS+GI0Z/HOPF+QER2em/vSYiK8MRpz+W8WK9xh/ndhHZLCJnRWKcw9qtFhGviFwXyvhOiGG87+m5ItLl/55uF5FvRGKc/jbn+mOsFJEXQx3jsDjG+57+x7Dv527/z0BGBMaZKiKPi8gO//c08DlXY8ysv2FNZh8CygAnsANYckKbHGA18F3gixEc5xlAuv/ry4BNERxrEm/NU60A9kVinMPaPYe1wOG6CP6engv8IxzxTTDONKwqA3P893MiNdYT2l8FPBeJcQJfBb7v/zobaAecgVxfewSWN8thGGPcwPFyGG8yxjQbY94AhsIRoF8gcb5mjOnw392ItTcjHAKJtdf4f2qBREbYTBgC48bp92ngYaA5lMGdINBYwy2QON8PPGKMqQXr9yvEMR430e/pDcADIYns7QKJ0wDJ/mKeSViJIKC62ZoILKOVuog0E43zZmB9UCMaXUCxisi7RGQf8E8gOKUVxzZunCJSCLwLuDuEcY0k0H//df7hgfUisjQ0ob1NIHEuANJF5AUR2SIiHw5ZdG8X8O+UiCQAl2J9IAi1QOK8C1iMtSl3F/AZY4yPAMye8wjGFlCpiwgQcJwich5WIgjLuDuBlw/5G/A3EXkH8G3gwmAHdoJA4vxf4EvGGK+E9zCKQGLdCsw1xvSKyOXAo1jVfUMpkDgdwGnABUA8sEFENhpjDgQ7uBNM5Hf/KuBVY0x7EOMZTSBxXgJsB84HyoFnRORlY8zJJ/acQHsElmgpdRFQnCKyArgXuMYYE65SqxP6nhpjXgLKRSQr2IGdIJA4K4AHRaQGuA74hYi8MyTRvd24sRpjuo0xvf6vnwBiIvR7Wg88aYzpM8a0Ai8B4VjYMJGf0+sJz7AQBBbnTVjDbcYYUwUcBhYFdPVwTNBE2g3r00k1UMpbEzFLR2l7O+GbLB43TmAO1k7tMyL9ewrM463J4lOBhuP3IynOE9r/nvBNFgfyPc0b9j1dA9RG4vcUawjjWX/bBGA3sCwSv6f+dqlYY+6JEfxv/0vgdv/Xuf7fp6xArq9DQwRWDkNE8oDNQArgE5HPYs3aj9vtCmWcwDeATKxPrQAeE4YKigHGei1WrakhYAB4n/H/FEdYnBEhwFivAz4hIh6s7+n1kfg9NcbsFZEngZ1Yp7ffa4zZHco4A43V3/RdwNPGmL5QxziBOL8N/F5EdmENJX3JWL2tcWmJCaWUmuV0jkAppWY5TQRKKTXLaSJQSqlZThOBUkrNcpoIlFJqltNEoNQ4RCRNRD7p/7pARP4a7piUmk66fFSpcYhICVZFz2XhjkWpYNANZUqN73+wyl9sBw4Ci40xy0TkRuCdWBt8lgE/xtr1+SFgELjcGNMuIuXAz7FKA/cDHzPG7Av1X0Kp0ejQkFLj+zJwyBizCviPE55bhlVSeQ3WWRX9xphTgA3A8Yqa9wCfNsacBnwR+EUoglYqUNojUGpqnjfG9AA9ItIFPO5/fBewQkSSsA4LemhY5dLY0Iep1Og0ESg1NYPDvvYNu+/D+v2yAZ3+3oRSEUmHhpQaXw+QPJkX+osSHhaR9wCIJWznSCs1Ek0ESo3DWGc6vCoiu4EfTuISHwBuFpEdQCWRebykmsV0+ahSSs1y2iNQSqlZThOBUkrNcpoIlFJqltNEoJRSs5wmAqWUmuU0ESil1CyniUAppWa5/w8QddTSfl9kqwAAAABJRU5ErkJggg==",
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABCN0lEQVR4nO3dd3zb5bX48c+RZHnvvRKP7A04CWGUvVdbaAudcGkpbentvLfj1/bSdXu7b4G2lNLFLYWWQim0hFH2SALZiTMdx/GIHe9tWZb0/P74KmASD3lo2ef9eukVS3r01Ylj5+hZ5xFjDEoppWYvW7gDUEopFV6aCJRSapbTRKCUUrOcJgKllJrlNBEopdQs5wh3ABOVlZVlSkpKwh2GUkpFlS1btrQaY7JHei7qEkFJSQmbN28OdxhKKRVVROTIaM/p0JBSSs1ymgiUUmqW00SglFKznCYCpZSa5TQRKKXULKeJQCmlZjlNBEopNcsFPRGIiF1EtonIP0Z4TkTkDhGpEpGdInJqsONRSin1dqHoEXwG2DvKc5cB8/23W4BfhiAepZRSwwR1Z7GIFAFXAN8FPj9Ck2uA+4x1Os5GEUkTkXxjTGMw41Iz00+fORBw289dtCCIkSgVXYLdI/hf4D8B3yjPFwJ1w+7X+x9TSikVIkFLBCJyJdBsjNkyVrMRHjvp7EwRuUVENovI5paWlmmLUSmlVHB7BGcCV4tIDfAgcL6I/PGENvVA8bD7RcDREy9kjLnHGFNhjKnIzh6xeJ5SSqlJCtocgTHmK8BXAETkXOCLxpgPntDsMeA2EXkQWAt06fyAGm4i4/5KqckJeRlqEbkVwBhzN/AEcDlQBfQDN4U6HjXz9A566Hd7sImQHOcg1mEPd0hKRbSQJAJjzAvAC/6v7x72uAE+FYoY1MzWO+hhe10nB4/10O3yvO25nORYFuYms6QghbgYTQpKnSjqDqZRajifMWyr7WRDdRs+n6EkK5FVxWkkx8XgM4aOPjc1bf28XNXKpsPtVJSkc8qctHCHrVRE0USgotagx8v63U0caeunLCuRs+dnkZbgPKnd2rJMWnoG2VjdxmuH2jh4rJd3ripkfm5yGKJWKvJorSEVlVxDXh7e0kBdez/nL8rhyhX5IyaB47KTY7lqZQFXrcind9DDlXe+wh9eq8EanVRqdtMegYo6bo+Pv21roL3fzdUrC5ibmRjwa8uyk/hAShz/2nuM/3qskgffqOX8RTk4bCN/JtIdyGo20B6BiirGGJ7e00RL7yCXL8+bUBI4LjHWwdUrCzi9NIO9jT08srWBfrdn/BcqNUNpIlBR5fWadg619PGO+dmUZSVN+joiwtqyTC5flkdLzyAPvlFHa+/gNEaqVPTQRKCiRmPXAJuq21mUl8zKotRpueb83GSuO60InzH8ZXMd1S2903JdpaKJJgIVFYa8Pp6qPEZynINzF2YjMlKZqsnJTYnj+oo5pCc4eXxnI2/UtOsksppVNBGoqLDpcDtdA0NctCQ3KDuFk+IcvOe0IhbkJvHaoTaerGzC4x2taK5SM4smAhXxWnsH2VbbwZL8FIrSE4L2Pg67jUuX5nFGeSYHjvXy0JZ6GrsGgvZ+SkUKTQQqohljeOlgC067jbPmZwX9/USE1SUZXLUin45+N1fc8QovHdDS52pm030EKqIdae+nrn2Ad8zPIj6EdYLKspO4fvUcNh1u4yO/e51PnzePz1y4ALtt9LkJPSFNRSvtEaiI5fMZXqlqJTU+hhVFaSF//4xEJ3//1Flce2oRdzxXxQfv3URzjyvkcSgVbJoIVMR6srKJtl43p5dljPlJPJjinXZ+9J6V/OC6FWyr6+CKO15hw6G2sMSiVLBoIlARyecz/OxfB8lIcLIgAorDvbeimEc/dSbJcQ4+cO9Gfv58FT6fLjFVM4POEaiwGG88/VBLL/uP9XDp0jxs07hnYCoW5aXw2G1n8ZVHdvHDp/bz+uF2fvq+VWQkjl7sTqlooD0CFZG21XaSHOdgfs7ky0gEQ1KsgzuuX8V33rmMDYfauOrOV3Q3sop6QUsEIhInIq+LyA4RqRSRb47Q5lwR6RKR7f7bN4IVj4oezT0uGjoHWFWUhi1McwNjERE+ePpc/vqJdbiGvLz3VxvY39QT7rCUmrRg9ggGgfONMSuBVcClInL6CO1eNsas8t++FcR4VJTYXtdJjF1YWpAS7lDGtKIojYduXYfdJnzoN5voGhgKd0hKTUrQEoGxHO8zx/hvOrumxtQ36OFAUy9L8lOIjYLzhcuyk/i/m9fi9vp4dFsDg0PecIek1IQFdY5AROwish1oBp4xxmwaodk6//DRehFZOsp1bhGRzSKyuaVFd3nOZLsauvAaw8ritHCHErAFucn8+sMVdLuGWF/ZpAXrVNQJaiIwxniNMauAImCNiCw7oclWYK5/+OhO4NFRrnOPMabCGFORnZ0dzJBVGHl9hp31XZRmJZI+xrGTkWh1SQbnLMjmSFs/22o7wx2OUhMSklVDxphO4AXg0hMe7z4+fGSMeQKIEZHgF5RREelwax8DQ15WFE7PWQOhtrwwlfLsRF471KaH3KioEsxVQ9kikub/Oh64ENh3Qps88ReWF5E1/nh02+Ystaexm8RYO3Myg1dhNJhEhPMX5eB02Hh2b7MOEamoEcweQT7wvIjsBN7AmiP4h4jcKiK3+ttcB+wWkR3AHcD1Rn97ZqW+QQ81bX0sykuJmA1kk5HgdHDW/Cyaul3saewOdzhKBSRoO4uNMTuBU0Z4/O5hX98F3BWsGFT02H+sB2NgcV74y0lM1eK8ZHY3dPFqVRvzc5JxOnTfpops+hOqws4Yw57GbnJTYslMig13OFMmIrxjfjYDQ16213WGOxylxqWJQIVdS+8gbb1uluRH9gayichLjaMsK5EttR24dG+BinCaCFTY7W3swS4SEVVGp9PpZZm4PT7tFaiIp9VHVVj5jGF/Uw+l2YnEReBO4omcOnai7ORYSrMS2VnfxWlz04mx6+cuFZn0J1OFVUPHAANDXhZEWJXR6XLanHQGhry6gkhFNE0EKqyqmntx2ISSrMRwhxIUBWlx5KXEsb22U/cVqIiliUCFjc8Yqlp6KclKnLHDJiLCyqJUOgeGqOsYCHc4So1oZv72qajQ2Omi3+2NuMNnptu8nCTiYmzsqu8KdyhKjUgTgQqbg8092G1CSebMHBY6zmG3sTQ/lUOtvfQNesIdjlIn0USgwsIcHxbKTJgVO2+XFqZgDOzTk8xUBJr5v4EqIjV2uegb9DJvhg8LHZee4CQvJY69jd06aawijiYCFRYHm3ux24TSGbpaaCSL85Np63PT0qMlqlVk0USgQs4YQ3VLL8Xp8cQ6Im8TWbAsyE3GLsJeHR5SEUYTgQq5g829dLs8lGXNjmGh4+JirLMWqpp7dXhIRRRNBCrknt3bDDCrhoWOm5eTRO+ghx26lFRFEE0EKuSe3XuMnORYkuJmX6mrsqxEbALrdzeGOxSl3qSJQIVUe5+brbUds7I3ANbwUHF6Aut3NenwkIoYwTyzOE5EXheRHSJSKSLfHKGNiMgdIlIlIjtF5NRgxaMiw/P7mvGZ2TksdNy8nCRq2/u1EJ2KGMHsmw8C5xtjekUkBnhFRNYbYzYOa3MZMN9/Wwv80v+nmgl23n7SQ89tKicnPomrbL/G1jf6SzcmfjpoYYVbeXYSLxxoYf2uJpYWpIY7HKWCemaxAXr9d2P8txP7wtcA9/nbbhSRNBHJN8boAOoM5PYKLx5N5aq5bdgi9Hz60/vunNLrA0lg8U47a0szeGJ3I1+4eAEiEfrNULNGUOcIRMQuItuBZuAZY8ymE5oUAnXD7tf7HzvxOreIyGYR2dzS0hK0eFVwvd6STK/HzgWFneEOZXxBHr+/bHk+1S19HGzuHb+xUkEW1GUbxhgvsEpE0oC/icgyY8zuYU1G+ih00m+gMeYe4B6AiooKnWGLUs83pOK0+Tgzr4ftteGO5u1SvbUscj3GEtcjJPjaicGFDztuSaLHlkebvZxO+1yMTM8GuEuW5PL1R3fzzJ5jM+6IThV9QrJ+zxjTKSIvAJcCwxNBPVA87H4RcDQUManQe7ExldNze4h3+IL6PhMZ3knwtTDHvZF03xEM0GfLpt1RzhAJCB7iTBdp3lqyvftxSyJ1MWtoti8GmVpn+v5NteQkx3L/xiO4PWN/Pz530YIpvZdS4wlaIhCRbGDInwTigQuB75/Q7DHgNhF5EGuSuEvnB2am+j4nVd3x3DAvMob27GaQOUMbyPPsxkOs/z/4JbhtI+x2Nj7SfEcoGtpCuft5cmx7qHJeiMuWPqUYSjITeaOmHdeQNyLPa1azRzB7BPnAH0TEjjUX8RdjzD9E5FYAY8zdwBPA5UAV0A/cFMR4VBi91GitjjknP/w7apO9R5nvfhqn6eOoYyX1MWvwSuzoLxAbnfZSOm0lZHoPUup+kRWuP1PlvJB2x7xJx1GSlcDrNe3Utvfr8JAKq2CuGtoJnDLC43cP+9oAnwpWDCpyvHg0lcLEQcpTXOELwhjyPDsoGXqVQUlhd+y19NrzAn+9CG2OBXTbCljoXs9C95McMWdwNGZy219yU+KIi7FR09qniUCF1ezb469CbsgnvNqUwtUlbYRrpaQYL6VDL5HrqaTNXsYh54V4xTmpaw3ZkqiMfRfz3M8yd+g17GaQupjTmehfzibC3IxEatr6McboMlIVNpoIVNBtbU2k12MP27CQw7hYMPgkqb566h2nTeo/7RMZcXDQeRFedwxFni14xcnRmNMmfJ2SrAT2H+vhWM8geSlxU4pJqcnSRKCC7sWjqTjExxl5oS+pEOfrYNHgP4g1PRx0XkirY9H0XVxsVDvPw+b2MHdoA0MSP+FLzM2wSm3UtPZpIlBho4lABd2Ljamcmt1Lckxwl42eKMVbz8LB9RiEPbHvpMdeMP1vIsIh5wXEDLoodz/PZV2fp8NROqFLPJ90Oq3NnVCWOf3xKRUArT6qgqp5wEFlRyLn5Ie2N5A7tJvFg4/hlkR2xb03OEnAz4id/bGX0mvLYYH7SRJ9zRN6/aq0Fqr7Uul3e4IUoVJj00SggurlEC8btRk3pe7nKRt6gS5bMbvjrmPQlhL09/WJk32xVzIkCSwcXI/DDAT82pVpLRiEuvbAX6PUdNJEoILqxcZUsuKGWJLeH/T3SvC1cF3njeR5Kql3nMa+2CsmvTJoMjwSz/7Yy4gx/cwffApMYENhpYndJNrd1LYH/3uk1Eg0Eaig8frg5aYU3pHfFfRqo+WDz/Ch9qvI8VRywHkxdc51Uy4DMRl9thyqneeS5qtnztCGgF5jE1iaam0s08NqVDhoIlBBs7sjkY7BmKAOCzl9PVzc/RWu7r6NHnsBf0p/mDZHeGvztDgW0+RYTqFnG+me6oBesyy1jd5BD50DQ0GOTqmT6aohFTQvHk1BMJwdhGWjNuNmxcCDrO3/BXGmi40Jn2BTwifxidM6EinMamLOIsnXRLn7WXbackauYTTMitRWAGrb+0lPCN1wllKgiUBNo58+c+Bt99cfTqQssYv9R49N23s4zAALXf9kTf/dpPnqqI05nZcT/5PmmKXT9h7TwYidg85LWOH6M/PcT7Mn9p1jDlXlxA6QEuegrr2flUVpIYtTKdBEoILENeTlYG8a7yo8NPWLGUOWdz9LXQ+zxPUocaabZvsiHkn9NUdizp7yLuFgcdnSqHaew3z3vyj0bKEhZvWobUWgOCOBg829+HwGW6Qe4aZmJE0EKijq2vsxCCvTWid3AeMj37ODeYNPM2/wGdJ8dXiJ4WDsxeyMu4GGmIqITQDDtdoXkmavpXjodbptRfTY80dtOycjgcqj3RzrcZGfOvFdykpNliYCFRQ1bf0k2t2UJ01gotgYknzHyPIeYLnrzyT5WvASQ61zHW84b+FQ7IUM2DKCF3QwiHDYeS7Jribmu59mR9z78MrIpSSK0xMAqGsf0ESgQkoTgZp2xhiOtPexPLUNu4y/HNJm3OR49pDv2Umc6caHnUPO8zkYewmHnefitkV3iWavODnovISlgw9T7n6eA85LR+zNxDvt5CTHUtvez5rSKEt4KqppIlDTrrXXTd+glxWFYw8L2Y2bfM928oZ2EMMg3bYC6mNW024v49WkL4Yo2tDotedSF7OWuUMbyPHupdmxZMR2xRkJbKvtwO3x4XTo6m4VGkH7SRORYhF5XkT2ikiliHxmhDbnikiXiGz3374RrHhU6NS09QGMOT+Q6TnIKtcfKR56nR57Abtir6Uy7t20OBaPfVpYFDvqOJUuWxEl7peI83WM2KY4PR6fgaOdWm5ChU4wewQe4AvGmK0ikgxsEZFnjDF7Tmj3sjHmyiDGoULsSFs/2cmxZDhPXtBvN25K3S+Q7T1Ary2Hfc4r6LPnhiHKMBChynkhK1wPMH/waXbHXYeRt59VXJgWj90m1Lb3U5KVGKZA1WwTtB6BMabRGLPV/3UPsBcoDNb7qcgwOOTlaNcAJZkJJz0X5+tkmeshsrwHqY1Zy67Y62ZPEvBz25I45LyAJNNC8dDGk5532G0UpMZR26F1h1TohGQQUkRKsM4v3jTC0+tEZIeIrBeRyNoVpCbMqpcDJZlv/zSb6D3GMtdDxJh+9sReY62pD0MtoEjQ4SijybGMQs82Ur11Jz0/JyOBtl43fYNallqFRtB/E0UkCXgY+Kwx5sRaA1uBucaYlcCdwKOjXOMWEdksIptbWlqCGq+ampq2fmIdtredtpXsPcrSwUfxSiy74t5Lt70ojBFGhiMxZ9Iv6cxzP3NSyeriDP8yUu0VqBAJaiIQkRisJHC/MeaRE583xnQbY3r9Xz8BxIhI1gjt7jHGVBhjKrKzs4MZspoCYww1bX3MzUh4c2dskreJxYOPMyhJVMa+m0FbapijjAw+ieFg7CU4jIty93MwrOpodnIscQ6blqVWIRPMVUMC/AbYa4z5ySht8vztEJE1/njaghWTCq6WnkH63d43JznjfB0sHnycIUlgT+w14xZem236bVkciTmDDO9hcj2733zcJkJxRgJ17QNallqFRECJQEQeFpErRCY0qHsm8CHg/GHLQy8XkVtF5FZ/m+uA3SKyA7gDuN7oT37UqmmzPsHOzUwg1tfJosF/YLCxJ/YahjQJjKjJsZIO2xxKhl4h03PwzcfnZCTQO+iho1/LUqvgC3T56C+Bm4A7ROQh4PfGmH1jvcAY8wowZjEYY8xdwF0BxqAiXE1bHznJsSTFGK7s+iyxpoc9se8KyVGRUUuEQ7EXsmLgAS7r/gIPpD+EV2LfnCfQ4SEVCgElAmPMv4B/iUgqcAPwjIjUAb8G/miM0Y8ts1xnv5umLherSzM4t/e7zBnaQJXzgjGLrI3l9L47pznCyDUkCRyKvYDFg//gnN7v8Vzy7aTGx5AaH6OJQIVEwEM9IpIJ3Ah8FNgG/Aw4FXgmKJGpqPLSwVYM8O60Z1npeoA34j9Ki2NxuMOKGp32EjbH38xK1wMscj0GQHFGPA0dAwx5Azv7WKnJCnSO4BHgZSABuMoYc7Ux5s/GmE8DOvireGFfM4sTj/Ju7/eoi1nDq4mfD3dIUeeVxM9TH7OaC3u+QabnIHMyEnB7feyo6wx3aGqGC7RHcK8xZokx5nvGmEYAEasgjDGmImjRqajg8xk2HGzgrjk/wCtO1if/6KTSCWp8Rhw8kfwT3JLIld2fpjzNWjfx8sFJnumgVIACTQTfGeGxDdMZiIpeOxu6+HjK3ZTHVPFU8vdnXdmI6dRnz+GJlJ+Q5j3Cla5vkJvi5NUqTQQquMZMBP51/qcB8SJyioic6r+dizVMpBTV2//MjVn/4PXYD3M49txwhxP16p1reTXxcywcXM/H855gW10nPS5dj6GCZ7xVQ5dgTRAXAcM3hfUAXw1STCqaDLZxTudXqaWMjcn/Ee5oZozN8R+lYGgbH+EunorPZWP1aVy0RHtaKjjG7BEYY/5gjDkPuNEYc96w29UjlYxQs0/fq58gxdbF1sKf4RVnuMOZOcTGk8k/oNNWzN1zv0flge3hjkjNYOMNDX3Q/2WJiHz+xFsI4lORrPYhEpse4mfHbuC0VeeHO5oZx21L5rHUX+K0G67u+iQM9YY7JDVDjTdZfLyWcBKQPMJNzVYDx+CNT1DlXcxztpve3Amrpleno5QXs+9krqMG10sfAKN7CtT0G3OOwBjzK/+f3wxNOCpqbPl3zFAPnzj0HS47Q88bCqbS5e/iu/e/yjfk17Drm7BCfx3V9Ap0Q9kPRCRFRGJE5FkRaR02bKRmm/rHoPYv7E7/DAddxVy0JC/cEc1oi/KSecz1HjZwNez+FtQ+FO6Q1AwT6D6Ci/2HylwJ1AMLAF0iMhu5u+CNT0Dacn5+7F3kp8axrFCLygWTiHDuwhw+XfUxTOYZsOEj0L4l3GGpGSTQRBDj//Ny4AFjTHuQ4lGRbvuXwNWE69R7eOFgJxctycV/pIQKovMW5tDaL2wv/Q3EZsOLV0N/Q7jDUjNEoIngcRHZB1QAz4pINuAKXlgqIjW/BFW/goWf5bmWubiGfFyyVIeFQuHsBVnYbcIz1cA5j8NQN7x0DXi0OqmauoASgTHmy8A6oMJfcroPuCaYgakI43XBpo9BYims+Bb/3NlIVpKTtaUZ4Y5sVkiJi6FibjrP7WuG9BVwxp+gfas1TKQridQUTeTEscXA+0Tkw1gni10cnJBURNr1Leg5AGvvod8Xy3P7mrl0WR4Oe1CPvVbDnL8oh31NPTR2DUDRVXDKD6Dur9ZKIqWmINBVQ/8H/Ag4C1jtv41ZdVREikXkeRHZKyKVIvKZEdqIiNwhIlUislNETp3E30EFW8d22PsDKLsR8i7k+X0tDAx5uWJ5Qbgjm1XOW5QDwPP7WqwHFn0Bym6yVhLVPBDGyFS0C/SoygpgyQTPE/YAXzDGbBWRZGCLiDxjjNkzrM1lwHz/bS3WkZhrJ/AeKth8Hth4M8RmwSk/BuCfu46SlRTLGh0WCqn5OUkUpsXz/P5m3r92DojA6ruh9xBsvAmSyiBLf33UxAXar98NTGhW0BjTaIzZ6v+6B9gLnLjz6BrgPmPZCKSJyOTONlTBsf8O6NgKp90BsRn0DXp4bl8zly/Pw27T1UKhJCKcvyiHV6taGfR4rQftTjjrYYgvsCaP++rCG6SKSoEmgixgj4g8JSKPHb8F+iYiUgKcAmw64alCYPhPbj0nJwtE5BYR2Swim1taWgJ9WzVVfXWw6xtQcAXMeQ8Az+1rxjXk44rlmq/D4bxF2fS7vWyqHraCOy4Lzv0HeAfgxau0JpGasECHhm6f7BuISBLwMPBZ/6a0tz09wktOGn4yxtwD3ANQUVExkeEpNRVbP2utSKm40xqGAP65s5Gc5FgqSnRYKBzOKM8iPsbO03uaeMeC7LeeOPIXKLgKjvwJnlwNc9775r/ZuFbcHpRYVfQIdPnoi0ANEOP/+g1g63ivE5EYrCRw/yhlq+uB4mH3i4CjgcSkgqzhCah7BJZ9HZJKAeh2DfH8/mYuX56vw0JhEhdj5/xFOTy5+xhe3wmfiZLnQf4l0LMPjj0XngBVVAqoRyAiHwNuATKAcqzhm7uBC8Z4jQC/AfYaY34ySrPHgNtE5EGsSeKu42ciqzDy9MPm2yBlMT9ruArf0QMA7D7axaDHx4Dby0+fORDmIGee0/vuHPmJnZlvu3tpSgb/7C1n83M/Zm3OCcNAGWvA1QKtr1gT/OkrgxStmkkCHRr6FLAG/xi/MeagiOSM85ozgQ8Bu0Rku/+xrwJz/Ne4G3gCq2xFFdAP3DSR4FWQVH4X+g7DBS+wZsOv3nz46fo15MfFco3910hf+MKb7c4r6CTW7mN9bcbJiUAECi4DdzscfRxiMyCheOQLKeUX6GTxoDHGffyOiDgYYSx/OGPMK8YYMcasMMas8t+eMMbc7U8C+FcLfcoYU26MWW6M2Tz5v4qaFl17Ye8PofTDkHvOmw+3DsaxtyeDs7KOBjz0rIIjKcbHOfldPFmXzomjQwCIHYrfA44Uq1KpRyeP1dgCTQQvishXsQ6xvwh4CHg8eGGpsDAG3vgkOJLglB++7alXWq3NY2dl6RROJLisuIOmASfb2xJHbuCItyaMvQNQ97CWoVBjCjQRfBloAXYBH8ca0vlasIJSYVLzR2h+AVZ9H+LeGvkzxkoEC5PbyYkbCF986k0XFHYSY7OGh0YVnwcFV0JfDRx7NmSxqegT0ByBMcYnIo8CjxpjdCF/NNl5e2DtPANw8C6IL4K++re97nBfCg0DSdxcWhmEANVkpDi9nJ3XzRN16Xz1lLrRh+vSV8JAPbS+BvGFkLokpHGq6DDe4fUiIreLSCuwD9gvIi0i8o3QhKdCpvk5axih4IqT1p+/3FqAQ3ycnqkLuiLJpcUdNPTFjj48dFzepVaCb/i7taJIqROMNzT0WazVP6uNMZnGmAysZZ5nisjngh2cChHXMevEq4wKazhhGI9PeLW1gFPTm0lyeMIUoBrJpcUdxNm9/LU6a+yGNru1M1wcUPcX8A6GJkAVNcZLBB8GbjDGHD7+gDGmGvig/zkV7YyBxifBHgc555309JaOHHo8Ts7Lrg9DcGosKU4vlxV38NiRDFyecZZyxaRA8XUw2AYNj1n/7kr5jZcIYowxrSc+6J8niBmhvYo23fusycSc86yVJid4vrmIDOcAK9JO+jFQEeA9Za30DDl4qj59/MZJpZB7AXTvgbaNwQ9ORY3xEoF7ks+paODzQNPTEJsDGaed9PTRPic7u7I4J7sBrSgRmU7P7aEocZCHxhseOi7rDEheBE3/gn7t5SnLeIlgpYh0j3DrAZaHIkAVRK2vwVCnVZ9GTv5R+Gt1JgbhnGw9JD1S2QSuLW3l1aYUGvqc479ABIquhphk63Qzry4HVuMkAmOM3RiTMsIt2RijQ0PRbKgbWl6BlEXWgSYn8Bl4qDqLpSlt5OregYh2XVkrBuHh6szxGwPY4635gqEeqNf5AjWxM4vVTNL0LOCDvJGPnt54LJm6vjjOy9Hhg0hXnOTmjNxuHqrOGrnkxEgSiiDvQqtS6YFRit2pWUMTwWzU3wBdOyFzHThHnmT8U1U2KTEeVmccC3FwajKuL2+hri+OF46mBv6izNMheQFs+yK0aZmv2UwTwWxjDBz7F9gTIPusEZu0DFirUK4ra8Vp0xo10eCyOR3kxrv53f7cwF8kAoXvhLg8eOW94O4MVngqwmkimG16D/mXi74D7LEjNvlLdTZDPhsfmK+7UKNFjM3wofnNvNyUysGuuMBf6IiHMx+E/lrY9FGdL5ilAj2PQM0Ex3sDMWmQXjFiE6/PGhZal9tNeYqLZt0+EHYbqtsCanfDPAd3VeZzz948fnh6TeBvkH0GrPxv2P4lOPgLWPCpyQWqopb2CGaTrl1WOYnc862yAyN4sTGVhr5YPji/OcTBqanKjPNwfXkrfzucGdhS0uEWfxEKLoetn4f2cU+hVTOMJoLZwuexzrGNy4fUZaM2u/9gDtlxbi4u6gxdbGrafGxxEwD37M0bp+UJxAan/wFis635gqHuIESnIlXQEoGI/FZEmkVk9yjPnysiXSKy3X/TiqbB1P4GDHVZSwZHqVlc3+fkuaOpXF/eSoxNx4qjUWGim+vK2nigKpu63gn2CuKyrPmCvsPWAUU6XzBrBHOO4PfAXcB9Y7R52RhzZRBjUABeF7S8DEnlI24eO+6BqmxE4Pp5OkkczT6zvIFHazL40Y4ifnZm9fgv2Hn72+9nvwNq7oehXkhfNfZrV9w+9vMqKgStR2CMeQloD9b11QS0vmqVEsi9YNQmbq/w50PZnF/QSWGilpGKZvkJQ9y86Bh/P5LJzraEiV8g+2xImAuNT1jVStWMF+45gnUiskNE1ovI0tEaicgtIrJZRDa3tOin1Qnx9EHbJkhdCvH5ozZ7uj6NVleMLhmdIW5d0khG7BD/va144iM8YoPid/vPL/irNb+kZrRwJoKtwFxjzErgTuDR0RoaY+4xxlQYYyqys7NDFd/M0Poa+IYg+5wxm913MIeixEHOye8KUWAqmJJjfHxm2VE2Nqfw/ER2Gx8XkwKF14CryVpyrGa0sCUCY0y3MabX//UTQIyIBFhLVwXE1Qxtb0DqcogbPYFWdsTzenMKH1lwTMtNzyA3zGuhLHmA27fMod8ziV/1lIWQucbqUXYfmP4AVcQIWyIQkTwRa/mKiKzxx6IDktNpzw/AeCBn7N7A7/bnkuDw8t5y3T02kzjthu+traG2N44f7yic3EVyL7JKUDQ8qktKZ7CgrRoSkQeAc4EsEakH/gv/qWbGmLuB64BPiIgHGACuN0bXq02bgSZrl2jaCogdvTxxq8vBYzWZvK+8hVSnN4QBqlBYm9PLB+c389v9uVwxt51Ts/qAwHcrA6wrvhYO3QP1f4OSD414doWKbkFLBMaYG8Z5/i6s5aUqGPZ8H3xuayngGO4/mIPbZ+PGhbqTeKb60qo6nmtI5UubSvnHpZXE2if4eSs2C/Ivh4a/W2dY5Iz9M6Wij9YamokGGqHqbij9EMRmnPT08U+DHp/wu33LWZXWQnNrg9YVmqGSY3x8d80RbnphAT/ZWchXTpnEGRNpK6G3GppfgMQSSJwz3WGqMNI+3kxU+T/WSqFlXx+z2Ya2PDqH4rg070iIAlPhcl5BFx+Y18yv9ubzwtGUiV9ABAquAGca1D+iR1zOMJoIZpr+Bqj6FZTdOOYuYmNgfVMJhfG9rEjVrsBs8PVTa1mU1s8XNpTR7h65BPmY7LFQdK11xGXD41qCYgbRRDDTVH4PjBeWfm3MZgd60zjcl8oleUdGKz2kZpg4h+GuMw/R77Hx86oVgR9rOVxCIeRdAN17oWPLtMeowkMTwUzSVweHfg3l/wZJJWM2Xd84l0T7EGdnHQ1NbCoizEt18c2KWvZ0Z/JI/bzJXSRznVW3qvEp6ByxpqSKMpoIZpLK/wYMLP1/YzY75orn9fY8zs+pI86uS0Znm/eUtXJOdj2PNJSzo3MSezhFoOid1lDRq+8DT/+0x6hCSxPBTNF3BKp/A2U3j7ui45+NJdjFcGm+ThLPRiJwU8keihJ6uatqBW2DEzja8jhHEhS9C7r2wNbPTX+QKqQ0EcwUld+z/lz61TGbtbocvNBcxNlZDWQ4B0MQmIpEsXYfn5u/Da+x8b8HV+HxTWKiKKkclnwJqu6B2oenP0gVMpoIZoK+I1D9Wyj/KCQWj9n0vgM5eIyNKwpqQhObilj58f18vGwXVb1p3F+7cHIXWfFtyKiA12+xVqypqKSJYCao/G9AYMlXxmzWN2TjDwdyOS29mcL4vtDEpiLa2sxjXJZXw5NNJWxom+DxlgC2GDjjfuvwo403gvFNe4wq+HRncbTrOwKHfgvzbhm3N/DgoWy63A6unh/AqVUq6kykftBw75+zn6reVO45tIy5CT0UTPRDQsoCOO2n8PrHYf/PYJHOGUQb7RFEu93ftYqALR27N+DyCr/em8eanG7mJ+uZA+otDpvhM/N3EGPz8b8HVuH2TeK/hfKPQeHVsP3L0Llr+oNUQaWJIJr11kD176xfwoSiMZv+5VA2TQNOPrNM9w2ok2XGuvhE+U7qBpJ5sHbBxC8gAmvvBWc6vPp+a6hIRQ1NBNGs8nhv4MtjNhv0Cr+ozGd1dg9n5PaEKDgVbU5Jb+Xi3COsbyphZ+fopctHFZcNp/8OunbD9rFXr6nIookgWvUehurfW3MD4/YGsqzewPKjWk5CjekDc/dTGN/LLw8tp3soZuIXKLgM5n8K9v8UGp+Z/gBVUGgiiFbHewNLAugN7CmgIruHM3P1hCk1NqfNx23zdtDjcfLbw0smd5FTfggpi2HjR2BQDx2MBsE8oey3wJVAszFm2QjPC/Az4HKgH7jRGLM1WPHMKL3VUP0HmH+rVQRsDH+pzqKx38kP1h7W3oAKSEliD9cWVfGXugW80d7IurIxypDsvH3kx7PPhup74Zl3QPF1jPrDt2KU16uQCmaP4PfApWM8fxkw33+7BfhlEGOZWXZ/F8Q+bm+g3+3hjl0FrM7u4aw87Q2owF2Vf5i5Cd389vASutz2iV8gPh9yzoPuPdClq4giXdASgTHmJaB9jCbXAPcZy0YgTUTygxXPjNG9Hw7/wT83MHZv4LevHKbF5eTLq+q1N6AmxGEzfLx8N91DTr6zdez9KaPKOgMSiuHoE+DWJcuRLJwbygqBumH36/2PNYYnnAi28/a3vq59yJobsDlH75YD7S4Hv3p+ORcV9XBadm+QA1QzUWliN1cW1PBQdRlXz23n7PwJ9irFZhWmq7obGh6Fkg+PPkSkwiqck8Uj/USMeFSGiNwiIptFZHNLS0uQw4pgA0etrnbmOqv64xh+XplPn8fOf66cxPm0SvldW1RFabKL/9o8h0HvJP4Td6ZD3iXQVwNtG6c9PjU9wpkI6oHhfc4iYMTdTsaYe4wxFcaYiuzs7JAEF5GOPQf2eMhaN2az+j4n/3cwh/eUtTI/VTf2qMlz2nz812lHqO6J5959k6hFBJB+CiQvgGPPgqt5egNU0yKcieAx4MNiOR3oMsbosNBoeg9D7yFrNYZ97Prx/7OtCJsYPrtcq0GqqTu3oJtLijq4a3c+DX3OiV9ABAqvAlss1P8NfHoYUqQJWiIQkQeADcBCEakXkZtF5FYRudXf5AmgGqgCfg18MlixRD1jrE9TjhTIWD1m09ebk/hHbSa3LmkiP2EoRAGqme7rp9ViYPITx44kKxm4mqD5hekMTU2DoE0WG2NuGOd5A3wqWO8/o3TthoEGKLwGbKP/k3l98M0tcyhIGOTji5tCGKCayY5XNb26IJ6/1C3gV1sPsSJt5I1i68rGKE2RsgjSVkHrq5A8f9yT9FTo6M7iSOfpg6Z/QVw+pK0cs+lD1VlUdiTylVPqiXdoXXg1va7MP0xeXB+/r1nC0GRONAPIvxRiUqH+UfDqCXmRQhNBpNv7I/B0Q/4lYy6963Lb+eGOIlZn93DlnLG2byg1OTE2w0dK9tLoSuSfjaWTu4g91jr4fqgDmp6e1vjU5GkiiGT99bDn+5CyBBLnjtn0RzsK6XA7+K/TanWptgqaVWmtVKQf49GGMtrdsZO7SOJcyDoTOrZC/ePTG6CaFD2hLJJt/woer5ed7goGxzh9qqo3lT8ezOGmhcdYltEfwgDVbPTBufv4jx1n8afahdw2b+fkLpJzLvRWwesfhazdVglrFTbaI4hUx16Emj+yJeHfGLSljNrMa4R7q5eS7hzk8yt0uagKvty4Aa7Ir+HV1gIO9KRN7iI2h7Xr2N1pHXxvRtxLqkJEewSRyDsIb3wcEkt5Pf5WKvrvHbXpU01zONKfwmfnb2NX3Szeda1C6prCal5sKeT3NYv5zrIN2CYzHBmXCyu/C9v+wzpbo/ym6Q5TBUh7BJGo8ntWcbnVv8Qj8aM2axuM5aG6+axKa2ZNxrEQBqhmuzi7lw/M3c/hvlRebBm7+OGYFn4Ocs6BLZ+xNk2qsNBEEGm69sGe78HcG6DgkjGb3ndkMT6Em0r26gSxCrkzMhtZmNzBg7UL6PNMcnDBZod1f7C+3vAR3XUcJpoIIonxWUNC9gQ49adjNt3Skc3r7Xm8u7CKnLiBEAWo1FtE4CMle+nxOHmkft7kL5Q4FyruhJaXYd9Ppi9AFTBNBJFk/8+g+SU49UcQnztqs36Pg98eXkpxQg9X5NeELj6lTlCa2M15OfU8dWwODQOJU7jQh6H43bDza9AxyZVIatI0EUSKzl2w/ctQeDWU/duYTf9Uu5AOdywfL9uFw6arLVR4va/4ALE2L/fVLJr84h8RWP0rq2z1hg/qruMQ00QQCbwueO0D1i/B2nvH3EFc2ZXBs83FXJFfQ3mSHj+pwi8lZohri6rY2ZXNsw2pk79QXBas/Y31oWjn16cvQDUuTQSRYMf/s3741/52zI01g14b91QvIy+uj+uKDoYwQKXGdnFuLYXxvXx76yQPsDmu8ArrGNa9P7KGSVVIaCIIt/rHrAmy+Z+EwsvHbPqX+vk0DybwsbLdxNq1qJyKHA6b4cNz93KkN27yB9gcd8qPIakcXn0/DOiy6FDQDWWhsvP2kx8bbIVDv4b4AohJH7HN6X1WaYmqnlTWN5ZwYU4tS1I6ghmpUpOyIq2NS4o6uHN3PlfPbaM4yT25C8UkwdkPwdNnwCvvgQueBVvM9Aar3kZ7BOHiHYTaP4M4oPi9Y54z4PEJv6peRrrTxQ1z9ocwSKUm5vaKIzgE/t8bJVOrGpG+ypova3kZtn5husJTo9BEEA7GWPXYB9ug+D3gHHuC7dGGcuoHkrm5dA8JDt1woyJXfsIQX1xZz0uNqTx2JGNqFyt5v7Xz+MCdUH3f9ASoRhTURCAil4rIfhGpEpEvj/D8uSLSJSLb/bdvBDOeiGAMND0FPfsg72JIKhmz+eG+FB49WsaZWUc5NV1rCanI96H5zazM7OVbW+bQOWif2sVO+QHknmdVKW16bnoCVCcJ5pnFduDnwGXAEuAGEVkyQtOXjTGr/LdvBSueiNH6GrRtgsy11m0MLq/w86rlpDjc3FiyN0QBKjU1dht8b00NnW4H39wyxeMobQ44+2FIXgAvvRM6dkxLjOrtgtkjWANUGWOqjTFu4EHgmiC+X+Tr2A7H/gWpyyBv7BPHAH64vYiGgWQ+Xr6bJIceRK+ix5L0AW5bepS/1WSxvjZ9ahdzpsN5T1pDqC9cBr010xKjekswE0EhUDfsfr3/sROtE5EdIrJeRJaOdCERuUVENovI5paWKB0e6doDDY9BYpl1CP04SeC1Y8n8Zn8eF+UeYWVaa4iCVGr63LaskRUZfXz1jbk0D0xxgWJCEZz7JHgG4LmLoK92eoJUQHCXj470P92J6wi2AnONMb0icjnwKDD/pBcZcw9wD0BFRUX01VSo+RPU/dX6YZ4z9gohgG63nS9uKKU02cX75xwIUZBKTa8Ym+En66q54smlfHlTKb855+DJn3923j6xixZfCzV/hCdWwKWbIXkKxe7Um4LZI6gHiofdLwKODm9gjOk2xvT6v34CiBGRrCDGFHrVf4DXPmhVWJz7Qevw7nHcvmUOxwac/GRdNXF2XSWkote8VBdfWlnPc0fT+P2BnKlfMKEISj8CPjc8czZ0Vk79miqoieANYL6IlIqIE7geeGx4AxHJE7E+I4jIGn88ox/OG02Mgb0/gY03Qd6FMPf9YHeO+7K/12TwyOEsPrX0KKdk9YUgUKWC68aFx7iwsIPvbi1mS0vS1C8Ynw+lN1pfP3MWHF0/9WvOckFLBMYYD3Ab8BSwF/iLMaZSRG4VkVv9za4DdovIDuAO4HpjZsDhpT4PbL4Ntn3B6sqe81hAOyOruuL4yuslrM7u4dPLGkMQqFLBZxP48brDFCa6+eQr5bRMdb4AIC4HLn4VEufAC1fArm9b53moSQnqPgJjzBPGmAXGmHJjzHf9j91tjLnb//VdxpilxpiVxpjTjTGvBTOekBjqsZa5HfwFLP5POOvPYI8b92X9HhuffKWceLuPO888RIyWl1YzSKrTyy/PrqLLbefTr5Yz5JuGI/WSyuDiDVDyAdj1DXjxKuivn/p1ZyHdWTydOivhqdXQ+CSsvhtO+T7I+N9iY+Brb8zlYFc8/3tGNXkJulRUzTxL0gf47zVH2Nicwpc2TbEExXGOBFh3H1TcBceeh38sgf136pGXE6SJYLocvh+eWgPuTjj/WZj/8YBfevfePB45nMVnlh/l7Hw9Y0DNXO8ubePzy+t55HAWP9hRND0XFYEFn4IrdkPWOtjy7/D06XD0KaYn28x8mgimaqgHNt5snaqUcRpctg1yzwn45f+sTef724u5ck4b/77s6PgvUCrKfXpZI++f18wv9+Tzu/3TsJLouKQya+PZGfeD6xi8cCk8cyYcfVLnD8ahiWAqml+GJ1bC4d/Dkq9Y5XLj8wN++ZaWRD73WhmnZfXwo3WHsU3DsKlSkU4Evl1xhIuLOvjmlrncu2/087kndfGS98NVB2H1L605gxcug8cXwp4fgitKN6QGmZ5HMBlDPbDzG9Zh80mlcOHLkH3GhC6xsy2Bm15YQH6Cm3veUUWcXbuwKvptqA5s9fe6skzuOvMQn32tjO9snUPfkJ1/X3Z0vA33gbPHwvxboewmqP0rVP0Ktv8n7Pgq5J4Pc66Fondaq48UEm2rNSsqKszmzZvD8+Y7b4fufda6ZU83ZFRA7kUB7Q8YbltrIh9+fgGpTi8PXLBvzAM8Av3FUiqarCvLBMDjgy+/XsJfq7P5yIJjfO3UuuCtmHO1QOcO6N4L7nZrIUf2WVbJl/xLIHXJuKVfopmIbDHGVIz0nPYIAtW1B2ruh94qiMuFOe+xdjlO0KbmJG5+YQGZcUP86YL9FCZO8hQnpWYAhw1+sLaGdKeXX+/LY29HAnedVUVOvGf63ywu29rcmXuB1SOoe9i6bfuCdUsoskrD518CeRdAbOb0xxChtEcwHlcL7Lrd6lqKA3LeYZWPlonVWTcG/u9gDt/aUsycpEHuv2A/+QEsE9UegZqJjvcIhvt7TQZf3lRCstPLT9Yd5qy8EK2gc3dB7yH/rRp8Luvx+ELr7OSkMogvApv/d37F7aGJa5ppj2AyXM2w90dw4OfgG4R5t4I9HhyJE77UgMfG196Yy8OHs7igoJOfnFFNqlPXOSs13DUl7SxMG+CTL5fzwecWcm1pK//vlDoy4oLQOxjOmQoZp1o344OBo1bPv+eQdVRmy0sgMdYu5sRSaN9qHaUZwB6haKE9ghO9LQG4YO4NsOzrkLJw4pUSgecbUvn65rnU98Xy2eUN/PuyoxNaHaQ9AjUTjdQjOM7lEe6sLOBXe/JIcXr55NJGPjCvhXhHGJaAel3QVwO9h6HvMAz6Vx05M6yT0/IusIaakudH/PzCWD2C2ZUIdt4++nOuFmjfZJ2AZLzW4TE574DYyRVDreqK48c7C1lfl8G8lAG+u6aGtTm9gP7nrlQgavuTuK9mMZXdmWTFDXHL4iauK20Nfg9hLEM9Vs/h2LPQ9C/o9x+5klBkJYTjiSGhIHwxjkITwXE7b3/7fWOsLmDbJmt8UOyQuhyyz5x0AtjdnsAvKvNZX5dOrN3Hp5c18rFFTTiHLQ/VRKBU4BxJc/nZ7gJeaUolxubjgsJO3lnSzhm53aSEc4jVGGv1Ud/ht3oM3gHrOWemNZSU4L8509/qMYRpjkHnCE401G0dG9m5w/qHdCRBznnWzuBJzAG0uRw8fiSDh6qzqOxIJDnGwyeXNvJvC4+RGc5PL0rNAKtzevnj+QfY2xHPw4ezeLQmkyfrMrCL4ZSsXtZk97A0o5+l6f0UJw5iD9XQvYi1sig201pKbgy4mvyJoQa69kLHNqutI8mfFIqhbbM1xzDOAVWhNHt6BF4XvPYh6x+mtxowkFgC6adCypK3VgQEwBjY3xXPsw1pPNeQytbWJAzC8ow+ritt5Z2lbWNOBmuPQKnAnTifMOQTtrYm8nJjKi83plDZkYDHWP/7O8RHQaKbokQ3RYmDFCS6yYt3kxs/RG6Cm7z4IdJjPaEZzjfGmlPor7WO1uyvg6FO6zlHorX6MHMtZK6GzDWQMNJJvtNHh4YADv0WNt0MMSmQtsrKyM7xD9XeUN2Gz0DDQBL7utPZ35PO3p502t3xAJQldnFKejNrMo4xJ6F34nEppaZkyCfUDyRT05fMMVcCLYPxb966hmIxJ5yaGyNe0p2D1i3GRYZzkHSni6zYAQri+8iL68dpe2tieqyJ7YkH2/32xOBqBvzv5UiylqzGF/hvheCIf/vrpzCspENDYG0AO/a8tfxrnGVfzQMOKtsTqexI4Lm6uRzoSaPXY+0eTotxsSi5gxVpVaxKayXdORiK6JVSo4ixGUoTuylNPHnfgccndA7F0u6Oo8Nt/dnujqVjyLpf05/Cts5YBn1v/VcoGLJjByiI76Ugvo8aH5SnuChPcU19qDcmxVqIkrrMuu8bsgrkDTRYy1b7G6Bn/7D26dbEc3wBxOWDuyOgD7ATNXsSQUyytTlkmG63neruOKp74jjUHUdlewKVHQm0uN4qGZEf18fq9GMsTOlgUXIHObEDkb5KTCnl57AZsmJdZMW6Rm1jDPR7HbQMxnN0IJGjA0kcdSVydCCRyq5Mnmh8a9g4PXbozaRQnuJiXsoA5SkuiiY7N2GLsVYcDa9S4HXBQOOw5FAPXf6zmWMz4bSfTOKNxhbURCAilwI/A+zAvcaY/znhefE/fznQD9xojNkazJi+taWY3e2JVPfE0ep66/hIuxjmpw5wdn43y9L7WZrRx5L0fnbXabVCpWYyEUh0eEh09FCS2PO253wG5uTmU9UVzyH/B8ZDXfH8qz6NPw++9f+H0+ajNNnFutwebq+onVpA9jirmGVS6VuPefpgoAnK/21q1x5F0BKBiNiBnwMXAfXAGyLymDFmz7BmlwHz/be1wC/9fwZNVXc8Xa4hlid3kpfTT0FcH/nxfeTG9uMYVuzK1wu7dchfqVnNJlDf3EgcsNQJS7MA/8rynqEYGl2JNBzvRQwk4vIGacmSIxGSyyFtWXAuH5SrWtYAVcaYagAReRC4BhieCK4B7vMfWL9RRNJEJN8YE7ST2+8774Cu2lFKTVlyzBDJMZ0sSO5887FpnVgOoWAmgkKgbtj9ek7+tD9Sm0LgbYlARG4BbvHf7RWR/YROFtAawvebimiJNVrihOiJNVrihOiJNQLj/OZoTwQS69zRnghmIhhpSvXEtaqBtMEYcw9wz3QENVEisnm0JVeRJlpijZY4IXpijZY4IXpijZY4YeqxBnMPXj1QPOx+EXDiobyBtFFKKRVEwUwEbwDzRaRURJzA9cBjJ7R5DPiwWE4HuoI5P6CUUupkQRsaMsZ4ROQ24Cms5aO/NcZUisit/ufvBp7AWjpahbV89KZgxTMFYRmSmqRoiTVa4oToiTVa4oToiTVa4oQpxhp1JSaUUkpNr5lzxI5SSqlJ0USglFKznCYCPxG5VET2i0iViHx5hOcXicgGERkUkS+GI0Z/HOPF+QER2em/vSYiK8MRpz+W8WK9xh/ndhHZLCJnRWKcw9qtFhGviFwXyvhOiGG87+m5ItLl/55uF5FvRGKc/jbn+mOsFJEXQx3jsDjG+57+x7Dv527/z0BGBMaZKiKPi8gO//c08DlXY8ysv2FNZh8CygAnsANYckKbHGA18F3gixEc5xlAuv/ry4BNERxrEm/NU60A9kVinMPaPYe1wOG6CP6engv8IxzxTTDONKwqA3P893MiNdYT2l8FPBeJcQJfBb7v/zobaAecgVxfewSWN8thGGPcwPFyGG8yxjQbY94AhsIRoF8gcb5mjOnw392ItTcjHAKJtdf4f2qBREbYTBgC48bp92ngYaA5lMGdINBYwy2QON8PPGKMqQXr9yvEMR430e/pDcADIYns7QKJ0wDJ/mKeSViJIKC62ZoILKOVuog0E43zZmB9UCMaXUCxisi7RGQf8E8gOKUVxzZunCJSCLwLuDuEcY0k0H//df7hgfUisjQ0ob1NIHEuANJF5AUR2SIiHw5ZdG8X8O+UiCQAl2J9IAi1QOK8C1iMtSl3F/AZY4yPAMye8wjGFlCpiwgQcJwich5WIgjLuDuBlw/5G/A3EXkH8G3gwmAHdoJA4vxf4EvGGK+E9zCKQGLdCsw1xvSKyOXAo1jVfUMpkDgdwGnABUA8sEFENhpjDgQ7uBNM5Hf/KuBVY0x7EOMZTSBxXgJsB84HyoFnRORlY8zJJ/acQHsElmgpdRFQnCKyArgXuMYYE65SqxP6nhpjXgLKRSQr2IGdIJA4K4AHRaQGuA74hYi8MyTRvd24sRpjuo0xvf6vnwBiIvR7Wg88aYzpM8a0Ai8B4VjYMJGf0+sJz7AQBBbnTVjDbcYYUwUcBhYFdPVwTNBE2g3r00k1UMpbEzFLR2l7O+GbLB43TmAO1k7tMyL9ewrM463J4lOBhuP3IynOE9r/nvBNFgfyPc0b9j1dA9RG4vcUawjjWX/bBGA3sCwSv6f+dqlYY+6JEfxv/0vgdv/Xuf7fp6xArq9DQwRWDkNE8oDNQArgE5HPYs3aj9vtCmWcwDeATKxPrQAeE4YKigHGei1WrakhYAB4n/H/FEdYnBEhwFivAz4hIh6s7+n1kfg9NcbsFZEngZ1Yp7ffa4zZHco4A43V3/RdwNPGmL5QxziBOL8N/F5EdmENJX3JWL2tcWmJCaWUmuV0jkAppWY5TQRKKTXLaSJQSqlZThOBUkrNcpoIlFJqltNEoNQ4RCRNRD7p/7pARP4a7piUmk66fFSpcYhICVZFz2XhjkWpYNANZUqN73+wyl9sBw4Ci40xy0TkRuCdWBt8lgE/xtr1+SFgELjcGNMuIuXAz7FKA/cDHzPG7Av1X0Kp0ejQkFLj+zJwyBizCviPE55bhlVSeQ3WWRX9xphTgA3A8Yqa9wCfNsacBnwR+EUoglYqUNojUGpqnjfG9AA9ItIFPO5/fBewQkSSsA4LemhY5dLY0Iep1Og0ESg1NYPDvvYNu+/D+v2yAZ3+3oRSEUmHhpQaXw+QPJkX+osSHhaR9wCIJWznSCs1Ek0ESo3DWGc6vCoiu4EfTuISHwBuFpEdQCWRebykmsV0+ahSSs1y2iNQSqlZThOBUkrNcpoIlFJqltNEoJRSs5wmAqWUmuU0ESil1CyniUAppWa5/w8QddTSfl9kqwAAAABJRU5ErkJggg==\n",
       "text/plain": [
        "<Figure size 432x288 with 1 Axes>"
       ]
      },
      "metadata": {
       "needs_background": "light"
-     }
+     },
+     "output_type": "display_data"
     }
    ],
-   "metadata": {}
+   "source": [
+    "# clustered distribution in blue\n",
+    "sns.histplot(\n",
+    "    pd.DataFrame(my_outcomes_clustered)[\"time\"], kde=True, stat=\"density\", linewidth=0\n",
+    ")\n",
+    "\n",
+    "# non-clustered distribution in orange\n",
+    "sns.histplot(\n",
+    "    pd.DataFrame(my_outcomes)[\"time\"],\n",
+    "    kde=True,\n",
+    "    stat=\"density\",\n",
+    "    linewidth=0,\n",
+    "    color=\"orange\",\n",
+    ")"
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "### Calculating standard errors and confidence intervals from existing outcomes"
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "Once we have the outcomes, we can compute confidence intervals and standard errors from them"
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
-   "source": [
-    "from estimagic.inference import bootstrap_from_outcomes\n",
-    "\n",
-    "my_results = bootstrap_from_outcomes(\n",
-    "    data=df, outcome=ols_fit, bootstrap_outcomes=my_outcomes, ci_method=\"t\"\n",
-    ")\n",
-    "my_results[\"summary\"]"
-   ],
+   "execution_count": 12,
+   "metadata": {},
    "outputs": [
     {
-     "output_type": "execute_result",
      "data": {
-      "text/html": [
-       "<div>\n",
-       "<style scoped>\n",
-       "    .dataframe tbody tr th:only-of-type {\n",
-       "        vertical-align: middle;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe tbody tr th {\n",
-       "        vertical-align: top;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe thead th {\n",
-       "        text-align: right;\n",
-       "    }\n",
-       "</style>\n",
-       "<table border=\"1\" class=\"dataframe\">\n",
-       "  <thead>\n",
-       "    <tr style=\"text-align: right;\">\n",
-       "      <th></th>\n",
-       "      <th>mean</th>\n",
-       "      <th>std</th>\n",
-       "      <th>lower_ci</th>\n",
-       "      <th>upper_ci</th>\n",
-       "    </tr>\n",
-       "  </thead>\n",
-       "  <tbody>\n",
-       "    <tr>\n",
-       "      <th>constant</th>\n",
-       "      <td>93.767044</td>\n",
-       "      <td>1.501594</td>\n",
-       "      <td>90.756341</td>\n",
-       "      <td>96.598375</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>time</th>\n",
-       "      <td>0.391806</td>\n",
-       "      <td>0.128111</td>\n",
-       "      <td>0.113098</td>\n",
-       "      <td>0.617271</td>\n",
-       "    </tr>\n",
-       "  </tbody>\n",
-       "</table>\n",
-       "</div>"
-      ],
       "text/plain": [
-       "               mean       std   lower_ci   upper_ci\n",
-       "constant  93.767044  1.501594  90.756341  96.598375\n",
-       "time       0.391806  0.128111   0.113098   0.617271"
+       "(constant    90.756341\n",
+       " time         0.122825\n",
+       " dtype: float64,\n",
+       " constant    96.637964\n",
+       " time         0.617271\n",
+       " dtype: float64)"
       ]
      },
+     "execution_count": 12,
      "metadata": {},
-     "execution_count": 10
+     "output_type": "execute_result"
     }
    ],
-   "metadata": {}
+   "source": [
+    "my_results = bootstrap(\n",
+    "    data=df,\n",
+    "    outcome=ols_fit,\n",
+    "    existing_outcomes=my_outcomes,\n",
+    ")\n",
+    "my_results.ci(ci_method=\"t\")"
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
-    "You can use this to calculate confidence intervals with several methods (e.g. \"bca\" and \"percentile\") without duplicated evaluations of the bootstrap outcome function."
-   ],
-   "metadata": {}
+    "You can use this to calculate confidence intervals with several methods (e.g. \"percentile\" and \"bc\") without duplicated evaluations of the bootstrap outcome function."
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "## Extending bootstrap results with more draws\n",
     "\n",
@@ -777,109 +543,69 @@
     "As an example we will take `my_outcomes` defined above and extend it with another 500 draws. \n",
     "\n",
     "Note: It is very important to use a different random seed when you calculate the additional outcomes!!!"
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [],
    "source": [
-    "more_outcomes = bootstrap(data=df, outcome=ols_fit, seed=5471, n_draws=500)[\"outcomes\"]\n",
-    "combined_outcomes = pd.concat([my_outcomes, more_outcomes])\n",
-    "combined_results = bootstrap_from_outcomes(\n",
-    "    data=df,\n",
-    "    outcome=ols_fit,\n",
-    "    bootstrap_outcomes=combined_outcomes,\n",
-    ")\n",
-    "combined_results[\"summary\"]"
-   ],
+    "more_outcomes = bootstrap(data=df, outcome=ols_fit, seed=5471, n_draws=500).outcomes\n",
+    "\n",
+    "combined_outcomes = my_outcomes + more_outcomes"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
    "outputs": [
     {
-     "output_type": "execute_result",
      "data": {
-      "text/html": [
-       "<div>\n",
-       "<style scoped>\n",
-       "    .dataframe tbody tr th:only-of-type {\n",
-       "        vertical-align: middle;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe tbody tr th {\n",
-       "        vertical-align: top;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe thead th {\n",
-       "        text-align: right;\n",
-       "    }\n",
-       "</style>\n",
-       "<table border=\"1\" class=\"dataframe\">\n",
-       "  <thead>\n",
-       "    <tr style=\"text-align: right;\">\n",
-       "      <th></th>\n",
-       "      <th>mean</th>\n",
-       "      <th>std</th>\n",
-       "      <th>lower_ci</th>\n",
-       "      <th>upper_ci</th>\n",
-       "    </tr>\n",
-       "  </thead>\n",
-       "  <tbody>\n",
-       "    <tr>\n",
-       "      <th>constant</th>\n",
-       "      <td>93.761361</td>\n",
-       "      <td>1.499120</td>\n",
-       "      <td>90.852839</td>\n",
-       "      <td>96.765793</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>time</th>\n",
-       "      <td>0.392384</td>\n",
-       "      <td>0.125592</td>\n",
-       "      <td>0.163723</td>\n",
-       "      <td>0.650190</td>\n",
-       "    </tr>\n",
-       "  </tbody>\n",
-       "</table>\n",
-       "</div>"
-      ],
       "text/plain": [
-       "               mean       std   lower_ci   upper_ci\n",
-       "constant  93.761361  1.499120  90.852839  96.765793\n",
-       "time       0.392384  0.125592   0.163723   0.650190"
+       "(constant    90.791794\n",
+       " time         0.124142\n",
+       " dtype: float64,\n",
+       " constant    96.688548\n",
+       " time         0.619296\n",
+       " dtype: float64)"
       ]
      },
+     "execution_count": 14,
      "metadata": {},
-     "execution_count": 12
+     "output_type": "execute_result"
     }
    ],
-   "metadata": {}
+   "source": [
+    "combined_results = bootstrap(\n",
+    "    data=df,\n",
+    "    outcome=ols_fit,\n",
+    "    existing_outcomes=combined_outcomes,\n",
+    ")\n",
+    "combined_results.ci(ci_method=\"t\")"
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "## Accessing the bootstrap samples"
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
-    "It is also possible to just access the bootstrap samples. You can use this for example, if you want to calculate your the bootstrap outcomes in parallel in a way that is not yet supported by estimagic (e.g. on a large cluster or supercomputer)."
-   ],
-   "metadata": {}
+    "It is also possible to just access the bootstrap samples. You may do so, for example, if you want to calculate your bootstrap outcomes in parallel in a way that is not yet supported by estimagic (e.g. on a large cluster or super-computer)."
+   ]
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
-   "source": [
-    "from estimagic.inference import get_bootstrap_samples\n",
-    "\n",
-    "my_samples = get_bootstrap_samples(data=df, seed=1234)\n",
-    "\n",
-    "my_samples[0]"
-   ],
+   "execution_count": 15,
+   "metadata": {},
    "outputs": [
     {
-     "output_type": "execute_result",
      "data": {
       "text/html": [
        "<div>\n",
@@ -1014,40 +740,44 @@
        "</div>"
       ],
       "text/plain": [
-       "    id     diet  pulse  time     kind  constant\n",
-       "47  16   no fat     89    30  walking         1\n",
-       "83  28   no fat    140    30  running         1\n",
-       "38  13  low fat     93    30  walking         1\n",
-       "53  18   no fat    101    30  walking         1\n",
-       "76  26   no fat    126    15  running         1\n",
-       "..  ..      ...    ...   ...      ...       ...\n",
-       "0    1  low fat     85     1     rest         1\n",
-       "33  12  low fat     93     1  walking         1\n",
-       "21   8   no fat     92     1     rest         1\n",
-       "71  24  low fat    120    30  running         1\n",
-       "68  23  low fat     99    30  running         1\n",
+       "    id     diet  pulse time     kind  constant\n",
+       "47  16   no fat     89   30  walking         1\n",
+       "83  28   no fat    140   30  running         1\n",
+       "38  13  low fat     93   30  walking         1\n",
+       "53  18   no fat    101   30  walking         1\n",
+       "76  26   no fat    126   15  running         1\n",
+       "..  ..      ...    ...  ...      ...       ...\n",
+       "0    1  low fat     85    1     rest         1\n",
+       "33  12  low fat     93    1  walking         1\n",
+       "21   8   no fat     92    1     rest         1\n",
+       "71  24  low fat    120   30  running         1\n",
+       "68  23  low fat     99   30  running         1\n",
        "\n",
        "[90 rows x 6 columns]"
       ]
      },
+     "execution_count": 15,
      "metadata": {},
-     "execution_count": 13
+     "output_type": "execute_result"
     }
    ],
-   "metadata": {}
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "source": [],
-   "outputs": [],
-   "metadata": {}
+   "source": [
+    "from estimagic.inference import get_bootstrap_samples\n",
+    "\n",
+    "my_samples = get_bootstrap_samples(data=df, seed=1234)\n",
+    "\n",
+    "my_samples[0]"
+   ]
   }
  ],
  "metadata": {
+  "interpreter": {
+   "hash": "5cdb9867252288f10687117449de6ad870b49795ca695c868016dc0022895cce"
+  },
   "kernelspec": {
-   "name": "python3",
-   "display_name": "Python 3.8.11 64-bit ('estimagic': conda)"
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -1059,10 +789,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.11"
-  },
-  "interpreter": {
-   "hash": "5cdb9867252288f10687117449de6ad870b49795ca695c868016dc0022895cce"
+   "version": "3.9.13"
   }
  },
  "nbformat": 4,
diff --git a/docs/source/reference_guides/index.rst b/docs/source/reference_guides/index.rst
index 9518c30dd..22b2b5923 100644
--- a/docs/source/reference_guides/index.rst
+++ b/docs/source/reference_guides/index.rst
@@ -135,10 +135,10 @@ Bootstrap
 
     .. autofunction:: bootstrap
 
-.. dropdown:: bootstrap_from_outcomes
-
-    .. autofunction:: bootstrap_from_outcomes
+.. dropdown::  BootstrapResult
 
+    .. autoclass:: BootstrapResult
+        :members:
 
 
 
diff --git a/src/estimagic/__init__.py b/src/estimagic/__init__.py
index 7d6c40d01..083afd688 100644
--- a/src/estimagic/__init__.py
+++ b/src/estimagic/__init__.py
@@ -9,7 +9,6 @@
 from estimagic.estimation.estimate_msm import MomentsResult
 from estimagic.estimation.msm_weighting import get_moments_cov
 from estimagic.inference.bootstrap import bootstrap
-from estimagic.inference.bootstrap import bootstrap_from_outcomes
 from estimagic.logging.read_log import OptimizeLogReader
 from estimagic.optimization.optimize import maximize
 from estimagic.optimization.optimize import minimize
diff --git a/src/estimagic/config.py b/src/estimagic/config.py
index dbe7caca0..37d4f46e4 100644
--- a/src/estimagic/config.py
+++ b/src/estimagic/config.py
@@ -82,13 +82,6 @@
 else:
     IS_JAX_INSTALLED = True
 
-try:
-    import seaborn  # noqa: F401
-except ImportError:
-    IS_SEABORN_INSTALLED = False
-else:
-    IS_SEABORN_INSTALLED = True
-
 
 # =================================================================================
 # Dashboard Defaults
diff --git a/src/estimagic/estimation/estimate_ml.py b/src/estimagic/estimation/estimate_ml.py
index 7a8ab7f10..f7883ec48 100644
--- a/src/estimagic/estimation/estimate_ml.py
+++ b/src/estimagic/estimation/estimate_ml.py
@@ -1,9 +1,13 @@
 from dataclasses import dataclass
+from dataclasses import field
+from functools import cached_property
 from typing import Any
+from typing import Dict
 from typing import Union
 
 import numpy as np
 import pandas as pd
+from estimagic.config import DEFAULT_SEED
 from estimagic.differentiation.derivatives import first_derivative
 from estimagic.differentiation.derivatives import second_derivative
 from estimagic.exceptions import InvalidFunctionError
@@ -14,18 +18,22 @@
 from estimagic.inference.ml_covs import cov_robust
 from estimagic.inference.ml_covs import cov_strata_robust
 from estimagic.inference.shared import calculate_ci
-from estimagic.inference.shared import calculate_inference_quantities
+from estimagic.inference.shared import calculate_estimation_summary
+from estimagic.inference.shared import calculate_free_estimates
 from estimagic.inference.shared import calculate_p_values
-from estimagic.inference.shared import check_is_optimized_and_derivative_case
-from estimagic.inference.shared import convert_flat_params_to_pytree
+from estimagic.inference.shared import calculate_summary_data_estimation
+from estimagic.inference.shared import FreeParams
 from estimagic.inference.shared import get_derivative_case
 from estimagic.inference.shared import transform_covariance
+from estimagic.inference.shared import transform_free_cov_to_cov
+from estimagic.inference.shared import transform_free_values_to_params_tree
 from estimagic.optimization.optimize import maximize
 from estimagic.optimization.optimize_result import OptimizeResult
 from estimagic.parameters.block_trees import block_tree_to_matrix
 from estimagic.parameters.block_trees import matrix_to_block_tree
 from estimagic.parameters.conversion import Converter
 from estimagic.parameters.conversion import get_converter
+from estimagic.parameters.space_conversion import InternalParams
 from estimagic.shared.check_option_dicts import check_numdiff_options
 from estimagic.shared.check_option_dicts import check_optimization_options
 from estimagic.utilities import to_pickle
@@ -128,7 +136,6 @@ def estimate_ml(
     # ==================================================================================
     # Check and process inputs
     # ==================================================================================
-
     is_optimized = optimize_options is False
 
     if not is_optimized:
@@ -144,9 +151,6 @@ def estimate_ml(
     jac_case = get_derivative_case(jacobian)
     hess_case = get_derivative_case(hessian)
 
-    check_is_optimized_and_derivative_case(is_optimized, jac_case)
-    check_is_optimized_and_derivative_case(is_optimized, hess_case)
-
     check_numdiff_options(numdiff_options, "estimate_ml")
     numdiff_options = {} if numdiff_options in (None, False) else numdiff_options
     loglike_kwargs = {} if loglike_kwargs is None else loglike_kwargs
@@ -230,8 +234,9 @@ def estimate_ml(
     # ==================================================================================
 
     if jac_case == "closed-form":
-        x = converter.params_to_internal(estimates)
-        int_jac = converter.derivative_to_internal(jacobian_eval, x)
+        int_jac = converter.derivative_to_internal(
+            jacobian_eval, internal_estimates.values
+        )
     elif jac_case == "numerical":
 
         def func(x):
@@ -325,6 +330,8 @@ def func(x):
     # create a LikelihoodResult object
     # ==================================================================================
 
+    free_estimates = calculate_free_estimates(estimates, internal_estimates)
+
     res = LikelihoodResult(
         _params=estimates,
         _converter=converter,
@@ -336,7 +343,8 @@ def func(x):
         _internal_jacobian=int_jac,
         _internal_hessian=int_hess,
         _design_info=design_info,
-        _flat_params=internal_estimates,
+        _internal_estimates=internal_estimates,
+        _free_estimates=free_estimates,
         _has_constraints=constraints not in [None, []],
     )
 
@@ -348,7 +356,8 @@ class LikelihoodResult:
     """Likelihood estimation results object."""
 
     _params: Any
-    _flat_params: Any
+    _internal_estimates: InternalParams
+    _free_estimates: FreeParams
     _converter: Converter
     _has_constraints: bool
     _optimize_result: Union[OptimizeResult, None] = None
@@ -359,6 +368,7 @@ class LikelihoodResult:
     _internal_jacobian: Union[np.ndarray, None] = None
     _internal_hessian: Union[np.ndarray, None] = None
     _design_info: Union[pd.DataFrame, None] = None
+    _cache: Dict = field(default_factory=dict)
 
     def __post_init__(self):
         if self._internal_jacobian is None and self._internal_hessian is None:
@@ -389,40 +399,27 @@ def _get_free_cov(
         seed,
     ):
         if method not in self._valid_methods:
-            raise ValueError()
-
-        internal_jac = self._internal_jacobian
-        internal_hess = self._internal_hessian
-        design_info = self._design_info
-        converter = self._converter
-        flat_params = self._flat_params
-
-        if method == "jacobian":
-            int_cov = cov_jacobian(internal_jac)
-        elif method == "hessian":
-            int_cov = cov_hessian(internal_hess)
-        elif method == "robust":
-            int_cov = cov_robust(jac=internal_jac, hess=internal_hess)
-        elif method == "cluster_robust":
-            int_cov = cov_cluster_robust(
-                jac=internal_jac, hess=internal_hess, design_info=design_info
-            )
-        elif method == "strata_robust":
-            int_cov = cov_strata_robust(
-                jac=internal_jac, hess=internal_hess, design_info=design_info
-            )
-        else:
-            raise ValueError(f"Invalid method: {method}")
-
-        np.random.seed(seed)
+            msg = f"Invalid method: {method}. Valid methods are {self._valid_methods}."
+            raise ValueError(msg)
+        args = (method, n_samples, bounds_handling, seed)
+        is_cached = args in self._cache
 
-        free_cov = transform_covariance(
-            flat_params=flat_params,
-            internal_cov=int_cov,
-            converter=converter,
-            n_samples=n_samples,
-            bounds_handling=bounds_handling,
-        )
+        if is_cached:
+            free_cov = self._cache[args]
+        else:
+            free_cov = _calculate_free_cov_ml(
+                method=method,
+                internal_estimates=self._internal_estimates,
+                converter=self._converter,
+                internal_jacobian=self._internal_jacobian,
+                internal_hessian=self._internal_hessian,
+                n_samples=n_samples,
+                design_info=self._design_info,
+                bounds_handling=bounds_handling,
+                seed=seed,
+            )
+            if seed is not None:
+                self._cache[args] = free_cov
 
         return free_cov
 
@@ -440,7 +437,6 @@ def jacobian(self):
             raise NotAvailableError(
                 f"No jacobian is available because {self._no_jacobian_reason}."
             )
-
         return self._jacobian
 
     @property
@@ -449,26 +445,25 @@ def hessian(self):
             raise NotAvailableError(
                 f"No hessian is available because {self._no_hessian_reason}."
             )
-
         return self._hessian
 
-    @property
+    @cached_property
     def _se(self):
         return self.se()
 
-    @property
+    @cached_property
     def _cov(self):
         return self.cov()
 
-    @property
+    @cached_property
     def _summary(self):
         return self.summary()
 
-    @property
+    @cached_property
     def _ci(self):
         return self.ci()
 
-    @property
+    @cached_property
     def _p_values(self):
         return self.p_values()
 
@@ -477,7 +472,7 @@ def se(
         method="jacobian",
         n_samples=10_000,
         bounds_handling="clip",
-        seed=None,
+        seed=DEFAULT_SEED,
     ):
         """Calculate standard errors.
 
@@ -513,9 +508,13 @@ def se(
         )
 
         free_se = np.sqrt(np.diagonal(free_cov))
-        out = convert_flat_params_to_pytree(free_se, self._flat_params, self._converter)
 
-        return out
+        se = transform_free_values_to_params_tree(
+            values=free_se,
+            free_params=self._free_estimates,
+            params=self._params,
+        )
+        return se
 
     def cov(
         self,
@@ -523,9 +522,9 @@ def cov(
         n_samples=10_000,
         bounds_handling="clip",
         return_type="pytree",
-        seed=None,
+        seed=DEFAULT_SEED,
     ):
-        """Calculate the variance-covariance matrix of the estimated parameters.
+        """Calculate the variance-covariance (matrix) of the estimated parameters.
 
         Args:
             method (str): One of "jacobian", "hessian", "robust", "cluster_robust",
@@ -552,8 +551,8 @@ def cov(
                 transforming constraints.
 
         Returns:
-            Any: The covariance matrix of the estimated parameters as block-pytree or
-                numpy array.
+            Any: The covariance matrix of the estimated parameters as block-pytree,
+                numpy.ndarray or pandas.DataFrame.
         """
         free_cov = self._get_free_cov(
             method=method,
@@ -561,20 +560,13 @@ def cov(
             bounds_handling=bounds_handling,
             seed=seed,
         )
-        if return_type == "array":
-            out = free_cov
-        elif return_type == "dataframe":
-            free_index = np.array(self._flat_params.names)[self._flat_params.free_mask]
-            out = pd.DataFrame(data=free_cov, columns=free_index, index=free_index)
-        elif return_type == "pytree":
-            if len(free_cov) != len(self._flat_params.values):
-                raise NotAvailableError(
-                    "Covariance matrices in block-pytree format are only available if "
-                    "there are no constraints that reduce the number of free "
-                    "parameters."
-                )
-            out = matrix_to_block_tree(free_cov, self._params, self._params)
-        return out
+        cov = transform_free_cov_to_cov(
+            free_cov=free_cov,
+            free_params=self._free_estimates,
+            params=self._params,
+            return_type=return_type,
+        )
+        return cov
 
     def summary(
         self,
@@ -582,7 +574,7 @@ def summary(
         n_samples=10_000,
         ci_level=0.95,
         bounds_handling="clip",
-        seed=None,
+        seed=DEFAULT_SEED,
     ):
         """Create a summary of estimation results.
 
@@ -611,18 +603,19 @@ def summary(
         Returns:
             Any: The estimation summary as pytree of DataFrames.
         """
-        free_cov = self._get_free_cov(
+        summary_data = calculate_summary_data_estimation(
+            self,
+            free_estimates=self._free_estimates,
             method=method,
+            ci_level=ci_level,
             n_samples=n_samples,
             bounds_handling=bounds_handling,
             seed=seed,
         )
-
-        summary = calculate_inference_quantities(
-            estimates=self._params,
-            internal_estimates=self._flat_params,
-            free_cov=free_cov,
-            ci_level=ci_level,
+        summary = calculate_estimation_summary(
+            summary_data=summary_data,
+            names=self._free_estimates.all_names,
+            free_names=self._free_estimates.free_names,
         )
         return summary
 
@@ -632,7 +625,7 @@ def ci(
         n_samples=10_000,
         ci_level=0.95,
         bounds_handling="clip",
-        seed=None,
+        seed=DEFAULT_SEED,
     ):
         """Calculate confidence intervals.
 
@@ -671,18 +664,18 @@ def ci(
             seed=seed,
         )
 
-        free_values = self._flat_params.values[self._flat_params.free_mask]
-        free_se = np.sqrt(np.diagonal(free_cov))
-
-        free_lower, free_upper = calculate_ci(free_values, free_se, ci_level)
-
-        lower = convert_flat_params_to_pytree(
-            free_lower, self._flat_params, self._converter
-        )
-        upper = convert_flat_params_to_pytree(
-            free_upper, self._flat_params, self._converter
+        free_lower, free_upper = calculate_ci(
+            free_values=self._free_estimates.values,
+            free_standard_errors=np.sqrt(np.diagonal(free_cov)),
+            ci_level=ci_level,
         )
 
+        lower, upper = (
+            transform_free_values_to_params_tree(
+                values, free_params=self._free_estimates, params=self._params
+            )
+            for values in (free_lower, free_upper)
+        )
         return lower, upper
 
     def p_values(
@@ -690,7 +683,7 @@ def p_values(
         method="jacobian",
         n_samples=10_000,
         bounds_handling="clip",
-        seed=None,
+        seed=DEFAULT_SEED,
     ):
         """Calculate p-values.
 
@@ -727,15 +720,15 @@ def p_values(
             seed=seed,
         )
 
-        free_values = self._flat_params.values[self._flat_params.free_mask]
-        free_se = np.sqrt(np.diagonal(free_cov))
-        free_p_values = calculate_p_values(free_values, free_se)
-
-        out = convert_flat_params_to_pytree(
-            free_p_values, self._flat_params, self._converter
+        free_p_values = calculate_p_values(
+            free_values=self._free_estimates.values,
+            free_standard_errors=np.sqrt(np.diagonal(free_cov)),
         )
 
-        return out
+        p_values = transform_free_values_to_params_tree(
+            free_p_values, free_params=self._free_estimates, params=self._params
+        )
+        return p_values
 
     def to_pickle(self, path):
         """Save the LikelihoodResult object to pickle.
@@ -744,3 +737,41 @@ def to_pickle(self, path):
             path (str, pathlib.Path): A str or pathlib.path ending in .pkl or .pickle.
         """
         to_pickle(self, path=path)
+
+
+def _calculate_free_cov_ml(
+    method,
+    internal_estimates,
+    converter,
+    internal_jacobian,
+    internal_hessian,
+    n_samples,
+    design_info,
+    bounds_handling,
+    seed,
+):
+    if method == "jacobian":
+        int_cov = cov_jacobian(internal_jacobian)
+    elif method == "hessian":
+        int_cov = cov_hessian(internal_hessian)
+    elif method == "robust":
+        int_cov = cov_robust(jac=internal_jacobian, hess=internal_hessian)
+    elif method == "cluster_robust":
+        int_cov = cov_cluster_robust(
+            jac=internal_jacobian, hess=internal_hessian, design_info=design_info
+        )
+    elif method == "strata_robust":
+        int_cov = cov_strata_robust(
+            jac=internal_jacobian, hess=internal_hessian, design_info=design_info
+        )
+
+    np.random.seed(seed)
+
+    free_cov = transform_covariance(
+        internal_params=internal_estimates,
+        internal_cov=int_cov,
+        converter=converter,
+        n_samples=n_samples,
+        bounds_handling=bounds_handling,
+    )
+    return free_cov
diff --git a/src/estimagic/estimation/estimate_msm.py b/src/estimagic/estimation/estimate_msm.py
index 0025ab055..19512f3a4 100644
--- a/src/estimagic/estimation/estimate_msm.py
+++ b/src/estimagic/estimation/estimate_msm.py
@@ -2,29 +2,36 @@
 import functools
 from collections.abc import Callable
 from dataclasses import dataclass
+from dataclasses import field
+from functools import cached_property
 from typing import Any
+from typing import Dict
 from typing import Union
 
 import numpy as np
 import pandas as pd
+from estimagic.config import DEFAULT_SEED
 from estimagic.differentiation.derivatives import first_derivative
 from estimagic.estimation.msm_weighting import get_weighting_matrix
 from estimagic.exceptions import InvalidFunctionError
-from estimagic.exceptions import NotAvailableError
 from estimagic.inference.msm_covs import cov_optimal
 from estimagic.inference.msm_covs import cov_robust
 from estimagic.inference.shared import calculate_ci
-from estimagic.inference.shared import calculate_inference_quantities
+from estimagic.inference.shared import calculate_estimation_summary
+from estimagic.inference.shared import calculate_free_estimates
 from estimagic.inference.shared import calculate_p_values
-from estimagic.inference.shared import check_is_optimized_and_derivative_case
-from estimagic.inference.shared import convert_flat_params_to_pytree
+from estimagic.inference.shared import calculate_summary_data_estimation
+from estimagic.inference.shared import FreeParams
 from estimagic.inference.shared import get_derivative_case
 from estimagic.inference.shared import transform_covariance
+from estimagic.inference.shared import transform_free_cov_to_cov
+from estimagic.inference.shared import transform_free_values_to_params_tree
 from estimagic.optimization.optimize import minimize
 from estimagic.parameters.block_trees import block_tree_to_matrix
 from estimagic.parameters.block_trees import matrix_to_block_tree
 from estimagic.parameters.conversion import Converter
 from estimagic.parameters.conversion import get_converter
+from estimagic.parameters.space_conversion import InternalParams
 from estimagic.parameters.tree_registry import get_registry
 from estimagic.sensitivity.msm_sensitivity import calculate_actual_sensitivity_to_noise
 from estimagic.sensitivity.msm_sensitivity import (
@@ -165,8 +172,6 @@ def estimate_msm(
 
     jac_case = get_derivative_case(jacobian)
 
-    check_is_optimized_and_derivative_case(is_optimized, jac_case)
-
     check_numdiff_options(numdiff_options, "estimate_msm")
 
     numdiff_options = {} if numdiff_options in (None, False) else numdiff_options.copy()
@@ -322,10 +327,11 @@ def func(x):
     # Create MomentsResult
     # ==================================================================================
 
+    free_estimates = calculate_free_estimates(estimates, internal_estimates)
+
     res = MomentsResult(
         _params=estimates,
         _weights=weights,
-        _flat_params=internal_estimates,
         _converter=converter,
         _internal_weights=internal_weights,
         _internal_moments_cov=internal_moments_cov,
@@ -333,6 +339,8 @@ def func(x):
         _jacobian=jacobian_eval,
         _no_jacobian_reason=_no_jac_reason,
         _empirical_moments=empirical_moments,
+        _internal_estimates=internal_estimates,
+        _free_estimates=free_estimates,
         _has_constraints=constraints not in [None, []],
     )
     return res
@@ -450,8 +458,9 @@ class MomentsResult:
     """Method of moments estimation results object."""
 
     _params: Any
+    _internal_estimates: InternalParams
+    _free_estimates: FreeParams
     _weights: Any
-    _flat_params: Any
     _converter: Converter
     _internal_moments_cov: np.ndarray
     _internal_weights: np.ndarray
@@ -460,6 +469,33 @@ class MomentsResult:
     _has_constraints: bool
     _jacobian: Any = None
     _no_jacobian_reason: Union[str, None] = None
+    _cache: Dict = field(default_factory=dict)
+
+    def _get_free_cov(self, method, n_samples, bounds_handling, seed):
+        if method not in {"optimal", "robust"}:
+            msg = f"Invalid method {method}. method must be in {'optimal', 'robust'}"
+            raise ValueError(msg)
+        args = (method, n_samples, bounds_handling, seed)
+        is_cached = args in self._cache
+
+        if is_cached:
+            free_cov = self._cache[args]
+        else:
+            free_cov = _calculate_free_cov_msm(
+                internal_estimates=self._internal_estimates,
+                internal_jacobian=self._internal_jacobian,
+                internal_moments_cov=self._internal_moments_cov,
+                internal_weights=self._internal_weights,
+                converter=self._converter,
+                method=method,
+                n_samples=n_samples,
+                bounds_handling=bounds_handling,
+                seed=seed,
+            )
+            if seed is not None:
+                self._cache[args] = free_cov
+
+        return free_cov
 
     @property
     def params(self):
@@ -473,48 +509,23 @@ def weights(self):
     def jacobian(self):
         return self._jacobian
 
-    @property
+    @cached_property
     def _se(self):
         return self.se()
 
-    @property
+    @cached_property
     def _cov(self):
         return self.cov()
 
-    @property
+    @cached_property
     def _summary(self):
         return self.summary()
 
-    @property
+    @cached_property
     def _ci(self):
         return self.ci()
 
-    def _get_free_cov(self, method, n_samples, bounds_handling, seed):
-
-        int_jac = self._internal_jacobian
-        weights = self._internal_weights
-        converter = self._converter
-        flat_params = self._flat_params
-        moments_cov = self._internal_moments_cov
-
-        if method == "optimal":
-            int_cov = cov_optimal(int_jac, weights)
-        else:
-            int_cov = cov_robust(int_jac, weights, moments_cov)
-
-        np.random.seed(seed)
-
-        free_cov = transform_covariance(
-            flat_params=flat_params,
-            internal_cov=int_cov,
-            converter=converter,
-            n_samples=n_samples,
-            bounds_handling=bounds_handling,
-        )
-
-        return free_cov
-
-    @property
+    @cached_property
     def _p_values(self):
         return self.p_values()
 
@@ -523,7 +534,7 @@ def se(
         method="robust",
         n_samples=10_000,
         bounds_handling="clip",
-        seed=None,
+        seed=DEFAULT_SEED,
     ):
         """Calculate standard errors.
 
@@ -560,9 +571,13 @@ def se(
         )
 
         free_se = np.sqrt(np.diagonal(free_cov))
-        out = convert_flat_params_to_pytree(free_se, self._flat_params, self._converter)
 
-        return out
+        se = transform_free_values_to_params_tree(
+            values=free_se,
+            free_params=self._free_estimates,
+            params=self._params,
+        )
+        return se
 
     def cov(
         self,
@@ -570,7 +585,7 @@ def cov(
         n_samples=10_000,
         bounds_handling="clip",
         return_type="pytree",
-        seed=None,
+        seed=DEFAULT_SEED,
     ):
         """Calculate the variance-covariance matrix of the estimated parameters.
 
@@ -609,21 +624,13 @@ def cov(
             bounds_handling=bounds_handling,
             seed=seed,
         )
-
-        if return_type == "array":
-            out = free_cov
-        elif return_type == "dataframe":
-            free_index = np.array(self._flat_params.names)[self._flat_params.free_mask]
-            out = pd.DataFrame(data=free_cov, columns=free_index, index=free_index)
-        elif return_type == "pytree":
-            if len(free_cov) != len(self._flat_params.values):
-                raise NotAvailableError(
-                    "Covariance matrices in block-pytree format are only available if "
-                    "there are no constraints that reduce the number of free "
-                    "parameters."
-                )
-            out = matrix_to_block_tree(free_cov, self._params, self._params)
-        return out
+        cov = transform_free_cov_to_cov(
+            free_cov=free_cov,
+            free_params=self._free_estimates,
+            params=self._params,
+            return_type=return_type,
+        )
+        return cov
 
     def summary(
         self,
@@ -631,7 +638,7 @@ def summary(
         n_samples=10_000,
         ci_level=0.95,
         bounds_handling="clip",
-        seed=None,
+        seed=DEFAULT_SEED,
     ):
         """Create a summary of estimation results.
 
@@ -660,18 +667,19 @@ def summary(
         Returns:
             Any: The estimation summary as pytree of DataFrames.
         """
-        free_cov = self._get_free_cov(
+        summary_data = calculate_summary_data_estimation(
+            self,
+            free_estimates=self._free_estimates,
             method=method,
+            ci_level=ci_level,
             n_samples=n_samples,
             bounds_handling=bounds_handling,
             seed=seed,
         )
-
-        summary = calculate_inference_quantities(
-            estimates=self._params,
-            internal_estimates=self._flat_params,
-            free_cov=free_cov,
-            ci_level=ci_level,
+        summary = calculate_estimation_summary(
+            summary_data=summary_data,
+            names=self._free_estimates.all_names,
+            free_names=self._free_estimates.free_names,
         )
         return summary
 
@@ -681,7 +689,7 @@ def ci(
         n_samples=10_000,
         ci_level=0.95,
         bounds_handling="clip",
-        seed=None,
+        seed=DEFAULT_SEED,
     ):
         """Calculate confidence intervals.
 
@@ -721,17 +729,18 @@ def ci(
             seed=seed,
         )
 
-        free_values = self._flat_params.values[self._flat_params.free_mask]
-        free_se = np.sqrt(np.diagonal(free_cov))
-        free_lower, free_upper = calculate_ci(free_values, free_se, ci_level)
-
-        lower = convert_flat_params_to_pytree(
-            free_lower, self._flat_params, self._converter
-        )
-        upper = convert_flat_params_to_pytree(
-            free_upper, self._flat_params, self._converter
+        free_lower, free_upper = calculate_ci(
+            free_values=self._free_estimates.values,
+            free_standard_errors=np.sqrt(np.diagonal(free_cov)),
+            ci_level=ci_level,
         )
 
+        lower, upper = (
+            transform_free_values_to_params_tree(
+                values, free_params=self._free_estimates, params=self._params
+            )
+            for values in (free_lower, free_upper)
+        )
         return lower, upper
 
     def p_values(
@@ -739,7 +748,7 @@ def p_values(
         method="robust",
         n_samples=10_000,
         bounds_handling="clip",
-        seed=None,
+        seed=DEFAULT_SEED,
     ):
         """Calculate p-values.
 
@@ -774,22 +783,22 @@ def p_values(
             seed=seed,
         )
 
-        free_values = self._flat_params.values[self._flat_params.free_mask]
-        free_se = np.sqrt(np.diagonal(free_cov))
-        free_p_values = calculate_p_values(free_values, free_se)
-
-        out = convert_flat_params_to_pytree(
-            free_p_values, self._flat_params, self._converter
+        free_p_values = calculate_p_values(
+            free_values=self._free_estimates.values,
+            free_standard_errors=np.sqrt(np.diagonal(free_cov)),
         )
 
-        return out
+        p_values = transform_free_values_to_params_tree(
+            free_p_values, free_params=self._free_estimates, params=self._params
+        )
+        return p_values
 
     def sensitivity(
         self,
         kind="bias",
         n_samples=10_000,
         bounds_handling="clip",
-        seed=None,
+        seed=DEFAULT_SEED,
         return_type="pytree",
     ):
         """Calculate sensitivity measures for moments estimates.
@@ -924,7 +933,7 @@ def sensitivity(
             )
         elif return_type == "dataframe":
             registry = get_registry(extended=True)
-            row_names = self._flat_params.names
+            row_names = self._internal_estimates.names
             col_names = leaf_names(self._empirical_moments, registry=registry)
             out = pd.DataFrame(
                 data=raw,
@@ -946,3 +955,34 @@ def to_pickle(self, path):
             path (str, pathlib.Path): A str or pathlib.path ending in .pkl or .pickle.
         """
         to_pickle(self, path=path)
+
+
+def _calculate_free_cov_msm(
+    internal_estimates,
+    internal_jacobian,
+    internal_moments_cov,
+    internal_weights,
+    converter,
+    method,
+    n_samples,
+    bounds_handling,
+    seed,
+):
+
+    if method == "optimal":
+        internal_cov = cov_optimal(internal_jacobian, internal_weights)
+    else:
+        internal_cov = cov_robust(
+            internal_jacobian, internal_weights, internal_moments_cov
+        )
+
+    np.random.seed(seed)
+
+    free_cov = transform_covariance(
+        internal_params=internal_estimates,
+        internal_cov=internal_cov,
+        converter=converter,
+        n_samples=n_samples,
+        bounds_handling=bounds_handling,
+    )
+    return free_cov
diff --git a/src/estimagic/inference/__init__.py b/src/estimagic/inference/__init__.py
index 9b6beff88..6b440b43c 100644
--- a/src/estimagic/inference/__init__.py
+++ b/src/estimagic/inference/__init__.py
@@ -1,5 +1,4 @@
 from estimagic.inference.bootstrap import bootstrap  # noqa: F401
-from estimagic.inference.bootstrap import bootstrap_from_outcomes  # noqa: F401
 from estimagic.inference.bootstrap_outcomes import get_bootstrap_outcomes  # noqa: F401
 from estimagic.inference.bootstrap_samples import get_bootstrap_indices  # noqa: F401
 from estimagic.inference.bootstrap_samples import get_bootstrap_samples  # noqa: F401
diff --git a/src/estimagic/inference/bootstrap.py b/src/estimagic/inference/bootstrap.py
index f712cadac..c9d8f4037 100644
--- a/src/estimagic/inference/bootstrap.py
+++ b/src/estimagic/inference/bootstrap.py
@@ -1,13 +1,15 @@
 import functools
 from dataclasses import dataclass
+from functools import cached_property
 from typing import Any
 
 import numpy as np
 import pandas as pd
 from estimagic.batch_evaluators import joblib_batch_evaluator
-from estimagic.inference.bootstrap_ci import compute_ci
+from estimagic.inference.bootstrap_ci import calculate_ci
 from estimagic.inference.bootstrap_helpers import check_inputs
 from estimagic.inference.bootstrap_outcomes import get_bootstrap_outcomes
+from estimagic.inference.shared import calculate_estimation_summary
 from estimagic.parameters.block_trees import matrix_to_block_tree
 from estimagic.parameters.tree_registry import get_registry
 from pybaum import leaf_names
@@ -17,9 +19,10 @@
 
 
 def bootstrap(
-    data,
     outcome,
     *,
+    data=None,
+    existing_outcomes=None,
     outcome_kwargs=None,
     n_draws=1_000,
     cluster_by=None,
@@ -28,15 +31,20 @@ def bootstrap(
     error_handling="continue",
     batch_evaluator=joblib_batch_evaluator,
 ):
-    """Calculate bootstrap estimates, standard errors and confidence intervals
-    for statistic of interest in given original sample.
+    """Use the bootstrap to calculate inference quantities.
 
     Args:
-        data (pandas.DataFrame): Original dataset.
-        outcome (callable): Function of the data calculating statistic of interest.
-            Returns a general pytree (e.g. pandas Series, dict, numpy array, etc.).
+        outcome (Union[callable, Any]): Either a function that computes the statistic
+            of interest, or an evaluation of that function. If it is an evaluation,
+            existing_outcomes must be passed, and the evaluation step is skipped.
+        data (pd.DataFrame): Dataset. Default None.
+        existing_outcomes (Union[BootstrapResult, List[Any]]): Evaluations of the
+            outcome function. If None, then outcome must be callable and new outcomes
+            are generated. If not None and outcome is callable, the new evaluations are
+            appended to the existing ones. Default None.
         outcome_kwargs (dict): Additional keyword arguments for outcome.
-        n_draws (int): number of bootstrap samples to draw.
+        n_draws (int): Number of bootstrap samples to draw. If len(existing_outcomes) >=
+            n_draws, a random subset of existing_outcomes is used.
         cluster_by (str): Column name of variable to cluster by or None.
         seeds (numpy.array): Array of seeds for bootstrap samples, default is none.
         n_cores (int): number of jobs for parallelization.
@@ -51,88 +59,85 @@ def bootstrap(
         BootstrapResult: A BootstrapResult object storing information on summary
             statistics, the covariance matrix, and estimated boostrap outcomes.
     """
-    check_inputs(data, cluster_by)
-
-    if outcome_kwargs is not None:
-        outcome = functools.partial(outcome, **outcome_kwargs)
-
-    bootstrap_outcomes = get_bootstrap_outcomes(
-        data=data,
-        outcome=outcome,
-        cluster_by=cluster_by,
-        seed=seed,
-        n_draws=n_draws,
-        n_cores=n_cores,
-        error_handling=error_handling,
-        batch_evaluator=batch_evaluator,
-    )
-
-    base_outcome = outcome(data)
-    out = bootstrap_from_outcomes(base_outcome, bootstrap_outcomes)
-
-    return out
-
-
-def bootstrap_from_outcomes(base_outcome, bootstrap_outcomes):
-    """Create BootstrapResults object.
+    if existing_outcomes is None:
+        existing_outcomes = []
+    elif isinstance(existing_outcomes, BootstrapResult):
+        existing_outcomes = existing_outcomes.outcomes
+    elif not isinstance(existing_outcomes, list):
+        raise ValueError("existing_outcomes must be a list or BootstrapResult.")
+
+    n_existing = len(existing_outcomes)
+
+    if callable(outcome) and n_draws > n_existing:
+
+        check_inputs(data=data, cluster_by=cluster_by)
+
+        if outcome_kwargs is not None:
+            outcome = functools.partial(outcome, **outcome_kwargs)
+
+        new_outcomes = get_bootstrap_outcomes(
+            data=data,
+            outcome=outcome,
+            cluster_by=cluster_by,
+            seed=seed,
+            n_draws=n_draws - n_existing,
+            n_cores=n_cores,
+            error_handling=error_handling,
+            batch_evaluator=batch_evaluator,
+        )
+        base_outcome = outcome(data)
+    else:
+        new_outcomes = []
+        base_outcome = outcome
 
-    Args:
-        base_outcome (pytree): Pytree of base outcomes, i.e. the outcome
-            statistic(s) evaluated on the original data set.
-        bootstrap_outcomes (list): List of pytrees of estimated
-            bootstrap outcomes.
+    if n_draws <= n_existing:
+        np.random.seed(seed)
+        random_indices = np.random.choice(n_existing, n_draws, replace=False)
+        all_outcomes = [existing_outcomes[k] for k in random_indices]
+    else:
+        all_outcomes = existing_outcomes + new_outcomes
 
-    Returns:
-        BootstrapResult: A BootstrapResult object storing information on summary
-            statistics, the covariance matrix, and the estimated boostrap outcomes.
-    """
-    if isinstance(bootstrap_outcomes, list):
-        registry = get_registry(extended=True)
+    # ==================================================================================
+    # Process results
+    # ==================================================================================
 
-        flat_outcomes = [
-            tree_just_flatten(est, registry=registry) for est in bootstrap_outcomes
-        ]
-        internal_outcomes = np.array(flat_outcomes)
-    else:
-        raise TypeError(
-            "bootstrap_outcomes must be a list of pytrees, "
-            f"not {type(bootstrap_outcomes)}."
-        )
+    registry = get_registry(extended=True)
+    flat_outcomes = [
+        tree_just_flatten(_outcome, registry=registry) for _outcome in all_outcomes
+    ]
+    internal_outcomes = np.array(flat_outcomes)
 
-    out = BootstrapResult(
+    result = BootstrapResult(
         _base_outcome=base_outcome,
         _internal_outcomes=internal_outcomes,
+        _internal_cov=np.cov(internal_outcomes, rowvar=False),
     )
-
-    return out
+    return result
 
 
 @dataclass
 class BootstrapResult:
     _base_outcome: Any
     _internal_outcomes: np.ndarray
+    _internal_cov: np.ndarray
 
-    @property
-    def _outcomes(self):
-        return self.outcomes()
-
-    @property
+    @cached_property
     def _se(self):
         return self.se()
 
-    @property
+    @cached_property
     def _cov(self):
         return self.cov()
 
-    @property
+    @cached_property
     def _ci(self):
         return self.ci()
 
-    @property
+    @cached_property
     def _p_values(self):
         return self.p_values()
 
-    @property
+    @cached_property
     def _summary(self):
         return self.summary()
 
@@ -146,11 +151,12 @@ def base_outcome(self):
         """
         return self._base_outcome
 
+    @cached_property
     def outcomes(self):
         """Returns the estimated bootstrap outcomes.
 
         Returns:
-            Any: The boostrap outcomes as a list of pytrees.
+            List[Any]: The boostrap outcomes as a list of pytrees.
         """
         registry = get_registry(extended=True)
         _, treedef = tree_flatten(self._base_outcome, registry=registry)
@@ -159,25 +165,23 @@ def outcomes(self):
             tree_unflatten(treedef, out, registry=registry)
             for out in self._internal_outcomes
         ]
-
         return outcomes
 
     def se(self):
         """Calculate standard errors.
 
         Returns:
-            list: A list of pytrees containing standard errors for the bootstrapped
-                statistic.
+            Any: The standard errors of the estimated parameters as a block-pytree,
+                numpy.ndarray, or pandas.DataFrame.
         """
-        cov = np.cov(self._internal_outcomes, rowvar=False)
+        cov = self._internal_cov
         se = np.sqrt(np.diagonal(cov))
 
         registry = get_registry(extended=True)
         _, treedef = tree_flatten(self._base_outcome, registry=registry)
 
-        out = tree_unflatten(treedef, se, registry=registry)
-
-        return out
+        se = tree_unflatten(treedef, se, registry=registry)
+        return se
 
     def cov(self, return_type="pytree"):
         """Calculate the variance-covariance matrix of the estimated parameters.
@@ -191,26 +195,22 @@ def cov(self, return_type="pytree"):
 
         Returns:
             Any: The covariance matrix of the estimated parameters as a block-pytree,
-                numpy array, or pandas DataFrame.
+                numpy.ndarray, or pandas.DataFrame.
         """
-        cov = np.cov(self._internal_outcomes, rowvar=False)
+        cov = self._internal_cov
 
-        if return_type == "array":
-            out = cov
-        elif return_type == "dataframe":
+        if return_type == "dataframe":
             registry = get_registry(extended=True)
-            leafnames = leaf_names(self._base_outcome, registry=registry)
-            free_index = np.array(leafnames)
-            out = pd.DataFrame(data=cov, columns=free_index, index=free_index)
+            names = np.array(leaf_names(self._base_outcome, registry=registry))
+            cov = pd.DataFrame(cov, columns=names, index=names)
         elif return_type == "pytree":
-            out = matrix_to_block_tree(cov, self._base_outcome, self._base_outcome)
-        else:
-            raise TypeError(
+            cov = matrix_to_block_tree(cov, self._base_outcome, self._base_outcome)
+        elif return_type != "array":
+            raise ValueError(
                 "return_type must be one of pytree, array, or dataframe, "
                 f"not {return_type}."
             )
-
-        return out
+        return cov
 
     def ci(self, ci_method="percentile", ci_level=0.95):
         """Calculate confidence intervals.
@@ -230,13 +230,12 @@ def ci(self, ci_method="percentile", ci_level=0.95):
         registry = get_registry(extended=True)
         base_outcome_flat, treedef = tree_flatten(self._base_outcome, registry=registry)
 
-        lower_flat, upper_flat = compute_ci(
+        lower_flat, upper_flat = calculate_ci(
             base_outcome_flat, self._internal_outcomes, ci_method, ci_level
         )
 
         lower = tree_unflatten(treedef, lower_flat, registry=registry)
         upper = tree_unflatten(treedef, upper_flat, registry=registry)
-
         return lower, upper
 
     def p_values(self):
@@ -246,7 +245,8 @@ def p_values(self):
             Any: A pytree with the same structure as base_outcome containing p-values
                 for the parameter estimates.
         """
-        raise NotImplementedError("Bootstrapped p-values are not implemented yet.")
+        msg = "Bootstrap p_values are not yet implemented."
+        raise NotImplementedError(msg)
 
     def summary(self, ci_method="percentile", ci_level=0.95):
         """Create a summary of bootstrap results.
@@ -262,4 +262,26 @@ def summary(self, ci_method="percentile", ci_level=0.95):
                 on the mean, standard errors, as well as the confidence intervals.
                 Soon this will be a pytree.
         """
-        raise NotImplementedError("summary is not implemented yet.")
+        registry = get_registry(extended=True)
+        names = leaf_names(self.base_outcome, registry=registry)
+        summary_data = _calulcate_summary_data_bootstrap(
+            self, ci_method=ci_method, ci_level=ci_level
+        )
+        summary = calculate_estimation_summary(
+            summary_data=summary_data,
+            names=names,
+            free_names=names,
+        )
+        return summary
+
+
+def _calulcate_summary_data_bootstrap(bootstrap_result, ci_method, ci_level):
+    lower, upper = bootstrap_result.ci(ci_method=ci_method, ci_level=ci_level)
+    summary_data = {
+        "value": bootstrap_result.base_outcome,
+        "standard_error": bootstrap_result.se(),
+        "ci_lower": lower,
+        "ci_upper": upper,
+        "p_value": np.full(len(lower), np.nan),  # p-values are not implemented yet
+    }
+    return summary_data
diff --git a/src/estimagic/inference/bootstrap_ci.py b/src/estimagic/inference/bootstrap_ci.py
index a8a4f4da4..2a107f7bd 100644
--- a/src/estimagic/inference/bootstrap_ci.py
+++ b/src/estimagic/inference/bootstrap_ci.py
@@ -3,7 +3,7 @@
 from scipy.stats import norm
 
 
-def compute_ci(
+def calculate_ci(
     base_outcome_flat,
     estimates,
     ci_method="percentile",
@@ -31,7 +31,7 @@ def compute_ci(
         np.ndarray: 1d array of the upper confidence interval, where the k'th entry
             contains the upper confidence interval for the k'th parameter.
     """
-    check_inputs(ci_method=ci_method, ci_level=ci_level)
+    check_inputs(ci_method=ci_method, ci_level=ci_level, skipdata=True)
 
     alpha = 1 - ci_level
 
diff --git a/src/estimagic/inference/bootstrap_helpers.py b/src/estimagic/inference/bootstrap_helpers.py
index 06a20c86b..e0a761a2f 100644
--- a/src/estimagic/inference/bootstrap_helpers.py
+++ b/src/estimagic/inference/bootstrap_helpers.py
@@ -1,7 +1,9 @@
 import pandas as pd
 
 
-def check_inputs(data=None, cluster_by=None, ci_method="percentile", ci_level=0.95):
+def check_inputs(
+    data=None, cluster_by=None, ci_method="percentile", ci_level=0.95, skipdata=False
+):
     """Check validity of inputs.
 
     Args:
@@ -9,24 +11,26 @@ def check_inputs(data=None, cluster_by=None, ci_method="percentile", ci_level=0.
         cluster_by (str): column name of variable to cluster by.
         ci_method (str): method of choice for confidence interval computation.
         alpha (float): significance level of choice.
-    """
+        skipdata (bool): Whether to skip all checks on the data argument.
 
+    """
     ci_method_list = ["percentile", "bc", "t", "normal", "basic"]
 
-    if not isinstance(data, pd.DataFrame) and data is not None:
-        raise ValueError("Input 'data' must be DataFrame.")
-
-    elif (cluster_by is not None) and (cluster_by not in data.columns.tolist()):
-        raise ValueError(
-            "Input 'cluster_by' must be None or a column name of DataFrame."
+    if not skipdata:
+        if data is None:
+            raise ValueError("Data cannot be None if outcome is callable.")
+        elif not isinstance(data, pd.DataFrame):
+            raise TypeError("Data must be a pandas.DataFrame.")
+        elif (cluster_by is not None) and (cluster_by not in data.columns.tolist()):
+            raise ValueError(
+                "Input 'cluster_by' must be None or a column name of 'data'."
+            )
+
+    if ci_method not in ci_method_list:
+        msg = (
+            "ci_method must be 'percentile', 'bc', 't', 'basic' or 'normal', "
+            f"'{ci_method}' was supplied"
         )
-
-    elif ci_method not in ci_method_list:
-        raise ValueError(
-            "ci_method must be 'percentile', 'bc',"
-            " 't', 'basic' or 'normal', '{method}'"
-            " was supplied".format(method=ci_method)
-        )
-
-    elif ci_level > 1 or ci_level < 0:
+        raise ValueError(msg)
+    if ci_level > 1 or ci_level < 0:
         raise ValueError("Input 'ci_level' must be in [0,1].")
diff --git a/src/estimagic/inference/shared.py b/src/estimagic/inference/shared.py
index c9a337dcc..fd3b3210f 100644
--- a/src/estimagic/inference/shared.py
+++ b/src/estimagic/inference/shared.py
@@ -1,13 +1,16 @@
+from typing import NamedTuple
+
 import numpy as np
 import pandas as pd
 import scipy
+from estimagic.parameters.block_trees import matrix_to_block_tree
 from estimagic.parameters.tree_registry import get_registry
 from pybaum import tree_just_flatten
 from pybaum import tree_unflatten
 
 
 def transform_covariance(
-    flat_params,
+    internal_params,
     internal_cov,
     converter,
     n_samples,
@@ -16,8 +19,15 @@ def transform_covariance(
     """Transform the internal covariance matrix to an external one, given constraints.
 
     Args:
-        flat_params (FlatParams): NamedTuple with internal parameter values and names,
-            lower_bounds and upper_bounds, and free_mask.
+        internal_params (InternalParams): NamedTuple with entries:
+            - value (np.ndarray): Internal parameter values.
+            - lower_bounds (np.ndarray): Lower bounds on the internal params.
+            - upper_bounds (np.ndarray): Upper bounds on the internal params.
+            - soft_lower_bounds (np.ndarray): Soft lower bounds on the internal params.
+            - soft_upper_bounds (np.ndarray): Soft upper bounds on the internal params.
+            - name (list): List of names of the external parameters.
+            - free_mask (np.ndarray): Boolean mask representing which external parameter
+              is free.
         internal_cov (np.ndarray or pandas.DataFrame) with a covariance matrix of the
             internal parameter vector. For background information about internal and
             external params see :ref:`implementation_of_constraints`.
@@ -38,18 +48,15 @@ def transform_covariance(
             the index is a subset of params.index. The columns are the same as the
             index.
     """
-    if isinstance(internal_cov, pd.DataFrame):
-        internal_cov = internal_cov.to_numpy()
-
     if converter.has_transforming_constraints:
         _from_internal = converter.params_from_internal
 
-        is_free = flat_params.free_mask
-        lower_bounds = flat_params.lower_bounds
-        upper_bounds = flat_params.upper_bounds
+        is_free = internal_params.free_mask
+        lower_bounds = internal_params.lower_bounds
+        upper_bounds = internal_params.upper_bounds
 
         sample = np.random.multivariate_normal(
-            mean=flat_params.values,
+            mean=internal_params.values,
             cov=internal_cov,
             size=n_samples,
         )
@@ -77,18 +84,51 @@ def transform_covariance(
     return free_cov
 
 
-def calculate_inference_quantities(estimates, internal_estimates, free_cov, ci_level):
-    """Add standard errors, pvalues and confidence intervals to params.
+def calculate_summary_data_estimation(
+    estimation_result,
+    free_estimates,
+    ci_level,
+    method,
+    n_samples,
+    bounds_handling,
+    seed,
+):
+    se = estimation_result.se(
+        method=method, n_samples=n_samples, bounds_handling=bounds_handling, seed=seed
+    )
+    lower, upper = estimation_result.ci(
+        method=method,
+        n_samples=n_samples,
+        ci_level=ci_level,
+        bounds_handling=bounds_handling,
+        seed=seed,
+    )
+    p_values = estimation_result.p_values(
+        method=method, n_samples=n_samples, bounds_handling=bounds_handling, seed=seed
+    )
+    summary_data = {
+        "value": estimation_result.params,
+        "standard_error": se,
+        "ci_lower": lower,
+        "ci_upper": upper,
+        "p_value": p_values,
+        "free": free_estimates.free_mask,
+    }
+    return summary_data
+
+
+def calculate_estimation_summary(
+    summary_data,
+    names,
+    free_names,
+):
+    """Create estimation summary using pre-calculated results.
 
-    Args
-        params (pytree): The input parameter pytree.
-        internal_estimates (FlatParams): NamedTuple with internal estimated parameter
-            values and names, lower_bounds and upper_bounds, and free_mask.
-        free_cov (pd.DataFrame): Quadratic DataFrame containing the covariance matrix
-            of the free parameters. If parameters were fixed (explicitly or by other
-            constraints) the index is a subset of params.index. The columns are the same
-            as the index.
-        ci_level (float): Confidence level for the calculation of confidence intervals.
+    Args:
+        summary_data (dict): Dictionary with entries ['params', 'p_value', 'ci_lower',
+        'ci_upper', 'standard_error'].
+        names (List[str]): List of parameter names, corresponding to result_object.
+        free_names (List[str]): List of parameter names for free parameters.
 
     Returns:
         pytree: A pytree with the same structure as params. Each leaf in the params
@@ -99,34 +139,20 @@ def calculate_inference_quantities(estimates, internal_estimates, free_cov, ci_l
             reproduced for convenience.
 
     """
-    if not isinstance(free_cov, pd.DataFrame):
-        free_index = np.array(internal_estimates.names)[internal_estimates.free_mask]
-        free_cov = pd.DataFrame(data=free_cov, columns=free_index, index=free_index)
     # ==================================================================================
-    # Construct summary data frame for flat estimates
+    # Flatten summary and construct data frame for flat estimates
     # ==================================================================================
-    registry = get_registry(extended=True)
-
-    df = pd.DataFrame(index=internal_estimates.names)
-    df["value"] = tree_just_flatten(estimates, registry=registry)
-    df["free"] = internal_estimates.free_mask
-    df.loc[free_cov.index, "standard_error"] = np.sqrt(np.diag(free_cov))
 
-    df["p_value"] = calculate_p_values(
-        df["value"].to_numpy(),
-        df["standard_error"].to_numpy(),
-    )
+    registry = get_registry(extended=True)
+    flat_data = {
+        key: tree_just_flatten(val, registry=registry)
+        for key, val in summary_data.items()
+    }
 
-    lower, upper = calculate_ci(
-        df["value"].to_numpy(),
-        df["standard_error"].to_numpy(),
-        ci_level=ci_level,
-    )
-    df["ci_lower"] = lower
-    df["ci_upper"] = upper
+    df = pd.DataFrame(flat_data, index=names)
 
-    df.loc[free_cov.index, "stars"] = pd.cut(
-        df.loc[free_cov.index, "p_value"],
+    df.loc[free_names, "stars"] = pd.cut(
+        df.loc[free_names, "p_value"],
         bins=[-1, 0.01, 0.05, 0.1, 2],
         labels=["***", "**", "*", ""],
     )
@@ -136,9 +162,9 @@ def calculate_inference_quantities(estimates, internal_estimates, free_cov, ci_l
     # ==================================================================================
 
     # create tree with values corresponding to indices of df
-    indices = tree_unflatten(estimates, internal_estimates.names, registry=registry)
+    indices = tree_unflatten(summary_data["value"], names, registry=registry)
 
-    estimates_flat = tree_just_flatten(estimates)
+    estimates_flat = tree_just_flatten(summary_data["value"])
     indices_flat = tree_just_flatten(indices)
 
     # use index chunks in indices_flat to access the corresponding sub data frame of df,
@@ -151,20 +177,20 @@ def calculate_inference_quantities(estimates, internal_estimates, free_cov, ci_l
             loc = [index_leaf]
             index = [0]
         elif isinstance(params_leaf, pd.DataFrame) and "value" in params_leaf:
-            loc = index_leaf["value"].values.flatten()
+            loc = index_leaf["value"].to_numpy().flatten()
             index = params_leaf.index
         elif isinstance(params_leaf, pd.DataFrame):
-            loc = index_leaf.values.flatten()
+            loc = index_leaf.to_numpy().flatten()
             # use product of existing index and columns for regular pd.DataFrame
             index = pd.MultiIndex.from_tuples(
                 [
-                    (*row, col)
+                    (*row, col) if isinstance(row, tuple) else (row, col)
                     for row in params_leaf.index
                     for col in params_leaf.columns
                 ]
             )
         elif isinstance(params_leaf, pd.Series):
-            loc = index_leaf.values.flatten()
+            loc = index_leaf.to_numpy().flatten()
             index = params_leaf.index
         else:
             # array case (numpy or jax)
@@ -181,7 +207,7 @@ def calculate_inference_quantities(estimates, internal_estimates, free_cov, ci_l
 
         summary_flat.append(df_chunk)
 
-    summary = tree_unflatten(estimates, summary_flat)
+    summary = tree_unflatten(summary_data["value"], summary_flat)
     return summary
 
 
@@ -240,7 +266,7 @@ def _to_numpy(df_or_array, name):
     elif isinstance(df_or_array, np.ndarray):
         arr = df_or_array
     else:
-        raise ValueError(
+        raise TypeError(
             f"{name} must be a DataFrame or numpy array, not {type(df_or_array)}."
         )
     return arr
@@ -259,9 +285,7 @@ def _check_names_coincide(name_dict):
 
 def get_derivative_case(derivative):
     """Determine which kind of derivative should be used."""
-    if isinstance(derivative, (pd.DataFrame, np.ndarray)):
-        case = "pre-calculated"
-    elif callable(derivative):
+    if callable(derivative):
         case = "closed-form"
     elif derivative is False:
         case = "skip"
@@ -270,33 +294,66 @@ def get_derivative_case(derivative):
     return case
 
 
-def check_is_optimized_and_derivative_case(is_minimized, derivative_case):
-    if (not is_minimized) and derivative_case == "pre-calculated":
-        raise ValueError(
-            "Providing a pre-calculated derivative is only possible if the "
-            "optimization was done outside of the estimate_function, i.e. if "
-            "optimize_options=False."
-        )
-
-
-def calculate_ci(flat_values, flat_standard_errors, ci_level):
-
+def calculate_ci(free_values, free_standard_errors, ci_level):
     alpha = 1 - ci_level
     scale = scipy.stats.norm.ppf(1 - alpha / 2)
-    lower = flat_values - scale * flat_standard_errors
-    upper = flat_values + scale * flat_standard_errors
+    lower = free_values - scale * free_standard_errors
+    upper = free_values + scale * free_standard_errors
     return lower, upper
 
 
-def calculate_p_values(flat_values, flat_standard_error):
-    tvalues = flat_values / np.clip(flat_standard_error, 1e-300, np.inf)
+def calculate_p_values(free_values, free_standard_errors):
+    tvalues = free_values / np.clip(free_standard_errors, 1e-300, np.inf)
     pvalues = 2 * scipy.stats.norm.sf(np.abs(tvalues))
     return pvalues
 
 
-def convert_flat_params_to_pytree(free_statistic, flat_params, converter):
-    helper = np.full(len(flat_params.values), np.nan)
-    helper[flat_params.free_mask] = free_statistic
-    out = converter.params_from_internal(helper)
+def calculate_free_estimates(estimates, internal_estimates):
+    mask = internal_estimates.free_mask
+    names = internal_estimates.names
+
+    registry = get_registry(extended=True)
+    external_flat = np.array(tree_just_flatten(estimates, registry=registry))
+
+    free_estimates = FreeParams(
+        values=external_flat[mask],
+        free_mask=mask,
+        all_names=names,
+        free_names=np.array(names)[mask].tolist(),
+    )
+    return free_estimates
+
+
+def transform_free_cov_to_cov(free_cov, free_params, params, return_type):
+    """Fill non-free values and project to params block-tree."""
+    mask = free_params.free_mask
+    cov = np.full((len(mask), len(mask)), np.nan)
+    cov[np.ix_(mask, mask)] = free_cov
+    if return_type == "dataframe":
+        names = free_params.all_names
+        cov = pd.DataFrame(cov, columns=names, index=names)
+    elif return_type == "pytree":
+        cov = matrix_to_block_tree(cov, params, params)
+    elif return_type != "array":
+        raise ValueError(
+            "return_type must be one of pytree, array, or dataframe, "
+            f"not {return_type}."
+        )
+    return cov
+
+
+def transform_free_values_to_params_tree(values, free_params, params):
+    """Fill non-free values and project to params tree structure."""
+    mask = free_params.free_mask
+    flat = np.full(len(mask), np.nan)
+    flat[np.ix_(mask)] = values
+    registry = get_registry(extended=True)
+    pytree = tree_unflatten(params, flat, registry=registry)
+    return pytree
+
 
-    return out
+class FreeParams(NamedTuple):
+    values: np.ndarray  # free external parameter values
+    free_mask: np.ndarray  # boolean mask to filter free params from external params
+    free_names: list  # names of free external parameters
+    all_names: list  # names of all external parameters
diff --git a/src/estimagic/optimization/convergence_report.py b/src/estimagic/optimization/convergence_report.py
index 374d98c95..4685066ed 100644
--- a/src/estimagic/optimization/convergence_report.py
+++ b/src/estimagic/optimization/convergence_report.py
@@ -2,12 +2,11 @@
 from estimagic.optimization.history_tools import get_history_arrays
 
 
-def get_convergence_report(history, direction, converter=None):
+def get_convergence_report(history, direction):
 
     history_arrs = get_history_arrays(
         history=history,
         direction=direction,
-        converter=converter,
     )
 
     critvals = history_arrs["criterion"][history_arrs["is_accepted"]]
diff --git a/src/estimagic/optimization/history_tools.py b/src/estimagic/optimization/history_tools.py
index 2abd8b4af..3d0c31eef 100644
--- a/src/estimagic/optimization/history_tools.py
+++ b/src/estimagic/optimization/history_tools.py
@@ -5,10 +5,13 @@
 from pybaum import tree_just_flatten
 
 
-def get_history_arrays(history, direction, converter=None):
-
-    if converter is not None:
-        to_internal = converter.params_to_internal
+def get_history_arrays(history, direction):
+    parhist = history["params"]
+    is_flat = (
+        len(parhist) > 0 and isinstance(parhist[0], np.ndarray) and parhist[0].ndim == 1
+    )
+    if is_flat:
+        to_internal = lambda x: x.tolist()
     else:
         registry = get_registry(extended=True)
         to_internal = partial(tree_just_flatten, registry=registry)
diff --git a/src/estimagic/optimization/process_results.py b/src/estimagic/optimization/process_results.py
index df3b280bd..01ff982ce 100644
--- a/src/estimagic/optimization/process_results.py
+++ b/src/estimagic/optimization/process_results.py
@@ -45,7 +45,6 @@ def process_internal_optimizer_result(
             conv_report = get_convergence_report(
                 history=hist,
                 direction=fixed_kwargs["direction"],
-                converter=converter,
             )
 
             res.convergence_report = conv_report
@@ -101,7 +100,6 @@ def _process_one_result(res, converter, primary_key, fixed_kwargs, skip_checks):
         conv_report = get_convergence_report(
             history=res["history"],
             direction=fixed_kwargs["direction"],
-            converter=converter,
         )
 
     else:
diff --git a/src/estimagic/parameters/consolidate_constraints.py b/src/estimagic/parameters/consolidate_constraints.py
index faa2a71e0..b5d9b0b1b 100644
--- a/src/estimagic/parameters/consolidate_constraints.py
+++ b/src/estimagic/parameters/consolidate_constraints.py
@@ -289,7 +289,7 @@ def simplify_covariance_and_sdcorr_constraints(
 
         if uncorrelated:
             lower[diag_indices] = np.maximum(0, lower[diag_indices])
-        elif dim <= 2:
+        elif dim <= 2 and constr["type"] == "sdcorr":
             lower[diag_indices] = np.maximum(0, lower[diag_indices])
             lower[off_indices] = -1
             upper[off_indices] = 1
diff --git a/src/estimagic/parameters/constraint_tools.py b/src/estimagic/parameters/constraint_tools.py
index 6713ea9d1..9eaa33704 100644
--- a/src/estimagic/parameters/constraint_tools.py
+++ b/src/estimagic/parameters/constraint_tools.py
@@ -15,7 +15,7 @@ def count_free_params(params, constraints=None, lower_bounds=None, upper_bounds=
         int: Number of (free) parameters
 
     """
-    _, flat_params = get_converter(
+    _, internal_params = get_converter(
         params=params,
         constraints=constraints,
         lower_bounds=lower_bounds,
@@ -26,7 +26,7 @@ def count_free_params(params, constraints=None, lower_bounds=None, upper_bounds=
         scaling_options={},
     )
 
-    return int(flat_params.free_mask.sum())
+    return int(internal_params.free_mask.sum())
 
 
 def check_constraints(params, constraints, lower_bounds=None, upper_bounds=None):
diff --git a/src/estimagic/parameters/conversion.py b/src/estimagic/parameters/conversion.py
index b82dd47f8..ca19bcb4d 100644
--- a/src/estimagic/parameters/conversion.py
+++ b/src/estimagic/parameters/conversion.py
@@ -5,7 +5,7 @@
 from estimagic.parameters.process_selectors import process_selectors
 from estimagic.parameters.scale_conversion import get_scale_converter
 from estimagic.parameters.space_conversion import get_space_converter
-from estimagic.parameters.tree_conversion import FlatParams
+from estimagic.parameters.space_conversion import InternalParams
 from estimagic.parameters.tree_conversion import get_tree_converter
 
 
@@ -51,12 +51,13 @@ def get_converter(
             func at params. Used for consistency checks.
         soft_lower_bounds (pytree): As lower_bounds
         soft_upper_bounds (pytree): As upper_bounds
-        add_soft_bounds (bool): Whether soft bounds should be added to the flat_params
+        add_soft_bounds (bool): Whether soft bounds should be added to the
+            internal_params
 
     Returns:
         Converter: NamedTuple with methods to convert between internal and external
             parameters, derivatives and function outputs.
-        FlatParams: NamedTuple with internal parameter values, lower_bounds and
+        InternalParams: NamedTuple with internal parameter values, lower_bounds and
             upper_bounds.
 
     """
@@ -77,7 +78,7 @@ def get_converter(
             primary_key=primary_key,
         )
 
-    tree_converter, flat_params = get_tree_converter(
+    tree_converter, internal_params = get_tree_converter(
         params=params,
         lower_bounds=lower_bounds,
         upper_bounds=upper_bounds,
@@ -93,15 +94,15 @@ def get_converter(
         constraints=constraints,
         params=params,
         tree_converter=tree_converter,
-        param_names=flat_params.names,
+        param_names=internal_params.names,
     )
 
     space_converter, internal_params = get_space_converter(
-        flat_params=flat_params, flat_constraints=flat_constraints
+        internal_params=internal_params, internal_constraints=flat_constraints
     )
 
     scale_converter, scaled_params = get_scale_converter(
-        flat_params=internal_params,
+        internal_params=internal_params,
         scaling=scaling,
         scaling_options=scaling_options,
     )
@@ -148,7 +149,7 @@ def _derivative_to_internal(derivative_eval, x, jac_is_flat=False):
     def _func_to_internal(func_eval):
         return tree_converter.func_flatten(func_eval)
 
-    flat_params = scaled_params._replace(free_mask=internal_params.free_mask)
+    internal_params = scaled_params._replace(free_mask=internal_params.free_mask)
 
     converter = Converter(
         params_to_internal=_params_to_internal,
@@ -158,7 +159,7 @@ def _func_to_internal(func_eval):
         has_transforming_constraints=space_converter.has_transforming_constraints,
     )
 
-    return converter, flat_params
+    return converter, internal_params
 
 
 class Converter(NamedTuple):
@@ -237,14 +238,14 @@ def _fast_derivative_to_internal(derivative_eval, x, jac_is_flat=True):
     else:
         upper_bounds = upper_bounds.astype(float)
 
-    flat_params = FlatParams(
+    internal_params = InternalParams(
         values=params.astype(float),
         lower_bounds=lower_bounds,
         upper_bounds=upper_bounds,
         free_mask=np.full(len(params), True),
         names=[str(i) for i in range(len(params))],
     )
-    return converter, flat_params
+    return converter, internal_params
 
 
 def _is_fast_path(
diff --git a/src/estimagic/parameters/scale_conversion.py b/src/estimagic/parameters/scale_conversion.py
index c92314662..fa0c3b0b9 100644
--- a/src/estimagic/parameters/scale_conversion.py
+++ b/src/estimagic/parameters/scale_conversion.py
@@ -2,19 +2,19 @@
 from typing import NamedTuple
 
 import numpy as np
-from estimagic.parameters.tree_conversion import FlatParams
+from estimagic.parameters.space_conversion import InternalParams
 
 
 def get_scale_converter(
-    flat_params,
+    internal_params,
     scaling,
     scaling_options,
 ):
     """Get a converter between scaled and unscaled parameters.
 
     Args:
-        flat_params (FlatParams): NamedTuple of flattened and possibly reparametrized
-            but not yet scaled parameter values and bounds.
+        internal_params (InternalParams): NamedTuple of internal and possibly
+            reparametrized but not yet scaled parameter values and bounds.
         func (callable): The criterion function. Possibly used to calculate a scaling
             factor.
         scaling (bool): Whether scaling should be done.
@@ -23,20 +23,27 @@ def get_scale_converter(
     Returns:
         ScaleConverter: NamedTuple with methods to convert between scaled and unscaled
             internal parameters and derivatives.
-        FlatParams: NamedTuple of 1d numpy array with flat, internal and scaled
-            parameters and bounds.
+        InternalParams: NamedTuple with entries:
+            - value (np.ndarray): Internal parameter values.
+            - lower_bounds (np.ndarray): Lower bounds on the internal params.
+            - upper_bounds (np.ndarray): Upper bounds on the internal params.
+            - soft_lower_bounds (np.ndarray): Soft lower bounds on the internal params.
+            - soft_upper_bounds (np.ndarray): Soft upper bounds on the internal params.
+            - name (list): List of names of the external parameters.
+            - free_mask (np.ndarray): Boolean mask representing which external parameter
+              is free.
 
     """
     # fast path
     if not scaling:
-        return _fast_path_scale_converter(), flat_params
+        return _fast_path_scale_converter(), internal_params
 
     scaling_options = {} if scaling_options is None else scaling_options
     valid_keys = {"method", "clipping_value", "magnitude"}
     scaling_options = {k: v for k, v in scaling_options.items() if k in valid_keys}
 
     factor, offset = calculate_scaling_factor_and_offset(
-        flat_params=flat_params, **scaling_options
+        internal_params=internal_params, **scaling_options
     )
 
     _params_to_internal = partial(
@@ -64,21 +71,21 @@ def _derivative_from_internal(derivative):
         derivative_from_internal=_derivative_from_internal,
     )
 
-    if flat_params.soft_lower_bounds is not None:
-        _soft_lower = converter.params_to_internal(flat_params.soft_lower_bounds)
+    if internal_params.soft_lower_bounds is not None:
+        _soft_lower = converter.params_to_internal(internal_params.soft_lower_bounds)
     else:
         _soft_lower = None
 
-    if flat_params.soft_upper_bounds is not None:
-        _soft_upper = converter.params_to_internal(flat_params.soft_upper_bounds)
+    if internal_params.soft_upper_bounds is not None:
+        _soft_upper = converter.params_to_internal(internal_params.soft_upper_bounds)
     else:
         _soft_upper = None
 
-    params = FlatParams(
-        values=converter.params_to_internal(flat_params.values),
-        lower_bounds=converter.params_to_internal(flat_params.lower_bounds),
-        upper_bounds=converter.params_to_internal(flat_params.upper_bounds),
-        names=flat_params.names,
+    params = InternalParams(
+        values=converter.params_to_internal(internal_params.values),
+        lower_bounds=converter.params_to_internal(internal_params.lower_bounds),
+        upper_bounds=converter.params_to_internal(internal_params.upper_bounds),
+        names=internal_params.names,
         soft_lower_bounds=_soft_lower,
         soft_upper_bounds=_soft_upper,
     )
@@ -104,14 +111,14 @@ def _fast_path_scale_converter():
 
 
 def calculate_scaling_factor_and_offset(
-    flat_params,
+    internal_params,
     method="start_values",
     clipping_value=0.1,
     magnitude=1,
 ):
-    x = flat_params.values
-    lower_bounds = flat_params.lower_bounds
-    upper_bounds = flat_params.upper_bounds
+    x = internal_params.values
+    lower_bounds = internal_params.lower_bounds
+    upper_bounds = internal_params.upper_bounds
 
     if method == "start_values":
         raw_factor = np.clip(np.abs(x), clipping_value, np.inf)
diff --git a/src/estimagic/parameters/space_conversion.py b/src/estimagic/parameters/space_conversion.py
index 841617ddc..203aac9e5 100644
--- a/src/estimagic/parameters/space_conversion.py
+++ b/src/estimagic/parameters/space_conversion.py
@@ -38,12 +38,11 @@
 import estimagic.parameters.kernel_transformations as kt
 import numpy as np
 from estimagic.parameters.process_constraints import process_constraints
-from estimagic.parameters.tree_conversion import FlatParams
 
 
 def get_space_converter(
-    flat_params,
-    flat_constraints,
+    internal_params,
+    internal_constraints,
 ):
     """Get functions to convert between in-/external space of params and derivatives.
 
@@ -51,21 +50,29 @@ def get_space_converter(
     for bounds.
 
     Args:
-        flat_params (FlatParams): NamedTuple with flattened parameter values and bounds.
-        flat_constraints (list): List of constraints with processed selector fields.
+        internal_params (InternalParams): NamedTuple with internal parameter values and
+            bounds.
+        internal_constraints (list): List of constraints with processed selector fields.
 
     Returns:
-        SpaceConverter
-        FlatParams: NamedTuple of 1d numpy array with flat and internal params and
-            bounds.
+        SpaceConverter: The space converter.
+        InternalParams: NamedTuple with entries:
+            - value (np.ndarray): Internal parameter values.
+            - lower_bounds (np.ndarray): Lower bounds on the internal params.
+            - upper_bounds (np.ndarray): Upper bounds on the internal params.
+            - soft_lower_bounds (np.ndarray): Soft lower bounds on the internal params.
+            - soft_upper_bounds (np.ndarray): Soft upper bounds on the internal params.
+            - name (list): List of names of the external parameters.
+            - free_mask (np.ndarray): Boolean mask representing which external parameter
+              is free.
 
     """
     transformations, constr_info = process_constraints(
-        constraints=flat_constraints,
-        params_vec=flat_params.values,
-        lower_bounds=flat_params.lower_bounds,
-        upper_bounds=flat_params.upper_bounds,
-        param_names=flat_params.names,
+        constraints=internal_constraints,
+        params_vec=internal_params.values,
+        lower_bounds=internal_params.lower_bounds,
+        upper_bounds=internal_params.upper_bounds,
+        param_names=internal_params.names,
     )
     _params_to_internal = partial(
         reparametrize_to_internal,
@@ -110,27 +117,33 @@ def get_space_converter(
     )
 
     free_mask = constr_info["internal_free"]
-    if flat_params.soft_lower_bounds is not None and not _has_transforming_constraints:
-        _soft_lower = flat_params.soft_lower_bounds[free_mask]
+    if (
+        internal_params.soft_lower_bounds is not None
+        and not _has_transforming_constraints
+    ):
+        _soft_lower = internal_params.soft_lower_bounds[free_mask]
     else:
         _soft_lower = None
 
-    if flat_params.soft_upper_bounds is not None and not _has_transforming_constraints:
-        _soft_upper = flat_params.soft_upper_bounds[free_mask]
+    if (
+        internal_params.soft_upper_bounds is not None
+        and not _has_transforming_constraints
+    ):
+        _soft_upper = internal_params.soft_upper_bounds[free_mask]
     else:
         _soft_upper = None
 
-    internal_params = FlatParams(
-        values=converter.params_to_internal(flat_params.values),
+    params = InternalParams(
+        values=converter.params_to_internal(internal_params.values),
         lower_bounds=constr_info["lower_bounds"],
         upper_bounds=constr_info["upper_bounds"],
-        names=flat_params.names,
+        names=internal_params.names,
         free_mask=free_mask,
         soft_lower_bounds=_soft_lower,
         soft_upper_bounds=_soft_upper,
     )
 
-    return converter, internal_params
+    return converter, params
 
 
 class SpaceConverter(NamedTuple):
@@ -494,3 +507,13 @@ def post_replace_jacobian(post_replacements):
     jacobian[positions_out, :] *= 0
     jacobian[positions_out, positions_in] = 1
     return jacobian
+
+
+class InternalParams(NamedTuple):
+    values: np.ndarray
+    lower_bounds: np.ndarray
+    upper_bounds: np.ndarray
+    soft_lower_bounds: np.ndarray = None
+    soft_upper_bounds: np.ndarray = None
+    names: list = None
+    free_mask: np.ndarray = None
diff --git a/src/estimagic/visualization/slice_plot.py b/src/estimagic/visualization/slice_plot.py
index bea9512c9..23b1bd242 100644
--- a/src/estimagic/visualization/slice_plot.py
+++ b/src/estimagic/visualization/slice_plot.py
@@ -98,7 +98,7 @@ def slice_plot(
 
     func_eval = func(params)
 
-    converter, flat_params = get_converter(
+    converter, internal_params = get_converter(
         params=params,
         constraints=None,
         lower_bounds=lower_bounds,
@@ -109,7 +109,7 @@ def slice_plot(
         scaling_options=None,
     )
 
-    n_params = len(flat_params.values)
+    n_params = len(internal_params.values)
 
     selected = np.arange(n_params, dtype=int)
     if selector is not None:
@@ -119,26 +119,26 @@ def slice_plot(
             tree_just_flatten(selector(helper), registry=registry), dtype=int
         )
 
-    if not np.isfinite(flat_params.lower_bounds[selected]).all():
+    if not np.isfinite(internal_params.lower_bounds[selected]).all():
         raise ValueError("All selected parameters must have finite lower bounds.")
 
-    if not np.isfinite(flat_params.upper_bounds[selected]).all():
+    if not np.isfinite(internal_params.upper_bounds[selected]).all():
         raise ValueError("All selected parameters must have finite upper bounds.")
 
     evaluation_points, metadata = [], []
     for pos in selected:
-        lb = flat_params.lower_bounds[pos]
-        ub = flat_params.upper_bounds[pos]
+        lb = internal_params.lower_bounds[pos]
+        ub = internal_params.upper_bounds[pos]
         grid = np.linspace(lb, ub, n_gridpoints)
-        name = flat_params.names[pos]
+        name = internal_params.names[pos]
         for param_value in grid:
-            if param_value != flat_params.values[pos]:
+            if param_value != internal_params.values[pos]:
                 meta = {
                     "name": name,
                     "Parameter Value": param_value,
                 }
 
-                x = flat_params.values.copy()
+                x = internal_params.values.copy()
                 x[pos] = param_value
                 point = converter.params_from_internal(x)
                 evaluation_points.append(point)
@@ -162,8 +162,8 @@ def slice_plot(
     func_values += [converter.func_to_internal(func_eval)] * len(selected)
     for pos in selected:
         meta = {
-            "name": flat_params.names[pos],
-            "Parameter Value": flat_params.values[pos],
+            "name": internal_params.names[pos],
+            "Parameter Value": internal_params.values[pos],
         }
         metadata.append(meta)
 
@@ -188,7 +188,7 @@ def slice_plot(
 
     plots_dict = {}
     for pos in selected:
-        par_name = flat_params.names[pos]
+        par_name = internal_params.names[pos]
         if param_names is not None and par_name in param_names:
             par_name = param_names[par_name]
 
@@ -201,7 +201,7 @@ def slice_plot(
         )
         subfig.add_trace(
             go.Scatter(
-                x=[flat_params.values[pos]],
+                x=[internal_params.values[pos]],
                 y=[converter.func_to_internal(func_eval)],
                 marker={"color": color},
             )


Test patch:
diff --git a/tests/estimation/test_estimate_ml.py b/tests/estimation/test_estimate_ml.py
index f00575c34..1ffa852b3 100644
--- a/tests/estimation/test_estimate_ml.py
+++ b/tests/estimation/test_estimate_ml.py
@@ -10,6 +10,7 @@
 from estimagic.examples.logit import logit_hessian
 from estimagic.examples.logit import logit_loglike
 from estimagic.examples.logit import logit_loglike_and_derivative as llad
+from numpy.testing import assert_array_equal
 from scipy.stats import multivariate_normal
 from statsmodels.base.model import GenericLikelihoodModel
 
@@ -90,6 +91,7 @@ def logit_np_inputs():
     spector_data = sm.datasets.spector.load_pandas()
     spector_data.exog = sm.add_constant(spector_data.exog)
     x_df = sm.add_constant(spector_data.exog)
+
     out = {
         "y": spector_data.endog,
         "x": x_df.to_numpy(),
@@ -197,7 +199,8 @@ def test_estimate_ml_with_logit_no_constraints(
 
         # compare estimated standard errors
         exp_se = getattr(exp, f"bse{statsmodels_suffix_map[method]}")
-        aaae(got.se(method=method), exp_se, decimal=3)
+        got_se = got.se(method=method)
+        aaae(got_se, exp_se, decimal=3)
 
         # compare estimated confidence interval
         if method == "hessian":
@@ -231,10 +234,104 @@ def test_estimate_ml_with_logit_no_constraints(
         aaae(got._p_values, got.p_values())
 
 
-def test_estimate_ml_optimize_options_false(fitted_logit_model, logit_inputs):
+test_cases_constr = list(
+    itertools.product(
+        [None, logit_jacobian],  # jacobian
+        [
+            {"loc": [1, 2, 3], "type": "covariance"},
+            {"loc": [0, 1], "type": "linear", "lower_bound": -20, "weights": 1},
+            {"loc": [0, 1], "type": "increasing"},
+        ],
+    )
+)
+
+
+@pytest.mark.parametrize("jacobian, constraints", test_cases_constr)
+def test_estimate_ml_with_logit_constraints(
+    fitted_logit_model,
+    logit_np_inputs,
+    jacobian,
+    constraints,
+):
+    """
+    Test that estimate_ml computes correct params and standard errors under different
+    scenarios with constraints.
+    """
+    seed = 1234
+
+    # ==================================================================================
+    # estimate
+    # ==================================================================================
+
+    kwargs = {"y": logit_np_inputs["y"], "x": logit_np_inputs["x"]}
+
+    optimize_options = {
+        "algorithm": "scipy_lbfgsb",
+        "algo_options": {"convergence.relative_criterion_tolerance": 1e-12},
+    }
+
+    if "criterion_and_derivative" in optimize_options:
+        optimize_options["criterion_and_derivative_kwargs"] = kwargs
+
+    got = estimate_ml(
+        loglike=logit_loglike,
+        params=logit_np_inputs["params"],
+        loglike_kwargs=kwargs,
+        optimize_options=optimize_options,
+        jacobian=jacobian,
+        jacobian_kwargs=kwargs,
+        constraints=constraints,
+    )
+
+    # ==================================================================================
+    # test
+    # ==================================================================================
+
+    exp = fitted_logit_model
+
+    methods = ["jacobian", "hessian", "robust"]
+
+    statsmodels_suffix_map = {
+        "jacobian": "jac",
+        "hessian": "",
+        "robust": "jhj",
+    }
+
+    # compare estimated parameters
+    aaae(got.params, exp.params, decimal=3)
+
+    for method in methods:
+
+        # compare estimated standard errors
+        exp_se = getattr(exp, f"bse{statsmodels_suffix_map[method]}")
+        got_se = got.se(method=method, seed=seed)
+        corr = np.corrcoef(got_se, exp_se)
+        aaae(corr, np.ones_like(corr), decimal=4)
+
+        # compare estimated confidence interval
+        if method == "hessian":
+            lower, upper = got.ci(method=method, seed=seed)
+            exp_lower = exp.conf_int().T[0]
+            exp_upper = exp.conf_int().T[1]
+            corr_lower = np.corrcoef(lower, exp_lower)
+            corr_upper = np.corrcoef(upper, exp_upper)
+            aaae(corr_lower, np.ones_like(corr), decimal=4)
+            aaae(corr_upper, np.ones_like(corr), decimal=4)
+
+        summary = got.summary(method=method, seed=seed)
+
+        aaae(summary["value"], exp.params, decimal=3)
+        aaae(summary["standard_error"], got.se(method=method, seed=seed))
+        lower, upper = got.ci(method=method, seed=seed)
+        aaae(summary["ci_lower"], lower)
+        aaae(summary["ci_upper"], upper)
+        aaae(summary["p_value"], got.p_values(method=method, seed=seed))
+
+
+def test_estimate_ml_optimize_options_false(fitted_logit_model, logit_np_inputs):
     """Test that estimate_ml computes correct covariances given correct params."""
 
-    kwargs = {"y": logit_inputs["y"], "x": logit_inputs["x"]}
+    kwargs = {"y": logit_np_inputs["y"], "x": logit_np_inputs["x"]}
 
     params = pd.DataFrame({"value": fitted_logit_model.params})
 
@@ -328,3 +425,23 @@ def test_to_pickle(normal_inputs, tmp_path):
     )
 
     got.to_pickle(tmp_path / "bla.pkl")
+
+
+def test_caching(normal_inputs):
+    kwargs = {"y": normal_inputs["y"]}
+
+    start_params = {"mean": 5, "sd": 3}
+
+    got = estimate_ml(
+        loglike=normal_loglike,
+        params=start_params,
+        loglike_kwargs=kwargs,
+        optimize_options="scipy_lbfgsb",
+        lower_bounds={"sd": 0.0001},
+        jacobian_kwargs=kwargs,
+        constraints=[{"selector": lambda p: p["sd"], "type": "sdcorr"}],
+    )
+
+    assert got._cache == {}
+    cov = got.cov(method="robust", return_type="array")
+    assert_array_equal(list(got._cache.values())[0], cov)
diff --git a/tests/estimation/test_estimate_msm.py b/tests/estimation/test_estimate_msm.py
index 9a052dfb1..8a5b5bc39 100644
--- a/tests/estimation/test_estimate_msm.py
+++ b/tests/estimation/test_estimate_msm.py
@@ -8,6 +8,7 @@
 from estimagic.shared.check_option_dicts import check_numdiff_options
 from estimagic.shared.check_option_dicts import check_optimization_options
 from numpy.testing import assert_array_almost_equal as aaae
+from numpy.testing import assert_array_equal
 
 
 def _sim_pd(params):
@@ -176,3 +177,24 @@ def test_to_pickle(tmp_path):
     )
 
     calculated.to_pickle(tmp_path / "bla.pkl")
+
+
+def test_caching():
+    start_params = np.array([3, 2, 1])
+
+    # abuse simulate_moments to get empirical moments in correct format
+    empirical_moments = _sim_np(np.zeros(3))
+    if isinstance(empirical_moments, dict):
+        empirical_moments = empirical_moments["simulated_moments"]
+
+    got = estimate_msm(
+        simulate_moments=_sim_np,
+        empirical_moments=empirical_moments,
+        moments_cov=cov_np,
+        params=start_params,
+        optimize_options="scipy_lbfgsb",
+    )
+
+    assert got._cache == {}
+    cov = got.cov(method="robust", return_type="array")
+    assert_array_equal(list(got._cache.values())[0], cov)
diff --git a/tests/inference/test_bootstrap.py b/tests/inference/test_bootstrap.py
index be3d9de2f..2c78af01f 100644
--- a/tests/inference/test_bootstrap.py
+++ b/tests/inference/test_bootstrap.py
@@ -2,9 +2,8 @@
 import pandas as pd
 import pytest
 from estimagic.inference.bootstrap import bootstrap
-from estimagic.inference.bootstrap import bootstrap_from_outcomes
-from pandas.testing import assert_frame_equal as afe
-from pandas.testing import assert_series_equal as ase
+from pandas.testing import assert_frame_equal
+from pandas.testing import assert_series_equal
 
 
 @pytest.fixture
@@ -76,113 +75,119 @@ def test_bootstrap_with_outcome_kwargs(shift, setup):
     )
 
     expected = pd.Series([2.5, 7.0], index=["x1", "x2"])
-    ase(result.base_outcome, expected + shift)
+    assert_series_equal(result.base_outcome, expected + shift)
+
+
+def test_bootstrap_existing_outcomes(setup):
+    result = bootstrap(
+        data=setup["df"],
+        outcome=_outcome_func,
+        n_draws=2,
+    )
+    assert len(result.outcomes) == 2
+    result = bootstrap(
+        data=setup["df"],
+        outcome=_outcome_func,
+        existing_outcomes=result,
+        n_draws=1,
+    )
+    assert len(result.outcomes) == 1
 
 
 def test_bootstrap_from_outcomes(setup, expected):
 
-    result = bootstrap_from_outcomes(
-        base_outcome=_outcome_func(setup["df"]),
-        bootstrap_outcomes=setup["estimates_pytree"],
+    result = bootstrap(
+        outcome=_outcome_func(setup["df"]),
+        existing_outcomes=setup["estimates_pytree"],
     )
 
-    outcomes = result.outcomes()
+    outcomes = result.outcomes
     lower_ci, upper_ci = result.ci()
     covariance = result.cov()
     standard_errors = result.se()
 
-    with pytest.raises(NotImplementedError) as error:
-        assert result._p_values()
-    assert str(error.value) == "Bootstrapped p-values are not implemented yet."
-
-    with pytest.raises(NotImplementedError) as error:
-        assert result._summary()
-    assert str(error.value) == "summary is not implemented yet."
+    with pytest.raises(NotImplementedError):
+        assert result._p_values
 
     # Use rounding to adjust precision and ensure reproducibility accross all
     # supported pandas versions.
     for i in range(len(outcomes)):
-        ase(outcomes[i], setup["estimates_pytree"][i])
+        assert_series_equal(outcomes[i], setup["estimates_pytree"][i])
 
-    ase(lower_ci.round(2), expected["lower_ci"].round(2))
-    ase(upper_ci.round(2), expected["upper_ci"].round(2))
-    afe(covariance.round(2), expected["cov"].round(2))
-    ase(standard_errors.round(2), expected["se"].round(2))
+    assert_series_equal(lower_ci.round(2), expected["lower_ci"].round(2))
+    assert_series_equal(upper_ci.round(2), expected["upper_ci"].round(2))
+    assert_frame_equal(covariance.round(2), expected["cov"].round(2))
+    assert_series_equal(standard_errors.round(2), expected["se"].round(2))
 
 
 def test_bootstrap_from_outcomes_private_methods(setup, expected):
 
-    result = bootstrap_from_outcomes(
-        base_outcome=_outcome_func(setup["df"]),
-        bootstrap_outcomes=setup["estimates_pytree"],
+    result = bootstrap(
+        outcome=_outcome_func(setup["df"]),
+        existing_outcomes=setup["estimates_pytree"],
     )
 
-    outcomes = result._outcomes
+    outcomes = result.outcomes
     lower_ci, upper_ci = result._ci
     covariance = result._cov
     standard_errors = result._se
 
-    with pytest.raises(NotImplementedError) as error:
-        assert result._p_values()
-    assert str(error.value) == "Bootstrapped p-values are not implemented yet."
-
-    with pytest.raises(NotImplementedError) as error:
-        assert result._summary()
-    assert str(error.value) == "summary is not implemented yet."
+    with pytest.raises(NotImplementedError):
+        assert result._p_values
 
     for i in range(len(outcomes)):
-        ase(outcomes[i], setup["estimates_pytree"][i])
+        assert_series_equal(outcomes[i], setup["estimates_pytree"][i])
 
-    ase(lower_ci.round(2), expected["lower_ci"].round(2))
-    ase(upper_ci.round(2), expected["upper_ci"].round(2))
-    afe(covariance.round(2), expected["cov"].round(2))
-    ase(standard_errors.round(2), expected["se"].round(2))
+    assert_series_equal(lower_ci.round(2), expected["lower_ci"].round(2))
+    assert_series_equal(upper_ci.round(2), expected["upper_ci"].round(2))
+    assert_frame_equal(covariance.round(2), expected["cov"].round(2))
+    assert_series_equal(standard_errors.round(2), expected["se"].round(2))
 
 
 def test_bootstrap_from_outcomes_single_outcome(setup, expected):
 
-    result = bootstrap_from_outcomes(
-        base_outcome=_outcome_func(setup["df"]["x1"]),
-        bootstrap_outcomes=setup["estimates_pytree_x1"],
+    result = bootstrap(
+        outcome=_outcome_func(setup["df"]["x1"]),
+        existing_outcomes=setup["estimates_pytree_x1"],
     )
 
-    outcomes = result.outcomes()
+    outcomes = result.outcomes
     lower_ci, upper_ci = result.ci()
 
     for i in range(len(outcomes)):
-        ase(outcomes[i], setup["estimates_pytree_x1"][i])
+        assert_series_equal(outcomes[i], setup["estimates_pytree_x1"][i])
 
-    ase(lower_ci.round(2), expected["lower_ci_x1"].round(2))
-    ase(upper_ci.round(2), expected["upper_ci_x1"].round(2))
+    assert_series_equal(lower_ci.round(2), expected["lower_ci_x1"].round(2))
+    assert_series_equal(upper_ci.round(2), expected["upper_ci_x1"].round(2))
 
 
 @pytest.mark.parametrize("input_type", ["arr", "df", "dict"])
 def test_wrong_input_types(input_type, setup):
-    with pytest.raises(TypeError):
-        assert bootstrap_from_outcomes(
-            base_outcome=_outcome_func(setup["df"]),
-            bootstrap_outcomes=setup["estimates_" + input_type],
+    with pytest.raises(ValueError):
+        assert bootstrap(
+            outcome=_outcome_func(setup["df"]),
+            existing_outcomes=setup["estimates_" + input_type],
         )
 
 
 @pytest.mark.parametrize("return_type", ["array", "dataframe", "pytree"])
 def test_cov_correct_return_type(return_type, setup):
-    result = bootstrap_from_outcomes(
-        base_outcome=_outcome_func(setup["df"]),
-        bootstrap_outcomes=setup["estimates_pytree"],
+    result = bootstrap(
+        outcome=_outcome_func(setup["df"]),
+        existing_outcomes=setup["estimates_pytree"],
     )
     _ = result.cov(return_type=return_type)
 
 
 def test_cov_wrong_return_type(setup):
-    result = bootstrap_from_outcomes(
-        base_outcome=_outcome_func(setup["df"]),
-        bootstrap_outcomes=setup["estimates_pytree"],
+    result = bootstrap(
+        outcome=_outcome_func(setup["df"]),
+        existing_outcomes=setup["estimates_pytree"],
     )
 
     expected_msg = "return_type must be one of pytree, array, or dataframe, not dict."
 
-    with pytest.raises(TypeError) as error:
+    with pytest.raises(ValueError) as error:
         assert result.cov(return_type="dict")
 
     assert str(error.value) == expected_msg
diff --git a/tests/inference/test_bootstrap_ci.py b/tests/inference/test_bootstrap_ci.py
index acbb1efc8..4380eaec4 100644
--- a/tests/inference/test_bootstrap_ci.py
+++ b/tests/inference/test_bootstrap_ci.py
@@ -3,7 +3,7 @@
 import numpy as np
 import pandas as pd
 import pytest
-from estimagic.inference.bootstrap_ci import compute_ci
+from estimagic.inference.bootstrap_ci import calculate_ci
 from estimagic.inference.bootstrap_helpers import check_inputs
 from estimagic.parameters.tree_registry import get_registry
 from numpy.testing import assert_array_almost_equal as aaae
@@ -72,7 +72,7 @@ def outcome_flat(data):
         return tree_just_flatten(outcome(data), registry=registry)
 
     base_outcome = outcome_flat(setup["df"])
-    lower, upper = compute_ci(base_outcome, setup["estimates"], ci_method=method)
+    lower, upper = calculate_ci(base_outcome, setup["estimates"], ci_method=method)
 
     aaae(lower, expected[method + "_ci"][:, 0])
     aaae(upper, expected[method + "_ci"][:, 1])
@@ -80,16 +80,16 @@ def outcome_flat(data):
 
 def test_check_inputs_data():
     data = "this is not a data frame"
-    with pytest.raises(ValueError) as excinfo:
+    with pytest.raises(TypeError) as excinfo:
         check_inputs(data=data)
-    assert "Input 'data' must be DataFrame." == str(excinfo.value)
+    assert "Data must be a pandas.DataFrame." == str(excinfo.value)
 
 
 def test_check_inputs_cluster_by(setup):
     cluster_by = "this is not a column name of df"
     with pytest.raises(ValueError) as excinfo:
         check_inputs(data=setup["df"], cluster_by=cluster_by)
-    assert "Input 'cluster_by' must be None or a column name of DataFrame." == str(
+    assert "Input 'cluster_by' must be None or a column name of 'data'." == str(
         excinfo.value
     )
 
diff --git a/tests/inference/test_shared.py b/tests/inference/test_shared.py
index 5d9988096..2de522983 100644
--- a/tests/inference/test_shared.py
+++ b/tests/inference/test_shared.py
@@ -1,15 +1,19 @@
+from typing import NamedTuple
+
 import numpy as np
 import pandas as pd
 import pytest
-from estimagic.inference.shared import calculate_inference_quantities
+from estimagic.inference.shared import _to_numpy
+from estimagic.inference.shared import calculate_estimation_summary
+from estimagic.inference.shared import get_derivative_case
 from estimagic.inference.shared import process_pandas_arguments
-from estimagic.parameters.tree_conversion import FlatParams
+from estimagic.inference.shared import transform_covariance
+from estimagic.inference.shared import transform_free_cov_to_cov
+from estimagic.inference.shared import transform_free_values_to_params_tree
 from estimagic.parameters.tree_registry import get_registry
-from pandas.testing import assert_frame_equal
+from numpy.testing import assert_array_almost_equal as aaae
 from pybaum import leaf_names
 from pybaum import tree_equal
-from pybaum import tree_just_flatten
-from pybaum import tree_map
 
 
 @pytest.fixture
@@ -40,109 +44,230 @@ def test_process_pandas_arguments_incompatible_names(inputs):
         process_pandas_arguments(**inputs)
 
 
-def test_calculate_inference_quantities():
-    """Test calculate inference quantities for many relevant cases."""
-
-    # ==================================================================================
-    # create test inputs
-    # ==================================================================================
-
-    estimates = [
-        0,
-        {
-            "a": 1,
-            "b": (2, 3),
-            "c": np.array([4, 5]),
-            "d": pd.Series([6, 7], index=["a", "b"]),
-            "e": pd.DataFrame(
-                {"col1": [8, 9, 10], "col2": [11, 12, 13]},
-                index=pd.MultiIndex.from_tuples([("a", "b"), ("c", "d"), ("e", "f")]),
-            ),
-            "f": np.array([[15, 16], [17, 18]]),
-            "g": pd.DataFrame({"value": [19, 20], "col2": [-1, -1]}, index=["g", "h"]),
-        },
-    ]
+def _from_internal(x, return_type="flat"):
+    return x
 
-    registry = get_registry(extended=True)
-    estimates_flat = tree_just_flatten(estimates, registry=registry)
-    names = leaf_names(estimates, registry=registry)
 
-    flat_estimates = FlatParams(
-        values=estimates_flat, lower_bounds=None, upper_bounds=None, names=names
+class FakeConverter(NamedTuple):
+    has_transforming_constraints: bool = True
+    params_from_internal: callable = _from_internal
+
+
+class FakeInternalParams(NamedTuple):
+    values: np.ndarray = np.arange(2)
+    lower_bounds: np.ndarray = np.full(2, -np.inf)
+    upper_bounds: np.ndarray = np.full(2, np.inf)
+    free_mask: np.ndarray = np.array([True, True])
+
+
+def test_transform_covariance_no_bounds():
+    internal_cov = np.eye(2)
+
+    converter = FakeConverter()
+    internal_params = FakeInternalParams()
+
+    np.random.seed(1234)
+
+    got = transform_covariance(
+        internal_params=internal_params,
+        internal_cov=internal_cov,
+        converter=converter,
+        n_samples=100,
+        bounds_handling="ignore",
+    )
+
+    np.random.seed(1234)
+    expected_sample = np.random.multivariate_normal(np.arange(2), np.eye(2), 100)
+    expected = np.cov(expected_sample, rowvar=False)
+
+    aaae(got, expected)
+
+
+def test_transform_covariance_with_clipping():
+    internal_cov = np.eye(2)
+
+    converter = FakeConverter()
+    internal_params = FakeInternalParams(
+        lower_bounds=np.ones(2), upper_bounds=np.ones(2)
     )
 
-    free_cov = pd.DataFrame(
-        np.diag(np.arange(len(estimates_flat))), index=names, columns=names
+    np.random.seed(1234)
+
+    got = transform_covariance(
+        internal_params=internal_params,
+        internal_cov=internal_cov,
+        converter=converter,
+        n_samples=100,
+        bounds_handling="clip",
     )
 
-    ci_level = 0  # implies that ci_lower = value = ci_upper
+    expected = np.zeros((2, 2))
 
-    # ==================================================================================
-    # create expected output
-    # ==================================================================================
+    aaae(got, expected)
 
-    df = pd.DataFrame(
-        {
-            "value": flat_estimates.values,
-            "standard_error": np.sqrt(np.diag(free_cov)),
-            "ci_lower": flat_estimates.values,
-            "ci_upper": flat_estimates.values,
-        },
-        index=names,
+
+def test_transform_covariance_invalid_bounds():
+    internal_cov = np.eye(2)
+
+    converter = FakeConverter()
+    internal_params = FakeInternalParams(
+        lower_bounds=np.ones(2), upper_bounds=np.ones(2)
     )
-    df[["ci_lower", "ci_upper"]] = df[["ci_lower", "ci_upper"]].astype(float)
-
-    df_c = df.loc[["1_c_0", "1_c_1"]]
-    df_c.index = pd.RangeIndex(stop=2)
-
-    df_e = df.loc[
-        [
-            "1_e_a_b_col1",
-            "1_e_a_b_col2",
-            "1_e_c_d_col1",
-            "1_e_c_d_col2",
-            "1_e_e_f_col1",
-            "1_e_e_f_col2",
-        ]
-    ]
-    df_e.index = pd.MultiIndex.from_tuples(
-        [
-            ("a", "b", "col1"),
-            ("a", "b", "col2"),
-            ("c", "d", "col1"),
-            ("c", "d", "col2"),
-            ("e", "f", "col1"),
-            ("e", "f", "col2"),
-        ]
+
+    np.random.seed(1234)
+
+    with pytest.raises(ValueError):
+        transform_covariance(
+            internal_params=internal_params,
+            internal_cov=internal_cov,
+            converter=converter,
+            n_samples=10,
+            bounds_handling="raise",
+        )
+
+
+class FakeFreeParams(NamedTuple):
+    free_mask: np.ndarray = np.array([True, False, True])
+    all_names: list = ["a", "b", "c"]
+    free_names: list = ["a", "c"]
+
+
+def test_transform_free_cov_to_cov_pytree():
+    got = transform_free_cov_to_cov(
+        free_cov=np.eye(2),
+        free_params=FakeFreeParams(),
+        params={"a": 1, "b": 2, "c": 3},
+        return_type="pytree",
     )
 
-    df_f = df.loc[["1_f_0_0", "1_f_0_1", "1_f_1_0", "1_f_1_1"]]
-    df_f.index = pd.MultiIndex.from_tuples([(0, 0), (0, 1), (1, 0), (1, 1)])
-
-    expected = [
-        df.loc[["0"]].set_axis([0]),
-        {
-            "a": df.loc[["1_a"]].set_axis([0]),
-            "b": (df.loc[["1_b_0"]].set_axis([0]), df.loc[["1_b_1"]].set_axis([0])),
-            "c": df_c,
-            "d": df.loc[["1_d_a", "1_d_b"]].set_axis(["a", "b"]),
-            "e": df_e,
-            "f": df_f,
-            "g": df.loc[["1_g_g", "1_g_h"]].set_axis(["g", "h"]),
-        },
-    ]
+    assert got["a"]["a"] == 1
+    assert got["c"]["c"] == 1
+    assert got["a"]["c"] == 0
+    assert got["c"]["a"] == 0
+    assert np.isnan(got["a"]["b"])
 
-    # ==================================================================================
-    # compute and compare
-    # ==================================================================================
 
-    got = calculate_inference_quantities(estimates, flat_estimates, free_cov, ci_level)
+def test_transform_free_cov_to_cov_array():
+    got = transform_free_cov_to_cov(
+        free_cov=np.eye(2),
+        free_params=FakeFreeParams(),
+        params={"a": 1, "b": 2, "c": 3},
+        return_type="array",
+    )
+
+    expected = np.array([[1, np.nan, 0], [np.nan, np.nan, np.nan], [0, np.nan, 1]])
+
+    assert np.array_equal(got, expected, equal_nan=True)
+
+
+def test_transform_free_cov_to_cov_dataframe():
+    got = transform_free_cov_to_cov(
+        free_cov=np.eye(2),
+        free_params=FakeFreeParams(),
+        params={"a": 1, "b": 2, "c": 3},
+        return_type="dataframe",
+    )
+
+    expected = np.array([[1, np.nan, 0], [np.nan, np.nan, np.nan], [0, np.nan, 1]])
+
+    assert np.array_equal(got.to_numpy(), expected, equal_nan=True)
+    assert isinstance(got, pd.DataFrame)
+    assert list(got.columns) == list("abc")
+    assert list(got.index) == list("abc")
+
+
+def test_transform_free_cov_to_cov_invalid():
+    with pytest.raises(ValueError):
+        transform_free_cov_to_cov(
+            free_cov=np.eye(2),
+            free_params=FakeFreeParams(),
+            params={"a": 1, "b": 2, "c": 3},
+            return_type="bla",
+        )
+
+
+def test_transform_free_values_to_params_tree():
+    got = transform_free_values_to_params_tree(
+        values=np.array([10, 11]),
+        free_params=FakeFreeParams(),
+        params={"a": 1, "b": 2, "c": 3},
+    )
+
+    assert got["a"] == 10
+    assert got["c"] == 11
+    assert np.isnan(got["b"])
+
+
+def test_get_derivative_case():
+    assert get_derivative_case(lambda x: True) == "closed-form"
+    assert get_derivative_case(False) == "skip"
+    assert get_derivative_case(None) == "numerical"
+
+
+def test_to_numpy_invalid():
+    with pytest.raises(TypeError):
+        _to_numpy(15)
+
+
+def test_calculate_estimation_summary():
+    # input data
+    summary_data = {
+        "value": {
+            "a": pd.Series([0], index=["i"]),
+            "b": pd.DataFrame({"c1": [1], "c2": [2]}),
+        },
+        "standard_error": {
+            "a": pd.Series([0.1], index=["i"]),
+            "b": pd.DataFrame({"c1": [0.2], "c2": [0.3]}),
+        },
+        "ci_lower": {
+            "a": pd.Series([-0.2], index=["i"]),
+            "b": pd.DataFrame({"c1": [-0.4], "c2": [-0.6]}),
+        },
+        "ci_upper": {
+            "a": pd.Series([0.2], index=["i"]),
+            "b": pd.DataFrame({"c1": [0.4], "c2": [0.6]}),
+        },
+        "p_value": {
+            "a": pd.Series([0.001], index=["i"]),
+            "b": pd.DataFrame({"c1": [0.2], "c2": [0.07]}),
+        },
+        "free": np.array([True, True, True]),
+    }
+
+    registry = get_registry(extended=True)
+    names = leaf_names(summary_data["value"], registry=registry)
+    free_names = names
 
-    # drop irrelevant columns
-    got = tree_map(lambda df: df.drop(columns=["stars", "p_value", "free"]), got)
+    # function call
+    summary = calculate_estimation_summary(summary_data, names, free_names)
 
-    # for debugging purposes we first compare each leaf
-    for got_leaf, exp_leaf in zip(tree_just_flatten(got), tree_just_flatten(expected)):
-        assert_frame_equal(got_leaf, exp_leaf)
+    # expectations
+    expectation = {
+        "a": pd.DataFrame(
+            {
+                "value": 0,
+                "standard_error": 0.1,
+                "ci_lower": -0.2,
+                "ci_upper": 0.2,
+                "p_value": 0.001,
+                "free": True,
+                "stars": "***",
+            },
+            index=["i"],
+        ),
+        "b": pd.DataFrame(
+            {
+                "value": [1, 2],
+                "standard_error": [0.2, 0.3],
+                "ci_lower": [-0.4, -0.6],
+                "ci_upper": [0.4, 0.6],
+                "p_value": [0.2, 0.7],
+                "free": [True, True],
+                "stars": ["", "*"],
+            },
+            index=pd.MultiIndex.from_tuples([(0, "c1"), (0, "c2")]),
+        ),
+    }
 
-    assert tree_equal(expected, got)
+    tree_equal(summary, expectation)
diff --git a/tests/optimization/test_with_constraints.py b/tests/optimization/test_with_constraints.py
index e4dd86861..1fca2cbe7 100644
--- a/tests/optimization/test_with_constraints.py
+++ b/tests/optimization/test_with_constraints.py
@@ -11,6 +11,7 @@
 import numpy as np
 import pandas as pd
 import pytest
+import statsmodels.api as sm
 from estimagic.examples.criterion_functions import rosenbrock_dict_criterion
 from estimagic.examples.criterion_functions import rosenbrock_gradient
 from estimagic.examples.criterion_functions import (
@@ -23,13 +24,13 @@
 from estimagic.examples.criterion_functions import sos_ls_jacobian
 from estimagic.examples.criterion_functions import trid_gradient
 from estimagic.examples.criterion_functions import trid_scalar_criterion
+from estimagic.examples.logit import logit_loglike
 from estimagic.exceptions import InvalidConstraintError
 from estimagic.exceptions import InvalidParamsError
 from estimagic.optimization.optimize import maximize
 from estimagic.optimization.optimize import minimize
 from numpy.testing import assert_array_almost_equal as aaae
 
-
 FUNC_INFO = {
     "sos": {
         "criterion": sos_dict_criterion,
@@ -149,6 +150,7 @@ def test_constrained_minimization(
         algorithm=algorithm,
         derivative=derivative,
         constraints=constraints,
+        algo_options={"convergence.relative_criterion_tolerance": 1e-12},
     )
 
     if params_type == "pandas":
@@ -189,11 +191,11 @@ def test_three_independent_constraints():
         params=params,
         algorithm="scipy_lbfgsb",
         constraints=constraints,
+        algo_options={"convergence.relative_criterion_tolerance": 1e-12},
     )
-
     expected = np.array([0] * 4 + [4, 5] + [0] + [7.5] * 2 + [0])
 
-    aaae(res.params, expected, decimal=5)
+    aaae(res.params, expected, decimal=4)
 
 
 INVALID_CONSTRAINT_COMBIS = [
@@ -219,7 +221,7 @@ def test_incompatible_constraints_raise_errors(constraints):
 
 
 def test_bug_from_copenhagen_presentation():
-    # make sure maximum of work hours is optimal
+    # Make sure maximum of work hours is optimal
     def u(params):
         return params["work"]["hours"] ** 2
 
@@ -282,3 +284,23 @@ def criterion(x):
             algorithm="scipy_lbfgsb",
             constraints=[{"loc": [1, 2], "type": "probability"}],
         )
+
+
+def test_covariance_constraint_in_2_by_2_case():
+    spector_data = sm.datasets.spector.load_pandas()
+    spector_data.exog = sm.add_constant(spector_data.exog)
+    x_df = sm.add_constant(spector_data.exog)
+
+    start_params = np.array([-10, 2, 0.2, 2])
+    kwargs = {"y": spector_data.endog, "x": x_df.to_numpy()}
+
+    result = maximize(
+        criterion=logit_loglike,
+        criterion_kwargs=kwargs,
+        params=start_params,
+        algorithm="scipy_lbfgsb",
+        constraints={"loc": [1, 2, 3], "type": "covariance"},
+    )
+
+    expected = np.array([-13.0213351, 2.82611417, 0.09515704, 2.37867869])
+    aaae(result.params, expected, decimal=4)
diff --git a/tests/parameters/test_scale_conversion.py b/tests/parameters/test_scale_conversion.py
index a513c363d..e2b281de2 100644
--- a/tests/parameters/test_scale_conversion.py
+++ b/tests/parameters/test_scale_conversion.py
@@ -1,20 +1,20 @@
 import numpy as np
 import pytest
 from estimagic import first_derivative
+from estimagic.parameters.conversion import InternalParams
 from estimagic.parameters.scale_conversion import get_scale_converter
-from estimagic.parameters.tree_conversion import FlatParams
 from numpy.testing import assert_array_almost_equal as aaae
 from numpy.testing import assert_array_equal as aae
 
 
 TEST_CASES = {
-    "start_values": FlatParams(
+    "start_values": InternalParams(
         values=np.array([0, 1, 1, 1, 1, 1]),
         lower_bounds=np.array([-2, 0, 0.5, 2 / 3, 3 / 4, 4 / 5]),
         upper_bounds=np.array([2, 2, 1.5, 4 / 3, 5 / 4, 6 / 5]),
         names=None,
     ),
-    "bounds": FlatParams(
+    "bounds": InternalParams(
         values=np.full(6, 0.5),
         lower_bounds=np.zeros(6),
         upper_bounds=np.ones(6),
@@ -28,7 +28,7 @@
 
 @pytest.mark.parametrize("method, expected", PARAMETRIZATION, ids=IDS)
 def test_get_scale_converter_active(method, expected):
-    params = FlatParams(
+    params = InternalParams(
         values=np.arange(6),
         lower_bounds=np.arange(6) - 1,
         upper_bounds=np.arange(6) + 1,
@@ -41,7 +41,7 @@ def test_get_scale_converter_active(method, expected):
     }
 
     converter, scaled = get_scale_converter(
-        flat_params=params,
+        internal_params=params,
         scaling=True,
         scaling_options=scaling_options,
     )
@@ -63,7 +63,7 @@ def test_get_scale_converter_active(method, expected):
 
 
 def test_scale_conversion_fast_path():
-    params = FlatParams(
+    params = InternalParams(
         values=np.arange(6),
         lower_bounds=np.arange(6) - 1,
         upper_bounds=np.arange(6) + 1,
@@ -71,7 +71,7 @@ def test_scale_conversion_fast_path():
     )
 
     converter, scaled = get_scale_converter(
-        flat_params=params,
+        internal_params=params,
         scaling=False,
         scaling_options=None,
     )
diff --git a/tests/parameters/test_space_conversion.py b/tests/parameters/test_space_conversion.py
index 15d755664..b7afa7cc0 100644
--- a/tests/parameters/test_space_conversion.py
+++ b/tests/parameters/test_space_conversion.py
@@ -4,13 +4,13 @@
 from estimagic.parameters.space_conversion import _multiply_from_left
 from estimagic.parameters.space_conversion import _multiply_from_right
 from estimagic.parameters.space_conversion import get_space_converter
-from estimagic.parameters.tree_conversion import FlatParams
+from estimagic.parameters.space_conversion import InternalParams
 from numpy.testing import assert_array_almost_equal as aaae
 
 
 def _get_test_case_no_constraint():
     n_params = 10
-    fp = FlatParams(
+    fp = InternalParams(
         values=np.arange(n_params),
         lower_bounds=np.full(n_params, -1),
         upper_bounds=np.full(n_params, 11),
@@ -22,7 +22,7 @@ def _get_test_case_no_constraint():
 
 
 def _get_test_case_fixed(with_value):
-    fp = FlatParams(
+    fp = InternalParams(
         values=np.arange(5),
         lower_bounds=np.full(5, -np.inf),
         upper_bounds=np.full(5, np.inf),
@@ -33,7 +33,7 @@ def _get_test_case_fixed(with_value):
     else:
         constraints = [{"index": [0, 2, 4], "type": "fixed"}]
 
-    internal = FlatParams(
+    internal = InternalParams(
         values=np.array([1, 3]),
         lower_bounds=np.full(2, -np.inf),
         upper_bounds=np.full(2, np.inf),
@@ -44,14 +44,14 @@ def _get_test_case_fixed(with_value):
 
 
 def _get_test_case_increasing(as_one):
-    fp = FlatParams(
+    fp = InternalParams(
         values=np.array([0.1, 2.2, 2.3, 10.1, -1]),
         lower_bounds=np.full(5, -np.inf),
         upper_bounds=np.full(5, np.inf),
         names=list("abcde"),
     )
 
-    internal = FlatParams(
+    internal = InternalParams(
         values=np.array([0.1, -2.1, -0.1, -7.8, -1]),
         lower_bounds=np.full(5, -np.inf),
         upper_bounds=np.array([np.inf, 0, 0, 0, np.inf]),
@@ -70,14 +70,14 @@ def _get_test_case_increasing(as_one):
 
 
 def _get_test_case_decreasing(as_one):
-    fp = FlatParams(
+    fp = InternalParams(
         values=np.array([0.1, 2.2, 2.3, 10.1, -1]),
         lower_bounds=np.full(5, -np.inf),
         upper_bounds=np.full(5, np.inf),
         names=list("abcde"),
     )
 
-    internal = FlatParams(
+    internal = InternalParams(
         values=np.array([0.1, -2.1, -0.1, -7.8, -1]),
         lower_bounds=np.full(5, -np.inf),
         upper_bounds=np.array([np.inf, 0, 0, 0, np.inf]),
@@ -96,14 +96,14 @@ def _get_test_case_decreasing(as_one):
 
 
 def _get_test_case_equality(as_one):
-    fp = FlatParams(
+    fp = InternalParams(
         values=np.array([0, 1.5, 1.5, 0, 1.5, 1]),
         lower_bounds=np.array([-10, 1, 0.9, -np.inf, -np.inf, -10]),
         upper_bounds=np.full(6, np.inf),
         names=list("abcdef"),
     )
 
-    internal = FlatParams(
+    internal = InternalParams(
         values=np.array([0, 1.5, 0, 1]),
         lower_bounds=np.array([-10, 1, -np.inf, -10]),
         upper_bounds=np.full(4, np.inf),
@@ -122,14 +122,14 @@ def _get_test_case_equality(as_one):
 
 
 def _get_test_case_probability():
-    fp = FlatParams(
+    fp = InternalParams(
         values=np.array([0.1, 0.2, 0.2, 0.5, 10]),
         lower_bounds=np.full(5, -np.inf),
         upper_bounds=np.full(5, np.inf),
         names=list("abcde"),
     )
 
-    internal = FlatParams(
+    internal = InternalParams(
         values=np.array([0.2, 0.4, 0.4, 10]),
         lower_bounds=np.array([0, 0, 0, -np.inf]),
         upper_bounds=np.full(4, np.inf),
@@ -142,14 +142,14 @@ def _get_test_case_probability():
 
 
 def _get_test_case_uncorrelated_covariance():
-    fp = FlatParams(
+    fp = InternalParams(
         values=np.array([1, 0, 4, 0, 0, 9, 10]),
         lower_bounds=np.full(7, -np.inf),
         upper_bounds=np.full(7, np.inf),
         names=list("abcdefg"),
     )
 
-    internal = FlatParams(
+    internal = InternalParams(
         values=np.array([1, 4, 9, 10]),
         lower_bounds=np.array([0, 0, 0, -np.inf]),
         upper_bounds=np.full(4, np.inf),
@@ -165,14 +165,14 @@ def _get_test_case_uncorrelated_covariance():
 
 
 def _get_test_case_covariance():
-    fp = FlatParams(
+    fp = InternalParams(
         values=np.array([1, -0.2, 1.2, -0.2, 0.1, 1.3, 0.1, -0.05, 0.2, 1, 10]),
         lower_bounds=np.full(11, -np.inf),
         upper_bounds=np.full(11, np.inf),
         names=list("abcdefghijk"),
     )
 
-    internal = FlatParams(
+    internal = InternalParams(
         values=np.array(
             [
                 1,
@@ -201,14 +201,14 @@ def _get_test_case_covariance():
 
 
 def _get_test_case_normalized_covariance():
-    fp = FlatParams(
+    fp = InternalParams(
         values=np.array([4, 0.1, 2, 0.2, 0.3, 3, 10]),
         lower_bounds=np.full(7, -np.inf),
         upper_bounds=np.full(7, np.inf),
         names=list("abcdefg"),
     )
 
-    internal = FlatParams(
+    internal = InternalParams(
         values=np.array([0.05, 1.4133294025, 0.1, 0.2087269956, 1.7165177078, 10]),
         lower_bounds=[-np.inf, 0, -np.inf, -np.inf, 0, -np.inf],
         upper_bounds=np.full(6, np.inf),
@@ -224,14 +224,14 @@ def _get_test_case_normalized_covariance():
 
 
 def _get_test_case_sdcorr():
-    fp = FlatParams(
+    fp = InternalParams(
         values=np.array([2, 1.5, 3, 0.2, 0.15, 0.33, 10]),
         lower_bounds=np.full(7, -np.inf),
         upper_bounds=np.full(7, np.inf),
         names=list("abcdefg"),
     )
 
-    internal = FlatParams(
+    internal = InternalParams(
         values=np.array([2, 0.3, 1.46969385, 0.45, 0.91855865, 2.82023935, 10]),
         lower_bounds=np.array([0, -np.inf, 0, -np.inf, -np.inf, 0, -np.inf]),
         upper_bounds=np.full(7, np.inf),
@@ -270,8 +270,8 @@ def _get_test_case_sdcorr():
 )
 def test_space_converter_with_params(constraints, params, expected_internal):
     converter, internal = get_space_converter(
-        flat_params=params,
-        flat_constraints=constraints,
+        internal_params=params,
+        internal_constraints=constraints,
     )
 
     aaae(internal.values, expected_internal.values)