# Haskell Quiz/Geodesic Dome Faces/Solution Jkramar

### From HaskellWiki

< Haskell Quiz | Geodesic Dome Faces(Difference between revisions)

m |
m |
||

(15 intermediate revisions by one user not shown) | |||

Line 1: | Line 1: | ||

− | This problem seems to be strongly IO-bound, so actually computing the geodesic faces is not too time-sensitive. Hence there is time for computing the faces from the vertices by trying each triple of vertices to check if it's a face. |
+ | [[Category:Haskell Quiz solutions|Geodesic Dome Faces]] |

+ | This problem seems to be strongly IO-bound, so actually computing the geodesic faces is not too time-sensitive. Hence there is time for computing the faces from the vertices by trying each triple of vertices to check if it's a face. (In the problem statement, the data for each polyhedron includes a list of the faces, not just the vertices.) |
||

− | Originally I wrote this program representing vectors just by ordinary tuples, and in some places that was definitely more understandable. But of course it's more fun this way. |
+ | Originally this program represented vectors with ordinary tuples, and in some places that was definitely more understandable. But of course it's more fun this way. |

<haskell> |
<haskell> |
||

− | import Prelude hiding (maximum, foldr, foldr1, concat, sequence_, sum) |
+ | import Prelude hiding (sum, foldr1, foldr, maximum, mapM_) |

− | import Control.Monad hiding (sequence_) |
||

import Control.Applicative |
import Control.Applicative |
||

import Data.Foldable |
import Data.Foldable |
||

Line 26: | Line 26: | ||

dot :: (Num a) => V3 a -> V3 a -> a |
dot :: (Num a) => V3 a -> V3 a -> a |
||

− | dot a b = sum$liftA2 (*) a b |
+ | dot = (sum.).liftA2 (*) |

cross :: (Num a) => V3 a -> V3 a -> V3 a |
cross :: (Num a) => V3 a -> V3 a -> V3 a |
||

Line 36: | Line 36: | ||

-- matrix multiplication |
-- matrix multiplication |
||

mmul :: (Num a) => V3 (V3 a) -> V3 (V3 a) -> V3 (V3 a) |
mmul :: (Num a) => V3 (V3 a) -> V3 (V3 a) -> V3 (V3 a) |
||

− | mmul t s = foldr1 (.+) <$> (flip (.*)<$>t<*>) <$> s |
+ | mmul t = fmap$foldr1 (.+).(<*>t).fmap (.*) |

− | |||

− | normize :: (Floating a) => V3 a -> V3 a |
||

− | normize v = (1/sqrt (dot v v)).*v |
||

-- chooses xs!!n gives the combinations of n elements from xs |
-- chooses xs!!n gives the combinations of n elements from xs |
||

− | chooses :: (Foldable t, Alternative f) => t a -> [f [a]] |
+ | chooses :: [a] -> [[[a]]] |

− | chooses = foldr consider$pure []:repeat empty where |
+ | chooses = foldr consider$[[]]:repeat [] where |

− | consider x cs = zipWith (flip (<|>)) cs$fmap (x:)<$>empty:cs |
+ | consider x cs = zipWith (flip (++)) cs$map (x:)<$>[]:cs |

− | |||

− | tris :: (Foldable t, Alternative f) => t (V3 a) -> f (Tri a) |
||

− | tris xs = (\[a,b,c]->V a b c) <$> (chooses xs!!3) |
||

− | |||

− | orient :: (Num a) => Tri a -> Tri a |
||

− | orient t@(V a b c) = if (==1)$signum$det t then t else V a c b |
||

− | |||

− | farth :: (Real a, Foldable t) => V3 a -> t (V3 a) -> a |
||

− | farth dir = maximum.map (dot dir).toList |
||

-- inefficient function to generate all the positively-oriented faces of the |
-- inefficient function to generate all the positively-oriented faces of the |
||

-- triangle-faced polyhedron with vertices at vs |
-- triangle-faced polyhedron with vertices at vs |
||

− | faces :: (Real a, Foldable t, MonadPlus m) => t (V3 a) -> m (Tri a) |
+ | faces :: (Real a) => [V3 a] -> [Tri a] |

− | faces vs = do |
+ | faces vs = filter isFace $ orient <$> tris where |

− | t@(V a b c) <- unwrapMonad $ orient <$> tris vs |
+ | orient t@(V a b c) = if (signum$det t)==1 then t else V a c b |

