Log.hpp

/*Copyright 2009,2010 Alex Graves

This file is part of RNNLIB.

RNNLIB is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

RNNLIB 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 General Public License for more details.

You should have received a copy of the GNU General Public License
along with RNNLIB.  If not, see <http://www.gnu.org/licenses/>.*/

#ifndef _INCLUDED_Log_h  
#define _INCLUDED_Log_h

#include <cmath>
#include <iostream>
#include <limits>

using namespace std;

template<class T> class Log
{
	//data
	T expVal;
	T logVal;
	
public:

	//static data
	static const T expMax;
	static const T expMin;
	static const T expLimit;
	static const T logZero;
	static const T logInfinity;
		
	//static functions
	static T safe_exp(T x)
	{
		if (x == logZero)
		{
			return 0;
		}
		if (x >= expLimit)
		{
			return expMax;
		}
		return std::exp(x);
	}
	static T safe_log(T x)
	{
		if (x <= expMin)
		{
			return logZero;
		}
		return std::log(x);
	}
	static T log_add(T x, T y)
	{
		if (x == logZero)
		{
			return y;
		}
		if (y == logZero)
		{
			return x;
		}
		if (x < y) 
		{
			swap(x, y);
		}
		return x + std::log(1.0 + safe_exp(y - x));
	}
	static T log_subtract(T x, T y)
	{
		if (y == logZero)
		{
			return x;
		}	
		if (y >= x)
		{
			return logZero;
		}
		return x + std::log(1.0 - safe_exp(y - x));
	}
	static T log_multiply(T x, T y)
	{
		if (x == logZero || y == logZero)
		{
			return logZero;
		}
		return x + y;
	}
	static T log_divide(T x, T y)
	{
		if (x == logZero)
		{
			return logZero;
		}
		if (y == logZero)
		{
			return logInfinity;
		}
		return x - y;
	}
	
	//functions
	Log(T v = 0, bool logScale = false):
		expVal(logScale ? -1 : v),
		logVal(logScale ? v : safe_log(v))
	{}
	Log<T>& operator =(const Log<T>& l)
	{
		logVal = l.logVal;
		expVal = l.expVal;
		return *this;
	}
	Log<T>& operator +=(const Log<T>& l)
	{
		logVal = log_add(logVal, l.logVal);
		expVal = -1;
		return *this;
	}
	Log<T>& operator -=(const Log<T>& l)
	{
		logVal = log_subtract(logVal, l.logVal);
		expVal = -1;
		return *this;
	}
	Log<T>& operator *=(const Log<T>& l)
	{
		logVal = log_multiply(logVal, l.logVal);
		expVal = -1;
		return *this;
	}
	Log<T>& operator /=(const Log<T>& l)
	{
		logVal = log_divide(logVal, l.logVal);
		expVal = -1;
		return *this;
	}
	T exp()
	{
		if (expVal < 0)
		{
			expVal = safe_exp(logVal);
		}
		return expVal; 
	}
	T log() const
	{
		return logVal;
	}
};

//helper functions
template<class T> Log<T> operator +(Log<T> log1, Log<T> log2)
{
	return Log<T>(Log<T>::log_add(log1.log(), log2.log()), true);
}
template<class T> Log<T> operator -(Log<T> log1, Log<T> log2)
{
	return Log<T>(Log<T>::log_subtract(log1.log(), log2.log()), true);
}
template<class T> Log<T> operator *(Log<T> log1, Log<T> log2)
{
	return Log<T>(Log<T>::log_multiply(log1.log(), log2.log()), true);
}
template<class T> Log<T> operator /(Log<T> log1, Log<T> log2)
{
	return Log<T>(Log<T>::log_divide(log1.log(), log2.log()), true);
}
template<class T> bool operator >(Log<T> log1, Log<T> log2)
{
	return (log1.log() > log2.log());
}
template<class T> bool operator <(Log<T> log1, Log<T> log2)
{
	return (log1.log() < log2.log());
}
template<class T> bool operator ==(Log<T> log1, Log<T> log2)
{
	return (log1.log() == log2.log());
}
template<class T> bool operator <=(Log<T> log1, Log<T> log2)
{
	return (log1.log() <= log2.log());
}
template<class T> bool operator >=(Log<T> log1, Log<T> log2)
{
	return (log1.log() >= log2.log());
}
template<class T> ostream& operator <<(ostream& out, const Log<T>& l)
{
	out << l.log();
	return out;
}
template<class T> istream& operator >>(istream& in, Log<T>& l)
{
	T d;
	in >> d;
	l = Log<T>(d, true);
	return in;
}

template <class T> const T Log<T>::expMax = numeric_limits<T>::max();
template <class T> const T Log<T>::expMin = numeric_limits<T>::min();
template <class T> const T Log<T>::expLimit = std::log(expMax);
template <class T> const T Log<T>::logInfinity = 1e100;
template <class T> const T Log<T>::logZero = -Log<T>::logInfinity;

#endif