prob_util.h

Go to the documentation of this file.

/* $Id: prob_util.h 18931 2015-04-22 21:50:53Z cdawson $ */
/* ======================================================================== *
 |                                                                          |
 | Copyright (c) 2007-2010, by members of University of Arizona Computer    |
 | Vision group (the authors) including                                     |
 |       Kobus Barnard, Kyle Simek, Ernesto Brau.                           |
 |                                                                          |
 | For use outside the University of Arizona Computer Vision group please   |
 | contact Kobus Barnard.                                                   |
 |                                                                          |
 * ======================================================================== */


#ifndef PROB_UTIL_H_INCLUDED
#define PROB_UTIL_H_INCLUDED

#include <iterator>
#include <algorithm>
#include <cmath>

namespace kjb{

    //forward declarations
    class Vector;
    class Matrix;
    
#if 0
NOTE: I moved this to sample_cpp to resolve circular dependency between prob_util.h and prob_distribution.h.  Besides, this is basically sampling, even
  though its name isnt sample()...
  Kyle May 19, 2012
template<class Iterator, class distance_type>
inline
Iterator element_uar(Iterator first, distance_type size)
{
    int dist = kjb::sample(kjb::Categorical_distribution<size_t>(0, size - 1, 1));
    Iterator p = first;
    std::advance(p, dist);

    return p;
}
#endif


inline double log_sum(double n1, double n2)
{
/*
 * the implementation below relies on this identity:
 *  
 * log{ x1 + x2 + x3 + ...} =
 *    pi + log(sum_i exp( log xi - pi))
 *
 * where pi = max_i (log(xi))
 */
    if(n1 < n2)
        std::swap(n1, n2);

    return n1 + log(1.0 + exp(n2-n1));
}

inline double log_sum(double n1, double n2, double n3)
{
/*
 * the implementation below relies on this identity:
 *  
 * log{ x1 + x2 + x3 + ...} =
 *    pi + log(sum_i exp( log xi - pi))
 *
 * where pi = max_i (log(xi))
 */
    if(n1 < n2)
        std::swap(n1, n2);
    if(n1 < n3)
        std::swap(n1, n3);

    return n1 + log(1.0 + exp(n2-n1) + exp(n3 - n1));
}

template<class Iterator>
double log_sum(Iterator first, Iterator last)
{
/*
 * the implementation below relies on this identity:
 *  
 * log{ x1 + x2 + x3 + ...} =
 *    pi + log(sum_i exp( log xi - pi))
 *
 * where pi = max_i (log(xi))
 *
 */
    double pi = *std::max_element(first, last);

    double accum = 0;

    for(; first != last; first++)
    {
        accum += exp(*first - pi);
    }
    return pi + log(accum);
}

inline double log_diff(double n1, double n2)
{
/*
 * the implementation below relies on this identity:
 *  
 * log{ x1 - x2} = log(x2) + log(exp(log(x2/x1) - 1))
 *
 */
    if(n1 < n2)
        std::swap(n1, n2);

    return n2 + log(exp(n2-n1) - 1);
}
    
Vector log_normalize(const Vector& vec);
    
Vector exponentiate(const Vector& vec);
    
Vector log_normalize_and_exponentiate(const Vector& vec);
    
Matrix log_normalize_rows(const Matrix& mat);
    
Matrix log_normalize_and_exponentiate(const Matrix& mat);

Vector log_marginalize_over_cols(const Matrix& mat);
    
Vector log_marginalize_over_rows(const Matrix& mat);
    
}; // namespace kjb

#endif /*PROB_UTIL_H_INCLUDED */