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