logfloat (empty) → 0.8.1
raw patch · 4 files changed
+438/−0 lines, 4 filesdep +basebuild-type:Customsetup-changed
Dependencies added: base
Files
- Data/Number/LogFloat.lhs +329/−0
- LICENSE +33/−0
- Setup.hs +46/−0
- logfloat.cabal +30/−0
+ Data/Number/LogFloat.lhs view
@@ -0,0 +1,329 @@+%% 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 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.01+|+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 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 non-negative numbers+for efficiency's sake, see the forthcoming "Data.Number.LogFloat.Signed"+for doing signed log-domain calculations.+----------------------------------------------------------------++> module Data.Number.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.+>+> 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 = 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.
+ LICENSE view
@@ -0,0 +1,33 @@+Copyright (c) 2007--2008, wren ng thornton.+ALL RIGHTS RESERVED.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of the copyright holders nor the names of+ other contributors may be used to endorse or promote products+ derived from this software without specific prior written+ permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS+FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE+COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,+BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE+POSSIBILITY OF SUCH DAMAGE.+
+ Setup.hs view
@@ -0,0 +1,46 @@+#!/usr/bin/env runhaskell++module Main (main) where++-- <http://www.haskell.org/ghc/docs/latest/html/libraries/Cabal/Distribution-Simple.html>+import Distribution.Simple+import Distribution.Simple.Setup (CleanFlags, HaddockFlags)+import Distribution.Simple.LocalBuildInfo (LocalBuildInfo)+import Distribution.PackageDescription (HookedBuildInfo+ , emptyHookedBuildInfo+ , PackageDescription+ )+import System.Cmd (system)++main :: IO ()+main = defaultMainWithHooks defaultUserHooks+ { preHaddock = preHaddockScript+ , postClean = postCleanScript+ }+++preHaddockScript :: Args -> HaddockFlags -> IO HookedBuildInfo+preHaddockScript _ _ = do + putStrLn "Building lhs2hs..."+ system "ghc --make lhs2hs.hs -o lhs2hs"+ putStrLn "Illiterating Data.Number.LogFloat for Haddock..."+ system "cat Data/Number/LogFloat.lhs | ./lhs2hs > Data/Number/LogFloat.hs"+ return emptyHookedBuildInfo+++postCleanScript :: Args+ -> CleanFlags+ -> PackageDescription+ -> Maybe LocalBuildInfo+ -> IO ()+postCleanScript _ _ _ _ = do + putStrLn $ "removing files: " ++ commafy files+ removeAll files+ where+ files = ["Data/Number/LogFloat.hs", "lhs2hs", "lhs2hs.hi", "lhs2hs.o"]+ + removeAll = sequence_ . map (system . (++) "rm -f ")+ + commafy [] = ""+ commafy [x] = x+ commafy (x:xs@(_:_)) = x++", "++commafy xs
+ logfloat.cabal view
@@ -0,0 +1,30 @@+----------------------------------------------------------------+-- wren ng thornton <wren@cpan.org> ~ 2008.08.01+----------------------------------------------------------------++Name: logfloat+Version: 0.8.1+Cabal-Version: >= 1.2+Build-Type: Custom+Stability: stable+Copyright: Copyright (c) 2007--2008 wren ng thornton+License: BSD3+License-File: LICENSE+Author: wren ng thornton+Maintainer: wren@community.haskell.org+Homepage: http://code.haskell.org/~wren/+Category: Data, Math, Natural Language Processing+Synopsis: Log-domain floating point numbers+Description: 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+ probabilities as is done in Hidden Markov Models.+ It is also helpful for preventing overflow.++Library+ Exposed-Modules: Data.Number.LogFloat+ Build-Depends: base+ Ghc-Options: -O2 -fvia-C -optc-O3++----------------------------------------------------------------+----------------------------------------------------------- fin.