{-# LANGUAGE TypeFamilies, FlexibleInstances #-}

-- |
-- Module     : Simulation.Aivika.Branch.Generator
-- Copyright  : Copyright (c) 2016-2017, David Sorokin <david.sorokin@gmail.com>
-- License    : BSD3
-- Maintainer : David Sorokin <david.sorokin@gmail.com>
-- Stability  : experimental
-- Tested with: GHC 7.10.3
--
-- Here is defined a random number generator,
-- where 'BR' can be an instance of 'MonadGenerator'.
--
module Simulation.Aivika.Branch.Generator () where

import Control.Monad
import Control.Monad.Trans

import qualified System.Random.MWC as MWC
import System.Random

import Data.IORef

import Simulation.Aivika.Trans
import Simulation.Aivika.Trans.Generator.Primitive
import Simulation.Aivika.Branch.Internal.BR

instance MonadGenerator (BR IO) where

  data Generator (BR IO) =
    Generator { Generator (BR IO) -> BR IO Double
generator01 :: BR IO Double,
                -- ^ the generator of uniform numbers from 0 to 1
                Generator (BR IO) -> BR IO Double
generatorNormal01 :: BR IO Double,
                -- ^ the generator of normal numbers with mean 0 and variance 1
                Generator (BR IO) -> BR IO Int
generatorSequenceNo :: BR IO Int
                -- ^ the generator of sequence numbers
              }

  generateUniform :: Generator (BR IO) -> Double -> Double -> BR IO Double
generateUniform = BR IO Double -> Double -> Double -> BR IO Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateUniform01 (BR IO Double -> Double -> Double -> BR IO Double)
-> (Generator (BR IO) -> BR IO Double)
-> Generator (BR IO)
-> Double
-> Double
-> BR IO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (BR IO) -> BR IO Double
generator01

  generateUniformInt :: Generator (BR IO) -> Int -> Int -> BR IO Int
generateUniformInt = BR IO Double -> Int -> Int -> BR IO Int
forall (m :: * -> *). Monad m => m Double -> Int -> Int -> m Int
generateUniformInt01 (BR IO Double -> Int -> Int -> BR IO Int)
-> (Generator (BR IO) -> BR IO Double)
-> Generator (BR IO)
-> Int
-> Int
-> BR IO Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (BR IO) -> BR IO Double
generator01

  generateTriangular :: Generator (BR IO) -> Double -> Double -> Double -> BR IO Double
generateTriangular = BR IO Double -> Double -> Double -> Double -> BR IO Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> Double -> m Double
generateTriangular01 (BR IO Double -> Double -> Double -> Double -> BR IO Double)
-> (Generator (BR IO) -> BR IO Double)
-> Generator (BR IO)
-> Double
-> Double
-> Double
-> BR IO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (BR IO) -> BR IO Double
generator01

  generateNormal :: Generator (BR IO) -> Double -> Double -> BR IO Double
generateNormal = BR IO Double -> Double -> Double -> BR IO Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateNormal01 (BR IO Double -> Double -> Double -> BR IO Double)
-> (Generator (BR IO) -> BR IO Double)
-> Generator (BR IO)
-> Double
-> Double
-> BR IO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (BR IO) -> BR IO Double
generatorNormal01

  generateLogNormal :: Generator (BR IO) -> Double -> Double -> BR IO Double
generateLogNormal = BR IO Double -> Double -> Double -> BR IO Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateLogNormal01 (BR IO Double -> Double -> Double -> BR IO Double)
-> (Generator (BR IO) -> BR IO Double)
-> Generator (BR IO)
-> Double
-> Double
-> BR IO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (BR IO) -> BR IO Double
generatorNormal01

  generateExponential :: Generator (BR IO) -> Double -> BR IO Double
generateExponential = BR IO Double -> Double -> BR IO Double
forall (m :: * -> *). Monad m => m Double -> Double -> m Double
generateExponential01 (BR IO Double -> Double -> BR IO Double)
-> (Generator (BR IO) -> BR IO Double)
-> Generator (BR IO)
-> Double
-> BR IO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (BR IO) -> BR IO Double
generator01

  generateErlang :: Generator (BR IO) -> Double -> Int -> BR IO Double
generateErlang = BR IO Double -> Double -> Int -> BR IO Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Int -> m Double
generateErlang01 (BR IO Double -> Double -> Int -> BR IO Double)
-> (Generator (BR IO) -> BR IO Double)
-> Generator (BR IO)
-> Double
-> Int
-> BR IO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (BR IO) -> BR IO Double
generator01

  generatePoisson :: Generator (BR IO) -> Double -> BR IO Int
generatePoisson = BR IO Double -> Double -> BR IO Int
forall (m :: * -> *). Monad m => m Double -> Double -> m Int
generatePoisson01 (BR IO Double -> Double -> BR IO Int)
-> (Generator (BR IO) -> BR IO Double)
-> Generator (BR IO)
-> Double
-> BR IO Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (BR IO) -> BR IO Double
generator01

  generateBinomial :: Generator (BR IO) -> Double -> Int -> BR IO Int
generateBinomial = BR IO Double -> Double -> Int -> BR IO Int
forall (m :: * -> *). Monad m => m Double -> Double -> Int -> m Int
generateBinomial01 (BR IO Double -> Double -> Int -> BR IO Int)
-> (Generator (BR IO) -> BR IO Double)
-> Generator (BR IO)
-> Double
-> Int
-> BR IO Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (BR IO) -> BR IO Double
generator01

  generateGamma :: Generator (BR IO) -> Double -> Double -> BR IO Double
generateGamma Generator (BR IO)
g = BR IO Double -> BR IO Double -> Double -> Double -> BR IO Double
forall (m :: * -> *).
Monad m =>
m Double -> m Double -> Double -> Double -> m Double
generateGamma01 (Generator (BR IO) -> BR IO Double
generatorNormal01 Generator (BR IO)
g) (Generator (BR IO) -> BR IO Double
generator01 Generator (BR IO)
g)

  generateBeta :: Generator (BR IO) -> Double -> Double -> BR IO Double
generateBeta Generator (BR IO)
g = BR IO Double -> BR IO Double -> Double -> Double -> BR IO Double
forall (m :: * -> *).
Monad m =>
m Double -> m Double -> Double -> Double -> m Double
generateBeta01 (Generator (BR IO) -> BR IO Double
generatorNormal01 Generator (BR IO)
g) (Generator (BR IO) -> BR IO Double
generator01 Generator (BR IO)
g)

  generateWeibull :: Generator (BR IO) -> Double -> Double -> BR IO Double
generateWeibull = BR IO Double -> Double -> Double -> BR IO Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateWeibull01 (BR IO Double -> Double -> Double -> BR IO Double)
-> (Generator (BR IO) -> BR IO Double)
-> Generator (BR IO)
-> Double
-> Double
-> BR IO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (BR IO) -> BR IO Double
generator01

  generateDiscrete :: Generator (BR IO) -> DiscretePDF a -> BR IO a
generateDiscrete = BR IO Double -> DiscretePDF a -> BR IO a
forall (m :: * -> *) a. Monad m => m Double -> DiscretePDF a -> m a
generateDiscrete01 (BR IO Double -> DiscretePDF a -> BR IO a)
-> (Generator (BR IO) -> BR IO Double)
-> Generator (BR IO)
-> DiscretePDF a
-> BR IO a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (BR IO) -> BR IO Double
generator01

  generateSequenceNo :: Generator (BR IO) -> BR IO Int
generateSequenceNo = Generator (BR IO) -> BR IO Int
generatorSequenceNo

  newGenerator :: GeneratorType (BR IO) -> BR IO (Generator (BR IO))
newGenerator GeneratorType (BR IO)
tp =
    case GeneratorType (BR IO)
tp of
      GeneratorType (BR IO)
SimpleGenerator ->
        do let g :: IO (IO Double)
g = Gen (PrimState IO) -> IO Double
forall a (m :: * -> *).
(Variate a, PrimMonad m) =>
Gen (PrimState m) -> m a
MWC.uniform (Gen (PrimState IO) -> IO Double)
-> IO (Gen (PrimState IO)) -> IO (IO Double)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
                   IO (Gen (PrimState IO))
MWC.createSystemRandom
           IO Double
g' <- IO (IO Double) -> BR IO (IO Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO IO (IO Double)
g
           BR IO Double -> BR IO (Generator (BR IO))
forall (m :: * -> *).
MonadGenerator m =>
m Double -> m (Generator m)
newRandomGenerator01 (IO Double -> BR IO Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO IO Double
g')
      SimpleGeneratorWithSeed Word32
x ->
        [Char] -> BR IO (Generator (BR IO))
forall a. HasCallStack => [Char] -> a
error [Char]
"Unsupported generator type SimpleGeneratorWithSeed: newGenerator"
      CustomGenerator BR IO (Generator (BR IO))
g ->
        BR IO (Generator (BR IO))
g
      CustomGenerator01 BR IO Double
g ->
        BR IO Double -> BR IO (Generator (BR IO))
forall (m :: * -> *).
MonadGenerator m =>
m Double -> m (Generator m)
newRandomGenerator01 BR IO Double
g

  newRandomGenerator :: g -> BR IO (Generator (BR IO))
newRandomGenerator g
g = 
    do IORef g
r <- IO (IORef g) -> BR IO (IORef g)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef g) -> BR IO (IORef g))
-> IO (IORef g) -> BR IO (IORef g)
forall a b. (a -> b) -> a -> b
$ g -> IO (IORef g)
forall a. a -> IO (IORef a)
newIORef g
g
       let g01 :: BR IO Double
g01 = do g
g <- IO g -> BR IO g
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO g -> BR IO g) -> IO g -> BR IO g
forall a b. (a -> b) -> a -> b
$ IORef g -> IO g
forall a. IORef a -> IO a
readIORef IORef g
r
                    let (Double
x, g
g') = g -> (Double, g)
forall a g. (Random a, RandomGen g) => g -> (a, g)
random g
g
                    IO () -> BR IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> BR IO ()) -> IO () -> BR IO ()
forall a b. (a -> b) -> a -> b
$ IORef g -> g -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef g
r g
g'
                    Double -> BR IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
x
       BR IO Double -> BR IO (Generator (BR IO))
forall (m :: * -> *).
MonadGenerator m =>
m Double -> m (Generator m)
newRandomGenerator01 BR IO Double
g01

  newRandomGenerator01 :: BR IO Double -> BR IO (Generator (BR IO))
newRandomGenerator01 BR IO Double
g01 =
    do BR IO Double
gNormal01 <- BR IO Double -> BR IO (BR IO Double)
newNormalGenerator01 BR IO Double
g01
       IORef Int
gSeqNoRef <- IO (IORef Int) -> BR IO (IORef Int)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Int) -> BR IO (IORef Int))
-> IO (IORef Int) -> BR IO (IORef Int)
forall a b. (a -> b) -> a -> b
$ Int -> IO (IORef Int)
forall a. a -> IO (IORef a)
newIORef Int
0
       let gSeqNo :: BR IO Int
gSeqNo =
             do Int
x <- IO Int -> BR IO Int
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Int -> BR IO Int) -> IO Int -> BR IO Int
forall a b. (a -> b) -> a -> b
$ IORef Int -> IO Int
forall a. IORef a -> IO a
readIORef IORef Int
gSeqNoRef
                IO () -> BR IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> BR IO ()) -> IO () -> BR IO ()
forall a b. (a -> b) -> a -> b
$ IORef Int -> (Int -> Int) -> IO ()
forall a. IORef a -> (a -> a) -> IO ()
modifyIORef' IORef Int
gSeqNoRef (Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
                Int -> BR IO Int
forall (m :: * -> *) a. Monad m => a -> m a
return Int
x
       Generator (BR IO) -> BR IO (Generator (BR IO))
forall (m :: * -> *) a. Monad m => a -> m a
return Generator :: BR IO Double -> BR IO Double -> BR IO Int -> Generator (BR IO)
Generator { generator01 :: BR IO Double
generator01 = BR IO Double
g01,
                          generatorNormal01 :: BR IO Double
generatorNormal01 = BR IO Double
gNormal01,
                          generatorSequenceNo :: BR IO Int
generatorSequenceNo = BR IO Int
gSeqNo }

-- | Create a normal random number generator with mean 0 and variance 1
-- by the specified generator of uniform random numbers from 0 to 1.
newNormalGenerator01 :: BR IO Double
                        -- ^ the generator
                        -> BR IO (BR IO Double)
newNormalGenerator01 :: BR IO Double -> BR IO (BR IO Double)
newNormalGenerator01 BR IO Double
g =
  do IORef Double
nextRef <- IO (IORef Double) -> BR IO (IORef Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Double) -> BR IO (IORef Double))
-> IO (IORef Double) -> BR IO (IORef Double)
forall a b. (a -> b) -> a -> b
$ Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     IORef Bool
flagRef <- IO (IORef Bool) -> BR IO (IORef Bool)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Bool) -> BR IO (IORef Bool))
-> IO (IORef Bool) -> BR IO (IORef Bool)
forall a b. (a -> b) -> a -> b
$ Bool -> IO (IORef Bool)
forall a. a -> IO (IORef a)
newIORef Bool
False
     IORef Double
xi1Ref  <- IO (IORef Double) -> BR IO (IORef Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Double) -> BR IO (IORef Double))
-> IO (IORef Double) -> BR IO (IORef Double)
forall a b. (a -> b) -> a -> b
$ Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     IORef Double
xi2Ref  <- IO (IORef Double) -> BR IO (IORef Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Double) -> BR IO (IORef Double))
-> IO (IORef Double) -> BR IO (IORef Double)
forall a b. (a -> b) -> a -> b
$ Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     IORef Double
psiRef  <- IO (IORef Double) -> BR IO (IORef Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Double) -> BR IO (IORef Double))
-> IO (IORef Double) -> BR IO (IORef Double)
forall a b. (a -> b) -> a -> b
$ Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     let loop :: BR IO ()
loop =
           do Double
psi <- IO Double -> BR IO Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> BR IO Double) -> IO Double -> BR IO Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
psiRef
              if (Double
psi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
1.0) Bool -> Bool -> Bool
|| (Double
psi Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0.0)
                then do Double
g1 <- BR IO Double
g
                        Double
g2 <- BR IO Double
g
                        let xi1 :: Double
xi1 = Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
g1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1.0
                            xi2 :: Double
xi2 = Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
g2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1.0
                            psi :: Double
psi = Double
xi1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
xi1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
xi2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
xi2
                        IO () -> BR IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> BR IO ()) -> IO () -> BR IO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi1Ref Double
xi1
                        IO () -> BR IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> BR IO ()) -> IO () -> BR IO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi2Ref Double
xi2
                        IO () -> BR IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> BR IO ()) -> IO () -> BR IO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef Double
psi
                        BR IO ()
loop
                else IO () -> BR IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> BR IO ()) -> IO () -> BR IO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef (Double -> IO ()) -> Double -> IO ()
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
sqrt (- Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
psi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
psi)
     BR IO Double -> BR IO (BR IO Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (BR IO Double -> BR IO (BR IO Double))
-> BR IO Double -> BR IO (BR IO Double)
forall a b. (a -> b) -> a -> b
$
       do Bool
flag <- IO Bool -> BR IO Bool
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Bool -> BR IO Bool) -> IO Bool -> BR IO Bool
forall a b. (a -> b) -> a -> b
$ IORef Bool -> IO Bool
forall a. IORef a -> IO a
readIORef IORef Bool
flagRef
          if Bool
flag
            then do IO () -> BR IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> BR IO ()) -> IO () -> BR IO ()
forall a b. (a -> b) -> a -> b
$ IORef Bool -> Bool -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Bool
flagRef Bool
False
                    IO Double -> BR IO Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> BR IO Double) -> IO Double -> BR IO Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
nextRef
            else do IO () -> BR IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> BR IO ()) -> IO () -> BR IO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi1Ref Double
0.0
                    IO () -> BR IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> BR IO ()) -> IO () -> BR IO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi2Ref Double
0.0
                    IO () -> BR IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> BR IO ()) -> IO () -> BR IO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef Double
0.0
                    BR IO ()
loop
                    Double
xi1 <- IO Double -> BR IO Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> BR IO Double) -> IO Double -> BR IO Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
xi1Ref
                    Double
xi2 <- IO Double -> BR IO Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> BR IO Double) -> IO Double -> BR IO Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
xi2Ref
                    Double
psi <- IO Double -> BR IO Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> BR IO Double) -> IO Double -> BR IO Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
psiRef
                    IO () -> BR IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> BR IO ()) -> IO () -> BR IO ()
forall a b. (a -> b) -> a -> b
$ IORef Bool -> Bool -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Bool
flagRef Bool
True
                    IO () -> BR IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> BR IO ()) -> IO () -> BR IO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
nextRef (Double -> IO ()) -> Double -> IO ()
forall a b. (a -> b) -> a -> b
$ Double
xi2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
psi
                    Double -> BR IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> BR IO Double) -> Double -> BR IO Double
forall a b. (a -> b) -> a -> b
$ Double
xi1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
psi