test_renyi_entropy.cpp

#include "estimation.h"

#include "tim/core/renyi_entropy.h"
#include "tim/core/signal_tools.h"

#include <pastel/sys/random.h>
#include <pastel/sys/string/string_algorithms.h>

#include <pastel/math/matrix/cholesky_decomposition.h>
#include <pastel/math/matrix/matrix_determinant.h>
#include <pastel/math/matrix/random_matrix.h>

using namespace Tim;

namespace
{

    void testRenyiEntropyCase(
        const std::string& name,
        const Signal& signal,
        dreal q,
        dreal correct)
    {
        //const integer kNearest = 1;
        //const integer timeWindowRadius = signal->samples() / 10;

        Signal signalSet[] = {signal};

        const dreal estimate = renyiEntropyLps(
            range(signalSet), q);

        /*
       const dreal estimate = renyiEntropyKl(
           range(signalSet), q, 
           kNearest, 
           norm);
       */

        /*
       std::vector<dreal> entropySet;
       entropySet.reserve(signal->samples());
       temporalRenyiEntropyKl(
           range(signalSet), signal->samples() / 10, std::back_inserter(entropySet), 
           kNearest);
       const dreal estimate = std::accumulate(entropySet.begin(), entropySet.end(), (dreal)0) /
           entropySet.size();
       */

        /*
       log() << name << ": " << estimate << ", correct: " 
           << correct << logNewLine;
       */
        log() << name << ": " <<
            realToString(estimate, 4) << " ("
            << realToString(100 * relativeError<dreal>(estimate, correct), 2) << "%)" << logNewLine;
    }

    void testRenyiEntropy()
    {
        //Timer timer;
        //timer.setStart();

        log() << "Computing Renyi entropies using Leonenko-Pronzato-Savani estimator..." << logNewLine;
        log() << "Relative errors to correct analytic results shown in brackets." << logNewLine;

        const integer dimension = 5;
        const integer samples = 10000;
        const dreal q = 2;

        testRenyiEntropyCase(
            "Gaussian(0, 1)",
            generateGaussian(dimension, samples),
            q,
            gaussianRenyiEntropy(q, dimension, 1));

        if (dimension > 1)
        {
            for (integer i = 1;i <= 32;i *= 2)
            {
                Matrix<dreal> covariance(dimension, dimension);
                const dreal det = (dreal)i;
                const dreal cond = 1;
                setRandomSymmetricPositiveDefinite(
                    det, cond, covariance);

                const CholeskyDecompositionInplace<dreal> cholesky(
                    covariance);

                testRenyiEntropyCase(
                    "Cor.G.(" + realToString(determinant(covariance), 2) + ", " + 
                    realToString(cond) + ")",
                    generateCorrelatedGaussian(dimension, samples, cholesky),
                    q,
                    gaussianRenyiEntropy(q, dimension, determinant(cholesky)));
            }

            for (integer i = 1;i <= 32;i *= 2)
            {
                Matrix<dreal> covariance(dimension, dimension);
                const dreal det = 1;
                const dreal cond = (dreal)i;
                setRandomSymmetricPositiveDefinite(
                    det, cond, covariance);

                const CholeskyDecompositionInplace<dreal> cholesky(
                    covariance);

                testRenyiEntropyCase(
                    "Cor.G.(" + realToString(determinant(covariance), 2) + ", " + 
                    realToString(cond) + ")",
                    generateCorrelatedGaussian(dimension, samples, cholesky),
                    q,
                    gaussianRenyiEntropy(q, dimension, determinant(cholesky)));
            }
        }

        testRenyiEntropyCase(
            "Uniform(-1, 1)",
            generateUniform(dimension, samples),
            q,
            uniformRenyiEntropy(std::pow((dreal)2, (dreal)dimension)));

        //timer.store();
        //log() << "Finished in " << timer.seconds() << " seconds." << logNewLine;
    }

    void testAdd()
    {
        timTestList().add("renyi_entropy", testRenyiEntropy);
    }

    CallFunction run(testAdd);

}