Updated to add integrated intensity

File now knows how to compute the integrated intensity of the generated
process, and this is written to a file. Next up we need to add in some
proper commandline arguments.

src/HawkesGen.hs

@@ -3,21 +3,28 @@
 module Main
   where
 
+import qualified Data.Text as T
 import qualified Data.Vector.Unboxed as U
 
 import Control.Arrow ((&&&))
-import Control.Monad.Loops (unfoldrM)
+import Data.Double.Conversion.Text (toShortest)
+import Data.Text.IO (writeFile, hPutStrLn)
 import Data.List (foldl')
+import System.IO (openFile, IOMode(..), hClose)
 import System.Random.Mersenne 
   ( MTGen
   , newMTGen
   , random
   )
 
-main :: IO () 
+main :: IO ()
 main = do gen <- newMTGen Nothing
-          evs <- generateHawkes gen hc 1000
-          print evs
+          evs <- generateHawkes gen hc 10000
+          let vec = U.fromList evs
+          let tis = U.zip vec (U.tail vec)
+          h <- openFile "iis.txt" WriteMode
+          U.mapM_ (hPutStrLn h . toShortest) . U.map (integratedIntensity hc vec) $ tis
+          hClose h
   where hc = HC 1.2 0.6 0.8 
 
 -- | This is the context that a user can pass around in order to
@@ -46,7 +53,7 @@ generateHawkes g hc t =
   do let lstar = hcLambda hc
      ev <- (\x -> -(1 / lstar * log x)) `fmap` random g
      let hs = HS g hc t [ev] lstar
-     if ev <= t then hsEvents `fmap` untilM generateHawkes' hs
+     if ev <= t then (reverse . hsEvents) `fmap` untilM generateHawkes' hs
                 else return []
 
 
@@ -83,3 +90,11 @@ intensity :: HawkesContext -- ^ Process Context
 intensity (HC l a b) ev evs = foldl' binop l evs
   where binop z x = z + a * (exp $ -b * (ev - x))
 
+integratedIntensity :: HawkesContext   -- ^ Process context
+                    -> U.Vector Double -- ^ Events
+                    -> (Double,Double) -- ^ Events
+                    -> Double          -- ^ Integrated intenstity
+integratedIntensity (HC l a b) evs (t1,t2) = U.foldl' binop seed . U.filter (<= t1) $ evs
+  where seed =  (t2 - t1) * l
+        ab = a / b
+        binop z x = z + ab * ((exp $ -b * (t1 - x)) - (exp $ -b * (t2 - x)))