{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Line-of-sight sampling\n", "This notebook starts from a posterior distribution in the angular diameter distances infered from a strong gravitational lens time-delay analysis (TDSL) that has not yet been corrected for the line-of-sight contribution. The notebook reads-in the files and post-processes the external convergences and then saves those files as the final posteriors." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import os\n", "import pickle\n", "\n", "import matplotlib.pyplot as plt\n", "from matplotlib import rc\n", "import matplotlib.patches as mpatches\n", "\n", "%matplotlib inline\n", "\n", "cwd = os.getcwd()\n", "base_path, _ = os.path.split(cwd)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEf9JREFUeJzt3X+MpdVdx/H3R2gx1sYusq4ItLNt1j9oolRHaKw/UJRfVamxaSDarpVkmwpJjf7h1sZgapqg8Uc0USraVZqoFG21m3YVt9j64w9aBkRgqciU0rCbLWwLtmgNBvz6xz2Ll+nMzp25M/femfN+JTf3uec5z51z5t6Zzz3PeZ7npqqQJPXna6bdAEnSdBgAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE6dPu0GnMpZZ51Vc3Nz026GJG0pd9999xeqaudq9WY6AObm5lhYWJh2MyRpS0nyuVHquQtIkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdWrVAEhyXpKPJ3kwyZEk72jlv5LkWJJ72+3KoW3emWQxyUNJLhsqv7yVLSbZvzldkiSNYpTzAJ4FfqGq7knyUuDuJIfbut+uqt8YrpzkfOBq4NXAtwAfS/KtbfXvAT8MHAXuSnKwqh7ciI5IktZm1QCoquPA8bb8dJJPA+ecYpOrgFur6hngs0kWgQvbusWqegQgya2trgEgSVOwpjOBk8wBrwE+CbwOuD7JW4AFBqOEpxiEw51Dmx3l/wPjsSXlFy3zM/YB+wBe/vKXr6V5mpC5/R99fvnRG18/xZZIGsfIk8BJvh74IPBzVfVl4CbgVcAFDEYIv7kRDaqqm6tqvqrmd+5c9VIWkqR1GmkEkORFDP75/2lVfQigqh4fWv+HwEfaw2PAeUObn9vKOEW5JGnCRjkKKMD7gE9X1W8NlZ89VO3HgQfa8kHg6iRnJNkN7AE+BdwF7EmyO8mLGUwUH9yYbkiS1mqUEcDrgDcD9ye5t5X9EnBNkguAAh4F3gZQVUeS3MZgcvdZ4Lqqeg4gyfXA7cBpwIGqOrKBfdEmGt7vL2l7GOUooH8GssyqQ6fY5j3Ae5YpP3Sq7SRJkzPT3weg6fJTv7S9GQDaMB4eKm0tXgtIkjrlCEBjcTeRtHU5ApCkThkAktQpA0CSOmUASFKnDABJ6pRHAWlTeE6ANPscAUhSpwwASeqUu4D0Ap7YJfXDEYAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpU54Ipk3ndYGk2eQIQJI6ZQBIUqcMAEnqlHMAnfPib1K/HAFIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkTq0aAEnOS/LxJA8mOZLkHa38zCSHkzzc7ne08iT53SSLSe5L8h1Dz7W31X84yd7N65YkaTWjjACeBX6hqs4HXgtcl+R8YD9wR1XtAe5ojwGuAPa02z7gJhgEBnADcBFwIXDDydCQJE3eqgFQVcer6p62/DTwaeAc4CrgllbtFuANbfkq4P01cCfwsiRnA5cBh6vqyap6CjgMXL6hvZEkjWxNcwBJ5oDXAJ8EdlXV8bbq88CutnwO8NjQZkdb2UrlkqQpGDkAknw98EHg56rqy8PrqqqA2ogGJdmXZCHJwokTJzbiKSVJyxgpAJK8iME//z+tqg+14sfbrh3a/ROt/Bhw3tDm57aylcpfoKpurqr5qprfuXPnWvoiSVqDUY4CCvA+4NNV9VtDqw4CJ4/k2Qt8eKj8Le1ooNcCX2q7im4HLk2yo03+XtrK1JG5/R99/iZpuka5GujrgDcD9ye5t5X9EnAjcFuSa4HPAW9q6w4BVwKLwFeAtwJU1ZNJfhW4q9V7d1U9uSG9kCSt2aoBUFX/DGSF1ZcsU7+A61Z4rgPAgbU0UJK0OTwTWJI6ZQBIUqcMAEnqlF8J2SGPwJEEjgAkqVsGgCR1ygCQpE4ZAJLUKSeBNTXDk9GP3vj6KbZE6pMjAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSp/xCGM0cvyhGmgxHAJLUKQNAkjplAEhSpwwASeqUk8CaCcMTv5ImwxGAJHVq1QBIciDJE0keGCr7lSTHktzbblcOrXtnksUkDyW5bKj88la2mGT/xndFkrQWo4wA/gS4fJny366qC9rtEECS84GrgVe3bX4/yWlJTgN+D7gCOB+4ptWVJE3JqnMAVfWPSeZGfL6rgFur6hngs0kWgQvbusWqegQgya2t7oNrbrHWxX3skpYaZw7g+iT3tV1EO1rZOcBjQ3WOtrKVyiVJU7LeALgJeBVwAXAc+M2NalCSfUkWkiycOHFio55WkrTEugKgqh6vqueq6n+BP+T/d/McA84bqnpuK1upfLnnvrmq5qtqfufOnetpniRpBOsKgCRnDz38ceDkEUIHgauTnJFkN7AH+BRwF7Anye4kL2YwUXxw/c2WJI1r1UngJH8OXAycleQocANwcZILgAIeBd4GUFVHktzGYHL3WeC6qnquPc/1wO3AacCBqjqy4b2RJI1slKOArlmm+H2nqP8e4D3LlB8CDq2pdZKkTeOZwJLUKQNAkjrlxeA00/x2MGnzOAKQpE4ZAJLUKQNAkjplAEhSpwwASeqURwFpy1h6SWuPCpLG4whAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcrDQLexpYdNStIwRwCS1CkDQJI6ZQBIUqcMAEnqlJPA2rL8tjBpPI4AJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktSpVQMgyYEkTyR5YKjszCSHkzzc7ne08iT53SSLSe5L8h1D2+xt9R9OsndzuiNJGtUoI4A/AS5fUrYfuKOq9gB3tMcAVwB72m0fcBMMAgO4AbgIuBC44WRoSBthbv9Hn79JGs2qAVBV/wg8uaT4KuCWtnwL8Iah8vfXwJ3Ay5KcDVwGHK6qJ6vqKeAwXx0qkqQJWu8cwK6qOt6WPw/sasvnAI8N1TvaylYq/ypJ9iVZSLJw4sSJdTZPkrSasSeBq6qA2oC2nHy+m6tqvqrmd+7cuVFPK0laYr0B8HjbtUO7f6KVHwPOG6p3bitbqVySNCXrDYCDwMkjefYCHx4qf0s7Gui1wJfarqLbgUuT7GiTv5e2MmnDOSEsjWbV7wRO8ufAxcBZSY4yOJrnRuC2JNcCnwPe1KofAq4EFoGvAG8FqKonk/wqcFer9+6qWjqxLEmaoFUDoKquWWHVJcvULeC6FZ7nAHBgTa2TJG0azwSWpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTq54HoK3Fs18ljcoRgCR1ygCQpE4ZAJLUKecAtK0Nz4k8euPrp9gSafY4ApCkThkAktQpA0CSOuUcgLrhfID0Qo4AJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlOcBqEueEyA5ApCkbhkAktQpA0CSOmUASFKnnATeBvwieEnr4QhAkjplAEhSpwwASerUWHMASR4FngaeA56tqvkkZwIfAOaAR4E3VdVTSQL8DnAl8BXgp6vqnnF+vrTRPEFMPdmISeAfqKovDD3eD9xRVTcm2d8e/yJwBbCn3S4Cbmr30lQ5ia5ebcYuoKuAW9ryLcAbhsrfXwN3Ai9LcvYm/HxJ0gjGDYAC/i7J3Un2tbJdVXW8LX8e2NWWzwEeG9r2aCt7gST7kiwkWThx4sSYzZMkrWTcXUDfU1XHknwTcDjJvw2vrKpKUmt5wqq6GbgZYH5+fk3bShvJ+QBtd2ONAKrqWLt/Avgr4ELg8ZO7dtr9E636MeC8oc3PbWWSpClYdwAkeUmSl55cBi4FHgAOAntbtb3Ah9vyQeAtGXgt8KWhXUWSpAkbZxfQLuCvBkd3cjrwZ1X1t0nuAm5Lci3wOeBNrf4hBoeALjI4DPStY/xsSdKY1h0AVfUI8O3LlH8RuGSZ8gKuW+/Pk6bJ+QBtR54JLEmdMgAkqVMGgCR1yu8DkMbg3IC2MkcAktQpRwDSGnnxOG0XjgAkqVOOALYoP4VKGpcjAEnqlAEgSZ1yF5C0QU61W85DRDWLHAFIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkTnkYqDQBXjVUs8gAkCbMMNCsMAC2EK//s70ZDJo0A0CaolFC3WDQZjEApBm0UjAYBtpIHgUkSZ0yACSpU+4CkrYZdxNpVAaAtEV5VJjGZQBI25ijAZ2KASB1YumIwUCQASB1aqVdSAZDPwyAGed+Xs0KdydtPwaApBcY50OHIbG1GACS1mycf/SGxOyYeAAkuRz4HeA04I+q6sZJt2HWudtHW8lar2c0Sp2lwWBobI5U1eR+WHIa8O/ADwNHgbuAa6rqweXqz8/P18LCwsTaNysMAGntVgqGHsMjyd1VNb9avUmPAC4EFqvqEYAktwJXAcsGgCSNaqNGIuMaJ2QmHVaTDoBzgMeGHh8FLppwG2aGn/Sl7Wcr/V3P3CRwkn3AvvbwP5M8NMbTnQV8YfxWTd126QfYl1m1XfqyXfpBfm2svrxilEqTDoBjwHlDj89tZc+rqpuBmzfihyVZGGU/2KzbLv0A+zKrtktftks/YDJ9mfTloO8C9iTZneTFwNXAwQm3QZLEhEcAVfVskuuB2xkcBnqgqo5Msg2SpIGJzwFU1SHg0IR+3IbsSpoB26UfYF9m1Xbpy3bpB0ygLxM9D0CSNDv8SkhJ6tSWD4AkZyY5nOThdr9jhXp/m+Q/knxkSfnuJJ9MspjkA21yeuLW0I+9rc7DSfYOlX8iyUNJ7m23b5pc659vw+WtDYtJ9i+z/oz2O15sv/O5oXXvbOUPJblsku1ear39SDKX5L+HXoP3TrrtS43Ql+9Lck+SZ5O8ccm6Zd9r0zJmX54bel2mfuDJCH35+SQPJrkvyR1JXjG0buNel6ra0jfg14H9bXk/8Gsr1LsE+FHgI0vKbwOubsvvBd4+q/0AzgQeafc72vKOtu4TwPwUX4fTgM8ArwReDPwrcP6SOj8LvLctXw18oC2f3+qfAexuz3PaFuzHHPDAtF6DdfZlDvg24P3AG0d5r221vrR1/znt12ONffkB4Ova8tuH3mMb+rps+REAg0tJ3NKWbwHesFylqroDeHq4LEmAHwT+crXtJ2CUflwGHK6qJ6vqKeAwcPmE2rea5y/zUVX/A5y8zMew4T7+JXBJew2uAm6tqmeq6rPAYnu+aRinH7Nm1b5U1aNVdR/wv0u2nbX32jh9mTWj9OXjVfWV9vBOBudMwQa/LtshAHZV1fG2/Hlg1xq2/UbgP6rq2fb4KIPLVUzDKP1Y7lIaw+394zbE/eUp/ENarW0vqNN+519i8BqMsu2kjNMPgN1J/iXJPyT53s1u7CrG+b3O0msC47fna5MsJLkzybQ+5J201r5cC/zNOrc9pZm7FMRyknwM+OZlVr1r+EFVVZKZPaxpk/vxk1V1LMlLgQ8Cb2YwFNbkHAdeXlVfTPKdwF8neXVVfXnaDROvaH8frwT+Psn9VfWZaTdqNUl+CpgHvn8znn9LBEBV/dBK65I8nuTsqjqe5GzgiTU89ReBlyU5vX2S+6pLU2ykDejHMeDiocfnMtj3T1Uda/dPJ/kzBsPMSQbAqpf5GKpzNMnpwDcweA1G2XZS1t2PGuykfQagqu5O8hngW4FpXdN8nN/riu+1KRnrPTL09/FIkk8Ar2GwH34aRupLkh9i8OHw+6vqmaFtL16y7SfW25DtsAvoIHByJnwv8OFRN2x/sB8HTh4xsKbtN9go/bgduDTJjnaU0KXA7UlOT3IWQJIXAT8CPDCBNg8b5TIfw318I/D37TU4CFzdjq7ZDewBPjWhdi+17n4k2ZnBd17QPmnuYTBJNy3jXHpl2ffaJrVzFOvuS+vDGW35LOB1TPcS9Kv2JclrgD8Afqyqhj8MbuzrMu0Z8XFvDPa93gE8DHwMOLOVzzP4xrGT9f4JOAH8N4P9Zpe18lcy+GezCPwFcMaM9+NnWlsXgbe2spcAdwP3AUdo37g2hT5cyeALfz4DvKuVvbu9iQG+tv2OF9vv/JVD276rbfcQcMWU31Pr6gfwE+33fy9wD/Cj0+zHiH35rvb38F8MRmNHTvVe24p9Ab4buJ/B0Tb3A9dugb58DHi8vZfuBQ5uxuvimcCS1KntsAtIkrQOBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ36P6djMi1CrmE+AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "-0.003989534930386585 0.005210819641515942 0.041154649764267845\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuYXFWZ7/HvW1XdnXuApAmXBDoh4RJA0NMGdRDRjBpEiZegyaMjahwOapxnBh2McwQZxnGAMyOORzyYEUaEo4CZM9pKNB4IkhEhpoFwSSLYuUBukE7nnk6nu6re80ftdCqV6tTurl237t/nefKwa+9VVe/uhF+tXnvXWubuiIjI0BCrdAEiIlI+Cn0RkSFEoS8iMoQo9EVEhhCFvojIEKLQFxEZQhT6IiJDiEJfRGQIUeiLiAwhiUoXkGv8+PHe1NRU6TJERGrK008/vcPdGwu1q7rQb2pqorW1tdJliIjUFDN7JUw7De+IiAwhCn0RkSFEoS8iMoSECn0zm2VmL5lZm5ktzHO8wcweDI6vMLOmYP/HzWxV1p+0mV0c7SmIiEhYBUPfzOLAncAVwHRgnplNz2k2H9jl7lOBO4DbANz9/7j7xe5+MfAXwAZ3XxXlCYiISHhhevozgDZ3X+/u3cADwOycNrOBe4PtxcBMM7OcNvOC54qISIWECf3TgU1ZjzcH+/K2cfcksAcYl9PmY8BPBlamiIhEoSwXcs3sEqDT3V/s4/i1ZtZqZq3t7e3lKKnfupNpHmrdhJaXFJFaFib0twCTsh5PDPblbWNmCWAs0JF1fC7H6eW7+yJ3b3b35sbGgl8oq4jvPPonblj8PEteeK3SpYiIDFiY0F8JTDOzyWZWTybAW3LatADXBNtzgGUedInNLAZ8lBoez0+nne8+1gbA/kM9Fa5GRGTgCoZ+MEa/AFgKrAUecvfVZnaLmV0VNLsbGGdmbcD1QPZtnZcBm9x9fbSll8+W3Qd7t/d1JdlzUMEvIrXJqm2Murm52att7p2tuw/ytluXHbVv461XVqgaEZFjmdnT7t5cqJ2+kRtCPJZ796mISG1S6IdwzDcORERqlEI/hHS60hWIiERDoR/CN5esrXQJIiKRUOiH8Pzm3ZUuQUQkEgr9EHpS1XWHk4jIQCn0Q+hJaVBfRAYHhX4IybR6+iIyOCj0Q1BPX0QGC4V+CEmN6YvIIKHQDyGpG/VFZJBQ6Iegu3dEZLBQ6A/Qy6/vq3QJIiL9ptAfoPfcsbzSJYiI9JtCX0RkCFHoi4gMIQp9EZEhRKFfgCZbE5HBRKFfwHX3PV3pEkREIqPQLyCmpRJFZBAJFfpmNsvMXjKzNjNbmOd4g5k9GBxfYWZNWcfeYGZPmtlqM3vBzIZFV37pJRT6IjKIFAx9M4sDdwJXANOBeWY2PafZfGCXu08F7gBuC56bAO4HrnP384HLgZ7Iqi+DRFy/DInI4BEm0WYAbe6+3t27gQeA2TltZgP3BtuLgZlmZsB7gOfd/TkAd+9w91Q0pZeHevoiMpiECf3TgU1ZjzcH+/K2cfcksAcYB5wNuJktNbNnzOyGfG9gZteaWauZtba3t/f3HErqj69pugURGTxKPXaRAC4FPh7890NmNjO3kbsvcvdmd29ubGwscUkD15A4+seV0uIqIlJjwoT+FmBS1uOJwb68bYJx/LFAB5nfCpa7+w537wSWAG8qtuhKmTfjDDbeemXvY026JiK1JkzorwSmmdlkM6sH5gItOW1agGuC7TnAMnd3YClwoZmNCD4M3gGsiab08jpnwmjmXzr5qH1pV09fRGpLolADd0+a2QIyAR4H7nH31WZ2C9Dq7i3A3cB9ZtYG7CTzwYC77zKzb5H54HBgibs/XKJzKYlxI+uZdcEp/OOHLjzmmDJfRGpNwdAHcPclZIZmsvfdlLXdBVzdx3PvJ3PbZk1Kufd5B88Pf7+Rf776ojJXJCIycLoJvYBUyvv8Vu7PV+Ve2hARqW4K/QJS7sQtf+jHdQ+/iNQYhf5xPL95N53dKeLx/OGeiOnHJyK1Ral1HFd99wmAPnv6Y4fXlbMcEZGiKfRD6OtC7rwZk/LuFxGpVgr9EDS9sogMFgr9EPoa3nmodXOZKxERKY5CP4QHWzcd9fjEEZmx/Fd3dlaiHBGRAVPoh3DGSSOOevzrv76sQpWIiBRHoR/C7XPecNTjCWOOLP6lmTZFpJYo9ENoHN3Q57FFy9eXsRIRkeIo9ENoSMT7PPa9x9rKWImISHEU+kXadyhZ6RJEREILNcvmUJWIGW8648RKlyEiEhn19I9j3Kh6pjSOrHQZIiKRUej3wd3Z35XE+vhilohILVLo9+E7j7ZxoDvFC1t2V7oUEZHIKPT7cMcjLwPw4pa9Fa5ERCQ6Cn0RkSEkVOib2Swze8nM2sxsYZ7jDWb2YHB8hZk1BfubzOygma0K/twVbfnVIa1v5YpIjSh4y6aZxYE7gXcDm4GVZtbi7muyms0Hdrn7VDObC9wGfCw4ts7dL4647pKLGYTN8pQ7MXTBV0SqX5ie/gygzd3Xu3s38AAwO6fNbODeYHsxMNNq/LaXWD/K1/w7IlIrwoT+6UD23MKbg31527h7EtgDjAuOTTazZ83scTN7e5H1lk2hGL/s7Mbe7bQr9EWkNpT6Qu424Ax3fyNwPfBjMxuT28jMrjWzVjNrbW9vL3FJ4Rzuvc+bcUbe47d/5MjMm0n19EWkRoQJ/S1A9mKwE4N9eduYWQIYC3S4+yF37wBw96eBdcDZuW/g7ovcvdndmxsbG3MPV9SEMfln2KxPHPnR6UKuiNSKMKG/EphmZpPNrB6YC7TktGkBrgm25wDL3N3NrDG4EIyZTQGmATUxF/GH35gZwfrgxbkjWRnZI/4a0xeRWlHw7h13T5rZAmApEAfucffVZnYL0OruLcDdwH1m1gbsJPPBAHAZcIuZ9QBp4Dp331mKE4naGeMyq2WdOW5E3uPZMZ/SmL6I1IhQs2y6+xJgSc6+m7K2u4Cr8zzvP4D/KLLGijjcee/rJiTPCvp0uhwViYgUT9/I7YO7EzvOXZuJ+JEfnXr6IlIrFPp9SLsf9179scPrmHnuyZm2GtMXkRqh0O9D2gt/QevKN5wKQE9K4zsiUhsU+n1Iu1PoS7lbdh0E4F3/8ngZKhIRKZ5Cvw8eoqffcaC7TNWIiERDod+HdPr4F3IBhtXFy1OMiEhEFPp9CDOmH9dPT0RqjGKrD2l3Cs2WrDl3RKTWKPT7kHYnUWB8J5VS6ItIbVHo92H5y+3sP5Q8bptrL5sCwHmnjmHTzs5ylCUiUhSFfh5dPSk2dnTSU6Anf/KYYbxlykms3baXt9/+GH/YUBPTConIEKbQz6M/sypkX+xd176/BNWIiERHoZ9Hf1bC+v26jt7teG2vECkiQ4BCP4+BLn8YK3Rjv4hIhSn08xjonZi6b19Eqp1iKo+BzppZ6MtcIiKVptDPY6DDO3EN74hIlVPo5zHQ4R319EWk2in08zjc0y/0jVyAJ7/6rt5t9fRFpNqFCn0zm2VmL5lZm5ktzHO8wcweDI6vMLOmnONnmNl+M/tyNGWX1v1PvQLANz54QcG2wxJHZtqsT+gzVESqW8GUMrM4cCdwBTAdmGdm03OazQd2uftU4A7gtpzj3wJ+VXy55fG/lrUB4YZrGuoU9CJSO8Ik1gygzd3Xu3s38AAwO6fNbODeYHsxMNMsk5hm9kFgA7A6mpLLJ8wQfUNWT9+1QLqIVLkwoX86sCnr8eZgX9427p4E9gDjzGwU8BXg74svtfzCjNFnt0lrqVwRqXKlHpu4GbjD3Y87KY2ZXWtmrWbW2t7eXuKSwuvv3Tg/fXoTO7WEoohUsTChvwWYlPV4YrAvbxszSwBjgQ7gEuB2M9sI/DXwd2a2IPcN3H2Ruze7e3NjY2O/T6JU+nsH5tLVr/PFnzxTmmJERCIQJvRXAtPMbLKZ1QNzgZacNi3ANcH2HGCZZ7zd3ZvcvQn4NvBNd/9uRLWX3KFkuPGaH3/2kt7t7XsPlaocEZGiJQo1cPdk0DtfCsSBe9x9tZndArS6ewtwN3CfmbUBO8l8MNS8VMhvaTWObihxJSIi0SgY+gDuvgRYkrPvpqztLuDqAq9x8wDqq6iwN+NY1jiQvpQrItVMN5kfR108XILrm7giUisU+scxPuSwzfC6I/fqG/oAEJHqpdA/jsvPDncn0ciGeOFGIiJVQKF/HBZygH5k/ZFLIxrTF5FqptCPgJZJFJFaodCPyNkTRlW6BBGRghT6EdECKiJSCxT6EVHoi0gtCPXlrKHm9BOGc9Gksf16ju7VF5FaoJ5+HrHY0fPkh6GOvojUAoV+Hpt2HqQn1b/J8ZX5IlILFPo5tuw+CMAvn9/WvycGXf2w9/aLiFSCQj/Hns6eAT1PUS8itUChn2OgHfXDz1P4i0g1U+hHJBHcvbNm295+Xw8QESkXhX5EsmfXVOiLSLVS6EcklbXiStgVt0REyk2hn2OgY/o3vPec3u20OvoiUqUU+jkGugjKJVPG9W4nlfoiUqUU+hH69J81ARreEZHqFSr0zWyWmb1kZm1mtjDP8QYzezA4vsLMmoL9M8xsVfDnOTP7ULTlR88ZeGCfM2E0AH/avj+qckREIlUw9M0sDtwJXAFMB+aZ2fScZvOBXe4+FbgDuC3Y/yLQ7O4XA7OA75tZVU/yVkwv/fCkax//wYqoyhERiVSYnv4MoM3d17t7N/AAMDunzWzg3mB7MTDTzMzdO909GewfBkV0o8vEi6hQM22KSLULE/qnA5uyHm8O9uVtE4T8HmAcgJldYmargReA67I+BKpSFD19EZFqVfILue6+wt3PB94MfNXMhuW2MbNrzazVzFrb29tLXdJxpYvo6hfzXBGRcggT+luASVmPJwb78rYJxuzHAh3ZDdx9LbAfuCD3Ddx9kbs3u3tzY2Nj+OpLoJgbb3pSCn0RqW5hQn8lMM3MJptZPTAXaMlp0wJcE2zPAZa5uwfPSQCY2ZnAucDGSCovkcO99Xs+1dzv59bHdQesiFS3gnfSuHvSzBYAS4E4cI+7rzazW4BWd28B7gbuM7M2YCeZDwaAS4GFZtYDpIHPu/uOUpxIVNJBV7+/K2cBXDixf0ssioiUW6jbJ919CbAkZ99NWdtdwNV5nncfcF+RNZZVMgj9ugH02s9qHMXsi09j5YadUZclIhIJjUfk6A5myEzEB3Ynzoj6BFv3dPHNJWujLEtEJBIK/RzJ4GLsQMfnN+/qBGDR8vWR1SQiEhWFfo5kkT397qQmWxOR6qXQz9G+/xAwsDH9Yp4nIlIOSqgcN/18NXBk+cP+0gybIlLNFPoRy55L/7cvba9gJSIix1Lo5/jwGzPTCp05buSAnp/9rdxfPLctkppERKKi0M9RF49xyphjpgcKLbunX8zc/CIipaDQz9GTSlOXGPhsmRdPOuHIA2W+iFQZhX6OnrRTFxv4j+XG9x9ZX0aZLyLVRqGfI5lKD/gefcjM2TNhTAOgqZZFpPoo9HP0pNIkiujpA8Qt86GhzBeRaqPQz7Grs4cTRtQV9RrxIn5TEBEpJYV+jl0HujlxZH1Rr6Fv5YpItVI65ehJp2koMrS//oHzAWh5bitrtu6NoiwRkUgo9HOk0xArcoHzd5x9ZMnH36+r6jVjRGSIUejnSKbTA553J594hK8lIlIshX6OVAQ9/WxRfoCIiBRLoZ8jFXFPP8oPEBGRYoUKfTObZWYvmVmbmS3Mc7zBzB4Mjq8ws6Zg/7vN7GkzeyH477uiLT96qbQTs+iCWvfqi0g1KRj6ZhYH7gSuAKYD88xsek6z+cAud58K3AHcFuzfAXzA3S8ErqEGFklPpT3Snr4yX0SqSZie/gygzd3Xu3s38AAwO6fNbODeYHsxMNPMzN2fdfetwf7VwHAza4ii8FJJuUd68fXk0VV9uiIyxIQJ/dOBTVmPNwf78rZx9ySwBxiX0+YjwDPufmhgpZZHKh1t6G/f2xXZa4mIFKssF3LN7HwyQz7/vY/j15pZq5m1tre3l6OkPvWkogn93/zNZQDcGCy/KCJSDcKE/hZgUtbjicG+vG3MLAGMBTqCxxOB/wQ+6e7r8r2Buy9y92Z3b25sbMzXpCw27ewEYHdnT9GvpbVyRaQahQn9lcA0M5tsZvXAXKAlp00LmQu1AHOAZe7uZnYC8DCw0N2fiKroUtnblQn7t52VOzLVf4eS6cKNRETKrGDoB2P0C4ClwFrgIXdfbWa3mNlVQbO7gXFm1gZcDxy+rXMBMBW4ycxWBX9OjvwsItIdBPWwunjRrzV6WKJ3e8f+qr6MISJDSKJwE3D3JcCSnH03ZW13AVfned43gG8UWWPZHF7UPIpZMs9qHNW7/fLr+xg/SnfxiEjl6Ru5WQ739OsT0fxYzpkwGoCO/d2RvJ6ISLEU+ln2H8qM6UcV+j/8zJsBzbQpItVDoZ/luvufAaA+okVQTh07nJNG1tPVo4u6IlIdFPp5RNXTB4LQT7F198HIXlNEZKAU+nk0RBj6DYkYv3rxNd5267Le7wGIiFSKQj+PKHv62R8gr2tKBhGpMIV+HlGN6QOcMnZY77a+oysilabQz6Muwp7+hDFZoa/UF5EKU+jnMbK++G/kHrZq0+7IXktEpFgK/TwswpWzNu3UXTsiUj0U+jn+aua0SF/vO/Mu7t3evEt374hIZSn0A+lgKuSo1zF/65QjM3Ze/9Bz0b64iEg/KfQDqeAqazzCoR2IdqhIRKRYCv1AOgj9WNRdfeCP/zAr8tcUERkIhX7gkTXbAXhg5auRv3aU9/2LiBRDaRQ4fJG1FHfbZP/2cLA7Ffnri4iEpdAPJMrUG1/ywrayvI+ISD4K/cBJI+sAOO/UMSV9ny/9VHfwiEjlKPQDY4dnQv+W2eeX5PXnvnlS7/b/W/M67fu0bq6IlJ9CP5AM1scdHsGi6PncMOvc3u2//FErb/7HR0ryPiIixxMq9M1slpm9ZGZtZrYwz/EGM3swOL7CzJqC/ePM7DEz229m34229Gilgi9nJeKlua9+zLBQa9CLiJRUwdA3szhwJ3AFMB2YZ2bTc5rNB3a5+1TgDuC2YH8XcCPw5cgqLpGew6Ffgvv0IXOh+KqLTivJa4uIhBWmpz8DaHP39e7eDTwAzM5pMxu4N9heDMw0M3P3A+7+OzLhX9VS6cw6tvFY6Ua8PnPp5KMeP7W+o2TvJSKST5iEOx3YlPV4c7Avbxt3TwJ7gHGEZGbXmlmrmbW2t7eHfVqkepKZnn5diYZ3AE4cUXfU45UbdpbsvURE8qmKC7nuvsjdm929ubGxsSI17D+UBGBUQ+nG3k8ePeyox92pdMneS0QknzChvwWYlPV4YrAvbxszSwBjgZoauzgQhP7IEob+8JzFWQ4lFfoiUl5hQn8lMM3MJptZPTAXaMlp0wJcE2zPAZa519bigPu7k9QnYtSV+Ju5X//AkWvgnd3Jkr6XiEiuggkXjNEvAJYCa4GH3H21md1iZlcFze4GxplZG3A90Htbp5ltBL4FfMrMNue586fi3J3vP76ecnxOffrPjlzMvf+pV3tvFRURKYdQYxnuvgRYkrPvpqztLuDqPp7bVER9ZfHMq7sA6EmVJ4Anjx/Jhh0HAPjXR17m+vecU5b3FRGpigu5lZYsU9gfln1/0HeWtfVeRIbMbx0a9hGRUlHoA3WJ8v4Y/j5nfp//fPbIdfHvLmtj+k1L2XWgu6w1icjQoNCn/IucvH1a37el/mxV5gOg44AmZBOR6Cn0gcPXb6eMH1m29/zFgkt7t18JxvdFREpNoQ8kgykYbvxA+W4sunDi2N7tH/xuA3u7egDQvTwiUkoKfbJm2CzRZGt92Xjrlb3b/7z0Jbbv7cpK/fLWIiJDg+b7BZJB6MfLHPrZfvTkK/zoyVcq9v4iMjSopw/8aft+AOJW/tD/wjvPyrs/XVtfaBaRGqHQB2782YsAvLa3/DNA/+17z827f9ueqp+NWkRqkEIfOKsxc9fOjMknVeT9/3Xuxcfsu+aeP/DH1/b2XuAVEYmCQh+YeOII3jBxLKeOHV6R95998eksvu6tx+yf9e3/4jP/vrICFYnIYKXQB/Yc7GHs8LrCDUuouekkvvXRi47Z3/rKrgpUIyKDlUKf6gh9gA+/aWKlSxCRQU6hD+zu7OaEEZUPfYAbZh074+YPn9jAPo3ti0gEhnzoJ1NpdnX2kCjhguj98fnLp/K1K887at/Nv1jDhTf/hkPJVIWqEpHBojqSroKe27wHoCwLqIT12bdPYf0333fM/nO+9mu6ehT8IjJwQz70Dw+bfOCi0ypcydFiMWPjrVdy1yf+21H7z73x11z7o1aaFj7MgUNJHn5+W1V9YIlIdRvyob/46c0AnDSyvsKV5DfrglP4aPPRF3h/s+Z1AM7/+lK+8ONnWLFhZyVKE5EaNOTn3unY381JI+uZXMZplfvr9jkXcfuci/j1i69x3f1PH3N87qKnADjv1DHcMOscRjUkeHTtdr4y6xysj6kl5v9wJRt2HGDZly8vZekiUmVC9fTNbJaZvWRmbWa2MM/xBjN7MDi+wsyaso59Ndj/kpm9N7rSi5dOO2tf28u7zj25z3CsJrMuOIUN//Q+fvvly/niu6Yec3zttr18+t9XcvVdT3LX4+uY/NUl/GHDTq655w80LXyYe363obfto3/cznrN4y8y5BQMfTOLA3cCVwDTgXlmljvx/Hxgl7tPBe4AbgueOx2YC5wPzAK+F7xeVfi/z25hd2cPb50yrtKlhGZmNI0fyZfecw7P3vhuHrn+HZw6dhinn5D/28Qf/f6TPP5yOwC3/HINTQsfpmnhw73H12zdy12Pr+PVjk56Uune/as27ebD33uCF7fs4d+Wr2fWt5cDhdfwfbWjkx/813pdZxCpUlbof04zeytws7u/N3j8VQB3/6esNkuDNk+aWQJ4DWgEFma3zW7X1/s1Nzd7a2trUScVxqsdnXz4fz/Bjv3dPHPju6t2TL8/lrywjYdf2Mb4kfWs2ryH5zbtjvT1P/nWM/ld2w7Wtx/g3FNG8/Zp43l0beY3hq/MOpfdnd18f/l6AG58/3ROGF7HU+s7eOe5J1Mfj9GdSnP+aWPYeaCbKY2jGFkfJx4znlzfwZhhdUw/dQyxrOmtz/q7JbzznJP5wTXNbNxxgHXt+zmrcRRNWUNx7k5XTxozaEjE+P26Dt521jjMjFTacXcS8RgHu1OYwbC6gfU5kqk0+7qSnNiPfyfuftzfIPd29VAfjw24poE4Xk3JVJpEmZcOleiY2dPu3lywXYjQnwPMcvfPBo//ArjE3RdktXkxaLM5eLwOuAS4GXjK3e8P9t8N/MrdF/f1fgMN/T++tpcFP34Wd8c9sxZJunfbCRbH6t13eEbNz11+Fl+ZlX+my8Hg9b1dbN97iIknDmf11r088+oufvn8VjZ2dNKdTDNhTAOv782sx2t2ZOnIconHrHcRm3jMepeOqU/E6OzO3J6aiFnvmgcAoxsSmEFdPEZHHwvI1ydi9KTSxM0YP6qh9+97/Kh66uIxYv0cztuy+yAAw+vivR2Ew+svxGPGwe4Uw+vjGLCrs5vO7hQxMyaMaegzZDcEw2tTggn/jlo2zcIto9PVk/mHXRc39nUlj/mSYfZ77+vqYc/BHk7L81vh3oM97NjfTdO4EQMO/u5kmkTMjvrglv65/OxGvvb+ga3gFzb0q+JCrpldC1wLcMYZZwzoNYYl4pw9YRRmmeAwM2LB/zgxy2zEeo/BhDHDOP+0sbz3/AmRnUc1mjBmGBPGDAPg0mnjuXTaeP5q5rSCzzv84dmdSnMomeZQMsW67QeoT8R4decBXti8l9NOGMYpY4eRSjtPrusgFjOaxo1gR3BxfH37fhoSceoTMXZ1dtOQiHPmuBEsf7mdU8YMY/zoBhIxY9Wm3Qyri3NW40jq4jHM4GB3mkfWvo4Z/Pl5E3ho5Sb2HUrSfOaJTJswinjMSDu8vqeLtvb9jGpIMPXkUfx81VY+1jyJeNxoe30/ibgx8cThLH95B9v3dfGOs08mZv1flnL5y+1s33eId0+fQF08hrv3rnmQcqiPxziUTGFmbNt9kNNOGM6wulhvKB/z8wW2Bu2mnzqmd7+ZZX722Q2Pk6EHDiVJxIyGRDwznUh26Oec5K7ObjbsOHDU+x22ryvJ4y+3M/20MdgAV207FIR+JRcjqnWn9jFMG6UhO7wjIjKYhO3ph/k9biUwzcwmm1k9mQuzLTltWoBrgu05wDLPfJq0AHODu3smA9OAP4Q9CRERiVbB4R13T5rZAmApEAfucffVZnYL0OruLcDdwH1m1gbsJPPBQNDuIWANkAS+4O6aR0BEpEIKDu+Um4Z3RET6L8rhHRERGSQU+iIiQ4hCX0RkCFHoi4gMIQp9EZEhpOru3jGzduCVStcR0nhgR6WLKKHBfH46t9o1mM+vmHM7090bCzWqutCvJWbWGuYWqVo1mM9P51a7BvP5lePcNLwjIjKEKPRFRIYQhX5xFlW6gBIbzOenc6tdg/n8Sn5uGtMXERlC1NMXERlCFPohFLMwfLULcW7Xm9kaM3vezB41szMrUedAFTq/rHYfMTM3s5q5KyTMuZnZR4O/v9Vm9uNy1zhQIf5dnmFmj5nZs8G/zfdVos6BMLN7zGx7sOJgvuNmZt8Jzv15M3tTpAVkVkjSn77+kJlOeh0wBagHngOm57T5PHBXsD0XeLDSdUd4bu8ERgTbn6uVcwt7fkG70cBy4CmgudJ1R/h3Nw14FjgxeHxypeuO8NwWAZ8LtqcDGytddz/O7zLgTcCLfRx/H/ArMmumvQVYEeX7q6df2Aygzd3Xu3s38AAwO6fNbODeYHsxMNOOtyJ29Sh4bu7+mLt3Bg+fAiaWucZihPm7A/gH4Dagq5zFFSnMuf0lcKe77wJw9+1lrnGgwpybA4fXfRwLbC1jfUVx9+Vk1h3py2zgR57xFHCCmZ0a1fsr9As7HdiU9XhzsC9vG3dPAnuAcWWprjhhzi3bfDI9kFpR8PyCX50nufvD5SwsAmH+7s4GzjazJ8zsKTObVbbqihPm3G4GPmFmm4ElwBfLU1pZ9Pf/y35fNsAEAAABoklEQVSpioXRpfqZ2SeAZuAdla4lKmYWA74FfKrCpZRKgswQz+VkfkNbbmYXuvvuilYVjXnAD939X4J1vO8zswvcPf9K9NJLPf3CtgCTsh5PDPblbRMsDD8W6ChLdcUJc26Y2Z8D/wO4yt0Plam2KBQ6v9HABcBvzWwjmfHTlhq5mBvm724z0OLuPe6+AXiZzIdAtQtzbvOBhwDc/UlgGJl5awaDUP9fDpRCv7BiFoavdgXPzczeCHyfTODXypjwYcc9P3ff4+7j3b3J3ZvIXLO4yt1rYb3OMP8uf0aml4+ZjScz3LO+nEUOUJhzexWYCWBm55EJ/fayVlk6LcAng7t43gLscfdtUb24hncK8CIWhq92Ic/tfwKjgJ8G16ZfdferKlZ0P4Q8v5oU8tyWAu8xszVACvhbd6/630BDntuXgH8zs78hc1H3UzXS0cLMfkLmw3h8cE3i60AdgLvfReYaxfuANqAT+HSk718jPycREYmAhndERIYQhb6IyBCi0BcRGUIU+iIiQ4hCX0RkCFHoi4gMIQp9EZEhRKEvIjKE/H/5HNVF6ZGiZQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "0.00536026591197984\n" ] } ], "source": [ "# ====================================\n", "# LENS REDSHIFTS AND DEFAULT COSMOLOGY\n", "# ====================================\n", "\n", "# redshifts\n", "z_d = 0.745\n", "z_s = 1.789\n", "\n", "\n", "# =========================\n", "# LINE-OF-SIGHT CONVERGENCE\n", "# =========================\n", "\n", "# this is the convergence distribution used in Birrer et al. 2018\n", "\n", "kappa_file = 'kappahist_J1206_measured_5innermask_nobeta_zgap-1.0_-1.0_fiducial_120_gal_120_zoverr_45_gal_45_zoverr_24_med_increments4_4_4_4.cat'\n", "path2kappa = os.path.join(base_path, 'DATA', 'LOS', kappa_file)\n", "output = np.loadtxt(path2kappa)\n", "kappa_list = np.linspace(-0.1, 1, len(output))\n", "pdf_list = output\n", "from lenstronomy.Util.prob_density import Approx\n", "pdf_approx = Approx(kappa_list, pdf_list)\n", "pdf_draw = pdf_approx.draw(n=50000)\n", "plt.hist(pdf_draw, bins=np.linspace(-0.1, 0.2, 100))\n", "plt.show()\n", "print(np.median(pdf_draw), np.mean(pdf_draw), np.std(pdf_draw))\n", "\n", "plt.plot(kappa_list, output)\n", "plt.show()\n", "\n", "mean_kappa = np.sum(kappa_list*output) / np.sum(output)\n", "print(mean_kappa)\n", "std_kappa = 'test'\n", "\n", "def draw_kappa(n=1):\n", " return pdf_approx.draw(n=n)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "33848\n", "1230.47625397\n", "5854.735165912209\n", "1804.2437530184063\n", "67696\n", "-4.547473508864641e-13\n" ] } ], "source": [ "# we import posterior samples from the analysis pre-LOS\n", "input_name = 'test.txt' # pickle file with [Dd_ds_dds_samples, D_d_samples, kappa_pert] samples from combined analysis of imaging, kinematic and time delays\n", "# and assign the post-processed posterior a name\n", "processed_name = 'test'\n", "\n", "# we saved the following distributions (you can process more if you like):\n", "# - 'angular_diameter_pre_LOS.txt': the final BIC selected sample, processed name: 'final'\n", "# - 'angular_diameter_pre_LOS_all.txt': combination of all models regardless of BIC values, processed name: 'marginalize_all'\n", "# - 'angular_diameter_pre_LOS_power_law.txt': BIC selection of the power-law profiles, processed name: 'final_power_law'\n", "# - 'angular_diameter_pre_LOS_composite.txt': BIC selection of the composite profiles, processed name: 'final_composite'\n", "\n", "file_name = os.path.join(base_path, 'Posteriors', input_name)\n", "f = open(file_name)\n", "[Dd_ds_dds_samples, D_d_samples, kappa_pert] = pickle.load(f)\n", "f.close()\n", "n_enhance = 2 # draws from each posterior n_enhance numbers from the kappa_ext distribution, this ensures a smooth final posterior distribution for a good KDE likelihood computation\n", "Dd_ds_dds_samples = np.repeat(Dd_ds_dds_samples, n_enhance)\n", "print(len(D_d_samples))\n", "D_d_samples = np.repeat(D_d_samples, n_enhance)\n", "kappa_pert = np.repeat(kappa_pert, n_enhance)\n", "\n", "# set kappa_pert = 0 if you do not agree with the subtraction\n", "#kappa_pert = 0\n", "\n", "kappa_ext = pdf_approx.draw(n = len(Dd_ds_dds_samples)) - kappa_pert\n", "Dd_ds_dds_samples /= (1 - kappa_ext)\n", "D_dt_samples = Dd_ds_dds_samples * (1 + z_d)\n", " \n", "np.random.seed(1206)\n", "\n", "D_dt_samples_save = D_dt_samples[(D_dt_samples >= 0) & (D_d_samples >= 0)]\n", "D_d_samples_save = D_d_samples[(D_dt_samples >= 0) & (D_d_samples >= 0)]\n", "np.save(name+'_D_dt.npy', D_dt_samples_save)\n", "np.save(name+'_D_d.npy', D_d_samples_save)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.6" } }, "nbformat": 4, "nbformat_minor": 2 }