Python scipy.interpolate.interp1d() Examples

def _advect_surface(self, dt):

        if self.top:
            # Extract top surface
            x = self.model.mesh.data[self.top.data][:, 0]
            y = self.model.mesh.data[self.top.data][:, 1]

            # Extract velocities from top
            vx = self.model.velocityField.data[self.top.data][:, 0]
            vy = self.model.velocityField.data[self.top.data][:, 1]

            # Advect top surface
            x2 = x + vx * nd(dt)
            y2 = y + vy * nd(dt)

            # Spline top surface
            f = interp1d(x2, y2, kind='cubic', fill_value='extrapolate')

            self.TField.data[self.top.data, 0] = f(x)
        comm.Barrier()
        self.TField.syncronise() 

def calclogf(self):
        """
        # wavelength dependence, from Table 19 in Leinert et al 1998
        # interpolated w/ a quadratic in log-log space
        Returns:
            interpolant (object):
                a 1D quadratic interpolant of intensity vs wavelength

        """
        self.zodi_lam = np.array([0.2, 0.3, 0.4, 0.5, 0.7, 0.9, 1.0, 1.2, 2.2, 3.5,
                4.8, 12, 25, 60, 100, 140]) # um
        self.zodi_Blam = np.array([2.5e-8, 5.3e-7, 2.2e-6, 2.6e-6, 2.0e-6, 1.3e-6,
                1.2e-6, 8.1e-7, 1.7e-7, 5.2e-8, 1.2e-7, 7.5e-7, 3.2e-7, 1.8e-8,
                3.2e-9, 6.9e-10]) # W/m2/sr/um
        x = np.log10(self.zodi_lam)
        y = np.log10(self.zodi_Blam)
        return interp1d(x, y, kind='quadratic') 

def test_cubic(self):
        """ Check the actual implementation of spline interpolation.
        """

        interp10 = interp1d(self.x10, self.y10, kind='cubic')
        assert_array_almost_equal(
            interp10(self.x10),
            self.y10,
        )
        assert_array_almost_equal(
            interp10(1.2),
            np.array([1.2]),
        )
        assert_array_almost_equal(
            interp10([2.4, 5.6, 6.0]),
            np.array([2.4, 5.6, 6.0]),
        ) 

def test_linear(self):
        """ Check the actual implementation of linear interpolation.
        """

        interp10 = interp1d(self.x10, self.y10)
        assert_array_almost_equal(
            interp10(self.x10),
            self.y10,
        )
        assert_array_almost_equal(
            interp10(1.2),
            np.array([1.2]),
        )
        assert_array_almost_equal(
            interp10([2.4, 5.6, 6.0]),
            np.array([2.4, 5.6, 6.0]),
        ) 

def _resample_regressor(hr_regressor, hr_frame_times, frame_times):
    """ this function sub-samples the regressors at frame times

    Parameters
    ----------
    hr_regressor : array of shape(n_samples),
        the regressor time course sampled at high temporal resolution

    hr_frame_times : array of shape(n_samples),
        the corresponding time stamps

    frame_times: array of shape(n_scans),
         the desired time stamps

    Returns
    -------
    regressor: array of shape(n_scans)
         the resampled regressor
    """
    from scipy.interpolate import interp1d
    f = interp1d(hr_frame_times, hr_regressor)
    return f(frame_times).T 

def monotone_fn_inverter(fn, x, vectorized=True, **keywords):
    """
    Given a monotone function fn (no checking is done to verify monotonicity)
    and a set of x values, return an linearly interpolated approximation
    to its inverse from its values on x.
    """
    x = np.asarray(x)
    if vectorized:
        y = fn(x, **keywords)
    else:
        y = []
        for _x in x:
            y.append(fn(_x, **keywords))
        y = np.array(y)

    a = np.argsort(y)

    return interp1d(y[a], x[a]) 

def apply_sync(sync_file, times, forward=True):
    """
    :param sync_file: probe sync file (usually of the form _iblrig_ephysData.raw.imec1.sync.npy)
    :param times: times in seconds to interpolate
    :param forward: if True goes from probe time to session time, from session time to probe time
    otherwise
    :return: interpolated times
    """
    sync_points = np.load(sync_file)
    if forward:
        fcn = interp1d(sync_points[:, 0],
                       sync_points[:, 1], fill_value='extrapolate')
    else:
        fcn = interp1d(sync_points[:, 1],
                       sync_points[:, 0], fill_value='extrapolate')
    return fcn(times) 

