src/Grid.hs

   1 {-# LANGUAGE BangPatterns #-}
   2 -- | The Grid module contains the Grid type, its tests, and the 'zoom'
   3 --   function used to build the interpolation.
   4 module Grid (
   5   cube_at,
   6   grid_tests,
   7   slow_tests,
   8   zoom
   9   )
  10 where
  11
  12 import qualified Data.Array.Repa as R
  13 import qualified Data.Array.Repa.Operators.Traversal as R (unsafeTraverse)
  14 import Test.HUnit (Assertion, assertEqual)
  15 import Test.Framework (Test, testGroup)
  16 import Test.Framework.Providers.HUnit (testCase)
  17 import Test.Framework.Providers.QuickCheck2 (testProperty)
  18 import Test.QuickCheck ((==>),
  19                         Arbitrary(..),
  20                         Gen,
  21                         Property,
  22                         choose,
  23                         vectorOf)
  24 import Assertions (assertAlmostEqual, assertTrue)
  25 import Comparisons ((~=))
  26 import Cube (Cube(Cube),
  27              find_containing_tetrahedron,
  28              tetrahedra,
  29              tetrahedron)
  30 import Examples (trilinear, trilinear9x9x9, zeros)
  31 import FunctionValues (make_values, value_at)
  32 import Point (Point(..))
  33 import ScaleFactor (ScaleFactor)
  34 import Tetrahedron (
  35   Tetrahedron(v0,v1,v2,v3),
  36   c,
  37   polynomial,
  38   )
  39 import Values (Values3D, dims, empty3d, zoom_shape)
  40
  41
  42 -- | Our problem is defined on a Grid. The grid size is given by the
  43 --   positive number h, which we have defined to always be 1 for
  44 --   performance reasons (and simplicity). The function values are the
  45 --   values of the function at the grid points, which are distance h=1
  46 --   from one another in each direction (x,y,z).
  47 data Grid = Grid { function_values :: Values3D }
  48           deriving (Show)
  49
  50
  51 instance Arbitrary Grid where
  52     arbitrary = do
  53       x_dim <- choose (1, 27)
  54       y_dim <- choose (1, 27)
  55       z_dim <- choose (1, 27)
  56       elements <- vectorOf (x_dim * y_dim * z_dim) (arbitrary :: Gen Double)
  57       let new_shape = (R.Z R.:. x_dim R.:. y_dim R.:. z_dim)
  58       let fvs = R.fromListUnboxed new_shape elements
  59       return $ Grid fvs
  60
  61
  62
  63 -- | Takes a grid and a position as an argument and returns the cube
  64 --   centered on that position. If there is no cube there, well, you
  65 --   shouldn't have done that. The omitted "otherwise" case actually
  66 --   does improve performance.
  67 cube_at :: Grid -> Int -> Int -> Int -> Cube
  68 cube_at !g !i !j !k =
  69    Cube i j k fvs' tet_vol
  70    where
  71      fvs = function_values g
  72      fvs' = make_values fvs i j k
  73      tet_vol = 1/24
  74
  75
  76 --   The first cube along any axis covers (-1/2, 1/2). The second
  77 --   covers (1/2, 3/2).  The third, (3/2, 5/2), and so on.
  78 --
  79 --   We translate the (x,y,z) coordinates forward by 1/2 so that the
  80 --   first covers (0, 1), the second covers (1, 2), etc. This makes
  81 --   it easy to figure out which cube contains the given point.
  82 calculate_containing_cube_coordinate :: Grid -> Double -> Int
  83 calculate_containing_cube_coordinate g coord
  84     -- Don't use a cube on the boundary if we can help it. This
  85     -- returns cube #1 if we would have returned cube #0 and cube #1
  86     -- exists.
  87     | coord < offset = 0
  88     | coord == offset && (xsize > 1 && ysize > 1 && zsize > 1) = 1
  89     | otherwise = (ceiling (coord + offset)) - 1
  90     where
  91       (xsize, ysize, zsize) = dims (function_values g)
  92       offset = 1/2
  93
  94
  95 -- | Takes a 'Grid', and returns a 'Cube' containing the given 'Point'.
  96 --   Since our grid is rectangular, we can figure this out without having
  97 --   to check every cube.
  98 find_containing_cube :: Grid -> Point -> Cube
  99 find_containing_cube g (Point x y z) =
 100     cube_at g i j k
 101     where
 102       i = calculate_containing_cube_coordinate g x
 103       j = calculate_containing_cube_coordinate g y
 104       k = calculate_containing_cube_coordinate g z
 105
 106
 107 zoom_lookup :: Values3D -> ScaleFactor -> a -> (R.DIM3 -> Double)
 108 zoom_lookup v3d scale_factor _ =
 109     zoom_result v3d scale_factor
 110
 111
 112 zoom_result :: Values3D -> ScaleFactor -> R.DIM3 -> Double
 113 zoom_result v3d (sfx, sfy, sfz) (R.Z R.:. m R.:. n R.:. o) =
 114   f p
 115   where
 116     g = Grid v3d
 117     offset = 1/2
 118     m' = (fromIntegral m) / (fromIntegral sfx) - offset
 119     n' = (fromIntegral n) / (fromIntegral sfy) - offset
 120     o' = (fromIntegral o) / (fromIntegral sfz) - offset
 121     p  = Point m' n' o'
 122     cube = find_containing_cube g p
 123     t = find_containing_tetrahedron cube p
 124     f = polynomial t
 125
 126
 127 --
 128 -- Instead of IO, we could get away with a generic monad 'm'
 129 -- here. However, /we/ only call this function from within IO.
 130 --
 131 zoom :: Values3D -> ScaleFactor -> IO Values3D
 132 zoom v3d scale_factor
 133     | xsize == 0 || ysize == 0 || zsize == 0 = return empty3d
 134     | otherwise =
 135         R.computeUnboxedP $ R.unsafeTraverse v3d transExtent f
 136         where
 137           (xsize, ysize, zsize) = dims v3d
 138           transExtent = zoom_shape scale_factor
 139           f = zoom_lookup v3d scale_factor
 140
 141
 142 -- | Check all coefficients of tetrahedron0 belonging to the cube
 143 --   centered on (1,1,1) with a grid constructed from the trilinear
 144 --   values. See example one in the paper.
 145 --
 146 --   We also verify that the four vertices on face0 of the cube are
 147 --   in the correct location.
 148 --
 149 trilinear_c0_t0_tests :: Test.Framework.Test
 150 trilinear_c0_t0_tests =
 151   testGroup "trilinear c0 t0"
 152     [testGroup "coefficients"
 153       [testCase "c0030 is correct" test_trilinear_c0030,
 154        testCase "c0003 is correct" test_trilinear_c0003,
 155        testCase "c0021 is correct" test_trilinear_c0021,
 156        testCase "c0012 is correct" test_trilinear_c0012,
 157        testCase "c0120 is correct" test_trilinear_c0120,
 158        testCase "c0102 is correct" test_trilinear_c0102,
 159        testCase "c0111 is correct" test_trilinear_c0111,
 160        testCase "c0210 is correct" test_trilinear_c0210,
 161        testCase "c0201 is correct" test_trilinear_c0201,
 162        testCase "c0300 is correct" test_trilinear_c0300,
 163        testCase "c1020 is correct" test_trilinear_c1020,
 164        testCase "c1002 is correct" test_trilinear_c1002,
 165        testCase "c1011 is correct" test_trilinear_c1011,
 166        testCase "c1110 is correct" test_trilinear_c1110,
 167        testCase "c1101 is correct" test_trilinear_c1101,
 168        testCase "c1200 is correct" test_trilinear_c1200,
 169        testCase "c2010 is correct" test_trilinear_c2010,
 170        testCase "c2001 is correct" test_trilinear_c2001,
 171        testCase "c2100 is correct" test_trilinear_c2100,
 172        testCase "c3000 is correct" test_trilinear_c3000],
 173
 174     testGroup "face0 vertices"
 175       [testCase "v0 is correct" test_trilinear_f0_t0_v0,
 176        testCase "v1 is correct" test_trilinear_f0_t0_v1,
 177        testCase "v2 is correct" test_trilinear_f0_t0_v2,
 178        testCase "v3 is correct" test_trilinear_f0_t0_v3]
 179     ]
 180   where
 181     g = Grid trilinear
 182     cube = cube_at g 1 1 1
 183     t = tetrahedron cube 0
 184
 185     test_trilinear_c0030 :: Assertion
 186     test_trilinear_c0030 =
 187       assertAlmostEqual "c0030 is correct" (c t 0 0 3 0) (17/8)
 188
 189     test_trilinear_c0003 :: Assertion
 190     test_trilinear_c0003 =
 191       assertAlmostEqual "c0003 is correct" (c t 0 0 0 3) (27/8)
 192
 193     test_trilinear_c0021 :: Assertion
 194     test_trilinear_c0021 =
 195       assertAlmostEqual "c0021 is correct" (c t 0 0 2 1) (61/24)
 196
 197     test_trilinear_c0012 :: Assertion
 198     test_trilinear_c0012 =
 199       assertAlmostEqual "c0012 is correct" (c t 0 0 1 2) (71/24)
 200
 201     test_trilinear_c0120 :: Assertion
 202     test_trilinear_c0120 =
 203       assertAlmostEqual "c0120 is correct" (c t 0 1 2 0) (55/24)
 204
 205     test_trilinear_c0102 :: Assertion
 206     test_trilinear_c0102 =
 207       assertAlmostEqual "c0102 is correct" (c t 0 1 0 2) (73/24)
 208
 209     test_trilinear_c0111 :: Assertion
 210     test_trilinear_c0111 =
 211       assertAlmostEqual "c0111 is correct" (c t 0 1 1 1) (8/3)
 212
 213     test_trilinear_c0210 :: Assertion
 214     test_trilinear_c0210 =
 215       assertAlmostEqual "c0210 is correct" (c t 0 2 1 0) (29/12)
 216
 217     test_trilinear_c0201 :: Assertion
 218     test_trilinear_c0201 =
 219       assertAlmostEqual "c0201 is correct" (c t 0 2 0 1) (11/4)
 220
 221     test_trilinear_c0300 :: Assertion
 222     test_trilinear_c0300 =
 223       assertAlmostEqual "c0300 is correct" (c t 0 3 0 0) (5/2)
 224
 225     test_trilinear_c1020 :: Assertion
 226     test_trilinear_c1020 =
 227       assertAlmostEqual "c1020 is correct" (c t 1 0 2 0) (8/3)
 228
 229     test_trilinear_c1002 :: Assertion
 230     test_trilinear_c1002 =
 231       assertAlmostEqual "c1002 is correct" (c t 1 0 0 2) (23/6)
 232
 233     test_trilinear_c1011 :: Assertion
 234     test_trilinear_c1011 =
 235       assertAlmostEqual "c1011 is correct" (c t 1 0 1 1) (13/4)
 236
 237     test_trilinear_c1110 :: Assertion
 238     test_trilinear_c1110 =
 239       assertAlmostEqual "c1110 is correct" (c t 1 1 1 0) (23/8)
 240
 241     test_trilinear_c1101 :: Assertion
 242     test_trilinear_c1101 =
 243       assertAlmostEqual "c1101 is correct" (c t 1 1 0 1) (27/8)
 244
 245     test_trilinear_c1200 :: Assertion
 246     test_trilinear_c1200 =
 247       assertAlmostEqual "c1200 is correct" (c t 1 2 0 0) 3
 248
 249     test_trilinear_c2010 :: Assertion
 250     test_trilinear_c2010 =
 251       assertAlmostEqual "c2010 is correct" (c t 2 0 1 0) (10/3)
 252
 253     test_trilinear_c2001 :: Assertion
 254     test_trilinear_c2001 =
 255       assertAlmostEqual "c2001 is correct" (c t 2 0 0 1) 4
 256
 257     test_trilinear_c2100 :: Assertion
 258     test_trilinear_c2100 =
 259       assertAlmostEqual "c2100 is correct" (c t 2 1 0 0) (7/2)
 260
 261     test_trilinear_c3000 :: Assertion
 262     test_trilinear_c3000 =
 263       assertAlmostEqual "c3000 is correct" (c t 3 0 0 0) 4
 264
 265     test_trilinear_f0_t0_v0 :: Assertion
 266     test_trilinear_f0_t0_v0 =
 267       assertEqual "v0 is correct" (v0 t) (Point 1 1 1)
 268
 269     test_trilinear_f0_t0_v1 :: Assertion
 270     test_trilinear_f0_t0_v1 =
 271       assertEqual "v1 is correct" (v1 t) (Point 0.5 1 1)
 272
 273     test_trilinear_f0_t0_v2 :: Assertion
 274     test_trilinear_f0_t0_v2 =
 275       assertEqual "v2 is correct" (v2 t) (Point 0.5 0.5 1.5)
 276
 277     test_trilinear_f0_t0_v3 :: Assertion
 278     test_trilinear_f0_t0_v3 =
 279       assertEqual "v3 is correct" (v3 t) (Point 0.5 1.5 1.5)
 280
 281
 282 test_trilinear_reproduced :: Assertion
 283 test_trilinear_reproduced =
 284     assertTrue "trilinears are reproduced correctly" $
 285              and [p (Point i' j' k') ~= value_at trilinear i j k
 286                     | i <- [0..2],
 287                       j <- [0..2],
 288                       k <- [0..2],
 289                       c0 <- cs,
 290                       t <- tetrahedra c0,
 291                       let p = polynomial t,
 292                       let i' = fromIntegral i,
 293                       let j' = fromIntegral j,
 294                       let k' = fromIntegral k]
 295     where
 296       g = Grid trilinear
 297       cs = [ cube_at g ci cj ck | ci <- [0..2], cj <- [0..2], ck <- [0..2] ]
 298
 299
 300 test_zeros_reproduced :: Assertion
 301 test_zeros_reproduced =
 302     assertTrue "the zero function is reproduced correctly" $
 303              and [p (Point i' j' k') ~= value_at zeros i j k
 304                     | i <- [0..2],
 305                       j <- [0..2],
 306                       k <- [0..2],
 307                       let i' = fromIntegral i,
 308                       let j' = fromIntegral j,
 309                       let k' = fromIntegral k,
 310                       c0 <- cs,
 311                       t0 <- tetrahedra c0,
 312                       let p = polynomial t0 ]
 313     where
 314       g = Grid zeros
 315       cs = [ cube_at g ci cj ck | ci <- [0..2], cj <- [0..2], ck <- [0..2] ]
 316
 317
 318 -- | Make sure we can reproduce a 9x9x9 trilinear from the 3x3x3 one.
 319 test_trilinear9x9x9_reproduced :: Assertion
 320 test_trilinear9x9x9_reproduced =
 321     assertTrue "trilinear 9x9x9 is reproduced correctly" $
 322       and [p (Point i' j' k') ~= value_at trilinear9x9x9 i j k
 323             | i <- [0..8],
 324               j <- [0..8],
 325               k <- [0..8],
 326               t <- tetrahedra c0,
 327               let p = polynomial t,
 328               let i' = (fromIntegral i) * 0.5,
 329               let j' = (fromIntegral j) * 0.5,
 330               let k' = (fromIntegral k) * 0.5]
 331     where
 332       g = Grid trilinear
 333       c0 = cube_at g 1 1 1
 334
 335
 336
 337 prop_cube_indices_never_go_out_of_bounds :: Grid -> Gen Bool
 338 prop_cube_indices_never_go_out_of_bounds g =
 339   do
 340     let coordmin = negate (1/2)
 341
 342     let (xsize, ysize, zsize) = dims $ function_values g
 343     let xmax = (fromIntegral xsize) - (1/2)
 344     let ymax = (fromIntegral ysize) - (1/2)
 345     let zmax = (fromIntegral zsize) - (1/2)
 346
 347     x <- choose (coordmin, xmax)
 348     y <- choose (coordmin, ymax)
 349     z <- choose (coordmin, zmax)
 350
 351     let idx_x = calculate_containing_cube_coordinate g x
 352     let idx_y = calculate_containing_cube_coordinate g y
 353     let idx_z = calculate_containing_cube_coordinate g z
 354
 355     return $
 356       idx_x >= 0 &&
 357       idx_x <= xsize - 1 &&
 358       idx_y >= 0 &&
 359       idx_y <= ysize - 1 &&
 360       idx_z >= 0 &&
 361       idx_z <= zsize - 1
 362
 363
 364 -- | Given in Sorokina and Zeilfelder, p. 80, (2.9). Note that the
 365 --   third and fourth indices of c-t10 have been switched. This is
 366 --   because we store the triangles oriented such that their volume is
 367 --   positive. If T and T-tilde share \<v1,v2,v3\> and v0,v0-tilde point
 368 --   in opposite directions, one of them has to have negative volume!
 369 prop_c0120_identity :: Grid -> Property
 370 prop_c0120_identity g =
 371   xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
 372     c t0 0 1 2 0 ~= (c t0 1 0 2 0 + c t10 1 0 0 2) / 2
 373   where
 374     fvs = function_values g
 375     (xsize, ysize, zsize) = dims fvs
 376     cube0 = cube_at g 1 1 1
 377     cube1 = cube_at g 0 1 1
 378     t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
 379     t10 = tetrahedron cube1 10
 380
 381
 382 -- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
 383 --   'prop_c0120_identity'.
 384 prop_c0111_identity :: Grid -> Property
 385 prop_c0111_identity g =
 386   xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
 387     c t0 0 1 1 1 ~= (c t0 1 0 1 1 + c t10 1 0 1 1) / 2
 388   where
 389     fvs = function_values g
 390     (xsize, ysize, zsize) = dims fvs
 391     cube0 = cube_at g 1 1 1
 392     cube1 = cube_at g 0 1 1
 393     t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
 394     t10 = tetrahedron cube1 10
 395
 396
 397 -- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
 398 --   'prop_c0120_identity'.
 399 prop_c0201_identity :: Grid -> Property
 400 prop_c0201_identity g =
 401   xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
 402     c t0 0 2 0 1 ~= (c t0 1 1 0 1 + c t10 1 1 1 0) / 2
 403   where
 404     fvs = function_values g
 405     (xsize, ysize, zsize) = dims fvs
 406     cube0 = cube_at g 1 1 1
 407     cube1 = cube_at g 0 1 1
 408     t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
 409     t10 = tetrahedron cube1 10
 410
 411
 412 -- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
 413 --   'prop_c0120_identity'.
 414 prop_c0102_identity :: Grid -> Property
 415 prop_c0102_identity g =
 416   xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
 417     c t0 0 1 0 2 ~= (c t0 1 0 0 2 + c t10 1 0 2 0) / 2
 418   where
 419     fvs = function_values g
 420     (xsize, ysize, zsize) = dims fvs
 421     cube0 = cube_at g 1 1 1
 422     cube1 = cube_at g 0 1 1
 423     t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
 424     t10 = tetrahedron cube1 10
 425
 426
 427 -- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
 428 --   'prop_c0120_identity'.
 429 prop_c0210_identity :: Grid -> Property
 430 prop_c0210_identity g =
 431   xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
 432     c t0 0 2 1 0 ~= (c t0 1 1 1 0 + c t10 1 1 0 1) / 2
 433   where
 434     fvs = function_values g
 435     (xsize, ysize, zsize) = dims fvs
 436     cube0 = cube_at g 1 1 1
 437     cube1 = cube_at g 0 1 1
 438     t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
 439     t10 = tetrahedron cube1 10
 440
 441
 442 -- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
 443 --   'prop_c0120_identity'.
 444 prop_c0300_identity :: Grid -> Property
 445 prop_c0300_identity g =
 446   xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
 447     c t0 0 3 0 0 ~= (c t0 1 2 0 0 + c t10 1 2 0 0) / 2
 448   where
 449     fvs = function_values g
 450     (xsize, ysize, zsize) = dims fvs
 451     cube0 = cube_at g 1 1 1
 452     cube1 = cube_at g 0 1 1
 453     t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
 454     t10 = tetrahedron cube1 10
 455
 456
 457 -- | All of the properties from Section (2.9), p. 80. These require a
 458 --   grid since they refer to two adjacent cubes.
 459 p80_29_properties :: Test.Framework.Test
 460 p80_29_properties =
 461   testGroup "p. 80, Section (2.9) Properties" [
 462     testProperty "c0120 identity" prop_c0120_identity,
 463     testProperty "c0111 identity" prop_c0111_identity,
 464     testProperty "c0201 identity" prop_c0201_identity,
 465     testProperty "c0102 identity" prop_c0102_identity,
 466     testProperty "c0210 identity" prop_c0210_identity,
 467     testProperty "c0300 identity" prop_c0300_identity ]
 468
 469
 470 grid_tests :: Test.Framework.Test
 471 grid_tests =
 472     testGroup "Grid Tests" [
 473       trilinear_c0_t0_tests,
 474       p80_29_properties,
 475       testProperty "cube indices within bounds"
 476         prop_cube_indices_never_go_out_of_bounds ]
 477
 478
 479 -- Do the slow tests last so we can stop paying attention.
 480 slow_tests :: Test.Framework.Test
 481 slow_tests =
 482     testGroup "Slow Tests" [
 483       testCase "trilinear reproduced" test_trilinear_reproduced,
 484       testCase "trilinear9x9x9 reproduced" test_trilinear9x9x9_reproduced,
 485       testCase "zeros reproduced" test_zeros_reproduced ]