finished implementing ProductDomain

Author: ihincks

src/qinfer/domains.py

@@ -34,9 +34,9 @@
 from operator import mul
 from scipy.special import binom
 from math import factorial
-from itertools import combinations_with_replacement
+from itertools import combinations_with_replacement, product
 import numpy as np
-from numpy.lib.recfunctions import merge_arrays
+from .utils import join_struct_arrays, separate_struct_array
 
 import abc
 
@@ -152,14 +152,24 @@ class ProductDomain(Domain):
     """
     A domain made from the cartesian product of other domains.
 
-    :param list domains: A list of domains.
+    :param Domain domains: ``Domain`` objects as separate arguments, 
+        or as a singe list of ``Domain``s.
     """
-    def __init__(self, domains):
-        super(ProductDomain, self).__init__()
+    def __init__(self, *domains):
+        
+        if len(domains) == 1:
+            try:
+                self._domains = list(domains[0])
+            except:
+                self._domains = domains
+        else:
+            self._domains = domains
+        
         self._domains = domains
-        self._example_point = merge_arrays(
-            [domain.example_point for domain in self._domains],
-            flatten = True, usemask = False)
+        self._dtypes = [domain.example_point.dtype for domain in self._domains]
+        self._example_point = join_struct_arrays(
+            [np.array(domain.example_point) for domain in self._domains]
+        )
         self._dtype = self._example_point.dtype
 
     @property
@@ -224,12 +234,49 @@ def values(self):
         :rtype: `np.ndarray`
         """
         if self.is_finite:
-            raise NotImplemented()
+            separate_values = [domain.values for domain in self._domains]
+            return np.concatenate([
+                join_struct_arrays(map(np.array, value)) 
+                for value in product(*separate_values)
+            ])
         else:
             return self.example_point
 
+    ## METHODS ##
+    
+    def _mytype(self, array):
+        # astype does weird stuff with struct names, and possibly 
+        # depends on numpy version; hopefully 
+        # the following is a bit more predictable since it passes through 
+        # uint8
+        return separate_struct_array(array, self.dtype)[0]
+    
+    def to_regular_arrays(self, array):
+        """
+        Expands from an array of type `self.dtype` into a list of
+        arrays with dtypes corresponding to the factor domains.
 
-    ## ABSTRACT METHODS ##
+        :param np.ndarray array: An `np.array` of type `self.dtype`.
+
+        :rtype: ``list``
+        """
+        return separate_struct_array(self._mytype(array), self._dtypes)
+
+    def from_regular_arrays(self, arrays):
+        """
+        Merges a list of arrays (of the same shape) of dtypes 
+        corresponding to the factor domains into a single array 
+        with the dtype of the ``ProductDomain``.
+
+        :param list array: A list of ``np.ndarray``s
+
+        :rtype: `np.ndarray`
+        """
+        return self._mytype(join_struct_arrays([
+            array.astype(dtype)
+            for dtype, array in zip(self._dtypes, arrays)
+        ]))
+        
 
     def in_domain(self, points):
         """
@@ -240,7 +287,11 @@ def in_domain(self, points):
 
         :rtype: `bool`
         """
-        raise NotImplemented()
+        return all([
+            domain.in_domain(array)
+            for domain, array in
+                zip(self._domains, separate_struct_array(points, self._dtypes))
+        ])
 
 
 ## CLASSES ###################################################################

src/qinfer/utils.py

@@ -454,6 +454,48 @@ def pretty_time(secs, force_h=False, force_m=False):
 def safe_shape(arr, idx=0, default=1):
     shape = np.shape(arr)
     return shape[idx] if idx < len(shape) else default
+    
+def join_struct_arrays(arrays):
+    """
+    Takes a list of possibly structured arrays, concatenates their
+    dtypes, and returns one big array with that dtype. Does the 
+    inverse of ``separate_struct_array``.
+    
+    :param list arrays: List of ``np.ndarray``s
+    """
+    # taken from http://stackoverflow.com/questions/5355744/numpy-joining-structured-arrays
+    sizes = np.array([a.itemsize for a in arrays])
+    offsets = np.r_[0, sizes.cumsum()]
+    shape = arrays[0].shape
+    joint = np.empty(shape + (offsets[-1],), dtype=np.uint8)
+    for a, size, offset in zip(arrays, sizes, offsets):
+        joint[...,offset:offset+size] = np.atleast_1d(a).view(np.uint8).reshape(shape + (size,))
+    dtype = sum((a.dtype.descr for a in arrays), [])
+    return joint.ravel().view(dtype)
+    
+def separate_struct_array(array, dtypes):
+    """
+    Takes an array with a structured dtype, and separates it out into 
+    a list of arrays with dtypes coming from the input ``dtypes``.
+    Does the inverse of ``join_struct_arrays``.
+    
+    :param np.ndarray array: Structured array.
+    :param dtypes: List of ``np.dtype``, or just a ``np.dtype`` and the number of
+        them is figured out automatically by counting bytes.
+    """
+    try:
+        offsets = np.cumsum([np.dtype(dtype).itemsize for dtype in dtypes])
+    except TypeError:
+        dtype_size = np.dtype(dtypes).itemsize
+        num_fields = int(array.nbytes / (array.size * dtype_size))
+        offsets = np.cumsum([dtype_size] * num_fields)
+        dtypes = [dtypes] * num_fields
+    offsets = np.concatenate([[0], offsets]).astype(int)
+    uint_array = array.view(np.uint8).reshape(array.shape + (-1,))
+    return [
+        uint_array[..., offsets[idx]:offsets[idx+1]].flatten().view(dtype)
+        for idx, dtype in enumerate(dtypes)
+    ]
 
 
 #==============================================================================