prob.h

Go to the documentation of this file.

/*  This file is part of libDAI - http://www.libdai.org/
 *
 *  Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
 *
 *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
 */




#ifndef __defined_libdai_prob_h
#define __defined_libdai_prob_h


#include <cmath>
#include <vector>
#include <ostream>
#include <algorithm>
#include <numeric>
#include <functional>
#include <dai/util.h>
#include <dai/exceptions.h>


namespace dai {


template<typename T> struct fo_id : public std::unary_function<T, T> {
    T operator()( const T &x ) const {
        return x;
    }
};


template<typename T> struct fo_abs : public std::unary_function<T, T> {
    T operator()( const T &x ) const {
        if( x < (T)0 )
            return -x;
        else
            return x;
    }
};


template<typename T> struct fo_exp : public std::unary_function<T, T> {
    T operator()( const T &x ) const {
        return exp( x );
    }
};


template<typename T> struct fo_log : public std::unary_function<T, T> {
    T operator()( const T &x ) const {
        return log( x );
    }
};


template<typename T> struct fo_log0 : public std::unary_function<T, T> {
    T operator()( const T &x ) const {
        if( x )
            return log( x );
        else
            return 0;
    }
};


template<typename T> struct fo_inv : public std::unary_function<T, T> {
    T operator()( const T &x ) const {
        return 1 / x;
    }
};


template<typename T> struct fo_inv0 : public std::unary_function<T, T> {
    T operator()( const T &x ) const {
        if( x )
            return 1 / x;
        else
            return 0;
    }
};


template<typename T> struct fo_plog0p : public std::unary_function<T, T> {
    T operator()( const T &p ) const {
        return p * dai::log0(p);
    }
};


template<typename T> struct fo_divides0 : public std::binary_function<T, T, T> {
    T operator()( const T &x, const T &y ) const {
        if( y == (T)0 )
            return (T)0;
        else
            return x / y;
    }
};


template<typename T> struct fo_KL : public std::binary_function<T, T, T> {
    T operator()( const T &p, const T &q ) const {
        if( p == (T)0 )
            return (T)0;
        else
            return p * (log(p) - log(q));
    }
};


template<typename T> struct fo_Hellinger : public std::binary_function<T, T, T> {
    T operator()( const T &p, const T &q ) const {
        T x = sqrt(p) - sqrt(q);
        return x * x;
    }
};


template<typename T> struct fo_pow : public std::binary_function<T, T, T> {
    T operator()( const T &x, const T &y ) const {
        if( y != 1 )
            return pow( x, y );
        else
            return x;
    }
};


template<typename T> struct fo_max : public std::binary_function<T, T, T> {
    T operator()( const T &x, const T &y ) const {
        return (x > y) ? x : y;
    }
};


template<typename T> struct fo_min : public std::binary_function<T, T, T> {
    T operator()( const T &x, const T &y ) const {
        return (x > y) ? y : x;
    }
};


template<typename T> struct fo_absdiff : public std::binary_function<T, T, T> {
    T operator()( const T &x, const T &y ) const {
        return dai::abs( x - y );
    }
};



template <typename T>
class TProb {
    public:
        typedef std::vector<T> container_type;

        typedef TProb<T> this_type;

    private:
        container_type _p;

    public:

        TProb() : _p() {}

        explicit TProb( size_t n ) : _p( n, (T)1 / n ) {}

        explicit TProb( size_t n, T p ) : _p( n, p ) {}


        template <typename TIterator>
        TProb( TIterator begin, TIterator end, size_t sizeHint ) : _p() {
            _p.reserve( sizeHint );
            _p.insert( _p.begin(), begin, end );
        }


        template <typename S>
        TProb( const std::vector<S> &v ) : _p() {
            _p.reserve( v.size() );
            _p.insert( _p.begin(), v.begin(), v.end() );
        }

        typedef typename container_type::const_iterator const_iterator;
        typedef typename container_type::iterator iterator;
        typedef typename container_type::const_reverse_iterator const_reverse_iterator;
        typedef typename container_type::reverse_iterator reverse_iterator;


        iterator begin() { return _p.begin(); }
        const_iterator begin() const { return _p.begin(); }

        iterator end() { return _p.end(); }
        const_iterator end() const { return _p.end(); }

        reverse_iterator rbegin() { return _p.rbegin(); }
        const_reverse_iterator rbegin() const { return _p.rbegin(); }

        reverse_iterator rend() { return _p.rend(); }
        const_reverse_iterator rend() const { return _p.rend(); }


        void resize( size_t sz ) {
            _p.resize( sz );
        }


        T get( size_t i ) const { 
#ifdef DAI_DEBUG
            return _p.at(i);
#else
            return _p[i];
#endif
        }

        void set( size_t i, T val ) {
            DAI_DEBASSERT( i < _p.size() );
            _p[i] = val;
        }


        const container_type& p() const { return _p; }

        container_type& p() { return _p; }

        T operator[]( size_t i ) const { return get(i); }

