IncrementalMoments.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_INCREMENTALMOMENTS_H_
#define _MIRA_INCREMENTALMOMENTS_H_

#include <math/Eigen.h>

namespace mira {


template <typename T, int D, int Order>
class IncrementalMoments
{
    static_assert(sizeof(T)==0, "There is no implementation for the moment of the given dimension D and Order");
};


// 1D cases
template <typename T>
class IncrementalMoments<T,1,1>
{
public:

    IncrementalMoments() {
        reset();
    }

    void reset()
    {
        k = T(0);
        m1 = T(0);
    }

    void insert(const T& x, T weight = T(1))
    {
        Ctx ctx;
        _insert(x, weight, ctx);
    }

public:

    T sumWeight() const { return k; }

    T mean() const { return m1; }

protected:

    struct Ctx {
        T k;
        T nk;
        T alpha;
        T beta;
    };

    T k;
    T m1;

protected:

    void _insert(const T& x, T weight, Ctx& ctx)
    {
        ctx.k = k;
        ctx.nk = k+weight;
        ctx.alpha = k/ctx.nk;      // = k / (k+weight);
        ctx.beta  = weight/ctx.nk; // = weight / (k+weight);

        // update mean
        m1 = ctx.alpha*m1 + ctx.beta*x;
        k = ctx.nk;
    }
};

template <typename T>
class IncrementalMoments<T,1,2> : public IncrementalMoments<T,1,1>
{
    typedef IncrementalMoments<T,1,1> Base;
    typedef typename Base::Ctx Ctx;

public:

    IncrementalMoments() {
        reset();
    }

    void reset()
    {
        Base::reset();
        m2 = T(0);
    }

    void insert(const T& x, T weight = T(1))
    {
        Ctx ctx;
        _insert(x, weight, ctx);
    }

    T var() const { return m2; }

protected:

    T m2;

protected:

    void _insert(const T& x, T weight, Ctx& ctx)
    {
        Base::_insert(x, weight, ctx);
        float alpha_beta = ctx.alpha*ctx.beta; // = (k*weight) / (k+weight);
        // update covariance
        T d = this->m1 - x;
        m2 = (ctx.alpha*m2) + alpha_beta * (d*d);
    }
};


// n-D cases
template <typename T, int D>
class IncrementalMoments<T,D,1>
{
public:
    typedef Eigen::Matrix<T,D,1> Vector;
    typedef Eigen::Matrix<T,D,D> Matrix;

public:

    IncrementalMoments() {
        reset();
    }

    void reset()
    {
        k = T(0);
        m1.setZero();
    }

    void insert(const Vector& x, T weight = T(1))
    {
        Ctx ctx;
        _insert(x, weight, ctx);
    }

public:

    T sumWeight() const { return k; }

    const Vector& mean() const { return m1; }

protected:

    struct Ctx {
        T k;
        T nk;
        T alpha;
        T beta;
    };

    T k;
    Vector m1;

protected:

    void _insert(const Vector& x, T weight, Ctx& ctx)
    {
        ctx.k = k;
        ctx.nk = k+weight;
        ctx.alpha = k/ctx.nk;      // = k / (k+weight);
        ctx.beta  = weight/ctx.nk; // = weight / (k+weight);

        // update mean
        m1 = ctx.alpha*m1 + ctx.beta*x;
        k = ctx.nk;
    }
};

template <typename T, int D>
class IncrementalMoments<T,D,2> : public IncrementalMoments<T,D,1>
{
    typedef IncrementalMoments<T,D,1> Base;
    typedef typename Base::Vector Vector;
    typedef typename Base::Matrix Matrix;
    typedef typename Base::Ctx Ctx;

public:

    IncrementalMoments() {
        reset();
    }

    void reset()
    {
        Base::reset();
        m2.setZero();
    }

    void insert(const Vector& x, T weight = T(1))
    {
        Ctx ctx;
        _insert(x, weight, ctx);
    }

    Matrix cov() const { return m2; }

protected:

    Matrix m2;

protected:

    void _insert(const Vector& x, T weight, Ctx& ctx)
    {
        Base::_insert(x, weight, ctx);
        float alpha_beta = ctx.alpha*ctx.beta; // = (k*weight) / (k+weight);
        // update covariance
        Vector d = this->m1 - x;
        m2 = (ctx.alpha*m2) + alpha_beta * (d*d.transpose());
    }
};




} // namespace

#endif