src/FunctionValues.hs

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