EntropyTrainer-checkpoint.ipynb 134 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Ufb5Oc9LHkOd"
   },
   "source": [
    "# Test NN for Entropy Closure #\n",
    "\n",
    "Test Data is generated from RTSN code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Zs9ZGENoItzq"
   },
   "source": [
    "Load all needed modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "import pandas as pd\n",
    "import random                        # for generating random numbers\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "from tensorflow import keras\n",
    "import matplotlib.pyplot as plt      # MATLAB like plotting routines\n",
    "from sklearn.preprocessing import normalize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Z9bPtTIHIxnR"
   },
   "source": [
    "Load and preprocess Training Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 78,
     "resources": {
      "http://localhost:8080/nbextensions/google.colab/files.js": {
       "data": "Ly8gQ29weXJpZ2h0IDIwMTcgR29vZ2xlIExMQwovLwovLyBMaWNlbnNlZCB1bmRlciB0aGUgQXBhY2hlIExpY2Vuc2UsIFZlcnNpb24gMi4wICh0aGUgIkxpY2Vuc2UiKTsKLy8geW91IG1heSBub3QgdXNlIHRoaXMgZmlsZSBleGNlcHQgaW4gY29tcGxpYW5jZSB3aXRoIHRoZSBMaWNlbnNlLgovLyBZb3UgbWF5IG9idGFpbiBhIGNvcHkgb2YgdGhlIExpY2Vuc2UgYXQKLy8KLy8gICAgICBodHRwOi8vd3d3LmFwYWNoZS5vcmcvbGljZW5zZXMvTElDRU5TRS0yLjAKLy8KLy8gVW5sZXNzIHJlcXVpcmVkIGJ5IGFwcGxpY2FibGUgbGF3IG9yIGFncmVlZCB0byBpbiB3cml0aW5nLCBzb2Z0d2FyZQovLyBkaXN0cmlidXRlZCB1bmRlciB0aGUgTGljZW5zZSBpcyBkaXN0cmlidXRlZCBvbiBhbiAiQVMgSVMiIEJBU0lTLAovLyBXSVRIT1VUIFdBUlJBTlRJRVMgT1IgQ09ORElUSU9OUyBPRiBBTlkgS0lORCwgZWl0aGVyIGV4cHJlc3Mgb3IgaW1wbGllZC4KLy8gU2VlIHRoZSBMaWNlbnNlIGZvciB0aGUgc3BlY2lmaWMgbGFuZ3VhZ2UgZ292ZXJuaW5nIHBlcm1pc3Npb25zIGFuZAovLyBsaW1pdGF0aW9ucyB1bmRlciB0aGUgTGljZW5zZS4KCi8qKgogKiBAZmlsZW92ZXJ2aWV3IEhlbHBlcnMgZm9yIGdvb2dsZS5jb2xhYiBQeXRob24gbW9kdWxlLgogKi8KKGZ1bmN0aW9uKHNjb3BlKSB7CmZ1bmN0aW9uIHNwYW4odGV4dCwgc3R5bGVBdHRyaWJ1dGVzID0ge30pIHsKICBjb25zdCBlbGVtZW50ID0gZG9jdW1lbnQuY3JlYXRlRWxlbWVudCgnc3BhbicpOwogIGVsZW1lbnQudGV4dENvbnRlbnQgPSB0ZXh0OwogIGZvciAoY29uc3Qga2V5IG9mIE9iamVjdC5rZXlzKHN0eWxlQXR0cmlidXRlcykpIHsKICAgIGVsZW1lbnQuc3R5bGVba2V5XSA9IHN0eWxlQXR0cmlidXRlc1trZXldOwogIH0KICByZXR1cm4gZWxlbWVudDsKfQoKLy8gTWF4IG51bWJlciBvZiBieXRlcyB3aGljaCB3aWxsIGJlIHVwbG9hZGVkIGF0IGEgdGltZS4KY29uc3QgTUFYX1BBWUxPQURfU0laRSA9IDEwMCAqIDEwMjQ7CgpmdW5jdGlvbiBfdXBsb2FkRmlsZXMoaW5wdXRJZCwgb3V0cHV0SWQpIHsKICBjb25zdCBzdGVwcyA9IHVwbG9hZEZpbGVzU3RlcChpbnB1dElkLCBvdXRwdXRJZCk7CiAgY29uc3Qgb3V0cHV0RWxlbWVudCA9IGRvY3VtZW50LmdldEVsZW1lbnRCeUlkKG91dHB1dElkKTsKICAvLyBDYWNoZSBzdGVwcyBvbiB0aGUgb3V0cHV0RWxlbWVudCB0byBtYWtlIGl0IGF2YWlsYWJsZSBmb3IgdGhlIG5leHQgY2FsbAogIC8vIHRvIHVwbG9hZEZpbGVzQ29udGludWUgZnJvbSBQeXRob24uCiAgb3V0cHV0RWxlbWVudC5zdGVwcyA9IHN0ZXBzOwoKICByZXR1cm4gX3VwbG9hZEZpbGVzQ29udGludWUob3V0cHV0SWQpOwp9CgovLyBUaGlzIGlzIHJvdWdobHkgYW4gYXN5bmMgZ2VuZXJhdG9yIChub3Qgc3VwcG9ydGVkIGluIHRoZSBicm93c2VyIHlldCksCi8vIHdoZXJlIHRoZXJlIGFyZSBtdWx0aXBsZSBhc3luY2hyb25vdXMgc3RlcHMgYW5kIHRoZSBQeXRob24gc2lkZSBpcyBnb2luZwovLyB0byBwb2xsIGZvciBjb21wbGV0aW9uIG9mIGVhY2ggc3RlcC4KLy8gVGhpcyB1c2VzIGEgUHJvbWlzZSB0byBibG9jayB0aGUgcHl0aG9uIHNpZGUgb24gY29tcGxldGlvbiBvZiBlYWNoIHN0ZXAsCi8vIHRoZW4gcGFzc2VzIHRoZSByZXN1bHQgb2YgdGhlIHByZXZpb3VzIHN0ZXAgYXMgdGhlIGlucHV0IHRvIHRoZSBuZXh0IHN0ZXAuCmZ1bmN0aW9uIF91cGxvYWRGaWxlc0NvbnRpbnVlKG91dHB1dElkKSB7CiAgY29uc3Qgb3V0cHV0RWxlbWVudCA9IGRvY3VtZW50LmdldEVsZW1lbnRCeUlkKG91dHB1dElkKTsKICBjb25zdCBzdGVwcyA9IG91dHB1dEVsZW1lbnQuc3RlcHM7CgogIGNvbnN0IG5leHQgPSBzdGVwcy5uZXh0KG91dHB1dEVsZW1lbnQubGFzdFByb21pc2VWYWx1ZSk7CiAgcmV0dXJuIFByb21pc2UucmVzb2x2ZShuZXh0LnZhbHVlLnByb21pc2UpLnRoZW4oKHZhbHVlKSA9PiB7CiAgICAvLyBDYWNoZSB0aGUgbGFzdCBwcm9taXNlIHZhbHVlIHRvIG1ha2UgaXQgYXZhaWxhYmxlIHRvIHRoZSBuZXh0CiAgICAvLyBzdGVwIG9mIHRoZSBnZW5lcmF0b3IuCiAgICBvdXRwdXRFbGVtZW50Lmxhc3RQcm9taXNlVmFsdWUgPSB2YWx1ZTsKICAgIHJldHVybiBuZXh0LnZhbHVlLnJlc3BvbnNlOwogIH0pOwp9CgovKioKICogR2VuZXJhdG9yIGZ1bmN0aW9uIHdoaWNoIGlzIGNhbGxlZCBiZXR3ZWVuIGVhY2ggYXN5bmMgc3RlcCBvZiB0aGUgdXBsb2FkCiAqIHByb2Nlc3MuCiAqIEBwYXJhbSB7c3RyaW5nfSBpbnB1dElkIEVsZW1lbnQgSUQgb2YgdGhlIGlucHV0IGZpbGUgcGlja2VyIGVsZW1lbnQuCiAqIEBwYXJhbSB7c3RyaW5nfSBvdXRwdXRJZCBFbGVtZW50IElEIG9mIHRoZSBvdXRwdXQgZGlzcGxheS4KICogQHJldHVybiB7IUl0ZXJhYmxlPCFPYmplY3Q+fSBJdGVyYWJsZSBvZiBuZXh0IHN0ZXBzLgogKi8KZnVuY3Rpb24qIHVwbG9hZEZpbGVzU3RlcChpbnB1dElkLCBvdXRwdXRJZCkgewogIGNvbnN0IGlucHV0RWxlbWVudCA9IGRvY3VtZW50LmdldEVsZW1lbnRCeUlkKGlucHV0SWQpOwogIGlucHV0RWxlbWVudC5kaXNhYmxlZCA9IGZhbHNlOwoKICBjb25zdCBvdXRwdXRFbGVtZW50ID0gZG9jdW1lbnQuZ2V0RWxlbWVudEJ5SWQob3V0cHV0SWQpOwogIG91dHB1dEVsZW1lbnQuaW5uZXJIVE1MID0gJyc7CgogIGNvbnN0IHBpY2tlZFByb21pc2UgPSBuZXcgUHJvbWlzZSgocmVzb2x2ZSkgPT4gewogICAgaW5wdXRFbGVtZW50LmFkZEV2ZW50TGlzdGVuZXIoJ2NoYW5nZScsIChlKSA9PiB7CiAgICAgIHJlc29sdmUoZS50YXJnZXQuZmlsZXMpOwogICAgfSk7CiAgfSk7CgogIGNvbnN0IGNhbmNlbCA9IGRvY3VtZW50LmNyZWF0ZUVsZW1lbnQoJ2J1dHRvbicpOwogIGlucHV0RWxlbWVudC5wYXJlbnRFbGVtZW50LmFwcGVuZENoaWxkKGNhbmNlbCk7CiAgY2FuY2VsLnRleHRDb250ZW50ID0gJ0NhbmNlbCB1cGxvYWQnOwogIGNvbnN0IGNhbmNlbFByb21pc2UgPSBuZXcgUHJvbWlzZSgocmVzb2x2ZSkgPT4gewogICAgY2FuY2VsLm9uY2xpY2sgPSAoKSA9PiB7CiAgICAgIHJlc29sdmUobnVsbCk7CiAgICB9OwogIH0pOwoKICAvLyBXYWl0IGZvciB0aGUgdXNlciB0byBwaWNrIHRoZSBmaWxlcy4KICBjb25zdCBmaWxlcyA9IHlpZWxkIHsKICAgIHByb21pc2U6IFByb21pc2UucmFjZShbcGlja2VkUHJvbWlzZSwgY2FuY2VsUHJvbWlzZV0pLAogICAgcmVzcG9uc2U6IHsKICAgICAgYWN0aW9uOiAnc3RhcnRpbmcnLAogICAgfQogIH07CgogIGNhbmNlbC5yZW1vdmUoKTsKCiAgLy8gRGlzYWJsZSB0aGUgaW5wdXQgZWxlbWVudCBzaW5jZSBmdXJ0aGVyIHBpY2tzIGFyZSBub3QgYWxsb3dlZC4KICBpbnB1dEVsZW1lbnQuZGlzYWJsZWQgPSB0cnVlOwoKICBpZiAoIWZpbGVzKSB7CiAgICByZXR1cm4gewogICAgICByZXNwb25zZTogewogICAgICAgIGFjdGlvbjogJ2NvbXBsZXRlJywKICAgICAgfQogICAgfTsKICB9CgogIGZvciAoY29uc3QgZmlsZSBvZiBmaWxlcykgewogICAgY29uc3QgbGkgPSBkb2N1bWVudC5jcmVhdGVFbGVtZW50KCdsaScpOwogICAgbGkuYXBwZW5kKHNwYW4oZmlsZS5uYW1lLCB7Zm9udFdlaWdodDogJ2JvbGQnfSkpOwogICAgbGkuYXBwZW5kKHNwYW4oCiAgICAgICAgYCgke2ZpbGUudHlwZSB8fCAnbi9hJ30pIC0gJHtmaWxlLnNpemV9IGJ5dGVzLCBgICsKICAgICAgICBgbGFzdCBtb2RpZmllZDogJHsKICAgICAgICAgICAgZmlsZS5sYXN0TW9kaWZpZWREYXRlID8gZmlsZS5sYXN0TW9kaWZpZWREYXRlLnRvTG9jYWxlRGF0ZVN0cmluZygpIDoKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgJ24vYSd9IC0gYCkpOwogICAgY29uc3QgcGVyY2VudCA9IHNwYW4oJzAlIGRvbmUnKTsKICAgIGxpLmFwcGVuZENoaWxkKHBlcmNlbnQpOwoKICAgIG91dHB1dEVsZW1lbnQuYXBwZW5kQ2hpbGQobGkpOwoKICAgIGNvbnN0IGZpbGVEYXRhUHJvbWlzZSA9IG5ldyBQcm9taXNlKChyZXNvbHZlKSA9PiB7CiAgICAgIGNvbnN0IHJlYWRlciA9IG5ldyBGaWxlUmVhZGVyKCk7CiAgICAgIHJlYWRlci5vbmxvYWQgPSAoZSkgPT4gewogICAgICAgIHJlc29sdmUoZS50YXJnZXQucmVzdWx0KTsKICAgICAgfTsKICAgICAgcmVhZGVyLnJlYWRBc0FycmF5QnVmZmVyKGZpbGUpOwogICAgfSk7CiAgICAvLyBXYWl0IGZvciB0aGUgZGF0YSB0byBiZSByZWFkeS4KICAgIGxldCBmaWxlRGF0YSA9IHlpZWxkIHsKICAgICAgcHJvbWlzZTogZmlsZURhdGFQcm9taXNlLAogICAgICByZXNwb25zZTogewogICAgICAgIGFjdGlvbjogJ2NvbnRpbnVlJywKICAgICAgfQogICAgfTsKCiAgICAvLyBVc2UgYSBjaHVua2VkIHNlbmRpbmcgdG8gYXZvaWQgbWVzc2FnZSBzaXplIGxpbWl0cy4gU2VlIGIvNjIxMTU2NjAuCiAgICBsZXQgcG9zaXRpb24gPSAwOwogICAgd2hpbGUgKHBvc2l0aW9uIDwgZmlsZURhdGEuYnl0ZUxlbmd0aCkgewogICAgICBjb25zdCBsZW5ndGggPSBNYXRoLm1pbihmaWxlRGF0YS5ieXRlTGVuZ3RoIC0gcG9zaXRpb24sIE1BWF9QQVlMT0FEX1NJWkUpOwogICAgICBjb25zdCBjaHVuayA9IG5ldyBVaW50OEFycmF5KGZpbGVEYXRhLCBwb3NpdGlvbiwgbGVuZ3RoKTsKICAgICAgcG9zaXRpb24gKz0gbGVuZ3RoOwoKICAgICAgY29uc3QgYmFzZTY0ID0gYnRvYShTdHJpbmcuZnJvbUNoYXJDb2RlLmFwcGx5KG51bGwsIGNodW5rKSk7CiAgICAgIHlpZWxkIHsKICAgICAgICByZXNwb25zZTogewogICAgICAgICAgYWN0aW9uOiAnYXBwZW5kJywKICAgICAgICAgIGZpbGU6IGZpbGUubmFtZSwKICAgICAgICAgIGRhdGE6IGJhc2U2NCwKICAgICAgICB9LAogICAgICB9OwogICAgICBwZXJjZW50LnRleHRDb250ZW50ID0KICAgICAgICAgIGAke01hdGgucm91bmQoKHBvc2l0aW9uIC8gZmlsZURhdGEuYnl0ZUxlbmd0aCkgKiAxMDApfSUgZG9uZWA7CiAgICB9CiAgfQoKICAvLyBBbGwgZG9uZS4KICB5aWVsZCB7CiAgICByZXNwb25zZTogewogICAgICBhY3Rpb246ICdjb21wbGV0ZScsCiAgICB9CiAgfTsKfQoKc2NvcGUuZ29vZ2xlID0gc2NvcGUuZ29vZ2xlIHx8IHt9OwpzY29wZS5nb29nbGUuY29sYWIgPSBzY29wZS5nb29nbGUuY29sYWIgfHwge307CnNjb3BlLmdvb2dsZS5jb2xhYi5fZmlsZXMgPSB7CiAgX3VwbG9hZEZpbGVzLAogIF91cGxvYWRGaWxlc0NvbnRpbnVlLAp9Owp9KShzZWxmKTsK",
       "headers": [
        [
         "content-type",
         "application/javascript"
        ]
       ],
       "ok": true,
       "status": 200,
       "status_text": "OK"
      }
     }
    },
    "colab_type": "code",
    "id": "tbAk6wxEJt8J",
    "outputId": "43ee9205-39a1-4f74-dd6a-1a2ce5c89030"
   },
   "outputs": [],
   "source": [
    "def preprocess_data(filename):\n",
    "    # reading csv file\n",
    "    dataFrameInput = pd.read_csv(filename) #outputs a dataframe object\n",
    "\n",
    "    idx = 0\n",
    "    xDataList = list()\n",
    "    yDataList = list()\n",
    "    data = dataFrameInput.values  # numpy array\n",
    "\n",
    "    ## Somehow t = 0 has an odd number of elements, just cut it out\n",
    "\n",
    "    for row in data:\n",
    "        if (row[2] > 0):\n",
    "            if (idx % 2 == 0):\n",
    "                xDataList.append(row)\n",
    "            else:\n",
    "                yDataList.append(row)\n",
    "\n",
    "            idx = idx + 1\n",
    "\n",
    "    # merge the lists\n",
    "    DataList = list()\n",
    "    for rowX, rowY in zip(xDataList, yDataList):\n",
    "        DataList.append([rowX, rowY])\n",
    "\n",
    "    # Shuffle data\n",
    "\n",
    "    random.shuffle(DataList)\n",
    "\n",
    "    DataArray = np.asarray(DataList)\n",
    "    #print(DataArray.shape)\n",
    "\n",
    "    # Strip off header information, i.e. the first 3 cols\n",
    "    DataArraySlim = DataArray[:, :, 3:]\n",
    "    #print(DataArraySlim.shape)\n",
    "\n",
115
    "    '''\n",
116
117
118
119
120
121
122
123
124
125
126
    "    # split in train and test data (ratio 4:1)\n",
    "    DataTrain = DataArraySlim[:4 * int(DataArraySlim.shape[0] / 5)]\n",
    "    DataTest = DataArraySlim[4 * int(DataArraySlim.shape[0] / 5):]\n",
    "\n",
    "    # Split in x (input) and y (output) data\n",
    "    xDataTrain = DataTrain[:, 0, :]\n",
    "    yDataTrain = DataTrain[:, 1, :]\n",
    "    xDataTest = DataTest[:, 0, :]\n",
    "    yDataTest = DataTest[:, 1, :]\n",
    "\n",
    "    #Normalize Input\n",
127
128
129
130
131
132
133
134
135
136
137
138
    "    #xDataTrain = normalize(xDataTrain, axis=1, norm='l1')\n",
    "    #xDataTest = normalize(xDataTest, axis=1, norm='l1')\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    return DataArraySlim\n",
    "\n",
    "def split_data(Data, ratio): # 0<ratio<1 ratio of test data\n",
    "    \n",
    "    # split in train and test data (ratio 4:1)\n",
    "    DataTrain = Data[:int((1-ratio)* int(Data.shape[0] ))]\n",
    "    DataTest = Data[int((1-ratio)* int(Data.shape[0] )):]\n",
139
    "\n",
140
141
142
143
144
145
146
    "    # Split in x (input) and y (output) data\n",
    "    xDataTrain = DataTrain[:, 0, :]\n",
    "    yDataTrain = DataTrain[:, 1, :]\n",
    "    xDataTest = DataTest[:, 0, :]\n",
    "    yDataTest = DataTest[:, 1, :]\n",
    "    \n",
    "    return (xDataTrain,yDataTrain,xDataTest,yDataTest)\n"
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "N2NZasiqKFm_"
   },
   "source": [
    "Review Training Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the Data\n",
166
    "(xDataTrain,yDataTrain,xDataTest,yDataTest) = split_data(preprocess_data(\"trainNN.csv\"),0.2) \n"
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zmHBl9XYLa4q"
   },
   "source": [
    "## Start the Neural Network ##\n",
    "First, load the needed python modules\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Idh6eJsZUhZj"
   },
   "source": [
    "##Build the Network - Take one: Simple FCN (Fully connected Network) ###\n",
    "\n",
    "* Input: Solution of the Kinetic solver at time-step n and a cell, called u \n",
    "* Output: Legendre multiplier of the entropy closure, called alpha\n",
    "\n",
    "u and alpha have the same length, which is given by the size of the moment basis.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "_piriejwUY6I"
   },
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "from tensorflow import keras\n",
    "from tensorflow.keras import layers\n",
    "from tensorflow.keras.utils import plot_model\n",
    "# Build the network:\n",
    "def create_model():\n",
    "  model = tf.keras.models.Sequential([\n",
    "    keras.layers.Dense(256, activation='relu', input_shape=(4,)),\n",
    "    keras.layers.Dropout(0.2),\n",
    "    keras.layers.Dense(512, activation='relu', input_shape=(64,)),\n",
    "    keras.layers.Dropout(0.2),\n",
    "    keras.layers.Dense(256, activation='relu', input_shape=(256,)),\n",
    "    keras.layers.Dropout(0.2),\n",
    "    keras.layers.Dense(128, activation='relu', input_shape=(128,)),\n",
    "    keras.layers.Dropout(0.2),\n",
    "    keras.layers.Dense(4,)\n",
    "  ])\n",
    "\n",
    "  model.summary()\n",
    "  #plot_model(model, to_file='model_plot.png', show_shapes=True, show_layer_names=True)\n",
    "  '''\n",
    "  bfgs_opt = tfp.optimizer.bfgs_minimize(\n",
    "    value_and_gradients_function, initial_position, tolerance=1e-08, x_tolerance=0,\n",
    "    f_relative_tolerance=0, initial_inverse_hessian_estimate=None,\n",
    "    max_iterations=50, parallel_iterations=1, stopping_condition=None,\n",
    "    validate_args=True, max_line_search_iterations=50, name=None\n",
    "    )\n",
    "  '''\n",
    "  model.compile(loss=tf.keras.losses.MeanSquaredError(), optimizer='adam', metrics=['accuracy'])\n",
    "\n",
    "  return model\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "j7bZnd0QWdou"
   },
   "source": [
    "### Compile the model ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 495
    },
    "colab_type": "code",
    "id": "_UJ-kksOWdMN",
    "outputId": "cb9e61a8-9dc4-4d2b-871b-28722435db0b"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model: \"sequential\"\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "dense (Dense)                (None, 256)               1280      \n",
      "_________________________________________________________________\n",
      "dropout (Dropout)            (None, 256)               0         \n",
      "_________________________________________________________________\n",
      "dense_1 (Dense)              (None, 512)               131584    \n",
      "_________________________________________________________________\n",
      "dropout_1 (Dropout)          (None, 512)               0         \n",
      "_________________________________________________________________\n",
      "dense_2 (Dense)              (None, 256)               131328    \n",
      "_________________________________________________________________\n",
      "dropout_2 (Dropout)          (None, 256)               0         \n",
      "_________________________________________________________________\n",
      "dense_3 (Dense)              (None, 128)               32896     \n",
      "_________________________________________________________________\n",
      "dropout_3 (Dropout)          (None, 128)               0         \n",
      "_________________________________________________________________\n",
      "dense_4 (Dense)              (None, 4)                 516       \n",
      "=================================================================\n",
      "Total params: 297,604\n",
      "Trainable params: 297,604\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "model = create_model()\n",
    "\n",
    "# Let's use the Adam optimizer for learning\n",
    "#model.compile(loss=tf.keras.losses.MeanSquaredError(), optimizer='adam', metrics=['accuracy'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HSIghzZCW41n"
   },
   "source": [
    "## Train the model ##\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 862
    },
    "colab_type": "code",
    "id": "tA1P9EwUW9bh",
    "outputId": "76cd7455-2388-4cb6-8d3b-6d8a6cd8f6b4"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model: \"sequential_1\"\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "dense_5 (Dense)              (None, 256)               1280      \n",
      "_________________________________________________________________\n",
      "dropout_4 (Dropout)          (None, 256)               0         \n",
      "_________________________________________________________________\n",
      "dense_6 (Dense)              (None, 512)               131584    \n",
      "_________________________________________________________________\n",
      "dropout_5 (Dropout)          (None, 512)               0         \n",
      "_________________________________________________________________\n",
      "dense_7 (Dense)              (None, 256)               131328    \n",
      "_________________________________________________________________\n",
      "dropout_6 (Dropout)          (None, 256)               0         \n",
      "_________________________________________________________________\n",
      "dense_8 (Dense)              (None, 128)               32896     \n",
      "_________________________________________________________________\n",
      "dropout_7 (Dropout)          (None, 128)               0         \n",
      "_________________________________________________________________\n",
      "dense_9 (Dense)              (None, 4)                 516       \n",
      "=================================================================\n",
      "Total params: 297,604\n",
      "Trainable params: 297,604\n",
      "Non-trainable params: 0\n",
