NormalRandomGenerator.h

Go to the documentation of this file.

/*
 * Copyright (C) 2012 by
 *   MetraLabs GmbH (MLAB), GERMANY
 * and
 *   Neuroinformatics and Cognitive Robotics Labs (NICR) at TU Ilmenau, GERMANY
 * All rights reserved.
 *
 * Contact: info@mira-project.org
 *
 * Commercial Usage:
 *   Licensees holding valid commercial licenses may use this file in
 *   accordance with the commercial license agreement provided with the
 *   software or, alternatively, in accordance with the terms contained in
 *   a written agreement between you and MLAB or NICR.
 *
 * GNU General Public License Usage:
 *   Alternatively, this file may be used under the terms of the GNU
 *   General Public License version 3.0 as published by the Free Software
 *   Foundation and appearing in the file LICENSE.GPL3 included in the
 *   packaging of this file. Please review the following information to
 *   ensure the GNU General Public License version 3.0 requirements will be
 *   met: http://www.gnu.org/copyleft/gpl.html.
 *   Alternatively you may (at your option) use any later version of the GNU
 *   General Public License if such license has been publicly approved by
 *   MLAB and NICR (or its successors, if any).
 *
 * IN NO EVENT SHALL "MLAB" OR "NICR" BE LIABLE TO ANY PARTY FOR DIRECT,
 * INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF
 * THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF "MLAB" OR
 * "NICR" HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * "MLAB" AND "NICR" SPECIFICALLY DISCLAIM ANY WARRANTIES, INCLUDING,
 * BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
 * FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS
 * ON AN "AS IS" BASIS, AND "MLAB" AND "NICR" HAVE NO OBLIGATION TO
 * PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS OR MODIFICATIONS.
 */

#ifndef _MIRA_NORMALRANDOMGENERATOR_H_
#define _MIRA_NORMALRANDOMGENERATOR_H_

#include <math/Math.h>
#include <math/Eigen.h>

#include <utils/PParam.h>
#include <math/RandomGenerator.h>
#ifndef Q_MOC_RUN
#include <boost/random/uniform_01.hpp>
#endif

namespace mira {


namespace detail {
template <typename T>
class NormalRandomDistributionBase
{
public:

    NormalRandomDistributionBase() : first(true) {}

    template <typename Engine>
    T normal(Engine& eng)
    {
        if(first) {
            float u;
            u = boost::uniform_01<T>()(eng);
            v = boost::uniform_01<T>()(eng);
            sqrt_2ln_u = std::sqrt(-T(2) * std::log(T(1)-u));
            // using log(1-u) above instead of log(u), since u is element of [0,1)
            // i.e. u may become 0, where log(0) = -inf. When using log(1-u)
            // we can avoid the -inf.
        }
        first = !first;
        return sqrt_2ln_u * (first ? std::sin(T(2)*pi<T>()*v) :
                                     std::cos(T(2)*pi<T>()*v));
    }

private:
    bool first;
    T sqrt_2ln_u,v;
};
}


template <int D, typename T=float>
class NormalRandomDistribution : private detail::NormalRandomDistributionBase<T>
{
public:
    typedef T input_type;
    typedef Eigen::Matrix<T,D,1> result_type;
    typedef detail::NormalRandomDistributionBase<T> Base;

    typedef Eigen::Matrix<T,D,D> VarianceType;

public:

    NormalRandomDistribution() {
        mL = VarianceType::Identity();
    }

    NormalRandomDistribution(const Eigen::Matrix<T,D,D>& sigma) {
        setSigma(sigma);
    }

    bool setSigma(const Eigen::Matrix<T,D,D>& sigma) {
        // Cholesky decomposition: sigma = L*L'
        auto llt = sigma.llt();
        if(llt.info()!=Eigen::Success)
            return false;
        mL = llt.matrixL();
        return true;
    }

    template <typename Engine>
    Eigen::Matrix<T,D,1> operator()(Engine& eng) {
        result_type n;
        for(int i=0;i<D;++i)
            n(i) = this->normal(eng);
        return mL*n;
    }

private:
    Eigen::Matrix<T,D,D> mL;
};

// specialization for univariate distribution (D=1)
template <typename T>
class NormalRandomDistribution<1,T> : private detail::NormalRandomDistributionBase<T>
{
public:
    typedef T input_type;
    typedef T result_type;

    typedef T VarianceType;

public:

    NormalRandomDistribution() : mSigma(T(1)) {}

    NormalRandomDistribution(T sigma) : mSigma(sigma) {}

    void setSigma(T sigma) {
        mSigma = sigma;
    }


    template <typename Engine>
    T operator()(Engine& eng) {
        return mSigma * this->normal(eng);
    }

private:
    T mSigma;
};


template <int D, typename T=float>
class NormalRandomGenerator : public RandomGenerator<NormalRandomDistribution<D,T>>
{
    // if new C++11 typedefs are supported this could become a templated typedef.
public:
    NormalRandomGenerator() {}
    MIRA_RANDOM_GENERATOR_COMMON(NormalRandomGenerator, MIRA_PPARAM(NormalRandomDistribution<D,T>))
};





} // namespace

#endif