-- SciErr.hs (C) 2004 Hardcore Software (www.hcsw.org) and Doug Hoyte
--
-- This program is distributed under the terms of the GNU GPL.
-- See www.gnu.org for more information.
--
-- Haskell Scientific Error Module
--
-- This module introduces a data type and implements instances of
-- various type classes for the data type. The data type, SciErr, holds
-- 2 pieces of information: A "value" and the "absolute error" of that
-- value.
--
-- The module encapsulates and automates the conversion to and from
-- absolute and relative errors while permitting you to use familiar
-- arithmetic operators.
--
-- For example, let's say you are trying to determine how fast the tips of
-- the blades on a windmill are spinning:
--
-- \-\ /-/
-- \ \ / /
-- \ \/ /
-- \ * \
-- /-/\ \
-- / / \ \
-- /-/ \-\
--
-- Excuse the horrible ASCII art. Assume that the lengths of all the blades
-- are equal.
--
-- Let's say that you measure the length of a blade from the centre of the
-- windmill to the tip of the blade to be 10 meters, plus or minus
-- 0.1 meters (10 centimeters), and you use your wristwatch to time how long
-- it takes for a blade to complete a revolution. Let's say it turns out
-- to be 5 seconds, plus or minus 1 second.
--
-- The first step is to figure out the circumference of the circle being
-- traversed by the tip of the windmill blades. Recall that the formula
-- for circumference is
--
-- 2 * PI * radius
--
-- Notice that 2 and PI have no scientific errors as they are known exactly.
--
-- Using this scientific error module, we evaluate the expression:
--
-- (2 +/- 0) * (pi +/- 0) * (10 +/- 0.1)
--
-- It evaluates to 62.8319 +/- 0.628319
--
-- So we know that the circumference of the circle is about 63 meters, give
-- or take a meter.
--
-- Notice that we just as easily could've evaluated
--
-- (10 +/- 0.1) * ((2 * pi) +/- 0)
--
-- The next step is to divide this distance by the time it took the blade
-- to traverse the circle. So, we evaluate:
--
-- (2 +/- 0) * (pi +/- 0) * (10 +/- 0.1) / (5 +/- 1)
--
-- This evaluates to: 12.5664 +/- 2.63894
--
-- So, we deduce that the windmill is rotating at roughly 13 meters per second,
-- plus or minus around 3 meters per second.
--
--
-- This simplistic example illustrated basic usage of this module.
--
-- Although you may create a SciErr datatype by using the data definition directly,
-- it is recommended that you use the supplied operator (+/-), as in the
-- example above.
--
-- Given a SciErr data type, you can extract the value and the absolute error
-- from the data type with the sciErrToTuple function. It returns
--
-- (value, error)
--
-- The familiar operations you may perform on SciErr datatypes are:
--
-- (+), (-), (*), (/), (^), (^^), and show
--
-- As SciErr derives the Eq class, the (==) operator will return True if, and
-- only if, the values and errors are identical. Since this is unlikely to happen
-- in serious usage, better solutions are the supplied (+/-==) and (+/-/=) operators.
--
-- Assuming (and this may be a fairly shakey assumption) that the actual values
-- fall equally within their specified error ranges, (+/-==) will tell you the
-- probability of the 2 SciErr datatypes being equal. (+/-/=) is, as you would
-- expect, the complement of (+/-==).
--
-- Note that ghc will warn that signum is undefined. This is deliberate, as the
-- sign of the number isn't necessarily known. (Consider 0 +/- 1)
module SciErr (SciErr, (+/-), sciErrToTuple, (+/-==), (+/-/=)) where
data SciErr = SciErr (Double, Double)
deriving (Read, Eq)
(+/-) :: Double -> Double -> SciErr
(+/-) s1 s2 | s2 < 0 = error "Absolute error can't be less than 0"
(+/-) s1 s2 = SciErr (s1, s2)
sciErrToTuple :: SciErr -> (Double, Double)
sciErrToTuple (SciErr (n, e)) = (n, e)
(+/-==), (+/-/=) :: SciErr -> SciErr -> Double
(+/-==) = shared
(+/-/=) s1 s2 = 1 - (s1 +/-== s2)
-- Internal functions for computing (in)equalities
overlap, spread, shared :: SciErr -> SciErr -> Double
shared s1 s2 = if toverlap == 0 then
0
else
toverlap / tspread
where toverlap = overlap s1 s2
tspread = spread s1 s2
overlap s1 s2 = if (val s2 - abserr s2 < val s1 - abserr s1) then
overlap s2 s1
else if (val s1 + abserr s1 < val s2 - abserr s2) then
0
else if (val s1 + abserr s1 < val s2 + abserr s1) then
(val s1 + abserr s1) - (val s2 - abserr s2)
else
(2 * abserr s2)
spread s1 s2 = (max (val s1 + abserr s1) (val s2 + abserr s2)) - (min (val s1 - abserr s1) (val s2 - abserr s2))
-- Internal accessor functions
val, abserr, relerr :: SciErr -> Double
val (SciErr (n, e)) = n
abserr (SciErr (n, e)) = e
relerr (SciErr (n, e)) = e / n
-- Various type class instances
instance Show SciErr where
show (SciErr (a,b)) = show a ++ " +/- " ++ show b
instance Num SciErr where
fromInteger n = SciErr (fromInteger n, 0)
abs s1 = SciErr (abs $ val s1, abserr s1)
(+) s1 s2 = SciErr (val s1 + val s2, abserr s1 + abserr s2)
(-) s1 s2 = SciErr (val s1 - val s2, abserr s1 + abserr s2)
(*) s1 s2 = SciErr (v, v * (relerr s1 + relerr s2))
where v = val s1 * val s2
instance Fractional SciErr where
fromRational r = SciErr (fromRational r, 0)
(/) s1 s2 = SciErr (v, v * (relerr s1 + relerr s2))
where v = val s1 / val s2