− | let dir = ((a.-b) `cross` (a.-c)) |
+ | tris = (\[a,b,c]->V a b c) <$> chooses vs!!3 |

− | guard (farth dir [a,b,c]==farth dir vs) >> return t |
+ | farth (V a b c) = maximum.(dot ((a.-b) `cross` (a.-c))<$>) |

+ | isFace t = farth t (toList t) == farth t vs |
||

-- each triangle side is broken into n pieces, unlike in the problem statement, |
-- each triangle side is broken into n pieces, unlike in the problem statement, |
||

-- where they use n+1 for some reason |
-- where they use n+1 for some reason |
||

− | shatter :: (Integral a, Fractional b, MonadPlus m) => a -> Tri b -> m (Tri b) |
+ | shatter :: (Fractional a) => Int -> Tri a -> [Tri a] |

− | shatter n t = msum$return<$>mmul t'<$>coords where |
+ | shatter n t = mmul t'<$>coords where |

basis = V (V 1 0 0) (V 0 1 0) (V 0 0 1) |
basis = V (V 1 0 0) (V 0 1 0) (V 0 0 1) |
||

fwd = [(V k j (n-1-k-j).+) <$> basis|k<-[0..n-1], j<-[0..n-1-k]] |
fwd = [(V k j (n-1-k-j).+) <$> basis|k<-[0..n-1], j<-[0..n-1-k]] |
||

Line 61: | Line 61: | ||

coords = (fmap fromIntegral<$>)<$>fwd++bwd |
coords = (fmap fromIntegral<$>)<$>fwd++bwd |
||

− | geode :: (RealFloat b, Integral a, MonadPlus m, Foldable t) => |
+ | geode :: (RealFloat a) => [V3 a] -> Int -> [Tri a] |

− | t (V3 b) -> a -> m (Tri b) |
+ | geode vs n = fmap normize<$>(shatter n=<<faces vs) where |

− | geode vs n = return.fmap normize=<<shatter n=<<faces vs |
+ | normize v = (1/sqrt (dot v v)).*v |

cyc :: V3 a -> [V3 a] |
cyc :: V3 a -> [V3 a] |
||

Line 78: | Line 78: | ||

icosahed = cyc =<< [V x (y*(1+sqrt 5/2)) 0|x<-[-1, 1], y<-[-1, 1]] |
icosahed = cyc =<< [V x (y*(1+sqrt 5/2)) 0|x<-[-1, 1], y<-[-1, 1]] |
||

− | -- this is the test the Ruby Quiz people were doing to see how fast the code is |
+ | -- this is the test the Ruby Quiz people used |

+ | main :: IO () |
||

+ | main = mapM_ print (geode octahed 51::[Tri Double]) |
||

+ | </haskell> |
||

+ | |||

+ | Here is a reimplementation using the vector and matrix types from the hmatrix package: |
||

+ | |||

+ | <haskell> |
||

+ | import Numeric.LinearAlgebra |
||

+ | import Data.List |
||

+ | import Control.Applicative ((<$>)) |
||

+ | |||

+ | cross :: Vector Double -> Vector Double -> Vector Double |
||

+ | cross v w = fromList$map (det.fromColumns.(:[v,w]))$toColumns$ident 3 |
||

+ | |||

+ | -- chooses xs!!n gives the combinations of n elements from xs |
||

+ | chooses :: [a] -> [[[a]]] |
||

+ | chooses = foldr consider$[[]]:repeat [] where |
||

+ | consider x cs = zipWith (flip (++)) cs$map (x:)<$>[]:cs |
||

+ | |||

+ | cyc :: [a] -> [[a]] |
||

+ | cyc = map (take 3).take 3.tails.cycle |
||

+ | |||

+ | -- inefficient function to generate all the positively-oriented faces of the |
||

+ | -- triangle-faced polyhedron with vertices at vs |
||

+ | faces :: [Vector Double] -> [[Vector Double]] |
||

+ | faces vs = filter isFace $ orient <$> chooses vs!!3 where |
||

+ | orient t@[a,b,c] = if (signum$det$fromColumns t)==1 then t else [a,c,b] |
||

+ | farth = (maximum.).map.dot.sum.map (\(v:w:_)->cross v w).cyc |
||

