CRNMultivariateRandomTools.cpp

Go to the documentation of this file.

/* Copyright 2008-2016 INSA-Lyon, ENS-Lyon
 * 
 * This file is part of libcrn.
 * 
 * libcrn is free software: you can 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 3 of the License, or
 * (at your option) any later version.
 * 
 * libcrn is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with libcrn.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * file: MultivariateRandomTools.cpp
 * \author Jean DUONG, Yann LEYDIER
 */

#include <CRNMath/CRNMultivariateRandomTools.h>
#include <CRNMath/CRNMatrixInt.h>
#include <CRNMath/CRNMultivariateGaussianPDF.h>
#include <CRNMath/CRNMultivariateGaussianMixture.h>
#include <CRNMath/CRNMatrixDouble.h>
#include <CRNMath/CRNUnivariateRandomTools.h>
#include <CRNMath/CRNSquareMatrixDouble.h>

using namespace crn;

/*****************************************************************************/
MatrixDouble MultivariateRandomTools::NewGaussianSample(const MatrixDouble& Mu, const SquareMatrixDouble& Sigma, size_t n, size_t m, bool reseed)
{
    size_t dimension = Mu.GetRows();
    MatrixDouble M(n, dimension, 0.0);
    SquareMatrixDouble L = Sigma.MakeCholesky();    

    for (size_t v = 0; v < n; v++) 
    {
        MatrixDouble X = UnivariateRandomTools::NewGaussianSample(0.0, 1.0, dimension, m, reseed);
        MatrixDouble Y = makeMLTransform(Mu, L, X);
        
        for (size_t c = 0; c < dimension; c++)
        {
            M.At(v, c) = Y.At(c, 0);
        }
    }
    return M;
}

/*****************************************************************************/

MatrixDouble MultivariateRandomTools::NewGaussianSample(const MultivariateGaussianPDF& P, size_t n, size_t m, bool reseed)
{
    size_t dimension = P.GetDimension();
    
    MatrixDouble Mu = P.GetMean();
    SquareMatrixDouble Sigma = P.GetVariance();
    
    MatrixDouble Patterns(n, dimension, 0.0);
    SquareMatrixDouble L = Sigma.MakeCholesky();    
        
    for (size_t r = 0; r < n; r++) 
    {
        MatrixDouble X = UnivariateRandomTools::NewGaussianSample(0.0, 1.0, dimension, m, reseed);
        MatrixDouble Y = makeMLTransform(Mu, L, X);
        
        for (size_t c = 0; c < dimension; c++)
            Patterns.At(r, c) = Y.At(c, 0);
    }
    return Patterns;
}

MatrixDouble MultivariateRandomTools::makeMLTransform(const MatrixDouble& M, const SquareMatrixDouble& L, const MatrixDouble& X)
{
    size_t n = X.GetRows();
    
    MatrixDouble Y(n, 1, 0.0);

    for (size_t r = 0; r < n; r++) 
    {
        double Yr = M.At(r, 0);

        for (size_t c = 0; c <= r; c++) 
            Yr += L.At(r, c) * X.At(c, 0);

        Y.At(r, 0) = Yr;
    }

    return Y;
}

/*****************************************************************************/

MatrixDouble MultivariateRandomTools::NewGaussianMixtureSample(const MultivariateGaussianMixture& Mx, size_t n, size_t m, bool reseed)
{
    size_t dim = Mx.GetDimension();
    size_t nbPDF = Mx.GetNbMembers();
    
    MatrixDouble Patterns(n, dim, 0.0);
    
    MatrixDouble CumulWeights(nbPDF, 1, 0.0);
    MatrixDouble IndexeSample = UnivariateRandomTools::NewUniformSample(n, reseed);
    MatrixInt Pop(nbPDF, 1, 0);
    
    // Cumulative weights from the mixture
    
    double Mass = 0.0;
        
    for (size_t k = 0; k < nbPDF; k++)
    {
        Mass += Mx.GetWeight(k);
        
        CumulWeights.At(k, 0) = Mass;
    }
    
    CumulWeights *= 1.0 / Mass;
    
    // Build indexes to indicate which PDF of the mixture is used to generate each pattern
    
    for (size_t k = 0; k < n; k++)
    {
        double d = IndexeSample.At(k, 0);
        bool Continue = true;
        size_t Id = 0;
        
        while (Continue)
        {
            if ((CumulWeights.At(Id, 0) >= d) || (Id == nbPDF - 1))
            {
                Continue = false;
            }
            else
            {
                Id++;
            }
        }
        
        Pop.IncreaseElement(Id, 0, 1);
    }
    
    // Create patterns
    
    size_t PatternIndex = 0;
    
    for (size_t p = 0; p < nbPDF; p++)
    {
        size_t SubSamplePop = Pop.At(p, 0);
        
        MatrixDouble SubSample = NewGaussianSample(Mx.GetMember(p), SubSamplePop, m, reseed);
        
        for (size_t r = 0; r < SubSamplePop; r++)
        {
            for (size_t c = 0; c < dim; c++)
                Patterns.At(PatternIndex, c) = SubSample.At(r, c);
            
            PatternIndex++;
        }
    }
    return Patterns;
}