354
      "_________________________________________________________________\n"
355
356
357
358
359
     ]
    },
    {
     "data": {
      "text/plain": [
360
       "\"\\nhistory = model.fit(xDataTrain, yDataTrain,validation_split=0.33, epochs=1, batch_size=1000, verbose=1)\\n\\nprint(history.history.keys())\\nplt.plot(history.history['accuracy'])\\nplt.plot(history.history['val_accuracy'])\\nplt.title('model accuracy')\\nplt.ylabel('accuracy')\\nplt.xlabel('epoch')\\nplt.legend(['train', 'test'], loc='upper left')\\nplt.show()\\n# summarize history for loss\\nplt.plot(history.history['loss'])\\nplt.plot(history.history['val_loss'])\\nplt.title('model loss')\\nplt.ylabel('loss')\\nplt.xlabel('epoch')\\nplt.legend(['train', 'test'], loc='upper left')\\nplt.show()\\n\""
361
362
      ]
     },
363
364
365
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
366
367
368
369
370
371
372
373
374
375
376
    }
   ],
   "source": [
    "#from tensorflow.keras.callbacks import ReduceLROnPlateau\n",
    "\n",
    "#rlrop = ReduceLROnPlateau(monitor='val_loss', factor=0.1, patience=100)\n",
    "\n",
    "\n",
    "# Create the model\n",
    "model = create_model()\n",
    "\n",
377
378
    "'''\n",
    "history = model.fit(xDataTrain, yDataTrain,validation_split=0.33, epochs=1, batch_size=1000, verbose=1)\n",
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
    "\n",
    "print(history.history.keys())\n",
    "plt.plot(history.history['accuracy'])\n",
    "plt.plot(history.history['val_accuracy'])\n",
    "plt.title('model accuracy')\n",
    "plt.ylabel('accuracy')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'test'], loc='upper left')\n",
    "plt.show()\n",
    "# summarize history for loss\n",
    "plt.plot(history.history['loss'])\n",
    "plt.plot(history.history['val_loss'])\n",
    "plt.title('model loss')\n",
    "plt.ylabel('loss')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'test'], loc='upper left')\n",
395
396
    "plt.show()\n",
    "'''"
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Uw9gWsZm7gXr"
   },
   "source": [
    "## Evaluate the Model accuracy\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "FkW8CoZl7rNn"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
422
423
424
      "5857/5857 [==============================] - 5s 812us/step - loss: 117.9833 - accuracy: 0.7676\n",
      "Test score: 117.9832992553711\n",
      "Test accuracy: 0.767593502998352\n"
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
     ]
    }
   ],
   "source": [
    "score = model.evaluate(xDataTest, yDataTest)\n",
    "print('Test score:', score[0])\n",
    "print('Test accuracy:', score[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "G3z1HCkGQJP_"
   },
   "source": [
    "##Make a prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Q9EU3H2hQNgX"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 4)\n",
458
      "[[-2.2854481e-04 -6.9507652e-05 -5.4984037e-05  5.2336970e-04]]\n"
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
     ]
    }
   ],
   "source": [
    "print(np.asarray([xDataTest[0]]).shape)\n",
    "pred = model.predict(np.asarray([xDataTest[0]]))\n",
    "print(pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "33xsSsD5YeA0"
   },
   "source": [
    "## Save trained model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "iCgbAY22YdFo"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
491
      "2.2.0\n"
492
493
494
495
496
497
     ]
    }
   ],
   "source": [
    "print(tf.version.VERSION)\n",
    "\n",
498
499
    "#!mkdir -p saved_model\n",
    "#model.save('saved_model/my_model')\n"
500
501
502
503
504
505
506
507
508
509
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "4mJTKJVzwxE5"
   },
510
   "outputs": [],
511
512
513
   "source": [
    "# Prepare model for download\n",
    "\n",
514
    "#!zip -r model.zip saved_model"
515
516
517
518
   ]
  },
  {
   "cell_type": "code",
519
   "execution_count": 11,
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "zxdU5YCdxGMo"
   },
   "outputs": [],
   "source": [
    "# download model\n",
    "#from google.colab import files\n",
    "#files.download(\"model.zip\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "-EcApWoWf8xk"
   },
   "source": [
    "## Take two: A more sophisticated Loss Function ##\n",
    "\n",
    "The Loss is defined as the Objective function of the minimal entropy problem:\n",
    "\n",
    "$L(u,\\alpha, \\theta) = <\\eta_*(\\alpha(\\theta) \\cdot m(v))>_v - \\alpha(\\theta) \\cdot u$\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Easiest case \n",
    "$\\alpha, u\\in\\mathbb{R}$\n",
    "\n",
    "Entropy $\\eta(\\alpha\\cdot m) = exp(\\alpha\\cdot m)$ \n",
    "\n",
    "Spherical harmonics basis for $M_0$: $ m_0 = \\sqrt{\\left(\\frac{1}{4\\pi}\\right)}$"
   ]
  },
  {
   "cell_type": "code",
563
   "execution_count": 12,
564
565
   "metadata": {},
   "outputs": [
566
567
568
569
570
571
572
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(37.846815845621244, shape=(), dtype=float64)\n"
     ]
    },
