X-Git-Url: http://gitweb.michael.orlitzky.com/?p=spline3.git;a=blobdiff_plain;f=src%2FGrid.hs;h=269b37cec68d5c2f6211971dc1af32316b86c7ae;hp=893fe5c360d3f39c73a99ebab26d1320fd798815;hb=fd3d394c27e3a90de8238b98bd112c4fe3468ee0;hpb=89a8f8f3728c6e8730a37304708312b67bc23274

diff --git a/src/Grid.hs b/src/Grid.hs
index 893fe5c..269b37c 100644
--- a/src/Grid.hs
+++ b/src/Grid.hs
@@ -1,109 +1,83 @@
--- | The Grid module just contains the Grid type and two constructors
---   for it. We hide the main Grid constructor because we don't want
---   to allow instantiation of a grid with h <= 0.
+{-# LANGUAGE BangPatterns #-}
+-- | The Grid module contains the Grid type, its tests, and the 'zoom'
+--   function used to build the interpolation.
 module Grid (
   cube_at,
   grid_tests,
-  make_grid,
   slow_tests,
   zoom
   )
 where
 
-import Data.Array (Array, array, (!))
 import qualified Data.Array.Repa as R
-import Test.HUnit
+import qualified Data.Array.Repa.Operators.Traversal as R (unsafeTraverse)
+import Test.HUnit (Assertion, assertEqual)
 import Test.Framework (Test, testGroup)
 import Test.Framework.Providers.HUnit (testCase)
 import Test.Framework.Providers.QuickCheck2 (testProperty)
-import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
-
-import Assertions
-import Comparisons
+import Test.QuickCheck ((==>),
+                        Arbitrary(..),
+                        Gen,
+                        Property,
+                        choose,
+                        vectorOf)
+import Assertions (assertAlmostEqual, assertTrue)
+import Comparisons ((~=))
 import Cube (Cube(Cube),
              find_containing_tetrahedron,
              tetrahedra,
              tetrahedron)
-import Examples
-import FunctionValues
-import Point (Point)
-import PolynomialArray (PolynomialArray)
-import ScaleFactor
-import Tetrahedron (Tetrahedron, c, number, polynomial, v0, v1, v2, v3)
-import ThreeDimensional
+import Examples (trilinear, trilinear9x9x9, zeros)
+import FunctionValues (make_values, value_at)
+import Point (Point(..))
+import ScaleFactor (ScaleFactor)
+import Tetrahedron (
+  Tetrahedron(v0,v1,v2,v3),
+  c,
+  polynomial,
+  )
 import Values (Values3D, dims, empty3d, zoom_shape)
 
 
-type CubeGrid = Array (Int,Int,Int) Cube
-
-
 -- | Our problem is defined on a Grid. The grid size is given by the
---   positive number h. The function values are the values of the
---   function at the grid points, which are distance h from one
---   another in each direction (x,y,z).
-data Grid = Grid { h :: Double, -- MUST BE GREATER THAN ZERO!
-                   function_values :: Values3D,
-                   cube_grid :: CubeGrid }
-          deriving (Eq, Show)
+--   positive number h, which we have defined to always be 1 for
+--   performance reasons (and simplicity). The function values are the
+--   values of the function at the grid points, which are distance h=1
+--   from one another in each direction (x,y,z).
+data Grid = Grid { function_values :: Values3D }
+          deriving (Show)
 
 
 instance Arbitrary Grid where
     arbitrary = do
-      (Positive h') <- arbitrary :: Gen (Positive Double)
-      fvs <- arbitrary :: Gen Values3D
-      return (make_grid h' fvs)
-
-
--- | The constructor that we want people to use. If we're passed a
---   non-positive grid size, we throw an error.
-make_grid :: Double -> Values3D -> Grid
-make_grid grid_size values
-    | grid_size <= 0 = error "grid size must be positive"
-    | otherwise = Grid grid_size values (cubes grid_size values)
-
-
--- | Returns a three-dimensional array of cubes centered on the grid
---   points (h*i, h*j, h*k) with the appropriate 'FunctionValues'.
-cubes :: Double -> Values3D -> CubeGrid
-cubes delta fvs 
-  = array (lbounds, ubounds)
-           [ ((i,j,k), cube_ijk)
-                 | i <- [0..xmax],
-                   j <- [0..ymax],
-                   k <- [0..zmax],
-                   let tet_vol = (1/24)*(delta^(3::Int)),
-                   let cube_ijk =
-                         Cube delta i j k (make_values fvs i j k) tet_vol]
-     where
-       xmax = xsize - 1
-       ymax = ysize - 1
-       zmax = zsize - 1
-       lbounds = (0, 0, 0)
-       ubounds = (xmax, ymax, zmax)
-       (xsize, ysize, zsize) = dims fvs
+      x_dim <- choose (1, 27)
+      y_dim <- choose (1, 27)
+      z_dim <- choose (1, 27)
+      elements <- vectorOf (x_dim * y_dim * z_dim) (arbitrary :: Gen Double)
+      let new_shape = (R.Z R.:. x_dim R.:. y_dim R.:. z_dim)
+      let fvs = R.fromListUnboxed new_shape elements
+      return $ Grid fvs
+
 
 
 -- | Takes a grid and a position as an argument and returns the cube
---   centered on that position. If there is no cube there (i.e. the
---   position is outside of the grid), it will throw an error.
+--   centered on that position. If there is no cube there, well, you
+--   shouldn't have done that. The omitted "otherwise" case actually
+--   does improve performance.
 cube_at :: Grid -> Int -> Int -> Int -> Cube
-cube_at g i j k
-    | i < 0      = error "i < 0 in cube_at"
-    | i >= xsize = error "i >= xsize in cube_at"
-    | j < 0      = error "j < 0 in cube_at"
-    | j >= ysize = error "j >= ysize in cube_at"
-    | k < 0      = error "k < 0 in cube_at"
-    | k >= zsize = error "k >= zsize in cube_at"
-    | otherwise = (cube_grid g) ! (i,j,k)
-      where
-        fvs = function_values g
-        (xsize, ysize, zsize) = dims fvs
-
---   The first cube along any axis covers (-h/2, h/2). The second
---   covers (h/2, 3h/2).  The third, (3h/2, 5h/2), and so on.
+cube_at !g !i !j !k =
+   Cube i j k fvs' tet_vol
+   where
+     fvs = function_values g
+     fvs' = make_values fvs i j k
+     tet_vol = 1/24
+
+
+--   The first cube along any axis covers (-1/2, 1/2). The second
+--   covers (1/2, 3/2).  The third, (3/2, 5/2), and so on.
 --
---   We translate the (x,y,z) coordinates forward by 'h/2' so that the
---   first covers (0, h), the second covers (h, 2h), etc. This makes
+--   We translate the (x,y,z) coordinates forward by 1/2 so that the
+--   first covers (0, 1), the second covers (1, 2), etc. This makes
 --   it easy to figure out which cube contains the given point.
 calculate_containing_cube_coordinate :: Grid -> Double -> Int
 calculate_containing_cube_coordinate g coord
@@ -112,59 +86,57 @@ calculate_containing_cube_coordinate g coord
     -- exists.
     | coord < offset = 0
     | coord == offset && (xsize > 1 && ysize > 1 && zsize > 1) = 1
-    | otherwise = (ceiling ( (coord + offset) / cube_width )) - 1
+    | otherwise = (ceiling (coord + offset)) - 1
     where
       (xsize, ysize, zsize) = dims (function_values g)
-      cube_width = (h g)
-      offset = cube_width / 2
+      offset = 1/2
 
 
 -- | Takes a 'Grid', and returns a 'Cube' containing the given 'Point'.
 --   Since our grid is rectangular, we can figure this out without having
 --   to check every cube.
 find_containing_cube :: Grid -> Point -> Cube
-find_containing_cube g p =
+find_containing_cube g (Point x y z) =
     cube_at g i j k
     where
-      (x, y, z) = p
       i = calculate_containing_cube_coordinate g x
       j = calculate_containing_cube_coordinate g y
       k = calculate_containing_cube_coordinate g z
 
 
-{-# INLINE zoom_lookup #-}
-zoom_lookup :: Grid -> PolynomialArray -> ScaleFactor -> a -> (R.DIM3 -> Double)
-zoom_lookup g polynomials scale_factor _ =
-    zoom_result g polynomials scale_factor
+zoom_lookup :: Values3D -> ScaleFactor -> a -> (R.DIM3 -> Double)
+zoom_lookup v3d scale_factor _ =
+    zoom_result v3d scale_factor
 
 
-{-# INLINE zoom_result #-}
-zoom_result :: Grid -> PolynomialArray -> ScaleFactor -> R.DIM3 -> Double
-zoom_result g polynomials (sfx, sfy, sfz) (R.Z R.:. m R.:. n R.:. o) =
-  (polynomials ! (i, j, k, (number t))) p
+zoom_result :: Values3D -> ScaleFactor -> R.DIM3 -> Double
+zoom_result v3d (sfx, sfy, sfz) (R.Z R.:. m R.:. n R.:. o) =
+  f p
   where
-    offset = (h g)/2
+    g = Grid v3d
+    offset = 1/2
     m' = (fromIntegral m) / (fromIntegral sfx) - offset
     n' = (fromIntegral n) / (fromIntegral sfy) - offset
     o' = (fromIntegral o) / (fromIntegral sfz) - offset
-    p  = (m', n', o') :: Point
+    p  = Point m' n' o'
     cube = find_containing_cube g p
-    -- Figure out i,j,k without importing those functions.
-    Cube _ i j k _ _ = cube
     t = find_containing_tetrahedron cube p
-
-    
-zoom :: Grid -> PolynomialArray -> ScaleFactor -> Values3D
-zoom g polynomials scale_factor
-    | xsize == 0 || ysize == 0 || zsize == 0 = empty3d
-    | otherwise =
-        R.force $ R.traverse arr transExtent (zoom_lookup g polynomials scale_factor)
-          where
-            arr = function_values g
-            (xsize, ysize, zsize) = dims arr
-            transExtent = zoom_shape scale_factor
+    f = polynomial t
 
 
+--
+-- Instead of IO, we could get away with a generic monad 'm'
+-- here. However, /we/ only call this function from within IO.
+--
+zoom :: Values3D -> ScaleFactor -> IO Values3D
+zoom v3d scale_factor
+    | xsize == 0 || ysize == 0 || zsize == 0 = return empty3d
+    | otherwise =
+        R.computeUnboxedP $ R.unsafeTraverse v3d transExtent f
+        where
+          (xsize, ysize, zsize) = dims v3d
+          transExtent = zoom_shape scale_factor
+          f = zoom_lookup v3d scale_factor
 
 
 -- | Check all coefficients of tetrahedron0 belonging to the cube
@@ -206,7 +178,7 @@ trilinear_c0_t0_tests =
        testCase "v3 is correct" test_trilinear_f0_t0_v3]
     ]
   where
-    g = make_grid 1 trilinear
+    g = Grid trilinear
     cube = cube_at g 1 1 1
     t = tetrahedron cube 0
 
@@ -292,60 +264,62 @@ trilinear_c0_t0_tests =
 
     test_trilinear_f0_t0_v0 :: Assertion
     test_trilinear_f0_t0_v0 =
-      assertEqual "v0 is correct" (v0 t) (1, 1, 1)
+      assertEqual "v0 is correct" (v0 t) (Point 1 1 1)
 
     test_trilinear_f0_t0_v1 :: Assertion
     test_trilinear_f0_t0_v1 =
-      assertEqual "v1 is correct" (v1 t) (0.5, 1, 1)
+      assertEqual "v1 is correct" (v1 t) (Point 0.5 1 1)
 
     test_trilinear_f0_t0_v2 :: Assertion
     test_trilinear_f0_t0_v2 =
-      assertEqual "v2 is correct" (v2 t) (0.5, 0.5, 1.5)
+      assertEqual "v2 is correct" (v2 t) (Point 0.5 0.5 1.5)
 
     test_trilinear_f0_t0_v3 :: Assertion
     test_trilinear_f0_t0_v3 =
-      assertClose "v3 is correct" (v3 t) (0.5, 1.5, 1.5)
+      assertEqual "v3 is correct" (v3 t) (Point 0.5 1.5 1.5)
 
 
 test_trilinear_reproduced :: Assertion
 test_trilinear_reproduced =
     assertTrue "trilinears are reproduced correctly" $
-             and [p (i', j', k') ~= value_at trilinear i j k
+             and [p (Point i' j' k') ~= value_at trilinear i j k
                     | i <- [0..2],
                       j <- [0..2],
                       k <- [0..2],
+                      c0 <- cs,
                       t <- tetrahedra c0,
                       let p = polynomial t,
                       let i' = fromIntegral i,
                       let j' = fromIntegral j,
                       let k' = fromIntegral k]
     where
-      g = make_grid 1 trilinear
-      c0 = cube_at g 1 1 1
+      g = Grid trilinear
+      cs = [ cube_at g ci cj ck | ci <- [0..2], cj <- [0..2], ck <- [0..2] ]
 
 
 test_zeros_reproduced :: Assertion
 test_zeros_reproduced =
     assertTrue "the zero function is reproduced correctly" $
-             and [p (i', j', k') ~= value_at zeros i j k
+             and [p (Point i' j' k') ~= value_at zeros i j k
                     | i <- [0..2],
                       j <- [0..2],
                       k <- [0..2],
                       let i' = fromIntegral i,
                       let j' = fromIntegral j,
-                      let k' = fromIntegral k]
+                      let k' = fromIntegral k,
+                      c0 <- cs,
+                      t0 <- tetrahedra c0,
+                      let p = polynomial t0 ]
     where
-      g = make_grid 1 zeros
-      c0 = cube_at g 1 1 1
-      t0 = tetrahedron c0 0
-      p = polynomial t0
+      g = Grid zeros
+      cs = [ cube_at g ci cj ck | ci <- [0..2], cj <- [0..2], ck <- [0..2] ]
 
 
 -- | Make sure we can reproduce a 9x9x9 trilinear from the 3x3x3 one.
 test_trilinear9x9x9_reproduced :: Assertion
 test_trilinear9x9x9_reproduced =
     assertTrue "trilinear 9x9x9 is reproduced correctly" $
-      and [p (i', j', k') ~= value_at trilinear9x9x9 i j k
+      and [p (Point i' j' k') ~= value_at trilinear9x9x9 i j k
             | i <- [0..8],
               j <- [0..8],
               k <- [0..8],
@@ -355,40 +329,20 @@ test_trilinear9x9x9_reproduced =
               let j' = (fromIntegral j) * 0.5,
               let k' = (fromIntegral k) * 0.5]
     where
-      g = make_grid 1 trilinear
+      g = Grid trilinear
       c0 = cube_at g 1 1 1
 
 
--- | The point 'p' in this test lies on the boundary of tetrahedra 12 and 15.
---   However, the 'contains_point' test fails due to some numerical innacuracy.
---   This bug should have been fixed by setting a positive tolerance level.
---
---   Example from before the fix:
---
---   b1 (tetrahedron c 20) (0, 17.5, 0.5)
---   -0.0
---
-test_tetrahedra_collision_sensitivity :: Assertion
-test_tetrahedra_collision_sensitivity =
-  assertTrue "tetrahedron collision tests isn't too sensitive" $
-             contains_point t20 p
-  where
-    g = make_grid 1 naturals_1d
-    cube = cube_at g 0 18 0
-    p = (0, 17.5, 0.5) :: Point
-    t20 = tetrahedron cube 20
-
 
 prop_cube_indices_never_go_out_of_bounds :: Grid -> Gen Bool
 prop_cube_indices_never_go_out_of_bounds g =
   do
-    let delta = Grid.h g
-    let coordmin = negate (delta/2)
+    let coordmin = negate (1/2)
 
     let (xsize, ysize, zsize) = dims $ function_values g
-    let xmax = delta*(fromIntegral xsize) - (delta/2)
-    let ymax = delta*(fromIntegral ysize) - (delta/2)
-    let zmax = delta*(fromIntegral zsize) - (delta/2)
+    let xmax = (fromIntegral xsize) - (1/2)
+    let ymax = (fromIntegral ysize) - (1/2)
+    let zmax = (fromIntegral zsize) - (1/2)
 
     x <- choose (coordmin, xmax)
     y <- choose (coordmin, ymax)
@@ -407,21 +361,125 @@ prop_cube_indices_never_go_out_of_bounds g =
       idx_z <= zsize - 1
 
 
+-- | Given in Sorokina and Zeilfelder, p. 80, (2.9). Note that the
+--   third and fourth indices of c-t10 have been switched. This is
+--   because we store the triangles oriented such that their volume is
+--   positive. If T and T-tilde share \<v1,v2,v3\> and v0,v0-tilde point
+--   in opposite directions, one of them has to have negative volume!
+prop_c0120_identity :: Grid -> Property
+prop_c0120_identity g =
+  xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
+    c t0 0 1 2 0 ~= (c t0 1 0 2 0 + c t10 1 0 0 2) / 2
+  where
+    fvs = function_values g
+    (xsize, ysize, zsize) = dims fvs
+    cube0 = cube_at g 1 1 1
+    cube1 = cube_at g 0 1 1
+    t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
+    t10 = tetrahedron cube1 10
+
+
+-- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
+--   'prop_c0120_identity'.
+prop_c0111_identity :: Grid -> Property
+prop_c0111_identity g =
+  xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
+    c t0 0 1 1 1 ~= (c t0 1 0 1 1 + c t10 1 0 1 1) / 2
+  where
+    fvs = function_values g
+    (xsize, ysize, zsize) = dims fvs
+    cube0 = cube_at g 1 1 1
+    cube1 = cube_at g 0 1 1
+    t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
+    t10 = tetrahedron cube1 10
+
+
+-- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
+--   'prop_c0120_identity'.
+prop_c0201_identity :: Grid -> Property
+prop_c0201_identity g =
+  xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
+    c t0 0 2 0 1 ~= (c t0 1 1 0 1 + c t10 1 1 1 0) / 2
+  where
+    fvs = function_values g
+    (xsize, ysize, zsize) = dims fvs
+    cube0 = cube_at g 1 1 1
+    cube1 = cube_at g 0 1 1
+    t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
+    t10 = tetrahedron cube1 10
+
+
+-- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
+--   'prop_c0120_identity'.
+prop_c0102_identity :: Grid -> Property
+prop_c0102_identity g =
+  xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
+    c t0 0 1 0 2 ~= (c t0 1 0 0 2 + c t10 1 0 2 0) / 2
+  where
+    fvs = function_values g
+    (xsize, ysize, zsize) = dims fvs
+    cube0 = cube_at g 1 1 1
+    cube1 = cube_at g 0 1 1
+    t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
+    t10 = tetrahedron cube1 10
+
+
+-- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
+--   'prop_c0120_identity'.
+prop_c0210_identity :: Grid -> Property
+prop_c0210_identity g =
+  xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
+    c t0 0 2 1 0 ~= (c t0 1 1 1 0 + c t10 1 1 0 1) / 2
+  where
+    fvs = function_values g
+    (xsize, ysize, zsize) = dims fvs
+    cube0 = cube_at g 1 1 1
+    cube1 = cube_at g 0 1 1
+    t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
+    t10 = tetrahedron cube1 10
+
+
+-- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
+--   'prop_c0120_identity'.
+prop_c0300_identity :: Grid -> Property
+prop_c0300_identity g =
+  xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
+    c t0 0 3 0 0 ~= (c t0 1 2 0 0 + c t10 1 2 0 0) / 2
+  where
+    fvs = function_values g
+    (xsize, ysize, zsize) = dims fvs
+    cube0 = cube_at g 1 1 1
+    cube1 = cube_at g 0 1 1
+    t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
+    t10 = tetrahedron cube1 10
+
+
+-- | All of the properties from Section (2.9), p. 80. These require a
+--   grid since they refer to two adjacent cubes.
+p80_29_properties :: Test.Framework.Test
+p80_29_properties =
+  testGroup "p. 80, Section (2.9) Properties" [
+    testProperty "c0120 identity" prop_c0120_identity,
+    testProperty "c0111 identity" prop_c0111_identity,
+    testProperty "c0201 identity" prop_c0201_identity,
+    testProperty "c0102 identity" prop_c0102_identity,
+    testProperty "c0210 identity" prop_c0210_identity,
+    testProperty "c0300 identity" prop_c0300_identity ]
+
 
 grid_tests :: Test.Framework.Test
 grid_tests =
     testGroup "Grid Tests" [
       trilinear_c0_t0_tests,
-      testCase "tetrahedra collision test isn't too sensitive"
-         test_tetrahedra_collision_sensitivity,
-      testCase "trilinear reproduced" test_trilinear_reproduced,
-      testCase "zeros reproduced" test_zeros_reproduced ]
+      p80_29_properties,
+      testProperty "cube indices within bounds"
+        prop_cube_indices_never_go_out_of_bounds ]
 
 
 -- Do the slow tests last so we can stop paying attention.
 slow_tests :: Test.Framework.Test
 slow_tests =
     testGroup "Slow Tests" [
-      testProperty "cube indices within bounds"
-                   prop_cube_indices_never_go_out_of_bounds,
-      testCase "trilinear9x9x9 reproduced" test_trilinear9x9x9_reproduced ]
+      testCase "trilinear reproduced" test_trilinear_reproduced,
+      testCase "trilinear9x9x9 reproduced" test_trilinear9x9x9_reproduced,
+      testCase "zeros reproduced" test_zeros_reproduced ]