src/FunctionValues.hs

   1 {-# LANGUAGE BangPatterns #-}
   2 -- | The FunctionValues module contains the 'FunctionValues' type and
   3 --   the functions used to manipulate it.
   4 module FunctionValues (
   5   FunctionValues(..),
   6   empty_values,
   7   eval,
   8   make_values,
   9   rotate,
  10   function_values_tests,
  11   function_values_properties,
  12   value_at
  13   )
  14 where
  15
  16 import Prelude hiding (LT)
  17 import Test.HUnit (Assertion)
  18 import Test.Framework (Test, testGroup)
  19 import Test.Framework.Providers.HUnit (testCase)
  20 import Test.Framework.Providers.QuickCheck2 (testProperty)
  21 import Test.QuickCheck (Arbitrary(..), choose)
  22
  23 import Assertions (assertTrue)
  24 import Cardinal ( Cardinal(..), cwx, cwy, cwz )
  25 import Examples (trilinear)
  26 import Values (Values3D, dims, idx)
  27
  28 -- | The FunctionValues type represents the value of our function f at
  29 --   the 27 points surrounding (and including) the center of a
  30 --   cube. Each value of f can be accessed by the name of its
  31 --   direction.
  32 data FunctionValues =
  33     FunctionValues { front              :: !Double,
  34                      back               :: !Double,
  35                      left               :: !Double,
  36                      right              :: !Double,
  37                      top                :: !Double,
  38                      down               :: !Double,
  39                      front_left         :: !Double,
  40                      front_right        :: !Double,
  41                      front_down         :: !Double,
  42                      front_top          :: !Double,
  43                      back_left          :: !Double,
  44                      back_right         :: !Double,
  45                      back_down          :: !Double,
  46                      back_top           :: !Double,
  47                      left_down          :: !Double,
  48                      left_top           :: !Double,
  49                      right_down         :: !Double,
  50                      right_top          :: !Double,
  51                      front_left_down    :: !Double,
  52                      front_left_top     :: !Double,
  53                      front_right_down   :: !Double,
  54                      front_right_top    :: !Double,
  55                      back_left_down     :: !Double,
  56                      back_left_top      :: !Double,
  57                      back_right_down    :: !Double,
  58                      back_right_top     :: !Double,
  59                      interior           :: !Double }
  60       deriving (Eq, Show)
  61
  62
  63 instance Arbitrary FunctionValues where
  64     arbitrary = do
  65       front'  <- choose (min_double, max_double)
  66       back'   <- choose (min_double, max_double)
  67       left'   <- choose (min_double, max_double)
  68       right'  <- choose (min_double, max_double)
  69       top'    <- choose (min_double, max_double)
  70       down'   <- choose (min_double, max_double)
  71       front_left' <- choose (min_double, max_double)
  72       front_right' <- choose (min_double, max_double)
  73       front_top' <- choose (min_double, max_double)
  74       front_down' <- choose (min_double, max_double)
  75       back_left' <- choose (min_double, max_double)
  76       back_right' <- choose (min_double, max_double)
  77       back_top' <- choose (min_double, max_double)
  78       back_down' <- choose (min_double, max_double)
  79       left_top' <- choose (min_double, max_double)
  80       left_down' <- choose (min_double, max_double)
  81       right_top' <- choose (min_double, max_double)
  82       right_down' <- choose (min_double, max_double)
  83       front_left_top' <- choose (min_double, max_double)
  84       front_left_down' <- choose (min_double, max_double)
  85       front_right_top' <- choose (min_double, max_double)
  86       front_right_down' <- choose (min_double, max_double)
  87       back_left_top' <- choose (min_double, max_double)
  88       back_left_down' <- choose (min_double, max_double)
  89       back_right_top' <- choose (min_double, max_double)
  90       back_right_down' <- choose (min_double, max_double)
  91       interior' <- choose (min_double, max_double)
  92
  93       return empty_values { front = front',
  94                             back  = back',
  95                             left  = left',
  96                             right = right',
  97                             top   = top',
  98                             down  = down',
  99                             front_left = front_left',
 100                             front_right = front_right',
 101                             front_top = front_top',
 102                             front_down = front_down',
 103                             back_left = back_left',
 104                             back_right = back_right',
 105                             back_top = back_top',
 106                             back_down = back_down',
 107                             left_top = left_top',
 108                             left_down = left_down',
 109                             right_top = right_top',
 110                             right_down = right_down',
 111                             front_left_top = front_left_top',
 112                             front_left_down = front_left_down',
 113                             front_right_top = front_right_top',
 114                             front_right_down = front_right_down',
 115                             back_left_top = back_left_top',
 116                             back_left_down = back_left_down',
 117                             back_right_top = back_right_top',
 118                             back_right_down = back_right_down',
 119                             interior = interior' }
 120       where
 121         -- | We perform addition with the function values contained in a
 122         --   FunctionValues object. If we choose random doubles near the machine
 123         --   min/max, we risk overflowing or underflowing the 'Double'. This
 124         --   places a reasonably safe limit on the maximum size of our generated
 125         --   'Double' members.
 126         max_double :: Double
 127         max_double = 10000.0
 128
 129         -- | See 'max_double'.
 130         min_double :: Double
 131         min_double = (-1) * max_double
 132
 133
 134 -- | Return a 'FunctionValues' with a bunch of zeros for data points.
 135 empty_values :: FunctionValues
 136 empty_values =
 137     FunctionValues 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 138
 139
 140 -- | The eval function is where the magic happens for the
 141 --   FunctionValues type. Given a 'Cardinal' direction and a
 142 --   'FunctionValues' object, eval will return the value of the
 143 --   function f in that 'Cardinal' direction. Note that 'Cardinal' can
 144 --   be a composite type; eval is what performs the \"arithmetic\" on
 145 --   'Cardinal' directions.
 146 eval :: FunctionValues -> Cardinal -> Double
 147 eval f F = front f
 148 eval f B = back f
 149 eval f L = left f
 150 eval f R = right f
 151 eval f T = top f
 152 eval f D = down f
 153 eval f FL = front_left f
 154 eval f FR = front_right f
 155 eval f FD = front_down f
 156 eval f FT = front_top f
 157 eval f BL = back_left f
 158 eval f BR = back_right f
 159 eval f BD = back_down f
 160 eval f BT = back_top f
 161 eval f LD = left_down f
 162 eval f LT = left_top f
 163 eval f RD = right_down f
 164 eval f RT = right_top f
 165 eval f FLD = front_left_down f
 166 eval f FLT = front_left_top f
 167 eval f FRD = front_right_down f
 168 eval f FRT = front_right_top f
 169 eval f BLD = back_left_down f
 170 eval f BLT = back_left_top f
 171 eval f BRD = back_right_down f
 172 eval f BRT = back_right_top f
 173 eval f I   = interior f
 174 eval _ (Scalar x) = x
 175 eval f (Sum x y) = (eval f x) + (eval f y)
 176 eval f (Difference x y) = (eval f x) - (eval f y)
 177 eval f (Product x y) = (eval f x) * (eval f y)
 178 eval f (Quotient x y) = (eval f x) / (eval f y)
 179
 180
 181 -- | Takes a three-dimensional list of 'Double' and a set of 3D
 182 --   coordinates (i,j,k), and returns the value at (i,j,k) in the
 183 --   supplied list. If there is no such value, we calculate one
 184 --   according to Sorokina and Zeilfelder, remark 7.3, p. 99.
 185 --
 186 --   We specifically do not consider values more than one unit away
 187 --   from our grid.
 188 --
 189 --   Examples:
 190 --
 191 --   >>> value_at Examples.trilinear 0 0 0
 192 --   1.0
 193 --
 194 --   >>> value_at Examples.trilinear (-1) 0 0
 195 --   0.0
 196 --
 197 --   >>> value_at Examples.trilinear 0 0 4
 198 --   1.0
 199 --
 200 --   >>> value_at Examples.trilinear 1 3 0
 201 --   5.0
 202 --
 203 value_at :: Values3D -> Int -> Int -> Int -> Double
 204 value_at v3d !i !j !k
 205   -- Put the most common case first!
 206   | (valid_i i) && (valid_j j) && (valid_k k) =
 207       idx v3d i j k
 208
 209   -- The next three are from the first line in (7.3). Analogous cases
 210   -- have been added where the indices are one-too-big. These are the
 211   -- "one index is bad" cases.
 212   | not (valid_i i) =
 213       if (dim_i == 1)
 214       then
 215         -- We're one-dimensional in our first coordinate, so just
 216         -- return the data point that we do have. If we try to use
 217         -- the formula from remark 7.3, we go into an infinite loop.
 218         value_at v3d 0 j k
 219       else
 220         if (i == -1)
 221         then
 222           2*(value_at v3d 0 j k) - (value_at v3d 1 j k)
 223         else
 224           2*(value_at v3d (i-1) j k) - (value_at v3d (i-2) j k)
 225
 226   | not (valid_j j) =
 227       if (dim_j == 1)
 228       then
 229         -- We're one-dimensional in our second coordinate, so just
 230         -- return the data point that we do have. If we try to use
 231         -- the formula from remark 7.3, we go into an infinite loop.
 232         value_at v3d i 0 k
 233       else
 234         if (j == -1)
 235         then
 236           2*(value_at v3d i 0 k) - (value_at v3d i 1 k)
 237         else
 238           2*(value_at v3d i (j-1) k) - (value_at v3d i (j-2) k)
 239
 240   | not (valid_k k) =
 241       if (dim_k == 1)
 242       then
 243         -- We're one-dimensional in our third coordinate, so just
 244         -- return the data point that we do have. If we try to use
 245         -- the formula from remark 7.3, we go into an infinite loop.
 246         value_at v3d i j 0
 247       else
 248         if (k == -1)
 249         then
 250           2*(value_at v3d i j 0) - (value_at v3d i j 1)
 251         else
 252           2*(value_at v3d i j (k-1)) - (value_at v3d i j (k-2))
 253
 254   | otherwise =
 255       let istr = show i
 256           jstr = show j
 257           kstr = show k
 258           coordstr = "(" ++ istr ++ "," ++ jstr ++ "," ++ kstr ++ ")"
 259       in
 260         error $ "value_at called outside of domain: " ++ coordstr
 261   where
 262     (dim_i, dim_j, dim_k) = dims v3d
 263
 264     valid_i :: Int -> Bool
 265     valid_i i' = (i' >= 0) && (i' < dim_i)
 266
 267     valid_j :: Int -> Bool
 268     valid_j j' = (j' >= 0) && (j' < dim_j)
 269
 270     valid_k :: Int -> Bool
 271     valid_k k' = (k' >= 0) && (k' < dim_k)
 272
 273
 274
 275 -- | Given a three-dimensional list of 'Double' and a set of 3D
 276 --   coordinates (i,j,k), constructs and returns the 'FunctionValues'
 277 --   object centered at (i,j,k)
 278 make_values :: Values3D -> Int -> Int -> Int -> FunctionValues
 279 make_values values !i !j !k =
 280     empty_values { front  = value_at values (i-1) j k,
 281                    back   = value_at values (i+1) j k,
 282                    left   = value_at values i (j-1) k,
 283                    right  = value_at values i (j+1) k,
 284                    down   = value_at values i j (k-1),
 285                    top    = value_at values i j (k+1),
 286                    front_left = value_at values (i-1) (j-1) k,
 287                    front_right = value_at values (i-1) (j+1) k,
 288                    front_down =value_at values (i-1) j (k-1),
 289                    front_top = value_at values (i-1) j (k+1),
 290                    back_left = value_at values (i+1) (j-1) k,
 291                    back_right = value_at values (i+1) (j+1) k,
 292                    back_down = value_at values (i+1) j (k-1),
 293                    back_top = value_at values (i+1) j (k+1),
 294                    left_down = value_at values i (j-1) (k-1),
 295                    left_top = value_at values i (j-1) (k+1),
 296                    right_down = value_at values i (j+1) (k-1),
 297                    right_top = value_at values i (j+1) (k+1),
 298                    front_left_down = value_at values (i-1) (j-1) (k-1),
 299                    front_left_top = value_at values (i-1) (j-1) (k+1),
 300                    front_right_down = value_at values (i-1) (j+1) (k-1),
 301                    front_right_top = value_at values (i-1) (j+1) (k+1),
 302                    back_left_down = value_at values (i+1) (j-1) (k-1),
 303                    back_left_top = value_at values (i+1) (j-1) (k+1),
 304                    back_right_down = value_at values (i+1) (j+1) (k-1),
 305                    back_right_top = value_at values (i+1) (j+1) (k+1),
 306                    interior = value_at values i j k }
 307
 308 -- | Takes a 'FunctionValues' and a function that transforms one
 309 --   'Cardinal' to another (called a rotation). Then it applies the
 310 --   rotation to each element of the 'FunctionValues' object, and
 311 --   returns the result.
 312 rotate :: (Cardinal -> Cardinal) -> FunctionValues -> FunctionValues
 313 rotate rotation fv =
 314     FunctionValues { front  = eval fv (rotation F),
 315                      back   = eval fv (rotation B),
 316                      left   = eval fv (rotation L),
 317                      right  = eval fv (rotation R),
 318                      down   = eval fv (rotation D),
 319                      top    = eval fv (rotation T),
 320                      front_left = eval fv (rotation FL),
 321                      front_right = eval fv (rotation FR),
 322                      front_down = eval fv (rotation FD),
 323                      front_top = eval fv (rotation FT),
 324                      back_left = eval fv (rotation BL),
 325                      back_right = eval fv (rotation BR),
 326                      back_down = eval fv (rotation BD),
 327                      back_top = eval fv (rotation BT),
 328                      left_down = eval fv (rotation LD),
 329                      left_top = eval fv (rotation LT),
 330                      right_down = eval fv (rotation RD),
 331                      right_top = eval fv (rotation RT),
 332                      front_left_down = eval fv (rotation FLD),
 333                      front_left_top = eval fv (rotation FLT),
 334                      front_right_down = eval fv (rotation FRD),
 335                      front_right_top = eval fv (rotation FRT),
 336                      back_left_down = eval fv (rotation BLD),
 337                      back_left_top = eval fv (rotation BLT),
 338                      back_right_down = eval fv (rotation BRD),
 339                      back_right_top = eval fv (rotation BRT),
 340                      interior = interior fv }
 341
 342
 343
 344 -- | Ensure that the trilinear values wind up where we think they
 345 --   should.
 346 test_directions :: Assertion
 347 test_directions =
 348     assertTrue "all direction functions work" (and equalities)
 349         where
 350           fvs = make_values trilinear 1 1 1
 351           equalities = [ interior fvs == 4,
 352                          front fvs == 1,
 353                          back  fvs == 7,
 354                          left  fvs == 2,
 355                          right fvs == 6,
 356                          down  fvs == 3,
 357                          top   fvs == 5,
 358                          front_left fvs == 1,
 359                          front_right fvs == 1,
 360                          front_down fvs == 1,
 361                          front_top fvs == 1,
 362                          back_left fvs == 3,
 363                          back_right fvs == 11,
 364                          back_down fvs == 5,
 365                          back_top fvs == 9,
 366                          left_down fvs == 2,
 367                          left_top fvs == 2,
 368                          right_down fvs == 4,
 369                          right_top fvs == 8,
 370                          front_left_down fvs == 1,
 371                          front_left_top fvs == 1,
 372                          front_right_down fvs == 1,
 373                          front_right_top fvs == 1,
 374                          back_left_down fvs == 3,
 375                          back_left_top fvs == 3,
 376                          back_right_down fvs == 7,
 377                          back_right_top fvs == 15]
 378
 379
 380 function_values_tests :: Test.Framework.Test
 381 function_values_tests =
 382     testGroup "FunctionValues Tests"
 383                   [ testCase "test directions" test_directions ]
 384
 385
 386 prop_x_rotation_doesnt_affect_front :: FunctionValues -> Bool
 387 prop_x_rotation_doesnt_affect_front fv0 =
 388     expr1 == expr2
 389     where
 390       fv1 = rotate cwx fv0
 391       expr1 = front fv0
 392       expr2 = front fv1
 393
 394 prop_x_rotation_doesnt_affect_back :: FunctionValues -> Bool
 395 prop_x_rotation_doesnt_affect_back fv0 =
 396     expr1 == expr2
 397     where
 398       fv1 = rotate cwx fv0
 399       expr1 = back fv0
 400       expr2 = back fv1
 401
 402
 403 prop_y_rotation_doesnt_affect_left :: FunctionValues -> Bool
 404 prop_y_rotation_doesnt_affect_left fv0 =
 405     expr1 == expr2
 406     where
 407       fv1 = rotate cwy fv0
 408       expr1 = left fv0
 409       expr2 = left fv1
 410
 411 prop_y_rotation_doesnt_affect_right :: FunctionValues -> Bool
 412 prop_y_rotation_doesnt_affect_right fv0 =
 413     expr1 == expr2
 414     where
 415       fv1 = rotate cwy fv0
 416       expr1 = right fv0
 417       expr2 = right fv1
 418
 419
 420 prop_z_rotation_doesnt_affect_down :: FunctionValues -> Bool
 421 prop_z_rotation_doesnt_affect_down fv0 =
 422     expr1 == expr2
 423     where
 424       fv1 = rotate cwz fv0
 425       expr1 = down fv0
 426       expr2 = down fv1
 427
 428
 429 prop_z_rotation_doesnt_affect_top :: FunctionValues -> Bool
 430 prop_z_rotation_doesnt_affect_top fv0 =
 431     expr1 == expr2
 432     where
 433       fv1 = rotate cwz fv0
 434       expr1 = top fv0
 435       expr2 = top fv1
 436
 437
 438 function_values_properties :: Test.Framework.Test
 439 function_values_properties =
 440   let tp = testProperty
 441   in
 442     testGroup "FunctionValues Properties" [
 443       tp "x rotation doesn't affect front" prop_x_rotation_doesnt_affect_front,
 444       tp "x rotation doesn't affect back" prop_x_rotation_doesnt_affect_back,
 445       tp "y rotation doesn't affect left" prop_y_rotation_doesnt_affect_left,
 446       tp "y rotation doesn't affect right" prop_y_rotation_doesnt_affect_right,
 447       tp "z rotation doesn't affect top" prop_z_rotation_doesnt_affect_top,
 448       tp "z rotation doesn't affect down" prop_z_rotation_doesnt_affect_down ]