[Haskell-cafe] fast Eucl. dist. - Haskell vs C

[Kenneth Hoste]

On May 19, 2009, at 13:24 , Daniel Schüssler wrote:

> Hello!
>
> On Monday 18 May 2009 14:37:51 Kenneth Hoste wrote:
>> I'm mostly interested in the range 10D to 100D
>
> is the dimension known at compile-time? Then you could consider  
> Template
> Haskell.

In general, no. :-)

It will be known for some applications, but not for others.

I'm more and more amazed what comes into play just to implement
something simple like n-dim. Euclidean distance relatively fast using  
Haskell.

It seems to me that GHC is missing several critical optimizations
(yes, I know, patches welcome) to enable it to emit fast code
for HPC applications.

I'm still a big fan of Haskell, for a variety of reasons, but it seems
like it's not ready yet for the task I had in mind, which is a shame.

Just to be clear, this isn't a flame bait post or anything, just my 2  
cents.

K.

> I wrote up some code for generating the vector types and vector
> subtraction/inner product below, HTH. One problem is that I'm using a
> typeclass and apparently you can't make {-# SPECIALISE #-} pragmas  
> with TH,
> so let's hope it is automatically specialised by GHC.
>
> Greetings,
> Daniel
>
> TH.hs
> ------------------
>
>
> {-# LANGUAGE TemplateHaskell #-}
> {-# OPTIONS -fglasgow-exts #-}
>
> module TH where
>
> import Language.Haskell.TH
> import Control.Monad
>
> -- Non-TH stuff
> class InnerProductSpace v r | v -> r where
>    innerProduct :: v -> v -> r
>
> class AbGroup v where
>    minus :: v -> v -> v
>
> euclidean x y = case minus x y of
>                  z -> sqrt $! innerProduct z z
>
> -- TH
> noContext :: Q Cxt
> noContext = return []
> strict :: Q Type -> StrictTypeQ
> strict = liftM ((,) IsStrict)
>
> makeVectors :: Int -- ^ Dimension
>            -> Q Type -- ^ Component type, assumed to be a 'Num'
>            -> String -- ^ Name for the generated type
>            -> Q [Dec]
> makeVectors n ctyp name0 = do
>  -- let's assume ctyp = Double, name = Vector for the comments
>
>  -- generate names for the variables we will need
>  xs <- replicateM n (newName "x")
>  ys <- replicateM n (newName "y")
>
>  let
>      name = mkName name0
>
>      -- shorthands for arithmetic expressions; the first takes  
> expressions,
>      -- the others take variable names
>      sumE  e1 e2     = infixE (Just       e1)  [|(+)|] (Just       e2)
>      varDiffE e1 e2  = infixE (Just (varE e1)) [|(-)|] (Just (varE  
> e2))
>      varProdE e1 e2  = infixE (Just (varE e1)) [|(*)|] (Just (varE  
> e2))
>
>
>      conPat vars = conP name (fmap varP vars)
>
>      -- > data Vector = Vector !Double ... !Double
>      theDataD =
>          dataD noContext name [] -- no context, no params
>                    [normalC name (replicate n (strict ctyp))]
>                    [''Eq,''Ord,''Show] -- 'deriving' clause
>
>      innerProdD =
>          -- > instance InnerProductSpace Vector Double where ...
>          instanceD noContext ( conT ''InnerProductSpace
>                                `appT` conT name
>                                `appT` ctyp)
>                    -- > innerProduct = ...
>                    [valD
>                     (varP 'innerProduct)
>                     (normalB
>                      -- \(Vector x1 x2 ... xn) (Vector y1 y2 ... yn)  
> ->
>                      (lamE [conPat xs, conPat ys]
>                       -- x1*y1 + .... + xn*yn + 0
>                       (foldl sumE [|0|] $
>                              zipWith varProdE xs ys)
>                       ))
>
>                     [] -- no 'where' clause
>                    ]
>
>      abGroupD =
>          instanceD noContext ( conT ''AbGroup
>                                `appT` conT name)
>                    -- > minus = ...
>                    [valD
>                     (varP 'minus)
>                     (normalB
>                      -- \(Vector x1 x2 ... xn) (Vector y1 y2 ... yn)  
> ->
>                      (lamE [conPat xs, conPat ys]
>                       -- Vector (x1-y1) ... (xn-yn)
>                       (foldl appE (conE name) $
>                              zipWith varDiffE xs ys)
>                       ))
>
>                     [] -- no 'where' clause
>                    ]
>
>
>  sequence [theDataD,innerProdD,abGroupD]
>
>
>
> Main.hs
> ------------------
> {-# LANGUAGE TemplateHaskell #-}
> {-# LANGUAGE MultiParamTypeClasses #-}
>
> module Main where
>
> import TH
>
> $(makeVectors 3 [t|Double|] "Vec3")
>
> main = print $ euclidean (Vec3 1 1 1) (Vec3 0 0 0)
> _______________________________________________
> Haskell-Cafe mailing list
> Haskell-Cafe at haskell.org
> http://www.haskell.org/mailman/listinfo/haskell-cafe

-- 

Kenneth Hoste
Paris research group - ELIS - Ghent University, Belgium
email: kenneth.hoste at elis.ugent.be
website: http://www.elis.ugent.be/~kehoste
blog: http://boegel.kejo.be