Merge branch 'uncached_vector'

Conflicts:
	spline3/src/Main.hs

diff --git a/src/Grid.hs b/src/Grid.hs
index 893fe5c360d3f39c73a99ebab26d1320fd798815..742360f5d799993f47ce2cb38adf3f3cfc46ae6f 100644 (file)
--- a/src/Grid.hs
+++ b/src/Grid.hs
@@ -27,9 +27,8 @@ import Cube (Cube(Cube),
 import Examples
 import FunctionValues
 import Point (Point)
-import PolynomialArray (PolynomialArray)
 import ScaleFactor
-import Tetrahedron (Tetrahedron, c, number, polynomial, v0, v1, v2, v3)
+import Tetrahedron (Tetrahedron, c, polynomial, v0, v1, v2, v3)
 import ThreeDimensional
 import Values (Values3D, dims, empty3d, zoom_shape)
 
@@ -133,15 +132,15 @@ find_containing_cube g p =
 
 
 {-# 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 :: Grid -> ScaleFactor -> a -> (R.DIM3 -> Double)
+zoom_lookup g scale_factor _ =
+    zoom_result g 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 :: Grid -> ScaleFactor -> R.DIM3 -> Double
+zoom_result g (sfx, sfy, sfz) (R.Z R.:. m R.:. n R.:. o) =
+  f p
   where
     offset = (h g)/2
     m' = (fromIntegral m) / (fromIntegral sfx) - offset
@@ -149,16 +148,14 @@ zoom_result g polynomials (sfx, sfy, sfz) (R.Z R.:. m R.:. n R.:. o) =
     o' = (fromIntegral o) / (fromIntegral sfz) - offset
     p  = (m', n', o') :: Point
     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
-
+    f = polynomial t
     
-zoom :: Grid -> PolynomialArray -> ScaleFactor -> Values3D
-zoom g polynomials scale_factor
+zoom :: Grid -> ScaleFactor -> Values3D
+zoom g scale_factor
     | xsize == 0 || ysize == 0 || zsize == 0 = empty3d
     | otherwise =
-        R.force $ R.traverse arr transExtent (zoom_lookup g polynomials scale_factor)
+        R.force $ R.unsafeTraverse arr transExtent (zoom_lookup g scale_factor)
           where
             arr = function_values g
             (xsize, ysize, zsize) = dims arr

diff --git a/src/Main.hs b/src/Main.hs
index 963b31ed70448c8c88860ab039b9c1d85fbd401d..8a1efdb76de21e6b94421e63e6ebd07630ff5990 100644 (file)
--- a/src/Main.hs
+++ b/src/Main.hs
@@ -9,10 +9,7 @@ import Data.Array.Repa (
 
 import System.Environment (getArgs)
 
-import Cube (tetrahedron)
-import Grid (cube_at, make_grid, zoom)
-import PolynomialArray (make_polynomial_array)
-import Tetrahedron (polynomial)
+import Grid (make_grid, zoom)
 import Values (read_values_3d, write_values_1d)
 
 mri_shape :: DIM3
@@ -30,13 +27,5 @@ main = do
   mridata <- read_values_3d mri_shape in_file
 
   let g = make_grid 1 mridata
-  let polynomials = make_polynomial_array (255,255,0,23)
-        [ ((i,j,k,tet), polynomial t) | i <- [0..255],
-                                        j <- [0..255],
-                                        k <- [0],
-                                        tet <- [0..23],
-                                        let c = cube_at g i j k,
-                                        let t = tetrahedron c tet ]
-
-  let output = zoom g polynomials (4,4,1)
+  let output = zoom g (4,4,1)
   write_values_1d output out_file

diff --git a/src/PolynomialArray.hs b/src/PolynomialArray.hs
deleted file mode 100644 (file)
index 2cf96dd..0000000
--- a/src/PolynomialArray.hs
+++ /dev/null
@@ -1,26 +0,0 @@
--- | The PolynomialArray module contains one type, 'PolynomialArray',
---   and a smart constructor for it. The array is essentially a hash;
---   the unique keys can be thought of as a cube/tetrahedron pair, and
---   the associated value is the polynomial function over that
---   tetrahedron. Rather than define 'Ord' instances and a proper map,
---   we just use the indices i,j,k along with the tetrahedron number
---   as the keys.
-module PolynomialArray (
-  PolynomialArray,
-  make_polynomial_array
-  )
-where
-
-import RealFunction(RealFunction)
-import Point (Point)
-
-import Data.Array (Array, array)
-
-type PolynomialArray =
-    Array (Int,Int,Int,Int) (RealFunction Point)
-
-make_polynomial_array :: (Int, Int, Int, Int) ->
-                        [( (Int,Int,Int,Int), (RealFunction Point) )] ->
-                        PolynomialArray
-make_polynomial_array (max_i, max_j, max_k, max_tetrahedron_number) =
-  array ( (0,0,0,0), (max_i,max_j,max_k, max_tetrahedron_number) )

diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs
index 5b459fea2870de5c2fed88ee9f9e5b57c7e5cfc9..87332b94440afdb999b173c13f681095d48ac796 100644 (file)
--- a/src/Tetrahedron.hs
+++ b/src/Tetrahedron.hs
@@ -1,6 +1,11 @@
 module Tetrahedron
 where
 
+import qualified Data.Vector as V (
+  singleton,
+  snoc,
+  sum
+  )
 import Numeric.LinearAlgebra hiding (i, scale)
 import Prelude hiding (LT)
 import Test.QuickCheck (Arbitrary(..), Gen, choose)
@@ -86,11 +91,26 @@ instance ThreeDimensional Tetrahedron where
 
 polynomial :: Tetrahedron -> (RealFunction Point)
 polynomial t =
-    sum [ (c t i j k l) `cmult` (beta t i j k l) | i <- [0..3],
-                                                   j <- [0..3],
-                                                   k <- [0..3],
-                                                   l <- [0..3],
-                                                   i + j + k + l == 3]
+    V.sum $ V.singleton ((c t 0 0 0 3) `cmult` (beta t 0 0 0 3)) `V.snoc`
+            ((c t 0 0 1 2) `cmult` (beta t 0 0 1 2)) `V.snoc`
+            ((c t 0 0 2 1) `cmult` (beta t 0 0 2 1)) `V.snoc`
+            ((c t 0 0 3 0) `cmult` (beta t 0 0 3 0)) `V.snoc`
+            ((c t 0 1 0 2) `cmult` (beta t 0 1 0 2)) `V.snoc`
+            ((c t 0 1 1 1) `cmult` (beta t 0 1 1 1)) `V.snoc`
+            ((c t 0 1 2 0) `cmult` (beta t 0 1 2 0)) `V.snoc`
+            ((c t 0 2 0 1) `cmult` (beta t 0 2 0 1)) `V.snoc`
+            ((c t 0 2 1 0) `cmult` (beta t 0 2 1 0)) `V.snoc`
+            ((c t 0 3 0 0) `cmult` (beta t 0 3 0 0)) `V.snoc`
+            ((c t 1 0 0 2) `cmult` (beta t 1 0 0 2)) `V.snoc`
+            ((c t 1 0 1 1) `cmult` (beta t 1 0 1 1)) `V.snoc`
+            ((c t 1 0 2 0) `cmult` (beta t 1 0 2 0)) `V.snoc`
+            ((c t 1 1 0 1) `cmult` (beta t 1 1 0 1)) `V.snoc`
+            ((c t 1 1 1 0) `cmult` (beta t 1 1 1 0)) `V.snoc`
+            ((c t 1 2 0 0) `cmult` (beta t 1 2 0 0)) `V.snoc`
+            ((c t 2 0 0 1) `cmult` (beta t 2 0 0 1)) `V.snoc`
+            ((c t 2 0 1 0) `cmult` (beta t 2 0 1 0)) `V.snoc`
+            ((c t 2 1 0 0) `cmult` (beta t 2 1 0 0)) `V.snoc`
+            ((c t 3 0 0 0) `cmult` (beta t 3 0 0 0))
 
 
 -- | Returns the domain point of t with indices i,j,k,l.

diff --git a/src/Values.hs b/src/Values.hs
index 8e2b61e789e9cb4fc86fc4921055e1021c5cecc3..070abbbf2554ea5eec5276e5496fe1554fcaf3d3 100644 (file)
--- a/src/Values.hs
+++ b/src/Values.hs
@@ -12,7 +12,7 @@ import Data.Array.Repa (
   DIM3,
   extent,
   fromList,
-  index,
+  unsafeIndex,
   reshape,
   size
   )
@@ -72,7 +72,7 @@ dims v3d =
 
 idx :: Values3D -> Int -> Int -> Int -> Double
 idx v3d i j k =
-    index v3d shape
+    unsafeIndex v3d shape
     where
       shape :: DIM3
       shape = (Z :. i :. j :. k)