573
574
    {
     "data": {
575
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAa0klEQVR4nO3de5hddX3v8fdnZjJJJveQSYBcSBCIRFSEAcGKgnhBSku1jxXsOaJV46UeL6dHHtHzqKcen1MPbX1a20eLB0QrxUIVq0etQVqJPeU2aIAE5G5ukGQy9/v1e/7Ye2AY555Za+291+f1PPNkr7XX2r/vmj3zyW9+e63fUkRgZmb5UZV1AWZmli4Hv5lZzjj4zcxyxsFvZpYzDn4zs5xx8JuZ5YyD38wsZxz8ZimRdKqkPknfyroWyzcHv1l6/ha4L+sizBz8ZkWSQtIpY5ZvlPQ/5+m1rwDagDvm4/XMjoWD32wGJG2S1DbF1zum2Hc58KfAf02vYrPJ1WRdgFk5iIh9wMo57v554PqIOCBp/ooymyMHv9k8kvRj4ILi4vuBPcDrgVdkVpTZOA5+s+f1AHVjlo8HDkBhqAd4eIp93x8RN0XEm8eulPQxYDOwr9jbXwpUS9oWEWfNX+lmM+fgN3veLuAdkvYAbwBeCzTCc0M9S+fwmtcB3x6z/N8o/EfwwWMp1OxY+MNds+d9FPgdCmff/CHwvWN9wYjoiYhDo19AF9AXEU3H+tpmcyXfiMXMLF/c4zczyxkHv5lZzjj4zcxyxsFvZpYzZXE655o1a2Lz5s1Zl2FmVlbuv//+oxFRP359WQT/5s2baWxszLoMM7OyImnvROs91GNmljMOfjOznHHwm5nljIPfzCxnHPxmZjnj4DczyxkHv5lZzjj4zcxK0KH2Pv5ix6M81dQ176+dWPBLukHSEUm7x6z7nKSDknYVvy5Nqn0zs3K2r6WHL//rExxs6533106yx38jcMkE678UEWcWv36UYPtmZmWru38IgCUL53+ChcSCPyJ2Ai1Jvb6ZWSXrKgb/snIK/il8WNKDxaGgVRm0b2ZW8rrKscc/ia8ALwLOBJ4F/mKyDSVtl9QoqbGpybcnNbN8GR3qWbqozIM/Ig5HxHBEjABfA86dYtvrIqIhIhrq639jVlEzs4r2XI+/tsyDX9IJYxbfAuyebFszszzr6hti8YJqqqs076+d2Hz8km4GLgTWSDoAfBa4UNKZQAC/Bt6fVPtmZuWse2AokWEeSDD4I+LKCVZfn1R7ZmaVpKt/mKUJfLALvnLXzKwkdfUNOvjNzPKku3+YJQurE3ltB7+ZWQnq6h9yj9/MLE8c/GZmOdPdP5TIVbvg4DczK0md7vGbmeXH4PAIA0MjDn4zs7xIckpmcPCbmZWczr7kJmgDB7+ZWcnpHigGv3v8Zmb54KEeM7OceW6ox8FvZpYP3f3DgIPfzCw3nh/q8Vw9Zma50PncjdYXJPL6Dn4zsxLjHr+ZWc509A5SV1tNTXUyEe3gNzMrMZ19QyxL6OItcPCbmZWcjr5Bli9KZnwfEgx+STdIOiJp9wTP/YmkkLQmqfbNzMpVOff4bwQuGb9S0kbgjcC+BNs2MytbHX2DLF9chj3+iNgJtEzw1JeAq4FIqm0zs3LW0TvIsnIc6pmIpMuBgxHxwAy23S6pUVJjU1NTCtWZmZWGzr4hlpfpUM8LSKoDPgV8ZibbR8R1EdEQEQ319fXJFmdmViIigo6+yunxvwjYAjwg6dfABuAXko5PsQYzs5LWPzTC4HCwfHFyPf7kXnmciHgIWDu6XAz/hog4mlYNZmalrqN3EKA8e/ySbgbuArZKOiDpPUm1ZWZWKTqKUzInOcaf2CtHxJXTPL85qbbNzMpVR1+hx1+WF3CZmdnsjd6EJckxfge/mVkJKesxfjMzm73nevwOfjOzfBgd4y/XuXrMzGyWOvsGqa4SdbXJ3IQFHPxmZiWlo7cwM6ekxNpw8JuZlZDOhOfiBwe/mVlJ6Uh4Ln5w8JuZlZSOXvf4zcxyJem7b4GD38yspCQ9JTM4+M3MSkp77yAr6xz8Zma50D80TM/AMKsc/GZm+dDeU7hqd0VdbaLtOPjNzEpEW3GCtpWL3eM3M8uFtmKP32P8ZmY50dYzAMDKxR7qMTPLheeGetzjNzPLh+c/3C3T4Jd0g6QjknaPWfd5SQ9K2iVph6QTk2rfzKzctPUOUF0lli0s3yt3bwQuGbfu2oh4WUScCfxf4DMJtm9mVlbaegZZuXhBolMyQ4LBHxE7gZZx6zrGLC4BIqn2zczKTVvvYOLDPADJ/j0xAUlfAN4JtAMXTbHddmA7wKZNm9IpzswsQ209A4mfww8ZfLgbEZ+OiI3ATcCHp9juuohoiIiG+vr69Ao0M8tIW88gKxO+aheyPavnJuD3M2zfzKykjI7xJy3V4Jd06pjFy4Ffpdm+mVkpay/3MX5JNwMXAmskHQA+C1wqaSswAuwFPpBU+2Zm5WRweISu/iFWpTDUk1jwR8SVE6y+Pqn2zMzKWXtKV+2Cr9w1MysJoxO0rai0MX4zM5tYe29xgrYKP6vHzMyKWrvTmYsfHPxmZiWhpbvQ41+9xD1+M7NcaC4G/3FLHfxmZrnQ0t3PogVV1NUmP5OOg9/MrAQ0dw9w3JKFqbTl4DczKwHNXQOpjO+Dg9/MrCS0dDv4zcxypaV7gOMc/GZm+dHc3Z/KGT3g4Dczy1zPwBB9gyOs9oe7Zmb50NxVPIffQz1mZvnQnOJVu+DgNzPLXEt3PwCrPcZvZpYPHuoxM8uZNCdoAwe/mVnmWroHqK2pYunC5OfpAQe/mVnmmosXb0lKpb3Egl/SDZKOSNo9Zt21kn4l6UFJt0lamVT7ZmblormrP7VhHki2x38jcMm4dbcDZ0TEy4DHgGsSbN/MrCw0dfVTvyydi7cgweCPiJ1Ay7h1OyJiqLh4N7AhqfbNzMpFU2c/9UsrIPhn4I+AH0/2pKTtkholNTY1NaVYlplZeoZHgqNdA6xdXuHBL+nTwBBw02TbRMR1EdEQEQ319fXpFWdmlqLWngGGRyLVHv+Mzh2SdCrwv4BtwKLR9RFx8mwblPQu4DLg4oiI2e5vZlZJjnQUrtpdu3zRNFvOn5n2+L8OfIVCL/0i4JvAt2bbmKRLgKuB342Intnub2ZWaZq6isFfgh/uLo6IOwBFxN6I+Bzw21PtIOlm4C5gq6QDkt4D/A2wDLhd0i5JXz2G2s3Myt6Rjj6AVM/qmellYv2SqoDHJX0YOAgsnWqHiLhygtXXz7I+M7OKNtrjL8XTOT8K1AEfAc4G/jNwVVJFmZnlxZGOfpYurKGuNp3pGmCGPf6IuK/4sAt4d3LlmJnlS1NXf6rj+zDzs3pOAz4BnDR2n4h4XUJ1mZnlQlNHP2tKMfiBW4GvAl8DhpMrx8wsX5q6+nnJictTbXOmwT8UEV9JtBIzsxw60tHHRVvXptrmlMEvaXXx4Q8kfQi4DegffT4iWibc0czMptXdP0T3wHCqZ/TA9D3++4EARieJ/sSY5wKY9ZW7ZmZWcKQz/Yu3YJrgj4gtaRViZpY3z7b3AnDCivSma4CZn9WzCPgQ8GoKPf2fA1+NiL4EazMzq2jPthUi9ISVi1Ntd6Yf7n4T6AS+XFx+B/D3wNuSKMrMLA8OFadrOD7FCdpg5sF/RkRsG7P8b5IeTqIgM7O8eKatl1V1C1hcW51quzOdsuEXks4bXZD0SqAxmZLMzPLhUHsfx69Id5gHZt7jPxv4D0n7isubgEclPQRE8R66ZmY2C8+093Fiyh/swsyDf/xN083M7Bgdau/lrE0rU293phdwdU70vC/gMjObm96BYVp7Bjkx5TN6YPYXcI3eKlH4Ai4zsznL6owemMUFXMXe/6mMueeumZnNzXMXb60sseAfJem9FG7GsgHYBZwH/AdwcWKVmZlVsOcu3srgrJ7Z3IHrHGBvRFwEvAJon2oHSTdIOiJp95h1b5O0R9KIpIY5V21mVuZGh3rSnq4BZh78faPTM0haGBG/ArZOs8+N/ObZQLuBtwI7Z1OkmVmlOdDay3FLalm0IN2Lt2Dmp3MekLQS+B5wu6RWYO9UO0TETkmbx617BEDShPuYmeXFgdYeNqyuy6Ttmd5z9y3Fh5+T9G/ACuBfEqvKzKzC7W/p4Yz1KzJpe6ZDPc+JiDsj4vsRMZBEQaMkbZfUKKmxqakpyabMzFI1PBIcbOtlw6psevyzDv60RMR1EdEQEQ319fVZl2NmNm8Od/QxOBxsXJ3+GT1QwsFvZlap9rf0ALCx0nr8km4G7gK2Sjog6T2S3iLpAHA+8ENJP0mqfTOzUrW/tXDx1sZS/nB3LiLiykmeui2pNs3MysH+lh4kODGDq3bBQz1mZqnb39rD8csXsbAm/XP4wcFvZpa6Ay29mY3vg4PfzCx1+1t72JDRGT3g4DczS1Xf4DCHOvo4afWSzGpw8JuZpejXzd1EwJZ6B7+ZWS483dQNwMlrHPxmZrnw1NFC8G9x8JuZ5cPTR7tZt3whSxYmdhnVtBz8ZmYpevpod6a9fXDwm5ml6qmmLrasWZppDQ5+M7OUtHYP0NozmOkHu+DgNzNLzdPN2X+wCw5+M7PUPFU8lTPLc/jBwW9mlprHD3dSW13FSRlNxzzKwW9mlpJHD3fyorVLqanONnod/GZmKXnsUCdb12V7Rg84+M3MUtHRN8gz7X2cdvyyrEtx8JuZpeHxw50AbF3n4Dczy4XHDncBcFolB7+kGyQdkbR7zLrVkm6X9Hjx31VJtW9mVkoePdTJktpq1q/M7gYso5Ls8d8IXDJu3SeBOyLiVOCO4rKZWcV77HAnp6xbRlWVsi4lueCPiJ1Ay7jVlwPfKD7+BvB7SbVvZlYqIoI9z3Sw7YTlWZcCpD/Gvy4ini0+PgSsm2xDSdslNUpqbGpqSqc6M7MEHGjtpb13kJeuX5F1KUCGH+5GRAAxxfPXRURDRDTU19enWJmZ2fx66GA7AGesz2eP/7CkEwCK/x5JuX0zs9Q9dLCdmiqxtQTO4Yf0g//7wFXFx1cB/5xy+2Zmqdt9sJ3T1i1jYU111qUAyZ7OeTNwF7BV0gFJ7wH+DHiDpMeB1xeXzcwqVkSw+2B7yYzvAyR208eIuHKSpy5Oqk0zs1JzsK2X1p7BkhnfB1+5a2aWqF/uawPgZRtWZlrHWA5+M7ME3b+3lcULqtl2onv8Zma5cP/eVl6+cQULMp6Df6zSqcTMrML0DAzx8LMdNJy0OutSXsDBb2aWkAf2tzM8Epx9UmnNR+ngNzNLyP17C9OVnbXJwW9mlguNe1s5de1SVtQtyLqUF3Dwm5klYHB4hHufbuH8Fx2XdSm/wcFvZpaAX+5ro2dgmN86ZU3WpfwGB7+ZWQL+/YmjVAnOO9k9fjOzXPh/TxzlZRtWsmJxaY3vg4PfzGzedfYNsmt/G68uwWEecPCbmc27u55sZngkeNUppTfMAw5+M7N599NHDrNsUU3JXbE7ysFvZjaPhkeCOx45wkVb11JbU5oRW5pVmZmVqV/sa6W5e4A3vmRd1qVMysFvZjaPduw5xIJq8drT6rMuZVIOfjOzeRIR7Hj4MOe/aA3LFpXeaZyjHPxmZvNk1/429jb3cNnLTsi6lCllEvySPippt6Q9kj6WRQ1mZvPttl8eZGFNFW8+4/isS5lS6sEv6QzgfcC5wMuByySdknYdZmbzaXB4hB888Axv2LaupId5IJse/+nAPRHRExFDwJ3AWzOow8xs3tz5aBOtPYO85RXrsy5lWlkE/27gAknHSaoDLgU2jt9I0nZJjZIam5qaUi/SzGw2br53H2uW1vKaEj6bZ1TqwR8RjwBfBHYA/wLsAoYn2O66iGiIiIb6+tL/RppZfu1v6eFfHz3CFedsKqmbqk8mkwoj4vqIODsiXgO0Ao9lUYeZ2Xz41j17EfCOV27KupQZqcmiUUlrI+KIpE0UxvfPy6IOM7Nj1TswzC337ecN29Zx4srFWZczI5kEP/AdSccBg8AfR0RbRnWYmR2Tm+/dR2vPIO959clZlzJjmQR/RFyQRbtmZvOpf2iYv9v5JOduWc25W0pzJs6JlP6nEGZmJerWxgMc7ujnI687NetSZsXBb2Y2B939Q/zVHY9z9kmr+K0SveHKZBz8ZmZz8Hd3PklTZz+f/u3TkZR1ObPi4Dczm6UDrT1c9/On+J2Xn8hZm1ZlXc6sOfjNzGYhIrjmuw9RLfHJN78463LmxMFvZjYLt95/gJ8/fpRPXno668vkvP3xHPxmZjP09NFuPv+Dhzl3y2r+8NzyuEp3Ig5+M7MZ6B0Y5oPfup+aavGlt59JVVV5faA7VlZX7pqZlY2RkeDq7zzIo4c7+fq7zinbIZ5R7vGbmU0hIvjCjx7hBw88w9VvejEXbl2bdUnHzMFvZjaJiOAvb3+M6//9ad71qs184LXlMx/PVDzUY2Y2gZGR4H/8YA/fuGsvV5yzkc9ctq3sLtSajIPfzGyctp4BPvaPu/jZo02874ItfOrS8rs6dyoOfjOzMRp/3cLHb9nFofY+Pv97Z/CfXrmpokIfHPxmZgB09A3ylzse4xt3/ZoTVyzmH99/fllOxzATDn4zy7W+wWH+/q69/O3PnqC9d5B3nncSn7jkxSxdWLnxWLlHZmY2hWfbe7np7n38w737aOke4DWn1XP1m7ZyxvoVWZeWOAe/meVGc1c/P9lzmB8+9Ax3PdlMABe/eB3vvWAL551cXnPqHwsHv5lVpIjgcEc/Dx1s5+6nmrnryWYeOdRBBGxZs4QPXXgKf9CwkU3H1WVdauoyCX5JHwfeCwTwEPDuiOjLohYzK2+9A8McaO1hf2sP+5p72N/ay2OHO3n4mQ6auwcAqK2p4uxNq/j460/j9aev4/QTllXcmTqzkXrwS1oPfATYFhG9km4BrgBuTLsWM0vH8EgwODzC0EgwPBwMjowwNBwMjfl3cDjoGxymZ2D0a4iegWG6+4foHRima2CItu5BmrsHaOnup6V7gObuATr7hl7Q1qIFVZyydikXn76WbScs5yXrV/DS9StYtKA6o6MvPVkN9dQAiyUNAnXAM0k08td3PM73H5j8pSNi8ueme/FpNpjq6ananUnb0+xOTPEK0+477YFPtW8JH9cxtT3NcR1z23P/OZz+/czmZxwKYT80Ugj7Y/m5GlVbXcXKugWsXlLLcUtreemqlRy3pJY1S2vZuLqODavq2Lh6MfVLF+a6Nz8TqQd/RByU9OfAPqAX2BERO8ZvJ2k7sB1g06a5zXu9dtlCtq5bNvVGU/x8TPejM90P11TPTvdzmWTb0724ptlgqqanrzvBtqf9XS/T45q27bmHXJJ1VwlqqqtYUCWqq6qoqRYLqguPF1SLmuK6mio9t92iBdXU1VZTV1tD3cJqltTWsLi2sG5BtacWmy+arpc27w1Kq4DvAG8H2oBbgX+KiG9Ntk9DQ0M0NjamU6CZWYWQdH9ENIxfn8V/oa8Hno6IpogYBL4LvCqDOszMcimL4N8HnCepToW/US8GHsmgDjOzXEo9+CPiHuCfgF9QOJWzCrgu7TrMzPIqk7N6IuKzwGezaNvMLO/8MbmZWc44+M3McsbBb2aWMw5+M7OcSf0CrrmQ1ATsnePua4Cj81hOOfGx51Nejz2vxw2TH/tJEVE/fmVZBP+xkNQ40ZVreeBj97HnSV6PG2Z/7B7qMTPLGQe/mVnO5CH483xVsI89n/J67Hk9bpjlsVf8GL+Zmb1QHnr8ZmY2hoPfzCxnKir4Jb1N0h5JI5Iaxj13jaQnJD0q6U1j1l9SXPeEpE+mX/X8k3SmpLsl7ZLUKOnc4npJ+uvisT4o6aysa51vkv6LpF8Vfw7+95j1E77/lUbSn0gKSWuKy3l4z68tvucPSrpN0soxz1X8+z6nDIuIivkCTge2Aj8DGsas3wY8ACwEtgBPAtXFryeBk4Ha4jbbsj6Oefg+7ADeXHx8KfCzMY9/TOFufucB92Rd6zwf90XAT4GFxeW1U73/WdebwPFvBH5C4WLHNXl4z4vH+Eagpvj4i8AX8/K+zzXDKqrHHxGPRMSjEzx1OfDtiOiPiKeBJ4Bzi19PRMRTETEAfLu4bbkLYHnx8Qqev5n95cA3o+BuYKWkE7IoMCEfBP4sIvoBIuJIcf1k73+l+RJwNS+8D3qlv+dExI6IGCou3g1sKD7Ow/s+pwyrqOCfwnpg/5jlA8V1k60vdx8DrpW0H/hz4Jri+ko93lGnARdIukfSnZLOKa6v9ONG0uXAwYh4YNxTFX/s4/wRhb9wIB/HPqdjzORGLMdC0k+B4yd46tMR8c9p15OVqb4PFG5n+fGI+I6kPwCup3Cv47I3zXHXAKspDGmcA9wi6eQUy0vUNMf+KQpDHhVpJr/3kj4NDAE3pVlbOSq74I+IuQTYQQrjn6M2FNcxxfqSNtX3QdI3gY8WF28F/k/x8VTfh7IwzXF/EPhuFAY/75U0QmHyqrI/bpj82CW9lMIY9gOF21izAfhF8UP9ij72UZLeBVwGXFx8/6FCjn0aczrGvAz1fB+4QtJCSVuAU4F7gfuAUyVtkVQLXFHcttw9A7y2+Ph1wOPFx98H3lk80+M8oD0ins2iwIR8j8IHvEg6jcKHXUeZ/P2vCBHxUESsjYjNEbGZwp/7Z0XEISr/PUfSJRQ+2/jdiOgZ81RFv+9Fc8qwsuvxT0XSW4AvA/XADyXtiog3RcQeSbcAD1P4U/CPI2K4uM+HKZwJUQ3cEBF7Mip/Pr0P+CtJNUAfsL24/kcUzvJ4AugB3p1NeYm5AbhB0m5gALiq2Pub9P3PgUp/zwH+hsKZO7cX/+K5OyI+MNXvfaWIiKG5ZJinbDAzy5m8DPWYmVmRg9/MLGcc/GZmOePgNzPLGQe/mVnOOPjNpiCpK+sazOabg9/MLGcc/GYzULzy9VpJuyU9JOntxfUnSNpZvPfBbkkXSKqWdOOYbT+edf1mY1XUlbtmCXorcCbwcgrz/9wnaSfwDuAnEfEFSdVAXXG79RFxBsDYG4OYlQL3+M1m5tXAzRExHBGHgTspzAB6H/BuSZ8DXhoRncBTwMmSvlycR6Yjq6LNJuLgNzsGEbETeA2FGRFvlPTOiGil8JfBz4AP8PzsqGYlwcFvNjM/B95eHL+vpxD290o6CTgcEV+jEPBnFe93WxUR3wH+O1Bx97m18uYxfrOZuQ04n8I9TQO4OiIOSboK+ISkQaALeCeFOyB9XdJox+qaiV7QLCuendPMLGc81GNmljMOfjOznHHwm5nljIPfzCxnHPxmZjnj4DczyxkHv5lZzvx/I0k2ahz7yDoAAAAASUVORK5CYII=\n",
576
577
578
579
580
581
582
583
584
585
586
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
587
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcS0lEQVR4nO3de5RdZZnn8e+vrrmSIqQMMRcDEtS0XIwFgtLSkLEFdDrYSxFvpBl6MtrYQ6vdNmrPtDOrZ7WMY9ti22gUNXRrO4gC0WYUGkVlNEhAbuEyRExMQhISyL1SdS71zB/7rcMhJKlKqva5VP0+a9U6e797n32ezSH11Pu8795bEYGZmRlAS70DMDOzxuGkYGZmFU4KZmZW4aRgZmYVTgpmZlbhpGBmZhVOCmZmVuGkYJYjSR2SbpK0TlJI+r0DtkvSNZKeTT/XSFJ9ojVzUjCrhbuB9wJbDrJtGXAxcBpwKvDvgf9Us8jMDuCkYHYI6S/7k6rWvy7pb47kGBFRiIi/j4i7gfJBdlkKfCYiNkbEJuAzwB+NJG6zkXBSMDsKkuZJ2nmYn3cP81C/AzxYtf5gajOri7Z6B2DWjCLit0DXKBxqCrCran0XMEWSwjcmszpwT8GsvvYCx1StHwPsdUKwenFSMDu0XmBS1frxgwupfLT3MD/vGeZnrCEbZB50WmozqwuXj8wO7QHg3ZLWAG8CzgVWQ6V8NGU4B5HUCQxOM+2QNAHoT72BG4APS7oNCOAjwOdH8yTMjoR7CmaHdhXZFNGdwHuAW47yOE8A+4HZwA/T8svSti8B3wMeBh4B/jW1mdWFXLo0M7NB7imYmVmFk4KZmVU4KZiZWYWTgpmZVTT1lNQZM2bE/Pnz6x2GmVlTue+++7ZHRPfBtjV1Upg/fz6rV6+udxhmZk1F0vpDbXP5yMzMKnJNCpK60gNGHpf0mKSzJU2XdIekJ9PrsWlfSbpW0lpJD0lalGdsZmb2Ynn3FD4H/CAiXkl2T5fHgKuBOyNiAXBnWge4EFiQfpYB1+Ucm5mZHSC3pCBpGvBG4HqoPGxkJ7AEWJF2W0H21ClS+w2RWQV0SZqVV3xmZvZiefYUTgC2AV+T9CtJX5E0GZgZEZvTPluAmWl5NrCh6v0bU9sLSFomabWk1du2bcsxfDOz8SfPpNAGLAKui4jXAPt4vlQEQLpL5BHdfCkilkdET0T0dHcfdEaVmZkdpTyTwkZgY0Tck9ZvIksSWwfLQun1mbR9EzC36v1zUpuZmdVIbkkhIrYAGyS9IjUtBh4FVpI9rJz0emtaXglclmYhnQXsqiozmZkZsK+/xGduf4IHN+zM5fh5X7z2p8A3JHUATwGXkyWiGyVdAawHLkn73gZcBKwle+LV5TnHZmbWdPb2l/j8j9Zy/LQJnDa3a9SPn2tSiIgHgJ6DbFp8kH0DuDLPeMzMml2xPABAe0s+hR5f0Wxm1kRK5WxuTlurhtjz6DgpmJk1kdJA1lNoa3VPwcxs3CuUsp5Ch3sKZmZW6Sl4TMHMzIoeUzAzs0GlwdlHHlMwM7PSQNZTcFIwMzMKlZ6Cy0dmZuNeseTykZmZJYMDzU4KZmZWmZLq8pGZmVFw+cjMzAa5fGRmZhVFzz4yM7NBlaTQ5p6Cmdm4N1g+6nD5yMzMBnsKbS0uH5mZjXul8gAStDopmJlZoRy0t7YgOSmYmY17xfIA7Tn1EsBJwcysqRTLA7nNPAInBTOzplJM5aO8OCmYmTWRYnkgt+mo4KRgZtZUiuWB3B7FCTknBUnrJD0s6QFJq1PbdEl3SHoyvR6b2iXpWklrJT0kaVGesZmZNaNieaDpy0fnRcTpEdGT1q8G7oyIBcCdaR3gQmBB+lkGXFeD2MzMmspYHFNYAqxIyyuAi6vab4jMKqBL0qw6xGdm1rCynkKTlo+AAG6XdJ+kZaltZkRsTstbgJlpeTawoeq9G1ObmZkleZeP2nI7cuaciNgk6SXAHZIer94YESEpjuSAKbksA5g3b97oRWpm1gSy8lGT9hQiYlN6fQa4GTgT2DpYFkqvz6TdNwFzq94+J7UdeMzlEdETET3d3d15hm9m1nCadqBZ0mRJUweXgd8HHgFWAkvTbkuBW9PySuCyNAvpLGBXVZnJzMxo7vLRTODmdNOmNuCbEfEDSfcCN0q6AlgPXJL2vw24CFgL9AKX5xibmVlTKpbyLR/llhQi4ingtIO0PwssPkh7AFfmFY+Z2VhQHGjS8pGZmY0+3+bCzMwqiqVo3ttcmJnZ6Gra2UdmZjb6nBTMzKyiqS9eMzOz0eWegpmZARARlAbG3l1SzczsKBTL2a3iOvyMZjMzK5YHAGhr8ZiCmdm4N5gUXD4yM7NK+ajd5SMzMxvsKXR4SqqZmT0/puCegpnZuFcZU3D5yMzMKlNSXT4yMzOXj8zMrKJQcvnIzMySwaTQ6aRgZmb9g1NSnRTMzGywp+DHcZqZGf0pKUxod1IwMxv3nu8ptOb2GU4KZmZNopIUPKZgZmaFUhlwUjAzM6AwFmYfSWqV9CtJ30/rJ0i6R9JaSf9bUkdq70zra9P2+XnHZmbWTPqLY+M6hauAx6rWrwE+GxEnATuAK1L7FcCO1P7ZtJ+ZmSWF8gBSEz95TdIc4C3AV9K6gPOBm9IuK4CL0/KStE7avjjtb2ZmZAPNHa0t5PmrMe+ewt8DHwUG0vpxwM6IKKX1jcDstDwb2ACQtu9K+7+ApGWSVktavW3bthxDNzNrLP2lgVzHEyDHpCDprcAzEXHfaB43IpZHRE9E9HR3d4/moc3MGlp/aSDX8QSAthyP/QbgDyRdBEwAjgE+B3RJaku9gTnAprT/JmAusFFSGzANeDbH+MzMmkqhNEBnW34XrkGOPYWI+FhEzImI+cClwI8i4j3Aj4G3p92WArem5ZVpnbT9RxERecVnZtZsCuUmLh8dxl8CH5a0lmzM4PrUfj1wXGr/MHB1HWIzM2tYhVI515vhQb7lo4qIuAu4Ky0/BZx5kH36gHfUIh4zs2ZUaOaBZjMzG121GGh2UjAzaxLuKZiZWcVYHWg2M7OjMHhFc56cFMzMmoTLR2ZmVtHfzBevmZnZ6Grqex+ZmdnoKpTKnpJqZmYZzz4yMzMAIiIrH3n2kZmZlQaCiHwfxQlOCmZmTaFQyp5V5vKRmZk5KZiZ2fMKZScFMzNL+otZUvDFa2ZmRqFcBtxTMDMzsquZAU9JNTMz6EvlowntTgpmZuNefzErH01o95iCmdm411eqTVJoG85OkhYAfwssBCYMtkfEiTnFZWZmVRqtfPQ14DqgBJwH3AD8c15BmZnZC+0vpJ5Cg0xJnRgRdwKKiPUR8UngLfmFZWZm1QbLRxM7GqB8BPRLagGelPRBYBMwJb+wzMysWqV81CA9hauAScB/Bl4LvA9YmldQZmb2Qn1p9lFnzmMKw+opRMS9aXEvcPlw3iNpAvBToDN9zk0R8deSTgC+BRwH3Ae8LyIKkjrJxipeCzwLvDMi1h3BuZiZjVn9xTJSg9w6W9LJkr4s6XZJPxr8GeJt/cD5EXEacDpwgaSzgGuAz0bEScAO4Iq0/xXAjtT+2bSfmZkBfaUBOttakJTr5wx3TOHbwBeBLwPl4bwhIoKsZwHQnn4COB94d2pfAXySbGbTkrQMcBPwD5KUjmNmNq71Fcu5X6MAw08KpYi47kgPLqmVrER0EvAF4NfAzogopV02ArPT8mxgA0BElCTtIisxbT/gmMuAZQDz5s070pDMzJrS/kKZiTVICoctH0maLmk68D1JfyJp1mBbaj+siChHxOnAHOBM4JUjDTgilkdET0T0dHd3j/RwZmZNoa800BA9hfvISj6DRay/qNoWwLCuaI6InZJ+DJwNdElqS72FOWTTW0mvc4GNktqAaWQDzmZm415fsZz7IDMMkRQi4oSjPbCkbqCYEsJE4E1kg8c/Bt5ONgNpKXBresvKtP6LtP1HHk8wM8s01JhCml76J8A5ZD2EnwFfjIi+w7xtFrAijSu0ADdGxPclPQp8S9LfAL8Crk/7Xw/8k6S1wHPApUdzQmZmY1F/cSD3+x7B8AeabwD2AJ9P6+8G/gl4x6HeEBEPAa85SPtTZOMLB7b3He54ZmbjWV+pzPTJHbl/znCTwqsjYmHV+o/TX/xmZlYDfcUGmH1U5f504RkAkl4HrM4nJDMzO1BfsTFmHw16LfBzSb9N6/OAJyQ9THad2qm5RGdmZgDsL5YbakzhglyjMDOzw8qmpNa5p1B1gdqeg22PiOdGPSIzM3uR/gYpHx148drgdQPiCC5eMzOzo1ceCArlBpiSWn3xWuo1LKDqGc1mZpa//sGnrjVATwEASX9M9qCdOcADwFnAz4HFuUVmZmZA1VPXGmhK6lXAGcD6iDiP7KK0XblFZWZmFYNPXatF+Wi4n9A3eEsLSZ0R8TjwivzCMjOzQfsrSaFBykdkdy7tAm4B7pC0A1ifV1BmZva8/YUGSwoR8ba0+Ml0C+xpwA9yi8rMzCoGewqTO4b7d/zRO+JPiIif5BGImZkd3L7+7GGVkzobZ6DZzMzqpDeVjyZ1OCmYmY17g0mhFuUjJwUzswbXW0jlI/cUzMxsX/9g+cg9BTOzcW9/oYTUWBevmZlZnewrlJnc0YakoXceIScFM7MG11soMbEG4wngpGBm1vB6C2UmOymYmRlkA80TazDIDE4KZmYNr7dQck/BzMwyvYUykzqbvKcgaa6kH0t6VNIaSVel9umS7pD0ZHo9NrVL0rWS1kp6SNKivGIzM2smvYUSk2pwh1TIt6dQAj4SEQvJntR2paSFwNXAnRGxALgzrQNcSPa4zwXAMuC6HGMzM2sa+/rLNbkZHuSYFCJic0Tcn5b3AI8Bs4ElwIq02wrg4rS8BLghMquALkmz8orPzKxZZGMKTV4+qiZpPtkjPO8BZkbE5rRpCzAzLc8GNlS9bWNqMzMb13oL5Zrc9whqkBQkTQG+A/xZROyu3hYRAcQRHm+ZpNWSVm/btm0UIzUzazyl8gD9pYGa3PcIck4KktrJEsI3IuK7qXnrYFkovT6T2jcBc6vePie1vUBELI+Inojo6e7uzi94M7MG0Fus3bMUIN/ZRwKuBx6LiL+r2rQSWJqWlwK3VrVflmYhnQXsqiozmZmNS4PPZ67VQHOe/ZE3AO8DHpb0QGr7OPAp4EZJVwDrgUvSttuAi4C1QC9weY6xmZk1hcqjOGvUU8gtKUTE3cChbum3+CD7B3BlXvGYmTWjPX1ZUpja2V6Tz/MVzWZmDWxv6ilMnTAGBprNzGxk9vQVAZjipGBmZrtT+eiYCS4fmZmNe3v7XD4yM7NkcKB5crPfJdXMzEZuT1+Rie2ttLfW5te1k4KZWQPb21+qWekInBTMzBranr5SzWYegZOCmVlD29NfYmqNZh6Bk4KZWUPb01fkGPcUzMwMsvKRxxTMzAzIrlOYUqPpqOCkYGbW0Pb0FT2mYGZm2VPX9hXKLh+ZmRns2p/dDO/YSR01+0wnBTOzBrUzJYWuSS4fmZmNezt7CwB0uadgZmY7e1NPYaJ7CmZm414lKbh8ZGZmz48puHxkZjbu7eot0CKY6ovXzMxsR2+RaRPbaWlRzT7TScHMrEHt3F+saekInBTMzBrWzt5CTQeZwUnBzKxh7dqflY9qKbekIOmrkp6R9EhV23RJd0h6Mr0em9ol6VpJayU9JGlRXnGZmTWLZ/cWmD557JSPvg5ccEDb1cCdEbEAuDOtA1wILEg/y4DrcozLzKzhRQTb9vbTPaWzpp+bW1KIiJ8Czx3QvARYkZZXABdXtd8QmVVAl6RZecVmZtbo9vSXKJQGmDFWksIhzIyIzWl5CzAzLc8GNlTttzG1vYikZZJWS1q9bdu2/CI1M6uj7Xv6AZgxdeyUjw4rIgKIo3jf8ojoiYie7u7uHCIzM6u/7Xuzm+EdN3ls9xS2DpaF0uszqX0TMLdqvzmpzcxsXNq+N/UUxnj5aCWwNC0vBW6tar8szUI6C9hVVWYyMxt3KkmhxuWj3G6oIelfgN8DZkjaCPw18CngRklXAOuBS9LutwEXAWuBXuDyvOIyM2sG2/f0I8H0Gl/RnFtSiIh3HWLT4oPsG8CVecViZtZstu0tMH1SB22ttS3o+IpmM7MGtG1Pf83HE8BJwcysIW3etZ9ZXRNq/rlOCmZmDWjzrj5mTZtY8891UjAzazB9xTLP7Svw0mnuKZiZjXubd/UB8NIu9xTMzMa9zTv3A3hMwczMqnoKHlMwM7OnU0/heI8pmJnZumd7mXlMJxPaW2v+2U4KZmYNZt2z+5h/3OS6fHZut7loZLev2cItDxz+JqxChz/IyDYjHX6Pod8/ss8fjRiG/m8wxPFH/N8w588f4X+Ahj+/IT4/O8Zw9hpJDIfeQYL21hbaW0VbSwvtbaK9pYW2VtHW2kJHam9rFR2tLUzsaGVKZxuTO9sqr5M7W+lsq/1f2yO1bvs+3rRw5tA75mBcJoWdvUWe3Lr3kNuHeshDdqumw2wfKoAhdsj984EhDkEMcZQh33/ET8o48P0jO8fcz2+Enz/UEUb++TX4fyTnfwcDERQHgmJ5YET/P7W3imMnddA9tZMZUzorr7OPncj84yYx/7jJvLRrIq0tI0uAo2V3X5Fn9xWYP8M9hZq55Iy5XHLG3KF3NLOGUE7JoVgeoFQOigMDFMtBqZy9FkoD7C+W2ddfYl9/ib3pdV+hzJ6+Ejv2Fdi2t5/te/t5cusetu3tp1h+PtN0tLZwYvdkTp0zjVPndHH63C5eNeuYuiSKddv3Abh8ZGZ2KK0torWlddQGXiOCrbv7+c32fax7dh/rtu/j8S17uOPRrdy4eiMAXZPaOeekGZx7cjdvWjiTrhrdwvrX27IqxondTgpmZjUhieOnTeD4aRM4++XHVdojgo079nP/b3fwsye385P/t43vP7SZjtYWzntlN3+4aA6LX/mSXG9n/djmPXS0tXCCy0dmZvUlibnTJzF3+iSWnD6biOCRTbu55YFN3PrA0/xwzVbmHDuRK845gUt65jK5c/R/hT769G5eMXMq7TV+jsIgT0k1MzsESZwyZxr/5a0LWfWx8/nie1/LzGMm8N++9yjnfvouvnnPbymVB0bt8yKCRzfvZuGsY0btmEfKScHMbBjaWlu44NXH850PvJ6b3n8284+bxMdvfpgLP/czVq97blQ+Y8vuPp7bV2DhS50UzMyaRs/86Xz7/WfzxfcuordQ5h1f+gWfXLmG3kJpRMe9b/0OAE6dM200wjwqTgpmZkdBEhe8ehY//NAbed9ZL+PrP1/HW6+9m8e37D7qY6566lkmd7RyymwnBTOzpjSls43/vuTVfPOPX8fuvhIXf+H/ctN9G4/qWPc89Rw986fnOrtpKE4KZmaj4PUnzeC2q87h9Lld/Pm3H+TjNz9Mf6k87Pdv2rmfJ5/Zy+urpsjWg5OCmdkoecnUCfzzFa/j/ee+nG/e81suXb6Krbv7hvXeHzyyBYA3/87xeYY4JCcFM7NR1NbawtUXvpJ/fM8intiyh7dcezf3DjE7KSL47v0bWTjrmLrd82iQk4KZWQ4uOmUWt1z5BqZOaONdy1dxwy/WHfImgvf85jnWPL2b95w1r8ZRvlhDJQVJF0h6QtJaSVfXOx4zs5E4eeZUbrnyDZx7cjf/9dY1/Pm3H6Kv+MJxhlJ5gL+97TG6p3byttfMrlOkz2uYpCCpFfgCcCGwEHiXpIX1jcrMbGSmTWzny5f1cNXiBXzn/o28/Ys/55e/eY6IoL9U5hM3P8KDG3fxV295FZM66n/nofpH8LwzgbUR8RSApG8BS4BH6xqVmdkItbSID73pZE6dM42P3vQQl3zpFxw3uYP+0gB7+0t88LyTWHJ6/XsJ0FhJYTawoWp9I/C6A3eStAxYBjBvXv3rb2Zmw7X4VTO5+y/PZ+WDm7hv/Q462lq46JRZvP7lM+odWkUjJYVhiYjlwHKAnp6eET7fy8ystiZ2tPLOM+bxzjMa84/ahhlTADYB1Y9Dm5PazMysRhopKdwLLJB0gqQO4FJgZZ1jMjMbVxqmfBQRJUkfBH4ItAJfjYg1dQ7LzGxcaZikABARtwG31TsOM7PxqpHKR2ZmVmdOCmZmVuGkYGZmFU4KZmZWoUPdta8ZSNoGrD/Kt88Ato9iOI3G59fcfH7NrdHP72UR0X2wDU2dFEZC0uqI6Kl3HHnx+TU3n19za+bzc/nIzMwqnBTMzKxiPCeF5fUOIGc+v+bm82tuTXt+43ZMwczMXmw89xTMzOwATgpmZlYxLpKCpHdIWiNpQFLPAds+JmmtpCckvbmq/YLUtlbS1bWP+uhIOl3SKkkPSFot6czULknXpvN5SNKiesd6tCT9qaTH03f6P6vaD/pdNiNJH5EUkmak9ab//iR9On1vD0m6WVJX1bYx8d016++NF4iIMf8DvAp4BXAX0FPVvhB4EOgETgB+TXbb7ta0fCLQkfZZWO/zGOa53g5cmJYvAu6qWv4/gICzgHvqHetRnt95wL8BnWn9JYf7Lusd71Ge41yyW8ivB2aMle8P+H2gLS1fA1wzlr67Zv69Uf0zLnoKEfFYRDxxkE1LgG9FRH9E/AZYC5yZftZGxFMRUQC+lfZtBgEck5anAU+n5SXADZFZBXRJmlWPAEfoA8CnIqIfICKeSe2H+i6b0WeBj5J9l4Oa/vuLiNsjopRWV5E9XRHGznfXzL83KsZFUjiM2cCGqvWNqe1Q7c3gz4BPS9oA/C/gY6m9mc+p2snA70q6R9JPJJ2R2sfE+UlaAmyKiAcP2DQmzq/KfyDr+cDYObcxcR4N9ZCdkZD0b8DxB9n0iYi4tdbx5Olw5wosBj4UEd+RdAlwPfDvahnfSA1xfm3AdLISyhnAjZJOrGF4IzbE+X2crMzSlIbz71DSJ4AS8I1axmbDM2aSQkQczS++TWT120FzUhuHaa+7w52rpBuAq9Lqt4GvpOXDnWtDGeL8PgB8N7Ii7i8lDZDdfKzpz0/SKWQ19QclQXYO96fJAk1xfkP9O5T0R8BbgcXpO4QmObdhGBPnMd7LRyuBSyV1SjoBWAD8ErgXWCDpBEkdwKVp32bwNHBuWj4feDItrwQuS7NYzgJ2RcTmegQ4QreQDTYj6WSyAb3tHPq7bBoR8XBEvCQi5kfEfLLyw6KI2MIY+P4kXUA2VvIHEdFbtanpv7ukmX9vVIyZnsLhSHob8HmgG/hXSQ9ExJsjYo2kG4FHybqzV0ZEOb3ng2QzQFqBr0bEmjqFf6T+I/A5SW1AH7Astd9GNoNlLdALXF6f8Ebsq8BXJT0CFICl6S/OQ36XY8RY+P7+gWyG0R2pJ7QqIt5/uH+HzSQiSk38e6PCt7kwM7OK8V4+MjOzKk4KZmZW4aRgZmYVTgpmZlbhpGBmZhVOCmZHQdLeesdglgcnBTMzq3BSMBuBdIXxpyU9IulhSe9M7bMk/TQ91+IRSb8rqVXS16v2/VC94zc70Li4otksR38InA6cRnYPpnsl/RR4N/DDiPgfklqBSWm/2RHxaoDqh8yYNQr3FMxG5hzgXyKiHBFbgZ+Q3b31XuBySZ8ETomIPcBTwImSPp/uA7S7XkGbHYqTglkOIuKnwBvJ7pL5dUmXRcQOsh7FXcD7ef4OtmYNw0nBbGR+BrwzjRd0kyWCX0p6GbA1Ir5M9st/UXrecktEfAf4K6DpnrNsY5/HFMxG5mbgbLLn8Qbw0YjYImkp8BeSisBe4DKyp3B9TdLgH2MfO9gBzerJd0k1M7MKl4/MzKzCScHMzCqcFMzMrMJJwczMKpwUzMyswknBzMwqnBTMzKzi/wOsLEoq04ViOwAAAABJRU5ErkJggg==\n",
588
589
590
591
592
593
594
595
596
597
598
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
599
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWt0lEQVR4nO3df5RfdX3n8ec7mWSSEJIYM2AIkoQWqBRFNPyote0RtLK1LW1dW9Qq1XY5tbW1Ho8eXbtbztn2rF277lp3jyxaFarVrT/LVq1YK2KPFQgI8rsgEggQmEBIZjIz35nvzHv/+N4J4zQJX8Lce7/zvc/HOXNyv/d753vfN3fmlU8+997PJzITSVJzLKm7AElStQx+SWoYg1+SGsbgl6SGMfglqWEMfklqGINfkhrG4JeeoYg4LyLujIixiPhmRGw+zLb/JSJuiYh2RFxSYZnSAQa/9AxExAbgC8B/AtYD24H/e5hvuQd4F/Dl8quTDs7gV2NFREbEj895/YmI+NOn+TG/BtyWmZ/NzAngEuD0iPiJg22cmZdn5leBkSOtW3qmDH7pICLihIh44jBfrys2/Ung5tnvy8z9wA+K9VJPGqi7AKkXZeb9wLouNl0NDM9btxc4eqFrkhaKLX7pmRkF1sxbtwa7ctTDDH412Riwas7r58wuFF09o4f5en2x6W3A6XO+7yjgx4r1Uk8y+NVkNwGvi4ilEXE+8HOzb2Tm/Zm5+jBfnyo2/SJwWkS8OiJWAP8Z+H5m3nmwHUbEsmK7JcBARKyIiKWlHqU0j8GvJnsb8EvAE8DrgS893Q/IzGHg1cCfAXuAs4ELZ9+PiEsj4tI53/IRYBx4LfDeYvkNR1S9dITCiVgkqVls8UtSwxj8ktQwBr8kNYzBL0kNsyie3N2wYUNu2bKl7jIkaVG54YYbdmfm0Pz1iyL4t2zZwvbt2+suQ5IWlYjYcbD1dvVIUsMY/JLUMAa/JDWMwS9JDWPwS1LDGPyS1DAGvyQ1jMEvST3orl0jfOCqu9g92lrwzzb4JakH3blrH3/5T/ewd3xqwT/b4JekHtSamgFgcGDhY7q04I+Ij0XEoxFx65x16yPi6xFxd/Hns8ravyQtZq3pTvAvX0zBD3wCOH/euncD38jMk4BvFK8lSfNMtmdb/As/JXNpwZ+Z1wCPz1t9AXB5sXw58Ctl7V+SFrNWexpYZF09h3BsZj5cLO8Cjq14/5K0KMy2+JcvXfzBf0B2Znk/5EzvEXFxRGyPiO3Dw8MVViZJ9Wu1Z1i2NFiyJBb8s6sO/kciYiNA8eejh9owMy/LzG2ZuW1o6N/MIyBJfW2yPVNKax+qD/4rgYuK5YuAv6t4/5K0KEy2ZxhctvAXdqHc2zk/DfwLcEpE7IyI3wbeB7wiIu4GXl68liTN02pPl9biL23qxcx87SHeOq+sfUpSv5hsz5RyDz/45K4k9aTJ6ZlSbuUEg1+SelJryha/JDWKLX5Jahhb/JLUMK3pGZaXME4PGPyS1JMm23b1SFKjtNrTdvVIUpPY4pekhjH4JalhWn00SJskqQuLcpA2SdKRm5y2xS9JjdGenmF6Ju3jl6SmmJwupl00+CWpGVpTBr8kNcpsi3/QIRskqRkm27b4JalRWu1pAC/uSlJTtGzxS1Kz2NUjSQ0z2+If9AEuSWqGA8HvkA2S1AwTU52LuyuW2eKXpEZ4Mvht8UtSI8w+uWvwS1JDTBT38a/wrh5JaobZrh4v7kpSQ0zMdvXY4pekZpiYmmZgSTDgffyS1Ayt9kxpF3ahpuCPiLdHxG0RcWtEfDoiVtRRhyT1oomp6dLu4Ycagj8iNgF/CGzLzNOApcCFVdchSb1qYmqmtLH4ob6ungFgZUQMAKuAh2qqQ5J6zkS7z1r8mfkg8BfA/cDDwN7MvGr+dhFxcURsj4jtw8PDVZcpSbVpTU33Vx9/RDwLuADYChwHHBURvzl/u8y8LDO3Zea2oaGhqsuUpNp0unr6qMUPvBz4YWYOZ+YU8AXgJTXUIUk9aaLfWvx0unjOiYhVERHAecAdNdQhST2p727nzMxrgc8BNwK3FDVcVnUdktSryr6dc6C0Tz6MzPwT4E/q2Lck9bqJ9jQr+vB2TknSIUxMzZQ2QBsY/JLUc/ruyV1J0uG1+vTJXUnSQUzPJJPTM7b4JakpJtvlTrsIBr8k9ZQDE6332ZO7kqRDODDfri1+SWqGA9MuGvyS1AwHunq8uCtJzTAb/N7OKUkNMdvVM2iLX5KaoeXFXUlqlidv5zT4JakRxiY7wb9qucEvSY1g8EtSw4wXwb/S4JekZhifmm3xlzdPlsEvST1kbHKa5QNLWLokStuHwS9JPWR8sl1q/z4Y/JLUU8Ymp1lZ4j38YPBLUk8Zm5ou9cIuGPyS1FPGJ6ft6pGkJhmfnGbVsvLu6AGDX5J6il09ktQw3tUjSQ3jXT2S1DDjk3b1SFKjjHlXjyQ1R2YyPjXNyhLH6YGagj8i1kXE5yLizoi4IyJ+qo46JKmXzE67WHaLv9x/Vg7tg8A/ZOa/j4jlwKqa6pCknjE22QYo/eJu5cEfEWuBnwV+CyAzJ4HJquuQpF4zVsFY/FBPV89WYBj4eER8LyI+GhFH1VCHJPWUJ8fi77/gHwBeBHw4M88A9gPvnr9RRFwcEdsjYvvw8HDVNUpS5cYrmHYR6gn+ncDOzLy2eP05Ov8Q/IjMvCwzt2XmtqGhoUoLlKQ6HOjq6bexejJzF/BARJxSrDoPuL3qOiSp14xPdS7u9utdPX8AfKq4o+de4E011SFJPaOqi7u1BH9m3gRsq2PfktSrxlr928cvSTqI0Vanq2f1YLlt8q4+PSJOAv4rcCqwYnZ9Zp5YUl2S1Dj7i+A/quTg77bF/3Hgw0AbeBlwBfDJsoqSpCYabbUZHFjCsqXldsZ0++krM/MbQGTmjsy8BHhVeWVJUvOMttqld/NA9xd3WxGxBLg7It4KPAisLq8sSWqe/a126d080H2L/210BlL7Q+DFwBuAi8oqSpKaaLSi4O9qD5l5fbE4ivfcS1IpRlttju6V4I+Ik4F3Apvnfk9mnltSXZLUOPtb02xYvbz0/XT7T8tngUuBjwDT5ZUjSc21v9Vm87PLn56k2+BvZ+aHS61EkhpupBfu6omI9cXi/4uI3wO+CLRm38/Mx0usTZIaZX8vBD9wA5BAFK/fOee9BHxyV5IWwMxMMjY5Xf9dPZm5tfQKJEnsn6xmnB7o/q6eFcDvAS+l09L/NnBpZk6UWJskNcZoReP0QPcXd68ARoAPFa9fB/w18JoyipKkpnlygLZyh2SG7oP/tMw8dc7rb0aEs2ZJ0gIZLcbiP3pF7wzZcGNEnDP7IiLOBraXU5IkNc/oRNHiX947XT0vBr4TEfcXr08A7oqIW4DMzBeUUp0kNUQv9vGfX2oVktRw+yuafQu6f4Br5GDv+wCXJC2MXmrxz3+AK4s/Ax/gkqQFMzIxBcCalT30AFfR+j+JOXPuSpIWxt7xKVYsW8LgQI/czhkRv0NnMpbjgZuAc4DvAOeVVpkkNci+8TZrViyrZF9PZwauM4Edmfky4Axgb2lVSVLD7JuYYs3K3gr+idnhGSJiMDPvBE4pryxJapZ9E1OsqeDhLej+ds6dEbEO+BLw9YjYA+woqyhJapp94+1KZt+C7ufc/dVi8ZKI+CawFviH0qqSpIbZNzHFiUNHVbKvp/3/isz8VhmFSFKT7R2f6rmLu5KkkmQm+8anKrmHHwx+Sard/slpZpL+b/FHxNKI+F5E/H1dNUhSL9g3PvvUbp8HP51nA+6ocf+S1BP2FcM1rO3n4I+I44FXAR+tY/+S1Ev2jhUt/j7v6vmfwLuAmUNtEBEXR8T2iNg+PDxcWWGSVLV9xSQsfXtxNyJ+EXg0M2843HaZeVlmbsvMbUNDQxVVJ0nVO9DH38ct/p8Gfjki7gM+A5wbEZ+soQ5J6gn7Jvr84m5mviczj8/MLcCFwD9l5m9WXYck9Yonij7+KiZaB+/jl6TaPTE2yZoVAyxbWk0kV/PPyyFk5tXA1XXWIEl1e3xsivVHVTNAG9jil6Ta7dk/ybMMfklqjj1jkzxrlcEvSY2xZ7/BL0mN8vjYJOuPquZWTjD4JalW45PTTEzN2McvSU2xZ2wSgPV29UhSMzy+vxP8tvglqSFmW/xe3JWkhthTDNfgxV1Jaog9+23xS1KjzPbxVzX7Fhj8klSr3aMt1h+1nIGKBmgDg1+SajU80mJo9WCl+zT4JalGw6Mtho42+CWpMYZHWhxj8EtSM2Rmp6vH4JekZtg30abVnjH4JakphkdaAAa/JDXFgeD3rh5JaobhUVv8ktQodvVIUsM8OjLB8qVLKh2uAQx+SarNI3snGDp6kIiodL8GvyTV5KG9E2xat7Ly/Rr8klSTh54Y57h1Kyrfr8EvSTWYnkl27Z3gOFv8ktQMu0dbtGeSjQa/JDXDg0+MA7DJrh5JaoaHiuBvRFdPRDw3Ir4ZEbdHxG0R8baqa5Ckuj38xAQAG9dWH/wDle8R2sA7MvPGiDgauCEivp6Zt9dQiyTV4sEnxlk9OMCaFdXHcOUt/sx8ODNvLJZHgDuATVXXIUl12rlnnE3rVlb+8BbU3McfEVuAM4BrD/LexRGxPSK2Dw8PV16bJJVpx2P72fzsVbXsu7bgj4jVwOeBP8rMffPfz8zLMnNbZm4bGhqqvkBJKsnMTLLj8TG2bjiqlv3XEvwRsYxO6H8qM79QRw2SVJeH900w2Z5h87MbEvzR6dD6K+COzPxA1fuXpLrt2L0fgC0N6ur5aeANwLkRcVPx9Qs11CFJtbjvsTEANtfU1VP5fUSZ+c9A9ZexJalH7HhsP8sHlrBxTfVP7YJP7kpS5e7dvZ/N61exZEk9bWCDX5IqdvcjI5x07Ora9m/wS1KFxien2fH4GKccu6a2Ggx+SarQ3Y+OkAmnPMcWvyQ1wl27RgA4+dija6vB4JekCt21a4TBgSW1PbwFBr8kVequ4sLu0pru6AGDX5Iqk5l8f+denr9pba11GPySVJH7Hhtj7/gUpx+/rtY6DH5JqshND+wB4IUnrKu1DoNfkipy8wN7WbV8KScdU98dPWDwS1Jlvnf/Hk7btLbWC7tg8EtSJfZNTHHLg3s5Z+v6uksx+CWpCtfe+zgzCS/58Q11l2LwS1IVvvOD3axYtoQzar6wCwa/JFXin+/ezZlb1jM4sLTuUgx+SSrbD3fv5+5HRzn3J46puxTA4Jek0n3ttl0AvPInn1NzJR0GvySV7Ku37uIFx6/luHUr6y4FMPglqVT3PDrCzQ88wauev7HuUg4w+CWpRJ+57gEGlgSvfvHxdZdygMEvSSUZn5zm8zfu5BWnHsuG1YN1l3OAwS9JJfn0dfezZ2yKN790a92l/AiDX5JKMD45zf+55gecvXU9Z26pf5iGuQx+SSrBh6++h0f2tXjHz59Sdyn/hsEvSQvsjof3cek193LBC4/jrB4YlG0+g1+SFtBoq81b/+ZG1q5cxh+/6tS6yzmogboLkKR+MTE1zcVXbOe+x8b46zefxdDRvXMnz1wGvyQtgN2jLd7yyRu4/r49fODXT++J4ZcPxeCXpGcgM7ny5of40y/fwcjEFB967Rn80unH1V3WYdUS/BFxPvBBYCnw0cx8Xx11SNKRGh5pcdXtu7j8O/fxr4+M8oLj13L5m87i1OPW1F3aU6o8+CNiKfC/gVcAO4HrI+LKzLy96lqkXpaZXW7XxTYLvc+uP6+bz1q443w6uv280VabR/ZNMDzSYueeMW5/eIRbH9zLrQ/tJROet3EN//01p/MrZ2yqfS7dbtXR4j8LuCcz7wWIiM8AFwALHvy//zc38u1/HX7K7bo6/13+kNTxC9btD/BC/4J1/XtYwy//Qp6HhQw5LX7rVi3j1I1rePvLT+a85x3DqRvXELE4An9WHcG/CXhgzuudwNnzN4qIi4GLAU444YQj2tFLfuzZDHU5PkY35y3o7uR2+zPQ7Y9KV7V1udOufzy7PoaF+ztZyL+PzudVX1u3xdVxrHX8XHa2W7hQrOPcr1y2lGPXrOCYNYMct3Ylx64ZXHRBP1/PXtzNzMuAywC2bdt2RG2p15+9eUFrkqR+UMcDXA8Cz53z+vhinSSpAnUE//XASRGxNSKWAxcCV9ZQhyQ1UuVdPZnZjoi3Al+jczvnxzLztqrrkKSmqqWPPzO/Anyljn1LUtM5SJskNYzBL0kNY/BLUsMY/JLUMNHt0AF1iohhYMcRfvsGYPcCltPrmna80Lxj9nj730Id8+bMHJq/clEE/zMREdszc1vddVSlaccLzTtmj7f/lX3MdvVIUsMY/JLUME0I/svqLqBiTTteaN4xe7z9r9Rj7vs+fknSj2pCi1+SNIfBL0kN01fBHxGviYjbImImIrbNe+89EXFPRNwVEa+cs/78Yt09EfHu6qteGBHxwoj4bkTcFBHbI+KsYn1ExF8Wx/f9iHhR3bUulIj4g4i4szjn/23O+oOe634REe+IiIyIDcXrvjzHEfH+4vx+PyK+GBHr5rzXl+e4sjzKzL75Ap4HnAJcDWybs/5U4GZgENgK/IDOkNBLi+UTgeXFNqfWfRxHeOxXAf+uWP4F4Oo5y1+lM6PeOcC1dde6QMf7MuAfgcHi9TGHO9d117uAx/1cOkOa7wA29Pk5/nlgoFj+c+DP+/kcV5lHfdXiz8w7MvOug7x1AfCZzGxl5g+Be+hM+n5g4vfMnARmJ35fjBJYUyyvBR4qli8ArsiO7wLrImJjHQUusLcA78vMFkBmPlqsP9S57hf/A3gXPzoHfF+e48y8KjPbxcvv0pmtD/r3HFeWR30V/IdxsAneNx1m/WL0R8D7I+IB4C+A9xTr++kY5zoZ+JmIuDYivhURZxbr+/V4iYgLgAcz8+Z5b/XtMc/xZjr/q4H+Pd7KjqtnJ1s/lIj4R+A5B3nrvZn5d1XXU6XDHTtwHvD2zPx8RPw68FfAy6usb6E9xfEOAOvpdG2cCfxtRJxYYXmleIpj/o90uj/6Rje/zxHxXqANfKrK2vrZogv+zDySMDvcBO+LZuL3wx17RFwBvK14+Vngo8Xyop3c/imO9y3AF7LTOXpdRMzQGdhq0R4vHPqYI+L5dPqzb44I6BzXjcVF/EV7zE/1+xwRvwX8InBeca5hER/vU6jsuJrS1XMlcGFEDEbEVuAk4Dr6a+L3h4CfK5bPBe4ulq8E3ljc+XEOsDczH66jwAX2JToXeImIk+lcDNvNoc/1opaZt2TmMZm5JTO30OkGeFFm7qJPz3FEnE/nesYvZ+bYnLf68hxTYR4tuhb/4UTErwIfAoaAL0fETZn5ysy8LSL+Fridzn8Zfz8zp4vv6ZeJ3/8D8MGIGAAmgIuL9V+hc9fHPcAY8KZ6yltwHwM+FhG3ApPARUWL8JDnuo/16zn+X3Tu3Pl68b+c72bm7x7u93kxy8x2VXnkkA2S1DBN6eqRJBUMfklqGINfkhrG4JekhjH4JalhDH7pMCJitO4apIVm8EtSwxj8UheKp2LfHxG3RsQtEfEbxfqNEXFNMQ/CrRHxMxGxNCI+MWfbt9ddvzRXXz25K5Xo14AXAqfTGRPo+oi4Bngd8LXM/LOIWAqsKrbblJmnAcydQETqBbb4pe68FPh0Zk5n5iPAt+iMCno98KaIuAR4fmaOAPcCJ0bEh4rxZvbVVbR0MAa/9Axk5jXAz9IZRfETEfHGzNxD538GVwO/y5MjpUo9weCXuvNt4DeK/vshOmF/XURsBh7JzI/QCfgXFXPhLsnMzwN/DPTFHLjqH/bxS935IvBTdOZBTeBdmbkrIi4C3hkRU8Ao8EY6syZ9PCJmG1bvOdgHSnVxdE5Jahi7eiSpYQx+SWoYg1+SGsbgl6SGMfglqWEMfklqGINfkhrm/wPGrAXj7qAu/gAAAABJRU5ErkJggg==\n",
600
601
602
603
604
605
606
607
608
609
610
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
611
    "import math\n",
612
613
614
615
    "\n",
    "def loss1D(u, alpha):\n",
    "    return 4*np.pi*np.exp(alpha*np.sqrt(1/(4*np.pi))) - alpha*u\n",
    "\n",
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
    "# Custom Loss\n",
    "def custom_loss1d(u_input, alpha_pred): # (label,prediciton)\n",
    "     return (0.5*alpha_pred*alpha_pred -u*alpha_pred)\n",
    "    #return 4*math.pi*tf.math.exp(alpha_pred*np.sqrt(1/(4*np.pi))) - alpha_pred*u_input\n",
    "\n",
    "def custom_loss1dMB(u_input, alpha_pred): # (label,prediciton)\n",
    "    return 4*math.pi*tf.math.exp(alpha_pred*np.sqrt(1/(4*np.pi))) - alpha_pred*u_input\n",
    "  \n",
    "# Custom Loss\n",
    "def custom_loss1dMBPrime(u_input, alpha_pred): # (label,prediciton)\n",
    "\n",
    "    return 0.5*tf.square(4*math.pi*np.sqrt(1/(4*np.pi))*tf.math.exp(alpha_pred*np.sqrt(1/(4*np.pi))) - u_input)\n",
    "#tf.norm(4*math.pi*np.sqrt(1/(4*np.pi))*tf.math.exp(alpha_pred*np.sqrt(1/(4*np.pi))) - u_input, ord='euclidean')\n",
    "\n",
    "\n",
631
632
    "# not realizable moment\n",
    "u = -4\n",
633
634
635
    "print(custom_loss1dMBPrime(u,1))\n",
    "\n",
    "alpha = np.arange(-100.0, -3, 0.05)\n",
636
    "plt.figure()\n",
637
    "plt.plot(alpha,custom_loss1dMBPrime(u,alpha))\n",
638
639
640
641
642
643
644
    "plt.ylabel('alpha')\n",
    "plt.xlabel('loss')\n",
    "plt.title('u=-4')\n",
    "plt.show()\n",
    "\n",
    "#realizable moment\n",
    "u = 10\n",
645
    "alpha = np.arange(-100.0, 9.0, 0.05)\n",
646
    "plt.figure()\n",
647
    "plt.plot(alpha,custom_loss1dMBPrime(u,alpha))\n",
648
649
650
651
652
653
654
    "plt.ylabel('alpha')\n",
    "plt.xlabel('loss')\n",
    "plt.title('u=10')\n",
    "plt.show()\n",
    "\n",
    "# moment close to the boundary of the realizable set\n",
    "u = 0.1\n",
655
    "alpha = np.arange(-100.0, 1.0, 0.05)\n",
656
    "plt.figure()\n",
657
    "plt.plot(alpha,custom_loss1dMBPrime(u,alpha))\n",
658
659
660
661
662
663
664
665
    "plt.ylabel('alpha')\n",
    "plt.xlabel('loss')\n",
    "plt.title('u=0.1')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
666
   "execution_count": 66,
667
668
669
670
671
672
673
674
675
   "metadata": {},
   "outputs": [],
   "source": [
    "# build the network\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "from tensorflow import keras\n",
    "from tensorflow.keras import layers\n",
    "import math\n",
676
    "from tensorflow.keras.callbacks import EarlyStopping,ModelCheckpoint\n",
677
    "\n",
678
    "def custom_loss1dMB(u_input, alpha_pred): # (label,prediciton)\n",
679
680
    "    return 4*math.pi*tf.math.exp(alpha_pred*np.sqrt(1/(4*np.pi))) - alpha_pred*u_input\n",
    "  \n",
681
682
683
684
685
686
687
688
689
    "# Custom Loss\n",
    "def custom_loss1dMBPrime(u_input, alpha_pred): # (label,prediciton)\n",
    "    return 0.5*tf.square(4*math.pi*np.sqrt(1/(4*np.pi))*tf.math.exp(alpha_pred*np.sqrt(1/(4*np.pi))) - u_input)\n",
    "\n",
    "# Custom Loss\n",
    "def custom_loss1dQuad(u_input, alpha_pred): # (label,prediciton)\n",
    "    return (0.5*alpha_pred*alpha_pred -u*alpha_pred)\n",
    "    #return 4*math.pi*tf.math.exp(alpha_pred*np.sqrt(1/(4*np.pi))) - alpha_pred*u_input\n",
    "  \n",
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
    "\n",
    "\n",
    "# Build the network:\n",
    "def create_model():\n",
    "    \n",
    "    # Define the input\n",
    "    input_ = keras.Input(shape=(1,))\n",
    "    \n",
    "    # Hidden layers\n",
    "    hidden1 = layers.Dense(4, activation=\"tanh\")(input_)  \n",
    "    hidden2 = layers.Dense(8, activation=\"tanh\")(hidden1)  \n",
    "    hidden3 = layers.Dense(32, activation=\"tanh\")(hidden2)  \n",
    "    hidden4 = layers.Dense(8, activation=\"tanh\")(hidden3)  \n",
    "    hidden5 = layers.Dense(4, activation=\"tanh\")(hidden4)  \n",
    "    \n",
    "    #Define the ouput\n",
    "    output_ = layers.Dense(1)(hidden5)\n",
707
    "    \n",
708
709
710
711
712
713
714
    "    # Create the model\n",
    "    model = keras.Model(inputs=[input_], outputs=[output_] )\n",
    "    \n",
    "    \n",
    "    model.summary()\n",
    "    \n",
    "    # tf.keras.losses.MeanSquaredError()\n",
715
716
717
    "    #custom_loss1d\n",
    "    model.compile(loss=custom_loss1dMBPrime , optimizer='adam',\n",
    "              metrics=[custom_loss1dMB, custom_loss1dMBPrime])\n",
718
719
720
721
722
723
    "\n",
    "    return model\n"
   ]
  },
  {
   "cell_type": "code",
724
   "execution_count": 67,
725
726
727
728
729
730
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
731
      "Model: \"model_5\"\n",
732
733
734
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
735
      "input_6 (InputLayer)         [(None, 1)]               0         \n",
736
      "_________________________________________________________________\n",
737
      "dense_40 (Dense)             (None, 4)                 8         \n",
738
      "_________________________________________________________________\n",
739
      "dense_41 (Dense)             (None, 8)                 40        \n",
740
      "_________________________________________________________________\n",
741
      "dense_42 (Dense)             (None, 32)                288       \n",
742
      "_________________________________________________________________\n",
743
      "dense_43 (Dense)             (None, 8)                 264       \n",
744
      "_________________________________________________________________\n",
745
      "dense_44 (Dense)             (None, 4)                 36        \n",
746
      "_________________________________________________________________\n",
747
      "dense_45 (Dense)             (None, 1)                 5         \n",
748
749
750
751
752
753
754
755
756
757
758
759
760
761
      "=================================================================\n",
      "Total params: 641\n",
      "Trainable params: 641\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "model = create_model()\n"
   ]
  },
  {
   "cell_type": "code",
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.1      0.100001 0.100002 ... 9.999997 9.999998 9.999999]\n",
      "[6.842216 4.260433 9.455449 ... 3.610957 2.044917 0.218258]\n"
     ]
    }
   ],
   "source": [
    "# build training data and shuffe!\n",
    "uTrain = np.arange(0.1, 10,0.000001)\n",
    "print(uTrain)\n",
    "random.shuffle(uTrain)\n",
    "print(uTrain)"
   ]
  },
  {
   "cell_type": "markdown",
784
785
   "metadata": {},
   "source": [
786
    "### Study of Optimizers ###\n",
787
    "\n",
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
    "Use the dataset np.arange(0.1, 2,0.000001)\n",
    "* RMSPROP: minimal Loss  4.8240e-05  in 79 epochs\n",
    "* ADAM:    minimal Loss  1.2580e-06 in 68 epochs\n",
    "* SGD:  minimal Loss 6.9069e-05 in 167 epochs\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6.842216 4.260433 9.455449 ... 3.610957 2.044917 0.218258]\n",
      "Epoch 1/1500\n",
      "6930/6930 [==============================] - 9s 1ms/step - loss: 0.0828 - custom_loss1dMB: 8.3398 - custom_loss1dMBPrime: 0.0828 - val_loss: 0.0080 - val_custom_loss1dMB: 8.2817 - val_custom_loss1dMBPrime: 0.0080\n",
      "Epoch 2/1500\n",
      "6930/6930 [==============================] - 9s 1ms/step - loss: 0.0060 - custom_loss1dMB: 8.2853 - custom_loss1dMBPrime: 0.0060 - val_loss: 0.0065 - val_custom_loss1dMB: 8.2806 - val_custom_loss1dMBPrime: 0.0065\n",
      "Epoch 3/1500\n",
      "6930/6930 [==============================] - 9s 1ms/step - loss: 0.0035 - custom_loss1dMB: 8.2837 - custom_loss1dMBPrime: 0.0035 - val_loss: 0.0025 - val_custom_loss1dMB: 8.2781 - val_custom_loss1dMBPrime: 0.0025\n",
      "Epoch 4/1500\n",
      "6930/6930 [==============================] - 9s 1ms/step - loss: 0.0027 - custom_loss1dMB: 8.2831 - custom_loss1dMBPrime: 0.0027 - val_loss: 7.9423e-04 - val_custom_loss1dMB: 8.2772 - val_custom_loss1dMBPrime: 7.9423e-04\n",
      "Epoch 5/1500\n",
      "6930/6930 [==============================] - 9s 1ms/step - loss: 0.0021 - custom_loss1dMB: 8.2828 - custom_loss1dMBPrime: 0.0021 - val_loss: 0.0015 - val_custom_loss1dMB: 8.2775 - val_custom_loss1dMBPrime: 0.0015\n",
      "Epoch 6/1500\n",
      "6930/6930 [==============================] - 9s 1ms/step - loss: 0.0017 - custom_loss1dMB: 8.2825 - custom_loss1dMBPrime: 0.0017 - val_loss: 0.0025 - val_custom_loss1dMB: 8.2781 - val_custom_loss1dMBPrime: 0.0025\n",
      "Epoch 7/1500\n",
      "6930/6930 [==============================] - 9s 1ms/step - loss: 0.0015 - custom_loss1dMB: 8.2824 - custom_loss1dMBPrime: 0.0015 - val_loss: 6.3264e-04 - val_custom_loss1dMB: 8.2771 - val_custom_loss1dMBPrime: 6.3264e-04\n",
      "Epoch 8/1500\n",
      "6930/6930 [==============================] - 9s 1ms/step - loss: 0.0012 - custom_loss1dMB: 8.2823 - custom_loss1dMBPrime: 0.0012 - val_loss: 5.3574e-04 - val_custom_loss1dMB: 8.2770 - val_custom_loss1dMBPrime: 5.3574e-04\n",
      "Epoch 9/1500\n",
      "6930/6930 [==============================] - 9s 1ms/step - loss: 0.0011 - custom_loss1dMB: 8.2822 - custom_loss1dMBPrime: 0.0011 - val_loss: 5.7144e-04 - val_custom_loss1dMB: 8.2771 - val_custom_loss1dMBPrime: 5.7144e-04\n",
      "Epoch 10/1500\n",
      "6930/6930 [==============================] - 9s 1ms/step - loss: 9.2444e-04 - custom_loss1dMB: 8.2822 - custom_loss1dMBPrime: 9.2444e-04 - val_loss: 4.0372e-04 - val_custom_loss1dMB: 8.2769 - val_custom_loss1dMBPrime: 4.0372e-04\n",
      "Epoch 11/1500\n",
      " 558/6930 [=>............................] - ETA: 6s - loss: 9.3096e-04 - custom_loss1dMB: 8.2815 - custom_loss1dMBPrime: 9.3096e-04"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-74-0264eab8996e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     16\u001b[0m \u001b[0;31m#,tf.keras.callbacks.History() ,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     17\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 18\u001b[0;31m \u001b[0mhistory\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0muTrain\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0muTrain\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mvalidation_split\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepochs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1500\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbatch_size\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcallbacks\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mes\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mmc\u001b[0m\u001b[0;34m]\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     19\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     20\u001b[0m \u001b[0;31m# summarize history for loss\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/tensorflow/python/keras/engine/training.py\u001b[0m in \u001b[0;36m_method_wrapper\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m     64\u001b[0m   \u001b[0;32mdef\u001b[0m \u001b[0m_method_wrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     65\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_in_multi_worker_mode\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# pylint: disable=protected-access\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 66\u001b[0;31m       \u001b[0;32mreturn\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\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     67\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     68\u001b[0m     \u001b[0;31m# Running inside `run_distribute_coordinator` already.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/tensorflow/python/keras/engine/training.py\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, x, y, batch_size, epochs, verbose, callbacks, validation_split, validation_data, shuffle, class_weight, sample_weight, initial_epoch, steps_per_epoch, validation_steps, validation_batch_size, validation_freq, max_queue_size, workers, use_multiprocessing)\u001b[0m\n\u001b[1;32m    846\u001b[0m                 batch_size=batch_size):\n\u001b[1;32m    847\u001b[0m               \u001b[0mcallbacks\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mon_train_batch_begin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstep\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 848\u001b[0;31m               \u001b[0mtmp_logs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtrain_function\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0miterator\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    849\u001b[0m               \u001b[0;31m# Catch OutOfRangeError for Datasets of unknown size.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    850\u001b[0m               \u001b[0;31m# This blocks until the batch has finished executing.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/tensorflow/python/eager/def_function.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwds)\u001b[0m\n\u001b[1;32m    578\u001b[0m         \u001b[0mxla_context\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mExit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    579\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 580\u001b[0;31m       \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_call\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\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    581\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    582\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mtracing_count\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_tracing_count\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/tensorflow/python/eager/def_function.py\u001b[0m in \u001b[0;36m_call\u001b[0;34m(self, *args, **kwds)\u001b[0m\n\u001b[1;32m    609\u001b[0m       \u001b[0;31m# In this case we have created variables on the first call, so we run the\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    610\u001b[0m       \u001b[0;31m# defunned version which is guaranteed to never create variables.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 611\u001b[0;31m       \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_stateless_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# pylint: disable=not-callable\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    612\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_stateful_fn\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    613\u001b[0m       \u001b[0;31m# Release the lock early so that multiple threads can perform the call\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/tensorflow/python/eager/function.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m   2418\u001b[0m     \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_lock\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2419\u001b[0m       \u001b[0mgraph_function\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_maybe_define_function\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2420\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mgraph_function\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_filtered_call\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# pylint: disable=protected-access\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2421\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2422\u001b[0m   \u001b[0;34m@\u001b[0m\u001b[0mproperty\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/tensorflow/python/eager/function.py\u001b[0m in \u001b[0;36m_filtered_call\u001b[0;34m(self, args, kwargs)\u001b[0m\n\u001b[1;32m   1659\u001b[0m       \u001b[0;31m`\u001b[0m\u001b[0margs\u001b[0m\u001b[0;31m`\u001b[0m \u001b[0;32mand\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m`\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;31m`\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1660\u001b[0m     \"\"\"\n\u001b[0;32m-> 1661\u001b[0;31m     return self._call_flat(\n\u001b[0m\u001b[1;32m   1662\u001b[0m         (t for t in nest.flatten((args, kwargs), expand_composites=True)\n\u001b[1;32m   1663\u001b[0m          if isinstance(t, (ops.Tensor,\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/tensorflow/python/eager/function.py\u001b[0m in \u001b[0;36m_call_flat\u001b[0;34m(self, args, captured_inputs, cancellation_manager)\u001b[0m\n\u001b[1;32m   1743\u001b[0m         and executing_eagerly):\n\u001b[1;32m   1744\u001b[0m       \u001b[0;31m# No tape is watching; skip to running the function.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1745\u001b[0;31m       return self._build_call_outputs(self._inference_function.call(\n\u001b[0m\u001b[1;32m   1746\u001b[0m           ctx, args, cancellation_manager=cancellation_manager))\n\u001b[1;32m   1747\u001b[0m     forward_backward = self._select_forward_and_backward_functions(\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/tensorflow/python/eager/function.py\u001b[0m in \u001b[0;36mcall\u001b[0;34m(self, ctx, args, cancellation_manager)\u001b[0m\n\u001b[1;32m    591\u001b[0m       \u001b[0;32mwith\u001b[0m \u001b[0m_InterpolateFunctionError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    592\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcancellation_manager\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 593\u001b[0;31m           outputs = execute.execute(\n\u001b[0m\u001b[1;32m    594\u001b[0m               \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msignature\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    595\u001b[0m               \u001b[0mnum_outputs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_num_outputs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/tensorflow/python/eager/execute.py\u001b[0m in \u001b[0;36mquick_execute\u001b[0;34m(op_name, num_outputs, inputs, attrs, ctx, name)\u001b[0m\n\u001b[1;32m     57\u001b[0m   \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     58\u001b[0m     \u001b[0mctx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mensure_initialized\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 59\u001b[0;31m     tensors = pywrap_tfe.TFE_Py_Execute(ctx._handle, device_name, op_name,\n\u001b[0m\u001b[1;32m     60\u001b[0m                                         inputs, attrs, num_outputs)\n\u001b[1;32m     61\u001b[0m   \u001b[0;32mexcept\u001b[0m \u001b[0mcore\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_NotOkStatusException\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# build training data and shuffe!\n",
    "#uTrain = np.arange(0.1, 1, 0.000001)\n",
    "#random.shuffle(uTrain)\n",
    "\n",
    "# create training data, consiting of only one moment u = 1\n",
    "#uTrain = np.ones(1000000)\n",
    "print(uTrain)\n",
    "#uTrain2 = np.ones(1000000)*2\n",
    "#uTrain = np.concatenate(uTrain1,uTrain2)\n",
    "\n",
    "#Create Early Stopping callback\n",
    "es = EarlyStopping(monitor='loss', mode='min', min_delta =0.00005, patience=50, verbose=10) # loss == custom_loss1dMBPrime by model definition\n",
    "mc = ModelCheckpoint('best_model.h5', monitor='loss', mode='min', save_best_only=True)\n",
    "\n",
    "# use uTrain as labels, since the custom loss is prepared that way \n",
    "#,tf.keras.callbacks.History() ,\n",
    "\n",
    "history = model.fit(uTrain, uTrain,validation_split=0.3, epochs=1500, batch_size=1000, verbose=1, callbacks = [es,mc])\n",
868
869
870
871
872
873
874
875
    "\n",
    "# summarize history for loss\n",
    "plt.plot(history.history['loss'])\n",
    "plt.plot(history.history['val_loss'])\n",
    "plt.title('model loss')\n",
    "plt.ylabel('loss')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'test'], loc='upper left')\n",
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
    "plt.show()\n",
    "\n",
    "plt.plot(history.history['val_loss'])\n",
    "plt.title('model loss')\n",
    "plt.ylabel('loss')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['test'], loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-4.494337]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYKklEQVR4nO3de5QcZ3nn8e8zV2l0sZA0voAtywLbxEAAM4CJAa8NGJNl4bBZEm4BvLtRCDfDcuDEm5yz/LF7soEcSAIsHMdgAjgQApglxGAggI0DNoyNwXdjDL5hyWP5ostopqe7n/2ja6SxdjQz1kx1T5e+n3PG091VXe9bLul3Xj1V9VZkJpKk6unpdAckSeUw4CWpogx4SaooA16SKsqAl6SKMuAlqaIMeEmqKANeOkQRMRARX4qIX0dERsS/63SfpJkMeGlxrgTeAGzrdEekAxnwOiwVI+4nzXj/6Yj4n49lG5lZy8y/zswrgcaSd1JaJANeOkBEbIqIh+f4eV2n+ygtRF+nOyAtN5l5F7Cu0/2QFssRvCRVlAGvw9U4MDTj/dHTL4oSze45fl7f/u5Kj50lGh2urgNeFxE3Ai8BzgBGYV+JZvVCNhIRg0AUbwciYgUwmc7DrWXAEbwOV+cB/wF4GHg98NVD3M6twF7gCcBlxevjF989afHCgYYkVZMjeEmqKANekirKgJekijLgJamiltVlkhs3bszNmzd3uhuS1DWuueaaBzJzeLZlyyrgN2/ezOjoaKe7IUldIyLuPNgySzSSVFEGvCRVlAEvSRVlwEtSRRnwklRRBrwkVZQBL0kVZcBLUgd9+6btfOLyX5aybQNekjrou7ds55NX/qqUbRvwktRBtXoy0FtOFBvwktRBU40mA30GvCRVzlSjSX9vzL/iITDgJamDWgHvCF6SKqfWSANekqpoqt70JKskVdFUo0l/nzV4Saoca/CSVFHW4CWpoqYa1uAlqZJqda+Dl6RKsgYvSRXVuorGgJekyql163XwEbEuIr4UEbdExM0R8bwy25OkbjPVyNImG+srZav7/Q3wzcz8TxExAAyV3J4kdZUyJxsrLeAj4gjghcCbATKzBtTKak+Suk2zmdSb3Xkd/AnAGHBRRPw0Ii6MiFUHrhQRWyNiNCJGx8bGSuyOJC0vU80mQFcGfB9wKvDxzHwmsAf40wNXyswLMnMkM0eGh4dL7I4kLS9TjQToypOs9wD3ZObVxfsv0Qp8SRKtmSSB7rvRKTO3AXdHxMnFRy8CbiqrPUnqNlONIuC79CqadwAXF1fQ3AGcW3J7ktQ1ao1ya/ClBnxmXgeMlNmGJHWrbq7BS5LmMFXyCN6Al6QOqXXrSVZJ0txqJZ9kNeAlqUOmL5O0Bi9JFTN9ktUavCRVzP6TrNbgJalSpmvwZU0XbMBLUodMj+CtwUtSxXgdvCRV1FS9OMlqiUaSqqXmSVZJqiZr8JJUUdbgJamivNFJkirKycYkqaKmGk36e4MIA16SKqVWb5ZWngEDXpI6pjWCN+AlqXJqjTTgJamKphpNBko6wQoGvCR1zFSjWdo0BQB9pW0ZiIhfA7uABlDPzJEy25OkblKrN0u7ixVKDvjCmZn5QBvakaSuMllvMthvDV6SKqfsEXzZAZ/AtyLimojYOtsKEbE1IkYjYnRsbKzk7kjS8jFZbzDY11va9ssO+Odn5qnAy4C3RcQLD1whMy/IzJHMHBkeHi65O5K0fNTqzdIe1wclB3xm3lv8vh+4BHhOme1JUjeZrDcZ7MaAj4hVEbFm+jVwNnBDWe1JUrcpewRf5lU0RwGXFJPo9AH/kJnfLLE9SeoqrRF8eTX40gI+M+8Anl7W9iWp2012cw1eknRwtXqjO2vwkqS5eaOTJFVQZlJrNBns4hudJEmzmGokmTDY3703OkmSZlFrtJ7H2s1TFUiSZjE51QCwBi9JVeMIXpIqanKqFfCO4CWpYvaP4D3JKkmVsm8E741OklQttUbrJKtTFUhSxTiCl6SKmpyuwRvwklQt+0fwnmSVpEqpOYKXpGradyerAS9J1TI9gjfgJalirMFLUkVZg5ekipoewRvwklQxtUaDvp6gtydKa6P0gI+I3oj4aUR8vey2JKlbTE41Sx29Q3tG8OcBN7ehHUnqGhP1BitKfFwflBzwEXEs8O+BC8tsR5K6zcRUkxVdPoL/a+B9QPNgK0TE1ogYjYjRsbGxkrsjScvDxFSDFQNdOoKPiJcD92fmNXOtl5kXZOZIZo4MDw+X1R1JWlZaI/guDXjgdOAVEfFr4AvAWRHxuRLbk6SuMVlvsKLEx/VBiQGfmedn5rGZuRl4DfDdzHxDWe1JUjeZmOryk6ySpNntbUPA95W69UJmfh/4fjvakqRuMDHV7N4SjSTp4CamGl19klWSdBATU00GrcFLUvVMTjVYacBLUvW0TrJag5ekSqk3mtSbuTyuoomIE4G/AE4BVkx/nplbSuqXJFXWRL01e8tyGcFfBHwcqANnAp8BvCtVkg7BRPHA7eVyo9PKzPxXIDLzzsx8P61ZIiVJj9G+gC/5MsmF3ug0GRE9wC8i4u3AvcDq8rolSdU1Mf3A7WVSojkPGALeCTwL+EPgTWV1SpKqrF0lmgWN4DPzJ8XL3cC55XVHkqpvOuDLvg5+oVfRnAS8Fzh+5ncy86yS+iVJlTVdolkWI3jgn4BPAH8HNMrrjiRV3/4STbk1+IUGfD0zP15qTyTpMDFRXwY1+IhYX7z854h4K3AJMDm9PDMfLLFvklRJ+0o0Hb5M8hoggSjev3fGsgS8k1WSHqO9y6FEk5knlNq6JB2GJouAL3u64IVeRbMCeCvwfFoj9x8An8jMiRL7JkmVtLe2jC6TpDX3zC7gI8X71wGfBV5dRqckqcrGpxr09wYDfcvjKpqnZuYpM95/LyJuKqNDklR1e2vlP+wDFj5VwbURcdr0m4h4LjBaTpckqdrGa3WGBhY6vj50C23hWcAPI+Ku4v0m4NaIuB7IzPztA79Q1O2vAAaLdr6Umf9jCfosSV1tT63B0ED5I/iFBvw5h7DtSeCszNwdEf3AlRHxjcy86hC2JUmVsbfWYGiwwwE/40anXbMtn+tGp8xMWpOTAfQXP3kIfZSkShmv1Rnq73yJ5sAbnaYDOljAjU4R0Vts40nAxzLz6lnW2QpsBdi0adOCOy5J3WpvrcG6oYHS21nwjU7FaP5EZjyTdT6Z2QCeERHrgEsi4qmZecMB61wAXAAwMjLiCF9S5Y3XGjx+3TKpwUfEf6X10I9jgeuA04AfAi9ayPcz8+GI+B6tWv4N860vSVU2Xmuwsg0nWR/LE52eDdyZmWcCzwQemesLETFcjNyJiJXAS4BbDr2rklQN47U6q5bRZZITmTkREUTEYGbeEhEnz/OdY4C/L+rwPcAXM/Pri+qtJFXA+DK7TPKeYjT+VeDbEfEQcOdcX8jMn9Ma6UuSCo1mMllvtqVEs9Bnsr6qePn+opZ+BPDN0nolSRU1PVXwchrB75OZl5fREUk6HIzX6gBtmaqg3KnMJEmPMj7ZvhG8AS9JbTReM+AlqZL2TrVKNCst0UhStTiCl6SKMuAlqaL27gt4SzSSVCm7Jls1+FVtmA/egJekNtpTBPyawf7S2zLgJamNdk/U6e0JVvSXH78GvCS10e7JOqsH+4iI+VdeJANektpo10Qr4NvBgJekNto9OWXAS1IV7ZlssHqFAS9JlbNr0hKNJFXS7okpR/CSVEW7J+uscQQvSdWze6LOKgNekqql2Uz21BrW4CWpavYUj+tb0+01+Ig4LiK+FxE3RcSNEXFeWW1JUjfYXcxD064RfJmt1IH3ZOa1EbEGuCYivp2ZN5XYpiQtW7snpmeS7PIRfGbel5nXFq93ATcDTyirPUla7qanCq7UZZIRsRl4JnD1LMu2RsRoRIyOjY21ozuS1BHTI/jKXCYZEauBLwPvysydBy7PzAsycyQzR4aHh8vujiR1zCN7pwA4YmX5c8FDyQEfEf20wv3izPxKmW1J0nJXmYCP1mTHnwRuzswPldWOJHWL6YBf2+0BD5wO/CFwVkRcV/z8bontSdKytnPvFAN9PazoL/95rFDiZZKZeSVQ/iNLJKlLPLJ3qm3lGfBOVklqGwNekirKgJekijLgJamiDHhJqigDXpIqqNFMdk3U23YNPBjwktQWuybaexcrGPCS1BbtnqYADHhJaouHxlsB/7ghA16SKuXBPZMArF810LY2DXhJaoMdu2sAbFg12LY2DXhJaoMH97QCfv1qR/CSVCkP7qkx0NfDqoH2zCQJBrwktcUDu2tsXDVA61EZ7WHAS1IbPLhnsq3lGTDgJaktHtxTY30bT7CCAS9JbbFjT40NbbxEEgx4SWqL1gjegJekSpmYajBea7DBGrwkVcvYrtZdrButwUtStWzbOQHAUUesaGu7pQV8RHwqIu6PiBvKakOSusH2IuCPXluRgAc+DZxT4vYlqStse6QYwa+tSIkmM68AHixr+5LULbbvnGCwr6etc8HDMqjBR8TWiBiNiNGxsbFOd0eSltz2nZMctXZFW6cpgGUQ8Jl5QWaOZObI8PBwp7sjSUtu286JttffYRkEvCRV3f07JziyzfV3MOAlqVSZybadExxVpRF8RHwe+BFwckTcExH/pay2JGm52rGnxsRUk2Mft7LtbfeVteHMfG1Z25akbnHXg+MAHPe4oba3bYlGkkp0dxHwmzYY8JJUKfc8tBegIyUaA16SSnTXjnE2rh5gaKC0ivhBGfCSVKK7HxrnuPXtL8+AAS9JpbrrwfGOnGAFA16SSjNeq3Pvw3t54vDqjrRvwEtSSX55/x4y4aSjDHhJqpTbtu8C4MQOBXz7T+uW4Pc+/kMm642DLg/mnsFtvgne5p3/bZ4NzPf9xbY/3wx1cy2dv+3FdW7x+77IY7eI7S924r/FHJeFtF/mn4v525/727090NfbQ39PtH73Bn09PfT1Bv29PfT2xL5lfb3BQG8Pqwb7GBroZdVAH0ODrd+rBnsZGuhj3VB/R65CWazb7t9Ff29w/IZVHWm/+/6PzeLINYPU6s1Zl+U8382ce435vz/P8kW2P5/52z/4CvN+dxHbXtj351kn9/2ntPYP/t3O/rmYbwOLb//Q/98tZN+azWSq2aTeSOqNJlPN1u96Y8bnzcf2Z39lfy8bVg+wYdUAG1YPMrx6kGMft5JNG4Y4bv0Qm9YPsWHVQNun5J3LL7bvZsvG1fT3dqZYUomA//gbntXpLkh6jDJbId9oJpNTTcan6uyZbDBem/G71mDPZJ2Hx6fYsXuSHXtq7NhTY/vOCa6/95F9D7Oetm6on1OOWctTHr+WUx6/lmdvXs+xHbqCBeDG3zzCaVs2dKz9SgS8pO4TEfT3Bv29sKK/lyN47E872ltrcM9D49z90Dh37hjntu27uPE3O/n7H92571/1m9YP8bwtGzjj5GHOPPlIVg70LvWuzOq+R/ayfeckzzhuXVvam40BL6lrrRzo5cSj1nDiUWse9Xm90eS27bu5+lc7+NEvd/CNG+7jH0fvZmV/L2c+eZhXPP0JvPi3jqSvxNLJdXc9DGDAS9JS6uvt4ZSiTHPu6SfQaCY//tWDXHr9fXzjhm1cev02jl67gtc+ZxOvP20TG1cv/cM4rrv7YQaKfnRKLPYk31IaGRnJ0dHRTndDUoXVG02+e8v9fO7qu7jitjFW9PfwxudtZusLtyxp0L/q//wbAXzlracv2TZnExHXZObIbMscwUs6rPT19nD2U47m7KcczR1ju/no927nwh/cwWd/dCdvOeOJ/PEZW1jRv7g6/cPjNX5298O8/awTl6jXh8YbnSQdtrYMr+ZDv/8MvvPfzuCsJx/Jh79zGy/+0OVcduO2RW33ytsfoJlwxkkbl6inh8aAl3TY2zK8mo+9/lT+4Y+ey9BAL3/82Wt428XXsmP35PxfnsU3rt/G+lUDPP3YdUvb0cfIgJekwu88cSP/8s4X8N6Xnsy3b9rO2R++gm/ecN9j2sauiSm+c/N2Xv7bx5R6lc5CGPCSNEN/bw9vO/NJ/PM7ns8x61bwls9dyzs//1Me2lNb0Pf/8Sd3M1lv8nunHltyT+dnwEvSLE4+eg2XvPV03v3ik7j0+vt4yYev4Fvz1ObHa3Uu/MGveO4J63l6B69/n1ZqwEfEORFxa0TcHhF/WmZbkrTU+nt7OO/FJ/K1tz+fI9cMsvWz1/CuL/yUh8dnjOZ3bYOLXga7tvPBy25l284J3nP2yZ3r9AylBXxE9AIfA14GnAK8NiJOKas9SSrLKY9fy1ffdjrnvehEvv7z+3jBB77H3/7rL7h/1wRc/gHyrqv42cXnc9G//Zo3/85mnnPC+k53GSjxRqeIeB7w/sx8afH+fIDM/IuDfccbnSQtd7ds28lfXXYbH73jHFbE1P+3PPsGiT+/v239metGpzJLNE8A7p7x/p7is0eJiK0RMRoRo2NjYyV2R5IW78lHr+XCN43wmzddxa3D51CL1t2vjd4V5NNeTZx3fYd7uF/HT7Jm5gWZOZKZI8PDw53ujiQtyJYtT+LkTY9ngCnoW0Fvs0YMroU1R3W6a/uUOVXBvcBxM94fW3wmSdWw53541rkwci6MXgS7t3e6R49SZsD/BDgxIk6gFeyvAV5XYnuS1F6vuXj/65d/qHP9OIjSAj4z6xHxduAyoBf4VGbeWFZ7kqRHK3U2ycy8FLi0zDYkSbPr+ElWSVI5DHhJqigDXpIqyoCXpIpaVs9kjYgx4M5D/PpG4IEl7M5y5/5W3+G2z+7voTk+M2e9S3RZBfxiRMToweZjqCL3t/oOt312f5eeJRpJqigDXpIqqkoBf0GnO9Bm7m/1HW777P4uscrU4CVJj1alEbwkaQYDXpIqqusCPiJeHRE3RkQzIkYOWHZ+8YDvWyPipTM+r8zDvyPiGRFxVURcVzwJ6znF5xERf1vs488j4tRO93WpRMQ7IuKW4rh/YMbnsx7vKoiI90RERsTG4n2Vj+8Hi+P784i4JCLWzVhWyWPctkzKzK76AX4LOBn4PjAy4/NTgJ8Bg8AJwC9pTVPcW7zeAgwU65zS6f1YxP5/C3hZ8fp3ge/PeP0NIIDTgKs73dcl2t8zge8Ag8X7I+c63p3u7xLt83G0ptm+E9hY5eNb7NvZQF/x+i+Bv6zyMW5nJnXdCD4zb87MW2dZ9ErgC5k5mZm/Am4HnlP83J6Zd2RmDfhCsW63SmBt8foI4DfF61cCn8mWq4B1EXFMJzq4xP4E+N+ZOQmQmdNPMz7Y8a6CDwPvo3Wsp1X1+JKZ38rMevH2KlpPf4PqHuO2ZVLXBfwcDvaQ7wU9/LuLvAv4YETcDfwVcH7xedX2c9pJwAsi4uqIuDwinl18Xsn9jYhXAvdm5s8OWFTJ/Z3Ff6b1LxWo7j63bb9KfeDHoYqI7wBHz7LozzLz/7a7P+021/4DLwLenZlfjojfBz4JvLid/Vtq8+xvH7CeVlni2cAXI2JLG7u35ObZ3/9Oq2RRKQv5Ox0RfwbUgYtnWU+HYFkGfGYeSmDN9ZDvrnr491z7HxGfAc4r3v4TcGHxumsfcj7P/v4J8JVsFS9/HBFNWpM0VW5/I+JptGrNP4sIaO3TtcWJ9K7dX5j/73REvBl4OfCi4lhDl+/zHNq2X1Uq0XwNeE1EDBYP+j4R+DEzHv4dEQO0Hv79tQ72c7F+A5xRvD4L+EXx+mvAG4urLU4DHsnM+zrRwSX2VVonWomIk2idlHqAgx/vrpWZ12fmkZm5OTM30/qn+6mZuY3qHl8i4hxa5xxekZnjMxZV7hgX2pZJy3IEP5eIeBXwEWAY+JeIuC4zX5qZN0bEF4GbaP0z722Z2Si+U6WHf/8R8DcR0QdMAFuLzy+ldaXF7cA4cG5nurfkPgV8KiJuAGrAm4oR3kGPd0VV9fgCfJTWlTLfLv7lclVmvmWuv9PdLDPr7cokpyqQpIqqUolGkjSDAS9JFWXAS1JFGfCSVFEGvCRVlAGvw15E7O50H6QyGPCSVFEGvFQo7hL9YETcEBHXR8QfFJ8fExFXFHPw3xARL4iI3oj49Ix1393p/ksH6ro7WaUS/UfgGcDTac1385OIuAJ4HXBZZv6viOgFhor1npCZTwWY+ZAKablwBC/t93zg85nZyMztwOW0ZrD8CXBuRLwfeFpm7gLuALZExEeKuVR2dqrT0sEY8NI8MvMK4IW0Zvz7dES8MTMfojXS/z7wFvbP6iktGwa8tN8PgD8o6uvDtEL9xxFxPLA9M/+OVpCfWjwrtSczvwz8OVCZZ6SqOqzBS/tdAjyP1jMyE3hfZm6LiDcB742IKWA38EZaT+C5KCKmB0nnz7ZBqZOcTVKSKsoSjSRVlAEvSRVlwEtSRRnwklRRBrwkVZQBL0kVZcBLUkX9P+GmvKxN8MuZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-3.0447404]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZSklEQVR4nO3deZBdZZ3/8fe393RnI0knhITQYEIEwo+tBWQYNfEHAurgUiqgozI6lDpx0KJk3KpGq5yaRWcct5JfRhYRHJdRHJ2fC7gQDEugAwTZwhIIATTpkITe7/qdP+65nSZ0ujudfs699+nPq6qrz7333PN8T07lmyff85znMXdHRETiU1fpAEREJAwleBGRSCnBi4hESgleRCRSSvAiIpFSghcRiZQSvIhIpJTgRSbAzJrM7L/M7GkzczN73Tj732pmQ2bWl/xsSSdSkX2U4EUmbgPwHuBPE9x/rbvPTH5WBoxLZFRK8BK9pMe9fMTr68zsCwdzDHfPuvu/u/sGoDDlQYoEoAQv05qZLTOzvWP8XHIIh/9HM9tlZrePV9IRCaGh0gGIVJK7PwPMDXDovwMeBrLARcDPzOxkd38yQFsio1IPXiQAd9/o7r3unnH3bwO3AxdUOi6ZXpTgZToYAFpHvD68vJGUaPrG+Hn3FMXggE3RsUQmRCUamQ7uBy4xs4eAc4DXAl0wXKKZOZGDmFkz+5J0k5m1ABnfb85tM5sLnAGsB/LAu4DXAJcf6omIHAz14GU6uBx4M7AXeDfwk0keZwswCCwBfpVsHwVgZp82s18k+zUCXwC6gV3AR4G3uPtjk2xXZFJMC36IiMRJPXgRkUgpwYuIREoJXkQkUkrwIiKRqqphkgsWLPCOjo5KhyEiUjM2bdq0y93bR/usqhJ8R0cHXV1dlQ5DRKRmmNm2A32mEo2ISKSU4EVEIhW0RGNmTwO9lObPzrt7Z8j2RERknzRq8KvdfVcK7YiIyAgq0YiIRCp0gnfgZjPbZGaXjbaDmV1mZl1m1tXd3R04HBGR6SN0gj/b3U8Fzgf+xsxes/8O7r7O3TvdvbO9fdShnCIiMglBE7y7P5f83gncBJwesj0RkVpzy8M7uGp9mJUcgyV4M2szs1nlbeBc4MFQ7YmI1KLfPrqDqzc8FeTYIUfRLAJuMrNyO991918GbE9EpObkCk5jXZjVHIMleHffCpwU6vgiIjHIF4o01IcppmiYpIhIBeWKTkN9mB68EryISAXl8kWa1IMXEYlPXj14EZE45QpFGurUgxcRiU6+4DSqBy8iEp98UT14EZEo5QqqwYuIRClfLNKoUTQiIvHJF5yGQE+yKsGLiFRQrqAevIhIlHIaRSMiEifNRSMiEqlcUT14EZEo5fUkq4hInPIaBy8iEqecxsGLiMRJ4+BFRCLk7sl0werBi4hEJV90gGBrsirBi4hUSK5QBKCxQT14EZGo5AqlHrxq8CIikcmXe/CqwYuIxKVcg9c4eBGRyAzX4PUkq4hIXPIF9eBFRKKUL5Z68BoHLyISmfIoGo2DFxGJzL4SjXrwIiJRyQ4Pk1QPXkQkKhoHLyISqeFx8LVagzezejO7z8z+J3RbIiK1pDwOvpZr8JcDj6TQjohITSnfZK3JGryZLQXeCHwrZDsiIrVoeBx8jT7J+u/AlUDxQDuY2WVm1mVmXd3d3YHDERGpHrla7cGb2ZuAne6+aaz93H2du3e6e2d7e3uocEREqk4t1+D/DPgLM3sa+B6wxsxuCNieiEhNKSf4plpb8MPdP+XuS929A7gI+K27vydUeyIitSablGiaarAHLyIiY8jlkx58oATfEOSo+3H3W4Fb02hLRKRW7FuTtcZusoqIyNiyeU1VICISpeFRNLU6VYGIiIwuW3CaGuowU4IXEYlKrlAMdoMVlOBFRComVygGe4oVlOBFRCqmlODVgxcRiU4mrwQvIhKlXHKTNRQleBGRCsnldZNVRCRKuUIx2FOsoAQvIlIxWd1kFRGJU1Y3WUVE4qQHnUREIqVRNCIikdKTrCIikdJNVhGRSGU1Dl5EJE6ai0ZEJFK6ySoiEqmcxsGLiMQpq6kKRETi4+5k9aCTiEh8CkXHHZVoRERikys4gG6yiojEJlsoAurBi4hEJ5ck+CZNVSAiEpdsXj14EZEo5VSiERGJ03CJRjdZRUTiks2XRtGoBy8iEpl9PfgavMlqZi1mdreZbTazh8zs86HaEhGpNWkMk2wIdmTIAGvcvc/MGoENZvYLd78rYJsiIjUhl8IommAJ3t0d6EteNiY/Hqo9EZFakkkSfHOt3mQ1s3ozux/YCdzi7htH2ecyM+sys67u7u6Q4YiIVI1MvgBAc0N9sDaCJnh3L7j7ycBS4HQzWzXKPuvcvdPdO9vb20OGIyJSNYZ78I012oMvc/e9wO+A89JoT0Sk2tV0icbM2s1sbrI9AzgHeDRUeyIitaSc4EM+6BRyFM1i4NtmVk/pH5IfuPv/BGxPRKRmZHLha/AhR9E8AJwS6vgiIrWsPA6+Jks0IiJyYJmcEryISJQy+dJ6rGY1OFWBiIgcWCZfCNp7ByV4EZGKyOSLQcfAgxK8iEhFZPPFoCNoQAleRKQiMvli0DHwMMFhkma2AvhH4Higpfy+ux8TKC4RkahlctVTg78W+CaQB1YD1wM3hApKRCR2mXyxahL8DHf/DWDuvs3dPwe8MVxYIiJxS6MGP9EnWTNmVgc8bmZrgeeAmeHCEhGJWyZfoLUp5GwxE+/BXw60An8LnAb8JfC+UEGJiMQujRLNhP75cPd7ks0+4NJw4YiITA9pjIOf6CiaY4FPAEeN/I67rwkUl4hI1KqpBv9D4CrgP4BCuHBERKaHTL5AU8AFt2HiCT7v7t8MGomIyDRS8RKNmc1LNn9mZh8BbgIy5c/dfXfA2EREopXJVf4m6ybAgfJ8lp8Y8ZkDepJVRGQSsoUK1+Dd/eigrYuITEP5QpFC0SvegwfAzFqAjwBnU+q5/x64yt2HAsYmIhKlNBbchonfZL0e6AW+lry+BPgO8I4QQYmIxKyc4KuiBw+scvfjR7z+nZk9HCIgEZHYZcsJvrE65oO/18zOLL8wszOArjAhiYjELZMvPU5ULT3404A7zOyZ5PUyYIuZ/QFwd/8/QaITEYnQUK66avDnBY1CRGQaGcqVevCtTRUcJjniQafe0T7Xg04iIgdvMEnwLRWei2b/B508+W3oQScRkUkZTvCV7MGPfNAp6c2vYMSarCIicvAySYKfEXgUzUQfdPogpUU/lgL3A2cCdwCvDxaZiEikBlNK8AezotOrgG3uvho4BXgxWFQiIhEbzJZG0bRUSYIfKk9LYGbN7v4osDJcWCIi8UqrBz/RYZLPmtlc4CfALWa2B9gWKigRkZgNDd9krYJx8O7+1mTzc2b2O2AO8MuxvmNmR1Kaw2YRpRE369z9K4cQq4hIFIZyBeqMqlnRaZi7r5/grnngCne/18xmAZvM7BZ31xw2IjKtDWYLtDTWY2bj73wIgv3z4e5/dPd7k+1e4BFgSaj2RERqxVC+ELz+DgET/Ehm1kFp5M3GUT67zMy6zKyru7s7jXBERCpqMFsMPoIGUkjwZjYT+BHwMXfv2f9zd1/n7p3u3tne3h46HBGRihvKFZgR+ClWCJzgzayRUnK/0d1/HLItEZFaMZgr0NIYvoASrAUr3T24GnjE3f8tVDsiIrVmKFf7Nfg/A/4SWGNm9yc/FwRsT0SkJpR68OET/EEPk5wod9/AvlkoRUQkMZgt0D6zOXg7qYyiERGRfTL5SEbRiIjISw1ma78GLyIioxiMYZikiIi83FBKN1mV4EVEUlQselKDr+Fx8CIi8nJD+XTmggcleBGRVA1kkwSvGryISFwGMqUE39YU7DGkYUrwIiIp6svkAWhrVg9eRCQqA9lyglcPXkQkKuUefKtKNCIicSnfZFWJRkQkMv3lGrx68CIicRlO8KrBi4jEpT8p0bRqHLyISFwGsnka6ozmBk1VICISlf5MgdamekqrmoalBC8ikqL+TD6V+jsowYuIpGogW1CCFxGJUV8mT1sKN1hBCV5EJFUD2XwqT7GCEryISKr6MoVUnmIFJXgRkVQNZHWTVUQkSqVhkkrwIiLR6ddNVhGR+OQKRQZzBWa1NKbSnhK8iEhK+oZKE43NnqESjYhIVHqGcgDMTqkHn84/I4Ht6BnC/cCfjzflw7gzQoz7/bF3ONT2x5uzYvzvj9f+GDuM991aPrcJfH88h9J+1f/ZpTBXynTTM1juwSvBT9jrvngrg7lCpcMQmVYa642Gujoa6o2GOqOhvo7G5HdDvdFY/ix5f0ZTPW1NDbQ27/s9s6mB1uYGZjU3MH9mE/NnNjO/rYkFM5uZkdKNyDTt68Gnk3qjSPCfv/AECsXRu/Bj9ewBnLF3GP/74xjnAON9f9z2Ax5//O+Oe/aTbrvUfthrE7r9sb9b2+fm7uSKTr5QJFdw8sUi+YKTL7+X/M4XSvvl8kV6h/Ls6BmiP1NgIJunP1sgmy8esI3WpnoWz2nhqPltLJvXyrJ5rRzd3saqI+bQPqt57ACrVM9gkuDVg5+4d3YeWekQRGQScoUiA9kCPYM5dvdneaE/w66+LC/0ZdnVl+G5PYNs2z3Axq0vDC+UAbBodjOrjpjDaR2HcfbyBZxwxBzq66q/pFTuwc+q9R68mV0DvAnY6e6rQrUjIrWrsb6OOTPqmDOjkSPntR5wP3dnd3+Wx3f28dDzPTz43Is88OxefvPoTv6FLcyZ0chrj23nzScdwWuOXUBzQ3WWd3qH4qnBXwd8Hbg+YBsiMg2YWak+P7OZM4+ZP/z+zt4h7nzyBTY8votfP7KDn25+ntktDbz9tKW8/6wOjprfVsGoX65nMIcZzEzpSdZgrbj7bWbWEer4IiILZ7Vw4clLuPDkJeQKRTY8sYub7n2O79y5jevueJpzjlvEFeeuZOXhsyodKgA9Q3lmNTdQl1I5qeLj4M3sMjPrMrOu7u7uSocjIjWqsb6O1SsX8tWLT+H2T65h7erl3Ln1Bc77ym1c8YPNdPdmKh0iPYO51MozUAUJ3t3XuXunu3e2t7dXOhwRicCi2S1cce5KbvvEaj549tH8bPPznPPl9dx037OHPPrrUPQM5VJ7yAmqIMGLiIRyWFsTn3nj8fz88rM5ZkEbH//+ZtZ+9z76MvmKxNMzlE9tmgJQgheRaWD5wln88ENn8XfnvZJfPPhH3vKN29na3Zd6HD2DudQmGoOACd7M/hO4E1hpZs+a2QdCtSUiMp76OuPDr3sFN3zgDPb0Z3n7N+/ggWf3phpDz2AkJRp3v9jdF7t7o7svdferQ7UlIjJRZy1fwI8/chYzWxq4eN1d3LX1hdTa3j2QZf7MptTaU4lGRKado+a38V8fOosj5s7gA9fdw+bte4O3OZgtMJQrclirEryISFCLZrdwwwfPYN7MJt5/7d08sbM3aHsv9JeGac5ri6BEIyJS7RbNbuE7f3UG9XXGB77dxYsDuWBt7ekvHVs9eBGRlHQsaOOq95zG83sHufz79x1wZtpDtXsgC8C8NiV4EZHUdHbM4+/ffAK3bunm6799Ikgbe/qV4EVEKuLdZyzjracs4au/fZz7A9x03a0ELyJSGWbG5y88gcNnt/Dx79/PQHZqn3bd3Z+lztJbjxWU4EVEhs1uaeRL7ziJp1/o519vfmxKj717IMthrU2pzSQJSvAiIi/x6lfM55LTl3HdHU/z0PMvTtlxu3szqT7kBErwIiIvc+UbXsncGY189icPUpyiUTU7ezMsmt0yJceaKCV4EZH9zGlt5NMXHMd9z+zlh5u2T8kxd/YMsXCWEryISMW97dQlnLpsLv9682OHfMO1WHS6ezMsmt08RdFNjBK8iMgozIxPX3AcO3szXP37pw7pWHsGsuSLzsJZSvAiIlWhs2MebzhhEVetf/KQlvzb0VP6rmrwIiJV5MrzXslQvsg3fjf5J1x39g4BsFAlGhGR6vGK9pm8/dQlfPfuZ9jRMzSpY5S/p5usIiJVZu3qFRSKzlXrn5zU95/dM0h9nbF4jhK8iEhVWTa/lbeesoTvbnyGnZPoxW/fPcDiOS001KebcpXgRUQmYO3q5eSLzlXrtx70d7fvGWTpYTMCRDU2JXgRkQnoWNDGW05ewo0btw3fNJ2o7bsHOPKw1kCRHZgSvIjIBK1ds5xcoci3DmJc/FCuwM7eDEfOU4IXEalaRy9o480nHcENd20bnt99PNt3DwCwTAleRKS6rV29nMFcgas3TKwW/9iOPgCWL5wZMqxRKcGLiByEFYtmccGqxXz7jm0TWqR7y45e6kwJXkSkJqxds5y+TJ5rbh+/Fv/4jl6Omt9GS2N9CpG9lBK8iMhBOm7xbM45fhHX3v4UvUNj9+K37OhlRQV676AELyIyKX+7ZgU9Q3muv3PbAfd5cSDH1u5+TlwyJ8XI9lGCFxGZhBOXzmH1yna+9fut9GdGny/+3mf2AHBax2FphjZMCV5EZJI++voV7BnIcePG0XvxXdt2U19nnHzk3HQDSyjBi4hM0qnLDuPs5Qv4f+u38uLgy2vx6x/r5qSlc2htaqhAdErwIiKH5JPnv5I9A1m+9KstL3l/++4BHnyuh3NPOLxCkQVO8GZ2npltMbMnzOyTIdsSEamEVUvm8N5Xd3DDxm3c/dRu6P0TXHs+/73hPuoM3nji4orFFizBm1k98A3gfOB44GIzOz5UeyIilXLFucfSMb+Nj9x4L3t/+QV8253M3/Rlzlt1eEXmoCkLWRg6HXjC3bcCmNn3gAuBhwO2KSKSulktjfxm8J3UFTLwUOm9i+0WLn78FvhCM3x2Z0XiClmiWQJsH/H62eS9lzCzy8ysy8y6uru7A4YjIhJO3cceYHDl28jVldZdLTbMgBPfAZf/oXIxVazlhLuvc/dOd+9sb2+vdDgiIpMz63BmzJxLo+egoaXUm2+eDbMWVSykkCWa54AjR7xemrwnIhKn/p1w2qXQeSl0XQt9OyoaTsgEfw+wwsyOppTYLwIuCdieiEhlXXTjvu03/Vvl4kgES/DunjeztcCvgHrgGnd/KFR7IiLyUkEfr3L3nwM/D9mGiIiMruI3WUVEJAwleBGRSCnBi4hESgleRCRS5u6VjmGYmXUDB14eZWwLgF1TGE610/nGb7qds853co5y91GfEq2qBH8ozKzL3TsrHUdadL7xm27nrPOdeirRiIhESgleRCRSMSX4dZUOIGU63/hNt3PW+U6xaGrwIiLyUjH14EVEZAQleBGRSNVcgjezd5jZQ2ZWNLPO/T77VLLA9xYze8OI96NZ/NvMTjazu8zs/mQlrNOT983Mvpqc4wNmdmqlY50qZvZRM3s0ue7/MuL9Ua93DMzsCjNzM1uQvI75+n4xub4PmNlNZjZ3xGdRXuPUcpK719QPcBywErgV6Bzx/vHAZqAZOBp4ktI0xfXJ9jFAU7LP8ZU+j0M4/5uB85PtC4BbR2z/AjDgTGBjpWOdovNdDfwaaE5eLxzrelc63ik65yMpTbO9DVgQ8/VNzu1coCHZ/mfgn2O+xmnmpJrrwbv7I+6+ZZSPLgS+5+4Zd38KeILSwt/Di3+7exYoL/5dqxyYnWzPAZ5Pti8ErveSu4C5Zra4EgFOsQ8D/+TuGQB3L69efKDrHYMvA1dSutZlsV5f3P1md88nL++itPobxHuNU8tJNZfgx3CgRb4ntPh3DfkY8EUz2w58CfhU8n5s51l2LPDnZrbRzNab2auS96M8XzO7EHjO3Tfv91GU5zuKv6L0PxWI95xTO6+gC35Mlpn9Gjh8lI8+4+7/nXY8aRvr/IHXAx939x+Z2TuBq4H/m2Z8U22c820A5lEqS7wK+IGZHZNieFNunPP9NKWSRVQm8nfazD4D5IEbR9lPJqEqE7y7TyZhjbXId00t/j3W+ZvZ9cDlycsfAt9Ktmt2kfNxzvfDwI+9VLy828yKlCZpiu58zexESrXmzWYGpXO6N7mRXrPnC+P/nTaz9wNvAl6fXGuo8XMeQ2rnFVOJ5qfARWbWnCz0vQK4mxGLf5tZE6XFv39awTgP1fPAa5PtNcDjyfZPgfcmoy3OBF509z9WIsAp9hNKN1oxs2Mp3ZTaxYGvd81y9z+4+0J373D3Dkr/dT/V3f9EvNcXMzuP0j2Hv3D3gREfRXeNE6nlpKrswY/FzN4KfA1oB/6/md3v7m9w94fM7AfAw5T+m/c37l5IvhPT4t9/DXzFzBqAIeCy5P2fUxpp8QQwAFxamfCm3DXANWb2IJAF3pf08A54vSMV6/UF+DqlkTK3JP9zucvdPzTW3+la5u75tHKSpioQEYlUTCUaEREZQQleRCRSSvAiIpFSghcRiZQSvIhIpJTgZdozs75KxyASghK8iEiklOBFEslTol80swfN7A9m9q7k/cVmdlsyB/+DZvbnZlZvZteN2PfjlY5fZH819ySrSEBvA04GTqI03809ZnYbcAnwK3f/BzOrB1qT/Za4+yqAkYtUiFQL9eBF9jkb+E93L7j7DmA9pRks7wEuNbPPASe6ey+wFTjGzL6WzKXSU6mgRQ5ECV5kHO5+G/AaSjP+XWdm73X3PZR6+rcCH2LfrJ4iVUMJXmSf3wPvSurr7ZSS+t1mdhSww93/g1IiPzVZK7XO3X8EfBaIZo1UiYdq8CL73AS8mtIamQ5c6e5/MrP3AZ8wsxzQB7yX0go815pZuZP0qdEOKFJJmk1SRCRSKtGIiERKCV5EJFJK8CIikVKCFxGJlBK8iEiklOBFRCKlBC8iEqn/BTg8TFepw67RAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-6.951536]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXCElEQVR4nO3de5SkdX3n8fe3u6e7p+cKzOgwCAwaRAVFcVCMMVkka9DVsLnshkWjYS9s4poQT6InrrurfyQnMXqyazRHDvHCMWF1vYckJmKygok5IgNB7q6gchkG6GGYoXu6u6qr6rt/VNXQzJlLz0w/VV2/eb/O6UPVU089v+/DM/OZX/+eyy8yE0lSeYb6XYAkqRoGvCQVyoCXpEIZ8JJUKANekgplwEtSoQx4SSqUAS8tUkRcFBH3RsRMRHwjIk4/xLo/iojZiJju/Fzfy1olMOClRYmIDcCXgP8OnAhsA/7PYb72psxc3fl5XdU1Svsz4FW8iMiI+LEF76+JiN89ws38PHBXZn4+M+eA9wPnRsQLlrBUaUkZ8DquRcRpEbH7ED+XdVY9G/hu93uZuRe4v7P8YK6NiMmIuD4izq1wN6QDGul3AVI/ZeaDwPpFrLoamNxv2R5gzUHWfzNwKxDAlcDXIuIFmbn76CqVjpw9eGlxpoG1+y1bC0wdaOXM/FZmzmbmTGb+PrAbeE21JUrPZMDreDADTCx4v6n7ojNEM32Inzd3Vr0LOHfB91YBz+ssX4yk3ZuXesaA1/HgNuCyiBiOiIuBn+p+kJkPLrjS5UA/13ZW/TJwTkT8QkSMA/8DuD0z792/sc4/Gq+OiNGIGI+IdwEbgG9VvqfSAga8jgdXAm+iPUzyZuArR7qBzJwEfgH4PeBJ4JXApd3PI+KqiLiq83YN8LHOetuBi4HXZ+YTR70H0lEIJ/yQpDLZg5ekQhnwklQoA16SCmXAS1KhltWdrBs2bMgtW7b0uwxJGhi33HLLzszceKDPllXAb9myhW3btvW7DEkaGBHxwME+c4hGkgplwEtSoQx4SSqUAS9JhTLgJalQBrwkFcqAl6RCGfCS1Edfv/sxrrrx/kq2bcBLUh/933sf45P/+MNKtm3AS1If1RotVgxXE8UGvCT1Ub3RYmzEgJek4sw37cFLUpHqjRaj9uAlqTzzzWTFcFSybQNekvrIHrwkFaruGLwklcmraCSpUPPNAR2iiYh3RsRdEXFnRHwmIsarbE+SBs1ADtFExCnAbwBbM/McYBi4tKr2JGkQzTdajA5awHeMACsjYgSYAB6puD1JGij1ZosVgzZEk5nbgQ8BDwI7gD2Zef3+60XEFRGxLSK2TU5OVlWOJC1L9UHswUfECcAlwBnAZmBVRLxl//Uy8+rM3JqZWzdu3FhVOZK0LNUH9CTrTwM/zMzJzJwHvgT8eIXtSdLAGcgePO2hmQsiYiIiArgIuKfC9iRpoDRbSSsZvKtoMvMm4AvArcAdnbaurqo9SRo09UYLoLIhmpFKttqRme8D3ldlG5I0qOrNdsD7sDFJKky3B++jCiSpMPP7evAGvCQVpeoxeANekvqk24M34CWpMLWGQzSSVKTuVTSeZJWkwtTmHaKRpCLVGk0AxkaGK9m+AS9JfeJ18JJUqJoBL0llejrgHaKRpKLsG6JZYQ9ekory9ElWA16SilLzUQWSVKZ918F7J6sklaXebDIyFIwY8JJUltp8dRNugwEvSX1Ta7QqO8EKBrwk9U290arsGngw4CWpb2qNZmXXwIMBL0l9U2u0KruCBgx4SeqbWqNlD16SSuQYvCQVqtZoOkQjSSVyiEaSClX3OnhJKlPNMXhJKlNtvumjCiSpRD6qQJIK5WWSklSoWsOnSUpScVqtpN50iEaSilNvVjvhNhjwktQX3flYHYOXpMLUGk2gugm3wYCXpL7oTrjtGLwkFWbfGPygBnxErI+IL0TEvRFxT0S8qsr2JGlQ9KIHP1LZlts+DPxtZv5iRIwCExW3J0kDoTsGX+VJ1soCPiLWAT8J/ApAZtaBelXtSdIg6V5FM6gnWc8AJoFPRcQ/R8THI2LV/itFxBURsS0itk1OTlZYjiQtH3Pz7R78+IrBvExyBDgP+FhmvgzYC/zO/itl5tWZuTUzt27cuLHCciRp+ZjrjMGPD+iNTg8DD2fmTZ33X6Ad+JJ03Ov24FcOYg8+Mx8FHoqIszqLLgLurqo9SRokvRiiqfoqml8Hru1cQfMD4PKK25OkgTDbgx58pQGfmbcBW6tsQ5IGUXcMfuXoAA7RSJIOrtuDH9g7WSVJB1abbzK+YoiIqKwNA16S+mB2vlnpCVYw4CWpL+bmm5WeYAUDXpL6Yna+ZQ9ekko0W3eIRpKKVGs0WVnhYwrAgJekvrAHL0mFmmt4klWSimQPXpIKNedVNJJUprnOnaxVMuAlqQ+80UmSCpSZzM43K32SJBjwktRz882kldVO9gEGvCT13GwPZnMCA16Seq62L+A9ySpJRenFdH1gwEtSz3Wn63OIRpIKYw9ekgo1W+/Mx+oYvCSVZa5hD16SitTtwXujkyQVZqYT8KtGRyptZ1Fbj4gzgd8HXgSMd5dn5nMrqkuSijVTbwAwsUx68J8CPgY0gAuBTwN/XlVRklSyvbVOD36s2h78YgN+ZWb+PRCZ+UBmvh/4V9WVJUnlmqk3GAoYG6l2lHyx/3zUImII+H5EvAPYDqyurixJKtdMvcnE6AgRUWk7i/3n40pgAvgN4OXALwNvq6ooSSrZTL1R+fg7LLIHn5k3d15OA5dXV44klW9vrVn5+Dss/iqa5wPvAk5f+J3MfG1FdUlSsZZVDx74PHAV8KdAs7pyJKl8e2vNyq+Bh8UHfCMzP1ZpJZJ0nJiZb7J+5YrK2znkSdaIODEiTgT+MiLeHhEnd5d1lkuSjtBMrcGqsf4P0dwCJNC9luddCz5LwDtZJekIdS+TrNohW8jMMyqvQJKOM3vrDVYtl5OsETEOvB34Cdo9938ArsrMuQprk6QizdSaTPTgMsnF3uj0aeBs4CPARzuv/2wxX4yI4Yj454j4q6MrUZLKMd9sUW+2mKj4WfCw+KtozsnMFy14/42IuHuR370SuAdYe0SVSVKBuo8KXk49+Fsj4oLum4h4JbDtcF+KiOfQfijZx4+uPEkqS/dRwctmDJ7282f+KSIe7Lw/DfheRNwBZGa+5CDf+1/Au4E1B9twRFwBXAFw2mmnLbIcSRpM3UcF96IHv9gWLj7SDUfEG4HHM/OWiPgXB1svM68GrgbYunVrHmk7kjRIlk0PfsHNTFMH+jwzdx3i668GfjYi3kB7Fqi1EfHnmfmWo6pUkgrQ7cFXPR8rHPmNTt0ednCYG50y8z3AewA6PfjfNtwlHe9m57s9+GV0o1OnN38mC+ZklSQdmaen6+t/Dx6AiPiPtC93fA5wG3AB8E/ARYv5fmbeANxwNAVKUkmma50e/DK6TPJK4Hzggcy8EHgZsKeyqiSpUNNz7YBfM97np0kuMNd9LEFEjGXmvcBZ1ZUlSWWampsngmV1J+vDEbEe+Arw9Yh4EnigqqIkqVRPzTVYPTbC0FC1E27D4udk/bnOy/dHxDeAdcDfVlaVJBVqutZgTQ/G32HxPfh9MvPGKgqRpOPB1Nx8T8bfYfFj8JKkJTA112DNeG968Aa8JPXQdK3BagNeksrT7sE7RCNJxWmPwduDl6TiTM317ioaA16SeqTeaFFrtOzBS1Jpus+hcQxekgozNTcPwGqHaCSpLFP7HjRmwEtSUboB73XwklSY7hDNWsfgJaksDtFIUqE8ySpJhdo92w74dSsdopGkouyeaT+mYGS4N9FrwEtSj+yZnWf9RG9672DAS1LP7J6pc8LEaM/aM+AlqUeenJnv2fg7GPCS1DPtIRp78JJUnN0zddbbg5eksrRayZ7ZeU7wJKsklWWq1qCVsM4hGkkqy+6ZOoBDNJJUmt0z7btYvQ5ekgrTfUyBV9FIUmH2DdHYg5eksuzp9uAdg5eksuza2+7BeyerJBVm53SNEyZW9OxJkmDAS1JP7Jyqs2H1WE/bNOAlqQee2Fsz4CWpRDun65y0uneXSEKFAR8Rp0bENyLi7oi4KyKurKotSVrudk71vgdf5cyvDeC3MvPWiFgD3BIRX8/MuytsU5KWnbn5JlO1BhvXFDJEk5k7MvPWzusp4B7glKrak6Tl6onOJZInrSpkiGahiNgCvAy4qRftSdJysnOqBlDeSdaIWA18EfjNzHzqAJ9fERHbImLb5ORk1eVIUs/tnO4EfClDNAARsYJ2uF+bmV860DqZeXVmbs3MrRs3bqyyHEnqiyemCxuiiYgAPgHck5l/VFU7krTcTXZ68MWcZAVeDfwy8NqIuK3z84YK25OkZWlyqsaasRHGVwz3tN3KLpPMzH8EoqrtS9Kg2LFnlk3rxnverneySlLFduyZ4+T1K3vergEvSRXbsWeOk9fag5ekotQbLXZO1zh5vQEvSUV57Kk5MuFkx+AlqSyPPjUHwKZ1jsFLUlEe2T0LwGZ78JJUlkf3dHvwBrwkFWXHnjlWj42wZrx3k213GfCSVKEHd81w6okTfWnbgJekCv3oib1sOcmAl6SiNFvJw7tmOc2Al6Sy7NgzS73ZYstJq/rSvgEvSRV58IkZAE53DF6SyvKjbsBvsAcvSUV5YNdeRoeH2NSHB42BAS9Jlfnh5F5OO2mC4aH+TI1hwEtSRf7fY1Oc9ew1fWvfgJekCszUGzywa4bnG/CSVJb7Hp8mE87atLpvNRjwklSB7z06BWAPXpJKc++jU4yNDHF6n25yAgNekipx+8O7OXvz2r5dQQMGvCQtuUazxR3b9/CS56zvax0GvCQtse8/Ps3cfIuXnrq+r3UY8JK0xG57aDcA5xrwklSWm37wBBtWj/btOfBdBrwkLaHM5Fv3P8GrnreBiP6dYAUDXpKW1P2T00xO1fjx553U71IMeElaSjd8bxKAVz9vQ58rMeAlaUldf9djvGDTmr5N07eQAS9JS2RyqsbND+zi4nM29bsUwICXpCXzF7dtJxPe8OKT+10KYMBL0pLITD5780O89NT1fX3A2EIGvCQtgW/d9wT3PT7Nv3vFqf0uZR8DXpKOUWbyx3//fTatHedfv+yUfpezjwEvScfob+58lO/8aBdvv/B5jI0M97ucfQx4SToGO6drvO+6uzjnlLVc9orT+l3OM4z0u4Cl8B+uuZl6s3VM2zjcLcWLueH4cHclH24bx1rD4u6KPkwbx7wPi6ngGGs4bI2LKuJYPi7iz8tiVjrmY3WY5ociGBkORoaCkeEhRoaDFUNDz1w2FKxY8Nn46DCrRodZOTrMqtERJkaHmRgbYdXoMBOjI4yO9K7fOlNv8J//7Bam5ua55vLzGRleXn3mSgM+Ii4GPgwMAx/PzD+oop2ZepNao3nQz/Mw38/DrHC47y9mI8dew2G2v4gij3U/czGNHHMNx7afi6nwcPtx2G0sgxqO9f/jorZxjId7MX9empk0mkmjlTSaLeY7/20dQ9sTo8OcMDHKiatGOWHVKCetGuWEiVE2rhlj8/pxTlm/ks3rV/KsNWOLC+SpR+ELl8MvXgNrnr1v8UO7ZnjH/76VO7bv4aOXncfZm9cdfdEVqSzgI2IY+BPgXwIPAzdHxHWZefdSt/WZKy5Y6k1K6qNWqxP6rRbzzXboN1rJfLPF3HyTmXqTvbUmM/UGM/Wn/7u31mD3zDy79tbZNVPnyb11frhzml3TdfbWn9kJHB4KNq0dZ/P6cTZ3Qv+U7s8J7ferx0bgxj+EB78NN36A+dd/iNsf3s1ffncHn735QYYjuOotL+d1Zy+PG5v2V2UP/hXAfZn5A4CI+CxwCbDkAS+pLENDwehQMLqEpwmnaw127J5l++5ZHtk9xyO7Z3lkzyzbn5zl1gef5K9v30Fjwa8O9469DWL+6Q1s+wQrtn2Cs3MFlzY+zZvO3cxvv+4sNq9fuWQ1LrUqA/4U4KEF7x8GXrn/ShFxBXAFwGmnLa8TFJLKsXpshDOfvYYzD3ITUrOVTE7VOv8AzPK5x/6Kc+/5EC/c/U1Gs0Y9xrj/pAvZ/sr38p0XvYATVo32eA+OXN9Psmbm1cDVAFu3bj32QV5JOgrDQ8GmdeNsWjfOy08/AdgMtVPh1nkYGWe0WeeFW07hhee/pN+lLlqVAb8dWHhL13M6yyRpMOx9HF5+OWy9HLZ9CqYf63dFR6TKgL8ZODMizqAd7JcCl1XYniQtrUuvffr1G/+of3UcpcoCPjMbEfEO4Gu0L5P8ZGbeVVV7kqRnqnQMPjO/Cny1yjYkSQe2vG67kiQtGQNekgplwEtSoQx4SSpULMUDpJZKREwCDxzl1zcAO5ewnOXO/S3f8bbP7u/ROT0zNx7og2UV8MciIrZl5tZ+19Er7m/5jrd9dn+XnkM0klQoA16SClVSwF/d7wJ6zP0t3/G2z+7vEitmDF6S9Ewl9eAlSQsY8JJUqIEL+Ij4NxFxV0S0ImLrfp+9JyLui4jvRcTPLFh+cWfZfRHxO72veulExEsj4tsRcVtEbIuIV3SWR0T8cWcfb4+I8/pd61KJiF+PiHs7x/0PFyw/4PEuQUT8VkRkRGzovC/5+H6wc3xvj4gvR8T6BZ8VeYx7lkmZOVA/wAuBs4AbgK0Llr8I+C4wBpwB3E/7McXDndfPBUY767yo3/txDPt/PfD6zus3ADcseP03QAAXADf1u9Yl2t8Lgb8Dxjrvn3Wo493vepdon0+l/ZjtB4ANJR/fzr69DhjpvP4A8IGSj3EvM2ngevCZeU9mfu8AH10CfDYza5n5Q+A+2hN/75v8OzPrQHfy70GVwNrO63XAI53XlwCfzrZvA+sj4uR+FLjEfg34g8ysAWTm453lBzveJfifwLtpH+uuUo8vmXl9ZjY6b79Ne/Y3KPcY9yyTBi7gD+FAk3yfcojlg+o3gQ9GxEPAh4D3dJaXtp9dzwdeExE3RcSNEXF+Z3mR+xsRlwDbM/O7+31U5P4ewL+n/ZsKlLvPPduvvk+6fSAR8XfApgN89N7M/Ite19Nrh9p/4CLgnZn5xYj4t8AngJ/uZX1L7TD7OwKcSHtY4nzgcxHx3B6Wt+QOs7//lfaQRVEW83c6It4LNIBrD7CejsKyDPjMPJrAOtQk3wM1+feh9j8iPg1c2Xn7eeDjndcDO8n5Yfb314AvZXvw8jsR0aL9kKbi9jciXkx7rPm7EQHtfbq1cyJ9YPcXDv93OiJ+BXgjcFHnWMOA7/Mh9Gy/ShqiuQ64NCLGOhN9nwl8hwWTf0fEKO3Jv6/rY53H6hHgpzqvXwt8v/P6OuCtnastLgD2ZOaOfhS4xL5C+0QrEfF82ieldnLw4z2wMvOOzHxWZm7JzC20f3U/LzMfpdzjS0RcTPucw89m5syCj4o7xh09y6Rl2YM/lIj4OeAjwEbgryPitsz8mcy8KyI+B9xN+9e8/5KZzc53Spr8+z8BH46IEWAOuKKz/Ku0r7S4D5gBLu9PeUvuk8AnI+JOoA68rdPDO+jxLlSpxxfgo7SvlPl65zeXb2fmrx7q7/Qgy8xGrzLJRxVIUqFKGqKRJC1gwEtSoQx4SSqUAS9JhTLgJalQBryOexEx3e8apCoY8JJUKANe6ujcJfrBiLgzIu6IiF/qLD85Ir7ZeQb/nRHxmogYjohrFqz7zn7XL+1v4O5klSr088BLgXNpP+/m5oj4JnAZ8LXM/L2IGAYmOuudkpnnACycpEJaLuzBS0/7CeAzmdnMzMeAG2k/wfJm4PKIeD/w4sycAn4APDciPtJ5lspT/SpaOhgDXjqMzPwm8JO0n/h3TUS8NTOfpN3TvwH4VZ5+qqe0bBjw0tP+Afilzvj6Rtqh/p2IOB14LDP/lHaQn9eZK3UoM78I/DegmDlSVQ7H4KWnfRl4Fe05MhN4d2Y+GhFvA94VEfPANPBW2jPwfCoiup2k9xxog1I/+TRJSSqUQzSSVCgDXpIKZcBLUqEMeEkqlAEvSYUy4CWpUAa8JBXq/wO6rFMxCQOgVAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make prediction\n",
    "pred = model.predict([1])\n",
    "print(pred)\n",
    "# plot prediction to minimal entropy\n",
    "u = 1\n",
    "alpha = np.arange(-100.0, 1.0, 0.05)\n",
    "plt.figure()\n",
    "plt.plot(alpha,custom_loss1dMBPrime(u,alpha))\n",
    "plt.plot(pred,custom_loss1dMBPrime(u,pred), marker=\"*\")\n",
    "plt.ylabel('alpha')\n",
    "plt.xlabel('loss')\n",
    "plt.title('u=1')\n",
    "plt.show()\n",
    "# plot prediction to minimal entropy\n",
    "u = 1.5\n",
    "pred = model.predict([u])\n",
    "print(pred)\n",
    "alpha = np.arange(-100.0, 1.0, 0.05)\n",
    "plt.figure()\n",
    "plt.plot(alpha,custom_loss1dMBPrime(u,alpha))\n",
    "plt.plot(pred,custom_loss1dMBPrime(u,pred), marker=\"*\")\n",
    "plt.ylabel('alpha')\n",
    "plt.xlabel('loss')\n",
    "plt.title('u=1.5')\n",
    "plt.show()\n",
    "\n",
    "u = 0.5\n",
    "pred = model.predict([u])\n",
    "print(pred)\n",
    "alpha = np.arange(-100.0, 1.0, 0.05)\n",
    "plt.figure()\n",
    "plt.plot(alpha,custom_loss1dMBPrime(u,alpha))\n",
    "plt.plot(pred,custom_loss1dMBPrime(u,pred), marker=\"*\")\n",
    "plt.ylabel('alpha')\n",
    "plt.xlabel('loss')\n",
    "plt.title('u=0.5')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Test the model ####\n",
    "\n",
    "Use alphas computed with a Newton solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
For faster browsing, not all history is shown. View entire blame