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 --   4.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   | (i >= 0) && (j >= 0) && (k >= 0) =
 204       idx v3d i j k
 205
 206   -- The next three are from the first line in (7.3).
 207   | (i == -1) && (j >= 0) && (k >= 0) =
 208       2*(value_at v3d 0 j k) - (value_at v3d 1 j k)
 209
 210   | (i >= 0) && (j == -1) && (k >= 0) =
 211       2*(value_at v3d i 0 k) - (value_at v3d i 1 k)
 212
 213   | (i >= 0) && (j >= 0) && (k == -1) =
 214       2*(value_at v3d i j 0) - (value_at v3d i j 1)
 215
 216   -- The next two are from the second line in (7.3).
 217   | (i == -1) && (j == -1) && (k >= 0) =
 218       2*(value_at v3d i 0 k) - (value_at v3d i 1 k)
 219
 220   | (i == -1) && (j == ysize) && (k >= 0) =
 221       2*(value_at v3d i (ysize - 1) k) - (value_at v3d i (ysize - 2) k)
 222
 223   -- The next two are from the third line in (7.3).
 224   | (i == -1) && (j >= 0) && (k == -1) =
 225       2*(value_at v3d i j 0) - (value_at v3d i j 1)
 226
 227   | (i == -1) && (j >= 0) && (k == zsize) =
 228       2*(value_at v3d i j (zsize - 1)) - (value_at v3d i j (zsize - 2))
 229
 230   -- Repeat the above (j and k) cases for i >= 0.
 231   | (i >= 0) && (j == -1) && (k == -1) =
 232       2*(value_at v3d i j 0) - (value_at v3d i j 1)
 233
 234   | (i == xsize) && (j == -1) && (k >= 0) =
 235       2*(value_at v3d (xsize - 1) j k) - (value_at v3d (xsize - 2) j k)
 236
 237   -- These two cases I made up.
 238   | (i == -1) && (j == -1) && (k == -1) =
 239       2*(value_at v3d i j 0) - (value_at v3d i j 1)
 240
 241   | (i == xsize) && (j == ysize) && (k == zsize) =
 242       2*(value_at v3d i j (zsize - 1)) - (value_at v3d i j (zsize - 2))
 243
 244   | otherwise =
 245       let istr = show i
 246           jstr = show j
 247           kstr = show k
 248           coordstr = "(" ++ istr ++ "," ++ jstr ++ "," ++ kstr ++ ")"
 249       in
 250         error $ "value_at called outside of domain: " ++ coordstr
 251   where
 252     (xsize, ysize, zsize) = dims v3d
 253
 254
 255 -- | Given a three-dimensional list of 'Double' and a set of 3D
 256 --   coordinates (i,j,k), constructs and returns the 'FunctionValues'
 257 --   object centered at (i,j,k)
 258 make_values :: Values3D -> Int -> Int -> Int -> FunctionValues
 259 make_values values i j k =
 260     empty_values { front  = value_at values (i-1) j k,
 261                    back   = value_at values (i+1) j k,
 262                    left   = value_at values i (j-1) k,
 263                    right  = value_at values i (j+1) k,
 264                    down   = value_at values i j (k-1),
 265                    top    = value_at values i j (k+1),
 266                    front_left = value_at values (i-1) (j-1) k,
 267                    front_right = value_at values (i-1) (j+1) k,
 268                    front_down =value_at values (i-1) j (k-1),
 269                    front_top = value_at values (i-1) j (k+1),
 270                    back_left = value_at values (i+1) (j-1) k,
 271                    back_right = value_at values (i+1) (j+1) k,
 272                    back_down = value_at values (i+1) j (k-1),
 273                    back_top = value_at values (i+1) j (k+1),
 274                    left_down = value_at values i (j-1) (k-1),
 275                    left_top = value_at values i (j-1) (k+1),
 276                    right_down = value_at values i (j+1) (k-1),
 277                    right_top = value_at values i (j+1) (k+1),
 278                    front_left_down = value_at values (i-1) (j-1) (k-1),
 279                    front_left_top = value_at values (i-1) (j-1) (k+1),
 280                    front_right_down = value_at values (i-1) (j+1) (k-1),
 281                    front_right_top = value_at values (i-1) (j+1) (k+1),
 282                    back_left_down = value_at values (i+1) (j-1) (k-1),
 283                    back_left_top = value_at values (i+1) (j-1) (k+1),
 284                    back_right_down = value_at values (i+1) (j+1) (k-1),
 285                    back_right_top = value_at values (i+1) (j+1) (k+1),
 286                    interior = value_at values i j k }
 287
 288 -- | Takes a 'FunctionValues' and a function that transforms one
 289 --   'Cardinal' to another (called a rotation). Then it applies the
 290 --   rotation to each element of the 'FunctionValues' object, and
 291 --   returns the result.
 292 rotate :: (Cardinal -> Cardinal) -> FunctionValues -> FunctionValues
 293 rotate rotation fv =
 294     FunctionValues { front  = eval fv (rotation F),
 295                      back   = eval fv (rotation B),
 296                      left   = eval fv (rotation L),
 297                      right  = eval fv (rotation R),
 298                      down   = eval fv (rotation D),
 299                      top    = eval fv (rotation T),
 300                      front_left = eval fv (rotation FL),
 301                      front_right = eval fv (rotation FR),
 302                      front_down = eval fv (rotation FD),
 303                      front_top = eval fv (rotation FT),
 304                      back_left = eval fv (rotation BL),
 305                      back_right = eval fv (rotation BR),
 306                      back_down = eval fv (rotation BD),
 307                      back_top = eval fv (rotation BT),
 308                      left_down = eval fv (rotation LD),
 309                      left_top = eval fv (rotation LT),
 310                      right_down = eval fv (rotation RD),
 311                      right_top = eval fv (rotation RT),
 312                      front_left_down = eval fv (rotation FLD),
 313                      front_left_top = eval fv (rotation FLT),
 314                      front_right_down = eval fv (rotation FRD),
 315                      front_right_top = eval fv (rotation FRT),
 316                      back_left_down = eval fv (rotation BLD),
 317                      back_left_top = eval fv (rotation BLT),
 318                      back_right_down = eval fv (rotation BRD),
 319                      back_right_top = eval fv (rotation BRT),
 320                      interior = interior fv }
 321
 322
 323
 324 -- | Ensure that the trilinear values wind up where we think they
 325 --   should.
 326 test_directions :: Assertion
 327 test_directions =
 328     assertTrue "all direction functions work" (and equalities)
 329         where
 330           fvs = make_values trilinear 1 1 1
 331           equalities = [ interior fvs == 4,
 332                          front fvs == 1,
 333                          back  fvs == 7,
 334                          left  fvs == 2,
 335                          right fvs == 6,
 336                          down  fvs == 3,
 337                          top   fvs == 5,
 338                          front_left fvs == 1,
 339                          front_right fvs == 1,
 340                          front_down fvs == 1,
 341                          front_top fvs == 1,
 342                          back_left fvs == 3,
 343                          back_right fvs == 11,
 344                          back_down fvs == 5,
 345                          back_top fvs == 9,
 346                          left_down fvs == 2,
 347                          left_top fvs == 2,
 348                          right_down fvs == 4,
 349                          right_top fvs == 8,
 350                          front_left_down fvs == 1,
 351                          front_left_top fvs == 1,
 352                          front_right_down fvs == 1,
 353                          front_right_top fvs == 1,
 354                          back_left_down fvs == 3,
 355                          back_left_top fvs == 3,
 356                          back_right_down fvs == 7,
 357                          back_right_top fvs == 15]
 358
 359
 360 function_values_tests :: Test.Framework.Test
 361 function_values_tests =
 362     testGroup "FunctionValues Tests"
 363                   [ testCase "test directions" test_directions ]
 364
 365
 366 prop_x_rotation_doesnt_affect_front :: FunctionValues -> Bool
 367 prop_x_rotation_doesnt_affect_front fv0 =
 368     expr1 == expr2
 369     where
 370       fv1 = rotate cwx fv0
 371       expr1 = front fv0
 372       expr2 = front fv1
 373
 374 prop_x_rotation_doesnt_affect_back :: FunctionValues -> Bool
 375 prop_x_rotation_doesnt_affect_back fv0 =
 376     expr1 == expr2
 377     where
 378       fv1 = rotate cwx fv0
 379       expr1 = back fv0
 380       expr2 = back fv1
 381
 382
 383 prop_y_rotation_doesnt_affect_left :: FunctionValues -> Bool
 384 prop_y_rotation_doesnt_affect_left fv0 =
 385     expr1 == expr2
 386     where
 387       fv1 = rotate cwy fv0
 388       expr1 = left fv0
 389       expr2 = left fv1
 390
 391 prop_y_rotation_doesnt_affect_right :: FunctionValues -> Bool
 392 prop_y_rotation_doesnt_affect_right fv0 =
 393     expr1 == expr2
 394     where
 395       fv1 = rotate cwy fv0
 396       expr1 = right fv0
 397       expr2 = right fv1
 398
 399
 400 prop_z_rotation_doesnt_affect_down :: FunctionValues -> Bool
 401 prop_z_rotation_doesnt_affect_down fv0 =
 402     expr1 == expr2
 403     where
 404       fv1 = rotate cwz fv0
 405       expr1 = down fv0
 406       expr2 = down fv1
 407
 408
 409 prop_z_rotation_doesnt_affect_top :: FunctionValues -> Bool
 410 prop_z_rotation_doesnt_affect_top fv0 =
 411     expr1 == expr2
 412     where
 413       fv1 = rotate cwz fv0
 414       expr1 = top fv0
 415       expr2 = top fv1
 416
 417
 418 function_values_properties :: Test.Framework.Test
 419 function_values_properties =
 420   let tp = testProperty
 421   in
 422     testGroup "FunctionValues Properties" [
 423       tp "x rotation doesn't affect front" prop_x_rotation_doesnt_affect_front,
 424       tp "x rotation doesn't affect back" prop_x_rotation_doesnt_affect_back,
 425       tp "y rotation doesn't affect left" prop_y_rotation_doesnt_affect_left,
 426       tp "y rotation doesn't affect right" prop_y_rotation_doesnt_affect_right,
 427       tp "z rotation doesn't affect top" prop_z_rotation_doesnt_affect_top,
 428       tp "z rotation doesn't affect down" prop_z_rotation_doesnt_affect_down ]