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