--- title: "Neural-network imputation with deepImp" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Neural-network imputation with deepImp} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- > This vignette is **pre-computed**: its outputs were generated on a machine with a > working `torch` backend (see `vignettes/deepImp.Rmd.orig` and > `vignettes/precompute.R`). This keeps the package check fast and free of a Python > or `libtorch` dependency while still showing real results. ## Why neural-network imputation? Most real data sets contain missing values, and discarding incomplete rows wastes information and can bias an analysis. *Imputation* fills the gaps with plausible values so the data can be analysed as a whole. Classical imputation methods (mean imputation, $k$-nearest neighbours, regression/PMM) assume a fixed, usually linear, relationship between a variable and its predictors. `deepImp` instead learns that relationship with a **neural network**, which can capture non-linear and interaction structure without the user specifying it. It imputes one variable at a time, using all the others as predictors, and cycles through the variables repeatedly until the imputations stabilise (an *iterative, chained* scheme). The package handles two settings with the same engine: * **mixed-type data** -- numeric, count, semi-continuous, binary and nominal variables -- via `impNNet()`; * **compositional data with rounded zeros** -- values that fall below a detection limit -- via `impNNetCoDa()`. Networks run on a native [torch](https://torch.mlverse.org/) backend by default (no Python required), or optionally on [keras3](https://keras3.posit.co/). ``` r library(deepImp) ``` ## The architecture: `deepimp_arch()` The network is described by a `deepimp_arch()` object. Rather than hard-coding a topology, you choose the depth and width (`hidden`), the regularisation (`dropout`, `batchnorm`), the `activation`, and the `optimizer` and its `learning_rate`: ``` r arch <- deepimp_arch(hidden = c(128, 64, 32), dropout = 0.1, activation = "relu") arch #> #> hidden layers : 128 -> 64 -> 32 #> dropout : 0.1 #> activation : relu #> batchnorm : TRUE #> optimizer : adam (lr = 0.001 ) ``` `hidden = c(128, 64, 32)` means three hidden layers of 128, 64 and 32 units. For small data sets or quick runs, `deepimp_arch_small()` is a lighter preset: ``` r deepimp_arch_small() #> #> hidden layers : 64 -> 32 #> dropout : 0.1 #> activation : relu #> batchnorm : TRUE #> optimizer : adam (lr = 0.001 ) ``` The same architecture object is translated to whichever backend you select, so a configuration is portable across `torch` and `keras3`. ## Imputing mixed-type data with `impNNet()` We use the `sleep` data from the `VIM` package -- mammal sleep measurements that already contain genuine missing values. ``` r data(sleep, package = "VIM") str(sleep) #> 'data.frame': 62 obs. of 10 variables: #> $ BodyWgt : num 6654 1 3.38 0.92 2547 ... #> $ BrainWgt: num 5712 6.6 44.5 5.7 4603 ... #> $ NonD : num NA 6.3 NA NA 2.1 9.1 15.8 5.2 10.9 8.3 ... #> $ Dream : num NA 2 NA NA 1.8 0.7 3.9 1 3.6 1.4 ... #> $ Sleep : num 3.3 8.3 12.5 16.5 3.9 9.8 19.7 6.2 14.5 9.7 ... #> $ Span : num 38.6 4.5 14 NA 69 27 19 30.4 28 50 ... #> $ Gest : num 645 42 60 25 624 180 35 392 63 230 ... #> $ Pred : int 3 3 1 5 3 4 1 4 1 1 ... #> $ Exp : int 5 1 1 2 5 4 1 5 2 1 ... #> $ Danger : int 3 3 1 3 4 4 1 4 1 1 ... colSums(is.na(sleep)) #> BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred #> 0 0 14 12 4 4 4 0 #> Exp Danger #> 0 0 ``` `impNNet()` first infers the measurement scale of each column with `guessType()` (you can override this with the `vartypes` argument). This determines whether a column is modelled by regression (numeric/count) or classification (binary/nominal): ``` r guessType(sleep)$type #> [1] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" #> [8] "count" "count" "count" ``` Now we impute. We use the small architecture and a modest number of epochs to keep the example quick; in practice you would train longer. The `seed` makes the run reproducible (exactly so when `dropout = 0`): ``` r imp <- impNNet(sleep, arch = deepimp_arch_small(), epochs = 50, iterations = 3, seed = 1) imp #> #> observations : 62 #> variables : 10 #> imputations : 1 #> backend : torch #> architecture : #> #> hidden layers : 64 -> 32 #> dropout : 0.1 #> activation : relu #> batchnorm : TRUE #> optimizer : adam (lr = 0.001 ) ``` The result is a `"deepimp"` object. `summary()` reports how many cells were imputed per variable, and `getImputed()` (or `as.data.frame()`) returns the completed data frame: ``` r summary(imp) #> summary #> observations : 62 #> variables : 10 #> imputations : 1 #> missing per variable: #> BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred #> 0 0 14 12 4 4 4 0 #> Exp Danger #> 0 0 #> converged in : 3 iteration(s); last change = 0.4711 completed <- getImputed(imp) ``` To see the imputation at work, look at cells that were missing before. For example, the non-dreaming sleep variable `NonD` had missing entries; here are the values `deepImp` filled in: ``` r missing_rows <- which(is.na(sleep$NonD)) round(completed$NonD[missing_rows], 2) #> [1] 4.21 10.04 12.37 4.19 5.01 7.77 10.43 9.16 8.03 6.29 10.55 2.90 #> [13] 9.37 12.10 ``` `plot()` shows the convergence of the chained imputation -- the mean change in the imputed cells from one iteration to the next, which should decrease: ``` r plot(imp) ``` ![plot of chunk impnnet-plot](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAEDmlDQ1BrQ0dDb2xvclNwYWNlR2VuZXJpY1JHQgAAOI2NVV1oHFUUPpu5syskzoPUpqaSDv41lLRsUtGE2uj+ZbNt3CyTbLRBkMns3Z1pJjPj/KRpKT4UQRDBqOCT4P9bwSchaqvtiy2itFCiBIMo+ND6R6HSFwnruTOzu5O4a73L3PnmnO9+595z7t4LkLgsW5beJQIsGq4t5dPis8fmxMQ6dMF90A190C0rjpUqlSYBG+PCv9rt7yDG3tf2t/f/Z+uuUEcBiN2F2Kw4yiLiZQD+FcWyXYAEQfvICddi+AnEO2ycIOISw7UAVxieD/Cyz5mRMohfRSwoqoz+xNuIB+cj9loEB3Pw2448NaitKSLLRck2q5pOI9O9g/t/tkXda8Tbg0+PszB9FN8DuPaXKnKW4YcQn1Xk3HSIry5ps8UQ/2W5aQnxIwBdu7yFcgrxPsRjVXu8HOh0qao30cArp9SZZxDfg3h1wTzKxu5E/LUxX5wKdX5SnAzmDx4A4OIqLbB69yMesE1pKojLjVdoNsfyiPi45hZmAn3uLWdpOtfQOaVmikEs7ovj8hFWpz7EV6mel0L9Xy23FMYlPYZenAx0yDB1/PX6dledmQjikjkXCxqMJS9WtfFCyH9XtSekEF+2dH+P4tzITduTygGfv58a5VCTH5PtXD7EFZiNyUDBhHnsFTBgE0SQIA9pfFtgo6cKGuhooeilaKH41eDs38Ip+f4At1Rq/sjr6NEwQqb/I/DQqsLvaFUjvAx+eWirddAJZnAj1DFJL0mSg/gcIpPkMBkhoyCSJ8lTZIxk0TpKDjXHliJzZPO50dR5ASNSnzeLvIvod0HG/mdkmOC0z8VKnzcQ2M/Yz2vKldduXjp9bleLu0ZWn7vWc+l0JGcaai10yNrUnXLP/8Jf59ewX+c3Wgz+B34Df+vbVrc16zTMVgp9um9bxEfzPU5kPqUtVWxhs6OiWTVW+gIfywB9uXi7CGcGW/zk98k/kmvJ95IfJn/j3uQ+4c5zn3Kfcd+AyF3gLnJfcl9xH3OfR2rUee80a+6vo7EK5mmXUdyfQlrYLTwoZIU9wsPCZEtP6BWGhAlhL3p2N6sTjRdduwbHsG9kq32sgBepc+xurLPW4T9URpYGJ3ym4+8zA05u44QjST8ZIoVtu3qE7fWmdn5LPdqvgcZz8Ww8BWJ8X3w0PhQ/wnCDGd+LvlHs8dRy6bLLDuKMaZ20tZrqisPJ5ONiCq8yKhYM5cCgKOu66Lsc0aYOtZdo5QCwezI4wm9J/v0X23mlZXOfBjj8Jzv3WrY5D+CsA9D7aMs2gGfjve8ArD6mePZSeCfEYt8CONWDw8FXTxrPqx/r9Vt4biXeANh8vV7/+/16ffMD1N8AuKD/A/8leAvFY9bLAAAAOGVYSWZNTQAqAAAACAABh2kABAAAAAEAAAAaAAAAAAACoAIABAAAAAEAAAJYoAMABAAAAAEAAAGQAAAAANOlcRYAAEAASURBVHgB7J0J3Ez198ePfV8SsodsIUu2smTLTvZkKSqkQogQFbKFJAmhkDVbkohkCRWJZM2WLUt2KVnnfz/n/7vTPGPmmTuPO/PM8jmv1zxz597v/S7vO8/cc8/3fM9J4DBEKCRAAiRAAiRAAiRAArYRSGhbTayIBEiABEiABEiABEhACVDB4heBBEiABEiABEiABGwmQAXLZqCsjgRIgARIgARIgASoYPE7QAIkQAIkQAIkQAI2E6CCZTNQVkcCJEACJEACJEACVLD4HSABEiABEiABEiABmwlQwbIZKKsjARIgARIgARIgASpY/A6QAAmQAAmQAAmQgM0EqGDZDJTVkQAJkAAJkAAJkAAVLH4HSIAESIAESIAESMBmAlSwbAbK6kiABEiABEiABEiACha/AyRAAiRAAiRAAiRgMwEqWDYDZXUkQAIkQAIkQAIkQAWL3wESIAESIAESIAESsJkAFSybgbI6EiABEiABEiABEqCCxe8ACZAACZAACZAACdhMgAqWzUBZHQmQAAmQAAmQAAlQweJ3gARIgARIgARIgARsJkAFy2agrI4ESIAESIAESIAEqGDxO0ACJEACJEACJEACNhOggmUzUFZHAiRAAiRAAiRAAlSw+B0gARIgARIgARIgAZsJUMGyGSirIwESIAESIAESIAEqWPwOkAAJkAAJkAAJkIDNBKhg2QyU1ZEACZAACZAACZAAFSx+B0iABEiABEiABEjAZgJUsGwGyupIgARIgARIgARIgAoWvwMkQAIkQAIkQAIkYDMBKlg2A2V1JEACJEACJEACJEAFi98BEiABEiABEiABErCZABUsm4GyOhIgARIgARIgARKggsXvAAmQAAmQAAmQAAnYTIAKls1AWR0JkAAJkAAJkAAJJCYCEiCB4BNYu3atXLhwQZIkSSL169cPfgd8tLhjxw45cOCAlqpevbqkTZvWxxk8TAIkQAIk4EoggcMQ1x3cJgESCDyBMmXKyJYtWyRNmjRy+fLlwDfoZwvdunWT999/X8/69ddf5aGHHvKzBhYnARIggegmwCnC6L7+HD0JkAAJkAAJkEAACHCKMABQWSUJhDsBWLCeeuopHcYDDzwQ7sNh/0mABEgg6ASoYAUdORuMNgKbNm0S+FxdvHhRHn30UWnQoEGsCK5evSpLly6VrVu3SooUKaREiRJSu3ZtSZo06R3nWSn7119/yaJFi/TcRx55RHLmzCmrVq2SjRs36na9evUkT548Meo+c+aM/Pbbb7qvYMGCkjJlSjl58qSsXLlS99WqVUsSJ04sy5YtE/hrFStWTH3J7rnnHvn777+1/z///LMUKlRIatSooe2YDaxfv14OHTokyZMnlxYtWsiePXu03tOnT0u5cuWUT8KE1o3rJ06c0PHs3LlTEiVKJOjvk08+qX022zTfb926Jdu2bZPvvvtO0F7RokWlWrVqkj17drOIvl+7dk3mzp2r26VKldJxmMzQ7ypVqkiFChX0+NGjR2XNmjW6XaRIESldurRum39++eUX2b59u36EP1uOHDl028q1Q0F8d44cOaLTyfjuTJ06VY4dOyY1a9aUSpUqaV24xosXLxZM5+bNm1fq1KkjGTNmlIULF+rxsmXLyoMPPqjb5h+URd2HDx/W8T322GP6bh7HO767e/fu1WvdunVrLYvvAL4bGGuTJk0kffr0rqfotj/XxEo/7miAO0ggHAjAB4tCAiRgP4Hbt287evToAR/HGC9DoXEYSoDuM3ywYjRs3Gwcxo0wRnmcX7JkSYdxU4tT2f379zvrGzBggNbl2qfUqVM7Pvvssxh1v/LKK85z0CfIN99849w3cuRIh6EoOD+jPkPJchjKi6N48eIx9mfIkCFG340btR7PlCmT46OPPnIYylSM8obi4Dh//nyM/nj7MG3aNIfhgB/jfPQFdRs+bjFOAwf3PqMsxj9+/PgYZf/8809nna+//rrDUCidn012/fr103POnTvnMJQuPW4owzHqwQdDEdNjhrLsHJc/17lx48Z6vqEEO5577jlnPwwFR9tCXfny5XPuR//wvQJbs69jxoxx9stQMh3ouzt3Q2F2DBs2zIHvrSkvvfSS1oG+G0q6w1C0nXWi7vvvv99hLIYwi+u71WviTz9iNMAPJBAmBCRM+slukkDYETAsIM6bUYIECRyG9cFRvnx55z7zRmgO7N9//3UY1gfnccOy4mjYsKHzRogbKm5KEH/KuipYhoVH6zesFQ7DuuRsK1myZI7ff/9d68YfXwoW+g4lEeUKFy7srAf7DSuWo3379jHqR3ummAqWeYNH+UaNGjmMFZXOejp27GgW9/r+ww8/OMAVbeLdsJQ5DMuSk1fWrFkdhpVIzzcsZg7Dcues37AkOurWrRtDYTAsQ862XBUsMDOsNI62bds6mjdv7qwD+w3Lkp7Tpk0b535TIcUBV/ZPP/20lvXn2uEEU8Eyx5oqVSoHlKGhQ4fq98FYgOBs27DEqRIGpdMsDz6uCtbkyZOd5e+9915Hhw4dYiiersq2qWChLlwvwxLm6Nq1qyN37tzOOnC+Kf5cE3/6YdbPdxIIJwJUsMLparGvYUXAVfFYvny5s++uNxZXC9Z7773nvGkNGTLEWd7VEjFnzhzd709Z15s8brbG9KPWAWVt0KBBzjZdb5S+FCxYLowpT61n9+7dzjpwI163bp3uv3nzplqS0CYUFFNMBQv7q1at6jCmt/QQlJVcuXJpXa7Ki3me+7sxPaZlceP/8ccfnYe7dOmiygWUAGMKTPe3atXK2UdXZcOYnnQYU696DIqeaTlzVbCg+G3evNlZv6syZV5XY9rTWf+rr77qLNu/f3/nfpSB+HPtUN5UsMALCiSURmMKV/u6YMECZ/1QUo2pTZyi1wbWNJyDF9qE4HjmzJl1H67JlStXdP+NGzec7GFBNa1YpoKFOlC/Ka7fKWNa19ztsHpN/O2HswFukEAYEaCCFUYXi10NHwLXr193mNYiWAlMyxNGgO106dLpTc5VwTL8hpw3RNzQz549qy/czMy6TMuOP2Vdb4bZsmWLARE3Okz/4AZq+Bo5j/lSsKDEmII6zBt5lixZzN36boSj0GNQUkxxVbCgPLqKq8L35Zdfuh6KsQ0FwJyWc73BoxCUhkuXLsUoD2sW+ghLnanQmQVcp/9WrFihu10VLFgdXWXcuHHO8c6cOdN5CBZGtAEGUC7RR1NhhNJiij/XDue4KliGz5tZjb6/8cYbzr4YvlExjk2YMMF5zFSwXJXhpk2bOr9j+K516tTJWd7wodK6XBUs9/oNHy8tj+lJiD/XxN9+aAP8QwJhRoBO7sYvIoUE7CYAJ2RDkdJqK1euLIaVxdkEtuHobCgBzn3YMBQh52c4JXuSP/74Q3f7U9a1Hjh0uwoc59EX1AfHc6sCB2pTXJ3v3Z3F4RAOMX4XzeIx3uH07SquKxZj68/x48fFmGrTUw2l0bUKMabQYnw2pj7VQR87DQuQGD5XMY4j0KuhWOm+Xbt2qfO4awHDn8v1YwznefMao8ALL7wgxvSZnDp1Sp32DWVO4AAPMayD+o4/cb12ODd//vx4cwrGZgquo6u4f8Yx17bhAG86wbueh218zwylNMZubxxMBv5ck7vpR4xO8QMJhDABKlghfHHYtfAl4Br53Jh+iTEQw7rhvPG6HkBUd4gxzSbPPvtsDKXMLGeu9vOnrHku3k2lxHUfVv1BDKua6+5Yt8323Qu5jtv9mKfP7v0x+4KysfXHVUnyFajVsChpxHxcB1NBde0LFANTsArSXaAouYqrsuy63/Cxkt69ewtWBxp+TNomjuP8Z555xlnUZGflOjtP+t+G67ixy7CAOosY05vObWz89NNPMT7jg9k2trE61X3FI/ZDTMX4/z/9/19fHFz75uua3E0/XPvEbRIIZQJUsEL56rBvYUsAFh4oCLBSIVyBMX3iVJgMR2DBsnp3wfJ6Y2pQrT3GNJqGD0AZKB0IZYBl9qZ1xp+yru2gbWP60hnyAQoHLC4QV+uR6zmB3Db8tXS5v9kGWJkSW3+gCIGxMbWloQmMaUpVZHDuV199JS+//LKGYDCmvTR8BJQJKBwI5QDLGPiZsmTJEnPTY8R6KEJWBOEKEHbCWEUnxvSm83ojlIExTeysIq7XDhW4Wgvx2dWihZAZZugIHEOoD3dxHTeUM8Mf0FkE1jsoSca0pir5zgP/2/DFwZ9r4tpvf/vh3i9+JoFQJZAwVDvGfpFAuBMwfGd0CIgJ1LlzZ02JA4Vg8ODBHodmOGI79w8fPlxjDkExe/PNNwWpdXAjMpzctYw/ZZ2VGhtQqMy+IFWP4bCtyh/KQBEItgwcOFBjUhl+UzJ69GiZPn26dgGxurxZV8w+ItYVBDG7sI3xIHYYeCFuFBQt0yLlOjUK5QtTVOa1gPIKwVTuww8/rNtx/YNpQgisSagfYvjN6bv5J67XDucbvnhmNfqOukxr04gRI8RY6ajxu4zFA8oiRmHjA2KEmWPcsGGDWtowxQdeUM6MhQFq2YISHhexek0C3Y+49J3nkIDtBMLMZ4zdJYGwIWBMPcWI0QRHdcMKoEvsDeuMOggbSpNzPHASNgJE6n7jH10d201HYnxGmAdT/ClrKBPOOg0LhW5j5R3qNF9Y8YiVZKb4cnJHiABXMevBqkBXMVeVIayAKa5O7t76M2vWLLO413fEn4JDudm2+ztCKpiCRQfNmjXzWhYLEcDJFFcnd8MqZe7Wd8R5MtsyFMIYx/AB8cDM4wUKFLjjuD/XDie7OrmbKzddK4XTPb5XZpvmu+Fv5txnOrnjvNWrV8cIT4HvmPl9wHVyXZHp6uSOFZeuYihjWj/ic5nizzXxpx9m/XwngXAiQAuW8WtEIYFAEIDDNyJhG0FCtXpYCuA4jCkpd+duFMAUDKaW4MeD6S+UhxUE0zovvviiRvA2++lPWfMcvLds2VKMMAXO6TT4EyE6OKKrGzdX16JB2TaCl2r0dljqIJhKg/+Sq5XHW0eMAKaCKOmwFLr69CDqfM+ePZ3WMJyP44jMDraILm8K/IqM1XSC6TFjNZy5+67eTSsWKjHigd1RV1yv3R0V/W8HLHIYG6xWmJZGtgBMU3bv3t15iqv/FMp9//33ggj1uOb4jmHqERH3DcVWr4fzRD83/LkmgeyHn91mcRIICIEE0AYDUjMrJQEScBJAWhb4XflzE4fzNaZqXP1mnBW6bcRW1oi07fTVwXSVERpB4PCNFChGPCvx1zHdrWm/P2JaEjdyCKZPoXTCDwyO0fDNgQLir4DTvn37dLoMPkTuvkru9UGpMKxUYliY4kWxdO1PbNfOtZynbXDD9B6u43333ReDHaZb27Vrp6dhatnMLelaDxYZYLoU30ukZbJT/LkmgeyHnWNiXSTgD4HgP7L60zuWJYEIIYCbH17+iKdl9t7O96cs6oBFx4gA7q26oO/HSj+84ipQqJBX0KrAQohXKIi/1861z4eNPIKwWEHgZI9FAlDI4Wtn+rPhmGlFxbarwH8rUN8Df65JIPvhOl5uk0AwCVDBCiZttkUCJEACNhKAwzoUS6yORDJxrLzEYgjXVaqYboVTOYUESCC4BOiDFVzebI0Egk4AflaYBsTLXHEW9E64NIipKLM/cZkOdKkq6jdhJcJqScTgMv3QTOUKilavXr3EiOge9ZwIgATigwB9sOKDOtskARIgAZsJIIAtfP2MlXzq1+Yeed3m5lgdCZCADwJUsHwA4mESIAESIAESIAES8JcApwj9JcbyJEACJEACJEACJOCDABUsH4B4mARIgARIgARIgAT8JUAFy19iLE8CJEACJEACJEACPghQwfIBiIdJgARIgARIgARIwF8CVLD8JcbyJEACJEACJEACJOCDABUsH4B4mARIgARIgARIgAT8JUAFy19iLE8CJEACJEACJEACPghQwfIBiIdJgARIgARIgARIwF8CVLD8JcbyJEACJEACJEACJOCDABUsH4B4mARIgARIgARIgAT8JUAFy19iLE8CJEACJEACJEACPghQwfIBiIdJgARIgARIgARIwF8CVLD8JcbyJEACJEACJEACJOCDABUsH4B4mARIgARIgARIgAT8JUAFy19iLE8CJEACJEACJEACPghQwfIBiIdJgARIgARIgARIwF8CVLD8JcbyJEACJEACJEACJOCDABUsH4B4mARIgARIgARIgAT8JUAFy19iLE8CJEACJEACJEACPghQwfIBiIdJgARIgARIgARIwF8CVLD8JcbyJEACJEACJEACJOCDABUsH4B4mARIgARIgARIgAT8JUAFy19iLE8CJEACJEACJEACPghQwfIBiIdJgARIgARIgARIwF8CVLD8JcbyJEACJEACJEACJOCDABUsH4B4mARIgARIgARIgAT8JZDY6gnXr1+XQ4cOSaFChfSUSZMmyaZNm6RatWrSunVrq9VEZLlt27bJggULInJsHBQJkAAJkAAJhBKBZs2aScmSJUOpSx77YsmC9eeff0qePHmkZ8+eWsnYsWPlhRdekEWLFkmbNm1k5MiRHiuPlp1QrtauXRstw+U4SYAESIAESCBeCOBeGy4GDUsWrGHDhkmyZMmkf//+ChQKVuXKleWbb76RcePGycCBA1X5SpAgQbwAD4VGq1SpIkOGDAmFrrAPJEACJEACJBCRBPr16xc247Jkwdq+fbu0bdtWHnnkEdm7d68cPHhQnnrqKUmSJIm0aNFCLl26JIcPHw6bQbOjJEACJEACJEACJBBIApYsWP/++6+Y1qmvv/5a+1O7dm19P3bsmL5D2aKQAAmQAAmQAAmQgDuB3377TU6fPi358+eXrFmzuh+OyM+WFKyyZcvK3LlzBe8ffvihlChRQnLnzq1O72+99ZbkyJFDsmXLFpGAOCgSIAESIAESIIG4Efjll1/k2WeflZMnT6qusG/fPqlZs6ZgoVyGDBniVmmYnGVpirB37946nDp16iikiRMn6mc4ukMr/fjjjyVhQktVhQkWdpMESIAESIAESOBuCMB1qEaNGtKpUyfVHbZs2SInTpyQ9OnTS926deXWrVt3U33In2tJK4I5b8eOHfrCisJy5crpwIYOHao+WdBGKSRAAiRAAiRAAiRgEoCOAOsVjDGmm1Hq1Kll8uTJcvv27bBZDWiOx993r1OEgwYNkqtXr/qsD6EaIFhpSCEBEiABEiABEiABEFi3bp3MmjXrDhhQthDLCsexUC5SxauCNW3aNLlw4YLlcVPBsoyKBUmABEiABEiABCKcgFcFC1HbKSRAAiRAAiRAAiQQFwKIlzlv3jwpXbp0jNMdDodOD7766qsx9kfaB0s+WJE2aI6HBEiABEiABEggsARef/11mTp1qnz00UcCpQpy5coV6dChgy6MwzRhJItXC5ZVHywTDqcITRJ8JwESIAESIAESQDgnZHyBo7sZ0skM07Bs2TJJlChRREPyqmDRByuirzsHRwIkQAIkQAIBJ4C4mdu2bdOQTgw0+j/c9MEK+PeODZAACZAACZBAVBAoWLCg4BVN4tWC5QkCgoquWbNGcxEOHjxYEKG1VKlSkjixX9V4qtrvfdeuXdO4GkmTJo1x7hdffCGbN2/WY4jXVb9+fWf8jRgF+YEESIAESIAESMAWAvCx6t+/v5QvX17q1atnS53hXollzQgZrEeMGCE3b97UMb/55puC15kzZ2TGjBlSpEiRoLKoXr26PPDAAzJ9+nRt9/r169KgQQNZuXKlfkacDVxwXOwVK1YIgptRSIAESIAESIAE7CVw48YN9bPau3evdOvWzd7Kw7g2S6sI58+fr8rVO++8I2vXrnUOd/jw4QKwAwYMcO6Lrw042UO56tOnj+zevVuOHz8uo0ePVmsbVixQSIAESIAESIAE7CXw999/q3Hj1KlTOsOVKVMmexsI49osWbCQ6PnFF1+UHj16CFYAmFKyZElN/gzLERStJEmSmIeC/r506VJ59NFHY0SU7969uy4F7dWrl2BKMVmyZEHvFxskARIgARIggUgkcPbsWZ0OzJMnj3z66afi7rITiWP2Z0yWLFiwBmXPnt1jvTly5JDLly/L0aNHPR4P1k4od40aNbqjuebNm6vy9/vvv99xjDtIgARIgARIgAT8J3DkyBGpWLGi5iaePXs2lSsPCC0pWIULF1Y/Kw/ny5QpU9QyhHgXwRYs+fzjjz+0WQQsw1JQd1m+fLnugiJIIQESIAESIAESuDsCO3fulAoVKkirVq1k7NixOlN0dzVG5tmWpgi7du2qWiocxrEqD/L111/LwoULZfHixdK5c+egBwxLnjy5Oq9DccqSJYvkypVLfvrpJ00g2bRpUzl//ry88cYb8sknn2gySTq5R+YXmKMiARIgARIIHoENGzZI48aN1R2nffv2wWs4DFuypGDB1wqO7l26dBGsJoQ8+eSTqrV27NhRELIh2LJq1SrBtB9CRcByhXcoW6ZFC+bL8ePHS82aNdVPLNj9Y3skQAIkQAIkEEkElixZIu3atdP0Nw0bNoykoQVkLJYULLQMmLVq1ZIdO3YIgpCmT59eQzPE59QbHOvwgjZtipnvCAHNDhw4oKEczGN8JwESIAESIAES8J/Axx9/LL179xbEmqxUqZL/FUThGZYVLLDBtFyZMmX0hc///vsv3kJKEP8KkjJlSipXIXVl2BkSIAESIIFwJIAQTR988IGGYXjooYfCcQjx0mdLTu7o2bp16zRo5+3bt7WjFy5cUAVmwoQJukovXnrvpVE4vyMoKvpIIQESIAESIAES8J8AZoTggz116lTZuHGjULnyj6ElBWv9+vVSpUoVSZs2rcaTQhOIdwE/LJgMQy2QJ/yw0K9z5875R4OlSYAESIAESIAEBNlRWrZsKZs2bRI4tt9///2k4icBS1OEn3/+uWBl3oIFC5zVp0qVSt577z2dLmzTpo2aD9OkSeM8HuiNsmXLysWLFz02g6CikBo1ajiDn7oGSPV4EneSAAmQAAmQAAnIlStX1Lc5UaJEsnr1asH9nuI/AUsKFpIn16lTx2PtcHaDGfHkyZMSTAULuQ+nTZsmOXPm1JWCrp1DdFkEPi1durSkS5fO9RC3SYAESIAESIAEvBD4888/pW7dulKoUCGdGozPDC1euhg2uy0pWHBsxxxs37597wgoBiUHFwCr+YIp6A9CMLz00ks6bTlu3DinMrV161Zd6YD8hPny5Qtmt9gWCZAACZAACYQlAYQ+wn0VCtaYMWPEXDQWloMJgU5bUrBatGghEydOlCZNmggCiyE0w6VLlwS+WSNHjpROnTo5p+KCOSbMDyOa7DPPPKPOd9OnT5eqVava1gVMiQ4cONBnfYi5BW1/yJAhPsuyAAmQAAmQAAmEGoHt27frTFW3bt3ktddeC7XuhWV/LClYjzzyiHzzzTeCvH6IgWFK4sSJVbnCEs74EkRwxxzxqFGj9MsBixbS5tghiPuFKUhfAif/hAktrRfwVRWPkwAJkAAJkEBQCSBKAPyscR9FIFGKPQQsKVhoCkkd4WcF3yYEG4VvE4J5ZsqUyZ6e3EUtUG6gccOpvXXr1jJr1qy7qO2/U+FTVq5cuf92eNkKpu+Zly5wNwmQAAmQAAn4TWDRokU6MzVjxgypV6+e3+fzBO8ELCtYZhWwGOEVioKUPj///LP6iv3www+ahDoU+8k+kQAJkAAJkEB8E/joo4+kf//+8tVXX8mjjz4a392JuPb9VrBCnUCKFCnUOc/sJ/yj4C9WrFgxcxffSYAESIAESCCqCSCH8KRJkzSIeOHChaOaRaAGH/GOQ7169ZLixYsHih/rJQESIAESIIGwIYBsLPBVnj17tkZnp3IVuEsXcRYsd1RYAUnrlTsVfiYBEiABEog2AgjCjcDgyHaC6OwZMmSINgRBHa9XBQvTagggalXSp09vtWhQy2FlBIUESIAESIAEopnA5cuXpVGjRgI3mlWrVknKlCmjGUdQxu5VwULgUH+SJfujjNk5MrS7ZcsWOX78uJw6dUquXr2qORPz5s0rSKeTOnVqO5tjXSRAAiRAAiQQVgRwb0Q2FrjLTJkyRRBiiRJ4Al4pv//++87Ezjdv3pR+/frp6sHnnntOkz4iVxECjSKi+oABAwLfUw8tfPrppxoI9NChQx6OioaSQCBUBEdlRFqPiLiTBEiABEggggkcPHhQo7PDeoU4V7wXBu9ie1Wwnn76aWcv3n77bYEj3Jo1a2Jovq1atZJ77rlHJkyYIH369HGWD8YGHPQQEA39HD9+vMBihflkaObnz5/XmF0Iitq9e3fZv3+/jBgxIhjdYhskQAIkQAIkEBIEkDYOsa2w2KtHjx4h0ado6oRXBcsVAqK4Q/v1ZFZ8+eWXBTn/EA7h/vvvdz0toNtQqhDSf/To0Xe0gyComOIsX768ThN27dpVEG2emvsdqLiDBEiABEggAgl8++238uSTTwpmo+DYTgk+AUthGmAdmjt3rsfeIbAnFJps2bJ5PB6onTB7Vq9e3Wf11apVkxMnTgjKU0iABEiABEgg0gnMmzdPU9shqwmVq/i72pYULOT2gyM5kjpv3rxZl3ju3LlTA3q2bdtWkHQ5SZIkQR0F8gTiS+RLEP4/efLkkjt3bl9FeZwESIAESIAEwprAhx9+KF26dJGvv/5aateuHdZjCffOW5oirF+/vsyfP199nhBa31UwTQfHuWDL888/r6sikBcRvmCwsmXMmFEVPdMHa8mSJZoCYOjQoR6nN4PdZ7ZHAiRAAiRAAoEigAVn06ZNk++++05zBQeqHdZrjYAlBQtVwYpVt25dgeVq9+7dkiNHDilVqpQ6uVtryt5SlSpV0qTTiEjbu3dvQXRadylTpozMmTNH56Hdj/EzCZAACZAACUQCgVu3bsmLL74oyMG7ceNGyZ49eyQMK+zHYFnBwkgRmAyxpfAKBYEj+/Lly+Wvv/6SPXv2aM5BOLLDHwyvUA1+Ggrs2AcSIAESIIHwJ/Dvv/+qm87Zs2c1dBLve6FzTf1SsH777TcN1QCHcSSK/OWXX9SK5Wl1YTCHmCZNmpBR+oI5brZFAiRAAiQQvQQuXrwoTzzxhM4krVy5UqO0Ry+N0Bu5JSd3dBuBRosWLapmSPhcXb9+Xd58801VbHbt2hV6I2OPSIAESIAESCBCCWB1/GOPPSYFChSQRYsWUbkKwetsScGCgzsCdSKW1Nq1a53DGD58uNy4cSPeIrk7O8INEiABEiABEogSAvv27ZMKFSqoXzRS3yRKlChKRh5ew7SkYCEGFhzoEAk2a9aszhGWLFlSsCQUpkkoWhQSIAESIAESIIHAEfjpp58Ei7yQpQRGDkroErCkYCGRsrdVCVhNiCzdR48eDd1RsmckQAIkQAIkEOYEYMxAbCtEZ0eGEkpoE7CkYCEPIQJ2ehKYJ5MlS8ZAnp7gcB8JkAAJkAAJ2EAA+XcR1Puzzz6Tp556yoYaWUWgCVhaRQhNuVy5cprbD0FHIYgSu3DhQlm8eLF07tyZc8CBvlKsnwRIgARIICoJjBkzRqcDYcFC/ElKeBCwpGDB1wqO7gi/j9WEECSRTJgwoXTs2FFDNoTHcNlLEiABEiABEggfArjnImD2+vXrJX/+/OHTcfZULClYiHeFaUKsXEBqmkOHDmkQzyJFimgKmsmTJ8vLL7+sCheZkgAJkAAJkAAJ3B2BmzdvqgFj69atGp3ddYHZ3dXMs4NFwJIPFuJeYf4XSZORfqZFixaCZMtwcP/xxx/V2e7ChQvB6jPbIQESIAESIIGIJXD16lVp3LixGjPWrVsXY/V+xA46Agfm1YIFR7ohQ4bokI8dOyZJkyZVnytXBsh/dOTIEcmUKZPce++9roe4TQIkQAIkQAIk4CcBGCvg65wlSxZZsWKFLiLzswoWDxECXhWsJk2aaFqcf/75R65cuSJp06aVYsWKxeg2fLCqVasmbdu2jbGfH0iABEiABEiABPwjgJBICMNQsWJFGT9+PN1u/MMXcqW9KlhJkiSRiRMnaocxRYg4WFgiSiEBEiABEiABErCXwJ49e1S5evrpp7lwzF608VabVwXLtUc9e/Z0/chtEiABEiABEiABmwjAlxlJmwcMGCAvvfSSTbWymvgmYEnBwhQhkjvHJhkyZIjtMI+RAAmQAAmQAAm4EVi+fLm0bt1aJk2aJM2aNXM7yo/hTMCSgtWpUyeZNWtWrON0OByxHudBEiABEiABEiCB/wh8+umnmlMQQburVq363wFuRQQBSwpW+/btpUqVKjEGfP78edm+fbtGcp86dWqMY/xAAiRAAiRAAiTgncDo0aNl5MiRsmrVKkEwb0rkEbCkYEG5clewTBTDhg0TfFEQ2Z1CAiRAAiRAAiTgnQBme3r16qXGiY0bN0revHm9F+aRsCZgKdBobCPEnPGmTZvkzJkzsRXjMRIgARIgARKIagKIzo6wRqtXr9bo7FSuIvvrYMmCFRsCpNGB0AcrNko8RgIkQAIkEM0EEFMSBolr164JorOnSZMmmnFExdgtKViI6m4qUiYVrCqE1QrHypYtK5kzZzYPBfUdit2WLVsEAdpOnTolSDGAoKh4MkC/UqdOHdT+sDESIAESIAEScCVw7tw5qVevnuTKlUunBpEZhRL5BCwpWDBnzps3LwaNBAkSqAaO2B1w1IsPwQqMgQMHar4mT+2nS5dO+wYnffSXQgIkQAIkQALBJHD06FHN3Vu9enUZO3Yso7MHE348t2VJwfroo48Er1ASJJ9u166dIOotUgrAYoVYXIkTJxascDx58qR88cUXugR2//79MmLEiFDqPvtCAiRAAiQQ4QR27dql0dnxkP/WW29F+Gg5PHcClhQs86SdO3eqQ/vu3bslT548UqpUKSlXrly8aORQqrp166YrGM3+me+wXKF/5cuX12nCrl27yjvvvEMrlgmI7yRAAiRAAgElgBWCjRo1kqFDh0qHDh0C2hYrD00ClhQs+DlhKm7QoEE6invuuUcuXbokt27d0vD+sCalSpUqqCM8ePCg9O3b12ebSEZ94sQJQfl8+fL5LM8CJEACJEACJHA3BL788ktdLfjJJ5+oknU3dfHc8CVgKUzD9OnTZfjw4Zr8+e+//xY47GElBAKkwZqFZafBllq1at3hF+apDzNmzJDkyZNL7ty5PR3mPhIgARIgARKwjQCUqmeffVad2WHBokQvAUsWrKVLl0qfPn2kY8eOTlKJEiUSOO1BgYGVCKv3UqRI4Twe6I3nn39e6tSpIzt27JBWrVqpD1bGjBklSZIkTh+sJUuWyFdffaUmWvhmUUiABEiABEggUATg6/v+++9rnKtixYoFqhnWGyYELGkdFy9elBw5cngcEhzLoVzBsTx79uweywRiZ6VKlVS5Qubx3r17y+3bt+9opkyZMjJnzhxGmb+DDHeQAAmQAAnYRQBuNN27dxckbobvFWdM7CIb3vVYUrAaN24sY8aMkYYNG0qmTJlijBhThw8//HBQlSuzA3Bkxxf6r7/+kj179qhfGMIxZMuWTV/p06c3i/KdBEiABEiABGwncOPGDXWTwWr1DRs23HGPtL1BVhg2BCwpWHAOv3LlimTNmlVXDT7++ONqMUI02vXr16sTH1b0mYLpu4ceesj8GPB3RMRFUFFTTp8+LZMmTdKVG3DIp5AACZAACZCA3QRwX2zSpImuUF+zZg0DW9sNOMzrs6RgYUUEIrcjWvvvv/8ukydPdg4bShdyEeJlSu3atYOqYJntmu9//PGHThvii08Fy6TCdxIgARIgAbsIIJNJ3bp1JX/+/IKFYPD/pZCAKwFLCta4ceMEr1ASWKzgG+ZJsMIRUqNGDeeXft++fZ6Kch8JkAAJkAAJ+EXg8OHDUrNmTVWw3nvvPcZY9Ite9BS2pGCFIo4iRYrItGnTJGfOnPpFd+3j2bNnBekJSpcuLQg6SiEBEiABEiABOwj8+uuvuoIdC6z69etnR5WsI0IJWFKw4MSHcAeff/654MvlSbZt2+Zpd8D2TZ06VRUrfMlhsYKFzVSmtm7dqmlyhg0bxuCiAbsCrJgESIAEoovAd999pz5XyL+LWFcUEoiNgCUFa8GCBRprChYhzDeHylxzy5YtpUKFCvLMM8+ozxfmwatWrRrbeP06hlUhWKXoS5D3kCsWfVHicRIgARIIXwIwMGABF+4zDRo0CN+BsOdBI2BJwUIqnBYtWsjcuXOD1jGrDeXKlUuDuo0aNcpptm3WrJnV02Mth9WI3ix2riciun2wUwW5ts9tEiABEiCBwBHAqnRMB2LBFx7qKSRghYAlBStZsmRy3333WakvXsokTJhQXnvtNXVqb926tcyaNcuWflSsWFHw8iX8h/NFiMdJgARIIDwJDBkyRNPErV27VuD7SyEBqwQs5SKEWRTWq71791qtN17KlSxZUn7++WfB1CFWGUIxpJAACZAACZCAvwSQHaRz584yc+ZMjc5O5cpfgixvyYKFnH9IWvnggw9KwYIFpUCBAneQQ96/UBDkQ0TUeQoJkAAJkAAJxIUA4j62adNGV6MjOvu9994bl2p4TpQTsKRgITs45qBTpkypztzIOxgucuTIEU2hw8Sb4XLF2E8SIAESiD8CSL0GgwJmQFavXq33vfjrDVsOZwKWFKz58+erYx9SAYTKCkKr0Hv16iXoP5JxUkiABEiABEjAGwEsbMKMTdGiRQWGhcSJLd0ivVXH/VFOwNK3J1GiRFK5cuWwU65wbbH6kdarKP+Wc/gkQAIk4IPAwYMHpVatWvLEE0/Iu+++y+jsPnjxsG8Clpzcn3zySVm0aJEg91K4SdOmTaV///4h3+0LFy6og/7x48dDvq/sIAmQAAlEEgEEysaK8U6dOsno0aOpXEXSxY3HsViyYGXPnl1NpTCblitXThB7yt10Gl+O5Zj627Jli0AxOXXqlFy9elXSpk0refPm1ZWEqVOnjke8vpvGfD9WqiCYK7giaCk4T5kyRQoVKuS7ApYgARIgARKIMwG4viB2Iu5hTz/9dJzr4Ykk4E7AkoKFCLbnzp0TTBVCmcHLXeJDwfr0009l4MCBcujQIffu6GekzkFKg/bt24fkEwmWAderV0/zKR47dkwyZMggN2/elLFjx0qVKlXUogXllkICJEACJGA/ATzYduzYUWMnwveKQgJ2ErCkYCHPH16hJIgu365dO33iGD9+vFqsoKDAsoZVjrAEffHFF9K9e3dBypsRI0aEUve1L1BcYcGCogjlFYL+9+jRQ5cHDx06VD788EPdzz8kQAIkQAL2EcB9Y8CAAbJs2TJ55JFH7KuYNZHA/whYUrBCkRb+Obp166bz5e79g+UqT548Ur58eZ0m7Nq1q7zzzjshZ8XCEmA44ZvKles4WrVqJW3btnXdxW0SIAESIAEbCGDmA6sEkbyZrhg2AGUVHgl4VbD69Okj2bJlEygnM2bMkJ9++sljBeZOTGsFU7Dio2/fvj6brFatmpw4cUJQPl++fD7LB7PArVu37vBlM9uHJQvHKSRAAiRAAvYQgFvGiy++qJHZN27cKDly5LCnYtZCAh4IeFWwli5dqpo9FCwoV1hFGJsEW8HCctp58+apD1Ns/YJymDx5csmdO3dsxeLlGFatTJgwQXr27HlH++BtJQ/iHSdyBwmQAAmQwB0Erl27pmnUEOtq/fr1cs8999xRhjtIwE4CXhWsnTt3OtuB8hRsBcrZuJcN5EeEU+KOHTsE02lYNZgxY0aN1WX6YCF9z1dffSXwZXJf9eil2qDuxvQg4q288sor6iNm5k5E7isoXps2bQpqf9gYCZAACUQigUuXLknDhg0lTZo0smrVKkFKNQoJBJqAVwUr0A3fbf2VKlVS5eqll16S3r17C0y/7lKmTBmZM2eOII5XKAqi4n/99dfqrI/p2Iceekid27EfjpehNqUZigzZJxIgARKIjQAWPNWuXVtKlSqlKd9C8WE7tv7zWPgSCFsFC8jhyL58+XJdibdnzx7NOZggQQL1HYPCkj59+pC/Mvfdd5+OAT5i+/btk8yZM0vJkiUlYULPMWDxQ4GYLc8++6xH5/iQHzA7SAIkQAJBIoAV5HAnad68uS50ClKzbIYElEBYK1jmNYTZt2zZsubHsHx/4IEHBC9fgtyKWICA+F7Dhw/XpKRQKikkQAIkQAL/EUC8xvr16+vvJVacU0gg2AQ8m0mC3Qu2Z5lA9erVddHB4MGD5Y033lCz98qVKy2fz4IkQAIkEOkEvvnmG7VcvffeexrOJ9LHy/GFJgG/LFhXrlzRaThEG3eX+++/330XPweQAEzeyLOIIKVYdpwlSxZ1lK9QoUIAW2XVJEACJBDaBObOnSsvv/yy+t/WrFkztDvL3kU0AUsWrL///ltNrcjxh7ghCHng/opoSiE6OPhpIZr93r17dSUlViXiB8VXzLIQHQ67RQIkQAJ3RQCr3bEqG4uHqFzdFUqebAMBSxast99+Wx2x8cUtXrw4l7jaAN7OKrDqEE9sCF2BHxjkN6xcubKmgShSpIidTbEuEiABEghJAv3799ecgohxVaBAgZDsIzsVXQQsKVgbNmzQmzfmsymhSwABVV977TVB6Ao4wSOUxRNPPKG+WlYc6EN3ZOwZCZAACXgmgIwXSNgMy/33338vWbNm9VyQe0kgyAQsTREigCcDswX5ytxFc6lTpxbk2jpw4IAGX0X8F1i4/vjjj7uolaeSAAmQQGgRuHr1qjRp0kR/62C5onIVWtcn2ntjScGCbw8cBy9fvhztvMJq/BkyZJBRo0apjxYCsRYuXFgtXGfOnAmrcbCzJEACJOBO4MKFC+pnBV/UFStWSLp06dyL8DMJxCsBS1OEmTJl0icDzGtXqVJFEBzTPfbSmDFj4nUgbNw7AawwROodxM9CaAdcR8SFwYs/St658QgJkEBoEoA1HkGXH330Uf1tS5QoUWh2lL2KagKWFKzFixfLn3/+qdHFv/vuO4/AqGB5xBJSOxFKA2EdfvvtN3n99dc1sCmULvhspUyZMqT6ys6QAAmQgCcC+P3CCsE2bdrIkCFDPBXhPhIICQKWpgjHjRsnJ06ciPUVEqNhJywRKFiwoCxcuFCTniLxKRzgYeG6fv26pfNZiARIgATig8CmTZvkscceU1cHKlfxcQXYpj8ELClY/lTIsuFDoESJEhovBsrW7NmzdeoQFi6syqGQAAmQQCgRQGyrOnXqyAcffKCLdkKpb+wLCXgi4HWKEFNHSJjctWtXmTFjhs/glYi/RAlPAuXLlxeswEF6iX79+snQoUPV9I7VOe6+duE5QvaaBEggnAnMnDlTA4guWLBAqlWrFs5DYd+jiIBXBWvp0qVSqFAhVbAQX2TRokWxYqGCFSuesDhYo0YNwQs/YgMGDBDkOxw2bJg6k4bFANhJEiCBiCOA+IsjRozQB8CHH3444sbHAUUugQQOQyJ3eMEZmZn/b+PGjcFpMAitIKwDLJeI4o9Vo8OHD9fApUFomk2QAAmQgODW1Lt3b/UXRUJ7BkvmlwIEMMsCCQcfPPpg6aXiH3cCiC3Ttm1bjaHVunVrzXUI69bmzZvdi/IzCZAACdhK4ObNm5pnFW4LiM5O5cpWvKwsSASoYAUJdLg2kzhxYg3jsH//fl0aXb9+fWnWrJns3LkzXIfEfpMACYQwgX/++UcaNmwox44dk3Xr1qkFPYS7y66RgFcCVLC8ouEBVwLIc9irVy85dOiQFC1aVJdKw8KFdDwUEiABErCDwPnz56V69eoal2/58uWSNm1aO6plHSQQLwSoYMUL9vBtFHkO4QB/8OBByZw5s5QuXVotXMePHw/fQbHnJEAC8U4AFquKFStKyZIl5bPPPpNkyZLFe5/YARK4GwJhq2Bdu3bNY2DML774Qp3g+vfvL19++aU6St4NIJ7rmcA999wjI0eO1KjwKFGkSBG1cDHPoWde3EsCJOCdwO7duwXhYp588kkZP368Zg3xXppHSCA8CPilYCFFwcSJE/VGCgUHUXXhjBgfAjNyhw4dnE0jCnmtWrWkUaNGGscJsZyeeOIJfSK6cuWKsxw37CWAFYb4QdyxY4emU8qfP79auC5dumRvQ6yNBEggIgnAib1y5cqCh2JYxykkECkELCtYWBoJ35sXX3xRRo0apdajN998U8qWLSu7du2Kdx6I14SlvAiQiqchTFmNHj1ap7JcFbF472iEdiBXrlwyffp0XWUIB/i8efOqhQsOqxQSIAES8EQA8RaxcOajjz6SF154wVMR7iOBsCVgScGaP3++Bnp75513ZO3atc7BIjbSjRs3QuKpA/+oyKwORevBBx/UKPTdu3eXvn37ahwVWNwogSdQoEABDVS6evVqwQuKFixczHMYePZsgQTCicC0adM0FMPnn38uyBpBIYFII2BJwZo7d65arnr06CFZs2Z1MoAz4ocffqiWIyha8SlJkiTR6UH3PjRv3lyVwN9//939ED8HkEDx4sUFq4Dw44nvD6YOYeFinsMAQmfVJBAmBOC/+frrr8u3336r04Nh0m12kwT8ImBJwcJ0W/bs2T1WnCNHDrl8+bIcPXrU4/FA7jx9+rT88ccf2gRiM23btu2O5nCTh6CflOATgFXxu+++k48//lgtWYULF1YLFxMIBP9asEUSiG8C+L/Hg/rkyZMFmS/wIEYhgUglYEnBwk0RaVM8yZQpU3Q5be7cuT0dDtg+xGVasWKFKk6wqmFZL14LFy7UNhFP5eWXX5bOnTtLixYtBOEFKPFH4PHHH9dFEZjCHTRokC7FNpXf+OsVWyYBEggWAcxytGnTRh+4oFzlyZMnWE2zHRKIFwKJrbTatWtXKVeunC6jhUMi5Ouvv1ZlZvHixarEJEqUyEpVtpVZtWqVYNrvl19+UcsV3mGlMi1aR44cUYtJzZo1dRrTtoZZ0V0RgK8FVnrOmjVLE4lnypRJ8xw+9thjd1UvTyYBEghdAn///bc0bdpUXQTgx8sH3tC9VuyZfQQsJ3tGfKkuXbpo+gKzeeSr69ixoyDbOSxKoSAwQSdIkECweu3kyZNByWEVicmeg3EtEeIDFlAk7SxYsKC+Q5GnkAAJRA6Bs2fPSt26dfW3GH6YSZMmjZzBcSRBJxCRyZ6RG2rfvn26DB9Oy7BgwUo0YcKEkFGuzCuN+FzoF8IFxLfzvdknvt9JAHkOO3XqpOl26tSpo3HLYOFCTC0KCZBA+BPAPQIPoI888ojMnj2bylX4X1KOwA8ClnywYAlCkFH8syA3VIkSJQQ+VzD7Yj9SHMRHvKN58+ZJlSpV1J8HqxkhDRo00H/mnj176lQUogOfOnXKDyQsGmwCSInx6quvaswyfLdwTZ9++mnmOQz2hWB7JGAjATwoQbl65plnZOzYsTqzYGP1rIoEQp6AJQULSX4LFSrk9YUgk2nSpNEEwFDGgiGYsmzZsqVcuHBBn4rgzI7kw4gKjACj+/fv16lL+GThGCX0CcAvA8FrkUA6W7ZsmucQFi4o8BQSIIHwIbB+/XqpWrWqDBw4UFOXhU/P2VMSsI+AJQUL8UqyZMmi+eYQKR0R07FaD1YGOLdDoRkzZoxauGB9CEZQz3Hjxmn727dv19Vp7777rnz66afqE4YAo/ny5ZNu3brJG2+8IYsWLVJrm33YWFMgCSDPIYLawjqKaURkEICF688//wxks6ybBEjABgJ4+IVLySeffCLPP/+8DTWyChIITwKWFCw4JiI6+q+//qqR0mvUqKGRd6HQYJUe4hzBAR5xqA4dOiQ///xzwGls3bpV/4nNhmC9gtM9fHlcBSvW4Ph++PBh193cDgMCyHMIRRpTDQi7gSjxsHBdvHgxDHrPLpJA9BHAopX27dvLkiVL1Kcy+ghwxCTwHwFLChZiltSrV89jhvPnnnvOmT4nQ4YM8vDDD8tPP/30XwsB2kLsK1jSTLn33nvlq6++Ut8wcx/ezTKYxoyLXL16VXy9GJ08LmStn4NrN3XqVF1gsXfvXl2NBAsXfAApJEACoUFg6NChmjZtzZo1UrFixdDoFHtBAvFIwJKChSX0c+bM8dhNPKlg+tAURHR/4IEHzI8Be0cQUSQIheJ38OBBbad27dpy//336/aZM2fklVde0RQ/rVq1Uh8xfzuDm3rKlCl9vrBqEVHlKYElAAsWFjYgxyGspviewcIVjCnpwI6MtZNA+BK4ffu2zmBgRgMP45jSp5AACYhYUrCwdB7TfnAqR7R0rMqDf8yIESM0zxz2I5kvfJ4wlQPnxkDLiy++KEjwjDQ+npQbhGjAyhVMZ5orDP3t07PPPiv48fD1QjoYTGdRgkMA6TVgrUSQ2wULFmieQyjDiKtFIQESCB4B/O7jAXbz5s2yYcMG5wNu8HrAlkggdAlYDjQKywGmA12nZZBgGYoOAo3CTwbTg2+//bYm8QzWkE3lB87QrnLlyhU5d+5cUP7hGWjUlXzwt5EwFsHnoNwPHjxYkOAbwWYpJEACgSPw119/SePGjQX3ATzopEqVKnCNsWYS+B+BcAo0alnBwtjwD4WUNLAOZcyYUdPnmL5N8FOCMzmm1KJNqGCFxhX//PPP5a233lLlCv4gmD6mkAAJ2E8AK3qxoAiLn2A9hpJFIYFgEAgnBcvSFKEJDbGuKlWqpFYrWAlM5QrHU6RIEZLKFYKjYvUjJfIJ4GkaYTsQtw2hOjB1i7xnFBIgAfsIYKU4HioRkmfGjBlUruxDy5oijIBfCham3RC4E0qL+ytUueBmC58dSnQQwNRgmzZtZPfu3dKuXTuNIl29enX58ccfowMAR0kCASSAGQysEEQoBsQe5FR8AGGz6rAnENNxyctw4HfVokULWbZsmU4DeiqG6cFQFPS7WLFiodg19imABOCT98ILL6iShUUOCHyItEmILM3vQwDBs+qIJQBrcLNmzVSxQtxBCgmQQOwELClYcFxfvny5hj2ANQjTgeEiTZs2DZeusp8BIIA8hz169FBlCxkHsMIVviPw1cqfP38AWmSVJBB5BJC5o0OHDjolSN/GyLu+HFFgCFhSsLD8FikPsFow1ASWsy1btmi4BoSPgLM9ElLnzZtXypYtK8hvRyEBrHBC2iRkHBg+fLiUKVNGrbL9+/eXnDlzEhAJkIAXAhMnTtQMCgiNAr9GCgmQgDUClhQsrBgMRasVAtthygdOl54kXbp0MnLkSPUXoK+AJ0LRty99+vSqYMGqBcssgiLi4QE5NjNnzhx9QDhiEoiFwKBBgwTpb9atW6crBmMpykMkQAJuBCw5ucOPae7cuXL58mW30+Pv4+zZs9W/Bg6XX3/9tezbt0/Onj2reeqgcCGiMHxwsJqsd+/e8ddRthySBKBMffDBBxpyBLkNESUeFi7mOQzJy8VOBZkA4gsixiECS3///fdUroLMn81FBgFLFqxMmTIJcv/hJoSluYha7m4RGjNmTFCJjB8/XiPHw6/GXWC5ypMnjzo1Y5qwa9eugtx17n12P4+fo48Apgc/+eQT6du3r2C6EOl3sPIUU4kMnBh93weOWDT1VOvWreXEiROyfv16QY5ZCgmQgP8ELFmwkJIEgeUSJkyoOeDmz5+vOeEQ3d18+d/03Z2B/INYfu9LqlWrpj8UZr5CX+V5PDoJwOEdT+tIVAufQ/jwwcLFPIfR+X2I1lFjlgI5Xf/9919ZtWoVlato/SJw3LYQsKRgIaEunmZie9nSGz8qqVWrlip3vk5BILzkyZNL7ty5fRXlcRLQEA7IcYkk5osWLdKVhrBwMc8hvxyRTgCLhCpXrqzpxfBQHY1ZOSL9GnN8wSVgaYrQ7BLiYSHYKFbu4XXr1i3NTYg5eiRGDqbAMRnL7ZEDEclGYXGAMz5SNiAn3cmTJ/UmiZUvSJvinqswmH1lW+H7+3F7AAAqpElEQVRHoFy5cmrNWr16teY5xMpDOMU/+eSTnGoOv8vJHvsgcODAAcFDK8LajBgxwkdpHiYBErBCwFIuQjj+PvHEEzof763S+Ag0+vvvv8tLL70kK1euFDhluguW4vfs2VNviu7H7PzMXIR20gzNur744gtdqo7vORT2+vXrh2ZH2SsS8JPA1q1bpW7duroYCIuCKCQQygTCKRehJQsWnt4RawrKDBI+Yxsr9PA+a9asWBWvQF4oOLIjACr6tGfPHrl06ZJaF7JlyyZ4YUk+hQTsIIBI8HjImDNnjgYuHTx4sAwbNkwDl9pRP+sggfgg8O233+oD6NixYwWO7RQSIAH7CFhSsPCE06BBA0HKkU2bNun2K6+8or0oVKiQrtCDz0p8CZJQY7UghQQCSQCrUDEdjWnCqVOnapgQrDocMmQIAzAGEjzrDggBLOrAQzNC3mB6kEICJGAvAUtO7rAMmRF8H3zwQTlz5ozGnUJXnnvuOYGf0z///GNvz1gbCYQoAfjzIW0IYq/BqtW4cWPNdbh9+/YQ7TG7RQIxCWDhEsLXYAaAylVMNvxEAnYRsKRgYQm7efNAGhpMvyGQJwQ3G/ilIIwDhQSiiQDyHHbr1k0QAgQWVIQEgYULiheFBEKVAALqjho1SkPu0PIfqleJ/YoEApYULDyhz5w5Uzp37qxjfuyxxwQBPhEz6LXXXtPgjAyDEAlfB44hLgQQkBSOl1C08H+AmxYsXEePHo1LdTyHBAJCAKu+8b1ECAas/C5YsGBA2mGlJEAC/0/AkoKFpbvI1YaUNJA333xTLVZ4Ysf8val4/X+V/EsC0UkAiyqwwhAWLMQQeuihh9TCdfr06egEwlGHDAEEDsXv+N69e3VREmYhKCRAAoElYEnBQhcQA2j37t3aG/hhYfvzzz/Xp3TT4T2wXWXtJBAeBJDn8P3335ddu3Zp3DhYCmDhunDhQngMgL2MKAIIs1OzZk0d0zfffMPV1RF1dTmYUCZgWcHCIOBvdezYMZ0KwT8tntAReJRpaEL5ErNv8UUgR44cMmXKFA1nggTkWHEICxf+ZygkEAwCyL4Blw7kkV24cKFmtQhGu2yDBEhAxJKCdfXqVWnSpIkmv82VK5fky5fvjhdhkgAJeCaA/xfEz1q3bp38+OOPqmjBwsU8h555ca89BDBVXb58eQ2rA0U/UaJE9lTMWkiABCwRsBQHC3F+kJ+tZcuWmmAZTr0UEiAB/wjA4ot4cZs3b5a+ffvqSq633npL42kxlZN/LFk6dgI//fSTZhvo37+/dOnSJfbCPEoCJBAQApYUrA0bNkijRo1k+vTpAekEKyWBaCKAVYaIoI1VuPDNeuedd2TQoEHSokULSZjQklE5mnBxrH4SWLFihYYLGT9+vH6n/DydxUmABGwiYOnXHJHS4T9CIQESsI9A1apVdbn8u+++q0pW8eLF1cJlXwusKdoIYFU3YrEhSjsUdgoJkED8EbCkYCFH1bx58+TUqVPx11O2TAIRSgDR4Ldt2yavv/669OrVS8qVK6cWrggdLocVIALvvfee5slcuXKlPP744wFqhdWSAAlYJeB1ihA/9q7pby5fviwlSpRQp8mcOXNqUmXXRsaMGeP6kdskQAJ+EECeQ/g4mnkOkYIqb968uurQTFPlR3UsGmUE4NMHqxXcObCogkICJBD/BBIYaW4cnrpRuHBhQSgGq4LlwNEqFSpU0KGb6YOilQPHbR+B69evy4QJE2T48OFSpkwZ9dHCAw6FBFwJ3Lx5U6OzwwKKQNBZsmRxPcxtEog4AvBbhWDxXaiLVwuWGVQ0tgEgOnDy5MljK8JjJEACcSCQNGlSQQBfpDaBdRhTPjVq1JABAwYwxUkceEbiKQifA4sn4qohBEi6dOkicZgcEwmELQFLPlgYHf6BEVPl9u3bOlhEpYbjO56yb9y4EbYA2HESCGUCSLmD6XoE88WU4SOPPCLt27eXI0eOhHK32bcAEzh//rwq3VDEYbmichVg4KyeBOJAwJKCtX79eqlSpYqkTZvWGRwR/9h4eurdu7c+ZcehbZ5CAiRgkQBuoDCJI3hk6tSppVixYmrh4sITiwAjqNjx48elUqVKmklj/vz5kixZsggaHYdCApFDwJKChZyDSBSKJ6UUKVLo6BFsFKtWJk6cKJ9++qn89ddfkUOFIyGBECWQKVMmnTJEnkMsQilUqJBauJjnMEQvmM3d2rNnj84kNGvWTH97GTfNZsCsjgRsJGBJwULk6ZIlS3psFk9S8JM/efKkx+PcSQIkYD8B5DmcPHmy/Pzzz3L48GGdroeFi3kO7WcdKjUizVLlypU1C8DAgQNDpVvsBwmQgBcCXp3cXctjFdPUqVP1H9v9iWnatGmSJEkSyZMnj+spQduGcrdlyxaB2RzTJXD8xFQm/FUQMRvTKRQSiFQC8INEcMmdO3dqVHh87+Gz1alTJy5AiaCLvmzZMmnTpo0q1ZhNoJAACYQ+AUsKFiICYyoQCZ/hYIun50uXLgl8s0aOHKk/5lCygi2YmsST3KFDhzw2Db8V9A99RpwhCglEKoGiRYvKF198oXkOoWCNGjVKkOfw2WefFeY5DO+rjt+5Hj16yMKFCwXR/ykkQALhQcCSgoWVS9988400b95cf8TNoeGHG0/KyKUWbMFTe7t27eTpp58W5NzCk3uGDBn0ZoIVNpiyxA2ne/fusn//fhkxYkSwu8j2SCDoBGC1XbVqlaxduzZGnsOnnnqKeQ6DfjXuvkEoyqNHj9Zryjhod8+TNZBAMAl4DTTqrRNHjx6VHTt26LLgggULCpxu40MqVqyoU4D48YlN8NTXtWtXnUIMlBWLgUZjuwI8Fp8Eli5dKm+88YYgIOXgwYOlYcOG8dkdtm2RAFwfevbsqbkpkbwZD5AUEiAB0QdHcAiHQKOWnNxdL2quXLmkXr16AgUnvpQr9AdxgapXr+7aNY/b1apVE0SZR3kKCUQbgfr168vWrVulf//+GlLFtHBFG4dwGi/iCj7zzDNqhUR2CCpX4XT12FcS+I+A3wrWf6fG71atWrU0AbWvXsyYMUOdfXPnzu2rKI+TQEQSgOUWfpQI7YApffgkIq7d999/H5HjDedB/f3334Lk33BxwDRv5syZw3k47DsJRDUBSz5YoUjo+eeflzp16uh0ZatWrfQpL2PGjLqi0fTBWrJkiXz11VeaMJeOvqF4FdmnYBJIlCiRIIk0VqN99NFHGtuuVKlS8vbbb3sNwxLM/kV7W+fOndPZgfvvv1/wYIhgzhQSIIHwJRC2Chbib8EX7KWXXtKpDzOFj+ulQHiJOXPmaMR51/3cJoFoJoAbd5cuXQQPKWPHjtUch5hux4pcBC6lBJ8AfFthlUfOyffff58LEoJ/CdgiCdhOIGynCEECsbeWL18uFy9elE2bNsnKlSt1tSOmQhDZGgFSkc6HQgIkcCcB5Dns06eP+icWKFBAHn30UbVwMc/hnawCuQcxzJDntXXr1vLBBx9QuQokbNZNAkEk4JcFC1GiEf8KK5LcBWbt+JI0adLoisL4ap/tkkA4E0C8OEwTYrXtsGHDNM8hnKz79esnWbJkCeehhXzf4cTeqFEjdWPo0KFDyPeXHSQBErBOwJIFC46XWI2ECOkIMgqHcfeX9SYDWxLLm2HNevfddzUOFlbkUEiABHwTwKpghD3ZvXu3XL9+XacLYeGCTyPFfgLwEW3QoIFGZ6dyZT9f1kgC8U3AkgULT7eYinvllVekePHizoTP8d35efPmaZBRWNWwMurll1/WHyw4tptSunRp+fLLL/kkbgLhOwn4IJA9e3Z1gn/ttdc0Gny+fPk0YG+3bt0E1mLK3RP45JNPBHwRDBn+pBQSIIHII2BJwdqwYYM6xL733nshQwA/TC1bthSkCEmePLl07txZfa6w9BxP4XgyRJBFRHDHsQULFoRM39kREggHAshzOHPmTA3vgDha+AyLFhaW4H+OEjcCmIYdN26crFmzRh566KG4VcKzSIAEQp6ApSlChD9IkSJFSA0GP1BIk7N9+3bnlCBydnXs2FGftvHUjSduRLFetGiRYJqTQgIk4D+BIkWKyOeff64hT77++mvB/xbCPHD63T+WcF/ALMC0adMEvldUrvzjx9IkEG4ELClYCFI4d+5cuXz5csiMD9GpXdN+tG3bVlffIDaWq8CBFD9shw8fdt3NbRIgAT8JIOwJVurOmjVL4zQhpAMsXJ5CpPhZdcQXh08b4vX98MMPqlzBh5VCAiQQ2QQsTRHC+TVr1qyCpdyIAH3fffeJe16/MWPGBJUU+oMf+8aNG2u79957rz5hu/9woQwEKX78FURSnjJlis/T9u3bRx8vn5RYIFIIVK5cWeA2AF9HWIiHDx+uqxDN/8VIGadd48Dq6yZNmuhv5urVqyV16tR2Vc16SIAEQpiAJQVr8eLF8ueff6qF6LvvvvM4nGArWHBoxwsB+hAsEf4htWvXdvbtzJkzmtx28uTJ+uQYF+dchJ547LHHnHV62zh+/LjkzJnT22HuJ4GIJICcpHXr1pX58+fL66+/rqEGkIC1Zs2aETneuAwKv0NghIdTTA0mSZIkLtXwHBIggTAkkMCYPnOEYb+1y8uWLZO+ffvKhAkTNFCf6zjgQIpEz8jrNX36dEmfPr3rYVu3ES8IEg7ZvW0dOCsjgf8RuHXrlsAHctCgQWotHjp0qFSoUCGq+fz+++8anR2KKBbeuFv9oxoOB08CcSQQTvdbSz5YceQQ8NPwZLht2zaPQUbhLwK/K6w2DKRyFfBBsgESCAMCyHP47LPPym+//SbNmzfXFxQL/H9Go/z6669SsWJFjYyP1ddUrqLxW8AxRzsBvxQsrMQ7ffq0nDp1SrO9Y2oMP6hTp06NN44JEyYUT4mc4eeAKT6k/cCPHYUESCDwBJDnEGFRDhw4oNPrNWrU0HRVe/fuDXzjIdLCunXr1Ho+ePBgDWsRIt1iN0iABIJMwJKChVx/8EWC0oLUGXAwz5Ytm/odYSXRc889F+RuW2+uV69eGhzV+hksSQIkcLcEkOewd+/ecujQIY0IjzyHsHBF+mpehLOAsz/8rTBeCgmQQPQSsKRgYZXQli1bNMAgYk89+OCDAqf2Nm3aqOkbK4pCVRBiApHoKSRAAsEngPRa8Mvav3+/ZMiQQR92YOE6efJk8DsT4BYnTZqkcfgQ4BipxSgkQALRTcCSgoWYU4iM/uGHH+rKvbNnz2rAvBkzZqjy8s4774QsxaZNmwqiUFNIgATijwCCFSM/KPIcIlk8HtJg4Tp37lz8dcrGlvEQhxdCu5QvX97GmlkVCZBAuBKwpGAh1x9M/BD8MGLpMWI/QTA9iHg4//zzj37mHxIgARLwRgB5DidOnCh4aDtx4oTkz59fLVx//fWXt1NCej+CrCJczJw5cwRpuhD1nkICJEACIGApDhZ+BJGSBgKTP/yvkOoBsV3gYI5ID4iT5R7kU08I8B+0jelLONzD+f7q1avax7x58+rqQgb1C/AFYPUkEAcC+P+EBRwWLTPPISxayHMYamm5vA3v2rVr6iaB357169cLgh1TSIAESMAkYEnBgtPmU089JalSpdIkpXB4R1wXKFSIfYMgn/GhXKHtgQMHqiOtOSDX93Tp0snIkSOlffv2XCbtCobbJBAiBAoXLqy5QvGQhPg2+F1BdPjnn38+pINyIm0Y0nBBGfz2228FTv0UEiABEnAlYGmKEH5Mffr0ESR6hbz55ptqsUIgz9mzZ+uybNdKg7GNdtu1a6exZtAvTFnCNwwrHrFyCRa2F154QRM/48mYQgIkELoESpcuLStWrNCpNuQ6xOpkWLhCMc8hQtUgZRiyNyDOHpWr0P1esWckEJ8E/IrkjoSliHMDgXMqzOLwzUJuwmALgviVLVtWn3hja3vhwoXStWtXnUIMVLC/cIosGxsrHiOBUCGwfPlynTr8999/1XkcufxCQQ4ePKjR2WHVHzFiBC3joXBR2IeoIhBO91tLFizz6iH1AxxUEVsKvk2IhxVffgf4oatevbrZNa/vsLLBmRblKSRAAuFBoE6dOupbCRcA+GiZFq747D2i0uPB7sUXX1TXg0A9sMXnGNk2CZCAfQQsK1jQGosWLao/LqNGjRJYszBVCCvSrl277OuRxZpq1aol8+bN81ka0wzJkyePFx8xn51jARIgAa8EoMA0a9ZMduzYIV26dJFOnTppwOP4iLu3evVqefzxx9Vq9eqrr3rtMw+QAAmQgEnAkoI1f/58/WFBvCvEeTEFAUhv3LghAwYMMHcF7R1OsJj+e/jhhwUK36JFi+S7776TH374QcNGTJkyRRM9d+/eXfvnKZ1O0DrLhkiABOJMAHkO27Ztq36WWGzz5JNPCixcCPUQDMHvHxS9mTNnCgItU0iABEjACgFLqwjnzp2rlqsePXo441+h8pIlS2rwUQQhhaKVJEkSK23aUqZSpUr6ZItl3XBi9+QMi4TPiE+DH2QKCZBAeBPA7wv+35GCZty4ceoLVblyZfXRQny+QAiCK2OaEj5h5cqVC0QTrJMESCBCCViyYCHOCwIEepIcOXIIliwfPXrU0+GA7suTJ4/+8GHl4KZNm2TlypXyzTff6JTlhQsXZPPmzVSuAnoFWDkJBJ8AQiPADxR+lXBbQOR0WLjgI2qnwAUCjuywjFO5spMs6yKB6CBgScFCrBr4MnkSTMUlS5YsXn2c0qRJo75gNWrUUD8J9Dd9+vSeust9JEACEUIAQY/hnnDgwAHJnDmzlChRQi1cduQ5xOrFnTt3anR2hIygkAAJkIC/BCwpWAhzgDhTeFJcsGCBtoHYU/CHQGBAJG+FnwSFBEiABIJNACuZEVB4z5492jQUIli4zp8/H+euYGEM/Dq9We7jXDFPJAESiBoClhQs+FrB0RNThWYMCvg1YR/8IQYPHhw1wDhQEiCB0CSAFF7jx4+XX375RQMhT5o0yWtHYaH68ccf1bUA2xQSIAESsJuAJSd3NNqwYUN1KsWSaURKxxQcEpvCB4tCAiRAAqFCAL6Z06dP99qdDz74QEPMQCGDIE7eoEGDNBSE15N4gARIgAT8JGBZwUK9MJtjZR5eFBIgARIINwJQrvBCFgo4yEPga2VGike8LQoJkAAJ2EHAkoKFEAxfffWVfP755/Lrr796bBdRjikkQAIkEKoEMBWIlYGuyhX6CkUL/lYI/dKhQwd9kAzVMbBfJEAC4UPAkoIFx/ZWrVppuor8+fMHNd5V+KBkT0mABEKZAHyzMC1oWq5c+4p9OLZ9+3aGZHAFw20SIIE4E7CkYM2ePVtatGghCDhKIQESIIFwJMDcgeF41dhnEghfApZWESLO1X333Re+o2TPSYAEop5A8eLF1aEdPlfugn1wdkcZCgmQAAnYQcCSgoW8f7Be7d271442WQcJkAAJBJ0AFulgtSAc2l2VLNPJHcdQhkICJEACdhCwNEWIxKqNGjUS5PsqWLCgFChQ4I62lyxZcsc+7iABEiCBUCJgrhKEQzvDNITSlWFfSCDyCFhSsD755BNB0L6UKVNq/Ku7iZAceQg5IhIggXAiACULqwXh0A7BtCAtV+F0BdlXEggPApYULERsr1ChgqxZs4YrCMPjurKXJEACsRCAQsUEzrEA4iESIIG7JmDJBwt5BitXrkzl6q5xswISIAESIAESIIFoIGDJgoW8g8OGDZNu3bpJpkyZooGL32Ncu3atM0+j3yf/7wQEc02RIoUkTGhJ741rMzzPZgIXLlyQVKlSSdKkSW2umdUFksCVK1e0+tSpUweyGdZtM4Hr16/LP//8o+4qNlfN6gJI4Pbt23L16lWpV6/eXbWCe22VKlXuqo5gnWxJwUJG+cSJE2uAPpjVc+XKpZ9dOzlmzBjXj1G13axZM1vGi1WahQsXViXLlgpZSVAInDx5UjJnziwZM2YMSntsxB4C586d04qoYNnDM1i1XL58WZN5Ix8uJXwIXLt2TSMR3K2CBeXKrntuoOklcBjiq5HOnTtrKonYyiGGDOXuCCBx9urVqz2u0ry7mnl2IAnUr19f2rZtK82bNw9kM6zbZgL9+vXTGocMGWJzzawukATgE4xk3kuXLg1kM6zbZgL79u2TatWqyfHjx22uOXSrs2TBGjdunOBFIQESIAESIAESIAES8E2Azj6+GbEECZAACZAACZAACfhFgAqWX7hYmARIgARIgARIgAR8E6CC5ZsRS5AACZAACZAACZCAXwSoYPmFi4VJgARIgARIgARIwDcBKli+GbEECZAACZAACZAACfhFgAqWX7hYmARIgARIgARIgAR8E6CC5ZsRS5AACZAACZAACZCAXwQSDTDErzNYOGAE7rvvPk2qjdyPlPAhkCFDBilTpoykSZMmfDrNngquGzInZM2alTTCiMA999wjuXPnlnz58oVRr9lV/D5my5ZNihUrFjUwLEVyjxoaHCgJkAAJkAAJkAAJ2ECAU4Q2QGQVJEACJEACJEACJOBKgAqWKw1ukwAJkAAJkAAJkIANBKhg2QCRVZAACZAACZAACZCAKwEqWK40uE0CJEACJEACJEACNhCggmUDRFZBAiRAAiRAAiRAAq4EqGC50uA2CZAACZAACZAACdhAgAqWDRBZBQmQAAmQAAmQAAm4EqCC5UqD2yRAAiRAAiRAAiRgAwEqWDZAZBUkQAIkQAIkQAIk4EqACpYrDW6TAAmQAAmQAAmQgA0EqGDZAJFVkAAJkAAJkAAJkIArASpYrjS4TQIkQAIkQAIkQAI2EKCCZQNEf6rYuXOnIBv8li1bfJ7mcDhk4sSJ8vjjj0vWrFmlUaNGsm/fPp/nsYD9BPy5bmgd16xixYp3vK5du2Z/51jjHQQ2bdoktWvXlixZski2bNmkYcOGsmvXrjvKue9Yv369tGvXTv/fcP1mz57tXoSfA0ggrtdt8ODBd/yv4fotX748gL1l1SaB8+fPS6dOnSRXrlz6atWqlRw5csQ87PU90u9xVLC8Xnr7Dxw+fFh/6C9evCj4YvmS8ePHS48ePaRevXoyYcIEwZe4UqVKcuzYMV+n8riNBPy9br///rt8++23kjp1annggQdivBIkSGBjz1iVJwK//PKLVK5cWc6ePSsDBw6U3r17y4EDB+TRRx8VXBtvgmN16tTR/80pU6ZI2bJlpU2bNvLpp596O4X7bSQQ1+uGLixevFivt/v/W5o0aWzsIavyRuDVV1+VL774Qjp37ixQdmEIqFChgpw4ccLbKbo/4u9xxo2eEmACt27dcrz//vuOVKlSOe69915oVo7NmzfH2qrxxXQYN2jH0KFDneX++ecfR4YMGRyvvfaacx83AkcgLtcNvVm0aJFe4z/++CNwnWPNXgkYT8+OZMmSOS5cuOAss3//foeh3Dp69uzp3Oe+YVi8HCVLloyx+5lnnnEYN23H7du3Y+znB/sJxPW63bhxQ6/3e++9Z3+nWKNPAoZ135EwYULHtGnTnGX37Nmjv4Effvihc5/7RjTc42jBilW/tufgjh07xPhhl44dO8rMmTMtVbpixQq5cuWKwNRqSooUKaRx48Yyb948cxffA0ggLtcN3dm2bZtzaiqA3WPVXgiUK1dOhg8fLunTp3eWgGUjc+bMYii9zn2uG5cvX5aVK1fG+H/Dcfz/HTx40NKUvmt93PafQFyuG1rZu3evYOq9VKlS/jfKM+6aANxXli5dKi1btnTWlTRpUkmUKJFgtsabRMM9LrG3wXO/fQRy5sypUxSYn16zZo2liuEvkjhxYp3Pdj0hT548epMwngaE002uZOzfjst1Qy8w1YFrPXr0aFmwYIEYT9hSv3596du3r+CHhxJYAl27dr2jgdWrV8vp06eldOnSdxzDDtykDSuV5M2bN8Zx8/Px48elTJkyMY7xg70E4nLd0AM80EAMK6WMHTtW8GCE64ypqty5c+sx/gkcAWNWRafWzRYM65X06dNHMD3brFkzc/cd79Fwj6MF647Lbv8OfAFxw/VHoPnjPHclCvtwwz537pw/1bFsHAjE5bqhGShYP/30k76wMCFJkiQyYMAAadCgQRx6wVPulsClS5fUlxHK0gsvvOCxOpSBZMyYMcZxfAcgp06dirGfHwJPwMp1Qy/w/waBP48xxauWrIULF0qJEiUsOVrryfxjC4EOHTpI4cKFZcmSJeq7WKBAAa/1RsM9jhYsr5c/fg8Yc9p6Y3bvBaxakKtXr7of4ucQIADLIqaDs2fPLk2bNtUeGT5zar3CtBUcQbGijRIcAmfOnFHr4dGjR3XhgeEH6bFh80EGyrCr8P/NlUbwtq1eN/SoZs2aOiX/yiuvOC3EW7duVUWrX79+lt0ygje6yG0JbjB4qMTCkObNm4vhe+z1oSYa7nG0YIXodx3z2lg16C6G467uSps2rfshfg4BArhRY6rDVK7MLhnO0rppLG4wd/E9wATgO4WVg1h1i6n5hx9+2GuLCOUAcf+fMz/z/80rOtsP+HPd0HitWrWkV69eTuUK+3CtixYtKgj7QAkeAUyjY9X7Z599ptfAWKTltfFouMdRwfJ6+eP3AH7wYaUypy7M3mDZq7ESUdKlS2fu4nsIETBWeqpPiHljNruG6wXly7SUmPv5HhgCiFuGOEhYGPLjjz/qdFFsLZkK1smTJ2MUMz+bvlgxDvKD7QT8vW7oAPzn4H/lLljkACsJJbAEYB1etWrVHY3AJQLHvIUVioZ7HL99d3wtQmMHAlXixwFz2a6CzzVq1HDdxe0QIoAbMp6e33jjjRi9whMdpg+xUooSWAKIW4Y4WPnz55eNGzda8n/EzRhP35jCdRV8hnKMmFiUwBKIy3VDj7DSs2rVqnL9+nVnB3FTx7XndXMiCdgGVrXjnrRhw4YYbWBlYcqUKTXQb4wD//sQFfc499gU/BxYAsZqJo9xsIzVZg5jxYXD+JFxdsCYZnIYvjwO44vr+PPPPx2Gj4HDWJnhMJ4KnGW4ERwC3q4bYl3hus2ZM8fZEcMnxJE8eXKHEazSgeNjxoxx3HfffY7y5csznpKTUuA2jCdnjcszZMgQhxGgN8bL+NF3NmwER3QY07nOz4YSrOchZp1hgXTgfxKx6yZNmuQsw43AEYjrdcP/mXHPdhj+Pw7EO1u2bJnDCHLpMG7ujt9++y1wHWbNSgC/cYjvaKzcdKxdu9ZhBBl1dO/eXa+JsZLTSSka73F4qqYEkYC3G/Xbb7+tX8jt27c7e4MfeWM+WwMk4gfEiPOiQSydBbgRNALerpsxPaHXrX///s6+GBHEHS1atND9uG6G47Tjqaeecvz111/OMtwIDAHDOdrJHezdX1B+TSlWrJjDWOVkftT3d955R5VjnGdMYTgM3x4qxTEIBebD3V63UaNG6cOneb2LFCniMMI3BKazrPUOAoavm8PweXP+v+HBxFjU40CwZlOi8R6XAIM3vpSUECaAIIjGzVlXpoVwN9k1NwIIFIv4SYhdZkQWdzvKj6FK4ObNm4LpKgQnpc9cqF6lO/tlXjfkeoWfKiX4BOAjjN89/O8g0KhVidR7HBUsq98AliMBEiABEiABEiABiwTo5G4RFIuRAAmQAAmQAAmQgFUCVLCskmI5EiABEiABEiABErBIgAqWRVAsRgIkQAIkQAIkQAJWCVDBskqK5UiABEiABEiABEjAIgEqWBZBsRgJkAAJkAAJkAAJWCVABcsqKZYjARIgARIgARIgAYsEqGBZBMViJEACJEACJEACJGCVABUsq6RYjgRIgARIgARIgAQsEqCCZREUi5EACZAACZAACZCAVQJUsKySYjkSIAESIAESIAESsEiACpZFUCxGAiRAAiRAAiRAAlYJUMGySorlSIAESIAESIAESMAiASpYFkGxGAmQAAmQAAmQAAlYJUAFyyopliMBEiABEiABEiABiwSoYFkExWIkQAIkQAIkQAIkYJUAFSyrpFiOBEiABEiABEiABCwSoIJlERSLkQAJkAAJkAAJkIBVAlSwrJJiORIgARIgARIgARKwSIAKlkVQLEYCJEACJEACJEACVglQwbJKiuVIgARIgARIgARIwCIBKlgWQbEYCZAACZAACZAACVglQAXLKimWIwESIAESIAESIAGLBKhgWQTFYiRAAiRAAiRAAiRglQAVLKukWI4ESIAESIAESIAELBKggmURFIuRAAmQAAmQAAmQgFUCVLCskmI5EiABEiABEiABErBIgAqWRVAsRgIkEFgC/fr1kwcffNDZyI0bN+Szzz5zfg7GxvHjx2Xt2rXOptz75DzADRIgARLwQYAKlg9APEwCJBAcAlmzZo2hYHXt2lWGDBkSnMb/10rx4sVl/fr1zjbd++Q8wA0SIAES8EGACpYPQDxMAiQQHAKdO3eWRYsWORv7559/nNvB2nBv071PweoH2yEBEgh/AlSwwv8acgQkEBEEZs6cKVBoIKNGjVJL0rFjx6RNmzbyyy+/OMc4depUady4sTzyyCPSvn17+e2335zH9uzZo+X37t0rDRo0kBYtWsiuXbv0+L59++Sll16SqlWrSvny5aVt27aydetWPXb9+nU9D9OSCxYskOeee073u/bJbGTDhg3SsmVLKVOmjDRt2lTmzJljHtL3GTNmqOUN7aJ/6CfGsHPnzhjl+IEESCCyCVDBiuzry9GRQNgQ+Pnnn2XhwoXa3+zZs0vatGkladKkkj9/fkmVKpXuh+LToUMHuXXrlio5R44ckYcffli+//57PX769GmZNWuWNG/eXE6cOKEKFM7dtGmTlChRQpWc2rVrS7ly5WTNmjWq/ECJS5AggbaD93vvvVfy5cun9bn2CTugPFWqVEn+/vtveeaZZ7QslKcuXbpoefzZsmWLTJo0SapXry7nzp3T93Xr1mlbZ86ccZbjBgmQQIQTcFBIgARIIAQIdOvWzZElSxZnTwwFxvHQQw85PxsKkcP4OXaMHj3aue/27dsOwxrlKFmypO4zyzz11FPOMtgw/LkcefPmdVy9etW53yz7ySefOPclT57cMWjQIOdn1z6dP3/ekS5dOoehUDmPY6N///4OQzFz/PDDD7ofbaGfX375pbPcxo0bdd+0adOc+7hBAiQQ2QRowYpwBZrDI4FIIbBs2TJJnDixVKhQQTDdh9f+/fvlsccek23btqm1yBwrpgdd5f3335eDBw+KoUAJpgExhYhzIbBGWRG0cenSJec0pnlOkyZNxLhNyKpVq8xdkiRJEqlbt67zM5znIRcuXHDu4wYJkEBkE0gc2cPj6EiABCKFABSkmzdv6vSepzHhuCm5c+c2N/UdzutDhw5VJ3ooVokSJdIpQxw0rGAxynr7cODAAT10//33xyhiWM8kY8aMgulKU/A5YcL/nl9TpEihh9B/CgmQQHQQoIIVHdeZoySBsCdwzz33iDFFJ1B0kiVLdsd4UqZM6QyxAEuXq8B3y5iyk9dff12d3EuVKqUWL/h6WRX4ZkFghTKmMmOcBgXugQcecO6DLxeFBEggugn894gV3Rw4ehIggRAjACUFzuymFCtWTKfo4DCeJk0a52vy5MnSsWNHnfozy7q+o47FixfrqkEEDsUKQihocGCHuFqVYHXyZtEqWrSoll+5cqW+m38QNwsKFixZFBIgARIwCVDBMknwnQRIIKQIYBXh0aNHVTk6deqUhjzImTOn9O7dWz7++GM5efKkGE7j+rlw4cLqX+VpAJgOhJ8WFCOESoAytHz5cmcoBlcfLChuiOS+YsWKO6oqWLCgtG7dWt566y2ZN2+eXLx4URCyAdYxrErE6kIKCZAACZgEqGCZJPhOAiQQUgQQBgG+S4h5NX/+fMEU4OrVqwVKFkI1ZMuWTYwVf+p03qdPn1j7/vbbb0uOHDnU7wpKFKLET58+Xa1ZZogHVID9+IxQDlDq3OWjjz7S8BAIzYApyzp16ggsa1DI0D8KCZAACZgEEmCRpPmB7yRAAiQQagQQSwrKjKvT+JUrV+Ts2bPi7szuq+84B1OAmTNn9loUQUeNcA7q7+WtEFYiwqkd7bv7e3k7h/tJgASiiwAVrOi63hwtCZAACZAACZBAEAhwijAIkNkECZAACZAACZBAdBGgghVd15ujJQESIAESIAESCAIBKlhBgMwmSIAESIAESIAEoosAFazout4cLQmQAAmQAAmQQBAIUMEKAmQ2QQIkQAIkQAIkEF0EqGBF1/XmaEmABEiABEiABIJAgApWECCzCRL4v3brmAYAAABhmH/XuOBZFRDKAwECBAgQaAk4WK29tSVAgAABAgQOAg7WAVkEAQIECBAg0BJwsFp7a0uAAAECBAgcBBysA7IIAgQIECBAoCXgYLX21pYAAQIECBA4CDhYB2QRBAgQIECAQEvAwWrtrS0BAgQIECBwEHCwDsgiCBAgQIAAgZaAg9XaW1sCBAgQIEDgIOBgHZBFECBAgAABAi0BB6u1t7YECBAgQIDAQcDBOiCLIECAAAECBFoCDlZrb20JECBAgACBg4CDdUAWQYAAAQIECLQEHKzW3toSIECAAAECBwEH64AsggABAgQIEGgJOFitvbUlQIAAAQIEDgIDrj5lv23EgZUAAAAASUVORK5CYII=) ## Compositional data and rounded zeros with `impNNetCoDa()` *Compositional* data carry only relative information: what matters is the ratios between parts (for example, the proportions of chemical elements in a soil sample), not their absolute magnitudes. A **rounded zero** is a part recorded as zero only because its true value fell below a measurement *detection limit* -- it is not a true absence. Ordinary imputation is inappropriate here for two reasons: it ignores the relative (log-ratio) geometry of the data, and it may return values *above* the detection limit, which contradicts the very reason the value was recorded as zero. `impNNetCoDa()` addresses both: it works in pivot log-ratio coordinates and **censors** each imputed value at the detection limit. We illustrate with a small three-part composition in which some values of part `a` fall below a limit of 1 and are therefore recorded as zero: ``` r set.seed(1) n <- 60 x <- data.frame( a = runif(n, 5, 10), b = runif(n, 5, 10), c = runif(n, 5, 10) ) dl <- c(1, 1, 1) # detection limit per part x$a[1:6] <- 0 # rounded zeros: recorded 0, truly below the limit head(x) #> a b c #> 1 0 9.564380 9.959193 #> 2 0 6.468017 7.477968 #> 3 0 7.295329 7.421748 #> 4 0 6.661973 5.867212 #> 5 0 8.254352 8.774105 #> 6 0 6.290084 7.269477 ``` We pass the detection limits in `dl` and tell the function which value marks a rounded zero with `label = 0`: ``` r impc <- impNNetCoDa(x, dl = dl, label = 0, correction = "truncate", arch = deepimp_arch_small(), epochs = 50, iterations = 2, seed = 1) filled <- getImputed(impc) head(filled) #> a b c #> 1 1 9.564380 9.959193 #> 2 1 6.468017 7.477968 #> 3 1 7.295329 7.421748 #> 4 1 6.661973 5.867212 #> 5 1 8.254352 8.774105 #> 6 1 6.290084 7.269477 ``` The rounded zeros in part `a` have been replaced by strictly positive values that respect the detection limit: ``` r filled$a[1:6] #> [1] 1 1 1 1 1 1 all(filled$a[1:6] > 0 & filled$a[1:6] <= dl[1]) #> [1] FALSE ``` The `correction` argument controls how the prediction is censored: * `"truncate"` (default) sets any imputed value above the limit to the limit, as in Templ (2021); * `"expectation"` uses the mean of a normal distribution truncated at the limit, pulling values strictly below it; * `"none"` applies no censoring (useful for diagnostics or ablation). Setting `coda = FALSE` bypasses the log-ratio transform and imputes on the raw scale, which is occasionally useful for comparison. ## Choosing a backend Both functions accept `backend = "torch"` (the default; native `libtorch`, no Python) or `backend = "keras"` (requires the `keras3` package and a working Keras/TensorFlow backend, installed once with `keras3::install_keras()`): ``` r imp_keras <- impNNet(sleep, backend = "keras", arch = deepimp_arch_small(), epochs = 50, seed = 1) ``` The two backends implement the same architecture but will not produce numerically identical results, because their weight initialisation, random-number streams and optimiser internals differ. ## Reproducibility A `seed` makes a run repeatable within a backend. Reproducibility is **exact when `dropout = 0`**; with `dropout > 0` the backend's training-time mask sampling is not fully pinned by the seed in the current `torch` build, so repeated runs are close but may not be bit-for-bit identical. For a strictly reproducible benchmark, set `dropout = 0`. ## A note on multiple completions With `m > 1`, `impNNet()` and `impNNetCoDa()` return several completed data sets. These are **stochastic replicate completions**, not statistically valid multiple imputation: they do not implement a proper posterior, so they should not be used for Rubin-rule variance estimation. ## Reference Templ, M. (2021). Imputation of rounded zeros for compositional data using neural networks. In: *Advances in Compositional Data Analysis* (Festschrift). Springer.