5-Plots-Benchmark.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using ValueShapes\n",
    "using ArraysOfArrays\n",
    "using StatsBase \n",
    "using LinearAlgebra\n",
    "using Statistics\n",
    "using Distributions \n",
    "using IntervalSets\n",
    "\n",
    "using HCubature\n",
    "using JLD\n",
    "using CPUTime\n",
    "\n",
    "using Plots\n",
    "using Colors\n",
    "using ColorSchemes\n",
    "using LaTeXStrings\n",
    "\n",
    "using SpecialFunctions\n",
    "\n",
    "pyplot(size=(750,500))\n",
    "line_colors = ColorSchemes.tab20b;\n",
    "\n",
    "import PyPlot\n",
    "\n",
    "plt = PyPlot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SMALL_SIZE = 12\n",
    "MEDIUM_SIZE = 13\n",
    "BIGGER_SIZE = 13\n",
    "\n",
    "plt.rc(\"font\", size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc(\"axes\", titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "plt.rc(\"axes\", labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc(\"xtick\", labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc(\"ytick\", labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc(\"legend\", fontsize=SMALL_SIZE)    # legend fontsize\n",
    "plt.rc(\"figure\", titlesize=BIGGER_SIZE)  # fontsize of the figure title"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PATH_1 = \"../../AHMI_publication/NormalDistributionData/normal_dist-ffcor.jld\"\n",
    "# PATH_2 = \"../../AHMI_publication/NormalDistributionData/normal_dist_4-iid.jld\"\n",
    "# PATH_3 = \"../../AHMI_publication/NormalDistributionData/normal_dist-ffcor-2.jld\" \n",
    "\n",
    "PATH_1 = \"../../AHMI_publication/NormalDistributionData/normal_dist-ffcor.jld\"\n",
    "PATH_2 = \"../../AHMI_publication/NormalDistributionData/normal_dist-ffcor-3.jld\"\n",
    "PATH_3 = \"../../AHMI_publication/NormalDistributionData/normal_dist_4-iid.jld\"\n",
    "\n",
    "TITLE = \"AHMI Tests (Multivariate Normal Distribution)\"\n",
    "\n",
    "\n",
    "saved_1 = load(PATH_1);\n",
    "saved_2 = load(PATH_2);\n",
    "saved_3 = load(PATH_3);\n",
    "\n",
    "# Data #1 \n",
    "dim_array = saved_1[\"dim_array\"]\n",
    "integrals_ahmi = saved_1[\"integrals_ahmi_array\"]\n",
    "uns_ahmi = saved_1[\"uns_ahmi_array\"]\n",
    "integrals_true = saved_1[\"integrals_true_array\"]\n",
    "\n",
    "dim_range_1 = sort(unique(dim_array));\n",
    "integrals_1 = exp.(integrals_ahmi .- integrals_true)\n",
    "unsertanity_ahmi_1 = exp.(uns_ahmi .- integrals_true);\n",
    "unsertanity_scatter_1 = std(integrals_1, dims=1)'\n",
    "\n",
    "mean_integrals_1 = mean(integrals_1, dims=1)'\n",
    "unsertanity_ahmi_1 = mean(unsertanity_ahmi_1, dims=1)';\n",
    "\n",
    "# Data #2 \n",
    "dim_array = saved_2[\"dim_array\"]\n",
    "integrals_ahmi = saved_2[\"integrals_ahmi_array\"]\n",
    "uns_ahmi = saved_2[\"uns_ahmi_array\"]\n",
    "integrals_true = saved_2[\"integrals_true_array\"]\n",
    "\n",
    "dim_range_2 = sort(unique(dim_array));\n",
    "integrals_2 = exp.(integrals_ahmi .- integrals_true)\n",
    "unsertanity_ahmi_2 = exp.(uns_ahmi .- integrals_true);\n",
    "unsertanity_scatter_2 = std(integrals_2, dims=1)'\n",
    "\n",
    "mean_integrals_2 = mean(integrals_2, dims=1)'\n",
    "unsertanity_ahmi_2 = mean(unsertanity_ahmi_2, dims=1)';\n",
    "\n",
    "# Data #3\n",
    "dim_array = saved_3[\"dim_array\"]\n",
    "integrals_ahmi = saved_3[\"integrals_ahmi_array\"]\n",
    "uns_ahmi = saved_3[\"uns_ahmi_array\"]\n",
    "integrals_true = saved_3[\"integrals_true_array\"]\n",
    "\n",
    "dim_range_3 = sort(unique(dim_array));\n",
    "integrals_3 = exp.(integrals_ahmi .- integrals_true)\n",
    "unsertanity_ahmi_3 = exp.(uns_ahmi .- integrals_true);\n",
    "unsertanity_scatter_3 = std(integrals_3, dims=1)'\n",
    "\n",
    "mean_integrals_3 = mean(integrals_3, dims=1)'\n",
    "unsertanity_ahmi_3 = mean(unsertanity_ahmi_3, dims=1)';\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_minimum = 0.94\n",
    "y_maximum = 1.06\n",
    "\n",
    "fig, ax = plt.subplots(3,1, figsize=(13, 7), sharex=true)\n",
    "\n",
    "fig.subplots_adjust(hspace=0.0)\n",
    "\n",
    "ax[1].axhline(1, c=\"red\", label=\"Truth\")\n",
    "ax[2].axhline(1, c=\"red\", label=\"Truth\")\n",
    "ax[3].axhline(1, c=\"red\", label=\"Truth\")\n",
    "\n",
    "ax[1].plot(dim_range_1, mean_integrals_1[:,1], label=L\"$\\langle\\hat{I}\\rangle$\", color=\"k\", alpha=1)\n",
    "ax[1].fill_between(dim_range_1, mean_integrals_1[:,1] .- unsertanity_scatter_1[:,1], mean_integrals_1[:,1] .+ unsertanity_scatter_1[:,1], alpha=1, color=plt.cm.Blues(0.3), label=L\"$\\sigma(\\hat{I})$\")\n",
    "ax[1].plot(dim_range_1, mean_integrals_1[:,1] .- unsertanity_ahmi_1[:,1], alpha=1, ls=\"--\", color=\"k\",  label=L\"$\\langle\\sigma\\rangle$\")\n",
    "ax[1].plot(dim_range_1, mean_integrals_1[:,1] .+ unsertanity_ahmi_1[:,1], alpha=1, ls=\"--\", color=\"k\")\n",
    "ax[1].set_ylim(y_minimum, y_maximum)\n",
    "ax[1].legend(loc=\"upper left\", frameon=true, framealpha=0.8, ncol=5)\n",
    "\n",
    "ind = 1:20\n",
    "ax[2].plot(dim_range_2[ind], mean_integrals_2[ind,1], label=L\"t=100\", color=\"k\", alpha=1)\n",
    "ax[2].fill_between(dim_range_2[ind], mean_integrals_2[ind,1] .- unsertanity_scatter_2[ind,1], mean_integrals_2[ind,1] .+ unsertanity_scatter_2[ind,1], alpha=1, color=plt.cm.Blues(0.3), )\n",
    "ax[2].plot(dim_range_2[ind], mean_integrals_2[ind,1] .- unsertanity_ahmi_2[ind,1], alpha=1, ls=\"--\", color=\"k\",  label=L\"MCMC: $\\langle\\sigma\\rangle$\",)\n",
    "ax[2].plot(dim_range_2[ind], mean_integrals_2[ind,1] .+ unsertanity_ahmi_2[ind,1], alpha=1, ls=\"--\", color=\"k\")\n",
    "ax[2].set_ylim(y_minimum, y_maximum)\n",
    "\n",
    "ax[1].text(2.5, 0.96, \"t=1000\", fontsize=14, verticalalignment=\"top\")\n",
    "ax[2].text(2.5, 0.96, \"t=500\", fontsize=14, verticalalignment=\"top\")\n",
    "ax[3].text(2.5, 0.96, \"t=100\", fontsize=14, verticalalignment=\"top\")\n",
    "\n",
    "ax[3].plot(dim_range_3, mean_integrals_3[:,1], label=L\"t=10\", color=\"k\", alpha=1)\n",
    "ax[3].fill_between(dim_range_3, mean_integrals_3[:,1] .- unsertanity_scatter_3[:,1], mean_integrals_3[:,1] .+ unsertanity_scatter_3[:,1], alpha=1, color=plt.cm.Blues(0.3), )\n",
    "ax[3].plot(dim_range_3, mean_integrals_3[:,1] .- unsertanity_ahmi_3[:,1], alpha=1, ls=\"--\", color=\"k\",  label=L\"MCMC: $\\langle\\sigma\\rangle$\",)\n",
    "ax[3].plot(dim_range_3, mean_integrals_3[:,1] .+ unsertanity_ahmi_3[:,1], alpha=1, ls=\"--\", color=\"k\")\n",
    "ax[3].set_ylim(y_minimum, y_maximum)\n",
    "\n",
    "ax[3].set_xlabel(\"# Dimensions\")\n",
    "# ax[3].set_ylabel(\"Ratio to Truth\")\n",
    "ax[2].set_ylabel(\"Ratio to Truth\")\n",
    "# ax[1].set_ylabel(\"Ratio to Truth\")\n",
    "\n",
    "ax[1].set_xlim(2, 25)\n",
    "\n",
    "# fig.savefig(\"/home/iwsatlas1/vhafych/MPP-Project/AHMI_publication/NormalDistributionData/normal_benchmarking-b.pdf\", bbox_inches = \"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tot_volumes_accepted1_array = saved_2[\"tot_volumes_accepted1_array\"]\n",
    "tot_volumes_accepted2_array = saved_2[\"tot_volumes_accepted2_array\"]\n",
    "\n",
    "ahmi_time_array = saved_2[\"ahmi_time_array\"]\n",
    "\n",
    "tot_volumes_accepted1_array_mean = mean(tot_volumes_accepted1_array, dims=1)'\n",
    "tot_volumes_accepted2_array_mean = mean(tot_volumes_accepted2_array, dims=1)'\n",
    "\n",
    "tot_volumes_accepted1_array_std = std(tot_volumes_accepted1_array, dims=1)'\n",
    "tot_volumes_accepted2_array_std = std(tot_volumes_accepted2_array, dims=1)'\n",
    "\n",
    "volumes_sum = tot_volumes_accepted1_array_mean .+ tot_volumes_accepted2_array_mean\n",
    "volumes_std = sqrt.( tot_volumes_accepted1_array_std.^2 .+ tot_volumes_accepted2_array_std.^2)\n",
    "\n",
    "mean_time = mean(ahmi_time_array, dims=1)'\n",
    "std_time = std(ahmi_time_array, dims=1)'\n",
    "\n",
    "fig, ax = plt.subplots(1,2, figsize=(13, 4))\n",
    "\n",
    "ax[1].plot(dim_range_2[ind], mean_time[ind], color=\"k\")\n",
    "ax[1].fill_between(dim_range_2[ind], mean_time[ind,1] .- std_time[ind,1], mean_time[ind,1] .+ std_time[ind,1], alpha=1, color=plt.cm.Blues(0.3),)\n",
    "ax[1].set_xlabel(\"# Dimensions\")\n",
    "ax[1].set_ylabel(\"Integration Time, [sec.]\")\n",
    "ax[1].set_xticks(dim_range_2[1:2:end])\n",
    "ax[1].set_xlim(2, 21)\n",
    "\n",
    "ax[2].plot(dim_range_2[ind], volumes_sum[ind], color=\"k\")\n",
    "ax[2].fill_between(dim_range_2[ind], volumes_sum[ind,1] .- volumes_std[ind,1], volumes_sum[ind,1] .+ volumes_std[ind,1], alpha=1, color=plt.cm.Blues(0.3),)\n",
    "ax[2].set_xlabel(\"# Dimensions\")\n",
    "ax[2].set_ylabel(\"# Hyper Rectangles\")\n",
    "ax[2].set_xticks(dim_range_2[1:2:end])\n",
    "ax[2].set_xlim(2, 21)\n",
    "\n",
    "# fig.savefig(\"/home/iwsatlas1/vhafych/MPP-Project/AHMI_publication/NormalDistributionData/time_cubes.pdf\", bbox_inches = \"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cauchy Distribution: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# PATH_1 = \"../../AHMI_publication/CaushyDistributionData/mmod-caushy_dist_4-ffcor-2.jld\"\n",
    "# PATH_2 = \"../../AHMI_publication/CaushyDistributionData/mmod-caushy_dist_4-ffcor.jld\"\n",
    "\n",
    "# saved_1 = load(PATH_1);\n",
    "# saved_2 = load(PATH_2);\n",
    "\n",
    "# # Data #1 \n",
    "# dim_array = saved_1[\"dim_array\"]\n",
    "# integrals_ahmi = saved_1[\"integrals_ahmi_array\"]\n",
    "# uns_ahmi = saved_1[\"uns_ahmi_array\"]\n",
    "# integrals_true = saved_1[\"integrals_true_array\"]\n",
    "\n",
    "# dim_range_1 = sort(unique(dim_array));\n",
    "# integrals_1 = exp.(integrals_ahmi .- integrals_true)\n",
    "# unsertanity_ahmi_1 = exp.(uns_ahmi .- integrals_true);\n",
    "# unsertanity_scatter_1 = std(integrals_1, dims=1)'\n",
    "\n",
    "# mean_integrals_1 = mean(integrals_1, dims=1)'\n",
    "# unsertanity_ahmi_1 = mean(unsertanity_ahmi_1, dims=1)';\n",
    "\n",
    "# # Data #2 \n",
    "# dim_array = saved_2[\"dim_array\"]\n",
    "# integrals_ahmi = saved_2[\"integrals_ahmi_array\"]\n",
    "# uns_ahmi = saved_2[\"uns_ahmi_array\"]\n",
    "# integrals_true = saved_2[\"integrals_true_array\"]\n",
    "\n",
    "# dim_range_2 = sort(unique(dim_array));\n",
    "# integrals_2 = exp.(integrals_ahmi .- integrals_true)\n",
    "# unsertanity_ahmi_2 = exp.(uns_ahmi .- integrals_true);\n",
    "# unsertanity_scatter_2 = std(integrals_2, dims=1)'\n",
    "\n",
    "# mean_integrals_2 = mean(integrals_2, dims=1)'\n",
    "# unsertanity_ahmi_2 = mean(unsertanity_ahmi_2, dims=1)';\n",
    "\n",
    "\n",
    "# fig, ax = plt.subplots(1,1, figsize=(13, 6))\n",
    "\n",
    "# ax.axhline(1, c=\"red\", label=\"Truth\")\n",
    "\n",
    "# ax.plot(dim_range_1, mean_integrals_1[:,1], label=L\"$\\langle\\hat{I}\\rangle$\", color=\"k\", alpha=1)\n",
    "# ax.fill_between(dim_range_1, mean_integrals_1[:,1] .- unsertanity_scatter_1[:,1], mean_integrals_1[:,1] .+ unsertanity_scatter_1[:,1], alpha=1, color=plt.cm.Blues(0.3), label=L\"$\\sigma(\\hat{I})$\" )\n",
    "# ax.plot(dim_range_1, mean_integrals_1[:,1] .- unsertanity_ahmi_1[:,1], alpha=1, ls=\"--\", color=\"k\",  label=L\"$\\langle\\sigma\\rangle$\",)\n",
    "# ax.plot(dim_range_1, mean_integrals_1[:,1] .+ unsertanity_ahmi_1[:,1], alpha=1, ls=\"--\", color=\"k\")\n",
    "\n",
    "# # ax.plot(dim_range_2, mean_integrals_2[:,1], label=L\"$\\langle\\hat{I}\\rangle$\", color=\"k\", alpha=1)\n",
    "# # ax.fill_between(dim_range_2, mean_integrals_2[:,1] .- unsertanity_scatter_2[:,1], mean_integrals_2[:,1] .+ unsertanity_scatter_2[:,1], alpha=1, color=plt.cm.Blues(0.3), label=L\"$\\sigma(\\hat{I})$\" )\n",
    "# # ax.plot(dim_range_2, mean_integrals_2[:,1] .- unsertanity_ahmi_2[:,1], alpha=1, ls=\"--\", color=\"k\",  label=L\"$\\langle\\sigma\\rangle$\",)\n",
    "# # ax.plot(dim_range_2, mean_integrals_2[:,1] .+ unsertanity_ahmi_2[:,1], alpha=1, ls=\"--\", color=\"k\")\n",
    "\n",
    "# ax.set_xlabel(\"# Dimensions\")\n",
    "# ax.set_ylabel(\"Ratio to Truth\")\n",
    "\n",
    "# # ins = ax.inset_axes([0.2,0.4,0.3,0.4])\n",
    "# # ins.hist(integrals_1[:, 13], bins=5, color=plt.cm.Blues(0.3), alpha=0.4)\n",
    "# # ins.hist(integrals_2[:, 13], bins=5, color=plt.cm.Oranges(0.3), alpha=0.4)\n",
    "# # ins.axvline(1, c=\"red\")\n",
    "# # ins.set_xlabel(L\"Integral Estimate $\\hat{I}$\")\n",
    "# # ins.get_yaxis().set_visible(false)\n",
    "# # ins.set_title(\"# Dimensions = 13\")\n",
    "\n",
    "# ax.set_ylim(0.98, 1.02)\n",
    "# ax.set_xlim(2, 7)\n",
    "# ax.set_xticks(dim_range_1)\n",
    "\n",
    "# ax.legend(loc=\"upper left\", frameon=true, framealpha=0.8, ncol=5)\n",
    "\n",
    "# # fig.savefig(\"/home/iwsatlas1/vhafych/MPP-Project/AHMI_publication/CaushyDistributionData/caushy_benchmarking-b-3.pdf\", bbox_inches = \"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PATH_1 = \"../../AHMI_publication/CaushyDistributionData/mmod-caushy_dist_4-ffcor-2.jld\"\n",
    "# PATH_2 = \"../../AHMI_publication/CaushyDistributionData/mmod-caushy_dist_4-ffcor.jld\"\n",
    "\n",
    "PATH_1 = \"../../AHMI_publication/CaushyDistributionData/mmod-caushy_dist_4-ffcor-3.jld\"\n",
    "PATH_2 = \"../../AHMI_publication/CaushyDistributionData/mmod-caushy_dist_4-ffcor-4.jld\"\n",
    "\n",
    "saved_1 = load(PATH_1);\n",
    "saved_2 = load(PATH_2);\n",
    "\n",
    "dim_array = vcat(saved_1[\"dim_array\"][:,1:6], saved_2[\"dim_array\"][:,1:6]) \n",
    "integrals_ahmi = vcat(saved_1[\"integrals_ahmi_array\"][:,1:6], saved_2[\"integrals_ahmi_array\"][:,1:6])  \n",
    "uns_ahmi = vcat(saved_1[\"uns_ahmi_array\"][:,1:6], saved_2[\"uns_ahmi_array\"][:,1:6]) \n",
    "integrals_true = vcat(saved_1[\"integrals_true_array\"][:,1:6], saved_2[\"integrals_true_array\"][:,1:6])\n",
    "\n",
    "dim_range_1 = sort(unique(dim_array));\n",
    "integrals_1 = exp.(integrals_ahmi .- integrals_true)\n",
    "unsertanity_ahmi_1 = exp.(uns_ahmi .- integrals_true);\n",
    "unsertanity_scatter_1 = std(integrals_1, dims=1)'\n",
    "\n",
    "mean_integrals_1 = mean(integrals_1, dims=1)'\n",
    "unsertanity_ahmi_1 = mean(unsertanity_ahmi_1, dims=1)';\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(13, 6))\n",
    "\n",
    "ax.axhline(1, c=\"red\", label=\"Truth\")\n",
    "\n",
    "ax.plot(dim_range_1, mean_integrals_1[:,1], label=L\"$\\langle\\hat{I}\\rangle$\", color=\"k\", alpha=1)\n",
    "ax.fill_between(dim_range_1, mean_integrals_1[:,1] .- unsertanity_scatter_1[:,1], mean_integrals_1[:,1] .+ unsertanity_scatter_1[:,1], alpha=1, color=plt.cm.Blues(0.3), label=L\"$\\sigma(\\hat{I})$\" )\n",
    "ax.plot(dim_range_1, mean_integrals_1[:,1] .- unsertanity_ahmi_1[:,1], alpha=1, ls=\"--\", color=\"k\",  label=L\"$\\langle\\sigma\\rangle$\",)\n",
    "ax.plot(dim_range_1, mean_integrals_1[:,1] .+ unsertanity_ahmi_1[:,1], alpha=1, ls=\"--\", color=\"k\")\n",
    "\n",
    "ax.set_xlabel(\"# Dimensions\")\n",
    "ax.set_ylabel(\"Ratio to Truth\")\n",
    "\n",
    "ax.set_ylim(0.95, 1.05)\n",
    "ax.set_xlim(2, 7)\n",
    "ax.set_xticks(dim_range_1)\n",
    "\n",
    "ax.text(2.1, 0.958, \"t=500\", fontsize=14, verticalalignment=\"top\")\n",
    "\n",
    "ax.legend(loc=\"upper left\", frameon=true, framealpha=0.8, ncol=5)\n",
    "\n",
    "# fig.savefig(\"/home/iwsatlas1/vhafych/MPP-Project/AHMI_publication/CaushyDistributionData/caushy_benchmarking-b-3.pdf\", bbox_inches = \"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Funnel: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PATH_1 = \"../../AHMI_publication/FunnelDistributionData/funnel_dist_2-ffcor-3.jld\"\n",
    "PATH_2 = \"../../AHMI_publication/FunnelDistributionData/funnel_dist_2-ffcor-4.jld\"\n",
    "\n",
    "saved_1 = load(PATH_1);\n",
    "saved_2 = load(PATH_2);\n",
    "\n",
    "dim_array = vcat(saved_1[\"dim_array\"], saved_2[\"dim_array\"]) \n",
    "integrals_ahmi = vcat(saved_1[\"integrals_ahmi_array\"], saved_2[\"integrals_ahmi_array\"])  \n",
    "uns_ahmi = vcat(saved_1[\"uns_ahmi_array\"], saved_2[\"uns_ahmi_array\"]) \n",
    "integrals_true = vcat(saved_1[\"integrals_true_array\"], saved_2[\"integrals_true_array\"])\n",
    "\n",
    "dim_range_1 = sort(unique(dim_array));\n",
    "integrals_1 = exp.(integrals_ahmi .- integrals_true)\n",
    "unsertanity_ahmi_1 = exp.(uns_ahmi .- integrals_true);\n",
    "unsertanity_scatter_1 = std(integrals_1, dims=1)'\n",
    "\n",
    "mean_integrals_1 = mean(integrals_1, dims=1)'\n",
    "unsertanity_ahmi_1 = mean(unsertanity_ahmi_1, dims=1)';\n",
    "\n",
    "# Plots:\n",
    "\n",
    "fig, ax = plt.subplots(1,1, figsize=(13, 6))\n",
    "\n",
    "ax.axhline(1, c=\"red\", label=\"Truth\")\n",
    "\n",
    "ax.plot(dim_range_1, mean_integrals_1[:,1], label=L\"$\\langle\\hat{I}\\rangle$\", color=\"k\", alpha=1)\n",
    "ax.fill_between(dim_range_1, mean_integrals_1[:,1] .- unsertanity_scatter_1[:,1], mean_integrals_1[:,1] .+ unsertanity_scatter_1[:,1], alpha=1, color=plt.cm.Blues(0.3), label=L\"$\\sigma(\\hat{I})$\" )\n",
    "ax.plot(dim_range_1, mean_integrals_1[:,1] .- unsertanity_ahmi_1[:,1], alpha=1, ls=\"--\", color=\"k\",  label=L\"$\\langle\\sigma\\rangle$\",)\n",
    "ax.plot(dim_range_1, mean_integrals_1[:,1] .+ unsertanity_ahmi_1[:,1], alpha=1, ls=\"--\", color=\"k\")\n",
    "\n",
    "ax.set_xlabel(\"# Dimensions\")\n",
    "ax.set_ylabel(\"Ratio to Truth\")\n",
    "\n",
    "ax.set_ylim(0.95, 1.05)\n",
    "ax.set_xlim(2, 6)\n",
    "ax.set_xticks(dim_range_1)\n",
    "\n",
    "ax.text(2.1, 0.958, \"t=500\", fontsize=14, verticalalignment=\"top\")\n",
    "\n",
    "ax.legend(loc=\"upper left\", frameon=true, framealpha=0.8, ncol=5)\n",
    "\n",
    "# fig.savefig(\"/home/iwsatlas1/vhafych/MPP-Project/AHMI_publication/FunnelDistributionData/funnel_benchmarking-b.pdf\", bbox_inches = \"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Shell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PATH_1 = \"../../AHMI_publication/GaussShellDistributionData/gaus_shell_dist_5-ffcor-4.jld\"\n",
    "PATH_2 = \"../../AHMI_publication/GaussShellDistributionData/gaus_shell_dist_5-ffcor-4.jld\"\n",
    "\n",
    "saved_1 = load(PATH_1);\n",
    "saved_2 = load(PATH_2);\n",
    "\n",
    "# Data #1 \n",
    "dim_array = saved_1[\"dim_array\"]\n",
    "integrals_ahmi = saved_1[\"integrals_ahmi_array\"]\n",
    "uns_ahmi = saved_1[\"uns_ahmi_array\"]\n",
    "integrals_true = saved_1[\"integrals_true_array\"]\n",
    "\n",
    "dim_range_1 = sort(unique(dim_array));\n",
    "integrals_1 = exp.(integrals_ahmi .- integrals_true)\n",
    "unsertanity_ahmi_1 = exp.(uns_ahmi .- integrals_true);\n",
    "unsertanity_scatter_1 = std(integrals_1, dims=1)'\n",
    "\n",
    "mean_integrals_1 = mean(integrals_1, dims=1)'\n",
    "unsertanity_ahmi_1 = mean(unsertanity_ahmi_1, dims=1)';\n",
    "\n",
    "# Data #2 \n",
    "dim_array = saved_2[\"dim_array\"]\n",
    "integrals_ahmi = saved_2[\"integrals_ahmi_array\"]\n",
    "uns_ahmi = saved_2[\"uns_ahmi_array\"]\n",
    "integrals_true = saved_2[\"integrals_true_array\"]\n",
    "\n",
    "dim_range_2 = sort(unique(dim_array));\n",
    "integrals_2 = exp.(integrals_ahmi .- integrals_true)\n",
    "unsertanity_ahmi_2 = exp.(uns_ahmi .- integrals_true);\n",
    "unsertanity_scatter_2 = std(integrals_2, dims=1)'\n",
    "\n",
    "mean_integrals_2 = mean(integrals_2, dims=1)'\n",
    "unsertanity_ahmi_2 = mean(unsertanity_ahmi_2, dims=1)';\n",
    "\n",
    "# Make plot: \n",
    "\n",
    "fig, ax = plt.subplots(1,1, figsize=(13, 6))\n",
    "\n",
    "ax.axhline(1, c=\"red\", label=\"Truth\")\n",
    "\n",
    "ax.plot(dim_range_2, mean_integrals_2[:,1], label=L\"$\\langle\\hat{I}\\rangle$\", color=\"k\", alpha=1)\n",
    "ax.fill_between(dim_range_2, mean_integrals_2[:,1] .- unsertanity_scatter_2[:,1], mean_integrals_2[:,1] .+ unsertanity_scatter_2[:,1], alpha=1, color=plt.cm.Blues(0.3), label=L\"$\\sigma(\\hat{I})$\" )\n",
    "ax.plot(dim_range_2, mean_integrals_2[:,1] .- unsertanity_ahmi_2[:,1], alpha=1, ls=\"--\", color=\"k\",  label=L\"$\\langle\\sigma\\rangle$\",)\n",
    "ax.plot(dim_range_2, mean_integrals_2[:,1] .+ unsertanity_ahmi_2[:,1], alpha=1, ls=\"--\", color=\"k\")\n",
    "\n",
    "ax.set_xlabel(\"# Dimensions\")\n",
    "ax.set_ylabel(\"Ratio to Truth\")\n",
    "\n",
    "ax.set_ylim(0.95, 1.05)\n",
    "ax.set_xlim(2, 17)\n",
    "\n",
    "ax.set_xticks(dim_range_2)\n",
    "\n",
    "ax.text(2.3, 0.958, \"t=500\", fontsize=14, verticalalignment=\"top\")\n",
    "\n",
    "ax.legend(loc=\"upper left\", frameon=true, framealpha=0.8, ncol=5)\n",
    "\n",
    "# fig.savefig(\"/home/iwsatlas1/vhafych/MPP-Project/AHMI_publication/GaussShellDistributionData/shell_benchmarking-b.pdf\", bbox_inches = \"tight\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.4.0-rc1",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}