diff --git a/py/extract/ex2d.py b/py/extract/ex2d.py index e5ca710..31a31ac 100644 --- a/py/extract/ex2d.py +++ b/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 diff --git a/py/util/cachedict.py b/py/util/cachedict.py index a12b6b1..236e813 100644 --- a/py/util/cachedict.py +++ b/py/util/cachedict.py @@ -3,15 +3,35 @@ """ 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]] @@ -19,4 +39,3 @@ def __setitem__(self, key, value): self._keys[i] = key dict.__setitem__(self, key, value) - \ No newline at end of file