diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c)2011, Qi Qi
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Qi Qi nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,4 @@
+import Numeric.LinearAlgebra
+import Dynamic
+import System.Environment
+main = defaultMain
diff --git a/hVOIDP.cabal b/hVOIDP.cabal
new file mode 100644
--- /dev/null
+++ b/hVOIDP.cabal
@@ -0,0 +1,35 @@
+Name:	hVOIDP
+Version: 1.0
+Synopsis: Optimal variable selection in chain graphical model.
+Description: 
+   An implementation of the VOIDP algorithm that selects the optimal observations in chain graphical models.
+
+    .
+
+    References:
+    .
+
+    * Andreas Krause and Carlos Guestrin. 2009. Optimal value of information in graphical models. J. Artif. Int. Res. 35, 1 (July 2009), 557-591. 
+
+    .
+    * Qi Qi; Yi Shang; Hongchi Shi; , "An improved algorithm for optimal subset selection in chain graphical models," Evolutionary Computation (CEC), 2010 IEEE Congress on , vol., no., pp.1-6, 18-23 July 2010
+   
+   
+License:     GPL
+License-file: LICENSE
+Author:      Qi Qi <qiqi789@gmail.com>
+Maintainer:  Qi Qi <qiqi789@gmail.com>
+Homepage:    http://dslsrv4.cs.missouri.edu/~qqbm9
+Category:    AI
+Build-Type:  Simple
+cabal-version:  >=1.6
+
+Source-Repository head
+  Type:     git
+  Location: git://github.com/qiqi789/VOIDP--haskell-implementation-for-the-optimal-value-of-information-in-chain-models.git
+
+Executable hVOIDP
+  Main-is:    voidp.hs
+  Build-Depends:   array,hmatrix,base >=3 && <5
+  extra-libraries: lapack blas
+  extra-lib-dirs:  /usr/lib 
diff --git a/voidp.hs b/voidp.hs
new file mode 100644
--- /dev/null
+++ b/voidp.hs
@@ -0,0 +1,175 @@
+import Numeric.LinearAlgebra
+
+import Dynamic
+
+import System.Environment
+
+-- compile: ghc --make -O2 voidp.hs -rtsopts -optl -llapack -lblas -fforce-recomp
+-- run: voidp +RTS -Ksize 
+-- run: voidp +RTS -H    ; make it run a little faster
+
+disp :: Matrix Double -> IO ( )
+disp = putStr . dispf 6
+
+dispv :: Vector Double -> IO ( )
+dispv = putStr . vecdisp (dispf 4)
+
+
+qt :: Int -> Matrix Double -> Matrix Double -> Matrix Double -> Matrix Double
+qt n prior transmat obsmat
+    | n==1 = prior
+    | otherwise = fromColumns $ [fromList ( map (foldVector (+) 0.0) ( toColumns qtM) )]
+                  where qtM = qt' <> ((1><4) [1,1,1,1::Double]) * transmat             
+                        qt' = qt (n-1) prior transmat obsmat
+
+mQT :: Int -> Matrix Double -> Matrix Double -> Matrix Double -> Matrix Double
+mQT n prior transmat obsmat = fromBlocks [[qt i1 prior transmat obsmat | i1<-[1..n]]]
+
+
+ot :: Int -> Matrix Double -> Matrix Double -> Matrix Double -> Matrix Double
+ot n prior transmat obsmat  = fromColumns $ [fromList ( map (foldVector (+) 0.0) ( toColumns otM) )]
+                                 where otM = ot' <> ones * obsmat             
+                                       ot' = qt n prior transmat obsmat
+                                       ones = repmat (ident 1) 1 10
+
+mOT :: Int -> Matrix Double -> Matrix Double -> Matrix Double -> Matrix Double
+mOT n prior transmat obsmat = fromBlocks [[ot i2 prior transmat obsmat | i2<-[1..n]]]
+
+
+entropy :: Vector Double -> Double
+entropy v = vec' <.> logvec'
+              where vec' = fromList $ filter (>0) (toList v)
+                    logvec' = mapVector (logBase 2) vec'
+
+rew_base_ary :: Int -> Matrix Double -> Matrix Double -> Matrix Double -> Vector Double
+rew_base_ary n prior transmat obsmat = fromList $ [entropy v | v <- toColumns (mOT n prior transmat obsmat)]
+           
+                                       
+exp_rew :: Vector Double -> Matrix Double -> Vector Double -> Double
+exp_rew vq obsmat vo = sumElements $ zipVectorWith (*) vqo vo'
+                       where vqm = ( fromColumns [vq] <> (repmat (ident 1) 1 (cols obsmat ) ) ) * obsmat
+                             idx = find (>0) vo
+                             vo' = fromList $ map (atIndex vo) idx
+                             vqm' = trans (extractRows idx $ trans vqm)
+                             vo'' = mapVector recip vo'
+                             vo''m = (repmat (ident 1) (dim vq )  1) <> (fromRows [vo''])
+                             mqo = vqm'  * vo''m
+                             vqo = fromList $ [entropy v' | v <- toColumns mqo, let vm =  (fromColumns [v]) <> (repmat (ident 1) 1 (cols obsmat) ) * obsmat, let v' = fromList $ (map sumElements (toColumns vm)) ]
+                             
+            
+exp_rew_ary :: Int -> Matrix Double -> Matrix Double -> Matrix Double -> Vector Double
+exp_rew_ary n mq obsmat mo = fromList $ [exp_rew (vqlist!!ii) obsmat (volist!!ii)  | ii<-[0..n-1], let vqlist=toColumns mq, let volist =toColumns mo] 
+
+
+type Entry = (Double,Int)
+type Coord = (Int,Int,Int)
+
+
+onePoint :: Int -> Int -> Int -> Int -> Matrix Double -> Matrix Double -> Matrix Double -> Matrix Double -> Vector Double
+onePoint a b k n mot mqt obsmat transmat  | n == a+1 = qt1' / otvec
+                                          | otherwise = (qt2' . qt2) $ (onePoint a b k (n-1) mot mqt obsmat transmat)
+                                            where  qt1 = (asColumn ((toColumns mqt) !! (a-1)) ) <> (asRow $ constant 1.0 4) 
+                                                   ok = (asColumn ((toColumns obsmat) !! (k-1)) ) <>  (asRow $ constant 1.0 4)
+                                                   qt1' = fromList (map sumElements ( toColumns (qt1 * ok * transmat)))
+                                                   otvec =  if (mot @@> (k-1,a-1))/=0 
+                                                            then constant (mot @@> (k-1,a-1)) 4
+                                                            else  error "zero error"
+                                                   qt2 vec = (asColumn vec) <> (asRow $ constant 1.0 4) * transmat
+                                                   qt2' mat = fromList (map sumElements (toColumns mat))
+
+
+abPoints :: Int -> Int -> Int -> Matrix Double -> Matrix Double -> Matrix Double -> Matrix Double -> [Vector Double]
+abPoints a b k mot mqt obsmat transmat = take (b-a-1) restqt
+                                          where qt1 = (asColumn ((toColumns mqt) !! (a-1)) ) <> (asRow $ constant 1.0 4) 
+                                                ok = (asColumn ((toColumns obsmat) !! (k-1)) ) <>  (asRow $ constant 1.0 4)
+                                                qt1' = fromList (map sumElements ( toColumns (qt1 * ok * transmat)))
+                                                otvec =  if (mot @@> (k-1,a-1))/=0 
+                                                         then constant (mot @@> (k-1,a-1)) 4
+                                                         else  error "zero error"
+                                                qt1'' = (qt1' / otvec)
+                                                qt2 vec = (asColumn vec) <> (asRow $ constant 1.0 4) * transmat
+                                                qt2' mat = fromList (map sumElements (toColumns mat))
+                                                restqt = iterate (qt2' . qt2) qt1''
+       
+
+calc_Lab0 :: Int -> Int -> Matrix Double -> Matrix Double -> Matrix Double -> Matrix Double -> Vector Double -> Entry
+calc_Lab0 a b mot mqt transmat obsmat rewBases 
+  | a >= (b-1) = (0.0, 0)
+  | a == 0 = let a0 = foldVector (+) (0.0 :: Double) $ subVector 0 (b-a-1) rewBases
+             in (a0, 0)
+  | otherwise = ((sumElements $ fromRows [oneCase ii |ii<-[1..10] ]), 0)
+    where dif_n = b-a-1
+          oneCase k | mot @@> (k-1,a-1) == 0 = constant 0.0 dif_n
+--                    | otherwise = fromList [ cal_entropy (onePoint a b k t mot mqt obsmat transmat ) | t <- [a+1..b-1]]
+                    | otherwise = fromList [ (cal_entropy v) * (mot @@> (k-1,a-1 ) ) | v <- (abPoints a b k mot mqt obsmat transmat) ]
+                        where cal_entropy vec = let vm =  (asColumn vec) <> (asRow $ constant 1.0 (cols obsmat) ) * obsmat 
+                                                    v' = fromList $ map sumElements (toColumns vm)
+                                                in entropy v'
+
+
+calc_Labk :: Int -> Int -> Int -> Vector Double -> Table Entry Coord -> Entry
+calc_Labk a b k expRews t 
+  | a >= (b-1) = (0.0, 0)
+  | otherwise = if max_selj > sel0 then (max_selj, max_selj_i) else (sel0, -1)
+                where selj j = (expRews @> (j-1)) - 0.0 + (fst $ findTable t (0,a,j)) + (fst $ findTable t (k-1,j,b)) 
+                      selj_ary = fromList [selj jj |jj <- [a+1..b-1]]
+                      max_selj = maxElement selj_ary
+                      max_selj_i = maxIndex selj_ary + (a+1)
+                      sel0 = fst $ findTable t (0,a,b)
+                             
+                               
+                            
+computeTbl :: Int -> Matrix Double -> Matrix Double -> Matrix Double -> Matrix Double -> Vector Double -> Vector Double -> Table Entry Coord -> Coord -> Entry
+computeTbl budget mot mqt transmat obsmat rewBases expRews table (k,ii,j)
+  | k==0 = calc_Lab0 ii j mot mqt transmat obsmat rewBases
+  | k>0 && k<=budget =  calc_Labk ii j k expRews table
+  | otherwise = (-100.0, -1)
+    
+    
+bndsTbl :: Int -> Int -> ((Int,Int,Int),(Int,Int,Int))  
+bndsTbl n budget = ((0,0,0),(budget,n+1,n+1))  -- n-1,n+1
+
+selsList :: Int -> Int -> Int-> Table Entry Coord  -> [Int]
+selsList k a b table =
+      let sel = snd $ findTable table (k,a,b)
+          a' = sel
+          k' = k -1
+      in sel: if k'>0 && sel>0 
+              then selsList k' a' b table
+              else []
+           
+                   
+
+selectAlgo :: Int -> Int -> Matrix Double -> Matrix Double -> Matrix Double -> Matrix Double -> Vector Double -> Vector Double -> [Int]
+selectAlgo n budget mot mqt transmat obsmat rewBases expRews = selsList budget 0 (n+1) t
+  where t = dynamic (computeTbl budget mot mqt transmat obsmat rewBases expRews) (bndsTbl n budget)
+        
+
+
+main :: IO ( )
+main = do
+     args <- getArgs
+     prior <- loadMatrix $ args !! 0  --"hmm_prior.txt"
+     transmat <- loadMatrix $ args!!1 --"hmm_transmat.txt"
+     obsmat <- loadMatrix $ args!!2 --"hmm_obsmat.txt"
+
+--     disp $ mQT 24 prior transmat obsmat
+--     disp $ mOT 24 prior transmat obsmat
+  
+--     dispv $ rew_base_ary 24 prior transmat obsmat
+--     print $ sumElements $ rew_base_ary 24 prior transmat obsmat
+     let n = read $ args!!3 --24 -- 24
+         mqt = mQT n prior transmat obsmat
+         mot = mOT n prior transmat obsmat
+         expRews =  exp_rew_ary n mqt obsmat mot
+         rewBases = rew_base_ary n prior transmat obsmat
+         budget = read $ args!!4
+         sels = selectAlgo n budget mot mqt transmat obsmat rewBases expRews
+         --entry = calc_Lab0 2 4 mot mqt transmat obsmat rewBases
+         --val = onePoint 1 4 3 mot mqt obsmat transmat
+         --t = dynamic (computeTbl budget mot mqt transmat obsmat rewBases expRews) (bndsTbl n budget)
+         --entry = calc_Labk 2 7 5 expRews t
+     print $ sels
+     --print $ val
+     --dispv $ val
+