Assemble a set of unique magnitude tuples for the ROI
975510
def mask_roi_unique(self):
        
        # There is no good inherent way in numpy to do this...
        # http://stackoverflow.com/q/16970982/

        # Also possible and simple:
        #return np.unique(zip(self.mask_1.mask_roi_sparse,self.mask_2.mask_roi_sparse))

        A = np.vstack([self.mask_1.mask_roi_sparse,self.mask_2.mask_roi_sparse]).T
        B = A[np.lexsort(A.T[::-1])]
        return B[np.concatenate(([True],np.any(B[1:]!=B[:-1],axis=1)))]
1076302
def posterior_to_xarray(self):
        """"""
        data = {}
        for idx, var_name in enumerate(self.var_names):
            # Use emcee3 syntax, else use emcee2
            data[var_name] = (
                self.sampler.get_chain()[(..., idx)].T
                if hasattr(self.sampler, "get_chain")
                else self.sampler.chain[(..., idx)]
            )
        return dict_to_dataset(data, library=self.emcee, coords=self.coords, dims=self.dims)
44742
def _filtdim(items, shape, dim, nsl):
    """"""
    normshape = tuple(stop - start for start, stop in shape)
    nsl_type = type(nsl)
    newitems = list()
    # Number of groups
    num = reduce(operator.mul, normshape[:dim+1])
    # Size of each group
    size = len(items) // num
    # Size of the dimension
    n = normshape[dim]
    if nsl_type is int:
        for i in range(num):
            if i % n == nsl:
                newitems += items[size*i:size*(i+1)]
        # Collapse dimension
        newshape = shape[:dim] + shape[dim+1:]
    elif nsl_type is slice:
        for i in range(num):
            if nsl.start <= (i % n) < nsl.stop:
                newitems += items[size*i:size*(i+1)]
        # Reshape dimension
        offset = shape[dim][0]
        redim = (offset + nsl.start, offset + nsl.stop)
        newshape = shape[:dim] + (redim, ) + shape[dim+1:]
    # farray
    else:
        if nsl.size < clog2(n):
            fstr = "expected dim {} select to have >= {} bits, got {}"
            raise ValueError(fstr.format(dim, clog2(n), nsl.size))
        groups = [list() for _ in range(n)]
        for i in range(num):
            groups[i % n] += items[size*i:size*(i+1)]
        for muxins in zip(*groups):
            it = boolfunc.iter_terms(nsl._items)
            xs = [reduce(operator.and_, (muxin, ) + next(it))
                  for muxin in muxins]
            newitems.append(reduce(operator.or_, xs))
        # Collapse dimension
        newshape = shape[:dim] + shape[dim+1:]
    return newitems, newshape
980011
def make_model_rois(self, components, name_factory):
        
        ret_dict = {}

        # Figure out which sources need to be split by components
        master_roi_source_info = {}
        sub_comp_sources = {}
        for comp_name, model_comp in self.model_components.items():
            comp_info = model_comp.info
            if comp_info.components is None:
                master_roi_source_info[comp_name] = model_comp
            else:
                sub_comp_sources[comp_name] = model_comp

        # Build the xml for the master
        master_roi = SourceFactory.make_roi(master_roi_source_info)
        master_xml_mdl = name_factory.master_srcmdl_xml(
            modelkey=self.model_name)
        print("Writing master ROI model to %s" % master_xml_mdl)
        master_roi.write_xml(master_xml_mdl)
        ret_dict['master'] = master_roi

        # Now deal with the components
        for comp in components:
            zcut = "zmax%i" % comp.zmax
            compkey = "%s_%s" % (zcut, comp.make_key(
                '{ebin_name}_{evtype_name}'))
            # name_keys = dict(zcut=zcut,
            #                 modelkey=self.model_name,
            #                 component=compkey)
            comp_roi_source_info = {}
            for comp_name, model_comp in sub_comp_sources.items():
                comp_info = model_comp.info
                if comp_info.selection_dependent:
                    key = comp.make_key('{ebin_name}_{evtype_name}')
                elif comp_info.moving:
                    key = zcut
                info_clone = comp_info.components[key].clone_and_merge_sub(key)
                comp_roi_source_info[comp_name] =\
                    ModelComponent(info=info_clone,
                                   spectrum=model_comp.spectrum)

            # Build the xml for the component
            comp_roi = SourceFactory.make_roi(comp_roi_source_info)
            comp_xml_mdl = name_factory.comp_srcmdl_xml(modelkey=self.model_name,
                                                        component=compkey)
            print("Writing component ROI model to %s" % comp_xml_mdl)
            comp_roi.write_xml(comp_xml_mdl)
            ret_dict[compkey] = comp_roi

        return ret_dict
