<commit_before>// ---------------------------------------------------------------------
//
// Copyright (C) 2019 by the SampleFlow authors.
//
// This file is part of the SampleFlow library.
//
// The SampleFlow library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of SampleFlow.
//
// ---------------------------------------------------------------------

#ifndef SAMPLEFLOW_CONSUMERS_AUTOCOVARIANCETRACE_H
#define SAMPLEFLOW_CONSUMERS_AUTOCOVARIANCETRACE_H

#include <sampleflow/consumer.h>
#include <sampleflow/types.h>
#include <sampleflow/element_access.h>
#include <mutex>
#include <deque>

#include <boost/numeric/ublas/matrix.hpp>


namespace SampleFlow
{
  namespace Consumers
  {
    /**
     * This is a Consumer class that implements computing the trace of the
     * running sample auto-covariance matrix function:
     * @f{align*}{
     *   \hat\gamma(l)
     *   &=
     *   \frac{1}{n-l-1} \sum_{t=1}^{n-l}{(x_{t+l}-\bar{x})^T(x_{t}-\bar{x})}
     *   \\
     *   &=
     *   \frac{1}{n-l-1} \text{trace}\left[
     *      \sum_{t=1}^{n-l}{(x_{t+l}-\bar{x}) (x_{t}-\bar{x})^T}
     *   \right]
     *   \\
     *   &=
     *   \text{trace}\; \gamma(l).
     * @f}
     * In other words, it calculates the trace of the (auto-)covariance matrix
     * $\gamma(l)$ of samples $x_{t+l}$ and $x_t$
     * with a lag between zero and $k$. The fraction in front is $n-l-1$ because
     * we are comparing $n-l$ pairs of samples, of which $n-l-1$ are statistically
     * independent. This is consistent with the definition of the covariance matrix
     * used by the CovarianceMatrix class.
     *
     * This class updates $\hat\gamma(l), l=0,1,2,3,\ldots,k$ for each new, incoming
     * sample. The value of the maximum lag $k$ is set in the constructor. The
     * auto-covariances $\hat\gamma(l)$ are only defined if the class has seen
     * $n>k+1$ samples.
     *
     * @note This class only calculates a quantity derived from the actual
     *   auto-covariance of the samples (which is a matrix), namely the
     *   trace. The trace is an indication for the size of a matrix and in the
     *   case of the auto-covariance, is used to see how quickly the
     *   correlation of samples in a Markov chain decreases -- in other words,
     *   how quickly $\hat \gamma(l)$ goes to zero as $l\to\infty$. As such,
     *   the quantity computed by this class is a cheap way to assess the decay
     *   of the correlation of successive samples. At
     *   the same time, the trace is not a *measure* on the auto-covariance
     *   matrices since matrices may be large even if their trace is small.
     *   For example, one could imagine processes in which the diagonal
     *   elements of $\gamma(l)$ are small but the off-diagonal entries
     *   are large. In those cases, $\hat\gamma(l)=\text{trace}\;\hat\gamma(l)$
     *   will be small even though the actual auto-covariance matrices are
     *   large. A more rigorous approach is therefore to use the
     *   AutoCovarianceMatrix class to compute the matrices $\gamma(l)$
     *   and to compute some kind of norm of these matrices; the downside is,
     *   of course, that computing $\gamma(l)$ is substantially more expensive
     *   than computing the quantity $\hat\gamma(l)$ used here.
     *
     * @note In the case of scalar data types, the (auto-)covariance matrices
     *   are also just numbers. In those cases, the results of the current
     *   AutoCovarianceTrace and the AutoCovarianceMatrix classes are the
     *   same, though the latter returns $1\times 1$ matrices whereas the
     *   current class just returns numbers for each lag $l$.
     *
     *
     * <h3> Algorithm </h3>
     *
     * The approach to deriving an update formula for $\gamma(l)$ is the
     * same as for many other classes, and similar to the one used in the
     * MeanValue class.
     * In the following, let us only consider the case when we have already seen
     * $n\ge k$ samples, so that all of the $\hat\gamma(l)$ can actually be computed.
     *
     * Let us expand the formula above and then denote some of its parts as $\hat\alpha$
     * and $\hat\beta$:
     * @f{align*}{
     *   \hat\gamma(l)
     *   &=
     *   \frac{1}{n-l-1}\sum_{t=1}^{n-l}{(x_{t+l}-\bar{x}_n)^T(x_{t}-\bar{x}_n)}
     * \\&=
     *   \underbrace{\frac{1}{n-l-1}\sum_{t=1}^{n-l}{x_{t+l}^T x_{t}}}_{\hat\alpha_n(l)}
     *   -
     *   \bar{x}_n^T
     *   \underbrace{\left[ \frac{1}{n-l-1}\sum_{t=1}^{n-l}(x_{t+l}+x_{t}) \right]}_{\hat\beta_n(l)}
     *   +
     *   \frac{n-l}{n-l-1}\bar{x}_n^T \bar{x}_n
     * \\&=
     *   \hat\alpha_n(l)-\bar{x}_n^T \hat \beta_n(l)+\frac{n-l}{n-l-1} \bar{x}_n^T \bar{x}_n
     * \\&=
     *   \hat\alpha_n(l)-\bar{x}_n^T \hat \beta_n(l)+\left(1+\frac{1}{n-l-1}\right) \bar{x}_n^T \bar{x}_n.
     * @f}
     *
     * For each new sample, we then need to update the scalars $\hat\alpha_{n+1}(l)$,
     * vectors $\hat\beta_{n+1}(l)$, and sample mean $\bar{x}_{n+1}$. The
     * principle of updating $\bar x_{n+1}$ is equivalent to what the MeanValue
     * class does. For $\hat\alpha$, we use that
     * @f{align*}{
     *  \hat\alpha_{0}(l) = \cdots = \hat\alpha_{l+1}(l)
     *  &= 0,
     *  \\
     *  \hat\alpha_{l+2}(l)
     *  &=
     *  \sum_{t=1}^2 x_{t+l} x_t
     *  \\
     *  &=
     *  x_{1+l} x_1 + x_{2+l} x_2,
     *  \\
     *  \hat\alpha_{n+1}(l)
     *  &=
     *  \frac{1}{n-l}\sum_{t=1}^{n+1-l}{x_{t+l}^T x_{t}}
     *  \qquad\qquad\qquad (\text{for}\, n\ge l+2)
     *  \\
     *  &=
     *  \frac{1}{n-l}\left[\sum_{t=1}^{n-l}{x_{t+l}^T x_{t}} + x_{n+1}^T x_{n+1-l} \right]
     *  \\
     *  &=
     *  \frac{1}{n-l}\left[(n-l-1)\hat\alpha_n(l) + x_{n+1}^T x_{n+1-l} \right]
     *  \\
     *  &=
     *  \frac{n-l-1}{n-l} \alpha_n(l) + \frac{1}{n-l} x_{n+1}^T x_{n+1-l}.
     *  \\
     *  &=
     *  \hat\alpha_n(l) - \frac{1}{n-l} \hat\alpha_n(l) + \frac{1}{n-l} x_{n+1}^T x_{n+1-l}.
     * @f}
     * Similarly, for $\hat\beta$ we can use that
     * @f{align*}
     *  \hat\beta_{0}(l) = \cdots = \hat\beta_{l+1}(l)
     *  &= 0,
     *  \\
     *  \hat\beta_{l+2}(l)
     *  &=
     *  \frac{1}{l+2-l-1}\sum_{t=1}^{l+2-l}(x_{t+l}+x_{t})
     *  \\
     *  &=
     *  \sum_{t=1}^{2}(x_{t+l}+x_{t})
     *  \\
     *  &=
     *  (x_{l+1}+x_{1}) + (x_{l+2}+x_{2}),
     *  \\
     *   \hat\beta_{n+1}(l)
     *   &= \frac{1}{n-l}\sum_{t=1}^{n+1-l}(x_{t+l}+x_{t})
     *   \qquad\qquad\qquad (\text{for}\, n\ge l+2)
     * \\
     *   &= \frac{1}{n-l} \left[\sum_{t=1}^{n-l}(x_{t+l}+x_{t}) + (x_{n+1}+x_{n+1-l})\right]
     * \\
     *   &= \frac{1}{n-l} \left[(n-l-1)\beta_n(l) + (x_{n+1}+x_{n+1-l})\right]
     * \\
     *   &= \beta_n(l) - \frac{1}{n-l}\beta_n(l) + \frac{1}{n-l} (x_{n+1}+x_{n+1-l}).
     * @f}
     *
     *<commit_msg>Fix up a couple of formulas.<commit_after>// ---------------------------------------------------------------------
//
// Copyright (C) 2019 by the SampleFlow authors.
//
// This file is part of the SampleFlow library.
//
// The SampleFlow library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of SampleFlow.
//
// ---------------------------------------------------------------------

#ifndef SAMPLEFLOW_CONSUMERS_AUTOCOVARIANCETRACE_H
#define SAMPLEFLOW_CONSUMERS_AUTOCOVARIANCETRACE_H

#include <sampleflow/consumer.h>
#include <sampleflow/types.h>
#include <sampleflow/element_access.h>
#include <mutex>
#include <deque>

#include <boost/numeric/ublas/matrix.hpp>


namespace SampleFlow
{
  namespace Consumers
  {
    /**
     * This is a Consumer class that implements computing the trace of the
     * running sample auto-covariance matrix function:
     * @f{align*}{
     *   \hat\gamma(l)
     *   &=
     *   \frac{1}{n-l-1} \sum_{t=1}^{n-l}{(x_{t+l}-\bar{x})^T(x_{t}-\bar{x})}
     *   \\
     *   &=
     *   \frac{1}{n-l-1} \text{trace}\left[
     *      \sum_{t=1}^{n-l}{(x_{t+l}-\bar{x}) (x_{t}-\bar{x})^T}
     *   \right]
     *   \\
     *   &=
     *   \text{trace}\; \gamma(l).
     * @f}
     * In other words, it calculates the trace of the (auto-)covariance matrix
     * $\gamma(l)$ of samples $x_{t+l}$ and $x_t$
     * with a lag between zero and $k$. The fraction in front is $n-l-1$ because
     * we are comparing $n-l$ pairs of samples, of which $n-l-1$ are statistically
     * independent. This is consistent with the definition of the covariance matrix
     * used by the CovarianceMatrix class.
     *
     * This class updates $\hat\gamma(l), l=0,1,2,3,\ldots,k$ for each new, incoming
     * sample. The value of the maximum lag $k$ is set in the constructor. The
     * auto-covariances $\hat\gamma(l)$ are only defined if the class has seen
     * $n>k+1$ samples.
     *
     * @note This class only calculates a quantity derived from the actual
     *   auto-covariance of the samples (which is a matrix), namely the
     *   trace. The trace is an indication for the size of a matrix and in the
     *   case of the auto-covariance, is used to see how quickly the
     *   correlation of samples in a Markov chain decreases -- in other words,
     *   how quickly $\hat \gamma(l)$ goes to zero as $l\to\infty$. As such,
     *   the quantity computed by this class is a cheap way to assess the decay
     *   of the correlation of successive samples. At
     *   the same time, the trace is not a *measure* on the auto-covariance
     *   matrices since matrices may be large even if their trace is small.
     *   For example, one could imagine processes in which the diagonal
     *   elements of $\gamma(l)$ are small but the off-diagonal entries
     *   are large. In those cases, $\hat\gamma(l)=\text{trace}\;\hat\gamma(l)$
     *   will be small even though the actual auto-covariance matrices are
     *   large. A more rigorous approach is therefore to use the
     *   AutoCovarianceMatrix class to compute the matrices $\gamma(l)$
     *   and to compute some kind of norm of these matrices; the downside is,
     *   of course, that computing $\gamma(l)$ is substantially more expensive
     *   than computing the quantity $\hat\gamma(l)$ used here.
     *
     * @note In the case of scalar data types, the (auto-)covariance matrices
     *   are also just numbers. In those cases, the results of the current
     *   AutoCovarianceTrace and the AutoCovarianceMatrix classes are the
     *   same, though the latter returns $1\times 1$ matrices whereas the
     *   current class just returns numbers for each lag $l$.
     *
     *
     * <h3> Algorithm </h3>
     *
     * The approach to deriving an update formula for $\gamma(l)$ is the
     * same as for many other classes, and similar to the one used in the
     * MeanValue class.
     * In the following, let us only consider the case when we have already seen
     * $n\ge k$ samples, so that all of the $\hat\gamma(l)$ can actually be computed.
     *
     * Let us expand the formula above and then denote some of its parts as $\hat\alpha$
     * and $\hat\beta$:
     * @f{align*}{
     *   \hat\gamma(l)
     *   &=
     *   \frac{1}{n-l-1}\sum_{t=1}^{n-l}{(x_{t+l}-\bar{x}_n)^T(x_{t}-\bar{x}_n)}
     * \\&=
     *   \underbrace{\frac{1}{n-l-1}\sum_{t=1}^{n-l}{x_{t+l}^T x_{t}}}_{\hat\alpha_n(l)}
     *   -
     *   \bar{x}_n^T
     *   \underbrace{\left[ \frac{1}{n-l-1}\sum_{t=1}^{n-l}(x_{t+l}+x_{t}) \right]}_{\hat\beta_n(l)}
     *   +
     *   \frac{n-l}{n-l-1}\bar{x}_n^T \bar{x}_n
     * \\&=
     *   \hat\alpha_n(l)-\bar{x}_n^T \hat \beta_n(l)+\frac{n-l}{n-l-1} \bar{x}_n^T \bar{x}_n
     * \\&=
     *   \hat\alpha_n(l)-\bar{x}_n^T \hat \beta_n(l)+\left(1+\frac{1}{n-l-1}\right) \bar{x}_n^T \bar{x}_n.
     * @f}
     *
     * For each new sample, we then need to update the scalars $\hat\alpha_{n+1}(l)$,
     * vectors $\hat\beta_{n+1}(l)$, and sample mean $\bar{x}_{n+1}$. The
     * principle of updating $\bar x_{n+1}$ is equivalent to what the MeanValue
     * class does. For $\hat\alpha$, we use that
     * @f{align*}{
     *  \hat\alpha_{1}(l) = \cdots = \hat\alpha_{l+1}(l)
     *  &= 0,
     *  \\
     *  \hat\alpha_{l+2}(l)
     *  &=
     *  \sum_{t=1}^2 x_{t+l}^T x_t
     *  \\
     *  &=
     *  x_{1+l}^T x_1 + x_{2+l}^T x_2,
     *  \\
     *  \hat\alpha_{n+1}(l)
     *  &=
     *  \frac{1}{n-l}\sum_{t=1}^{n+1-l}{x_{t+l}^T x_{t}}
     *  \qquad\qquad\qquad (\text{for}\, n\ge l+2)
     *  \\
     *  &=
     *  \frac{1}{n-l}\left[\sum_{t=1}^{n-l}{x_{t+l}^T x_{t}} + x_{n+1}^T x_{n+1-l} \right]
     *  \\
     *  &=
     *  \frac{1}{n-l}\left[(n-l-1)\hat\alpha_n(l) + x_{n+1}^T x_{n+1-l} \right]
     *  \\
     *  &=
     *  \frac{n-l-1}{n-l} \alpha_n(l) + \frac{1}{n-l} x_{n+1}^T x_{n+1-l}.
     *  \\
     *  &=
     *  \hat\alpha_n(l) - \frac{1}{n-l} \hat\alpha_n(l) + \frac{1}{n-l} x_{n+1}^T x_{n+1-l}.
     * @f}
     * Similarly, for $\hat\beta$ we can use that
     * @f{align*}
     *  \hat\beta_{1}(l) = \cdots = \hat\beta_{l+1}(l)
     *  &= 0,
     *  \\
     *  \hat\beta_{l+2}(l)
     *  &=
     *  \frac{1}{l+2-l-1}\sum_{t=1}^{l+2-l}(x_{t+l}+x_{t})
     *  \\
     *  &=
     *  \sum_{t=1}^{2}(x_{t+l}+x_{t})
     *  \\
     *  &=
     *  (x_{l+1}+x_{1}) + (x_{l+2}+x_{2}),
     *  \\
     *   \hat\beta_{n+1}(l)
     *   &= \frac{1}{n-l}\sum_{t=1}^{n+1-l}(x_{t+l}+x_{t})
     *   \qquad\qquad\qquad (\text{for}\, n\ge l+2)
     * \\
     *   &= \frac{1}{n-l} \left[\sum_{t=1}^{n-l}(x_{t+l}+x_{t}) + (x_{n+1}+x_{n+1-l})\right]
     * \\
     *   &= \frac{1}{n-l} \left[(n-l-1)\beta_n(l) + (x_{n+1}+x_{n+1-l})\right]
     * \\
     *   &= \beta_n(l) - \frac{1}{n-l}\beta_n(l) + \frac{1}{n-l} (x_{n+1}+x_{n+1-l}).
     * @f}
     *
     *
{'fraction_non_alphanumeric': 0.1988872455806747, 'fraction_numerical': 0.014648918937953377, 'mean_word_length': 2.1675217488289094, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3133807', 'n_tokens_mistral': 5380, 'n_tokens_neox': 5298, 'n_words': 1357}