Python numpy.histogram2d() Examples

def pred_table(self, threshold=.5):
        """
        Prediction table

        Parameters
        ----------
        threshold : scalar
            Number between 0 and 1. Threshold above which a prediction is
            considered 1 and below which a prediction is considered 0.

        Notes
        ------
        pred_table[i,j] refers to the number of times "i" was observed and
        the model predicted "j". Correct predictions are along the diagonal.
        """
        model = self.model
        actual = model.endog
        pred = np.array(self.predict() > threshold, dtype=float)
        bins = np.array([0, 0.5, 1])
        return np.histogram2d(actual, pred, bins=bins)[0] 

def _mutual_information_varoquaux(x, y, bins=256, sigma=1, normalized=True):
    """Based on Gael Varoquaux's implementation: https://gist.github.com/GaelVaroquaux/ead9898bd3c973c40429."""
    jh = np.histogram2d(x, y, bins=bins)[0]

    # smooth the jh with a gaussian filter of given sigma
    scipy.ndimage.gaussian_filter(jh, sigma=sigma, mode="constant", output=jh)

    # compute marginal histograms
    jh = jh + np.finfo(float).eps
    sh = np.sum(jh)
    jh = jh / sh
    s1 = np.sum(jh, axis=0).reshape((-1, jh.shape[0]))
    s2 = np.sum(jh, axis=1).reshape((jh.shape[1], -1))

    if normalized:
        mi = ((np.sum(s1 * np.log(s1)) + np.sum(s2 * np.log(s2))) / np.sum(jh * np.log(jh))) - 1
    else:
        mi = np.sum(jh * np.log(jh)) - np.sum(s1 * np.log(s1)) - np.sum(s2 * np.log(s2))

    return mi 

def hist(self, nplan, xedges, yedges):
        """Returns completeness histogram for Monte Carlo simulation
        
        This function uses the inherited Planet Population module.
        
        Args:
            nplan (float):
                number of planets used
            xedges (float ndarray):
                x edge of 2d histogram (separation)
            yedges (float ndarray):
                y edge of 2d histogram (dMag)
        
        Returns:
            h (ndarray):
                2D numpy ndarray containing completeness histogram
        
        """
        
        s, dMag = self.genplans(nplan)
        # get histogram
        h, yedges, xedges = np.histogram2d(dMag, s.to('AU').value, bins=1000,
                range=[[yedges.min(), yedges.max()], [xedges.min(), xedges.max()]])
        
        return h, xedges, yedges 

def fit(self, magnitude, time, dt_bins, dm_bins):
        def delta_calc(idx):
            t0 = time[idx]
            m0 = magnitude[idx]
            deltat = time[idx + 1 :] - t0
            deltam = magnitude[idx + 1 :] - m0

            deltat[np.where(deltat < 0)] *= -1
            deltam[np.where(deltat < 0)] *= -1

            return np.column_stack((deltat, deltam))

        lc_len = len(time)
        n_vals = int(0.5 * lc_len * (lc_len - 1))

        deltas = np.vstack(tuple(delta_calc(idx) for idx in range(lc_len - 1)))

        deltat = deltas[:, 0]
        deltam = deltas[:, 1]

        bins = [dt_bins, dm_bins]
        counts = np.histogram2d(deltat, deltam, bins=bins, normed=False)[0]
        result = np.fix(255.0 * counts / n_vals + 0.999).astype(int)

        return {"DMDT": result} 

def gridded_mean(lat, lon, data, grid):
    """Grid data along latitudes and longitudes

    Args:
        lat: Grid points along latitudes as 1xN dimensional array.
        lon: Grid points along longitudes as 1xM dimensional array.
        data: The data as NxM numpy array.
        grid: A tuple with two numpy arrays, consisting of latitude and
            longitude grid points.

    Returns:
        Two matrices in grid form: the mean and the number of points of `data`.
    """
    grid_sum, _, _ = np.histogram2d(lat, lon, grid, weights=data)
    grid_number, _, _ = np.histogram2d(lat, lon, grid)

    return grid_sum / grid_number, grid_number 

def morista_index(points):
    # Morisita Index of Dispersion

    N = points.shape[1]

    ims = []
    for i in range(1, N):
        bins, _, _ = np.histogram2d(points[0], points[1], i)

        # I_M  = Q * (\sum_{k=1}^{Q}{n_k * (n_k - 1)})/(N * (N _ 1))
        Q = len(bins)  # num_quadrants

        # Eqn 1.
        I_M = Q * np.sum(np.ravel(bins) * (np.ravel(bins) - 1)) / (N * (N - 1))
        ims.append([i, I_M])


    return np.array(ims).T[1].max() 

def pred_table(self, threshold=.5):
        """
        Prediction table

        Parameters
        ----------
        threshold : scalar
            Number between 0 and 1. Threshold above which a prediction is
            considered 1 and below which a prediction is considered 0.

        Notes
        ------
        pred_table[i,j] refers to the number of times "i" was observed and
        the model predicted "j". Correct predictions are along the diagonal.
        """
        model = self.model
        actual = model.endog
        pred = np.array(self.predict() > threshold, dtype=float)
        bins = np.array([0, 0.5, 1])
        return np.histogram2d(actual, pred, bins=bins)[0] 

def pix_histogram(self, mask, lbl):
    '''
      get individual mask and label and create 2d hist
    '''
    # flatten mask and cast
    flat_mask = mask.flatten().astype(np.uint32)
    # flatten label and cast
    flat_label = lbl.flatten().astype(np.uint32)
    # get the histogram
    histrange = np.array([[-0.5, self.num_classes - 0.5],
                          [-0.5, self.num_classes - 0.5]], dtype='float64')
    h_now, _, _ = np.histogram2d(np.array(flat_mask),
                                 np.array(flat_label),
                                 bins=self.num_classes,
                                 range=histrange)
    return h_now 

def test_displaced_two_mode_state_hafnian(self, sample_func):
        """Test the sampling routines by comparing the photon number frequencies and the exact
        probability distribution of a two mode coherent state
        """
        n_samples = 1000
        n_cut = 6
        sigma = np.identity(4)
        mean = 5 * np.array([0.1, 0.25, 0.1, 0.25])
        samples = sample_func(sigma, samples=n_samples, mean=mean, cutoff=n_cut)
        # samples = hafnian_sample_classical_state(sigma, mean = mean, samples = n_samples)
        probs = np.real_if_close(
            np.array(
                [
                    [density_matrix_element(mean, sigma, [i, j], [i, j]) for i in range(n_cut)]
                    for j in range(n_cut)
                ]
            )
        )
        freq, _, _ = np.histogram2d(samples[:, 1], samples[:, 0], bins=np.arange(0, n_cut + 1))
        rel_freq = freq / n_samples

        assert np.allclose(
            rel_freq, probs, rtol=rel_tol / np.sqrt(n_samples), atol=rel_tol / np.sqrt(n_samples)
        ) 

def hist2d(self, x, y, weights, bins):   #bins[nx,ny]
        ### generates a 2d hist from input 1d axis for x,y. repeats them to match shape of weights X*Y (data points)
        ### x will be 0-100%
        freqs = np.repeat(np.array([y], dtype=np.float64), len(x), axis=0)
        throts = np.repeat(np.array([x], dtype=np.float64), len(y), axis=0).transpose()
        throt_hist_avr, throt_scale_avr = np.histogram(x, 101, [0, 100])

        hist2d = np.histogram2d(throts.flatten(), freqs.flatten(),
                                range=[[0, 100], [y[0], y[-1]]],
                                bins=bins, weights=weights.flatten(), normed=False)[0].transpose()

        hist2d = np.array(abs(hist2d), dtype=np.float64)
        hist2d_norm = np.copy(hist2d)
        hist2d_norm /=  (throt_hist_avr + 1e-9)

        return {'hist2d_norm':hist2d_norm, 'hist2d':hist2d, 'throt_hist':throt_hist_avr,'throt_scale':throt_scale_avr} 

def test_gate_fraction_2_error_large_gate_fraction(self):
        bins = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5]
        self.assertRaises(
            ValueError,
            FlowCal.gate.density2d,
            self.pyramid,
            bins=bins,
            gate_fraction=1.1,
            sigma=0.0,
            )

    # Test implicit gating (when values exist outside specified bins)
    #
    # The expected behavior is that density2d() should mimic
    # np.histogram2d()'s behavior: values outside the specified bins are
    # ignored (in the context of a gate function, this means they are
    # implicitly gated out). 

def test_sub_bin_1_mask(self):
        bins = [0.5, 2.5, 4.5]
        np.testing.assert_array_equal(
            FlowCal.gate.density2d(
                self.slope, bins=bins, gate_fraction=13.0/30, sigma=0.0,
                full_output=True).mask,
            np.array([0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,0,0,0,1,1,0,1,1,
                1,1,1,1,1], dtype=bool)
            )

    # Test bins edge case (when bin edges = values)
    #
    # Again, the expected behavior is that density2d() should mimic
    # np.histogram2d()'s behavior which will group the right-most two values
    # together in the same bin (since the last bins interval is fully closed,
    # as opposed to all other bins intervals which are half-open). 

def spectral_distribution(points, cumulative=True):
    """ 
    Returns the distribution of complex values (in r,theta-space).
    """

    points = np.array([(np.abs(z), np.angle(z)) for z in points])
    r, theta = np.split(points, 2, axis=1)

    r = np.array([logr(x, r.max()) for x in r])

    Z, R, THETA = np.histogram2d(
        x=r[:, 0],
        y=theta[:, 0],
        bins=(np.linspace(0, 1, 101), np.linspace(0, np.pi, 101)),
    )

    if cumulative:
        Z = Z.cumsum(axis=0).cumsum(axis=1)
        Z = Z / Z.max()

    return Z.flatten() 

def mutual_information(x, y, bins=8):
    """Mutual information score with set number of bins

    Helper function for `sklearn.metrics.mutual_info_score` that builds a
    contingency table over a set number of bins.
    Credit: `Warran Weckesser <https://stackoverflow.com/a/20505476/3996580>`_.


    Parameters
    ----------
    x : array-like, shape=[n_samples]
        Input data (feature 1)
    y : array-like, shape=[n_samples]
        Input data (feature 2)
    bins : int or array-like, (default: 8)
        Passed to np.histogram2d to calculate a contingency table.

    Returns
    -------
    mi : float
        Mutual information between x and y.

    Examples
    --------
    >>> import scprep
    >>> data = scprep.io.load_csv("my_data.csv")
    >>> mi = scprep.stats.mutual_information(data['GENE1'], data['GENE2'])
    """
    x, y = _vector_coerce_two_dense(x, y)
    c_xy = np.histogram2d(x, y, bins)[0]
    mi = metrics.mutual_info_score(None, None, contingency=c_xy)
    return mi 

def hist(self, nplan, xedges, yedges):
        """Returns completeness histogram for Monte Carlo simulation
        
        This function uses the inherited Planet Population module.
        
        Args:
            nplan (float):
                number of planets used
            xedges (float ndarray):
                x edge of 2d histogram (separation)
            yedges (float ndarray):
                y edge of 2d histogram (dMag)
        
        Returns:
            float ndarray:
                2D numpy ndarray containing completeness frequencies
        
        """
        
        s, dMag = self.genplans(nplan)
        # get histogram
        h, yedges, xedges = np.histogram2d(dMag, s.to('AU').value, bins=1000,
                range=[[yedges.min(), yedges.max()], [xedges.min(), xedges.max()]])
        
        return h, xedges, yedges 

def mutual_information(x, y, nbins=32, normalized=False):
    """
    Compute mutual information
    :param x: 1D numpy.array : flatten data from an image
    :param y: 1D numpy.array : flatten data from an image
    :param nbins: number of bins to compute the contingency matrix (only used if normalized=False)
    :return: float non negative value : mutual information
    """
    import sklearn.metrics
    if normalized:
        mi = sklearn.metrics.normalized_mutual_info_score(x, y)
    else:
        c_xy = np.histogram2d(x, y, nbins)[0]
        mi = sklearn.metrics.mutual_info_score(None, None, contingency=c_xy)
    # mi = adjusted_mutual_info_score(None, None, contingency=c_xy)
    return mi 

