Fix inverse predictions for replicate measurements

For details, see NEWS.md

7 files changed, 148 insertions, 87 deletions

diff --git a/vignettes/chemCal.R b/vignettes/chemCal.R
index 701db7b..d9015e9 100644
--- a/vignettes/chemCal.R
+++ b/vignettes/chemCal.R
@@ -1,36 +1,23 @@
-### R code from vignette source '/home/jranke/git/chemCal/vignettes/chemCal.Rnw'
-
-###################################################
-### code chunk number 1: chemCal.Rnw:38-42
-###################################################
+## ------------------------------------------------------------------------
 library(chemCal)
-data(massart97ex3)
 m0 <- lm(y ~ x, data = massart97ex3)
 calplot(m0)
 
+## ------------------------------------------------------------------------
+plot(m0, which=3)
 
-###################################################
-### code chunk number 2: chemCal.Rnw:49-50
-###################################################
-plot(m0,which=3)
-
+## ---- message = FALSE, echo = TRUE---------------------------------------
+weights <- with(massart97ex3, {
+  yx <- split(y, x)
+  ybar <- sapply(yx, mean)
+  s <- round(sapply(yx, sd), digits = 2)
+  w <- round(1 / (s^2), digits = 3)
+})
+massart97ex3.means <- aggregate(y ~ x, massart97ex3, mean)
 
-###################################################
-### code chunk number 3: chemCal.Rnw:56-63
-###################################################
-attach(massart97ex3)
-yx <- split(y, x)
-ybar <- sapply(yx, mean)
-s <- round(sapply(yx, sd), digits = 2)
-w <- round(1 / (s^2), digits = 3)
-weights <- w[factor(x)]
-m <- lm(y ~ x, w = weights)
+m <- lm(y ~ x, w = weights, data = massart97ex3.means)
 
-
-###################################################
-### code chunk number 4: chemCal.Rnw:69-71
-###################################################
+## ------------------------------------------------------------------------
 inverse.predict(m, 15, ws=1.67)
 inverse.predict(m, 90, ws = 0.145)
 
-
diff --git a/vignettes/chemCal.Rmd b/vignettes/chemCal.Rmd
index 2515abb..ccbbdc7 100644
--- a/vignettes/chemCal.Rmd
+++ b/vignettes/chemCal.Rmd
@@ -8,19 +8,14 @@ output:
     toc_float: true
     code_folding: show
     fig_retina: null
-bibliography: refs.bib
+bibliography: references.bib
 vignette: >
   %\VignetteEngine{knitr::rmarkdown}
   %\VignetteIndexEntry{Introduction to chemCal}
+  %\VignetteEncoding{UTF-8}
 ---
 
