[Haskell-cafe] performance question http://www.haskell.org/pipermail/haskell-cafe/2013-February/106269.html

.gitignore

*.o
*.hi
*.hp
*.aux
*.prof
/printMatrixDecay.ps
/printMatrixDecay
/dist
/cabal-dev

createMatrixDump.py

#!/usr/bin/env python

import argparse
import random

parser = argparse.ArgumentParser()

parser.add_argument("-N",
    help = "Produce a matrix of size NxN",
    default = 1)

options = parser.parse_args()

N = int(options.N)

A = [ [ 0 for i in range(N) ] for j in range(N) ]

for i in range(N):
  for j in range(N):
    A[i][j] = random.random()

print("N = %i" % (N))
for i in range(N):
  for j in range(N):
    print("matrix(%i,%i) = %e" % (i, j, A[i][j]))

LICENSE

Portions by Nicolas Bock are of unknown license.

The revisions made by Bob Ippolito are in the public domain.

printMatrixDecay.cabal

-- Initial printMatrixDecay.cabal generated by cabal init.  For further 
-- documentation, see http://haskell.org/cabal/users-guide/

name:                printMatrixDecay
version:             0.1.0.0
-- synopsis:            
-- description:         
homepage:            https://gist.github.com/etrepum/4747507
-- license:             
license-file:        LICENSE
author:              Nicolas Bock and Bob Ippolito
maintainer:          bob@redivi.com
-- copyright:           
category:            Data
build-type:          Simple
cabal-version:       >=1.8

executable printMatrixDecay
  main-is:             printMatrixDecay.hs
  ghc-options:         -O2
                       -rtsopts
  ghc-prof-options:    -prof
                       -fprof-auto
  -- other-modules:       
  build-depends:       base ==4.5.*,
                       attoparsec,
                       bytestring,
                       vector,
                       regex-posix
  

printMatrixDecay.hs

{-# OPTIONS -funbox-strict-fields #-}
{-# LANGUAGE ForeignFunctionInterface #-}
{-# LANGUAGE OverloadedStrings #-}
import Text.Printf (printf)
import Data.Either (rights)
import qualified Data.ByteString.Char8 as B
import qualified Data.Vector.Unboxed as V
import Data.Attoparsec (Parser, parseOnly, skipMany, endOfInput, string)
import Data.Attoparsec.ByteString.Char8 (double, skipSpace, notChar)
import Control.Applicative ((*>), (<*))

data Strata = Strata
              { lowerBound :: !Double
              , upperBound :: !Double
              , count      :: !Int
              }
            deriving (Show, Eq)

strataBounds :: [Double]
strataBounds = [ 0.0, 1.0e-8, 1.0e-7, 1.0e-6, 1.0e-5, 1.0e-4
               , 1.0e-3, 1.0e-2, 1.0e-1, 1.0, 2.0 ]

strataVec :: V.Vector Double
strataVec = V.fromList strataBounds

-- Binary search over strataVec for the index
strataIndex :: Double -> Int
strataIndex d = go 0 lastIndex
  where vec = strataVec
        lastIndex = V.length vec - 1
        go lo hi
          | hi >= lo = case compare d pivot of
            EQ -> if mid == lastIndex
                  then outOfBounds
                  else mid
            LT -> go lo (mid - 1)
            GT -> go (mid + 1) hi
          | otherwise = if hi < 0 then outOfBounds else hi
          where mid = lo + (hi - lo) `div` 2
                pivot = V.unsafeIndex vec mid
                outOfBounds = error $ "can not stratify aij = " ++ show d

main :: IO ()
main = do
  l <- B.getContents
  printStrataCounts $ stratify $ rights $ map readMatrixLine $ B.lines l

printStrataCounts :: ([Strata], Int) -> IO ()
printStrataCounts (ss0, n) = go ss0 0
  where
    go :: [Strata] -> Int -> IO ()
    go (Strata lb ub c:ss) total = do
      let total' = total + c
      _ <- printf "[%1.2e, %1.2e) = %i (%1.2f%%) %i\n"
           lb ub c (100.0 * fromIntegral c / fromIntegral n :: Double) total'
      go ss total'
    go _ss _total = return ()

readMatrixLine :: B.ByteString -> Either String Int
readMatrixLine = fmap (strataIndex . abs) . parseOnly matrixParser

matrixParser :: Parser Double
matrixParser = skipSpace *>
               string "matrix(" *>
               skipMany (notChar '=') *>
               string "= " *>
               double <* endOfInput

{-
regex = makeRegex $ B.pack "^\\s*matrix\\(.*= ([0-9eE.+-]+)$"
-}

stratify :: [Int] -> ([Strata], Int)
stratify = total . V.accum (+) acc0 . flip zip (repeat 1)
  where
    acc0 = V.replicate (V.length strataVec - 1) 0
    wrap = map toStrata . zip3 strataBounds (tail strataBounds) . V.toList
    total v = (wrap v, V.foldl1' (+) v)
    toStrata (lb, ub, c) = Strata lb ub c

printMatrixDecay.py

#!/usr/bin/env python
#
# vim: tw=0

import argparse
import math
import re
import sys

strataBounds = [ 0.0, 1.0e-8, 1.0e-7, 1.0e-6, 1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2, 1.0e-1, 1.0, 2.0 ]

parser = argparse.ArgumentParser()

parser.add_argument("FILE",
    help = "Read metric from FILE",
    nargs = "+")

options = parser.parse_args()

for f in options.FILE:
  if f == "-":
    fd = sys.stdin
  else:
    fd = open(f, "r")

  N = 0
  strataCounts = [ 0 ] * (len(strataBounds)-1)
  for line in fd:
    result = re.compile("matrix.*= ([0-9.eE+-]+)$").search(line)
    if result:
      N += 1
      aij = math.fabs(float(result.group(1)))
      foundMatrixElement = False
      for i in range(len(strataBounds)-1):
        if (aij >= strataBounds[i]) and (aij < strataBounds[i+1]):
          strataCounts[i] += 1
          foundMatrixElement = True
          break
      if not foundMatrixElement:
        print("could not place \"%s\" anywhere" % (line.strip()))
        sys.exit(1)

  fd.close()

  print("(%s) read %i matrix elements (%ix%i = %i)" % (f, N,
    int(math.sqrt(N)), int(math.sqrt(N)), int(math.sqrt(N))**2))

  total = 0
  for i in range(len(strataBounds)-1):
    total += strataCounts[i]
    print("[%1.2e, %1.2e) = %i (%1.2f%%) %i" % (strataBounds[i], strataBounds[i+1],
      strataCounts[i], 100*(strataCounts[i]/N), total))

Setup.hs

import Distribution.Simple
main = defaultMain