def find_joint_probability_distribution(TS, rang, nbins):
    """
    Assign each node j to a vector of length nbins where each element is the product of its own
    individual_probability_distribution and its neighbors'. P(x) * P(y) <-- as opposed to P(x,y)

    Parameters
    ----------
    TS (np.ndarray): Array consisting of $L$ observations from $N$ sensors.
    rang (list): list of the minimum and maximum value in the time series
    nbins (int): number of bins for the pre-processing step (to yield a discrete probability distribution)

    Returns
    -------
    JointP (dict): a dictionary where the keys are nodes in the graph and the
                         are nbins x nbins arrays corresponding to joint probability vectors
    """

    N, L = TS.shape  # N nodes and L length
    JointP = dict()  # create a dict to put the joint prob distributions into

    for l in range(N):  # for each node
        for j in range(l):  # for each possible edge
            P, _, _ = np.histogram2d(
                TS[j], TS[l], bins=nbins, range=np.array([rang, rang])
            )
            JointP[(j, l)] = P / L

    return JointP 

def test_missing_range(self):
        x = np.array([1, 2, 3, 4, 5])
        y = np.array([5, 7, 1, 5, 9])
        with self.assertWarns(PrivacyLeakWarning):
            res = histogram2d(x, y, epsilon=1, range=[(0, 10), None])
        self.assertIsNotNone(res) 

def weighted_mode_avr(self, values, weights, vertrange, vertbins):
        ### finds the most common trace and std
        threshold = 0.5  # threshold for std calculation
        filt_width = 7  # width of gaussian smoothing for hist data

        resp_y = np.linspace(vertrange[0], vertrange[-1], vertbins, dtype=np.float64)
        times = np.repeat(np.array([self.time_resp],dtype=np.float64), len(values), axis=0)
        weights = np.repeat(weights, len(values[0]))

        hist2d = np.histogram2d(times.flatten(), values.flatten(),
                                range=[[self.time_resp[0], self.time_resp[-1]], vertrange],
                                bins=[len(times[0]), vertbins], weights=weights.flatten())[0].transpose()
        ### shift outer edges by +-1e-5 (10us) bacause of dtype32. Otherwise different precisions lead to artefacting.
        ### solution to this --> somethings strage here. In outer most edges some bins are doubled, some are empty.
        ### Hence sometimes produces "divide by 0 error" in "/=" operation.

        if hist2d.sum():
            hist2d_sm = gaussian_filter1d(hist2d, filt_width, axis=0, mode='constant')
            hist2d_sm /= np.max(hist2d_sm, 0)


            pixelpos = np.repeat(resp_y.reshape(len(resp_y), 1), len(times[0]), axis=1)
            avr = np.average(pixelpos, 0, weights=hist2d_sm * hist2d_sm)
        else:
            hist2d_sm = hist2d
            avr = np.zeros_like(self.time_resp)
        # only used for monochrome error width
        hist2d[hist2d <= threshold] = 0.
        hist2d[hist2d > threshold] = 0.5 / (vertbins / (vertrange[-1] - vertrange[0]))

        std = np.sum(hist2d, 0)

        return avr, std, [self.time_resp, resp_y, hist2d_sm]

    ### calculates weighted avverage and resulting errors 

def test_f32_rounding(self):
        # gh-4799, check that the rounding of the edges works with float32
        x = np.array([276.318359, -69.593948, 21.329449], dtype=np.float32)
        y = np.array([5005.689453, 4481.327637, 6010.369629], dtype=np.float32)
        counts_hist, xedges, yedges = np.histogram2d(x, y, bins=100)
        assert_equal(counts_hist.sum(), 3.) 

def _get_smoothed_histogram2d(self, chain, param1, param2):  # pragma: no cover
        # No test coverage here because
        smooth = chain.config["smooth"]
        bins = chain.config["bins"]
        kde = chain.config["kde"]

        x = chain.get_data(param1)
        y = chain.get_data(param2)
        w = chain.weights

        if chain.grid:
            binsx = get_grid_bins(x)
            binsy = get_grid_bins(y)
            hist, x_bins, y_bins = np.histogram2d(x, y, bins=[binsx, binsy], weights=w)
        else:
            binsx, smooth = get_smoothed_bins(smooth, bins, x, w, marginalised=False)
            binsy, _ = get_smoothed_bins(smooth, bins, y, w, marginalised=False)
            hist, x_bins, y_bins = np.histogram2d(x, y, bins=[binsx, binsy], weights=w)

        if chain.power is not None:
            hist = hist ** chain.power

        x_centers = 0.5 * (x_bins[:-1] + x_bins[1:])
        y_centers = 0.5 * (y_bins[:-1] + y_bins[1:])

        if kde:
            nn = x_centers.size * 2  # Double samples for KDE because smooth
            x_centers = np.linspace(x_bins.min(), x_bins.max(), nn)
            y_centers = np.linspace(y_bins.min(), y_bins.max(), nn)
            xx, yy = meshgrid(x_centers, y_centers, indexing="ij")
            coords = np.vstack((xx.flatten(), yy.flatten())).T
            data = np.vstack((x, y)).T
            hist = MegKDE(data, w, kde).evaluate(coords).reshape((nn, nn))
            if chain.power is not None:
                hist = hist ** chain.power
        elif smooth:
            hist = gaussian_filter(hist, smooth, mode=self.parent._gauss_mode)

        return hist, x_centers, y_centers 

def make_density_map(paths, map_config):
    xs = np.linspace(map_config.xs[0], map_config.xs[1], num=map_config.xres+1)
    ys = np.linspace(map_config.ys[0], map_config.ys[1], num=map_config.yres+1)
    y = paths[:,0]
    x = paths[:,1]
    H, xedges, yedges = np.histogram2d(y, x, bins=(xs, ys))
    H = H.astype(np.float)
    H = H/np.max(H)
    return H.T 

def convertMovie(
    runner,
    photondist,
    structures,
    imagesize,
    frames,
    psf,
    photonrate,
    background,
    noise,
    mode3Dstate,
    cx,
    cy,
):
    edges = range(0, imagesize + 1)

    photonposframe = distphotonsxy(
        runner, photondist, structures, psf, mode3Dstate, cx, cy
    )

    if len(photonposframe) == 0:
        simframe = _np.zeros((imagesize, imagesize))
    else:
        x = photonposframe[:, 0]
        y = photonposframe[:, 1]
        simframe, xedges, yedges = _np.histogram2d(y, x, bins=(edges, edges))
        simframe = _np.flipud(simframe)  # to be consistent with render

    return simframe 

def calc_score(labels, y_pred):
    if y_pred.sum() == 0 and labels.sum() == 0:
        return 1
    if labels.sum() == 0 and y_pred.sum() > 0 or y_pred.sum() == 0 and labels.sum() > 0:
        return 0

    true_objects = len(np.unique(labels))
    pred_objects = len(np.unique(y_pred))

    intersection = np.histogram2d(labels.flatten(), y_pred.flatten(), bins=(true_objects, pred_objects))[0]

    # Compute areas (needed for finding the union between all objects)
    area_true = np.histogram(labels, bins=true_objects)[0]
    area_pred = np.histogram(y_pred, bins=pred_objects)[0]
    area_true = np.expand_dims(area_true, -1)
    area_pred = np.expand_dims(area_pred, 0)

    # Compute union
    union = area_true + area_pred - intersection

    # Exclude background from the analysis
    intersection = intersection[1:, 1:]
    union = union[1:, 1:]
    union[union == 0] = 1e-9

    # Compute the intersection over union
    iou = intersection / union
    return precision_at(0.5, iou)


# Precision helper function 

