src/FunctionValues.hs

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