Sample the interpolator of a root 2d hist.

    Root's hist2d has a weird internal interpolation routine,
    also using neighbouring bins.
961771
def interpol_hist2d(h2d, oversamp_factor=10):
    
    from rootpy import ROOTError

    xlim = h2d.bins(axis=0)
    ylim = h2d.bins(axis=1)
    xn = h2d.nbins(0)
    yn = h2d.nbins(1)
    x = np.linspace(xlim[0], xlim[1], xn * oversamp_factor)
    y = np.linspace(ylim[0], ylim[1], yn * oversamp_factor)
    mat = np.zeros((xn, yn))
    for xi in range(xn):
        for yi in range(yn):
            try:
                mat[xi, yi] = h2d.interpolate(x[xi], y[yi])
            except ROOTError:
                continue
    return mat, x, y
1014533
def hist2d(x, y, label=None, opacity=None):
    
    x = np.atleast_1d(x)
    y = np.atleast_1d(y)
    data = [go.Histogram2d(x=x, y=y, opacity=opacity, name=label)]
    return Chart(data=data)
919597
def save_root(histogram: HistogramBase, path: str, name: Optional[str] = None):
    
    if name is None:
        name = histogram.name
    if os.path.isfile(path):
        # TODO: Not supported currently
        hfile = uproot.write.TFile.TFileUpdate(path)
    else:
        hfile = uproot.write.TFile.TFileCreate(path)
    write_root(histogram, hfile, name)
    hfile.close()
724656
protected void computeHistogramBorder(T image, RectangleRotate_F32 region) {
		for( int i = 0; i < samplePts.size(); i++ ) {
			Point2D_F32 p = samplePts.get(i);

			squareToImageSample(p.x, p.y, region);

			// make sure its inside the image
			if( !BoofMiscOps.checkInside(image, imageX, imageY)) {
				sampleHistIndex[ i ] = -1;
			} else {
				// use the slower interpolation which can handle the border
				interpolate.get(imageX, imageY, value);

				int indexHistogram = computeHistogramBin(value);

				sampleHistIndex[ i ] = indexHistogram;
				histogram[indexHistogram] += weights[i];
			}
		}
	}
1088467
def _sample_item(self, **kwargs):
        
        t = self.t_
        if 'fixed_stratum' in kwargs:
            stratum_idx = kwargs['fixed_stratum']
        else:
            stratum_idx = None
        if stratum_idx is not None:
            # Sample in given stratum
            loc = self.strata._sample_in_stratum(stratum_idx, replace=False)
            # Record instrumental distribution
            if self.record_inst_hist:
                self.inst_pmf_[stratum_idx,t] = 1
        else:
            # Choose stratum based on instrumental distribution
            self._calc_inst_pmf()

            if self.record_inst_hist:
                inst_pmf = self.inst_pmf_[:,t]
            else:
                inst_pmf = self.inst_pmf_

            loc, stratum_idx = self.strata.sample(pmf = inst_pmf, replace=False)

        return loc, 1, {'stratum': stratum_idx}
1057310
def representative_sample(X, num_samples, save=False):
    
    X = X.values if hasattr(X, 'values') else np.array(X)
    N, M = X.shape
    rownums = np.arange(N)
    np.random.shuffle(rownums)

    idx = AnnoyIndex(M)
    for i, row in enumerate(X):
        idx.add_item(i, row)
    idx.build(int(np.log2(N)) + 1)

    if save:
        if isinstance(save, basestring):
            idxfilename = save
        else:
            idxfile = tempfile.NamedTemporaryFile(delete=False)
            idxfile.close()
            idxfilename = idxfile.name
        idx.save(idxfilename)
        idx = AnnoyIndex(M)
        idx.load(idxfile.name)

    samples = -1 * np.ones(shape=(num_samples,), dtype=int)
    samples[0] = rownums[0]
    # FIXME: some integer determined by N and num_samples and distribution
    j, num_nns = 0, min(1000, int(num_samples / 2. + 1))
    for i in rownums:
        if i in samples:
            continue
        nns = idx.get_nns_by_item(i, num_nns)
        # FIXME: pick vector furthest from past K (K > 1) points or outside of a hypercube
        #        (sized to uniformly fill the space) around the last sample
        samples[j + 1] = np.setdiff1d(nns, samples)[-1]
        if len(num_nns) < num_samples / 3.:
            num_nns = min(N, 1.3 * num_nns)
        j += 1
    return samples
