Personal tools

Haskell Quiz/Geodesic Dome Faces/Solution Jkramar

From HaskellWiki

< Haskell Quiz | Geodesic Dome Faces(Difference between revisions)
Jump to: navigation, search
m
Line 80: Line 80:
 
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 used
 
main :: IO ()
 
main :: IO ()
 
main = mapM_ print (geode octahed 51::[Tri Double])
 
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 :: [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 :: [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).take 3.tails.cycle
  +
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 = toColumns.(fromColumns t*/n<>).fromColumns<$>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*/pnorm PNorm2 v
  +
  +
cyc :: [a] -> [[a]]
  +
cyc = map (take 3).take 3.tails.cycle
  +
  +
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
 
</haskell>
 
</haskell>

Revision as of 07:00, 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 :: [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 :: [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).take 3.tails.cycle
  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 = toColumns.(fromColumns t*/n<>).fromColumns<$>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*/pnorm PNorm2 v
 
cyc :: [a] -> [[a]]
cyc = map (take 3).take 3.tails.cycle
 
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