914743
def squeeze(attrs, inputs, proto_obj):
    """"""
    new_attrs = translation_utils._fix_attribute_names(attrs,
                                                       {'axes' : 'axis'})
    return 'squeeze', new_attrs, inputs
320114
def mask_roi_digi(self):
        
        # http://stackoverflow.com/q/24205045/#24206440
        A = np.vstack([self.mask_1.mask_roi_sparse,self.mask_2.mask_roi_sparse]).T
        B = self.mask_roi_unique

        AA = np.ascontiguousarray(A)
        BB = np.ascontiguousarray(B)
         
        dt = np.dtype((np.void, AA.dtype.itemsize * AA.shape[1]))
        a = AA.view(dt).ravel()
        b = BB.view(dt).ravel()
         
        idx = np.argsort(b)
        indices = np.searchsorted(b[idx],a)
        return idx[indices]
105619
def _bg_combine(self, bgs):
        """"""
        out = np.ones(self.h5["raw"].shape, dtype=float)
        # bg is an h5py.DataSet
        for bg in bgs:
            out *= bg[:]
        return out
1303058
def fan_triangulate(indices):
    """"""
    if len(indices[0]) != 4:
        raise ValueError("Assumes working with a sequence of quad indices")
    new_indices = []
    for face in indices:
        new_indices.extend([face[:-1], face[1:]])
    return np.array(new_indices)
982125
def _make_clusters(self, matrix, num_clusters_per_roi, metric):
        """"""

        from scipy.cluster.hierarchy import fclusterdata

        # maxclust needed to ensure t is interpreted as # clusters in heirarchical clustering
        group_ids = fclusterdata(matrix, metric=metric, t=num_clusters_per_roi,
                                 criterion='maxclust')
        group_set = np.unique(group_ids)
        clusters = [
            self._summary_func(matrix[group_ids == group, :], axis=0, keepdims=True)
            for group in group_set]

        return np.vstack(clusters).squeeze()
305879
def posterior_to_xarray(self):
        """"""
        data = {}
        for i, var_name in enumerate(self.var_names):
            data[var_name] = np.expand_dims(self.posterior[i], axis=0)
        return dict_to_dataset(data, library=self.tfp, coords=self.coords, dims=self.dims)
198867
def build_quantity(orig_text, text, item, values, unit, surface, span, uncert):
    """"""
    # Discard irrelevant txt2float extractions, cardinal numbers, codes etc.
    if surface.lower() in ['a', 'an', 'one'] or \
            re.search(r'1st|2nd|3rd|[04-9]th', surface) or \
            re.search(r'\d+[A-Z]+\d+', surface) or \
            re.search(r'\ba second\b', surface, re.IGNORECASE):
        logging.debug(u'\tMeaningless quantity ("%s"), discard', surface)
        return

    # Usually "$3T" does not stand for "dollar tesla"
    elif unit.entity.dimensions and \
            unit.entity.dimensions[0]['base'] == 'currency':
        if len(unit.dimensions) > 1:
            try:
                suffix = re.findall(r'\d(K|M|B|T)\b(.*?)$', surface)[0]
                values = [i * r.SUFFIXES[suffix[0]] for i in values]
                unit = l.UNITS[unit.dimensions[0]['base']][0]
                if suffix[1]:
                    surface = surface[:surface.find(suffix[1])]
                    span = (span[0], span[1] - len(suffix[1]))
                logging.debug(u'\tCorrect for "$3T" pattern')
            except IndexError:
                pass
        else:
            try:
                suffix = re.findall(r'%s(K|M|B|T)\b' % re.escape(surface),
                                    orig_text)[0]
                surface += suffix
                span = (span[0], span[1] + 1)
                values = [i * r.SUFFIXES[suffix] for i in values]
                logging.debug(u'\tCorrect for "$3T" pattern')
            except IndexError:
                pass

    # Usually "1990s" stands for the decade, not the amount of seconds
    elif re.match(r'[1-2]\d\d0s', surface):
        unit = l.NAMES['dimensionless']
        surface = surface[:-1]
        span = (span[0], span[1] - 1)
        logging.debug(u'\tCorrect for decade pattern')

    # Usually "in" stands for the preposition, not inches
    elif unit.dimensions[-1]['base'] == 'inch' and \
            re.search(r' in$', surface) and '/' not in surface:
        if len(unit.dimensions) > 1:
            unit = get_unit_from_dimensions(unit.dimensions[:-1], orig_text)
        else:
            unit = l.NAMES['dimensionless']
        surface = surface[:-3]
        span = (span[0], span[1] - 3)
        logging.debug(u'\tCorrect for "in" pattern')

    elif is_quote_artifact(text, item.span()):
        if len(unit.dimensions) > 1:
            unit = get_unit_from_dimensions(unit.dimensions[:-1], orig_text)
        else:
            unit = l.NAMES['dimensionless']
        surface = surface[:-1]
        span = (span[0], span[1] - 1)
        logging.debug(u'\tCorrect for quotes')

    elif re.search(r' time$', surface) and len(unit.dimensions) > 1 and \
            unit.dimensions[-1]['base'] == 'count':
        unit = get_unit_from_dimensions(unit.dimensions[:-1], orig_text)
        surface = surface[:-5]
        span = (span[0], span[1] - 5)
        logging.debug(u'\tCorrect for "time"')

    objs = []
    for value in values:
        obj = c.Quantity(value=value,
                         unit=unit,
                         surface=surface,
                         span=span,
                         uncertainty=uncert)
        objs.append(obj)

    return objs
