Thursday, June 14, 2012

Getting Started with Simulations

To start numerical simulations, we first need the ability to set the initial condition. It will also be a good excercise for Paraiso building blocks, so let's give it a try. Please look at the example.

table :: Named (StaticValue TArray Double)
table = "table" `isNameOf` StaticValue TArray undefined

createBuilder :: Builder Vec2 Int Annotation ()
createBuilder = do
  x01 <- bind $ (cast $ loadIndex (Axis 0)) 
              / (cast $ broadcast $ loadSize (Axis 0))
  y01 <- bind $ (cast $ loadIndex (Axis 1)) 
              / (cast $ broadcast $ loadSize (Axis 1))
  x <- bind $ 4 * (x01-0.5)
  y <- bind $ 5 * (y01-0.5)
  z <- bind $ atan((1- x^2 - (y-(x^2)**(1/3))^2)*10)
  store table z

The Generator.hs makes a bit tedious use of casts and broadcasts we have just learned, to calculate x and y coordinates in real number inbetween 0 and 1 (they are namely x01 and y01.) See how input and output types for casts and broadcasts are inferred. The index and size are of type Int (as specified in the type annotation for createBuilder) while x01, y01, x, y and z are all of type (Builder of) Double because table has that type and we cannot convert type without an explicit cast. The careful reader also might have noticed that 1/3 is treated as Double instead of just being truncated to 0.

Now, Paraiso and the author give hearts to the readers! Thank you!

Thursday, June 7, 2012

Life in Paraiso

Now, let's try Paraiso in a real-world problem. Simulation of artificial life will be a good starting point.

The sample program, like every other Haskell programs, starts with some language extension and imports.