def calc_iou(gt_masks, predicted_masks, height=768, width=768):
    true_objects = gt_masks.shape[2]
    pred_objects = predicted_masks.shape[2]
    labels = np.zeros((height, width), np.uint16)
    
    for index in range(0, true_objects):
        labels[gt_masks[:, :, index] > 0] = index + 1
    y_true = labels.flatten()
    labels_pred = np.zeros((height, width), np.uint16)
    
    for index in range(0, pred_objects):
        if sum(predicted_masks[:, :, index].shape) == height + width:
            labels_pred[predicted_masks[:, :, index] > 0] = index + 1
    y_pred = labels_pred.flatten()
    
    intersection = np.histogram2d(y_true, y_pred, bins=(true_objects + 1, pred_objects + 1))[0]
    
    area_true = np.histogram(labels, bins=true_objects + 1)[0]
    area_pred = np.histogram(y_pred, bins=pred_objects + 1)[0]
    area_true = np.expand_dims(area_true, -1)
    area_pred = np.expand_dims(area_pred, 0)
    
    # Compute union
    union = area_true + area_pred - intersection
    
    # Exclude background from the analysis
    intersection = intersection[1:, 1:]
    union = union[1:, 1:]
    union[union == 0] = 1e-9
    # print(union)
    # print(intersection)
    iou = intersection / union
    return iou 

def test_f32_rounding(self):
        # gh-4799, check that the rounding of the edges works with float32
        x = np.array([276.318359  , -69.593948  , 21.329449], dtype=np.float32)
        y = np.array([5005.689453, 4481.327637, 6010.369629], dtype=np.float32)
        counts_hist, xedges, yedges = np.histogram2d(x, y, bins=100)
        assert_equal(counts_hist.sum(), 3.) 

def test_no_params(self):
        x = np.array([1, 2, 3, 4, 5])
        y = np.array([5, 7, 1, 5, 9])
        with self.assertWarns(PrivacyLeakWarning):
            res = histogram2d(x, y)
        self.assertIsNotNone(res) 

def calc_MI(X, Y, bins):
    c_XY = np.histogram2d(X, Y, bins)[0]
    c_X = np.histogram(X, bins)[0]
    c_Y = np.histogram(Y, bins)[0]
    H_X = shan_entropy(c_X)
    H_Y = shan_entropy(c_Y)
    H_XY = shan_entropy(c_XY)
    MI = H_X + H_Y - H_XY
    return MI 

def __call__(self):
        source1 = self.request.GET['source1']
        source2 = self.request.GET['source2']
        data1, max1, name1 = self.getData(source1)
        data2, max2, name2 = self.getData(source2)
        student_ids =  set(data1.keys()).intersection(list(data2.keys()))
        data = np.array([(float(data1[s_id]), float(data2[s_id])) for s_id in student_ids])
        if len(data):
            x = data[:,0]
            y = data[:,1]
        else: x,y = [], []
        bins1 = self.getBins(max1)
        bins2 = self.getBins(max2)

        try:
            hist,xedges,yedges = np.histogram2d(x,y,bins=[bins1, bins2])
        except ValueError:
            # x,y = [] raises ValueErrors in old numpy
            hist = np.array([[0]])
            xedges = [0,1]
            yedges = xedges

        color_stepsize = int(math.ceil(hist.max()/10.0)) or 1
        color_ticks = list(range(0, int(hist.max()) or 1, color_stepsize))
        color_bins = [-color_stepsize/2.0]
        color_bins.extend([t+color_stepsize/2.0 for t in color_ticks])

        norm = matplotlib.colors.BoundaryNorm(color_bins, len(color_ticks))

        extent = [xedges[0], xedges[-1], yedges[0], yedges[-1] ]
        ax = self.fig.add_subplot(111)
        im = ax.imshow(hist.T,extent=extent,interpolation='nearest',
                origin='lower',
                cmap = pyplot.cm.gray_r,
                aspect='auto')
        ax.set_xlabel(name1)
        ax.set_ylabel(name2)
        self.fig.colorbar(im, norm=norm, boundaries=color_bins, ticks=color_ticks)
        corrcoef = np.corrcoef(x, y)[0,1] if len(x)>0 else 0
        ax.set_title("Korrelation = %.2f" % corrcoef)
        return self.createResponse() 

def test_no_range(self):
        x = np.array([1, 2, 3, 4, 5])
        y = np.array([5, 7, 1, 5, 9])
        with self.assertWarns(PrivacyLeakWarning):
            res = histogram2d(x, y, epsilon=1)
        self.assertIsNotNone(res)