Repository: h2r/pomdp-py (language: python)

Issue:
(Redundant) code in POUCT implementation
In the _simulate function of `po_uct.pyx`:
```
        root[action].num_visits += 1
        root.value = root.value + (total_reward - root.value) / (root.num_visits)
        root[action].value = root[action].value + (total_reward - root[action].value) / (root[action].num_visits)
```
Both the value of `root` (VNode) and `root[action]` (QNode) are updated based on `total_reward`. However, in fact, the algorithm in [the paper](https://papers.nips.cc/paper/2010/file/edfbe1afcf9246bb0d40eb4d8027d90f-Paper.pdf) only requires updating the value of the QNode, i.e. `root[action]`.

I also noticed in the source code of the [original author](https://www.davidsilver.uk/applications/) the expected discounted cumulative value is also not maintained in both the VNode and the QNode. 

Also in the current POUCT implementation in pomdp_py, commenting out `root.value = ...` and stick to only updating the QNode's value according to the paper, does not change the output behavior of the planner, since it eventually outputs an action based on the values of the QNodes that are immediate children of the root node. So we should remove this redundant line because it causes confusion.


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index a5c59c3f..4f425a1a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -130,3 +130,7 @@ Thumbs.db
 # Editor backup files #
 #######################
 
+# pomdp-solve generated files
+#############################
+*.alpha
+*.pg
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 61b92444..66a38b46 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,17 @@
+### Version 1.2.2 (pending)
+* Resolved [issue #10](https://github.com/h2r/pomdp-py/issues/10).
+  Set value in V-Node to be the max among its children Q-Nodes.
+  Initial V-Node value set to negative infinity.
+* Avoid search tree building during rollout (thanks Jason)
+* Documentation clarification about `.sample` and `.argmax` functions in the `RewardModel`.
+* Small pomdps (with enumerable state, action, observation spaces)
+  defined in `pomdp_py` can be converted to `.pomdp` and `.pomdpx` file formats.
+* Added interfacing with [pomdp_solve](https://www.pomdp.org/code/) and tested.
+* Added interfacing with [sarsop](https://github.com/AdaCompNUS/sarsop) and tested.
+* Added `utils/templates.py` that contains some convenient implementations of the POMDP interface.
+* Bug fixes (in `histogram.pyx`)
+
+
 ### Version 1.2.1
 
 * Fixed preferred rollout and action prior implementation; Previously the initial visits and values were not applied.
diff --git a/pomdp_py/__init__.py b/pomdp_py/__init__.py
index 5f9ebd4d..28a712f5 100644
--- a/pomdp_py/__init__.py
+++ b/pomdp_py/__init__.py
@@ -1,10 +1,20 @@
 from pomdp_py.utils import util
+
+# Framework
 from pomdp_py.framework.basics import *
 from pomdp_py.framework.oopomdp import *
 from pomdp_py.framework.planner import *
+
+# Representations
 from pomdp_py.representations.distribution.particles import Particles
 from pomdp_py.representations.distribution.histogram import Histogram
 from pomdp_py.representations.distribution.gaussian import Gaussian
+from pomdp_py.representations.belief.histogram import update_histogram_belief
+from pomdp_py.representations.belief.particles import update_particles_belief
+from pomdp_py.utils.interfaces.conversion import to_pomdp_file, to_pomdpx_file, AlphaVectorPolicy, PolicyGraph
+from pomdp_py.utils.interfaces.solvers import vi_pruning, sarsop
+
+# Algorithms
 from pomdp_py.algorithms.value_iteration import ValueIteration  # Cython compiled
 from pomdp_py.algorithms.value_function import value, qvalue
 from pomdp_py.algorithms.pomcp import POMCP
@@ -14,5 +24,6 @@
     tree_stats, ActionPrior
 from pomdp_py.algorithms.bsp.blqr import BLQR
 from pomdp_py.algorithms.visual import visual
-from pomdp_py.representations.belief.histogram import update_histogram_belief
-from pomdp_py.representations.belief.particles import update_particles_belief
\ No newline at end of file
+
+# Templates
+from pomdp_py.utils.templates import *
diff --git a/pomdp_py/algorithms/po_uct.pxd b/pomdp_py/algorithms/po_uct.pxd
index dda45f19..36a6af0a 100644
--- a/pomdp_py/algorithms/po_uct.pxd
+++ b/pomdp_py/algorithms/po_uct.pxd
@@ -10,7 +10,7 @@ cdef class QNode(TreeNode):
     pass
 
 cdef class VNode(TreeNode):
-    pass
+    cpdef argmax(VNode self)
 
 cdef class RootVNode(VNode):
     cdef public tuple history
@@ -27,7 +27,7 @@ cdef class POUCT(Planner):
     cdef RolloutPolicy _rollout_policy
     cdef Agent _agent
     cdef int _last_num_sims
-    cdef float _last_planning_time    
+    cdef float _last_planning_time
 
     cpdef _search(self)
     cpdef _simulate(POUCT self,
@@ -44,7 +44,6 @@ cdef class RolloutPolicy(PolicyModel):
 
 cdef class RandomRollout(RolloutPolicy):
     pass
-    
+
 cdef class ActionPrior:
     cpdef get_preferred_actions(ActionPrior self, State state, tuple history)
-    
diff --git a/pomdp_py/algorithms/po_uct.pyx b/pomdp_py/algorithms/po_uct.pyx
index ab04c6bd..0822f640 100644
--- a/pomdp_py/algorithms/po_uct.pyx
+++ b/pomdp_py/algorithms/po_uct.pyx
@@ -2,7 +2,7 @@
 presented in the POMCP paper :cite:`silver2010monte` as an extension to the UCT
 algorithm :cite:`kocsis2006bandit` that combines MCTS and UCB1
 for action selection.
-    
+
 In other words, this is just POMCP without particle belief,
 but arbitrary belief representation.
 
@@ -20,7 +20,7 @@ A note on the exploration constant. Quote from :cite:`gusmao2012towards`:
 
     "This constant should reflect the agent’s prior knowledge regarding
     the amount of exploration required."
-  
+
 In the POMCP paper, they set this constant following:
 
     "The exploration constant for POMCP was set to c = Rhi - Rlo,
@@ -87,7 +87,19 @@ cdef class VNode(TreeNode):
     def print_children_value(self):
         for action in self.children:
             print("   action %s: %.3f" % (str(action), self[action].value))
-    
+
+    cpdef argmax(VNode self):
+        """argmax(VNode self)
+        Returns the action of the child with highest value"""
+        cdef Action action, best_action
+        cdef float best_value = float("-inf")
+        best_action = None
+        for action in self.children:
+            if self[action].value > best_value:
+                best_action = action
+                best_value = self[action].value
+        return best_action
+
 
 cdef class RootVNode(VNode):
     def __init__(self, num_visits, value, history):
@@ -168,7 +180,7 @@ def tree_stats_helper(root, depth, stats, max_depth=None):
             stats['total_qnodes'] += 1
             stats['total_qnodes_children'] += len(root.children)
             stats['max_qnodes_children'] = max(stats['max_qnodes_children'], len(root.children))
-        
+
         for c in root.children:
             tree_stats_helper(root[c], depth+1, stats, max_depth=max_depth)
 
@@ -193,18 +205,18 @@ cdef class ActionPrior:
         is possible that only a subset of all actions is legal;
         This is particularly true in the RockSample problem."""
         raise NotImplementedError
-    
-    
+
+
 cdef class RolloutPolicy(PolicyModel):
     cpdef Action rollout(self, State state, tuple history):
         """rollout(self, State state, tuple history=None)"""
         pass
-    
+
 cdef class RandomRollout(RolloutPolicy):
     cpdef Action rollout(self, State state, tuple history):
         """rollout(self, State state, tuple history=None)"""
         return random.sample(self.get_all_actions(state=state, history=history), 1)[0]
-    
+
 cdef class POUCT(Planner):
 
     """ POUCT (Partially Observable UCT) :cite:`silver2010monte` is presented in the POMCP
@@ -253,7 +265,7 @@ cdef class POUCT(Planner):
         # to simplify function calls; plan only for one agent at a time
         self._agent = None
         self._last_num_sims = -1
-        self._last_planning_time = -1        
+        self._last_planning_time = -1
 
     @property
     def updates_agent_belief(self):
@@ -273,14 +285,14 @@ cdef class POUCT(Planner):
         cdef Action action
         cdef float time_taken
         cdef int sims_count
-    
+
         self._agent = agent   # switch focus on planning for the given agent
         if not hasattr(self._agent, "tree"):
             self._agent.add_attr("tree", None)
         action, time_taken, sims_count = self._search()
         self._last_num_sims = sims_count
         self._last_planning_time = time_taken
-        return action            
+        return action
 
     cpdef public update(self, Agent agent, Action real_action, Observation real_observation):
         """
@@ -312,7 +324,7 @@ cdef class POUCT(Planner):
         cdef tuple preference
         cdef int num_visits_init
         cdef float value_init
-        
+
         for action in self._agent.valid_actions(state=state, history=history):
             if vnode[action] is None:
                 history_action_node = QNode(self._num_visits_init,
@@ -351,15 +363,10 @@ cdef class POUCT(Planner):
             if self._num_sims > 0\
                and sims_count >= self._num_sims:
                 break
-            
-        best_action, best_value = None, float('-inf')            
-        for action in self._agent.tree.children:
-            if self._agent.tree[action].value > best_value:
-                best_value = self._agent.tree[action].value
-                best_action = action
-                
-        return best_action, time_taken, sims_count
 
+        best_action = self._agent.tree.argmax()
+        self._agent.tree.value = self._agent.tree[best_action].value
+        return best_action, time_taken, sims_count
 
     cpdef _simulate(POUCT self,
                     State state, tuple history, VNode root, QNode parent,
@@ -397,7 +404,6 @@ cdef class POUCT(Planner):
                                                                                depth+nsteps)
         root.num_visits += 1
         root[action].num_visits += 1
-        root.value = root.value + (total_reward - root.value) / (root.num_visits)
         root[action].value = root[action].value + (total_reward - root[action].value) / (root[action].num_visits)
         return total_reward
 
@@ -408,22 +414,15 @@ cdef class POUCT(Planner):
         cdef State next_state
         cdef Observation observation
         cdef float reward
-        
+
         while depth < self._max_depth:
             action = self._rollout_policy.rollout(state, history)
             next_state, observation, reward, nsteps = sample_generative_model(self._agent, state, action)
-            if root[action] is None:
-                history_action_node = QNode(self._num_visits_init, self._value_init)
-                root[action] = history_action_node
-            if observation not in root[action]:
-                root[action][observation] = self._VNode()
-                self._expand_vnode(root[action][observation], history, state=next_state)
             history = history + ((action, observation),)
             depth += nsteps
             total_discounted_reward += reward * discount
             discount *= (self._discount_factor**nsteps)
             state = next_state
-            root = root[action][observation]
         return total_discounted_reward
 
     cpdef Action _ucb(self, VNode root):
@@ -448,8 +447,8 @@ cdef class POUCT(Planner):
         '''
         cdef State next_state
         cdef Observation observation
-        cdef float reward        
-        
+        cdef float reward
+
         if self._agent.transition_model is None:
             next_state, observation, reward = self._agent.generative_model.sample(state, action)
         else:
@@ -463,9 +462,8 @@ cdef class POUCT(Planner):
         that this function is about creating a VNode object."""
         if root:
             return RootVNode(self._num_visits_init,
-                             self._value_init,
+                             float("-inf"),
                              self._agent.history)
         else:
             return VNode(self._num_visits_init,
-                         self._value_init)
-    
+                         float("-inf"))
diff --git a/pomdp_py/framework/basics.pyx b/pomdp_py/framework/basics.pyx
index 3ec028ca..c1db04c4 100644
--- a/pomdp_py/framework/basics.pyx
+++ b/pomdp_py/framework/basics.pyx
@@ -97,13 +97,13 @@ cdef class ObservationModel:
         """
         get_all_observations(self)
         Returns a set of all possible observations, if feasible."""
-        raise NotImplementedError        
-    
+        raise NotImplementedError
+
 cdef class TransitionModel:
     """
     A TransitionModel models the distribution :math:`T(s,a,s')=\Pr(s'|s,a)`.
     """
-    
+
     def probability(self, next_state, state, action, **kwargs):
         """
         probability(self, next_state, state, action, **kwargs)
@@ -115,9 +115,9 @@ cdef class TransitionModel:
             action (~pomdp_py.framework.basics.Action): the action :math:`a`
         Returns:
             float: the probability :math:`\Pr(s'|s,a)`
-        """        
+        """
         raise NotImplementedError
-        
+
     def sample(self, state, action, **kwargs):
         """sample(self, state, action, **kwargs)
         Returns next state randomly sampled according to the
@@ -162,13 +162,14 @@ cdef class RewardModel:
             next_state (State): the next state :math:`s'`
         Returns:
             float: the probability :math:`\Pr(r|s,a,s')`
-        """        
+        """
         raise NotImplementedError
-        
+
     def sample(self, state, action, next_state, **kwargs):
         """sample(self, state, action, next_state, **kwargs)
         Returns reward randomly sampled according to the
-        distribution of this reward model.
+        distribution of this reward model. This is required,
+        i.e. assumed to be implemented for a reward model.
 
         Args:
             state (~pomdp_py.framework.basics.State): the next state :math:`s`
@@ -178,28 +179,28 @@ cdef class RewardModel:
             float: the reward :math:`r`
         """
         raise NotImplementedError
-    
+
     def argmax(self, state, action, next_state, **kwargs):
         """
         argmax(self, state, action, next_state, **kwargs)
-        Returns the most likely reward"""
+        Returns the most likely reward. This is optional."""
         raise NotImplementedError
     def get_distribution(self, state, action, next_state, **kwargs):
         """get_distribution(self, state, action, next_state, **kwargs)
         Returns the underlying distribution of the model"""
-        raise NotImplementedError    
+        raise NotImplementedError
 
 cdef class BlackboxModel:
     """
     A BlackboxModel is the generative distribution :math:`G(s,a)`
     which can generate samples where each is a tuploe :math:`(s',o,r)`.
-    """    
+    """
     def sample(self, state, action, **kwargs):
         """
         sample(self, state, action, **kwargs)
         Sample (s',o,r) ~ G(s',o,r)"""
         raise NotImplementedError
-    
+
     def argmax(self, state, action, **kwargs):
         """
         argmax(self, state, action, **kwargs)
@@ -215,7 +216,7 @@ cdef class PolicyModel:
     The reason to have a policy model is to accommodate problems
     with very large action spaces, and the available actions may vary
     depending on the state (that is, certain actions have probabilty=0)"""
-    
+
     def probability(self, action, state, **kwargs):
         """
         probability(self, action, state, **kwargs)
@@ -228,7 +229,7 @@ cdef class PolicyModel:
             float: the probability :math:`\pi(a|s)`
         """
         raise NotImplementedError
-        
+
     def sample(self, state, **kwargs):
         """sample(self, state, **kwargs)
         Returns action randomly sampled according to the
@@ -241,7 +242,7 @@ cdef class PolicyModel:
             Action: the action :math:`a`
         """
         raise NotImplementedError
-    
+
     def argmax(self, state, **kwargs):
         """
         argmax(self, state, **kwargs)
@@ -262,7 +263,7 @@ cdef class PolicyModel:
         update(self, state, next_state, action, **kwargs)
         Policy model may be updated given a (s,a,s') pair."""
         pass
-    
+
 # Belief distribution is just a distribution. There's nothing special,
 # except that the update/abstraction function can be performed over them.
 # But it would make the class hierarchy a lot more complicated if belief
@@ -294,7 +295,7 @@ cdef class Action:
     def __eq__(self, other):
         raise NotImplementedError
     def __hash__(self):
-        raise NotImplementedError        
+        raise NotImplementedError
 cdef class State:
     """
     The State class. State must be `hashable`.
@@ -302,7 +303,7 @@ cdef class State:
     def __eq__(self, other):
         raise NotImplementedError
     def __hash__(self):
-        raise NotImplementedError        
+        raise NotImplementedError
 
 cdef class Observation:
     """
@@ -311,7 +312,7 @@ cdef class Observation:
     def __eq__(self, other):
         raise NotImplementedError
     def __hash__(self):
-        raise NotImplementedError        
+        raise NotImplementedError
 
 cdef class Agent:
     """ An Agent operates in an environment by taking actions, receiving
@@ -319,7 +320,7 @@ cdef class Agent:
     planner (:class:`Planner`), and the belief update is the job taken care of
     by the belief representation or the planner. But, the Agent supplies the
     :class:`TransitionModel`, :class:`ObservationModel`, :class:`RewardModel`,
-    OR :class:`BlackboxModel` to the planner or the belief update algorithm. 
+    OR :class:`BlackboxModel` to the planner or the belief update algorithm.
 
     __init__(self, init_belief,
              policy_model,
@@ -336,7 +337,7 @@ cdef class Agent:
                  blackbox_model=None):
         self._init_belief = init_belief
         self._policy_model = policy_model
-        
+
         self._transition_model = transition_model
         self._observation_model = observation_model
         self._reward_model = reward_model
@@ -490,7 +491,7 @@ cdef class Environment:
     def reward_model(self):
         """The :class:`RewardModel` underlying the environment"""
         return self._reward_model
-    
+
     def state_transition(self, action, execute=True, **kwargs):
         """
         state_transition(self, action, execute=True, **kwargs)
@@ -536,8 +537,8 @@ cdef class Environment:
             Observation: an observation sampled from :math:`\Pr(o|s',a)`.
         """
         return observation_model.sample(self.state, action, **kwargs)
-    
-    
+
+
 cdef class Option(Action):
     """An option is a temporally abstracted action that
     is defined by (I, pi, B), where I is a initiation set,
@@ -570,10 +571,10 @@ cdef class Option(Action):
         Convenience function; Can be overriden if don't
         feel like writing a PolicyModel class"""
         return self.policy.sample(state, **kwargs)
-    
+
     def __eq__(self, other):
         raise NotImplementedError
-    
+
     def __hash__(self):
         raise NotImplementedError
 
@@ -589,7 +590,7 @@ cpdef sample_generative_model(Agent agent, State state, Action action, float dis
     """
     sample_generative_model(Agent agent, State state, Action action, float discount_factor=1.0)
     :math:`(s', o, r) \sim G(s, a)`
-    
+
     If the agent has transition/observation models, a `black box` will be created
     based on these models (i.e. :math:`s'` and :math:`o` will be sampled according
     to these models).
@@ -666,4 +667,3 @@ cpdef sample_explict_models(TransitionModel T, ObservationModel O, RewardModel R
         return next_state, observation, reward, nsteps
     else:
         return next_state, reward, nsteps
-    
diff --git a/pomdp_py/representations/distribution/histogram.pyx b/pomdp_py/representations/distribution/histogram.pyx
index a85b21ee..069b9060 100644
--- a/pomdp_py/representations/distribution/histogram.pyx
+++ b/pomdp_py/representations/distribution/histogram.pyx
@@ -61,11 +61,7 @@ cdef class Histogram(GenerativeDistribution):
         if not isinstance(other, Histogram):
             return False
         else:
-            for s in self._histogram:
-                if s not in other._histogram:
-                    return False
-                if self[s] != other[s]:
-                    return False
+            return self.histogram == other.histogram
 
     def __iter__(self):
         return iter(self._histogram)
@@ -95,10 +91,9 @@ cdef class Histogram(GenerativeDistribution):
         """get_histogram(self)
         Returns a dictionary from value to probability of the histogram"""
         return self._histogram
-    
+
     # Deprecated; it's assuming non-log probabilities
     def is_normalized(self, epsilon=1e-9):
         """Returns true if this distribution is normalized"""
         prob_sum = sum(self._histogram[state] for state in self._histogram)
         return abs(1.0-prob_sum) < epsilon
-
diff --git a/pomdp_py/utils/interfaces/conversion.py b/pomdp_py/utils/interfaces/conversion.py
new file mode 100644
index 00000000..7b7636a2
--- /dev/null
+++ b/pomdp_py/utils/interfaces/conversion.py
@@ -0,0 +1,374 @@
+"""
+Provides utility to convert a POMDP written in pomdp_py to
+a POMDP file format (.pomdp or .pomdpx). Output to a file.
+"""
+import subprocess
+import os
+import pomdp_py
+import numpy as np
+import xml.etree.ElementTree as ET
+
+def to_pomdp_file(agent, output_path=None,
+                  discount_factor=0.95):
+    """
+    Pass in an Agent, and use its components to generate
+    a .pomdp file to `output_path`.
+
+    The .pomdp file format is specified at:
+    http://www.pomdp.org/code/pomdp-file-spec.html
+
+    Note:
+    * It is assumed that the reward is independent of the observation.
+    * The state, action, and observations of the agent must be
+      explicitly enumerable.
+    * The state, action and observations of the agent must be
+      convertable to a string that does not contain any blank space.
+
+    Args:
+        agent (~pomdp_py.framework.basics.Agent): The agent
+        output_path (str): The path of the output file to write in. Optional.
+                           Default None.
+        discount_factor (float): The discount factor
+    Returns:
+        (list, list, list): The list of states, actions, observations that
+           are ordered in the same way as they are in the .pomdp file.
+    """
+    # Preamble
+    try:
+        all_states = list(agent.all_states)
+        all_actions = list(agent.all_actions)
+        all_observations = list(agent.all_observations)
+    except NotImplementedError:
+        print("S, A, O must be enumerable for a given agent to convert to .pomdp format")
+
+    content = "discount: %f\n" % discount_factor
+    content += "values: reward\n" # We only consider reward, not cost.
+
+    list_of_states = " ".join(str(s) for s in all_states)
+    assert len(list_of_states.split(" ")) == len(all_states),\
+        "states must be convertable to strings without blank spaces"
+    content += "states: %s\n" % list_of_states
+
+    list_of_actions = " ".join(str(a) for a in all_actions)
+    assert len(list_of_actions.split(" ")) == len(all_actions),\
+        "actions must be convertable to strings without blank spaces"
+    content += "actions: %s\n" % list_of_actions
+
+    list_of_observations = " ".join(str(o) for o in all_observations)
+    assert len(list_of_observations.split(" ")) == len(all_observations),\
+        "observations must be convertable to strings without blank spaces"
+    content += "observations: %s\n" % list_of_observations
+
+    # Starting belief state - they need to be normalized
+    total_belief = sum(agent.belief[s] for s in all_states)
+    content += "start: %s\n" % (" ".join(["%f" % (agent.belief[s]/total_belief)
+                                          for s in all_states]))
+
+    # State transition probabilities - they need to be normalized
+    for s in all_states:
+        for a in all_actions:
+            probs = []
+            for s_next in all_states:
+                prob = agent.transition_model.probability(s_next, s, a)
+                probs.append(prob)
+            total_prob = sum(probs)
+            for i, s_next in enumerate(all_states):
+                prob_norm = probs[i] / total_prob
+                content += 'T : %s : %s : %s %f\n' % (a, s, s_next, prob_norm)
+
+    # Observation probabilities - they need to be normalized
+    for s_next in all_states:
+        for a in all_actions:
+            probs = []
+            for o in all_observations:
+                prob = agent.observation_model.probability(o, s_next, a)
+                probs.append(prob)
+            total_prob = sum(probs)
+            assert total_prob > 0.0,\
+                "No observation is probable under state={} action={}"\
+                .format(s_next, a)
+            for i, o in enumerate(all_observations):
+                prob_norm = probs[i] / total_prob
+                content += 'O : %s : %s : %s %f\n' % (a, s_next, o, prob_norm)
+
+    # Immediate rewards
+    for s in all_states:
+        for a in all_actions:
+            for s_next in all_states:
+                # We will take the argmax reward, which works for deterministic rewards.
+                r = agent.reward_model.sample(s, a, s_next)
+                content += 'R : %s : %s : %s : *  %f\n' % (a, s, s_next, r)
+
+    if output_path is not None:
+        with open(output_path, "w") as f:
+            f.write(content)
+    return all_states, all_actions, all_observations
+
+
+
+def to_pomdpx_file(agent, pomdpconvert_path,
+                   output_path=None,
+                   discount_factor=0.95):
+    """
+    Converts an agent to a pomdpx file. Requires
+    the usage of `pomdpconvert` from github://AdaCompNUS/sarsop
+
+    Follow the instructions at https://github.com/AdaCompNUS/sarsop
+    to download and build sarsop (I tested on Ubuntu 18.04, gcc version 7.5.0)
+
+    See documentation for pomdpx at:
+    https://bigbird.comp.nus.edu.sg/pmwiki/farm/appl/index.php?n=Main.PomdpXDocumentation
+
+    First converts the agent into .pomdp, then convert it to pomdpx.
+    """
+    pomdp_path = "./temp-pomdp.pomdp"
+    to_pomdp_file(agent, pomdp_path,
+                  discount_factor=discount_factor)
+    proc = subprocess.Popen([pomdpconvert_path, pomdp_path])
+    proc.wait()
+
+    pomdpx_path = pomdp_path + "x"
+    assert os.path.exists(pomdpx_path), "POMDPx conversion failed."
+
+    with open(pomdpx_path, 'r') as f:
+        content = f.read()
+
+    if output_path is not None:
+        os.rename(pomdpx_path, output_path)
+
+    # Delete temporary files
+    os.remove(pomdp_path)
+    if os.path.exists(pomdpx_path):
+        os.remove(pomdpx_path)
+
+
+def parse_pomdp_solve_output(alpha_path, pg_path=None):
+    """Parse the output of pomdp_solve, given
+    by an .alpha file and a .pg file.
+
+    Given a path to a .alpha file, read and interpret its contents.
+    The file formats are specified at:
+    https://www.pomdp.org/code/alpha-file-spec.html
+    https://www.pomdp.org/code/pg-file-spec.html
+
+    Note on policy graph (from the official website): To use this first requires
+    knowing which of the policy graph states to start in. This can be achieved
+    by finding the alpha vector with the maximal dot product with the initial
+    starting state.
+
+    Note: Parsing the .alpha file is required. The .pg path is optional (this
+        is because I noticed some errors in the .pg file produced)
+
+    Returns:
+        alphas: [(alpha_vector, action_number) ...]
+        policy_graph: a mapping from node number to (action_number, edges)
+    """
+    alphas = []  # (alpha_vector, action_number) tuples
+    with open(alpha_path, 'r') as f:
+        action_number = None
+        alpha_vector = None
+        mode = "action"
+        for line in f:
+            line = line.rstrip()
+            if len(line) == 0:
+                continue
+
+            if mode == "action":
+                action_number = int(line)
+                mode = "alpha"
+            elif mode == "alpha":
+                alpha_vector = tuple(map(float, line.split(" ")))
+                alphas.append((alpha_vector, action_number))
+                mode = "action"
+                action_number = None
+                alpha_vector = None
+
+    policy_graph = {}  # a mapping from node number to (action_number, edges)
+    if pg_path is None:
+        return alphas
+    else:
+        with open(pg_path, 'r') as f:
+            for line in f:
+                line = line.rstrip()
+                if len(line) == 0:
+                    continue
+                parts = list(map(int, line.split()))  # Splits on whitespace
+                assert parts[0] not in policy_graph,\
+                    "The node id %d already exists. Something wrong" % parts[0]
+                policy_graph[parts[0]] = (parts[1], parts[2:])
+        return alphas, policy_graph
+
+
+
+class AlphaVectorPolicy(pomdp_py.Planner):
+    """
+    An offline POMDP policy is specified by a collection
+    of alpha vectors, each associated with an action. When
+    planning is needed, the dot product of these alpha vectors
+    and the agent's belief vector is computed and the alpha
+    vector leading to the maximum is the 'dominating' alpha
+    vector and we return its corresponding action.
+
+    An offline policy can be optionally represented as
+    a policy graph. In this case, the agent can plan without
+    actively maintaining a belief because the policy graph
+    is a finite state machine that transitions by observations.
+
+    This can be constructed using .policy file created by sarsop.
+    """
+    def __init__(self, alphas, states):
+        """
+        Args:
+            alphas (list): A list of (alpha_vector, action) tuples.
+                An alpha_vector is a list of floats
+            states (list): List of states, ordered as in .pomdp file
+        """
+        self.alphas = alphas
+        self.states = states
+
+    def plan(self, agent):
+        b = [agent.belief[s] for s in self.states]
+        _, action = max(self.alphas,
+                 key=lambda va: np.dot(b, va[0]))
+        return action
+
+    def value(self, belief):
+        """
+        Returns the value V(b) under this alpha vector policy.
+
+        V(b) = max_{alpha} alpha dot b
+        """
+        b = [belief[s] for s in self.states]
+        alpha_vector, _ = max(self.alphas,
+                              key=lambda va: np.dot(b, va[0]))
+        return np.dot(b, alpha_vector)
+
+    @classmethod
+    def construct(cls, policy_path,
+                  states, actions, solver="pomdpsol"):
+        """
+        Returns an AlphaVectorPolicy, given `alphas`,
+        which are the output of parse_appl_policy_file.
+
+        Args:
+            alphas (list): List of ( [V1, V2, ... VN], A ) tuples. A is an action number
+            actions (list): List of actions, ordered as in .pomdp file
+        """
+        if solver == "pomdp-solve" or solver == "vi":
+            return cls.construct_from_pomdp_solve(policy_path, states, actions)
+        elif solver == "pomdpsol" or solver == "sarsop":
+            alphas = []
+            root = ET.parse(policy_path).getroot()
+            for vector in root.find("AlphaVector").iter("Vector"):
+                action = actions[int(vector.attrib["action"])]
+                alpha_vector = tuple(map(float, vector.text.split()))
+                alphas.append((alpha_vector, action))
+            return AlphaVectorPolicy(alphas,
+                                     states)
+
+    @classmethod
+    def construct_from_pomdp_solve(cls, alpha_path,
+                                   states, actions):
+        alphas_with_action_numbers = parse_pomdp_solve_output(alpha_path)
+        alphas = [(alpha_vector, actions[action_number])
+                  for alpha_vector, action_number in alphas_with_action_numbers]
+        return AlphaVectorPolicy(alphas, states)
+
+
+class PGNode:
+    """A node on the policy graph"""
+    def __init__(self, node_id, alpha_vector, action):
+        self.node_id = node_id
+        self.alpha_vector = alpha_vector
+        self.action = action
+    def __eq__(self, other):
+        if isinstance(other, PolicyNode):
+            return self.node_id == other.node_id
+        return False
+    def __hash__(self):
+        return hash(self.node_id)
+    def __str__(self):
+        return repr(self)
+    def __repr__(self):
+        return "NodeID(%d)::AlphaVector(%s)::Action(%s)\n"\
+            % (self.node_id, str(self.alpha_vector), self.action)
+
+
+class PolicyGraph(pomdp_py.Planner):
+    """A PolicyGraph encodes a POMDP plan. It
+    can be constructed from the alphas and policy graph
+    format output by Cassandra's pomdp-solver."""
+
+    def __init__(self, nodes, edges, states):
+        """
+        Args:
+            nodes (list): A list of PGNodes
+            edges (dict): Mapping from node_id to a dictionary {observation -> node_id}
+            states (list): List of states, ordered as in .pomdp file
+        """
+        self.nodes = {n.node_id:n for n in nodes}
+        self.edges = edges
+        self.states = states
+        self._current_node = None
+
+    @classmethod
+    def construct(cls, alpha_path, pg_path,
+                  states, actions, observations):
+        """
+        See parse_pomdp_solve_output for detailed definitions of
+        alphas and pg.
+
+        Args:
+            alpha_path (str): Path to .alpha file
+            pg_path (str): Path to .pg file
+            states (list): List of states, ordered as in .pomdp file
+            actions (list): List of actions, ordered as in .pomdp file
+            observations (list): List of observations, ordered as in .pomdp file
+        """
+        # alphas (list): List of ( [V1, V2, ... VN], A ) tuples
+        # pg (dict): {node_id -> (A, edges)}
+        alphas, pg = parse_pomdp_solve_output(alpha_path, pg_path)
+
+        nodes = []
+        for node_id, (alpha_vector, action_number) in enumerate(alphas):
+            node = PGNode(node_id, alpha_vector, actions[action_number])
+            nodes.append(node)
+
+        edges = {}
+        for node_id in pg:
+            assert 0 <= node_id < len(nodes), "Invalid node id in policy graph"
+
+            action_number, o_links = pg[node_id]
+            assert actions[action_number] == nodes[node_id].action,\
+                "Inconsistent action mapping"
+
+            edges[node_id] = {}
+
+            for o_index, next_node_id in enumerate(o_links):
+                observation = observations[o_index]
+                edges[node_id][observation] = next_node_id
+        return PolicyGraph(nodes, edges, states)
+
+    def plan(self, agent):
+        if self._current_node is None:
+            self._current_node = self._find_node(agent)
+        return self._current_node.action
+
+    def _find_node(self, agent):
+        """Locate the node in the policy graph corresponding to the agent's current
+        belief state. """
+        b = [agent.belief[s] for s in self.states]
+        nid = max(self.nodes,
+                   key=lambda nid: np.dot(b, self.nodes[nid].alpha_vector))
+        return self.nodes[nid]
+
+    def update(self, agent, action, observation):
+        """
+        Updates the planner based on real action and observation.
+        Basically sets the current node pointer based on the incoming
+        observation."""
+        if self._current_node is None:
+            # Find out the node number using agent current belief
+            self._current_node = self._find_node(agent)
+        # Transition the current node following the graph
+        self._current_node = self.nodes[self.edges[self._current_node.node_id][observation]]
diff --git a/pomdp_py/utils/interfaces/solvers.py b/pomdp_py/utils/interfaces/solvers.py
new file mode 100644
index 00000000..a3a8950e
--- /dev/null
+++ b/pomdp_py/utils/interfaces/solvers.py
@@ -0,0 +1,117 @@
+"""
+Provides function calls to use external solvers,
+given a POMDP defined using pomdp_py interfaces.
+
+Currently, we interface with:
+* pomdp-solve: http://www.pomdp.org/code/index.html
+* SARSOP
+
+We plan to interface with:
+* POMDP.jl
+"""
+import pomdp_py
+from pomdp_py.utils.interfaces.conversion\
+    import to_pomdp_file, PolicyGraph, AlphaVectorPolicy, parse_pomdp_solve_output
+import subprocess
+import os
+
+def vi_pruning(agent, pomdp_solve_path,
+               discount_factor=0.95,
+               options=[],
+               pomdp_name="temp-pomdp",
+               remove_generated_files=False,
+               return_policy_graph=False):
+    """
+    Value Iteration with pruning, using the software pomdp-solve
+    https://www.pomdp.org/code/ developed by Anthony R. Cassandra.
+
+    Args:
+        agent (pomdp_py.Agent): The agent that contains the POMDP definition
+        pomdp_solve_path (str): Path to the `pomdp_solve` binary generated after
+            compiling the pomdp-solve library.
+        options (list): Additional options to pass in to the command line interface.
+             The options should be a list of strings, such as ["-stop_criteria", "weak", ...]
+             Some useful options are:
+                 -horizon <int>
+                 -time_limit <int>
+        pomdp_name (str): The name used to create the .pomdp file.
+        remove_generated_files (bool): True if after policy is computed,
+            the .pomdp, .alpha, .pg files are removed. Default is False.
+        return_policy_graph (bool): True if return the policy as a PolicyGraph.
+            By default is False, in which case an AlphaVectorPolicy is returned.
+    """
+    try:
+        all_states = list(agent.all_states)
+        all_actions = list(agent.all_actions)
+        all_observations = list(agent.all_observations)
+    except NotImplementedError:
+        print("S, A, O must be enumerable for a given agent to convert to .pomdp format")
+
+    pomdp_path = "./%s.pomdp" % pomdp_name
+    to_pomdp_file(agent, pomdp_path, discount_factor=discount_factor)
+    proc = subprocess.Popen([pomdp_solve_path,
+                             "-pomdp", pomdp_path,
+                             "-o", pomdp_name] + list(map(str,options)))
+    proc.wait()
+
+    # Read the value and policy graph files
+    alpha_path = "%s.alpha" % pomdp_name
+    pg_path = "%s.pg" % pomdp_name
+    if return_policy_graph:
+        policy = PolicyGraph.construct(alpha_path, pg_path,
+                                       all_states, all_actions, all_observations)
+    else:
+        policy = AlphaVectorPolicy.construct(
+            alpha_path, all_states, all_actions, solver="pomdp-solve")
+
+    # Remove temporary files
+    if remove_generated_files:
+        os.remove(pomdp_path)
+        os.remove(alpha_path)
+        os.remove(pg_path)
+    return policy
+
+
+def sarsop(agent, pomdpsol_path,
+           discount_factor=0.95,
+           timeout=30, memory=100,
+           precision=0.5,
+           pomdp_name="temp-pomdp",
+           remove_generated_files=False):
+    """
+    SARSOP, using the binary from https://github.com/AdaCompNUS/sarsop
+    This is an anytime POMDP planning algorithm
+
+    Args:
+        agent (pomdp_py.Agent): The agent that defines the POMDP models
+        pomdpsol_path (str): Path to the `pomdpsol` binary
+        timeout (int): The time limit (seconds) to run the algorithm until termination
+        memory (int): The memory size (mb) to run the algorithm until termination
+        precision (float): solver runs until regret is less than `precision`
+    """
+    try:
+        all_states = list(agent.all_states)
+        all_actions = list(agent.all_actions)
+        all_observations = list(agent.all_observations)
+    except NotImplementedError:
+        print("S, A, O must be enumerable for a given agent to convert to .pomdpx format")
+
+    pomdp_path = "./%s.pomdp" % pomdp_name
+    to_pomdp_file(agent, pomdp_path, discount_factor=discount_factor)
+    proc = subprocess.Popen([pomdpsol_path,
+                             "--timeout", str(timeout),
+                             "--memory", str(memory),
+                             "--precision", str(precision),
+                             "--output", "%s.policy" % pomdp_name,
+                             pomdp_path])
+    proc.wait()
+
+    policy_path = "%s.policy" % pomdp_name
+    policy = AlphaVectorPolicy.construct(policy_path,
+                                         all_states, all_actions)
+
+    # Remove temporary files
+    if remove_generated_files:
+        os.remove(pomdp_path)
+        os.remove(policy_path)
+    return policy
diff --git a/pomdp_py/utils/templates.py b/pomdp_py/utils/templates.py
new file mode 100644
index 00000000..a5ff3392
--- /dev/null
+++ b/pomdp_py/utils/templates.py
@@ -0,0 +1,118 @@
+"""
+Some particular implementations of the interface for convenience
+"""
+import pomdp_py
+import random
+
+class SimpleState(pomdp_py.State):
+    """A SimpleState is a state that stores
+    one piece of hashable data and the equality
+    of two states of this kind depends just on
+    this data"""
+    def __init__(self, data):
+        self.data = data
+    def __hash__(self):
+        return hash(self.data)
+    def __eq__(self, other):
+        if isinstance(other, SimpleState):
+            return self.data == other.data
+        return False
+    def __str__(self):
+        return str(self.data)
+    def __repr__(self):
+        return "SimpleState({})".format(self.data)
+
+class SimpleAction(pomdp_py.Action):
+    """A SimpleAction is an action defined by a string name"""
+    def __init__(self, name):
+        self.name = name
+    def __hash__(self):
+        return hash(self.name)
+    def __eq__(self, other):
+        if isinstance(other, SimpleAction):
+            return self.name == other.name
+        return False
+    def __str__(self):
+        return self.name
+    def __repr__(self):
+        return "SimpleAction(%s)" % self.name
+
+class SimpleObservation(pomdp_py.Observation):
+    """A SimpleObservation is an observation
+    with a piece of hashable data that defines
+    the equality."""
+    def __init__(self, data):
+        self.data = data
+    def __hash__(self):
+        return hash(self.data)
+    def __eq__(self, other):
+        if isinstance(other, SimpleObservation):
+            return self.data == other.data
+        return False
+    def __str__(self):
+        return str(self.data)
+    def __repr__(self):
+        return "SimpleObservation(%s)" % self.data
+
+class DetTransitionModel(pomdp_py.TransitionModel):
+    """A DetTransitionModel is a deterministic transition model.
+    A probability of 1 - epsilon is given for correct transition,
+    and epsilon is given for incorrect transition."""
+    def __init__(self, epsilon=1e-12):
+        self.epsilon = epsilon
+
+    def probability(self, next_state, state, action):
+        """According to problem spec, the world resets once
+        action is open-left/open-right. Otherwise, stays the same"""
+        if self.sample(state, action) == next_state:
+            return 1.0 - self.epsilon
+        else:
+            return self.epsilon
+
+    def sample(self, state, action):
+        raise NotImplementedError
+
+
+class DetObservationModel(pomdp_py.ObservationModel):
+    """A DetTransitionModel is a deterministic transition model.
+    A probability of 1 - epsilon is given for correct transition,
+    and epsilon is given for incorrect transition."""
+    def __init__(self, epsilon=1e-12):
+        self.epsilon = epsilon
+
+    def probability(self, observation, next_state, action):
+        """According to problem spec, the world resets once
+        action is open-left/open-right. Otherwise, stays the same"""
+        if self.sample(next_state, action) == observation:
+            return 1.0 - self.epsilon
+        else:
+            return self.epsilon
+
+    def sample(self, next_state, action):
+        raise NotImplementedError
+
+
+class DetRewardModel(pomdp_py.RewardModel):
+    """A DetRewardModel is a deterministic reward model (the most typical kind)."""
+    def reward_func(self, state, action, next_state):
+        raise NotImplementedError
+
+    def sample(self, state, action, next_state):
+        # deterministic
+        return self.reward_func(state, action, next_state)
+
+    def argmax(self, state, action, next_state):
+        return self.sample(state, action, next_state)
+
+class UniformPolicyModel(pomdp_py.RolloutPolicy):
+    def __init__(self, actions):
+        self.actions = actions
+
+    def sample(self, state, **kwargs):
+        return random.sample(self.actions, 1)[0]
+
+    def get_all_actions(self, state=None, history=None):
+        return self.actions
+
+    def rollout(self, state, history=None):
+        return random.sample(self.actions, 1)[0]


Test patch:
diff --git a/pomdp_py/utils/interfaces/tests/README.md b/pomdp_py/utils/interfaces/tests/README.md
new file mode 100644
index 00000000..b8f71193
--- /dev/null
+++ b/pomdp_py/utils/interfaces/tests/README.md
@@ -0,0 +1,52 @@
+# Test output
+
+01/07/2021
+
+Ubuntu 18.04, gcc 7.5.0
+
+```
+[testing] test_pomdp_file_conversion
+[testing] Running pomdp-solve on generated ./test_tiger.POMDP
+
+** Warning **
+        lp_solve reported 2 LPs with numerical instability.
+Pass.
+[testing] test_pomdpx_file_conversion
+[testing] converting to .pomdpx file
+input file   : ./temp-pomdp.pomdp
+Pass.
+
+... (solver output)
+++++++++++++++++++++++++++++++++++++++++
+User time = 0 hrs., 0 mins, 2.55 secs. (= 2.55 secs)
+System time = 0 hrs., 0 mins, 0.01 secs. (= 0.01 secs)
+Total execution time = 0 hrs., 0 mins, 2.56 secs. (= 2.56 secs)
+
+** Warning **
+        lp_solve reported 2 LPs with numerical instability.
+[testing] running tiger simulation with computed policy graph (step=0)
+[testing] running tiger simulation with computed policy graph (step=1)
+[testing] running tiger simulation with computed policy graph (step=2)
+[testing] running tiger simulation with computed policy graph (step=3)
+[testing] running tiger simulation with computed policy graph (step=4)
+[testing] running tiger simulation with computed policy graph (step=5)
+[testing] running tiger simulation with computed policy graph (step=6)
+[testing] running tiger simulation with computed policy graph (step=7)
+[testing] running tiger simulation with computed policy graph (step=8)
+[testing] running tiger simulation with computed policy graph (step=9)
+Pass.
+```
+
+
+## Download pomdp_solve
+
+```
+wget https://www.pomdp.org/code/pomdp-solve-5.4.tar.gz
+```
+Then, build the package
+according to [Installing from Source section in the documentation](https://www.pomdp.org/code/).
+
+
+## Download sarsop
+
+Follow the instructions on its [github repo](https://github.com/AdaCompNUS/sarsop).
diff --git a/pomdp_py/utils/interfaces/tests/test_conversion_pomdp-solve.py b/pomdp_py/utils/interfaces/tests/test_conversion_pomdp-solve.py
new file mode 100644
index 00000000..7e22a25d
--- /dev/null
+++ b/pomdp_py/utils/interfaces/tests/test_conversion_pomdp-solve.py
@@ -0,0 +1,43 @@
+# Tests the conversion to .pomdp file format code
+
+import sys
+import pomdp_py
+import subprocess
+from pomdp_py.utils.test_utils import *
+from pomdp_py.utils.interfaces.conversion import to_pomdp_file
+import os
+import io
+
+def test_pomdp_file_conversion(pomdp_solve_path):
+    print("[testing] test_pomdp_file_conversion")
+    tiger = make_tiger()
+
+    # Generate a .pomdp file
+    filename = "./test_tiger.POMDP"
+    to_pomdp_file(tiger.agent, filename, discount_factor=0.95)
+    assert os.path.exists(filename)
+
+    print("[testing] Running pomdp-solve on generated %s" % filename)
+    proc = subprocess.Popen([pomdp_solve_path, "-pomdp", filename],
+                            stdout=subprocess.PIPE)
+    solution_found = False
+    for line in io.TextIOWrapper(proc.stdout, encoding="utf-8"):
+        if "Solution found" in line:
+            solution_found = True
+    assert solution_found, "Something wrong - generated tiger POMDP could not be solved"
+    print("Pass.")
+
+    # Delete the generated pomdp file
+    os.remove(filename)
+
+    # Delete pomdp-solve generated files
+    remove_files("./*.pg")
+    remove_files("./*.alpha")
+
+if __name__ == "__main__":
+    if len(sys.argv) < 2:
+        print("To run test, do: %s <pomdp-solver-path>" % sys.argv[0])
+        print("Download pomdp-solve from https://www.pomdp.org/code/")
+        exit(1)
+    solver_path = sys.argv[1]
+    test_pomdp_file_conversion(solver_path)
diff --git a/pomdp_py/utils/interfaces/tests/test_conversion_pomdpx.py b/pomdp_py/utils/interfaces/tests/test_conversion_pomdpx.py
new file mode 100644
index 00000000..ed79a9db
--- /dev/null
+++ b/pomdp_py/utils/interfaces/tests/test_conversion_pomdpx.py
@@ -0,0 +1,36 @@
+# Tests the conversion to .pomdpx file format code
+
+import sys
+from pomdp_py.utils.test_utils import *
+from pomdp_py.utils.interfaces.conversion import to_pomdpx_file
+import os
+
+
+def test_pomdpx_file_conversion(pomdpconvert_path):
+    """
+    test pomdpx file conversion.
+    If pomdp file conversion works, and .pomdpx file is created,
+    then we assume the .pomdpx file is correct (because it is
+    converted using an external software).
+    """
+    print("[testing] test_pomdpx_file_conversion")
+    tiger = make_tiger()
+
+    filename = "./test_tiger.POMDPX"
+    print("[testing] converting to .pomdpx file")
+    to_pomdpx_file(tiger.agent, pomdpconvert_path,
+                   output_path=filename,
+                   discount_factor=0.95)
+    assert os.path.exists(filename), ".pomdpx file not created."
+    print("Pass.")
+
+    # Remove file
+    os.remove(filename)
+
+if __name__ == "__main__":
+    if len(sys.argv) < 2:
+        print("To run test, do: %s <pomdpconvert-path>" % sys.argv[0])
+        print("Download pomdpconvert from https://github.com/AdaCompNUS/sarsop")
+        exit(1)
+    converter_path = sys.argv[1]
+    test_pomdpx_file_conversion(converter_path)
diff --git a/pomdp_py/utils/interfaces/tests/test_local.sh b/pomdp_py/utils/interfaces/tests/test_local.sh
new file mode 100755
index 00000000..061a2a14
--- /dev/null
+++ b/pomdp_py/utils/interfaces/tests/test_local.sh
@@ -0,0 +1,7 @@
+python test_conversion_pomdp-solve.py ~/software/pomdp-solve-5.4/src/pomdp-solve
+
+python test_conversion_pomdpx.py ~/software/sarsop/src/pomdpconvert
+
+python test_vi_pruning.py ~/software/pomdp-solve-5.4/src/pomdp-solve
+
+python test_sarsop.py ~/software/sarsop/src/pomdpsol
diff --git a/pomdp_py/utils/interfaces/tests/test_sarsop.py b/pomdp_py/utils/interfaces/tests/test_sarsop.py
new file mode 100644
index 00000000..677e405f
--- /dev/null
+++ b/pomdp_py/utils/interfaces/tests/test_sarsop.py
@@ -0,0 +1,55 @@
+import sys
+import pomdp_py
+import subprocess
+from pomdp_py.utils.test_utils import *
+from pomdp_py.utils.interfaces.conversion import to_pomdp_file
+from pomdp_py.utils.interfaces.solvers import sarsop
+import os
+import io
+
+
+def test_sarsop(pomdpsol_path):
+    print("[testing] test_sarsop")
+    tiger = make_tiger()
+
+    # Building a policy graph
+    print("[testing] solving the tiger problem...")
+    policy = sarsop(tiger.agent, pomdpsol_path, discount_factor=0.95,
+                    timeout=10, memory=20, precision=0.000001,
+                    remove_generated_files=True)
+
+    assert str(policy.plan(tiger.agent)) == "listen",\
+        "Bad solution. Test failed."
+
+    # Plan with the graph for several steps. So we should get high rewards
+    # eventually in the tiger domain.
+    got_high_reward = False
+    for step in range(10):
+        true_state = tiger.env.state
+        action = policy.plan(tiger.agent)
+        observation = tiger.agent.observation_model.sample(true_state, action)
+        reward = tiger.env.reward_model.sample(true_state, action, None)
+        print("[testing] running computed policy graph"\
+              "(step=%d, action=%s, observation=%s, reward=%d)" % (step, action, observation, reward))
+
+        # belief update
+        new_belief = pomdp_py.update_histogram_belief(tiger.agent.cur_belief,
+                                                      action, observation,
+                                                      tiger.agent.observation_model,
+                                                      tiger.agent.transition_model)
+        tiger.agent.set_belief(new_belief)
+
+        assert reward == -1 or reward == 10, "Reward is negative. Failed."
+        if reward == 10:
+            got_high_reward = True
+    assert got_high_reward, "Should have gotten high reward. Failed."
+    print("Pass.")
+
+
+if __name__ == "__main__":
+    if len(sys.argv) < 2:
+        print("To run test, do: %s <pomdp-solver-path>" % sys.argv[0])
+        print("Download pomdp-solve from https://www.pomdp.org/code/")
+        exit(1)
+    solver_path = sys.argv[1]
+    test_sarsop(solver_path)
diff --git a/pomdp_py/utils/interfaces/tests/test_vi_pruning.py b/pomdp_py/utils/interfaces/tests/test_vi_pruning.py
new file mode 100644
index 00000000..534774d9
--- /dev/null
+++ b/pomdp_py/utils/interfaces/tests/test_vi_pruning.py
@@ -0,0 +1,60 @@
+import sys
+import pomdp_py
+import subprocess
+from pomdp_py.utils.test_utils import *
+from pomdp_py.utils.interfaces.conversion import to_pomdp_file
+from pomdp_py.utils.interfaces.solvers import vi_pruning
+import os
+import io
+
+
+def test_vi_pruning(pomdp_solve_path, return_policy_graph=True):
+    print("[testing] test_vi_pruning")
+    tiger = make_tiger()
+
+    # Building a policy graph
+    print("[testing] solving the tiger problem...")
+    policy = vi_pruning(tiger.agent, pomdp_solve_path, discount_factor=0.95,
+                        options=["-horizon", "100"],
+                        remove_generated_files=False,
+                        return_policy_graph=return_policy_graph)
+
+    assert str(policy.plan(tiger.agent)) == "listen",\
+        "Bad solution. Test failed."
+
+    # Plan with the graph for several steps. So we should get high rewards
+    # eventually in the tiger domain.
+    got_high_reward = False
+    for step in range(10):
+        true_state = tiger.env.state
+        action = policy.plan(tiger.agent)
+        observation = tiger.agent.observation_model.sample(true_state, action)
+        reward = tiger.env.reward_model.sample(true_state, action, None)
+        print("[testing] simulating computed policy graph"\
+              "(step=%d, action=%s, observation=%s, reward=%d)" % (step, action, observation, reward))
+        # No belief update needed. Just update the policy graph
+        if return_policy_graph:
+            # We use policy graph. No belief update is needed. Just update the policy.
+            policy.update(tiger.agent, action, observation)
+        else:
+            # belief update is needed
+            new_belief = pomdp_py.update_histogram_belief(tiger.agent.cur_belief,
+                                                          action, observation,
+                                                          tiger.agent.observation_model,
+                                                          tiger.agent.transition_model)
+            tiger.agent.set_belief(new_belief)
+
+        assert reward == -1 or reward == 10, "Reward is negative. Failed."
+        if reward == 10:
+            got_high_reward = True
+    assert got_high_reward, "Should have gotten high reward. Failed."
+    print("Pass.")
+
+if __name__ == "__main__":
+    if len(sys.argv) < 2:
+        print("To run test, do: %s <pomdp-solver-path>" % sys.argv[0])
+        print("Download pomdp-solve from https://www.pomdp.org/code/")
+        exit(1)
+    solver_path = sys.argv[1]
+    # test_vi_pruning(solver_path, return_policy_graph=True)
+    test_vi_pruning(solver_path, return_policy_graph=False)
diff --git a/pomdp_py/utils/test_utils.py b/pomdp_py/utils/test_utils.py
new file mode 100644
index 00000000..97051c73
--- /dev/null
+++ b/pomdp_py/utils/test_utils.py
@@ -0,0 +1,18 @@
+import pomdp_py
+from pomdp_problems.tiger.tiger_problem import TigerProblem
+from pomdp_problems.tiger.tiger_problem import State as TigerState
+import os
+import glob
+
+def remove_files(pattern):
+    file_list = glob.glob(pattern)
+    for file_path in file_list:
+        os.remove(file_path)
+
+def make_tiger(noise=0.15, init_state="tiger-left", init_belief=[0.5, 0.5]):
+    """Convenient function to quickly build a tiger domain.
+    Useful for testing"""
+    tiger = TigerProblem(noise, TigerState(init_state),
+                         pomdp_py.Histogram({TigerState("tiger-left"): init_belief[0],
+                                             TigerState("tiger-right"): init_belief[1]}))
+    return tiger