content: {-# LANGUAGE FlexibleContexts, ParallelListComp #-}
  {-# OPTIONS_GHC -fno-warn-unused-binds #-}
  
  module Validate (Similar(..), validate, validate') where
  
  import Data.Int
  import Data.Word
  import Data.Array.IArray
  import Foreign.C.Types
  import Foreign.Storable
  import Control.Exception                (assert)
  import Unsafe.Coerce
  
  class Similar a where
    sim :: a -> a -> Bool
  
  instance Similar Int     where sim = (==)
  instance Similar Int8    where sim = (==)
  instance Similar Int16   where sim = (==)
  instance Similar Int32   where sim = (==)
  instance Similar Int64   where sim = (==)
  instance Similar Word    where sim = (==)
  instance Similar Word8   where sim = (==)
  instance Similar Word16  where sim = (==)
  instance Similar Word32  where sim = (==)
  instance Similar Word64  where sim = (==)
  instance Similar CShort  where sim = (==)
  instance Similar CUShort where sim = (==)
  instance Similar CInt    where sim = (==)
  instance Similar CUInt   where sim = (==)
  instance Similar CLong   where sim = (==)
  instance Similar CULong  where sim = (==)
  instance Similar CLLong  where sim = (==)
  instance Similar CULLong where sim = (==)
  
  instance Similar Bool    where sim = (==)
  instance Similar Char    where sim = (==)
  instance Similar CChar   where sim = (==)
  instance Similar CSChar  where sim = (==)
  instance Similar CUChar  where sim = (==)
  
  instance Similar Float   where sim = absoluteOrRelative
  instance Similar CFloat  where sim = absoluteOrRelative
  instance Similar Double  where sim = absoluteOrRelative
  instance Similar CDouble where sim = absoluteOrRelative
  
  instance (Similar a, Similar b) => Similar (a,b) where
    (x,y) `sim` (u,v) = x `sim` u && y `sim` v
  
  --
  -- http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
  --
  
  absoluteOrRelative :: (Fractional a, Ord a) => a -> a -> Bool
  absoluteOrRelative u v
    | abs (u-v) < epsilonAbs = True
    | abs u > abs v          = abs ((u-v) / u) < epsilonRel
    | otherwise              = abs ((v-u) / v) < epsilonRel
    where
      epsilonRel = 0.001
      epsilonAbs = 0.00001
  
  
  -- Comparisons using lexicographically ordered floating-point numbers
  -- reinterpreted as twos-complement integers.
  --
  lexicographic32 :: (Num a, Storable a) => Int -> a -> a -> Bool
  lexicographic32 maxUlps a b
    = assert (sizeOf a == 4 && maxUlps > 0 && maxUlps < 4 * 1024 * 1024)
    $ intDiff < fromIntegral maxUlps
    where
      intDiff = abs (toInt a - toInt b)
      toInt x | x' < 0    = 0x80000000 - x'
              | otherwise = x'
              where x'    = unsafeCoerce x :: Int32
  
  
  lexicographic64 :: (Num a, Storable a) => Int -> a -> a -> Bool
  lexicographic64 maxUlps a b
    = assert (sizeOf a == 8 && maxUlps > 0 && maxUlps < 8 * 1024 * 1024)
    $ intDiff < fromIntegral maxUlps
    where
      intDiff = abs (toInt a - toInt b)
      toInt x | x' < 0    = 0x8000000000000000 - x'
              | otherwise = x'
              where x'    = unsafeCoerce x :: Int64
  
  
  -- Compare two vectors element-wise for equality, for a given measure of
  -- similarity. The index and values are returned for pairs that fail.
  --
  validate
    :: (IArray array e, Ix ix, Similar e)
    => array ix e
    -> array ix e
    -> [(ix,(e,e))]
  validate ref arr = validate' (assocs ref) (elems arr)
  
  validate' :: (Ix ix, Similar e) => [(ix,e)] -> [e] -> [(ix,(e,e))]
  validate' ref arr =
    filter (not . uncurry sim . snd) [ (i,(x,y)) | (i,x) <- ref | y <- arr ]
  
  