Let us first explore an example that falls under novelty detection. Here, we train a model on data with some distribution and no outliers. The test data, has some "novel" subset of data that does not follow that distribution.
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm
%matplotlib inline
Use the np.random module to generate a normal distribution of 1,000 data points in two dimensions (e.g. x, y) - choose whatever mean and sigma^2 you like. Generate another 1,000 data points with a normal distribution in two dimensions that are well separated from the first set. You now have two "clusters". Concatenate them so you have 2,000 data points in two dimensions. Plot the points. This will be the training set.
In [ ]:
mu,sigma=3,0.1
x=np.random.normal(mu,sigma,1000)
y=np.random.normal(mu,sigma,1000)
x_0=np.random.normal(2,sigma,1000)
y_0=np.random.normal(2,sigma,1000)
X_train_normal=np.ndarray(shape=(2000,2))
for i in range(0,2000):
if (i<1000):
X_train_normal[i]=[x[i],y[i]]
else:
X_train_normal[i]=[x_0[i-1001],y_0[i-1001]]
print(xy)
In [ ]:
p_x=X_train_normal[:,0]
p_y=X_train_normal[:,1]
print(p_y)
plt.scatter(p_x,p_y)
plt.show()
Generate 100 data points with the same distribution as your first random normal 2-d set, and 100 data points with the same distribution as your second random normal 2-d set. This will be the test set labeled X_test_normal.
In [ ]:
X_test_normal=np.concatenate((0.1*np.random.randn(100,2)+2,0.1*np.random.randn(100,2)+3))
plt.scatter(X_test_normal[:,0],X_test_normal[:,1])
Generate 100 data points with a random uniform distribution. This will be the test set labeled X_test_uniform.
In [ ]:
X_test_uniform=np.random.rand(100,2)+3
plt.scatter(X_test_uniform[:,0],X_test_uniform[:,1])
Define a model classifier with the svm.OneClassSVM
In [ ]:
model = svm.OneClassSVM()
Fit the model to the training data.
In [ ]:
model.fit(X_train_normal)
Use the trained model to predict whether X_test_normal data point are in the same distributions. Calculate the fraction of "false" predictions.
In [ ]:
predicted=model.predict(X_test_normal)-1
print(np.count_nonzero(predicted))
Use the trained model to predict whether X_test_uniform is in the same distribution. Calculate the fraction of "false" predictions.
In [ ]:
uniform=model.predict(X_test_uniform)-1
print(np.count_nonzero(uniform))
Use the trained model to see how well it recovers the training data. (Predict on the training data, and calculate the fraction of "false" predictions.)
In [ ]:
trained=model.predict(X_train_normal)-1
print(np.count_nonzero(trained))
Create another instance of the model classifier, but change the kwarg value for nu. Hint: Use help to figure out what the kwargs are.
In [ ]:
new_model=svm.OneClassSVM(nu=0.1)
new_model.fit(X_train_normal)
Redo the prediction on the training set, prediction on X_test_random, and prediction on X_test.
In [ ]:
new_predicted=new_model.predict(X_test_normal)-1
new_uniform=new_model.predict(X_test_uniform)-1
new_trained=new_model.predict(X_train_normal)-1
print(np.count_nonzero(new_trained))
print(np.count_nonzero(new_predicted))
print(np.count_nonzero(new_uniform))
Plot in scatter points the X_train in blue, X_test_normal in red, and X_test_uniform in black. Overplot the trained model decision function boundary for the first instance of the model classifier.
In [193]:
plt.scatter(X_train_normal[:,0],X_train_normal[:,1],color='blue')
plt.scatter(X_test_uniform[:,0],X_test_uniform[:,1],color='black')
plt.scatter(X_test_normal[:,0],X_test_normal[:,1],color='red')
xx1, yy1 = np.meshgrid(np.linspace(1.5, 4, 1000), np.linspace(1.5, 4,1000))
Z1 =model.decision_function(np.c_[xx1.ravel(), yy1.ravel()])
Z1 = Z1.reshape(xx1.shape)
plt.contour(xx1, yy1, Z1, levels=[0],
linewidths=2)
Out[193]:
Do the same for the second instance of the model classifier.
In [194]:
plt.scatter(X_train_normal[:,0],X_train_normal[:,1],color='blue')
plt.scatter(X_test_uniform[:,0],X_test_uniform[:,1],color='black')
plt.scatter(X_test_normal[:,0],X_test_normal[:,1],color='red')
xx1, yy1 = np.meshgrid(np.linspace(1.5, 4, 1000), np.linspace(1.5, 4,1000))
Z1 =new_model.decision_function(np.c_[xx1.ravel(), yy1.ravel()])
Z1 = Z1.reshape(xx1.shape)
plt.contour(xx1, yy1, Z1, levels=[0],
linewidths=2)
Out[194]:
In [ ]:
from sklearn.covariance import EllipticEnvelope
Test how well EllipticEnvelope predicts the outliers when you concatenate the training data with the X_test_uniform data.
In [ ]:
train_uniform=np.concatenate((X_train_normal,X_test_uniform))
envelope=EllipticEnvelope()
envelope.fit(train_uniform)
envelope.predict(train_uniform)
Compute and plot the mahanalobis distances of X_test, X_train_normal, X_train_uniform
In [195]:
print(range(100))
In [197]:
plt.scatter(range(100),envelope.mahalanobis(X_test_uniform),color='black')
plt.scatter(range(2000),envelope.mahalanobis(X_train_normal),color='blue')
plt.scatter(range(200),envelope.mahalanobis(X_test_normal),color='red')
plt.show()
In [ ]:
data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAAEp0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMC4wcmMyKzI5MjAuZzExNWJhZGUsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy9bT2XBAAAgAElEQVR4nO3df5AkZ33f8fd353aB2cOWNHdRXSR2RsTEKcWVgNhScIGphAMMFwcpiaNwLGHBlLdYxS5RVMoIb1XKSdW6AFewZWKOWiHIWTMGORhKKkIAmeDE/gOZPRA/hJB16HZOUkm6HxKg4wBJu9/80T27s7PdMz0z3fOj5/Oq6tqZ3p6eZ57u/s4z3366H3N3RERk/E0NuwAiIpIOBXQRkZxQQBcRyQkFdBGRnFBAFxHJCQV0EZGcUEAXEckJBXQRkZxQQBcRyYl9SRYys0uAjwG/BDjwG8ADwB1ABdgAbnD3p9qt58CBA16pVHovrYjIBDpx4sQ5dz/YaTlLcum/mR0H/trdP2ZmM0AR+F3gSXd/v5ndDFzq7u9tt575+XlfX19P9glERAQAMzvh7vOdluuYcjGznwdeDdwG4O7PuPsPgOuA4+Fix4Hrey+uiIj0K0kO/SrgLPAJM/uGmX3MzGaBy939sXCZx4HLo15sZktmtm5m62fPnk2n1CIiskeSgL4PuAY45u4vA34M3Ny8gAd5m8jcjbuvufu8u88fPNgxBSQiIj1KEtAfAR5x93vC558mCPBPmNkhgPDvmWyKKCIiSXQM6O7+OPCwmf1iOOsw8F3gLmAxnLcI3JlJCUVEJJGk/dB/G6iZ2beAlwK/D7wfeJ2ZPQi8Nnw+MLVajUqlwtTUFJVKhVqtNsi3FxEZOYn6obv7vUBUl5nD6RYnmVqtxtLSEhcvXgSgXq+ztLQEwMLCwjCKJCIydGN5pejKysp2MG+4ePEiKysrQyqRiMjwjWVAP336dFfzRUQmwVgG9Lm5ua7mi4hMgrEM6KurqxSLxV3zisUiq6urQyqRiMjwjWVAX1hYYG1tjXK5jJlRLpdZW1vTCVERmWiJbs6VFt2cS0Ske6ndnEtERMaDArqISE4ooIuI5IQCuohITiigi4jkhAK6iEhOKKCLiOSEArqISE4ooIuI5IQCuohITiigi4jkhAK6iEhOKKCLiOSEArqISE4ooIuI5IQCuohITiigi4jkhAK6iEhOKKCLiOTEviQLmdkG8DSwCTzn7vNmdhlwB1ABNoAb3P2pbIopIiKddNNC/xfu/tKmgUpvBr7s7i8Bvhw+FxGRIekn5XIdcDx8fBy4vv/iiIhIr5IGdAe+ZGYnzGwpnHe5uz8WPn4cuDzqhWa2ZGbrZrZ+9uzZPosrIiJxEuXQgVe5+6Nm9veAu83se83/dHc3M496obuvAWsA8/PzkcuIiEj/ErXQ3f3R8O8Z4LPAtcATZnYIIPx7JqtCiohIZx0DupnNmtkLG4+B1wPfAe4CFsPFFoE7syqkiIh0liTlcjnwWTNrLP9n7v4FM/sa8Odm9k6gDtyQXTFFRKSTjgHd3R8C/mnE/PPA4SwKJSIi3dOVoiIiOTG2Ab1Wq1GpVJiamqJSqVCr1YZdJBGRoUrabXGk1Go1lpaWuHjxIgD1ep2lpaB7/MLCwjCLJiIyNGPZQl9ZWdkO5g0XL15kZWVlSCUSERm+sQzop0+f7mq+iMgkGMuAPjc319V8EZFJMJYBfXV1lWKxuGtesVhkdXV1SCUSERm+sQzoCwsLrK2tUS6XMTPK5TJra2s6ISoiE83cB3e/rPn5eV9fXx/Y+4mI5IGZnWgaiyLWWLbQRURkLwV0EZGcUEAXEckJBXQRkZxQQBcRyQkFdBGRnFBAFxHJCQV0EZGcUEAXEckJBXQRkZxQQBcRyQkFdBGRnMhlQNd4oyIyicZyTNF2NN6oiEyq3LXQNd6oiEyqxAHdzApm9g0z+1z4/Cozu8fMTprZHWY2k10x47WmV+r1euRyGm9URPKumxb6TcD9Tc8/APyhu/8C8BTwzjQLlkQjvVKv13F36vU6Zha5rMYbFZG8SxTQzexK4F8CHwufG/Aa4NPhIseB67MoYKvmFvni4uKe9Iq77wnqGm9URCZB0hb6HwG/A2yFz0vAD9z9ufD5I8AVKZdtj9YW+ebmZuRy7q7xRkVk4nTs5WJmvwaccfcTZvbPu30DM1sClqD/tEfUCc8o5XKZjY2Nvt5LRGTcJGmhvxJ4k5ltAJ8iSLXcAlxiZo0vhCuBR6Ne7O5r7j7v7vMHDx7sq7BJTmwqvSIik6pjQHf397n7le5eAd4M/B93XwC+Avx6uNgicGdmpQzFtfALhYLSKyIy8frph/5e4D1mdpIgp35bOkWKt7q6SrFY3DWvWCxy/Phxtra22NjYUDAXkYnV1ZWi7v5XwF+Fjx8Crk2/SPEawXplZYXTp08zNzfH6uqqgriICDm8UlREZFKN1b1cdJ8WEZF4Y9VC131aRETijVVAj+u2qPu0iIiMWUCP67ao+7SIiIxZQI/rtqgLiURExiygLywssLa2pvu0iIhEMHcf2JvNz8/7+vr6wN5PRCQPzOyEu893Wm6sWugiIhJPAV1EJCcU0EVEckIBXUQkJxTQRURyQgFdRCQnFNBFRHJCAV1EJCcU0EVEckIBXUQkJxTQRURyYuwDeq1Wo1KpMDU1RaVSoVarDbtIIiJDMVZD0LXSkHQiIjvGuoWuIelERHaMdUDXkHQiIjvGOqBrSDoRkR1jHdA1JJ2IyI6OAd3Mnm9mf2tm3zSz+8zsv4TzrzKze8zspJndYWYz2Rd3Nw1JJyKyo+MQdGZmwKy7XzCzaeBvgJuA9wCfcfdPmdlHgW+6+7F269IQdCIi3UttCDoPXAifToeTA68BPh3OPw5c32NZRUQkBYly6GZWMLN7gTPA3cD3gR+4+3PhIo8AV2RTRBERSSJRQHf3TXd/KXAlcC3wj5K+gZktmdm6ma2fPXu2x2LupqtDRUT26qqXi7v/APgK8MvAJWbWuNL0SuDRmNesufu8u88fPHiwr8LCztWh9Xodd9++OlRBXUQmXZJeLgfN7JLw8QuA1wH3EwT2Xw8XWwTuzKqQzXR1qIhItCQt9EPAV8zsW8DXgLvd/XPAe4H3mNlJoATcll0xd8RdBVqv19VKF5GJlqSXy7fc/WXu/k/c/Zfc/b+G8x9y92vd/Rfc/d+5+8+yL277q0CbUy/Ks4vIpBm7uy0eOXKEY8eiu7s3p150F0YRmTQdLyxKUxoXFlUqFer1euz/zYy5ubnIZcrlMhsbG329v4jIoKV2YdGo6XQnxbm5OeXZRWQijV1Ab5dDb9yYK2meXUQkT8YuoEfdYRGgVCpt35grbhlQF0cRya+xC+hRd1isVqucO3du+4RnY5k4GgBDRPJo7E6KdiPuBKpOjorIOMntSdFuaAAMEZkkuQ7oGgBDRCZJrlMuIiJ5oJRLE90GQEQmwdhd+t+txu12dRsAEcm73LfQdbtdEZkUuQ/o7W4DICKSJ7kP6HG3ATAz5dJFJFdyH9CPHDkSOd/dlXYRkVzJfUD//Oc/H/s/3QJARPIk9wG9XdBud1dGEZFxk/uA3i6HrlsAiEie5D6gR93Pxcx417vepX7oIpIruQ/oUfdzuf322/nIRz4y7KKJiKRK93IRERlxupeLiMiEUUAXEcmJsQ7ououiiMiOjgHdzF5kZl8xs++a2X1mdlM4/zIzu9vMHgz/Xpp9cXc07qJYr9dx9+27KCqoi8ik6nhS1MwOAYfc/etm9kLgBHA98HbgSXd/v5ndDFzq7u9tt640T4pqvFARmRSpnRR198fc/evh46eB+4ErgOuA4+FixwmC/MDEXQGqy/lFZFJ1lUM3swrwMuAe4HJ3fyz81+PA5amWrIO4K0CnpqaUUxeRiZQ4oJvZfuAvgHe7+4+a/+dB3iYyd2NmS2a2bmbrZ8+e7auwEOTODxw4EHs/883NTeXURWQiJQroZjZNEMxr7v6ZcPYTYX69kWc/E/Vad19z93l3nz948GBfha3VarzjHe/g/Pnze/43NbX3o2hkIhGZJEl6uRhwG3C/u3+o6V93AYvh40XgzvSLt9vKygrPPvts5P+2trYi5yunLiKTIkkL/ZXAfwBeY2b3htMR4P3A68zsQeC14fNM9RKcdYtcmVS1GlQqMDUV/FX2Mf/2dVrA3f8GsJh/H063OO3Nzc3F5s5LpRI/+clPdg0IXSwWdYtcmUi1GiwtQeNwqNeD5wC6yWh+jdWVoqurq0xPT++ZPzMzwy233MLi4iKFQgGAQqHA4uKibpErE2llZSeYN1y8GMyX/Bq7uy3WajVuuumm7ROjU1NTbG1tUSqVePrpp3nmmWe2ly0Wi6ytrSmoy8SZmoKoQ9sMYk43yQjL7d0WFxYWuOWWWyiVSsDOydDz58/vCuagXi4yueJOHemUUr6NXUBv3MMlqutilHq9rguNZOKsrkLLQF0Ui8F8ya+xC+grKyu7TnwmoQuNZNIsLMDaGpTLQZqlXA6eK/uYb2MX0PvpV64UTL6oW157CwuwsRHkzDc2FMwnwdgF9Hb9yqenpymVSgTXQkWr1+vceOONWRRNBqjRLa9eD07+NbrlKahLEnltDIxdQF9dXaXYmhwk6If+iU98gnPnzrG1tUW5XI5dx7FjxxTUx5y65fUmr4GsG50aA2NdR+4+sOnlL3+596tarXqpVGrcDMxLpZIvLy97uVx2M/NyuezVatWr1aoXi8Xt5VqnQqHQd1lkeMzcg8Nx92Q27JKNrmrVvVjcXV/FYjB/kpTL0ftOuTy6dQSse4IYO1b90Bs9XJpPik5PT2Nmkf3PAd761rfGrm+Qn13SVakELatW5XKQL5a9VGeBdn305+ZGs46S9kMfq4AeN0pRlMbIRfv27WNzc3PP/wuFAs8991zPZZHhar20HYJueerJEU8XGwXafbGdPj2adZTLC4u66eHSCPxLjRtYtIibL+NB3fK6p4uNAu366I99HSXJy6Q19ZtDL5fLsTnx1snMvBomvpaXl71QKGznzpeXl/sqRxaq1SCHZ7aTyxNJ06jmh4ch7ngb1ToiYQ59rAJ61InO6enp2KBeKpX2nCwdRaO6E0n+qOHQ2SjWUdKAPlY5dAhOjK6srHD69Gnm5uZYXV1te+Kz2ajerEsnq0SknVzm0CG4OdfGxgZbW1tsbEe7twCngM3w79HI147qlaJxpwY02JKIdGPsAnqzG2+8kbe+9X8Ba0CF4ONUgFuJC+qjOCTd2J+IEZGRMJYBvVarceDAAY4dOwb8PjDbssRsOH+vURySTnfGE5E0jF1A33v73LgAvXf+zMxMJkPS9XupsLrgiUgaxu6k6N6Li04RpFlabQBX7ZpTKpU4d+5cX+/fShe4iEymWi24d9Dp00F6dHU1u2M+tydF9+bAfxf4ccu8H4fzdzt//nzqA11keZOosb5JkHSk7Tu+RvZun0n6NqY1pXFzruiLi446nHLYDP8ebfu/YrGYWp/0rG4Spb7p+abtO97a3eArC+TxwiL36IuL4qejDhdaKv2Cw1Evp1TzpVL0hi2V+lvvoHcYGSxt396NwoU/g77bZ9KAPnYpl4WFBdbW1iiXy5hZ2/uet+sBM4rdF5upb3q+afv2ZlRSHaPa1XjsAjrsvbgoPqjH94D5rcsuSyWB+eST3c1PalR3GEmHtu9uSc8njMrAJiPb1bhTEx74OHAG+E7TvMuAu4EHw7+XJvk5kEbKJUp8GuZU5M+io/yhPzsz03cCs1p1LxSy+emsHGu+afvu6KYuekl1ZJWiGWTqh7Ry6MCrgWtaAvoHgZvDxzcDH0jyZmncnKtxs61SqbQ9clHznRQJ77RImxx63fb3HYWjdsK0D8xRyBXKbmluk1HYvqNQhm7OJ3R77qGfL85RqJuG1AJ6sC4qLQH9AeBQ+PgQ8ECS9fQT0JOeDN1798XdvVymp9/uW1F7RJuv+agNG7djFQoKvHmVt1b1MD5P1LHUTau7WnVv/XE9MxNf5l5PPneqm0EH+6wD+g+aHlvz84jXLgHrwPrc3FzPH6ibe6F3mh7uIk8St2G7/E6QHMhbz5RBf564Yymup1hUOapV9+np3ctNT8cH1F57o4zauKMDC+jh86eSrKefFvpOGqX/6Sj4hdYtFbNFymX3o1T9FGXfxPwUZT9KNbPcuYyuvA1MPejPExckS6XkAbLbL6Fev7Ta1c0wvtiTBvRee7k8YWaHAMK/Z3pcT2Jp3lTrk8BvAo8UCh1vnvLKeo1bWaJCnSmcCnVuZYkbNmujeZZbMtNrz5RhXBGa5D0H3dMmrkvmk08mv5dRN909azW4cGHv/CTHabu6Gekup0miPntb6H/A7pOiH0yynkHk0AGfmZlJtJyZdUyGPVwoR34dP1woj9RJk6xk/RnTXn+W5e3lp/aw8tRJ3nPQZUujZZt0HXGdFmZng18EnfaPdnXTTRnS2hdJsZfLJ4HHgGeBR4B3AiXgywTdFv8SuCzJm6XZy2V2djYyQC8vL3u1Wt1O0RwFPwW+Gf492rT8byf4rbdF9G+vLcb0d3YX4nbq5eV0dtSo9Td+6vay3qj1TU8nO4C7eY+4sSij5g/j53k37znIRkkaXyDt1tH8WeJSoq1T44RqVD30M+5o2l+WqQX0NKc0A3q5XPbl5eU9zxtdGRvBvDVXvhlOdTP/yf747ovVahAITlGOXObpUrnHz7B7J0krOGYhLjC05hd73VHj1t/rejutL6sWaLuDN408dTdBt1qN/+yDyPV3Kuvy8k6wLRSC592uKy74tuus0G6ane2+4dLpc6b9RZ67gB6Vcmm0whvBvLXL4qleti5B67txJv0oVb/A7q19gaK/fbqaSgsy6qBrt5MPUlwwSmtHTbL+btabtLyFQrpfoO0O3k69JToF6qStwcb7tKuDrE/YJ+nql7TV2m0LN8mXebdTPw2XtE845y6g99JtcbPHLflwobxrVlQvl14OkKQ7XSO1P2zdHCS97KhJ1t+83ubA1WjlJUlvtJvSaLG3O3jbpa3apYdKpfjufM37XtKW6SD6y3dqlWZxAVG7nHYWU9JjXi30DuK6LbbLkZ/qZYvNzPhbwoCdZhBr91M4zQ3frXatxHY57qjydpuPjepT3EvgagTCduXLsq7bBZS4tFoaQajTezdPg2ggdGqVdtNq7fQlmeQXSdL9YWoqeZ0nPeaVQ+8gqoUelSO/wE5Q/zDEXxUaM/2U6e0WeGPqpYXeHNxKpb1Xt6W143R633YnBLv5Od+c80/a4kzSA6RdvTS/vttfC5C83lvrut0XU1QOeHm5feCIqtMud8s9U6EQrKvXtFUWJ0MH0UI3S34stbtoqTFNT0fvu+0aLknrcCR7uaQ59ZtDT5ojP9Xh/5222inK20+7yaF303Lo5SDsXEedr2JtrLs5MLVOne7lHrWj9pIr7tSybc699lp/SXo+tB6kcV9Ys7PRr9+3L3mZooJHr1OnOozarnEnENPovZQkh97psv20jqFG+TudU4g70Xr4cPRrGue3or7Es0xr5TKgt/Yvj8uRb3b4/xZBnrz964OnSXu5dGqldTPF9Rnu1DsmzVxitztmpxZq1OdL8hO8194LjXV0+qXUehI6ac+eXqe01gNBubtJi0HwpRaXYkgjQHVK4bW7bL+fnirNU/MXdK+57E4NlCTvnabcBfSolMupmFo9xRUOR/1MzP/PUHKID9abmB+l6kepJrqRV6eWQC8HarMkO3parb7WHTdpiy3uAGjXKk5ysPX6JdUob2u9RAWz5sCVZh1mPTW+jJrTCo1WeFpf7kkCVKOLb2sZku4jnVIy3Uy99gdv3dfb1XncL7XG/7OQu4AedVI0OodeDHPcm36G6CjXCOhHqfpmzIVDZyjtSbU0T+f2lzv+lE/rIEq6o6ddjnYHS9LcersdP6512fyF1umLsl2ru5sAkeQn+iCmxjZsPv/RyzZK83N06v8dd2K7NbD3etK0223YKuq8UmtdR+1D7bZPuzJkIXcBPa7bYtDL5Yo9JyzBY4P1Jta0TPSWaXcydedLI72DJmoHb+hmR0+7pR530MS9z/79uw/4dmOuxgXc5i+OJGO2RuU7s66HLKfWFuQguuW128dmZuK/sDuddGzk5pN0wezncyZpGXeb0unlC6a5MaKTom1Uq1WP73O+uatSG71S4oJy80nPuLRL3Gu3INNg3ryDNyTd0VvTJN10x+pmJ09SnkarPar1Vigk66lQLsf/vG0E9Hat0ax/sWQ5zc7ublVm/X5XX905957F+7aeqO+1lR51Erz5V8QgvhQb09RU0MhQt8UOGpf17+17fut2pUX1SmmeWlvXcb1YzhB9FDV/GWQ19ZJDb7SiWtMg7Q7ExsHTTZfKxvq7OVijdvg0D552/4+6eCfr7ZfFlMWX86hNSa6ijntdp18BozL1mpLJZUCvVqv+FvbmzbfAf8j+7ZZ5XG1ugX+Y5T3/iupnHhfos26dN3bQdl0i46bWwJx2EOh0t7lRnZpbgWn0Adc0vG3YOjUaLMM+99HN5+hFLgO6u/u5Nlvup8x0vNz/h8w6xF8s1Dz/DCX/IbO+RfBlcIbSQAI6xH+TD2vHbT7Btbw8nDL0M+3b1zkvr2m8pnEJ4s2TWujNqm26EYZTkv9/gcN7TpheoOgfZnlPq7x1fYNqpTd22FJpdx55GDtxc756WK3zNH5tNFro45p20TTe0yBy6BYsOxjz8/O+vr7e8+svHDjA/vPn+y6HEwyE2uo5Cuxjs+PrNyhzFRt9l2NcVKvB36UluHix9/VMTcHWVjpl6lWhAJudN7FIKsrlYCSjublglKSoUZiSMLMT7j7fabl9va1+OIopBHOIDuYAhQTBHGCOURhrajBKpWAnrFT6C+Yw/GAOCuZRCgV4wQuih2uT/qQRzLvR65iiQ3Ge2WEXAYDTZDTo4ogpFuGGG4JgXq8PuzSShVIJrrxSwTwr7sGxs7Q0mLFkxyag33gjwPMyW39cGiZquc9xJLNyjJJKBW67TcE8z86f1/YdhIsX4W1vyz6oj01AX1uDEk9ltv4kwbyx3L/nzzMrxyj57nfhmWeGXQqRfNjagsXFbIP62AT0zU34WYYt9G4cIJ1cvohMls1NuOmm7NY/NgH9i7yW5/PTYRdDRKQvKfXtiDQ2Af11fDlxWiRrg+voKSKS3NgE9FEyKl8sIjJ+LMMAooDeo6MMoA+SiOROltdy9hXQzewNZvaAmZ00s5vTKlSUn/L8LFffFQNuIcMzGyIiPeg5oJtZAfgT4I3A1cBRM7s6rYK1eh4/y2rVPVFPFxEZNf200K8FTrr7Q+7+DPAp4Lp0irXXpFydKSL5dnVmzd7+AvoVwMNNzx8J52XiQ6VVfsp0Vqvv2jlKwy6CiIwZM7jvvuzWn/lJUTNbMrN1M1s/e/Zsz+v5Z7cssDz9Cc5S2h57bovOXQiTLBf3v7j5P2OGm7ilwztL3pVKcPhwcHOrKDMzwR0mZaeusuzhMeqmp+H22zN+kyT32I2agF8Gvtj0/H3A+9q9pv8Ri/aOF/jXy1U/R2l7EIofMuvnrORbmJ+e2hm44i3bA1cEQ9c1lt+0KT//96/eNW+LYCDpL3DY6xaMTbo5VfAtgiHofnO2GjvW5f79O+MYdjOQQqm0e6T65vuel0rBGIVp3gu9MWBFr+M4Nj5nY7s0f9bGZ+lnoOao+p2aCtbdPMReY39oXX7//r1D8rWOLxk3PF+hELw2zUF+m/fhXveLduVoLWvca6Lev3nwkqj6afwtlXbXf9yoWp3qK26fa8xL+n6dPnvrcdT6WaNee/hw5/2ik6g6HPlBogluvfsQcBUwA3wT+MftXpPGiEVZyuIAHpTmss/O7gwIUSgEO+kwPleSHbtdnY/z9pD2tG27kzSg9zXAhZkdAf4IKAAfd/fVdsv3O8CFiMgkGsgAF+7+eeDz/axDRETSoVM2IiI5oYAuIpITCugiIjmhgC4ikhMK6CIiOaGALiKSEwroIiI50deFRV2/mdlZoJ7Cqg4A51JYT9pGsVyjWCZQuboximUClasb/Zap7O4HOy000ICeFjNbT3LV1KCNYrlGsUygcnVjFMsEKlc3BlUmpVxERHJCAV1EJCfGNaCvDbsAMUaxXKNYJlC5ujGKZQKVqxsDKdNY5tBFRGSvcW2hi4hIi7EL6Gb2BjN7wMxOmtnNA3zfF5nZV8zsu2Z2n5ndFM7/PTN71MzuDacjTa95X1jOB8zsVzMs24aZfTt8//Vw3mVmdreZPRj+vTScb2b2x2G5vmVm12RQnl9sqo97zexHZvbuYdSVmX3czM6Y2Xea5nVdN2a2GC7/oJktZlSuPzCz74Xv/VkzuyScXzGznzTV20ebXvPycNufDMve1yBvMeXqeruleZzGlOmOpvJsmNm94fxB1lVcTBje/pVkFIxRmQgG0vg+8GJ2Rkm6ekDvfQi4Jnz8QuDvgKuB3wP+U8TyV4flex7BqE7fBwoZlW0DONAy74PAzeHjm4EPhI+PAP8bMOAVwD0D2GaPA+Vh1BXwauAa4Du91g1wGcHoXJcBl4aPL82gXK8H9oWPP9BUrkrzci3r+duwrBaW/Y0ZlKur7Zb2cRpVppb//zfgPw+hruJiwtD2r3FroV8LnHT3h9z9GeBTwHWDeGN3f8zdvx4+fhq4H7iizUuuAz7l7j9z91PASYLyD8p1wPHw8XHg+qb5f+qBrwKXmNmhDMtxGPi+u7e7oCyzunL3/wc8GfF+3dTNrwJ3u/uT7v4UcDfwhrTL5e5fcvfnwqdfBa5st46wbD/n7l/1IDL8adNnSa1cbcRtt1SP03ZlClvZNwCfbLeOjOoqLiYMbf8at4B+BfBw0/NHaB9UM2FmFeBlwD3hrN8Kf0J9vPHzisGW1YEvmdkJM1sK513u7o+Fjx8HLh9CuQDezO6Dbdh1Bd3XzTD2u98gaM01XGVm3zCz/2tmvxLOuyIsyyDK1c12G2R9/QrwhLs/2DRv4HXVEhOGtn+NW0AfOjPbD/wF8G53/xFwDPgHwEuBxwh+/g3aq9z9GuCNwH80s1c3/zNskQy8O5OZzQBvAv5nOGsU6lm6fWsAAAJGSURBVGqXYdVNO2a2AjwH1MJZjwFz7v4y4D3An5nZzw2wSCO33ZocZXeDYeB1FRETtg16/xq3gP4o8KKm51eG8wbCzKYJNlzN3T8D4O5PuPumu28Bt7KTKhhYWd390fDvGeCzYRmeaKRSwr9nBl0ugi+Yr7v7E2H5hl5XoW7rZmDlM7O3A78GLITBgDClcT58fIIgP/0PwzI0p2UyKVcP220g9WVm+4B/A9zRVNaB1lVUTGCI+9e4BfSvAS8xs6vC1t+bgbsG8cZhru424H53/1DT/Ob8878GGmfi7wLebGbPM7OrgJcQnJRJu1yzZvbCxmOCE2vfCd+/cbZ8EbizqVxvC8+4vwL4YdPPw7Ttaj0Nu66adFs3XwReb2aXhumG14fzUmVmbwB+B3iTu19smn/QzArh4xcT1M9DYdl+ZGavCPfPtzV9ljTL1e12G9Rx+lrge+6+nUoZZF3FxQSGuX/1c5Z3GBPBmeK/I/jmXRng+76K4KfTt4B7w+kIcDvw7XD+XcChpteshOV8gD7PqLcp14sJehF8E7ivUSdACfgy8CDwl8Bl4XwD/iQs17eB+YzKNQucB36+ad7A64rgC+Ux4FmC3OQ7e6kbgpz2yXB6R0blOkmQS23sXx8Nl/234ba9F/g68K+a1jNPEGC/D/x3wosFUy5X19stzeM0qkzh/P8BvKtl2UHWVVxMGNr+pStFRURyYtxSLiIiEkMBXUQkJxTQRURyQgFdRCQnFNBFRHJCAV1EJCcU0EVEckIBXUQkJ/4/88JKSM/h1cwAAAAASUVORK5CYII=
In [ ]:
In [ ]: