Merge branch 'master' of https://github.com/sbailey/specter

py/extract/ex2d.py

@@ -5,6 +5,7 @@
 import sys
 import numpy as np
 import scipy.sparse
+import scipy.linalg
 from scipy.sparse import spdiags, issparse
 from scipy.sparse.linalg import spsolve
 
@@ -103,8 +104,11 @@ def ex2d(image, ivar, psf, specrange, wavelengths, xyrange=None,
         print "Dumping {} for debugging".format(outfile)
         import fitsio
         fitsio.write(outfile, image, clobber=True)
-        fitsio.write(outfile, ivar)
-        fitsio.write(outfile, A.data)
+        fitsio.write(outfile, ivar, extname='IVAR')
+        fitsio.write(outfile, A.data, extname='ADATA') 
+        fitsio.write(outfile, A.indices, extname='AINDICES')
+        fitsio.write(outfile, A.indptr, extname='AINDPTR')
+        fitsio.write(outfile, iCov.toarray(), extname='ICOV')
         raise err
 
     #- Convolve with Resolution matrix to decorrelate errors
@@ -133,8 +137,8 @@ def sym_sqrt(a):
 
     WRITTEN: Adam S. Bolton, U. of Utah, 2009
     """
-
-    w, v = np.linalg.eigh(a)
+    ### w, v = np.linalg.eigh(a)  #- np.linalg.eigh is single precision !?!
+    w, v = scipy.linalg.eigh(a)
 
     #- Trim meaningless eigenvalues below machine precision
     ibad = w < w.max()*sys.float_info.epsilon
@@ -166,9 +170,12 @@ def resolution_from_icov(icov):
 
     WRITTEN: Adam S. Bolton, U. of Utah, 2009
     """
+    #- force symmetry since due to rounding it might not be exactly symmetric
+    icov = 0.5*(icov + icov.T)
+
     if issparse(icov):
         icov = icov.toarray()
-        
+
     sqrt_icov = sym_sqrt(icov)
     norm_vector = np.sum(sqrt_icov, axis=1)
     R = np.outer(norm_vector**(-1), np.ones(norm_vector.size)) * sqrt_icov

py/util/cachedict.py

@@ -3,20 +3,39 @@
 """
 
 class CacheDict(dict):
-    def __init__(self, n):
+    """
+    A dictionary that keeps only the last n items added
+    """
+    def __init__(self, n, d=None):
+        """
+        Create CacheDict with n keys cached.
+
+        If d is a dict, use that to initialize the key/value pairs
+        """
         self._keys = [None,]*n
         self._current = -1
         self._n = n
+
+        if type(d) == dict:
+            for key, value in d.items():
+                self[key] = value
+
+    #- Needed by pickle to reconstruct the object.
+    #- Pickle tries to reconstruct the dictionary via __setitem__ before
+    #- filling in _keys, _current, _n.  This allows it to create a new
+    #- object first to get those right, before calling __setitem__
+    def __reduce__(self):
+        return type(self), (self._n, dict(self))
 
     def __setitem__(self, key, value):
+        """Sets the key/value, possibly dropping an earlier cached key/value"""
         if key in self:
             return
-            
+
         i = self._current = (self._current + 1) % self._n
         if self._keys[i] is not None:
             del self[self._keys[i]]
 
         self._keys[i] = key
         dict.__setitem__(self, key, value)
 
-