781075
def get_inputs_from_cm(index, cm):
    """"""
    return tuple(i for i in range(cm.shape[0]) if cm[i][index])
20994
def make_roi_plots(self, gta, mcube_tot, **kwargs):
        

        fmt = kwargs.get('format', self.config['format'])
        figsize = kwargs.get('figsize', self.config['figsize'])
        prefix = kwargs.get('prefix', '')
        loge_bounds = kwargs.get('loge_bounds', None)
        weighted = kwargs.get('weighted', False)

        roi_kwargs = {}
        roi_kwargs.setdefault('loge_bounds', loge_bounds)
        roi_kwargs.setdefault(
            'graticule_radii', self.config['graticule_radii'])
        roi_kwargs.setdefault('label_ts_threshold',
                              self.config['label_ts_threshold'])
        roi_kwargs.setdefault('cmap', self.config['cmap'])
        roi_kwargs.setdefault('catalogs', self._catalogs)

        if loge_bounds is None:
            loge_bounds = (gta.log_energies[0], gta.log_energies[-1])
        esuffix = '_%.3f_%.3f' % (loge_bounds[0], loge_bounds[1])

        mcube_diffuse = gta.model_counts_map('diffuse')
        counts_map = gta.counts_map()

        if weighted:
            wmap = gta.weight_map()
            counts_map = copy.deepcopy(counts_map)
            mcube_tot = copy.deepcopy(mcube_tot)
            counts_map.data *= wmap.data
            mcube_tot.data *= wmap.data
            mcube_diffuse.data *= wmap.data

        # colors = ['k', 'b', 'g', 'r']
        data_style = {'marker': 's', 'linestyle': 'None'}

        fig = plt.figure(figsize=figsize)

        if gta.projtype == "WCS":
            xmin = -1
            xmax = 1
        elif gta.projtype == "HPX":
            hpx2wcs = counts_map.make_wcs_mapping(proj='CAR', oversample=2)
            counts_map = counts_map.to_wcs(hpx2wcs=hpx2wcs)
            mcube_tot = mcube_tot.to_wcs(hpx2wcs=hpx2wcs)
            mcube_diffuse = mcube_diffuse.to_wcs(hpx2wcs=hpx2wcs)
            xmin = None
            xmax = None

        fig = plt.figure(figsize=figsize)
        rp = ROIPlotter(mcube_tot, roi=gta.roi, **roi_kwargs)
        rp.plot(cb_label='Counts', zscale='pow', gamma=1. / 3.)
        plt.savefig(os.path.join(gta.config['fileio']['workdir'],
                                 '%s_model_map%s.%s' % (
                                     prefix, esuffix, fmt)))
        plt.close(fig)

        rp = ROIPlotter(counts_map, roi=gta.roi, **roi_kwargs)

        rp.plot(cb_label='Counts', zscale='sqrt')
        plt.savefig(os.path.join(gta.config['fileio']['workdir'],
                                 '%s_counts_map%s.%s' % (
                                     prefix, esuffix, fmt)))
        plt.close(fig)

        for iaxis, xlabel, psuffix in zip([0, 1],
                                          ['LON Offset [deg]', 'LAT Offset [deg]'],
                                          ['xproj', 'yproj']):

            fig = plt.figure(figsize=figsize)
            rp.plot_projection(iaxis, label='Data', color='k',
                               xmin=xmin, xmax=xmax, **data_style)
            rp.plot_projection(iaxis, data=mcube_tot, label='Model', xmin=xmin, xmax=xmax,
                               noerror=True)
            rp.plot_projection(iaxis, data=mcube_diffuse, label='Diffuse', xmin=xmin, xmax=xmax,
                               noerror=True)
            plt.gca().set_ylabel('Counts')
            plt.gca().set_xlabel(xlabel)
            plt.gca().legend(frameon=False)
            annotate(loge_bounds=loge_bounds)
            plt.savefig(os.path.join(gta.config['fileio']['workdir'],
                                     '%s_counts_map_%s%s.%s' % (prefix, psuffix,
                                                                esuffix, fmt)))
            plt.close(fig)