        size_t size() const { return _p.size(); }


        template<typename unOp> T accumulateSum( T init, unOp op ) const {
            T t = op(init);
            for( const_iterator it = begin(); it != end(); it++ )
                t += op(*it);
            return t;
        }


        template<typename unOp> T accumulateMax( T init, unOp op, bool minimize ) const {
            T t = op(init);
            if( minimize ) {
                for( const_iterator it = begin(); it != end(); it++ )
                    t = std::min( t, op(*it) );
            } else {
                for( const_iterator it = begin(); it != end(); it++ )
                    t = std::max( t, op(*it) );
            }
            return t;
        }

        T entropy() const { return -accumulateSum( (T)0, fo_plog0p<T>() ); }

        T max() const { return accumulateMax( (T)(-INFINITY), fo_id<T>(), false ); }

        T min() const { return accumulateMax( (T)INFINITY, fo_id<T>(), true ); }

        T sum() const { return accumulateSum( (T)0, fo_id<T>() ); }

        T sumAbs() const { return accumulateSum( (T)0, fo_abs<T>() ); }

        T maxAbs() const { return accumulateMax( (T)0, fo_abs<T>(), false ); }

        bool hasNaNs() const {
            bool foundnan = false;
            for( const_iterator x = _p.begin(); x != _p.end(); x++ )
                if( dai::isnan( *x ) ) {
                    foundnan = true;
                    break;
                }
            return foundnan;
        }