51756
def histogram1d(data, bins=None, *args, **kwargs):
    
    import dask
    if not hasattr(data, "dask"):
        data = dask.array.from_array(data, chunks=int(data.shape[0] / options["chunk_split"]))

    if not kwargs.get("adaptive", True):
        raise RuntimeError("Only adaptive histograms supported for dask (currently).")
    kwargs["adaptive"] = True

    def block_hist(array):
        return original_h1(array, bins, *args, **kwargs)

    return _run_dask(
        name="dask_adaptive1d",
        data=data,
        compute=kwargs.pop("compute", True),
        method=kwargs.pop("dask_method", "threaded"),
        func=block_hist)
455723
def inverseHistogram(hist, bin_range):
    
    data = hist.astype(float) / np.min(hist[np.nonzero(hist)])
    new_data = np.empty(shape=np.sum(data, dtype=int))
    i = 0
    xvals = np.linspace(bin_range[0], bin_range[1], len(data))
    for d, x in zip(data, xvals):
        new_data[i:i + d] = x
        i += int(d)
    return new_data
854983
def interpolation_points(self, N):
        
        if N == 1:
            return np.array([0.])
        return np.cos(np.arange(N)*np.pi/(N-1))
470664
def sample_given_context(self, c, c_dims):
        '''
        
        '''
        index = self.discrete_progress.sample_given_context(c, c_dims, self.space)
        return self.space.rand_value(index).flatten()[list(set(range(len(self.space.cardinalities))) - set(c_dims))]
344245
def histogram2d(data1, data2, bins=None, *args, **kwargs):
    """"""
    # TODO: currently very unoptimized! for non-dasks
    import dask
    if "axis_names" not in kwargs:
        if hasattr(data1, "name") and hasattr(data2, "name"):
            kwargs["axis_names"] = [data1.name, data2.name]
    if not hasattr(data1, "dask"):
        data1 = dask.array.from_array(data1, chunks=data1.size() / 100)
    if not hasattr(data2, "dask"):
        data2 = dask.array.from_array(data2, chunks=data2.size() / 100)

    data = dask.array.stack([data1, data2], axis=1)
    kwargs["dim"] = 2
    return histogramdd(data, bins, *args, **kwargs)
1115297
def interpolate(self):
        
        nans, x = utils.nan_idx(self)
        self[nans] = np.interp(x(nans), x(~nans), self[~nans])
        return self
1127424
def interpolateCircular2dStructuredIDW(grid, mask, kernel=15, power=2, 
                               fr=1, fphi=1, cx=0, cy=0):
    '''
     
    '''
    gx = grid.shape[0]
    gy = grid.shape[0]

    #FOR EVERY PIXEL
    for i in range(gx):
        for j in range(gy):
            
            if mask[i,j]:
                
                xmn = i-kernel
                if xmn < 0:
                    xmn = 0
                xmx = i+kernel
                if xmx > gx:
                    xmx = gx
                    
                ymn = j-kernel
                if ymn < 0:
                    ymn = 0
                ymx = j+kernel
                if ymx > gx:
                    ymx = gy
                    
                sumWi = 0.0
                value = 0.0 

                #radius and radian to polar center:
                R = ((i-cx)**2+(j-cy)**2)**0.5
                PHI = atan2(j-cy, i-cx)
                
                #FOR EVERY NEIGHBOUR IN KERNEL              
                for xi in range(xmn,xmx):
                    for yi in range(ymn,ymx):
                        if  (xi != i or yi != j) and not mask[xi,yi]:
                            nR = ((xi-cx)**2+(yi-cy)**2)**0.5
                            dr =  R - nR
                            #average radius between both p:
                            midR = 0.5*(R+nR)
                            #radian of neighbour p:
                            nphi = atan2(yi-cy, xi-cx)
                            #relative angle between both points:
                            dphi = min((2*np.pi) - abs(PHI - nphi), 
                                       abs(PHI - nphi))    
                            dphi*=midR       
                            
                            dist = ((fr*dr)**2+(fphi*dphi)**2)**2

                            wi = 1 / dist**(0.5*power)
                            sumWi += wi
                            value += wi * grid[xi,yi]  
                if sumWi:
                    grid[i,j] = value / sumWi               

    return grid
