%% This module should be run through lhs2hs.pl 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 QuickCheckness, 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/rewriterules.html
http://www.randomhacks.net/articles/2007/02/10/mapfusionandhaskellperformance
>
>
>
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.01
|
Module : Data.Number.LogFloat
Copyright : Copyright (c) 2007
License : BSD3
Maintainer : wren@community.haskell.org
Stability : stable
Portability : portable
This module presents a class for storing numbers in the logdomain.
The main reason for doing this is to prevent underflow when multiplying
many probabilities as is done in Hidden Markov Models. It is also
helpful for preventing overflow. In certain rare cases it may speed
up computations (addition is faster than multiplication, but
logarithms are really slow), but the primary goal is to improve
accuracy of results.
The 'LogFloat' of this module is restricted to nonnegative numbers
for efficiency's sake, see the forthcoming "Data.Number.LogFloat.Signed"
for doing signed logdomain calculations.
> module Data.Number.LogFloat
> (
>
>
>
>
>
> infinity, negativeInfinity, notANumber
>
>
> , log, toFractional
>
>
> , LogFloat
> , logFloat, logToLogFloat
> , fromLogFloat, logFromLogFloat
> ) where
>
> import Prelude hiding (log)
> import qualified Prelude (log)
>
>
>
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.
>
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.
>
The type signature is necessary for them not to default to Double.
> infinity, negativeInfinity, notANumber :: (Fractional a) => a
> infinity = 1 / 0
>
> negativeInfinity = negate infinity
> notANumber = infinity infinity
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.
>
> 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
builtin numeric types are 'Real', though 'Int' and 'Integer' aren't
'Fractional'.
>
>
> toFractional :: (Real a, Fractional b) => a -> b
> toFractional = fromRational . toRational
>
>
>
>
>
>
>
| 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 normaldomain, and so will underflow
at that point. This behavior may change in the future.
Performing operations in the logdomain is cheap, prevents underflow,
and is otherwise very nice for dealing with miniscule probabilities.
However, crossing into and out of the logdomain 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 normaldomain multiplication won't underflow, because that
way you enter the logdomain 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 logdomain 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 normaldomain
to logdomain.
>
> 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 logdomain.
>
> logToLogFloat :: (Real a) => a -> LogFloat
> logToLogFloat = LogFloat . guardIsANumber "logToLogFloat" . toFractional
| Return our logdomain value back into normaldomain. Beware of
overflow/underflow.
>
> fromLogFloat :: (Floating a) => LogFloat -> a
> fromLogFloat (LogFloat x) = toFractional (exp x)
| Return the logdomain value itself without costly conversion
>
> logFromLogFloat :: (Floating a) => LogFloat -> a
> logFromLogFloat (LogFloat x) = toFractional x
These are our modulespecific versions of "log/exp" and "exp/log";
They do the same things but also have a @LogFloat@ in between the
logarithm and exponentiation.
>
To show it, we want to show the normaldomain value rather than the
logdomain 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 logdomain 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 floatingpoint 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)))
>
>
> () (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"
>
>
>
>
>
> negate _ = errorOutOfRange "negate"
>
> abs = id
>
> fromInteger = LogFloat . log
> . guardNonNegative "fromInteger" . fromInteger
>
>
> instance Fractional LogFloat where
>
> (/) (LogFloat x) (LogFloat y)
> | x == negativeInfinity
> && y == negativeInfinity = errorOutOfRange "(/)"
> | otherwise = LogFloat (xy)
>
> fromRational = LogFloat . log
> . guardNonNegative "fromRational" . fromRational
>
>
>
>
> instance Real LogFloat where
> toRational (LogFloat x) = toRational (exp x)