        bool hasNegatives() const {
            return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<T>(), (T)0 ) ) != _p.end());
        }

        std::pair<size_t,T> argmax() const {
            T max = _p[0];
            size_t arg = 0;
            for( size_t i = 1; i < size(); i++ ) {
              if( _p[i] > max ) {
                max = _p[i];
                arg = i;
              }
            }
            return std::make_pair( arg, max );
        }

        size_t draw() {
            Real x = rnd_uniform() * sum();
            T s = 0;
            for( size_t i = 0; i < size(); i++ ) {
                s += get(i);
                if( s > x )
                    return i;
            }
            return( size() - 1 );
        }


        bool operator<( const this_type& q ) const {
            DAI_DEBASSERT( size() == q.size() );
            return lexicographical_compare( begin(), end(), q.begin(), q.end() );
        }

        bool operator==( const this_type& q ) const {
            if( size() != q.size() )
                return false;
            return p() == q.p();
        }

        std::string toString() const {
            std::stringstream ss;
            ss << *this;
            return ss.str();
        }


        template<typename unaryOp> this_type pwUnaryTr( unaryOp op ) const {
            this_type r;
            r._p.reserve( size() );
            std::transform( _p.begin(), _p.end(), back_inserter( r._p ), op );
            return r;
        }

        this_type operator- () const { return pwUnaryTr( std::negate<T>() ); }

        this_type abs() const { return pwUnaryTr( fo_abs<T>() ); }

        this_type exp() const { return pwUnaryTr( fo_exp<T>() ); }


        this_type log(bool zero=false) const {
            if( zero )
                return pwUnaryTr( fo_log0<T>() );
            else
                return pwUnaryTr( fo_log<T>() );
        }


        this_type inverse(bool zero=true) const {
            if( zero )
                return pwUnaryTr( fo_inv0<T>() );
            else
                return pwUnaryTr( fo_inv<T>() );
        }


        this_type normalized( ProbNormType norm = dai::NORMPROB ) const {
            T Z = 0;
            if( norm == dai::NORMPROB )
                Z = sum();
            else if( norm == dai::NORMLINF )
                Z = maxAbs();
            if( Z == (T)0 ) {
                DAI_THROW(NOT_NORMALIZABLE);
                return *this;
            } else
                return pwUnaryTr( std::bind2nd( std::divides<T>(), Z ) );
        }


        template<typename unaryOp> this_type& pwUnaryOp( unaryOp op ) {
            std::transform( _p.begin(), _p.end(), _p.begin(), op );
            return *this;
        }

        this_type& randomize() {
            std::generate( _p.begin(), _p.end(), rnd_uniform );
            return *this;
        }

        this_type& setUniform () {
            fill( (T)1 / size() );
            return *this;
        }

        this_type& takeAbs() { return pwUnaryOp( fo_abs<T>() ); }

        this_type& takeExp() { return pwUnaryOp( fo_exp<T>() ); }


        this_type& takeLog(bool zero=false) {
            if( zero ) {
                return pwUnaryOp( fo_log0<T>() );
            } else
                return pwUnaryOp( fo_log<T>() );
        }


        T normalize( ProbNormType norm=dai::NORMPROB ) {
            T Z = 0;
            if( norm == dai::NORMPROB )
                Z = sum();
            else if( norm == dai::NORMLINF )
                Z = maxAbs();
            if( Z == (T)0 )
                DAI_THROW(NOT_NORMALIZABLE);
            else
                *this /= Z;
            return Z;
        }


        this_type& fill( T x ) {
            std::fill( _p.begin(), _p.end(), x );
            return *this;
        }

        this_type& operator+= (T x) {
            if( x != 0 )
                return pwUnaryOp( std::bind2nd( std::plus<T>(), x ) );
            else
                return *this;
        }

        this_type& operator-= (T x) {
            if( x != 0 )
                return pwUnaryOp( std::bind2nd( std::minus<T>(), x ) );
            else
                return *this;
        }

        this_type& operator*= (T x) {
            if( x != 1 )
                return pwUnaryOp( std::bind2nd( std::multiplies<T>(), x ) );
            else
                return *this;
        }

        this_type& operator/= (T x) {
            if( x != 1 )
                return pwUnaryOp( std::bind2nd( fo_divides0<T>(), x ) );
            else
                return *this;
        }

        this_type& operator^= (T x) {
            if( x != (T)1 )
                return pwUnaryOp( std::bind2nd( fo_pow<T>(), x) );
            else
                return *this;
        }


        this_type operator+ (T x) const { return pwUnaryTr( std::bind2nd( std::plus<T>(), x ) ); }

        this_type operator- (T x) const { return pwUnaryTr( std::bind2nd( std::minus<T>(), x ) ); }

        this_type operator* (T x) const { return pwUnaryTr( std::bind2nd( std::multiplies<T>(), x ) ); }

        this_type operator/ (T x) const { return pwUnaryTr( std::bind2nd( fo_divides0<T>(), x ) ); }

        this_type operator^ (T x) const { return pwUnaryTr( std::bind2nd( fo_pow<T>(), x ) ); }



        template<typename binaryOp> this_type& pwBinaryOp( const this_type &q, binaryOp op ) {
            DAI_DEBASSERT( size() == q.size() );
            std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), op );
            return *this;
        }


        this_type& operator+= (const this_type & q) { return pwBinaryOp( q, std::plus<T>() ); }


        this_type& operator-= (const this_type & q) { return pwBinaryOp( q, std::minus<T>() ); }


        this_type& operator*= (const this_type & q) { return pwBinaryOp( q, std::multiplies<T>() ); }


        this_type& operator/= (const this_type & q) { return pwBinaryOp( q, fo_divides0<T>() ); }


        this_type& divide (const this_type & q) { return pwBinaryOp( q, std::divides<T>() ); }


        this_type& operator^= (const this_type & q) { return pwBinaryOp( q, fo_pow<T>() ); }



        template<typename binaryOp> this_type pwBinaryTr( const this_type &q, binaryOp op ) const {
            DAI_DEBASSERT( size() == q.size() );
            TProb<T> r;
            r._p.reserve( size() );
            std::transform( _p.begin(), _p.end(), q._p.begin(), back_inserter( r._p ), op );
            return r;
        }


        this_type operator+ ( const this_type& q ) const { return pwBinaryTr( q, std::plus<T>() ); }


        this_type operator- ( const this_type& q ) const { return pwBinaryTr( q, std::minus<T>() ); }


        this_type operator* ( const this_type &q ) const { return pwBinaryTr( q, std::multiplies<T>() ); }


        this_type operator/ ( const this_type &q ) const { return pwBinaryTr( q, fo_divides0<T>() ); }


        this_type divided_by( const this_type &q ) const { return pwBinaryTr( q, std::divides<T>() ); }


        this_type operator^ ( const this_type &q ) const { return pwBinaryTr( q, fo_pow<T>() ); }


        template<typename binOp1, typename binOp2> T innerProduct( const this_type &q, T init, binOp1 binaryOp1, binOp2 binaryOp2 ) const {
            DAI_DEBASSERT( size() == q.size() );
            return std::inner_product( begin(), end(), q.begin(), init, binaryOp1, binaryOp2 );
        }
};



template<typename T> T dist( const TProb<T> &p, const TProb<T> &q, ProbDistType dt ) {
    switch( dt ) {
        case DISTL1:
            return p.innerProduct( q, (T)0, std::plus<T>(), fo_absdiff<T>() );
        case DISTLINF:
            return p.innerProduct( q, (T)0, fo_max<T>(), fo_absdiff<T>() );
        case DISTTV:
            return p.innerProduct( q, (T)0, std::plus<T>(), fo_absdiff<T>() ) / 2;
        case DISTKL:
            return p.innerProduct( q, (T)0, std::plus<T>(), fo_KL<T>() );
        case DISTHEL:
            return p.innerProduct( q, (T)0, std::plus<T>(), fo_Hellinger<T>() ) / 2;
        default:
            DAI_THROW(UNKNOWN_ENUM_VALUE);
            return INFINITY;
    }
}



template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& p) {
    os << "(";
    for( size_t i = 0; i < p.size(); i++ )
        os << ((i != 0) ? ", " : "") << p.get(i);
    os << ")";
    return os;
}



template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
    return a.pwBinaryTr( b, fo_min<T>() );
}



template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
    return a.pwBinaryTr( b, fo_max<T>() );
}


typedef TProb<Real> Prob;


} // end of namespace dai


#endif