#!/usr/bin/env runhaskell
{-# LANGUAGE NoImplicitPrelude, OverloadedStrings #-}
{-# OPTIONS -Wall #-}

import           Control.Monad
import           Data.Tensor.TypeLevel
import           Language.Paraiso.Annotation (Annotation)
import qualified Language.Paraiso.Annotation.Boundary as Boundary
import           Language.Paraiso.Generator (generateIO)
import qualified Language.Paraiso.Generator.Native as Native
import           Language.Paraiso.Name
import           Language.Paraiso.OM
import           Language.Paraiso.OM.Builder
import           Language.Paraiso.OM.Builder.Boolean (select,eq,ge,le)
import           Language.Paraiso.OM.DynValue as DVal
import           Language.Paraiso.OM.Realm 
import qualified Language.Paraiso.OM.Reduce as Reduce
import           Language.Paraiso.OM.Value (StaticValue(..))
import           Language.Paraiso.Optimization
import           Language.Paraiso.Prelude
import           NumericPrelude hiding ((||),(&&))

-- the main program
main :: IO ()
main = do
  _ <- generateIO mySetup myOM
  return ()

Let's say our computation region to be 80x48 and with cyclic boundary conditions. The new Orthotope Machine will have the name "Life".

-- the code generation setup
mySetup :: Native.Setup Vec2 Int
mySetup = 
  (Native.defaultSetup $ Vec :~ 80 :~ 48)
  { = "./dist/" ,
    Native.boundary = compose $ const Boundary.Cyclic

-- the orthotope machine to be generated
myOM :: OM Vec2 Int Annotation
myOM = optimize O3 $
  makeOM (mkName "Life") [] myVars myKernels

We use an array variable cell for the cell states, population for counting the number of alive cells and generation for keeping track of the time.

-- the variables we use
cell :: Named (StaticValue TArray Int)
cell = "cell" `isNameOf` StaticValue TArray  undefined

population :: Named (StaticValue TScalar Int)
population = "population" `isNameOf` StaticValue TScalar undefined

generation :: Named (StaticValue TScalar Int)
generation = "generation" `isNameOf` StaticValue TScalar undefined

myVars :: [Named DynValue]
myVars = [f2d cell, f2d population, f2d generation]

Then, let's make two kernels, one for the initialization of the state and the other for updating the state for one generation.

-- our kernel
myKernels :: [Named (Builder Vec2 Int Annotation ())]
myKernels = ["init" `isNameOf` initBuilder,
             "proceed" `isNameOf` proceedBuilder]

Initialization, is just trivial.

initBuilder :: Builder Vec2 Int Annotation ()
initBuilder = do 
  -- store the initial states.
  store cell 0
  store population 0
  store generation 0

Now, how shall we write the rule of Conway's game of Life? Let's first define the adjacency. In Conway's game of Life, one cell has eight neighbours:

adjVecs :: [Vec2 Int]
adjVecs = zipWith (\x y -> Vec :~ x :~ y)

A timestep of the game begins by loading a Array variable called "cell" as the old state of the simulation.

proceedBuilder :: Builder Vec2 Int Annotation ()
proceedBuilder = do 
  oldCell <- bind $ load cell

We also load a Scalar variable called "generation."

  gen  <- bind $ load generation

shiftedCell <- bind $ shift v oldCell is an expression that yields a cell pattern shifted by amount v from the old cell. Then mapping it over every v in adjVecs creates the neighbour list.

  neighbours <- forM adjVecs $
    \v -> bind $ shift v oldCell

Count the neighbour cells that are alive.

  num <- bind $ sum neighbours

Apply the rule of Conway's game of Life.

  isAlive <- bind $
             (oldCell `eq` 0) && (num `eq` 3) ||
             (oldCell `eq` 1) && (num `ge` 2) && (num `le` 3) 

Update the cell according to the rule. The expression select isAlive 1 0, means isAlive ? 1 : 0 in C.

  newCell <- bind $ select isAlive 1 0

Count the number of alive cells, increment the generation, and store the new cell state.

  store population $ reduce Reduce.Sum newCell
  store generation $ gen + 1
  store cell $ newCell

Wednesday, June 6, 2012

Boundary Condition in Paraiso

The shift instruction of Orthotope Machine is displace the array by a constant vector. Because every array has a finite size, the users must set the boundary condition to specify how the out-of-bound access are treated.
Current Paraiso supports two boundary conditions; open boundary and cyclic boundary.
You can choose boundary condition for every dimension independently
by setting a value of the type vector of Condition in boundary field of the Setup data.
Of the two, the cyclic boundary condition is simple. Out-of-bound accesses see the other side of the array.
The open boundary condition is rather complicated. First, the size of margin regions are determined so that any Paraiso kernel can cover the region you said you wanted to calculate. You have specified the region 0 <= i < localSize in the Setup. Next, each kernel start the calculation from the maximum region including the margins, and calculates as large region as possible. The array values out of that region are undefined. Also, the reduction operators always reduce over the entire margined region. You need to manually mask the margins if you want to.
Please look at the sample program to understand how Paraiso deals with the boundary conditions. This program has three kernels:
myKernels :: [Named (Builder Vec1 Int Annotation ())]
myKernels =
  ["init" `isNameOf` do
      store table $ loadIndex (Axis 0)

  ,"increment" `isNameOf` do
      store table $ 1 + load table

  ,"calculate" `isNameOf` do
      center <- bind $ load table
      right  <- bind $ shift (Vec :~ (-1)) center
      left   <- bind $ shift (Vec :~ ( 1)) center
      ret    <- bind $ 10000 * left + 100 * center + right
      store table ret
      store total $ reduce Reduce.Sum ret
The first kernel initializes the array elements with their respective indices. The second adds 1 to every element. The third reads adjacent cells and perform some calculations. Let's give it a go:
  cout << "total: " << << endl;
  cout << endl;
Say make and two program will be generated. They are the same but for the boundary conditions. Under the cyclic boundary condition, you can read across the boundary, as follows:
$ ./main-cyclic.out
index:      0      1      2      3      4      5      6      7
value:      0      1      2      3      4      5      6      7
index:      0      1      2      3      4      5      6      7
value:      1      2      3      4      5      6      7      8
index:      0      1      2      3      4      5      6      7
value:  80102  10203  20304  30405  40506  50607  60708  70801
total: 363636
On the other hand if you use the open boundary condition, Paraiso adds the margin regions for you, and you can access indices beyond 0 and localSize.
$ ./main-open.out
index:     -1      0      1      2      3      4      5      6      7      8
value:     -1      0      1      2      3      4      5      6      7      8
value:      0      1      2      3      4      5      6      7      8      9
value:      0    102  10203  20304  30405  40506  50607  60708  70809      0
total: 283644

Tuesday, June 5, 2012

Builder Monad

--   type level data
  Value rea con = 
  -- | data obtained from the data-flow graph.
  -- insignificant and can be 'undefined'.
  -- | data obtained as an immediate value.

bind trick

bind = fmap return
Like this:
  x <- bind $ loadIndex (Axis 0) 
  y <- bind $ loadIndex (Axis 1) 
  z <- bind $ x*y
  w <- bind $ x+y
The above program originally looked like this:
  x <- loadIndex (Axis 0) 
  y <- loadIndex (Axis 1) 
  z <- return x * return y
  w <- return x + return y

OM instructions

Here is the instruction set:
data Inst vector gauge 
  = Load StaticIdx
  | Store StaticIdx
  | Reduce R.Operator 
  | Broadcast 
  | LoadIndex (Axis vector) 
  | LoadSize (Axis vector) 
  | Shift (vector gauge) 
  | Imm Dynamic 
  | Arith A.Operator 
and the corresponding monad combinators:
--           options                    input nodes            output nodes
imm       :: c                                              -> B (Value r c)
In short,
  • Imm is to introduce an immediate value.
  • Arith is for arithmetic operations.
We will gradually see the detail via examples.

