1 -- The "value_at" function pattern matches on some integers, but
2 -- doesn't handle the "otherwise" case, for performance reasons.
3 {-# OPTIONS_GHC -Wno-incomplete-patterns #-}
4 {-# LANGUAGE BangPatterns #-}
6 -- | The FunctionValues module contains the 'FunctionValues' type and
7 -- the functions used to manipulate it.
9 module FunctionValues (
15 function_values_tests,
16 function_values_properties,
33 import Test.Tasty ( TestTree, testGroup )
34 import Test.Tasty.HUnit ( Assertion, testCase )
35 import Test.Tasty.QuickCheck ( Arbitrary( arbitrary ), choose, testProperty )
37 import Assertions ( assertTrue )
39 Cardinal(F, B, L, R, D, T, FL, FR, FD, FT, BL, BR, BD, BT, LD, LT, RD,
40 RT, FLD, FLT, FRD, FRT, BLD, BLT, BRD, BRT, I, Scalar, Sum,
41 Difference, Product, Quotient ),
45 import Examples ( trilinear )
46 import Values ( Values3D, dims, idx )
48 -- | The FunctionValues type represents the value of our function f at
49 -- the 27 points surrounding (and including) the center of a
50 -- cube. Each value of f can be accessed by the name of its
54 FunctionValues { front :: !Double,
60 front_left :: !Double,
61 front_right :: !Double,
62 front_down :: !Double,
65 back_right :: !Double,
70 right_down :: !Double,
72 front_left_down :: !Double,
73 front_left_top :: !Double,
74 front_right_down :: !Double,
75 front_right_top :: !Double,
76 back_left_down :: !Double,
77 back_left_top :: !Double,
78 back_right_down :: !Double,
79 back_right_top :: !Double,
84 instance Arbitrary FunctionValues where
86 front' <- choose (min_double, max_double)
87 back' <- choose (min_double, max_double)
88 left' <- choose (min_double, max_double)
89 right' <- choose (min_double, max_double)
90 top' <- choose (min_double, max_double)
91 down' <- choose (min_double, max_double)
92 front_left' <- choose (min_double, max_double)
93 front_right' <- choose (min_double, max_double)
94 front_top' <- choose (min_double, max_double)
95 front_down' <- choose (min_double, max_double)
96 back_left' <- choose (min_double, max_double)
97 back_right' <- choose (min_double, max_double)
98 back_top' <- choose (min_double, max_double)
99 back_down' <- choose (min_double, max_double)
100 left_top' <- choose (min_double, max_double)
101 left_down' <- choose (min_double, max_double)
102 right_top' <- choose (min_double, max_double)
103 right_down' <- choose (min_double, max_double)
104 front_left_top' <- choose (min_double, max_double)
105 front_left_down' <- choose (min_double, max_double)
106 front_right_top' <- choose (min_double, max_double)
107 front_right_down' <- choose (min_double, max_double)
108 back_left_top' <- choose (min_double, max_double)
109 back_left_down' <- choose (min_double, max_double)
110 back_right_top' <- choose (min_double, max_double)
111 back_right_down' <- choose (min_double, max_double)
112 interior' <- choose (min_double, max_double)
114 return empty_values { front = front',
120 front_left = front_left',
121 front_right = front_right',
122 front_top = front_top',
123 front_down = front_down',
124 back_left = back_left',
125 back_right = back_right',
126 back_top = back_top',
127 back_down = back_down',
128 left_top = left_top',
129 left_down = left_down',
130 right_top = right_top',
131 right_down = right_down',
132 front_left_top = front_left_top',
133 front_left_down = front_left_down',
134 front_right_top = front_right_top',
135 front_right_down = front_right_down',
136 back_left_top = back_left_top',
137 back_left_down = back_left_down',
138 back_right_top = back_right_top',
139 back_right_down = back_right_down',
140 interior = interior' }
142 -- | We perform addition with the function values contained in a
143 -- FunctionValues object. If we choose random doubles near the machine
144 -- min/max, we risk overflowing or underflowing the 'Double'. This
145 -- places a reasonably safe limit on the maximum size of our generated
150 -- | See 'max_double'.
152 min_double = (-1) * max_double
155 -- | Return a 'FunctionValues' with a bunch of zeros for data points.
156 empty_values :: FunctionValues
158 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
161 -- | The eval function is where the magic happens for the
162 -- FunctionValues type. Given a 'Cardinal' direction and a
163 -- 'FunctionValues' object, eval will return the value of the
164 -- function f in that 'Cardinal' direction. Note that 'Cardinal' can
165 -- be a composite type; eval is what performs the \"arithmetic\" on
166 -- 'Cardinal' directions.
167 eval :: FunctionValues -> Cardinal -> Double
174 eval f FL = front_left f
175 eval f FR = front_right f
176 eval f FD = front_down f
177 eval f FT = front_top f
178 eval f BL = back_left f
179 eval f BR = back_right f
180 eval f BD = back_down f
181 eval f BT = back_top f
182 eval f LD = left_down f
183 eval f LT = left_top f
184 eval f RD = right_down f
185 eval f RT = right_top f
186 eval f FLD = front_left_down f
187 eval f FLT = front_left_top f
188 eval f FRD = front_right_down f
189 eval f FRT = front_right_top f
190 eval f BLD = back_left_down f
191 eval f BLT = back_left_top f
192 eval f BRD = back_right_down f
193 eval f BRT = back_right_top f
194 eval f I = interior f
195 eval _ (Scalar x) = x
196 eval f (Sum x y) = (eval f x) + (eval f y)
197 eval f (Difference x y) = (eval f x) - (eval f y)
198 eval f (Product x y) = (eval f x) * (eval f y)
199 eval f (Quotient x y) = (eval f x) / (eval f y)
202 -- | Takes a three-dimensional list of 'Double' and a set of 3D
203 -- coordinates (i,j,k), and returns the value at (i,j,k) in the
204 -- supplied list. If there is no such value, we calculate one
205 -- according to Sorokina and Zeilfelder, remark 7.3, p. 99.
207 -- We specifically do not consider values more than one unit away
212 -- >>> value_at Examples.trilinear 0 0 0
215 -- >>> value_at Examples.trilinear (-1) 0 0
218 -- >>> value_at Examples.trilinear 0 0 4
221 -- >>> value_at Examples.trilinear 1 3 0
224 value_at :: Values3D -> Int -> Int -> Int -> Double
225 value_at v3d !i !j !k
226 -- Put the most common case first!
227 | (valid_i i) && (valid_j j) && (valid_k k) =
230 -- The next three are from the first line in (7.3). Analogous cases
231 -- have been added where the indices are one-too-big. These are the
232 -- "one index is bad" cases.
236 -- We're one-dimensional in our first coordinate, so just
237 -- return the data point that we do have. If we try to use
238 -- the formula from remark 7.3, we go into an infinite loop.
243 2*(value_at v3d 0 j k) - (value_at v3d 1 j k)
245 2*(value_at v3d (i-1) j k) - (value_at v3d (i-2) j k)
250 -- We're one-dimensional in our second coordinate, so just
251 -- return the data point that we do have. If we try to use
252 -- the formula from remark 7.3, we go into an infinite loop.
257 2*(value_at v3d i 0 k) - (value_at v3d i 1 k)
259 2*(value_at v3d i (j-1) k) - (value_at v3d i (j-2) k)
264 -- We're one-dimensional in our third coordinate, so just
265 -- return the data point that we do have. If we try to use
266 -- the formula from remark 7.3, we go into an infinite loop.
271 2*(value_at v3d i j 0) - (value_at v3d i j 1)
273 2*(value_at v3d i j (k-1)) - (value_at v3d i j (k-2))
275 (dim_i, dim_j, dim_k) = dims v3d
277 valid_i :: Int -> Bool
278 valid_i i' = (i' >= 0) && (i' < dim_i)
280 valid_j :: Int -> Bool
281 valid_j j' = (j' >= 0) && (j' < dim_j)
283 valid_k :: Int -> Bool
284 valid_k k' = (k' >= 0) && (k' < dim_k)
288 -- | Given a three-dimensional list of 'Double' and a set of 3D
289 -- coordinates (i,j,k), constructs and returns the 'FunctionValues'
290 -- object centered at (i,j,k)
291 make_values :: Values3D -> Int -> Int -> Int -> FunctionValues
292 make_values values !i !j !k =
293 empty_values { front = value_at values (i - 1) j k,
294 back = value_at values (i + 1) j k,
295 left = value_at values i (j - 1) k,
296 right = value_at values i (j + 1) k,
297 down = value_at values i j (k - 1),
298 top = value_at values i j (k + 1),
299 front_left = value_at values (i - 1) (j - 1) k,
300 front_right = value_at values (i - 1) (j + 1) k,
301 front_down =value_at values (i - 1) j (k - 1),
302 front_top = value_at values (i - 1) j (k + 1),
303 back_left = value_at values (i + 1) (j - 1) k,
304 back_right = value_at values (i + 1) (j + 1) k,
305 back_down = value_at values (i + 1) j (k - 1),
306 back_top = value_at values (i + 1) j (k + 1),
307 left_down = value_at values i (j - 1) (k - 1),
308 left_top = value_at values i (j - 1) (k + 1),
309 right_down = value_at values i (j + 1) (k - 1),
310 right_top = value_at values i (j + 1) (k + 1),
311 front_left_down = value_at values (i - 1) (j - 1) (k - 1),
312 front_left_top = value_at values (i - 1) (j - 1) (k + 1),
313 front_right_down = value_at values (i - 1) (j + 1) (k - 1),
314 front_right_top = value_at values (i - 1) (j + 1) (k + 1),
315 back_left_down = value_at values (i + 1) (j - 1) (k - 1),
316 back_left_top = value_at values (i + 1) (j - 1) (k + 1),
317 back_right_down = value_at values (i + 1) (j + 1) (k - 1),
318 back_right_top = value_at values (i + 1) (j + 1) (k + 1),
319 interior = value_at values i j k }
321 -- | Takes a 'FunctionValues' and a function that transforms one
322 -- 'Cardinal' to another (called a rotation). Then it applies the
323 -- rotation to each element of the 'FunctionValues' object, and
324 -- returns the result.
325 rotate :: (Cardinal -> Cardinal) -> FunctionValues -> FunctionValues
327 FunctionValues { front = eval fv (rotation F),
328 back = eval fv (rotation B),
329 left = eval fv (rotation L),
330 right = eval fv (rotation R),
331 down = eval fv (rotation D),
332 top = eval fv (rotation T),
333 front_left = eval fv (rotation FL),
334 front_right = eval fv (rotation FR),
335 front_down = eval fv (rotation FD),
336 front_top = eval fv (rotation FT),
337 back_left = eval fv (rotation BL),
338 back_right = eval fv (rotation BR),
339 back_down = eval fv (rotation BD),
340 back_top = eval fv (rotation BT),
341 left_down = eval fv (rotation LD),
342 left_top = eval fv (rotation LT),
343 right_down = eval fv (rotation RD),
344 right_top = eval fv (rotation RT),
345 front_left_down = eval fv (rotation FLD),
346 front_left_top = eval fv (rotation FLT),
347 front_right_down = eval fv (rotation FRD),
348 front_right_top = eval fv (rotation FRT),
349 back_left_down = eval fv (rotation BLD),
350 back_left_top = eval fv (rotation BLT),
351 back_right_down = eval fv (rotation BRD),
352 back_right_top = eval fv (rotation BRT),
353 interior = interior fv }
357 -- | Ensure that the trilinear values wind up where we think they
359 test_directions :: Assertion
361 assertTrue "all direction functions work" (and equalities)
363 fvs = make_values trilinear 1 1 1
364 equalities = [ interior fvs == 4,
372 front_right fvs == 1,
376 back_right fvs == 11,
383 front_left_down fvs == 1,
384 front_left_top fvs == 1,
385 front_right_down fvs == 1,
386 front_right_top fvs == 1,
387 back_left_down fvs == 3,
388 back_left_top fvs == 3,
389 back_right_down fvs == 7,
390 back_right_top fvs == 15]
393 function_values_tests :: TestTree
394 function_values_tests =
395 testGroup "FunctionValues tests"
396 [ testCase "test directions" test_directions ]
399 prop_x_rotation_doesnt_affect_front :: FunctionValues -> Bool
400 prop_x_rotation_doesnt_affect_front fv0 =
407 prop_x_rotation_doesnt_affect_back :: FunctionValues -> Bool
408 prop_x_rotation_doesnt_affect_back fv0 =
416 prop_y_rotation_doesnt_affect_left :: FunctionValues -> Bool
417 prop_y_rotation_doesnt_affect_left fv0 =
424 prop_y_rotation_doesnt_affect_right :: FunctionValues -> Bool
425 prop_y_rotation_doesnt_affect_right fv0 =
433 prop_z_rotation_doesnt_affect_down :: FunctionValues -> Bool
434 prop_z_rotation_doesnt_affect_down fv0 =
442 prop_z_rotation_doesnt_affect_top :: FunctionValues -> Bool
443 prop_z_rotation_doesnt_affect_top fv0 =
451 function_values_properties :: TestTree
452 function_values_properties =
453 testGroup "FunctionValues properties" [
455 "x rotation doesn't affect front"
456 prop_x_rotation_doesnt_affect_front,
458 "x rotation doesn't affect back"
459 prop_x_rotation_doesnt_affect_back,
461 "y rotation doesn't affect left"
462 prop_y_rotation_doesnt_affect_left,
464 "y rotation doesn't affect right"
465 prop_y_rotation_doesnt_affect_right,
467 "z rotation doesn't affect top"
468 prop_z_rotation_doesnt_affect_top,
470 "z rotation doesn't affect down"
471 prop_z_rotation_doesnt_affect_down ]