+ | isFace ps = farth ps ps == farth ps vs |
||

+ | |||

+ | -- each triangle side is broken into n pieces, unlike in the problem statement, |
||

+ | -- where they use n+1 for some reason |
||

+ | shatter :: Int -> [Vector Double] -> [[Vector Double]] |
||

+ | shatter n' t = map (fromColumns t*/n<>)<$>fwd++bwd where |
||

+ | n = fromIntegral n'; basis = toColumns$ident 3 |
||

+ | fwd = [(fromList [k, j, n-1-k-j]+) <$> basis|k<-[0..n-1], j<-[0..n-1-k]] |
||

+ | bwd = [(fromList [k, j, n+1-k-j]-) <$> basis|k<-[1..n], j<-[1..n-k]] |
||

+ | |||

+ | geode :: [Vector Double] -> Int -> [[Vector Double]] |
||

+ | geode vs n = map proj<$>(shatter n=<<faces vs) where proj v = v*/sqrt (v<.>v) |
||

+ | |||

+ | tetrahed :: [Vector Double] |
||

+ | tetrahed = map fromList$[[x*sqrt 1.5, -sqrt 2/3, -1/3]|x<-[-1, 1]]++ |
||

+ | [[0, 2*sqrt 2/3, -1/3], [0, 0, 1]] |
||

+ | |||

+ | octahed :: [Vector Double] |
||

+ | octahed = map fromList$cyc =<< [[x, 0, 0]|x<-[-1, 1]] |
||

+ | |||

+ | icosahed :: [Vector Double] |
||

+ | icosahed = map fromList$cyc =<< [[x, y*(1+sqrt 5/2), 0]|x<-[-1, 1], y<-[-1, 1]] |
||

+ | |||

+ | -- this is the test the Ruby Quiz people used |
||

main :: IO () |
main :: IO () |
||

− | main = sequence_$map print$(geode octahed::Int->[Tri Double]) 51 |
+ | main = mapM_ print$geode octahed 51 |

</haskell> |
</haskell> |

## Latest revision as of 15:33, 18 November 2008

This problem seems to be strongly IO-bound, so actually computing the geodesic faces is not too time-sensitive. Hence there is time for computing the faces from the vertices by trying each triple of vertices to check if it's a face. (In the problem statement, the data for each polyhedron includes a list of the faces, not just the vertices.)

Originally this program represented vectors with ordinary tuples, and in some places that was definitely more understandable. But of course it's more fun this way.

