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
  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 -- | The eval function is where the magic happens for the
 139 --   FunctionValues type. Given a 'Cardinal' direction and a
 140 --   'FunctionValues' object, eval will return the value of the
 141 --   function f in that 'Cardinal' direction. Note that 'Cardinal' can
 142 --   be a composite type; eval is what performs the \"arithmetic\" on
 143 --   'Cardinal' directions.
 144 eval :: FunctionValues -> Cardinal -> Double
 145 eval f F = front f
 146 eval f B = back f
 147 eval f L = left f
 148 eval f R = right f
 149 eval f T = top f
 150 eval f D = down f
 151 eval f FL = front_left f
 152 eval f FR = front_right f
 153 eval f FD = front_down f
 154 eval f FT = front_top f
 155 eval f BL = back_left f
 156 eval f BR = back_right f
 157 eval f BD = back_down f
 158 eval f BT = back_top f
 159 eval f LD = left_down f
 160 eval f LT = left_top f
 161 eval f RD = right_down f
 162 eval f RT = right_top f
 163 eval f FLD = front_left_down f
 164 eval f FLT = front_left_top f
 165 eval f FRD = front_right_down f
 166 eval f FRT = front_right_top f
 167 eval f BLD = back_left_down f
 168 eval f BLT = back_left_top f
 169 eval f BRD = back_right_down f
 170 eval f BRT = back_right_top f
 171 eval f I   = interior f
 172 eval _ (Scalar x) = x
 173 eval f (Sum x y) = (eval f x) + (eval f y)
 174 eval f (Difference x y) = (eval f x) - (eval f y)
 175 eval f (Product x y) = (eval f x) * (eval f y)
 176 eval f (Quotient x y) = (eval f x) / (eval f y)
 177
 178 -- | Takes a three-dimensional list of 'Double' and a set of 3D
 179 --   coordinates (i,j,k), and returns the value at (i,j,k) in the
 180 --   supplied list. If there is no such value, we calculate one
 181 --   according to Sorokina and Zeilfelder, remark 7.3, p. 99.
 182 --
 183 --   We specifically do not consider values more than one unit away
 184 --   from our grid.
 185 --
 186 --   Examples:
 187 --
 188 --   >>> value_at Examples.trilinear 0 0 0
 189 --   1.0
 190 --
 191 --   >>> value_at Examples.trilinear (-1) 0 0
 192 --   0.0
 193 --
 194 --   >>> value_at Examples.trilinear 0 0 4
 195 --   1.0
 196 --
 197 --   >>> value_at Examples.trilinear 1 3 0
 198 --   5.0
 199 --
 200 value_at :: Values3D -> Int -> Int -> Int -> Double
 201 value_at v3d i j k
 202   -- Put the most common case first!
 203   | (valid_i i) && (valid_j j) && (valid_k k) =
 204       idx v3d i j k
 205
 206   -- The next three are from the first line in (7.3). Analogous cases
 207   -- have been added where the indices are one-too-big. These are the
 208   -- "one index is bad" cases.
 209   | not (valid_i i) =
 210       if (i == -1)
 211       then
 212         2*(value_at v3d 0 j k) - (value_at v3d 1 j k)
 213       else
 214         2*(value_at v3d (i-1) j k) - (value_at v3d (i-2) j k)
 215
 216   | not (valid_j j) =
 217       if (j == -1)
 218       then
 219         2*(value_at v3d i 0 k) - (value_at v3d i 1 k)
 220       else
 221         2*(value_at v3d i (j-1) k) - (value_at v3d i (j-2) k)
 222
 223   | not (valid_k k) =
 224       if (k == -1)
 225       then
 226         2*(value_at v3d i j 0) - (value_at v3d i j 1)
 227       else
 228         2*(value_at v3d i j (k-1)) - (value_at v3d i j (k-2))
 229
 230   | otherwise =
 231       let istr = show i
 232           jstr = show j
 233           kstr = show k
 234           coordstr = "(" ++ istr ++ "," ++ jstr ++ "," ++ kstr ++ ")"
 235       in
 236         error $ "value_at called outside of domain: " ++ coordstr
 237   where
 238     (xsize, ysize, zsize) = dims v3d
 239
 240     valid_i :: Int -> Bool
 241     valid_i i' = (i' >= 0) && (i' < xsize)
 242
 243     valid_j :: Int -> Bool
 244     valid_j j' = (j' >= 0) && (j' < ysize)
 245
 246     valid_k :: Int -> Bool
 247     valid_k k' = (k' >= 0) && (k' < zsize)
 248
 249
 250
 251 -- | Given a three-dimensional list of 'Double' and a set of 3D
 252 --   coordinates (i,j,k), constructs and returns the 'FunctionValues'
 253 --   object centered at (i,j,k)
 254 make_values :: Values3D -> Int -> Int -> Int -> FunctionValues
 255 make_values values i j k =
 256     empty_values { front  = value_at values (i-1) j k,
 257                    back   = value_at values (i+1) j k,
 258                    left   = value_at values i (j-1) k,
 259                    right  = value_at values i (j+1) k,
 260                    down   = value_at values i j (k-1),
 261                    top    = value_at values i j (k+1),
 262                    front_left = value_at values (i-1) (j-1) k,
 263                    front_right = value_at values (i-1) (j+1) k,
 264                    front_down =value_at values (i-1) j (k-1),
 265                    front_top = value_at values (i-1) j (k+1),
 266                    back_left = value_at values (i+1) (j-1) k,
 267                    back_right = value_at values (i+1) (j+1) k,
 268                    back_down = value_at values (i+1) j (k-1),
 269                    back_top = value_at values (i+1) j (k+1),
 270                    left_down = value_at values i (j-1) (k-1),
 271                    left_top = value_at values i (j-1) (k+1),
 272                    right_down = value_at values i (j+1) (k-1),
 273                    right_top = value_at values i (j+1) (k+1),
 274                    front_left_down = value_at values (i-1) (j-1) (k-1),
 275                    front_left_top = value_at values (i-1) (j-1) (k+1),
 276                    front_right_down = value_at values (i-1) (j+1) (k-1),
 277                    front_right_top = value_at values (i-1) (j+1) (k+1),
 278                    back_left_down = value_at values (i+1) (j-1) (k-1),
 279                    back_left_top = value_at values (i+1) (j-1) (k+1),
 280                    back_right_down = value_at values (i+1) (j+1) (k-1),
 281                    back_right_top = value_at values (i+1) (j+1) (k+1),
 282                    interior = value_at values i j k }
 283
 284 -- | Takes a 'FunctionValues' and a function that transforms one
 285 --   'Cardinal' to another (called a rotation). Then it applies the
 286 --   rotation to each element of the 'FunctionValues' object, and
 287 --   returns the result.
 288 rotate :: (Cardinal -> Cardinal) -> FunctionValues -> FunctionValues
 289 rotate rotation fv =
 290     FunctionValues { front  = eval fv (rotation F),
 291                      back   = eval fv (rotation B),
 292                      left   = eval fv (rotation L),
 293                      right  = eval fv (rotation R),
 294                      down   = eval fv (rotation D),
 295                      top    = eval fv (rotation T),
 296                      front_left = eval fv (rotation FL),
 297                      front_right = eval fv (rotation FR),
 298                      front_down = eval fv (rotation FD),
 299                      front_top = eval fv (rotation FT),
 300                      back_left = eval fv (rotation BL),
 301                      back_right = eval fv (rotation BR),
 302                      back_down = eval fv (rotation BD),
 303                      back_top = eval fv (rotation BT),
 304                      left_down = eval fv (rotation LD),
 305                      left_top = eval fv (rotation LT),
 306                      right_down = eval fv (rotation RD),
 307                      right_top = eval fv (rotation RT),
 308                      front_left_down = eval fv (rotation FLD),
 309                      front_left_top = eval fv (rotation FLT),
 310                      front_right_down = eval fv (rotation FRD),
 311                      front_right_top = eval fv (rotation FRT),
 312                      back_left_down = eval fv (rotation BLD),
 313                      back_left_top = eval fv (rotation BLT),
 314                      back_right_down = eval fv (rotation BRD),
 315                      back_right_top = eval fv (rotation BRT),
 316                      interior = interior fv }
 317
 318
 319
 320 -- | Ensure that the trilinear values wind up where we think they
 321 --   should.
 322 test_directions :: Assertion
 323 test_directions =
 324     assertTrue "all direction functions work" (and equalities)
 325         where
 326           fvs = make_values trilinear 1 1 1
 327           equalities = [ interior fvs == 4,
 328                          front fvs == 1,
 329                          back  fvs == 7,
 330                          left  fvs == 2,
 331                          right fvs == 6,
 332                          down  fvs == 3,
 333                          top   fvs == 5,
 334                          front_left fvs == 1,
 335                          front_right fvs == 1,
 336                          front_down fvs == 1,
 337                          front_top fvs == 1,
 338                          back_left fvs == 3,
 339                          back_right fvs == 11,
 340                          back_down fvs == 5,
 341                          back_top fvs == 9,
 342                          left_down fvs == 2,
 343                          left_top fvs == 2,
 344                          right_down fvs == 4,
 345                          right_top fvs == 8,
 346                          front_left_down fvs == 1,
 347                          front_left_top fvs == 1,
 348                          front_right_down fvs == 1,
 349                          front_right_top fvs == 1,
 350                          back_left_down fvs == 3,
 351                          back_left_top fvs == 3,
 352                          back_right_down fvs == 7,
 353                          back_right_top fvs == 15]
 354
 355
 356 function_values_tests :: Test.Framework.Test
 357 function_values_tests =
 358     testGroup "FunctionValues Tests"
 359                   [ testCase "test directions" test_directions ]
 360
 361
 362 prop_x_rotation_doesnt_affect_front :: FunctionValues -> Bool
 363 prop_x_rotation_doesnt_affect_front fv0 =
 364     expr1 == expr2
 365     where
 366       fv1 = rotate cwx fv0
 367       expr1 = front fv0
 368       expr2 = front fv1
 369
 370 prop_x_rotation_doesnt_affect_back :: FunctionValues -> Bool
 371 prop_x_rotation_doesnt_affect_back fv0 =
 372     expr1 == expr2
 373     where
 374       fv1 = rotate cwx fv0
 375       expr1 = back fv0
 376       expr2 = back fv1
 377
 378
 379 prop_y_rotation_doesnt_affect_left :: FunctionValues -> Bool
 380 prop_y_rotation_doesnt_affect_left fv0 =
 381     expr1 == expr2
 382     where
 383       fv1 = rotate cwy fv0
 384       expr1 = left fv0
 385       expr2 = left fv1
 386
 387 prop_y_rotation_doesnt_affect_right :: FunctionValues -> Bool
 388 prop_y_rotation_doesnt_affect_right fv0 =
 389     expr1 == expr2
 390     where
 391       fv1 = rotate cwy fv0
 392       expr1 = right fv0
 393       expr2 = right fv1
 394
 395
 396 prop_z_rotation_doesnt_affect_down :: FunctionValues -> Bool
 397 prop_z_rotation_doesnt_affect_down fv0 =
 398     expr1 == expr2
 399     where
 400       fv1 = rotate cwz fv0
 401       expr1 = down fv0
 402       expr2 = down fv1
 403
 404
 405 prop_z_rotation_doesnt_affect_top :: FunctionValues -> Bool
 406 prop_z_rotation_doesnt_affect_top fv0 =
 407     expr1 == expr2
 408     where
 409       fv1 = rotate cwz fv0
 410       expr1 = top fv0
 411       expr2 = top fv1
 412
 413
 414 function_values_properties :: Test.Framework.Test
 415 function_values_properties =
 416   let tp = testProperty
 417   in
 418     testGroup "FunctionValues Properties" [
 419       tp "x rotation doesn't affect front" prop_x_rotation_doesnt_affect_front,
 420       tp "x rotation doesn't affect back" prop_x_rotation_doesnt_affect_back,
 421       tp "y rotation doesn't affect left" prop_y_rotation_doesnt_affect_left,
 422       tp "y rotation doesn't affect right" prop_y_rotation_doesnt_affect_right,
 423       tp "z rotation doesn't affect top" prop_z_rotation_doesnt_affect_top,
 424       tp "z rotation doesn't affect down" prop_z_rotation_doesnt_affect_down ]