Calculate the spatial signal probability for each catalog object.

        Parameters:
        -----------
        None

        Returns:
        --------
        u_spatial : array of spatial probabilities per object
345730
def calc_signal_spatial(self):
        
        # Calculate the surface intensity
        self.surface_intensity_sparse = self.calc_surface_intensity()

        # Calculate the probability per object-by-object level
        self.surface_intensity_object = self.kernel.pdf(self.catalog.lon,
                                                        self.catalog.lat)
        
        # Spatial component of signal probability
        u_spatial = self.surface_intensity_object
        return u_spatial
1239636
def _photometricErrors(self, catalog=None, n_per_bin=100):
        

        if catalog is None:
            # Simple proxy for photometric errors
            release = self.config['data']['release']
            band_1 = self.config['catalog'].get('mag_1_band')
            if not band_1: band_1 = self.config['isochrone']['mag_1_field']
            band_2 = self.config['catalog'].get('mag_2_band')
            if not band_2: band_2 = self.config['isochrone']['mag_2_field']
            
            DELMIN = 0.0
            pars_1 = MAGERR_PARAMS[release][band_1]
            
            def photo_err_1(delta):
                p = pars_1
                return np.clip(np.exp(p[0]*delta+p[1])+p[2], 0, np.exp(p[0]*(DELMIN)+p[1])+p[2])

            pars_2 = MAGERR_PARAMS[release][band_2]
            def photo_err_2(delta):
                p = pars_2
                return np.clip(np.exp(p[0]*delta+p[1])+p[2], 0, np.exp(p[0]*(DELMIN)+p[1])+p[2])

        else:
            catalog.spatialBin(self.roi)

            if len(catalog.mag_1) < n_per_bin:
                logger.warning("Catalog contains fewer objects than requested to calculate errors.")
                #n_per_bin = int(len(catalog.mag_1) / 3)
                return self._photometricErrors(catalog=None)
             
            # Band 1
            mag_1_thresh = self.mask_1.mask_roi_sparse[catalog.pixel_roi_index] - catalog.mag_1
            sorting_indices = np.argsort(mag_1_thresh)
            mag_1_thresh_sort = mag_1_thresh[sorting_indices]
            mag_err_1_sort = catalog.mag_err_1[sorting_indices]
             
            # ADW: Can't this be done with np.median(axis=?)
            mag_1_thresh_medians = []
            mag_err_1_medians = []
            for i in range(0, int(len(mag_1_thresh) / float(n_per_bin))):
                mag_1_thresh_medians.append(np.median(mag_1_thresh_sort[n_per_bin * i: n_per_bin * (i + 1)]))
                mag_err_1_medians.append(np.median(mag_err_1_sort[n_per_bin * i: n_per_bin * (i + 1)]))
             
            if mag_1_thresh_medians[0] > 0.:
                mag_1_thresh_medians = np.insert(mag_1_thresh_medians, 0, -99.)
                mag_err_1_medians = np.insert(mag_err_1_medians, 0, mag_err_1_medians[0])
             
            photo_err_1 = scipy.interpolate.interp1d(mag_1_thresh_medians, mag_err_1_medians,
                                                     bounds_error=False, fill_value=mag_err_1_medians[-1])
             
            # Band 2
            mag_2_thresh = self.mask_2.mask_roi_sparse[catalog.pixel_roi_index] - catalog.mag_2
            sorting_indices = np.argsort(mag_2_thresh)
            mag_2_thresh_sort = mag_2_thresh[sorting_indices]
            mag_err_2_sort = catalog.mag_err_2[sorting_indices]
             
            mag_2_thresh_medians = []
            mag_err_2_medians = []
            for i in range(0, int(len(mag_2_thresh) / float(n_per_bin))):
                mag_2_thresh_medians.append(np.median(mag_2_thresh_sort[n_per_bin * i: n_per_bin * (i + 1)]))
                mag_err_2_medians.append(np.median(mag_err_2_sort[n_per_bin * i: n_per_bin * (i + 1)]))
             
            if mag_2_thresh_medians[0] > 0.:
                mag_2_thresh_medians = np.insert(mag_2_thresh_medians, 0, -99.)
                mag_err_2_medians = np.insert(mag_err_2_medians, 0, mag_err_2_medians[0])
             
            photo_err_2 = scipy.interpolate.interp1d(mag_2_thresh_medians, mag_err_2_medians,
                                                     bounds_error=False, fill_value=mag_err_2_medians[-1])
             
        self.photo_err_1=photo_err_1
        self.photo_err_2=photo_err_2

        return self.photo_err_1, self.photo_err_2
767366
def __d2f(disc_coords, v):
        
        return np.array([disc_coords[0] / v.shape[0],
                         disc_coords[1] / v.shape[1],
                         disc_coords[2] / v.shape[2]])
518376
def inference_likelihood(self, kern, X, Z, likelihood, Y):
        

        num_data, output_dim = Y.shape
        input_dim = Z.shape[0]
        if self.mpi_comm != None:
            from mpi4py import MPI
            num_data_all = np.array(num_data,dtype=np.int32)
            self.mpi_comm.Allreduce([np.int32(num_data), MPI.INT], [num_data_all, MPI.INT])
            num_data = num_data_all

        if isinstance(X, VariationalPosterior):
            uncertain_inputs = True
        else:
            uncertain_inputs = False

        #see whether we've got a different noise variance for each datum
        beta = 1./np.fmax(likelihood.variance, 1e-6)
        het_noise = beta.size > 1
        if het_noise:
            self.batchsize = 1

        psi0_full, psi1Y_full, psi2_full, YRY_full = self.gatherPsiStat(kern, X, Z, Y, beta, uncertain_inputs)

        #======================================================================
        # Compute Common Components
        #======================================================================

        Kmm = kern.K(Z).copy()
        diag.add(Kmm, self.const_jitter)
        if not np.isfinite(Kmm).all():
            print(Kmm)
        Lm = jitchol(Kmm)
        LmInv = dtrtri(Lm)

        LmInvPsi2LmInvT = LmInv.dot(psi2_full.dot(LmInv.T))
        Lambda = np.eye(Kmm.shape[0])+LmInvPsi2LmInvT
        LL = jitchol(Lambda)
        LLInv = dtrtri(LL)
        logdet_L = 2.*np.sum(np.log(np.diag(LL)))
        LmLLInv = LLInv.dot(LmInv)
        
        b  = psi1Y_full.dot(LmLLInv.T)
        bbt = np.square(b).sum()
        v = b.dot(LmLLInv).T
        LLinvPsi1TYYTPsi1LLinvT = tdot(b.T)
        
        tmp = -LLInv.T.dot(LLinvPsi1TYYTPsi1LLinvT+output_dim*np.eye(input_dim)).dot(LLInv)
        dL_dpsi2R = LmInv.T.dot(tmp+output_dim*np.eye(input_dim)).dot(LmInv)/2.
        
        # Cache intermediate results
        self.midRes['dL_dpsi2R'] = dL_dpsi2R
        self.midRes['v'] = v

        #======================================================================
        # Compute log-likelihood
        #======================================================================
        if het_noise:
            logL_R = -np.log(beta).sum()
        else:
            logL_R = -num_data*np.log(beta)
        logL = -(output_dim*(num_data*log_2_pi+logL_R+psi0_full-np.trace(LmInvPsi2LmInvT))+YRY_full-bbt)/2.-output_dim*logdet_L/2.

        #======================================================================
        # Compute dL_dKmm
        #======================================================================

        dL_dKmm =  dL_dpsi2R - output_dim*LmInv.T.dot(LmInvPsi2LmInvT).dot(LmInv)/2.

        #======================================================================
        # Compute the Posterior distribution of inducing points p(u|Y)
        #======================================================================

        if not self.Y_speedup or het_noise:
            wd_inv = backsub_both_sides(Lm, np.eye(input_dim)- backsub_both_sides(LL, np.identity(input_dim), transpose='left'), transpose='left')
            post = Posterior(woodbury_inv=wd_inv, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=Lm)
        else:
            post = None

        #======================================================================
        # Compute dL_dthetaL for uncertian input and non-heter noise
        #======================================================================

        if not het_noise:
            dL_dthetaL = (YRY_full*beta + beta*output_dim*psi0_full - num_data*output_dim*beta)/2. - beta*(dL_dpsi2R*psi2_full).sum() - beta*(v.T*psi1Y_full).sum()
            self.midRes['dL_dthetaL'] = dL_dthetaL

        return logL, dL_dKmm, post
764701
def correlate(self, signal):
        
        s = asarray(signal)

        if s.ndim == 1:
            if size(s) != self.shape[-1]:
                raise ValueError("Length of signal '%g' does not match record length '%g'"
                                 % (size(s), self.shape[-1]))

            return self.map(lambda x: corrcoef(x, s)[0, 1], index=[1])

        elif s.ndim == 2:
            if s.shape[1] != self.shape[-1]:
                raise ValueError("Length of signal '%g' does not match record length '%g'"
                                 % (s.shape[1], self.shape[-1]))
            newindex = arange(0, s.shape[0])
            return self.map(lambda x: array([corrcoef(x, y)[0, 1] for y in s]), index=newindex)

        else:
            raise Exception('Signal to correlate with must have 1 or 2 dimensions')
387780
def make_from_catalogue(cls, catalogue, spacing, dilate):
        '''
        
        '''
        new = cls()
        cat_bbox = get_catalogue_bounding_polygon(catalogue)

        if dilate > 0:
            cat_bbox = cat_bbox.dilate(dilate)

        # Define Grid spacing
        new.update({'xmin': np.min(cat_bbox.lons),
                    'xmax': np.max(cat_bbox.lons),
                    'xspc': spacing,
                    'ymin': np.min(cat_bbox.lats),
                    'ymax': np.max(cat_bbox.lats),
                    'yspc': spacing,
                    'zmin': 0.,
                    'zmax': np.max(catalogue.data['depth']),
                    'zspc': np.max(catalogue.data['depth'])})

        if new['zmin'] == new['zmax'] == new['zspc'] == 0:
            new['zmax'] = new['zspc'] = 1

        return new
986781
def _axis_in_detector(geometry):
    """"""
    du, dv = geometry.det_axes_init
    axis = geometry.axis
    c = np.array([np.vdot(axis, du), np.vdot(axis, dv)])
    cnorm = np.linalg.norm(c)

    # Check for numerical errors
    assert cnorm != 0

    return c / cnorm
117706
def spatialDomainNoGrid(self):
    
    self.w = np.zeros(self.xw.shape)

    if self.Debug:
      print("w = ")
      print(self.w.shape)
    
    for i in range(len(self.q)):
      # More efficient if we have created some 0-load points
      # (e.g., for where we want output)
      if self.q[i] != 0:
        dist = np.abs(self.xw - self.x[i])
        self.w -= self.q[i] * self.coeff * np.exp(-dist/self.alpha) * \
          ( np.cos(dist/self.alpha) + np.sin(dist/self.alpha) )
408219
def get_per_object_threshold(method, image, threshold, mask=None, labels=None,
                             threshold_range_min = None,
                             threshold_range_max = None,
                             **kwargs):
    
    if labels is None:
        labels = np.ones(image.shape,int)
        if not mask is None:
            labels[np.logical_not(mask)] = 0 
    label_extents = scipy.ndimage.find_objects(labels,np.max(labels))
    local_threshold = np.ones(image.shape,image.dtype)
    for i, extent in enumerate(label_extents, start=1):
        label_mask = labels[extent]==i
        if not mask is None:
            label_mask = np.logical_and(mask[extent], label_mask)
        values = image[extent]
        per_object_threshold = get_global_threshold(
            method, values, mask = label_mask, **kwargs)
        local_threshold[extent][label_mask] = per_object_threshold
    return local_threshold
865639
def sample_number_of_occurrences(self, n=1):
        
        # compute cdf from pmf
        cdf = numpy.cumsum(self.probs_occur)
        n_occ = numpy.digitize(numpy.random.random(n), cdf)
        return n_occ
