logfloat-0.8.3: Data/Number/LogFloat.hs
-- %% This module should be run through lhs2hs before running through
-- %% Haddock. (N.B. rember to include a copy in the cabalized)
-- %%
-- %% This module was originally translated from my Perl module
-- %% Math::LogFloat (version 0.3; revision 2007.12.20)
-- %%
-- %% N.B. Can't have `#' in the first column in GHC, not even if lhs
--
-- TODO: Add QuickCheck-ness, though beware of the fuzz.
-- TODO: Make sure rewrite rules really fire
-- TODO: profile to make sure we don't waste too much time constructing dictionaries
--
-- To turn on optimizations and look at the optimization records, cf:
-- http://www.haskell.org/ghc/docs/latest/html/users_guide/rewrite-rules.html
-- http://www.randomhacks.net/articles/2007/02/10/map-fusion-and-haskell-performance
-- {-# OPTIONS_GHC -ddump-simpl-stats #-}
{-# OPTIONS_GHC -O2 -fvia-C -optc-O3 #-}
-- Version History
-- (v0.8) Did a bunch of tweaking. Things should be decent now
-- (v0.7) Haddockified
-- (v0.6) Fixed monomorphism.
-- (v0.5) Added optimization rules.
-- (v0.4) Translated to Haskell at revision 2007.12.20.
-- (v0.3) Converted extensive comments to POD format.
-- (v0.2) Did a bunch of profiling, optimizing, and debugging.
-- (v0.1) Initial version created for hw5 for NLP with Jason Eisner.
--
----------------------------------------------------------------
-- ~ 2008.08.15
-- |
-- Module : Data.Number.LogFloat
-- Copyright : Copyright (c) 2007--2008 wren ng thornton
-- License : BSD3
-- Maintainer : wren@community.haskell.org
-- Stability : stable
-- Portability : portable
--
-- This module presents a class for storing numbers in the log-domain.
-- The main reason for doing this is to prevent underflow when multiplying
-- many small probabilities as is done in Hidden Markov Models and
-- other statistical models often used for natural language processing.
-- The log-domain also helps prevent overflow when multiplying many
-- large numbers. In rare cases it can speed up numerical computation
-- (since addition is faster than multiplication, though logarithms
-- are exceptionally slow), but the primary goal is to improve accuracy
-- of results. A secondary goal has been to maximize efficiency since
-- these computations are frequently done within a /O(n^3)/ loop.
--
-- The 'LogFloat' of this module is restricted to non-negative numbers
-- for efficiency's sake, see the forthcoming "Data.Number.LogFloat.Signed"
-- for doing signed log-domain calculations.
----------------------------------------------------------------
module Data.Number.LogFloat
(
-- * Documentation Note
-- | If you see no module description above, then the @lhs2hs@
-- script was not run correctly. Please rebuild the documentation
-- or see:
-- <http://code.haskell.org/~wren/logfloat/dist/doc/html/logfloat/>
-- * IEEE floating-point special values
-- | "GHC.Real" defines 'infinity' and 'notANumber' as
-- 'Rational'. We export variants which are polymorphic because
-- that can be more helpful at times.
--
-- BUG: At present these constants are broken for @Ratio@
-- types including 'Rational', since @Ratio@ types do not
-- typically permit a zero denominator. In GHC (6.8.2) the
-- result for 'infinity' is a rational with a numerator
-- sufficiently large that 'fromRational' will yield infinity
-- for @Float@ and @Double@. In Hugs (September 2006) it
-- yields an arithmetic overflow error. For GHC, our 'notANumber'
-- yields @0%1@ rather than @0%0@ as "GHC.Real" does.
infinity, negativeInfinity, notANumber
-- * Basic functions
, log, toFractional
-- * @LogFloat@ data type and conversion functions
, LogFloat
, logFloat, logToLogFloat
, fromLogFloat, logFromLogFloat
) where
import Prelude hiding (log)
import qualified Prelude (log)
-- Not portable, and we can do it ourselves.
-- import qualified GHC.Real (infinity, notANumber)
----------------------------------------------------------------
--
-- Try to add in some optimizations. Why the first few need to be down
-- here and localized to the module, I don't know. We don't do anything
-- foolish like this, but our clients might or they might be generated
-- by other code transformations.
{-# RULES
"log/exp" forall x. log (exp x) = x
"log.exp" log . exp = id
"exp/log" forall x. exp (log x) = x
"exp.log" exp . log = id
#-}
-- These are general rule versions of our operators for 'LogFloat'. I
-- had some issues inducing 'Ord' on @x@ and @y@, even though they're
-- 'Num' so I can't do "(+)/log" and "(-)/log" so easily.
{-# RULES
"(*)/log" forall x y. log x * log y = log (x + y)
"(/)/log" forall x y. log x / log y = log (x - y)
#-}
----------------------------------------------------------------
--
-- The type signature is necessary for them not to default to Double.
infinity, negativeInfinity, notANumber :: (Fractional a) => a
infinity = toFractional (1/0) -- == fromRational GHC.Real.infinity
{-# SPECIALIZE negativeInfinity :: Double #-}
negativeInfinity = negate infinity
notANumber = infinity - infinity -- == fromRational GHC.Real.notANumber
-- The dictionaries for these are really ugly in core.
-- TODO: be sure to check that these don't give eggregious performance hits
--
----------------------------------------------------------------
--
-- | Since the normal 'Prelude.log' throws an error on zero, we have
-- to redefine it in order for things to work right. Arguing from
-- limits it's obvious that @log 0 == negativeInfinity@.
--
-- If you're using some 'Floating' type that's not built in, verify
-- this equation holds for your @0@ and @negativeInfinity@. If it
-- doesn't, then you should avoid importing our 'log' and will probably
-- want converters to handle the discrepency.
{-# SPECIALIZE log :: Double -> Double #-}
log :: (Floating a) => a -> a
log 0 = negativeInfinity
log x = Prelude.log x
-- | The most generic numeric converter I can come up with. All the
-- built-in numeric types are 'Real', though 'Int' and 'Integer' aren't
-- 'Fractional'.
{-# SPECIALIZE toFractional :: (Real a) => a -> Double #-}
{-# SPECIALIZE toFractional :: (Fractional b) => Double -> b #-}
toFractional :: (Real a, Fractional b) => a -> b
toFractional = fromRational . toRational
-- This should only fire when it's type-safe
{-# RULES "toFractional/id" toFractional = id #-}
-- This should happen already, but who knows
-- TODO: see if it ever fires
{-# RULES
"toFractional/toFractional" forall x.
toFractional (toFractional x) = toFractional x
"toFractional.toFractional" toFractional . toFractional = toFractional
#-}
----------------------------------------------------------------
--
-- | Reduce the number of constant string literals we need to store.
errorOutOfRange :: String -> a
errorOutOfRange fun = error $ "Data.Number.LogFloat."++fun
++ ": argument out of range"
-- | We need these guards in order to ensure some invariants.
guardNonNegative :: String -> Double -> Double
guardNonNegative fun x | x >= 0 = x
| otherwise = errorOutOfRange fun
-- | It's unfortunate that notANumber is not equal to itself, but we
-- can hack around that. Is there any efficiency difference between
-- these two tests? If not, then we could use @log . guardNonNegative
-- fun = guardIsANumber fun . log@ in order to remove guardNonNegative.
guardIsANumber :: String -> Double -> Double
guardIsANumber fun x | x >= negativeInfinity = x
| otherwise = errorOutOfRange fun
----------------------------------------------------------------
--
-- | A @LogFloat@ is just a 'Double' with a special interpretation.
-- The 'logFloat' function is presented instead of the constructor,
-- in order to ensure semantic conversion. At present the 'Show'
-- instance will convert back to the normal-domain, and so will underflow
-- at that point. This behavior may change in the future.
--
-- Performing operations in the log-domain is cheap, prevents underflow,
-- and is otherwise very nice for dealing with miniscule probabilities.
-- However, crossing into and out of the log-domain is expensive and
-- should be avoided as much as possible. In particular, if you're
-- doing a series of multiplications as in @lp * logFloat q * logFloat
-- r@ it's faster to do @lp * logFloat (q * r)@ if you're reasonably
-- sure the normal-domain multiplication won't underflow, because that
-- way you enter the log-domain only once, instead of twice.
--
-- Even more particularly, you should /avoid addition/ whenever possible.
-- Addition is provided because it's necessary at times and the proper
-- implementation is not immediately transparent. However, between two
-- @LogFloat@s addition requires crossing the exp/log boundary twice;
-- with a @LogFloat@ and a regular number it's three times since the
-- regular number needs to enter the log-domain first. This makes addition
-- incredibly slow. Again, if you can parenthesize to do plain operations
-- first, do it!
newtype LogFloat = LogFloat Double
deriving (Eq, Ord)
-- | A constructor which does semantic conversion from normal-domain
-- to log-domain.
{-# SPECIALIZE logFloat :: Double -> LogFloat #-}
logFloat :: (Real a) => a -> LogFloat
logFloat = LogFloat . log . guardNonNegative "logFloat" . toFractional
-- This is simply a polymorphic version of the 'LogFloat' data
-- constructor. We present it mainly because we hide the constructor
-- in order to make the type a bit more opaque. If the polymorphism
-- turns out to be a performance liability because the rewrite rules
-- can't remove it, then we need to rethink all four constructors/destructors.
--
-- | Constructor which assumes the argument is already in the log-domain.
{-# SPECIALIZE logToLogFloat :: Double -> LogFloat #-}
logToLogFloat :: (Real a) => a -> LogFloat
logToLogFloat = LogFloat . guardIsANumber "logToLogFloat" . toFractional
-- | Return our log-domain value back into normal-domain. Beware of
-- overflow/underflow.
{-# SPECIALIZE fromLogFloat :: LogFloat -> Double #-}
fromLogFloat :: (Floating a) => LogFloat -> a
fromLogFloat (LogFloat x) = toFractional (exp x)
-- | Return the log-domain value itself without costly conversion
{-# SPECIALIZE logFromLogFloat :: LogFloat -> Double #-}
logFromLogFloat :: (Floating a) => LogFloat -> a
logFromLogFloat (LogFloat x) = toFractional x
-- These are our module-specific versions of "log/exp" and "exp/log";
-- They do the same things but also have a @LogFloat@ in between the
-- logarithm and exponentiation.
{-# RULES
-- Out of log-domain and back in
"log/fromLogFloat" forall x. log (fromLogFloat x) = logFromLogFloat x
"log.fromLogFloat" log . fromLogFloat = logFromLogFloat
"logFloat/fromLogFloat" forall x. logFloat (fromLogFloat x) = x
"logFloat.fromLogFloat" logFloat . fromLogFloat = id
-- Into log-domain and back out
"fromLogFloat/logFloat" forall x. fromLogFloat (logFloat x) = x
"fromLogFloat.logFloat" fromLogFloat . logFloat = id
#-}
----------------------------------------------------------------
-- To show it, we want to show the normal-domain value rather than the
-- log-domain value. Also, if someone managed to break our invariants
-- (e.g. by passing in a negative and noone's pulled on the thunk yet)
-- then we want to crash before printing the constructor, rather than
-- after. N.B. This means the show will underflow/overflow in the
-- same places as normal doubles since we underflow at the exp. Perhaps
-- this means we should show the log-domain value instead.
instance Show LogFloat where
show (LogFloat x) = let y = exp x
in y `seq` "LogFloat "++show y
----------------------------------------------------------------
-- These all work without causing underflow. However, do note that
-- they tend to induce more of the floating-point fuzz than using
-- regular floating numbers because @exp . log@ doesn't really equal
-- @id@. In any case, our main aim is for preventing underflow when
-- multiplying many small numbers (and preventing overflow for multiplying
-- many large numbers) so we're not too worried about +/- 4e-16.
instance Num LogFloat where
(*) (LogFloat x) (LogFloat y) = LogFloat (x+y)
(+) (LogFloat x) (LogFloat y)
| x >= y = LogFloat (x + log (1 + exp (y - x)))
| otherwise = LogFloat (y + log (1 + exp (x - y)))
-- Without the guard this would return NaN instead of error
(-) (LogFloat x) (LogFloat y)
| x >= y = LogFloat (x + log (1 - exp (y - x)))
| otherwise = errorOutOfRange "(-)"
signum (LogFloat x)
| x == negativeInfinity = 0
| x > negativeInfinity = 1
| otherwise = errorOutOfRange "signum"
-- The extra guard protects against NaN, in case someone
-- broke the invariant. That shouldn't be possible and
-- so noone else bothers to check, but we check here just
-- in case.
negate _ = errorOutOfRange "negate"
abs = id
fromInteger = LogFloat . log
. guardNonNegative "fromInteger" . fromInteger
instance Fractional LogFloat where
-- n/0 is handled seamlessly for us; we must catch 0/0 though
(/) (LogFloat x) (LogFloat y)
| x == negativeInfinity
&& y == negativeInfinity = errorOutOfRange "(/)" -- protect vs NaN
| otherwise = LogFloat (x-y)
fromRational = LogFloat . log
. guardNonNegative "fromRational" . fromRational
-- Just for fun. The more coersion functions the better. Though
-- it can underflow...
instance Real LogFloat where
toRational (LogFloat x) = toRational (exp x)
----------------------------------------------------------------
-- ----------------------------------------------------------- fin.