1234926
def construct_zernike_lookuptable(zernike_indexes):
    
    n_max = np.max(zernike_indexes[:,0])
    factorial = np.ones((1 + n_max,), dtype=float)
    factorial[1:] = np.cumproduct(np.arange(1, 1 + n_max, dtype=float))
    width = int(n_max//2 + 1)
    lut = np.zeros((zernike_indexes.shape[0],width), dtype=float)
    for idx, (n, m) in enumerate(zernike_indexes):
        alt = 1
        npmh = (n+m)//2
        nmmh = (n-m)//2
        for k in range(0,nmmh+1):
            lut[idx,k] = \
                (alt * factorial[n-k] /
                 (factorial[k]*factorial[npmh-k]*factorial[nmmh-k]))
            alt = -alt
    return lut
845747
def get_image(roi_rec, short, max_size, mean, std):
    
    im = imdecode(roi_rec['image'])
    if roi_rec["flipped"]:
        im = im[:, ::-1, :]
    im, im_scale = resize(im, short, max_size)
    height, width = im.shape[:2]
    im_info = np.array([height, width, im_scale], dtype=np.float32)
    im_tensor = transform(im, mean, std)

    # gt boxes: (x1, y1, x2, y2, cls)
    if roi_rec['gt_classes'].size > 0:
        gt_inds = np.where(roi_rec['gt_classes'] != 0)[0]
        gt_boxes = np.empty((len(gt_inds), 5), dtype=np.float32)
        gt_boxes[:, 0:4] = roi_rec['boxes'][gt_inds, :]
        gt_boxes[:, 4] = roi_rec['gt_classes'][gt_inds]
        # scale gt_boxes
        gt_boxes[:, 0:4] *= im_scale
    else:
        gt_boxes = np.empty((0, 5), dtype=np.float32)

    return im_tensor, im_info, gt_boxes
17997
def get_bin_indices(self, values):
        """"""
        return tuple([self.get_axis_bin_index(values[ax_i], ax_i)
                      for ax_i in range(self.dimensions)])
178210
def tdot_blas(mat, out=None):
    """"""
    if (mat.dtype != 'float64') or (len(mat.shape) != 2):
        return np.dot(mat, mat.T)
    nn = mat.shape[0]
    if out is None:
        out = np.zeros((nn, nn))
    else:
        assert(out.dtype == 'float64')
        assert(out.shape == (nn, nn))
        # FIXME: should allow non-contiguous out, and copy output into it:
        assert(8 in out.strides)
        # zeroing needed because of dumb way I copy across triangular answer
        out[:] = 0.0

    # # Call to DSYRK from BLAS
    mat = np.asfortranarray(mat)
    out = blas.dsyrk(alpha=1.0, a=mat, beta=0.0, c=out, overwrite_c=1,
                     trans=0, lower=0)

    symmetrify(out, upper=True)
    return np.ascontiguousarray(out)