1316875
def skimage_sinogram_space(geometry, volume_space, sinogram_space):
    """"""
    padded_size = int(np.ceil(volume_space.shape[0] * np.sqrt(2)))
    det_width = volume_space.domain.extent[0] * np.sqrt(2)
    skimage_detector_part = uniform_partition(-det_width / 2.0,
                                              det_width / 2.0,
                                              padded_size)

    skimage_range_part = geometry.motion_partition.insert(
        1, skimage_detector_part)

    skimage_range = uniform_discr_frompartition(skimage_range_part,
                                                interp=sinogram_space.interp,
                                                dtype=sinogram_space.dtype)

    return skimage_range
512377
def uvw(self, context):
        

        # Shape (ntime, na, 3)
        (lt, ut), (la, ua), (l, u) = context.array_extents(context.name)

        # Create empty UVW coordinates
        data = np.empty(context.shape, context.dtype)
        data[:,:,0] = np.arange(la+1, ua+1)    # U = antenna index
        data[:,:,1] = 0                        # V = 0
        data[:,:,2] = 0                        # W = 0

        return data
1135394
def getSparsity(self):
    
    sparsity = numpy.sum(self._unionSDR) / self._numInputs
    return sparsity
118030
def getCatalogPixels(self):
        
        filenames = self.config.getFilenames()

        nside_catalog = self.config.params['coords']['nside_catalog']
        nside_pixel = self.config.params['coords']['nside_pixel']
        # All possible catalog pixels spanned by the ROI
        superpix = ugali.utils.skymap.superpixel(self.pixels,nside_pixel,nside_catalog)
        superpix = np.unique(superpix)
        # Only catalog pixels that exist in catalog files
        pixels = np.intersect1d(superpix, filenames['pix'].compressed())
        return pixels
836704
def spatial_firing_rate(self,label_x='x',label_y='y',bins=None,resolution=1.0,geometry=None,weight_function=None,normalize_time=True,normalize_n=False,start_units_with_0=True):
        
        if bool(self):
            if type(bins) is tuple:
                if type(bins[0]) is int:
                    bins_x = self.spike_times.get_label(label_y).linspace_bins(bins=bins[0])
                else:
                    bins_x = bins[0]
                if type(bins[1]) is int:
                    bins_y = self.spike_times.get_label(label_x).linspace_bins(bins=bins[1])
                else:
                    bins_y = bins[1]
            else:
                if type(bins) is int:
                    bins_x = self.spike_times.get_label(label_y).linspace_bins(bins=bins)
                    bins_y = self.spike_times.get_label(label_x).linspace_bins(bins=bins)
                else:
                    bins_x = bins
                    bins_y = bins
            if weight_function is None:
                H,xed,yed = np.histogram2d(self.spike_times[label_x],self.spike_times[label_y],bins=(bins_x,bins_y))
            else:
                if start_units_with_0:
                    ## TODO : reintroduce a neuron number thing
                    spike_numbers = np.transpose([self.spikes[:,0]-self.n_range[0],self.spikes[:,1]])
                    weights = weight_function(spike_numbers)
                else:
                    weights = weight_function(self.spike_times[label_x],self.spike_times[label_y])
                H,xed,yed = np.histogram2d(self.spike_times[label_x],self.spike_times[label_y],bins=(bins_x,bins_y),weights=weights)
            if normalize_time:
                H = H/(self.spike_times.labels[0].convert('s').len())
            if normalize_n:
                H = H/cells_per_bin
            return H,xed,yed
        else:
            return (np.array([[0]]),np.array([0]),np.array([0]))
198313
def entropy_calc(item, POP):
    
    try:
        result = 0
        for i in item.keys():
            likelihood = item[i] / POP[i]
            if likelihood != 0:
                result += likelihood * math.log(likelihood, 2)
        return -result
    except Exception:
        return "None"
630747
def _get_magnitude_vector_properties(catalogue, config):
    ''''''

    mmin = config.get('input_mmin', np.min(catalogue['magnitude']))
    neq = np.float(np.sum(catalogue['magnitude'] >= mmin - 1.E-7))
    return neq, mmin