{-# LANGUAGE BangPatterns #-} -- | The FunctionValues module contains the 'FunctionValues' type and -- the functions used to manipulate it. -- module FunctionValues ( FunctionValues(..), empty_values, eval, make_values, rotate, function_values_tests, function_values_properties, value_at ) where import Prelude hiding ( LT ) import Test.Tasty ( TestTree, testGroup ) import Test.Tasty.HUnit ( Assertion, testCase ) import Test.Tasty.QuickCheck ( Arbitrary(..), choose, testProperty ) import Assertions ( assertTrue ) import Cardinal ( Cardinal(..), cwx, cwy, cwz ) import Examples ( trilinear ) 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 -- cube. Each value of f can be accessed by the name of its -- direction. -- data FunctionValues = FunctionValues { front :: !Double, back :: !Double, left :: !Double, right :: !Double, top :: !Double, down :: !Double, front_left :: !Double, front_right :: !Double, front_down :: !Double, front_top :: !Double, back_left :: !Double, back_right :: !Double, back_down :: !Double, back_top :: !Double, left_down :: !Double, left_top :: !Double, right_down :: !Double, right_top :: !Double, front_left_down :: !Double, front_left_top :: !Double, front_right_down :: !Double, front_right_top :: !Double, back_left_down :: !Double, back_left_top :: !Double, back_right_down :: !Double, back_right_top :: !Double, interior :: !Double } deriving (Eq, Show) instance Arbitrary FunctionValues where arbitrary = do front' <- choose (min_double, max_double) back' <- choose (min_double, max_double) left' <- choose (min_double, max_double) right' <- choose (min_double, max_double) top' <- choose (min_double, max_double) down' <- choose (min_double, max_double) front_left' <- choose (min_double, max_double) front_right' <- choose (min_double, max_double) front_top' <- choose (min_double, max_double) front_down' <- choose (min_double, max_double) back_left' <- choose (min_double, max_double) back_right' <- choose (min_double, max_double) back_top' <- choose (min_double, max_double) back_down' <- choose (min_double, max_double) left_top' <- choose (min_double, max_double) left_down' <- choose (min_double, max_double) right_top' <- choose (min_double, max_double) right_down' <- choose (min_double, max_double) front_left_top' <- choose (min_double, max_double) front_left_down' <- choose (min_double, max_double) front_right_top' <- choose (min_double, max_double) front_right_down' <- choose (min_double, max_double) back_left_top' <- choose (min_double, max_double) back_left_down' <- choose (min_double, max_double) back_right_top' <- choose (min_double, max_double) back_right_down' <- choose (min_double, max_double) interior' <- choose (min_double, max_double) return empty_values { front = front', back = back', left = left', right = right', top = top', down = down', front_left = front_left', front_right = front_right', front_top = front_top', front_down = front_down', back_left = back_left', back_right = back_right', back_top = back_top', back_down = back_down', left_top = left_top', left_down = left_down', right_top = right_top', right_down = right_down', front_left_top = front_left_top', front_left_down = front_left_down', front_right_top = front_right_top', front_right_down = front_right_down', back_left_top = back_left_top', back_left_down = back_left_down', back_right_top = back_right_top', back_right_down = back_right_down', interior = interior' } where -- | We perform addition with the function values contained in a -- FunctionValues object. If we choose random doubles near the machine -- min/max, we risk overflowing or underflowing the 'Double'. This -- places a reasonably safe limit on the maximum size of our generated -- 'Double' members. max_double :: Double max_double = 10000.0 -- | See 'max_double'. min_double :: Double min_double = (-1) * max_double -- | Return a 'FunctionValues' with a bunch of zeros for data points. empty_values :: FunctionValues empty_values = FunctionValues 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -- | The eval function is where the magic happens for the -- FunctionValues type. Given a 'Cardinal' direction and a -- 'FunctionValues' object, eval will return the value of the -- function f in that 'Cardinal' direction. Note that 'Cardinal' can -- be a composite type; eval is what performs the \"arithmetic\" on -- 'Cardinal' directions. eval :: FunctionValues -> Cardinal -> Double eval f F = front f eval f B = back f eval f L = left f eval f R = right f eval f T = top f eval f D = down f eval f FL = front_left f eval f FR = front_right f eval f FD = front_down f eval f FT = front_top f eval f BL = back_left f eval f BR = back_right f eval f BD = back_down f eval f BT = back_top f eval f LD = left_down f eval f LT = left_top f eval f RD = right_down f eval f RT = right_top f eval f FLD = front_left_down f eval f FLT = front_left_top f eval f FRD = front_right_down f eval f FRT = front_right_top f eval f BLD = back_left_down f eval f BLT = back_left_top f eval f BRD = back_right_down f eval f BRT = back_right_top f eval f I = interior f eval _ (Scalar x) = x eval f (Sum x y) = (eval f x) + (eval f y) eval f (Difference x y) = (eval f x) - (eval f y) eval f (Product x y) = (eval f x) * (eval f y) 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, we calculate one -- according to Sorokina and Zeilfelder, remark 7.3, p. 99. -- -- We specifically do not consider values more than one unit away -- from our grid. -- -- Examples: -- -- >>> value_at Examples.trilinear 0 0 0 -- 1.0 -- -- >>> value_at Examples.trilinear (-1) 0 0 -- 0.0 -- -- >>> value_at Examples.trilinear 0 0 4 -- 1.0 -- -- >>> value_at Examples.trilinear 1 3 0 -- 5.0 -- value_at :: Values3D -> Int -> Int -> Int -> Double value_at v3d !i !j !k -- Put the most common case first! | (valid_i i) && (valid_j j) && (valid_k k) = idx v3d i j k -- The next three are from the first line in (7.3). Analogous cases -- have been added where the indices are one-too-big. These are the -- "one index is bad" cases. | not (valid_i i) = if (dim_i == 1) then -- We're one-dimensional in our first coordinate, so just -- return the data point that we do have. If we try to use -- the formula from remark 7.3, we go into an infinite loop. value_at v3d 0 j k else if (i == -1) then 2*(value_at v3d 0 j k) - (value_at v3d 1 j k) else 2*(value_at v3d (i-1) j k) - (value_at v3d (i-2) j k) | not (valid_j j) = if (dim_j == 1) then -- We're one-dimensional in our second coordinate, so just -- return the data point that we do have. If we try to use -- the formula from remark 7.3, we go into an infinite loop. value_at v3d i 0 k else if (j == -1) then 2*(value_at v3d i 0 k) - (value_at v3d i 1 k) else 2*(value_at v3d i (j-1) k) - (value_at v3d i (j-2) k) | not (valid_k k) = if (dim_k == 1) then -- We're one-dimensional in our third coordinate, so just -- return the data point that we do have. If we try to use -- the formula from remark 7.3, we go into an infinite loop. value_at v3d i j 0 else if (k == -1) then 2*(value_at v3d i j 0) - (value_at v3d i j 1) else 2*(value_at v3d i j (k-1)) - (value_at v3d i j (k-2)) where (dim_i, dim_j, dim_k) = dims v3d valid_i :: Int -> Bool valid_i i' = (i' >= 0) && (i' < dim_i) valid_j :: Int -> Bool valid_j j' = (j' >= 0) && (j' < dim_j) valid_k :: Int -> Bool valid_k k' = (k' >= 0) && (k' < dim_k) -- | 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 :: 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, left = value_at values i (j-1) k, right = value_at values i (j+1) k, down = value_at values i j (k-1), top = value_at values i j (k+1), front_left = value_at values (i-1) (j-1) k, front_right = value_at values (i-1) (j+1) k, front_down =value_at values (i-1) j (k-1), front_top = value_at values (i-1) j (k+1), back_left = value_at values (i+1) (j-1) k, back_right = value_at values (i+1) (j+1) k, back_down = value_at values (i+1) j (k-1), back_top = value_at values (i+1) j (k+1), left_down = value_at values i (j-1) (k-1), left_top = value_at values i (j-1) (k+1), right_down = value_at values i (j+1) (k-1), right_top = value_at values i (j+1) (k+1), front_left_down = value_at values (i-1) (j-1) (k-1), front_left_top = value_at values (i-1) (j-1) (k+1), front_right_down = value_at values (i-1) (j+1) (k-1), front_right_top = value_at values (i-1) (j+1) (k+1), back_left_down = value_at values (i+1) (j-1) (k-1), back_left_top = value_at values (i+1) (j-1) (k+1), back_right_down = value_at values (i+1) (j+1) (k-1), back_right_top = value_at values (i+1) (j+1) (k+1), interior = value_at values i j k } -- | Takes a 'FunctionValues' and a function that transforms one -- 'Cardinal' to another (called a rotation). Then it applies the -- rotation to each element of the 'FunctionValues' object, and -- returns the result. rotate :: (Cardinal -> Cardinal) -> FunctionValues -> FunctionValues rotate rotation fv = FunctionValues { front = eval fv (rotation F), back = eval fv (rotation B), left = eval fv (rotation L), right = eval fv (rotation R), down = eval fv (rotation D), top = eval fv (rotation T), front_left = eval fv (rotation FL), front_right = eval fv (rotation FR), front_down = eval fv (rotation FD), front_top = eval fv (rotation FT), back_left = eval fv (rotation BL), back_right = eval fv (rotation BR), back_down = eval fv (rotation BD), back_top = eval fv (rotation BT), left_down = eval fv (rotation LD), left_top = eval fv (rotation LT), right_down = eval fv (rotation RD), right_top = eval fv (rotation RT), front_left_down = eval fv (rotation FLD), front_left_top = eval fv (rotation FLT), front_right_down = eval fv (rotation FRD), front_right_top = eval fv (rotation FRT), back_left_down = eval fv (rotation BLD), back_left_top = eval fv (rotation BLT), back_right_down = eval fv (rotation BRD), back_right_top = eval fv (rotation BRT), interior = interior fv } -- | Ensure that the trilinear values wind up where we think they -- should. test_directions :: Assertion test_directions = assertTrue "all direction functions work" (and equalities) where fvs = make_values trilinear 1 1 1 equalities = [ interior fvs == 4, front fvs == 1, back fvs == 7, left fvs == 2, right fvs == 6, down fvs == 3, top fvs == 5, front_left fvs == 1, front_right fvs == 1, front_down fvs == 1, front_top fvs == 1, back_left fvs == 3, back_right fvs == 11, back_down fvs == 5, back_top fvs == 9, left_down fvs == 2, left_top fvs == 2, right_down fvs == 4, right_top fvs == 8, front_left_down fvs == 1, front_left_top fvs == 1, front_right_down fvs == 1, front_right_top fvs == 1, back_left_down fvs == 3, back_left_top fvs == 3, back_right_down fvs == 7, back_right_top fvs == 15] function_values_tests :: TestTree function_values_tests = testGroup "FunctionValues tests" [ testCase "test directions" test_directions ] prop_x_rotation_doesnt_affect_front :: FunctionValues -> Bool prop_x_rotation_doesnt_affect_front fv0 = expr1 == expr2 where fv1 = rotate cwx fv0 expr1 = front fv0 expr2 = front fv1 prop_x_rotation_doesnt_affect_back :: FunctionValues -> Bool prop_x_rotation_doesnt_affect_back fv0 = expr1 == expr2 where fv1 = rotate cwx fv0 expr1 = back fv0 expr2 = back fv1 prop_y_rotation_doesnt_affect_left :: FunctionValues -> Bool prop_y_rotation_doesnt_affect_left fv0 = expr1 == expr2 where fv1 = rotate cwy fv0 expr1 = left fv0 expr2 = left fv1 prop_y_rotation_doesnt_affect_right :: FunctionValues -> Bool prop_y_rotation_doesnt_affect_right fv0 = expr1 == expr2 where fv1 = rotate cwy fv0 expr1 = right fv0 expr2 = right fv1 prop_z_rotation_doesnt_affect_down :: FunctionValues -> Bool prop_z_rotation_doesnt_affect_down fv0 = expr1 == expr2 where fv1 = rotate cwz fv0 expr1 = down fv0 expr2 = down fv1 prop_z_rotation_doesnt_affect_top :: FunctionValues -> Bool prop_z_rotation_doesnt_affect_top fv0 = expr1 == expr2 where fv1 = rotate cwz fv0 expr1 = top fv0 expr2 = top fv1 function_values_properties :: TestTree function_values_properties = let tp = testProperty in testGroup "FunctionValues properties" [ tp "x rotation doesn't affect front" prop_x_rotation_doesnt_affect_front, tp "x rotation doesn't affect back" prop_x_rotation_doesnt_affect_back, tp "y rotation doesn't affect left" prop_y_rotation_doesnt_affect_left, tp "y rotation doesn't affect right" prop_y_rotation_doesnt_affect_right, tp "z rotation doesn't affect top" prop_z_rotation_doesnt_affect_top, tp "z rotation doesn't affect down" prop_z_rotation_doesnt_affect_down ]