Difference between revisions of "Haskell Quiz/Geodesic Dome Faces/Solution Jkramar"

From HaskellWiki
Jump to navigation Jump to search
m
m
Line 5: Line 5:
   
 
<haskell>
 
<haskell>
import Prelude hiding (sum, foldr1, foldr, maximum, sequence_)
+
import Prelude hiding (sum, foldr1, foldr, maximum, mapM_)
 
import Control.Applicative
 
import Control.Applicative
 
import Data.Foldable
 
import Data.Foldable
Line 38: Line 38:
 
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 = fmap$foldr1 (.+).(<*>t).fmap (.*)
 
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
Line 46: Line 43:
 
chooses = foldr consider$[[]]:repeat [] where
 
chooses = foldr consider$[[]]:repeat [] where
 
consider x cs = zipWith (flip (++)) cs$map (x:)<$>[]:cs
 
consider x cs = zipWith (flip (++)) cs$map (x:)<$>[]:cs
 
tris :: [V3 a] -> [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 (signum$det t)==1 then t else V a c b
 
   
 
-- inefficient function to generate all the positively-oriented faces of the
 
-- inefficient function to generate all the positively-oriented faces of the
Line 57: Line 48:
 
faces :: (Real a) => [V3 a] -> [Tri a]
 
faces :: (Real a) => [V3 a] -> [Tri a]
 
faces vs = filter isFace $ orient <$> tris vs where
 
faces vs = filter isFace $ orient <$> tris vs where
 
orient t@(V a b c) = if (signum$det t)==1 then t else V a c b
 
tris xs = (\[a,b,c]->V a b c) <$> (chooses xs!!3)
 
farth (V a b c) = maximum.(dot ((a.-b) `cross` (a.-c))<$>)
 
farth (V a b c) = maximum.(dot ((a.-b) `cross` (a.-c))<$>)
 
isFace t@(V a b c) = farth t [a,b,c] == farth t vs
 
isFace t@(V a b c) = farth t [a,b,c] == farth t vs
Line 62: Line 55:
 
-- 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) => a -> Tri b -> [Tri b]
+
shatter :: (Fractional a) => Int -> Tri a -> [Tri a]
 
shatter n t = 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)
Line 70: Line 63:
 
coords = (fmap fromIntegral<$>)<$>fwd++bwd
 
coords = (fmap fromIntegral<$>)<$>fwd++bwd
   
geode :: (Integral a, RealFloat b) => [V3 b] -> a -> [Tri b]
+
geode :: (RealFloat a) => [V3 a] -> Int -> [Tri a]
geode vs n = fmap normize<$>(shatter n=<<faces vs)
+
geode vs n = fmap normize<$>(shatter n=<<faces vs) where
 
normize v = (1/sqrt (dot v v)).*v
   
 
cyc :: V3 a -> [V3 a]
 
cyc :: V3 a -> [V3 a]
Line 87: Line 81:
   
 
main :: IO ()
 
main :: IO ()
main = mapM_ print (geode octahed (51::Int)::[Tri Double])
+
main = mapM_ print (geode octahed 51::[Tri Double])
 
</haskell>
 
</haskell>

Revision as of 20:05, 17 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 vs where
  orient t@(V a b c) = if (signum$det t)==1 then t else V a c b
  tris xs = (\[a,b,c]->V a b c) <$> (chooses xs!!3)
  farth (V a b c) = maximum.(dot ((a.-b) `cross` (a.-c))<$>)
  isFace t@(V a b c) = farth t [a,b,c] == 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]]

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