src/Piecewise.hs

   1 {-# LANGUAGE ScopedTypeVariables #-}
   2
   3 module Piecewise (
   4   Piecewise(..),
   5   evaluate,
   6   evaluate',
   7   from_intervals,
   8   list_plot' )
   9 where
  10
  11 import qualified Algebra.Additive as Additive ( C )
  12 import qualified Algebra.Field as Field ( C )
  13 import Control.Arrow ( first )
  14 import Misc ( partition )
  15 import NumericPrelude
  16 import qualified Prelude as P
  17
  18 -- | A predicate is basically a function that returns True or
  19 --   False. In this case, the predicate is used to determine which
  20 --   piece of a 'Piecewise' function we want to use.
  21 type Predicate a = (a -> Bool)
  22
  23
  24 -- | One piece of a piecewise function. The predicate determines
  25 --   whether or not the point lies within this piece, and the second
  26 --   component (a function) is what we'll evaluate if it does.
  27 type Piece a = (Predicate a, (a -> a))
  28
  29 -- | A representation of piecewise functions of one variable. A first
  30 --   approach might use a list of intervals paired with the function's
  31 --   value on those intervals. Rather than limit ourselves to a list
  32 --   of intervals, we allow the use of any predicate. The simple
  33 --   intervals case can be implemented with the predicate, \"is x
  34 --   between x1 and x2?\".
  35 data Piecewise a =
  36   Piecewise { pieces :: [Piece a] }
  37
  38
  39 -- | Evaluate a piecewise function at a point. If the point is within
  40 --   the domain of the function, the function is evaluated at the
  41 --   point and the result is returned wrapped in a 'Just'. If the
  42 --   point is outside of the domain, 'Nothing' is returned.
  43 evaluate :: Piecewise a -> a -> Maybe a
  44 evaluate (Piecewise ps) x =
  45   foldl f Nothing ps
  46   where
  47     f (Just y) _ = Just y
  48     f Nothing (p,g) = if p x then Just (g x) else Nothing
  49
  50
  51 -- | Evaluate a piecewise function at a point. If the point is within
  52 --   the domain of the function, the function is evaluated at the
  53 --   point and the result is returned. If the point is outside of the
  54 --   domain, zero is returned.
  55 evaluate' :: (Additive.C a) => Piecewise a -> a -> a
  56 evaluate' pw x =
  57   case evaluate pw x of
  58     Nothing -> zero
  59     Just result -> result
  60
  61
  62 -- | Construct a piecewise function from a list of intervals paired
  63 --   with functions. This is a convenience function that automatically
  64 --   converts intervals to \"contained in\" predicates.
  65 --
  66 --   Examples:
  67 --
  68 --   >>> let p1 = ((-1,0), \x -> -x) :: ((Double,Double), Double->Double)
  69 --   >>> let p2 = ((0,1), \x -> x) :: ((Double,Double), Double->Double)
  70 --   >>> let ivs = [p1, p2]
  71 --   >>> let pw = from_intervals ivs
  72 --   >>> evaluate' pw (-0.5)
  73 --   0.5
  74 --   >>> evaluate' pw 0.5
  75 --   0.5
  76 --
  77 from_intervals :: forall a. (Ord a) => [((a,a), (a -> a))] -> Piecewise a
  78 from_intervals =
  79   Piecewise . map (first f)
  80   where
  81     f :: (a,a) -> Predicate a
  82     f (x1,x2) x = x1 <= x && x <= x2
  83
  84
  85 -- | Plot a piecewise function over an interval @(x0,x1)@ using
  86 --   'evaluate''. Return the result a list of @(x,y)@ tuples.
  87 --
  88 --   Examples:
  89 --
  90 --   >>> let p1 = ((-1,0), \x -> -x) :: ((Double,Double), Double->Double)
  91 --   >>> let p2 = ((0,1), \x -> x) :: ((Double,Double), Double->Double)
  92 --   >>> let ivs = [p1, p2]
  93 --   >>> let pw = from_intervals ivs
  94 --   >>> list_plot' pw (-2) 2 4
  95 --   [(-2.0,0.0),(-1.0,1.0),(0.0,-0.0),(1.0,1.0),(2.0,0.0)]
  96 --
  97 list_plot' :: (Additive.C a, Field.C a)
  98            => Piecewise a
  99            -> a
 100            -> a
 101            -> Int
 102            -> [(a,a)]
 103 list_plot' pw x0 x1 subintervals =
 104   [ (x, evaluate' pw x) | x <- nodes ]
 105   where
 106     nodes = x0 : (map snd $ partition subintervals x0 x1)