From 883ce9d78072c492000de94478189095032b6615 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Tue, 23 Aug 2011 14:40:35 -0400 Subject: [PATCH] Begin updating everything to use Repa arrays (Values3D). --- src/FunctionValues.hs | 15 +++++++------ src/Grid.hs | 28 +++++++++++------------- src/Main.hs | 11 ++++------ src/Values.hs | 50 +++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 76 insertions(+), 28 deletions(-) diff --git a/src/FunctionValues.hs b/src/FunctionValues.hs index 1fbc044..1d722b1 100644 --- a/src/FunctionValues.hs +++ b/src/FunctionValues.hs @@ -7,6 +7,7 @@ import Prelude hiding (LT) import Test.QuickCheck (Arbitrary(..), choose) import Cardinal +import Values (Values3D, dims, idx) -- | The FunctionValues type represents the value of our function f at -- the 27 points surrounding (and including) the center of a @@ -162,21 +163,23 @@ eval f (Quotient x y) = (eval f x) / (eval f y) -- | Takes a three-dimensional list of 'Double' and a set of 3D -- coordinates (i,j,k), and returns the value at (i,j,k) in the -- supplied list. If there is no such value, zero is returned. -value_at :: [[[Double]]] -> Int -> Int -> Int -> Double +value_at :: Values3D -> Int -> Int -> Int -> Double value_at values i j k | i < 0 = 0 | j < 0 = 0 | k < 0 = 0 - | length values <= k = 0 - | length (values !! k) <= j = 0 - | length ((values !! k) !! j) <= i = 0 - | otherwise = ((values !! k) !! j) !! i + | xsize <= i = 0 + | ysize <= j = 0 + | zsize <= k = 0 + | otherwise = idx values i j k + where + (xsize, ysize, zsize) = dims values -- | Given a three-dimensional list of 'Double' and a set of 3D -- coordinates (i,j,k), constructs and returns the 'FunctionValues' -- object centered at (i,j,k) -make_values :: [[[Double]]] -> Int -> Int -> Int -> FunctionValues +make_values :: Values3D -> Int -> Int -> Int -> FunctionValues make_values values i j k = empty_values { front = value_at values (i-1) j k, back = value_at values (i+1) j k, diff --git a/src/Grid.hs b/src/Grid.hs index 63b2cfa..1a01ab7 100644 --- a/src/Grid.hs +++ b/src/Grid.hs @@ -12,27 +12,27 @@ import Misc (flatten) import Point (Point) import Tetrahedron (polynomial) import ThreeDimensional (contains_point) - +import Values (Values3D, dims, empty3d) -- | 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 :: [[[Double]]] } + function_values :: Values3D } deriving (Eq, Show) instance Arbitrary Grid where arbitrary = do (Positive h') <- arbitrary :: Gen (Positive Double) - fvs <- arbitrary :: Gen [[[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 -> [[[Double]]] -> Grid +make_grid :: Double -> Values3D -> Grid make_grid grid_size values | grid_size <= 0 = error "grid size must be positive" | otherwise = Grid grid_size values @@ -40,24 +40,21 @@ make_grid grid_size values -- | Creates an empty grid with grid size 1. empty_grid :: Grid -empty_grid = Grid 1 [[[]]] +empty_grid = Grid 1 empty3d -- | Returns a three-dimensional list of cubes centered on the grid -- points of g with the appropriate 'FunctionValues'. cubes :: Grid -> [[[Cube]]] cubes g - | fvs == [[[]]] = [[[]]] - | head fvs == [[]] = [[[]]] + | xsize == 0 || ysize == 0 || zsize == 0 = [[[]]] | otherwise = [[[ Cube (h g) i j k (make_values fvs i j k) | i <- [0..xsize]] | j <- [0..ysize]] | k <- [0..zsize]] where fvs = function_values g - zsize = (length fvs) - 1 - ysize = length (head fvs) - 1 - xsize = length (head $ head fvs) - 1 + (xsize, ysize, zsize) = dims fvs -- | Takes a grid and a position as an argument and returns the cube @@ -87,8 +84,7 @@ find_containing_cubes g p = zoom :: Grid -> Int -> [[[Double]]] zoom g scale_factor - | fvs == [[[]]] = [] - | head fvs == [[]] = [] + | xsize == 0 || ysize == 0 || zsize == 0 = [] | otherwise = [[[f p | i <- [0..scaled_zsize], let i' = scale_dimension i, @@ -105,6 +101,8 @@ zoom g scale_factor scale_dimension x = (fromIntegral x) / (fromIntegral scale_factor) fvs = function_values g - scaled_zsize = ((length fvs) - 1) * scale_factor - scaled_ysize = (length (head fvs) - 1) * scale_factor - scaled_xsize = (length (head $ head fvs) - 1) * scale_factor + (xsize, ysize, zsize) = dims fvs + scaled_xsize = xsize * scale_factor + scaled_ysize = ysize * scale_factor + scaled_zsize = zsize * scale_factor + diff --git a/src/Main.hs b/src/Main.hs index bc63265..7b50de9 100644 --- a/src/Main.hs +++ b/src/Main.hs @@ -5,11 +5,10 @@ import Data.Array.Repa ( DIM3, Z(..), (:.)(..), - toList ) import Values ---import Grid(make_grid, zoom) +import Grid(make_grid, zoom) mri_shape :: DIM3 mri_shape = (Z :. 256 :. 256 :. 109) @@ -17,8 +16,6 @@ mri_shape = (Z :. 256 :. 256 :. 109) main :: IO () main = do mridata <- read_values_3d mri_shape "./data/mridata.txt" - let tmp = Data.Array.Repa.toList mridata - print tmp --- let g = make_grid 1 (Data.Array.Repa.toList mridata2) --- let output = zoom g 2 --- print "Hello, world." + let g = make_grid 1 mridata + let output = zoom g 1 + print output diff --git a/src/Values.hs b/src/Values.hs index 10d7fb9..bd7c4d6 100644 --- a/src/Values.hs +++ b/src/Values.hs @@ -1,24 +1,74 @@ +{-# LANGUAGE TypeSynonymInstances #-} + module Values where import Data.Array.Repa ( Array, + Z(..), + (:.)(..), DIM1, + DIM2, DIM3, + extent, + fromList, + index, reshape, ) import Data.Array.Repa.IO.Vector (readVectorFromTextFile) import System.FilePath () +import Test.QuickCheck (Arbitrary(..), Gen) type Values1D = Array DIM1 Double +type Values2D = Array DIM2 Double type Values3D = Array DIM3 Double + +instance Arbitrary Values3D where + arbitrary = do + x_dim <- arbitrary :: Gen Int + y_dim <- arbitrary :: Gen Int + z_dim <- arbitrary :: Gen Int + one_d <- arbitrary :: Gen Values1D + let new_shape = (Z :. x_dim :. y_dim :. z_dim) + let three_d = reshape new_shape one_d + return three_d + + +instance Arbitrary Values1D where + arbitrary = do + x <- arbitrary :: Gen [Double] + let shape = (Z :. (length x)) + let one_d = Data.Array.Repa.fromList shape x + return one_d + + read_values_1d :: FilePath -> IO Values1D read_values_1d path = readVectorFromTextFile path + read_values_3d :: DIM3 -> FilePath -> IO Values3D read_values_3d sh path = do one_d <- read_values_1d path return $ reshape sh one_d + + +empty3d :: Values3D +empty3d = Data.Array.Repa.fromList (Z :. 0 :. 0 :. 0) [] + + +dims :: Values3D -> (Int, Int, Int) +dims v3d = + let (Z :. x :. y :. z) = extent v3d + in + (x,y,z) + + +idx :: Values3D -> Int -> Int -> Int -> Double +idx v3d i j k = + index v3d shape + where + shape :: DIM3 + shape = (Z :. i :. j :. k) -- 2.43.2