Commit 784d5977 authored by Maximilian Schanner's avatar Maximilian Schanner
Browse files

Reran the tests so the figures and outputs appear.

parent b1b777e6
......@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
"execution_count": null,
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
......@@ -25,7 +25,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
......@@ -45,7 +45,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
......@@ -81,7 +81,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
......@@ -111,7 +111,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
......@@ -171,9 +171,22 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 6,
"metadata": {},
"outputs": [],
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAE9CAYAAAA8gKerAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAYX0lEQVR4nO3dfYxfV33n8ffHM+PHxE6aGDuJbRxYAwW2ATrNggpVEAkNESLAwhK0UlkqrTdVEe3uVlu0SFC26j6r2qUggleNCiyFdMUmhCU8kzTJqiGxswnNI3WchBi7seMkdvycsb/7h3+GwczY45nfmTszeb+k0fzuved3z9dXkf3JOfeem6pCkiRJbczrugBJkqS5zLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDQ12XUAL5557bq1du7brMiRJ0gvEpk2bnqqq5WMdm5Nha+3atWzcuLHrMiRJ0gtEksfHO9bpNGKSa5PsSHLfOMcvSbI7yT29n49Nd42SJElT0fXI1l8AnwI+f5I2t1XV26enHEmSpP7qdGSrqm4Fnu6yBkmSpJZmw9OIb0hyb5JvJHlV18VIkiSdjq6nEU/lbuDFVbU3yRXADcC6sRomWQ+sB1izZs30VShJknQSM3pkq6r2VNXe3uebgKEk547TdkNVDVfV8PLlYz55KUmSNO1mdNhKsjJJep8v5li9u7qtSpIkaeI6nUZM8iXgEuDcJFuBjwNDAFV1DfAe4HeSjAAHgKuqqjoqV5Ik6bR1Graq6v2nOP4pji0NIUmSNCvN6GlESZKk2c6wJUmS1JBhS5IkqSHDliRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1JBhS5IkqSHDliRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKmhTsNWkmuT7Ehy3zjHk+STSTYn+WGS1013jZIkSVPR9cjWXwCXn+T424B1vZ/1wGemoSZJkqS+6TRsVdWtwNMnaXIl8Pk65g7grCTnTU91kiRJU9f1yNapXAA8MWp7a2/fL0iyPsnGJBt37tw5LcVJkiSdykwPWxljX43VsKo2VNVwVQ0vX768cVmSJEkTM9PD1lZg9ajtVcC2jmqRJEk6bTM9bN0I/FbvqcTXA7uranvXRUmSJE3UYJedJ/kScAlwbpKtwMeBIYCquga4CbgC2AzsBz7YTaWSJEmT02nYqqr3n+J4Ab87TeVIkiT13UyfRpQkSZrVDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1JBhS5IkqSHDliRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDXUatpJcnuThJJuTfGSM45ck2Z3knt7Px7qoU5IkabIGu+o4yQDwaeAyYCtwV5Ibq+qBE5reVlVvn/YCJUmS+qDLka2Lgc1VtaWqDgNfBq7ssB5JkqS+6zJsXQA8MWp7a2/fid6Q5N4k30jyqukpTZIkqT86m0YEMsa+OmH7buDFVbU3yRXADcC6MU+WrAfWA6xZs6afdUqSJE1alyNbW4HVo7ZXAdtGN6iqPVW1t/f5JmAoybljnayqNlTVcFUNL1++vFXNkiRJp6XLsHUXsC7JhUnmA1cBN45ukGRlkvQ+X8yxendNe6WSJEmT1Nk0YlWNJPkQ8C1gALi2qu5PcnXv+DXAe4DfSTICHACuqqoTpxolSZJmrMzF7DI8PFwbN27sugxJkvQCkWRTVQ2PdcwV5CVJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1JBhS5IkqSHDliRJUkOGLUmSpIYGT3YwyZ5TfD/A9qp6Wf9KkiRJmnke3bmX2x/ZxZO7D7Ji2ULe+NJzuHD5Gaf83qlGth6pqqUn+TkT2NeXP0EfPbnnIF+443Ee3bm361IkSdIc8OjOvVy3cSv7Do6wctlC9h0c4bqNWyeUNU4Vtv7xBPqfSJtpNTQw77QugiRJ0snc/sguzlo0xNJFQ8xLODxyhC079/JHX3uAL9zxOBmcv2C87550GrGqtpyq84m06cLxi/BvvnIvK5YuYtnCIV5x/tIJD/lJkiQd9+TugwzOgzsf28VD23fzzP4RFgzMY/H8eYwcOcLgshVrx/vuScPWySTZUFXrJ/v9lnbsOcgnv/93HD0KR4HFg3s4/+yFLJk/wHXPHOB9w6sMXJIkacIG5sHND+/gqb2H2bX3MM8fgQPPH+XZg7B9z9NkcGjxeN+ddNgCPjuF7zb1/NFi5OjPtvePHGXzzv1s2fkoq85axIFDz/OJd/7D7gqUJEmzSoDd+0fYe2iEw0egRh079nneuLdmTXrph6raNNnvduUo8ONnD3Ddxie47eEdXZcjSZJmiZGjMH8g7D4w8nNB66cy/ndPtfTD12DscwJU1TsmWOOMcnCkuPqLm/jnb3oJV77mAqcUJUnSuG57eAdfu3crP3764KS+f6ppxP/a+/1uYCXwP3vb7wcem1SPM8T+w0f5my1P8+Rzh1j/ppcYuCRJ0i+47eEd/PtvPMRPnj04/ujTKZzqacS/Bkjyx1X1G6MOfS3JrZPsc0Yojg0H7tp3mNsf2WXYkiRJP3V8AdPPfH8zO/Ye/Ll7wU/XRG+QX57kJceXeUhyIbB88t3ODDueO8TKeeHJ3ZMbFpQkSXPP8QVMjx45OuWgBRO/Qf5fArckuSXJLcDNwO9PrWtIcnmSh5NsTvKRMY4nySd7x3+Y5HVT7XO0J/ccZMHQACuWLeznaSVJ0ix2wz0/YcvOvdz8o51TDlowwZGtqvpmknXAK3q7HqqqQ1PpOMkA8GngMmArcFeSG6vqgVHN3gas6/38I+Azvd99sefACOcsmc8bX3pOv04pSZJmsUd37uW2Hz3FuUsXsO/g4b6c83TW2VoHvBxYCFyUhKr6/BT6vhjYPGpq8svAlcDosHUl8PmqKuCOJGclOa+qtk+h3586Crxu1TLv15IkScCx1/Kcu3QB8wh7Dh3pyzknNI2Y5OPAn/V+3gz8Z2Cqyz5cADwxantrb9/ptpmSL/v+REmS1PPk7oO86rwz2X/4CIeen8awBbwHeAvw91X1QeAiYNwXLk7QWMt/nfhU5UTaHGuYrE+yMcnGI/t3T7iIv999kBvu+cmE20uSpLlrxbKFLBgcZHjtWSddqPR0TDRsHaiqo8BIkqXADuAlU+x7K7B61PYqYNsk2gBQVRuqariqhgcWLzuNMop7nnj2NNpLkqS56o0vPYdnDzzP/IEBFg5O+kU7P2eiZ9mY5CzgfwCbgLuBO6fY913AuiQXJpkPXAXceEKbG4Hf6j2V+Hpgd7/u1zpu36ER+hZdJUnSrHbh8jN43/AqliwcZMmCqbxC+mdOeZYkAf5DVT0LXJPkm8DSqvrhVDquqpEkHwK+BQwA11bV/Umu7h2/BrgJuALYDOwHPjiVPk80wLGXVL929emMhEmSpLnswuVncOHyM3jquYP89+9tnvL5Thm2qqqS3AD8am/7sSn3+rNz38SxQDV63zWj+wZ+t1/9nWhoAOYPDnDla/p6z70kSZoDhteczfx5cHiaFjW9I8mvTa2rmWUAmD80yFtesdylHyRJ0s95dOde/u+Wpzl78fwpn2uik5FvBv5FkseBfRy7yamq6lemXEFHzj97IRetOosPv+VlXZciSZJmmNsf2cVZi4ZYtmQ+T+07zJHJvoWaiYett02+i5ln5dIFXLT6bP71ZS9zVEuSJP2CJ3cfZOWyhZy9eIgFg2H/85NPWxN9Xc/jk+5hhjl70SBX/dpqRo5i0JIkSWNasWwhew+OsGLpIp466xBbnznAoZHJBa6T3rOV5O5TnWAibWaKMxfM47UvPpuFQ4O+fFqSJI3r+HpbK89cwNmL57No/gDz58HQeKtFnSSHnWpk65eTnGyJhwAzft2EhUPzWDw0jyULhvgHy8/g2QPP85uvWtF1WZIkaYY6vt7W7Y/sYt/hI8wLPLj9OQYH5jEycoR9h49y/CHFwXlQR0cOjXeuU4WtV0ygnv68OKiPBhLOWTLEgsF5LJo/wL5DR1g0NMDyMxew+pwlvPGl5ziFKEmSTur4elvH3fbwDj5726PsO/Q8u/YeZt/hEapg0fx5/Jgad2zrpGFrtt6rdcbCQc5aPJ+lCwdZdfZizlu6kHkD83jf8CpDliRJmpQ3vfxFrPqlxdz+yC5+tH0PTzyzn32HjrB4wSA/OLj36fG+lxo/iM1aw8PD9b++cQu3P7KLJ3cfZMWyhY5mSZKkZpJsqqrhsY7156U/M9CJQ3+SJEld6M/rrCVJkjSmk45sJXmOsR9mPL6C/NImVUmSJM0Rp7pB/szpKkSSJGkuchpRkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1JBhS5IkqSHDliRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaGuyi0yS/BFwHrAUeA/5JVT0zRrvHgOeAI8BIVQ1PX5WSJElT19XI1keA71XVOuB7ve3xvLmqXmPQkiRJs1FXYetK4HO9z58D3tlRHZIkSU11FbZWVNV2gN7vF43TroBvJ9mUZP20VSdJktQnze7ZSvJdYOUYhz56Gqf59araluRFwHeSPFRVt47T33pgPcCaNWtOu15JkqQWmoWtqrp0vGNJnkxyXlVtT3IesGOcc2zr/d6R5HrgYmDMsFVVG4ANAMPDwzXV+iVJkvqhq2nEG4EP9D5/APjqiQ2SLEly5vHPwFuB+6atQkmSpD7oKmz9R+CyJH8HXNbbJsn5SW7qtVkB3J7kXuBO4OtV9c1OqpUkSZqkTtbZqqpdwFvG2L8NuKL3eQtw0TSXJkmS1FeuIC9JktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1JBhS5IkqSHDliRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ52ErSTvTXJ/kqNJhk/S7vIkDyfZnOQj01mjJElSP3Q1snUf8G7g1vEaJBkAPg28DXgl8P4kr5ye8iRJkvpjsItOq+pBgCQna3YxsLmqtvTafhm4EnigeYGSJEl9MpPv2boAeGLU9tbevjElWZ9kY5KNO3fubF6cJEnSRDQb2UryXWDlGIc+WlVfncgpxthX4zWuqg3ABoDh4eFx20mSJE2nZmGrqi6d4im2AqtHba8Ctk3xnJIkSdNqJk8j3gWsS3JhkvnAVcCNHdckSZJ0Wrpa+uFdSbYCbwC+nuRbvf3nJ7kJoKpGgA8B3wIeBP6qqu7vol5JkqTJ6uppxOuB68fYvw24YtT2TcBN01iaJElSX83kaURJkqRZz7AlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1JBhS5IkqSHDliRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1FAnYSvJe5Pcn+RokuGTtHssyd8muSfJxumsUZIkqR8GO+r3PuDdwGcn0PbNVfVU43okSZKa6CRsVdWDAEm66F6SJGnazPR7tgr4dpJNSdZ3XYwkSdLpajayleS7wMoxDn20qr46wdP8elVtS/Ii4DtJHqqqW8fpbz2wHmDNmjWTqlmSJKnfmoWtqrq0D+fY1vu9I8n1wMXAmGGrqjYAGwCGh4drqn1LkiT1w4ydRkyyJMmZxz8Db+XYjfWSJEmzRldLP7wryVbgDcDXk3yrt//8JDf1mq0Abk9yL3An8PWq+mYX9UqSJE1WV08jXg9cP8b+bcAVvc9bgIumuTRJkqS+mrHTiJIkSXNBqubeveRJdgKPn8ZXzgVcOLX/vK7teG3b8Lq247Vtw+vazule2xdX1fKxDszJsHW6kmysqnFfG6TJ8bq247Vtw+vajte2Da9rO/28tk4jSpIkNWTYkiRJasiwdcyGrguYo7yu7Xht2/C6tuO1bcPr2k7frq33bEmSJDXkyJYkSVJDhq2eJH+c5IdJ7kny7STnd13TXJDkvyR5qHdtr09yVtc1zRVJ3pvk/iRHk/g00hQluTzJw0k2J/lI1/XMFUmuTbIjia9b66Mkq5PcnOTB3t8Dv9d1TXNBkoVJ7kxyb++6fqIv53Ua8ZgkS6tqT+/zh4FXVtXVHZc16yV5K/D9qhpJ8p8AquoPOy5rTkjyy8BR4LPAH1TVxo5LmrWSDAA/Ai4DtgJ3Ae+vqgc6LWwOSPIbwF7g81X16q7rmSuSnAecV1V3994jvAl4p//NTk2SAEuqam+SIeB24Peq6o6pnNeRrZ7jQatnCWAK7YOq+nZVjfQ27wBWdVnPXFJVD1bVw13XMUdcDGyuqi1VdRj4MnBlxzXNCVV1K/B013XMNVW1varu7n1+DngQuKDbqma/OmZvb3Oo9zPlPGDYGiXJnyR5AvinwMe6rmcO+m3gG10XIY3hAuCJUdtb8R8uzRJJ1gKvBX7QbSVzQ5KBJPcAO4DvVNWUr+sLKmwl+W6S+8b4uRKgqj5aVauBLwIf6rba2eNU17XX5qPACMeurSZoItdWfZEx9jm6rRkvyRnAV4DfP2GGRpNUVUeq6jUcm4m5OMmUp78Hp17W7FFVl06w6V8CXwc+3rCcOeNU1zXJB4C3A28pbxI8Lafx36ymZiuwetT2KmBbR7VIE9K7p+grwBer6n93Xc9cU1XPJrkFuByY0gMeL6iRrZNJsm7U5juAh7qqZS5Jcjnwh8A7qmp/1/VI47gLWJfkwiTzgauAGzuuSRpX70buPwcerKo/7bqeuSLJ8uNPzSdZBFxKH/KATyP2JPkK8HKOPd31OHB1Vf2k26pmvySbgQXArt6uO3zKsz+SvAv4M2A58CxwT1X9ZrdVzV5JrgD+GzAAXFtVf9JxSXNCki8BlwDnAk8CH6+qP++0qDkgyRuB24C/5di/WwD/tqpu6q6q2S/JrwCf49jfA/OAv6qqfzfl8xq2JEmS2nEaUZIkqSHDliRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JI06yU5kuSeJPcnuTfJv0ryC3+/Jbkkye4kU1qLKMk/S/KpMfa/KckDSaa02rSkucWwJWkuOFBVr6mqVwGXAVcw/uu2bquqK07cmWRgqkVU1W29viXppwxbkuaUqtoBrAc+1Hulybh6I103J/lLjq3ETZIbkmzqjZKtH9X2g0l+lOSvgV9v+WeQNLe8oF5ELemFoaq29KYRX8SxV8SczMXAq6vq0d72b1fV0733ot3Ve5XXfOATwK8Cu4Gbgf/XpnpJc41hS9JcddJRrVHuHBW0AD7ce+8kwGpgHbASuKWqdgIkuQ54Wd8qlTSnGbYkzTlJXgIcAXZMoPm+Ud+7BLgUeENV7U9yC7Cwd9gXyUqaFO/ZkjSnJFkOXAN8qqpONyAtA57pBa1XAK/v7f8BcEmSc5IMAe/tX8WS5jpHtiTNBYuS3AMMASPAF4A/ncR5vglcneSHwMPAHQBVtT3JHwF/A2wH7gYGAJK8Axiuqo9N9Q8haW7K6f+PnyTNTr1pwj+oqrc37GMt8H+q6tWt+pA0uziNKOmF5DDw6qkuajqeJG8CvgY81eL8kmYnR7YkSZIacmRLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGvr/3bJOeq1kQ4MAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 720x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"SF = sample_Fisher(1000, mu=(-1, 0, 0), kappa=kappa);\n",
"fig, ax = plt.subplots(1, 1, figsize=(10, 5))\n",
......@@ -199,7 +212,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
......@@ -218,7 +231,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
......@@ -266,7 +279,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
......@@ -302,7 +315,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
......@@ -329,7 +342,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
......@@ -345,7 +358,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
......@@ -361,7 +374,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
......@@ -416,7 +429,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
......@@ -476,7 +489,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.9"
"version": "3.7.3"
}
},
"nbformat": 4,
......
%% Cell type:code id: tags:
``` python
# Imports
import sys
import os
# relative import
sys.path.append(os.path.abspath('') + '/../')
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.stats import uniform, gamma
import pyfield
from corbass.utils import load, nez2dif
```
%% Cell type:code id: tags:
``` python
# seed for reproducability
seed = 161
np.random.seed(seed)
```
%% Cell type:markdown id: tags:
# Generate synthetic data for tests
This notebook serves the purpose of generating synthetic data for testing the `CORBASS` algorithm. We therefore generate data from a given set of Gauss coefficients and add synthetic errors from the Fisher-von Mises and the gamma distribution.
%% Cell type:code id: tags:
``` python
# some basic parameters
# the name of the output file
out = '../dat/synth_data_clean_complete.csv'
# switch for using locations and incompleteness structure from the example data file
real_locs = False
# the number of records to be generated, is only used if real_locs is False
n_points = 412
# the fraction of incomplete records, works as a switch if real_locs is False
r_inc = 0.
# switch for corrupting the data by noise
noise = False
# the error levels to be stored
ddec = 4.5
dinc = 4.5
dint = 8250
# the average concentration parameter from GEOMAGIA for the interval [750, today] is 650
kappa = 650
header = f"# This file was produced using the notebook Gen_Data.ipynb with the following parameters:\n" \
+ f"# real_locs={real_locs}, n_points={n_points:d}, r_inc={r_inc:.2f}, noise={noise}, ddec={ddec:.2f}, " \
+ f"dinc={dinc:.2f}, dint={dint:d}, kappa={kappa:.1f}, seed={seed:d}\n"
```
%% Cell type:markdown id: tags:
## Sampling from the Fisher-von Mises distribution
%% Cell type:code id: tags:
``` python
# The sampling process involves some rotations, thus we first define some convenience functions
def angles(vec):
return np.arctan2(vec[1], vec[0]), \
np.pi/2 - np.arctan2(vec[2], np.sqrt(vec[0]**2 + vec[1]**2))
def rot_z(ang):
return np.array([[np.cos(ang), np.sin(ang), 0],
[-np.sin(ang), np.cos(ang), 0],
[0, 0, 1]])
def rot_y(ang):
return np.array([[np.cos(ang), 0, np.sin(ang)],
[0, 1, 0],
[-np.sin(ang), 0, np.cos(ang)]])
def rotator(vec):
vec = np.asarray(vec)
p, t = angles(vec)
return np.dot(rot_y(t).T, rot_z(p))
```
%% Cell type:code id: tags:
``` python
# This cell contains the sampling procedure
def sample_Fisher(n, mu=(0, 0, 1), kappa=20):
""" Generate samples from the Fisher distribution
Parameters:
-----------
n : int
The number of samples to be generated
mu : array-like of length 3, optional
A vector pointing towards the center of the distribution. Its length is ignored.
kappa : float, optional
The concentration parameter.
Returns:
--------
numpy array of shape (3, n) containing the sampled vectors
Reference:
----------
[1]: W. Jakob, "Numerically stable sampling of the von Mises
Fisher distribution on S^2 (and other tricks)",
http://www.mitsuba-renderer.org/~wenzel/files/vmf.pdf,
2015
"""
if kappa <= 0:
raise ValueError(f"The concentration parameter has to be positive, but kappa={kappa} was given.\n"
f"For kappa=0 use a uniform sampler on the sphere.")
trafo_mat = rotator(mu)
# sample from the uniform circle, V in [1]
angles = uniform.rvs(scale=2*np.pi, size=n)
vs = np.array([np.cos(angles),
np.sin(angles)])
# sample W in [1] via inverse cdf sampling
def inv_cdf(x):
return 1 + np.log(x + (1-x)*np.exp(-2*kappa))/kappa
unis = uniform.rvs(size=n)
ws = inv_cdf(unis)
ret = np.sqrt(1-ws**2)*vs
res = np.einsum("i...,ij->j...", np.array([ret[0], ret[1], ws]), trafo_mat)
if n == 1:
return res.flatten()
else:
return res
```
%% Cell type:markdown id: tags:
We draw some samples and plot the result, to eye-check whether the procedure works:
%% Cell type:code id: tags:
``` python
SF = sample_Fisher(1000, mu=(-1, 0, 0), kappa=kappa);
fig, ax = plt.subplots(1, 1, figsize=(10, 5))
D = np.arctan2(SF[1], SF[0])
I = np.arctan2(SF[2], np.sqrt(SF[0]**2 + SF[1]**2))
ax.set_aspect(1)
ax.set_xlim((-np.pi, np.pi));
ax.set_xlabel('D [rad.]')
ax.set_ylim((-np.pi/2., np.pi/2.));
ax.set_ylabel('I [rad.]')
ax.scatter(D, I, alpha=0.4);
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
## Generate synthetic data
We first need a set of coefficients for the field. At the time of writing this, IGRF-13 [1] had just been released and we directly download the reported coefficients as a reference model.
%% Cell type:code id: tags:
``` python
IGRF = pd.read_csv('https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf13coeffs.txt', header=0, delim_whitespace=True, skiprows=3)
coeffs = IGRF[['2020.0']].to_numpy().flatten()
# retrieve the maximal degree using pyfield and the index of the last entry in coeffs
l_max = pyfield.i2lm_l(len(coeffs)-1)
```
%% Cell type:markdown id: tags:
Then we generate random points as points of observation:
%% Cell type:code id: tags:
``` python
def ran_sph(npoints, r=1., ndim=3, nocluster=True):
""" Sample random points on the ndim-sphere.
Parameters:
-----------
npoints : int
The number of points to sample
r : float, optional
The radius of the sphere to be sampled on
ndim : int, optional
The dimension of the sphere
nocluster : bool, optional
Whether toexclude points sampled outside of the ball, i.e. sample
uniformly distributed
Returns:
--------
Array of shape (ndim, npoints) including the sampled points.
References:
-----------
[1] http://mathworld.wolfram.com/SpherePointPicking.html
[2] http://mathworld.wolfram.com/HyperspherePointPicking.html
"""
vec = np.random.randn(ndim, npoints)
if nocluster:
n = np.linalg.norm(vec, axis=0)
vec = vec[:, n <= 1.]
while vec.shape[1] < npoints:
ad = np.random.randn(ndim, npoints)
n = np.linalg.norm(ad, axis=0)
ad = ad[:, n <= 1.]
vec = np.concatenate((vec, ad), axis=1)
if vec.shape[1] > npoints:
vec = vec[:, 0:npoints]
vec /= np.linalg.norm(vec, axis=0)
if npoints == 1:
vec = vec.flatten()
return vec*r
```
%% Cell type:code id: tags:
``` python
if real_locs:
# load reference data
pars = load('../examples/Example_Parfile.py')
# ungroup the data
ref_data = pars.data.filter(lambda x: True)
ref_data.reset_index(inplace=True)
n_points = len(ref_data)
x_obs = np.zeros((3, n_points), order='F')
x_obs[0] = ref_data['co-lat']
x_obs[1] = ref_data['lon']
x_obs[2] = ref_data['rad']
else:
x_obs = ran_sph(n_points, r=pyfield.REARTH)
# transform x_obs from cartesian coordinates to co-lat, lon, rad
pyfield.mapLoc(fromSys=pyfield.SYS_GEO,
fromForm=pyfield.COOR_CAR,
toSys=pyfield.SYS_GEO,
toForm=pyfield.COOR_CLR,
t=0,
x=np.asfortranarray(x_obs))
```
%% Cell type:markdown id: tags:
Next we use the `pyfield` library to get the basis functions and generate a field from the coefficients.
%% Cell type:code id: tags:
``` python
dspharm = np.empty((len(coeffs), 3*n_points), order='F')
pyfield.dspharm(src=pyfield.SOURCE_INTERNAL,
gSys=pyfield.SYS_GEO,
atSys=pyfield.SYS_GEO,
atForm=pyfield.COOR_CLR,
bSys=pyfield.SYS_GEO,
bForm=pyfield.FIELD_NED,
lmax=l_max,
R=pyfield.REARTH,
t=0.,
at=x_obs,
B=dspharm)
```
%% Cell type:markdown id: tags:
The fields $N, E, Z$ components are then easily calculated using a dot product with the coefficients.
%% Cell type:code id: tags:
``` python
field_obs = np.dot(coeffs, dspharm).reshape(-1, 3).T
```
%% Cell type:markdown id: tags:
Using a utility function from `CORBASS`, we can transform these to $D,I,F$.
%% Cell type:code id: tags:
``` python
decs, incs, ints = nez2dif(*field_obs)
```
%% Cell type:markdown id: tags:
Now we set up error levels for the components and add some synthetic noise.
%% Cell type:code id: tags:
``` python
if noise:
mus = np.array([np.cos(np.deg2rad(incs))*np.cos(np.deg2rad(decs)),
np.cos(np.deg2rad(incs))*np.sin(np.deg2rad(decs)),
np.sin(np.deg2rad(incs))])
# the angular components retrieve an error from the Fisher-vonMises distribution
# the intensities retrieve an error from the gamma distribution
for it, mu, d, i, f in zip(np.arange(n_points), mus.T, decs, incs, ints):
samp = sample_Fisher(1, mu=mu, kappa=kappa)
decs[it] = np.rad2deg(np.arctan2(samp[1], samp[0]))
incs[it] = np.rad2deg(np.arctan2(samp[2], np.sqrt(samp[0]**2 + samp[1]**2)))
b = ints[it]/dint**2
a = ints[it] * b
ints[it] = gamma.rvs(a=a, scale=1./b)
# mimic some incompleteness in the data
if r_inc != 0.:
if real_locs:
# incompleteness of reference data
decs[ref_data.query('not (D==D)').index] = np.nan
incs[ref_data.query('not (I==I)').index] = np.nan
ints[ref_data.query('not (F==F)').index] = np.nan
else:
if r_inc != 0.:
# generate indices of missing points
ind = np.arange(n_points)
np.random.shuffle(ind)
mnum = np.int(r_inc*n_points)
mind = ind[0:mnum]
# generate a boolean array of triples of True and False for missing
# values at each point
mdist = np.zeros((3, n_points), dtype=bool)
bools = [True, True, False, False]
for j in mind:
np.random.shuffle(bools)
mdist[:, j] = bools[0:3]
# set the missing values to nan
decs[mdist[0]] = np.nan
incs[mdist[1]] = np.nan
ints[mdist[2]] = np.nan
```
%% Cell type:markdown id: tags:
Finally we build a pandas `DataFrame` from the synthetic data and store it.
%% Cell type:code id: tags:
``` python
synth_data = pd.DataFrame({'co-lat': x_obs[0],
'lat': 90-x_obs[0],
'lon': x_obs[1],
'rad': x_obs[2],
't': 2020,
'dt': 0,
'D': decs,
'I': incs,
'F': ints,
'dD': ddec,
'dI': dinc,
'dF': dint})
out_frame = synth_data.to_csv()
outfile = open(out, "w")
outfile.write(header + out_frame)
outfile.close()
```
%% Cell type:markdown id: tags:
## References
[1] P. Alken, E. Thebault, C. Beggan, H. Amit, J. Aubert, J. Baerenzung, T.N. Bondar,
W. Brown, S. Cali, A. Chambodut, A. Chulliat, G. Cox, C. C. Finlay, A. Fournier,
N. Gillet, A. Grayver, M. Hammer, M. Holschneider, L. Huder, G. Hulot, T. Jager,
C. Kloss, M. Korte, W. Kuang, A. Kuvshinov, B. Langlais, J.-M. Leger, V. Lesur,
P. W. Livermore, F. J. Lowes, S. Macmillan, W. Magnes, M. Mandea, S. Marsal,
J. Matzka, M. C. Metman, T. Minami, A. Morschhauser, J. E. Mound, M. Nair,
S. Nakano, N. Olsen, F. J. Pavon-Carrasco, V. G. Petrov, G. Ropp, M. Rother,
T. J. Sabaka, S. Sanchez, D. Saturnino, N. Schnepf, X. Shen, C. Stolle,
A. Tangborn, L. Tner-Clausen, H. Toh, J. M. Torta, J. Varner, P. Vigneron,
F. Vervelidou, I. Wardinski, J. Wicht, A. Woods, Y. Yang, Z. Zeren and B. Zhou,
"International Geomagnetic Reference Field: the thirteenth generation",
submitted to Earth, Planets and Space, see also:
https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
......
This diff is collapsed.
This diff is collapsed.
Supports Markdown
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