Gibbs_Sampling_Jupiter_testing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import all nessesary libarys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "#import cupy as np\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import binom\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Definitions of all interchangeable variables of the script"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Definitions\n",
    "\n",
    "# Matrix size\n",
    "\n",
    "matrix_width = 50\n",
    "matrix_height = 50\n",
    "\n",
    "# Grayscale value range\n",
    "\n",
    "gray_value_min = 0\n",
    "gray_value_max = 2\n",
    "\n",
    "# values --> used for würfeln.. has to be merged with grayscale_value_range\n",
    "\n",
    "#values = [-1,1]\n",
    "\n",
    "# dampening factor\n",
    "\n",
    "beta = 2\n",
    "\n",
    "#itterations\n",
    "\n",
    "itteration = 10\n",
    "\n",
    "# definition Mask\n",
    "\n",
    "mask = np.array([\n",
    "        [-0.81,0,0,0,-0.81],\n",
    "        [0,0.91,0.4,0.91,0],\n",
    "        [0,0.4,0,0.4,0],\n",
    "        [0,0.91,0.4,0.91,0],\n",
    "        [-0.81,0,0,0,-0.81]\n",
    "       ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "applys the mask to the selected data point of the matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "def apply_mask(height,width, matrix, mask):\n",
    "    # create a copy of the matrix in the area of the mask\n",
    "    matrix_slice = matrix.take(range((height - (mask.shape[0]-1)//2),(height + (mask.shape[0]+1)//2)),mode='wrap', axis=0).take(range((width - (mask.shape[1]-1)//2),(width + (mask.shape[1]+1)//2)),mode='wrap',axis=1)\n",
    "    #print (matrix_slice)\n",
    "    # calculate the probability of a data point to flip\n",
    "    xk = np.arange(gray_value_max)\n",
    "    print(xk)\n",
    "    pk = [0,0]\n",
    "    #print (matrix[height][width])\n",
    "    #for i in xk:\n",
    "    pk[matrix[height][width]]= 1 / (1 + np.exp(-beta * np.sum(mask * np.absolute(matrix[height][width] - matrix_slice))))\n",
    "    pk = np.where(pk == 0, 1-pk[matrix[height][width]], pk)\n",
    "    print(pk)\n",
    "    #use np.random.choice instead\n",
    "    np.random.choice(a= xk, p= pk)\n",
    "    #if binom.rvs(1,p=probability): \n",
    "     #   matrix[height][width] = -matrix[height][width]\n",
    "    return matrix\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "generate Matrix of white noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 33 ms, sys: 833 µs, total: 33.8 ms\n",
      "Wall time: 28.6 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 1, ..., 1, 0, 1],\n",
       "       [0, 1, 1, ..., 1, 1, 0],\n",
       "       [1, 0, 0, ..., 1, 1, 1],\n",
       "       ...,\n",
       "       [1, 1, 0, ..., 0, 0, 1],\n",
       "       [1, 1, 1, ..., 0, 0, 1],\n",
       "       [1, 0, 0, ..., 0, 0, 0]])"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD6CAYAAABnLjEDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAASA0lEQVR4nO3dX4hc130H8O+3slQH0iA7XYTQmtqlpsEPjc0MIcZ9CEpNVcfEegglpgQVBHppwCEtiZJCSaBQ5yVOHtoEEZvsQ4id2AELQymqKgiBInvWdlzbIpESHCIjW2sckfglrZJfHuau2R3f1T1z5vyb/r4fWLQzO/ee387OT3fOmXN+h2YGEfn/7/dqByAiZSjZRZxQsos4oWQXcULJLuKEkl3EiYWSneQhkj8ieYHk8VRBiUh6jP2cneQuAD8GcDeAiwCeAXC/mb18jWO2NTYajQbbWV9f33Y75JhcZmPpMxtfyDExQtoZeq5ifp8YIbHFPCbX6yfF3znX63ToeXrllVfwxhtvsO/YRZL9TgBfMLO/7G5/DgDM7F+uccy2xkLaJrfHXXMS0GwsfWbjCzkmRkg7Q89VzO8TIyS2mMfkev2k+Dvnep0OPU/j8RiTyaT3F1jkbfwBAD/fcvtid5+INOi63A2QPAbgWO52ROTaFkn2VwHctOX2anffNmZ2AsAJ4J1v40Okeqs29NYs5m1XSDu53tbnentaK5aY5zJEirfTMeeI7boMHRPz2t60yNv4ZwDcSvIWknsAfBzAyQXOJyIZRV/ZzewqyU8C+A8AuwA8YmYvJYtMRJKKHo2PaixiND7wvNtu53obn+KtZ6nR+T4xbad4XmJH/Wt9ElNrBD/kmJC2zSz5aLyILJHso/FbjUYjTCaTt2/n+t8w1WfOKT6nzvWZba6r9FA7fecYOm9s/CmuhLneSeWy6CDkeDze8XG6sos4oWQXcULJLuJE1dH4ECnmgO8Qy+A5Ykbwc4zkxrYTM0pea/ygZixDcr02YtoOHLfQaLyIZ0p2ESeU7CJOFO2zj8djG/qcvdQsoxApPmePkaI/G9NOiFJr4EPbmrftXK+5XLP9Ip8D9dlFPFOyizihZBdxQsku4kTRhTCzSlYeaWWBRK4iibkG20LaGTpPqolQMUoNXIZI8Zobik0LYUREyS7ihZJdxImqffZSRSb6tDyhpOZGGENinv/YhTylnssSfenQY2JiKVFdVkSWiJJdxAklu4gTSnYRJ6oO0KWSa+AsRT30WoM/fW3XGsgMGXwrtYtrn5YGhYd+55BdXHeiK7uIE0p2ESeU7CJONLfXW46tlUPkqlobc85UVU9KVTNNoaVJTqUq92aMRZVqRDxTsos4oWQXcaL5z9lz7fUdIkd/NVchh5bGHFK00/eYUv36XMVCZpUu6qEru4gTSnYRJ5TsIk4MJjvJR0heJvnilvtuJHmK5Pnu3xvyhikiiwq5sn8TwKGZ+44DOG1mtwI43d0eNBqNYGZvf4Ugue2rz9Zz7nTekMfUOmb2d5w9x+zPYwcOUzyXIbGkiDU2/lmzv0/f17y/T86JRjljGUx2M/s+gDdn7r4PwFr3/RqAw8EtikgVsX32fWZ2qfv+NQD7dnogyWMkJyQnGxsbkc2JyKIWHqCz6Xu9Hd+vmtkJMxub2XhlZWXR5kQkUuykmtdJ7jezSyT3A7icMqh55doud0ipyRep2sm1HXaKdmLk2p0m1wSsXAuTchevOAngSPf9EQBPRp5HRAoJ+ejt2wD+G8CfkrxI8iiABwHcTfI8gL/obotIwwbfxpvZ/Tv86MOJYxGRjIouhFlfXy+y0CKmL1qqEELrscTItWAlRSHImOIPIXGUHIeY5+fXoumyIk4o2UWcULKLOKFkF3Gi6ADdaDTCZDJZ6By1Bkb65Ko0kqsCba4KqLl+51KVbWeleJ5iJ+/kHCzUlV3ECSW7iBNKdhEnqlaXLdlnzNH/a2lHmz4p+n8t9fNT9HlbWgiT6vkPpSu7iBNKdhEnlOwiTlTdxbVPrd08c32WGnNMzGfoqcYkYuKf9xyh5xk6b6r+d4p+cc3fuecY7eIq4pmSXcQJJbuIE0p2ESeqLoRJVZEzxeBIya2Ih46pWZ01xTG5Ku3EnKPUwFnMMTkGVXNUlxWRJaNkF3FCyS7iRNWFMKkKLpSarBBz3lxKVaQNkavgQo4JVn3njHkuU0x8Kl2wQ1d2ESeU7CJOKNlFnKjaZw9R6zPQ2LZLxZtCrkIOpaRa5BKzC0uthTzaEUZEBinZRZxQsos4oWQXcaK5SjWzalWO7VNqIk6pSiklB57mPUeIUgN0LQ1SBlakVaUaEc+U7CJODCY7yZtIniH5MsmXSD7Q3X8jyVMkz3f/3pA/XBGJNdhnJ7kfwH4ze5bkHwBYB3AYwN8CeNPMHiR5HMANZvbZa51rPB7bUPGKXP2/UhVbc0yKqDn5pWbl1ZjnMscCldar2Cbrs5vZJTN7tvv+VwDOATgA4D4Aa93D1jD9D0BEGjVXn53kzQDuAHAWwD4zu9T96DUA+9KGJiIpBSc7yXcDeALAp8zsl1t/ZtP3Fr3vl0geIzkhOdnY2FgoWBGJF5TsJHdjmujfMrPvdXe/3vXnN/v1l/uONbMTZjY2s/HKykqKmEUkwuCqN057/w8DOGdmX97yo5MAjgB4sPv3yRQBpah6kmpb4aFztFQdptT2vzETcUL+HqnOG/OYmFVvMY8JOSbF77yTkCWudwH4BID/Ifl8d9/nMU3y75A8CuBnAP46OgoRyW4w2c3sBwB2+u/kw2nDEZFcNINOxImilWrW19fn7nPkmuAQ0ufNsZVyrsUnIWJiifmda+4ak6PSbUtbZg/RjjAiomQX8ULJLuJEc9Vlcy3wyLG7aqpjYrS0U+qQ2L9hrvGCec8RO08gRs7dhnVlF3FCyS7ihJJdxAklu4gTzVWXrTkgNK9Sky1yTVLJVWm11GSjks/LvOcIkfH5V3VZEc+U7CJOKNlFnCg6qWY0GmFrddkYuQoh1No2ua/tVAt5UhwzK9fOLSFthRQpmfecqaRa8DR0nqH4tRBGRJTsIl4o2UWcqFq8ItfnsTGFKFJ95pmrqESKdmJiCWm31K43JccLhrRcvHMnurKLOKFkF3FCyS7ihJJdxInmJ9XkqgCSayJLK+30STGJJtXWyrNKDaTFHJMrtlKDept0ZRdxQsku4oSSXcSJ5qrL5th5NERMnzdmgkPNHWFC2m25v1qqcm+IFK/TVK+FULqyizihZBdxQsku4kRzffZanwXHnLfm58e5CjbG9EVrfebfJ1dfOoUU4yGLvJ50ZRdxQsku4oSSXcSJwWQneT3Jp0n+kORLJL/Y3X8LybMkL5B8jOSe/OGKSKyQAbpfAzhoZm+R3A3gByT/HcCnATxkZo+S/DqAowC+tmhAtSqllNpyN6TtmCqqpSrtlJ4Icq12YioSpWpnSOzrKeekpsEru0291d3c3X0ZgIMAHu/uXwNwOFlUIpJcUJ+d5C6SzwO4DOAUgJ8AuGJmV7uHXARwYIdjj5GckJxsbGwkCFlEYgQlu5n9xsxuB7AK4AMA3hfagJmdMLOxmY1XVlbiohSRhc01qcbMrpA8A+BOAHtJXtdd3VcBvJojwFmx/ddZpSY45OqD1dp5NEaqiSwxk6dyTbiKmbwTs3NRyolbIaPxKyT3dt+/C8DdAM4BOAPgY93DjgB4MqhFEaki5Mq+H8AayV2Y/ufwHTN7iuTLAB4l+c8AngPwcMY4RWRBg8luZi8AuKPn/p9i2n8XkSWgGXQiTjS//VOqAYqYFWw5KqOkGogqvWJqU8xElpYryvQ9pmYs855j9jzasllElOwiXijZRZyoWqkm1+KNVH3cmFhizjnUl6vZ542xDJVWr6XUhKVUu+toRxgR2UbJLuKEkl3Eieaqy86q+TlprUIauT6nbmkXmVpzFkKUmsOQK/6d6Mou4oSSXcQJJbuIE0p2ESeqDtDFDLa1vn1PjkUhtSaXAHETfFJs09R3X65KsTmqCcW+BhetgKOFMCKiZBfxQsku4gRLLrIgua2xXFsel9pKOdU2yUPnyFWFtGb8IXL8jq1X9030uuw9ia7sIk4o2UWcULKLONFc8YpS583V98y5C+e8sdRqp9SYSapFOSXOkfu4TfqcXUSU7CJeKNlFnFCyizhRdIBuNBphMpnMdUyqSjUxC1RyVPrMtZAnVaXemMfM227fY0oNgpWayJVr8s4idGUXcULJLuKEkl3Eieaqy5YqHFBrUk1IPz9XRd1lK75RanfYUotnZpXeBUdXdhEnlOwiTgQnO8ldJJ8j+VR3+xaSZ0leIPkYyT35whSRRc3TZ38AwDkA7+lufwnAQ2b2KMmvAzgK4GuJ43uHUosfQtuOOc/QMaUKIcQUxYj5PD9V8Ypcx8TEn+K5Lb27TtCVneQqgI8A+EZ3mwAOAni8e8gagMNBLYpIFaFv478C4DMAftvdfi+AK2Z2tbt9EcCBvgNJHiM5ITnZ2NhYJFYRWcBgspO8F8BlM1uPacDMTpjZ2MzGKysrMacQkQRC+ux3AfgoyXsAXI9pn/2rAPaSvK67uq8CeDVfmCKyqLmqy5L8EIB/MLN7SX4XwBNbBuheMLN/Gzh+7pGaUtVYcw2WpBhgSTUpJdV5h9optWtPropEIUpNLop8/SevLvtZAJ8meQHTPvzDC5xLRDKrWjc+hK7surLvRFf2/vOqbryIc83t4jqk5I6mpSalDMk1KaXm1anU7jS5lCpwkXLxj67sIk4o2UWcULKLOFG14GSpUdu+87S8qCXF7qUhx5UazY4tbJmjiERLnxSEtD3v4hntCCMiSnYRL5TsIk4o2UWcaH5STa5qJSkmOIScN5dSk1JqbjmdaqBy3lhSTGRJVeE45e+sK7uIE0p2ESeU7CJONLcjzKyau3PEHJNi8UxM27nGKYbOEXuekHOk2Hk3pu0UFXVLLkvWpBoR2UbJLuKEkl3ECSW7iBNVB+hy1UDL1VaKLZtjzhsyQJRrgDHFJKHYySOl6gbGxBZzTK4tp0Ppyi7ihJJdxAklu4gTzS2ESbH4IVVfKMVCkhyTPEIfM287IeepuRgoxeKfVO3kOCY3XdlFnFCyizihZBdxormFMLkquqb4LDvFDh4xsfVJ8Xl+zeIhs3LNYQhpJ0d/O/bvnLMfryu7iBNKdhEnlOwiTijZRZxYuoUwuRZVpBJTXSXFNs+5tl/Otfgk1/bFpSeqbKq5TZmqy4rINkp2ESeU7CJOsOQCepIbAH4G4A8BvFGs4cUsU6zAcsW7TLECyxHvH5nZSt8Piib7242SEzPbueZtQ5YpVmC54l2mWIHli3eW3saLOKFkF3GiVrKfqNRujGWKFViueJcpVmD54t2mSp9dRMrT23gRJ4omO8lDJH9E8gLJ4yXbDkHyEZKXSb645b4bSZ4ieb7794aaMW4ieRPJMyRfJvkSyQe6+1uN93qST5P8YRfvF7v7byF5tntNPEZyT+1YN5HcRfI5kk91t5uNNUSxZCe5C8C/AvgrALcBuJ/kbaXaD/RNAIdm7jsO4LSZ3QrgdHe7BVcB/L2Z3QbggwD+rns+W4331wAOmtn7AdwO4BDJDwL4EoCHzOxPAPwCwNF6Ib7DAwDObbndcqyDSl7ZPwDggpn91Mz+F8CjAO4r2P4gM/s+gDdn7r4PwFr3/RqAwyVj2omZXTKzZ7vvf4Xpi/IA2o3XzOyt7ubu7ssAHATweHd/M/GSXAXwEQDf6G4TjcYaqmSyHwDw8y23L3b3tW6fmV3qvn8NwL6awfQheTOAOwCcRcPxdm+LnwdwGcApAD8BcMXMrnYPaek18RUAnwHw2+72e9FurEE0QDcHm3500dTHFyTfDeAJAJ8ys19u/Vlr8ZrZb8zsdgCrmL7Te1/diPqRvBfAZTNbrx1LSiXXs78K4KYtt1e7+1r3Osn9ZnaJ5H5Mr0pNILkb00T/lpl9r7u72Xg3mdkVkmcA3AlgL8nruitmK6+JuwB8lOQ9AK4H8B4AX0WbsQYreWV/BsCt3YjmHgAfB3CyYPuxTgI40n1/BMCTFWN5W9eHfBjAOTP78pYftRrvCsm93ffvAnA3puMMZwB8rHtYE/Ga2efMbNXMbsb0dfpfZvY3aDDWuZhZsS8A9wD4MaZ9tX8s2XZgfN8GcAnA/2HaJzuKaV/tNIDzAP4TwI214+xi/XNM36K/AOD57uuehuP9MwDPdfG+COCfuvv/GMDTAC4A+C6A368d60zcHwLw1DLEOvSlGXQiTmiATsQJJbuIE0p2ESeU7CJOKNlFnFCyizihZBdxQsku4sTvALsHgWjhiSUrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%time\n",
    "# generate Matrix of white noise\n",
    "matrix = np.random.randint(gray_value_min,gray_value_max,size=(matrix_height,matrix_width))\n",
    "# replaces 0 in the matrix with -1s (nessesary for a few masks)\n",
    "#matrix = np.where(matrix == 0, -1, matrix)\n",
    "# Print matrix as Image\n",
    "plt.imshow(matrix,cmap='gray',vmin=gray_value_min, vmax=gray_value_max-1)\n",
    "# Return (and print) matrix\n",
    "matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix = np.load('matrix.npy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "apply mask to whole matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1]\n",
      "[0.72710822 0.        ]\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "probabilities do not sum to 1",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<timed exec>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n",
      "\u001b[0;32m/tmp/ipykernel_67179/115657037.py\u001b[0m in \u001b[0;36mapply_mask\u001b[0;34m(height, width, matrix, mask)\u001b[0m\n\u001b[1;32m     13\u001b[0m     \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m     \u001b[0;31m#use np.random.choice instead\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 15\u001b[0;31m     \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mchoice\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m=\u001b[0m \u001b[0mxk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m=\u001b[0m \u001b[0mpk\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     16\u001b[0m     \u001b[0;31m#if binom.rvs(1,p=probability):\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     17\u001b[0m      \u001b[0;31m#   matrix[height][width] = -matrix[height][width]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mmtrand.pyx\u001b[0m in \u001b[0;36mnumpy.random.mtrand.RandomState.choice\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: probabilities do not sum to 1"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# itterate thrue the whole matrix and appy mask\n",
    "counter = 0\n",
    "np.random.seed(0)\n",
    "while counter < itteration:\n",
    "    counter += 1\n",
    "    for h in range(0, matrix.shape[0]):\n",
    "        for w in range(0, matrix.shape[1]):\n",
    "            matrix = apply_mask(h, w, matrix, mask) #faster with inplace ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "plot new matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f7f721d9040>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD6CAYAAABnLjEDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAARyUlEQVR4nO3dX6il1XnH8e+vo9bQENTkMAwztloihLlolXMQxV4Ei9SaEL2QYAhlLoS5ScGQQKottA0U2tzkz0VpGaJkLkLUGqkigWInE0Ju1H2iJuqQOBElyuiMJJJ403aSpxfnNZzZ7jPvu9dea73vcf0+cJiz/7xrPfs9+5m119prrVcRgZm99/3e2AGYWR1OdrNGONnNGuFkN2uEk92sEU52s0aslOySbpb0E0knJd2dKygzy0+p37NL2gP8FLgJeBV4CvhURLxwnmOWrmx9fT0pvhw2Nzezlznk9QypN8d5WVRPX7kpx+QyX3fKuSx1TC6r1v3yyy/z5ptvatFjqyT79cA/RsRfdLfvAYiIfz7PMUtXNuakH2nhOVvJkNczpN4c52VRPX3lphyTy3zdKeey1DG5rFr3xsYGs9ls4RtolY/x+4Gfb7v9anefmU3QBaUrkHQYOFy6HjM7v1WS/TXg8m23D3T3nSMijgBHIO1j/JCPNTmeU+Ij+yK16hkiV5ei7zm56hlyTMpH7r66a3ZdSnYZVvkY/xRwlaQrJV0E3AE8micsM8stuWWPiLOS/hr4L2APcF9EPJ8tMjPLaqU+e0R8B/hOpljMrCDPoDNrRPHR+NxSB4xKDMjl+s621ABXynfmfWWMOZCWUteUB2JznYOhr8ktu1kjnOxmjXCymzVi8n32XH2ulH5xjjnSQ47J8RpTxgKGqDW2kctYaylqjWWs8vdwy27WCCe7WSOc7GaNmHyfPdf3lSnfv/Y9Z6x13EPrrrU2vdTYRonzXWo9e63v8/ti2djY2PExt+xmjXCymzXCyW7WCCe7WSOqDtCtr68zm81qVrmjHBNXau2UkkuOAa+aA5m1BkjH3GAyBy+EMbNzONnNGuFkN2vEe2JSTQ6ldkCtpdRCjBy7y6aObUzlb19qM5QSYyaeVGNmTnazVjjZzRrhZDdrxKgDdFMe8IIy8U15dd3UlHjdpc5/qQlXQ8rwpBozO4eT3awRTnazRlTts29ubp7Tvyg1sSLXJZtrXUWm1i62Q6TEkmtCz7yxruZSS65+vSfVmNk5nOxmjXCymzVi8ptX1Opb59itdYgp9Tun3tcuca5SvqcutWCo1OYhO3HLbtYIJ7tZI5zsZo3oTXZJ90k6Lem5bfddJulxSS92/15aNkwzW9WQlv0bwM1z990NHIuIq4Bj3e3JiIjenxSSen9yGFLPkHpzxDrkmPlzm1JPqb9ZX2ySeusZEktK/EP+Jsv+nTc3N3c8F73JHhHfB34xd/etwNHu96PAbX3lmNm4UvvseyPiVPf768DenZ4o6bCkmaTZmTNnEqszs1WtPEAXW59PdvyMFRFHImIjIjbW1tZWrc7MEqVOqnlD0r6IOCVpH3B6yEHzC2FSpC4kqTVBo1a586+n1mSjnfqV5zsmVYlFOaUWIuUcrykltWV/FDjU/X4IeCRPOGZWigbsmf0t4KPAh4A3gH8A/hN4EPhD4BXgkxExP4i3qKyV/8vf7S17ynTTXOWUatn7jkmVo2UvdUxfGWOKiIXB9CZ7Tk52J/synOxp3jPJPiU1N98oYcqxLVLzP5q+ekttrJFydZoFZSx8kqfLmjXCyW7WCCe7WSOc7GaNmNxONX0jo6kDIWNdYrfUjialRuNTBohKDdqVmqxTo96xJlx5d1kzc7KbtcLJbtaIUa/iukhfP3OVjSdylLNqvTXLyLEoJIfUfn6t3Xzn65lSHz0nt+xmjXCymzXCyW7WiFH77EP6cqVWgKXEkmLId91TWlSRMmZSYpOJoeUsW26t+Mc8Bztxy27WCCe7WSOc7GaNcLKbNWLX7VRTalJNrUGxIcbakWVR3TkGGHNN3qm1DdhY25wtkriYyTvVmLXMyW7WCCe7WSMmvxBm3pibV/SVmWoqE0FSj+l7zpjbLI951Z4+Oc7tMtyymzXCyW7WCCe7WSMm12efV+rqmCnfH+eod4iUjR1SFhVN7JJF77qvxMYftTbNGKLUOMtO3LKbNcLJbtYIJ7tZI5zsZo2Y/ADdvClPkoDxBpVyqXVt8lIDZaUmDs0rdd36koOmbtnNGuFkN2tEb7JLulzScUkvSHpe0l3d/ZdJelzSi92/l5YP18xSDWnZzwKfj4iDwHXAZyQdBO4GjkXEVcCx7vZ5ra+vExHn/ZnX93guknp/hhyTI94hZfTVu9MklWWPySFXPX1/j5QycpXT976OiCz1LrK9jvX19R2f15vsEXEqIn7Y/f5r4ASwH7gVONo97Shw26pBm1k5S/XZJV0BXAM8AeyNiFPdQ68De/OGZmY5DU52Se8Hvg18NiJ+tf2x2PpctvCzmaTDkmaSZmfOnFkpWDNLNyjZJV3IVqJ/MyIe7u5+Q9K+7vF9wOlFx0bEkYjYiIiNtbW1HDGbWYLeSTXaGkm4FzgREV/e9tCjwCHgX7p/HykS4bvjSTquxAq2MS/Tm+OywrlW+uVYTTfWZJgh5aasyKt1ObFlDJlBdwPwV8CPJT3T3fe3bCX5g5LuBF4BPpkchZkV15vsEfEDYKf/xv48bzhmVopn0Jk1oupCmM3NzaV3Uq21U00ppfqzOfrfY+5am6LU4pOUelLObQ6rjAW4ZTdrhJPdrBFOdrNGTH7zilx9xrG+E8/RR0/d6GFKu8emSLm66rwpj0HkmicwtBy37GaNcLKbNcLJbtYIJ7tZI6oO0K2vrzObzbKXm2OgJtcEkxxq7apaq9yUS1MtOi7H4NpuG7RcdnB2Y2Njx+e5ZTdrhJPdrBFOdrNGjLoQpuZVWPrqSllskiv+WvXkiCWljFyxlDLW+a89Ecotu1kjnOxmjXCymzVi1IUwKQvxay56mcrGiamvOUf/r8RGnanljtmXzqHU2NJQbtnNGuFkN2uEk92sEU52s0ZMfqeaebkGosacLJKj3vn4x3w9y+4YvOiYWhNMUmKZ0uIZ71RjZr2c7GaNcLKbNWLUPnupTRrG7GOVmoRSq54SV1DJ9fdIuVJLjnOZ0k9O7VuvurOwN68wMye7WSuc7GaNmPxCmFL97xxX6kyRuvlijrpSNnmcV3NjihJ96d02TpGTW3azRjjZzRrhZDdrRG+yS7pY0pOSnpX0vKQvdvdfKekJSSclPSDpovLhmlkqDdgZQ8AfRMTbki4EfgDcBXwOeDgi7pf078CzEfFvPWUVWUmSY7LFmLvbzBtzp5Q+pRaS7Lbzn2Mgs9QlmyNiYcG9LXtsebu7eWH3E8CNwEPd/UeB24YGbGb1DeqzS9oj6RngNPA48DPgrYg42z3lVWD/DsceljSTlP8ib2Y22KBkj4jfRMTVwAHgWuAjQyuIiCMRsRERO0/aNbPilppUExFvSToOXA9cIumCrnU/ALyWI6CU/l6O56T0K6d05dRaO62Wmgg15vnPcbWglGNybdiRbSGMpDVJl3S/vw+4CTgBHAdu7552CHhk6SjNrJohLfs+4KikPWz95/BgRDwm6QXgfkn/BDwN3FswTjNbUW+yR8SPgGsW3P8SW/13M9sFPIPOrBFVV72tr68zm632DVypyQpDjLUjai6lLsdcanfWEoOOOVa4QdqKwpS6hqymG1qXW3azRjjZzRrhZDdrxOR3qinVx0otp6/MkpfcXVWtyS5DHh9rUU6KUpNqcvEVYczsHE52s0Y42c0aUbXPvrm5eU7/olbf+r1orN1xS/WTS10dKEWOhVWl3tt9f2dfEcbMnOxmrXCymzXCyW7WiFEXwpQajKu1wKPUQFSuS/uW2AG11ESimpNQ+uTYHTfX4p+cg5Bu2c0a4WQ3a4ST3awRoy6EGfOSu7XKTVFrp9hcx0x5w46UCT+lxkNy7ZSceoxbdrNGONnNGuFkN2vEqH32RXJcbbVUH6vmd/x9Zaach1qvJ0Wt8ZGUDVPGHM9Zth4vhDEzJ7tZK5zsZo1wsps1YnIDdPNyDYSUuNJGrssX5xgEKzkZo4aUc1lq19pcC5Hm5Vh4tEo+uGU3a4ST3awRTnazRkyuzz6lDS2mPKlmSleUHevqsEPLOV+ZO91XQ64xn6Hcsps1wslu1ojByS5pj6SnJT3W3b5S0hOSTkp6QNJF5cI0s1Ut07LfBZzYdvtLwFci4sPAL4E7l61c0rt+UgwpIyLO+Rmi75j5x4c+p++Y+dczpIwh53JIOcvGOuS85Sq3lLFiGfJ3zGlQsks6AHwM+Hp3W8CNwEPdU44Ct2WNzMyyGtqyfxX4AvDb7vYHgbci4mx3+1Vg/6IDJR2WNJM0O3PmzCqxmtkKepNd0seB0xGxmVJBRByJiI2I2FhbW0spwswyGPI9+w3AJyTdAlwMfAD4GnCJpAu61v0A8Fq5MM1sVb0te0TcExEHIuIK4A7guxHxaeA4cHv3tEPAI8Wi7FFqEClHLCkDLimDb0MMKSNlkG/Zeoe+nhKDVymDqkPKSZHyd13lvbDK9+x/A3xO0km2+vD3rlCWmRW21HTZiPge8L3u95eAa/OHZGYleAadWSMmtxBmXsmFAWMotUhk6HGryrV4o9QxJV5zrjJzvHdXicUtu1kjnOxmjXCymzVi8n32nBvuLVNPal0lxhRW+V59u5R+cqkxkpTNN2pt2DH2wpztco5TuGU3a4ST3awRTnazRjjZzRpRdYBuc3Nz5QGUVRZALFvOWBN4hi4uKVHPWFeWGave1HL73k9TmojzDrfsZo1wsps1wslu1oiqffb19XVms9l5n9PXR6nVn526UguEdtPCoylPfpliPW7ZzRrhZDdrhJPdrBFOdrNGTG7VW47JL0Mmi5Ra0ZYy2WLKK7VylJu6UqvWisdaxn7NbtnNGuFkN2uEk92sEaP22ac+ESTHTi9D+lylFlWklFOiz1hqIlSp909K3bnObckFNm7ZzRrhZDdrhJPdrBGq+d2lpCKV5eiz5+pzTeX7/FJq7u7bJ2XOQql5DkOUep8ueM8tLNgtu1kjnOxmjXCymzXCyW7WiMkthOmTMkllTCkDdqUGC4fUU2siUY5jSl2yudSloUvtwrs9lo2NjR2Pdctu1ggnu1kjnOxmjag9qeYM8ArwIeDNahWvZjfFCrsr3t0UK+yOeP8oItYWPVA12X9XqTSLiJ1HEiZkN8UKuyve3RQr7L545/ljvFkjnOxmjRgr2Y+MVG+K3RQr7K54d1OssPviPccofXYzq88f480aUTXZJd0s6SeSTkq6u2bdQ0i6T9JpSc9tu+8ySY9LerH799IxY3yHpMslHZf0gqTnJd3V3T/VeC+W9KSkZ7t4v9jdf6WkJ7r3xAOSLho71ndI2iPpaUmPdbcnG+sQ1ZJd0h7gX4G/BA4Cn5J0sFb9A30DuHnuvruBYxFxFXCsuz0FZ4HPR8RB4DrgM935nGq8/wPcGBF/ClwN3CzpOuBLwFci4sPAL4E7xwvxXe4CTmy7PeVYe9Vs2a8FTkbESxHxv8D9wK0V6+8VEd8HfjF3963A0e73o8BtNWPaSUSciogfdr//mq035X6mG29ExNvdzQu7nwBuBB7q7p9MvJIOAB8Dvt7dFhONdaiayb4f+Pm22692903d3og41f3+OrB3zGAWkXQFcA3wBBOOt/tY/AxwGngc+BnwVkSc7Z4ypffEV4EvAL/tbn+Q6cY6iAfolhBbX11M6usLSe8Hvg18NiJ+tf2xqcUbEb+JiKuBA2x90vvIuBEtJunjwOmI2Bw7lpxqrmd/Dbh82+0D3X1T94akfRFxStI+tlqlSZB0IVuJ/s2IeLi7e7LxviMi3pJ0HLgeuETSBV2LOZX3xA3AJyTdAlwMfAD4GtOMdbCaLftTwFXdiOZFwB3AoxXrT/UocKj7/RDwyIix/E7Xh7wXOBERX9720FTjXZN0Sff7+4Cb2BpnOA7c3j1tEvFGxD0RcSAirmDrffrdiPg0E4x1KRFR7Qe4BfgpW321v6tZ98D4vgWcAv6PrT7ZnWz11Y4BLwL/DVw2dpxdrH/G1kf0HwHPdD+3TDjePwGe7uJ9Dvj77v4/Bp4ETgL/Afz+2LHOxf1R4LHdEGvfj2fQmTXCA3RmjXCymzXCyW7WCCe7WSOc7GaNcLKbNcLJbtYIJ7tZI/4fRBd+V9+QV90AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(matrix,cmap='gray',vmin=gray_value_min, vmax=gray_value_max-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "safe new matrix to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "np.save(\"matrix\", matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fragen für Mittwoch:\n",
    "Formel (was liegt mit in der summe? O.o)\n",
    "        e^-2 abhängigkeit zu anz graustufen\n",
    "        mehrere pixel zeitgleich angucken (maske gerade anzahl (8x8 oder so mit 4 mittleren pixeln)\n",
    "        exponenten ersetzen.. wie genau soll ich das umsetzen? O.o\n",
    "        und wo genau brauch ich fakultäten"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}