def percentileofscore(data, score):
    """Return the percentile-position of score relative to data.

    score: Array of scores at which the percentile is computed.

    Return percentiles (0-100).

    Example
            r = randn(50)
        x = linspace(-2,2,100)
        percentileofscore(r,x)

    Raise an error if the score is outside the range of data.
    """
    cdf = empiricalcdf(data)
    interpolator = interpolate.interp1d(np.sort(data), np.sort(cdf))
    return interpolator(score)*100. 

def load_spectrum(filename, grid):
    """
    Load a single spectrum
    """
    file_in = pyfits.open(filename)
    wl = np.array(file_in[0].data[2])
    flux = np.array(file_in[0].data[0])
    ivar = np.array((file_in[0].data[1]))
    # correct for radial velocity of star
    redshift = file_in[0].header['Z']
    wl_shifted = wl - redshift * wl
    # resample
    flux_rs = (interpolate.interp1d(wl_shifted, flux))(grid)
    ivar_rs = (interpolate.interp1d(wl_shifted, ivar))(grid)
    ivar_rs[ivar_rs < 0] = 0. # in interpolating you can end up with neg
    return flux_rs, ivar_rs 

def compute_precision_score_mapping(thresh, prec, score):
  ind = np.argsort(thresh);
  thresh = thresh[ind];
  prec = prec[ind];
  for i in xrange(1, len(prec)):
    prec[i] = max(prec[i], prec[i-1]);
  
  indexes = np.unique(thresh, return_index=True)[1]
  indexes = np.sort(indexes);
  thresh = thresh[indexes]
  prec = prec[indexes]
  
  thresh = np.vstack((min(-1000, min(thresh)-1), thresh[:, np.newaxis], max(1000, max(thresh)+1)));
  prec = np.vstack((prec[0], prec[:, np.newaxis], prec[-1]));
  
  f = interp1d(thresh[:,0], prec[:,0])
  val = f(score)
  return val 

def resampling_by_interpolate(self, x):
        """Resampling base on 1st order interpolation

        Parameters
        ---------
        x : array, shape ('int(len(x) * f0rate)')
            array of wsolaed waveform

        Returns
        ---------
        wsolaed: array, shape (`len(x)`)
            Array of resampled (F0 transformed) waveform sequence

        """

        # interpolate
        wedlen = len(x)
        intpfunc = interp1d(np.arange(wedlen), x, kind=1)
        x_new = np.arange(0.0, wedlen - 1, self.f0rate)
        resampled = intpfunc(x_new)

        return resampled 

def get_tukeyQcrit(k, df, alpha=0.05):
    '''
    return critical values for Tukey's HSD (Q)

    Parameters
    ----------
    k : int in {2, ..., 10}
        number of tests
    df : int
        degrees of freedom of error term
    alpha : {0.05, 0.01}
        type 1 error, 1-confidence level



    not enough error checking for limitations
    '''
    if alpha == 0.05:
        intp = interpolate.interp1d(crows, cv005[:,k-2])
    elif alpha == 0.01:
        intp = interpolate.interp1d(crows, cv001[:,k-2])
    else:
        raise ValueError('only implemented for alpha equal to 0.01 and 0.05')
    return intp(df) 

def highres(y,kind='cubic',res=100):
    '''
    Interpolate data onto a higher resolution grid by a factor of *res*

    Args:
        y (1d array/list): signal to be interpolated
        kind (str): order of interpolation (see docs for scipy.interpolate.interp1d)
        res (int): factor to increase resolution of data via linear interpolation
    
    Returns:
        shift (float): offset between target and reference signal 
    '''
    y = np.array(y)
    x = np.arange(0, y.shape[0])
    f = interp1d(x, y,kind='cubic')
    xnew = np.linspace(0, x.shape[0]-1, x.shape[0]*res)
    ynew = f(xnew)
    return xnew,ynew 

