{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "

# Nonnegative Matrix Factorization (NMF)

\n", "
\n", "\n", "
\n", "\n", "

\n", "Following Section 8.3.1 [Müller, FMP, Springer 2015], we cover in this notebook a technique known as nonnegative matrix factorization (NMF). The formal definition of the NMF problem and an iterative learning procedure for computing a factorization in practice was first described by Lee and Seoung. Implementations and applications for NMF can also be found in the NMF Toolbox. \n", "\n", "

\n", "

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Basic Idea\n", "\n", "**Nonnegative matrix factorization** (NMF) is a technique where a matrix $V$ with nonnegative entries is factored into two matrices $W$ and $H$ that also have only nonnegative entries:\n", "\n", "\n", "\n", "Typically, the matrices $W$ and $H$ are required to have a much lower rank than the original matrix $V$. This factorization is interpreted as follows: The columns of $V$ are regarded as data vectors. The underlying assumption is that these vectors can be represented as a weighted superposition of a relatively small number of **template** vectors. The columns of $W$ correspond to these templates. Furthermore, the rows of $H$—called **activations**—indicate where these templates occur in $V$. The nonnegativity constraints often lead to a decomposition that allows for a semantically meaningful interpretation of the coefficients. However, in most cases, the resulting factorization problem has no exact solution, thus requiring optimization procedures for finding suitable numerical approximations. In the following, to build up some intuition, we give a factorization example form the music domain. Subsequently, we formally introduce the NMF problem and discuss an iterative learning procedure. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example: Spectrogam Factorization \n", "\n", "As illustration of NMF, we show how this technique can be used to factorize a spectrogram of a music recording into musically meaningful components. Let us consider the first measures of the Prélude Op. 28, No. 4 by Frédéric Chopin. The following figure shows the musical score as well as a piano-roll representation of the score synchronized to an audio recording of a performance. For illustration purposes, all information related to the note number $p=71$ is highlighted by the red rectangular frames.\n", "\n", "\n", "\n", "
\n", "\n", "\n", "\n", "As for the original data matrix $V$, we use the [magnitude STFT](../C2/C2_STFT-Basic.html), which is a sequence of spectral vectors. Using NMF, this matrix can be factorized into a product of two nonnegative matrices $W$ and $H$. In the ideal case, the first matrix $W$ represents the spectral patterns of the notes' pitches that occur in the piece of music, while the second matrix $H$ exhibits the time positions where these spectral patterns are active in the audio recording. The following figure shows such a factorization for our Chopin example.\n", "\n", "\n", "\n", "In this case, each template specified by the matrix $W$ reflects how a note of a certain pitch is spectrally realized in $V$, and the activation matrix $H$ looks similar to the piano-roll representation of the score. In practice, however, it is often hard to predict which of the signal's properties are ultimately captured by the learned factors. \n", "To better control this factorization, we will show how additional score information can be used to [constrain NMF](../C8/C8S3_NMFAudioDecomp.html) to yield a musically meaningful decomposition. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Formal Definition of NMF\n", "\n", "A matrix with real-valued coefficients is called **nonnegative** if all the coefficients are either zero or positive. Let $V \\in \\mathbb{R}_{\\ge 0}^{K \\times N}$ be such a nonnegative matrix having $K\\in\\mathbb{N}$ rows and $N\\in\\mathbb{N}$ columns. The dimensions $K$ and $N$ of the matrix $V$ are usually thought to be large. Given a number $R\\in\\mathbb{N}$ smaller than both $K$ and $N$, the goal of NMF is to find two nonnegative matrices $W \\in \\mathbb{R}_{\\ge 0}^{K \\times R}$ and $H \\in \\mathbb{R}_{\\ge 0}^{R \\times N}$ such that \n", "\n", "\$$\n", " V \\approx W \\cdot H.\n", "\$$\n", "\n", "As said above, the columns of $V$ are regarded as $K$-dimensional data vectors, where $N$ is the number of data vectors. This matrix is then approximately factorized into a $(K \\times R)$ matrix $W$ and an $(R \\times N)$ matrix $H$. The parameter $R$, which is referred to as the **rank** of the factorization, is usually chosen to be much smaller than $K$ and $N$. Therefore, the number of coefficients in $W$ and $H$ is typically much smaller than the total number in $V$ (i.e., $KR+RN \\ll KN$), and the product $WH$ can be thought of as a compressed version of the original matrix $V$. As already mentioned before, the column vectors of $W$ are also referred to as **template vectors**, whereas the weights specified by $H$ are called **activations**. As opposed to arbitrary linear combinations as known from linear algebra, the linear combinations occurring in the NMF context only involve nonnegative weights of nonnegative template vectors. As a result, there are no effects such as destructive interferences, where a (positive) component can be canceled out by adding a kind of inverse (negative) component. Instead, the data vectors need to be explained in a purely constructive fashion only involving positive components.\n", "\n", "To find an approximate factorization $V \\approx W \\cdot H$, we need to specify a distance function that quantifies the \n", "quality of the approximation. There are many ways for defining such a distance function, leading to different NMF variants. In the following, we only consider one of these variants, which is based on the **Euclidean distance**. Let $A,B\\in\\mathbb{R}^{K \\times N}$ be two matrices with coefficients $A_{kn}$ and $B_{kn}$\n", "for $k\\in[1:K]$ and $n\\in[1:N]$. Then, the square of the Euclidean distance between $A$ and $B$ is defined by\n", "\n", "\$$\n", " \\|A-B\\|^2:= \\sum_{k=1}^{K}\\sum_{n=1}^{N}(A_{kn}-B_{kn})^2.\n", "\$$\n", "\n", "Based on this distance measure, we can formalize our NMF problem as follows: Given a nonnegative matrix $V\\in\\mathbb{R}_{\\ge 0}^{K \\times N}$ and a rank parameter $R$, minimize\n", "\n", "\$$\n", " \\|V-WH\\|^2\n", "\$$\n", "\n", "with respect to $W \\in \\mathbb{R}_{\\ge 0}^{K \\times R}$ and $H \\in \\mathbb{R}_{\\ge 0}^{R \\times N}$. In other words, regarding $\\|V-WH\\|^2$ as a joint function of $W$ and $H$, the objective is to find a minimum under the nonnegativity constraint for $W$ and $H$. For general matrices, this is a hard computational problem due to several reasons. First, it is in general difficult to enforce hard constraints such as nonnegativity on the solution of an optimization problem. Second, the joint optimization over both matrices $W$ and $H$ leads to computational challenges. In fact, when regarding $\\|V-WH\\|^2$ as a function of $W$ only or $H$ only, one can show that the resulting functions satisfy a strong property referred to as **convexity**. This property, which implies that any local minimum must be a global minimum, makes it possible to apply powerful tools from the field of convex analysis. However, $\\|V-WH\\|^2$ is **not convex** in both matrices together. Therefore, it is unrealistic to expect an algorithm that can solve this problem in the sense of finding a global minimum. However, there are many techniques from numerical optimization that can be applied to find local minima." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Multiplicative Update Rules\n", "\n", "Lee and Seung introduced in their seminal paper a powerful algorithm based on multiplicative update rules to iteratively learn a nonnegative matrix factorization. In the following, we only summarize their algorithm and refer to their paper or Section 8.3.1 of [Müller, FMP, Springer 2015] for a proof. The idea is based on the standard gradient descent approach, which is applied to our problem of minimizing $\\|V-WH\\|^2$ as a function of $W$ and $H$. Since the joint optimization is a very hard problem, one idea is to first fix the factor $W$ and to optimize with regard to $H$, and then to fix the learned factor $H$ and to optimize with regard to $W$. This process is then iterated, where the role of $W$ and $H$ is interchanged after each step. In standard gradient descent the update rules are additive, where a parameter needs to be chosen to control the **step size** towards the direction of the negative gradient. The main trick in the NMF optimization algorithm is that this step size parameter can be set in a specific way so that the **additive** update rules become **multiplicative** update rules. The following table shows the iterative algorithm for learning an NMF decomposition., where the multiplicative update rules are given in matrix notation. The operator $\\odot$ denotes pointwise multiplication and the operator $\\oslash$ pointwise division. \n", "\n", "\n", "\n", "The multiplicative update rules and their properties have a number of remarkable implications.\n", "\n", "* The first implication is that the matrix sequences $W^{(0)},W^{(1)},W^{(2)},\\ldots$ and $H^{(0)},H^{(1)},H^{(2)},\\ldots$ converge. Denoting the limit matrices by $W^{(\\infty)}$ and $H^{(\\infty)}$, the stationarity property implies that these matrices form a local minimum of the distance function $\\|V-WH\\|^2$.\n", "* Second, the multiplicative update rules are extremely easy to implement. \n", "* Third, in practice, the convergence turns out to be relatively fast in comparison with many other methods.\n", "* Fourth, one major benefit of using multiplicative update rules is that the nonnegativity constraints are enforced automatically.\n", "\n", "Indeed, starting with nonnegative matrices $V$, $W^{(0)}$, and $H^{(0)}$, all entries in the equations (1) and (2) of the above table are also nonnegative. Since all operations are multiplicative (either multiplication or division), also the matrices $W^{(\\ell)}$ and $H^{(\\ell)}$ remain nonnegative throughout the iteration.\n", "\n", "In practice, the iteration is performed until a specified stop criterion is fulfilled. For example, one may perform a certain number of iterations $\\ell=0,1,2,\\ldots, L$ for some user-specified parameter $L\\in\\mathbb{N}$. As another stop criterion, one may look at the distances between two subsequently computed template matrices and activation matrices. Specifying a threshold $\\varepsilon>0$, the iteration may be stopped when $\\|H^{(\\ell+1)}-H^{(\\ell)}\\|^2\\leq \\varepsilon$ and $\\|W^{(\\ell+1)}-W^{(\\ell)}\\|^2\\leq \\varepsilon$. Note that this stop criterion does not say anything about the quality of the approximation $V\\approx WH$ achieved by the procedure. Even in the limit case and even when converging to the global minimum (and not to a local one), the distance $\\|V-WH\\|^2$ may still be large. In particular, this may happen if the rank parameter $R$ is chosen too small." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementation\n", "\n", "In the following code cell, we provide a basic implementation that closely follows the NMF procedure described above. Here are some important notes:\n", "\n", "* All operation are performed using matrix-based operations.\n", "* To avoid devision by zero, a small value (machine epsilon) is added to the denominators in the multiplicative update rules.\n", "* In the factorization $V\\approx WH$, there is a degree of freedom with regard to columns and rows of $W$ and $H$, respectively. For example, one may divide a column of $W$ by a factor $\\alpha\\in\\mathbb{R}$ and multiply the correspond row of $H$ by the same factor without changing the product $WH$. In the following code, setting the parameter norm=True, we use this fact to enforce that all columns of the final template matrix $W$ are normalized with regard to the [maximum norm](../C3/C3S1_FeatureNormalization.html).\n", "* If not specified, the matrices $W$ and $H$ are initialized randomly. The power of NMF lies in the possibility of guiding the decomposition process by [initializing the matrices in a musically informed way](../C8/C8S3_NMFAudioDecomp.html)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import os, sys\n", "import numpy as np\n", "import scipy\n", "#import scipy.spatial\n", "import matplotlib.pyplot as plt\n", "import IPython.display as ipd\n", "import librosa, librosa.display\n", "#from matplotlib import gridspec\n", "from numba import jit\n", "\n", "sys.path.append('..')\n", "import LibFMP.B\n", "import LibFMP.C2\n", "import LibFMP.C6\n", "%matplotlib inline\n", "\n", "@jit(nopython=True)\n", "def NMF(V, R, thresh=0.001, L=1000, W=None, H=None, norm=False, report=False):\n", " \"\"\"NMF algorithm with Euclidean distance\n", " \n", " Notebook: C8/C8S3_NMFbasic.ipynb\n", " \n", " Args: \n", " V: Nonnegative matrix of size K x N\n", " R: Rank parameter\n", " thresh: Threshold used as stop criterion\n", " L: Maximal number of iteration\n", " W: Nonnegative matrix of size K x R used for initialization\n", " H: Nonnegative matrix of size R x N used for initialization\n", " norm (bool): Applies max-normalization of columns of final W\n", " report (bool): Reports errors during runtime\n", " \n", " Returns: \n", " W: Nonnegative matrix of size K x R\n", " H: Nonnegative matrix of size R x N\n", " V_approx: Nonnegative matrix W*H of size K x N\n", " V_approx_err: Error between V and V_approx\n", " H_W_error: History of errors of subsequent H and W matrices\n", " \"\"\" \n", " K = V.shape[0]\n", " N = V.shape[1]\n", " if W is None:\n", " W = np.random.rand(K,R)\n", " if H is None:\n", " H = np.random.rand(R,N)\n", " V = V.astype(np.float64)\n", " W = W.astype(np.float64)\n", " H = H.astype(np.float64) \n", " H_W_error = np.zeros((2,L))\n", " ell = 1\n", " below_thresh = False\n", " eps_machine = np.finfo(np.float32).eps\n", " while not below_thresh and ell <= L:\n", " H_ell = H\n", " W_ell = W \n", " H = H * ( W.transpose().dot(V) / (W.transpose().dot(W).dot(H)+ eps_machine) )\n", " W = W * ( V.dot(H.transpose()) / (W.dot(H).dot(H.transpose())+ eps_machine) )\n", " #H = np.multiply( H, np.divide( np.matmul(np.transpose(W), V), np.matmul(np.matmul(np.transpose(W), W), H))) #H+1 = H *p ((W^T * V) /p (W^T * W * H))\n", " #W = np.multiply( W, np.divide( np.matmul(V, np.transpose(H)), np.matmul(np.matmul(W, H), np.transpose(H)))) # W+1 = W *p ((V * H^T) /p (W * H * H^T)) \n", " H_error = np.linalg.norm(H-H_ell, ord=2)\n", " W_error = np.linalg.norm(W - W_ell, ord=2) \n", " H_W_error[:, ell-1] = [H_error, W_error]\n", " if report:\n", " print('Iteration: ',ell,', H_error: ',H_error,', W_error: ',W_error) \n", " if H_error < thresh and W_error < thresh:\n", " below_thresh = True\n", " H_W_error = H_W_error[:,0:ell]\n", " ell += 1 \n", " if norm: \n", " for r in range(R):\n", " v_max = np.max(W[:,r])\n", " if v_max > 0:\n", " W[:,r] = W[:,r] / v_max\n", " H[r,:] = H[r,:] * v_max \n", " V_approx = W.dot(H)\n", " V_approx_err = np.linalg.norm(V-V_approx, ord=2)\n", " return W, H, V_approx, V_approx_err, H_W_error " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Toy Example\n", "\n", "We illustrate the functioning of the NMF procedure by means of a small toy example $V \\in \\mathbb{R}_{\\ge 0}^{K \\times N}$ with $K=4$ and $N=8$. The rank parameter is set to $R=2$." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Matrix V and randomly initialized matrices W and H\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAz4AAACNCAYAAAB2fa+eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df7RdZX3n8fcnNwSQEAOTYIEEAiVFBCtI0LYsKUqVH9LQmbo6pAP+7GI6ioXWGRXtwi7smmHaLgs1/ooQwIEFdlCZ1AmlLIGqVZAEKBgiTYoI16Ak5UcSRDKB7/yx970crufcs/e5+5z9nH0+r7XOInuffZ77vZfnu8/z7L2f51FEYGZmZmZm1mSz6g7AzMzMzMys39zxMTMzMzOzxnPHx8zMzMzMGs8dHzMzMzMzazx3fMzMzMzMrPHc8TEzMzMzs8Zzx8fMzMzMzBrPHR8zMzMzM2s8d3zMzMzMLCmS/oekC+uOo0kkfU/S0XXHUSd3fAZE0pik5yQd0+a96yVdVUdcZoMg6Z8l/YeW7ddJCknvbNl3tKQdkubXE6WZmRVR8py+X/7eYVPK+LSkL3YofyHwTuAL/fodBknS+ZLWSXpe0tVl359y7FGSbpP0jKTNkv59h+OWSvq5pGtbdv8VcMlMfpdh547PgETEC8APgJf1tCUtA84EPl5HXGYD8jTwypbtC4Bnpuz7IPCliHh6kIGZmVlphc/pwLHATuCRKWUcA9zXofx3A2sj4rmygUmaXWRf2TJmaAvw58DqHt9vjev/AF8H9gfOA66V9CttDv8McPeUfWuAN0s6sHjozeKOz2B9H3jNlH1/CfxVRGypIR6zQXkGmAeTV/J+G7iK/Esyv8vzn4BP1xWgmZkVVuacfizwYETElDKOBu7tUP7pwD+27pB0kKSvSNoq6YeS/qjlvUckfUTS/cCzkmZ32HeUpDskPS1pg6Tl05XR819nioj4akTcBPxbL++3eDVwEPDXEfFCRNwG/BNwbutBks4m65x+Y8rP+TmwHnhbT79IA7jjM1gbaLnjI+m3gaVknR+zJpv8kgT+M/Bl4NGWfe8F/ikiflBDbGZmVk6Zc/pxZBd+J0l6FbAAuL9D+a8FHmo5fhbwd8A/AwcDpwAXSjq15TMrgLcD8yNi99R9gPIy/gE4gOyO1HWSjuxSRmvcX887Te1eX+/wu1RJHfZNDqOQNI/scbYPdShjI/C66kMbDu74DNbkHR9JY8ClwMcj4me1RmXWf88A8yTtAfwh2VXA7cAr8y+0DwCX1xifmZkVV+acfiywQtK2iRfZo/+bImJnh/LnAztatk8AFkbEJRGxKyIeBr4InN1yzN9ExGNTHo9r3fdrwFzg0ryM28geGVvRpYxJEXFmRMzv8Dqzy9+sCj8AngD+m6Q9JL0N+E3gFS3HfBK4MiIe61DGDrK/70iq+hlGm94GYGl+ongP8HOy51/Nmu4Zsqt77wD+OSIekvSrZI9FvB3YDfx9jfGZmVlxhc7pkvYEjgJ+j5ePN3k/8MvTlP8UsG/L9qHAQZJax4COAd9q2W7X0G/ddxDwWES82LLvR2R3kKYrIxkR8f8k/Q5ZR/MjwDrgb4HnASQdC/wW2V22TvYlewxuJLnjM1iPkHV2jgP+DFjR5plXsyZ6BjicbADsxfm+7WSPRXwQ+LRzwcxsaBQ6p+cz2c4Cbm29uyNpKdlYk07uB36FlzpLjwE/jIil03ym3XdI674twGJJs1o6P4cA/9KljEmSbgbe1OHtb0XE6dN9vgoRcT/ZXZ6JmL4DXJNvngwsAR6VBNkdrjFJr4mI1+fHHAW0zvQ2Uvyo2wDlDbuNwOeBuyLiH7t8xKwpniH7spgXEf+Q79tO9gjECcDVNcVlZmblFT2nH0c2scHUR9pOoPOMbgBraWncA98DtueTD+ydLxFyjKQTSsR8F/As8OH8MbGTySZluKFoARFxekTM7fDq2OnJJ1bYi+wu1ZikvVonT+j2/pSyfjV//xWS/itwIC/9vVeR3Uk7Nn99Hvi/wKn5Z/cEjgduLfo7N407PoP3fbJBaB+uOxCzAXqG7HGC1nE828lP2NM8521mZukpek4/lilTKuezwB3K9B2fLwFnSNobJpcE+e28vB8C24ArePn02dOKiF3AcrIZ47YBnwXeOaBJdf4UeA74KHBO/u8/LfK+pJslfazl2HOBx8nG+pwCvDUingeIiJ9FxE8mXmTTiP88Irbmn10O3DHKMwnLT5eYmZmZWUok/XfgiYi4rO5YmkLSXcD7IuL7XQ9uKHd8zMzMzMys8fyom5mNNEmrJT0hqe0VMGX+RtJmSfdLen3Le++StCl/vWtwUZuZmVlZ7viY2ai7GjhtmvdPJ1toeClwHvA5AEn7A58A3gi8AfiEpP36GqmZmZn1zB0fMxtpEfFN4MlpDjkL+FJk7gTmSzqQbJacWyPiyYh4imyWnOk6UGZmZlYjd3zMzKZ3MC9f1G4839dpv5mZmSWoLwuYSqpkxoQ999yzimIqLWuvvfaqpByoLqYq/05V/X5VxfTII4+wbds2FTm2Tb27JSJ8Bb5Cs2fPjjlz5sy4nCOOOKKCaDIPPPDAtohYCB3PPRvIFg6esCoiVpX4Ee3qX0yz38w66JCjPlcP0D777BP777//jMvZtWtXBdFUa+7cuZWVtXv37krK2XfffSspB2D79u2VlPPiiy92P6iAp556imeffXao2mh96fhUZfHixZWVtXTpdIv9Fldlg80xdbds2bJSx8+e/VKV3r1794JKgrBJc+bM4cgjj5xxOTfddFMF0WSWLFnyo9bt1joAsHv37p9HRLmK9HLjQOvJaBHZCuDjZKtkt+6/YwY/x2wktMlRn6sHaP/99+fCCy+ccTmPPvpoBdFU6zd+4zcqK+vJJ6d7Arq4k08+uZJyAL7xjW9UUs7OndUsnbdy5cpSx6fQRvOjbtYYkpg1a9bky0bP1DpQUT1YA7wzn93t14BnIuJx4BbgbZL2yyc1eFu+z8w66FOOmlniUmmjJX3Hx6yssbGxukOwmpWtA5KuJ7tzs0DSONlMbXsARMTngbXAGcBm4GfAe/L3npT0SV5alfySiKjmEqFZg/k8bTaaUsh9d3ysMSQlkVRWn17qQESs6PJ+AB/o8N5qYHWpH2g2wnyeNhtNqeR+4Y6PpDFgHfDjiDizfyGZ9SaVpBo2Tcpt1wGztDlH+6NJ53FrplRyv8wdnwuAjcC8PsViNmN+Xrwnjcpt1wGztDlH+6JR53FrphRyv1AEkhYBbweu6G84Zr2buJow8bLumpbbU+uA64FZWpyj1WvaedyaKZU2WtE7PpcBHwaqm4zcrA/8JVpa43LbdcAsbc7RyjXuPG7NlELud73jI+lM4ImIWN/luPMkrZO0rrLozEro5WqCpPmSbpT0A0kbJf16n8NMRi+5XdWCbv3iq8lmaXOOVquX8/izzz47oOjMXpLKHZ8ij7qdCCyX9AhwA/AWSddOPSgiVkXEshkuFGg2Iz0k1eXA30fEq4HXkT0jPSpK5/bUhQdT5EaVtZOvtxSSDpuy/9OSvlhXXKPIOVqp0ufxffbZZ9AxmgE9tdEq17XjExEXRcSiiFgCnA3cFhHn9D0ys5LKLo4laR5wEnAlQETsioin+xxmMpqY214c0aZxLLATeGTK/mOA+wYezYhyjlariedxa6ZUFjD1GccapeTVhMOBrcBVku6VdIUkXwobcr6abB0cCzyYr8vU6mjg3hriGVnOUbPRlELel+r4RMQdnh/eUtXm+dEFE88056/zpnxkNvB64HMRcRzwLPDRAYedhKbktscP2DSOA77fukPSq4AFwP21RDSCnKP905TzuDVTKmN80n9g36ygNotjbesy5mwcGI+Iu/LtGxnRjk9TpLJAmiXpWOAISWe17BsDNkXEzppiGjnOUbPRlEruu+NjjVLmudGI+ImkxyQdGREPAacAD/YtOBsIjxmwqSTtCRwF/B5wd8tb7wd+uZagRphz1Gw0pZD77vhYY/R4NeGDwHWS5gAPA++pPDAbmFSuKFlyjiF7tPvW1rs7kpYC004DbNVyjpqNplRy3x0fa5SySRUR9wGegr1BUjixWnKOI5vYYOojbSeQz+pog+McNRtNKeR+0h2fTZs21R2CDZFUriY02ezZs1m4cOGMy1myZMnMg2nDdcA6OJaXP+KGpIXAoXgq64HqJUclnUa25toYcEVEXDrl/b8G3pxvvgI4ICLm5++9ADyQv/doRCyfQfiNsGXLFi655JIZl3PKKadUEE3ma1/7WiXlVLnW3Nq1aysp58EHq3uC/rOf/Wwl5axYsaKScp577rnCx6by/Zx0x8esrBSeH7V69VIH3LBqtog4v82+rYBqCGfklclRSWPAZ4C3kk1Ic7ekNREx2ZqMiD9uOf6DZHf4JjwXEcfOOGgzm7EU2mju+FhjpHI1werT49VkN6zMBqSHHH0DsDkiHs4/fwNwFp0nolkBfGJGQZpZ5VJpo9Xf9TKrUApzxFu9elgjZLJhFRG7gImGVScrgOsrCNVsJJXM0YOBx1q2x/N9v0DSocBhwG0tu/fK13G7U9LvzCRuM5uZFNpovuNjjZHK1QSrT491oF3D6o0dyu/YsAJ2A5dGxE1lAzAbFR1ydEGeQxNWRcSqiY+0KSY6FH82cGNEvNCy75CI2CLpcOA2SQ9ExL/2FLyZ9SyVNpo7PtYYqSSV1aeHRhW4YWU2MB1ydLrFpseBxS3bi4AtHY49G/hA646I2JL/92FJd5A9pur8NBuwVNpo7vhYo6QwcM7q1aYOTNeoAjeszAaq5Hn6bmCppMOAH5Pl4O9PPUjSkcB+wHdb9u0H/Cwinpe0ADgR+IsZhG5mM5BCG80dH2uMVK4mWH16rAPJN6zGxsai6rrdjynFx8fHKy/zgAMOqLzMbdu2VV4m9GeNiu3bt1de5mGHHVZpeVu3bmXHjh2FZsgrm6MRsVvS+cAtZLMuro6IDZIuAdZFxJr80BXADRHRerf2KOALkl4kG9N8aeukJWY2OGVzX9Jq4EzgiYg4ps37IpuN9QzgZ8C7I+KebuW642ON4o6P9bCIbfINq7GxMX7pl36p0jJXrlxZaXkAH/nIRyov8/3vf3/lZa5evbryMgH222+/ysu8+eabKy+zijVcWl188cWlju8hR9cCa6fsu3jK9p+1+dx3gNeW+mFm1jclc/9qYCXwpQ7vnw4szV9vBD5Hh/G5rdzxscbwHR/rtQ64YWU2GD5Pm42mHu72flPSkmkOOQv4Un4x8k5J8yUdGBGPT1euOz7WKCk8P2r1ch0wS5tz1Gw0Tcn9bhMPddNpqnt3fGw0+EqiuQ6Ypc05ajaa2uR+t4mHuhbZZl+nGVknueNjjeEvVHMdMEubc9RsNPUh98vMyDrJ95utUVJYFdjqVXJVeDMbMOeo2WiqOO/XAO9U5teAZ7qN7wHf8bEG8ZVEG5Y60G2aTrOmGpYcNbNq9TCd9fXAyWRjgcaBTwB7AETE58kmJDoD2Ew2nfV7ipTrjo81igfN2pDUgauZfppOs8Yakhw1s4qVyf2IWNHl/WDKguJFuONjjeEriTYsdaDANJ1mjTQsOWpm1Uol993xsUZJIamsXk2pA5LOA86D5vxOZuD6bDaqUsj9vnR8jj/+eNatW9f9QLMK9Xo1QdIYsA74cUScWXlgDXLQQQeVXqW9ne9+97sVRJPZuXPn5L9TuaJUhXw9g1UAc+bM6TpFp9kwaFKODqt58+bx5je/ecblrFgx7ZNIpcyeXU1z9IQTTqikHIAdO3ZUUs7ChQsrKQfgpJNOqqSc0047rZJyvvOd7xQ+NpXc9x0fa5Qenx2/ANgIzKs2GquDxw+Ypc05ajaaUsj9+iMwq8jE1YQyUyVKWgS8Hbiir8HZQEytAylcXTKzlzhHzUZTL220fnDHxxqjTVItkLSu5XVem49dBnwYeHGgwVpfDEujKp+m87vAkZLGJb2v7pjMBmFYctTMqpVKx8ePulmjTEmmbRGxrNOxkibWUVkv6eR+x2aDMQwNqW7TdJo12TDkqJlVL4Xcd8fHGqOHgXMnAsslnQHsBcyTdG1EnNOXAK3vUhk8aWbtOUfNRlMqud/1UTdJiyXdLmmjpA2SLhhEYGa9mDVr1uSrm4i4KCIWRcQS4GzgtlHr9DQxv1vrQNGBlJJOk/SQpM2SPtrm/XdL2irpvvz1By3vvUvSpvz1rgp/FbNGKpujzs/OmngOt+Yq+93cD0Xu+OwGPhQR90jaF1gv6daIeLDPsZmVksrVhCHTqPzupQ7k05l/BngrMA7cLWlNm7/BlyPi/Cmf3R/4BLAMCLK/35qIeKrX38GsycrmqPOzq0adw625Ummjde1yRcTjEXFP/u8dZNP+HtzvwMx60evAuYi4YxTX8GlifvcwcPoNwOaIeDgidgE3AGcV/HGnArdGxJN5Y+pWoJoFEswaqmSOOj+n0cRzuDVXCpMblLrXJGkJcBxwVz+CMZuJVGYMGVZNyO8eZ4w6GHisZXuc9g2H35V0v6QbJS0u+Vkzo2OOTjcDp/OzoCacw625UmmjFZ7cQNJc4CvAhRGxvc375wHnARxyyCGVBWhWRgqLYw2j6fK7Nbdf9apX1RBdOW3qwAJJ61q2V0XEqpZttSkmpmz/HXB9RDwv6Q+Ba4C3FPzsjO2zzz4cf/zxlZZZ1crdrQ444IDKy3zwweqf2Hn1q19deZkAF154YeVlHnjggZWXefnll1da3k9/+tNSx7fJ0elm4Ew+P1NQpo229957Dzg6s0wKbbRCEUjagyyhrouIr7Y7JiJWRcSyiFi2cOHCKmM0KySVqwnDplt+t+b2/PnzBx9gCR2uJm+biD9/rZrysXFgccv2ImBL6wER8W8R8Xy++UXg+KKfNbOX9HBX1vnZRdk22p577jnYAM1Ip41WZFY3AVcCGyPiU/0Pyaw3qSTVMGlafvf4qNvdwFJJh0maQzbD35op5bZedl9O9hw9wC3A2yTtJ2k/4G35PjNro4ccdX5Oo2nncGuuVNpoRR51OxE4F3hA0n35vo9FxNr+hWXWG3d4SmtcfvcwscVuSeeTNYjGgNURsUHSJcC6iFgD/JGk5WQzKD0JvDv/7JOSPknWOAO4JCKerOY3MWumMjnq/Oyqcedwa64U2mhdOz4R8W3aPydrlpRUpkocJk3L717rQN5IWDtl38Ut/74IuKjDZ1cDq0v/ULMR1EuOOj87a9o53JorlTZa4ckNzIZBCgPnrF6uA2Zpc46ajaYUct8dH2uMVK4mWH1cB8zS5hw1G02p5H79XS+zCqUwcM7q1cPkBrWQdJqkhyRtlvTRuuMxG5RhyVEzq1YKee87PtYYqVxNsPoMSx2QNAZ8Bngr2ZS7d0taExHVL1pjlpBhyVEzq1YquZ90x2fz5s2VlbVp06ZKyqkypqocccQRlZW1dOnSSsqpMqYyUnh+tMkeeugh3vSmN824nA0bNlQQTeboo49+2faQ1IE3AJsj4mEASTcAZwHu+FjjDUmONtacOXMqWWj+mmuuqSCazIsvvlhJOVu3bq2kHKhuwe4zzjijknIAnnvuuUrKufPOOyspZ+fOnaWOTyH3k+74mJWRytUEq88Q1YGDgcdatseBN7Ye4JXWrYmGKEfNrEKp5H79XS+ziqSyOJbVp8cFTOvQbvrZeNmGV1q3BhqiHDWzCvXSRus2FlbSuyVtlXRf/vqDbmX6jo81ir9EbUjqwDiwuGV7EbClpljMBmpIctTMKlYm90uMhf1yRJxftFzf8bHGkMSsWbMmXzZ6ptaBhOvB3cBSSYdJmgOcDaypOSazvhuiHDWzCvXQRpscCxsRu4CJsbAz4js+1ii+kmjDUAciYrek84FbgDFgdURUN+ODWcKGIUfNrHpTcn+BpHUt26siYlXLdtexsLnflXQS8C/AH0fEY22OmeSOjzVGKgPnrD7DVAciYi2wtu44zAZpmHLUzKrTJve3RcSy6T7SZl9M2f474PqIeF7SHwLXAG+ZLg7fY7bGKDtwTtJiSbdL2ihpg6QLBhCm9ZEHTpulzTlqNpp6mNyg61jYiPi3iHg+3/wicHy3Qn3Hxxql5JfobuBDEXGPpH2B9ZJu9SKSw80NKbO0OUfNRlPJ3J8cCwv8mGws7O+3HiDpwIh4PN9cDmzsVqg7PtYYEwPnisqT5fH83zskbSR7ptQdnyFVtg6Y2WA5R81GUw9ttLZjYSVdAqyLiDXAH0laTnYh+0ng3d3KdcfHGqXkwLlJkpYAxwF39S04G4heriZLOg24nOzkekVEXDrl/T8B/oDs5LoVeG9E/Ch/7wXggfzQRyNiee/RmzVf2Rx1fpo1Q9ncbzcWNiIubvn3RcBFZcp0x8cao4eBcxOfmwt8BbgwIrb3Kz7rv14GThdcK+BeYFlE/EzSfwH+AviP+XvPRcSxM4/erPnK5qjz06wZUpnYxB0fa5QeGr17kHV6rouIr/YlKBuoHk6sk2sFAEiaWCtgsmEVEbe3HH8ncM4Mwyzl6aef3nbTTTf9qODhC4Bt/Yynk5/+9KdlDi8U52WXXdZzPBUp/Pe8+uqr+xvJ9Gr7/w4cWubgkjmafH6aWTHu+JhVqOzzo5IEXAlsjIhP9S0wG5gexw8UXStgwvuAm1u298ofqdwNXBoRN5UNoJuIWFj0WEnritzprJvjrNYQxdkuR6d7LDn5/DSz7lIZ3+eOjzVKyasJJwLnAg9Iui/f97H8mVIbUm3qQLexXkXWCsgOlM4BlgG/2bL7kIjYIulw4DZJD0TEv/YQutlIaJOj0z2W7Pw0awjf8TGrUNnnRyPi27T/UrUh1aEOdBvr1XWtgLzs3wI+Dvxmy7oBRMSW/L8PS7qDbJIMN6zM2ujhOX/np1kDpDLGp/57TmYV6WFxLGuYHhdHnFwrQNIcsrUC1kwp9zjgC8DyiHiiZf9+kvbM/72A7C5i3dOht525MEGOs1pDEWcPOdq0/DQbSam00XzHxxrFHR7rYbrMImsF/CUwF/jf2dCwyWlxjwK+IOlFsgtJl9a9AG6nKdtT4zirNSxxQrkcbVp+mo2yFNpo7vjUaNOmTZWUc/PNN3c/qKCqYtq8eXMl5ZSRysC5Jlu6dCkrV66ccTnnn39+BdH8ol7rQIG1An6rw+e+A7y29A80G1G95Kjzs1oRwQsvvDDjct773vdWEE3m9ttv735QAYceWmqCwWlt3bq1knLuvffeSsqB0jNndnTGGWdUUk6Z9mcqbbT6IzCrUAq3Ua1ePTzq1giSTpP0kKTNkj5adzztSFos6XZJGyVtkHRB3TFNR9KYpHslfb3uWKYjab6kGyX9IP/b/nrdMU1nVHPUbNSlkPe+42ONkcrAOavPqNaBgos8pmA38KGIuEfSvsB6SbcmGOeEC4CNwLy6A+nicuDvI+Id+TiYV9QdUCejmqNmoy6V3PcdH2uUFK4mWL1G9Gry5CKPEbELmFjkMSkR8XhE3JP/ewdZp+LgeqNqT9Ii4O3AFXXHMh1J84CTyNYkIyJ2RcTT9UY1vRHNUbORl0Leu+NjjTHx/OjEy0bP1DowQvWg3SKPSXYoJkhaQja18F31RtLRZcCHgRfrDqSLw4GtwFX5Y3lXSNqn7qA6GeEcNRtpqbTRfMaxRknhaoLVa0SvJhde5DEFkuYCXwEujIjtdcczlaQzgSciYn3dsRQwG3g98LmIOA54FkhyjNeEEc1Rs5GXQt4X6vgMw6BZs1TmiB82TcrvHtfxaYJCizymQNIeZJ2e6yLiq3XH08GJwHJJj5A9NvgWSdfWG1JH48B4REzcObuRrCOUpBHO0b5p0jncmiuVNlrXjk/LoNnTgdcAKyS9pt+BmZWVSlINk6bl9wg3qrou8pgCZYusXAlsjIhP1R1PJxFxUUQsioglZH/L2yLinJrDaisifgI8JunIfNcpJLxI5wjnaF807RxuzZVKG63IrG6Tg2YBJE0Mmk32xGqjy1+ipTUuv0exDnRa5LHmsNo5ETgXeEDSffm+j+XrtFjvPghcl3d6HwbeU3M80xrFHO2jxp3DrblSyP0iHZ92g2bf2J9wzHqXyuJYQ6ZR+T3KdaDdIo+piYhv0348UrIi4g7gjprDmFZE3AcsqzuOIkY5R/ukUedwa65Ucr9Ix6fQoFlJ5wHnARxyyCEzDMusNylcTRgyXfO7NbcPOOCAQcQ0I64DZmlzjlaqdBtt7ty5/Y7JrK0Ucr9I16vQoNmIWBURyyJi2cKFC6uKz6ywVJ4fHTJd87s1t1/5ylcONLiyPH7ALG3O0cqVbqPtvffeAwvObEIqbbQiHZ+hGDRrBmlMlThkGpffblSZpc05WqnGncOtuVLI+66Pug3RoFkbcak8PzpMmpbfrgNmaXOOVqtp53BrrlRyv8gYn6EYNGs2cRvVymlSfrsOmKXNOVq9Jp3DrblSyf36u15mFSp7G9ULvzVPL4/RdKsHkvaU9OX8/bskLWl576J8/0OSTq3sFzFrqLI56vw0a4ZBfjd34o6PNUbZgXNe+K15ehk4XbAevA94KiKOAP4a+J/5Z19D9kz90cBpwGfz8sysjbI56vw0a4Y+tdHa5v503PGxRil5NWFy4beI2AVMLPxmQ6yHOz5F6sFZwDX5v28ETpGkfP8NEfF8RPwQ2JyXZ2YdlMxR56dZQwzwu7kjd3ysMSYGzk28Cmi38NvBfQnOBmJqHaiwHkweExG7gWeAf1fws2aW6yFHnZ9mDdCnNlqn3O+o0OQGZa1fv36bpB91OWwBsK0fP38GHFMxg4zp0KIHrl+//paxsbEFLbv2krSuZXtVRKxq2S608Ju9ZNOmTdtOPfXUbrkN1dWRIuVM1pE2dQCqqQedjnEdMiuhhxx1flZs69at21auXDnjNtrKlSurC6qYoYypBl1juuqqq6r6WXW30UrneF86PhHRdQVTSesiYlk/fn6vHFMxKcYEEBGnlfxIoYXf7CVFchuqqyNly+mhDkCxejBxzLik2cArgScLftbMcn06Tzs/S3AbrTqOqbgB535HftTNRpkXfjMoVg/WAO/K//0O4LaIiHz/2fnMMocBS4HvDShus1Hg/DQbTTPJ/Y76csfHbJLbG2EAAAM6SURBVBh44TeDzvVA0iXAuohYA1wJ/C9Jm8muJp2df3aDpL8FHgR2Ax+IiBdq+UXMGsj5aTaaZpL701GXjlHfSDpvyrN8tXNMxaQYk6WlqjriumZmNngpnnsdUzEpxpSS2jo+ZmZmZmZmg+IxPmZmZmZm1ngD7/hIOk3SQ5I2S/rooH9+m3gWS7pd0kZJGyRdUHdMEySNSbpX0tfrjgVA0nxJN0r6Qf73+vW6Y7K0VJXfklZLekLS96uMz8zMOnMbrTi30YbTQB91kzQG/AvwVrIp6O4GVkTEgwML4hdjOhA4MCLukbQvsB74nTpjmiDpT4BlwLyIODOBeK4BvhURV+QzbLwiIp6uOy5LQ5X5LekkYCfwpYg4ptJAzczsF7iNVo7baMNp0Hd83gBsjoiHI2IXcANw1oBjeJmIeDwi7sn/vQPYSAIrO0taBLwduKLuWAAkzQNOIptBg4jY5YSyKSrL74j4Jl3m4jczs0q5jVaQ22jDa9Adn4OBx1q2x0mgAk+QtAQ4Drir3kgAuAz4MPBi3YHkDge2Alflt3avkLRP3UFZUpLObzMzm1bS53C30ablNlpBg+74qM2+JKaVkzQX+ApwYURsrzmWM4EnImJ9nXFMMRt4PfC5iDgOeBao/flfS0qy+W1mZl0lew53G60rt9EKGnTHZxxY3LK9CNgy4Bh+gaQ9yBLquoj4at3xACcCyyU9Qnar+S2Srq03JMaB8YiYuNJyI1mSmU1IMr/NzKyQJM/hbqMV4jZaQYPu+NwNLJV0WD7w6mxgzYBjeBlJInsmcmNEfKrOWCZExEURsSgilpD9jW6LiHNqjuknwGOSjsx3nUK2GrbZhOTy28zMCkvuHO42WuGY3EYraPYgf1hE7JZ0PnALMAasjogNg4yhjROBc4EHJN2X7/tYRKytMaZUfRC4Lj8hPgy8p+Z4LCFV5rek64GTgQWSxoFPRMSVlQVrZmYv4zba0HMbrYCBTmdtZmZmZmZWh4EvYGpmZmZmZjZo7viYmZmZmVnjueNjZmZmZmaN546PmZmZmZk1njs+ZmZmZmbWeO74mJmZmZlZ47njY2ZmZmZmjeeOj5mZmZmZNd7/B/YQdxVVxHvwAAAAAElFTkSuQmCC\n", "text/plain": [ "