-- 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