def _resample_regressor(hr_regressor, hr_frame_times, frame_times):
    """this function sub-samples the regressors at frame times

    Parameters
    ----------
    hr_regressor : array of shape(n_samples),
        the regressor time course sampled at high temporal resolution
    hr_frame_times : array of shape(n_samples),
        the corresponding time stamps
    frame_times: array of shape(n_scans),
         the desired time stamps

    Returns
    -------
    regressor: array of shape(n_scans)
         the resampled regressor
    """
    from scipy.interpolate import interp1d
    f = interp1d(hr_frame_times, hr_regressor)
    return f(frame_times).T 

def _build_f(self, X, y):
        """Build the f_ interp1d function."""

        # Handle the out_of_bounds argument by setting bounds_error
        if self.out_of_bounds not in ["raise", "nan", "clip"]:
            raise ValueError("The argument ``out_of_bounds`` must be in "
                             "'nan', 'clip', 'raise'; got {0}"
                             .format(self.out_of_bounds))

        bounds_error = self.out_of_bounds == "raise"
        if len(y) == 1:
            # single y, constant prediction
            self.f_ = lambda x: y.repeat(x.shape)
        else:
            self.f_ = interpolate.interp1d(X, y, kind='linear',
                                           bounds_error=bounds_error) 

def __init__(self, x, side='right'):
        step = True
        if step: #TODO: make this an arg and have a linear interpolation option?
            x = np.array(x, copy=True)
            x.sort()
            nobs = len(x)
            y = np.linspace(1./nobs,1,nobs)
            super(ECDF, self).__init__(x, y, side=side, sorted=True)
        else:
            return interp1d(x,y,drop_errors=False,fill_values=ival) 

def test_zero(self):
        """Check the actual implementation of zero-order spline interpolation.
        """
        interp10 = interp1d(self.x10, self.y10, kind='zero')
        assert_array_almost_equal(interp10(self.x10), self.y10)
        assert_array_almost_equal(interp10(1.2), np.array(1.))
        assert_array_almost_equal(interp10([2.4, 5.6, 6.0]),
                                  np.array([2., 6., 6.])) 