import Prelude hiding (sum, foldr1, foldr, maximum, mapM_) import Control.Applicative import Data.Foldable data V3 a = V !a !a !a deriving (Read,Show) instance Foldable V3 where foldr f x (V a b c) = f a$f b$f c x instance Functor V3 where fmap f (V a b c) = V (f a) (f b) (f c) instance Applicative V3 where pure a = V a a a V f g h <*> V a b c = V (f a) (g b) (h c) type Tri a = V3 (V3 a) (.+) :: (Num a) => V3 a -> V3 a -> V3 a (.+) = liftA2 (+) (.-) :: (Num a) => V3 a -> V3 a -> V3 a (.-) = liftA2 (-) (.*) :: (Num a) => a -> V3 a -> V3 a (.*) = fmap.(*) dot :: (Num a) => V3 a -> V3 a -> a dot = (sum.).liftA2 (*) cross :: (Num a) => V3 a -> V3 a -> V3 a cross (V x y z) (V x' y' z')=V (y*z'-z*y') (z*x'-x*z') (x*y'-y*x') det :: (Num a) => V3 (V3 a) -> a det (V a b c) = (a `cross` b) `dot` c -- matrix multiplication mmul :: (Num a) => V3 (V3 a) -> V3 (V3 a) -> V3 (V3 a) mmul t = fmap$foldr1 (.+).(<*>t).fmap (.*) -- chooses xs!!n gives the combinations of n elements from xs chooses :: [a] -> [[[a]]] chooses = foldr consider$[[]]:repeat [] where consider x cs = zipWith (flip (++)) cs$map (x:)<$>[]:cs -- inefficient function to generate all the positively-oriented faces of the -- triangle-faced polyhedron with vertices at vs faces :: (Real a) => [V3 a] -> [Tri a] faces vs = filter isFace $ orient <$> tris where orient t@(V a b c) = if (signum$det t)==1 then t else V a c b tris = (\[a,b,c]->V a b c) <$> chooses vs!!3 farth (V a b c) = maximum.(dot ((a.-b) `cross` (a.-c))<$>) isFace t = farth t (toList t) == farth t vs -- each triangle side is broken into n pieces, unlike in the problem statement, -- where they use n+1 for some reason shatter :: (Fractional a) => Int -> Tri a -> [Tri a] shatter n t = mmul t'<$>coords where basis = V (V 1 0 0) (V 0 1 0) (V 0 0 1) fwd = [(V k j (n-1-k-j).+) <$> basis|k<-[0..n-1], j<-[0..n-1-k]] bwd = [(V k j (n+1-k-j).-) <$> basis|k<-[1..n], j<-[1..n-k]] t' = ((1/fromIntegral n).*)<$>t coords = (fmap fromIntegral<$>)<$>fwd++bwd geode :: (RealFloat a) => [V3 a] -> Int -> [Tri a] geode vs n = fmap normize<$>(shatter n=<<faces vs) where normize v = (1/sqrt (dot v v)).*v cyc :: V3 a -> [V3 a] cyc (V a b c) = [V a b c, V b c a, V c a b] tetrahed :: (Floating a) => [V3 a] tetrahed = [V (x*sqrt 1.5) (-sqrt 2/3) (-1/3)|x<-[-1,1]]++ [V 0 (2*sqrt 2/3) (-1/3), V 0 0 1] octahed :: (Num a) => [V3 a] octahed = cyc =<< [V x 0 0|x<-[-1,1]] icosahed :: (Floating a) => [V3 a] icosahed = cyc =<< [V x (y*(1+sqrt 5/2)) 0|x<-[-1, 1], y<-[-1, 1]] -- this is the test the Ruby Quiz people used main :: IO () main = mapM_ print (geode octahed 51::[Tri Double])

Here is a reimplementation using the vector and matrix types from the hmatrix package:

import Numeric.LinearAlgebra import Data.List import Control.Applicative ((<$>)) cross :: Vector Double -> Vector Double -> Vector Double cross v w = fromList$map (det.fromColumns.(:[v,w]))$toColumns$ident 3 -- chooses xs!!n gives the combinations of n elements from xs chooses :: [a] -> [[[a]]] chooses = foldr consider$[[]]:repeat [] where consider x cs = zipWith (flip (++)) cs$map (x:)<$>[]:cs cyc :: [a] -> [[a]] cyc = map (take 3).take 3.tails.cycle -- inefficient function to generate all the positively-oriented faces of the -- triangle-faced polyhedron with vertices at vs faces :: [Vector Double] -> [[Vector Double]] faces vs = filter isFace $ orient <$> chooses vs!!3 where orient t@[a,b,c] = if (signum$det$fromColumns t)==1 then t else [a,c,b] farth = (maximum.).map.dot.sum.map (\(v:w:_)->cross v w).cyc isFace ps = farth ps ps == farth ps vs -- each triangle side is broken into n pieces, unlike in the problem statement, -- where they use n+1 for some reason shatter :: Int -> [Vector Double] -> [[Vector Double]] shatter n' t = map (fromColumns t*/n<>)<$>fwd++bwd where n = fromIntegral n'; basis = toColumns$ident 3 fwd = [(fromList [k, j, n-1-k-j]+) <$> basis|k<-[0..n-1], j<-[0..n-1-k]] bwd = [(fromList [k, j, n+1-k-j]-) <$> basis|k<-[1..n], j<-[1..n-k]] geode :: [Vector Double] -> Int -> [[Vector Double]] geode vs n = map proj<$>(shatter n=<<faces vs) where proj v = v*/sqrt (v<.>v) tetrahed :: [Vector Double] tetrahed = map fromList$[[x*sqrt 1.5, -sqrt 2/3, -1/3]|x<-[-1, 1]]++ [[0, 2*sqrt 2/3, -1/3], [0, 0, 1]] octahed :: [Vector Double] octahed = map fromList$cyc =<< [[x, 0, 0]|x<-[-1, 1]] icosahed :: [Vector Double] icosahed = map fromList$cyc =<< [[x, y*(1+sqrt 5/2), 0]|x<-[-1, 1], y<-[-1, 1]] -- this is the test the Ruby Quiz people used main :: IO () main = mapM_ print$geode octahed 51