src/FunctionValues.hs

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