-[Wissenschaftlicher Berater, Kronacher Str. 12, 79639 Grenzach-Wyhlen, Germany](http://www.jrwb.de)<br />
-
-```{r, include = FALSE}
-require(knitr)
-opts_chunk$set(engine='R', tidy=FALSE)
-```
-# Basic calibration functions for analytical chemistry
+# Basic calibration functions
 
 The `chemCal` package was first designed in the course of a lecture and lab
 course on "Analytics of Organic Trace Contaminants" at the University of Bremen
@@ -28,13 +23,17 @@ from October to December 2004. In the fall 2005, an email exchange with
 Ron Wehrens led to the belief that it would be desirable to implement the
 inverse prediction method given in @massart97 since it also covers the
 case of weighted regression. Studies of the IUPAC orange book and of DIN 32645
-as well as publications by Currie and the Analytical Method Committee of the
-Royal Society of Chemistry and a nice paper by Castillo and Castells provided
-further understanding of the matter.
-
-At the moment, the package consists of four functions, working on univariate
-linear models of class `lm` or `rlm`, plus two datasets for
-validation.
+(equivalent to ISO 11843), publications by @currie97 and the Analytical
+Method Committee of the Royal Society of Chemistry [@amc89] and a nice paper by
+Castells and Castillo [@castells00] provided some further understanding of the matter.
+
+At the moment, the package consists of four functions 
+([calplot](https://pkgdown.jrwb.de/chemCal/reference/calplot.lm.html),
+[lod](https://pkgdown.jrwb.de/chemCal/reference/lod.html),
+[loq](https://pkgdown.jrwb.de/chemCal/reference/loq.html) and
+[inverse.predict](https://pkgdown.jrwb.de/chemCal/reference/inverse.predict.html)),
+working on univariate linear models of class `lm` or `rlm`, plus several
+datasets for validation.
 
 A [bug report](http://bugs.r-project.org/bugzilla3/show_bug.cgi?id=8877)
 and the following e-mail exchange on the r-devel mailing list about
@@ -42,8 +41,21 @@ prediction intervals from weighted regression entailed some further studies
 on this subject. However, I did not encounter any proof or explanation of the
 formula cited below yet, so I can't really confirm that Massart's method is correct.
 
+In fact, in June 2018 I was made aware of the fact that the inverse prediction
+method implemented in chemCal version 0.1.37 and before did not take the
+variance of replicate calibration standards about their means into account, nor
+the number of replicates when calculating the degrees of freedom. Thanks to
+PhD student Anna Burniol Figols for reporting this issue!
+
+As a consequence, I rewrote `inverse.predict` not to automatically work with
+the mean responses for each calibration standard any more. The example
+calculations from @massart97 can still be reproduced when the regression model
+is calculated using the means of the calibration data as shown below.
+
+# Usage
+
 When calibrating an analytical method, the first task is to generate a suitable
-model. If we want to use the `chemCal` functions, we will have to restrict
+model. If we want to use the `chemCal` functions, we have to restrict
 ourselves to univariate, possibly weighted, linear regression so far.
 
 Once such a model has been created, the calibration can be graphically
@@ -64,16 +76,21 @@ plot(m0, which=3)
 ```
 
 Therefore, in Example 8 in @massart97, weighted regression
-is proposed which can be reproduced by
+is proposed which can be reproduced by the following code.
+Note that we are building the model on the mean values for
+each standard in order to be able to reproduce the results
+given in the book with the current version of chemCal.
 
 ```{r, message = FALSE, echo = TRUE}
-attach(massart97ex3)
-yx <- split(y, x)
-ybar <- sapply(yx, mean)
-s <- round(sapply(yx, sd), digits = 2)
-w <- round(1 / (s^2), digits = 3)
-weights <- w[factor(x)]
-m <- lm(y ~ x, w = weights)
+weights <- with(massart97ex3, {
+  yx <- split(y, x)
+  ybar <- sapply(yx, mean)
+  s <- round(sapply(yx, sd), digits = 2)
+  w <- round(1 / (s^2), digits = 3)
+})
+massart97ex3.means <- aggregate(y ~ x, massart97ex3, mean)
+
+m <- lm(y ~ x, w = weights, data = massart97ex3.means)
 ```
 
 If we now want to predict a new x value from measured y values,
@@ -88,7 +105,7 @@ The weight `ws` assigned to the measured y value has to be
 given by the user in the case of weighted regression, or alternatively,
 the approximate variance `var.s` at this location.
 
-# Some theory for `inverse.predict`
+# Background for `inverse.predict`
 
 Equation 8.28 in @massart97 gives a general equation for predicting the
 standard error $s_{\hat{x_s}}$ for an $x$ value predicted from measurements of
@@ -107,20 +124,32 @@ with
 s_e = \sqrt{ \frac{\sum w_i (y_i - \hat{y_i})^2}{n - 2}}
 \end{equation}
 
-where $w_i$ is the weight for calibration standard $i$, $y_i$ is the $y$
-value observed for standard $i$, $\hat{y_i}$ is the estimated value for
-standard $i$, $n$ is the number of calibration samples, $w_s$ is the weight
-attributed to the sample $s$, $m$ is the number of replicate measurements of
-sample $s$, $\bar{y_s}$ is the mean response for the sample,
-$\bar{y_w} = \frac{\sum{w_i y_i}}{\sum{w_i}}$ is the weighted mean of responses
-$y_i$, and $x_i$ is the given $x$ value for standard $i$.
+In chemCal version before 0.2, I interpreted $w_i$ to be the weight for
+calibration standard $i$, $y_i$ to be the mean value observed for standard $i$,
+and $n$ to be the number of calibration standards. With this implementation
+I was able to reproduce the examples given in the book. However, as noted above,
+I was made aware of the fact that this way of calculation does not take the
+variation of the y values about the means into account. Furthermore, I noticed
+that for the case of unweighted linear calibration with replicate standards,
+`inverse.predict` produced different results than `calibrate` from the
+`investr` package when using the Wald method.
+
+Both issues are now addressed in chemCal starting from version 0.2.1. Here,
+$y_i$ is calibration measurement $i$, $\hat{y_i}$ is the estimated value for
+calibration measurement $i$ and $n$ is the total number of calibration
+measurements.
+
+$w_s$ is the weight attributed to the sample $s$, $m$ is the number of
+replicate measurements of sample $s$, $\bar{y_s}$ is the mean response for the
+sample, $\bar{y_w} = \frac{\sum{w_i y_i}}{\sum{w_i}}$ is the weighted mean of
+responses $y_i$, and $x_i$ is the given $x$ value for standard $i$.
 
 The weight $w_s$ for the sample should be estimated or calculated in accordance
 to the weights used in the linear regression.
 
-I adjusted the above equation in order to be able to take a different
+I had also adjusted the above equation in order to be able to take a different
 precisions in standards and samples into account. In analogy to Equation 8.26
-from \cite{massart97} we get
+from \cite{massart97} I am using
 
 \begin{equation}
 s_{\hat{x_s}} = \frac{1}{b_1} \sqrt{\frac{{s_s}^2}{w_s m} +
diff --git a/vignettes/chemCal.html b/vignettes/chemCal.html
index 222ca8f..5779746 100644
--- a/vignettes/chemCal.html
+++ b/vignettes/chemCal.html
@@ -11,7 +11,7 @@
 
 <meta name="author" content="Johannes Ranke" />
 
-<meta name="date" content="2018-07-05" />
+<meta name="date" content="2018-07-17" />
 
 <title>Introduction to chemCal</title>
 
@@ -236,18 +236,22 @@ div.tocify {
 
 <h1 class="title toc-ignore">Introduction to chemCal</h1>
 <h4 class="author"><em>Johannes Ranke</em></h4>
-<h4 class="date"><em>2018-07-05</em></h4>
+<h4 class="date"><em>2018-07-17</em></h4>
 
 </div>
 
 
-<p><a href="http://www.jrwb.de">Wissenschaftlicher Berater, Kronacher Str. 12, 79639 Grenzach-Wyhlen, Germany</a><br /></p>
-<div id="basic-calibration-functions-for-analytical-chemistry" class="section level1">
-<h1>Basic calibration functions for analytical chemistry</h1>
-<p>The <code>chemCal</code> package was first designed in the course of a lecture and lab course on “Analytics of Organic Trace Contaminants” at the University of Bremen from October to December 2004. In the fall 2005, an email exchange with Ron Wehrens led to the belief that it would be desirable to implement the inverse prediction method given in <span class="citation">Massart et al. (1997)</span> since it also covers the case of weighted regression. Studies of the IUPAC orange book and of DIN 32645 as well as publications by Currie and the Analytical Method Committee of the Royal Society of Chemistry and a nice paper by Castillo and Castells provided further understanding of the matter.</p>
-<p>At the moment, the package consists of four functions, working on univariate linear models of class <code>lm</code> or <code>rlm</code>, plus two datasets for validation.</p>
+<div id="basic-calibration-functions" class="section level1">
+<h1>Basic calibration functions</h1>
+<p>The <code>chemCal</code> package was first designed in the course of a lecture and lab course on “Analytics of Organic Trace Contaminants” at the University of Bremen from October to December 2004. In the fall 2005, an email exchange with Ron Wehrens led to the belief that it would be desirable to implement the inverse prediction method given in <span class="citation">Massart et al. (1997)</span> since it also covers the case of weighted regression. Studies of the IUPAC orange book and of DIN 32645 (equivalent to ISO 11843), publications by <span class="citation">Currie (1997)</span> and the Analytical Method Committee of the Royal Society of Chemistry <span class="citation">(Analytical Methods Committee 1989)</span> and a nice paper by Castells and Castillo <span class="citation">(Castells and Castillo 2000)</span> provided some further understanding of the matter.</p>
+<p>At the moment, the package consists of four functions (<a href="https://pkgdown.jrwb.de/chemCal/reference/calplot.lm.html">calplot</a>, <a href="https://pkgdown.jrwb.de/chemCal/reference/lod.html">lod</a>, <a href="https://pkgdown.jrwb.de/chemCal/reference/loq.html">loq</a> and <a href="https://pkgdown.jrwb.de/chemCal/reference/inverse.predict.html">inverse.predict</a>), working on univariate linear models of class <code>lm</code> or <code>rlm</code>, plus several datasets for validation.</p>
 <p>A <a href="http://bugs.r-project.org/bugzilla3/show_bug.cgi?id=8877">bug report</a> and the following e-mail exchange on the r-devel mailing list about prediction intervals from weighted regression entailed some further studies on this subject. However, I did not encounter any proof or explanation of the formula cited below yet, so I can’t really confirm that Massart’s method is correct.</p>
-<p>When calibrating an analytical method, the first task is to generate a suitable model. If we want to use the <code>chemCal</code> functions, we will have to restrict ourselves to univariate, possibly weighted, linear regression so far.</p>
+<p>In fact, in June 2018 I was made aware of the fact that the inverse prediction method implemented in chemCal version 0.1.37 and before did not take the variance of replicate calibration standards about their means into account, nor the number of replicates when calculating the degrees of freedom. Thanks to PhD student Anna Burniol Figols for reporting this issue!</p>
+<p>As a consequence, I rewrote <code>inverse.predict</code> not to automatically work with the mean responses for each calibration standard any more. The example calculations from <span class="citation">Massart et al. (1997)</span> can still be reproduced when the regression model is calculated using the means of the calibration data as shown below.</p>
+</div>
+<div id="usage" class="section level1">
+<h1>Usage</h1>
+<p>When calibrating an analytical method, the first task is to generate a suitable model. If we want to use the <code>chemCal</code> functions, we have to restrict ourselves to univariate, possibly weighted, linear regression so far.</p>
 <p>Once such a model has been created, the calibration can be graphically shown by using the <code>calplot</code> function:</p>
 <pre class="r"><code>library(chemCal)
 m0 &lt;- lm(y ~ x, data = massart97ex3)
@@ -256,14 +260,16 @@ calplot(m0)</code></pre>
 <p>As we can see, the scatter increases with increasing x. This is also illustrated by one of the diagnostic plots for linear models provided by R:</p>
 <pre class="r"><code>plot(m0, which=3)</code></pre>
 <p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAqAAAAHgCAIAAAD17khjAAAACXBIWXMAAA7DAAAOwwHHb6hkAAAgAElEQVR4nOzdd1iTZ9sG8CtMUYSAiMheIgKCyBQQXhQRFNSqVG0RcUKtuK0LcW+qWK2rrbMt7gWOKs5aoILKUIZQVJYLGaISWfn+iF9KkeAieZJw/g6O9wh3xnOifT151n2zuFwuAQAAgHSRYToAAAAAtDwUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDyAGElNTf3iiy/09fXbtWtnZWW1cOHC8vLyz/lADQ2N7du3f+y7vLy83NzcPme7AMA4FDyAuLh165ajo+Pz589XrVp18ODBUaNG7dy5s3///jU1NUxHE5bAwMAhQ4YQ0cuXL1ksVnBwMNOJAKSHHNMBAOCtZcuWmZqaxsXFKSgoEJGfn1///v0dHBxOnz7Na0EpJisrO2zYMHt7e6aDAEgP7MEDiIuMjAwbGxteu/PY29t/8803cnJvfxGvra1dtGhR165d27dv7+Hh8ffff/PGKyoqJk2apK2traioaGxsvHz5ci6X2+Qmdu/ebWdnxzv+v3v37k8IWVtbu2zZsm7durVv397JyenUqVMNn/qoeM7Ozr/99tvJkydZLNabN28SExNra2vfuxVdXd3o6OgVK1ZYWVmpqqoOHz78+fPnn/CDAEg/LgCIh4EDB7Zt23bbtm2vXr1q8gVjx45VVVXdsmVLbGyst7c3m80uLi7mcrlBQUEdO3Zcu3btyZMnZ86cSUTR0dG8t3To0GHbtm28x1FRUXJycgsXLoyNjf3222+JaOvWrU1uqG/fvq6urk0+FRISoqiouHr16piYmLFjxxLR8ePHPy1eSUnJ0KFD+/fv/+jRo7q6Oh0dnY0bN753Kzo6OnZ2dlFRURwOJzs7u1OnTqGhoR//hw0g/VDwAOIiMzPTxMSEiNq3bz9o0KDNmzfn5OTwn83KymKxWIcPH+Z9W1lZ2bZt259++onL5Q4ZMmTv3r38V1pbW8+ZM4f3mF/wlZWVampqS5cu5b9s4sSJurq6TSYRVPAPHz6UlZX94Ycf+CMDBw60sbH55Hhff/314MGDeY/5Bd/MVngv69u3L/+pkJCQXr16NflTALRyOAcPIC7Mzc2zsrIuX758/vz5uLi4sLAwIhoyZMiePXtUVVUTExNlZGQGDx7Me7GysvLDhw8VFRWJ6Pjx40TE5XILCwsvXbqUmZnZr1+/Rh9+9+7dsrIyLy+vkpIS3kifPn1++uknDofTpk2bD0yYkpJSV1f35Zdf8ke+/PLL4OBg3gH2z4n3gVvhfaCrqyv/KWVl5fr6+g/MD9Cq4Bw8gFioqanhVW+/fv3Wr19/+/bt4uLiLVu2XLhwYdasWURUUFCgoaEhLy/Pf4uGhkb79u2J6Pbt2/3799fU1Ozevfvvv/+upqb27uc/ePCAiFxdXTv+v1GjRvE+Ni4ujvX/5s2b10zI4uJiFovVsWNH/kjnzp25XO7jx48/M94HboX37Xs/AQAIV9EDiImCggITE5M//vjD29ubN9K5c+dvv/02OTn58uXLRKSlpVVaWlpfXy8j8/b38vz8fHl5eSUlJRcXl8GDB586dcrJyUlGRsbR0fHdz+f1ZUFBga6ubqOntLW1MzMzeY/V1dWbCckr2pKSEk1NTd7IkydPiEhTU/Mz433gVpp/IwA0hD14ALFgaGjIP1/OH6ytrU1NTTU2NiYiOzu7mpqaM2fO8J6qrq52dHTct29fUlISh8OJjIzs1auXjIxMZWVldnb2u59vZWUlLy9/+vRp/khUVNSXX37J5XLbtWtn/v+aL1EbGxsZGZkjR47wRw4fPmxhYaGkpPSZ8T5wK82/EQAawh48gFiQkZHZtm3bl19+6ezsPGLECG1t7adPn0ZHR9+5c+fq1atEZGtr+8UXXwQHB69evdrIyGjnzp2vX78eMWIEb6d51apVgYGBpaWlK1eurKurS0hIyMrKMjc353++pqbm9OnTZ8+eXVlZaW1tHR8fv2LFitWrV7NYrCbzPH/+/MSJEw1HunTpYmlpOX78+NmzZ1dVVVlaWh47duzUqVNHjx795Hjy8vI5OTnJyck9evTgb8jQ0FDQVgDgIzB3fR8ANHbt2jVfX19dXV1FRUUTE5NRo0alpaXxn62qqpoxY4aRkZGysrK7u3tCQgJvfN++fSYmJm3btnVycjp27NilS5fs7Ox4V6Q3vE2urq5u/fr1vF1hc3PzH3/8sb6+vskYffv2ffffirlz53K53Orq6oiICDMzs3bt2jk4OJw4ceJz4l26dMnQ0FBZWbm8vLzhbXLNbKXhy7hc7qJFi7755pvP/XMHkEYsroAJMQAAAEBy4Rw8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFGpFM9nNmjWLNyMYAACAaCgpKR05cqRTp06i33QrKvj4+PgpU6Z0796d6SAAANBajBo16smTJyh4oTM3N7ezs2M6BQAAtBYMLpKEc/AAAABSCAUPAAAghVDwAAAAUggFDwDQ2rm7u7P+n5+fH9NxhIvL5To6OmZlZfFHjhw50qVLFzab7enpmZmZyWC2loWCBwBoRZ48eVJVVdVoMC8vLzc3t7KysrKy8vDhw4wEE42rV69OmjQpKSmJP/Lo0aNx48bt2bPn+fPnPj4+w4cPl5pV1FHwAADSr76+/scff9TS0rKxsenYsaOnp2dKSgrvKQ6HU15ebmxsrKysrKyszOBV3yKQmJgoKyvbpk2bhiP29vaurq6ysrIzZszIzMwsLy/nP3vlypX58+d/++23O3bsePcXIzGHggcAkH4LFy6Mjo6+cuXK48ePKyoqgoODfXx8MjIyiOj+/ftycnL29vZsNrtv3745OTlMhxWiuXPnbt++XVVVlT/Sr18//kGL69evGxoastlsIqqtrf3666+/+eYbFRUVS0vLS5cudevWLT09nZncn6R13QcPANAKlZaWbt++PTc3t0OHDkQkKys7ZsyYsrKyVatW/frrrxUVFY6Ojj/++KOBgUF4ePioUaOSk5OZjiw6vOMWXC735MmTISEhW7duZbFYRLRly5Znz56lpaXJy8sT0eTJk3/77bevvvoqLS2N9wLxhz14AAApl5KSYmtry2t3Pl9fX96paGdn5/Pnz3fp0kVBQWHp0qW3bt16/vw5Q0mZUVpaGhAQMG/evKNHjw4bNow3eODAgfDwcF6783z99dc1NTV37txhKOZHQ8EDAEg5GRmZ2traRoO1tbUyMjJE9Ndff126dIn/SllZWUVFRVFHZE51dbWPj4+mpmZaWpqbmxt//MmTJ/r6+o1ebGBg8OjRI9EG/HTiW/Dl5eU5OTl1dXVMBwEAkGw9e/a8c+dOYWFhw8Fjx4717t2biDgczogRIzIyMurq6tasWePp6amsrMxQUgacPHmyrq5uw4YN9fX1HA6Hw+HwrqLX09O7d+9ew1dyudycnJx3W19siUvBT5gw4cqVK7zHz549GzRokJqampmZmbq6+tatW6XmpgUAANFTUVFZsGCBj4/P5cuXa2pqSktL16xZs23btoULFxJR3759v/vuu759+2ppad2+fXvPnj1M5xWpmzdv3rp1S6mBiooKIhozZkxERMSLFy/4r9ywYYOWlpa5uTlzYT+OuFxk98svv9jb2//vf/8jokmTJv31118bN240Nzf/66+/pk+f3rFjx4CAAKYzAgBIqtmzZxsZGc2cOTM7O7tNmzY+Pj4JCQkGBga8Z+fMmTNnzhxmE4rS48eP+Y/XrFmzZs2ad18zbty4zMxMCwuLYcOGsdnsK1euvHjx4vjx4yKM+bnEpeD5SktLT5w4cerUKX9/fyLy8fGpq6uLjIxEwQMAfI5hw4YNGzastrZWTk7s/uUXQywWKzIycuzYsXFxcRUVFbNmzRo4cKCsrCzTuT6C2P01379/n4g8PT35I25ubj/88ANziQAApAfa/aNYWlpaWloyneITics5eD5DQ0M5ObmGF4Okp6c3nJQAAAAA3kuMfpVbvHjxsWPHzMzMjI2Np0+ffu7cOSKKi4vbvHmzj48P0+kAAAAkibgU/JEjR3Jzc3Nzc+/evfv69eurV6/yxv38/KysrFatWsVsPAAAAMkiLgXPnzyIh8Ph8B7ExsZ6enpK1nUNAAAAjBO7c/A8/KV+vLy8ampqXr58yWweAAAAySIue/DNCA4OPnjw4IfMdVNWVhYXFyfo2UePHkncYn8AAACfRgIK/sPnTXz48CF/1b93FRUVpaSkNLwBDwAAQFpJQMGHhISEhIR8yCt79Ohx6NAhQc+2b9++ffv2LZcLAABAfInpOXgAAAD4HCh4AAAAKYSCBwAAkELicg5+9uzZzb8gMjJSNEkAAACkgLgUfG1t7c6dO6uqqnR1dRUVFd99AQoeAADgw4lLwUdFRfXr18/Pz+/s2bNWVlZMxwEAAJBsYnQO3tfXt23btkynAAAAkAZiVPAyMjL79+/X0dFhOggAAIDEE5dD9DxDhw5lOgIAAIA0EKM9eAAAAGgpKHgAAAAphIIHAACQQih4AAAAKYSCBwAAkEIoeAAAACmEggcAAJBCKHgAAAAphIIHAACQQih4AAAAKYSCBwAAkEIoeAAAACmEggcAAJBCKHgAgNbuyJEjXbp0YbPZnp6emZmZvMHo6OguXbpoamoGBQVVVVUxmxA+AQoeAKBVe/To0bhx4/bs2fP8+XMfH5/hw4dzudz09PSwsLCjR48+fPiwrq5uyZIlTMeEj4aCBwBoFZKTk4cNG2Zqampvb79w4cLKykreeGJior29vaurq6ys7IwZMzIzM8vLy8+ePTtkyBBra2slJaWlS5dGR0czGx4+AQoeAED6/frrr0OGDBk4cOAff/yxY8eOkpKSnj17lpSUEFG/fv0OHz7Me9n169cNDQ3ZbHZNTQ3/vbKysoWFhbW1tcxEZ1plZWVxcTHTKT4FCh4AQMpxOJyZM2f+8ccf48aNMzExsbOz27Fjx8CBA1evXk1EysrKHTp04HK5J06cGDVq1Pr161ksVr9+/Y4fP56cnFxaWrpgwQIul/vixQumfw5RS05OdnV11dHRcXBw6NixY1RUVF1dHdOhPgIKHgBAyt28edPExMTS0rLhYHBw8IULF3iPS0tLAwIC5s2bd/To0WHDhhGRo6Pj+vXrg4KCevXq1bNnTzk5OTabzUB05qSlpfn5+YWFhVVUVBQVFcXHxx87dmzBggVM5/oIKHgAACn3+vVrZWXlRoPt27d//fo1EVVXV/v4+Ghqaqalpbm5ufGeffny5YABAzIyMrKzsz08PLp27Soj07r6YtWqVeHh4SNHjmSxWETUpUuX48eP79y5s7S0lOloH6p1/YUBALRClpaWt2/ffvXqVcPBa9euWVtbE9HJkyfr6uo2bNhQX1/P4XA4HA6Xy33w4EHXrl2zs7NfvXq1dOnSiRMnMpSdMTdv3uzfv3/DkQ4dOtja2qampjIV6WOh4AEApJy2tvagQYPGjBnz/Plz3si1a9fCw8PnzJlDRDdv3rx165ZSAxUVFVZWVitXrnR3d7ewsLC2tg4LC2P0J2CAnJxcwysNeaqrq2VlZRnJ8wlQ8AAA0m/r1q2GhoZdunRxcXHp1q3buHHjfv755169ehHRmjVruP/FO90+ZcqUJ0+ePHz4cPXq1a3t+DwRubm5HTlypOHIw4cPMzMz7ezsmIr0seSYDgAAAELXpk2byMjIRYsWZWRkqKmpmZqaysnh3//mhIeH9+rVS1ZWduLEiaqqqleuXJk2bdrixYvbtWvHdLQPhb9gAIDWQlVVlbfXDu9lYGBw48aNBQsWWFlZVVZW2tjYrF+/3t/fn+lcHwEFDwAA0ARdXd19+/YxneLTtbrTKgAAAK0BCh4AAEAKoeABAACkkAwRubm5/fPPP0wnAQAAgBYjR0SVlZW8GY5KSkoWLVpUX1/f/Ht69uwZEhIiinQAAADwSf5zFX27du3s7OzeW/BmZmbCjAQAAACf6z8Fr6SkNGHChCZfl5OTc+/ePUdHx44dO4okGAAAAHw6gRfZ5efne3t78+YfPnv2bLdu3fz8/MzMzJKSkkQYDwAAAD6FwIKfOnVqenq6u7s7Ea1cudLZ2Tk3N9fFxSUiIkKE8QAAAOBTCCz4a9euLVmyJCAgoLy8PD4+/ttvvzUxMRk9enRycrIo8wEAAMAnEFjw9fX1KioqRHT+/Hkul9unTx8ikpWVraqqEl06AAAA+CQC56J3dHTcsWOHgYHBunXrXFxcOnXq9OLFi/379+MSegAAAPEnsODXr1/fv39/V1dXJSWls2fPEpGrq2tWVtahQ4dEGA8AAAA+hcCCt7GxycvLy8jI0NfX19TUJKLFixdbWlp269ZNhPEAAKBllJWV7dy5My0tTUVFxcfHZ/DgwUwnAuFqbi76tm3b2tvb89qdiIYPH452BwCQRImJiZaWlvfv3x84cKC1tfXy5cv9/f2rq6uZzgVC1HgPfvv27e99T2hoqHDCAABAy6uvrw8MDNy5c6efnx9vJCQkZPDgwZs2bZozZw6z2UB4Ghd8eHj4e9+DggcAkCApKSlKSkr8diciGRmZuXPnzpo1CwUvxRoXfElJCSM5AABASJ49e6atrd1oUEdH59mzZ4zkAdH4uPXgMzMzo6KihBQFAACEwcjIKCsrq9FCYnfv3jUyMmIqEoiAwKvouVzugQMHbt++3fC/ifj4+JqamunTp4skGwAAtAAzMzNDQ8Ply5dHRESwWCwievLkyYIFCxYuXMh0NBAigQW/fPnyJUuWdO/ePTMzU01NzcDAgPcLYGxsrCjzAQDA54uOjh45cuSJEyfc3d3Ly8vPnTs3derUESNGMJ0LhEjgIfrdu3fPnDkzNTV1586dtra2N27cKCoqsrKyKi8vF2U+AAD4fNra2levXo2KijIxMenbt++NGzew+y71BO7BFxcX85aS69u374wZM+rr69u3b79gwYLVq1cPGTJEhAkBAKAFsFgsDw8PDw8PpoOAiAjcg1dXV8/OziYiXV1dWVlZ3jLwGhoa6enpoksHAAAAn0Rgwfv5+UVGRu7bt4/FYtnZ2UVFRRUWFu7fv//dey0AAACkVk0NVVQwHeJTCCz4yMhIHx+f48ePE9GaNWtOnDihp6f3008/RUREiDAeAAAAE8rK6PffaeRI6tSJ1q1jOs2nEHgOXlVVde/evbzHtra2jx49SkpKMjExMTY2FlU2AAAA0bp/ny5coJgY+vNPcnAgPz/asIEk89C1wIJvhM1m9+vXT6hRAAAAGFBfT7dvU0wMxcZSfj75+FBQEEVHk7Iy08k+i8CCbzhrcUMWFhbrJPNgBQAAwL+qqigujmJjKSaG1NTI35+iosjVlVgsppO1DIEFr6GhwX/M5XKfPn2akJDQqVOnoUOHiiQYAACAEDx7RmfPUmwsnT9PVlbk709XrpCZGdOxWp7Agt+zZ0+jkYqKii+++KK2tla4iQAAAFrc3btvd9azsqhPH/Lzo59+IlVVpmMJ0UcsNqOqqrpo0aItW7YILw0AAECLqa2l69dp2jTS16dBg6i4mJYsoUeP6NAhCgqS7nanD7/Ijuf169d5eXlCigIAANACSkvp4sW3F80ZG5OfHx0/TnZ2TMcSNYEF/+uvvzYaef78+ffff29vby/kSAAAAB8vL+9tqSclUe/e5O9P69aRlhbTsRgjsOAnTJjQaERGRsbGxmbbtm1CjgTQGuXn5/MmhHZ0dNTT02M6DoCEqKujlBSKiaHDh6mkhPr3p6lTydubFBWZTsY8gQXP4XBEmQOg1aqtrf3uu+9+//13Nzc3Ipo8eXJgYOC6detkZWWZjgYgrl6/posXKTaWTp0idXXy96cdO6TpDrcW8XHn4AGgxS1btiwzMzMrK4vNZhNRWVnZyJEjly9fvmTJEqajAYiZ/Hw6d+4/08wtWkS6ukzHElONC77h7e9NcnJyOn36tNDyALQu9fX127ZtS05O5rU7EampqW3fvt3Z2Xnx4sUs7I4AUIM73LKzydOTAgLot99IRYXpWOKuccGvWLGC94DD4YSHh2toaAwfPlxHR+fx48dHjx598+bN/PnzRR4SQGqVlZVxuVwDA4OGg0ZGRnV1dWVlZerq6kwFA2AYh0PXr1NMDB07RgoK5OdHS5bQ//5Hcjjw/KEa/0mFhobyHzg6Op47d05BQYE3snLlSj8/v/379/POFALA51NWVq6qquJwOG3atOEPVlVVVVVVKUv4PNgAn6KkhM6codhYunCBLC3J358uXCBzc6ZjSSSBE93ExMRMnjyZ3+5EJCcnFxoaGhMTI5JgAK2CoqKip6fn9u3bGw5u27atb9++Df/fB9AisrOzDx8+fOHChfLycqaz/FdeHm3aRP36kakpHT5MXl6UnU3Xr9PcuWj3TybwWAeXyy0qKmo0mJ+fL4fDIwAt6ocffvD09Lxz584XX3xBRMeOHbt48eLly5eZzgVSpaysLCQkJDEx0dnZuby8PC0tbcWKFe/eDi1SdXWUkECxsXTyJHE45O1NU6dS//6EX21biMC29vX1jYiIMDU1HTBgAIvF4nK5Z86ciYiIGDFihLAz1dbWlpaWqqmpycvLC3tbAIwzNja+e/fuDz/8wJtkwtXVddOmTTg+Dy0rKCjI0NAwJydHUVGRiO7du+fr66utrT1gwABRR3n1ii5dosOH6fRpMjIiPz/69Vfq2RN3uLU4gQW/cePGnJwcPz8/NTU1bW3t4uLisrIyDw+PyMhIIUVJSUnZsmVLbGzs06dPuVwui8XS1NQcMGDAtGnTbGxshLRRAHGgrKy8YMECplOA1Lp//35ycvKxY8f4e01mZmZr1qzZsmWL6Ar+wQM6f/4/d7itWUPa2iLaeqsksOBVVFSuXr168eLF5OTkx48f6+joODg4eHh4COm+nT///NPb29vAwCA4ONjQ0FBVVbWioqKwsDAmJsbZ2fnixYsuLi7C2C4AgNTLycmxtrZudEzU3t5+3rx5wt1wfT3dvv12+tj8fPLxoaAgio4myTlAdeXKlT/++OPFixfW1tZBQUFKSkpMJ/oIzZ1QZ7FYXl5eXl5eIsixYMECX1/fQ4cONTrHv2zZspCQkPnz51+9elUEMQAApI+qqmpJSUmjwWfPnvFnX2hhVVX0118UE0NHjhCbTf7+FBVFLi4k8xHrlzKutrZ2zJgxt27dCgoK0tPTu3Tp0urVq2NiYrp37850tA/VRMEHBwfb29s3s+779OnTWzxHamrqnj173r2CT0ZGJjg4mIGzRAAA0qJnz55Pnz69fv16w5ucN2/e7O/v35KbefaMzp6l2Fg6f56srMjfny5fJjOzltyECG3ZsuXZs2dpaWm8Ix+TJ0/+7bffvvrqq7S0NEmZgaqJgt+7d+/Lly8TExMFvUcYBW9iYpKQkDB06NB3n4qPjzc1NW3xLQIAtBLy8vI///zzsGHDvvnmG09Pz7Kysp07d5aWlja6P/MT8aeZy8qiPn3Iz49++kkKllo/cODAunXrGp7X+Prrr5cvX37nzh1J2YlvouArKytlZWVFfKYhLCxs/PjxRUVFw4YNMzQ0VFFRqaysLCgoiImJ2bVr1+7du0UZBgBAyvTv3z8pKWnDhg2LFy9WV1cfPHjw+PHjP/22Z/4dbsePU00N+fvTkiXk4UFSdOvTkydP9PX1Gw0aGBg8evRIggu+yftzcnJy7t275+jo2LFjR2HkGDdunLy8/Lp166KjoxuOW1lZ7du3LzAwUBgbBQBoPfT19aOioj7rI0pL6eLFtxfNGRuTnx/9/jvZ2bVQQPGip6d37949Q0ND/giXy7137967rS+2BP76lp+fP2HChK5du27evPns2bP+/v51dXVsNvv8+fMODg7CiDJ69OjAwMCCgoLi4uLS0lJ1dXVtbW09PT1JOdsBACCd8vLelnpSEvXuTf7+tHYtde7MdCzhCg4OXrRokZOTk+r/n27YsGFD586dzSVnZj2BBT916tT09PSJEycS0cqVK52dnffu3Tt16tSIiIizZ88KKQ2LxdLX12/0+xGHw6mtrcW8HwAAotPwDreCAurfn6ZOJW9vUlRkOpmIjB07NjMz08LCYtiwYWw2+8qVKy9fvjx27BjTuT6CwIK/du3a6tWrAwICysvL4+Pjf/vtNxMTk9GjR4eFhYkyHxEFBwcfPHiQy+W+95WJiYlhYWGCXllVVSV20y8DAIiV16/p4kWKjaVTp0hd/e0dbq6urXCaORaLtX79+uDg4EuXLpWXl8+ZM2fgwIEyEnWnn8CCr6+vV1FRIaLz589zudw+ffoQkaysbFVVlejSERGRp6fnB+6+29jYNHNRqLu7u6rkX9gJANDynj6lc+fo8OF/p5kLDyc9PaZjMc/S0tLS0pLpFJ9IYME7Ojru2LHDwMBg3bp1Li4unTp1evHixf79+81EflNjSEhISEjIh7xSSUnJTvDlHjIyMjidLxrFxcX79++/d++erq5uQECAlZUV04fqH58AACAASURBVIkAoCn8O9yys8nTkwIC6LffSEWF6VjQMgQebVi/fn1WVparq2tGRsaqVauIyNXV9ezZs4sWLRJhPJA8R44csbW1LS4udnV15XK53t7eK1asYDoUAPw/DofOnKHQUNLVpYAAKiujtWvp8WM6dIiCgtDu0kTgHryNjU1eXl5GRoa+vr6mpiYRLV682NLSslu3biKMBxLmyZMnkydPvnTpEn+vferUqY6Ojh4eHr1792Y2G0DrVV1NN27Q5ct0+TIlJZG9Pfn705UrhDnEpFpzsxy0bdvW3t6ew+Hk5uaampoOGzYMh7ihebGxsf379294TF5DQ2PKlCmHDh1CwQOIVF0dpaRQXBxdv05//UWdO5ObG33zDR07RkKagh7ETHMFv23btlWrVhUWFhIRl8sdMWKEo6PjrFmzhFHzs2fPbv4FwlumFlrQkydP9N65MMfAwCA+Pp6RPACtS10dZWXRX39RXBxduEDq6uTlRUFBtHcvqaszHQ5ETWDB79q1a/LkyRMnTuzTp8+oUaOIqFevXrNnz27Xrt0333zT4jlqa2t37txZVVWlq6ur2NR9lih4iaCvr3/ixIlGg1lZWQYGBozkAWgV8vIoLo7i4ujiRWKzycuLAgJo2zbq0IHpZMAoLpdrbW2dmprK/a9u3bpNmTKFy+U+e/aM9zIul7to0SJLS0uucMTGxhJRenq6kD5fWVn5p59+EtKHA09ZWZmWltaFCxf4I/fu3evcuXNKSgqDqQCk0D//cHfs4AYEcDt04BobcydN4u7dyy0oYDoWNNZkw4qGwD34Bw8e+Pr6Nhp0dXUV3p60r69v27ZthfThIBpsNvvQoUOBgYG2trY9evR48ODB2bNn161bZ2Njw3Q0AMnH21O/fp0uXyYFBXJ1JS8viowkyZkdHURJYMF37do1MTGx0ULsqampJiYmQooiIyOzf/9+HR0dIX0+iEbv3r0zMjJOnTp179693r17r127tlOnTkyHApBYeXlvr5I7d45qasjNjby8aNkyarAICkCTBBb8lClTQkND5eTk+vXrR0QFBQVnzpxZtmzZ0qVLhZemyfXgQeK0a9eOd90GAHyKR4/o+nWKi6Pz54nDod69ycuLpk4liZ1SDRghsODHjRtXWVm5bNmyxYsXE5G+vr6CgsL06dNnzJghwngAAK3D48f0559vj8CXlJCHB7m60qRJ1LNnK5wHHlpE0wX/5s2bR48ehYSETJgw4e7duw8ePOjYsaOVlRVvxhsAaFl5eXmLFy9OSEggIhcXlyVLlhgbGzMdCoTv6VO6evXtEfj8fHJyIjc3lDq0lKYLPi8vz8LCYt++faNHj3ZycnJychJxLIDW4++//x40aNB33323ZMkSIjp27FivXr1iY2MdHByYjgZCUFJCCQlvb1XPySFHR/Lyoh07yNaWJGqlMhB/TRe8ubl5nz59oqOjAwMDMXsdgFDNmDFjy5YtAQEBvG/nzJmjp6c3Y8aM69evMxtMzP3888/bt29/+vSpvr7+qlWr3N3dmU4kWGUl/f3321vV790jJyfy8qKoKHJ2JrnmZhsD+BxN/7fFYrHGjx+/bNkyBweH/v37d+jQoeEiuNOnTxdVPAAp9/r169TU1C+++KLh4PDhw8ePH19VVaWkpMRUMDHn5OR069YtHx+f3r1737hxw9PTc+LEic2sFs2AV68oIeHtOfU7d8jBgVxdac0acncnBQWmw0GrIPCXx++++46IXr58uXfv3kZPoeABWsqbN28UFRXl/rsbJycnp6CgwOFwUPBNWrt2bUpKSn5+fufOnXkj58+f9/X1DQ0N7dGjB5PJXr+m+Pi359STksjKitzcaMkS6t2bmpqgE0CoBBY8bwp6ABAqNputrKycmpracC6glJQUNputpqbGYDBxtm/fvkGDBvHbnYi8vb2NjIyioqL27Nkj6jRVVXTz5ttz6gkJZG5OXl40dy65uVGbNqIOA9AATv8AMInFYi1cuDAoKOjAgQO8tZjv3r0bFBS0cOFCpqOJr4qKCsN3pnnR1NR89OiRiBLU1lJq6ttz6g1L3dWVcNAFxAYKHoBhISEhMjIynp6eHTp0IKLnz5+vWLFiwoQJTOcSX507d75582ajwby8vJEjRwpxqw1LPTGRunZ9e5/6kSOkqirE7QJ8KhQ8APMmTpwYHBx87949FotlZmYmhyurmxUeHj506NAjR44MHz6cN/Ldd9+VlJQsWLCghbfEW1Kdd049Lu7tkuqTJtHhw1hSHcQf/h35UNXV1Qq49hWERl5e3hITkX6YwYMHT5s2bcSIEVpaWlpaWvfv33/16tWePXtaZiauJpdUDwig7duxpDpIFsyr8B55eXkBAQG8K56sra0PHTrEdCIAoA0bNuTm5o4YMUJPT2/q1KnPnz8PDAz8rE/My6OdO+nLL0lTkwYNops3KSCAcnLon39oxw4KCEC7g8RpvAevoaHR/BucnJxOnz4ttDzi5Z9//nFzc5szZ87PP//crl27+Pj4b7/99uHDh3PmzGE6GkBrZ2RktGHDhs/6CN7qq3FxdPkyqaiQlxf5+dGGDaSr20IZAZjUuOBXrFjBe8DhcMLDwzU0NIYPH66jo/P48eOjR4++efNm/vz5Ig/JmJUrV06dOnXmzJm8b93d3c+ePWttbR0aGtq+fXtmswHAp2i4pHpd3dvVV7GkOkijxgUfGhrKf+Do6Hju3Dn+ieeVK1f6+fnt37/fzc1NpBmZEx8f32hnXVdX18rK6vbt22I9LyYANFRc/PacesMl1ZcuJSMjppOBWKuqqvr+++8vXLhQXl5uY2Mzb948CwsLpkN9BIHn4GNiYiZPntzwsjI5ObnQ0NCYmBiRBBML9fX1srKyjQZlZGTq6uoYyQMAH+rRIzp8mEJCyNiYevakw4fJzo7OnKHiYjp0iCZNQrtD80pLS+3s7O7du7d8+fL9+/fb29v36dPn2LFjTOf6CAKvoudyuUVFRY0G8/PzW9UNPI6OjmfOnDEzM+OPPH36NC0tzdbWlsFUANC0J0/o2jUsqQ4tYvXq1X369NmyZQvvW2trazc3t4EDB/r5+UnKHVUC29rX1zciIsLU1HTAgAEsFovL5Z45cyYiImLEiBGizMesBQsW9OnTR1VVdfTo0XJycunp6RMnTgwLC2PjFlgAMfHkCV2+TFeu0OXLVFpKHh7k6UnTppFEHUoFMfTHH380WoqlZ8+eurq6t27dcnZ2ZirVRxFY8Bs3bszJyfHz81NTU9PW1i4uLi4rK/Pw8IiMjBRlPmZZWFicP39+5syZYWFhSkpKbdq0CQ8PnzRpEtO5AFqr168pI4PS0+nuXUpPpzt3iMOh3r3J05MmT6bu3bGnDi3l5cuXqu/MUchms1++fMlInk8gsOBVVFSuXr168eLF5OTkx48f6+joODg4eHh4tLbl4a2trePi4t68eVNZWfneewgBoCXV1lJ+Pt29SxkZdPcu3bxJeXlkbEyWlmRhQRMmkIUFWVig1EEYLCwsEhISjI2N+SNVVVUpKSkSdJ1dcyfUWSyWl5eXm5tbYWGhqakpl8ttbe3Op6ioqIjVHgGErbj4bZfz/vfWLWKzyc6OLC3Jy4umTSMrK6y7CqIxa9asMWPGmJub29nZEdGLFy+++eYbHx8fbW1tpqN9qOYKftu2batWreKtG8vlckeMGOHo6Dhr1qxWW/MA0JLKyv7t8owMSkmhujqysCA7O7Kzo9GjydaW2rVjOiW0Up6enhs2bBgyZEinTp3U1NRSUlKGDx/+/fffM53rIwgs+F27dk2ePHnixIl9+vQZNWoUEfXq1Wv27Nnt2rX75ptvRJgQAKTCmzeUm0s3b/7b6KWlZGLyttH9/cnamlpkMnmAFjJ8+HB/f//09PTy8nJra+uWWexAhAQWfGRk5JQpUzZv3lxSUsIbmTFjRllZ2Y8//oiCB4D3qKmhgoK3J855jd7w9PmkSWRpSUZGOH0OYk5RUdHe3p7pFJ9IYME/ePDA19e30aCrq2uruooeAD4U//Q5r9Gzs0lTkywsyNKS/Pxo8WIyN6d3po0CAOERWPBdu3ZNTEwcMGBAw8HU1FQTExPhpwLJxuFwtmzZkpKSYmJiEhoa2rlzZ6YTQUtrePr85k1KTaX27d/unfNWTLezIyUlplMCtGoCC37KlCmhoaFycnL9+vUjooKCgjNnzixbtmzp0qUijAeS5/Dhw0FBQfLy8np6en/88cfKlSvnz5+/fPlypnPBZ3jxgnJy/m305GR68+bt3rmFBfn7U48ehJtIAcSMwIIfN25cZWXlsmXLFi9eTET6+voKCgrTp0+fMWOGCOOBhHn69OmoUaOmTJkSFRXFG9m7d++4ceN69+7t7e3NbDb4UDU1dO/ef25X450+592uNmkSbdpEDW4OBgDxJLDgWSzW9OnTJ0yYcPfu3QcPHnTs2NHKykririEEEVu1apWGhga/3YlozJgxv/zyy+LFi1HwYqqujh4+/M/tallZ1KnT24vbAwJo8WLq1o1kBC5MBQDiSWDBl5SUqKqqKisrOzk5OTk58QYrKysvXrw4ZMgQUcUDCZOZmWn8zr6do6PjwYMHGckDTWg0mczt26Sq+vZgO28yGUtLatOG6ZQA8LkEFnzHjh1tbGwOHTrUcC21oqKiL774gsvliiQbSJ6OHTvev3+/0eDDhw9VVFQYyQNUXk537vzb6GlpVFPz7+nz0aOpRw9SVmY6JQC0vOZmsquurrazs9uxY8dXX30lskAg0cLCwlxcXK5du+bu7s4befLkSUxMDK7NFJHqasrJ+c9kMs+fk6np20b38iIHB9LSYjolAIhCcwUfHR199OjRwMDAy5cvb9q0qW3btiKLBRLKyckpMDDQ09PT09PT1dU1MzPzxIkTlpaWc+fOZTqaNOKvxdLkZDJBQViLBaA1a67gZWVlly1b5ubm9vXXXycmJh46dAiz0MN77d27d9SoUeHh4Tt27NDQ0Pj+++/DwsKYDiUtGk0mk5Hx71osfn40dy7WYgEAvuYKnsfb2zslJWXkyJH29va4Rw4+hI+Pj4+PD9MpJF+jyWTS0khZ+T+TyfTsSTiuBgACvL/giUhHR+fSpUvh4eErV64UdiCA1qWmhsrL//168IDu3KE7dyg9nerqqHt3srKi7t3pyy/J0pLYbKbjAoDEEFjw6enppqam/G/l5eXXrl3r6el5/fp1kQQDkEwNC7uigsrK/tPfDb8qKqi8nN68ITb73y99/bfH262sCFP8AsBnEFjwVlZW7w7i0Cu0RlVVVFb2ni8O5+3LHj8mRUVSU2viy8zs7QMlJWrT5u1jLS3MIQMAwtBEwQcHB9vb29fW1gp6z/Tp04UZCUDIPqqwnzwhBYWmC9vYuAUL+9mzZ7t27SKicePGdezYsaV/ZgBodZoo+L179758+TIxMVHQe1DwIF7eW9j8thZhYX8UOzu7W7duycjIENG8efPs7OySk5OFvVEAkG5NFHxlZaWsrKwSlnoEprRsYTdsa7E8JG5nZ5eSkrJ79+7g4GAi2rVr18SJEx0dHW/cuMF0NACQYE0UvDLmrYSW1coK+6PU1dXdunVr06ZNvHYnonHjxpWVlc2ePZvRXAAg8RoXvMb7FnV2cnI6ffq00PKAhCsooNu3336lplJJCVVVkaoqsdmkpvafy8XZbOrS5e1TDb/at2f6ZxCplJQUIpo6dWrDwVmzZs2ePTslJaVHjx4M5QIAide44FesWMF7wOFwwsPDNTQ0hg8frqOj8/jx46NHj75582b+/PkiDwniisul3Fy6fZtu3Xpb6jIyZGtLPXtSYCB9/z1para2wv5YvCWYX7582fDIWUVFBf8pAIBP07jgQ0ND+Q8cHR3PnTunoKDAG1m5cqWfn9/+/fvd3NxEmhHER20tZWf/O/N5QgIpKJCdHdnZ0ZQpZGlJ76wVC83T09OTlZUdOnTo+fPn+YNDhw6Vk5PT1tZmMBgASDqB98HHxMRs2rSJ3+5EJCcnFxoaOnny5B07dogkG4iBmhq6d49u3nz7xVs7nNfokybRL78Q9jI/2+LFiyMiIrp06TJr1iwiioyM/Oeff1atWsV0LgCQbAILnsvlFhUVNRrMz8+Xk/ug2W1BUlVWUmrqv/voKSlkYPC20QMCsHa4MCxatMja2nrMmDGTJ08mIhUVlZMnTw4aNIjpXAAg2QS2ta+vb0REhKmp6YABA1gsFpfLPXPmTERExIgRI0SZD4SuvJzu3Pl3H5233iiv0UePJjs7wg2Twjd48ODy8nKmUwCAVBFY8Bs3bszJyfHz81NTU9PW1i4uLi4rK/Pw8IiMjBRlPmh5xcX/7qDfvElFRWRlRXZ25OVF06ZR9+7U4LwMAABIqKYL/s2bN6WlpefOnYuPj09OTn78+LGOjo6Dg4OHhweWhJc8vEbnfSUlUXU1WVi8bfS5c8nCgvB3CtAKcDic3377LT09XVVV1dvb29XVlelEIFxNF3xeXp6FhcW+fftGjx7t5eUl4kzwWZq50H3SJNq0CRe6A7RCd+7cGTx4cI8ePdzd3SsqKoKDg93c3H7++WdZWVmmo4GwNF3w5ubmffr0iY6ODgwMxC67uGt0ofutW8Rm/9vou3YRVi4BaN24XO6oUaOWLFkyevRo3sjcuXN9fHy2b9/+7bffMpsNhKfpgmexWOPHj1+2bJmDg0P//v07dOgg02A2UCw2w7DmL3S3taV27ZiOCABiJD09vaamht/uRKSoqBgREbFo0SIUvBQTeJHdd999R0QvX77cu3dvo6daW8G/fPny0KFDeXl5gYGB5ubmDCRo/kJ3e3tq04aBVAAgIR49emRoaNho0MjI6N17oUGaCCz4wsJCUeYQW0FBQb/++isRycjIrFy5snPnzqmpqUJfrpt/WRxvH73hhe5z55K5OeG0GQB8MF1d3dzcXC6X2/CUa05Ojr6+PoOpQNg+bhmuzMzMqKgoIUURQ3PmzNm/f//06dPr6+tra2sTExPLysq6dOnS8lsqLqaYGFqyhPz9qVMnsrKitWuprIz8/OjQISoro+vXadMmCgoiS0u0OwB8FEtLSzabvXXrVv5IZWVlRETEmDFjGEwFwtbcTHYHDhy4fft2fX09fzA+Pr6mpqb1HKLfvHmzu7v7hg0beN86OTndvXvXxMTk0qVLffr0+fTP5V/ozttHv3WLFBX/vSxu507q3LllfgAAACIiOnDgwBdffHH8+PH//e9/5eXlBw8eHDZs2NixY5nOBUIksOCXL1++ZMmS7t27Z2ZmqqmpGRgYZGVl1dfXx8bGijIfs6qrq/mr7/AYGxvLy8sfPHjw4wq+0YXuqamkqfn2ZvSpU8nZGRe6A4BQmZqa3r59+8SJE6mpqZ07dz59+rS1tTXTocRdbW3tjh07Lly4UF5ebmNjM2vWLMk6qSGw4Hfv3j1z5szIyMg9e/YcOHDg3LlzlZWV/fr1a20Tar57EUpdXZ2qqup73vbiBaWl/buPnpWFC90BgFlycnLDhw8fPnw400EkQ2VlZd++fbW0tMaPH89ms69cueLg4LBv377+/fszHe1DCSz44uJid3d3Iurbt++MGTPq6+vbt2+/YMGC1atXDxkyRIQJmaShobF+/frZs2fzRzZu3FhfXz9nzpzGLy0rezvzK++rsJC6dHm7jz5pEi50BwCQLOvWrbO0tNy9ezfvWw8Pj759+3711Vd5eXmSsuiawJTq6urZ2dlEpKurKysrm5SU5OTkpKGhkZ6eLsJ4DDtx4oSbm5u6uvrs2bONjIy2bNkSHx/v7e3dsWPH/8z/mpFBz5/jQncAAKkRGxu7c+fOhiNubm4aGhq3b992cHBgKtVHEVjwfn5+kZGRnTp1CgoKsrOzi4qKWr9+/f79+7W1tUWZj1kuLi4ZGRleXl6LFi3qyuU6yslFurr2qqsjdXVq25ZsbcnWlsaPJ1tbMjBgOiwAALSYiooKDQ2NRoMaGhoSdJ5aYMFHRkZWV1cfP348KChozZo1Li4uBw4ckJWV3bNnjwjjMc/c3LywsJDS0mjECOre/W2p29pSp05MRwMAAGExMzO7efOmkZERf6S6ujotLa1r164MpvooAgteVVWVP4edra3to0ePkpKSTExMjFvnUiXW1pSZyXQIAAAQkbCwsGnTpllbW5uZmRERh8OZMWOGi4uLBF1I/6FXCrDZ7H79+gk1CgAAgJgYOHDg48eP3d3du3fvrqamlpCQ4Orqyr/mTiI0Lvh3Tzk04uTkdPr0aaHlAQAAEAvjx48fPnz4jRs3ysvLV65cKZRpTIWpccGvWLGC94DD4YSHh2toaAwfPlxHR+fx48dHjx598+bN/PnzRR4SAACAAaqqqpJ79LpxwfMnbgsNDXV0dDx37pyCggJvZOXKlX5+fvv373dzcxNBsvLy8mfPnhkbG8vifjMAAICPJHCxmZiYmMmTJ/PbnYjk5ORCQ0NjYmKEkWPChAlXrlzhPX727NmgQYPU1NTMzMzU1dW3bt3K5XKFsVEAAABpJbDguVzuu7O05ufnC2kGn19++SUrK4v3eNKkSX/99dfGjRvPnj07derU6dOnHzlyRBgbBQAAkFYC29rX1zciIsLU1HTAgAEsFovL5Z45cyYiImLEiBFCDVRaWnrixIlTp075+/sTkY+PT11dXWRkZEBAgFC3CwAAIE0E7sFv3LjRxsbGz8+vQ4cOVlZWHTp08PPz69GjR2RkpFAD3b9/n4g8PT35I25ubnfv3hXqRgEAAKSMwD14FRWVq1evXrx4MTk5+fHjxzo6Og4ODh4eHiwWS6iBDA0N5eTkCgsLzc3NeSPp6envX70NAAAAGmjuhDqLxfLy8vLy8hJNlMWLFx87dszMzMzY2Hj69Onnzp0jori4uM2bN/v4+IgmA7SI2tra3bt3Jycnm5mZjR07Vl1dnelEAACtjsCCLykpWbRoUWpqan19fcPxnj17bt26tcVzHDlyJDc3Nzc39+7du69fv7569Spv3M/Pz8rKatWqVS2+RRCS06dPjxgxor6+Xltb++DBg/PmzVuxYsXcuXOZzgUA0LoILPhJkybFx8cPGTJEWVm54Xjbtm2FkWPYsGENv+VwOLwHsbGxnp6euBVeUpSWlg4ZMiQoKOiXX37hjWzevHn69OnOzs4eHh7MZgMAaFUEFnxcXBzvEKso0/C1adOG98DLy4vD4VRVVTX6PQPE04oVK9hsNr/diSgsLCw6OnrBggV//fUXg8EAAFobgQWvp6dnIB5rnAcHBx88ePBD5rrhcDjNXG9fX1+PCXOE7e7du+9O1+zi4nLw4EFG8gAAtFoCC37kyJE7d+4MCwsTZZomeXp6fuDue0pKypQpUwQ9++bNm4qKipbLBU1QV1d/8OBBo8HCwkIcgAEAEDGBBW9qavr999+fPHnSzc2NzWY3fGr69OnCD/avkJCQkJCQD3mls7NzcnKyoGfbt2/f6AeBFvftt996eHgkJSU5ODjwRkpLS0+ePDlv3jxmgwEAtDYCC37OnDlE9PjxY/4MsnzCLvja2trS0lI1NTV5eXmhbghanJub29ChQ3v16jVw4EB3d/e7d+/+/vvvRkZGixcvZjoaAEDrInAmu0LBhBQlJSVlwoQJWlpaCgoKnTp1UlRU1NLSGjduXGpqqpC2CMJw+PDh6Ojo3NzclStXXrt2LTw8PDMzk+lQIIXOnz/fq1cvQ0PDvn37vntiCAA+buWYzMzMP/74Qxh78H/++ae3t7eBgUFwcLChoaGqqmpFRUVhYWFMTIyzs/PFixddXFxafKMgJAEBAVg7AITK1dU1Pj5eVVVVTU0tPj7eyMho6tSpmzZtYjoXgBgRWPBcLvfAgQO3b99uONFNfHx8TU2NMAp+wYIFvr6+hw4darRa3bJly0JCQubPn8+f+gYAWrlVq1bFx8f//vvvo0aN4o2EhIT88MMPISEhFhYWzGYDEB8CC3758uVLlizp3r17ZmammpqagYFBVlZWfX19bGysMHKkpqbu2bPn3bVoZWRkgoODBwwYIIyNAoAkioqKMjMz47c7Ee3YsWPfvn3Tpk27cOECg8EAxIrAc/C7d++eOXNmamrqzp07bW1tb9y4UVRUZGVlVV5eLowcJiYmCQkJTT4VHx9vamoqjI0CgCR6+fIlfzEqvg4dOjx8+JCRPADiSeAefHFxsbu7OxH17dt3xowZ9fX17du3X7BgwerVq4cMGdLiOcLCwsaPH19UVDRs2DBDQ0MVFZXKysqCgoKYmJhdu3bt3r27xbcIABJKSUkpLy+v0WB5ebmJiQkjeQDEk8CCV1dXz87OJiJdXV1ZWdmkpCQnJycNDY309HRh5Bg3bpy8vPy6deuio6MbjltZWe3bty8wMFAYGwUASTRmzJioqKi///7bycmJN7J27dpXr16tXr2a2WAAYkVgwfv5+UVGRnbq1CkoKMjOzi4qKmr9+vX79+/X1tYWUpTRo0cHBgYWFBQUFxeXlpaqq6tra2vr6ekJewV6AJAsGzZsOHXqlLOzs66urp6e3r17954/fz506FDcawPQkMCCj4yMrK6uPn78eFBQ0Jo1a1xcXA4cOCArK7tnzx7hpWGxWPr6+vr6+sLbBABIgdzc3M2bN2/YsCEjI0NbW/vo0aNYrhCgkeZuk9u9e7eMjAwR2draPnr0KCkpydDQUFNTU4TxAACaFhYWJg6LZQCILYFX0aupqT19+pT/LZvN7tevX3l5uZ6enkiCAQAAwKdrYg+ef5F8cHAwf112nnv37mloaIgiFwAAAHyGJgqev7Jnu3btlJSUGj7l5OQ0ZswYUeQCAACAz9BEwf/6669ElJubu3Pnzg4dOog8EgAAAHwugRfZJSYmijIHAAAAtKAmLrLj3R3Hn3O+oqJiwoQJJiYm3t7ef//9t2jjAQAAwKdoXPBlZWW9e/ceOnToxSIUhAAAGblJREFUmTNneCPBwcEHDhxwd3cvKSlxd3e/c+eOyEMCAADAx2l8iH7dunX3799PSkqyt7cnotzc3BMnTvz444+TJ0+uqalxd3dfvXr1b7/9xkRUAAAA+FCN9+BPnjw5Y8YMXrsT0dmzZ4lo6NChRCQvLz9y5MgbN26IOCIAAAB8rMYF/+DBA2tra/63Fy5csLW11dLS4n3bqVOnwsJC0aUDAACAT9K44NXU1Pgrvr9+/frKlSteXl78Z4uKijp27Ci6dAAAAPBJGhe8vb39jh07qquriWjfvn2VlZXe3t68p7hc7oEDBxru3wMAAIB4anyR3ZIlS5ydnW1tbW1sbE6ePGlqaurp6VlbW3vp0qWtW7cmJyf/+eefjAQFAACAD9d4D97W1jY+Pr5Lly7Jycl9+vSJjY2VlZUtLy/v37//9evXDx8+7ObmxkhQAAAA+HCN9+Dv37+vpqa2YcMG/kheXl5tbe2VK1d0dXVZLFZeXh6bzVZXVxdtTgAAAPgI/yn4wsJCHx+f2tra5t8zcODAH374QZipAAAA4LP8p+B1dXWzs7OZigIAAAAtRYaIWCzWL7/8curUqYcPHzKdBwAAAFqAHBFt2rTp3Llz27dvT0tLe/XqVY8ePaytrW1sbGxsbKysrBQVFZkOCQAAAB9Hjog8PDw8PDx435eXl9+5c+fmzZt///33rl27UlNTNTU1LSws7Ozs7Ozs7O3tO3fuzGhgAAAAeL/GV9Gz2Ww3Nzf+vXDV1dUZGRmpqampqalRUVEJCQmqqqpjx45dtWqVyKMCAADAh2pc8I0oKCj06NFDV1e3srIyISGBzWaPHDkyMDBQNOEAAADg07yn4J8/fz527Njr16/7+fktWbLEy8tLVlZWNMkAAADgk72n4BUVFW/cuHHq1ClMYAcAACBBGk9V24iysvLSpUuXLVsmmjQAAADQIt6zB09E48eP53A4IogCAAAALeU9e/BEJCcnN23aNBFEAQAAgJby/oIHAAAAiYOCBwAAkEIoeAAAACmEggcAAJBCKHgAAAAphIIHAACQQih4AAAAKYSCBwAAkEIoeAAAACmEggcAAJBCKHgAAAAphIIHAACQQih4AAAAKYSCBwAAkEIoeAAAACmEggcAAJBCckwHkADZ2dmXL18uLy+3sbHp37+/jAx+KwIAAHGHrmoOl8udP3++p6dnampqRUXF8uXLnZ2d8/Pzmc4FAADwHtiDb87evXvj4uIyMjLYbDZvZO3atV999dX169eZDQYAANA87ME3Z/fu3StWrOC3OxF99913RUVF2dnZDKYCAAB4LxR8c/Lz87t27dpwhMVimZmZPXz4kKlIAAAAHwIF3xxNTc2CgoJGg/n5+VpaWozkAQAA+EAo+OaMHDly5cqVtbW1/JEDBw7Iysp2796dwVQAAADvhYvsmjNlypTExMQePXqMGTNGXV390qVLf/75Z0xMDIvFYjoaAABAc7AH3xx5efmDBw9u3Ljx6dOnCQkJLi4umZmZNjY2TOcCAAB4D+zBv1+/fv369evHdAoAAICPgD14AAAAKYSCBwAAkEIoeAAAACmEggcAAJBCKHgAAAAphIIHAACQQih4AAAAKYSCBwAAkEIoeAAAACmEggcAAJBCKHgAAAAphIIHAACQQih4AAAAKSSOBV9bW/v06dOamhqmgwAAAEgqMSr4lJSUCRMmaGlpKSgodOrUSVFRUUtLa9y4campqUxHAwAAkDDish78n3/+6e3tbWBgEBwcbGhoqKqqWlFRUVhYGBMT4+zsfPHiRRcXF6YzAgAASAxxKfgFCxb4+voeOnRITu4/kZYtWxYSEjJ//vyrV68ylQ0AAEDiiMsh+tTU1MDAwEbtTkQyMjLBwcEpKSmMpAIAAJBQ4lLwJiYmCQkJTT4VHx9vamoq4jwAAAASTVwO0YeFhY0fP76oqGjYsGGGhoYqKiqVlZUFBQUxMTG7du3avXs30wEBAAAkibgU/Lhx4+Tl5detWxcdHd1w3MrKat++fYGBgUwFAwAAkETiUvBENHr06MDAwIKCguLi4tLSUnV1dW1tbT09PRaLxXQ0AAAACSNGBU9ELBZLX19fX1+/4SCHw6mtrVVWVmYqFQAAgMQRr4JvUnBw8MGDB7lc7ntfmZaWtmLFCkHPvnnz5uXLly0aDQAAQExJQMF7enp+4O67rq5uQECAoGf//vtvGxublssFAAAgvlgfsmcsHXr16rVx40ZnZ2emgwAAQGthY2Ozf/9+a2tr0W9aXO6DbwiLzQAAAHwmMSp4LDYDAADQUsTlHDwWmwEAAGhB4lLwWGwGAACgBYnLIXosNgMAANCCxKXgsdgMAABACxKXQ/RYbAYAAKAFiUvBY7EZAACAFiQuBU8iWWzm5MmTaWlpvMcPHjxIS0tTU1NrqQ9vDbhcbnFxsY6ODtNBJMyzZ89UVVUVFBSYDiJJqqqqqqqq1NXVmQ4iYYqKivD/0I/17NkzNzc3TU1NYXz48+fPhfGxH6IVzWT3888/JyUl8b9NTk6+f/++lpYWg5EkTl1dXW5ubteuXZkOImEKCgrU1dXbtWvHdBBJUlFR8erVK21tbaaDSBIul5udnW1ubs50EAlTVFTUo0cPIf25tWnTZvny5SoqKsL48Oa1ooJvZPv27ampqdu2bWM6iCQpLy83MjIqKytjOoiE8fX1nTZtmo+PD9NBJMmePXuuXr2K628+yps3b1RVVTmc/2vv3qNizv8/gL8/TZdpKtOFXGq6SlGRLSQkTikKXa1LkRTb5pLcrW9Y0iJrndU2Zy210emgm0IsG2dbFaVYR6XlsIZcK5VUamZ+f3zOb37zq2b0baN85vn46/N+fWben/e8jno2n/n4TEtfL+Qz4+vrGxQU5OPj09cL6WX95Sp6AAAA6EUIeAAAAAZCwAMAADAQAh4AAICBEPAAAAAMhIAHAABgIAQ8AAAAAyluwKuoqHT+8jqQT1lZWUVFpa9X8flB33oAP6E9oKSkhBsm9gBTf0IV90Y379+/b21t1dLS6uuFfGZqamr09PT6ehWfmbq6Oi6Xq6SkuH9P90BbW1tzc3Of3P/rs4af0B548+aNlpYWi8Xq64X0MsUNeAAAAAbDWwoAAAAGQsADAAAwEAIeAACAgRDwAAAADISABwAAYCAEPAAAAAMh4AEAABgIAQ8AAMBACHgAAAAGQsADAAAwEAIeAACAgRDwAAAADKSgAZ+YmGhvb8/lcqdOnfrnn3/29XL6o9bW1piYGCsrKw6HM2rUqH379r1//16yFw38oLa2tkmTJgUGBkoX0TdZiouLPTw8dHR0TE1N9+3bJ/0lWGhal9ra2vbs2WNhYaGhoTF27NjU1FTpvWhaZ5GRkevXr+9QlNMoJvRQrHiOHz9OCFm7dm16erq3tzebzS4rK+vrRfU7mzZtUlNTi4mJyc3N3b59u6qq6qpVq+hdaGB3rFu3jhCyaNEiSQV9k+XGjRsaGhoLFixIS0ujfwUfPHiQ3oWmybJx40Y1NbXvvvsuJycnIiKCEJKZmUnvQtM6EAqF58+f53A469atk67LaRQzeqhwAS8SiWxtbefPn08P29raLC0tQ0ND+3ZV/U17ezubzd6yZYuksnPnTmVl5Xfv3qGB3ZGVlaWtrW1kZCQJePRNDm9v7xkzZohEInoYFRXl4+MjRtNkE4lE+vr6UVFRksrEiRPnzp0rRtM6ycjI0NTUpN/QSge8nEYxpocKd4peIBDcuXPH39+fHiorK3t7e589e7ZvV9XfVFdXGxsbz5o1S1IxNTVtb29/+fIlGvhBDx8+DA4OPnbs2NChQyVF9E2W+vr67OzspUuXUhRFVw4cOJCRkUHQNLlYLBabzZYM2Ww2i8UiaFonLi4uhYWFd+7cMTQ0lK7LaRRjeqhwAV9dXU0IMTY2llRMTExevHjR3t7ed4vqd3g8XmVl5eTJk+nhu3fvjhw5Ym5ubmhoiAbK19raGhAQEBwc7OPjI11H32R5/PixSCRSUlLy8PDQ1tYePnz49u3bW1tbCZomG0VR69at+/HHH48dO3b9+vWdO3cWFRV99dVXBE3rREdHx8bGxsbGRk1NTboup1GM6aHCBXxtbS0hhMvlSipcLlcsFtN16Ky0tNTZ2bmsrCw5OZnFYqGB8kVFRbFYrL1793aoo2+y0L9Mw8PDHR0dU1JSQkJC4uLi1q5dS9A0ucLDw62srJYtW+bo6Lhjx47Fixe7uroSNK3b5DSKMT1UuIDX1dUlhDQ2NkoqDQ0NFEVpa2v33aL6qZqamqVLlzo4OPB4vDt37jg5ORE0UK709PTU1NRTp06pqqp22IW+yUK/r9qyZcuOHTs8PT23bt26bds2Pp/f1NSEpsnS2trq4OCgra199+7dpqam/Pz8vLy8L7/8UiwWo2ndJKdRjOmhwgX8sGHDCCECgUBSEQgE+vr6nX8jK7iqqipbW9vi4uKioqLMzEwTExO6jgbKkZ+fX1dXZ2JiQlEURVHXr19PSUmhKCorKwt9k8XAwIAQ4uDgIKmMHTtWLBYLBAI0TZbffvutoqLiyJEjo0aN4nA4kydP3rNnz+nTpx88eICmdZOcRjGmhwoX8Dwez9raOicnhx6KRKKcnBzpq8mAECISiXx8fEaPHl1cXDx+/HjpXWigHOHh4ZekWFlZTZ8+/dKlS05OTuibLKampjweLz8/X1IpLCxUU1MzMzND02RRUVEhhLx+/VpSefXqFSFET08PTesmOY1iTA+V+3oBnxpFUZs2bVqyZMnIkSMnTpyYlJRUVVWVnJzc1+vqXwoKCsrLy2fMmJGSkiJdX7hwIYfDQQNlsbS0tLS0lAy5XO7QoUPpT0YJIehbl5SVlTdu3BgVFSUUCh0dHYuKimJjY6Ojo+l3S2hal1xcXL744ouAgIBvvvnGyMiotLQ0JiYmJCRER0eHoGndIycLmBMTffYf9PpUYmKinZ2dlpbWlClTrl271tfL6Xf4fH6X/1qePXtGPwAN7I4JEyZI3+hGjL7JlpCQYGdnx+FwbGxsEhIShEKhZBea1qVXr16tXLnS1NSUzWbT95psaWmR7EXTOjM3N+9woxux3EYxoIeUWOqWkAAAAMAMCvcZPAAAgCJAwAMAADAQAh4AAICBEPAAAAAMhIAHAABgIAQ8AAAAAyHgAQAAGAgBDwAAwEAIeAAAAAZCwAMAADAQAh4AAICBEPAAAAAMhIAHAABgIAQ8AAAAAyHgAQAAGAgBDwAAwEAIeAAAAAZCwAMAADAQAh4AAICBEPAAAAAMhIAHAABgIAQ8AAAAAyHgAQAAGAgBDwAAwEAIeAAAAAZCwAP0d5GRkVRXLl++3OvH4vP5AwcO7JWpBg4cyOfze2UqAOgB5b5eAAB8mKamZmJiYoeira3t27dvtbS0lixZkpSURAgJDAx8+/ZtVlZWh20AUEAIeIDPgJqamr+/f+d6c3Ozn5+fg4PDp18SAPRzOEUP8BlTV1cvKipqb28nhDg6OqakpJw5c4aiKCsrK8n2mzdvCCGJiYn29vYaGho2NjbSJwMaGxvDwsIMDQ0NDQ3Dw8NbWlo6H8XHx2fs2LGSoUgkMjAw+Prrrwkh9fX1y5cvHzZsmJqampmZ2a5du8RicecZ2Gz2iRMnJMPg4GAvLy/JUNbabt++7eHhoaOjo6en5+vrKxAI/k2vABQNAh6AIc6dO+fr6+vu7v7s2bP8/HzJ9oABAw4dOrR8+fKZM2eeOnXKxcUlJCQkISGBECIWiz09PU+dOhUVFXXw4MHKyspt27Z1ntnf3//WrVuPHz+mh0VFRdXV1YGBgYSQ1atXZ2VlRUZGnj592sfHJzo6+uTJk//VsmWt7d27d+7u7i0tLfHx8TExMQUFBaGhof+2RwCKBKfoAT4DNTU1FEVJV3bt2tUhjPX09NTV1YVC4ZAhQwghku23b9/u3LnzP//5T3R0NCHE09Pz/fv3e/bsCQ8P//333/Pz83Nycuj303PmzBkxYkRTU1OHo3t5eamqqubk5ERERBBC0tPTTUxMJk6cSAhpaGiIi4tbvHgx/fTLly+XlpbOnz+/m69LztrKy8tfvHiRkZHh5OREv7o//vijJ70DUFQIeIDPQOeL7Kytrbv53Lt379bV1bm6ur5+/ZquTJ8+/ciRIy0tLSUlJTo6Op6ennRdTU0tICCAvl5PGpfLdXNzy87OjoiIEIvFGRkZgYGB9B8cmZmZhBCxWPzkyZO8vLyKigo3N7fuvy45azMyMmKz2atXr960aZObm1tAQEBAQED3ZwYABDzAZ0DWRXbd8ejRI0LIpEmTOtQFAsHz58+HDRsmfW7AwMCgy0n8/PxWrFhRX19///79R48eLVq0iK6XlZVt3ry5tLS0ra1twoQJOjo63VmS5HN6OWuzsLDIzc3dunXrvHnzKIqaPHny1q1bPTw8ujM/ABB8Bg/AeIMGDSKECAQC8f9nYWFhYGDw7Nkz6cvinj9/3uUkc+fOFYvFFy9eTE9Pt7e3t7KyIoS8efPGyclJR0cnOzu7trb24sWLxsbG3VnSy5cvP7g2QoiLi0tBQcHz589PnDghFAq9vLzu37//75oBoEAQ8AAMZ2Njo6Kicu7cOUnlhx9+mDdvnlgsHjduXG1tbW5uLl1vb29PT0/vchJdXd1p06adOXMmPT1d8va9uLi4paUlLi5u4sSJSkpKjY2N9+7d6/LpSkpKNTU19HZ9fX1hYeEH15aWlmZra9vQ0DB48OCFCxcmJycLhUJZ8wNAZzhFD8AcKioqf//9d0lJiZ2dnfR2ZGTk+vXrGxsbR48eXVBQsHv37tjYWIqipk6d6uzsvHDhwt27d/N4PD6f39raKmtyf3//VatWtbe3S66hMzc3V1JS2rNnT2BgYG1tbUxMjFAoLCwsrKyspN/iS4wZM+bAgQMjRoxQV1fftWuXvr4+XdfX15e1Nmtr64qKiqCgoKCgIKFQePz4cS6XO378+I/XPQCmEQNA/7ZmzRo9PT1Zew0MDA4ePEhv5+XlmZiYaGpqvnnzRnpbKBTu379/1KhR6urqVlZW8fHxIpGIfkpDQ0NoaKiBgcGQIUPCwsLy8/NnzZrV5YFevHihpKTk5uYmXUxOTjY3N+dwOBMmTMjIyMjLy7O3t6fXo6enl5CQQD+svLzc2dmZw+HY2Njw+fyjR49u2LCB3iVnbWlpaWPGjFFXV9fV1Z0xY8aNGzd63kQAxUOJu7orBQAAAHzW8Bk8AAAAAyHgAQAAGAgBDwAAwEAIeAAAAAZCwAMAADAQAh4AAICBEPAAjDJw4EA+n9+DJ0ZFRR09erTX19NZRUWFk5OTUCj8BMcCUGQIeAAgZWVlFy5cWLJkySc41siRIy0tLX/66adPcCwARYaABwASGxsbFhamrPyJ7l0dERGxd+/etra2T3M4AMWEgAdgLBMTk+Tk5IiICCMjIzMzs/j4+KdPn3p6eg4cONDY2Dg1NZV+WHV1dVZW1oIFCwghmzdv5vF4khtctrS0cLncffv29eDoly5doijqypUr9DAtLU1FReX27duEEHt7e01NzaysrF54kQAgAwIegMmio6M9PT0fPnwYHBy8cuVKV1fX7du3P3/+3NnZOSQkpLm5mRBy8eJFCwuLIUOGEEL8/PyePHly8+ZN+ukXLlxoaGiQfLvMf8XNzS0oKCg8PLy1tbWhoWHNmjWbNm0aM2YMIYSiKBcXl/Pnz/feCwWAjhDwAEw2bdq0WbNmsVisZcuWEUICAwPHjx+vrKy8dOnSlpYWgUBACLly5cq4cePoxzs4OPB4vMzMTHp48uTJKVOmGBkZdZj21q1bPB7Pzs4uLy+Prhw6dKiurq7Dww4cOPD69eu9e/dGR0draWlt27ZNsmvcuHFXr17t/RcMAP8LAQ/AZJKvbR08eLD0kP7CVpFIRAj5559/hg4dStcpivL19aVPnjc1NWVnZ0u+/V1aXFzcypUrV6xYMWfOnBUrVhQWFsbGxrLZ7A4PGzRo0Pfffx8TExMfH//LL79IP2DYsGGPHz/GtfQAHw8CHoDJWCyWnCGttrZ2wIABkqGvr295eXlVVdW5c+fev3/v7+/f+SnW1tahoaHh4eE3btx4+vSph4dHeHi4urp650fOnz9fU1Nz+PDhTk5O0nUulysSierr63v4wgDgQz7RRbMA0G/p6uo2NjZKhpMmTdLX18/KyioqKvLw8NDT0+v8lC1bttAbo0aNOnv2rJzJ9+/fr6qq+uDBg6SkpJCQEEm9oaFBSUmJy+X20osAgI4Q8ACKztjYuLq6WjJksVje3t4pKSn37t1LSkr6NzNXVlZ+++23qamppaWlGzZsmD179qBBg+hd1dXVhoaGXZ5RAIBegVP0AIrOxcWlpKREuuLn5/fXX38pKyvPnj27x9OKRKKwsDBXV1cfH5+tW7cOGDBg/fr1kr0lJSUuLi49nhwAPggBD6CgbG1ttbS0CCHu7u5VVVUvXryQ7HJxcdHU1PTx8dHQ0Ojx/Hw+v6Sk5PDhwxRFcTicw4cPJycn0/8tXiwWX7161dPT89+/CgCQhZLc0QIAFJa/v//kyZMjIyPp4cOHD83MzM6fPz9z5syPcbjS0lIvL69Hjx6pqqp+jPkBgCDgAYAQcvPmzcWLF9O3mWtubo6KisrNzX348KGKisrHONyyZctGjx69Zs2ajzE5ANBwkR0AEHt7e3d3919//XXq1KkWFhYURWVnZ3+kdK+srCwvL//5558/xuQAIIF38ADwf9ra2q5du2ZlZUXfuRYAPl8IeAAAAAbCVfQAAAAM9D8deZGBo3YifAAAAABJRU5ErkJggg==" /><!-- --></p>
-<p>Therefore, in Example 8 in <span class="citation">Massart et al. (1997)</span>, weighted regression is proposed which can be reproduced by</p>
-<pre class="r"><code>attach(massart97ex3)
-yx &lt;- split(y, x)
-ybar &lt;- sapply(yx, mean)
-s &lt;- round(sapply(yx, sd), digits = 2)
-w &lt;- round(1 / (s^2), digits = 3)
-weights &lt;- w[factor(x)]
-m &lt;- lm(y ~ x, w = weights)</code></pre>
+<p>Therefore, in Example 8 in <span class="citation">Massart et al. (1997)</span>, weighted regression is proposed which can be reproduced by the following code. Note that we are building the model on the mean values for each standard in order to be able to reproduce the results given in the book with the current version of chemCal.</p>
+<pre class="r"><code>weights &lt;- with(massart97ex3, {
+  yx &lt;- split(y, x)
+  ybar &lt;- sapply(yx, mean)
+  s &lt;- round(sapply(yx, sd), digits = 2)
+  w &lt;- round(1 / (s^2), digits = 3)
+})
+massart97ex3.means &lt;- aggregate(y ~ x, massart97ex3, mean)
+
+m &lt;- lm(y ~ x, w = weights, data = massart97ex3.means)</code></pre>
 <p>If we now want to predict a new x value from measured y values, we use the <code>inverse.predict</code> function:</p>
 <pre class="r"><code>inverse.predict(m, 15, ws=1.67)</code></pre>
 <pre><code>## $Prediction
@@ -291,8 +297,8 @@ m &lt;- lm(y ~ x, w = weights)</code></pre>
 ## [1] 36.20523 51.91526</code></pre>
 <p>The weight <code>ws</code> assigned to the measured y value has to be given by the user in the case of weighted regression, or alternatively, the approximate variance <code>var.s</code> at this location.</p>
 </div>
-<div id="some-theory-for-inverse.predict" class="section level1">
-<h1>Some theory for <code>inverse.predict</code></h1>
+<div id="background-for-inverse.predict" class="section level1">
+<h1>Background for <code>inverse.predict</code></h1>
 <p>Equation 8.28 in <span class="citation">Massart et al. (1997)</span> gives a general equation for predicting the standard error <span class="math inline">\(s_{\hat{x_s}}\)</span> for an <span class="math inline">\(x\)</span> value predicted from measurements of <span class="math inline">\(y\)</span> according to the linear calibration function <span class="math inline">\(y = b_0 + b_1 \cdot x\)</span>:</p>
 <span class="math display">\[\begin{equation}
 s_{\hat{x_s}} = \frac{s_e}{b_1} \sqrt{\frac{1}{w_s m} + \frac{1}{\sum{w_i}} +
@@ -304,9 +310,11 @@ s_{\hat{x_s}} = \frac{s_e}{b_1} \sqrt{\frac{1}{w_s m} + \frac{1}{\sum{w_i}} +
 <span class="math display">\[\begin{equation}
 s_e = \sqrt{ \frac{\sum w_i (y_i - \hat{y_i})^2}{n - 2}}
 \end{equation}\]</span>
-<p>where <span class="math inline">\(w_i\)</span> is the weight for calibration standard <span class="math inline">\(i\)</span>, <span class="math inline">\(y_i\)</span> is the mean <span class="math inline">\(y\)</span> value (!) observed for standard <span class="math inline">\(i\)</span>, <span class="math inline">\(\hat{y_i}\)</span> is the estimated value for standard <span class="math inline">\(i\)</span>, <span class="math inline">\(n\)</span> is the number calibration standards, <span class="math inline">\(w_s\)</span> is the weight attributed to the sample <span class="math inline">\(s\)</span>, <span class="math inline">\(m\)</span> is the number of replicate measurements of sample <span class="math inline">\(s\)</span>, <span class="math inline">\(\bar{y_s}\)</span> is the mean response for the sample, <span class="math inline">\(\bar{y_w} = \frac{\sum{w_i y_i}}{\sum{w_i}}\)</span> is the weighted mean of responses <span class="math inline">\(y_i\)</span>, and <span class="math inline">\(x_i\)</span> is the given <span class="math inline">\(x\)</span> value for standard <span class="math inline">\(i\)</span>.</p>
+<p>In chemCal version before 0.2, I interpreted <span class="math inline">\(w_i\)</span> to be the weight for calibration standard <span class="math inline">\(i\)</span>, <span class="math inline">\(y_i\)</span> to be the mean value observed for standard <span class="math inline">\(i\)</span>, and <span class="math inline">\(n\)</span> to be the number of calibration standards. With this implementation I was able to reproduce the examples given in the book. However, as noted above, I was made aware of the fact that this way of calculation does not take the variation of the y values about the means into account. Furthermore, I noticed that for the case of unweighted linear calibration with replicate standards, <code>inverse.predict</code> produced different results than <code>calibrate</code> from the <code>investr</code> package when using the Wald method.</p>
+<p>Both issues are now addressed in chemCal starting from version 0.2.1. Here, <span class="math inline">\(y_i\)</span> is calibration measurement <span class="math inline">\(i\)</span>, <span class="math inline">\(\hat{y_i}\)</span> is the estimated value for calibration measurement <span class="math inline">\(i\)</span> and <span class="math inline">\(n\)</span> is the total number of calibration measurements.</p>
+<p><span class="math inline">\(w_s\)</span> is the weight attributed to the sample <span class="math inline">\(s\)</span>, <span class="math inline">\(m\)</span> is the number of replicate measurements of sample <span class="math inline">\(s\)</span>, <span class="math inline">\(\bar{y_s}\)</span> is the mean response for the sample, <span class="math inline">\(\bar{y_w} = \frac{\sum{w_i y_i}}{\sum{w_i}}\)</span> is the weighted mean of responses <span class="math inline">\(y_i\)</span>, and <span class="math inline">\(x_i\)</span> is the given <span class="math inline">\(x\)</span> value for standard <span class="math inline">\(i\)</span>.</p>
 <p>The weight <span class="math inline">\(w_s\)</span> for the sample should be estimated or calculated in accordance to the weights used in the linear regression.</p>
-<p>I adjusted the above equation in order to be able to take a different precisions in standards and samples into account. In analogy to Equation 8.26 from  we get</p>
+<p>I had also adjusted the above equation in order to be able to take a different precisions in standards and samples into account. In analogy to Equation 8.26 from  I am using</p>
 <span class="math display">\[\begin{equation}
 s_{\hat{x_s}} = \frac{1}{b_1} \sqrt{\frac{{s_s}^2}{w_s m} +
     {s_e}^2 \left( \frac{1}{\sum{w_i}} +
@@ -315,6 +323,15 @@ s_{\hat{x_s}} = \frac{1}{b_1} \sqrt{\frac{{s_s}^2}{w_s m} +
 \end{equation}\]</span>
 <p>where I interpret <span class="math inline">\(\frac{{s_s}^2}{w_s}\)</span> as an estimator of the variance at location <span class="math inline">\(\hat{x_s}\)</span>, which can be replaced by a user-specified value using the argument <code>var.s</code> of the function <code>inverse.predict</code>.</p>
 <div id="refs" class="references">
+<div id="ref-amc89">
+<p>Analytical Methods Committee. 1989. “Robust Statistics — How Not to Reject Outliers. Part 1. Basic Concepts.” <em>The Analyst</em> 114: 1693–7.</p>
+</div>
+<div id="ref-castells00">
+<p>Castells, Reynaldo César, and Marcela Alejandra Castillo. 2000. “Systematic Errors: Detection and Correction by Means of Standard Calibration, Youden Calibration and Standard Additions Method in Conjunction with a Method Response Model.” <em>Analytica Chimica Acta</em> 423: 179–85.</p>
+</div>
+<div id="ref-currie97">
+<p>Currie, L. A. 1997. “Nomenclature in Evaluation of Analytical Methods Including Detection and Quantification Capabilities (IUPAC Recommendations 1995).” <em>Analytica Chimica Acta</em> 391: 105–26.</p>
+</div>
 <div id="ref-massart97">
 <p>Massart, D. L, B. G. M. Vandeginste, L. M. C. Buydens, S. De Jong, P. J. Lewi, and J Smeyers-Verbeke. 1997. <em>Handbook of Chemometrics and Qualimetrics: Part A</em>. Amsterdam: Elsevier.</p>
 </div>
diff --git a/vignettes/figure/unnamed-chunk-1-1.pdf b/vignettes/figure/unnamed-chunk-1-1.pdf
deleted file mode 100644
index c70b645..0000000
--- a/vignettes/figure/unnamed-chunk-1-1.pdf
+++ /dev/null
diff --git a/vignettes/figure/unnamed-chunk-2-1.pdf b/vignettes/figure/unnamed-chunk-2-1.pdf
deleted file mode 100644
index c1b934d..0000000
--- a/vignettes/figure/unnamed-chunk-2-1.pdf
+++ /dev/null
diff --git a/vignettes/references.bib b/vignettes/references.bib
new file mode 100644
index 0000000..a710662
--- /dev/null
+++ b/vignettes/references.bib
@@ -0,0 +1,35 @@
+% Encoding: UTF-8
+@book{massart97,
+  author = "Massart, D. L and Vandeginste, B. G. M. and Buydens, L. M. C. and De Jong, S. and Lewi, P. J. and Smeyers-Verbeke, J",
+  title = "Handbook of {C}hemometrics and {Q}ualimetrics: Part {A}",
+  publisher = "Elsevier", 
+  address = "Amsterdam",
+  year = "1997"
+}
+@article{currie97,
+  author = "Currie, L. A.",
+  year = "1997",
+  title = "Nomenclature in evaluation of analytical methods including
+  detection and quantification capabilities ({IUPAC Recommendations 1995})",
+  journal = "Analytica Chimica Acta",
+  volume = "391",
+  pages = "105 - 126"
+}
+@article{amc89,
+  Title                    = {Robust statistics --- how not to reject outliers. {Part} 1. {Basic} concepts},
+  Author                   = {{Analytical Methods Committee}},
+  Journal                  = {The Analyst},
+  Year                     = {1989},
+  Pages                    = {1693--1697},
+  Volume                   = {114}
+}
+@article{castells00,
+  author = {Reynaldo César Castells and Marcela Alejandra Castillo},
+  title = {Systematic errors: detection and correction by means of standard calibration, Youden calibration and standard additions method in conjunction with a method response model},
+  journal = "Analytica Chimica Acta",
+  volume = "423",
+  year = "2000",
+  pages = "179-185"
+}
+
+@Comment{jabref-meta: databaseType:bibtex;}
diff --git a/vignettes/refs.bib b/vignettes/refs.bib
deleted file mode 100644
index 514d76b..0000000
--- a/vignettes/refs.bib
+++ /dev/null
@@ -1,7 +0,0 @@
-@book{massart97,
-  author = "Massart, D. L and Vandeginste, B. G. M. and Buydens, L. M. C. and De Jong, S. and Lewi, P. J. and Smeyers-Verbeke, J",
-  title = "Handbook of {C}hemometrics and {Q}ualimetrics: Part {A}",
-  publisher = "Elsevier", 
-  address = "Amsterdam",
-  year = "1997"
-}