def get_log_spectrum(segment, fftlen, limit, logf, interp_logf):
    spectra = fft(segment, fftlen)
    # "fftlen is always even here."
    amplitude = np.abs(spectra[0:fftlen//2+1])
    # "ignore the zero frequency component"
    amplitude = amplitude[1:limit+2]
    interp_amplitude = interp1d(logf, amplitude)(interp_logf)
    interp_amplitude = interp_amplitude - min(interp_amplitude)
    return interp_amplitude


# ---- ComputeSHR ----- 

def _nd_check_shape(self, kind='linear'):
        # Check large ndim output shape
        a = [4, 5, 6, 7]
        y = np.arange(np.prod(a)).reshape(*a)
        for n, s in enumerate(a):
            x = np.arange(s)
            z = interp1d(x, y, axis=n, kind=kind)
            assert_array_almost_equal(z(x), y, err_msg=kind)

            x2 = np.arange(2*3*1).reshape((2,3,1)) / 12.
            b = list(a)
            b[n:n+1] = [2,3,1]
            assert_array_almost_equal(z(x2).shape, b, err_msg=kind) 

def invert_cdf(y, x):
    """
    Invert cumulative distribution function of the probability distribution

    Example:
    --------
    # analytical formula for the beam distribution
    f = lambda x: A * np.exp(-(x - mu) ** 2 / (2. * sigma ** 2))

    # we are interesting in range from -30 to 30 e.g. [um]
    x = np.linspace(-30, 30, num=100)

    # Inverted cumulative distribution function
    i_cdf = invert_cdf(y=f(x), x=x)

    # get beam distribution (200 000 coordinates)
    tau = i_cdf(np.random.rand(200000))


    :param y: array, [y0, y1, y2, ... yn] yi = y(xi)
    :param x: array, [x0, x1, x2, ... xn] xi
    :return: function
    """

    cum_int = integrate.cumtrapz(y, x, initial=0)
    #print("invert", np.max(cum_int), np.min(cum_int))
    cdf = cum_int / (np.max(cum_int) - np.min(cum_int))
    inv_cdf = interpolate.interp1d(cdf, x)
    return inv_cdf 

def test_circular_refs(self):
        # Test interp1d can be automatically garbage collected
        x = np.linspace(0, 1)
        y = np.linspace(0, 1)
        # Confirm interp can be released from memory after use
        with assert_deallocated(interp1d, x, y) as interp:
            new_y = interp([0.1, 0.2])
            del interp 

def interpolate_data(input_data, max_width):
    """
    Resample `input_data` to number of `max_width` counts
    """

    input_data = list(input_data)
    input_data_len = len(input_data)
    x = list(range(input_data_len))
    y = input_data
    xvals = np.linspace(0, input_data_len-1, max_width)
    yinterp = interp1d(x, y, kind='cubic')
    return yinterp(xvals) 

def multinterp(x, y, xquery, slow=False):
    """Multiple linear interpolations

    Parameters
    ----------
    x : array_like, shape=(N,)
        sorted array of x values
    y : array_like, shape=(N, M)
        array of y values corresponding to each x value
    xquery : array_like, shape=(M,)
        array of query values
    slow : boolean, default=False
        if True, use slow method (used mainly for unit testing)

    Returns
    -------
    yquery : ndarray, shape=(M,)
        The interpolated values corresponding to each x query.
    """
    x, y, xquery = map(np.asarray, (x, y, xquery))
    assert x.ndim == 1
    assert xquery.ndim == 1
    assert y.shape == x.shape + xquery.shape

    # make sure xmin < xquery < xmax in all cases
    xquery = np.clip(xquery, x.min(), x.max())

    if slow:
        from scipy.interpolate import interp1d
        return np.array([interp1d(x, y)(xq) for xq, y in zip(xquery, y.T)])
    elif len(x) == 3:
        # Most common case: use a faster approach
        yq_lower = y[0] + (xquery - x[0]) * (y[1] - y[0]) / (x[1] - x[0])
        yq_upper = y[1] + (xquery - x[1]) * (y[2] - y[1]) / (x[2] - x[1])
        return np.where(xquery < x[1], yq_lower, yq_upper)
    else:
        i = np.clip(np.searchsorted(x, xquery, side='right') - 1,
                    0, len(x) - 2)
        j = np.arange(len(xquery))
        return y[i, j] + ((xquery - x[i]) *
                          (y[i + 1, j] - y[i, j]) / (x[i + 1] - x[i])) 

def __call__(self, clips):
        if isinstance(clips, np.ndarray):
            H, W, _ = clips.shape
            # handle numpy array
            clips = clips.reshape((H,W,-1,self.dim)).transpose((3, 2, 0, 1))
            if self.modality == 'flow+mp4':
                if self._flow_ds_factor is not 0 or 1:
                    clips = np.transpose(clips, (1,0,2,3))
                    # downsample to make OF blocky
                    factor = self._flow_ds_factor
                    w_max = H
                    h_max = W
                    input_flow = block_reduce(clips[:,0:2, :, :], block_size=(1, 1, factor, factor), func=np.mean)
                    # resize to original size by repeating or interpolation
                    if self._upsample_interp is False:
                        input_flow = input_flow.repeat(factor, axis=2).repeat(factor, axis=3)
                    else:
                        # interpolate along certain dimension? only interp1d can do so
                        w_max_ds = input_flow.shape[2]
                        h_max_ds = input_flow.shape[3]
                        f_out = interpolate.interp1d(np.linspace(0, 1, w_max_ds), input_flow, kind='linear', axis=2)
                        input_flow = f_out(np.linspace(0, 1, w_max_ds * factor))
                        f_out = interpolate.interp1d(np.linspace(0, 1, h_max_ds), input_flow, kind='linear', axis=3)
                        input_flow = f_out(np.linspace(0, 1, h_max_ds * factor))
                    clips[:,0:2, :, :] = input_flow[:, :, :w_max, :h_max]
                clips = np.transpose(clips, (1,0,2,3))
                
            clips = torch.from_numpy(clips)
            #print(clips.shape)
            # backward compatibility
            return clips.float() / 255.0 

def _determine_particle_state_2D(self):

        known_xy = None
        known_z = None
        xs = None
        ys = None
        fact = dimensionalise(1.0, u.meter).magnitude
        if rank == 0:
            # points that we have known elevation for
            known_xy = self.badlands_model.recGrid.tinMesh['vertices'] / fact
            # elevation for those points
            known_z = self.badlands_model.elevation / fact
            xs = self.badlands_model.recGrid.regX / fact
            ys = self.badlands_model.recGrid.regY / fact

        known_xy = comm.bcast(known_xy, root=0)
        known_z = comm.bcast(known_z, root=0)
        xs = comm.bcast(xs, root=0)
        ys = comm.bcast(ys, root=0)

        comm.Barrier()

        grid_x, grid_y = np.meshgrid(xs, ys)
        interpolate_z = griddata(known_xy,
                                 known_z,
                                 (grid_x, grid_y),
                                 method='nearest').T
        interpolate_z = interpolate_z.mean(axis=1)

        f = interp1d(xs, interpolate_z)

        uw_surface = self.Model.swarm.particleCoordinates.data
        bdl_surface = f(uw_surface[:, 0])

        flags = uw_surface[:, 1] < bdl_surface

        return flags 

def extract_profile(field,
                    line,
                    nsamples=1000):
    """ Extract values along a line

    Parameters:
        field: The field to extract the data from
        line: list of (x,y, [z]) coordinates defining the sampling
              line.
        nsamples: number of sampling points
    """

    if size > 1:
        raise NotImplementedError("""The extract_profile function will not work
                                  in parallel""")

    coords = np.array([(nd(x), nd(y)) for (x, y) in line])

    x = np.linspace(coords[0, 0], coords[-1, 0], nsamples)

    if coords[0, 0] == coords[-1, 0]:
        y = np.linspace(coords[0, 1], coords[-1, 1], nsamples)
    else:
        from scipy import interpolate
        f = interpolate.interp1d(coords[:, 0], coords[:, 1])
        y = f(x)

    dx = np.diff(x)
    dy = np.diff(y)
    distances = np.zeros(x.shape)
    distances[1:] = np.sqrt(dx**2 + dy**2)
    distances = np.cumsum(distances)

    pts = np.array([x, y]).T
    values = field.evaluate(pts)

    return distances, values 

def createColorGradient(pal):
    inds = pal[0::2]
    cols = np.array(pal[1::2])
    red, grn, blu = cols >> 16, cols >> 8, cols
    cols = [c & 0xff for c in [red, grn, blu]]
    cols = [interpolate.interp1d(inds, c) for c in cols]
    img = np.arange(0, 255)
    img = np.dstack([fn(img) for fn in cols])
    return snowy.resize(img, 256, 32) 

def get_concentration_functions(composition_table_dict):

    meta = composition_table_dict["meta"]
    composition_table = Table.from_dict(composition_table_dict["data"])
    elements = [col for col in composition_table.columns if col not in meta]
    x = composition_table["X"].values
    y = composition_table["Y"].values
    cats = composition_table["X"].unique()
    concentration, conc, d, y_c, functions = {}, {}, {}, {}, RecursiveDict()

    for el in elements:
        concentration[el] = to_numeric(composition_table[el].values) / 100.0
        conc[el], d[el], y_c[el] = {}, {}, {}

        if meta["X"] == "category":
            for i in cats:
                k = "{:06.2f}".format(float(i))
                y_c[el][k] = to_numeric(y[where(x == i)])
                conc[el][k] = to_numeric(concentration[el][where(x == i)])
                d[el][k] = interp1d(y_c[el][k], conc[el][k])

            functions[el] = lambda a, b, el=el: d[el][a](b)

        else:
            functions[el] = interp2d(float(x), float(y), concentration[el])

    return functions 

def applyColorGradient(elevation_image, gradient_image):
    xvals = np.arange(256)
    yvals = gradient_image[0]
    apply_lut = interpolate.interp1d(xvals, yvals, axis=0)
    return apply_lut(snowy.unshape(np.clip(elevation_image, 0, 255))) 

def baseline_cumulative_hazard_function(self, params):
        """
        Returns a function that calculates the baseline cumulative
        hazard function for each stratum.

        Parameters
        ----------
        params : ndarray
            The model parameters.

        Returns
        -------
        A dict mapping stratum names to the estimated baseline
        cumulative hazard function.
        """

        from scipy.interpolate import interp1d
        surv = self.surv
        base = self.baseline_cumulative_hazard(params)

        cumhaz_f = {}
        for stx in range(surv.nstrat):
            time_h = base[stx][0]
            cumhaz = base[stx][1]
            time_h = np.r_[-np.inf, time_h, np.inf]
            cumhaz = np.r_[cumhaz[0], cumhaz, cumhaz[-1]]
            func = interp1d(time_h, cumhaz, kind='zero')
            cumhaz_f[self.surv.stratum_names[stx]] = func

        return cumhaz_f