1303067
def extrap1d(interpolator):
    
    
    xs = interpolator.x
    ys = interpolator.y

    def pointwise(x):
        if x < xs[0]:
            return ys[0] # +(x-xs[0])*(ys[1]-ys[0])/(xs[1]-xs[0])
        elif x > xs[-1]:
            return ys[-1] # +(x-xs[-1])*(ys[-1]-ys[-2])/(xs[-1]-xs[-2])
        else:
            return interpolator(x)

    def ufunclike(xs):
        return N.array(map(pointwise, N.array(xs)))

    return ufunclike
997394
def hierarchical_clusterization_2d(no_clusters = 0, distance_method = 0,
                                       vec = @values, debug = false)
      clusters = []

      if vec.length == 1
        hash = { vec[0] => 1 }
        cluster = PairCluster.new(hash)
        clusters.push(cluster)
        clusters
      end

      # Thresholds
      # threshold_distance = (0.25 * (vec.max-vec.min))
      threshold_density = (0.5 * vec.length).to_i

      # make a histogram from the input vector
      histogram = Hash[vec.group_by { |a| a }.map { |k, vs| [k, vs.length] }]

      # clusters = array of clusters
      # initially each length belongs to a different cluster
      histogram.each do |e|
        warn "pair (#{e[0].x} #{e[0].y}) appears #{e[1]} times" if debug
        hash = { e[0] => e[1] }
        cluster = PairCluster.new(hash)
        clusters.push(cluster)
      end

      clusters.each(&:print) if debug

      return clusters if clusters.length == 1

      # each iteration merge the closest two adiacent cluster
      # the loop stops according to the stop conditions
      iteration = 0
      loop do
        # stop condition 1
        break if no_clusters != 0 && clusters.length == no_clusters

        iteration += iteration
        warn "\nIteration #{iteration}" if debug

        min_distance = 100_000_000
        cluster1     = 0
        cluster2     = 0
        density      = 0

        [*(0..(clusters.length - 2))].each do |i|
          [*((i + 1)..(clusters.length - 1))].each do |j|
            dist = clusters[i].distance(clusters[j], distance_method)
            warn "distance between clusters #{i} and #{j} is #{dist}" if debug
            current_density = clusters[i].density + clusters[j].density
            if dist < min_distance
              min_distance = dist
              cluster1     = i
              cluster2     = j
              density      = current_density
            elsif dist == min_distance && density < current_density
              cluster1 = i
              cluster2 = j
              density  = current_density
            end
          end
        end

        # merge clusters 'cluster1' and 'cluster2'
        warn "clusters to merge #{cluster1} and #{cluster2}" if debug

        clusters[cluster1].add(clusters[cluster2])
        clusters.delete_at(cluster2)

        if debug
          clusters.each_with_index do |elem, i|
            warn "cluster #{i}"
            elem.print
          end
        end

        # stop condition 3
        # the density of the biggest clusters exceeds the threshold
        if no_clusters == 0 && clusters[cluster].density > threshold_density
          break
        end
      end

      @clusters = clusters
    end
155809
def _sampleFromEdges(self):
    
    dimensionsToChooseFrom = range(self.dimension)
    random.shuffle(dimensionsToChooseFrom)
    dim1, dim2 = dimensionsToChooseFrom[0], dimensionsToChooseFrom[1]
    coordinates = [random.uniform(-1, 1) * dim / 2. for dim in self.dimensions]
    coordinates[dim1] = self.dimensions[dim1] / 2. * random.choice([-1, 1])
    coordinates[dim2] = self.dimensions[dim2] / 2. * random.choice([-1, 1])
    return coordinates
544151
def sample(self, histogram_logits):
        
        histogram_probs = histogram_logits.exp()  # Batch size * actions * atoms
        atoms = self.support_atoms.view(1, 1, self.atoms)  # Need to introduce two new dimensions
        return (histogram_probs * atoms).sum(dim=-1).argmax(dim=1)