Commit d399dc90 authored by Maximilian Schanner's avatar Maximilian Schanner
Browse files

Further work on the examples.

parent d94247a4
This source diff could not be displayed because it is too large. You can view the blob instead.
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Exploration\n",
"The first step in the `CORBASS` algorithm is an exploration of the parameter space. This way the region of interest, i.e. the region where the probability mass is concentrated, is found. Here we give an example how to run the `corbass.exploration` module and how to interpret the output. First we import the function `explore` from the module:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import sys\n",
"import os\n",
"# relative import\n",
"sys.path.append(os.path.abspath('') + '/../')\n",
"from corbass.exploration import explore"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After the import we load the example parameter file:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from corbass.utils import load\n",
"pars = load('./Example_Parfile.py')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we use the `explore` function to explore the parameter space. This is done via grid evaluation. The model parameter bounds and the number of gridpoints are specified in the parameter file. Then for each point in the grid, i.e. each choice of model parameters, the posterior probability is calculated. `explore` then returns an instance of `ExplorationResult`, containing the posterior on the grid and a scale. We use the `write` member function to write the results to the disk, to use them in [the next example](http://localhost:8888/notebooks/Example_Integration.ipynb)."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
" 3% (1 of 27) | | Elapsed Time: 0:00:00 ETA: 0:00:03"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exploring 1450, this may take a while...\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100% (27 of 27) |########################| Elapsed Time: 0:00:04 Time: 0:00:04\n",
"N/A% (0 of 27) | | Elapsed Time: 0:00:00 ETA: --:--:--"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exploring 1550, this may take a while...\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100% (27 of 27) |########################| Elapsed Time: 0:00:04 Time: 0:00:04\n"
]
}
],
"source": [
"for bin_name, data in pars.data:\n",
" print(f\"Exploring {bin_name}, this may take a while...\")\n",
" rslt = explore(pars.bin_fname(bin_name), data, pars.n_ex, pars.bounds,\n",
" r_ref=pars.r_ref)\n",
" rslt.write()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can load the results for one of the bins and analyse them. The standard output name can be accessed from the `ExplorationResult` class, the output location is specified in the parameter file:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"from corbass.exploration import ExplorationResult\n",
"import numpy as np\n",
"\n",
"with np.load(f'{pars.bin_fname(1450)}{ExplorationResult.suffix}') as fh:\n",
" posterior = fh['posterior']\n",
" scale = fh['scale']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To demonstrate how we use this exploration, we focus on the marginal posterior for the error level and the residual scaling. First we get the corresponding parameter space arrays and build a meshgrid:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"from corbass.integration import bounds2grid\n",
"nd, err, res = bounds2grid(pars.bounds, pars.n_ex)\n",
"ER, RE = np.meshgrid(err[0], res[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we integrate out the non-dipole scaling, by means of a simple Riemann sum:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"marginal = posterior.sum(axis=0)*nd[1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can plot the marginal distribution:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.QuadMesh at 0x7f5054c56160>"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD4CAYAAAAD6PrjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAQVUlEQVR4nO3df6zV9X3H8eergEhsiZJWx7hstQnrhs1qlSCNSdPNZpJuGf4xE/5YJY0JmXFLmyxZtH9s6X/9q9n8QxfSdmLWzRDbTmJqHWVt9g+VQquliIy7uimBybqlLa4LFffeH+fjJyeXC/fce9F77uX5SL453/M+n8+5nw8fwovvj3NuqgpJkgDesdADkCSND0NBktQZCpKkzlCQJHWGgiSpW77QA5jJVVlZV3PNQg9DusCv/ebPF3oI0kUd/sG5H1fVe2bbb+xD4Wqu4bbcsdDDkC7wzDPPL/QQpItatvbEv8+ln6ePJEmdoSBJ6gwFSVJnKEiSOkNBktQZCpKkzlCQJHWGgiSpMxQkSZ2hIEnqDAVJUmcoSJI6Q0GS1BkKkqTOUJAkdSOFQpJrkzyR5MUkx5J8OMmaJPuSnGiP1w21fzDJZJLjSe4cqt+a5Eh77aEkeSsmJUmam1GPFP4K+EZV/TrwQeAY8ACwv6o2APvbc5JsBLYDNwFbgYeTLGvv8wiwE9jQtq2XaR6SpMtgxlBIshr4CPBFgKr6RVX9BNgG7G7NdgN3tf1twONVda6qXgImgc1J1gKrq+pAVRXw2FAfSdIYGOVI4X3AfwJ/k+T7Sb6Q5Brghqo6DdAer2/t1wGvDPU/2Wrr2v7U+gWS7ExyKMmh1zk3qwlJkuZulFBYDtwCPFJVHwL+h3aq6CKmu05Ql6hfWKzaVVWbqmrTClaOMERJ0uUwSiicBE5W1bPt+RMMQuLVdkqI9nhmqP36of4TwKlWn5imLkkaEzOGQlX9B/BKkve30h3AC8BeYEer7QCebPt7ge1JVia5kcEF5YPtFNPZJFvaXUf3DPWRJI2B5SO2+xPgy0muAn4EfJJBoOxJci/wMnA3QFUdTbKHQXCcB+6vqjfa+9wHPAqsAp5umyRpTGRwI9D4Wp01dVvuWOhhSBd45tTzCz0E6aKWrT1xuKo2zbafn2iWJHWGgiSpMxQkSZ2hIEnqDAVJUmcoSJI6Q0GS1BkKkqTOUJAkdYaCJKkzFCRJnaEgSeoMBUlSZyhIkjpDQZLUGQqSpM5QkCR1hoIkqTMUJEmdoSBJ6gwFSVJnKEiSOkNBktQZCpKkzlCQJHWGgiSpMxQkSZ2hIEnqDAVJUmcoSJI6Q0GS1I0UCkn+LcmRJM8lOdRqa5LsS3KiPV431P7BJJNJjie5c6h+a3ufySQPJcnln5Ikaa5mc6TwW1V1c1Vtas8fAPZX1QZgf3tOko3AduAmYCvwcJJlrc8jwE5gQ9u2zn8KkqTLZT6nj7YBu9v+buCuofrjVXWuql4CJoHNSdYCq6vqQFUV8NhQH0nSGBg1FAr4xySHk+xstRuq6jRAe7y+1dcBrwz1Pdlq69r+1PoFkuxMcijJodc5N+IQJUnztXzEdrdX1akk1wP7krx4ibbTXSeoS9QvLFbtAnYBrM6aadtIki6/kY4UqupUezwDfA3YDLzaTgnRHs+05ieB9UPdJ4BTrT4xTV2SNCZmDIUk1yR515v7wO8APwT2Ajtasx3Ak21/L7A9ycokNzK4oHywnWI6m2RLu+vonqE+kqQxMMrpoxuAr7W7R5cDf1dV30jyXWBPknuBl4G7AarqaJI9wAvAeeD+qnqjvdd9wKPAKuDptkmSxkQGNwKNr9VZU7fljoUehnSBZ049v9BDkC5q2doTh4c+QjAyP9EsSeoMBUlSZyhIkjpDQZLUGQqSpM5QkCR1hoIkqTMUJEmdoSBJ6gwFSVJnKEiSOkNBktQZCpKkzlCQJHWGgiSpMxQkSZ2hIEnqDAVJUmcoSJI6Q0GS1BkKkqTOUJAkdYaCJKkzFCRJnaEgSeoMBUlSZyhIkjpDQZLUGQqSpM5QkCR1I4dCkmVJvp/kqfZ8TZJ9SU60x+uG2j6YZDLJ8SR3DtVvTXKkvfZQklze6UiS5mM2RwqfAo4NPX8A2F9VG4D97TlJNgLbgZuArcDDSZa1Po8AO4ENbds6r9FLki6rkUIhyQTwu8AXhsrbgN1tfzdw11D98ao6V1UvAZPA5iRrgdVVdaCqCnhsqI8kaQyMeqTwl8CfAf83VLuhqk4DtMfrW30d8MpQu5Ottq7tT61fIMnOJIeSHHqdcyMOUZI0XzOGQpLfA85U1eER33O66wR1ifqFxapdVbWpqjatYOWIP1aSNF/LR2hzO/D7ST4OXA2sTvK3wKtJ1lbV6XZq6ExrfxJYP9R/AjjV6hPT1CVJY2LGI4WqerCqJqrqvQwuIP9TVf0hsBfY0ZrtAJ5s+3uB7UlWJrmRwQXlg+0U09kkW9pdR/cM9ZEkjYFRjhQu5nPAniT3Ai8DdwNU1dEke4AXgPPA/VX1RutzH/AosAp4um2SpDGRwY1A42t11tRtuWOhhyFd4JlTzy/0EKSLWrb2xOGq2jTbfn6iWZLUGQqSpM5QkCR1hoIkqTMUJEndfG5Jla5od/7yBxd6CNIlnJhTL48UJEmdoSBJ6gwFSVJnKEiSOkNBktQZCpKkzlCQJHWGgiSpMxQkSZ2hIEnqDAVJUmcoSJI6Q0GS1BkKkqTOUJAkdYaCJKkzFCRJnaEgSeoMBUlSZyhIkjpDQZLUGQqSpM5QkCR1hoIkqTMUJEndjKGQ5OokB5M8n+Roks+2+pok+5KcaI/XDfV5MMlkkuNJ7hyq35rkSHvtoSR5a6YlSZqLUY4UzgG/XVUfBG4GtibZAjwA7K+qDcD+9pwkG4HtwE3AVuDhJMvaez0C7AQ2tG3rZZyLJGmeZgyFGnitPV3RtgK2AbtbfTdwV9vfBjxeVeeq6iVgEticZC2wuqoOVFUBjw31kSSNgZGuKSRZluQ54Aywr6qeBW6oqtMA7fH61nwd8MpQ95Ottq7tT61P9/N2JjmU5NDrnJvNfCRJ8zBSKFTVG1V1MzDB4H/9H7hE8+muE9Ql6tP9vF1VtamqNq1g5ShDlCRdBrO6+6iqfgJ8m8G1gFfbKSHa45nW7CSwfqjbBHCq1SemqUuSxsQodx+9J8m1bX8V8DHgRWAvsKM12wE82fb3AtuTrExyI4MLygfbKaazSba0u47uGeojSRoDy0dosxbY3e4gegewp6qeSnIA2JPkXuBl4G6AqjqaZA/wAnAeuL+q3mjvdR/wKLAKeLptkqQxkcGNQONrddbUbbljoYchSYvKN+uJw1W1abb9/ESzJKkzFCRJnaEgSeoMBUlSZyhIkjpDQZLUGQqSpM5QkCR1hoIkqTMUJEmdoSBJ6gwFSVJnKEiSOkNBktQZCpKkzlCQJHWGgiSpMxQkSZ2hIEnqDAVJUmcoSJI6Q0GS1BkKkqTOUJAkdYaCJKkzFCRJnaEgSeoMBUlSZyhIkjpDQZLUzRgKSdYn+VaSY0mOJvlUq69Jsi/JifZ43VCfB5NMJjme5M6h+q1JjrTXHkqSt2ZakqS5GOVI4Tzwp1X1G8AW4P4kG4EHgP1VtQHY357TXtsO3ARsBR5Osqy91yPATmBD27ZexrlIkuZpxlCoqtNV9b22fxY4BqwDtgG7W7PdwF1tfxvweFWdq6qXgElgc5K1wOqqOlBVBTw21EeSNAZmdU0hyXuBDwHPAjdU1WkYBAdwfWu2DnhlqNvJVlvX9qfWp/s5O5McSnLodc7NZoiSpHkYORSSvBP4CvDpqvrZpZpOU6tL1C8sVu2qqk1VtWkFK0cdoiRpnkYKhSQrGATCl6vqq638ajslRHs80+ongfVD3SeAU60+MU1dkjQmRrn7KMAXgWNV9fmhl/YCO9r+DuDJofr2JCuT3MjggvLBdorpbJIt7T3vGeojSRoDy0doczvwCeBIkuda7TPA54A9Se4FXgbuBqiqo0n2AC8wuHPp/qp6o/W7D3gUWAU83TZJ0pjI4Eag8bU6a+q23LHQw5CkReWb9cThqto0235+olmS1BkKkqTOUJAkdYaCJKkzFCRJnaEgSeoMBUlSZyhIkjpDQZLUGQqSpM5QkCR1hoIkqTMUJEmdoSBJ6gwFSVJnKEiSOkNBktQZCpKkzlCQJHWGgiSpMxQkSZ2hIEnqDAVJUmcoSJI6Q0GS1BkKkqTOUJAkdYaCJKkzFCRJnaEgSeoMBUlSN2MoJPlSkjNJfjhUW5NkX5IT7fG6odceTDKZ5HiSO4fqtyY50l57KEku/3QkSfMxypHCo8DWKbUHgP1VtQHY356TZCOwHbip9Xk4ybLW5xFgJ7ChbVPfU5K0wGYMhar6Z+C/p5S3Abvb/m7grqH641V1rqpeAiaBzUnWAqur6kBVFfDYUB9J0phYPsd+N1TVaYCqOp3k+lZfB3xnqN3JVnu97U+tTyvJTgZHFQCvfbOeOD7HcS6kdwM/XuhBvA2uhHleCXME57nUvH8uneYaChcz3XWCukR9WlW1C9h1uQa1EJIcqqpNCz2Ot9qVMM8rYY7gPJeaJIfm0m+udx+92k4J0R7PtPpJYP1QuwngVKtPTFOXJI2RuYbCXmBH298BPDlU355kZZIbGVxQPthONZ1NsqXddXTPUB9J0piY8fRRkr8HPgq8O8lJ4C+AzwF7ktwLvAzcDVBVR5PsAV4AzgP3V9Ub7a3uY3An0yrg6bYtZYv69NcsXAnzvBLmCM5zqZnTPDO4GUiSJD/RLEkaYihIkjpDYZ6SbG1f6TGZ5IFpXv9okp8mea5tf74Q45yP6b7qZMrraV9dMpnkB0luebvHOF8jzHHRryNAkvVJvpXkWJKjST41TZulsJ6jzHNRr2mSq5McTPJ8m+Nnp2kz+7WsKrc5bsAy4F+B9wFXAc8DG6e0+Sjw1EKPdZ7z/AhwC/DDi7z+cQY3DgTYAjy70GN+C+a46NexzWMtcEvbfxfwL9P8nV0K6znKPBf1mrb1eWfbXwE8C2yZ71p6pDA/m4HJqvpRVf0CeJzBV30sKTX9V50M2wY8VgPfAa5983Msi8UIc1wSqup0VX2v7Z8FjnHhtwsshfUcZZ6LWluf19rTFW2beufQrNfSUJifdcArQ88v9vUdH26HeE8nuentGdrbatQ/h8VuSa1jkvcCH2LwP8xhS2o9LzFPWORrmmRZkucYfIB4X1XNey0v99dcXGlG+fqO7wG/WlWvJfk48A8MPtS3lMzqa0wWqSW1jkneCXwF+HRV/Wzqy9N0WZTrOcM8F/2a1uBzYDcnuRb4WpIPVNXwdbFZr6VHCvNzsa/16KrqZ28e4lXV14EVSd799g3xbTHjn8Nit5TWMckKBv9QfrmqvjpNkyWxnjPNcymtaVX9BPg2F/5KglmvpaEwP98FNiS5MclVDH6XxN7hBkl+6c1fKJRkM4M/8/9620f61toL3NPudNgC/LTat+guFUtlHdscvggcq6rPX6TZol/PUea52Nc0yXvaEQJJVgEfA16c0mzWa+npo3moqvNJ/hh4hsGdSF+qwVd9/FF7/a+BPwDuS3Ie+F9ge7XbAhaLTP9VJyugz/HrDO5ymAR+DnxyYUY6dyPMcdGvY3M78AngSDsXDfAZ4Fdg6awno81zsa/pWmB3Br/I7B3Anqp6asq/P7NeS7/mQpLUefpIktQZCpKkzlCQJHWGgiSpMxQkSZ2hIEnqDAVJUvf/U8Pf2nCoaLAAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"from matplotlib import pyplot as plt\n",
"plt.pcolormesh(ER, RE, marginal.T)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"By calculating the empirical mean and covariance of the individual components, we find the region where the probability mass is concentrated. Therefore the function `get_integ_bounds` from `corbass.integration` is available."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"ename": "ValueError",
"evalue": "Something went wrong exploring eps. eps_var == 0",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-8-d4c0edc03001>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mcorbass\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mintegration\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mget_integ_bounds\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m integ_bounds, moments = get_integ_bounds(pars.bounds, posterior,\n\u001b[0;32m----> 3\u001b[0;31m moments_out=True)\n\u001b[0m",
"\u001b[0;32m~/Documents/CORBASS/examples/../corbass/integration.py\u001b[0m in \u001b[0;36mget_integ_bounds\u001b[0;34m(expl_bounds, posterior, moments_out)\u001b[0m\n\u001b[1;32m 242\u001b[0m \u001b[0meps_var\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mepss\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mmarg_eps\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mepss\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0meps_mean\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 243\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0meps_var\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0.0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 244\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Something went wrong exploring eps. eps_var == 0'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 245\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 246\u001b[0m \u001b[0mmarg_rho\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mposterior\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mValueError\u001b[0m: Something went wrong exploring eps. eps_var == 0"
]
}
],
"source": [
"from corbass.integration import get_integ_bounds\n",
"integ_bounds, moments = get_integ_bounds(pars.bounds, posterior,\n",
" moments_out=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For illustration purposes, we show the marginal distribution for the error scaling and the residual, together with the Gaussian from the empirical mean and variance."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from scipy.stats import norm\n",
"fig, ax = plt.subplots(1, 2)\n",
"\n",
"ax[0].plot(err[0], marginal.sum(axis=1)*res[1])\n",
"ax[0].plot(err[0], norm.pdf(err[0], loc=moments[1, 0], scale=np.sqrt(moments[1, 1])))\n",
"ax[0].set_xlabel('err')\n",
"ax[0].set_ylabel('[Arb.]')\n",
"ax[0].set_yticks([])\n",
"\n",
"ax[1].plot(res[0], marginal.sum(axis=0)*err[1])\n",
"ax[1].plot(res[0], norm.pdf(res[0], loc=moments[2, 0], scale=np.sqrt(moments[2, 1])))\n",
"ax[1].set_xlabel('res')\n",
"ax[1].set_ylabel('[Arb.]')\n",
"ax[1].set_yticks([])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally we draw the found area of interest on top of the marginal distribution, to show wich area would be integrated over in the following step of `CORBASS`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from matplotlib.patches import Rectangle\n",
"fig, ax = plt.subplots()\n",
"ax.pcolormesh(ER, RE, marginal.T)\n",
"\n",
"rect = Rectangle((moments[1, 0] - 2*np.sqrt(moments[1, 1]),\n",
" moments[2, 0] - 2*np.sqrt(moments[2, 1])),\n",
" 4*np.sqrt(moments[1, 1]), 4*np.sqrt(moments[2, 1]), \n",
" fill=True, facecolor=(1, 1, 1, 0.2),\n",
" lw=2, ec='white')\n",
"ax.add_patch(rect)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
This source diff could not be displayed because it is too large. You can view the blob instead.
This diff is collapsed.
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -5,7 +5,7 @@ path = Path(__file__).parent # relative to the location of this very par file
# a prefix for this parameter file. Serves as identifier for output as well
# as a path to output locations
prefix = '../out/default'
prefix = '../out/example'
prefix = Path(path, prefix)
# specify the data file
......@@ -15,7 +15,7 @@ data_filename = Path(path, data_filename) # Relative to the parfile
# set up binning
# optional: specify initial and final year for the binning
initial = 1400
final = 1500
final = 1600
binlength = 100 # Specify (equal) bin width
# specify other parameters for the algorithm
......@@ -32,9 +32,9 @@ rho_min = 1000 # minimal residual
rho_max = 6500 # maximal residual
# number of gridpoints per dimension...
n_ex = 15 # ...in the exploration step
n_ex = 3 # ...in the exploration step
n_fine = 4 # ...in the integration step
# parameters for prediction
n_desi = 10 # number of field design points
n_desi = 1000 # number of field design points
l_max = 10 # maximal SH degree
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment